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

# In[51]:


import numpy as np
import scipy.stats as st
thetas = st.beta(8,2).rvs(15)
total = np.array([3, 5, 10, 15, 20, 100, 200, 300, 500, 1000, 1200, 1300, 1500, 2000, 3000])
data = np.random.binomial(n=total, p=thetas)


# In[52]:


data[0]/3


# In[53]:


thetas


# In[69]:


import pymc3 as pm
import scipy.optimize as opt
with pm.Model() as model:
    mu = pm.Beta('mu', alpha=1, beta=1)
    sd = pm.HalfCauchy('sd', 1)
    p = pm.Beta('p', mu=mu, sd=sd, shape=15)
    like = pm.Binomial('like', n=total, p=p, observed=data)

    start = pm.find_MAP(fmin=opt.fmin_l_bfgs_b)
    step = pm.NUTS(scaling=start)
    tr = pm.sample(50000, step, start=start)


# In[70]:


pm.forestplot(tr[-10000:], vars=['p'])


# In[34]:


list(zip(tr['p'][-30000:].mean(axis=0), thetas, data, total))


# In[16]:


import pymc3 as pm
import scipy.optimize as opt
with pm.Model() as model2:
    alpha = pm.Gamma('alpha', alpha=1, beta=0.01)
    beta = pm.Gamma('beta', alpha=1, beta=0.01)
    p = pm.Beta('p', alpha, beta, shape=15)
    like = pm.Binomial('like', n=total, p=p, observed=data)

    start = pm.find_MAP(fmin=opt.fmin_l_bfgs_b)
    step = pm.NUTS(scaling=start)
    tr2 = pm.sample(20000, step, start=start)


# In[17]:


list(zip(tr2['p'].mean(axis=0), thetas, data, total))


# In[54]:


import pymc

def beta_binomial():
    
    alpha = pymc.Gamma('alpha', alpha=1, beta=0.01)
    beta = pymc.Gamma('beta', alpha=1, beta=0.01)
    p = pymc.Beta('p',alpha, beta, value=[0.5]*15)
    like = pymc.Binomial('like', n=total, p=p, observed=True, value=data)
    
    return locals()



# In[60]:


M = pymc.MCMC(beta_binomial())


# In[61]:


M.sample(50000, 40000)
M.sample(50000, 40000)


# In[62]:


pymc.Matplot.summary_plot(M.p)


# In[ ]: