Zero-inflated Poisson (ZIP) family for GAMLSS.
Count data often has more zeros than a plain Poisson model can explain — for example, when some observations are structurally incapable of the event occurring at all (a “never-taker” always reports zero), in addition to the zeros that arise simply because the Poisson mean is low. ZeroInflatedPoisson models this as a mixture: with probability pi an observation is a structural zero, and with probability 1 - pi it is drawn from an ordinary Poisson(mu) distribution (which can itself still produce a zero). Both mu and pi are modeled as smooth functions of covariates through GAMLSS, so the excess-zero probability and the count intensity can each vary independently across the covariate space. If overdispersion remains even among the non-structural-zero counts, use ZeroInflatedNegativeBinomial instead.
Notes
Two distributional parameters are modeled, each with its own link:
g_{\mu}(\mu) = \log(\mu), \qquad g_{\pi}(\pi) = \log\!\left(\frac{\pi}{1-\pi}\right).
The probability mass function is a mixture of a point mass at zero and a Poisson distribution:
P(Y = 0) = \pi + (1-\pi) e^{-\mu}, \qquad
P(Y = k) = (1-\pi) \frac{\mu^{k} e^{-\mu}}{k!} \quad \text{for } k > 0.
Examples
Fit a GAMLSS for count data with excess zeros:
import numpy as np
import whittaker as wk
from scipy.special import expit
rng = np.random.default_rng(0)
n = 300
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(0.5 + 0.4 * np.sin(x))
pi = expit(-1.0 + 0.8 * np.cos(x))
is_structural_zero = rng.uniform(size=n) < pi
counts = rng.poisson(mu)
y = np.where(is_structural_zero, 0, counts).astype(float)
data = {"x": x, "y": y}
model = wk.GAMLSS(
formulas={"mu": "y ~ s(x)", "pi": "y ~ s(x)"},
family=wk.ZeroInflatedPoisson(),
)
model.fit(data)
print(model.summary())
GAMLSS fit summary
========================================
Family: ZeroInflatedPoisson(mu=log, pi=logit)
N obs: 300
Global deviance: 886.8940
AIC: 896.9039
BIC: 915.4410
Log-likelihood: -443.4470
Converged: True (6 iterations)
--- mu ---
EDF total: 2.00
Smooth 1: edf = 1.00
--- pi ---
EDF total: 3.00
Smooth 1: edf = 2.00
Attributes
|
Name
|
Description
|
|
parameter_names
|
Names of the distributional parameters modeled by this family.
|
parameter_names
Names of the distributional parameters modeled by this family.
parameter_names: tuple[str, …]
Methods
|
Name
|
Description
|
|
d2l_dtheta2()
|
Compute the (negative) second derivative of the log-likelihood.
|
|
dl_dtheta()
|
Compute the first derivative of the log-likelihood with respect to a parameter.
|
|
initialize()
|
Generate starting values for mu and pi before the first IRLS iteration.
|
|
link()
|
Map a distributional parameter from its natural scale to the link scale.
|
|
link_derivative()
|
Compute the derivative of the link function with respect to the natural parameter.
|
|
link_inverse()
|
Map a distributional parameter from the link scale back to its natural scale.
|
|
log_likelihood()
|
Compute the total log-likelihood under the zero-inflated Poisson model.
|
|
simulate()
|
Draw a random sample of counts from the fitted zero-inflated Poisson model.
|
d2l_dtheta2()
Compute the (negative) second derivative of the log-likelihood.
d2l_dtheta2(
param,
y,
params,
)
Implements the family-specific working weight used by GAMLSS’s Fisher-scoring updates. The returned values approximate -\partial^2 \ell / \partial \theta^2 for theta in {"mu", "pi"}, computed separately for zero and positive observations to reflect the zero-inflated mixture, and floored at _EPS to keep weights strictly positive.
Parameters
param: str-
Name of the parameter, either "mu" or "pi".
y: NDArray-
Observed response values.
params: dict of str to NDArray-
Current fitted values of all distributional parameters (
"mu" and "pi").
Returns
NDArray-
Elementwise working weight for
param.
dl_dtheta()
Compute the first derivative of the log-likelihood with respect to a parameter.
dl_dtheta(
param,
y,
params,
)
Implements the family-specific score function used by GAMLSS to build the working response at each IRLS iteration. For zero observations the derivative accounts for the mixture weight of the point mass versus the Poisson component at zero; for positive observations it reduces to the ordinary Poisson score. For mu this is \partial \ell / \partial \mu, and for pi it is \partial \ell / \partial \pi, both derived from the zero-inflated Poisson log-likelihood.
Parameters
param: str-
Name of the parameter to differentiate with respect to, either "mu" or "pi".
y: NDArray-
Observed response values.
params: dict of str to NDArray-
Current fitted values of all distributional parameters (
"mu" and "pi").
Returns
NDArray-
Elementwise first derivative of the log-likelihood with respect to
param.
initialize()
Generate starting values for mu and pi before the first IRLS iteration.
Sets mu to the mean of the strictly positive observations (or 1.0 if none are positive) for every observation, and sets pi to half the observed fraction of zeros, clipped to [0.01, 0.5], giving GAMLSS a reasonable starting point without any structural information about the covariates.
Parameters
y: NDArray-
Observed response values.
Returns
dict of str to NDArray-
Initial values for
"mu" and "pi", each broadcast to the shape of y.
link()
Map a distributional parameter from its natural scale to the link scale.
Applies the parameter-specific link function used internally by GAMLSS: the log link \eta = \log(\mu) for mu (clamped away from zero via _MU_FLOOR), and the logit link \eta = \log(\pi / (1-\pi)) for pi (clamped away from 0 and 1).
Parameters
param: str-
Name of the parameter to transform, either "mu" or "pi".
values: NDArray-
Values of
param on its natural scale.
Returns
NDArray-
Transformed values on the link (linear predictor) scale.
link_derivative()
Compute the derivative of the link function with respect to the natural parameter.
link_derivative(
param,
values,
)
Used by the IRLS fitting routine in GAMLSS to convert working responses between the link and natural scales. Returns d\eta/d\mu = 1/\mu for mu and d\eta/d\pi = 1/(\pi(1-\pi)) for pi, both evaluated with the same clamping used in link.
Parameters
param: str-
Name of the parameter, either "mu" or "pi".
values: NDArray-
Values of
param on its natural scale.
Returns
NDArray-
Derivative of the link function evaluated at values.
link_inverse()
Map a distributional parameter from the link scale back to its natural scale.
link_inverse(
param,
eta,
)
Implements the inverse of link: \mu = e^{\eta} (clipped to [-30, 30] before exponentiating to avoid overflow) for mu, and \pi = \mathrm{expit}(\eta) for pi.
Parameters
param: str-
Name of the parameter to transform, either "mu" or "pi".
eta: NDArray-
Linear predictor values on the link scale.
Returns
NDArray-
Values of
param on its natural scale.
log_likelihood()
Compute the total log-likelihood under the zero-inflated Poisson model.
log_likelihood(
y,
params,
)
Evaluates
\ell = \sum_{i: y_i = 0} \log\!\big(\pi_i + (1-\pi_i) e^{-\mu_i}\big)
+ \sum_{i: y_i > 0} \Big[\log(1-\pi_i) + y_i \log(\mu_i) - \mu_i - \log(y_i!)\Big],
summing the point-mass contribution for structural zeros with the Poisson contribution for observed counts.
Parameters
y: NDArray-
Observed response values.
params: dict of str to NDArray-
Fitted values of
"mu" and "pi" for each observation.
Returns
float-
Total log-likelihood summed over all observations.
simulate()
Draw a random sample of counts from the fitted zero-inflated Poisson model.
For each observation, flips a pi-weighted coin to decide whether it is a structural zero; observations that are not structural zeros are drawn from Poisson(mu), so a sampled value can still be zero even when it was not selected as a structural zero.
Parameters
params: dict of str to NDArray-
Fitted values of "mu" and "pi" for each observation to simulate.
rng: object-
Random number generator exposing
uniform and poisson methods, such as a numpy.random.Generator.
Returns
NDArray-
Simulated response values, one per row of
params.