Poisson Through MEstimator

Matching coefficients and robust covariance with explicit callbacks

This example fits the same unpenalized Poisson mean model with a native Poisson estimator and a callback-driven MEstimator. It compares sandwich covariance in both implementations. Comparing a score sandwich to Poisson’s default Fisher covariance would compare different inference assumptions.

Objective and scores

The objective and gradient must be derivatives of the same function. Clipping the linear predictor while retaining the unclipped gradient breaks that contract. This moderate-scale example uses the exact exponential throughout.

import numpy as np
import crabbymetrics as cm

rng = np.random.default_rng(42)
n = 700
x = rng.normal(size=(n, 3))
design = np.column_stack([np.ones(n), x])
theta_true = np.array([0.15, 0.2, 0.4, -0.6])
y = rng.poisson(np.exp(design @ theta_true)).astype(float)
data = {"X": design, "y": y, "n": n}

def sample(data):
    indices = data.get("indices", np.arange(data["n"]))
    return data["X"][indices], data["y"][indices]

def objective(theta, data):
    X, y = sample(data)
    eta = X @ theta
    mean = np.exp(eta)
    return float(np.mean(mean - y * eta)), X.T @ (mean - y) / len(y)

def scores(theta, data):
    X, y = sample(data)
    return X * (np.exp(X @ theta) - y)[:, None]

MEstimator does not add an intercept: it is the first column of design and the first element of theta. Both callbacks honor bootstrap indices. A score returns one row per observation, not a pre-aggregated gradient.

Fit and compare

theta0 = np.r_[np.log(y.mean()), np.zeros(x.shape[1])]
custom = cm.MEstimator(objective, scores, max_iterations=300, tolerance=1e-10)
custom.fit(data, theta0)
custom_result = custom.summary()

native = cm.Poisson(alpha=0.0, max_iterations=300, tolerance=1e-9)
native.fit(x, y)
native_result = native.summary(vcov="sandwich")
assert native_result["inference_available"]
native_theta = np.r_[native_result["intercept"], native_result["coef"]]
native_se = np.r_[native_result["intercept_se"], native_result["coef_se"]]

print("Custom coefficients:", custom_result["coef"])
print("Native coefficients:", native_theta)
print("Custom sandwich SE:", custom_result["se"])
print("Native sandwich SE:", native_se)
np.testing.assert_allclose(custom_result["coef"], native_theta, atol=1e-5)
np.testing.assert_allclose(custom_result["se"], native_se, rtol=1e-4, atol=1e-6)

Custom coefficients: [ 0.17295725  0.18071599  0.39032242 -0.59070055]
Native coefficients: [ 0.17295798  0.18071581  0.39032205 -0.59069979]
Custom sandwich SE: [0.03853367 0.03069967 0.02974408 0.0319522 ]
Native sandwich SE: [0.03853363 0.03069967 0.02974407 0.03195219]

In v0.9, fit() computes and stores the covariance before returning. It numerically differentiates the mean scores to obtain the bread; the meat is their uncentered outer product. Invalid scores or a singular bread reject the fit and clear its state. Repeated summaries use that stored covariance.

Bootstrap interface

draws = custom.bootstrap(20, seed=123)
print("Draw shape:", draws.shape)
print("Bootstrap coefficient mean:", draws.mean(axis=0))

Draw shape: (20, 4)
Bootstrap coefficient mean: [ 0.17146892  0.17443593  0.39456996 -0.60330473]

Twenty draws illustrate the interface, not a precise confidence interval. Keep callback data immutable for subsequent bootstrap or explicit covariance recomputation. Prefer the native estimator for this built-in model; callbacks are useful when the objective itself needs to change.