H34: Two-Gamma Mixture with Constant Slow Shape¶

An 8-parameter speed model. The mixture mean equals the scheduled speed exactly, and the slow component's mean is structurally bounded below the scheduled speed.

$$ f(V \mid v) = m(v)\;\mathrm{Gamma}(\alpha_1, \theta_1) + (1-m(v))\;\mathrm{Gamma}(\alpha_2, \theta_2) $$

with ensemble mean lock: $m \cdot E[\text{slow}] + (1-m) \cdot E[\text{fast}] = v$

Slow component¶

The slow component has a constant shape $\alpha_1$ and a mean that is a fraction of the scheduled speed:

$$ r(v) = \sigma(c_0 + c_1 v), \qquad \text{mean}_\text{slow} = r \cdot v, \qquad \theta_1 = \frac{r \cdot v}{\alpha_1} $$

Since $\sigma(\cdot) \in (0,1)$, the slow mean is always below $v$. The constant shape is a simplification of H33, which had a speed-dependent CV that turned out to be essentially flat ($d_1 \approx 0$).

Fast component (ensemble mean locked)¶

$$ \alpha_2 = \beta_0(v) \cdot v, \qquad \theta_2 = \frac{c}{\beta_0(v)}, \qquad c = \frac{1 - m \cdot r}{1 - m} $$

Mixture weight¶

$$\text{logit}(m) = e_0 + e_1 v$$

# Parameters Role
1--3 $\beta_0(v; s, e, k)$ Fast-component shape ramp
4--5 $c_0, c_1$ Slow mean ratio $r(v)$
6 $\alpha_1$ Slow shape (constant)
7--8 $e_0, e_1$ Mixture weight

Data and Fit¶

We use multi-span augmented AVL data (spans 1--5) to train. During training, each observation's log-likelihood is weighted by $1/\text{spans}$, so single-span observations (the production regime) count most and longer spans contribute less. This compensates for the fact that multi-span speeds are mean-reverted and tighter than the ballistic single-span distribution the model will face in production.

Parameters are loaded from a full-dataset weighted fit.

In [1]:
import json
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

from transit_speed.calibration.dataset import CalibrationDataset
from transit_speed.calibration.metrics import compute_calibration_metrics
from transit_speed.calibration.hypotheses.h34_2gamma_ratio_constcv import (
    H34_2GammaRatioConstCV, _slow_params, _ensemble_locked_params,
)
from transit_speed.speed_model import beta_0_vec

# Load pre-fitted parameters (weighted MLE: weight=1/spans)
with open("h34_weighted_params.json") as f:
    params = json.load(f)

a1 = np.exp(params["log_a1"])
print("H34 fitted parameters (full multi-span, weighted 1/spans):")
print(f"  Fast beta0:  start={params['start_b0']:.4f}  end={params['end_b0']:.4f}  kink={params['kink']:.1f}")
print(f"  Slow ratio:  logit(r) = {params['r_intercept']:.4f} + {params['r_slope']:.6f} * v")
print(f"  Slow shape:  alpha_1  = {a1:.4f}  (cv = {1/np.sqrt(a1):.2f})")
print(f"  Mix weight:  logit(m) = {params['mw_intercept']:.4f} + {params['mw_slope']:.6f} * v")

print("\nMixture at representative speeds:")
for v in [5, 15, 25, 40]:
    _a1, s1, r, ms = _slow_params(
        np.array([v]), params["r_intercept"], params["r_slope"], params["log_a1"])
    mw = 1.0 / (1.0 + np.exp(-(params["mw_intercept"] + params["mw_slope"] * v)))
    a2, s2, c = _ensemble_locked_params(
        np.array([v]), ms, np.array([mw]),
        params["start_b0"], params["end_b0"], params["kink"])
    print(f"  v={v:2d}: mw={mw:.3f} r={r[0]:.3f} c={c[0]:.3f}"
          f"  slow(mean={ms[0]:.1f})  fast(mean={a2[0]*s2[0]:.1f})")

H34 fitted parameters (full multi-span, weighted 1/spans):
  Fast beta0:  start=0.5714  end=0.1244  kink=21.9
  Slow ratio:  logit(r) = 0.0060 + 0.026456 * v
  Slow shape:  alpha_1  = 0.2664  (cv = 1.94)
  Mix weight:  logit(m) = -1.3474 + -0.018838 * v

Mixture at representative speeds:
  v= 5: mw=0.191 r=0.535 c=1.110  slow(mean=2.7)  fast(mean=5.6)
  v=15: mw=0.164 r=0.599 c=1.078  slow(mean=9.0)  fast(mean=16.2)
  v=25: mw=0.140 r=0.661 c=1.055  slow(mean=16.5)  fast(mean=26.4)
  v=40: mw=0.109 r=0.744 c=1.031  slow(mean=29.7)  fast(mean=41.3)

Parameter Variation with Scheduled Speed¶

The two key derived quantities -- slow mean ratio $r(v)$ and mixture weight $m(v)$ -- are smooth functions of scheduled speed. The slow shape $\alpha_1$ is a constant. Below we plot how each varies across the range of observed scheduled speeds.

In [2]:
v_range = np.linspace(1, 50, 200)

# Compute derived quantities across scheduled speed
a1_val = np.exp(params["log_a1"])
a1_arr, s1_arr, r_arr, ms_arr = _slow_params(
    v_range, params["r_intercept"], params["r_slope"], params["log_a1"])
mw_arr = 1.0 / (1.0 + np.exp(-(params["mw_intercept"] + params["mw_slope"] * v_range)))

fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))

# --- Slow mean ratio r(v) ---
ax = axes[0]
ax.plot(v_range, r_arr, "tab:orange", linewidth=2)
ax.set_xlabel("Scheduled speed (mph)")
ax.set_ylabel("$r(v)$")
ax.set_title("Slow Mean Ratio $r(v) = \\mathrm{mean_{slow}} / v$")
ax.set_ylim(0, 1)
ax.axhline(1.0, color="gray", linestyle=":", alpha=0.5)
ax.grid(True, alpha=0.3)

# Add secondary y-axis showing absolute slow mean
ax2 = ax.twinx()
ax2.plot(v_range, ms_arr, "tab:orange", linewidth=1, linestyle="--", alpha=0.5)
ax2.set_ylabel("Slow mean (mph)", color="tab:orange", alpha=0.5)
ax2.tick_params(axis="y", labelcolor="tab:orange", colors="tab:orange")

# --- Slow shape (constant) ---
ax = axes[1]
ax.axhline(a1_val, color="tab:red", linewidth=2, label=f"$\\alpha_1$ = {a1_val:.3f}")
ax.set_xlabel("Scheduled speed (mph)")
ax.set_ylabel("$\\alpha_1$")
ax.set_title(f"Slow Shape $\\alpha_1$ = {a1_val:.3f}  (CV = {1/np.sqrt(a1_val):.2f})")
ax.set_ylim(0, 0.5)
ax.legend()
ax.grid(True, alpha=0.3)

# --- Mixture weight m(v) ---
ax = axes[2]
ax.plot(v_range, mw_arr, "tab:blue", linewidth=2)
ax.set_xlabel("Scheduled speed (mph)")
ax.set_ylabel("$m(v)$")
ax.set_title("Slow Component Weight $m(v)$")
ax.set_ylim(0, 0.25)
ax.grid(True, alpha=0.3)

fig.suptitle("H34 Parameter Variation with Scheduled Speed", fontsize=13)
plt.tight_layout()
plt.show()
No description has been provided for this image

Distribution Shapes¶

