Distribution Families

Built-in distribution families.

Standard exponential family distributions: - GaussianFamily: Normal distribution for continuous data - BinomialFamily: Binomial distribution for binary/proportion data - PoissonFamily: Poisson distribution for count data - GammaFamily: Gamma distribution for positive continuous data

Additional distributions (Phase 5 Milestone 1): - BetaFamily: Beta distribution for proportions in (0, 1) - InverseGaussianFamily: Inverse Gaussian for positive durations

Heavy-tailed and robust distributions (Phase 5.5): - StudentTFamily: Student’s t for robust regression - CauchyFamily: Cauchy distribution (t with df=1) - NegativeBinomialFamily: Negative Binomial for overdispersed counts - TweedieFamily: Tweedie for zero-inflated continuous data

class aurora.distributions.families.GaussianFamily(variance=None, link=None)

Bases: Family

Gaussian (normal) distribution family for GLM.

The Gaussian family is the default for continuous response variables. It uses the identity link function by default and has constant variance V(mu) = variance.

Parameters:

• variance (float or None, default=None) – Known variance parameter. Set to a known value for weighted least squares; estimated from data when None (the default).

• link (LinkFunction, optional) – Link function. Defaults to IdentityLink(), which is the canonical link for the Gaussian family.

Variables:

default_link (LinkFunction) – The link function for this family (IdentityLink by default).

Examples

>>> import numpy as np
>>> from aurora.distributions.families import GaussianFamily
>>> from aurora.models.glm import fit_glm
>>> family = GaussianFamily()
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.variance(mu)
array([1., 1., 1.])
>>> # Use with fit_glm
>>> X = np.random.randn(100, 2)
>>> y = X @ np.array([1.5, -0.5]) + np.random.randn(100)
>>> result = fit_glm(X, y, family=family)

See also

PoissonFamily

For count data.

GammaFamily

For positive continuous data.

BinomialFamily

For binary/proportion data.

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.BinomialFamily(n=1.0, link=None)

Bases: Family

Binomial distribution family for binary and proportion data.

This family follows the R binomial() convention: the response is a proportion in [0, 1] and the mean parameter μ is a probability. The variance function is V(μ) = μ(1 − μ) and the number of trials enters as a prior weight (McCullagh & Nelder 1989, §2.3).

Usage conventions

• Binary (Bernoulli) data: pass 0/1 responses; no weights needed.

• Grouped data: pass observed proportions y = successes / trials and the trial counts as weights. With fit_glm use fit_glm(X, proportions, family=BinomialFamily(), weights=trials). A constant number of trials can also be set via BinomialFamily(n=k), which fit_glm treats as weights=k when no explicit weights are given. When trial counts are supplied (either way), the log-likelihood includes the binomial coefficient term log C(n_i, n_i·y_i) so that log_likelihood_/aic_ match R and statsmodels exactly.

param n:

Constant number of trials per observation. Used as the prior weight in deviance/log_likelihood (and by fit_glm as weights=n) when no explicit weights are passed. Leave at the default 1.0 for binary data or when passing per-observation trial counts through weights.

type n:

float, default=1.0

param link:

Link function. Defaults to LogitLink(), which is the canonical link for the Binomial family. Alternatives include ProbitLink and CLogLogLink.

type link:

LinkFunction, optional

ivar default_link:

The link function for this family (LogitLink by default).

vartype default_link:

LinkFunction

ivar n_trials:

The constant trial count n given at construction.

vartype n_trials:

float

Examples

>>> import numpy as np
>>> from aurora.distributions.families import BinomialFamily
>>> from aurora.models.glm import fit_glm
>>> family = BinomialFamily()
>>> mu = np.array([0.2, 0.5, 0.8])
>>> family.variance(mu)  # V(mu) = mu * (1 - mu)
array([0.16, 0.25, 0.16])
>>> # Use with fit_glm for logistic regression
>>> X = np.random.randn(200, 2)
>>> p = 1.0 / (1.0 + np.exp(-(X @ np.array([1.0, -0.5]))))
>>> y = np.random.binomial(1, p)
>>> result = fit_glm(X, y, family=family)

See also

GaussianFamily

For continuous data.

PoissonFamily

For count data.

BetaFamily

For continuous proportions in (0, 1).

property n_trials: float

Constant number of trials per observation (default 1.0).

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

Compute binomial deviance on the probability scale.

Unit deviance d_i = 2 [y_i log(y_i/μ_i) + (1−y_i) log((1−y_i)/(1−μ_i))] with the convention 0·log 0 = 0; the total deviance is Σ w_i d_i (R convention, McCullagh & Nelder §2.3).

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

Parameters:

• n (float)

• link (LinkFunction | None)

class aurora.distributions.families.PoissonFamily(link=None)

Bases: Family

Poisson distribution family for count data.

The Poisson family models count response variables where the variance equals the mean: V(mu) = mu. The canonical link is the log function.

Parameters:

link (LinkFunction, optional) – Link function. Defaults to LogLink(), which is the canonical link for the Poisson family.

Variables:

default_link (LinkFunction) – The link function for this family (LogLink by default).

Examples

>>> import numpy as np
>>> from aurora.distributions.families import PoissonFamily
>>> from aurora.models.glm import fit_glm
>>> family = PoissonFamily()
>>> mu = np.array([1.0, 2.0, 5.0])
>>> family.variance(mu)  # V(mu) = mu
array([1., 2., 5.])
>>> # Use with fit_glm
>>> X = np.random.randn(200, 2)
>>> y = np.random.poisson(np.exp(X @ np.array([0.5, -0.3])))
>>> result = fit_glm(X, y, family=family)

See also

GaussianFamily

For continuous data.

NegativeBinomialFamily

For overdispersed count data.

BinomialFamily

For binary/proportion data.

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.GammaFamily(shape=1.0, link=None)

Bases: Family

Gamma distribution family for positive continuous data.

The Gamma family is suitable for modeling strictly positive continuous responses such as durations, costs, or sizes. The variance function is V(mu) = mu^2 / shape. The canonical link is the inverse function.

Parameters:

• shape (float, default=1.0) – Shape parameter (also called the dispersion parameter). Must be positive. Controls the relationship between mean and variance: Var(Y) = mu^2 / shape.

• link (LinkFunction, optional) – Link function. Defaults to InverseLink(), which is the canonical link for the Gamma family. Common alternatives include LogLink and IdentityLink.

Variables:

default_link (LinkFunction) – The link function for this family (InverseLink by default).

Raises:

ValueError – If shape is not positive.

Examples

>>> import numpy as np
>>> from aurora.distributions.families import GammaFamily
>>> from aurora.models.glm import fit_glm
>>> family = GammaFamily(shape=2.0)
>>> mu = np.array([1.0, 2.0, 4.0])
>>> family.variance(mu)  # mu^2 / shape
array([0.5, 2. , 8. ])
>>> # Use with fit_glm
>>> X = np.random.randn(200, 2)
>>> y = np.exp(X @ np.array([0.5, -0.3]) + np.random.randn(200))
>>> result = fit_glm(X, y, family=GammaFamily(link=None))

See also

GaussianFamily

For continuous data.

InverseGaussianFamily

For positive skewed data with V(mu) = mu^3.

PoissonFamily

For count data.

property shape: float

Shape parameter (inverse dispersion).

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.BetaFamily(phi=1.0, link=None)

Bases: Family

Beta distribution for modeling continuous proportions in (0, 1).

Suitable for modeling success rates, conversion rates, market shares, and any continuous outcome bounded in (0, 1). The distribution uses a mean-precision parameterization where the precision parameter phi controls dispersion independently of the mean.

Parameters:

• phi (float or 'estimate', default=1.0) – Precision parameter (inverse dispersion). - phi large: Low variance, concentrated around mean - phi small: High variance, dispersed - ‘estimate’: Estimate from data via method-of-moments

• link (LinkFunction or None, default=None) – Link function. If None, uses LogitLink (canonical). Also supports ProbitLink, CLogLogLink.

Variables:

• phi (float or str) – Precision parameter or ‘estimate’ for automatic estimation.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.BinomialFamily

Binomial for discrete proportions

aurora.distributions.links.LogitLink

Canonical link for Beta regression

aurora.distributions.links.ProbitLink

Alternative link for Beta regression

Examples

>>> import numpy as np
>>> from aurora.distributions.families import BetaFamily
>>> family = BetaFamily(phi=10.0)
>>> mu = np.array([0.3, 0.5, 0.7])
>>> family.variance(mu)  # μ(1-μ)/(φ+1)
array([0.01909091, 0.02272727, 0.01909091])

>>> # Method-of-moments phi estimation
>>> y = np.random.beta(3, 7, size=100)  # True mean ≈ 0.3
>>> family = BetaFamily(phi='estimate')
>>> mu_init = family.initialize(y)

Notes

The precision phi controls the dispersion: - phi = 1: Uniform(0,1) when mu = 0.5 - phi -> inf: Concentration at mu (degenerate) - phi -> 0: Bimodal distribution at 0 and 1

The variance function is V(mu) = mu(1-mu), same as binomial, but the actual variance also depends on phi: Var(Y) = V(mu)/(phi+1).

__init__(phi=1.0, link=None)

Initialize Beta family.

Parameters:

• phi (float or 'estimate', default=1.0) – Precision parameter. Use ‘estimate’ for method-of-moments.

• link (LinkFunction or None) – Link function (default: LogitLink).

Return type:

None

property phi: float | str

Return the precision parameter.

log_likelihood(y, mu, **params)

Compute log-likelihood for Beta distribution.

Parameters:

• y (array) – Observed proportions in (0, 1).

• mu (array) – Mean parameter in (0, 1).

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Sum of log-likelihoods.

Return type:

scalar

Notes

Log-likelihood:

ℓ = log Γ(φ) - log Γ(μφ) - log Γ((1-μ)φ)

• (μφ - 1) log(y) + ((1-μ)φ - 1) log(1-y)

deviance(y, mu, **params)

Compute deviance for Beta distribution.

For Beta regression, the deviance is computed as:

D = 2 × φ × Σ[y log(y/μ) + (1-y) log((1-y)/(1-μ))]

This is the weighted KL-divergence scaled by 2φ.

Parameters:

• y (array) – Observed proportions.

• mu (array) – Fitted mean values.

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Deviance value.

Return type:

scalar

Notes

For the Beta distribution, deviance is always non-negative. When y = μ, the deviance is zero.

variance(mu, **params)

Compute variance function V(μ) = μ(1-μ).

Note: Actual variance is V(μ)/(φ+1), but for GLM/IRLS we return the variance function without dispersion.

Parameters:

• mu (array) – Mean parameter in (0, 1).

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Variance function values.

Return type:

array

Notes

The true variance is Var(Y) = μ(1-μ)/(φ+1), but in the GLM framework we separate the variance function V(μ) from the dispersion parameter. Here we return V(μ) = μ(1-μ).

full_variance(mu, **params)

Compute full variance: Var(Y) = μ(1-μ)/(φ+1).

Unlike variance(), this includes the precision parameter.

Parameters:

• mu (array) – Mean parameter.

• **params (dict) – Optional ‘phi’ parameter.

Returns:

Full variance values.

Return type:

array

initialize(y)

Initialize mean parameter from data.

Parameters:

y (array) – Observed proportions.

Returns:

Initial mean values, clamped to (0.01, 0.99).

Return type:

array

Notes

Uses sample mean clamped away from boundaries to avoid numerical issues with link function.

estimate_phi(y, mu=None)

Estimate precision parameter from data.

Uses method-of-moments estimation based on the sample variance:

φ̂ = ȳ(1 − ȳ) / Var(y) − 1

Note

Limitation: this is the marginal MoM estimator — it uses the sample mean/variance of y and ignores the fitted means mu (kept for API compatibility only). With heterogeneous means it is biased; a proper conditional estimator would solve Σ (yᵢ−μ̂ᵢ)² / [μ̂ᵢ(1−μ̂ᵢ)] = n / (φ+1) style equations.

Parameters:

• y (array) – Observed proportions.

• mu (array, optional) – Fitted means. Currently unused (see note above).

Returns:

Estimated precision parameter φ.

Return type:

float

property default_link: LinkFunction

Return the default link function (Logit).

class aurora.distributions.families.InverseGaussianFamily(lambda_=1.0, link=None)

Bases: Family

Inverse Gaussian (Wald) distribution for positive continuous data.

Suitable for modeling time-to-event data (failure times), reaction times, right-skewed positive continuous data, and insurance claim amounts. The distribution arises naturally as the first passage time of Brownian motion with positive drift.

Parameters:

• lambda (float or 'estimate', default=1.0) – Shape parameter (controls dispersion). - lambda large: Low variance, concentrated - lambda small: High variance, dispersed - ‘estimate’: Estimate from data via method-of-moments

• link (LinkFunction or None, default=None) – Link function. If None, uses InverseSquareLink (canonical). Common alternatives: LogLink, InverseLink, IdentityLink.

• lambda_ (float | str)

Variables:

• lambda (float or str) – Shape parameter or ‘estimate’ for automatic estimation.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.GammaFamily

Gamma for positive continuous data

aurora.distributions.families.TweedieFamily

Tweedie with power=3 is inverse Gaussian

aurora.distributions.links.InverseSquareLink

Canonical link for inverse Gaussian

Examples

>>> import numpy as np
>>> from aurora.distributions.families import InverseGaussianFamily
>>> from aurora.distributions.links import LogLink

>>> # Basic usage with canonical link
>>> family = InverseGaussianFamily(lambda_=2.0)
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.variance(mu)  # V(mu) = mu^3
array([ 1.,  8., 27.])

>>> # With log link (more common in practice)
>>> family = InverseGaussianFamily(lambda_=2.0, link=LogLink())

Notes

The variance function is V(mu) = mu^3, which means variance increases rapidly with the mean. This makes the inverse Gaussian suitable for data where larger values are more variable.

The canonical link (inverse square) is rarely used in practice; the log or inverse link are more common choices.

For failure time data, the inverse Gaussian arises naturally as the distribution of the first passage time for Brownian motion with drift.

__init__(lambda_=1.0, link=None)

Initialize Inverse Gaussian family.

Parameters:

• lambda (float or 'estimate', default=1.0) – Shape parameter. Use ‘estimate’ for method-of-moments.

• link (LinkFunction or None) – Link function (default: InverseSquareLink).

• lambda_ (float | str)

Return type:

None

property lambda_: float | str

Return the shape parameter.

log_likelihood(y, mu, **params)

Compute log-likelihood for Inverse Gaussian distribution.

Parameters:

• y (array) – Observed positive values.

• mu (array) – Mean parameter (positive).

• **params (dict) – Optional ‘lambda_’ parameter override.

Returns:

Sum of log-likelihoods.

Return type:

scalar

Notes

Log-likelihood:

ℓ = 0.5 × [log(λ) - log(2π) - 3log(y)] - λ(y-μ)² / (2μ²y)

deviance(y, mu, **params)

Compute deviance for Inverse Gaussian distribution.

The unit deviance is:

d(y, μ) = (y - μ)² / (μ² y)

Total deviance: D = Σ d(y, μ)

Parameters:

• y (array) – Observed values.

• mu (array) – Fitted mean values.

• **params (dict) – Optional parameters (not used for deviance).

Returns:

Deviance value.

Return type:

scalar

Notes

This is the “scaled deviance” without λ. The full deviance would be λ × D.

variance(mu, **params)

Compute variance function V(μ) = μ³.

Parameters:

• mu (array) – Mean parameter (positive).

• **params (dict) – Optional parameters (not used).

Returns:

Variance function values.

Return type:

array

Notes

The true variance is Var(Y) = μ³/λ, but the variance function (without dispersion) is V(μ) = μ³.

full_variance(mu, **params)

Compute full variance: Var(Y) = μ³/λ.

Unlike variance(), this includes the shape parameter.

Parameters:

• mu (array) – Mean parameter.

• **params (dict) – Optional ‘lambda_’ parameter.

Returns:

Full variance values.

Return type:

array

initialize(y)

Initialize mean parameter from data.

Parameters:

y (array) – Observed positive values.

Returns:

Initial mean values (sample mean, broadcast to array shape).

Return type:

array

Notes

Uses the sample mean as initial value, which is a natural choice for positive data.

estimate_lambda(y, mu)

Estimate shape parameter from data and fitted values.

Parameters:

• y (array) – Observed values.

• mu (array) – Fitted mean values.

Returns:

Estimated shape parameter λ.

Return type:

float

property default_link: LinkFunction

Return the default link function (InverseSquare).

aurora.distributions.families.WaldFamily

alias of InverseGaussianFamily

class aurora.distributions.families.StudentTFamily(df=5.0, link='identity')

Bases: Family

Student’s t distribution family for robust regression.

Provides robust regression that is less sensitive to outliers compared to Gaussian regression. The degree of robustness is controlled by the degrees of freedom parameter – lower values produce heavier tails and greater downweighting of extreme observations.

Parameters:

• df (float, default=5.0) – Degrees of freedom (nu). Controls tail heaviness: - df=1: Cauchy distribution (very heavy tails, no moments) - df=3-5: Heavy tails, robust to outliers - df=10: Moderate tails - df -> inf: Converges to Gaussian

• link (str, default='identity') – Link function: ‘identity’ or ‘log’ (for positive responses).

Variables:

• df (float) – Degrees of freedom parameter controlling tail heaviness.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.CauchyFamily

Special case with df=1

aurora.distributions.families.GaussianFamily

Limiting case as df -> inf

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The t-distribution is NOT in the exponential family, so standard IRLS does not apply directly. We use iteratively reweighted least squares with adaptive weights that downweight outliers:

w_i = (nu + 1) / (nu + z_i^2)

where z_i = (y_i - mu_i) / sigma is the standardized residual.

The variance is nu / (nu - 2) for nu > 2. For nu <= 2, variance is infinite.

Note

Backend support: this family is NumPy-only — its methods call numpy/scipy directly and do not use the backend namespace abstraction. Passing PyTorch tensors or JAX arrays is not supported.

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import StudentTFamily
>>> # Robust regression with df=5
>>> result = fit_glm(X, y, family=StudentTFamily(df=5))

>>> # Very robust (heavy tails)
>>> result = fit_glm(X, y, family=StudentTFamily(df=3))

References

[1]

Lange et al. (1989). Robust statistical modeling using the t distribution.

name = 'student_t'

valid_y_range = (-inf, inf)

__init__(df=5.0, link='identity')

Initialize Student’s t family.

Parameters:

• df (float) – Degrees of freedom. Must be > 0.

• link (str) – Link function name.

property default_link

Default link function.

variance(mu, **params)

Variance function.

For Student’s t with df > 2:

Var(Y) = σ² × df/(df-2)

Parameters:

• mu (ndarray) – Mean values (not used, variance is constant for location-scale t)

• **params (dict) – May contain ‘scale’ parameter

Returns:

Variance at each observation

Return type:

ndarray

initialize(y)

Initialize mean using median (robust).

Parameters:

y (ndarray) – Response values

Returns:

Initial mean estimates

Return type:

ndarray

log_likelihood(y, mu, **params)

Log-likelihood for t-distribution.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean (location) parameter

• **params (dict) – Must contain ‘scale’ or uses default

Returns:

Total log-likelihood

Return type:

float

deviance(y, mu, **params)

Deviance for t-distribution.

The deviance is -2 × (log-likelihood - saturated log-likelihood). For the t-distribution, we use the saturated model with μ = y.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted mean values

• **params (dict) – Additional parameters

Returns:

Total deviance

Return type:

float

weights(y, mu, **params)

Compute IRLS weights for robust regression.

The t-distribution uses adaptive weights that downweight outliers:

w_i = (ν + 1) / (ν + z_i²)

where z_i = (y_i - μ_i) / σ is the standardized residual.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Current mean estimates

• **params (dict) – Must contain ‘scale’

Returns:

Weights for each observation

Return type:

ndarray

estimate_scale(y, mu, ddof=1)

Estimate scale parameter robustly.

Uses the Median Absolute Deviation (MAD) for robustness.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted mean values

• ddof (int) – Degrees of freedom adjustment

Returns:

Estimated scale parameter

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of log-likelihood w.r.t. μ.

For Laplace approximation and gradient-based optimization.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Gradient of log-likelihood

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of log-likelihood w.r.t. μ.

For Hessian computation in Laplace approximation.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Negative Hessian (Fisher information approximation)

Return type:

ndarray

__repr__()

String representation.

Return type:

str

class aurora.distributions.families.CauchyFamily(link='identity')

Bases: StudentTFamily

Cauchy distribution (Student’s t with df=1).

The Cauchy distribution has extremely heavy tails and no finite moments. It is maximally robust to outliers but may have convergence issues. Use with caution – suitable only when extreme robustness is needed.

Parameters:

link (str, default='identity') – Link function: ‘identity’ or ‘log’.

Variables:

• df (float) – Fixed at 1.0 for Cauchy.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.StudentTFamily

General t-distribution family

Examples

>>> from aurora.models.glm import fit_glm
>>> result = fit_glm(X, y, family='cauchy')

name = 'cauchy'

__init__(link='identity')

Initialize Cauchy family (t with df=1).

Parameters:

link (str)

variance(mu, **params)

Variance is undefined (infinite) for Cauchy.

Parameters:

mu (NDArray)

Return type:

NDArray

__repr__()

String representation.

Return type:

str

class aurora.distributions.families.NegativeBinomialFamily(theta=1.0, link='log')

Bases: Family

Negative Binomial (NB2) distribution for overdispersed count data.

Models count data where the variance exceeds the mean: Var(Y) = mu + mu^2 / theta. This is an extension of Poisson regression for overdispersed data, commonly arising from unobserved heterogeneity, clustered or hierarchical data, and temporal or spatial correlation.

Parameters:

• theta (float or str, default=1.0) – Dispersion parameter (also called ‘size’ or ‘k’). - theta -> inf: Converges to Poisson - theta small: High overdispersion - ‘estimate’: Estimate from data during fitting

• link (str, default='log') – Link function: ‘log’ (default), ‘identity’, or ‘sqrt’.

Variables:

• theta (float) – Dispersion parameter (accessed via property; raises ValueError if not yet estimated when theta=’estimate’ was specified).

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.PoissonFamily

Poisson for equidispersed counts

aurora.distributions.families.ZeroInflatedPoisson

ZIP for excess zeros

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The NB2 parameterization is used:

P(Y=y) = Gamma(y+theta) / [Gamma(theta) y!] x (theta/(theta+mu))^theta x (mu/(theta+mu))^y

Properties:

E[Y] = mu Var(Y) = mu + mu^2 / theta

The negative binomial is a member of the exponential family when theta is fixed.

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import NegativeBinomialFamily
>>> # Fixed dispersion
>>> result = fit_glm(X, y, family=NegativeBinomialFamily(theta=2.0))

>>> # Estimate dispersion from data (glm.nb two-step flow)
>>> result = fit_glm(X, y, family=NegativeBinomialFamily(theta='estimate'))
>>> result.theta_  # estimated dispersion

References

[1]

Hilbe (2011). Negative Binomial Regression.

[2]

Cameron & Trivedi (2013). Regression Analysis of Count Data.

name = 'negative_binomial'

valid_y_range = (0, inf)

__init__(theta=1.0, link='log')

Initialize Negative Binomial family.

Parameters:

• theta (float or 'estimate') – Dispersion parameter. If ‘estimate’, will be estimated from data.

• link (str) – Link function name.

property theta: float

Dispersion parameter.

property default_link

Default link function (log).

variance(mu, **params)

Variance function: V(μ) = μ + μ²/θ.

Parameters:

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘theta’ override

Returns:

Variance at each observation

Return type:

