#!/usr/bin/env python
# coding: utf-8

# In[11]:


get_ipython().system('date')
import numpy as np, pymc as pm, matplotlib.pyplot as plt, seaborn as sns
get_ipython().run_line_magic('matplotlib', 'inline')


# In[12]:


# set random seed for reproducibility
np.random.seed(12345)


# In[13]:


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


# In[14]:


pm.Matplot.plot(m)


# In[15]:


print 'Pr[R2 <= R1|data] =', np.round(np.mean(R2.trace() <= R1.trace()), 4)


# Expecting 0.9574

# In[ ]: