开始使用简单的 pymc3 示例时遇到麻烦
trouble getting started with simple pymc3 example
我是 PyMC3 包的新手,我只是想从我正在学习的测量不确定性课程中实施一个例子。 (请注意,这是一门通过工作选修的员工教育课程,而不是分级 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,但由于零梯度错误而完全失败。 (相关的)课程幻灯片作为图像附加。最后,这是可比较的 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 采样为何不收敛以及我应该采取哪些不同做法的任何帮助或见解,我们将不胜感激。谢谢!
我也无法从代码中生成的数据和模型中获取有用的信息。在我看来,假数据中的噪声水平同样可以用模型中不同的方差来源来解释。这可能导致后验参数高度相关的情况。再加上极端的规模不平衡,那么这会有抽样问题是有道理的。
然而,看看 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
) 的先验相对较弱,但所有其他参数的强先验似乎在 [=12] 上传播了一个严格的不确定性界限=]后。最终,这是兴趣的后验,所以这似乎与练习的意图一致。另外,看到 mu
出现在后验中,正如所描述的分布,即 N(50.000215, 5.8e-6)
.
跟踪图
森林图
对图
然而,在这里,您可以看到核心问题仍然存在。 length
和 dT
之间存在很强的相关性,再加上标准误差之间有 4 或 5 个数量级的比例差异。在我真正相信结果之前,我肯定会做很长时间 运行。
带有 NUTS 的替代模型
为了获得这个 运行 NUTS,您必须解决缩放问题。也就是说,我们需要以某种方式重新参数化以使所有 tau
值更接近 1。然后,您将 运行 采样器并转换回您感兴趣的单位。不幸的是,我不知道现在有时间玩这个(我也必须弄清楚),但也许这是你可以开始自己探索的东西。
我是 PyMC3 包的新手,我只是想从我正在学习的测量不确定性课程中实施一个例子。 (请注意,这是一门通过工作选修的员工教育课程,而不是分级 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,但由于零梯度错误而完全失败。 (相关的)课程幻灯片作为图像附加。最后,这是可比较的 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 采样为何不收敛以及我应该采取哪些不同做法的任何帮助或见解,我们将不胜感激。谢谢!
我也无法从代码中生成的数据和模型中获取有用的信息。在我看来,假数据中的噪声水平同样可以用模型中不同的方差来源来解释。这可能导致后验参数高度相关的情况。再加上极端的规模不平衡,那么这会有抽样问题是有道理的。
然而,看看 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
) 的先验相对较弱,但所有其他参数的强先验似乎在 [=12] 上传播了一个严格的不确定性界限=]后。最终,这是兴趣的后验,所以这似乎与练习的意图一致。另外,看到 mu
出现在后验中,正如所描述的分布,即 N(50.000215, 5.8e-6)
.
跟踪图
森林图
对图
然而,在这里,您可以看到核心问题仍然存在。 length
和 dT
之间存在很强的相关性,再加上标准误差之间有 4 或 5 个数量级的比例差异。在我真正相信结果之前,我肯定会做很长时间 运行。
带有 NUTS 的替代模型
为了获得这个 运行 NUTS,您必须解决缩放问题。也就是说,我们需要以某种方式重新参数化以使所有 tau
值更接近 1。然后,您将 运行 采样器并转换回您感兴趣的单位。不幸的是,我不知道现在有时间玩这个(我也必须弄清楚),但也许这是你可以开始自己探索的东西。