RuntimeWarning: overflow encountered in _beta_ppf

[twiecki]

Code:

import numpy as np
import arviz as a
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

# PyMC Imports
import pymc as pm4 # PyMC 4.0

# Aesara and Theano imports
import aesara.tensor as at # used by PyMC 4.0
import aesara

data = pd.read_csv(pm4.get_data("radon.csv"))
county_names = data.county.unique()

data["log_radon"] = data["log_radon"].astype(aasara.config.floatX)

county_idx, counties = pd.factorize(data.county)
coords = {
    "county": counties,
    "obs_id": np.arange(len(county_idx)),
}

def build_model(pm, **kwargs):
    with pm.Model(coords=coords) as hierarchical_model:
        mu_a = pm.Normal("mu_a", mu=0.0, sigma=10)
        sigma_a = pm.HalfNormal("sigma_a", 5.0)
        mu_b = pm.Normal("mu_b", mu=0.0, sigma=10)
        sigma_b = pm.HalfNormal("sigma_b", 5.0)
        a = pm.Normal("a", dims="county") * sigma_a + mu_a
        b = pm.Normal("b", dims="county") * sigma_b + mu_b
        eps = pm.HalfNormal("eps", 5.0)
        radon_est = a[county_idx] + b[county_idx] * data.floor.values
        radon_like = pm.Normal(
            "radon_like", mu=radon_est, sigma=eps, observed=data.log_radon, 
            dims="obs_id"
        )
        
    return hierarchical_model

model_pymc4 = build_model(pm4)
with model_pymc4:
    idata_pymc4 = pm4.sample()

Samples fine but displays warning:

/Users/twiecki/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf
  return _boost._beta_ppf(q, a, b)

[ricardoV94]

Can you narrow when is the scipy function being called?

[twiecki]

I tried grepping for beta_ppf in the pymc code base but couldn't find anything. Any other ideas?

[ricardoV94]

I imagine it comes from Aesara, some Op that calls a scipy function that calls that. Maybe a RandomVariable by the looks of it.

[ricardoV94]

Do you get it with compute_initial_point?

[twiecki]

Only at the end of sampling.

[ricardoV94]

Hmm perhaps when computing the log-likelihood then

[twiecki]

beta_ppf is also not referenced in aesara anywhere.

[twiecki]

File ~/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/pymc/step_methods/step_sizes.py:73, in DualAverageAdaptation.warnings(self)
     71 n_bound = min(100, len(accept))
     72 n_good, n_bad = mean_accept * n_bound, (1 - mean_accept) * n_bound
---> 73 lower, upper = stats.beta(n_good + 1, n_bad + 1).interval(0.95)
     74 if target_accept < lower or target_accept > upper:
     75     msg = (
     76         f"The acceptance probability does not match the target. "
     77         f"It is {mean_accept:0.4g}, but should be close to {target_accept:0.4g}. "
     78         f"Try to increase the number of tuning steps."
     79     )

File ~/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:511, in rv_frozen.interval(self, alpha)
    510 def interval(self, alpha):
--> 511     return self.dist.interval(alpha, *self.args, **self.kwds)

File ~/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:1429, in rv_generic.interval(self, alpha, *args, **kwds)
   1427 q1 = (1.0-alpha)/2
   1428 q2 = (1.0+alpha)/2
-> 1429 a = self.ppf(q1, *args, **kwds)
   1430 b = self.ppf(q2, *args, **kwds)
   1431 return a, b

File ~/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py:2136, in rv_continuous.ppf(self, q, *args, **kwds)
   2134     goodargs = argsreduce(cond, *((q,)+args+(scale, loc)))
   2135     scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
-> 2136     place(output, cond, self._ppf(*goodargs) * scale + loc)
   2137 if output.ndim == 0:
   2138     return output[()]

File ~/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py:624, in beta_gen._ppf(self, q, a, b)
    623 def _ppf(self, q, a, b):
--> 624     return _boost._beta_ppf(q, a, b)

RuntimeWarning: overflow encountered in _beta_ppf

Some debugging:

> /Users/twiecki/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_continuous_distns.py(624)_ppf()
    622 
    623     def _ppf(self, q, a, b):
--> 624         return _boost._beta_ppf(q, a, b)
    625 
    626     def _stats(self, a, b):

ipdb> print(q)
[0.03]
ipdb> print(a)
[82.64]
ipdb> print(b)
[19.36]
ipdb> up
> /Users/twiecki/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py(2136)ppf()
   2134             goodargs = argsreduce(cond, *((q,)+args+(scale, loc)))
   2135             scale, loc, goodargs = goodargs[-2], goodargs[-1], goodargs[:-2]
-> 2136             place(output, cond, self._ppf(*goodargs) * scale + loc)
   2137         if output.ndim == 0:
   2138             return output[()]

ipdb> up
> /Users/twiecki/miniforge3/envs/pymc4b5/lib/python3.10/site-packages/scipy/stats/_distn_infrastructure.py(1429)interval()
   1427         q1 = (1.0-alpha)/2
   1428         q2 = (1.0+alpha)/2
-> 1429         a = self.ppf(q1, *args, **kwds)
   1430         b = self.ppf(q2, *args, **kwds)
   1431         return a, b

ipdb> q1
0.025000000000000022
ipdb> *args
*** SyntaxError: can't use starred expression here
ipdb> args
self = <scipy.stats._continuous_distns.beta_gen object at 0x176690af0>
alpha = array(0.95)
args = (82.63918461537455, 19.36081538462545)
kwds = {}

[fbarfi]

I get this error when using pm.sample with pymc3.11
The problem disappeared when I shifted to pymc and using pymc.sampling_jax.sample_numpyro_nuts

[ricardoV94]

This is not an error, it's a warning we get from calling a scipy function

[fbarfi]

I know I misspoke, sorry. It is just annoying. thanks.

[ricardoV94]

We need to check if the warning is harmless, in which case we can suppress it

[fbarfi]

Based on my limited knowledge, I did not notice any difference in the output as far the sampling is concerned (that is, when I get the warning and when I do not get it). thanks.

[twiecki]

Known scipy issue: scipy/scipy#14901

[fonnesbeck]

I'm seeing this now when running the multilevel notebook. Strange that they are just appearing now. I guess we can just filter them when they appear.

[twiecki]

I didn't manage to filter them.