ndarray

initialize(y)

Initialize mean with sample mean + small constant.

Parameters:

y (ndarray) – Response values (counts)

Returns:

Initial mean estimates

Return type:

ndarray

log_likelihood(y, mu, **params)

Log-likelihood for Negative Binomial.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean (fitted) values

• **params (dict) – May contain ‘theta’

Returns:

Total log-likelihood

Return type:

float

deviance(y, mu, **params)

Deviance for Negative Binomial.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Fitted mean values

• **params (dict) – May contain ‘theta’

Returns:

Total deviance

Return type:

float

estimate_theta(y, mu, method='ml')

Estimate dispersion parameter theta.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Fitted means

• method (str, default='ml') – Estimation method: - ‘moments’: Method of moments (fast) - ‘ml’: Maximum likelihood (more accurate)

Returns:

Estimated theta

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean values

Returns:

Gradient

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean values

Returns:

Negative Hessian diagonal

Return type:

ndarray

__repr__()

String representation.

Return type:

str

aurora.distributions.families.NegBinFamily

alias of NegativeBinomialFamily

class aurora.distributions.families.TweedieFamily(power=1.5, link='log', phi=1.0)

Bases: Family

Tweedie distribution for zero-inflated positive continuous data.

A compound Poisson-Gamma distribution that naturally handles data with a point mass at zero (no events) and continuous positive values (sum of event magnitudes). The variance function is Var(Y) = phi * mu^p, where the power parameter p controls the shape.

Special cases of the Tweedie family (for reference): - p = 0: Gaussian - p = 1: Poisson - 1 < p < 2: Compound Poisson-Gamma (Tweedie proper) - p = 2: Gamma - p = 3: Inverse Gaussian

This class implements only the compound Poisson-Gamma case 1 < p < 2; for the boundary cases use the dedicated families (GaussianFamily, PoissonFamily, GammaFamily, InverseGaussianFamily).

Parameters:

• power (float, default=1.5) – Variance power parameter p. For compound Poisson-Gamma, must be in (1, 2). Common choices: - p = 1.5: Balanced, works well for many applications - p = 1.6-1.7: More right-skewed, common for insurance claims - p closer to 1: More mass at zero - p closer to 2: Less mass at zero, more Gamma-like

• link (str, default='log') – Link function: ‘log’ (default), ‘identity’, or ‘power’.

• phi (float, default=1.0) – Dispersion parameter (scale).

Variables:

• power (float) – Variance power parameter in (1, 2).

• phi (float) – Dispersion (scale) parameter.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.CompoundPoissonGammaFamily

Explicit compound Poisson-Gamma alias for TweedieFamily

aurora.distributions.families.GammaFamily

Tweedie with power=2

aurora.distributions.families.InverseGaussianFamily

Tweedie with power=3

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The Tweedie distribution can be viewed as: - N ~ Poisson(lambda) events occur - Each event has magnitude X_i ~ Gamma(alpha, beta) - Total Y = X_1 + … + X_N

This naturally creates exact zeros (N=0) and continuous positives (N>0).

The probability of zero is: P(Y=0) = exp(-lambda)

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import TweedieFamily
>>> # Insurance claims (many zeros, heavy right tail)
>>> result = fit_glm(X, claims, family=TweedieFamily(power=1.7))

>>> # Rainfall amounts (zeros for dry days)
>>> result = fit_glm(X, rainfall, family=TweedieFamily(power=1.5))

References

[1]

Jorgensen (1987). Exponential dispersion models.

[2]

Smyth & Jorgensen (2002). Fitting Tweedie’s compound Poisson model.

name = 'tweedie'

valid_y_range = (0, inf)

__init__(power=1.5, link='log', phi=1.0)

Initialize Tweedie family.

Parameters:

• power (float) – Variance power parameter. Must be in (1, 2) for proper Tweedie.

• link (str) – Link function name.

• phi (float) – Dispersion parameter (scale).

property default_link

Default link function (log).

variance(mu, **params)

Variance function: V(μ) = μ^p.

Parameters:

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘phi’ (dispersion)

Returns:

Variance at each observation (without phi factor)

Return type:

ndarray

initialize(y)

Initialize mean using positive values.

Parameters:

y (ndarray) – Response values

Returns:

Initial mean estimates

Return type:

ndarray

deviance(y, mu, **params)

Deviance for Tweedie distribution.

The unit deviance for Tweedie with power p is:

For p ∈ (1, 2):

d(y, μ) = 2 × [y^(2-p)/((1-p)(2-p)) - y×μ^(1-p)/(1-p) + μ^(2-p)/(2-p)]

For y = 0:

d(0, μ) = 2 × μ^(2-p) / (2-p)

Parameters:

• y (ndarray) – Observed values (including zeros)

• mu (ndarray) – Fitted mean values

• **params (dict) – Additional parameters

Returns:

Total deviance

Return type:

float

log_likelihood(y, mu, **params)

Approximate (quasi-)log-likelihood for Tweedie.

Warning

The Tweedie density has no closed form. This method returns the quasi-likelihood -D/(2φ), which drops the normalization constant c(y, φ, p) of the full density (the series evaluation of Dunn & Smyth 2005 is not implemented). As a consequence, the returned value — and any AIC/BIC derived from it — is not comparable across families, nor across different values of the power parameter p or the dispersion φ. It is only meaningful for comparing models that share the same family, power and dispersion (e.g. nested mean models).

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘phi’

Returns:

Approximate (quasi-)log-likelihood

Return type:

float

weights(y, mu, **params)

IRLS weights for Tweedie.

Weights are 1/V(μ) = 1/μ^p.

Parameters:

• y (ndarray) – Observed values (not used directly)

• mu (ndarray) – Current mean estimates

Returns:

IRLS weights

Return type:

ndarray

estimate_power(y, mu, power_range=(1.1, 1.9), n_grid=20)

Estimate the power parameter via a deviance grid search.

Note

This is a heuristic: it minimizes the unit deviance over a grid of powers. It is not a profile likelihood in p — the full Tweedie density (Dunn & Smyth 2005 series evaluation) is not implemented, and the quasi-likelihood ignores the normalization constant that depends on p, so the criterion is not a valid profile likelihood for the power parameter. Treat the result as a rough starting value, not an MLE.

The family’s power attribute is always restored to its original value before returning.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted means (from current model)

• power_range (tuple) – Range of powers to search

• n_grid (int) – Number of grid points

Returns:

Estimated power parameter

Return type:

float

estimate_phi(y, mu, ddof=1)

Estimate dispersion parameter.

Uses Pearson estimator: φ = (1/n) Σ (y-μ)² / V(μ)

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted means

• ddof (int) – Degrees of freedom adjustment

Returns:

Estimated dispersion

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of quasi-log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Gradient

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of quasi-log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Negative Hessian diagonal (Fisher information)

Return type:

ndarray

probability_zero(mu, **params)

Probability of observing exactly zero.

For Tweedie with 1 < p < 2:

P(Y=0) = exp(-λ)

where λ = μ^(2-p) / [φ(2-p)]

Parameters:

mu (ndarray) – Mean values

Returns:

Probability of zero at each observation

Return type:

ndarray

__repr__()

String representation.

Return type:

str

class aurora.distributions.families.CompoundPoissonGammaFamily(power=1.5, link='log', phi=1.0)

Bases: TweedieFamily

Compound Poisson-Gamma distribution, an alias for TweedieFamily.

Identical to TweedieFamily but with more descriptive naming for users familiar with the compound Poisson-Gamma formulation.

The model assumes: - N ~ Poisson(lambda) events occur - Each event has magnitude G_i ~ Gamma(alpha, beta) - Total Y = G_1 + … + G_N (or 0 if N=0)

Parameters are related to Tweedie by: - lambda = mu^(2-p) / [phi(2-p)] (Poisson rate) - alpha = (2-p)/(p-1) (Gamma shape) - beta = phi(p-1)mu^(p-1) (Gamma rate)

Parameters:

• power (float, default=1.5) – Variance power parameter p. Must be in (1, 2).

• link (str, default='log') – Link function name.

• phi (float, default=1.0) – Dispersion parameter (scale).

Variables:

• power (float) – Variance power parameter in (1, 2).

• phi (float) – Dispersion (scale) parameter.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.TweedieFamily

Parent class with identical behavior

Examples

>>> from aurora.distributions.families import CompoundPoissonGammaFamily
>>> family = CompoundPoissonGammaFamily(power=1.5, phi=1.0)
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.get_poisson_rate(mu)       # lambda
array([...])
>>> family.get_gamma_shape()           # alpha
1.0

name = 'compound_poisson_gamma'

get_poisson_rate(mu, **params)

Get the underlying Poisson rate λ.

Parameters:

mu (ndarray) – Mean values

Returns:

Poisson rate for event count

Return type:

ndarray

get_gamma_shape()

Get the Gamma shape parameter α.

Return type:

float

get_gamma_rate(mu, **params)

Get the Gamma rate parameter β.

Parameters:

mu (ndarray) – Mean values

Returns:

Gamma rate for event magnitude

Return type:

ndarray

class aurora.distributions.families.gaussian.GaussianFamily(variance=None, link=None)

Bases: Family

Gaussian (normal) distribution family for GLM.

The Gaussian family is the default for continuous response variables. It uses the identity link function by default and has constant variance V(mu) = variance.

Parameters:

• variance (float or None, default=None) – Known variance parameter. Set to a known value for weighted least squares; estimated from data when None (the default).

• link (LinkFunction, optional) – Link function. Defaults to IdentityLink(), which is the canonical link for the Gaussian family.

Variables:

default_link (LinkFunction) – The link function for this family (IdentityLink by default).

Examples

>>> import numpy as np
>>> from aurora.distributions.families import GaussianFamily
>>> from aurora.models.glm import fit_glm
>>> family = GaussianFamily()
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.variance(mu)
array([1., 1., 1.])
>>> # Use with fit_glm
>>> X = np.random.randn(100, 2)
>>> y = X @ np.array([1.5, -0.5]) + np.random.randn(100)
>>> result = fit_glm(X, y, family=family)

See also

PoissonFamily

For count data.

GammaFamily

For positive continuous data.

BinomialFamily

For binary/proportion data.

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.poisson.PoissonFamily(link=None)

Bases: Family

Poisson distribution family for count data.

The Poisson family models count response variables where the variance equals the mean: V(mu) = mu. The canonical link is the log function.

Parameters:

link (LinkFunction, optional) – Link function. Defaults to LogLink(), which is the canonical link for the Poisson family.

Variables:

default_link (LinkFunction) – The link function for this family (LogLink by default).

Examples

>>> import numpy as np
>>> from aurora.distributions.families import PoissonFamily
>>> from aurora.models.glm import fit_glm
>>> family = PoissonFamily()
>>> mu = np.array([1.0, 2.0, 5.0])
>>> family.variance(mu)  # V(mu) = mu
array([1., 2., 5.])
>>> # Use with fit_glm
>>> X = np.random.randn(200, 2)
>>> y = np.random.poisson(np.exp(X @ np.array([0.5, -0.3])))
>>> result = fit_glm(X, y, family=family)

See also

GaussianFamily

For continuous data.

NegativeBinomialFamily

For overdispersed count data.

BinomialFamily

For binary/proportion data.

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.binomial.BinomialFamily(n=1.0, link=None)

Bases: Family

Binomial distribution family for binary and proportion data.

This family follows the R binomial() convention: the response is a proportion in [0, 1] and the mean parameter μ is a probability. The variance function is V(μ) = μ(1 − μ) and the number of trials enters as a prior weight (McCullagh & Nelder 1989, §2.3).

Usage conventions

• Binary (Bernoulli) data: pass 0/1 responses; no weights needed.

• Grouped data: pass observed proportions y = successes / trials and the trial counts as weights. With fit_glm use fit_glm(X, proportions, family=BinomialFamily(), weights=trials). A constant number of trials can also be set via BinomialFamily(n=k), which fit_glm treats as weights=k when no explicit weights are given. When trial counts are supplied (either way), the log-likelihood includes the binomial coefficient term log C(n_i, n_i·y_i) so that log_likelihood_/aic_ match R and statsmodels exactly.

param n:

Constant number of trials per observation. Used as the prior weight in deviance/log_likelihood (and by fit_glm as weights=n) when no explicit weights are passed. Leave at the default 1.0 for binary data or when passing per-observation trial counts through weights.

type n:

float, default=1.0

param link:

Link function. Defaults to LogitLink(), which is the canonical link for the Binomial family. Alternatives include ProbitLink and CLogLogLink.

type link:

LinkFunction, optional

ivar default_link:

The link function for this family (LogitLink by default).

vartype default_link:

LinkFunction

ivar n_trials:

The constant trial count n given at construction.

vartype n_trials:

float

Examples

>>> import numpy as np
>>> from aurora.distributions.families import BinomialFamily
>>> from aurora.models.glm import fit_glm
>>> family = BinomialFamily()
>>> mu = np.array([0.2, 0.5, 0.8])
>>> family.variance(mu)  # V(mu) = mu * (1 - mu)
array([0.16, 0.25, 0.16])
>>> # Use with fit_glm for logistic regression
>>> X = np.random.randn(200, 2)
>>> p = 1.0 / (1.0 + np.exp(-(X @ np.array([1.0, -0.5]))))
>>> y = np.random.binomial(1, p)
>>> result = fit_glm(X, y, family=family)

See also

GaussianFamily

For continuous data.

PoissonFamily

For count data.

BetaFamily

For continuous proportions in (0, 1).

property n_trials: float

Constant number of trials per observation (default 1.0).

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

Compute binomial deviance on the probability scale.

Unit deviance d_i = 2 [y_i log(y_i/μ_i) + (1−y_i) log((1−y_i)/(1−μ_i))] with the convention 0·log 0 = 0; the total deviance is Σ w_i d_i (R convention, McCullagh & Nelder §2.3).

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

Parameters:

• n (float)

• link (LinkFunction | None)

class aurora.distributions.families.gamma.GammaFamily(shape=1.0, link=None)

Bases: Family

Gamma distribution family for positive continuous data.

The Gamma family is suitable for modeling strictly positive continuous responses such as durations, costs, or sizes. The variance function is V(mu) = mu^2 / shape. The canonical link is the inverse function.

Parameters:

• shape (float, default=1.0) – Shape parameter (also called the dispersion parameter). Must be positive. Controls the relationship between mean and variance: Var(Y) = mu^2 / shape.

• link (LinkFunction, optional) – Link function. Defaults to InverseLink(), which is the canonical link for the Gamma family. Common alternatives include LogLink and IdentityLink.

Variables:

default_link (LinkFunction) – The link function for this family (InverseLink by default).

Raises:

ValueError – If shape is not positive.

Examples

>>> import numpy as np
>>> from aurora.distributions.families import GammaFamily
>>> from aurora.models.glm import fit_glm
>>> family = GammaFamily(shape=2.0)
>>> mu = np.array([1.0, 2.0, 4.0])
>>> family.variance(mu)  # mu^2 / shape
array([0.5, 2. , 8. ])
>>> # Use with fit_glm
>>> X = np.random.randn(200, 2)
>>> y = np.exp(X @ np.array([0.5, -0.3]) + np.random.randn(200))
>>> result = fit_glm(X, y, family=GammaFamily(link=None))

See also

GaussianFamily

For continuous data.

InverseGaussianFamily

For positive skewed data with V(mu) = mu^3.

PoissonFamily

For count data.

property shape: float

Shape parameter (inverse dispersion).

log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction

class aurora.distributions.families.beta.BetaFamily(phi=1.0, link=None)

Bases: Family

Beta distribution for modeling continuous proportions in (0, 1).

Suitable for modeling success rates, conversion rates, market shares, and any continuous outcome bounded in (0, 1). The distribution uses a mean-precision parameterization where the precision parameter phi controls dispersion independently of the mean.

