Portierung von pyMC2 Bayes'sches A/B-Testbeispiel auf pyMC3

Question

Ich arbeite daran, pyMC 3 zu lernen und habe Probleme. Da es begrenzte Tutorials für pyMC3 gibt, arbeite ich von Bayesian Methods for Hackers. Ich versuche, den pyMC 2 Code zu pyMC 3 im Bayesian A/B testing Beispiel zu portieren, ohne Erfolg. Soweit ich sehen kann, berücksichtigt das Modell die Beobachtungen überhaupt nicht.

Ich habe ein paar Änderungen aus dem Beispiel machen hatte, wie pyMC 3 ist ganz anders, so was sollte wie folgt aussehen: Import pymc als 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() 

stattdessen sieht aus wie

:

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() 

Entschuldigung für die lange Post, aber in meiner Anpassung gab es eine Reihe von kleinen Änderungen, z manuelles Erzeugen der Beobachtungen, da pm.rbernoulli nicht mehr existiert. Ich bin mir auch nicht sicher, ob ich den Start vor dem Ausführen der Ablaufverfolgung finden sollte. Wie sollte ich meine Implementierung ändern, damit sie korrekt ausgeführt wird?

Answer 1

Sie waren tatsächlich in der Nähe. Allerdings ist diese Zeile:

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

ist falsch, wie Sie gerade einen konstanten Wert für p (p_true == 0,05) einstellen. Daher ist Ihre oben definierte Zufallsvariable p, um ein einheitliches Prior zu haben, nicht durch die Wahrscheinlichkeit beschränkt, und Ihr Diagramm zeigt, dass Sie nur vom Vorherigen abtasten. Wenn Sie in Ihrem Code p_true durch p ersetzen, sollte es funktionieren. Hier ist die feste Version:

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); 

Answer 2

Sie waren sehr nah dran - Sie müssen nur die letzten zwei Zeilen, die das Trace Plot erzeugen, entfernen. Sie können sich das Trace Plot als eine Diagnose vorstellen, die nach dem Ende der Probenahme erfolgen soll. Folgendes funktioniert für mich:

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() 

Answer 3

Das funktionierte für mich. Ich habe die Beobachtungen generiert, bevor ich das Modell einleitete.

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)