# ZeroInflatedPoisson


Zero-inflated Poisson (ZIP) family for GAMLSS.


Usage

``` python
ZeroInflatedPoisson()
```


Count data often has more zeros than a plain [Poisson](Poisson.md#whittaker.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](ZeroInflatedPoisson.md#whittaker.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](GAMLSS.md#whittaker.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](ZeroInflatedNegativeBinomial.md#whittaker.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:


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

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](#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()](#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 Poisson model. |
| [simulate()](#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.


Usage

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


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


Usage

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


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


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]`, giving [GAMLSS](GAMLSS.md#whittaker.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.


Usage

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


Applies the parameter-specific link function used internally by [GAMLSS](GAMLSS.md#whittaker.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.


Usage

``` python
link_derivative(
    param,
    values,
)
```


Used by the IRLS fitting routine in [GAMLSS](GAMLSS.md#whittaker.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](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 Poisson model.


Usage

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


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