typestar

sm_regression_report.py en Python

Una regresión bien hecha: ajustar, diagnosticar, corregir, comparar, pronosticar.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson, jarque_bera

EFECTO_LENGUAJE = {"python": 9.0, "rust": 5.0, "sql": 0.0}


def crear_tabla(filas=240, semilla=31):
    """Historial de práctica: sesiones, horas, lenguaje, modo editor."""
    rng = np.random.default_rng(semilla)
    tabla = pd.DataFrame({
        "sesiones": rng.integers(1, 80, size=filas),
        "horas": rng.uniform(0.25, 3.0, size=filas),
        "lenguaje": rng.choice(list(EFECTO_LENGUAJE), size=filas),
        "editor": rng.integers(0, 2, size=filas),
    })
    senal = (72
              + 0.35 * tabla["sesiones"]
              + 3.0 * tabla["horas"]
              + tabla["lenguaje"].map(EFECTO_LENGUAJE)
              + 4.0 * tabla["editor"])
    # el ruido crece con las sesiones: errores heterocedásticos
    tabla["tpm"] = senal + rng.normal(0, 1 + tabla["sesiones"] / 25)
    return tabla


def diagnosticar(modelo, tabla):
    print("diagnósticos")
    _, bp_p, _, _ = het_breuschpagan(modelo.resid, modelo.model.exog)
    print(f"  breusch-pagan p {bp_p:.2e}"
          f" -> {'heterocedástico' if bp_p < 0.05 else 'constante'}")
    print(f"  durbin-watson   {durbin_watson(modelo.resid):.3f}")
    print(f"  jarque-bera p   {jarque_bera(modelo.resid)[1]:.4f}")
    return bp_p < 0.05


def principal():
    tabla = crear_tabla()

    base = smf.ols("tpm ~ sesiones + horas", data=tabla).fit()
    completo = smf.ols("tpm ~ sesiones + horas + C(lenguaje) + editor",
                   data=tabla).fit()

    print(f"base r2 {base.rsquared:.4f}   completo r2 {completo.rsquared:.4f}")
    f_stat, f_p, _ = completo.compare_f_test(base)
    print(f"sumar lenguaje y editor: F={f_stat:.2f} p={f_p:.2e}")

    usar_robusto = diagnosticar(completo, tabla)
    res = completo.get_robustcov_results("HC3") if usar_robusto else completo
    etiqueta = "robustos HC3" if usar_robusto else "clásicos"

    print(f"\ncoeficientes (errores estándar {etiqueta})")
    for nombre, coef, err in zip(completo.params.index, res.params, res.bse):
        print(f"  {nombre[:34]:<36}{coef:>8.3f} +/- {err:.3f}")

    print("\nanova (tipo 2)")
    print(anova_lm(completo, typ=2).round(3).to_string())

    futuro = pd.DataFrame({"sesiones": [40, 80], "horas": [1.5, 2.5],
                           "lenguaje": ["python", "rust"], "editor": [1, 0]})
    predicho = completo.get_prediction(futuro).summary_frame(alpha=0.05)
    print("\npronóstico")
    for fila, (media, bajo, alto) in enumerate(zip(
            predicho["mean"], predicho["obs_ci_lower"],
            predicho["obs_ci_upper"])):
        print(f"  fila {fila}: {media:.1f} tpm (95% obs {bajo:.1f}-{alto:.1f})")


if __name__ == "__main__":
    principal()

Cómo funciona

  1. La API de fórmulas mantiene el modelo legible junto a los datos.
  2. Los diagnósticos deciden si hacen falta errores robustos.
  3. Una comparación de modelos anidados justifica los términos extra.

Palabras clave y builtins usados aquí

El intento, en números

Líneas
73
Caracteres a escribir
2563
Tokens
712
Ritmo de tres estrellas
115 tpm

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

Escribe este fragmento

Paso 1 de 1 en Bis; paso 19 de 19 en Estadística con statsmodels.

← Anterior