Tweedie Power Parameter Estimation¶

The Tweedie family Tw_p(mu, phi) has variance function V(mu) = mu^p, where p controls the distribution shape:

p	Distribution
1	Poisson
1 < p < 2	Compound Poisson-Gamma (mass at zero + continuous)
2	Gamma

For insurance severity modelling, p typically falls in (1.1, 1.9) — a mix of exact zeros (no claim) and positive claim amounts.

Estimating p from data is a profile likelihood problem. For each candidate p, estimate_p():

1. Fits the full GLM via PIRLS (warm-started from the previous p)
2. Profiles phi by exact maximum likelihood
3. Evaluates the exact Tweedie likelihood and joint p/phi derivatives, with a diagnosed saddlepoint fallback when exact terms are not certifiable
4. Safeguarded Newton steps follow the joint likelihood valley; exact neighboring profiles certify the winner, with automatic Brent fallback when needed

The Pearson plug-in is explicit (phi_method="pearson") and faster, but it is not a likelihood profile. The zero mass, adjusted for prior weight w_i, is P(Y_i=0) = exp(-w_i*mu_i^(2-p) / ((2-p)*phi)); without zeros, p and phi are weakly separated.

This notebook demonstrates the full pipeline on synthetic data with known p and phi.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from superglm import Spline, SuperGLM, Tweedie, generate_tweedie_cpg

import matplotlib.pyplot as plt import numpy as np import pandas as pd from superglm import Spline, SuperGLM, Tweedie, generate_tweedie_cpg

1. Simulate Tweedie data¶

We generate data from a Compound Poisson-Gamma distribution with known parameters. The true effect is a smooth curve in x with varying exposure.

Here exposure plays two roles consistently:

• it scales the mean through mu_total = mu_rate * exposure
• it scales the effective dispersion through phi_eff = phi / exposure

That means the fitted model must include offset=np.log(exposure) so the mean specification matches the DGP, while sample_weight=exposure carries the exponential-dispersion-model prior-weight convention, not replication counts. Zero-weight observations must be removed consistently from X, y, weights, and the offset before profiling.

Key parameters:

• p = 1.6 — between Poisson (1) and Gamma (2)
• phi = 2.0 — baseline dispersion (scaled by 1/exposure per obs)

In [ ]:

rng = np.random.default_rng(42)
n = 50_000
TRUE_P = 1.6
TRUE_PHI = 2.0

# Covariates
x = rng.uniform(0, 1, n)
exposure = rng.uniform(0.5, 2.0, n)

# True mean: rate per unit exposure × exposure
mu_rate = np.exp(np.log(5.0) + 0.5 * np.sin(2 * np.pi * x))
mu_total = mu_rate * exposure

# Generate Tweedie responses with matching prior-weight convention.
# The fitted model will use offset=log(exposure) and sample_weight=exposure.
y = generate_tweedie_cpg(n, mu=mu_total, phi=TRUE_PHI / exposure, p=TRUE_P, rng=rng)

X = pd.DataFrame({"x": x})

offset = np.log(exposure)
zero_rate = np.mean(y == 0)
print(f"n = {n:,}")
print(f"True p = {TRUE_P}, True phi = {TRUE_PHI}")
print(f"Zero rate: {zero_rate:.1%}")
print(f"y range: [{y.min():.1f}, {y.max():.1f}]")
print(f"y mean: {y.mean():.1f}, y[y>0] mean: {y[y > 0].mean():.1f}")

rng = np.random.default_rng(42) n = 50_000 TRUE_P = 1.6 TRUE_PHI = 2.0 # Covariates x = rng.uniform(0, 1, n) exposure = rng.uniform(0.5, 2.0, n) # True mean: rate per unit exposure × exposure mu_rate = np.exp(np.log(5.0) + 0.5 * np.sin(2 * np.pi * x)) mu_total = mu_rate * exposure # Generate Tweedie responses with matching prior-weight convention. # The fitted model will use offset=log(exposure) and sample_weight=exposure. y = generate_tweedie_cpg(n, mu=mu_total, phi=TRUE_PHI / exposure, p=TRUE_P, rng=rng) X = pd.DataFrame({"x": x}) offset = np.log(exposure) zero_rate = np.mean(y == 0) print(f"n = {n:,}") print(f"True p = {TRUE_P}, True phi = {TRUE_PHI}") print(f"Zero rate: {zero_rate:.1%}") print(f"y range: [{y.min():.1f}, {y.max():.1f}]") print(f"y mean: {y.mean():.1f}, y[y>0] mean: {y[y > 0].mean():.1f}")

2. Estimate p¶

We start with a deliberately wrong initial guess (p=1.2) and let estimate_p() profile over the likelihood. Under the hood, each exact joint evaluation:

1. Fits a full PIRLS at the candidate p, warm-started where possible
2. Profiles phi by exact Newton updates in log(phi)
3. Evaluates the exact Tweedie log-likelihood and full p/phi curvature in one compiled series sweep, recording any bounded saddlepoint fallback

