4

我是贝叶斯统计和 MCMC 的绝对新手，所以我正在研究 John Kruschke 的“做贝叶斯数据分析：R 和 BUGS 教程”。为了测试我的理解，我正在尝试将他的示例从 BUGS 转换为 PyMC。

英寸。在图 8 中，他观察了两个（可能有偏差的）硬币中的每一个的一系列翻转，并试图估计它们的偏差的差异。低于每个硬币的偏差是theta，观察到的翻转是y：

import numpy as np
import pymc
priorA = 3
priorB = 3
theta1 = pymc.Beta('theta1', alpha=priorA, beta=priorB)
theta2 = pymc.Beta('theta2', alpha=priorA, beta=priorB)
y1 = pymc.Bernoulli('y1', p=theta1, value=np.array([1,1,1,1,1,0,0]), observed=True)
y2 = pymc.Bernoulli('y2', p=theta2, value=np.array([1,1,0,0,0,0,0]), observed=True)

当然，这两个硬币是独立的。如果我分别模拟它们，然后查看 thetas 的差异，我会得到与本书相同的答案。（它也与网格上的分析解决方案和集成一致。）

m1 = pymc.MCMC([y1, theta1])
m1.sample(iter=10000, burn=1000, thin=2)
m2 = pymc.MCMC([y2, theta2])
m2.sample(iter=10000, burn=1000, thin=2)
dt = m1.trace('theta1')[:] - m2.trace('theta2')[:]
print np.average(dt), pymc.utils.hpd(dt, alpha=0.05)
# 0.23098692189 [-0.15303663  0.55686801]

另一方面，如果我尝试同时模拟两个硬币，作为同一模型的一部分，我会得到一个完全不同的答案，与任何其他方法都不一致。

model = pymc.MCMC([y1, theta1, y2, theta2])
model.sample(iter=10000, burn=1000, thin=2)
dt = model.trace('theta1')[:] - model.trace('theta2')[:]
print np.average(dt), pymc.utils.hpd(dt, alpha=0.05)
# 0.330061191424 [-0.09999254  0.7190196 ]

所以我的第一个问题是，为什么？根据我对该理论的了解很少，我的猜测是 MCMC 很难充分探索二参数模型。（但是，BUGS 似乎处理得很好。）

真正奇怪的是，我一直在 iPython 笔记本中做这一切，而且 PyMC 中似乎存在一个错误。如果我运行独立硬币模型，重新启动我的内核（Kernel | Restart or File | Close and Halt），然后运行联合硬币模型，联合硬币将产生与独立硬币相似（但不相同）的答案（平均dtheta ~ 0.23)。如果我以相反的顺序运行模型（在两者之间重新启动内核），它们都会从联合硬币模型中产生（不正确的）答案，平均 dtheta ~ 0.33。如果我在两者之间完全关闭我的 iPython 笔记本服务器，我只能让这两个模型产生不同的答案。由于这也会从内存中卸载所有共享库，我希望这意味着 PyMC 的 Fortran/C 部分正在将这些模型中的某些内容缓存在一个内存位置 s 在 Python 解释器实例之间共享。版本是 Numpy 1.8.2、PyMC 2.3.3、Python 2.7.8、iPython 2.1.0、Anaconda 2.0.0。

对于理解这里发生的事情，我真的很感激任何帮助。我意识到这些都是愚蠢的、微不足道的模型，但 PyMC 的奇怪行为目前并不能激发信心！

4

0 回答 0