# GAMLSS Generalized Additive Model for Location, Scale, and Shape. Usage ``` python 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](QuantileFamily.md#whittaker.QuantileFamily.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](GAMLSS.md#whittaker.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](GAMLSSFamily.md#whittaker.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()](GaussianLS.md#whittaker.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 ``` python 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](#aic) | Akaike information criterion of the fitted model. | | [bic](#bic) | Bayesian information criterion of the fitted model. | | [converged](#converged) | Whether the RS algorithm converged before hitting `max_outer` iterations. | | [family](#family) | The [GAMLSSFamily](GAMLSSFamily.md#whittaker.GAMLSSFamily) used to build this model. | | [global_deviance](#global_deviance) | Global deviance of the fitted model. | | [is_fitted](#is_fitted) | Whether `fit()` has been called successfully. | | [log_likelihood](#log_likelihood) | Log-likelihood of the fitted model at convergence. | | [n_iter](#n_iter) | Number of outer RS iterations actually performed. | ------------------------------------------------------------------------ ### aic Akaike information criterion of the fitted model. `aic: float` 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. `bic: float` 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. `converged: bool` ------------------------------------------------------------------------ ### family The [GAMLSSFamily](GAMLSSFamily.md#whittaker.GAMLSSFamily) used to build this model. `family: GAMLSSFamily` ------------------------------------------------------------------------ ### global_deviance Global deviance of the fitted model. `global_deviance: float` 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. `is_fitted: bool` ------------------------------------------------------------------------ ### log_likelihood Log-likelihood of the fitted model at convergence. `log_likelihood: float` ------------------------------------------------------------------------ ### n_iter Number of outer RS iterations actually performed. `n_iter: int` ## Methods | Name | Description | |----|----| | [coefficients()](#coefficients) | Fitted basis coefficients for one distributional parameter. | | [edf()](#edf) | Effective degrees of freedom per smooth term for one distributional parameter. | | [fit()](#fit) | Fit the GAMLSS model via the RS algorithm. | | [fitted_values()](#fitted_values) | Fitted values on the response scale for one or all distributional parameters. | | [predict()](#predict) | Predict distributional parameters for new data. | | [simulate()](#simulate) | Simulate responses from the fitted model. | | [smoothing_params()](#smoothing_params) | Fitted smoothing parameters for one distributional parameter's smooth terms. | | [summary()](#summary) | Return a text summary of the fitted model. | ------------------------------------------------------------------------ ### coefficients() Fitted basis coefficients for one distributional parameter. Usage ``` python coefficients(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. Usage ``` python edf(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. Usage ``` python 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. Usage ``` python 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. Usage ``` python 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](GAMLSSPrediction.md#whittaker.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](GAMLSSPrediction.md#whittaker.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. Usage ``` python 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. Usage ``` python 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. Usage ``` python summary() ``` 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.