BUG: numba backend samples Cauchy/HalfCauchy with the wrong loc and scale

[notluquis]

Description

pytensor/link/numba/dispatch/random.py:

@numba_core_rv_funcify.register(ptr.CauchyRV)
def numba_core_CauchyRV(op, node):
    @numba_basic.numba_njit
    def random(rng, loc, scale):
        return (loc + rng.standard_cauchy()) / scale

A location-scale draw must be loc + scale * z. Written this way the sampled variable has location loc/scale and scale 1/scale.

rng_fn_scipy is correct, and JAX routes CauchyRV through the generic jax_sample_fn_loc_scale. Only numba is wrong — and it is the default linker. logp is unaffected, so NUTS posteriors are fine, but everything that draws is corrupted: pymc.draw, sample_prior_predictive, sample_posterior_predictive. pm.HalfCauchy inherits it via abs(cauchy(0, beta)) — a prior-predictive check on a HalfCauchy prior is silently wrong, which is how we hit it.

Reproduce

import numpy as np, pytensor, scipy.stats as st
from pytensor.tensor.random.basic import cauchy
from pytensor.compile.mode import Mode

def iqr(x): return float(np.subtract(*np.percentile(np.asarray(x), [75, 25])))
g = cauchy(3.0, 5.0, size=200_000)

print("numba (default):", np.median(g.eval()), iqr(g.eval()))
print("py linker      :", np.median(pytensor.function([], g, mode=Mode(linker="py", optimizer=None))()))
print("scipy          :", st.cauchy(3, 5).median(), st.cauchy(3, 5).ppf(.75) - st.cauchy(3, 5).ppf(.25))

numba (default): 0.5998   0.4002     <-- loc/scale = 3/5, 1/scale = 1/5
py linker      : 2.9871  10.0272
scipy          : 3.0     10.0

Of ten continuous distributions checked the same way (normal, gamma, exponential, laplace, logistic, t, vonmises, pareto, weibull), only cauchy is affected.

Why it survived

The expression was written twice — 5ea74e973 (2022-03-03, Aesara, as a string template) and 14da898c6 (2024-05-10, the Generator port, which rewrote the function and carried it forward verbatim). Both are spanned by the single test case, which passes scale=1.0:

(ptr.cauchy,
 [(pt.dvector(), np.array([1.0, 2.0])),   # loc
  (pt.dscalar(), np.array(1.0))],         # scale
 (2,), "cauchy", lambda *args: args),

At scale = 1 the buggy expression is algebraically identical to the correct one, so both errors vanish together and the Cramér-von Mises check passes:

loc	scale	median	IQR	scipy	CvM p
1.0	1.0	0.992	2.005	1.0 / 2.0	9.2e-01	passes
1.0	2.0	0.506	1.007	1.0 / 4.0	1.0e-08	fails
3.0	5.0	0.599	0.399	3.0 / 10.0	3.4e-08	fails

Any scale != 1 in the existing parametrisation catches it.

Fix

-        return (loc + rng.standard_cauchy()) / scale
+        return loc + scale * rng.standard_cauchy()

[ricardoV94]

yeah I hate tests with scale=1, loc=0 precisely because of that. Want to open a PR?

[notluquis]

PR ready. #2309