Because exposure scales the true mean, the fit includes offset=np.log(exposure) in addition to sample_weight=exposure. Exact joint ML is the default (phi_method="mle", method="auto"); explicit Brent remains available. The Pearson plug-in is explicit (phi_method="pearson") and intended for faster exploratory work, not likelihood-ratio inference. With fit_mode="reml", REML selects spline smoothing penalties inside each candidate fit; p and phi are still conditional plug-in estimates, not a joint mgcv-style REML fit.

In [ ]:

model = SuperGLM(
    family=Tweedie(p=1.2),  # deliberately wrong initial guess
    features={"x": Spline(n_knots=10)},
)

result = model.estimate_p(
    X, y, sample_weight=exposure, offset=offset, ci_alpha=0.05
)

print(f"Estimated p:   {result.p_hat:.4f}  (true: {TRUE_P})")
print(f"Estimated phi: {result.phi_hat:.4f}  (true: {TRUE_PHI})")
print(f"Converged:     {result.converged}")
print(f"Evaluations:   {result.n_evaluations}")

model = SuperGLM( family=Tweedie(p=1.2), # deliberately wrong initial guess features={"x": Spline(n_knots=10)}, ) result = model.estimate_p( X, y, sample_weight=exposure, offset=offset, ci_alpha=0.05 ) print(f"Estimated p: {result.p_hat:.4f} (true: {TRUE_P})") print(f"Estimated phi: {result.phi_hat:.4f} (true: {TRUE_PHI})") print(f"Converged: {result.converged}") print(f"Evaluations: {result.n_evaluations}")

3. Profile likelihood plot¶

The first profile plot shows the objective difference as a function of p. Each point on this curve required a full GLM fit + phi estimation + density evaluation. The orange dots show where Brent actually evaluated (typically 8–12 points); the smooth curve is filled in afterwards on a fine grid. The confidence interval was requested explicitly with ci_alpha=0.05, so the cached chi-squared cutoff is available here.

In [ ]:

fig = result.profile_plot(alpha=0.05)

# Add true p as a green vertical line
ax = fig.get_axes()[0]
ax.axvline(TRUE_P, color="green", linestyle="--", linewidth=1.5, label=f"True p = {TRUE_P}")
ax.legend()
plt.show()

fig_zoom = result.profile_plot(alpha=0.05)
ax_zoom = fig_zoom.get_axes()[0]
ax_zoom.axvline(TRUE_P, color="green", linestyle="--", linewidth=1.5, label=f"True p = {TRUE_P}")

zoom_xmin, zoom_xmax = 1.55, 1.65
ax_zoom.set_xlim(zoom_xmin, zoom_xmax)

zoom_ymax = None
for line in ax_zoom.get_lines():
    x = np.asarray(line.get_xdata(), dtype=float)
    y = np.asarray(line.get_ydata(), dtype=float)
    if x.shape != y.shape or x.size == 0:
        continue
    mask = np.isfinite(x) & np.isfinite(y) & (x >= zoom_xmin) & (x <= zoom_xmax)
    if np.any(mask):
        line_ymax = float(np.max(y[mask]))
        zoom_ymax = line_ymax if zoom_ymax is None else max(zoom_ymax, line_ymax)

if zoom_ymax is not None:
    ax_zoom.set_ylim(0, zoom_ymax * 1.1 if zoom_ymax > 0 else 1.0)

ax_zoom.set_title("Tweedie p profile likelihood (zoomed)")
ax_zoom.legend()
plt.show()

fig = result.profile_plot(alpha=0.05) # Add true p as a green vertical line ax = fig.get_axes()[0] ax.axvline(TRUE_P, color="green", linestyle="--", linewidth=1.5, label=f"True p = {TRUE_P}") ax.legend() plt.show() fig_zoom = result.profile_plot(alpha=0.05) ax_zoom = fig_zoom.get_axes()[0] ax_zoom.axvline(TRUE_P, color="green", linestyle="--", linewidth=1.5, label=f"True p = {TRUE_P}") zoom_xmin, zoom_xmax = 1.55, 1.65 ax_zoom.set_xlim(zoom_xmin, zoom_xmax) zoom_ymax = None for line in ax_zoom.get_lines(): x = np.asarray(line.get_xdata(), dtype=float) y = np.asarray(line.get_ydata(), dtype=float) if x.shape != y.shape or x.size == 0: continue mask = np.isfinite(x) & np.isfinite(y) & (x >= zoom_xmin) & (x <= zoom_xmax) if np.any(mask): line_ymax = float(np.max(y[mask])) zoom_ymax = line_ymax if zoom_ymax is None else max(zoom_ymax, line_ymax) if zoom_ymax is not None: ax_zoom.set_ylim(0, zoom_ymax * 1.1 if zoom_ymax > 0 else 1.0) ax_zoom.set_title("Tweedie p profile likelihood (zoomed)") ax_zoom.legend() plt.show()

4. Confidence interval¶

ci_alpha=0.05 asks estimate_p() to invert the likelihood-ratio test for an MLE dispersion profile before publishing the fitted model. Each new endpoint probe can require a Brent search + GLM refit, so this remains an explicit and potentially expensive request. result.ci() reads the cache on the detached returned result. Calling it later can extend that result's cache, but it does not mutate the model's independently owned summary cache.

In [ ]:

ci_lo, ci_hi = result.ci(alpha=0.05)

print(f"95% Profile CI: [{ci_lo:.4f}, {ci_hi:.4f}]")
print(f"True p = {TRUE_P} {'inside' if ci_lo <= TRUE_P <= ci_hi else 'OUTSIDE'} CI")
print(f"CI width: {ci_hi - ci_lo:.4f}")

# The interval requested by estimate_p() is cached, so this adds no CI fits.
result.profile_plot(alpha=0.05)

ci_lo, ci_hi = result.ci(alpha=0.05) print(f"95% Profile CI: [{ci_lo:.4f}, {ci_hi:.4f}]") print(f"True p = {TRUE_P} {'inside' if ci_lo <= TRUE_P <= ci_hi else 'OUTSIDE'} CI") print(f"CI width: {ci_hi - ci_lo:.4f}") # The interval requested by estimate_p() is cached, so this adds no CI fits. result.profile_plot(alpha=0.05)

5. Model summary¶

After estimate_p(), the model is refitted at the estimated p. Omitting ci_alpha leaves the summary at CI not computed without doing extra work. Here estimate_p(ci_alpha=0.05) computed the interval before atomic model publication, so the matching summary displays its independently owned copy.

In [ ]:

model.summary(alpha=0.05)

model.summary(alpha=0.05)

6. Effect of sample size on recovery¶

Profile likelihood gets sharper with more data. Let's compare p recovery and CI width across sample sizes.

In [ ]:

sample_sizes = [5_000, 10_000, 20_000, 50_000]
results_table = []

for n_sub in sample_sizes:
    rng_sub = np.random.default_rng(123)
    x_sub = rng_sub.uniform(0, 1, n_sub)
    exp_sub = rng_sub.uniform(0.5, 2.0, n_sub)
    mu_rate_sub = np.exp(np.log(5.0) + 0.5 * np.sin(2 * np.pi * x_sub))
    mu_sub = mu_rate_sub * exp_sub
    y_sub = generate_tweedie_cpg(n_sub, mu=mu_sub, phi=TRUE_PHI / exp_sub, p=TRUE_P, rng=rng_sub)
    X_sub = pd.DataFrame({"x": x_sub})

    m = SuperGLM(family=Tweedie(p=1.5), features={"x": Spline(n_knots=8)})
    r = m.estimate_p(X_sub, y_sub, sample_weight=exp_sub, offset=np.log(exp_sub))
    ci = r.ci(alpha=0.05)
    results_table.append(
        {
            "n": n_sub,
            "p_hat": r.p_hat,
            "phi_hat": r.phi_hat,
            "ci_lo": ci[0],
            "ci_hi": ci[1],
            "ci_width": ci[1] - ci[0],
            "evals": r.n_evaluations,
        }
    )

df_results = pd.DataFrame(results_table)
df_results["error"] = df_results["p_hat"] - TRUE_P
df_results[["n", "p_hat", "error", "ci_lo", "ci_hi", "ci_width", "phi_hat", "evals"]]

sample_sizes = [5_000, 10_000, 20_000, 50_000] results_table = [] for n_sub in sample_sizes: rng_sub = np.random.default_rng(123) x_sub = rng_sub.uniform(0, 1, n_sub) exp_sub = rng_sub.uniform(0.5, 2.0, n_sub) mu_rate_sub = np.exp(np.log(5.0) + 0.5 * np.sin(2 * np.pi * x_sub)) mu_sub = mu_rate_sub * exp_sub y_sub = generate_tweedie_cpg(n_sub, mu=mu_sub, phi=TRUE_PHI / exp_sub, p=TRUE_P, rng=rng_sub) X_sub = pd.DataFrame({"x": x_sub}) m = SuperGLM(family=Tweedie(p=1.5), features={"x": Spline(n_knots=8)}) r = m.estimate_p(X_sub, y_sub, sample_weight=exp_sub, offset=np.log(exp_sub)) ci = r.ci(alpha=0.05) results_table.append( { "n": n_sub, "p_hat": r.p_hat, "phi_hat": r.phi_hat, "ci_lo": ci[0], "ci_hi": ci[1], "ci_width": ci[1] - ci[0], "evals": r.n_evaluations, } ) df_results = pd.DataFrame(results_table) df_results["error"] = df_results["p_hat"] - TRUE_P df_results[["n", "p_hat", "error", "ci_lo", "ci_hi", "ci_width", "phi_hat", "evals"]]

7. Key takeaways¶

• Nested MLE + Brent are the defaults; Pearson is an explicit fast plug-in
• Wright-Bessel density evaluation has a diagnosed saddlepoint fallback
• Warm starts often reduce later PIRLS work, but every candidate is a fit
• ci_alpha publishes a requested interval for matching model summaries; result.ci() otherwise updates only the detached returned result
• CI narrows with more data — compare the ci_width column above
• This simulation uses EDM prior weights (phi_eff = phi / exposure), not frequency/replication weights
• Inspect density and boundary diagnostics; profiles near p=1 and p=2, or at a configured search bound, are naturally unstable