PyMC Bernoulli 模型检查答案

【问题标题】：PyMC Bernoulli model checkingPyMC Bernoulli 模型检查
【发布时间】：2015-11-10 00:30:20
【问题描述】：

我目前正在尝试使用 PyMC 进行模型检查，其中我的模型是伯努利模型并且我之前有一个 Beta。我想做一个 (i) gof 图以及 (ii) 计算后验预测 p 值。

我的代码使用二项式模型运行，但我很难找到使伯努利模型工作的正确方法。不幸的是，在任何地方都没有我可以使用的示例。我的代码如下所示：

import pymc as mc
import numpy as np
alpha = 2
beta = 2
n = 13
yes = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,0,0])

p = mc.Beta('p',alpha,beta)
surv = mc.Bernoulli('surv',p=p,observed=True,value=yes)
surv_sim = mc.Bernoulli('surv_sim',p=p)

mc_est = mc.MCMC({'p':p,'surv':surv,'surv_sim':surv_sim})
mc_est.sample(10000,5000,2)

import matplotlib.pylab as plt
plt.hist(mc_est.surv_sim.trace(),bins=range(0,3),normed=True)
plt.figure()
plt.hist(mc_est.p.trace(),bins=100,normed=True)

mc.Matplot.gof_plot(mc_est.surv_sim.trace(), 10/13., name='surv')

#here I have issues
D = mc.discrepancy(yes, surv_sim, p.trace())
mc.Matplot.discrepancy_plot(D)

我遇到的主要问题是确定discrepancy 函数的预期值。仅使用 p.trace() 在这里不起作用，因为这些是概率。不知何故，我需要在这里合并样本量，但我正在努力以与二项式模型类似的方式来做到这一点。我也不太确定，如果我正确地执行gof_plot。

希望有人可以在这里帮助我！谢谢！

【问题讨论】：

标签： python bayesian pymc bernoulli-probability

【解决方案1】：

根据discrepancy函数doc string，参数为：

observed : Iterable of observed values (size=(n,))
simulated : Iterable of simulated values (size=(r,n))
expected : Iterable of expected values (size=(r,) or (r,n))

所以你需要纠正两件事：

1) 修改模拟结果的大小为 n（例如，在您的示例中为 13）：

surv_sim = mc.Bernoulli('surv_sim', p=p, size=n)

2) 使用bernoulli_expval 方法封装您的p.trace()：

D = mc.discrepancy(yes, surv_sim.trace(), mc.bernoulli_expval(p.trace()))

（bernoulli_expval 只是吐回p。）

通过这两个更改，我得到以下信息：

【讨论】：

非常感谢您的回复。 size 参数的提示很棒！关于期望值，我没有使用expval 函数，因为正如你提到的，它只会吐回p。因此，我再次不确定是否使用此解决方案，差异是否在做它应该做的事情。它现在所做的是将数据和模拟的二元结果与概率参数进行比较。
另一个想法：如果我将此模型扩展到分类模型，我更不确定如何处理它。
我对@987654337@ 函数不是很熟悉。在this thread 中，pymc 的作者之一建议使用gof_plot 生成后验预测图。我不知道为什么。