Parameters:

• phi (float or 'estimate', default=1.0) – Precision parameter (inverse dispersion). - phi large: Low variance, concentrated around mean - phi small: High variance, dispersed - ‘estimate’: Estimate from data via method-of-moments

• link (LinkFunction or None, default=None) – Link function. If None, uses LogitLink (canonical). Also supports ProbitLink, CLogLogLink.

Variables:

• phi (float or str) – Precision parameter or ‘estimate’ for automatic estimation.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.BinomialFamily

Binomial for discrete proportions

aurora.distributions.links.LogitLink

Canonical link for Beta regression

aurora.distributions.links.ProbitLink

Alternative link for Beta regression

Examples

>>> import numpy as np
>>> from aurora.distributions.families import BetaFamily
>>> family = BetaFamily(phi=10.0)
>>> mu = np.array([0.3, 0.5, 0.7])
>>> family.variance(mu)  # μ(1-μ)/(φ+1)
array([0.01909091, 0.02272727, 0.01909091])

>>> # Method-of-moments phi estimation
>>> y = np.random.beta(3, 7, size=100)  # True mean ≈ 0.3
>>> family = BetaFamily(phi='estimate')
>>> mu_init = family.initialize(y)

Notes

The precision phi controls the dispersion: - phi = 1: Uniform(0,1) when mu = 0.5 - phi -> inf: Concentration at mu (degenerate) - phi -> 0: Bimodal distribution at 0 and 1

The variance function is V(mu) = mu(1-mu), same as binomial, but the actual variance also depends on phi: Var(Y) = V(mu)/(phi+1).

__init__(phi=1.0, link=None)

Initialize Beta family.

Parameters:

• phi (float or 'estimate', default=1.0) – Precision parameter. Use ‘estimate’ for method-of-moments.

• link (LinkFunction or None) – Link function (default: LogitLink).

Return type:

None

property phi: float | str

Return the precision parameter.

log_likelihood(y, mu, **params)

Compute log-likelihood for Beta distribution.

Parameters:

• y (array) – Observed proportions in (0, 1).

• mu (array) – Mean parameter in (0, 1).

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Sum of log-likelihoods.

Return type:

scalar

Notes

Log-likelihood:

ℓ = log Γ(φ) - log Γ(μφ) - log Γ((1-μ)φ)

• (μφ - 1) log(y) + ((1-μ)φ - 1) log(1-y)

deviance(y, mu, **params)

Compute deviance for Beta distribution.

For Beta regression, the deviance is computed as:

D = 2 × φ × Σ[y log(y/μ) + (1-y) log((1-y)/(1-μ))]

This is the weighted KL-divergence scaled by 2φ.

Parameters:

• y (array) – Observed proportions.

• mu (array) – Fitted mean values.

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Deviance value.

Return type:

scalar

Notes

For the Beta distribution, deviance is always non-negative. When y = μ, the deviance is zero.

variance(mu, **params)

Compute variance function V(μ) = μ(1-μ).

Note: Actual variance is V(μ)/(φ+1), but for GLM/IRLS we return the variance function without dispersion.

Parameters:

• mu (array) – Mean parameter in (0, 1).

• **params (dict) – Optional ‘phi’ parameter override.

Returns:

Variance function values.

Return type:

array

Notes

The true variance is Var(Y) = μ(1-μ)/(φ+1), but in the GLM framework we separate the variance function V(μ) from the dispersion parameter. Here we return V(μ) = μ(1-μ).

full_variance(mu, **params)

Compute full variance: Var(Y) = μ(1-μ)/(φ+1).

Unlike variance(), this includes the precision parameter.

Parameters:

• mu (array) – Mean parameter.

• **params (dict) – Optional ‘phi’ parameter.

Returns:

Full variance values.

Return type:

array

initialize(y)

Initialize mean parameter from data.

Parameters:

y (array) – Observed proportions.

Returns:

Initial mean values, clamped to (0.01, 0.99).

Return type:

array

Notes

Uses sample mean clamped away from boundaries to avoid numerical issues with link function.

estimate_phi(y, mu=None)

Estimate precision parameter from data.

Uses method-of-moments estimation based on the sample variance:

φ̂ = ȳ(1 − ȳ) / Var(y) − 1

Note

Limitation: this is the marginal MoM estimator — it uses the sample mean/variance of y and ignores the fitted means mu (kept for API compatibility only). With heterogeneous means it is biased; a proper conditional estimator would solve Σ (yᵢ−μ̂ᵢ)² / [μ̂ᵢ(1−μ̂ᵢ)] = n / (φ+1) style equations.

Parameters:

• y (array) – Observed proportions.

• mu (array, optional) – Fitted means. Currently unused (see note above).

Returns:

Estimated precision parameter φ.

Return type:

float

property default_link: LinkFunction

Return the default link function (Logit).

class aurora.distributions.families.inverse_gaussian.InverseGaussianFamily(lambda_=1.0, link=None)

Bases: Family

Inverse Gaussian (Wald) distribution for positive continuous data.

Suitable for modeling time-to-event data (failure times), reaction times, right-skewed positive continuous data, and insurance claim amounts. The distribution arises naturally as the first passage time of Brownian motion with positive drift.

Parameters:

• lambda (float or 'estimate', default=1.0) – Shape parameter (controls dispersion). - lambda large: Low variance, concentrated - lambda small: High variance, dispersed - ‘estimate’: Estimate from data via method-of-moments

• link (LinkFunction or None, default=None) – Link function. If None, uses InverseSquareLink (canonical). Common alternatives: LogLink, InverseLink, IdentityLink.

• lambda_ (float | str)

Variables:

• lambda (float or str) – Shape parameter or ‘estimate’ for automatic estimation.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.GammaFamily

Gamma for positive continuous data

aurora.distributions.families.TweedieFamily

Tweedie with power=3 is inverse Gaussian

aurora.distributions.links.InverseSquareLink

Canonical link for inverse Gaussian

Examples

>>> import numpy as np
>>> from aurora.distributions.families import InverseGaussianFamily
>>> from aurora.distributions.links import LogLink

>>> # Basic usage with canonical link
>>> family = InverseGaussianFamily(lambda_=2.0)
>>> mu = np.array([1.0, 2.0, 3.0])
>>> family.variance(mu)  # V(mu) = mu^3
array([ 1.,  8., 27.])

>>> # With log link (more common in practice)
>>> family = InverseGaussianFamily(lambda_=2.0, link=LogLink())

Notes

The variance function is V(mu) = mu^3, which means variance increases rapidly with the mean. This makes the inverse Gaussian suitable for data where larger values are more variable.

The canonical link (inverse square) is rarely used in practice; the log or inverse link are more common choices.

For failure time data, the inverse Gaussian arises naturally as the distribution of the first passage time for Brownian motion with drift.

__init__(lambda_=1.0, link=None)

Initialize Inverse Gaussian family.

Parameters:

• lambda (float or 'estimate', default=1.0) – Shape parameter. Use ‘estimate’ for method-of-moments.

• link (LinkFunction or None) – Link function (default: InverseSquareLink).

• lambda_ (float | str)

Return type:

None

property lambda_: float | str

Return the shape parameter.

log_likelihood(y, mu, **params)

Compute log-likelihood for Inverse Gaussian distribution.

Parameters:

• y (array) – Observed positive values.

• mu (array) – Mean parameter (positive).

• **params (dict) – Optional ‘lambda_’ parameter override.

Returns:

Sum of log-likelihoods.

Return type:

scalar

Notes

Log-likelihood:

ℓ = 0.5 × [log(λ) - log(2π) - 3log(y)] - λ(y-μ)² / (2μ²y)

deviance(y, mu, **params)

Compute deviance for Inverse Gaussian distribution.

The unit deviance is:

d(y, μ) = (y - μ)² / (μ² y)

Total deviance: D = Σ d(y, μ)

Parameters:

• y (array) – Observed values.

• mu (array) – Fitted mean values.

• **params (dict) – Optional parameters (not used for deviance).

Returns:

Deviance value.

