Pyndamics provides a way to describe a dynamical system in terms of the differential equations, or the stock-flow formalism. It is a wrapper around the Scipy odeint function, with further functionality for time plots, phase plots, and vector fields. The MCMC component of this package uses pymc (version 2 right now): http://pymc-devs.github.io/pymc/.
Page for this package: https://code.google.com/p/pyndamics/
from pyndamics import Simulation
from pyndamics.mcmc import MCMCModel
Data from http://www.seattlecentral.edu/qelp/sets/009/009.html
t=array([0,2,4,6,8,10,12])
h=array([1.0,1.5,2.2,3.2,4.3,5.2,5.6])
plot(t,h,'-o')
xlabel('Days')
ylabel('Height [cm]')
<matplotlib.text.Text at 0x1095d1a10>
Here, the constant value ($a=1$) is hand-picked, and doesn't fit the data particularly well.
sim=Simulation()
sim.add("h'=a",1,plot=True)
sim.add_data(t=t,h=h,plot=True)
sim.params(a=1)
sim.run(0,12)
<matplotlib.figure.Figure at 0x109771e90>
Specifying the prior probability distribution for $a$ as uniform between -10 and 10.
model=MCMCModel(sim,{'a':[-10,10]})
model.fit(iter=25000)
[-----------------100%-----------------] 25000 of 25000 complete in 15.7 sec Time taken: 15.76 s
What is the best fit parameter value?
model.a
0.3917653066312086
sim.run(0,12)
<matplotlib.figure.Figure at 0x10979a0d0>
model.plot_distributions()
model.plot_distributions(show_normal=True)
model=MCMCModel(sim,{'a':[-10,10],'initial_h':[0,4]})
model.fit(iter=25000)
[-----------------100%-----------------] 25000 of 25000 complete in 32.1 sec Time taken: 32.32 s
sim.run(0,12)
model.plot_distributions(show_normal=True)
<matplotlib.figure.Figure at 0x10a4c90d0>
sim.noplots=True # turn off the simulation plots
for i in range(500):
model.draw()
sim.run(0,12)
plot(sim.t,sim.h,'g-',alpha=.1)
sim.noplots=False # gotta love a double-negative
plot(t,h,'bo') # plot the data
[<matplotlib.lines.Line2D at 0x10977a9d0>]
sim=Simulation()
sim.add("h'=a*h*(1-h/K)",1,plot=True)
sim.add_data(t=t,h=h,plot=True)
sim.params(a=1,K=10)
sim.run(0,12)
<matplotlib.figure.Figure at 0x10972d450>
model=MCMCModel(sim,{'a':[0.001,5],'K':[0.1,40],'initial_h':[0,4]})
model.fit(iter=25000)
print sim.a,sim.K
[-----------------100%-----------------] 25000 of 25000 complete in 197.4 sec Time taken: 3 m, 19.49 s 0.300930890184 6.66506873289
sim.run(0,12)
model.plot_distributions(show_normal=True)
<matplotlib.figure.Figure at 0x10befeed0>
model.plot_joint_distribution('a','K',show_prior=False)
...showing the joint, with the prior, to see how much the inference constrains the final best estimates of the parameters
model.plot_joint_distribution('a','K',show_prior=True)
sim.noplots=True
for i in range(500):
model.draw()
sim.run(0,12)
plot(sim.t,sim.h,'g-',alpha=.1)
sim.noplots=False # gotta love a double-negative
plot(t,h,'bo')
xlabel('Days')
ylabel('Height [cm]')
<matplotlib.text.Text at 0x112f92750>
plot(t,model['h_data'].value,'o')
plot(t,model['h'].stats()['mean'],'--')
plot(t,model['h'].stats()['95% HPD interval'], 'r:')
xlabel('Days')
ylabel('Height [cm]')
<matplotlib.text.Text at 0x10f701890>