!date
Wed Dec 18 13:54:39 PST 2013
import pymc as pm, theano.tensor as T
x = array([ 0, 1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55,
60, 65, 70, 75, 80, 85, 90, 95, 100, 105])
y = array([-5.31126213, -6.88284349, -7.28148079, -7.20912457, -6.06006241,
-5.69987917, -5.72478151, -5.62202549, -5.36570549, -4.96331167,
-4.50282001, -3.99181652, -3.44459009, -2.88168406, -2.35241652,
-1.82025242, -1.25903034, -0.66321015, 0.06458783, 0.87678754,
1.70916784, 2.56996703, 3.38351035])
with pm.Model() as m:
a = pm.Flat('a')
b = pm.Flat('b')
c = pm.Flat('c')
d = pm.Flat('d')
e = pm.Flat('e')
f = pm.Flat('f')
g = pm.Flat('g')
h = pm.Flat('h')
t1 = a**((x+b)**c)
t2 = d * T.exp(-e * T.log(x/f)**2)
t3 = g*h**x
y_pred = t1 + t2 + t3
y_obs = pm.Normal('y_obs', mu=y_pred/y, sd=1., observed=ones_like(y))
with m: trace = pm.sample(20000, pm.Metropolis())
[-----------------100%-----------------] 20000 of 20000 complete in 12.9 sec
with m: trace = pm.sample(20000, pm.NUTS())