!date
import numpy as np, pymc as pm, matplotlib.pyplot as plt, seaborn as sns
%matplotlib inline
Mon Feb 2 15:07:16 PST 2015
# set random seed for reproducibility
np.random.seed(12345)
@pm.stochastic
def R1(value=1500):
return 1./value
@pm.stochastic
def R2(value=1500):
return 1./value
@pm.observed
def d1(value=2237, R=R1, nu=31):
x = nu * value/R
return pm.chi2_like(x, nu)
@pm.observed
def d2(value=1347, R=R2, nu=61):
x = nu * value/R
return pm.chi2_like(x, nu)
m = pm.MCMC([R1, R2, d1, d2])
m.sample(iter=200000, burn=100000, thin=10)
[-----------------100%-----------------] 200000 of 200000 complete in 22.9 sec
pm.Matplot.plot(m)
Plotting R1 Plotting R2
print 'Pr[R2 <= R1|data] =', np.round(np.mean(R2.trace() <= R1.trace()), 4)
Pr[R2 <= R1|data] = 0.9693
Expecting 0.9574