使用新观察到的数据更新 PyMC3 上的模型

Posted 2023-02-19

【中文标题】使用新观察到的数据更新 PyMC3 上的模型

【英文标题】：Updating model on PyMC3 with new observed data

【发布时间】：2019-04-12 04:05:21

【问题描述】：

我去年测量了 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)

然后我使用samples 的plot_posterior，同时返回平均值和HPD。

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

【问题讨论】：

Incremental model update with PyMC3的可能重复

@merv 我试图弄清楚 y0 是否是新值

看看他们在答案中链接到的笔记本：github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/… 主要是你没有像你说的那样使用“最佳发行版”，而是为所有人提取基于 KDE 的发行版跟踪结果中的变量，然后在下一轮采样中使用这些后验分布作为新的先验。

还可能值得注意的是，如果您切换到 sd 上的 InverseGamma 先验（或 tau 上的 Gamma），那么您的模型将是共轭的，而精确的后验则有一个封闭的形式。在这种情况下，您可以使用任意数量的新观察值进行在线更新，而无需运行 MCMC。 Wikipedia actually has a nice reference table。 This CrossValidated question 也可能提供信息。

@merv 我想你的建议会让我走上新的道路。最后。毕竟，我只想找到一种简单的方法来减少我的样本量——当然还有时间。时间就是金钱。

【参考方案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)

从后验中提取新的先验

然后我们可以使用这个模型的结果来提取参数上的 KDE 后验，使用来自the referenced notebook 的以下代码：

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。

---

共轭模型

但是，在您的特定情况下，预期分布是高斯分布，并且这些分布具有established closed-form conjugate models，即有原则的解决方案而不是近似值。我强烈建议您使用Kevin Murphy's Conjugate Bayesian analysis of the Gaussian distribution。

正反伽马模型

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