Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def make_model():
import pickle
with open('reaction_kinetics_data.pickle', 'rb') as fd:
data = pickle.load(fd, encoding='latin1')
y_obs = data['y_obs']
# The priors for the reaction rates:
k1 = pymc.Lognormal('k1', mu=2, tau=1./(10. ** 2), value=5.)
k2 = pymc.Lognormal('k2', mu=4, tau=1./(10. ** 2), value=5.)
# The noise term
#sigma = pymc.Uninformative('sigma', value=1.)
sigma = pymc.Exponential('sigma', beta=1.)
# The forward model
re_solver = ReactionKineticsSolver()
@pymc.deterministic
def model_output(value=None, k1=k1, k2=k2):
return re_solver(k1, k2)
# The likelihood term
@pymc.stochastic(observed=True)
def output(value=y_obs, mod_out=model_output, sigma=sigma, gamma=1.):
return gamma * pymc.normal_like(y_obs, mu=mod_out, tau=1/sigma ** 2)
return locals()
def estimate_alpha_beta(data, n_iter=1000, plot=False):
# data = data.dropna()
alpha_var = pm.Exponential('alpha', .5)
beta_var = pm.Exponential('beta', .5)
observations = pm.Beta('observations', alpha_var, beta_var, value=data,
observed=True)
model = pm.Model([alpha_var, beta_var])
mcmc = pm.MCMC(model)
mcmc.sample(n_iter)
if plot:
from pymc.Matplot import plot
plot(mcmc)
sns.despine()
alphas = mcmc.trace('alpha')[:]
betas = mcmc.trace('beta')[:]
disasters_array = array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
n = len(disasters_array)
# Define data and stochastics
switchpoint = pm.DiscreteUniform('switchpoint',lower=0,upper=110)
early_mean = pm.Exponential('early_mean',beta=1.)
late_mean = pm.Exponential('late_mean',beta=1.)
@pm.stochastic(observed=True, dtype=int)
def disasters( value = disasters_array,
early_mean = early_mean,
late_mean = late_mean,
switchpoint = switchpoint):
"""Annual occurences of coal mining disasters."""
return pm.poisson_like(value[:switchpoint],early_mean) + pm.poisson_like(value[switchpoint:],late_mean)
@pm.deterministic
def disasters_sim(early_mean = early_mean,
late_mean = late_mean,
switchpoint = switchpoint):
"""Coal mining disasters sampled from the posterior predictive distribution"""
return concatenate( (pm.rpoisson(early_mean, size=switchpoint), pm.rpoisson(late_mean, size=n-switchpoint)))
e_r,c,t,a ~ N(0, (gamma * W_r,c,t,a)**2 + sigma_e**2)
"""
# covariates
K1 = count_covariates(data, 'x')
K2 = count_covariates(data, 'w')
X = pl.array([data['x%d'%i] for i in range(K1)])
W = pl.array([data['w%d'%i] for i in range(K2)])
# priors
beta = mc.Laplace('beta', mu=0., tau=1., value=pl.zeros(K1))
gamma = mc.Exponential('gamma', beta=1., value=pl.zeros(K2))
sigma_e = mc.Exponential('sigma_e', beta=1., value=1.)
# hyperpriors for GPs (These seem to really matter!)
sigma_f = mc.Exponential('sigma_f', beta=1., value=[1., 1., .1])
tau_f = mc.Truncnorm('tau_f', mu=25., tau=5.**-2, a=10, b=pl.inf, value=[25., 25., 25.])
diff_degree = [2., 2., 2.]
# fixed-effect predictions
@mc.deterministic
def mu(X=X, beta=beta):
""" mu_i,r,c,t,a = beta * X_i,r,c,t,a"""
return pl.dot(beta, X)
@mc.deterministic
def sigma_explained(W=W, gamma=gamma):
""" sigma_explained_i,r,c,t,a = gamma * W_i,r,c,t,a"""
return pl.dot(pl.atleast_1d(gamma), pl.atleast_2d(W))
# GP random effects
## make index dicts to convert from region/country/age to array index
"""
__all__ = ['s','e','l','r','D']
from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform
import numpy as np
disasters_array = np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
s = DiscreteUniform('s', lower=0, upper=110)
e = Exponential('e', beta=1)
l = Exponential('l', beta=1)
@deterministic(plot=False)
def r(s=s, e=e, l=l):
"""Concatenate Poisson means"""
out = np.empty(len(disasters_array))
out[:s] = e
out[s:] = l
return out
D = Poisson('D', mu=r, value=disasters_array, observed=True)
if __name__ == '__main__':
from pymc import MCMC, Metropolis
M = MCMC([s,e,l,D], db='hdf5')
M.use_step_method(Metropolis, e, tally=True)
def make_gp_submodel(suffix, mesh, africa_val=None, with_africa_covariate=False):
"""
A small function that creates the mean and covariance object
of the random field.
"""
# from duffy import cut_matern
# The partial sill.
amp = pm.Exponential('amp_%s'%suffix, .1, value=1.)
# The range parameter. Units are RADIANS.
# 1 radian = the radius of the earth, about 6378.1 km
# scale = pm.Exponential('scale', 1./.08, value=.08)
scale = pm.Exponential('scale_%s'%suffix, 1, value=.08)
scale_in_km = scale*6378.1
# The nugget variance. Lower-bounded to preserve mixing.
V = pm.Exponential('V_%s'%suffix, 1, value=1.)
@pm.potential
def V_bound(V=V):
if V<.1:
return -np.inf
else:
return 0
# Create the covariance & its evaluation at the data locations.
@pm.deterministic(trace=True,name='C_%s'%suffix)
def C(amp=amp, scale=scale):
return pm.gp.FullRankCovariance(pm.gp.exponential.geo_rad, amp=amp, scale=scale)
# Changed beta for e
pm.Matplot.autocorrelation(S.e)
### SETUP ###
disasters_array = numpy.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
switchpoint = pm.DiscreteUniform('s', lower=0, upper=110)
early_mean = pm.Exponential('e', beta=1)
late_mean = pm.Exponential('l', beta=1)
@pm.stochastic(observed=True, dtype=int)
def disasters( value = disasters_array,
early_mean = early_mean,
late_mean = late_mean,
switchpoint = switchpoint):
"""Annual occurences of coal mining disasters."""
return pm.poisson_like(value[:switchpoint],early_mean) + \
pm.poisson_like(value[switchpoint:],late_mean)
@pm.deterministic
def disasters_sim(early_mean = early_mean,
late_mean = late_mean,
switchpoint = switchpoint):
pass
# ========================================================
# = Map BUGS distributions to PyMC Stochastic subclasses =
# ========================================================
bugs_dists = {'bern': (pm.Bernoulli, 'p'),
'bin': (pm.Binomial, 'p', 'n'),
'cat': (pm.Categorical, 'p'),
# 'negbin': (pm.NegativeBinomial, '') Need to implement standard parameterization or else translate with a Deterministic.
'pois': (pm.Poisson, 'mu'),
'beta': (pm.Beta, 'alpha', 'beta'),
'chisqr': (pm.Chi2, 'nu'),
# 'dexp': Double exponential distribution not implemented.
'exp': (pm.Exponential, 'beta'),
'gamma': (pm.Gamma, 'alpha', 'beta'),
# 'gen.gamma': Not implemented
'lnorm': (pm.Lognormal, 'mu', 'tau'),
# 'logis': Logistic distribution not implemented
'norm': (pm.Normal, 'mu', 'tau'),
# 'par': Pareto distribution not implemented.
# 't': T distribution not implemented !?
'unif': (pm.Uniform, 'lower', 'upper'),
# 'weib': Uses different parameterization than we do.
'multi': (pm.Multinomial, 'p', 'n'),
# 'dirch': Need to apply CompletedDirichlet
'mnorm': (pm.MvNormal, 'mu', 'tau'),
# 'mt': Multivariate student's T not implemented
'wish': (pm.Wishart, 'T', 'n')}
def make_gp_submodel(suffix, mesh, africa_val=None, with_africa_covariate=False):
"""
A small function that creates the mean and covariance object
of the random field.
"""
# from duffy import cut_matern
# The partial sill.
amp = pm.Exponential('amp_%s'%suffix, .1, value=1.)
# The range parameter. Units are RADIANS.
# 1 radian = the radius of the earth, about 6378.1 km
# scale = pm.Exponential('scale', 1./.08, value=.08)
scale = pm.Exponential('scale_%s'%suffix, 1, value=.08)
scale_in_km = scale*6378.1
# The nugget variance. Lower-bounded to preserve mixing.
V = pm.Exponential('V_%s'%suffix, 1, value=1.)
@pm.potential
def V_bound(V=V):
if V<.1:
return -np.inf
else:
return 0