from _api_doc_utils import *Poisson
Poisson GLM for count outcomes
1 Where it fits
Group: Regression
Poisson fits
\[ \mathbb E[Y_i\mid X_i=x_i]=\exp(\alpha+x_i'\beta). \]
The alpha constructor argument is an L2 penalty, not the intercept. For alpha=0, summary(vcov='vanilla') reports Fisher-information standard errors and summary(vcov='sandwich') reports robust QMLE-style standard errors. Penalized fits mark inference unavailable and omit se/vcov.
2 Likelihood and solver
With \(\eta_i=\alpha_0+x_i'\beta\) and \(\mu_i=\exp(\eta_i)\), the implemented criterion is the Poisson negative log-likelihood up to the data-only \(\log(y_i!)\) term:
\[ Q(\alpha_0,\beta) = \sum_{i=1}^n\{\mu_i-y_i\eta_i\} +\frac{\lambda}{2}\|\beta\|_2^2. \]
The intercept is unpenalized. Newton-CG uses the exact gradient and Hessian
\[ \nabla^2 Q = \tilde X'\operatorname{diag}(\mu_i)\tilde X +\operatorname{diag}(0,\lambda,\ldots,\lambda), \]
with \(10^{-8}\) added to the Hessian diagonal for numerical stability and a More-Thuente line search. The intercept starts at \(\log\{\max(\bar y,10^{-12})\}\) and slopes start at zero. A fit is stored only after Argmin reports solver convergence or a target cost. Reaching max_iterations raises ValueError and clears any previous fitted state. Successful summaries include converged, iterations, termination_reason, and the final penalized objective \(Q(\hat\alpha_0,\hat\beta)\).
Warning
The fit requires finite \(y_i\geq0\), a finite nonnegative penalty, finite regressors, a positive tolerance, and a positive iteration budget. Outcomes need not be integers, which permits Poisson QMLE for nonnegative responses.
3 Inference
Analytic inference is disabled for \(\lambda>0\). For an unpenalized fit, the vanilla covariance is the inverse Fisher information
\[ \widehat V_{\mathrm{vanilla}} = (\tilde X'\operatorname{diag}(\hat\mu_i)\tilde X+10^{-8}I)^{-1}. \]
The sandwich, also accepted under the alias QMLE, is
\[ \widehat V_{\mathrm{QMLE}} = B^{-1} \left\{ \sum_i \tilde x_i\tilde x_i'(y_i-\hat\mu_i)^2 \right\} B^{-1}, \qquad B=\tilde X'\operatorname{diag}(\hat\mu_i)\tilde X+10^{-8}I. \]
No HC finite-sample or cluster correction is applied. Wald tests use either covariance. The pairs bootstrap refits the fixed penalty in every draw, and any nonconverged replicate aborts the bootstrap.
4 Performance and numerical behavior
Every Newton Hessian evaluation costs \(O(np^2)\) and the dense Newton system is up to \(O(p^3)\), so this is more expensive per iteration than first-order GLM solvers at large \(p\). The code evaluates \(\exp(\eta)\) without clipping; extreme indices can overflow and abort optimization. There is no exposure or offset argument, so users must encode those effects in the supplied design only if estimating their coefficient is appropriate. Covariance and bootstrap add dense inversion and repeated-fit costs.
5 Python API
Constructor: cm.Poisson
Use fit(x, y) with nonnegative count-like outcomes. predict_lin(x) returns the log mean, predict(x) returns fitted conditional means, and summary() includes fit diagnostics. The class also supports bootstrap(B, seed=None).
print(inspect.signature(cm.Poisson))(alpha=0.0, max_iterations=100, tolerance=0.0001)cls = cm.Poisson
display(HTML(html_table(["Public method"], public_methods(cls))))| Public method |
|---|
bootstrap(self, /, n_bootstrap, seed=None) |
fit(self, /, x, y) |
predict(self, /, x) |
predict_lin(self, /, x) |
summary(self, /, vcov='vanilla') |
wald_test(self, /, r, q=None, vcov='vanilla') |
6 Minimal example
rng = np.random.default_rng(7)
x = rng.normal(size=(250, 2))
mu = np.exp(0.2 + x @ np.array([0.4, -0.25]))
y = rng.poisson(mu).astype(float)
model = cm.Poisson(max_iterations=200, tolerance=1e-08)
model.fit(x, y)
fit = model.summary(vcov='vanilla')
print({key: fit[key] for key in ['converged', 'iterations', 'termination_reason', 'objective']})
print(fit['coef'])
print(model.summary(vcov='sandwich')['coef_se'])
print(model.predict(x[:3])){'converged': True, 'iterations': 9, 'termination_reason': 'Solver converged', 'objective': 223.22902707069227}
[ 0.26320223 -0.29275249]
[0.0578437 0.05367217]
[1.01940438 1.34302229 1.319153 ]7 summary() contract
The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.
rng = np.random.default_rng(107)
x = rng.normal(size=(100, 2))
y = rng.poisson(np.exp(0.2 + x @ np.array([0.4, -0.25]))).astype(float)
model = cm.Poisson(max_iterations=200)
model.fit(x, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))| summary() key | shape |
|---|---|
intercept |
() |
coef |
(2,) |
penalty |
() |
converged |
() |
iterations |
() |
termination_reason |
() |
objective |
() |
inference_available |
() |
intercept_se |
() |
coef_se |
(2,) |
vcov |
(3, 3) |
vcov_type |
() |