chisquared entropy/logpdf lose float64 precision — pt.log(2) folds at float32

[maresb]

Summary

chisquared.entropy and chisquared.logpdf use pt.log(2) with a Python int literal. PyTensor stores 2 as int8, and log upgrades int8/int16 inputs to float32, so the constant folds to 0.6931471824645996 instead of the float64 log(2) = 0.6931471805599453. That ~1.9e-9 error leaks into both functions.

Affected code

• chisquared.py:54 (entropy): ... + pt.log(2) + ...
• chisquared.py:90 (logpdf): ... - (nu * pt.log(2)) / 2

Reproduction

import numpy as np, pytensor.tensor as pt
from scipy import stats
from pytensor_distributions import chisquared as C

print(float(pt.log(2).eval()))          # 0.6931471824645996  (float64 log(2) = 0.6931471805599453)

for nu in [2.0, 10.0, 100.0]:
    e  = float(C.entropy(pt.constant(nu, "float64")).eval())
    lp = float(C.logpdf(pt.constant(nu, "float64"), pt.constant(nu, "float64")).eval())
    print(nu, abs(e - stats.chi2(nu).entropy()), abs(lp - stats.chi2(nu).logpdf(nu)))

Measured error vs scipy:

nu	entropy	logpdf(x=nu)
2	1.9e-9	1.9e-9
10	1.9e-9	9.5e-9
100	1.9e-9	9.5e-8

entropy is a flat ~1.9e-9; logpdf's scales with nu (the nu·log(2) term), so it grows with the degrees of freedom and accumulates across a summed log-likelihood.

Root cause (upstream, not really a chisquared bug)

This is the "aggressive downcast of literals" root cause (pymc-devs/pytensor#1073) hitting direct user code: pt.log(<int>) → int8 → float32. A maintainer has confirmed #1073 is the underlying issue (comment):

the underlying issue is the same old aggressive downcast of literals: #1073

The obvious pt.log(2.0) does not work either: the Python float 2.0 autocasts to float32, so it's still 0.6931471824645996. Only a NumPy float64 constant is correct.

The real fix belongs upstream (#1073) — once PyTensor stops downcasting Python literals, pt.log(2) will be exact and nothing here needs changing. That root fix (pymc-devs/pytensor#1917) is still an open draft, so a stopgap may be warranted in the meantime.

Workaround (stopgap only — not the real fix)

Until #1073 lands, sidestep the downcast by feeding log a NumPy float64 constant, matching the pattern already used in this repo (genpareto.median, and the merged asymmetriclaplace.median fix in #60):

import numpy as np
# entropy:  ... + np.log(2.0) + ...
# logpdf:   ... - (nu * np.log(2.0)) / 2

This is a local workaround, not a fix for the underlying downcast; it should be revisited (and ideally reverted) once #1073 is resolved.

Scope

The same pt.log(<int>) pattern appears in ~9 modules (inversegamma, moyal, wishart, halfcauchy, weibull, halfstudentt, skew_studentt, polyagamma, plus chisquared) — all the same upstream symptom, all sidesteppable the same way pending #1073. Worth a repo-wide sweep if a stopgap is desired.