Estimating K₁, K₂, T₁, T₂
Using Hierarchical Bayesian NLSSEstimación de K₁, K₂, T₁, T₂
Mediante NLSS Bayesiano Jerárquico

Banister's IR-Model (1975) describes athletic performance as the superposition of two exponential responses to training load: a slow positive one — fitness — and a fast negative one — fatigue. Skiba & Clarke formalized the core question: how much does a session performed on day \(\tau\) contribute to the projected performance on day \(t\)?El modelo IR de Banister (1975) describe el rendimiento deportivo como la superposición de dos respuestas exponenciales a la carga de entrenamiento: una positiva lenta — fitness — y una negativa rápida — fatiga. Skiba & Clarke formalizaron la pregunta central: ¿cuánto contribuye una sesión realizada el día \(\tau\) al rendimiento proyectado en el día \(t\)?

The curve is biphasic: fatigue dominates immediately (fast decay \(T_2\)); after a zero-crossing the adaptation emerges (slow decay \(T_1\)). Projected performance on day \(t\) is the cumulative convolution of the entire TSS history:La curva es bifásica: la fatiga domina inmediatamente (decaimiento rápido \(T_2\)); después de un cruce por cero emerge la adaptación (decaimiento lento \(T_1\)). El rendimiento proyectado en el día \(t\) es la convolución acumulada de todo el historial de TSS:

The approach proposed by Skiba & Clarke for estimating the four parameters is nonlinear least squares (NLLS): find the set \(\{K_1,K_2,T_1,T_2\}\) that minimizes the discrepancy between predicted performance \(\hat{p}(d;\theta)\) and the maximum power tests collected throughout the season.El enfoque propuesto por Skiba & Clarke para estimar los cuatro parámetros es mínimos cuadrados no lineales (NLLS): encontrar el conjunto \(\{K_1,K_2,T_1,T_2\}\) que minimiza la discrepancia entre el rendimiento predicho \(\hat{p}(d;\theta)\) y los tests de potencia máxima recolectados a lo largo de la temporada.

The proposed improvement requires no additional sensors. It uses the same daily TSS any athlete already records, but applies two fundamental changes to the classical estimator.La mejora propuesta no requiere sensores adicionales. Utiliza el mismo TSS diario que cualquier atleta ya registra, pero aplica dos cambios fundamentales al estimador clásico.

The population prior \(\mathcal{P}_\text{pop}=\mathcal{N}(\mu_\text{pop},\Sigma_\text{pop})\) is built from the 30 calibrated AGMT2 athletes (intervals are ±2σ):El prior poblacional \(\mathcal{P}_\text{pop}=\mathcal{N}(\mu_\text{pop},\Sigma_\text{pop})\) se construye a partir de los 30 atletas AGMT2 calibrados (los intervalos son ±2σ):

The complete algorithm, step by stepEl algoritmo completo, paso a paso

The method replaces the unconstrained NLLS minimization with a regularized optimization via analytic gradient, executed every week on a sliding window. Each component is described with sufficient precision for direct implementation.El método reemplaza la minimización NLLS sin restricciones con una optimización regularizada mediante gradiente analítico, ejecutada cada semana en una ventana deslizante. Cada componente se describe con precisión suficiente para implementación directa.

A · Regularized cost functionA · Función de costo regularizada

The Bayesian method adds a second term penalizing departure from the known cohort distribution:El método Bayesiano añade un segundo término que penaliza el alejamiento de la distribución de la cohorte conocida:

B · Adaptive regularization factor \(\lambda(n)\)B · Factor de regularización adaptativo \(\lambda(n)\)

\(\lambda\) decreases with the number of tests in the current window. With few data points the prior dominates; with many, the athlete progressively individualizes.\(\lambda\) disminuye con el número de pruebas en la ventana actual. Con pocos puntos de datos el prior domina; con muchos, el atleta se individualiza progresivamente.

C · Analytic gradient of \(\mathcal{L}(\theta)\)C · Gradiente analítico de \(\mathcal{L}(\theta)\)

