Gaussian location-scale family for GAMLSS.
GaussianLS extends the plain Gaussian family to allow both the mean mu and the standard deviation sigma to vary smoothly with covariates, rather than assuming constant variance. Use it when a continuous, approximately symmetric response shows heteroscedasticity — for example, when the spread of measurements grows or shrinks over the range of a predictor — and you want the model to capture that varying spread rather than average it away. Each parameter has its own additive predictor and its own link function: mu uses the identity link (as in Gaussian), and sigma uses the log link, which keeps the fitted standard deviation positive.
Notes
The two parameters use distinct link functions:
g_{\mu}(\mu) = \mu \qquad \text{(identity)}, \qquad
g_{\sigma}(\sigma) = \log(\sigma).
The response density is the ordinary Gaussian density evaluated at the fitted, observation- specific mu and sigma:
f(y \mid \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}}
\exp\!\left(-\frac{(y-\mu)^2}{2\sigma^2}\right),
so the log-likelihood contribution for a single observation is \ell_i = -\log\sigma_i - \tfrac{1}{2}\log(2\pi)
- \tfrac{1}{2}\left(\frac{y_i - \mu_i}{\sigma_i}\right)^2. Unlike Gaussian, there is no separate scale parameter to estimate: sigma itself is the quantity being modeled by its own smooth predictor.
Examples
Fit a GAMLSS with a smoothly varying mean and standard deviation:
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 = 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
parameter_names
Names of the distributional parameters modeled by GaussianLS.
parameter_names: tuple[str, …]
Methods
|
Name
|
Description
|
|
d2l_dtheta2()
|
Expected Fisher information for mu or sigma.
|
|
dl_dtheta()
|
First derivative of the Gaussian log-likelihood with respect to mu or sigma.
|
|
initialize()
|
Starting values for mu and sigma from the raw response.
|
|
link()
|
Apply the link function for mu or sigma.
|
|
link_derivative()
|
Derivative of the link function for mu or sigma.
|
|
link_inverse()
|
Apply the inverse link for mu or sigma.
|
|
log_likelihood()
|
Total Gaussian log-likelihood at the current mu and sigma.
|
|
simulate()
|
Simulate response values from Normal(mu, sigma).
|
d2l_dtheta2()
Expected Fisher information for mu or sigma.
d2l_dtheta2(
param,
y,
params,
)
For the Gaussian distribution these expected information terms do not depend on y:
-E\!\left[\frac{\partial^2 \ell}{\partial \mu^2}\right] = \frac{1}{\sigma^2}, \qquad
-E\!\left[\frac{\partial^2 \ell}{\partial \sigma^2}\right] = \frac{2}{\sigma^2}.
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,).
dl_dtheta()
First derivative of the Gaussian log-likelihood with respect to mu or sigma.
dl_dtheta(
param,
y,
params,
)
The score functions are
\frac{\partial \ell}{\partial \mu} = \frac{y - \mu}{\sigma^2}, \qquad
\frac{\partial \ell}{\partial \sigma} = -\frac{1}{\sigma} + \frac{(y-\mu)^2}{\sigma^3}.
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.
mu is initialized at the observed values themselves and sigma at the sample standard deviation of y, constant across observations.
Parameters
y: NDArray-
Observed response values, shape
(n,).
Returns
dict[str, NDArray]-
{"mu": y.copy(), "sigma": <constant array of the sample std>}.
link()
Apply the link function for mu or sigma.
mu uses the identity link and sigma uses the log link, so that sigma’s additive predictor is unconstrained while the fitted standard deviation stays positive after applying link_inverse.
Parameters
param: str-
Either "mu" or "sigma".
values: NDArray-
Parameter values on the natural scale, shape
(n,).
Returns
NDArray-
Linked values: values unchanged for
"mu", log(values) for "sigma".
link_derivative()
Derivative of the link function for mu or sigma.
link_derivative(
param,
values,
)
For the identity link, d\eta/d\mu = 1; for the log link, d\eta/d\sigma = 1/\sigma.
Parameters
param: str-
Either "mu" or "sigma".
values: NDArray-
Parameter values on the natural scale at which to evaluate the derivative, shape
(n,).
Returns
NDArray-
An array of ones for
"mu", or 1 / values for "sigma".
link_inverse()
Apply the inverse link for mu or sigma.
link_inverse(
param,
eta,
)
Maps the additive predictor eta back to the natural scale: unchanged for mu (identity link), exponentiated for sigma (log link) so it stays positive.
Parameters
param: str-
Either "mu" or "sigma".
eta: NDArray-
Additive predictor values, shape
(n,).
Returns
NDArray-
eta unchanged for "mu", exp(eta) for "sigma".
log_likelihood()
Total Gaussian log-likelihood at the current mu and sigma.
log_likelihood(
y,
params,
)
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 Normal(mu, sigma).
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,).