Binding Internals: Poisson

Native likelihood optimization and inference

Unlike the closed-form OLS path, Poisson defines an objective, gradient, and Hessian in Rust. Python supplies arrays, not solver callbacks. The code is in src/estimators/mle.rs; the Poisson reference documents the public API.

Objective and solver

With \(\eta_i=b_0+x_i'\beta\) and \(\mu_i=\exp(\eta_i)\), the objective is

\[ \sum_i\{\mu_i-y_i\eta_i\}+\frac{\alpha}{2}\|\beta\|^2. \]

The intercept is not penalized. alpha is the slope penalty, not an intercept. The data-only log-factorial term is unnecessary for optimization.

fit() clears previous state, checks finite nonnegative outcomes and valid solver controls, and initializes slopes at zero and the intercept near \(\log(\bar y)\). Newton-CG uses analytic derivatives and a More-Thuente line search. Shared convergence checks reject a budget-exhausted fit instead of presenting it as successful. Valid coefficients, training arrays, and solver diagnostics are stored in Rust.

Prediction and inference

predict_lin(x) returns \(\eta\); predict(x) returns \(\exp(\eta)\). Prediction validates feature width and finite inputs, but extreme linear predictors can still overflow on exponentiation.

For alpha=0, summary(vcov="vanilla") uses inverse Fisher information and summary(vcov="sandwich") uses the individual score outer products as the covariance meat. The latter supports Poisson QMLE when the conditional mean is correct but the conditional variance is not Poisson.

Inference is marked unavailable for penalized fits or rank-deficient unpenalized information. No artificial diagonal ridge is added to make those standard errors look available. Check inference_available before accessing the inference fields.

Code

import numpy as np
import crabbymetrics as cm

rng = np.random.default_rng(19)
x = rng.normal(size=(500, 2))
y = rng.poisson(np.exp(0.2 + x @ np.array([0.3, -0.2]))).astype(float)
model = cm.Poisson(max_iterations=300, tolerance=1e-7)
model.fit(x, y)
result = model.summary(vcov="sandwich")
assert result["converged"] and result["inference_available"]
print("Coefficients:", result["coef"])
print("Sandwich SE:", result["coef_se"])
print("Predicted means:", model.predict(x[:3]))
np.testing.assert_allclose(model.predict(x[:3]), np.exp(model.predict_lin(x[:3])))
print("Bootstrap shape:", model.bootstrap(20, seed=19).shape)

Coefficients: [ 0.23770838 -0.18850806]
Sandwich SE: [0.0406713  0.03861668]
Predicted means: [0.93664738 1.53002348 1.90221872]
Bootstrap shape: (20, 3)

Bootstrap refits stay in Rust and stream row-index draws. Solver work and training-data copies are still repeated for each replicate. The custom Poisson example shows the same model through MEstimator when a Python-defined objective is needed.