import pymc as pm
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [12, 8]
# conversion rate
p = pm.Uniform('p', lower=0, upper=1)
p_true = 0.05
N = 1500
occurences = pm.rbernoulli(p_true, size=N)
print(occurences.sum())
74
obs = pm.Bernoulli('obs', p, value=occurences, observed=True)
mcmc = pm.MCMC([obs, p])
mcmc.sample(20000, 10000)
[-----------------100%-----------------] 20000 of 20000 complete in 4.1 sec
plt.vlines(p_true, 0, 170, linestyles="--",
label='true $p_A$ (unknown)')
plt.hist(mcmc.trace('p')[:],
bins=35, histtype="stepfilled", normed=True)
plt.title('Posterior distribution of $p_A$, '
'the true effectiveness of site A')
plt.xlabel('Value of $p_A$')
plt.ylabel('Density')
plt.legend()
plt.show()
p_true_a = 0.05
p_true_b = 0.04
n_a = 1500
n_b = 750
observations_a = pm.rbernoulli(p_true_a, n_a)
observations_b = pm.rbernoulli(p_true_b, n_b)
p_a = pm.Uniform('p_a', lower=0, upper=1)
p_b = pm.Uniform('p_b', lower=0, upper=1)
rate_conv_a = pm.Bernoulli('obs_a', p_a, value=observations_a, observed=True)
rate_conv_b = pm.Bernoulli('obs_b', p_b, value=observations_b, observed=True)
# delta = pm.Bernoulli('delta', p_a - p_b)
@pm.deterministic
def delta(p_a=p_a, p_b=p_b):
return p_a - p_b
mcmc = pm.MCMC([p_a, p_b, rate_conv_a, rate_conv_b, delta])
mcmc.sample(40000, 10000)
[-----------------100%-----------------] 40000 of 40000 complete in 6.8 sec
plt.rcParams['figure.figsize']=[12.5, 10]
ax = plt.subplot(311)
plt.vlines(p_true_a, 0, 80, linestyles='--', label='Value of $p_A$ (unknown)')
plt.hist(mcmc.trace('p_a')[:], bins=35,
histtype='stepfilled', normed=True)
plt.xlim([0, .1])
plt.legend(loc='upper right')
# plt.subtitle('Posterior distribution of $p_A$')
ax = plt.subplot(312)
plt.vlines(p_true_b, 0, 80, linestyles='--', label='Value of $p_B$ (unknown)')
plt.hist(mcmc.trace('p_b')[:], bins=35,
histtype='stepfilled', normed=True)
plt.xlim([0, .1])
plt.legend(loc='upper right')
# plt.subtitle('Posterior distribution of $p_A$')
ax = plt.subplot(313)
plt.vlines(0.01, 0, 80, linestyles='--', label='Value of $\Delta$ (unknown)')
plt.hist(mcmc.trace('delta')[:], bins=35,
histtype='stepfilled', normed=True)
# plt.subtitle('Posterior distribution of $p_A$')
plt.legend(loc='upper right')
plt.show()
plt.hist(mcmc.trace('p_a')[:], bins=35,
histtype='stepfilled', normed=True,alpha=.35)
plt.hist(mcmc.trace('p_b')[:], bins=35,
histtype='stepfilled', normed=True,alpha=.35)
plt.xlim([0, .1])
plt.ylim([0, 80])
plt.xlabel('value')
plt.ylabel('density')
plt.legend(loc='upper right')
plt.title('Posterior distribution of $p_A$ and $p_B$')
plt.show()
No handles with labels found to put in legend.