Zero-inflated negative binomial (ZINB) family for GAMLSS.
ZeroInflatedNegativeBinomial(theta=1.0)
ZeroInflatedNegativeBinomial combines the two departures from plain Poisson counts that are most common in practice: overdispersion (variance exceeding the mean, as in NegativeBinomial) and structural excess zeros (as in ZeroInflatedPoisson). It is appropriate for count data where, even after allowing for a mixture of structural and sampling zeros, the remaining positive counts are still more variable than a Poisson model would predict — for example, healthcare utilization counts, insurance claim frequencies, or ecological abundance data with many true absences plus overdispersed non-zero counts. mu (the NB mean) and pi (the zero-inflation probability) are modeled as smooth functions of covariates through GAMLSS, while the overdispersion parameter theta is fixed at construction rather than estimated per observation.
Parameters
theta: float = 1.0-
Negative binomial size (overdispersion) parameter, must be positive. Larger values mean less overdispersion (
theta -> infinity recovers ZeroInflatedPoisson); smaller values mean more overdispersion among the non-structural-zero counts. Unlike mu and pi, theta is a single fixed value shared across all observations rather than modeled by a smooth predictor.
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 Negative Binomial distribution using the mean-size parameterization (Var(NB) = mu + mu^2/theta):
P(Y = 0) = \pi + (1-\pi)\, \mathrm{NB}(0 \mid \mu, \theta), \qquad
P(Y = k) = (1-\pi)\, \mathrm{NB}(k \mid \mu, \theta) \quad \text{for } k > 0,
where
\mathrm{NB}(k \mid \mu, \theta) = \binom{k+\theta-1}{k}
\left(\frac{\theta}{\theta+\mu}\right)^{\theta} \left(\frac{\mu}{\theta+\mu}\right)^{k}.
Examples
Fit a GAMLSS for overdispersed 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))
theta = 2.0
is_structural_zero = rng.uniform(size=n) < pi
counts = rng.negative_binomial(theta, theta / (theta + 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.ZeroInflatedNegativeBinomial(theta=theta),
)
model.fit(data)
print(model.summary())
GAMLSS fit summary
========================================
Family: ZeroInflatedNegativeBinomial(theta=2, mu=log, pi=logit)
N obs: 300
Global deviance: 911.2129
AIC: 921.0281
BIC: 939.2047
Log-likelihood: -455.6064
Converged: True (3 iterations)
--- mu ---
EDF total: 2.00
Smooth 1: edf = 1.00
--- pi ---
EDF total: 2.91
Smooth 1: edf = 1.91
Attributes
|
Name
|
Description
|
|
parameter_names
|
Names of the distributional parameters modeled by this family.
|
|
theta
|
Fixed negative binomial size (overdispersion) parameter.
|
parameter_names
Names of the distributional parameters modeled by this family.
parameter_names: tuple[str, …]
theta
Fixed negative binomial size (overdispersion) parameter.
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 negative binomial model.
|
|
simulate()
|
Draw a random sample of counts from the fitted zero-inflated NB 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 for the zero-inflated negative binomial model, computed separately for zero and positive observations 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
"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 for the zero-inflated negative binomial model, combining the derivative of the NB(0 | mu, theta) probability weighted by the mixture for zero observations with the ordinary negative binomial score \partial \ell / \partial \mu = (y - \mu) / (\mu (1 + \mu/\theta)) for positive observations. For pi the score reflects the same zero/positive split as in ZeroInflatedPoisson.dl_dtheta.
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
"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]. The fixed overdispersion parameter theta is not part of the returned dictionary since it is not estimated per observation.
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 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). The overdispersion parameter theta is fixed and has no link.
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,
)
Returns d\eta/d\mu = 1/\mu for mu and d\eta/d\pi = 1/(\pi(1-\pi)) for pi, using the same clamping applied in link, for use by GAMLSS’s IRLS fitting routine.
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 negative binomial model.
log_likelihood(
y,
params,
)
Evaluates
\ell = \sum_{i: y_i = 0}
\log\!\big(\pi_i + (1-\pi_i)\, \mathrm{NB}(0 \mid \mu_i, \theta)\big)
+ \sum_{i: y_i > 0} \Big[\log(1-\pi_i) + \log \mathrm{NB}(y_i \mid \mu_i, \theta)\Big],
where \mathrm{NB}(\cdot \mid \mu, \theta) is the negative binomial probability mass function in the mean-size parameterization with the fixed theta.
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 NB 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 a negative binomial distribution with the fixed theta and success probability theta / (mu + theta), 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 negative_binomial methods, such as a numpy.random.Generator.
Returns
NDArray-
Simulated response values, one per row of
params.