python - 开始使用简单的 pymc3 示例时遇到麻烦

Question

我是 PyMC3 包的新手，我只是想实现我正在学习的测量不确定性类(class)中的一个示例。 (请注意，这是一门通过工作选修的员工教育类(class)，而不是我不应该在网上找到答案的分级类(class))。该类(class)使用 R，但我发现 python 更可取。

(简单的)问题如下:

假设您有一个在室温下的实际(未知)长度 length 和测量长度 m 的端规。两者的关系是:

length = m / (1 + alpha*dT)

其中 alpha 是膨胀系数，dT 是与室温的偏差，m 是测量量。目标是找到 length 的后验分布以确定其期望值和标准差(即测量不确定性)

该问题指定了 alpha 和 dT(标准差较小的高斯分布)的先验分布以及length(标准差较大的高斯分布)的松散先验分布。问题指定 m 测量了 25 次，平均值为 50.000215，标准差为 5.8e-6。我们假设 m 的测量值服从 m 真实值的均值正态分布。

我遇到的一个问题是，可能性似乎无法仅根据 PyMC3 中的这些统计数据来指定，因此我生成了一些虚拟测量数据(我最终进行了 1000 次测量，而不是 25 次测量)。同样，问题是在 length 上获得后验分布(在此过程中，虽然不太感兴趣，但在 alpha 和 dT 上更新了后验分布).

这是我的代码，它不起作用并且存在收敛问题:

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as stats
import pymc3 as pm
import theano.tensor as tt

basic_model = pm.Model()

xdata = np.random.normal(50.000215,5.8e-6*np.sqrt(1000),1000)

with basic_model:
    #prior distributions
    theta = pm.Normal('theta',mu=-.1,sd=.04)
    alpha = pm.Normal('alpha',mu=.0000115,sd=.0000012)
    length = pm.Normal('length',mu=50,sd=1)

mumeas = length*(1+alpha*theta)


with basic_model:
    obs = pm.Normal('obs',mu=mumeas,sd=5.8e-6,observed=xdata)
    #yobs = Normal('yobs',)
    start = pm.find_MAP()
    #trace = pm.sample(2000, step=pm.Metropolis, start=start)
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=200000,step=step,start=start,njobs=1)


length_samples = trace['length']

fig,ax=plt.subplots()
plt.hist(length_samples, histtype='stepfilled', bins=30, alpha=0.85,
     label="posterior of $\lambda_1$", color="#A60628", normed=True)

对于为什么这不起作用的任何帮助，我将不胜感激。我已经尝试了一段时间，但它从未收敛到 R 代码给出的预期解决方案。我尝试了默认采样器(我认为是 NUTS)以及 Metropolis，但由于零梯度错误而完全失败。 (相关的)类(class)幻灯片作为图像附加。最后，这是可比较的 R 代码:

library(rjags)

#Data
jags_data <- list(xbar=50.000215) 

jags_code <- jags.model(file = "calibration.txt",
                    data = jags_data, 
                    n.chains = 1, 
                    n.adapt = 30000)


post_samples <- coda.samples(model = jags_code, 
                         variable.names = 
c("l","mu","alpha","theta"),#,"ypred"),
                         n.iter = 30000)


summary(post_samples)

mean(post_samples[[1]][,"l"])
sd(post_samples[[1]][,"l"])

plot(post_samples)

和 calibration.txt 模型:

model{
  l~dnorm(50,1.0)
  alpha~dnorm(0.0000115,694444444444)
  theta~dnorm(-0.1,625)
  mu<-l*(1+alpha*theta)
  xbar~dnorm(mu,29726516052)
}

(请注意，我认为 dnorm 分布采用 1/sigma^2，因此方差看起来很奇怪)

对于 PyMC3 采样为何不收敛以及我应该采取哪些不同做法的任何帮助或见解，我们将不胜感激。谢谢!

Answer 1

我也无法从代码中生成的数据和模型中获取任何有用信息。在我看来，假数据中的噪声水平同样可以用模型中不同的方差来源来解释。这可能导致后验参数高度相关的情况。再加上极端的规模不平衡，那么这会有采样问题是有道理的。

然而，看看 JAGS 模型，他们似乎真的只使用了那个一个输入观察。我以前从未见过这种技术(？)，即输入数据的汇总统计信息而不是原始数据本身。我想它在 JAGS 中对他们有用，所以我决定尝试运行完全相同的 MCMC，包括使用高斯的精度 (tau) 参数化。

带有大都会的原始模型

with pm.Model() as m0:
    # tau === precision parameterization
    dT = pm.Normal('dT', mu=-0.1, tau=625)
    alpha = pm.Normal('alpha', mu=0.0000115, tau=694444444444)
    length = pm.Normal('length', mu=50.0, tau=1.0)

    mu = pm.Deterministic('mu', length*(1+alpha*dT))

    # only one input observation; tau indicates the 5.8 nm sd
    obs = pm.Normal('obs', mu=mu, tau=29726516052, observed=[50.000215])

    trace = pm.sample(30000, tune=30000, chains=4, cores=4, step=pm.Metropolis())

虽然它在采样 length 和 dT 方面仍然不是很好，但它至少看起来总体上是收敛的:

我认为这里值得注意的是，尽管 length (sd=1) 的先验相对较弱，但所有其他参数的强先验似乎传播了一个严格的不确定性绑定(bind)在 length 后部。最终，这是兴趣的后验，所以这似乎与练习的意图一致。另外，看到 mu 出现在后部，正如所描述的分布，即 N(50.000215, 5.8e-6)。

轨迹图

森林图

配对图

然而，在这里，您可以看到核心问题仍然存在。 length 和dT 之间存在很强的相关性，并且标准误差之间存在 4 或 5 个数量级的比例差异。在我真正相信结果之前，我肯定会跑很长时间。

带有 NUTS 的替代模型

为了让这个与 NUTS 一起运行，您必须解决缩放问题。也就是说，我们需要以某种方式重新参数化以使所有 tau 值更接近 1。然后，您将运行采样器并转换回您感兴趣的单位。不幸的是，我不知道现在没时间玩这个(我也得想办法)，但也许这是你可以开始自己探索的东西。