GAMLSSFamily

Abstract base class for GAMLSS distributional families.

Usage

GAMLSSFamily()

Ordinary Family subclasses (used by GAM) model only the mean mu of the response as a function of covariates, treating any other distributional parameters (e.g. the variance or dispersion) as constant across observations. A GAMLSSFamily, used by GAMLSS instead of GAM, generalizes the single-parameter Family abstraction to the location-scale-shape setting: instead of one response parameter with a mean link, it defines a full response distribution with K named parameters \theta_1, \ldots, \theta_K (e.g. location mu and scale sigma), each with its own link function g_k and its own additive predictor \eta_k = g_k(\theta_k). This is useful whenever more than the mean of the response changes systematically with covariates — for example, when the spread (heteroscedasticity), skew, or zero-inflation probability also varies across the range of the predictors.

GAMLSS fits every parameter’s additive predictor jointly by alternating penalized IRLS updates across parameters (the RS algorithm), which relies on each family supplying the per-parameter score dl_dtheta, (expected) Fisher information d2l_dtheta2, link/inverse-link/link-derivative, log-likelihood, initial values, and a simulate method. Subclasses must implement all of the abstract methods below.

Whittaker ships with the following concrete GAMLSS families:

GaussianLS — location-scale Gaussian: identity link for the mean, log link for the standard deviation.
GammaLS — location-scale Gamma: log link for both the mean and the coefficient of variation.
BetaLS — mean-precision Beta: logit link for the mean, log link for the precision.
ZeroInflatedPoisson — Poisson mean plus a zero-inflation probability, for count data with excess zeros.
ZeroInflatedNegativeBinomial — overdispersed counts with excess zeros, combining NegativeBinomial-style overdispersion with zero-inflation.

Examples

GAMLSS families are passed to GAMLSS, and each distributional parameter gets its own formula:

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.sin(x)
sigma = 0.2 + 0.3 * np.abs(np.cos(x))
y = rng.normal(mu, sigma)

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

# Model both the mean and the standard deviation as smooth functions of x
model = wk.GAMLSS(
    formulas={"mu": "y ~ s(x)", "sigma": "y ~ s(x)"},
    family=wk.GaussianLS(),
)
model.fit(data)
print(model.summary())

GAMLSS fit summary
========================================
Family: GaussianLS(mu=identity, sigma=log)
N obs: 300
Global deviance: 271.2831
AIC: 297.3004
BIC: 345.4816
Log-likelihood: -135.6416
Converged: True (4 iterations)

--- mu ---
  EDF total: 6.63
  Smooth 1: edf = 5.63

--- sigma ---
  EDF total: 6.38
  Smooth 1: edf = 5.38

Attributes

Name	Description
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, …]

Every concrete GAMLSSFamily names its parameters explicitly, e.g. ('mu', 'sigma') for GaussianLS and GammaLS, or ('mu', 'phi') for BetaLS. GAMLSS uses this tuple to determine which formulas it expects (one per name) and the order in which it cycles through parameters during the RS fitting algorithm.

Methods

Name	Description
d2l_dtheta2()	Negative expected second derivative of the log-likelihood with respect to `param`.
dl_dtheta()	First derivative of the log-likelihood with respect to `param`.
initialize()	Starting values for all distributional parameters given `y`.
link()	Apply the link function for `param`: `eta = g(theta)`.
link_derivative()	Derivative of the link for `param`: `d(eta)/d(theta)`.
link_inverse()	Apply the inverse link for `param`: `theta = g^-1(eta)`.
log_likelihood()	Full log-likelihood evaluated at the given parameter values.
simulate()	Simulate response values from the distribution.

d2l_dtheta2()

Negative expected second derivative of the log-likelihood with respect to param.

Usage

Source

d2l_dtheta2(
    param,
    y,
    params,
)

Returns the (expected) Fisher information -E[\partial^2 \ell / \partial \theta^2] for param, evaluated elementwise at the current parameter estimates. Must return positive values, since it is used directly as a working weight during fitting: it both scales the working pseudo-response and forms the diagonal weight matrix for the penalized IRLS update of param’s additive predictor.

