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 (Dunn & Smyth 2005). For each candidate p, estimate_p():

1. Fits the full GLM via PIRLS (warm-started from the previous p)

2. Takes the maximum-likelihood phi at that fitted mean

3. Scores the candidate by the resulting exact Tweedie log-likelihood

A bounded Brent search over p_bounds picks the power, and the model is refitted there. 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.

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)

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 candidate:

1. Fits a full PIRLS at the candidate p, warm-started where possible

2. Solves for the maximum-likelihood phi by Newton steps in log(phi)

3. Scores the candidate by the exact Tweedie log-likelihood at that phi

Because exposure scales the true mean, the fit includes offset=np.log(exposure) in addition to sample_weight=exposure. 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.

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:   {len(result.evaluations)}")

3. Profile likelihood plot

The profile plot shows the likelihood-ratio statistic at every power the search evaluated (typically 8–12), plus the powers the interval evaluated. Each point is a full GLM fit with its maximum-likelihood phi; the plot performs no new fits. The dashed line is the chi-squared cutoff. The confidence interval was requested explicitly with ci_alpha=0.05, so it is shaded here.

ax = result.profile_plot(alpha=0.05)

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

ax_zoom = result.profile_plot(alpha=0.05)
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 before publishing the fitted model. Each endpoint probe is a GLM refit at that power, 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.

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.

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.

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": len(r.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

• A bounded Brent search over p; phi is the maximum-likelihood dispersion at each candidate’s fitted mean

• 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

• Check result.warnings: an estimate at a search bound, particularly near p=1 and p=2, may have its optimum beyond the bound