# NegativeBinomial


Negative Binomial family with log link (NB2 parameterization).


Usage

``` python
NegativeBinomial(theta=1.0)
```


The Negative Binomial family models count data that is overdispersed relative to the Poisson distribution, i.e. where the observed variance exceeds the mean. This commonly arises when counts are driven by unobserved heterogeneity across observations (e.g. some individuals or locations are systematically more prone to events than others). Whittaker uses the NB2 parameterization, where the variance function is `V(mu) = mu + mu^2/theta` and `theta` controls the degree of overdispersion: as `theta -> infinity` the distribution converges to [Poisson](Poisson.md#whittaker.Poisson), while smaller `theta` implies heavier overdispersion. The canonical log link is used, giving the same multiplicative interpretation of coefficients as [Poisson](Poisson.md#whittaker.Poisson).


## Parameters


`theta: float = ``1.0`  
Overdispersion (size) parameter, must be positive. Smaller values of `theta` imply greater overdispersion; larger values make the distribution approach [Poisson](Poisson.md#whittaker.Poisson). The value supplied at construction is used as the starting point and is refined during fitting via an outer iteration around P-IRLS unless explicitly held fixed by the caller.


## Notes

The canonical link is the natural logarithm:

 g(\mu) = \log(\mu) 

The variance function is

 V(\mu) = \mu + \frac{\mu^2}{\theta}, 

so the variance always exceeds the mean by the extra term \mu^2/\theta. The deviance is

 D(y, \hat\mu) = 2 \sum_i \left\[ y_i \log\\\left(\frac{y_i}{\hat\mu_i}\right) - (y_i + \theta) \log\\\left(\frac{y_i + \theta}{\hat\mu_i + \theta}\right) \right\] . 


## Examples

Fit a GAM to overdispersed count data:


``` python
import numpy as np
import whittaker as wk

rng = np.random.default_rng(0)
n = 200
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(0.5 * np.sin(x))
theta = 3.0
y = rng.negative_binomial(theta, theta / (theta + mu))

data = {"x": x, "y": y}

model = wk.GAM("y ~ s(x)", family=wk.NegativeBinomial(theta=theta))
model.fit(data, method="REML")
print(model.summary())
```


    GAM fit summary
    ============================================================
    Formula:    y ~ s(x)
    Family:     NegativeBinomial(theta=4.721, link='log')
    Inference:  REML
    Observations: 200
    Coefficients: 10

    Parametric coefficients:
      Term                       Estimate    Std.Err    z value    p-value
      ------------------------ ---------- ---------- ---------- ----------
      (Intercept)                 -0.0054     0.0792     -0.069     0.9451

    Approximate significance of smooth terms:
      Term                        EDF Ref.df     Chi.sq    p-value
      ------------------------ ------ ------ ---------- ----------
      s(x)                       3.54      4     15.493   0.003781

    Total EDF:  4.54
    Scale est:  1.000000
    Deviance:   210.5567
    Null dev:   229.3245
    Dev. expl:  8.2%
    GCV score:  1.102291
    AIC:        543.06
    BIC:        558.05


## Attributes

| Name | Description |
|----|----|
| [scale_known](#scale_known) | Whether the dispersion parameter is fixed. Always `True` for NegativeBinomial. |
| [theta](#theta) | Overdispersion (size) parameter `theta` of the NB2 distribution. |

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


### scale_known


Whether the dispersion parameter is fixed. Always `True` for NegativeBinomial.


`scale_known: bool`


The NB2 dispersion is governed entirely by `theta` rather than a separate scale parameter estimated during P-IRLS, so the scale is treated as fixed at `1`. `theta` itself is refined by an outer iteration around P-IRLS rather than the usual scale estimation.


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


### theta


Overdispersion (size) parameter `theta` of the NB2 distribution.


`theta: float`


Smaller values imply greater overdispersion relative to Poisson; larger values make the distribution approach Poisson as `theta -> infinity`. This value is used by `variance`, `unit_deviance`, `log_likelihood`, and `simulate`.


## Methods

| Name | Description |
|----|----|
| [deviance()](#deviance) | Total NB2 deviance, the (weighted) sum of `unit_deviance`. |
| [initialize()](#initialize) | Starting values for `mu`: `y` nudged away from zero. |
| [link()](#link) | Apply the log link: \eta = \log(\mu). |
| [link_derivative()](#link_derivative) | Derivative of the log link: g'(\mu) = 1/\mu. |
| [link_inverse()](#link_inverse) | Apply the inverse log link: \mu = e^{\eta}. |
| [log_likelihood()](#log_likelihood) | NB2 log-likelihood evaluated using the current `theta`. |
| [simulate()](#simulate) | Simulate NB2-distributed response values with mean `mu` and the current `theta`. |
| [unit_deviance()](#unit_deviance) | Per-observation NB2 deviance contributions. |
| [variance()](#variance) | NB2 variance function: V(\mu) = \mu + \mu^2/\theta. |

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


### deviance()


Total NB2 deviance, the (weighted) sum of `unit_deviance`.


Usage

``` python
deviance(
    y,
    mu,
    *,
    weights=None,
)
```


#### Parameters


`y: NDArray`  
Observed response values (counts), shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.

`weights: NDArray | None = None`  
Optional prior weights, shape `(n,)`.


#### Returns


`float`  
The total (weighted) deviance.


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


### initialize()


Starting values for `mu`: `y` nudged away from zero.


Usage

``` python
initialize(y)
```


Since the log link requires strictly positive `mu`, small or zero counts are pushed away from zero to avoid `log(0)` on the first P-IRLS iteration.


#### Parameters


`y: NDArray`  
Observed response values (counts), shape `(n,)`.


#### Returns


`NDArray`  
Starting values for `mu`, shape `(n,)`.


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


### link()


Apply the log link: \eta = \log(\mu).


Usage

``` python
link(mu)
```


#### Parameters


`mu: NDArray`  
Conditional mean values, shape `(n,)`. Must be positive.


#### Returns


`NDArray`  
Linear predictor values \eta = \log(\mu), shape `(n,)`.


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


### link_derivative()


Derivative of the log link: g'(\mu) = 1/\mu.


Usage

``` python
link_derivative(mu)
```


#### Parameters


`mu: NDArray`  
Conditional mean values, shape `(n,)`.


#### Returns


`NDArray`  
Derivative values 1/\mu, shape `(n,)`.


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


### link_inverse()


Apply the inverse log link: \mu = e^{\eta}.


Usage

``` python
link_inverse(eta)
```


The linear predictor is clipped to `[-30, 30]` before exponentiating to guard against overflow while fitting.


#### Parameters


`eta: NDArray`  
Linear predictor values, shape `(n,)`.


#### Returns


`NDArray`  
Conditional mean values \mu = e^{\eta}, shape `(n,)`.


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


### log_likelihood()


NB2 log-likelihood evaluated using the current `theta`.


Usage

``` python
log_likelihood(
    y,
    mu,
    scale,
    *,
    weights=None,
)
```


The `scale` argument is accepted for interface compatibility but ignored; the overdispersion is instead governed by `theta` (see `scale_known`).


#### Parameters


`y: NDArray`  
Observed response values (counts), shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.

`scale: float`  
Ignored.

`weights: NDArray | None = None`  
Optional prior weights, shape `(n,)`.


#### Returns


`float`  
The total log-likelihood.


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


### simulate()


Simulate NB2-distributed response values with mean `mu` and the current `theta`.


Usage

``` python
simulate(
    mu,
    scale,
    rng,
)
```


#### Parameters


`mu: NDArray`  
Mean (fitted values), shape `(n,)`.

`scale: float`  
Ignored; overdispersion is governed by `theta`.

`rng: np.random.Generator`  
A `numpy.random.Generator` instance.


#### Returns


`NDArray`  
Simulated response values, shape `(n,)`.


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


### unit_deviance()


Per-observation NB2 deviance contributions.


Usage

``` python
unit_deviance(
    y,
    mu,
)
```


Computes d_i = 2 \left\[ y_i \log(y_i/\hat\mu_i) - (y_i + \theta) \log\\\left(\frac{y_i+\theta}{\hat\mu_i+\theta}\right) \right\], with the usual convention that the y_i \log(\cdot) term vanishes when y_i = 0.


#### Parameters


`y: NDArray`  
Observed response values (counts), shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.


#### Returns


`NDArray`  
Per-observation deviance contributions, shape `(n,)`.


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


### variance()


NB2 variance function: V(\mu) = \mu + \mu^2/\theta.


Usage

``` python
variance(mu)
```


The variance always exceeds the mean by the extra term \mu^2/\theta, which vanishes as `theta -> infinity` (recovering the Poisson variance `V(mu) = mu`).


#### Parameters


`mu: NDArray`  
Conditional mean values, shape `(n,)`.


#### Returns


`NDArray`  
Variance values \mu + \mu^2/\theta, shape `(n,)`.