Skip to content

Approximate Bayesian Computation is broken #787

Open
@Gabriel-p

Description

@Gabriel-p

Notebook title: Approximate Bayesian Computation
Notebook url: https://www.pymc.io/projects/examples/en/latest/samplers/SMC-ABC_Lotka-Volterra_example.html

Issue description

The 'Old good Gaussian fit' example runs with no issues, the 'Lotka–Volterra' example fails with:

---------------------------------------------------------------------------
_RemoteTraceback                          Traceback (most recent call last)
_RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/usr/local/lib/python3.11/dist-packages/pytensor/compile/function/types.py", line 1037, in __call__
    outputs = vm() if output_subset is None else vm(output_subset=output_subset)
              ^^^^
  File "/usr/local/lib/python3.11/dist-packages/pytensor/graph/op.py", line 531, in rval
    r = p(n, [x[0] for x in i], o)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pytensor/tensor/random/op.py", line 402, in perform
    self.rng_fn(rng, *args, None if size is None else tuple(size)),
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pymc/distributions/simulator.py", line 53, in rng_fn
    return cls.fn(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "<ipython-input-12-05bfa58b7eb6>", line 27, in competition_model
    return odeint(dX_dt, y0=X0, t=t, rtol=0.01, args=(a, b, c, d))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-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 "/usr/lib/python3.11/concurrent/futures/process.py", line 261, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/sampling.py", line 316, in _sample_smc_int
    smc._initialize_kernel()
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 253, in _initialize_kernel
    likelihoods = [self.likelihood_logp_func(sample) for sample in self.tempered_posterior]
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 253, in <listcomp>
    likelihoods = [self.likelihood_logp_func(sample) for sample in self.tempered_posterior]
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pytensor/compile/function/types.py", line 1047, in __call__
    raise_with_op(
  File "/usr/local/lib/python3.11/dist-packages/pytensor/link/utils.py", line 526, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/usr/local/lib/python3.11/dist-packages/pytensor/compile/function/types.py", line 1037, in __call__
    outputs = vm() if output_subset is None else vm(output_subset=output_subset)
              ^^^^
  File "/usr/local/lib/python3.11/dist-packages/pytensor/graph/op.py", line 531, in rval
    r = p(n, [x[0] for x in i], o)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pytensor/tensor/random/op.py", line 402, in perform
    self.rng_fn(rng, *args, None if size is None else tuple(size)),
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-packages/pymc/distributions/simulator.py", line 53, in rng_fn
    return cls.fn(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "<ipython-input-12-05bfa58b7eb6>", line 27, in competition_model
    return odeint(dX_dt, y0=X0, t=t, rtol=0.01, args=(a, b, c, d))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.11/dist-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: 6
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 0x793866AA9EE0, array([100,   2]), array([[1.58880866]]), array([[2.39588877]])]
Outputs clients: [[output[1](Simulator_sim_rv{"(),()->()"}.0)], [Composite{sqr((0.1 * (i0 - i1)))}(sim{[[10.97461 ... 67732525]]}, Simulator_sim_rv{"(),()->()"}.out)]]

Backtrace when the node is created (use PyTensor flag traceback__limit=N to make it longer):
  File "/usr/lib/python3.11/multiprocessing/process.py", line 314, in _bootstrap
    self.run()
  File "/usr/lib/python3.11/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File "/usr/lib/python3.11/concurrent/futures/process.py", line 261, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/sampling.py", line 316, in _sample_smc_int
    smc._initialize_kernel()
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 248, in _initialize_kernel
    self.likelihood_logp_func = _logp_forw(
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 645, in _logp_forw
    out_list, inarray0 = join_nonshared_inputs(
  File "/usr/local/lib/python3.11/dist-packages/pymc/pytensorf.py", line 588, in join_nonshared_inputs
    new_outputs = [
  File "/usr/local/lib/python3.11/dist-packages/pymc/pytensorf.py", line 589, in <listcomp>
    pytensor.clone_replace(output, replace, rebuild_strict=False) for output in outputs

HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.
"""

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
[<ipython-input-16-264bb474e972>](https://localhost:8080/#) in <cell line: 0>()
      5     sim = pm.Simulator("sim", competition_model, params=(a, b), epsilon=10, observed=observed)
      6 
----> 7     idata_lv = pm.sample_smc()

4 frames
[/usr/lib/python3.11/concurrent/futures/_base.py](https://localhost:8080/#) in __get_result(self)
    399         if self._exception:
    400             try:
--> 401                 raise self._exception
    402             finally:
    403                 # Break a reference cycle with the exception in self._exception

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: 6
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 0x793866AA9EE0, array([100,   2]), array([[1.58880866]]), array([[2.39588877]])]
Outputs clients: [[output[1](Simulator_sim_rv{"(),()->()"}.0)], [Composite{sqr((0.1 * (i0 - i1)))}(sim{[[10.97461 ... 67732525]]}, Simulator_sim_rv{"(),()->()"}.out)]]

Backtrace when the node is created (use PyTensor flag traceback__limit=N to make it longer):
  File "/usr/lib/python3.11/multiprocessing/process.py", line 314, in _bootstrap
    self.run()
  File "/usr/lib/python3.11/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File "/usr/lib/python3.11/concurrent/futures/process.py", line 261, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/sampling.py", line 316, in _sample_smc_int
    smc._initialize_kernel()
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 248, in _initialize_kernel
    self.likelihood_logp_func = _logp_forw(
  File "/usr/local/lib/python3.11/dist-packages/pymc/smc/kernels.py", line 645, in _logp_forw
    out_list, inarray0 = join_nonshared_inputs(
  File "/usr/local/lib/python3.11/dist-packages/pymc/pytensorf.py", line 588, in join_nonshared_inputs
    new_outputs = [
  File "/usr/local/lib/python3.11/dist-packages/pymc/pytensorf.py", line 589, in <listcomp>
    pytensor.clone_replace(output, replace, rebuild_strict=False) for output in outputs

HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

Expected output

If applicable, describe what should happen instead.

Proposed solution

If applicable, explain possible solutions and workarounds.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions