Skip to content

poisson

poisson

Poisson distribution for ngboost-lightning.

Poisson 

Poisson(params: NDArray[floating])

Bases: Distribution

Poisson distribution with log-rate parameterization.

Internal parameter is [log_rate] where rate = exp(log_rate). The log-link ensures rate stays positive during unconstrained boosting.

This is a discrete distribution: score uses logPMF, not logPDF.

ATTRIBUTE DESCRIPTION
n_params

Always 1 (log_rate).
rate

Poisson rate (lambda) values, shape [n_samples].

TYPE: NDArray[floating]

Construct Poisson from internal parameters.

PARAMETER DESCRIPTION
params

Internal parameters, shape [n_samples, 1]. Column 0 is log(rate).

TYPE: NDArray[floating]

Source code in ngboost_lightning/distributions/poisson.py 
29
30
31
32
33
34
35
36
37
38
39
def __init__(self, params: NDArray[np.floating]) -> None:
    """Construct Poisson from internal parameters.

    Args:
        params: Internal parameters, shape ``[n_samples, 1]``.
            Column 0 is log(rate).
    """
    log_rate: NDArray[np.floating] = params[:, 0]
    self.rate: NDArray[np.floating] = np.exp(log_rate)
    self._dist = sp_poisson(mu=self.rate)
    self._params = params

fit staticmethod 

fit(
    y: NDArray[floating],
    sample_weight: NDArray[floating] | None = None,
) -> NDArray[floating]

Estimate initial log_rate from target data.

PARAMETER DESCRIPTION
y

Target values (counts), shape [n_samples].

TYPE: NDArray[floating]

sample_weight

Per-sample weights, shape [n_samples].

TYPE: NDArray[floating] | None DEFAULT: None

RETURNS DESCRIPTION
NDArray[floating]

Parameter vector [log(rate)], shape [1].
Source code in ngboost_lightning/distributions/poisson.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
@staticmethod
def fit(
    y: NDArray[np.floating],
    sample_weight: NDArray[np.floating] | None = None,
) -> NDArray[np.floating]:
    """Estimate initial log_rate from target data.

    Args:
        y: Target values (counts), shape ``[n_samples]``.
        sample_weight: Per-sample weights, shape ``[n_samples]``.

    Returns:
        Parameter vector ``[log(rate)]``, shape ``[1]``.
    """
    rate = max(float(np.average(y, weights=sample_weight)), 1e-6)
    return np.array([np.log(rate)])

score 

score(y: NDArray[floating]) -> NDArray[floating]

Per-sample negative log-likelihood (using logPMF).

PARAMETER DESCRIPTION
y

Observed count values, shape [n_samples].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

NLL values, shape [n_samples].
Source code in ngboost_lightning/distributions/poisson.py 
58
59
60
61
62
63
64
65
66
67
def score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Per-sample negative log-likelihood (using logPMF).

    Args:
        y: Observed count values, shape ``[n_samples]``.

    Returns:
        NLL values, shape ``[n_samples]``.
    """
    return -self._dist.logpmf(y)

d_score 

d_score(y: NDArray[floating]) -> NDArray[floating]

Analytical gradient of NLL w.r.t. [log_rate].

Derivation

NLL = -(y * log(rate) - rate - gammaln(y + 1)) d(NLL)/d(rate) = -(y / rate - 1) = 1 - y / rate d(NLL)/d(log_rate) = d(NLL)/d(rate) * d(rate)/d(log_rate) = (1 - y / rate) * rate = rate - y

PARAMETER DESCRIPTION
y

Observed count values, shape [n_samples].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

Gradient array, shape [n_samples, 1].
Source code in ngboost_lightning/distributions/poisson.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def d_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Analytical gradient of NLL w.r.t. [log_rate].

    Derivation:
        NLL = -(y * log(rate) - rate - gammaln(y + 1))
        d(NLL)/d(rate) = -(y / rate - 1) = 1 - y / rate
        d(NLL)/d(log_rate) = d(NLL)/d(rate) * d(rate)/d(log_rate)
                           = (1 - y / rate) * rate
                           = rate - y

    Args:
        y: Observed count values, shape ``[n_samples]``.

    Returns:
        Gradient array, shape ``[n_samples, 1]``.
    """
    n = len(y)
    grad = np.empty((n, 1))
    grad[:, 0] = self.rate - y
    return grad

metric 

metric() -> NDArray[floating]

Fisher Information for Poisson with log-rate parameterization.

For Poisson(rate), Var(Y) = rate, and the Fisher information w.r.t. rate is 1/rate. Applying the chain rule for the log-rate parameterization: FI(log_rate) = rate^2 * (1/rate) = rate.

RETURNS DESCRIPTION
NDArray[floating]

FI tensor, shape [n_samples, 1, 1].
Source code in ngboost_lightning/distributions/poisson.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def metric(self) -> NDArray[np.floating]:
    """Fisher Information for Poisson with log-rate parameterization.

    For Poisson(rate), Var(Y) = rate, and the Fisher information
    w.r.t. rate is 1/rate. Applying the chain rule for the log-rate
    parameterization: FI(log_rate) = rate^2 * (1/rate) = rate.

    Returns:
        FI tensor, shape ``[n_samples, 1, 1]``.
    """
    n = len(self.rate)
    fi = np.empty((n, 1, 1))
    fi[:, 0, 0] = self.rate
    return fi

natural_gradient 

natural_gradient(y: NDArray[floating]) -> NDArray[floating]

Natural gradient via scalar Fisher (fast path).

Since n_params=1, the natural gradient is simply d_score / rate: nat_grad[:, 0] = (rate - y) / rate = 1 - y / rate

PARAMETER DESCRIPTION
y

Observed count values, shape [n_samples].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

Natural gradient, shape [n_samples, 1].
Source code in ngboost_lightning/distributions/poisson.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def natural_gradient(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Natural gradient via scalar Fisher (fast path).

    Since n_params=1, the natural gradient is simply d_score / rate:
        nat_grad[:, 0] = (rate - y) / rate = 1 - y / rate

    Args:
        y: Observed count values, shape ``[n_samples]``.

    Returns:
        Natural gradient, shape ``[n_samples, 1]``.
    """
    n = len(y)
    nat_grad = np.empty((n, 1))
    nat_grad[:, 0] = 1.0 - y / self.rate
    return nat_grad

mean 

mean() -> NDArray[floating]

Conditional mean (point prediction).

RETURNS DESCRIPTION
NDArray[floating]

Rate values, shape [n_samples].
Source code in ngboost_lightning/distributions/poisson.py 
122
123
124
125
126
127
128
def mean(self) -> NDArray[np.floating]:
    """Conditional mean (point prediction).

    Returns:
        Rate values, shape ``[n_samples]``.
    """
    return self.rate

sample 

sample(n: int) -> NDArray[floating]

Draw n samples per distribution instance.

PARAMETER DESCRIPTION
n

Number of samples to draw.

TYPE: int

RETURNS DESCRIPTION
NDArray[floating]

Samples, shape [n, n_samples].
Source code in ngboost_lightning/distributions/poisson.py 
130
131
132
133
134
135
136
137
138
139
def sample(self, n: int) -> NDArray[np.floating]:
    """Draw n samples per distribution instance.

    Args:
        n: Number of samples to draw.

    Returns:
        Samples, shape ``[n, n_samples]``.
    """
    return self._dist.rvs(size=(n, len(self)))

cdf 

cdf(y: NDArray[floating]) -> NDArray[floating]

Cumulative distribution function P(X <= y).

PARAMETER DESCRIPTION
y

Values at which to evaluate the CDF.

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

CDF values, same shape as y.
Source code in ngboost_lightning/distributions/poisson.py 
141
142
143
144
145
146
147
148
149
150
def cdf(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Cumulative distribution function P(X <= y).

    Args:
        y: Values at which to evaluate the CDF.

    Returns:
        CDF values, same shape as ``y``.
    """
    return self._dist.cdf(y)

ppf 

ppf(q: NDArray[floating]) -> NDArray[floating]

Percent point function (inverse CDF / quantile function).

For discrete distributions, returns the smallest integer k such that CDF(k) >= q.

PARAMETER DESCRIPTION
q

Quantiles, values in [0, 1].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

Integer-valued quantiles, same shape as q.
Source code in ngboost_lightning/distributions/poisson.py 
152
153
154
155
156
157
158
159
160
161
162
163
164
def ppf(self, q: NDArray[np.floating]) -> NDArray[np.floating]:
    """Percent point function (inverse CDF / quantile function).

    For discrete distributions, returns the smallest integer k such
    that CDF(k) >= q.

    Args:
        q: Quantiles, values in [0, 1].

    Returns:
        Integer-valued quantiles, same shape as ``q``.
    """
    return self._dist.ppf(q)

logpdf 

logpdf(y: NDArray[floating]) -> NDArray[floating]

Log probability mass function.

Note: For ABC compatibility this method is named logpdf, but for this discrete distribution it returns the log-PMF.

PARAMETER DESCRIPTION
y

Values at which to evaluate.

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

Log-PMF values, same shape as y.
Source code in ngboost_lightning/distributions/poisson.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
def logpdf(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Log probability mass function.

    Note: For ABC compatibility this method is named ``logpdf``, but
    for this discrete distribution it returns the log-PMF.

    Args:
        y: Values at which to evaluate.

    Returns:
        Log-PMF values, same shape as ``y``.
    """
    return self._dist.logpmf(y)

crps_score 

crps_score(y: NDArray[floating]) -> NDArray[floating]

Per-sample CRPS for Poisson.

Closed form (Czado, Gneiting & Held 2009): CRPS = (y - lam)*(2*F(y; lam) - 1) + 2*lam*f(floor(y); lam) - lam*exp(-2*lam)*(I_0(2*lam) + I_1(2*lam))

where lam = rate, F is the Poisson CDF, f is the PMF, and I_0, I_1 are modified Bessel functions of the first kind.

Uses exponentially-scaled Bessel functions (i0e, i1e) to avoid overflow for large lam: exp(-2*lam)*I_k(2*lam) = i_ke(2*lam)

PARAMETER DESCRIPTION
y

Observed count values, shape [n_samples].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

CRPS values, shape [n_samples].
Source code in ngboost_lightning/distributions/poisson.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
def crps_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Per-sample CRPS for Poisson.

    Closed form (Czado, Gneiting & Held 2009):
        ``CRPS = (y - lam)*(2*F(y; lam) - 1)
                 + 2*lam*f(floor(y); lam)
                 - lam*exp(-2*lam)*(I_0(2*lam) + I_1(2*lam))``

    where ``lam = rate``, ``F`` is the Poisson CDF, ``f`` is the PMF,
    and ``I_0``, ``I_1`` are modified Bessel functions of the first kind.

    Uses exponentially-scaled Bessel functions (``i0e``, ``i1e``) to
    avoid overflow for large ``lam``:
        ``exp(-2*lam)*I_k(2*lam) = i_ke(2*lam)``

    Args:
        y: Observed count values, shape ``[n_samples]``.

    Returns:
        CRPS values, shape ``[n_samples]``.
    """
    lam = self.rate
    y_floor = np.floor(y)
    cdf_y = self._dist.cdf(y_floor)
    pmf_y = self._dist.pmf(y_floor)
    two_lam = 2.0 * lam
    # exp(-2*lam) * I_0(2*lam) + exp(-2*lam) * I_1(2*lam)
    bessel_term = i0e(two_lam) + i1e(two_lam)
    result: NDArray[np.floating] = (
        (y - lam) * (2.0 * cdf_y - 1.0) + 2.0 * lam * pmf_y - lam * bessel_term
    )
    return result

crps_d_score 

crps_d_score(y: NDArray[floating]) -> NDArray[floating]

Gradient of CRPS w.r.t. [log_rate].

Uses central finite differences on the internal parameter. The analytical gradient involves derivatives of both the Poisson CDF and modified Bessel functions, making finite differences the more robust approach.

PARAMETER DESCRIPTION
y

Observed count values, shape [n_samples].

TYPE: NDArray[floating]

RETURNS DESCRIPTION
NDArray[floating]

Gradient array, shape [n_samples, 1].
Source code in ngboost_lightning/distributions/poisson.py 
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
def crps_d_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Gradient of CRPS w.r.t. [log_rate].

    Uses central finite differences on the internal parameter.
    The analytical gradient involves derivatives of both the Poisson
    CDF and modified Bessel functions, making finite differences
    the more robust approach.

    Args:
        y: Observed count values, shape ``[n_samples]``.

    Returns:
        Gradient array, shape ``[n_samples, 1]``.
    """
    eps = 1e-5
    params = self._params.copy()
    params_plus = params.copy()
    params_plus[:, 0] += eps
    params_minus = params.copy()
    params_minus[:, 0] -= eps
    score_plus = type(self)(params_plus).crps_score(y)
    score_minus = type(self)(params_minus).crps_score(y)
    n = len(y)
    grad = np.empty((n, 1))
    grad[:, 0] = (score_plus - score_minus) / (2.0 * eps)
    return grad

crps_metric 

crps_metric() -> NDArray[floating]

Riemannian metric for CRPS natural gradient (MC estimate).

Uses a Monte Carlo estimate: sample from the distribution, compute CRPS gradients, and average the outer product.

RETURNS DESCRIPTION
NDArray[floating]

Metric tensor, shape [n_samples, 1, 1].
Source code in ngboost_lightning/distributions/poisson.py 
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def crps_metric(self) -> NDArray[np.floating]:
    """Riemannian metric for CRPS natural gradient (MC estimate).

    Uses a Monte Carlo estimate: sample from the distribution, compute
    CRPS gradients, and average the outer product.

    Returns:
        Metric tensor, shape ``[n_samples, 1, 1]``.
    """
    n_mc = 50
    n = len(self.rate)
    met = np.zeros((n, 1, 1))
    rng = np.random.default_rng(42)
    for _ in range(n_mc):
        y_sample = rng.poisson(self.rate).astype(np.float64)
        g = self.crps_d_score(y_sample)
        met[:, 0, 0] += g[:, 0] ** 2
    met /= n_mc
    return met