Optimization is performed with L-BFGS-B (Byrd et al., 1995), which requires \(\partial\mathcal{L}/\partial\theta\). The gradient has two summable components:La optimización se realiza con L-BFGS-B (Byrd et al., 1995), que requiere \(\partial\mathcal{L}/\partial\theta\). El gradiente tiene dos componentes sumables:

D · Domain constraints (box constraints)D · Restricciones de dominio (box constraints)

L-BFGS-B enforces per-parameter bounds, eliminating physiologically impossible solutions. Combined with the Mahalanobis prior, they remove virtually all spurious local minima of the unconstrained NLLS.L-BFGS-B impone límites por parámetro, eliminando soluciones fisiológicamente imposibles. Combinados con el prior de Mahalanobis, eliminan virtualmente todos los mínimos locales espurios del NLLS sin restricciones.

E · Sliding window protocol — Algorithm 1E · Protocolo de ventana deslizante — Algoritmo 1

The recalibration algorithm runs every 7 days. The warm start (initializing L-BFGS-B from the previous week's \(\theta\)) reduces iterations from ~200 to ~30 in stable weeks, making the algorithm viable on limited hardware.El algoritmo de recalibración se ejecuta cada 7 días. El warm start (inicializar L-BFGS-B desde el \(\theta\) de la semana anterior) reduce las iteraciones de ~200 a ~30 en semanas estables, haciendo el algoritmo viable en hardware limitado.

PARAMETERS:   L = 90 days (window) · tol = 1e-5 · max_iter = 500
              λ_max = 5.0 · n_half = 4

INITIALIZATION (season day 0):
  θ_current  ←  μ_pop = [2.87, 4.09, 40.64, 6.58]

WEEKLY LOOP  (every 7 days, on day t):

  1. UPDATE WINDOW
       W_t ← { power tests in range [t − L,  t] }
       n   ← |W_t|   // number of tests in window

  2. UPDATE λ
       λ ← λ_max · exp(−n / n_half)    // eq. (6)

  3. IF n == 0  // no tests in window
       θ_current ← μ_pop    // no data → pure prior
       CONTINUE to next cycle

  4. BUILD PREDICTION for each test d ∈ W_t:
       p̂(d;θ) ← p₀ + Σ_{τ<d} TSS(τ)·[K₁·exp(−(d−τ)/T₁) − K₂·exp(−(d−τ)/T₂)]

  5. OPTIMIZE with L-BFGS-B:
       θ_init ← θ_current         // warm start from previous week
       θ_new  ← argmin L(θ)  subject to box constraints
                using analytic gradient eqs. (7)–(9)

  6. CONVERGENCE CHECK & REJECTION:
       IF ‖θ_new − θ_current‖ / ‖θ_current‖ < tol:
         θ_current ← θ_new            // normal convergence
       ELSE IF ‖θ_new − μ_pop‖ > 3·diag(Σ_pop)^0.5:
         θ_current ← θ_current        // REJECT: outside 3σ of prior
       ELSE:
         θ_current ← θ_new

  7. STORE  {t, θ_current, λ, n}  in athlete history

F · Mid-season regime change detectionF · Detección de cambio de régimen intra-temporada

An abrupt jump in parameters between consecutive weeks may signal a real physiological state change (altitude block, illness, form peak). The algorithm detects this via the relative parameter shift:Un salto abrupto en los parámetros entre semanas consecutivas puede señalar un cambio real en el estado fisiológico (bloque de altitud, enfermedad, pico de forma). El algoritmo detecta esto mediante el cambio relativo de parámetros:

Progressive convergence with available testsConvergencia progresiva con tests disponibles

The following four points were raised in external peer review of this methodology. Each identifies a genuine limitation of the current formulation and proposes a concrete mathematical extension. These are open research directions for the next iteration of the model.Los siguientes cuatro puntos fueron planteados en revisión por pares externa de esta metodología. Cada uno identifica una limitación genuina de la formulación actual y propone una extensión matemática concreta. Estas son direcciones de investigación abiertas para la próxima iteración del modelo.

A · Prior representativeness — sub-priors by athlete profileA · Representatividad del prior — sub-priors por perfil de atleta

Critique.Crítica. The model's Bayesian correction depends critically on \(\mu_\text{pop}\) and \(\Sigma_\text{pop}\) from the 30-athlete AGMT2 cohort. If that cohort consists of elite cyclists and the user is an amateur triathlete, the prior may introduce systematic bias rather than correcting it. A single population prior assumes all athletes belong to the same physiological distribution — an assumption that fails across age, sex, sport modality, and training age.La corrección bayesiana del modelo depende críticamente de \(\mu_\text{pop}\) y \(\Sigma_\text{pop}\) de la cohorte AGMT2 de 30 atletas. Si esa cohorte consiste en ciclistas de élite y el usuario es un triatleta aficionado, el prior puede introducir sesgo sistemático en lugar de corregirlo. Un prior poblacional único asume que todos los atletas pertenecen a la misma distribución fisiológica — una suposición que falla en edad, sexo, modalidad deportiva y edad de entrenamiento.

Proposed extension: hierarchical sub-priors.Extensión propuesta: sub-priors jerárquicos. Partition the cohort into \(G\) groups by profile vector \(z_i = [\text{sport},\,\text{level},\,\text{age\_bin},\,\text{sex}]\). For each group \(g\), estimate a separate prior from the athletes in that group:Dividir la cohorte en \(G\) grupos por vector de perfil \(z_i = [\text{sport},\,\text{level},\,\text{age\_bin},\,\text{sex}]\). Para cada grupo \(g\), estimar un prior separado de los atletas en ese grupo:

At inference time, select the matching group prior for athlete \(i\) and fall back to the global prior when \(|g| < n_\text{min}\) (e.g., fewer than 5 athletes in the group):En tiempo de inferencia, seleccionar el prior de grupo coincidente para el atleta \(i\) y retroceder al prior global cuando \(|g| < n_\text{min}\) (por ejemplo, menos de 5 atletas en el grupo):

A further extension uses soft group membership: a mixture prior with weights \(\pi_g\) proportional to the similarity between the new athlete's profile and each group centroid, implemented as a convex combination of group Gaussians. This avoids hard cluster boundaries and handles athletes who span multiple profiles (e.g., elite amateur triathlete with cycling background).Una extensión adicional usa membresía de grupo suave: un prior de mezcla con pesos \(\pi_g\) proporcionales a la similitud entre el perfil del nuevo atleta y cada centroide de grupo, implementado como una combinación convexa de gaussianas de grupo. Esto evita límites de cluster duros y maneja atletas que abarcan múltiples perfiles (por ejemplo, triatleta aficionado de élite con experiencia en ciclismo).

B · Inherent limitations of TSS — intensity-domain separationB · Limitaciones inherentes del TSS — separación por dominio de intensidad

Critique.Crítica. TSS collapses an entire session into a single scalar. Two sessions with identical TSS can elicit very different adaptive responses: 4 hours of endurance base and a 60-minute VO₂max interval session may produce the same TSS yet drive distinct molecular signaling cascades (AMPK vs. mTOR pathways, different fiber recruitment patterns). The model's \(K_1\) and \(K_2\) implicitly average over these differences, potentially creating systematic bias in athletes with polarized or pyramidal periodization.TSS colapsa una sesión completa en un escalar único. Dos sesiones con TSS idéntico pueden provocar respuestas adaptativas muy diferentes: 4 horas de base de resistencia y una sesión de intervalos de 60 minutos VO₂max pueden producir el mismo TSS pero impulsar cascadas de señalización molecular distintas (vías AMPK vs. mTOR, patrones de reclutamiento de fibras diferentes). Los \(K_1\) y \(K_2\) del modelo promedian implícitamente estas diferencias, potencialmente creando sesgo sistemático en atletas con periodización polarizada o piramidal.

Proposed extension: domain-decomposed load signal.Extensión propuesta: señal de carga descompuesta por dominio. Replace the scalar TSS with a vector of intensity-domain loads \(\boldsymbol{s}(\tau) = [s_1(\tau),\, s_2(\tau),\, s_3(\tau)]^\top\) separating zone 1–2 (aerobic base), zone 3 (threshold), and zone 4–5 (VO₂max/anaerobic) contributions, following the three-zone intensity model of Seiler (2010):Reemplazar el TSS escalar con un vector de cargas de dominio de intensidad \(\boldsymbol{s}(\tau) = [s_1(\tau),\, s_2(\tau),\, s_3(\tau)]^\top\) separando zona 1–2 (base aeróbica), zona 3 (umbral), y contribuciones de zona 4–5 (VO₂max/anaeróbica), siguiendo el modelo de intensidad de tres zonas de Seiler (2010):

This expands the parameter space from 4 to 12. The NLLS problem becomes significantly harder, but the Bayesian prior plays an even more critical stabilizing role: zone-specific \(\mu^{(k)}_\text{pop}\) and \(\Sigma^{(k)}_\text{pop}\) can be estimated from athletes with complete intensity-distribution records, and transferred to athletes with only aggregate TSS via a zone-imputation model.Esto expande el espacio de parámetros de 4 a 12. El problema NLLS se vuelve significativamente más difícil, pero el prior bayesiano juega un papel estabilizador aún más crítico: \(\mu^{(k)}_\text{pop}\) y \(\Sigma^{(k)}_\text{pop}\) específicos de zona pueden estimarse de atletas con registros de distribución de intensidad completos, y transferirse a atletas con solo TSS agregado mediante un modelo de imputación de zona.

C · Robustness to outlier tests — Huber loss functionC · Robustez ante tests atípicos — función de pérdida de Huber

Critique.Crítica. Maximum power tests (20', 5', etc.) are inherently noisy: motivation, pre-test nutrition, ambient temperature, and pacing strategy introduce measurement error that can easily shift a test result by 3–5%. Under the current squared-error loss \(\mathcal{L}_\text{data}\), a single anomalous test exerts quadratically larger influence on the gradient — it can dominate the optimization and drive the parameters toward a locally spurious fit.Los tests de potencia máxima (20', 5', etc.) son inherentemente ruidosos: la motivación, nutrición pre-test, temperatura ambiente, y estrategia de ritmo introducen errores de medición que pueden fácilmente desplazar un resultado de test por 3–5%. Bajo la actual pérdida de error cuadrado \(\mathcal{L}_\text{data}\), un único test anómalo ejerce una influencia cuadráticamente mayor en el gradiente — puede dominar la optimización e impulsar los parámetros hacia un ajuste localmente espurio.

Proposed extension: Huber loss.Extensión propuesta: pérdida de Huber. Replace the quadratic loss in \(\mathcal{L}_\text{data}\) with the Huber loss \(\mathcal{H}_\delta\), which is quadratic near zero (preserving gradient sensitivity for small residuals) but linear in the tails (downweighting large residuals):Reemplazar la pérdida cuadrática en \(\mathcal{L}_\text{data}\) con la pérdida de Huber \(\mathcal{H}_\delta\), que es cuadrática cerca de cero (preservando sensibilidad de gradiente para residuos pequeños) pero lineal en las colas (downweighting de residuos grandes):

The robust data term becomes:El término de datos robusto se convierte en:

The gradient of the Huber loss is piecewise, requiring a small modification to eq. (8):El gradiente de la pérdida de Huber es por partes, requiriendo una pequeña modificación a eq. (8):

Calibrating \(\delta\) requires an estimate of the expected test noise. A principled choice is \(\delta = k_\sigma \cdot \hat{\sigma}_r\) where \(\hat{\sigma}_r\) is the empirical standard deviation of residuals from a pilot fit and \(k_\sigma \in [1.0, 1.5]\). For most power meter athletes, a threshold of \(\delta \approx 8\text{–}12\text{ W}\) on a 200–350 W baseline (approximately 3–5% CV) provides good robustness without sacrificing sensitivity to genuine fitness changes.Calibrar \(\delta\) requiere una estimación del ruido de test esperado. Una opción principiada es \(\delta = k_\sigma \cdot \hat{\sigma}_r\) donde \(\hat{\sigma}_r\) es la desviación estándar empírica de residuos de un ajuste piloto y \(k_\sigma \in [1.0, 1.5]\). Para la mayoría de atletas con medidor de potencia, un umbral de \(\delta \approx 8\text{–}12\text{ W}\) en una línea base de 200–350 W (aproximadamente 3–5% CV) proporciona buena robustez sin sacrificar sensibilidad a cambios genuinos de fitness.

D · Parameter covariance — full \(\Sigma_\text{pop}\) vs. diagonal approximationD · Covarianza de parámetros — \(\Sigma_\text{pop}\) completa vs. aproximación diagonal

Critique.Crítica. If \(\Sigma_\text{pop}\) is diagonal (i.e., parameters treated as independent), the Mahalanobis distance degenerates to a sum of independent penalties \(\sum_j ({\theta_j - \mu_j})^2/\sigma_j^2\). This ignores physiologically meaningful correlations between parameters — most notably the \((K_1, T_1)\) relationship: athletes who accumulate fitness rapidly (high \(K_1\)) tend to also lose it more rapidly (lower \(T_1\)), reflecting the faster turnover of their training stimulus.Si \(\Sigma_\text{pop}\) es diagonal (es decir, parámetros tratados como independientes), la distancia de Mahalanobis degenera a una suma de penalizaciones independientes \(\sum_j ({\theta_j - \mu_j})^2/\sigma_j^2\). Esto ignora correlaciones fisiológicamente significativas entre parámetros — más notablemente la relación \((K_1, T_1)\): atletas que acumulan fitness rápidamente (alto \(K_1\)) tienden a también perderlo más rápidamente (menor \(T_1\)), reflejando la rotación más rápida de su estímulo de entrenamiento.

Current implementation status.Estado de implementación actual. The AGMT2 implementation uses the full sample covariance matrix \(\hat{\Sigma}_\text{pop}\) estimated from the 30-athlete cohort:La implementación AGMT2 usa la matriz de covarianza de muestra completa \(\hat{\Sigma}_\text{pop}\) estimada de la cohorte de 30 atletas:

The off-diagonal structure of the current \(\hat{\Sigma}_\text{pop}\) estimated from the AGMT2 cohort shows the dominant correlations:La estructura fuera de la diagonal de la actual \(\hat{\Sigma}_\text{pop}\) estimada de la cohorte AGMT2 muestra las correlaciones dominantes:

Practical consequence of ignoring off-diagonals.Consecuencia práctica de ignorar fuera de la diagonal. Consider an athlete where only \(K_1\) is being updated by a new test. With a diagonal prior, \(K_2\) is completely unconstrained by this update. With the full \(\hat{\Sigma}_\text{pop}\), a high observed \(K_1\) shifts the posterior expectation of \(K_2\) upward by \(\hat\rho_{K_1 K_2} \cdot \hat\sigma_{K_2}/\hat\sigma_{K_1} \cdot (K_1 - \mu_{K_1})\) — exactly as Bayesian inference prescribes. Given the cohort's \(\rho(K_1,K_2) \approx +0.63\), this cross-parameter update is substantial and physiologically meaningful.Considerar un atleta donde solo \(K_1\) está siendo actualizado por un nuevo test. Con un prior diagonal, \(K_2\) es completamente sin restricciones por esta actualización. Con la completa \(\hat{\Sigma}_\text{pop}\), un \(K_1\) observado alto desplaza la expectativa posterior de \(K_2\) hacia arriba por \(\hat\rho_{K_1 K_2} \cdot \hat\sigma_{K_2}/\hat\sigma_{K_1} \cdot (K_1 - \mu_{K_1})\) — exactamente como la inferencia bayesiana prescribe. Dado el \(\rho(K_1,K_2) \approx +0.63\) de la cohorte, esta actualización entre parámetros es sustancial y fisiológicamente significativa.

Caveat: estimation stability with \(N=30\).Advertencia: estabilidad de estimación con \(N=30\). A \(4\times4\) covariance matrix has 10 free parameters. With \(N=30\) observations the condition number of \(\hat\Sigma_\text{pop}\) is sensitive to outliers in the cohort. A regularized estimator (shrinkage toward the diagonal) is recommended:Una matriz de covarianza \(4\times4\) tiene 10 parámetros libres. Con \(N=30\) observaciones el número de condición de \(\hat\Sigma_\text{pop}\) es sensible a outliers en la cohorte. Se recomienda un estimador regularizado (shrinkage hacia la diagonal):

Summary of extensions, implementation readiness, and roadmapResumen de extensiones, estado de implementación y hoja de ruta

The population prior \(\mathcal{P}_\text{pop}=\mathcal{N}(\mu_\text{pop},\Sigma_\text{pop})\) was calibrated from a cohort of 30 athletes monitored through the AGMT2/ednaLabs platform over a minimum of 6 months each. The following table summarizes the cohort demographics and monitoring conditions.La prior poblacional \(\mathcal{P}_\text{pop}=\mathcal{N}(\mu_\text{pop},\Sigma_\text{pop})\) fue calibrada a partir de una cohorte de 30 atletas monitoreados a través de la plataforma AGMT2/ednaLabs durante un mínimo de 6 meses cada uno. La siguiente tabla resume las características demográficas de la cohorte y las condiciones de monitoreo.

The following self-contained Python script implements Algorithm 1 end-to-end, including the full covariance prior, adaptive \(\lambda\), analytic gradient, and weekly recalibration loop. It maps directly to equations (2), (4)–(9), and (15) in the paper. Dependencies are limited to NumPy, SciPy, and pandas.El siguiente script Python autocontenido implementa el Algoritmo 1 de extremo a extremo, incluyendo la prior de covarianza completa, \(\lambda\) adaptativo, gradiente analítico, y el bucle de recalibración semanal. Se asigna directamente a las ecuaciones (2), (4)–(9), y (15) en el artículo. Las dependencias se limitan a NumPy, SciPy, y pandas.

Input formatFormato de entrada

The script expects a CSV file with three columns. Every calendar day must have a TSS value; best20min_watts is left empty on non-test days.El script espera un archivo CSV con tres columnas. Cada día del calendario debe tener un valor de TSS; best20min_watts se deja vacío en días sin tests.

nlss_reference.py

#!/usr/bin/env python3
"""
NLSS — Nonlinear Least Squares with Shrinkage
Reference implementation of the Bayesian K₁K₂T₁T₂ estimation algorithm.

Companion code for:
  "Estimating K₁, K₂, T₁, T₂ Using Hierarchical Bayesian NLSS"
  Gabriel Della Mattia · AGMT2 · ednaLabs · 2026

Input:  CSV with columns  date, tss, best20min_watts
Output: Weekly estimates of θ = [K₁, K₂, T₁, T₂] with λ and test count.
Dependencies: numpy, scipy, pandas
"""

import numpy as np
from scipy.optimize import minimize
import pandas as pd
import sys, os

# ═══════════════════════════════════════════════════════════
# §A — AGMT2 Population Prior  (Table 7 in paper)
# ═══════════════════════════════════════════════════════════

MU_POP = np.array([2.87, 4.09, 40.64, 6.584])     # [K₁, K₂, T₁, T₂]

SD_POP = np.array([0.301, 0.334, 3.094, 0.651])

# Full covariance matrix with empirical correlations
CORR = np.array([
    [ 1.00,  +0.63, -0.41, -0.12],
    [+0.63,  1.00,  +0.09, -0.28],
    [-0.41, +0.09,  1.00,  +0.19],
    [-0.12, -0.28, +0.19,  1.00],
])

SIGMA_POP = np.outer(SD_POP, SD_POP) * CORR     # eq. (15)
SIGMA_INV = np.linalg.inv(SIGMA_POP)               # precision matrix

# ═══════════════════════════════════════════════════════════
# §B — Box Constraints  (Table 2)
# ═══════════════════════════════════════════════════════════

BOUNDS = [
    (0.20,  8.00),    # K₁
    (0.20, 12.00),    # K₂
    (14.0,  80.0),    # T₁ (days)
    ( 3.0,  18.0),    # T₂ (days)
]

# ═══════════════════════════════════════════════════════════
# §C — Adaptive λ  (eq. 6)
# ═══════════════════════════════════════════════════════════

LAMBDA_MAX = 5.0
N_HALF     = 4

def adaptive_lambda(n_tests: int) -> float:
    """eq. (6):  λ(n) = λ_max · exp(−n / n_half)"""
    return LAMBDA_MAX * np.exp(-n_tests / N_HALF)

# ═══════════════════════════════════════════════════════════
# §D — Performance prediction  (eq. 2)
# ═══════════════════════════════════════════════════════════

def predict_performance(theta, tss_history, test_days, p0=0.0):
    """p̂(t) = p₀ + Σ TSS(τ)·[K₁·e^{-(t-τ)/T₁} − K₂·e^{-(t-τ)/T₂}]"""
    K1, K2, T1, T2 = theta
    p_hat = np.empty(len(test_days))
    for i, d in enumerate(test_days):
        s = 0.0
        for tau in range(d):
            delta = d - tau
            s += tss_history[tau] * (K1 * np.exp(-delta / T1) - K2 * np.exp(-delta / T2))
        p_hat[i] = p0 + s
    return p_hat

# ═══════════════════════════════════════════════════════════
# §E — Cost function + analytic gradient  (eqs. 5, 7–9)
# ═══════════════════════════════════════════════════════════

def cost_and_gradient(theta, tss_history, test_days, test_watts, lam, p0=0.0):
    """L(θ) = L_data + λ·L_prior   with   ∇L in a single pass."""
    K1, K2, T1, T2 = theta
    n = len(test_days)
    if n == 0:
        diff = theta - MU_POP
        L_prior = diff @ SIGMA_INV @ diff
        g_prior = 2.0 * SIGMA_INV @ diff
        return lam * L_prior, lam * g_prior

    # Forward pass: residuals + cached exponentials
    residuals = np.empty(n)
    cache_e1, cache_e2, cache_te1, cache_te2 = [], [], [], []

    for i, d in enumerate(test_days):
        s = se1 = se2 = ste1 = ste2 = 0.0
        for tau in range(d):
            delta = d - tau
            w = tss_history[tau]
            exp1 = np.exp(-delta / T1)
            exp2 = np.exp(-delta / T2)
            s   += w * (K1 * exp1 - K2 * exp2)
            se1 += w * exp1
            se2 += w * exp2
            ste1 += w * (delta / (T1*T1)) * exp1
            ste2 += w * (delta / (T2*T2)) * exp2
        residuals[i] = test_watts[i] - (p0 + s)
        cache_e1.append(se1);  cache_e2.append(se2)
        cache_te1.append(ste1); cache_te2.append(ste2)

    # L_data  (eq. 5a)
    L_data = np.sum(residuals**2) / n

    # ∇L_data  (eq. 8)
    g_data = np.zeros(4)
    for i in range(n):
        r = residuals[i]
        g_data[0] += -2.0/n * r * cache_e1[i]       # ∂L/∂K₁
        g_data[1] +=  2.0/n * r * cache_e2[i]       # ∂L/∂K₂
        g_data[2] += -2.0*K1/n * r * cache_te1[i]  # ∂L/∂T₁
        g_data[3] +=  2.0*K2/n * r * cache_te2[i]  # ∂L/∂T₂

    # L_prior + ∇L_prior  (eqs. 5b, 9)
    diff = theta - MU_POP
    L_prior = diff @ SIGMA_INV @ diff
    g_prior = 2.0 * SIGMA_INV @ diff

    return L_data + lam*L_prior, g_data + lam*g_prior

# ═══════════════════════════════════════════════════════════
# §F — Weekly recalibration loop  (Algorithm 1)
# ═══════════════════════════════════════════════════════════

WINDOW_L = 90
DELTA_REL_THRESHOLD = 0.15

def run_nlss(df, p0=0.0, verbose=True):
    """Algorithm 1: weekly recalibration with sliding window."""
    tss   = df['tss'].values.astype(float)
    watts = df['best20min_watts'].values.astype(float)
    n_days = len(df)
    theta_current = MU_POP.copy()
    history = []

    for t in range(6, n_days, 7):       # every 7 days
        w_start = max(0, t - WINDOW_L)
        test_mask = ~np.isnan(watts[w_start:t+1])
        test_idx  = np.where(test_mask)[0] + w_start
        test_w    = watts[w_start:t+1][test_mask]
        n_tests   = len(test_w)
        lam = adaptive_lambda(n_tests)          # eq. (6)

        if n_tests == 0:
            theta_current = MU_POP.copy()       # pure prior
            continue

        # L-BFGS-B with warm start from previous θ
        result = minimize(
            lambda th: cost_and_gradient(th, tss, test_idx, test_w, lam, p0),
            theta_current, method='L-BFGS-B',
            jac=True, bounds=BOUNDS,
            options={'maxiter': 500, 'ftol': 1e-10, 'gtol': 1e-7}
        )
        theta_new = result.x

        # Convergence check + regime detection (eq. 10)
        d_rel = np.linalg.norm(theta_new - theta_current) \
              / (np.linalg.norm(theta_current) + 1e-12)
        if d_rel > DELTA_REL_THRESHOLD:
            lam = adaptive_lambda(0)          # regime change → reset λ
        elif np.any(np.abs(theta_new - MU_POP) > 3*SD_POP):
            theta_new = theta_current        # reject (>3σ)

        theta_current = theta_new.copy()
        history.append({
            'day': t, 'n_tests': n_tests, 'lambda': lam,
            'K1': theta_current[0], 'K2': theta_current[1],
            'T1': theta_current[2], 'T2': theta_current[3],
        })
    return history

# ═══════════════════════════════════════════════════════════
# §G — Main: load CSV and run
# ═══════════════════════════════════════════════════════════

if __name__ == '__main__':
    csv_path = sys.argv[1] if len(sys.argv) > 1 else 'athlete_data.csv'
    df = pd.read_csv(csv_path, parse_dates=['date']).sort_values('date')
    df = df.reset_index(drop=True)

    first_test = df.dropna(subset=['best20min_watts']).iloc[0]['best20min_watts']
    p0 = first_test * 0.5

    history = run_nlss(df, p0=p0)
    pd.DataFrame(history).to_csv(
        csv_path.replace('.csv', '_nlss_results.csv'), index=False
    )

Banister EW, Calvert TW, Savage MV, Bach T. A systems model of training for athletic performance. Australian Journal of Sports Medicine. 1975;7:57–61.

Busso T. Variable dose-response relationship between exercise training and performance. Medicine & Science in Sports & Exercise. 2003;35(7):1188–1195. doi:10.1249/01.MSS.0000074465.13621.37

Byrd RH, Lu P, Nocedal J, Zhu C. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing. 1995;16(5):1190–1208. doi:10.1137/0916069

De Pauw K, Roelands B, Cheung SS, de Geus B, Rietjens G, Meeusen R. Guidelines to classify subject groups in sport-science research. International Journal of Sports Physiology and Performance. 2013;8(2):111–122. doi:10.1123/ijspp.8.2.111

Ledoit O, Wolf M. A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis. 2004;88(2):365–411. doi:10.1016/S0047-259X(03)00096-4

Seiler S. What is best practice for training intensity and duration distribution in endurance athletes? International Journal of Sports Physiology and Performance. 2010;5(3):276–291. doi:10.1123/ijspp.5.3.276

Skiba PF, Clarke DC. The W'bal model: Mathematical background and practical applications. In: Science and Application of High-Intensity Interval Training. Human Kinetics; 2021. (Note: The influence curve formalization used here follows the Banister framework as extended by Skiba & Clarke in multiple applied works.)

Huber PJ. Robust estimation of a location parameter. The Annals of Mathematical Statistics. 1964;35(1):73–101.