Truncation errors (falling back to rejection sampling despite logcdf present)

I’ve been working with a truncated BetaBinomial distribution, and while trying to simulate data (via pm.draw() or pm.sample_prior_predictive), I often encounter Truncation Errors (TruncationError: Truncation did not converge in 10000 steps). From what I see in the documentation and the code behind pm.Truncated, this error should only arise if rejection sampling is used, and it should only be used when the truncated distribution does not implement a logcdf function. But pm.BetaBinomial has logcdf method, and so I am not sure why this is happening. Am I doing something wrong and are there ways to avoid / minimize the risk of such errors?

Sample code:

from pymc.distributions.mixture import _hurdle_mixture
class HurdleBetaBinomial:

    def __new__(cls, name, psi, n, alpha, beta, **kwargs):
        return _hurdle_mixture(
            name=name, nonzero_p=psi, nonzero_dist=pm.BetaBinomial.dist(n=n, alpha=alpha, beta=beta), dtype="int", **kwargs
        )


    @classmethod
    def dist(cls, psi, n, alpha, beta, **kwargs):
        return _hurdle_mixture(
            name=None, nonzero_p=psi, nonzero_dist=pm.BetaBinomial.dist(n=n, alpha=alpha, beta=beta), dtype="int", **kwargs
        )

N = 1
with pm.Model() as m: 
    treatment = pm.Bernoulli("treatment", p=0.5, shape=N)
    baseline = pm.Normal("inv_p", 0.2, 1, shape=N)        
    β = pm.Gamma("β", alpha=1, beta=4, shape=N)
    α = pm.Gamma("α", alpha=1, beta=2, shape=N)   

    ψ_treatment = pm.Normal("inv_ψ_effect", 0, 0.2)
    α_treatment = pm.Normal("α_effect", 0, 0.02)
    β_treatment = pm.Normal("β_effect", 0, 0.02)

    ψ = pm.Deterministic("ψ", pm.invlogit(baseline + ψ_treatment * treatment))
    
    days = HurdleBetaBinomial(
        "days", 
        n=7, 
        psi=ψ, 
        alpha=α + α_treatment * treatment, 
        beta=β + β_treatment * treatment, 
        shape=N, dims="tracks"
    )    
sim_model = pm.do(m, {
    "inv_ψ_effect": 0,
    "α_effect": 0.1,
    "β_effect": 0
})

#try 10 times
for i in range(10):
    with sim_model:
        simulate = pm.sample_prior_predictive(samples=10_000)

Logcdf isn’t enough, icdf is also needed to avoid rejection sampling. According to this issue, that method is still missing for BetaBinomial: Add icdf functions for distributions · Issue #6612 · pymc-devs/pymc · GitHub

In the meantime you can try to increase the number of steps for convergence, which is called max_n_steps

1 Like

Got it, thank you!