BUG: ZeroDivisionError when running OrderedLogistic example from the documentation

[LucHeuff]

Describe the issue:

Running the example code for pymc.OrderedLogistic results in a ZeroDivisionError during sampling.

Reproduceable code example:

# this is simply the example from the documentation

# Generate data for a simple 1 dimensional example problem
n1_c = 300; n2_c = 300; n3_c = 300
cluster1 = np.random.randn(n1_c) + -1
cluster2 = np.random.randn(n2_c) + 0
cluster3 = np.random.randn(n3_c) + 2

x = np.concatenate((cluster1, cluster2, cluster3))
y = np.concatenate((1*np.ones(n1_c),
                    2*np.ones(n2_c),
                    3*np.ones(n3_c))) - 1

# Ordered logistic regression
with pm.Model() as model:
    cutpoints = pm.Normal("cutpoints", mu=[-1,1], sigma=10, shape=2,
                          transform=pm.distributions.transforms.ordered)
    y_ = pm.OrderedLogistic("y", cutpoints=cutpoints, eta=x, observed=y)
    idata = pm.sample()

Error message:

pymc.sampling.parallel.RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pytensor/link/basic.py", line 6
65, in thunk
    outputs = fgraph_jit(*(x[0] for x in thunk_inputs))
ZeroDivisionError: division by zero

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/parallel.py", lin
e 154, in run
    self._start_loop()
    ~~~~~~~~~~~~~~~~^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/parallel.py", lin
e 211, in _start_loop
    point, stats = self._step_method.step(self._point)
                   ~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/arraystep.py"
, line 116, in step
    apoint, stats = self.astep(q)
                    ~~~~~~~~~~^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/base_hmc.
py", line 232, in astep
    hmc_step = self._hamiltonian_step(start, p0, step_size)
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/nuts.py",
 line 216, in _hamiltonian_step
    divergence_info, turning = tree.extend(direction)
                               ~~~~~~~~~~~^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/nuts.py",
 line 362, in extend
    tree, diverging, turning = self._build_subtree(
                               ~~~~~~~~~~~~~~~~~~~^
        self.left, self.depth, np.asarray(-self.step_size, dtype=self.floatX)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/nuts.py",
 line 450, in _build_subtree
    return self._single_step(left, epsilon)
           ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/nuts.py",
 line 406, in _single_step
    right = self.integrator.step(epsilon, left)
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/integrati
on.py", line 94, in step
    return self._step(epsilon, state)
           ~~~~~~~~~~^^^^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/integrati
on.py", line 127, in _step
    logp, q_new_grad = self._logp_dlogp_func(q_new, q.point_map_info)
                       ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/step_methods/hmc/integrati
on.py", line 50, in func
    return pytensor_function(q)
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pytensor/compile/function/types
.py", line 1039, in __call__
    outputs = vm() if output_subset is None else vm(output_subset=output_subset)
              ~~^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pytensor/link/basic.py", line 6
69, in thunk
    raise_with_op(self.fgraph, output_nodes[0], thunk)
    ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pytensor/link/utils.py", line 5
26, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pytensor/link/basic.py", line 6
65, in thunk
    outputs = fgraph_jit(*(x[0] for x in thunk_inputs))
ZeroDivisionError: division by zero
Apply node that caused the error: Add(Sum{axes=None}.0, Subtensor{i}.0, Sum{axes=None}.0)
Toposort index: 45
Inputs types: [TensorType(float64, shape=()), TensorType(float64, shape=()), TensorType(float64, shape=())]
Inputs shapes: [(2,)]
Inputs strides: [(8,)]
Inputs values: [array([-13.88371585, -48.70597313])]
Outputs clients: [[output[0](Add.0)]]

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was c
reated. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor op
timizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this App
ly node.
"""

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

ZeroDivisionError: division by zero
Apply node that caused the error: Add(Sum{axes=None}.0, Subtensor{i}.0, Sum{axes=None}.0)
Toposort index: 45
Inputs types: [TensorType(float64, shape=()), TensorType(float64, shape=()), TensorType(float64, shape=())]
Inputs shapes: [(2,)]
Inputs strides: [(8,)]
Inputs values: [array([-13.88371585, -48.70597313])]
Outputs clients: [[output[0](Add.0)]]

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was c
reated. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor op
timizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this App
ly node.

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

Traceback (most recent call last):
  File "/home/wsl/python/statistical-rethinking/chapter_12.py", line 631, in <module>
    test()
    ~~~~^^
  File "/home/wsl/python/statistical-rethinking/chapter_12.py", line 625, in test
    samples = pm.sample()
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/mcmc.py", line 93
5, in sample
    _mp_sample(**sample_args, **parallel_args)
    ~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/mcmc.py", line 14
11, in _mp_sample
    for draw in sampler:
                ^^^^^^^
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/parallel.py", lin
e 510, in __iter__
    draw = ProcessAdapter.recv_draw(self._active)
  File "/home/wsl/python/statistical-rethinking/.venv/lib/python3.13/site-packages/pymc/sampling/parallel.py", lin
e 381, in recv_draw
    raise error from old_error
pymc.sampling.parallel.ParallelSamplingError: Chain 2 failed with: division by zero
Apply node that caused the error: Add(Sum{axes=None}.0, Subtensor{i}.0, Sum{axes=None}.0)
Toposort index: 45
Inputs types: [TensorType(float64, shape=()), TensorType(float64, shape=()), TensorType(float64, shape=())]
Inputs shapes: [(2,)]
Inputs strides: [(8,)]
Inputs values: [array([-13.88371585, -48.70597313])]
Outputs clients: [[output[0](Add.0)]]

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was c
reated. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor op
timizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this App
ly node.

PyMC version information:

PyMC Version: PyMC v5.23.0 PyTensor Version v2.31.4 Operation system: Ubuntu 22.04 through WSL 2.0 How did you install PyMC: uv

Context for the issue:

I'm teaching myself how to use pymc for Bayesians statistics, so for me this is mostly an inconvenience, but I'm not sure this is a wider issue with OrderedLogistic models.

[welcome]

]
🎉 Welcome to PyMC! 🎉 We're really excited to have your input into the project! 💖

If you haven't done so already, please make sure you check out our Contributing Guidelines and Code of Conduct.

[ricardoV94]

Are you using numba mode? The c backend doesn't raise zero division error at runtime so that's perhaps why we didn't see it in these models

[LucHeuff]

I was sampling using the default pymc sampler (I think that's what you mean with numba mode?).

I just tried the nutpie sampler and that runs without issues.

[ricardoV94]

From your output it doesn't look like the default backend (backend is different than sampler, you can use Python/C/Numba/JAX backend with the default pymc sampler). Nutpie uses numba backend by default while PyMC uses C or Python if C isn't available.

Do you get a warning about not having a C compiler installed or something?

[LucHeuff]

I don't, but on earlier versions of pymc I used to get warnings (something about falling back to NumPy-C compiler from the top of my head).

I also tested removing numba from my venv, at which point the default sampler (so pm.sample()) fails on an ImportError for numba.

I have been able to succesfully sample with nutpie, numpyro and blackjax so that works for me as a work-around.

[ricardoV94]

I also tested removing numba from my venv, at which point the default sampler (so pm.sample()) fails on an ImportError for numba.

That's very odd because PyTensor doesn't use numba by default. You seem to have overridden pytensor.config.mode or be passing compile_kwargs to pm.sample

[LucHeuff]

Ah, I had a .pytensorrc hiding in my root setting PyTensor to Numba mode.
Sorry about the confusion.