Dimension error while replicating the Lotka Volterra model example

Questions v5
1

I’m trying to replicate the Lotka Volterra example in the Approximate Bayesian Computation section, and find that I’m getting an dimension mismatch error.

Specifically, while running the sample() method, I get

RuntimeError: The array return by func must be one-dimensional, but got ndim=3.
Apply node that caused the error: Simulator_sim_rv{"(),()->()"}(simulator_rng, [100   2], Exp.0, Exp.0)
Toposort index: 8
Inputs types: [RandomGeneratorType, TensorType(int64, shape=(2,)), TensorType(float64, shape=(1, 1)), TensorType(float64, shape=(1, 1))]
Inputs shapes: ['No shapes', (2,), (1, 1), (1, 1)]
Inputs strides: ['No strides', (8,), (8, 8), (8, 8)]
Inputs values: [Generator(PCG64) at 0x159748200, array([100,   2]), array([[0.1087194]]), array([[0.87183442]])]
Outputs clients: [[output[1](Simulator_sim_rv{"(),()->()"}.0)], [Composite{sqr((0.1 * (i0 - i1)))}(sim{[[10.12490 ... 04511178]]}, Simulator_sim_rv{"(),()->()"}.out)]]

I believe it is because it expects the sum_stat to return a one-dimensional output, but if the sum_stat is the identity function, then it should be able to compute the error with the simulated and observed data, right?

This is the version of PyMC that I’m currently using 5.19.1.

Thanks for your help!

2

Can you share the full code?

3

Yes, it is exactly as in the example. I’m pasting it here for your reference.

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm

print(f"Running on PyMC v{pm.__version__}")

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):
    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)

    idata_lv = pm.sample_smc()

The exact error message is as pasted below.

Traceback (most recent call last):
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/compile/function/types.py", line 992, in __call__
    outputs = vm() if output_subset is None else vm(output_subset=output_subset)
              ^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/graph/op.py", line 531, in rval
    r = p(n, [x[0] for x in i], o)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/tensor/random/op.py", line 402, in perform
    self.rng_fn(rng, *args, None if size is None else tuple(size)),
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pymc/distributions/simulator.py", line 53, in rng_fn
    return cls.fn(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "/var/folders/b9/hmz_c9_95kb0_r42zp984m7m0000gq/T/ipykernel_27808/1881729530.py", line 27, in competition_model
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/scipy/integrate/_odepack_py.py", line 247, in odeint
    output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: The array return by func must be one-dimensional, but got ndim=3.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/concurrent/futures/process.py", line 264, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pymc/smc/sampling.py", line 306, in _sample_smc_int
    smc._initialize_kernel()
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pymc/smc/kernels.py", line 249, in _initialize_kernel
    likelihoods = [self.likelihood_logp_func(sample) for sample in self.tempered_posterior]
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/compile/function/types.py", line 1002, in __call__
    raise_with_op(
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/link/utils.py", line 526, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/compile/function/types.py", line 992, in __call__
    outputs = vm() if output_subset is None else vm(output_subset=output_subset)
              ^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/graph/op.py", line 531, in rval
    r = p(n, [x[0] for x in i], o)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pytensor/tensor/random/op.py", line 402, in perform
    self.rng_fn(rng, *args, None if size is None else tuple(size)),
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/pymc/distributions/simulator.py", line 53, in rng_fn
    return cls.fn(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "/var/folders/b9/hmz_c9_95kb0_r42zp984m7m0000gq/T/ipykernel_27808/1881729530.py", line 27, in competition_model
  File "/Users/maerad/mamba/envs/pymc_env/lib/python3.12/site-packages/scipy/integrate/_odepack_py.py", line 247, in odeint
    output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: The array return by func must be one-dimensional, but got ndim=3.
Apply node that caused the error: Simulator_sim_rv{"(),()->()"}(simulator_rng, [100   2], Exp.0, Exp.0)
Toposort index: 8
Inputs types: [RandomGeneratorType, TensorType(int64, shape=(2,)), TensorType(float64, shape=(1, 1)), TensorType(float64, shape=(1, 1))]
Inputs shapes: ['No shapes', (2,), (1, 1), (1, 1)]
Inputs strides: ['No strides', (8,), (8, 8), (8, 8)]
Inputs values: [Generator(PCG64) at 0x162E0CBA0, array([100,   2]), array([[0.78711716]]), array([[0.64662201]])]
Outputs clients: [[output[1](Simulator_sim_rv{"(),()->()"}.0)], [Composite{sqr((0.1 * (i0 - i1)))}(sim{[[ 8.89815 ... 51544454]]}, Simulator_sim_rv{"(),()->()"}.out)]]

The scipy version I’m using is 1.14.1