import numpy as np
from crabbymetrics import MEstimator, Poisson
np.set_printoptions(precision=4, suppress=True)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
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")============================================================
Poisson Regression Comparison
============================================================
True intercept: 0.15
True coefficients: [ 0.2 0.4 -0.6]
1. Using MEstimator (general M-estimation framework)
------------------------------------------------------------
Call 1: theta=[0.4474 0.1 0.1 0.1 ], obj=624.59, |grad|=781.46
Initial objective: 624.5944
Initial gradient norm: 781.4648
Call 2: theta=[0.4474 0.1 0.1 0.1 ], obj=624.59, |grad|=781.46
Call 3: theta=[0.4474 0.1 0.1 0.1 ], obj=624.59, |grad|=781.46
Call 4: theta=[ 1.0535 79.3389 279.85 -725.2606], obj=169820124863.99, |grad|=223113598379.23
Call 5: theta=[ 1.0535 79.3389 279.85 -725.2606], obj=169820124863.99, |grad|=223113598379.23
Call 20: theta=[ 0.2757 0.1558 0.3222 -0.4638], obj=343.69, |grad|=142.34
Call 40: theta=[ 0.173 0.1807 0.3903 -0.5907], obj=332.98, |grad|=0.00
Estimated coefficients: [ 0.173 0.1807 0.3903 -0.5907]
Standard errors: [0.0374 0.0307 0.0318 0.03 ]
2. Using built-in Poisson estimator
------------------------------------------------------------
Estimated intercept: 0.1731653097473674
Estimated coefficients: [ 0.1807 0.3902 -0.5905]
Standard errors: [0.0379 0.0305 0.0307 0.0309]
3. Comparison
------------------------------------------------------------
Max coefficient difference: 2.35e-04
Max SE difference: 1.14e-03
Max relative coef error: 0.12%
Max relative SE error: 3.73%
\nSUCCESS: MEstimator matches built-in Poisson!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)}")
4. Testing bootstrap (5 iterations)
------------------------------------------------------------
Call 60: theta=[ 0.1817 0.1899 0.4001 -0.5876], obj=307.03, |grad|=3.10
Call 80: theta=[ 0.186 0.1894 0.3988 -0.5841], obj=307.02, |grad|=0.00
Call 100: theta=[ 0.1933 0.1674 0.3994 -0.576 ], obj=336.76, |grad|=5.15
Call 120: theta=[ 0.1996 0.1683 0.3946 -0.5726], obj=336.74, |grad|=0.00
Call 140: theta=[ 0.1821 0.2107 0.3564 -0.624 ], obj=290.03, |grad|=2.90
Call 160: theta=[ 0.1788 0.2125 0.3579 -0.6255], obj=290.03, |grad|=0.00
Call 180: theta=[ 0.1788 0.2125 0.3579 -0.6255], obj=290.03, |grad|=0.00
Call 200: theta=[ 0.1788 0.2125 0.3579 -0.6255], obj=290.03, |grad|=0.00
Call 220: theta=[ 1.0819 0.6115 -0.421 -0.7755], obj=1967.77, |grad|=4099.39
Call 240: theta=[ 0.2156 0.1933 0.3563 -0.5809], obj=379.07, |grad|=0.01
Call 260: theta=[ 0.127 0.1761 0.3481 -0.5477], obj=405.09, |grad|=60.24
Call 280: theta=[ 0.1443 0.1801 0.3494 -0.5741], obj=404.13, |grad|=0.00
Bootstrap draws shape: (5, 4)
Bootstrap mean: [ 0.1848 0.1887 0.3714 -0.5875]
Bootstrap std: [0.0238 0.0147 0.0209 0.0195]