from pyndamics import Simulation
from pyndamics.mcmc import MCMCModel

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]')

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)

model=MCMCModel(sim,{'a':[-10,10]})
model.fit(iter=25000)

model.a

sim.run(0,12)

model.plot_distributions()

model.plot_distributions(show_normal=True)

model=MCMCModel(sim,{'a':[-10,10],'initial_h':[0,4]})
model.fit(iter=25000)

sim.run(0,12)
model.plot_distributions(show_normal=True)

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

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)

model=MCMCModel(sim,{'a':[0.001,5],'K':[0.1,40],'initial_h':[0,4]})
model.fit(iter=25000)
print sim.a,sim.K

sim.run(0,12)
model.plot_distributions(show_normal=True)

model.plot_joint_distribution('a','K',show_prior=False)

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]')

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]')