In [3]:
def h34_pdf(v, params, x):
    """Return (pdf_mix, pdf_slow, pdf_fast, mw, c) at scheduled speed v."""
    a1, s1, r, ms = _slow_params(
        np.atleast_1d(v), params["r_intercept"], params["r_slope"], params["log_a1"])
    mw = 1.0 / (1.0 + np.exp(-(params["mw_intercept"] + params["mw_slope"] * v)))
    a2, s2, c = _ensemble_locked_params(
        np.atleast_1d(v), ms, np.atleast_1d(mw),
        params["start_b0"], params["end_b0"], params["kink"])
    # a1 is scalar (constant shape); s1, a2, s2, c are arrays
    a1 = float(a1)
    s1, a2, s2, c = s1[0], a2[0], s2[0], c[0]
    pdf_slow = stats.gamma.pdf(x, a=a1, scale=s1)
    pdf_fast = stats.gamma.pdf(x, a=a2, scale=s2)
    pdf_mix = mw * pdf_slow + (1 - mw) * pdf_fast
    return pdf_mix, pdf_slow, pdf_fast, mw, c

sched_values = [5, 15, 25, 40]
fig, axes = plt.subplots(1, 4, figsize=(20, 4))

for ax, v in zip(axes, sched_values):
    x = np.linspace(0.01, v * 3.5, 500)
    pdf_mix, pdf_slow, pdf_fast, mw, c = h34_pdf(v, params, x)

    ax.plot(x, pdf_mix, "k-", linewidth=2, label="Mixture")
    ax.plot(x, mw * pdf_slow, "--", color="tab:orange", linewidth=1.2,
            label=f"Slow (m={mw:.2f})")
    ax.plot(x, (1 - mw) * pdf_fast, "--", color="tab:blue", linewidth=1.2,
            label=f"Fast (1-m={1-mw:.2f})")
    ax.axvline(v, color="tab:red", linestyle=":", alpha=0.7, label="$v_{sched}$ = E[V]")

    ax.set_title(f"$v$ = {v} mph  (c={c:.3f})", fontsize=10)
    ax.set_xlabel("Speed (mph)")
    ax.set_xlim(0, v * 3.5)
    ax.grid(True, alpha=0.2)
    if ax is axes[0]:
        ax.set_ylabel("Density")
        ax.legend(fontsize=8)

fig.suptitle("H34 Mixture PDF (ensemble mean = $v_{sched}$)", fontsize=13, y=1.03)
plt.tight_layout()
plt.show()
No description has been provided for this image

Comparison to Data¶

In [4]:
# Load data and build multi-span augmentation
dataset = CalibrationDataset(
    csv_path="data/avl_log_20260310_030816.csv",
    gtfs_dir="data/gtfs",
)
ms = dataset.build_multispans(max_spans=5)

print(f"Single-span: {len(dataset.full_df):,}")
print(f"Multi-span:  {len(ms.full_df):,}")
print(f"  Train: {len(ms.train_df):,}  Test: {len(ms.test_df):,}")
print(f"\nSpan counts in test set:")
print(ms.test_df["spans"].value_counts().sort_index().to_string())

Single-span: 631,400
Multi-span:  3,018,109
  Train: 2,147,198  Test: 870,911

Span counts in test set:
spans
1    182517
2    176558
3    173582
4    170610
5    167644
In [5]:
speed_bands = [(4, 5), (7, 8), (8, 14), (14, 22), (22, 35), (35, 55)]
band_labels = ["4-5", "7-8", "8-14", "14-22", "22-35", "35-55"]

# Use single-span test data for comparison to model PDFs
test_1span = ms.test_df[ms.test_df["spans"] == 1]

fig, axes = plt.subplots(len(speed_bands), 1, figsize=(10, 3 * len(speed_bands)))

