MEstimator Poisson Example

This page mirrors examples/mestimator_poisson_example.py and shows how to match the built-in Poisson estimator with custom objective and score callbacks.

1 Define The Objective And Scores

import numpy as np
from crabbymetrics import MEstimator, Poisson

np.set_printoptions(precision=4, suppress=True)

call_counter = {"n": 0}


def poisson_objective(theta, data):
    call_counter["n"] += 1

    X = data["X"]
    y = data["y"]
    indices = data.get("indices", np.arange(len(y)))

    X_sample = X[indices]
    y_sample = y[indices]

    eta = X_sample @ theta
    eta = np.clip(eta, -20, 20)
    mu = np.exp(eta)

    obj = np.sum(mu - y_sample * eta)
    grad = X_sample.T @ (mu - y_sample)

    if call_counter["n"] <= 5 or call_counter["n"] % 20 == 0:
        print(
            f"    Call {call_counter['n']}: theta={theta}, "
            f"obj={obj:.2f}, |grad|={np.linalg.norm(grad):.2f}"
        )

    return obj, grad


def poisson_scores(theta, data):
    X = data["X"]
    y = data["y"]

    eta = X @ theta
    eta = np.clip(eta, -20, 20)
    mu = np.exp(eta)

    return X * (mu - y)[:, np.newaxis]

2 Compare MEstimator To Built-In Poisson

rng = np.random.default_rng(42)
n = 700
k = 3
intercept_true = 0.15
beta_true = np.array([0.2, 0.4, -0.6])

X_raw = rng.normal(size=(n, k))
X = np.column_stack([np.ones(n), X_raw])
theta_true = np.concatenate([[intercept_true], beta_true])
eta = X @ theta_true
mu = np.exp(eta)
y = rng.poisson(mu).astype(float)

print("=" * 60)
print("Poisson Regression Comparison")
print("=" * 60)
print(f"True intercept:    {intercept_true}")
print(f"True coefficients: {beta_true}")
print()

print("1. Using MEstimator (general M-estimation framework)")
print("-" * 60)

data = {"X": X, "y": y, "n": n}
theta0 = np.concatenate([[np.log(np.maximum(y.mean(), 1e-12))], np.full(k, 0.1)])

obj_init, grad_init = poisson_objective(theta0, data)
print(f"  Initial objective: {obj_init:.4f}")
print(f"  Initial gradient norm: {np.linalg.norm(grad_init):.4f}")

mestim = MEstimator(
    objective_fn=poisson_objective,
    score_fn=poisson_scores,
    max_iterations=200,
    tolerance=1e-6,
)

mestim.fit(data, theta0)
summary_m = mestim.summary()

print(f"Estimated coefficients: {summary_m['coef']}")
print(f"Standard errors:        {summary_m['se']}")
print()

print("2. Using built-in Poisson estimator")
print("-" * 60)

poisson_model = Poisson(alpha=0.0, max_iterations=200)
poisson_model.fit(X_raw, y)
summary_p = poisson_model.summary()

coef_p = np.concatenate([[summary_p["intercept"]], summary_p["coef"]])
se_p = np.concatenate([[summary_p["intercept_se"]], summary_p["coef_se"]])

print(f"Estimated intercept:    {summary_p['intercept']}")
print(f"Estimated coefficients: {summary_p['coef']}")
print(f"Standard errors:        {se_p}")
print()

print("3. Comparison")
print("-" * 60)
coef_diff = np.abs(summary_m["coef"] - coef_p)
se_diff = np.abs(summary_m["se"] - se_p)

print(f"Max coefficient difference: {np.max(coef_diff):.2e}")
print(f"Max SE difference:          {np.max(se_diff):.2e}")
print(
    f"Max relative coef error:    "
    f"{np.max(coef_diff / np.maximum(np.abs(coef_p), 1e-12)):.2%}"
)
print(f"Max relative SE error:      {np.max(se_diff / se_p):.2%}")

if np.max(coef_diff) < 1e-3 and np.max(se_diff / se_p) < 0.05:
    print("\\nSUCCESS: MEstimator matches built-in Poisson!")
else:
    print("\\nWARNING: Results differ more than expected")

3 Bootstrap The Custom Estimator

print()
print("4. Testing bootstrap (5 iterations)")
print("-" * 60)
boot_draws = mestim.bootstrap(5, seed=123)
print(f"Bootstrap draws shape: {boot_draws.shape}")
print(f"Bootstrap mean: {boot_draws.mean(axis=0)}")
print(f"Bootstrap std:  {boot_draws.std(axis=0)}")