"""Student-t observation noise (ADR-0058) -- the robust-regression likelihood."""
import numpy as np
from scipy.special import digamma, loggamma
from ..printing import PybnfError
from .base import NoiseModel
from .location import MEDIAN
from .scale import LINEAR
[docs]
class StudentT(NoiseModel):
"""Student-t (heavy-tailed) observation noise -- the outlier-robust likelihood, a
Gaussian with a tail-heaviness knob. It is the **first two-parameter** noise family
(ADR-0058): a location-scale family (like Gaussian/Laplace) carrying both a scale
``sigma`` and a shape ``df`` (degrees of freedom, nu). Small ``df`` gives fat tails
that downweight outliers (robust regression); ``df -> inf`` recovers the Gaussian.
This is Stan's/PyMC's ``student_t(nu, mu, sigma)``, with ``mu`` pinned to the
prediction as for every family.
Both noise parameters are **independently sourced** by the objective -- each may be
``fix_at`` a constant or ``fit`` a free parameter -- so a fit estimates 0, 1, or 2
noise parameters. ``df`` is the one parameter with a default (``DEFAULT_DF``, a
fixed 4): omit the ``df`` field and ``noise_param_defaults`` fills it, giving the
common fixed-nu robust recipe. Estimating ``df`` is statistically weakly identified
(the likelihood is nearly flat in large nu), so pair ``df = fit nu__FREE`` with a
positive prior (gamma / half_*, ADR-0057) -- not enforced here.
With ``z = (mu - forward(obs)) / sigma`` the per-point NLL splits (oracle:
``scipy.stats.t(df=nu, loc=mu, scale=sigma).logpdf``):
- ``data_fit`` (always summed): ``(nu+1)/2 * log(1 + z**2/nu)`` -- the
parameter-dependent core (depends on sigma via z, and on nu).
- the ``sigma`` normalizer ``log sigma`` -- summed iff sigma is estimated.
- the ``df`` normalizer ``-logGamma((nu+1)/2) + logGamma(nu/2) + 0.5*log(nu*pi)`` --
summed iff df is estimated. When df is fixed this whole block is a constant the
sampler drops; when df is free it is the term that keeps the fit honest. Either
way ``log_density`` (LOO/WAIC) includes it, so student_t needs no
``_density_constant`` -- the "constant when fixed" the Gaussian carries as
``0.5*log(2*pi)`` is, for student_t, this estimated-gated normalizer (ADR-0058).
Configured by the same two axes as Gaussian -- the scale its noise is additive on
and the location interpretation -- but exposed on the **linear** scale only (no
``log_student_t`` token). On the linear scale t is symmetric, so mean = median = mu
trivially. On a log scale ``base**StudentT`` has **no finite mean** (its tails are
too heavy for the MGF to exist, for any nu), so ``location = mean`` on a log scale
raises (only median centering is safe there) -- the Laplace log-scale-mean guard
(#419) taken to its limit.
"""
DEFAULT_DF = 4.0
noise_params = ('sigma', 'df')
noise_param_defaults = {'df': DEFAULT_DF}
def __init__(self, additive_on=LINEAR, location=MEDIAN):
self.additive_on = additive_on
self.location = location
[docs]
def with_location(self, location):
return type(self)(additive_on=self.additive_on, location=location)
[docs]
def mean_offset(self, noise):
"""0 on the linear scale (t is symmetric: mean = median). On **any** log scale
the original-space mean does not exist (the t-distribution's tails are heavier
than any exponential, so ``E[base**T]`` diverges for every nu), so mean-centering
is undefined -- raise, directing the user to ``location = median`` (ADR-0058)."""
t = self.additive_on.ln_base
if t == 0.0:
return 0.0
raise PybnfError(
"log-Student-t has no finite mean (its tails are too heavy for the mean to "
"exist on a log scale, for any df), so mean-centering is undefined. Use "
"location = median (the only safe centering for a Student-t on a log scale).")
def _mu(self, prediction, sigma):
"""The additive-space location parameter for ``prediction``."""
return self.additive_on.forward(prediction) - self.location.offset(self, sigma)
[docs]
def data_fit(self, prediction, observation, noise, extra=None):
sigma, nu = noise, extra['df']
z = (self._mu(prediction, sigma) - self.additive_on.forward(observation)) / sigma
return (nu + 1.) / 2. * np.log1p(z * z / nu)
[docs]
def d_data_fit_d_prediction(self, prediction, observation, noise, extra=None):
"""``d(data_fit)/d(prediction) = (nu+1) z/(nu + z**2) * forward'(pred)/sigma`` (#454),
with ``z = (mu - forward(obs))/sigma``. The factor ``(nu+1)/(nu + z**2)`` is the IRLS
weight ``w(z)`` that downweights an outlier (large ``z``); ``w(z)*z`` is the
robustified residual. PyBNF emits **only** this scalar data-fit gradient for Student-t
(no least-squares residual; the IRLS pseudo-residual is a possible later trust-region
refinement). The location offset is prediction-independent, so MEAN and MEDIAN agree
on ``d/d pred``; ``forward'(pred) = 1`` on the linear scale."""
sigma, nu = noise, extra['df']
z = (self._mu(prediction, sigma) - self.additive_on.forward(observation)) / sigma
return (nu + 1.) * z / (nu + z * z) * self.additive_on.dforward(prediction) / sigma
[docs]
def residual(self, prediction, observation, noise, extra=None):
"""The exact **square-root-loss residual** ``r = sign(z) * sqrt(2 * data_fit) =
sign(z) * sqrt((nu+1) * log1p(z**2/nu))`` with ``z = (mu - forward(obs))/sigma`` (#459) --
the least-squares residual #386's LM/TRF solver minimizes.
Unlike the IRLS pseudo-residual ``sqrt(w(z)) z``, this satisfies **both** invariants:
``1/2 r**2 == data_fit`` (so ``scipy.least_squares`` minimizes the *true* Student-t loss,
not a frozen-weight reweighted surrogate) **and** ``r * d_residual_d_prediction ==
d_data_fit_d_prediction`` (so its residual-Jacobian reproduces the objective gradient).
It is **smooth through z=0**: ``r ~ sqrt((nu+1)/nu) z`` near the origin (an odd, C-infinity
function of ``z``), behaving like a Gaussian residual at the center and downweighting the
tails as ``z`` grows -- which is why a fixed-scale Student-t fit is ``least_squares_exact``,
the Gaussian's exact-least-squares status recovered for the robust family (#459). (Laplace
has no such clean form -- ``sqrt(2*data_fit) ~ sqrt|z|`` is a cusp -- so it stays
scalar-only.) ``sign(0) == 0`` makes the residual exactly 0 at ``z=0``."""
sigma, nu = noise, extra['df']
z = (self._mu(prediction, sigma) - self.additive_on.forward(observation)) / sigma
return np.sign(z) * np.sqrt((nu + 1.) * np.log1p(z * z / nu))
[docs]
def d_residual_d_prediction(self, prediction, observation, noise, extra=None):
"""``d(r)/d(prediction) = (d r/d z) * forward'(pred)/sigma`` (#459), the residual-Jacobian
seam paired with :meth:`residual`. With ``u = z**2/nu``::
d r/d z = sqrt((nu+1)/nu) * sqrt(u / log1p(u)) / (1 + u)
the closed form of ``d_data_fit_d_prediction / r`` that stays finite at the residual zero:
as ``z -> 0`` the factor ``sqrt(u/log1p(u)) -> 1``, so ``d r/d z -> sqrt((nu+1)/nu)`` (the
slope of the linear ``r ~ sqrt((nu+1)/nu) z`` core) -- the smooth ``z->0`` limit, where the
naive ``(d data_fit/d pred)/r`` is ``0/0``. The tiny-``u`` cusp of ``u/log1p(u)`` is
series-guarded (``1 + u/2``). ``forward'(pred) = 1`` on the linear scale, the scale's chain
factor on a log scale; the offset is prediction-independent, so MEAN and MEDIAN agree."""
sigma, nu = noise, extra['df']
z = (self._mu(prediction, sigma) - self.additive_on.forward(observation)) / sigma
u = z * z / nu
# sqrt(u / log1p(u)): the 0/0 limit at z=0 is 1 (log1p(u) ~ u); series-guard the tiny-u cusp.
ratio = u / np.log1p(u) if u > 1e-8 else 1. + u / 2.
d_r_d_z = np.sqrt((nu + 1.) / nu) * np.sqrt(ratio) / (1. + u)
return d_r_d_z * self.additive_on.dforward(prediction) / sigma
[docs]
def d_nll_d_noise_params(self, prediction, observation, noise, extra=None):
"""The estimated-scale gradient columns ``{'sigma': ..., 'df': ...}`` (#451/#454/#385) --
the first **multi-parameter** estimated-noise gradient (ADR-0058). Each is the derivative
of ``data_fit`` plus that parameter's own normalizer (``log sigma`` for sigma, the df-block
for df); the normalizers depend only on their own parameter, so the cross terms vanish.
With ``z = (mu - forward(obs))/sigma``::
d loss/d sigma = nu (1 - z**2) / (sigma (nu + z**2))
d loss/d nu = 1/2 log1p(z**2/nu) - (nu+1) z**2 / (2 nu (nu + z**2))
+ 1/2 (digamma(nu/2) - digamma((nu+1)/2) + 1/nu)
The sigma column folds ``d(log sigma)/d sigma = 1/sigma`` into ``d(data_fit)/d sigma``
and reduces to Gaussian's ``(1 - z**2)/sigma`` as ``nu -> inf``; the df column folds the
data fit's nu-dependence with the df-block's digamma derivative -- the term that keeps a
free df honest. The location offset is noise-independent wherever Student-t's mean centering
is defined -- the linear scale, where the offset is 0; on a log scale Student-t has no finite
mean, so mean-centering raises in ``mean_offset`` regardless (no ``d_mean_offset_d_noise``
coupling term is needed, unlike Gaussian/Laplace on a log scale, #385)."""
sigma, nu = noise, extra['df']
z = (self._mu(prediction, sigma) - self.additive_on.forward(observation)) / sigma
d_sigma = nu * (1. - z * z) / (sigma * (nu + z * z))
d_data_fit_d_nu = 0.5 * np.log1p(z * z / nu) - (nu + 1.) * z * z / (2. * nu * (nu + z * z))
d_block_d_nu = 0.5 * (digamma(nu / 2.) - digamma((nu + 1.) / 2.) + 1. / nu)
return {'sigma': d_sigma, 'df': d_data_fit_d_nu + d_block_d_nu}
[docs]
def param_normalizers(self, noise, extra=None):
sigma, nu = noise, extra['df']
df_block = -loggamma((nu + 1.) / 2.) + loggamma(nu / 2.) + 0.5 * np.log(nu * np.pi)
return {'sigma': np.log(sigma), 'df': df_block}
