使用新的观察数据更新 PyMC3 上的模型
Updating model on PyMC3 with new observed data
我去年测量了 80 个水果的直径，在检查了值的最佳分布后，我创建了一个 PyMC3 模型
with Model() as diam_model:
mu = Normal('mu',mu=57,sd=5.42)
sigma = Uniform('sigma',0,10)
之后，据我所知，我 "trained" 模型与我的先前数据（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 越来越小？
内核密度估计更新先验
使用建议的另一个答案作为副本，可以使用 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-base 近似值在多大程度上引入了错误，以及它们在重复使用时如何在模型中传播。这是一个巧妙的技巧，但在没有进一步验证其稳健性的情况下将其投入生产时应该谨慎。
特别是，我会非常关注后验分布具有强相关结构的情况。此处提供的代码仅使用每个潜在变量的边缘生成一个“先验”分布。这对于像这样的简单模型来说似乎很好，而且不可否认，初始先验也缺乏相关性，所以这在这里可能不是一个大问题。然而，一般来说，对边缘进行汇总涉及丢弃有关变量如何相关的信息，而在其他情况下，这可能相当重要。例如，Beta 分布的默认参数化总是导致后验相关参数，因此上述技术是不合适的。相反，需要推断一个包含所有潜在变量的 multi-dimensional KDE。
共轭模型
但是，在您的特定情况下，预期分布是高斯分布，并且这些分布具有 established closed-form conjugate models, i.e., principled solutions rather than approximations. I strongly recommend working through Kevin Murphy's Conjugate Bayesian analysis of the Gaussian distribution.
Normal-Inverse伽玛模型
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
然后您可以绘制先验的 log-likelihood，如下所示：
# 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
为了说明模型的变化，让我们从稍微不同的分布中生成一些数据，然后绘制生成的后验结果 log-likelihood:
np.random.seed(53211277)
Y1 = np.random.normal(loc=62, scale=7.0, size=20)
产生
在这里，20 个观测值不足以完全移动到我提供的新位置和比例，但两个参数似乎都朝那个方向移动了。
我去年测量了 80 个水果的直径，在检查了值的最佳分布后，我创建了一个 PyMC3 模型
with Model() as diam_model:
mu = Normal('mu',mu=57,sd=5.42)
sigma = Uniform('sigma',0,10)
之后，据我所知，我 "trained" 模型与我的先前数据（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 越来越小？
内核密度估计更新先验
使用建议的另一个答案作为副本，可以使用 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-base 近似值在多大程度上引入了错误，以及它们在重复使用时如何在模型中传播。这是一个巧妙的技巧，但在没有进一步验证其稳健性的情况下将其投入生产时应该谨慎。
特别是，我会非常关注后验分布具有强相关结构的情况。此处提供的代码仅使用每个潜在变量的边缘生成一个“先验”分布。这对于像这样的简单模型来说似乎很好，而且不可否认，初始先验也缺乏相关性，所以这在这里可能不是一个大问题。然而，一般来说，对边缘进行汇总涉及丢弃有关变量如何相关的信息，而在其他情况下，这可能相当重要。例如，Beta 分布的默认参数化总是导致后验相关参数，因此上述技术是不合适的。相反，需要推断一个包含所有潜在变量的 multi-dimensional KDE。
共轭模型
但是，在您的特定情况下，预期分布是高斯分布，并且这些分布具有 established closed-form conjugate models, i.e., principled solutions rather than approximations. I strongly recommend working through Kevin Murphy's Conjugate Bayesian analysis of the Gaussian distribution.
Normal-Inverse伽玛模型
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
然后您可以绘制先验的 log-likelihood，如下所示：
# 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
为了说明模型的变化，让我们从稍微不同的分布中生成一些数据，然后绘制生成的后验结果 log-likelihood:
np.random.seed(53211277)
Y1 = np.random.normal(loc=62, scale=7.0, size=20)
产生
在这里，20 个观测值不足以完全移动到我提供的新位置和比例，但两个参数似乎都朝那个方向移动了。