Approximate Bayesian Computation is broken

[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.

[fonnesbeck]

Yep, looks like the competition model is getting vector-valued parameters for some reason.