GammaLS

Gamma location-scale family for GAMLSS.

Usage

Source

GammaLS()

GammaLS extends the plain Gamma family by letting both the mean mu and the coefficient of variation sigma vary smoothly with covariates, instead of assuming a fixed shape parameter. Use it for strictly positive, right-skewed responses where not only the typical magnitude but also the relative spread (coefficient of variation) changes systematically across the range of the predictors — for example, cost or duration data whose relative volatility grows with the covariates rather than staying proportional to mu alone. Both parameters use the log link, keeping mu > 0 and sigma > 0.

Notes

GammaLS parameterizes the Gamma distribution by its mean mu > 0 and its coefficient of variation sigma > 0, where sigma = 1 / sqrt(shape) and shape = alpha = 1 / sigma^2. Both parameters use the log link:

g_{\mu}(\mu) = \log(\mu), \qquad g_{\sigma}(\sigma) = \log(\sigma).

The response density is the Gamma density with shape alpha = 1/sigma^2 and rate alpha/mu:

f(y \mid \mu, \sigma) = \frac{(\alpha/\mu)^{\alpha}}{\Gamma(\alpha)}\, y^{\alpha - 1} \exp\!\left(-\frac{\alpha y}{\mu}\right), \qquad \alpha = \frac{1}{\sigma^{2}}.

Because sigma is the coefficient of variation, \operatorname{Var}(Y) = \sigma^2 \mu^2, so this is a direct generalization of Gamma’s variance function V(mu) = mu^2 in which the proportionality constant sigma^2 is itself allowed to depend on covariates.

Examples

Fit a GAMLSS where both the mean and relative spread of a positive response vary smoothly:

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.4 * np.sin(x))
sigma = 0.2 + 0.15 * np.abs(np.cos(x))
shape = 1.0 / sigma**2
y = rng.gamma(shape, mu / shape)

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

model = wk.GAMLSS(
    formulas={"mu": "y ~ s(x)", "sigma": "y ~ s(x)"},
    family=wk.GammaLS(),
)
model.fit(data)
print(model.summary())
GAMLSS fit summary
========================================
Family: GammaLS(mu=log, sigma=log)
N obs: 300
Global deviance: 666.8366
AIC: 690.3882
BIC: 734.0031
Log-likelihood: -333.4183
Converged: True (5 iterations)

--- mu ---
  EDF total: 5.93
  Smooth 1: edf = 4.93

--- sigma ---
  EDF total: 5.84
  Smooth 1: edf = 4.84

Attributes

Name Description
parameter_names Names of the distributional parameters modeled by GammaLS.

parameter_names

Names of the distributional parameters modeled by GammaLS.

parameter_names: tuple[str, …]

Methods

Name Description
d2l_dtheta2() Expected Fisher information for mu or sigma.
dl_dtheta() First derivative of the Gamma log-likelihood with respect to mu or sigma.
initialize() Starting values for mu and sigma from the raw response.
link() Apply the log link for mu or sigma.
link_derivative() Derivative of the log link for mu or sigma.
link_inverse() Apply the inverse log link for mu or sigma.
log_likelihood() Total Gamma log-likelihood at the current mu and sigma.
simulate() Simulate response values from a Gamma distribution with the given mu, sigma.

d2l_dtheta2()

Expected Fisher information for mu or sigma.

Usage

Source

d2l_dtheta2(
    param,
    y,
    params,
)

With shape alpha = 1/sigma^2 and trigamma function \psi_1:

-E\!\left[\frac{\partial^2 \ell}{\partial \mu^2}\right] = \frac{\alpha}{\mu^2}, \qquad -E\!\left[\frac{\partial^2 \ell}{\partial \sigma^2}\right] = \frac{4}{\sigma^4}\big(\alpha\,\psi_1(\alpha) - 1\big).

The sigma term is additionally floored at _SIGMA_FLOOR to keep the working weight strictly positive.

Parameters

param: str

Either "mu" or "sigma".

y: NDArray

Observed response values, shape (n,) (unused, kept for interface consistency).

params: dict[str, NDArray]
Current estimates with keys "mu" and "sigma", each of shape (n,).

Returns

NDArray
Elementwise working weights, shape (n,), guaranteed positive.

dl_dtheta()

First derivative of the Gamma log-likelihood with respect to mu or sigma.

Usage

Source

dl_dtheta(
    param,
    y,
    params,
)

With shape alpha = 1/sigma^2, the score functions are

\frac{\partial \ell}{\partial \mu} = \frac{\alpha}{\mu}\left(\frac{y}{\mu} - 1\right), \qquad \frac{\partial \ell}{\partial \sigma} = \frac{2}{\sigma^3}\left( \psi(\alpha) - \log\alpha - \log\frac{y}{\mu} - 1 + \frac{y}{\mu}\right),

where \psi is the digamma function.

Parameters

param: str

Either "mu" or "sigma".

y: NDArray

Observed response values, shape (n,).

params: dict[str, NDArray]
Current estimates with keys "mu" and "sigma", each of shape (n,).

Returns

NDArray
Elementwise first derivatives, shape (n,).

initialize()

Starting values for mu and sigma from the raw response.

Usage

Source

initialize(y)

mu is initialized at the (floored) observed values themselves, and sigma at the sample coefficient of variation of y (clamped to at least 0.1), constant across observations.

Parameters

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

Returns

dict[str, NDArray]
{"mu": y_safe.copy(), "sigma": <constant array of the sample CV>}.

log_likelihood()

Total Gamma log-likelihood at the current mu and sigma.

Usage

Source

log_likelihood(
    y,
    params,
)

Evaluated using the shape/rate parameterization alpha = 1/sigma^2, rate = alpha/mu, with y and mu floored at _MU_FLOOR and sigma floored at _SIGMA_FLOOR for numerical stability.

Parameters

y: NDArray

Observed response values, shape (n,).

params: dict[str, NDArray]
Current estimates with keys "mu" and "sigma", each of shape (n,).

Returns

float
The log-likelihood summed over all observations.

simulate()

Simulate response values from a Gamma distribution with the given mu, sigma.

Usage

Source

simulate(
    params,
    rng,
)

Converts to the shape/scale parameterization (shape = 1/sigma^2, scale = mu * sigma^2) expected by rng.gamma.

Parameters

params: dict[str, NDArray]

Current estimates with keys "mu" and "sigma", each of shape (n,).

rng: np.random.Generator
A NumPy random generator used to draw the simulated values.

Returns

NDArray
Simulated response values, shape (n,).