时间戳

首页 > 解决方案 > 开始使用简单的 pymc3 示例时遇到问题

python - 开始使用简单的 pymc3 示例时遇到问题

问题描述

我是使用 PyMC3 包的新手，我只是想从我正在学习的测量不确定性课程中实现一个示例。（请注意，这是通过工作进行的可选员工教育课程，而不是我不应该在网上找到答案的分级课程）。该课程使用 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，但完全失败，梯度误差为零。（相关）课程幻灯片作为图像附加。最后，这是可比较的 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 采样为什么不收敛以及我应该做些什么不同的帮助或见解将不胜感激。谢谢！

标签： pythonrbayesianpymc3jags

解决方案

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

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

Metropolis 原创模型

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) 上的先验相对较弱，但所有其他参数的强先验似乎在length后验上传播了严格的不确定性界限。最终，这是感兴趣的后验，因此这似乎与练习的意图一致。此外，请参阅mu后验结果，正如所描述的分布一样，即N(50.000215, 5.8e-6)。

跟踪图

在此处输入图像描述

森林图

在此处输入图像描述

配对图

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

在此处输入图像描述

带 NUTS 的替代模型

为了使用 NUTS 运行它，您必须解决扩展问题。也就是说，我们需要以某种方式重新参数化以使所有tau值更接近 1。然后，您将运行采样器并转换回您感兴趣的单位。不幸的是，我没有时间玩这个现在（我也必须弄清楚），但也许你可以自己开始探索。

推荐阅读