python – 将pyMC2贝叶斯A / B测试示例移植到pyMC3

我正在努力学习pyMC 3并遇到一些麻烦.由于pyMC3的教程有限,我从 Bayesian Methods for Hackers开始工作.我试图在 Bayesian A/B testing示例中将pyMC 2代码移植到pyMC 3,但没有成功.从我所看到的,模型没有考虑到观察结果.

我必须从示例中做一些更改,因为pyMC 3是完全不同的,所以应该是这样的：
导入pymc为pm

# The parameters are the bounds of the Uniform.
p = pm.Uniform('p',lower=0,upper=1)

# set constants
p_true = 0.05  # remember,this is unkNown.
N = 1500

# sample N Bernoulli random variables from Ber(0.05).
# each random variable has a 0.05 chance of being a 1.
# this is the data-generation step
occurrences = pm.rbernoulli(p_true,N)

print occurrences  # Remember: Python treats True == 1,and False == 0
print occurrences.sum()

# Occurrences.mean is equal to n/N.
print "What is the observed frequency in Group A? %.4f" % occurrences.mean()
print "Does this equal the true frequency? %s" % (occurrences.mean() == p_true)

# include the observations,which are Bernoulli
obs = pm.Bernoulli("obs",p,value=occurrences,observed=True)

# To be explained in chapter 3
mcmc = pm.MCMC([p,obs])
mcmc.sample(18000,1000)

figsize(12.5,4)
plt.title("Posterior distribution of $p_A$,the true effectiveness of site A")
plt.vlines(p_true,90,linestyle="--",label="true $p_A$(unkNown)")
plt.hist(mcmc.trace("p")[:],bins=25,histtype="stepfilled",normed=True)
plt.legend()

相反看起来像：

import pymc as pm

import random
import numpy as np
import matplotlib.pyplot as plt

with pm.Model() as model:
    # Prior is uniform: all cases are equally likely
    p = pm.Uniform('p',upper=1)

    # set constants
    p_true = 0.05  # remember,this is unkNown.
    N = 1500

    # sample N Bernoulli random variables from Ber(0.05).
    # each random variable has a 0.05 chance of being a 1.
    # this is the data-generation step
    occurrences = []  # pm.rbernoulli(p_true,N)
    for i in xrange(N):
        occurrences.append((random.uniform(0.0,1.0) <= p_true))
    occurrences = np.array(occurrences)
    obs = pm.Bernoulli('obs',p_true,observed=occurrences)

    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(18000,step,start)
    pm.traceplot(trace);
    plt.show()

对于冗长的帖子表示道歉,但在我的改编中,有一些小的变化,例如：手动生成观察结果,因为pm.rbernoulli不再存在.我也不确定在运行跟踪之前是否应该找到开始.我该如何更改我的实现以正确运行？

解决方法

你确实很亲密.但是,这一行：

obs = pm.Bernoulli('obs',observed=occurrences)

这是错误的,因为你只是为p设置一个常量值(p_true == 0.05).因此,上面定义的具有统一先验的随机变量p不受可能性的约束,并且您的图表显示您只是从先前采样.如果在代码中用p替换p_true,它应该可以工作.这是固定版本：

import pymc as pm

import random
import numpy as np
import matplotlib.pyplot as plt

with pm.Model() as model:
    # Prior is uniform: all cases are equally likely
    p = pm.Uniform('p',start)

pm.traceplot(trace);