Tweedie family with log link.
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). The Poisson (p = 1) and Gamma (p = 2) boundary cases are excluded here — use Poisson or 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 and Gamma, which have simpler, exact deviance and likelihood formulas. If p is unknown, use 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-like behavior (p near 1) and Gamma-like behavior (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:
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
|
Variance power p in V(\mu) = \mu^p.
|
|
scale_known
|
Whether the dispersion parameter is fixed. Always False for Tweedie.
|
p
Variance power p in V(\mu) = \mu^p.
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.
The dispersion phi is estimated from the data during fitting rather than fixed, independently of whether the variance power p itself is fixed or estimated (see TweedieEstimated).
Methods
|
Name
|
Description
|
|
deviance()
|
Total Tweedie deviance, the (weighted) sum of unit_deviance.
|
|
initialize()
|
Starting values for mu: y nudged away from zero.
|
|
link()
|
Apply the log link: \eta = \log(\mu).
|
|
link_derivative()
|
Derivative of the log link: g'(\mu) = 1/\mu.
|
|
link_inverse()
|
Apply the inverse log link: \mu = e^{\eta}.
|
|
log_likelihood()
|
Tweedie log-likelihood via the Dunn & Smyth (2005) saddlepoint approximation.
|
|
simulate()
|
Simulate Tweedie-distributed response values with mean mu and dispersion scale.
|
|
unit_deviance()
|
Per-observation Tweedie deviance contributions.
|
|
variance()
|
Tweedie variance function: V(\mu) = \mu^p for the current variance power p.
|
deviance()
Total Tweedie deviance, the (weighted) sum of unit_deviance.
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.
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).
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.
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}.
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.
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.
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.
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, 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.
Parameters
mu: NDArray-
Conditional mean values, shape
(n,).
Returns
NDArray-
Variance values \mu^p, shape
(n,).