6

我正在努力学习 pyMC 3 并且遇到了一些麻烦。由于 pyMC3 的教程有限，我正在使用Bayesian Methods for Hackers。我试图在贝叶斯 A/B 测试示例中将 pyMC 2 代码移植到 pyMC 3，但没有成功。据我所知，该模型根本没有考虑到观察结果。

我不得不对示例进行一些更改，因为 pyMC 3 完全不同，所以应该如下所示： import pymc as 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, 0, 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', 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)
    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 不再存在。我也不确定是否应该在运行跟踪之前找到开始。我应该如何更改我的实现以正确运行？

4

3 回答 3

4

你确实很亲近。但是，这一行：

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

是错误的，因为您只是为 p 设置了一个常量值（p_true == 0.05）。因此，您在上面定义为具有统一先验的随机变量 p 不受可能性的限制，并且您的绘图表明您只是从先验中采样。如果您在代码中将 p_true 替换为 p 它应该可以工作。这是固定版本：

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', 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)
    for i in xrange(N):
        occurrences.append((random.uniform(0.0, 1.0) <= p_true))
    occurrences = np.array(occurrences)
    obs = pm.Bernoulli('obs', p, observed=occurrences)

    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(18000, step, start)

pm.traceplot(trace);
于 2014-03-30T09:25:51.163 回答
1

这对我有用。我在启动模型之前生成了观察结果。

true_p_A = 0.05
true_p_B = 0.04
N_A = 1500
N_B = 750

obs_A = np.random.binomial(1, true_p_A, size=N_A)
obs_B = np.random.binomial(1, true_p_B, size=N_B)

with pm.Model() as ab_model:
    p_A = pm.Uniform('p_A', 0, 1)
    p_B = pm.Uniform('p_B', 0, 1)
    delta = pm.Deterministic('delta',p_A - p_B)
    obs_A = pm.Bernoulli('obs_A', p_A, observed=obs_A)
    osb_B = pm.Bernoulli('obs_B', p_B, observed=obs_B)

with ab_model:
    trace = pm.sample(2000)

pm.traceplot(trace)
于 2016-08-14T15:31:35.977 回答
0

你非常接近 - 你只需要取消最后两行，它们会产生跟踪图。您可以将绘制跟踪图视为完成采样后应进行的诊断。以下对我有用：

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', 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)
    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)

#Now plot
pm.traceplot(trace)
plt.show()
于 2014-03-29T02:04:55.663 回答