typestar

sp_fit_report.py en Python

Ajustar una curva de aprendizaje, bootstrapear sus parámetros y probar los residuos.

import numpy as np
from scipy import optimize, stats

CEILING, RATE, NOISE = 118.0, 0.11, 3.5


def learning(x, ceiling, rate):
    """Tokens por minuto tras x sesiones: acercamiento exponencial."""
    return ceiling * (1.0 - np.exp(-rate * x))


def make_observations(sesiones, rng):
    return learning(sesiones, CEILING, RATE) + rng.normal(0, NOISE,
                                                          size=sesiones.size)


def fit(sesiones, observadas):
    params, covariance = optimize.curve_fit(
        learning, sesiones, observadas, p0=[100.0, 0.05])
    return params, np.sqrt(np.diag(covariance))


def bootstrap_params(sesiones, observadas, rng, draws=300):
    keep = np.empty((draws, 2))
    n = sesiones.size
    for i in range(draws):
        idx = rng.integers(0, n, size=n)
        try:
            keep[i] = fit(sesiones[idx], observadas[idx])[0]
        except RuntimeError:
            keep[i] = np.nan
    return keep[~np.isnan(keep).any(axis=1)]


def principal():
    rng = np.random.default_rng(21)
    sesiones = np.arange(1, 61, dtype=float)
    observadas = make_observations(sesiones, rng)

    params, errores = fit(sesiones, observadas)
    print("curve_fit")
    print(f"  techo {params[0]:7.3f} +/- {errores[0]:.3f}  (real {CEILING})")
    print(f"  tasa  {params[1]:7.4f} +/- {errores[1]:.4f}  (real {RATE})")

    draws = bootstrap_params(sesiones, observadas, rng)
    low, high = np.percentile(draws[:, 0], [2.5, 97.5])
    print(f"intervalo 95% del techo (bootstrap): {low:.2f} a {high:.2f}"
          f"  ({draws.shape[0]} tiradas útiles)")

    residuos = observadas - learning(sesiones, *params)
    print("residuos")
    print(f"  de {residuos.std(ddof=2):.3f} (el ruido era {NOISE})")
    print(f"  shapiro p {stats.shapiro(residuos).pvalue:.4f}")
    print(f"  correlación {stats.pearsonr(sesiones, residuos)[1]:.4f}")

    at_90 = optimize.root_scalar(
        lambda x: learning(x, *params) - 0.9 * params[0],
        bracket=[1, 500], method="brentq")
    print(f"90% del techo se alcanza en la sesión {at_90.root:.1f}")

    total = float(np.trapezoid(learning(sesiones, *params), sesiones))
    print(f"área bajo la curva ajustada: {total:.0f}")


if __name__ == "__main__":
    principal()

Cómo funciona

  1. curve_fit da los parámetros; la covarianza da su error.
  2. Un bootstrap sobre filas remuestreadas verifica ese error empíricamente.
  3. Después los residuos se prueban por normalidad y estructura.

Palabras clave y builtins usados aquí

El intento, en números

Líneas
66
Caracteres a escribir
1991
Tokens
576
Ritmo de tres estrellas
115 tpm

Al ritmo de tres estrellas de 115 tokens por minuto, este intento toma unos 301 segundos.

Escribe este fragmento

Paso 1 de 1 en Bis; paso 23 de 23 en Computación científica con SciPy.

← Anterior