Parameters

param: str: Name of the distributional parameter to differentiate with respect to, one of parameter_names.
y: NDArray: Observed response values, shape (n,).
params: dict[str, NDArray]: Current estimates of all distributional parameters, keyed by name, each of shape (n,).

Returns

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

dl_dtheta()

First derivative of the log-likelihood with respect to param.

Usage

Source

dl_dtheta(
    param,
    y,
    params,
)

Returns \partial \ell / \partial \theta evaluated elementwise at the current parameter estimates. This score forms the basis of the working pseudo-response used by GAMLSS’s RS algorithm when updating the additive predictor for param, holding the other distributional parameters fixed.

Parameters

param: str: Name of the distributional parameter to differentiate with respect to, one of parameter_names.
y: NDArray: Observed response values, shape (n,).
params: dict[str, NDArray]: Current estimates of all distributional parameters, keyed by name, each of shape (n,).

Returns

NDArray: Elementwise first derivatives, shape (n,).

initialize()

Starting values for all distributional parameters given y.

Usage

Source

initialize(y)

Produces a reasonable initial guess for every parameter in parameter_names from the raw response alone, before any covariate information is used. GAMLSS uses these as the starting point for the first RS iteration.

Parameters

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

Returns

dict[str, NDArray]: Initial values for each parameter, keyed by name, each of shape (n,).

link()

Apply the link function for param: eta = g(theta).

Usage

Source

link(
    param,
    values,
)

Maps a distributional parameter from its natural scale (e.g. mu on the real line, sigma > 0) onto the unconstrained scale of the additive predictor eta, so that GAMLSS can model it with an ordinary smooth linear predictor.

Parameters

param: str: Name of the distributional parameter being linked, one of parameter_names.
values: NDArray: Parameter values on the natural scale, shape (n,).

Returns

NDArray: Linked values eta = g(theta), shape (n,).

link_derivative()

Derivative of the link for param: d(eta)/d(theta).

Usage

Source

link_derivative(
    param,
    values,
)

Used by the RS algorithm to convert the score dl_dtheta into a working pseudo-response on the scale of the additive predictor, via the chain rule \partial \ell / \partial \eta = (\partial \ell / \partial \theta)\, (d\eta/d\theta)^{-1}.

Parameters

param: str: Name of the distributional parameter being linked, one of parameter_names.
values: NDArray: Parameter values on the natural scale at which to evaluate the derivative, shape (n,).

Returns

NDArray: Elementwise derivatives d(eta)/d(theta), shape (n,).

link_inverse()

Apply the inverse link for param: theta = g^-1(eta).

Usage

Source

link_inverse(
    param,
    eta,
)

Maps the additive predictor back onto the natural scale of the parameter, enforcing any constraints (e.g. positivity for a scale parameter, or the unit interval for a probability) implied by the chosen link.

Parameters

param: str: Name of the distributional parameter being linked, one of parameter_names.
eta: NDArray: Additive predictor values, shape (n,).

Returns

NDArray: Parameter values on the natural scale, shape (n,).

log_likelihood()

Full log-likelihood evaluated at the given parameter values.

Usage

Source

log_likelihood(
    y,
    params,
)

Sums the per-observation log-density of the response distribution across all observations, using the current estimate of every distributional parameter. Used for convergence checks during fitting and for computing information criteria such as AIC/BIC in model summaries.

Parameters

y: NDArray: Observed response values, shape (n,).
params: dict[str, NDArray]: Current estimates of all distributional parameters, keyed by name, each of shape (n,).

Returns

float: The total log-likelihood summed over all observations.

simulate()

Simulate response values from the distribution.

Usage

Source

simulate(
    params,
    rng,
)

Draws one random sample per observation from the response distribution at the given per-observation parameter values. Used for posterior/parametric-bootstrap style simulation from a fitted GAMLSS model.

Parameters

params: dict[str, NDArray]: Distributional parameter values to simulate from, keyed by name, each of shape (n,).
rng: np.random.Generator: A NumPy random generator (e.g. numpy.random.default_rng()) used to draw the simulated values.

Returns

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