# ZeroInflatedNegativeBinomial


Zero-inflated negative binomial (ZINB) family for GAMLSS.


Usage

``` python
ZeroInflatedNegativeBinomial(theta=1.0)
```


[ZeroInflatedNegativeBinomial](ZeroInflatedNegativeBinomial.md#whittaker.ZeroInflatedNegativeBinomial) combines the two departures from plain [Poisson](Poisson.md#whittaker.Poisson) counts that are most common in practice: overdispersion (variance exceeding the mean, as in [NegativeBinomial](NegativeBinomial.md#whittaker.NegativeBinomial)) and structural excess zeros (as in [ZeroInflatedPoisson](ZeroInflatedPoisson.md#whittaker.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](GAMLSS.md#whittaker.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](ZeroInflatedPoisson.md#whittaker.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:


``` python
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](#parameter_names) | Names of the distributional parameters modeled by this family. |
| [theta](#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.


`theta: float`


## Methods

| Name | Description |
|----|----|
| [d2l_dtheta2()](#d2l_dtheta2) | Compute the (negative) second derivative of the log-likelihood. |
| [dl_dtheta()](#dl_dtheta) | Compute the first derivative of the log-likelihood with respect to a parameter. |
| [initialize()](#initialize) | Generate starting values for `mu` and `pi` before the first IRLS iteration. |
| [link()](#link) | Map a distributional parameter from its natural scale to the link scale. |
| [link_derivative()](#link_derivative) | Compute the derivative of the link function with respect to the natural parameter. |
| [link_inverse()](#link_inverse) | Map a distributional parameter from the link scale back to its natural scale. |
| [log_likelihood()](#log_likelihood) | Compute the total log-likelihood under the zero-inflated negative binomial model. |
| [simulate()](#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.


Usage

``` python
d2l_dtheta2(
    param,
    y,
    params,
)
```


Implements the family-specific working weight used by [GAMLSS](GAMLSS.md#whittaker.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.


Usage

``` python
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](ZeroInflatedPoisson.md#whittaker.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.


Usage

``` python
initialize(y)
```


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.


Usage

``` python
link(
    param,
    values,
)
```


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.


Usage

``` python
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](GAMLSS.md#whittaker.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](TermsPredictionResult.md#whittaker.TermsPredictionResult.values).


------------------------------------------------------------------------


### link_inverse()


Map a distributional parameter from the link scale back to its natural scale.


Usage

``` python
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.


Usage

``` python
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.


Usage

``` python
simulate(
    params,
    rng,
)
```


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`.