# Tweedie


Tweedie family with log link.


Usage

``` python
Tweedie(p=1.5)
```


The Tweedie family is an exponential dispersion model whose variance function is a power of the mean, `V(mu) = mu^p`. It is most widely used for `1 < p < 2`, the compound Poisson-Gamma case: a distribution with a point mass at zero and a continuous, right-skewed density on the positive reals. This shape is a natural fit for aggregate insurance claims (many policyholders have zero claims, and the rest have positive, Gamma-like claim amounts), precipitation totals, and biomass or catch data with structural zeros. Whittaker also supports `p > 2` for purely positive, heavy-tailed continuous data (including the inverse Gaussian case at `p = 3`; see [InverseGaussian](InverseGaussian.md#whittaker.InverseGaussian)). The Poisson (`p = 1`) and Gamma (`p = 2`) boundary cases are excluded here -- use [Poisson](Poisson.md#whittaker.Poisson) or [Gamma](Gamma.md#whittaker.Gamma) directly for exact likelihood computations at those values.


## Parameters


`p: float = ``1.5`  
Variance power. Must satisfy `1 < p < 2` (compound Poisson-Gamma, the typical choice for insurance-type data with zeros) or `p > 2` (positive continuous, heavy-tailed data). Values of exactly `1` or `2` are rejected because they correspond to [Poisson](Poisson.md#whittaker.Poisson) and [Gamma](Gamma.md#whittaker.Gamma), which have simpler, exact deviance and likelihood formulas. If [p](Tweedie.md#whittaker.Tweedie.p) is unknown, use [tw()](tw.md#whittaker.tw) to estimate it from the data instead of fixing it here.


## Notes

The link is the natural logarithm:

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

The variance function is a power of the mean:

 V(\mu) = \mu^{p}, 

which interpolates between [Poisson](Poisson.md#whittaker.Poisson)-like behavior ([p](Tweedie.md#whittaker.Tweedie.p) near 1) and [Gamma](Gamma.md#whittaker.Gamma)-like behavior ([p](Tweedie.md#whittaker.Tweedie.p) near 2), or heavier-tailed behavior for `p > 2`. The unit deviance is

 d(y, \hat\mu) = 2 \left\[ \frac{y^{2-p}}{(1-p)(2-p)} - \frac{y\\\hat\mu^{1-p}}{1-p} + \frac{\hat\mu^{2-p}}{2-p} \right\], 

summed over observations to give the total deviance. Because the Tweedie density has no closed form for `1 < p < 2`, the log-likelihood is evaluated using the Dunn & Smyth (2005) saddlepoint approximation.


## Examples

Fit a GAM to compound Poisson-Gamma data with a point mass at zero:


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

rng = np.random.default_rng(0)
n = 300
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(1.0 + 0.5 * np.sin(x))

p = 1.5
scale = 1.0
lam = mu ** (2 - p) / (scale * (2 - p))
alpha = (2 - p) / (p - 1)
gamma_scale = scale * (p - 1) * mu ** (p - 1)
n_claims = rng.poisson(lam)
y = np.array(
    [rng.gamma(alpha, gamma_scale[i], size=n_claims[i]).sum() for i in range(n)]
)

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

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


    GAM fit summary
    ============================================================
    Formula:    y ~ s(x)
    Family:     Tweedie(p=1.5, link='log')
    Inference:  REML
    Observations: 300
    Coefficients: 10

    Parametric coefficients:
      Term                       Estimate    Std.Err    t value    p-value
      ------------------------ ---------- ---------- ---------- ----------
      (Intercept)                  1.0177     0.0492     20.686    < 1e-16

    Approximate significance of smooth terms:
      Term                        EDF Ref.df     Chi.sq    p-value
      ------------------------ ------ ------ ---------- ----------
      s(x)                       4.44      5     69.085  1.589e-13

    Total EDF:  5.44
    Scale est:  1.186925
    Deviance:   349.6177
    Null dev:   436.9088
    Dev. expl:  20.0%
    GCV score:  1.208856
    AIC:        1236.44
    BIC:        1256.59


## Attributes

| Name | Description |
|----|----|
| [p](#p) | Variance power p in V(\mu) = \mu^p. |
| [scale_known](#scale_known) | Whether the dispersion parameter is fixed. Always `False` for Tweedie. |

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


### p


Variance power p in V(\mu) = \mu^p.


`p: float`


Determines the shape of the Tweedie distribution: `1 < p < 2` gives the compound Poisson-Gamma case with a point mass at zero, `p = 3` corresponds to the Inverse Gaussian, and other `p > 2` give heavier-tailed positive continuous distributions.


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


### scale_known


Whether the dispersion parameter is fixed. Always `False` for Tweedie.


`scale_known: bool`


The dispersion `phi` is estimated from the data during fitting rather than fixed, independently of whether the variance power [p](Tweedie.md#whittaker.Tweedie.p) itself is fixed or estimated (see [TweedieEstimated](TweedieEstimated.md#whittaker.TweedieEstimated)).


## Methods

| Name | Description |
|----|----|
| [deviance()](#deviance) | Total Tweedie 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) | Tweedie log-likelihood via the Dunn & Smyth (2005) saddlepoint approximation. |
| [simulate()](#simulate) | Simulate Tweedie-distributed response values with mean `mu` and dispersion `scale`. |
| [unit_deviance()](#unit_deviance) | Per-observation Tweedie deviance contributions. |
| [variance()](#variance) | Tweedie variance function: V(\mu) = \mu^p for the current variance power [p](Tweedie.md#whittaker.Tweedie.p). |

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


### deviance()


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


Usage

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


#### Parameters


`y: NDArray`  
Observed response values (non-negative), 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`, values (including structural zeros typical of compound Poisson-Gamma data) are pushed away from zero to avoid `log(0)` on the first P-IRLS iteration.


#### Parameters


`y: NDArray`  
Observed response values (non-negative), 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()


Tweedie log-likelihood via the Dunn & Smyth (2005) saddlepoint approximation.


Usage

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


Since the Tweedie density has no closed form for `1 < p < 2`, the log-density at each observation is approximated by `_saddlepoint_log_density` and summed (optionally weighted) over observations.


#### Parameters


`y: NDArray`  
Observed response values (non-negative), shape `(n,)`.

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

`scale: float`  
Dispersion parameter \phi.

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


#### Returns


`float`  
The total (approximate) log-likelihood.


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


### simulate()


Simulate Tweedie-distributed response values with mean `mu` and dispersion `scale`.


Usage

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


For `1 < p < 2` (compound Poisson-Gamma), a Poisson number of claims is drawn for each observation and each claim's size is drawn from a Gamma distribution, matching the point mass at zero characteristic of this regime. For `p > 2`, values are drawn directly from a Gamma distribution matched to the Tweedie mean and variance.


#### Parameters


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

`scale: float`  
Dispersion parameter \phi.

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


#### Returns


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


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


### unit_deviance()


Per-observation Tweedie deviance contributions.


Usage

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


Computes d_i = 2 \left\[ \frac{y_i^{2-p}}{(1-p)(2-p)} - \frac{y_i\\\hat\mu_i^{1-p}}{1-p} + \frac{\hat\mu_i^{2-p}}{2-p} \right\] for the current variance power [p](Tweedie.md#whittaker.Tweedie.p), with the convention that the first term vanishes when y_i = 0.


#### Parameters


`y: NDArray`  
Observed response values (non-negative), shape `(n,)`.

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


#### Returns


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


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


### variance()


Tweedie variance function: V(\mu) = \mu^p for the current variance power [p](Tweedie.md#whittaker.Tweedie.p).


Usage

``` python
variance(mu)
```


#### Parameters


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


#### Returns


`NDArray`  
Variance values \mu^p, shape `(n,)`.