用新的观测数据更新PyMC3模型

Question

去年我测量了80种水果的直径，在检查了这些值的最佳分布之后，我创建了一个PyMC3模型

with Model() as diam_model:
    mu = Normal('mu',mu=57,sd=5.42)
    sigma = Uniform('sigma',0,10)

之后，据我所知，我已经用我以前的数据（80个值）对模型进行了“训练”

^{pr2}$

然后我使用了samples的plot_posterior，还返回了平均值和HPD。在

我的想法是今年再次使用贝叶斯更新来减少样本量。如何添加单个值，并更新后面的值，期望HPD越来越小？在

Answer 1

核密度估计更新的先验知识

使用另一个建议的重复答案，可以使用this 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)

从后验中提取新的先验信息

然后，我们可以使用此模型的结果从the referenced notebook提取参数的KDE后验点，代码如下：

^{pr2}$

第二轮

现在，如果我们有更多的数据，我们可以使用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基近似在多大程度上引入了错误，以及当重复使用时它们如何在模型中传播。这是一个很好的技巧，但是在没有进一步验证其健壮性的情况下，应该谨慎地将其投入生产。在

但是，在您的例子中，期望的分布是高斯分布，并且这些分布有established closed-form conjugate models。我强烈建议您完成Kevin Murphy's Conjugate Bayesian analysis of the Gaussian distribution。在

正态反γ模型

正态反伽马模型估计观测到的正态随机变量的均值和方差。平均值用正态先验建模；方差用反伽玛函数建模。此模型使用了四个参数：

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个观测值不足以完全移动到我提供的新位置和比例，但两个参数似乎都朝着这个方向移动。在