Return type:

scalar

Notes

This is the “scaled deviance” without λ. The full deviance would be λ × D.

variance(mu, **params)

Compute variance function V(μ) = μ³.

Parameters:

• mu (array) – Mean parameter (positive).

• **params (dict) – Optional parameters (not used).

Returns:

Variance function values.

Return type:

array

Notes

The true variance is Var(Y) = μ³/λ, but the variance function (without dispersion) is V(μ) = μ³.

full_variance(mu, **params)

Compute full variance: Var(Y) = μ³/λ.

Unlike variance(), this includes the shape parameter.

Parameters:

• mu (array) – Mean parameter.

• **params (dict) – Optional ‘lambda_’ parameter.

Returns:

Full variance values.

Return type:

array

initialize(y)

Initialize mean parameter from data.

Parameters:

y (array) – Observed positive values.

Returns:

Initial mean values (sample mean, broadcast to array shape).

Return type:

array

Notes

Uses the sample mean as initial value, which is a natural choice for positive data.

estimate_lambda(y, mu)

Estimate shape parameter from data and fitted values.

Parameters:

• y (array) – Observed values.

• mu (array) – Fitted mean values.

Returns:

Estimated shape parameter λ.

Return type:

float

property default_link: LinkFunction

Return the default link function (InverseSquare).

class aurora.distributions.families.negative_binomial.NegativeBinomialFamily(theta=1.0, link='log')

Bases: Family

Negative Binomial (NB2) distribution for overdispersed count data.

Models count data where the variance exceeds the mean: Var(Y) = mu + mu^2 / theta. This is an extension of Poisson regression for overdispersed data, commonly arising from unobserved heterogeneity, clustered or hierarchical data, and temporal or spatial correlation.

Parameters:

• theta (float or str, default=1.0) – Dispersion parameter (also called ‘size’ or ‘k’). - theta -> inf: Converges to Poisson - theta small: High overdispersion - ‘estimate’: Estimate from data during fitting

• link (str, default='log') – Link function: ‘log’ (default), ‘identity’, or ‘sqrt’.

Variables:

• theta (float) – Dispersion parameter (accessed via property; raises ValueError if not yet estimated when theta=’estimate’ was specified).

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.PoissonFamily

Poisson for equidispersed counts

aurora.distributions.families.ZeroInflatedPoisson

ZIP for excess zeros

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The NB2 parameterization is used:

P(Y=y) = Gamma(y+theta) / [Gamma(theta) y!] x (theta/(theta+mu))^theta x (mu/(theta+mu))^y

Properties:

E[Y] = mu Var(Y) = mu + mu^2 / theta

The negative binomial is a member of the exponential family when theta is fixed.

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import NegativeBinomialFamily
>>> # Fixed dispersion
>>> result = fit_glm(X, y, family=NegativeBinomialFamily(theta=2.0))

>>> # Estimate dispersion from data (glm.nb two-step flow)
>>> result = fit_glm(X, y, family=NegativeBinomialFamily(theta='estimate'))
>>> result.theta_  # estimated dispersion

References

[1]

Hilbe (2011). Negative Binomial Regression.

[2]

Cameron & Trivedi (2013). Regression Analysis of Count Data.

name = 'negative_binomial'

valid_y_range = (0, inf)

__init__(theta=1.0, link='log')

Initialize Negative Binomial family.

Parameters:

• theta (float or 'estimate') – Dispersion parameter. If ‘estimate’, will be estimated from data.

• link (str) – Link function name.

property theta: float

Dispersion parameter.

property default_link

Default link function (log).

variance(mu, **params)

Variance function: V(μ) = μ + μ²/θ.

Parameters:

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘theta’ override

Returns:

Variance at each observation

Return type:

ndarray

initialize(y)

Initialize mean with sample mean + small constant.

Parameters:

y (ndarray) – Response values (counts)

Returns:

Initial mean estimates

Return type:

ndarray

log_likelihood(y, mu, **params)

Log-likelihood for Negative Binomial.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean (fitted) values

• **params (dict) – May contain ‘theta’

Returns:

Total log-likelihood

Return type:

float

deviance(y, mu, **params)

Deviance for Negative Binomial.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Fitted mean values

• **params (dict) – May contain ‘theta’

Returns:

Total deviance

Return type:

float

estimate_theta(y, mu, method='ml')

Estimate dispersion parameter theta.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Fitted means

• method (str, default='ml') – Estimation method: - ‘moments’: Method of moments (fast) - ‘ml’: Maximum likelihood (more accurate)

Returns:

Estimated theta

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean values

Returns:

Gradient

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed counts

• mu (ndarray) – Mean values

Returns:

Negative Hessian diagonal

Return type:

ndarray

__repr__()

String representation.

Return type:

str

class aurora.distributions.families.student_t.StudentTFamily(df=5.0, link='identity')

Bases: Family

Student’s t distribution family for robust regression.

Provides robust regression that is less sensitive to outliers compared to Gaussian regression. The degree of robustness is controlled by the degrees of freedom parameter – lower values produce heavier tails and greater downweighting of extreme observations.

Parameters:

• df (float, default=5.0) – Degrees of freedom (nu). Controls tail heaviness: - df=1: Cauchy distribution (very heavy tails, no moments) - df=3-5: Heavy tails, robust to outliers - df=10: Moderate tails - df -> inf: Converges to Gaussian

• link (str, default='identity') – Link function: ‘identity’ or ‘log’ (for positive responses).

Variables:

• df (float) – Degrees of freedom parameter controlling tail heaviness.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.CauchyFamily

Special case with df=1

aurora.distributions.families.GaussianFamily

Limiting case as df -> inf

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The t-distribution is NOT in the exponential family, so standard IRLS does not apply directly. We use iteratively reweighted least squares with adaptive weights that downweight outliers:

w_i = (nu + 1) / (nu + z_i^2)

where z_i = (y_i - mu_i) / sigma is the standardized residual.

The variance is nu / (nu - 2) for nu > 2. For nu <= 2, variance is infinite.

Note

Backend support: this family is NumPy-only — its methods call numpy/scipy directly and do not use the backend namespace abstraction. Passing PyTorch tensors or JAX arrays is not supported.

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import StudentTFamily
>>> # Robust regression with df=5
>>> result = fit_glm(X, y, family=StudentTFamily(df=5))

>>> # Very robust (heavy tails)
>>> result = fit_glm(X, y, family=StudentTFamily(df=3))

References

[1]

Lange et al. (1989). Robust statistical modeling using the t distribution.

name = 'student_t'

valid_y_range = (-inf, inf)

__init__(df=5.0, link='identity')

Initialize Student’s t family.

Parameters:

• df (float) – Degrees of freedom. Must be > 0.

• link (str) – Link function name.

property default_link

Default link function.

variance(mu, **params)

Variance function.

For Student’s t with df > 2:

Var(Y) = σ² × df/(df-2)

Parameters:

• mu (ndarray) – Mean values (not used, variance is constant for location-scale t)

• **params (dict) – May contain ‘scale’ parameter

Returns:

Variance at each observation

Return type:

ndarray

initialize(y)

Initialize mean using median (robust).

Parameters:

y (ndarray) – Response values

Returns:

Initial mean estimates

Return type:

ndarray

log_likelihood(y, mu, **params)

Log-likelihood for t-distribution.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean (location) parameter

• **params (dict) – Must contain ‘scale’ or uses default

Returns:

Total log-likelihood

Return type:

float

deviance(y, mu, **params)

Deviance for t-distribution.

The deviance is -2 × (log-likelihood - saturated log-likelihood). For the t-distribution, we use the saturated model with μ = y.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted mean values

• **params (dict) – Additional parameters

Returns:

Total deviance

Return type:

float

weights(y, mu, **params)

Compute IRLS weights for robust regression.

The t-distribution uses adaptive weights that downweight outliers:

w_i = (ν + 1) / (ν + z_i²)

where z_i = (y_i - μ_i) / σ is the standardized residual.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Current mean estimates

• **params (dict) – Must contain ‘scale’

