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

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

最後の 2 行をインデント解除するだけでトレースプロットが生成されます。traceplot をプロットすることは、サンプリングの終了後に発生する診断と考えることができます。以下は私にとってはうまくいきます：

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 に答える