python - 使用新观察到的数据更新 PyMC3 上的模型
问题描述
去年我测量了 80 个水果的直径,在检查了值的最佳分布后,我创建了一个 PyMC3 模型
with Model() as diam_model:
mu = Normal('mu',mu=57,sd=5.42)
sigma = Uniform('sigma',0,10)
之后,据我了解,我已经用我之前的数据(80 个值)“训练”了模型
with diam_model:
dist = Normal('dist',mu=mu,sd=sigma, observed=prior_data.values)
with diam_model:
samples=fit().sample(1000)
然后我使用 的plot_posterior
,samples
还返回均值和 HPD。
我的想法是今年再次使用贝叶斯更新来测量以减少样本量。如何添加单个值并更新后验,期望 HPD 变得越来越小?
解决方案
核密度估计更新先验
使用建议的另一个答案作为副本,可以使用此 Jupyter notebook中的代码提取先验的近似版本。
第一回合
我假设我们有来自第一轮抽样的数据,我们可以将平均值设为 57.0,标准差设为 5.42。
import numpy as np
import pymc3 as pm
from sklearn.preprocessing import scale
from scipy import stats
# generate data forced to match distribution indicated
Y0 = 57.0 + scale(np.random.normal(size=80))*5.42
with pm.Model() as m0:
# let's place an informed, but broad prior on the mean
mu = pm.Normal('mu', mu=50, sd=10)
sigma = pm.Uniform('sigma', 0, 10)
y = pm.Normal('y', mu=mu, sd=sigma, observed=Y0)
trace0 = pm.sample(5000, tune=5000)
从后验中提取新的先验
然后,我们可以使用该模型的结果从参考笔记本中使用以下代码提取参数的 KDE 后验:
def from_posterior(param, samples, k=100):
smin, smax = np.min(samples), np.max(samples)
width = smax - smin
x = np.linspace(smin, smax, k)
y = stats.gaussian_kde(samples)(x)
# what was never sampled should have a small probability but not 0,
# so we'll extend the domain and use linear approximation of density on it
x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
y = np.concatenate([[0], y, [0]])
return pm.Interpolated(param, x, y)
第二轮
现在,如果我们有更多数据,我们可以使用 KDE 更新的先验来运行一个新模型:
Y1 = np.random.normal(loc=57, scale=5.42, size=100)
with pm.Model() as m1:
mu = from_posterior('mu', trace0['mu'])
sigma = from_posterior('sigma', trace0['sigma'])
y = pm.Normal('y', mu=mu, sd=sigma, observed=Y1)
trace1 = pm.sample(5000, tune=5000)
同样,人们可以使用此跟踪来提取更新的后验估计,以用于未来的输入数据轮次。
买者自负
上述方法产生了对真实更新先验的近似值,并且在共轭先验不可能的情况下最有用。还应该注意的是,我不确定这种基于 KDE 的近似值在多大程度上引入了错误,以及它们在重复使用时如何在模型中传播。这是一个巧妙的技巧,但在未进一步验证其稳健性的情况下将其投入生产应谨慎。
特别是,我会非常关注后验分布具有强相关结构的情况。此处提供的代码仅使用每个潜在变量的边缘生成“先验”分布。这对于像这样的简单模型来说似乎很好,而且不可否认,初始先验也缺乏相关性,所以这可能不是一个大问题。然而,一般来说,对边缘进行总结涉及丢弃有关变量如何相关的信息,而在其他情况下,这可能相当重要。例如,Beta 分布的默认参数化总是导致后验中的相关参数,因此上述技术是不合适的。相反,需要推断出包含所有潜在变量的多维 KDE。
共轭模型
但是,在您的特定情况下,预期分布是高斯分布,并且这些分布已经建立了封闭形式的共轭模型,即原则解决方案而不是近似值。我强烈建议通过Kevin Murphy 的高斯分布共轭贝叶斯分析。
正反伽马模型
Normal-Inverse Gamma 模型估计观察到的正态随机变量的均值和方差。均值是用正态先验建模的;具有反伽马的方差。该模型使用四个参数作为先验:
mu_0 = prior mean
nu = number of observations used to estimate the mean
alpha = half the number of obs used to estimate variance
beta = half the sum of squared deviations
给定您的初始模型,我们可以使用这些值
mu_0 = 57.0
nu = 80
alpha = 40
beta = alpha*5.42**2
然后,您可以绘制先验的对数似然,如下所示:
# points to compute likelihood at
mu_grid, sd_grid = np.meshgrid(np.linspace(47, 67, 101),
np.linspace(4, 8, 101))
# normal ~ N(X | mu_0, sigma/sqrt(nu))
logN = stats.norm.logpdf(x=mu_grid, loc=mu_0, scale=sd_grid/np.sqrt(nu))
# inv-gamma ~ IG(sigma^2 | alpha, beta)
logIG = stats.invgamma.logpdf(x=sd_grid**2, a=alpha, scale=beta)
# full log-likelihood
logNIG = logN + logIG
# actually, we'll plot the -log(-log(likelihood)) to get nicer contour
plt.figure(figsize=(8,8))
plt.contourf(mu_grid, sd_grid, -np.log(-logNIG))
plt.xlabel("$\mu$")
plt.ylabel("$\sigma$")
plt.show()
更新参数
给定新数据,Y1
,更新参数如下:
# precompute some helpful values
n = Y1.shape[0]
mu_y = Y1.mean()
# updated NIG parameters
mu_n = (nu*mu_0 + n*mu_y)/(nu + n)
nu_n = nu + n
alpha_n = alpha + n/2
beta_n = beta + 0.5*np.square(Y1 - mu_y).sum() + 0.5*(n*nu/nu_n)*(mu_y - mu_0)**2
为了说明模型的变化,让我们从稍微不同的分布中生成一些数据,然后绘制生成的后验对数似然图:
np.random.seed(53211277)
Y1 = np.random.normal(loc=62, scale=7.0, size=20)
产生
在这里,20 个观测值不足以完全移动到我提供的新位置和尺度,但两个参数似乎都朝那个方向移动。
推荐阅读
- javascript - JavaScript 检查是否在 HTML 表单中选择了多个单选输入中的任何一个
- jquery - jQuery 解析 xml 无效的 XML:HTTP/1.1 200 OK
- php - 如何在 HTML 文件中执行 PHP 代码?
- c# - 实例化对象不实例化其他对象
- sql - 如何按群组显示从一周到另一周的登录变化
- ios - UITableViewCell 被格式化为默认大小
- ruby-on-rails-4 - Rails FFMPEG - 转换为 webm 为编码器“libvpx-vp9”提供错误
- mongodb - 将记录插入嵌套数组
- c++ - 为什么我在声明本身中收到“未在范围内声明”错误?
- python - 以印度逗号样式输出的 Pandas 数据框