In [19]:
%matplotlib inline
from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
import pickle

In [25]:
with open("/Users/philngo/Downloads/may_air.pkl") as f:
    models = pickle.load(f)["air1"]

In [51]:
means = []
covars = []
transmats = []
for model in models:
    means.append(models[model]._means_)
    covars.append([covar[0] for covar in models[model]._covars_])
    transmats.append(models[model].transmat_)
means = np.array(means)
covars = np.array(covars)
transmats = np.array(transmats)
print means[0]
print covars[0]
print transmats[0]


[[  7.93843724e-08]
 [  2.02002087e+00]
 [  2.37664225e+00]]
[[  2.64914591e-07]
 [  1.22243029e-01]
 [  5.19399436e-03]]
[[  9.95738711e-01   3.83217220e-03   4.29116343e-04]
 [  3.74575688e-02   9.52497364e-01   1.00450677e-02]
 [  1.38188267e-03   2.27407863e-02   9.75877331e-01]]

In [75]:
means = np.array([1,2,3,4,5,5,4,3,3,4,3,4,1,2,3,6,7,9,8,5,4,6,5,10,5])
alpha,loc,beta = gamma.fit(means)

print alpha, loc, beta

n = 300
generated_data = gamma.rvs(alpha,loc=loc,scale=beta,size=n)
print generated_data[0]


5.10671791388 -0.624047478884 0.999476690546
6.91942393258

In [78]:
fig, ax = plt.subplots(1, 1)

x = np.linspace(gamma.ppf(0.01,alpha,loc=loc,scale=beta),
                gamma.ppf(0.99,alpha,loc=loc,scale=beta), 100)

ax.plot(x, gamma.pdf(x, alpha,loc=loc,scale=beta),'r-', lw=5, alpha=0.6, label='gamma pdf')

ax.hist(means, normed=True, histtype='stepfilled', alpha=0.2)
ax.hist(generated_data, normed=True, histtype='stepfilled', alpha=0.2)

ax.legend(loc='best', frameon=False)
plt.show()


Make a real model. - see mcmc version which uses PyMC


In [ ]: