"""Negative-binomial observation noise (ADR-0011, ADR-0031)."""
import numpy as np
from scipy.optimize import brentq
from scipy.special import betainc, betaln, digamma, loggamma
from ..printing import PybnfError
from .base import NoiseModel
from .location import MEDIAN
#: The clip :meth:`NegBinomial.data_fit` puts on ``prob = r/(r+mean)`` to keep a
#: degenerate point's ``-logpmf`` finite.
_PROB_CLIP = 1e-10
#: The mean-side image of that clip: ``prob <= 1 - _PROB_CLIP`` iff
#: ``mean >= r * _PROB_CLIP/(1 - _PROB_CLIP)``. Capped at ``_MEAN_FLOOR_MAX``, it is the floor
#: (:func:`_mean_floor`) the score and every derivative that divides by the mean put on it, so a
#: zero prediction yields a finite score and slope rather than ``nan``/``-inf``.
_MEAN_FLOOR_REL = _PROB_CLIP / (1 - _PROB_CLIP)
#: The cap on that floor. ``r * _MEAN_FLOOR_REL`` grows with the dispersion, and from r ~ 1e7 on
#: it reached the counts being scored, so every expected count below it was scored as the floor
#: (#987). Capped at 1e-7, the floor is unchanged for r up to ~1e3 and keeps a zero prediction
#: finite at any r.
_MEAN_FLOOR_MAX = 1e-7
#: Where ``d(data_fit)/dr`` switches to its large-r series: ``r >= _SERIES_R_FACTOR * max(k, m, 1)``.
_SERIES_R_FACTOR = 1e4
def _d_fit_d_r_series(r, k, m):
"""``d(data_fit)/dr`` at a fixed mean for ``r`` far above ``k`` and ``m``. The closed form,
``psi(r) - psi(r + k) + log1p(m/r) + (k - m)/(r + m)``, has terms of order ``max(k, m)/r``
that cancel to order ``1/r**2``, and the digamma difference is held to absolute precision
only, so at a large ``r`` it is noise: at ``r = 1e12`` it summed to -0.096 for +4.5e-10, the
wrong sign (#987 review). The expansion in ``1/r`` instead: ``((k - m)**2 - k) / (2 r**2) +
(k**2/2 - k/6 - k**3/3 + k m**2 - 2 m**3/3) / r**3 + O(r**-4)``, relative error of order
``(max(k, m)/r)**2`` past the switch."""
x = 1.0 / r
second = ((k - m) ** 2 - k) / 2.0
third = k * k / 2.0 - k / 6.0 - k ** 3 / 3.0 + k * m * m - 2.0 * m ** 3 / 3.0
return x * x * (second + x * third)
[docs]
def draw_count(rng, mean, r):
"""One negative binomial count of ``mean`` and dispersion ``r`` from the numpy Generator
``rng``. Up to r ~ 1e3 it is the draw PyBNF has always made, ``nbinom.rvs`` at the clipped
``prob = r/(r + mean)``, so a seeded draw is unchanged to the bit. A larger ``r`` took the
clip's floor of ``r * 1e-10`` as the mean too (#987: at r = 1e12 every mean below 100 was
drawn as ~100), and ``prob`` rounds next to 1 there, so it draws gamma-Poisson with the
floor capped, in the same order of the same stream as numpy's own negative binomial."""
from scipy import stats
if r * _MEAN_FLOOR_REL <= _MEAN_FLOOR_MAX:
prob = np.clip(r / (r + mean), _PROB_CLIP, 1 - _PROB_CLIP)
return stats.nbinom.rvs(n=r, p=prob, size=1, random_state=rng)[0]
return rng.poisson(rng.gamma(r, max(mean, _mean_floor(r)) / r))
def _mean_floor(r):
"""The smallest mean the score and its derivatives take for dispersion ``r``."""
return min(r * _MEAN_FLOOR_REL, _MEAN_FLOOR_MAX)
def _d_betainc_d_b(a, b, x, rel=1e-6):
"""``d/db`` of the regularized incomplete beta ``betainc(a, b, x) = I_x(a, b)`` w.r.t. its
**second** parameter, by a central finite difference (issue #458).
The median centering puts the prediction in this second parameter (``b = target + 1``, the
continuous CDF ``I_p(r, pred + 1)``), so the implicit derivative ``d mean/d pred`` needs
``dI_x/db`` -- which is **not elementary**: it brings in the digamma function and a
non-elementary ``int_0^x t^(a-1)(1-t)^(b-1) ln(1-t) dt``. #458 takes the issue-sanctioned
*numerically-evaluated* parameter derivative -- a central difference in ``b`` with ``x`` held
fixed (``x = r/(r+mean)`` does not depend on ``pred``). ``b = target + 1 >= 1`` here, so the
relative step keeps ``b - db > 0``; at the median ``x`` is bounded away from 0 and 1, so
``betainc`` is smooth and the difference is accurate to ~1e-10 (validated against the integral
representation ``int .../B(a,b) - I_x(a,b)(psi(b) - psi(a+b))``)."""
db = rel * max(abs(b), 1.0)
return (betainc(a, b + db, x) - betainc(a, b - db, x)) / (2.0 * db)
def _d_betainc_d_a(a, b, x, rel=1e-6):
"""``d/da`` of the regularized incomplete beta ``betainc(a, b, x) = I_x(a, b)`` w.r.t. its
**first** parameter, by a central finite difference -- the sibling of :func:`_d_betainc_d_b`
(issue #458). A free dispersion under MEDIAN centering makes the median's mean depend on ``r``
through the first beta parameter ``a = r`` (the CDF is ``I_p(r, target+1)``), so the implicit
``d mean/d r`` needs ``dI_x/da`` -- as non-elementary as the second-parameter derivative.
``a = r > 0``, so the relative step keeps ``a - da > 0``; at the median ``x`` is bounded away
from 0 and 1, so the difference is accurate to ~1e-10 (validated against the integral
representation, like :func:`_d_betainc_d_b`)."""
da = rel * max(abs(a), 1.0)
return (betainc(a + da, b, x) - betainc(a - da, b, x)) / (2.0 * da)
def _mean_for_median(prediction, r):
"""Solve for the mean ``mu`` of ``NB(mean=mu, dispersion=r)`` whose **continuous**
0.5-quantile equals ``prediction`` -- the negative-binomial median realization
(issue #419, ADR-0031's "every means every").
The continuous CDF is ``F(x; mu, r) = I_p(r, x + 1)`` with ``p = r / (r + mu)``
(``scipy.special.betainc``, the regularized incomplete beta), which is exactly
``scipy.stats.nbinom.cdf(k, r, p)`` at integer ``k`` but smooth in its second
argument -- so we use it instead of the discrete ``nbinom.ppf`` step to keep the
objective continuous in the prediction (the optimizers need it). ``F`` is strictly
decreasing in ``mu`` (larger mean shifts the distribution right, lowering the mass
at or below the prediction), so there is a unique ``mu`` placing the median at the
prediction, found by a bounded root-find.
A prediction is a count median, so it is clamped to ``>= 0``; ``mu = 0`` gives
``p = 1`` and ``F = 1 > 0.5``, while ``mu -> inf`` gives ``F -> 0``, bracketing the
root in ``[0, hi]``.
"""
target = max(prediction, 0.0)
def gap(mu):
p = r / (r + mu)
return betainc(r, target + 1.0, p) - 0.5
# gap(0) == 0.5 > 0; grow the upper bound until the median exceeds the target.
hi = max(target, 1.0)
while gap(hi) > 0.0:
hi *= 2.0
if hi > 1e300:
# Below r ~ 1e-3 a median of 0 needs a mean beyond any double. Returning the bound
# (1e15, as this did) scored that mean instead, with no word (#987 review).
raise PybnfError(
'The negative binomial median cannot be placed at %r with dispersion %r'
% (prediction, r),
'At dispersion %r the mean that puts the median at %r is beyond 1e300. Use '
'location = mean, or a larger dispersion.' % (r, prediction))
return brentq(gap, 0.0, hi)
[docs]
class NegBinomial(NoiseModel):
"""Negative-binomial observation noise for count data. The dispersion ``r`` is the
noise parameter (``neg_bin`` reads it from the config constant ``neg_bin_r``;
``neg_bin_dynamic`` from the ``r__FREE`` free parameter).
The **location** axis (ADR-0011/0031) sets which distributional summary the
prediction is taken to be. The default is ``MEDIAN`` -- median is the universal
prediction-centering default for *every* noise family (ADR-0031, "every means
every"), true in code at the constructor like Gaussian/Laplace. ``MEDIAN``
interprets the prediction as the 0.5-quantile and solves for the mean placing the
continuous median there (issue #419). ``MEAN`` is the native parameterization (the
prediction *is* the mean) -- the legacy ``neg_bin`` objfuncs pin it explicitly to
stay frozen-mean. Unlike Gaussian/Laplace, the count family is **not additive on a
scale**, so it owns this realization directly rather than going through
``location.py``'s additive-offset abstraction -- it reuses the ``MEAN``/``MEDIAN``
markers, not the ``offset`` math.
A negative observed count contributes nothing to the *cost* (the count-domain guard on
``data_fit`` and every derivative) and is not a scored point at all for the *density*
(:meth:`observation_in_domain`, #523). A PMF is self-normalizing, so there is no
separable normalizer (``log_normalizer`` stays 0) and the full ``-logpmf`` lives in
``data_fit``.
"""
noise_params = ('dispersion',)
#: The count domain, for the objective's exclusion warning (#523).
observation_domain = 'a non-negative count'
#: A count family: a normalized column is refused for it when the configuration loads (#997).
scores_counts = True
[docs]
def observation_in_domain(self, observation):
"""A count is non-negative, so a **negative** observation is outside the negative
binomial's domain -- it has no probability under the family (#523).
Real surveillance data contains such values: a jurisdiction revising its cumulative
total downward produces a negative daily increment. The cost-path guards
(:meth:`data_fit` and the derivatives) score such a point as 0, which is right for
the fit -- it contributes nothing. It is emphatically *not* right for the density,
where a zero ``data_fit`` on a self-normalizing PMF reads as ``log p = 0``, a
probability of 1: better than the family assigns any real count, including one the
model predicts exactly. Saying the observation is out of domain routes it to the
same treatment as a NaN -- excluded from ``n`` for AIC/BIC and from the LOO/WAIC
observation axis -- instead of counted as a perfectly-fit point."""
return observation >= 0
def __init__(self, location=MEDIAN):
self.location = location
[docs]
def with_location(self, location):
return NegBinomial(location=location)
def _mean(self, prediction, noise):
"""The distribution mean for ``prediction`` under the location interpretation:
the prediction itself for ``MEAN``, the median inversion for ``MEDIAN``. A negative
prediction is the mean 0, since unclamped ``prob = r/(r + mean)`` leaves ``(0, 1]``."""
if self.location is MEDIAN:
return _mean_for_median(prediction, noise)
return max(prediction, 0.0)
def _mean_for_slope(self, prediction, noise):
""":meth:`_mean`, floored away from zero for the expressions that divide by it.
A MEAN-centered prediction of **exactly** zero is not pathological -- it is what
any model whose output is gated off over part of the fit window predicts there (an
epidemic model before its start time ``t0``, a stimulus-driven readout before the
stimulus). :meth:`data_fit` scores such a point finitely because it floors the mean
(:func:`_mean_floor`); every *derivative* below instead
divides by the mean (``r (mean - obs) / (mean (r + mean))``, ``r / (mean (r+mean))``),
which at ``mean == 0`` is ``0/0 -> nan`` for an observed zero and ``-inf`` otherwise.
That ``nan`` propagates into the gradient and out to the optimizer, which then tries
to assign a free parameter ``nan``.
So take the value path's floor here too: ``r _PROB_CLIP / (1 - _PROB_CLIP)``, capped at
``_MEAN_FLOOR_MAX`` (#987). The slope is then the finite value belonging to the objective
actually being scored -- large (order ``obs`` over the floor) where the model predicts
zero against a positive
count, which is correct: the objective there really is nearly flat-then-cliff, and a
large finite push toward a positive prediction is the right search direction.
MEDIAN centering never reaches the floor (``_mean_for_median`` returns a strictly
positive mean for any prediction), so it is a no-op there."""
return max(self._mean(prediction, noise), _mean_floor(noise))
[docs]
def data_fit(self, prediction, observation, noise, extra=None):
if observation < 0:
return 0
mean = self._mean(prediction, noise)
assert isinstance(noise, float)
prob = noise / (noise + mean)
# log of the negative-binomial PMF P(observation | r=noise, prob)
# == scipy.stats.nbinom.logpmf(observation, noise, prob).
if noise * _MEAN_FLOOR_REL <= _MEAN_FLOOR_MAX and prob >= _PROB_CLIP:
# r up to ~1e3, and a mean below r * 1e10: the clip binds only for a mean on the
# floor (a zero prediction), and the score is the one it has always been, to the bit.
prob = min(prob, 1 - _PROB_CLIP)
log_pmf = loggamma(observation + noise) - loggamma(observation + 1) - loggamma(noise) \
+ noise * np.log(prob) + observation * np.log(1 - prob)
else:
# A larger r, or a mean beyond r * 1e10 (#987): the clip would score the mean as
# another, and the form above loses digits to cancellation (2e-7 relative at
# r = 1e8, 0.1% at 1e12). The PMF of the mean itself, kept above the capped floor:
# log p and log(1 - p) without forming p, and the binomial coefficient as
# -log k - betaln(r, k) rather than loggamma(k + r) - loggamma(r).
m = max(mean, _mean_floor(noise))
log_r, log_m = np.log(noise), np.log(m)
log_sum = np.logaddexp(log_r, log_m)
# m / r overflows once m > r * 1.8e308, so a mean above r takes the difference of logs.
log_p = -np.log1p(m / noise) if m <= noise else log_r - log_sum
log_q = log_m - log_sum
log_coefficient = (-np.log(observation) - betaln(noise, observation)
if observation > 0 else 0.0)
log_pmf = log_coefficient + noise * log_p + observation * log_q
# A PMF is <= 1, so log_pmf <= 0; PyBNF minimizes the negative
# log-likelihood -log_pmf >= 0.
return -log_pmf
[docs]
def d_data_fit_d_prediction(self, prediction, observation, noise, extra=None):
"""``d(data_fit)/d(prediction)`` for one count point (#458, the deferred layer-G follow-up
of #385). The count family carries no least-squares residual (its data fit is a ``-logpmf``,
not a sum of squares), so PyBNF emits **only** this scalar data-fit gradient.
The chain ``d(data_fit)/d(mean) * d(mean)/d(prediction)``. The mean-slope is the closed-form
negative-binomial score ``r (mean - obs) / (mean (r + mean))`` (``r = noise`` the
dispersion); the mean-factor depends on the location interpretation:
* **MEAN**: the prediction *is* the mean, so ``d mean/d pred = 1`` and the slope is that
closed form directly -- the clean case (``neg_bin`` / ``neg_bin_dynamic``).
* **MEDIAN**: the mean is the CDF inversion ``_mean_for_median`` placing the continuous
median at the prediction, so ``d mean/d pred`` is the **implicit derivative** of that
root-find, ``-(dG/d pred) / (dG/d mean)`` with ``G(mean, pred) = betainc(r, target+1, p)
- 0.5`` and ``p = r/(r+mean)``. ``dG/d mean = -beta_pdf(p; r, target+1) * r/(r+mean)**2``
is the smooth beta-density chain rule (``F`` strictly decreasing in the mean, so this is
negative); ``dG/d pred`` puts the prediction in the beta's *second* parameter, so it is
the non-elementary ``d betainc/d b`` (:func:`_d_betainc_d_b`), times ``d target/d pred =
1``. The result is positive (a larger predicted median needs a larger mean).
A negative observation contributes nothing (the count-domain guard, mirroring
:meth:`data_fit`). A prediction clamped to the count floor (``pred <= 0``, where ``target``
floors at 0 and the median stops moving) has slope 0 -- a kink at ``pred == 0`` where PyBNF
takes the floor-side subgradient, like the Laplace kink (#454). The mean-slope divides by
the mean, so it reads it through :meth:`_mean_for_slope` -- the value path's ``prob`` clip
expressed on the mean -- which keeps a MEAN-centered zero prediction finite instead of
``nan``; away from that floor the form is the un-clipped analytic one."""
if observation < 0:
return 0.0
r = noise
mean = self._mean_for_slope(prediction, r)
d_fit_d_mean = r * (mean - observation) / (mean * (r + mean))
if self.location is not MEDIAN:
return d_fit_d_mean # MEAN: d mean/d pred = 1
if prediction <= 0.0:
return 0.0 # clamped to the count floor: the median stops moving
target = prediction
p = r / (r + mean)
# beta_pdf(p; r, target+1) = p**(r-1) (1-p)**target / B(r, target+1) -- via betaln for range.
beta_pdf = np.exp((r - 1.0) * np.log(p) + target * np.log1p(-p) - betaln(r, target + 1.0))
dG_d_mean = -beta_pdf * r / (r + mean) ** 2.
dG_d_pred = _d_betainc_d_b(r, target + 1.0, p)
d_mean_d_pred = -dG_d_pred / dG_d_mean
return d_fit_d_mean * d_mean_d_pred
[docs]
def d_nll_d_noise_params(self, prediction, observation, noise, extra=None):
"""``{'dispersion': d(data_fit)/d r}`` -- the estimated-dispersion gradient column for a free
``r`` (``neg_bin_dynamic``'s ``r__FREE``; #458, generalizing layer D/G of #385).
The negative-binomial PMF is **self-normalizing** (``log_normalizer == 0``), so -- unlike a
Gaussian's ``+log sigma`` or a Laplace's ``log(2 b)`` -- there is no separable normalizer:
the whole dispersion gradient lives in the data fit. With ``mean`` the distribution mean and
``prob = r/(r+mean)`` the partial holding the mean fixed is the negative-binomial dispersion
score::
d(data_fit)/d r |_mean = psi(r) - psi(obs + r) - log(prob) - 1 + (r + obs)/(r + mean)
(``psi`` the digamma; the ``-logpmf``'s ``r``-dependence through its gamma terms and ``prob``).
* **MEAN**: the prediction *is* the mean, r-independent, so this partial is the whole
derivative.
* **MEDIAN**: the mean is solved from ``r`` (the CDF inversion ``_mean_for_median``), so add
the coupling ``d(data_fit)/d mean * d mean/d r``. ``d mean/d r`` is the implicit derivative
of ``G(mean, r) = betainc(r, target+1, p) - 0.5 = 0`` (``p = r/(r+mean)``) holding the
prediction fixed, ``-(dG/d r)/(dG/d mean)``; ``dG/d r`` brings in the betainc **first**-
parameter derivative (``a = r``, :func:`_d_betainc_d_a`) plus ``p``'s own ``r``-dependence,
and ``dG/d mean`` is the same beta-density chain rule as the prediction gradient. The count
floor (``target = max(pred, 0)``) does not zero this -- the mean still moves with ``r`` even
at a clamped prediction. (Validated FD-first: the partial-only form is off 5-15%, sometimes
far more, on the median.)
A negative observation contributes nothing (the count-domain guard). The mean is read
through :meth:`_mean_for_slope` for consistency with :meth:`d_data_fit_d_prediction` --
it matters only for the MEDIAN coupling term, which divides by the mean; the MEAN
partial above is finite at a zero prediction either way, and the floor shifts it by
``O(_PROB_CLIP)``. MEDIAN centering never reaches the floor."""
if observation < 0:
return {'dispersion': 0.0}
r = noise
mean = self._mean_for_slope(prediction, r)
prob = r / (r + mean)
if r >= _SERIES_R_FACTOR * max(observation, mean, 1.0):
d_fit_d_r = _d_fit_d_r_series(r, observation, mean)
else:
d_fit_d_r = (digamma(r) - digamma(observation + r) - np.log(prob) - 1.0
+ (r + observation) / (r + mean))
if self.location is MEDIAN:
target = max(prediction, 0.0)
d_fit_d_mean = r * (mean - observation) / (mean * (r + mean))
beta_pdf = np.exp((r - 1.0) * np.log(prob) + target * np.log1p(-prob)
- betaln(r, target + 1.0))
dG_d_mean = -beta_pdf * r / (r + mean) ** 2.
dG_d_r = _d_betainc_d_a(r, target + 1.0, prob) + beta_pdf * mean / (r + mean) ** 2.
d_fit_d_r += d_fit_d_mean * (-dG_d_r / dG_d_mean)
return {'dispersion': d_fit_d_r}
[docs]
def location_fisher(self, prediction, observation, noise, extra=None):
"""``kappa = I_mean = r / (mean (r + mean))`` -- the negative-binomial **mean's** Fisher
information, the Gauss-Newton curvature the EFIM trust-region path (``job_type = gntr``,
#481) weights ``d(prediction)/d(theta)`` by. It is the variance of the mean score
``d(-logpmf)/d mean = r (mean - obs) / (mean (r + mean))`` (``Var(obs) = mean (r + mean)/r``,
so ``I_mean = [r/(mean(r+mean))]**2 * Var(obs) = r/(mean(r+mean))``).
* **MEAN**: the prediction *is* the mean (``d mean/d pred = 1``), so ``kappa`` is this
directly -- the case this cut supports (``neg_bin`` / ``neg_bin_dynamic`` pinned to MEAN).
* **MEDIAN**: the mean sits behind the betainc CDF inversion (``_mean_for_median``), whose
``d mean/d pred`` is the non-elementary implicit derivative; its location Fisher is out of
scope for this cut -- refused, pointing at ``job_type = lbfgs`` (which fits it via the
scalar data-fit gradient). A negative observation contributes no curvature (the
count-domain guard, mirroring :meth:`data_fit`).
A **zero mean** deliberately reports **no curvature** rather than the
:meth:`_mean_for_slope` floor the *gradient* uses. The two want opposite things: the
true ``I_mean`` diverges as ``1/mean``, so flooring it would put a ``~1/(r
_PROB_CLIP)`` weight on every gated-off row (an epidemic model's whole pre-``t0``
stretch) and let those rows dominate the Fisher matrix, crushing the trust-region
step. Zero says "this row carries no information about the mean", which is the stable
reading and leaves the gradient -- which does floor -- to supply the descent
direction."""
if self.location is MEDIAN:
from ..gradient.errors import GradientNotSupported
raise GradientNotSupported(
"A MEDIAN-centered negative-binomial has its mean behind a betainc CDF inversion, "
"whose location Fisher the EFIM trust-region path (job_type = gntr, #481) does not "
"assemble this cut; use job_type = lbfgs.")
mean = self._mean(prediction, noise)
if observation < 0 or mean <= 0.0:
return 0.0
return noise / (mean * (noise + mean))
