import pymc3 as pm
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%qtconsole --colors=linux
plt.style.use('ggplot')
from matplotlib import gridspec
from theano import tensor as tt
from scipy import stats
# pledger data:
s1 = 424
s2 = 5416
n1 = 777
n2 = 9072
# two-sided p-value = 0.005848:
pledger = np.array([[s1,s2], [n1-s1, n2-s2]])
chi2, p, dof, ex =stats.chi2_contingency(pledger)
print('The two-sided p-value is %.5f' % (p))
# Analytical Bayes factor:
from scipy.special import gammaln
def lchoose(N, k):
return gammaln(N+1) - gammaln(N-k+1) - gammaln(k+1)
logBF01 = [lchoose(n1, s1) + lchoose(n2, s2) + np.log(n1 + 1) + np.log(n2 + 1)
- lchoose(n1+n2, s1+s2) - np.log(n1 + n2 + 1)]
BF01 = np.exp(logBF01[0])
print('The analytical Bayes factor is %.5f' % (BF01))
The two-sided p-value is 0.00585 The analytical Bayes factor is 0.44974
with pm.Model() as model1:
theta1 = pm.Beta('theta1', alpha=1, beta=1)
theta2 = pm.Beta('theta2', alpha=1, beta=1)
delta = pm.Deterministic('delta', theta1-theta2)
s1 = pm.Binomial('s1', p=theta1, n=n1, observed=s1)
s2 = pm.Binomial('s2', p=theta2, n=n2, observed=s2)
trace1=pm.sample(3e3, model=model1)
burnin = 0
pm.traceplot(trace1[burnin:], varnames=['delta']);
plt.show()
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 100%|██████████| 3500/3500.0 [00:01<00:00, 1990.16it/s]
# BFs based on density estimation (using kernel smoothing instead of spline)
from scipy.stats.kde import gaussian_kde
from scipy.stats import cauchy
pm.summary(trace1, varnames=['delta'])
tmp = pm.df_summary(trace1, varnames=['delta'])
# 95% confidence interval:
x0 = tmp.values[0, 3]
x1 = tmp.values[0, 4]
fig = plt.figure(figsize=(12, 6))
gs = gridspec.GridSpec(1, 2)
ax1 = plt.subplot(gs[0])
t_delt = trace1['delta'][:]
my_pdf = gaussian_kde(t_delt)
x = np.linspace(-3, 3, 100)
ax1.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
ax1.plot(x, cauchy.pdf(x), 'r-', lw=2.5, alpha=0.6, label='Prior')
posterior = my_pdf(0) # this gives the pdf at point delta = 0
prior = cauchy.pdf(0) # height of order-restricted prior at delta = 0
BF01 = posterior/prior
print ('the Bayes Factor is %.5f' % (BF01))
ax1.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
ax2 = plt.subplot(gs[1])
x = np.linspace(-.015, .01, 20)
ax2.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
ax2.plot(x, cauchy.pdf(x), 'r-', lw=2.5, alpha=0.6, label='Prior')
ax2.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
plt.xlim([-.015, .01])
plt.show()
delta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- -0.051 0.018 0.000 [-0.088, -0.016] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| -0.087 -0.063 -0.051 -0.039 -0.015 the Bayes Factor is 1.97820
THe order-restriction of the parameter $\delta < 0$ requires $ \theta_{2} \gt \theta_{1}$. It was address previously in the book by using theta2 ~ dunif(0,1)
and theta1 ~ dunif(0,theta2)
. However, doing so introduce a bias in the prior:
def phi(x):
#'Cumulative distribution function for the standard normal distribution'
return 0.5 + 0.5 * pm.math.erf(x/pm.math.sqrt(2))
def normcdf1(thetap1, thetap2):
angle = 45*np.pi/180
return phi((np.cos(angle) * thetap1) - (np.sin(angle) * tt.abs_(thetap2)))
def normcdf2(thetap1, thetap2):
angle = 45*np.pi/180
return phi((np.sin(angle) * thetap1) + (np.cos(angle) * tt.abs_(thetap2)))
with pm.Model() as modelAE:
# the Approximate method
theta2a = pm.Uniform('theta2a', lower=0, upper=1)
theta1a = pm.Uniform('theta1a', lower=0, upper=theta2a)
deltaa = pm.Deterministic('deltaa', theta1a-theta2a)
## the Exact method
# Adaptation of the exact method as in the book using joint samples from a
# bivariate standard Gaussian then rotating them by 45 degree. The rotated
# sample is transform tinto rates that lie in the unit square
thetap = pm.MvNormal('thetap', mu=[0, 0], tau=pm.math.constant(np.eye(2)), shape=2)
theta1e = pm.Deterministic('theta1e', normcdf1(thetap[0], thetap[1]))
theta2e = pm.Deterministic('theta2e', normcdf2(thetap[0], thetap[1]))
deltae = pm.Deterministic('deltae', theta1e-theta2e)
traceAE=pm.sample(3e3, njobs=2)
fig = plt.figure(figsize=(15, 5))
gs = gridspec.GridSpec(1, 3)
ax1 = plt.subplot(gs[0])
ax1.plot(traceAE['theta1a'][:], traceAE['theta2a'][:], '.', alpha=.1)
plt.xlabel('theta1')
plt.ylabel('theta2')
plt.title('Approximate')
ax2 = plt.subplot(gs[1])
ax2.plot(traceAE['theta1e'][:], traceAE['theta2e'][:], '.', alpha=.1)
plt.xlabel('theta1')
plt.ylabel('theta2')
plt.title('Exact')
from scipy.interpolate import UnivariateSpline
ax2 = plt.subplot(gs[2])
y, binEdges = np.histogram(traceAE['deltaa'][:], bins=20, normed=True)
bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
s = UnivariateSpline(bincenters, y, s=1)
ax2.plot(bincenters,s(bincenters), 'g--', label='Approximate')
y2,binEdges2=np.histogram(traceAE['deltae'][:], bins=20, normed=True)
bincenters2 = 0.5*(binEdges2[1:]+binEdges2[:-1])
s2 = UnivariateSpline(bincenters2, y2, s=1)
ax2.plot(bincenters2, s2(bincenters2), 'r-', label='Exact')
plt.xlabel('theta1-theta2')
plt.ylabel('Density')
plt.title('Bias in the prior')
plt.legend(loc='upper left')
plt.show()
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 100%|██████████| 3500/3500.0 [00:02<00:00, 1335.68it/s]
# pledger data:
s1 = 424
s2 = 5416
n1 = 777
n2 = 9072
with pm.Model() as model2:
## the Exact method
# Adaptation of the exact method as in the book using joint samples from a
# bivariate standard Gaussian then rotating them by 45 degree. The rotated
# sample is transform tinto rates that lie in the unit square
thetap = pm.MvNormal('thetap', mu=[0, 0], tau= pm.math.constant(np.eye(2)), shape=2)
theta1 = pm.Deterministic('theta1', normcdf1(thetap[0], thetap[1]))
theta2 = pm.Deterministic('theta2', normcdf2(thetap[0], thetap[1]))
## the Approximate method
# theta2 = pm.Uniform("theta2",lower=0,upper=1)
# theta1 = pm.Uniform("theta1",lower=0,upper=theta2)
delta = pm.Deterministic('delta', theta1-theta2)
so1 = pm.Binomial('so1', p=theta1, n=n1, observed=s1)
so2 = pm.Binomial('so2', p=theta2, n=n2, observed=s2)
trace2 = pm.sample(3e3, njobs=2)
burnin=0
pm.traceplot(trace2[burnin:], varnames=['delta', 'theta1', 'theta2']);
plt.show()
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 98%|█████████▊| 3428/3500.0 [00:02<00:00, 1218.01it/s]/home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:451: UserWarning: The acceptance probability in chain 0 does not match the target. It is 0.91146674364, but should be close to 0.8. Try to increase the number of tuning steps. % (self._chain_id, mean_accept, target_accept)) 100%|██████████| 3500/3500.0 [00:03<00:00, 1165.58it/s]
pm.summary(trace2, varnames=['delta'])
tmp = pm.df_summary(trace2, varnames=['delta'])
# 95% confidence interval:
x0 = tmp.values[0, 3]
x1 = tmp.values[0, 4]
fig = plt.figure(figsize=(12, 6))
gs = gridspec.GridSpec(1, 2)
ax1 = plt.subplot(gs[0])
t_delt = trace1['delta'][:]
my_pdf = gaussian_kde(t_delt)
x = np.linspace(-3, 3, 100)
ax1.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
ax1.plot(x, cauchy.pdf(x), 'r-', lw=2.5, alpha=0.6, label='Prior')
posterior = my_pdf(0) # this gives the pdf at point delta = 0
prior = cauchy.pdf(0) # height of order-restricted prior at delta = 0
BF01 = posterior/prior
print ('the Bayes Factor is %.5f' %(BF01))
ax1.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
ax2 = plt.subplot(gs[1])
x = np.linspace(-.015, .01, 20)
ax2.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
ax2.plot(x, cauchy.pdf(x), 'r-', lw=2.5, alpha=0.6, label='Prior')
ax2.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
plt.xlim([-.015, .01])
plt.show()
delta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- -0.051 0.018 0.000 [-0.087, -0.015] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| -0.088 -0.063 -0.050 -0.038 -0.016 the Bayes Factor is 1.97820
### Zeelenberg data:
# Study Both:
sb = np.array([15,11,15,14,15,18,16,16,18,16,15,13,18,12,11,13,17,18,16,11,17,18,
12,18,18,14,21,18,17,10,11,12,16,18,17,15,19,12,21,15,16,20,15,19,
16,16,14,18,16,19,17,11,19,18,16,16,11,19,18,12,15,18,20, 8,12,19,
16,16,16,12,18,17,11,20])
nb = 21
# Study Neither:
sn = np.array([15,12,14,15,13,14,10,17,13,16,16,10,15,15,10,14,17,18,19,12,19,18,
10,18,16,13,15,20,13,15,13,14,19,19,19,18,13,12,19,16,14,17,15,16,
15,16,13,15,14,19,12,11,17,13,18,13,13,19,18,13,13,16,18,14,14,17,
12,12,16,14,16,18,13,13])
nn = 21
ns = len(sb)
# two-sided p-value = .03
stats.ttest_rel(sb,sn)
Ttest_relResult(statistic=2.1857463052684696, pvalue=0.032042184096639302)
with pm.Model() as model3:
mu = pm.HalfNormal('mu',sd=1)# standard Gaussian distribution prior. It is known as
# the "unit information prior", as it carries as much
# information as a single observation (Kass & Wasserman, 1995)
sigma = pm.Uniform('sigma', lower=0, upper=10)
delta = pm.HalfNormal('delta', sd=1)
sigma_alpha = pm.Uniform('sigma_alpha', lower=0, upper=10)
mu_alpha = delta*sigma_alpha
alpha_i = pm.Normal('alpha_i', mu=mu_alpha, sd=sigma_alpha, shape=ns)
phin = pm.Normal('phin', mu=mu, sd=sigma, shape=ns)
phib = phin + alpha_i
thetan = pm.Deterministic('thetan', phi(phin))
thetab = pm.Deterministic('thetab', phi(phib))
sno = pm.Binomial('sno', p=thetan, n=nn, observed=sn)
sbo = pm.Binomial('sbo', p=thetab, n=nb, observed=sb)
trace3 = pm.sample(1e4, njobs=2)
pm.traceplot(trace3, varnames=['delta']);
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 54%|█████▍ | 5693/10500.0 [00:33<00:33, 142.46it/s]/home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:451: UserWarning: The acceptance probability in chain 1 does not match the target. It is 0.676967302417, but should be close to 0.8. Try to increase the number of tuning steps. % (self._chain_id, mean_accept, target_accept)) /home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:467: UserWarning: Chain 1 contains 451 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize. % (self._chain_id, n_diverging)) 100%|██████████| 10500/10500.0 [01:01<00:00, 160.11it/s]/home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:451: UserWarning: The acceptance probability in chain 0 does not match the target. It is 0.702263637032, but should be close to 0.8. Try to increase the number of tuning steps. % (self._chain_id, mean_accept, target_accept)) /home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:467: UserWarning: Chain 0 contains 979 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize. % (self._chain_id, n_diverging))
burnin = 0
pm.summary(trace3[burnin:], varnames=['delta'])
tmp = pm.df_summary(trace3[burnin:], varnames=['delta'])
t_delt = trace3['delta'][burnin:]
# 95% confidence interval:
x0 = tmp.values[0, 3]
x1 = tmp.values[0, 4]
fig = plt.figure(figsize=(10, 6))
my_pdf = gaussian_kde(t_delt)
x = np.linspace(0, 4, 100)
plt.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
plt.plot(x, stats.norm.pdf(x)*2, 'r-', lw=2.5, alpha=0.6, label='Prior')
posterior = my_pdf(0) # this gives the pdf at point delta = 0
prior = stats.norm.pdf(0)*2 # height of order-restricted prior at delta = 0
BF01 = posterior/prior
print ('the Bayes Factor is %.5f' %(BF01))
plt.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
plt.show()
delta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- 0.775 0.441 0.023 [0.000, 1.650] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| 0.110 0.457 0.701 1.014 1.875 the Bayes Factor is 0.14434
### Geurts data:
# Normal Controls:
numerrors1 = np.array([15,10,61,11,60,44,63,70,57,11,67,21,89,12,63,11,96,10,37,19,44,
18,78,27,60,14])
nc = np.array([89,74,128,87,128,121,128,128,128,78,128,106,128,83,128,100,128,
73,128,86,128,86,128,100,128,79])
kc = nc - numerrors1
nsc = len(kc)
# ADHD:
numerrors2 = np.array([88,50,58,17,40,18,21,50,21,69,19,29,11,76,46,36,37,72,27,92,13,
39,53,31,49,57,17,10,12,21,39,43,49,17,39,13,68,24,21,27,48,54,
41,75,38,76,21,41,61,24,28,21])
na = np.array([128,128,128,86,128,117,89,128,110,128,93,107,87,128,128,113,128,
128,98,128,93,116,128,116,128,128,93,86,86,96,128,128,128,86,128,
78,128,111,100,95,128,128,128,128,128,128,98,127,128,93,110,96])
ka = na - numerrors2
nsa = len(ka)
# two-sided p-value = .72
stats.ttest_ind(kc/nc,ka/na,equal_var=False)
Ttest_indResult(statistic=-0.36711531382233059, pvalue=0.71545831925014036)
with pm.Model() as model4:
delta = pm.Normal('delta', mu=0, sd=1)
mu_ = pm.Normal('mu', mu=0, sd=1)
sigma = pm.Uniform('sigma', lower=0, upper=10)
alpha = delta*sigma
phic = pm.Normal('phic', mu=mu_+alpha/2, tau=1/sigma**2, shape=nsc)
phia = pm.Normal('phia', mu=mu_-alpha/2, tau=1/sigma**2, shape=nsa)
thetac = pm.Deterministic('thetac', phi(phic))
thetaa = pm.Deterministic('thetaa', phi(phia))
kco = pm.Binomial('kco', p=thetac, n=nc, observed=kc)
kao = pm.Binomial('kao', p=thetaa, n=na, observed=ka)
trace4 = pm.sample(3e3, njobs=2)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 100%|██████████| 3500/3500.0 [00:05<00:00, 697.02it/s]
pm.traceplot(trace4[:], varnames=['delta']);
plt.show()
pm.summary(trace4, varnames=['delta'])
tmp = pm.df_summary(trace4, varnames=['delta'])
# 95% confidence interval:
x0 = tmp.values[0, 3]
x1 = tmp.values[0, 4]
fig = plt.figure(figsize=(10, 6))
t_delt = trace4['delta'][:]
my_pdf = gaussian_kde(t_delt)
x = np.linspace(-3, 3, 100)
plt.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
plt.plot(x, stats.norm.pdf(x), 'r-', lw=2.5, alpha=0.6, label='Prior')
posterior = my_pdf(0) # this gives the pdf at point delta = 0
prior = stats.norm.pdf(0) # height of order-restricted prior at delta = 0
BF01 = posterior/prior
print ('the Bayes Factor is %.5f' %(BF01))
plt.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
plt.show()
delta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- -0.058 0.249 0.003 [-0.548, 0.416] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| -0.548 -0.226 -0.059 0.115 0.416 the Bayes Factor is 3.73140
with pm.Model() as model5:
delta = pm.HalfNormal('delta', sd=1)
mu_ = pm.Normal('mu', mu=0, sd=1)
sigma = pm.Uniform('sigma', lower=0, upper=10)
alpha = delta*sigma
phic = pm.Normal('phic', mu=mu_+alpha/2, sd=sigma, shape=nsc)
phia = pm.Normal('phia', mu=mu_-alpha/2, sd=sigma, shape=nsa)
thetac = pm.Deterministic('thetac', phi(phic))
thetaa = pm.Deterministic('thetaa', phi(phia))
kco = pm.Binomial('kco', p=thetac, n=nc, observed=kc)
kao = pm.Binomial('kao', p=thetaa, n=na, observed=ka)
trace5 = pm.sample(3e3, njobs=2)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... 98%|█████████▊| 3444/3500.0 [00:05<00:00, 755.97it/s]/home/laoj/Documents/Github/pymc3/pymc3/step_methods/hmc/nuts.py:467: UserWarning: Chain 0 contains 2 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize. % (self._chain_id, n_diverging)) 100%|██████████| 3500/3500.0 [00:05<00:00, 662.24it/s]
pm.traceplot(trace5[:], varnames=['delta']);
plt.show()
pm.summary(trace5, varnames=['delta'])
tmp = pm.df_summary(trace5, varnames=['delta'])
# 95% confidence interval:
x0 = tmp.values[0, 3]
x1 = tmp.values[0, 4]
fig = plt.figure(figsize=(10, 6))
t_delt = trace5['delta'][:]
my_pdf = gaussian_kde(t_delt)
x = np.linspace(0, 3, 100)
plt.plot(x, my_pdf(x), '--', lw=2.5, alpha=0.6, label='Posterior') # distribution function
plt.plot(x, stats.norm.pdf(x)*2,'r-', lw=2.5, alpha=0.6, label='Prior')
posterior = my_pdf(0) # this gives the pdf at point delta = 0
prior = stats.norm.pdf(0)*2 # height of order-restricted prior at delta = 0
BF01 = posterior/prior
print ('the Bayes Factor is %.5f' %(BF01))
plt.plot([0, 0], [posterior, prior], 'k-',
[0, 0], [posterior, prior], 'ko',
lw=1.5, alpha=1)
plt.xlabel('Delta')
plt.ylabel('Density')
plt.legend(loc='upper left')
plt.ylim([0,5])
plt.show()
delta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------- 0.173 0.138 0.002 [0.000, 0.445] Posterior quantiles: 2.5 25 50 75 97.5 |--------------|==============|==============|--------------| 0.006 0.064 0.142 0.249 0.514 the Bayes Factor is 2.50601