for ax, (lo, hi), label in zip(axes, speed_bands, band_labels):
    mask = (test_1span["scheduled_speed_mph"] >= lo) & (test_1span["scheduled_speed_mph"] < hi)
    obs = test_1span.loc[mask, "speed_mph"].values
    if len(obs) < 50:
        continue

    v_mid = (lo + hi) / 2
    x = np.linspace(0.01, hi * 3, 500)
    pdf_mix, pdf_slow, pdf_fast, mw, c = h34_pdf(v_mid, params, x)

    bins = np.linspace(0, 80, 150)
    ax.hist(obs, bins=bins, density=True, alpha=0.4, color="gray", edgecolor="none",
            label=f"Data (n={len(obs):,})")
    ax.plot(x, pdf_mix, "k-", linewidth=2, label="H34")
    ax.plot(x, mw * pdf_slow, "--", color="tab:orange", linewidth=1)
    ax.plot(x, (1 - mw) * pdf_fast, "--", color="tab:blue", linewidth=1)

    ax.set_title(f"Scheduled {label} mph  (n={len(obs):,})", fontsize=10)
    ax.set_xlim(0, 80)
    ax.set_ylim(0, 0.15)
    ax.set_ylabel("Density")
    ax.grid(True, alpha=0.2)
    if ax is axes[0]:
        ax.legend(fontsize=8)

axes[-1].set_xlabel("Speed (mph)")
plt.tight_layout()
plt.show()
No description has been provided for this image

Calibration by Span Count¶

The model is trained on multi-span augmented data (spans 1--5), but it does not condition on $\Delta t$ or span count. A key question is whether a single set of parameters produces well-calibrated predictions across different observation timescales:

  • 1-span (~25--30 s): raw polling intervals
  • 3-span (~75--90 s): moderate extrapolation horizon
  • 5-span (~125--150 s): the stalest observations

If calibration degrades at longer spans, it would suggest the model needs $\Delta t$-dependence. If it holds, the two-gamma mixture is capturing the speed distribution well regardless of averaging window.

In [6]:
# Evaluate PIT values on full test set
h34 = H34_2GammaRatioConstCV()
pit_all, logpdf_all = h34._evaluate(ms.test_df, params)

# Split by span count
span_groups = {1: "1-span (~30 s)", 3: "3-span (~90 s)", 5: "5-span (~150 s)"}
span_metrics = {}

for span, label in span_groups.items():
    mask = ms.test_df["spans"].values == span
    pit_span = pit_all[mask]
    pit_clean = pit_span[np.isfinite(pit_span)]
    m = compute_calibration_metrics(pit_clean)
    span_metrics[span] = m
    print(f"{label}:  n={m['n']:>10,d}  KS={m['ks_stat']:.4f}  mean_PIT={m['mean_pit']:.4f}")

# Also compute overall
pit_clean_all = pit_all[np.isfinite(pit_all)]
m_all = compute_calibration_metrics(pit_clean_all)
print(f"\n{'All spans':18s}:  n={m_all['n']:>10,d}  KS={m_all['ks_stat']:.4f}  mean_PIT={m_all['mean_pit']:.4f}")

1-span (~30 s):  n=   182,517  KS=0.0922  mean_PIT=0.4893
3-span (~90 s):  n=   173,582  KS=0.0882  mean_PIT=0.5255
5-span (~150 s):  n=   167,644  KS=0.1426  mean_PIT=0.5348

All spans         :  n=   870,911  KS=0.0636  mean_PIT=0.5188
In [7]:
fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)

for ax, (span, label) in zip(axes, span_groups.items()):
    mask = ms.test_df["spans"].values == span
    pit_span = pit_all[mask]
    pit_clean = pit_span[np.isfinite(pit_span)]
    m = span_metrics[span]

    ax.hist(pit_clean, bins=50, density=True, alpha=0.7, edgecolor="white")
    ax.axhline(1.0, color="tab:red", linestyle="--", linewidth=1.5)
    ax.set_xlabel("PIT value")
    ax.set_title(f"{label}\nKS={m['ks_stat']:.4f}  mean={m['mean_pit']:.4f}  n={m['n']:,}")
    ax.set_xlim(0, 1)
    ax.grid(True, alpha=0.3)

axes[0].set_ylabel("Density")
fig.suptitle("H34 PIT Histograms by Span Count", fontsize=13)
plt.tight_layout()
plt.show()
No description has been provided for this image

