The Zero Trick and Custom Likelihoods*

import arviz as az
import pymc as pm
from pymc.math import log, sqr
import numpy as np

rng = np.random.default_rng(6420)

10. The Zero Trick and Custom Likelihoods*#

Zero-trick Jeremy#

This introduces the “zero trick”, which is a method for specifying custom likelihoods in BUGS. For a more detailed treatment of these methods, see Ntzoufras [2009], page 276, which is where I got this explanation. These tricks are unnecessary in PyMC since we can just define custom distributions directly, but they do seem to work.

Adapted from Unit 6: zerotrickjeremy.odc.

Here’s the model we’re using:

\[\begin{split} \begin{align*} X_i \mid \theta &\sim N(\theta, 80) \\ \theta &\sim N(110, 120) \end{align*} \end{split}\]

Of course, BUGS can handle this model just fine. But, let’s pretend for a moment that there is no built-in normal distribution. The zero trick takes advantages of some properties of the Poisson distribution to recreate an arbitrary likelihood.

Given a log-likelihood of the form \(l_i = \log f(y; \theta)\),

\[\begin{split} \begin{align*} f(y | \theta) &= \prod_{i=1}^{n} e^{l_i} && \text{Log-likelihood to likelihood.}\\ &= \prod_{i=1}^{n} \frac{e^{-(-l_i)}(-l_i)^0}{0!} && \text{Matching the form of the Poisson PMF.} \\ &= \prod_{i=1}^{n} f_P(0; -l_i) && \text{Poisson evaluated at zero with mean $-l_i$.} \end{align*} \end{split}\]

–{cite:t} ntzoufras2009bayesian page 276.

But the rate, \(\lambda\), can’t be negative. So we need to add a constant, C, to keep that from happening.

The normal log-likelihood is:

\[l_i = -0.5 \log(2\pi) - \log(\sigma) - \frac{(y_i - \mu_i)^2}{2\sigma^2}\]

But the constant terms won’t affect the posterior.

Here’s the model in PyMC:

y = 98
μ = 110
σ = 80**0.5
τ = 120**0.5
C = 10000

inits = {"θ": 100}

with pm.Model() as m:
    θ = pm.Flat("θ")

    λ1 = pm.Deterministic("λ1", log(σ) + 0.5 * sqr(((y - θ) / σ)) + C)
    λ2 = pm.Deterministic("λ2", log(τ) + 0.5 * sqr(((θ - μ) / τ)) + C)

    z1 = pm.Poisson("z1", λ1, observed=0)
    z2 = pm.Poisson("z2", λ2, observed=0)

    trace = pm.sample(5000, tune=1000, initvals=inits, target_accept=0.88)

Hide code cell output

 
NUTS[nutpie]: [θ]













az.summary(trace, ci_prob=0.95, var_names="θ", kind="stats")










  
    

      
      mean
      sd
      eti95_lb
      eti95_ub
    

  
  
    

      θ
      100
      6.9
      89
      120
    

  









Custom likelihoods in PyMC#


PyMC is a lot more flexible. For one, many more built-in distributions are available, including mixtures and zero-inflated models, which is where I’ve seen the zero trick used most often (Karlis and Ntzoufras [2008], Ntzoufras [2009] page 288).


If you need a custom likelihood, take a look at the pm.Potential or pm.CustomDist classes. pm.Potential is used to adjust the likelihood, as we can see in the censored model ahead in Unit 8.6. For now, let’s take a look at how to use pm.CustomDist.




pm.CustomDist Motivating Example: Survival Model#


We can alternatively build the survival model in in Unit 6.9 with the custom distribution based on this post (but with equations corrected):



\[\begin{split}f(c,t) = \begin{cases} 
      \lambda \exp(-\lambda t) & \text{if } c = 1 \\
      \exp(-\lambda t) & \text{if } c = 0 \\
   \end{cases}\end{split}\]


c is an indicator variable. The data point is censored if c=0, and a failure or non-survival is denoted by c=1. When the point is not censored, the probability distribution is exponential. When the point is censored, the probability is the exponential distribution’s complementary CDF, which accounts for the probability to the right of the censored point.


To build this as a custom distribution, we create a logp or log-probability function, as PyMC prefers to work with log-probabilities. One tricky part for this distribution is that both observed data arrays (t and c) need to be concatenated into a matrix as input into the observed argument. We can see the summary results match those on Unit 6.9, which uses the different pm.Censored function.






# gamma dist parameters
α = 0.01
β = 0.1

# observed life within experiment dates
t = np.array([2, 72, 51, 60, 33, 27, 14, 24, 4, 21])
# censored indicator
c = np.array([1, 0, 1, 0, 1, 1, 1, 1, 1, 0])
# CustomDist requires all data to be together in a matrix
val = np.concat([[t], [c]])


# log-probability of custom distribution
def logp(value, lam):
    t = value[0, :]
    c = value[1, :]
    return (c * log(lam) - lam * t).sum()


with pm.Model() as m:
    λ = pm.Gamma("λ", α, β)
    μ = pm.Deterministic("μ", 1 / λ)

    exp_surv = pm.CustomDist("exp_surv", λ, logp=logp, observed=val)

    trace = pm.sample(5000, random_seed=rng)

az.summary(trace)









NUTS[nutpie]: [λ]








  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      λ
      0.0228
      0.0086
      0.011
      0.038
      8529
      8802
      1.00
      8.9e-05
      7.5e-05
    

    

      μ
      51.2
      23
      26
      92
      8529
      8802
      1.00
      0.26
      0.4
    

  









%load_ext watermark
%watermark -n -u -v -iv -p pytensor









Last updated: Sun, 18 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz: 1.3.0
numpy: 2.4.6
pymc : 6.3.1