Generalized Additive Model for Location, Scale, and Shape.
GAMLSS(
formulas,
family=None,
)
A GAMLSS extends the ordinary GAM by allowing every parameter of the response distribution, not just its mean, to depend on covariates through its own smooth additive predictor. For a distribution with parameters theta_1, ..., theta_K (e.g. location mu, scale sigma, and possibly shape parameters nu, tau), each parameter has its own link function g_k and its own formula:
g_k(\theta_k) = \eta_k = X_k \beta_k, \quad k = 1, \dots, K
This makes it possible to model, for example, both the mean and the variance of y as smooth functions of x, which ordinary (mean-only) GAMs cannot do. Fitting uses the RS (“Rigby and Stasinopoulos”) algorithm, which cycles through the parameters, holding all but one fixed, updating it via penalized IRLS, and repeating until the penalized log-likelihood converges.
Use GAMLSS when the assumption of a fixed dispersion (constant variance, constant shape) is implausible, such as heteroscedastic regression, or regression with distributions like the negative binomial or beta that have separate location and shape parameters.
Parameters
formulas: dict[str, str]-
Dict mapping parameter names to formula strings. All formulas must share the same response variable. Example: {"mu": "y ~ s(x1)", "sigma": "y ~ s(x2)"}. The set of keys must match family.parameter_names exactly.
family: GAMLSSFamily | None = None-
A GAMLSSFamily specifying the distributional model, including the number and names of parameters, their link functions, and log-likelihood derivatives used by the RS algorithm. Defaults to GaussianLS() (Gaussian location-scale, i.e.
mu and sigma both modeled).
Notes
Rigby & Stasinopoulos (2005) formulate GAMLSS fitting as penalized maximum likelihood. Within each outer RS iteration, and for each parameter theta_k in turn, a working response and weight are formed from the score and Fisher information of the log-likelihood with respect to theta_k:
z_k = \eta_k + \frac{\partial \ell / \partial \theta_k}{\partial^2 \ell / \partial \theta_k^2}
\cdot g_k'(\theta_k), \qquad w_k = -\frac{\partial^2 \ell}{\partial \theta_k^2}
\Big/ g_k'(\theta_k)^2
and a penalized weighted least squares problem is solved for beta_k, with all other parameters held at their current fitted values. Smoothing parameters for each parameter’s smooth terms can be selected by GCV, REML, or ML at every inner iteration. The algorithm alternates over parameters until the global deviance -2 * log_likelihood stops improving.
Examples
import numpy as np
from whittaker.gamlss import GAMLSS
from whittaker.families.gaussian_ls import GaussianLS
rng = np.random.default_rng(0)
n = 500
x = rng.uniform(0, 1, n)
mu = np.sin(2 * np.pi * x)
sigma = np.exp(-1 + 2 * x)
y = rng.normal(mu, sigma)
model = GAMLSS(
formulas={"mu": "y ~ s(x)", "sigma": "y ~ s(x)"},
family=GaussianLS(),
)
model.fit({"x": x, "y": y}, method="REML")
pred = model.predict({"x": x[:5]})
print(pred.values)
{'mu': array([-0.92639392, 0.91684274, 0.3068777 , 0.24081722, -1.08274642]), 'sigma': array([1.22798721, 0.57377032, 0.35711457, 0.33947392, 1.76958276])}
Attributes
|
Name
|
Description
|
|
aic
|
Akaike information criterion of the fitted model.
|
|
bic
|
Bayesian information criterion of the fitted model.
|
|
converged
|
Whether the RS algorithm converged before hitting max_outer iterations.
|
|
family
|
The GAMLSSFamily used to build this model.
|
|
global_deviance
|
Global deviance of the fitted model.
|
|
is_fitted
|
Whether fit() has been called successfully.
|
|
log_likelihood
|
Log-likelihood of the fitted model at convergence.
|
|
n_iter
|
Number of outer RS iterations actually performed.
|
aic
Akaike information criterion of the fitted model.
Computed from the global deviance and the total effective degrees of freedom summed across all distributional parameters’ smooth and parametric terms, and can be used to compare GAMLSS models fit with different formulas or families on the same data.
bic
Bayesian information criterion of the fitted model.
Like aic, but penalizes model complexity more heavily (using log(n) instead of 2 as the multiplier on total effective degrees of freedom), so it tends to favor simpler models when comparing fits.
converged
Whether the RS algorithm converged before hitting max_outer iterations.
global_deviance
Global deviance of the fitted model.
The RS algorithm minimizes this quantity at each outer iteration; it is defined as -2 * log_likelihood and is the primary measure of fit used to check convergence.
is_fitted
Whether fit() has been called successfully.
log_likelihood
Log-likelihood of the fitted model at convergence.
n_iter
Number of outer RS iterations actually performed.
Methods
|
Name
|
Description
|
|
coefficients()
|
Fitted basis coefficients for one distributional parameter.
|
|
edf()
|
Effective degrees of freedom per smooth term for one distributional parameter.
|
|
fit()
|
Fit the GAMLSS model via the RS algorithm.
|
|
fitted_values()
|
Fitted values on the response scale for one or all distributional parameters.
|
|
predict()
|
Predict distributional parameters for new data.
|
|
simulate()
|
Simulate responses from the fitted model.
|
|
smoothing_params()
|
Fitted smoothing parameters for one distributional parameter’s smooth terms.
|
|
summary()
|
Return a text summary of the fitted model.
|
coefficients()
Fitted basis coefficients for one distributional parameter.
Parameters
parameter: str-
Name of the distributional parameter (must be one of
family.parameter_names), e.g. "mu" or "sigma".
Returns
NDArray-
Coefficient vector
beta_k for that parameter’s linear predictor eta_k = X_k @ beta_k, in the order of its model matrix’s columns.
edf()
Effective degrees of freedom per smooth term for one distributional parameter.
Parameters
parameter: str-
Name of the distributional parameter, e.g.
"mu" or "sigma".
Returns
list[float]-
EDF of each smooth term in that parameter’s formula, reflecting how much smoothing was applied (lower values indicate heavier penalization toward linearity).
fit()
Fit the GAMLSS model via the RS algorithm.
fit(
data,
*,
method="GCV",
max_outer=50,
max_inner=20,
tol=1e-06,
select=False,
)
Parses each parameter’s formula, builds its model matrix (basis functions, penalties, and any parametric terms), and then runs the outer Rigby & Stasinopoulos loop: for each parameter in turn, forms the working response and IRLS weights from the family’s likelihood derivatives, selects smoothing parameters (if the formula includes smooth terms), and solves a penalized weighted least squares problem, holding the other parameters fixed. This repeats until the global deviance stabilizes or max_outer iterations are reached.
Parameters
data: InputData-
Column-oriented data dict (or any type accepted by prepare_data) containing the shared response column and all covariates referenced in formulas.
method: str = "GCV"-
Smoothing parameter selection method applied to every parameter’s smooth terms: "GCV" (generalized cross-validation, the default), "REML" (restricted maximum likelihood), or "ML" (maximum likelihood).
max_outer: int = 50-
Maximum number of outer RS iterations (full passes over all parameters).
max_inner: int = 20-
Maximum inner iterations per parameter within a single outer step, used to converge the penalized IRLS fit for that parameter before moving to the next one.
tol: float = 1e-06-
Relative convergence tolerance, applied both to the coefficient update within a parameter’s inner loop and to the change in overall log-likelihood between outer iterations.
select: bool = False-
If
True, add an extra shrinkage penalty to each smooth’s null space so that terms can be shrunk essentially to zero, enabling automatic term selection.
Returns
GAMLSS-
Returns
self, with is_fitted now True and per-parameter results accessible via coefficients(), fitted_values(), edf(), and related accessors.
fitted_values()
Fitted values on the response scale for one or all distributional parameters.
fitted_values(parameter=None)
Parameters
parameter: str | None = None-
If given, return only this parameter’s fitted values as a single array. If
None (default), return a dict of fitted values for every distributional parameter.
Returns
dict[str, NDArray] or NDArray-
theta_k = g_k^{-1}(eta_k) for the requested parameter(s), evaluated at the training covariate values used in fit().
predict()
Predict distributional parameters for new data.
predict(
new_data,
*,
parameter=None,
se=False,
)
For each parameter theta_k, builds the prediction design matrix from the fitted basis (using the same knots/constraints as training), forms the linear predictor eta_k = X_new @ beta_k (plus any offset), and maps it back to the response scale via the parameter’s inverse link, theta_k = g_k^{-1}(eta_k). When se=True, the standard error of eta_k is also computed from the Bayesian posterior covariance of beta_k.
Parameters
new_data: InputData-
Column-oriented new data dict containing all covariates used in the fitted formulas.
parameter: str | None = None-
If given, return only this parameter’s predicted values as an array on the response scale. Otherwise return a GAMLSSPrediction with all parameters.
se: bool = False-
If
True, compute standard errors on the linear predictor scale for each parameter, using the Bayesian covariance V_beta_k = (X_k' W_k X_k + S_k)^{-1} implied by the final IRLS weights and smoothing parameters for that parameter.
Returns
GAMLSSPrediction or NDArray-
A GAMLSSPrediction with all parameters’ values, linear predictors, and (optionally) standard errors, or a single
NDArray of predicted values if parameter is given.
simulate()
Simulate responses from the fitted model.
simulate(
n_sim=1,
*,
seed=None,
)
Draws n_sim independent replicate response vectors from the fitted distribution, using each observation’s fitted parameter values (mu, sigma, and any shape parameters) as the distribution’s parameters. This is useful for posterior-predictive checks, e.g. comparing the distribution of simulated data against the observed response.
Parameters
n_sim: int = 1-
Number of independent simulated replicates to draw per observation.
seed: int | None = None-
Random seed for reproducibility.
Returns
NDArray-
Simulated values, shape
(n, n_sim).
smoothing_params()
Fitted smoothing parameters for one distributional parameter’s smooth terms.
smoothing_params(parameter)
Parameters
parameter: str-
Name of the distributional parameter, e.g.
"mu" or "sigma".
Returns
list[float]-
One smoothing parameter (
lambda) per smooth term in that parameter’s formula, in the order the terms were declared. Empty if the formula has no smooth terms.
summary()
Return a text summary of the fitted model.
Includes global fit statistics (global deviance, AIC, BIC, log-likelihood, convergence status and iteration count) followed by a per-parameter section listing the total effective degrees of freedom and, if present, the EDF of each individual smooth term.
Returns
str-
Multi-line, human-readable summary suitable for printing.