PIT ECDFs overlaid¶

A direct overlay of the PIT ECDFs makes it easier to see how calibration shifts across span counts.

In [8]:
fig, ax = plt.subplots(figsize=(8, 5))

colors = {"1-span (~30 s)": "tab:blue", "3-span (~90 s)": "tab:orange", "5-span (~150 s)": "tab:green"}

for span, label in span_groups.items():
    mask = ms.test_df["spans"].values == span
    pit_span = pit_all[mask]
    pit_clean = np.sort(pit_span[np.isfinite(pit_span)])
    ecdf = np.arange(1, len(pit_clean) + 1) / len(pit_clean)
    # Subsample for plotting performance
    step = max(1, len(pit_clean) // 5000)
    ax.plot(pit_clean[::step], ecdf[::step], linewidth=1.2, color=colors[label],
            label=f"{label} (KS={span_metrics[span]['ks_stat']:.4f})")

ax.plot([0, 1], [0, 1], "k--", linewidth=1, alpha=0.5, label="Ideal (uniform)")
ax.set_xlabel("PIT value")
ax.set_ylabel("Cumulative probability")
ax.set_title("PIT ECDFs by Span Count")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
No description has been provided for this image

Empirical speed distributions by span count¶

The model does not condition on span count, yet the observed speed distribution changes: longer spans average out momentary stops and accelerations, producing a tighter distribution. Below we compare the empirical speed distributions at 1, 3, and 5 spans against the model PDF for a few representative scheduled speeds.

In [9]:
sched_bands = [(8, 14), (14, 22), (22, 35)]
span_list = [1, 3, 5]
span_colors = {1: "tab:blue", 3: "tab:orange", 5: "tab:green"}

fig, axes = plt.subplots(1, len(sched_bands), figsize=(18, 5))

for ax, (lo, hi) in zip(axes, sched_bands):
    v_mid = (lo + hi) / 2
    x = np.linspace(0.01, hi * 3, 500)
    pdf_mix, _, _, _, _ = h34_pdf(v_mid, params, x)

    ax.plot(x, pdf_mix, "k-", linewidth=2.5, label="H34 model", zorder=10)

    bins = np.linspace(0, hi * 3, 100)
    for span in span_list:
        sub = ms.test_df[ms.test_df["spans"] == span]
        mask = (sub["scheduled_speed_mph"] >= lo) & (sub["scheduled_speed_mph"] < hi)
        obs = sub.loc[mask, "speed_mph"].values
        if len(obs) < 50:
            continue
        ax.hist(obs, bins=bins, density=True, alpha=0.25, color=span_colors[span],
                edgecolor="none", label=f"{span}-span (n={len(obs):,})")

    ax.set_title(f"Scheduled {lo}--{hi} mph")
    ax.set_xlabel("Speed (mph)")
    ax.set_xlim(0, hi * 3)
    ax.set_ylim(0, 0.12)
    ax.grid(True, alpha=0.2)
    ax.legend(fontsize=8)

axes[0].set_ylabel("Density")
fig.suptitle("Empirical Speed Distributions by Span Count vs. H34 Model", fontsize=13)
plt.tight_layout()
plt.show()
No description has been provided for this image

Summary¶

H34 is an 8-parameter two-gamma mixture speed model with:

  1. Ensemble mean lock: the mixture mean equals the scheduled speed exactly, providing an unbiased position extrapolation.

  2. Mean-ratio slow component: the slow component's mean is $r(v) \cdot v$ where $r \in (0, 1)$ by construction (sigmoid), guaranteeing the slow mean stays below the scheduled speed at all $v$.

  3. Constant slow shape: $\alpha_1$ is a single fitted constant (~0.27), simplifying H33's speed-dependent CV which turned out to be flat. This removes one parameter with no loss in fit quality.

  4. Span-weighted training: observations are weighted by $1/\text{spans}$ during MLE, emphasizing the single-span (production) regime over the mean-reverted multi-span regime.

In [ ]: