# SPDX-License-Identifier: MIT
# Copyright (c) 2025 Lucy Eduardo Arias
"""Abstract base classes for GLM distribution families and link functions.
Mathematical Framework
----------------------
This module defines the core abstractions for exponential family distributions
and link functions used in Generalized Linear Models (GLM).
Exponential Family of Distributions
------------------------------------
A distribution belongs to the exponential family if its density/mass function
can be written as:
f(y; θ, φ) = exp{[y θ - b(θ)] / a(φ) + c(y, φ)}
where:
- y is the response variable
- θ is the canonical (natural) parameter
- φ is the dispersion parameter
- b(θ) is the cumulant function
- a(φ) is a known function (typically a(φ) = φ/w for weight w)
- c(y, φ) is a normalization term
**Properties**:
1. **Mean-variance relationship**:
- E[Y] = μ = b'(θ)
- Var(Y) = a(φ) b''(θ) = a(φ) V(μ)
where V(μ) is the variance function (family-specific).
2. **Log-likelihood**:
ℓ(θ, φ; y) = [y θ - b(θ)] / a(φ) + c(y, φ)
3. **Score function**:
∂ℓ/∂θ = [y - μ] / a(φ)
Link Functions
--------------
A link function g(·) relates the mean μ to the linear predictor η:
g(μ) = η = X^T β
**Requirements**:
- Monotonic: g is strictly increasing or decreasing
- Differentiable: g'(μ) exists and is continuous
- Domain: g maps (support of μ) → ℝ
**Canonical link**:
The canonical link is g(μ) = θ, i.e., when the linear predictor equals
the natural parameter. This leads to simplifications:
- Fisher information matrix is particularly simple
- IRLS = Fisher scoring
**Common links**:
| Family | Canonical Link | Other Links |
|-----------|----------------|-----------------------|
| Gaussian | Identity | Log, Inverse |
| Poisson | Log | Identity, Sqrt |
| Binomial | Logit | Probit, CLogLog |
| Gamma | Inverse | Log, Identity |
Variance Functions
------------------
The variance function V(μ) characterizes the mean-variance relationship:
Var(Y) = φ V(μ)
where φ is the dispersion parameter.
**Standard families**:
| Family | Support | V(μ) | φ |
|-----------|--------------|--------------|------------------|
| Gaussian | (-∞, ∞) | 1 | σ² |
| Poisson | {0,1,2,...} | μ | 1 (fixed) |
| Binomial | {0,1,...,n} | μ(n-μ)/n | 1 (fixed) |
| Gamma | (0, ∞) | μ² | shape^{-1} |
Deviance
--------
The deviance measures goodness-of-fit:
D(y; μ) = 2 φ [ℓ(y; y) - ℓ(μ; y)]
where ℓ(y; y) is the saturated model log-likelihood.
**Properties**:
- D ≥ 0 with equality iff μ = y
- Asymptotically χ² distributed under H₀: correct model
- Used for model comparison via likelihood ratio tests
Initialization
--------------
Starting values μ⁽⁰⁾ for IRLS must satisfy:
1. μ⁽⁰⁾ ∈ support of the distribution
2. V(μ⁽⁰⁾) > 0 (non-degenerate)
3. Avoid boundary of parameter space
**Common strategies**:
- Gaussian: μ⁽⁰⁾ = y
- Poisson: μ⁽⁰⁾ = y + 0.1 (avoid zeros)
- Binomial: μ⁽⁰⁾ = (y + 0.5) / (n + 1) (avoid 0 and 1)
- Gamma: μ⁽⁰⁾ = y (already positive)
Implementation Notes
--------------------
**Multi-backend support**:
All implementations must work with NumPy, PyTorch, and JAX via the
array namespace abstraction. Use the `namespace()` function to get
the appropriate backend.
**Numerical stability**:
- Avoid log(0), exp(large), division by zero
- Use log-sum-exp trick where appropriate
- Clip μ to valid range with small epsilon margin
**Type annotations**:
- Array: Generic array type (works with all backends)
- Scalar: Single numerical value (float or int)
Extending with Custom Families
-------------------------------
To add a new distribution family:
1. Inherit from `Family`:
```python
from aurora.distributions.base import Family, LinkFunction
class MyFamily(Family):
def log_likelihood(self, y, mu, **params):
# Implement log p(y | μ)
...
def deviance(self, y, mu, **params):
# Implement 2[ℓ(y; y) - ℓ(μ; y)]
...
def variance(self, mu, **params):
# Implement V(μ)
...
def initialize(self, y):
# Return starting μ⁽⁰⁾
...
@property
def default_link(self):
# Return canonical LinkFunction
return MyLink()
```
2. Implement the corresponding link function:
```python
class MyLink(LinkFunction):
def link(self, mu):
return ... # g(μ)
def inverse(self, eta):
return ... # g⁻¹(η)
def derivative(self, mu):
return ... # dg/dμ
```
3. Register in the family registry (optional):
```python
_FAMILY_REGISTRY['myfamily'] = MyFamily
```
References
----------
**Exponential family theory**:
- Barndorff-Nielsen, O. (2014). *Information and Exponential Families in
Statistical Theory*. John Wiley & Sons.
- Brown, L. D. (1986). *Fundamentals of Statistical Exponential Families*.
Institute of Mathematical Statistics.
**GLM families**:
- McCullagh, P., & Nelder, J. A. (1989). *Generalized Linear Models* (2nd ed.).
Chapman and Hall/CRC. Chapter 2: Model Specification.
https://doi.org/10.1007/978-1-4899-3242-6
- Dobson, A. J., & Barnett, A. G. (2018). *An Introduction to Generalized
Linear Models* (4th ed.). CRC Press.
**Variance functions and quasi-likelihood**:
- Wedderburn, R. W. M. (1974). "Quasi-likelihood functions, generalized linear
models, and the Gauss-Newton method." *Biometrika*, 61(3), 439-447.
https://doi.org/10.1093/biomet/61.3.439
**Deviance and model selection**:
- Hastie, T., Tibshirani, R., & Friedman, J. (2009). *The Elements of
Statistical Learning: Data Mining, Inference, and Prediction* (2nd ed.).
Springer. https://doi.org/10.1007/978-0-387-84858-7 (Chapter 7)
See Also
--------
aurora.distributions.families : Concrete family implementations
aurora.distributions.links : Link function implementations
aurora.models.glm.fitting : GLM fitting using families
Notes
-----
For detailed mathematical derivations, see REFERENCES.md in the repository root.
The exponential family framework unifies many common distributions and provides
a consistent theory for inference (MLE, Fisher information, asymptotic normality).
"""
from __future__ import annotations
from abc import ABC, abstractmethod
from ..core.types import Array, Scalar
[docs]
class LinkFunction(ABC):
"""Abstract base class for GLM link functions.
A link function g() relates the mean mu to the linear predictor eta:
g(mu) = eta = X^T beta
Subclasses must implement the forward link, its inverse, and its
derivative. The link must be monotonic and differentiable.
See Also
--------
aurora.distributions.links.LogitLink : Canonical link for binomial/beta
aurora.distributions.links.LogLink : Canonical link for Poisson/Gamma
aurora.distributions.links.IdentityLink : Canonical link for Gaussian
Examples
--------
>>> from aurora.distributions.links import LogitLink
>>> link = LogitLink()
>>> mu = np.array([0.2, 0.5, 0.8])
>>> eta = link.link(mu) # logit(mu)
>>> link.inverse(eta) # recovers mu
>>> link.derivative(mu) # d(mu)/d(eta)
"""
[docs]
@abstractmethod
def link(self, mu: Array) -> Array:
"""Apply the link function g(mu) to map mean to linear predictor.
Parameters
----------
mu : array
Mean parameter values.
Returns
-------
array
Linear predictor eta = g(mu).
"""
[docs]
@abstractmethod
def inverse(self, eta: Array) -> Array:
"""Apply the inverse link g^{-1}(eta) to map linear predictor to mean.
Parameters
----------
eta : array
Linear predictor values.
Returns
-------
array
Mean parameter mu = g^{-1}(eta).
"""
[docs]
@abstractmethod
def derivative(self, mu: Array) -> Array:
"""Compute the derivative dg/dmu of the link function.
Parameters
----------
mu : array
Mean parameter values at which to evaluate the derivative.
Returns
-------
array
Derivative dg/dmu evaluated at mu.
"""
[docs]
class Family(ABC):
"""Abstract base class for exponential family distributions.
Defines the interface that all distribution families must implement
for use in GLM, GAM, and GAMM fitting. Each subclass provides the
log-likelihood, deviance, variance function, and initialization
logic required by the IRLS fitting algorithm.
Subclasses include GaussianFamily, PoissonFamily, BinomialFamily,
GammaFamily, BetaFamily, InverseGaussianFamily, NegativeBinomialFamily,
StudentTFamily, TweedieFamily, and others.
See Also
--------
aurora.distributions.families : Concrete family implementations
aurora.distributions.links : Link function implementations
aurora.models.glm.fit_glm : GLM fitting using families
Examples
--------
>>> from aurora.distributions.families import GaussianFamily
>>> family = GaussianFamily()
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.variance(mu)
array([1., 1., 1.])
>>> family.log_likelihood(np.array([1.1, 1.9, 3.1]), mu)
np.float64(...)
"""
[docs]
@abstractmethod
def log_likelihood(self, y: Array, mu: Array, **params) -> Scalar:
"""Compute the total log-likelihood of observations given the mean.
Parameters
----------
y : array
Observed response values.
mu : array
Fitted mean values.
**params : dict
Additional family-specific parameters (e.g., dispersion).
Returns
-------
scalar
Sum of log-likelihoods across all observations.
"""
[docs]
@abstractmethod
def deviance(self, y: Array, mu: Array, **params) -> Scalar:
"""Compute the deviance between observed and fitted values.
The deviance measures goodness-of-fit:
D(y; mu) = 2 phi [l(y; y) - l(mu; y)]
Parameters
----------
y : array
Observed response values.
mu : array
Fitted mean values.
**params : dict
Additional family-specific parameters.
Returns
-------
scalar
Total deviance (non-negative, zero when mu = y).
"""
[docs]
@abstractmethod
def variance(self, mu: Array, **params) -> Array:
"""Evaluate the variance function V(mu) at given mean values.
The variance function characterizes the mean-variance relationship:
Var(Y) = phi * V(mu)
Parameters
----------
mu : array
Mean parameter values.
**params : dict
Additional family-specific parameters.
Returns
-------
array
Variance function values (same shape as mu).
"""
[docs]
@abstractmethod
def initialize(self, y: Array) -> Array:
"""Compute starting values for the mean parameter from data.
Used by IRLS to obtain initial mu^(0) before the first iteration.
The returned values must be in the valid range of the distribution.
Parameters
----------
y : array
Observed response values.
Returns
-------
array
Initial mean estimates (same shape as y).
"""
@property
@abstractmethod
def default_link(self) -> LinkFunction:
"""Return the default (canonical) link function for this family.
Returns
-------
LinkFunction
The canonical link function instance.
"""
__all__ = ["Family", "LinkFunction"]