-
Notifications
You must be signed in to change notification settings - Fork 11
/
reaction_kinetics_model.py
44 lines (38 loc) · 1.3 KB
/
reaction_kinetics_model.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
"""
The definition of the posterior of the reaction kinetics problem.
"""
import pymc
import warnings
warnings.simplefilter('ignore')
from reaction_kinetics_solver import *
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()
if __name__ == '__main__':
# Generate the observations
import pickle
re_solver = ReactionKineticsSolver()
k1 = 2
k2 = 4
data = {}
data['y_obs'] = re_solver(k1, k2)
with open('reaction_kinetics_data.pickle', 'wb') as fd:
pickle.dump(data, fd)