Returns:

Weights for each observation

Return type:

ndarray

estimate_scale(y, mu, ddof=1)

Estimate scale parameter robustly.

Uses the Median Absolute Deviation (MAD) for robustness.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted mean values

• ddof (int) – Degrees of freedom adjustment

Returns:

Estimated scale parameter

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of log-likelihood w.r.t. μ.

For Laplace approximation and gradient-based optimization.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Gradient of log-likelihood

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of log-likelihood w.r.t. μ.

For Hessian computation in Laplace approximation.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Negative Hessian (Fisher information approximation)

Return type:

ndarray

__repr__()

String representation.

Return type:

str

class aurora.distributions.families.tweedie.TweedieFamily(power=1.5, link='log', phi=1.0)

Bases: Family

Tweedie distribution for zero-inflated positive continuous data.

A compound Poisson-Gamma distribution that naturally handles data with a point mass at zero (no events) and continuous positive values (sum of event magnitudes). The variance function is Var(Y) = phi * mu^p, where the power parameter p controls the shape.

Special cases of the Tweedie family (for reference): - p = 0: Gaussian - p = 1: Poisson - 1 < p < 2: Compound Poisson-Gamma (Tweedie proper) - p = 2: Gamma - p = 3: Inverse Gaussian

This class implements only the compound Poisson-Gamma case 1 < p < 2; for the boundary cases use the dedicated families (GaussianFamily, PoissonFamily, GammaFamily, InverseGaussianFamily).

Parameters:

• power (float, default=1.5) – Variance power parameter p. For compound Poisson-Gamma, must be in (1, 2). Common choices: - p = 1.5: Balanced, works well for many applications - p = 1.6-1.7: More right-skewed, common for insurance claims - p closer to 1: More mass at zero - p closer to 2: Less mass at zero, more Gamma-like

• link (str, default='log') – Link function: ‘log’ (default), ‘identity’, or ‘power’.

• phi (float, default=1.0) – Dispersion parameter (scale).

Variables:

• power (float) – Variance power parameter in (1, 2).

• phi (float) – Dispersion (scale) parameter.

• default_link (LinkFunction) – The link function used to map the mean to the linear predictor.

See also

aurora.distributions.families.CompoundPoissonGammaFamily

Explicit compound Poisson-Gamma alias for TweedieFamily

aurora.distributions.families.GammaFamily

Tweedie with power=2

aurora.distributions.families.InverseGaussianFamily

Tweedie with power=3

aurora.models.glm.fit_glm

Fit a GLM with this family

Notes

The Tweedie distribution can be viewed as: - N ~ Poisson(lambda) events occur - Each event has magnitude X_i ~ Gamma(alpha, beta) - Total Y = X_1 + … + X_N

This naturally creates exact zeros (N=0) and continuous positives (N>0).

The probability of zero is: P(Y=0) = exp(-lambda)

Examples

>>> from aurora.models.glm import fit_glm
>>> from aurora.distributions.families import TweedieFamily
>>> # Insurance claims (many zeros, heavy right tail)
>>> result = fit_glm(X, claims, family=TweedieFamily(power=1.7))

>>> # Rainfall amounts (zeros for dry days)
>>> result = fit_glm(X, rainfall, family=TweedieFamily(power=1.5))

References

[1]

Jorgensen (1987). Exponential dispersion models.

[2]

Smyth & Jorgensen (2002). Fitting Tweedie’s compound Poisson model.

name = 'tweedie'

valid_y_range = (0, inf)

__init__(power=1.5, link='log', phi=1.0)

Initialize Tweedie family.

Parameters:

• power (float) – Variance power parameter. Must be in (1, 2) for proper Tweedie.

• link (str) – Link function name.

• phi (float) – Dispersion parameter (scale).

property default_link

Default link function (log).

variance(mu, **params)

Variance function: V(μ) = μ^p.

Parameters:

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘phi’ (dispersion)

Returns:

Variance at each observation (without phi factor)

Return type:

ndarray

initialize(y)

Initialize mean using positive values.

Parameters:

y (ndarray) – Response values

Returns:

Initial mean estimates

Return type:

ndarray

deviance(y, mu, **params)

Deviance for Tweedie distribution.

The unit deviance for Tweedie with power p is:

For p ∈ (1, 2):

d(y, μ) = 2 × [y^(2-p)/((1-p)(2-p)) - y×μ^(1-p)/(1-p) + μ^(2-p)/(2-p)]

For y = 0:

d(0, μ) = 2 × μ^(2-p) / (2-p)

Parameters:

• y (ndarray) – Observed values (including zeros)

• mu (ndarray) – Fitted mean values

• **params (dict) – Additional parameters

Returns:

Total deviance

Return type:

float

log_likelihood(y, mu, **params)

Approximate (quasi-)log-likelihood for Tweedie.

Warning

The Tweedie density has no closed form. This method returns the quasi-likelihood -D/(2φ), which drops the normalization constant c(y, φ, p) of the full density (the series evaluation of Dunn & Smyth 2005 is not implemented). As a consequence, the returned value — and any AIC/BIC derived from it — is not comparable across families, nor across different values of the power parameter p or the dispersion φ. It is only meaningful for comparing models that share the same family, power and dispersion (e.g. nested mean models).

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

• **params (dict) – May contain ‘phi’

Returns:

Approximate (quasi-)log-likelihood

Return type:

float

weights(y, mu, **params)

IRLS weights for Tweedie.

Weights are 1/V(μ) = 1/μ^p.

Parameters:

• y (ndarray) – Observed values (not used directly)

• mu (ndarray) – Current mean estimates

Returns:

IRLS weights

Return type:

ndarray

estimate_power(y, mu, power_range=(1.1, 1.9), n_grid=20)

Estimate the power parameter via a deviance grid search.

Note

This is a heuristic: it minimizes the unit deviance over a grid of powers. It is not a profile likelihood in p — the full Tweedie density (Dunn & Smyth 2005 series evaluation) is not implemented, and the quasi-likelihood ignores the normalization constant that depends on p, so the criterion is not a valid profile likelihood for the power parameter. Treat the result as a rough starting value, not an MLE.

The family’s power attribute is always restored to its original value before returning.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted means (from current model)

• power_range (tuple) – Range of powers to search

• n_grid (int) – Number of grid points

Returns:

Estimated power parameter

Return type:

float

estimate_phi(y, mu, ddof=1)

Estimate dispersion parameter.

Uses Pearson estimator: φ = (1/n) Σ (y-μ)² / V(μ)

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Fitted means

• ddof (int) – Degrees of freedom adjustment

Returns:

Estimated dispersion

Return type:

float

d_log_likelihood(y, mu, **params)

First derivative of quasi-log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Gradient

Return type:

ndarray

d2_log_likelihood(y, mu, **params)

Second derivative of quasi-log-likelihood w.r.t. μ.

Parameters:

• y (ndarray) – Observed values

• mu (ndarray) – Mean values

Returns:

Negative Hessian diagonal (Fisher information)

Return type:

ndarray

probability_zero(mu, **params)

Probability of observing exactly zero.

For Tweedie with 1 < p < 2:

P(Y=0) = exp(-λ)

where λ = μ^(2-p) / [φ(2-p)]

Parameters:

mu (ndarray) – Mean values

Returns:

Probability of zero at each observation

Return type:

ndarray

__repr__()

String representation.

Return type:

str

Base Class

class aurora.distributions.base.Family

Bases: 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(...)

abstractmethod log_likelihood(y, mu, **params)

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:

Sum of log-likelihoods across all observations.

Return type:

scalar

abstractmethod deviance(y, mu, **params)

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:

Total deviance (non-negative, zero when mu = y).

Return type:

scalar

abstractmethod variance(mu, **params)

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:

Variance function values (same shape as mu).

Return type:

array

abstractmethod initialize(y)

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:

Initial mean estimates (same shape as y).

Return type:

array

abstract property default_link: LinkFunction

Return the default (canonical) link function for this family.

Returns:

The canonical link function instance.

Return type:

LinkFunction