Hello,
Well the example in the page also does not work for me but it seems to be a sizing issue. The simulator class seems to resize everything to two dimensional arrays so for instance constants a,b become 1x1 arrays which breaks the dimension of odeint returning 1 x 1 x N sized arrays and an error. If I slightly modify it, it works for me and produces roughly the same results as in the page. I wonder if this resizing is an intended feature because it seems to make life generally difficult.
The page needs to be updated though.
import numpy as np
import pymc as pm
import arviz as az
from scipy.integrate import odeint
# Definition of parameters
a = 1.0
b = 0.1
c = 1.5
d = 0.75
# initial population of rabbits and foxes
X0 = [10.0, 5.0]
# size of data
size = 100
# time lapse
time = 15
t = np.linspace(0, time, size)
# Lotka - Volterra equation
def dX_dt(X, t, a, b, c, d):
"""Return the growth rate of fox and rabbit populations."""
return np.array([a * X[0] - b * X[0] * X[1], -c * X[1] + d * b * X[0] * X[1]])
# simulator function
def competition_model(rng, a, b, size=None):
if isinstance(a, np.ndarray) and a.ndim==2:
a = a[0,0]
if isinstance(b, np.ndarray) and b.ndim==2:
b = b[0,0]
return odeint(dX_dt, y0=X0, t=t, rtol=0.01, args=(a, b, c, d))
# function for generating noisy data to be used as observed data.
def add_noise(a, b):
noise = np.random.normal(size=(size, 2))
simulated = competition_model(None, a, b) + noise
return simulated
observed = add_noise(a, b)
with pm.Model() as model_lv:
a = pm.HalfNormal("a", 1.0)
b = pm.HalfNormal("b", 1.0)
sim = pm.Simulator("sim", competition_model, params=(a, b), epsilon=10, observed=observed)
sim.eval()
idata_lv = pm.sample_smc()
az.plot_trace(idata_lv, kind="rank_vlines")