Regresión lineal y descenso de gradiente

import numpy as np
import matplotlib.pyplot as plt

# Semilla fija: cada vez que corras esto saldrán los mismos datos
rng = np.random.default_rng(42)

n = 30                                # cantidad de estudiantes
# Horas de estudio en la semana del examen, entre 0 y 8
horas = rng.uniform(0, 8, size=n)
# Todo lo que las horas no explican (sueño, nervios, suerte)
ruido = rng.normal(0, 6, size=n)
# La regla «secreta»: 7 puntos por hora, partiendo de 30
nota = 30 + 7 * horas + ruido

print("Horas:", np.round(horas[:5], 1))
print("Notas:", np.round(nota[:5], 1))

plt.scatter(horas, nota)
plt.xlabel("Horas de estudio")
plt.ylabel("Nota del examen (sobre 100)")
plt.title("Treinta estudiantes (datos sintéticos)")
plt.show()

def predecir(x, w, b):
    """La predicción del modelo: una recta."""
    return w * x + b

# Una primera recta elegida «a ojo»
w, b = 5.0, 40.0
prediccion = predecir(horas, w, b)
residuos = nota - prediccion          # real menos predicho

for i in range(3):
    print(f"Estudiante {i}: estudió {horas[i]:.1f} h, "
          f"sacó {nota[i]:.1f}, la recta dice {prediccion[i]:.1f}, "
          f"residuo {residuos[i]:+.1f}")

# Dos puntos bastan para dibujar una recta
x_recta = np.array([0, 8])
plt.scatter(horas, nota)
plt.plot(x_recta, predecir(x_recta, w, b), color="C0",
         label="Recta a ojo: 5·x + 40")
plt.vlines(horas, prediccion, nota, color="C2", linewidth=1.2,
           label="Residuos")
plt.xlabel("Horas de estudio")
plt.ylabel("Nota del examen (sobre 100)")
plt.title("Cada residuo es un error de la recta")
plt.legend()
plt.show()

def mse(y, y_pred):
    """Error cuadrático medio."""
    return np.mean((y - y_pred) ** 2)

def mae(y, y_pred):
    """Error absoluto medio."""
    return np.mean(np.abs(y - y_pred))

for w, b in [(5, 40), (7, 30), (9, 20)]:
    p = predecir(horas, w, b)
    error = mse(nota, p)
    print(f"w = {w}, b = {b}:   MSE = {error:4.1f}   "
          f"RMSE = {np.sqrt(error):3.1f}   MAE = {mae(nota, p):3.1f}")

# 200 pendientes posibles entre 0 y 14
ws = np.linspace(0, 14, 200)
perdidas = [mse(nota, predecir(horas, w, 30)) for w in ws]
mejor = ws[np.argmin(perdidas)]       # la pendiente con menor pérdida

plt.plot(ws, perdidas)
plt.scatter([mejor], [min(perdidas)], color="C1", s=140, zorder=3,
            label="El punto más bajo")
plt.xlabel("Pendiente w (con b = 30 fijo)")
plt.ylabel("MSE")
plt.title("Un corte del valle: la pérdida según w")
plt.legend()
plt.show()

print(f"Con b = 30, la mejor pendiente está cerca de w = {mejor:.2f}")

def gradiente(x, y, w, b):
    """Derivadas parciales del MSE respecto de w y de b."""
    error = predecir(x, w, b) - y     # predicho menos real
    dw = 2 * np.mean(error * x)
    db = 2 * np.mean(error)
    return dw, db

def perdida_en(w, b):
    """El MSE de la recta (w, b) con los datos de las notas."""
    return mse(nota, predecir(horas, w, b))

# Comprobación con la derivada numérica (diferencia centrada),
# en la recta a ojo
h = 1e-6
w, b = 5.0, 40.0
dw_num = (perdida_en(w + h, b) - perdida_en(w - h, b)) / (2 * h)
db_num = (perdida_en(w, b + h) - perdida_en(w, b - h)) / (2 * h)
print("Con la fórmula:   ", np.round(gradiente(horas, nota, w, b), 4))
print("Derivada numérica:", np.round([dw_num, db_num], 4))

def descenso_gradiente(x, y, tasa, pasos, w=0.0, b=0.0):
    """Ajusta y ≈ w·x + b bajando por el MSE.

    Devuelve los w y b finales y el historial de (w, b, MSE).
    """
    historial = []                    # (w, b, MSE) antes de cada paso
    for paso in range(pasos):
        historial.append((w, b, mse(y, predecir(x, w, b))))
        dw, db = gradiente(x, y, w, b)  # hacia dónde sube el valle
        w = w - tasa * dw             # un paso en contra
        b = b - tasa * db
    historial.append((w, b, mse(y, predecir(x, w, b))))
    return w, b, np.array(historial)

w_final, b_final, historial = descenso_gradiente(
    horas, nota, tasa=0.02, pasos=2000)

for paso in [0, 1, 10, 100, 500, 2000]:
    w_p, b_p, perdida = historial[paso]
    print(f"Paso {paso:4d}:   w = {w_p:5.2f}   b = {b_p:5.2f}"
          f"   MSE = {perdida:7.1f}")

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for eje, paso in zip(ejes.flat, [0, 10, 100, 2000]):
    w_p, b_p, perdida = historial[paso]
    # Solo la recta final va resaltada (color "C1")
    color = "C1" if paso == 2000 else "C2"
    eje.scatter(horas, nota, s=14)
    eje.plot(x_recta, predecir(x_recta, w_p, b_p), color=color)
    # Título en dos líneas para que quepa en cada panel
    eje.set_title(f"Paso {paso}\nMSE = {perdida:.0f}")
for eje in ejes[1]:
    eje.set_xlabel("Horas")
for eje in ejes[:, 0]:
    eje.set_ylabel("Nota")
plt.show()

# Una cuadrícula de 200 × 200 valores de (w, b)
W, B = np.meshgrid(np.linspace(-2, 14, 200), np.linspace(-10, 60, 200))
# Para cada (w, b) de la cuadrícula, el MSE sobre los 30 estudiantes
Z = np.mean((W[..., None] * horas + B[..., None] - nota) ** 2, axis=-1)

niveles = np.geomspace(Z.min() + 1, Z.max(), 14)
plt.contour(W, B, Z, levels=niveles, colors="C2", linewidths=1)
plt.plot(historial[:, 0], historial[:, 1], color="C0",
         label="Camino del descenso")
algunos = [0, 1, 2, 3, 10, 100]       # algunos pasos, marcados
plt.scatter(historial[algunos, 0], historial[algunos, 1],
            color="C0", s=30, zorder=3)
plt.scatter([w_final], [b_final], color="C1", s=160, zorder=4,
            label="Después de 2000 pasos")
plt.xlabel("w (pendiente)")
plt.ylabel("b (intercepto)")
plt.title("El valle del MSE visto desde arriba")
plt.legend(loc="upper right", frameon=True)   # caja blanca sobre las curvas
plt.show()

tasas = [(0.002, "C0", "lenta"),
         (0.02, "C1", "adecuada"),
         (0.04, "C2", "demasiado grande")]

plt.figure()
for tasa, color, nombre in tasas:
    _, _, h = descenso_gradiente(horas, nota, tasa=tasa, pasos=600)
    plt.plot(h[:, 2], color=color, label=f"tasa {tasa}: {nombre}")
    print(f"tasa {tasa}:  MSE después de 600 pasos = {h[-1, 2]:.4g}")

# Escala logarítmica: cada marca del eje es 10 veces la anterior
plt.yscale("log")
plt.ylim(10, 1e5)
plt.xlabel("Paso")
plt.ylabel("MSE (escala logarítmica)")
plt.title("La misma bajada con tres tasas")
plt.legend()
plt.show()

from sklearn.linear_model import LinearRegression

modelo = LinearRegression()
# X tiene que ser una tabla: 30 filas y 1 columna
modelo.fit(horas.reshape(-1, 1), nota)

print(f"scikit-learn:      w = {modelo.coef_[0]:.4f}"
      f"   b = {modelo.intercept_:.4f}")
print(f"Nuestro descenso:  w = {w_final:.4f}   b = {b_final:.4f}")
print("Nota predicha para 6 horas:", modelo.predict([[6]]).round(1))

from sklearn.datasets import load_diabetes
import pandas as pd

diabetes = load_diabetes(as_frame=True)
X = diabetes.data      # tabla: 442 pacientes × 10 características
y = diabetes.target    # avance de la enfermedad un año después

print("Forma de X:", X.shape)
print("Columnas:", list(X.columns))
print("Objetivo: de", y.min(), "a", y.max(),
      "con media", round(y.mean(), 1))

# Estandarizar: cada columna queda con media 0 y desviación estándar 1
X_est = ((X - X.mean()) / X.std(ddof=0)).to_numpy()
y_arr = y.to_numpy()

def descenso_gradiente_varias(X, y, tasa, pasos):
    """Descenso de gradiente para y ≈ X @ w + b (un peso por columna)."""
    n, d = X.shape
    w = np.zeros(d)                   # un peso por característica
    b = 0.0
    for paso in range(pasos):
        error = X @ w + b - y         # los n errores de una vez
        w = w - tasa * 2 * X.T @ error / n
        b = b - tasa * 2 * np.mean(error)
    return w, b

w_gd, b_gd = descenso_gradiente_varias(X_est, y_arr, tasa=0.1,
                                       pasos=20000)
modelo_diabetes = LinearRegression().fit(X_est, y_arr)

comparacion = pd.DataFrame(
    {"descenso": w_gd, "scikit-learn": modelo_diabetes.coef_},
    index=X.columns)
print(comparacion.round(2))
print(f"b:  descenso = {b_gd:.2f}"
      f"   scikit-learn = {modelo_diabetes.intercept_:.2f}")

from sklearn.metrics import r2_score

pred = modelo_diabetes.predict(X_est)
residuos = y_arr - pred
r2 = 1 - np.sum(residuos ** 2) / np.sum((y_arr - y_arr.mean()) ** 2)
rmse = np.sqrt(np.mean(residuos ** 2))
error_abs = np.mean(np.abs(residuos))

print(f"RMSE = {rmse:.1f}   MAE = {error_abs:.1f}   R² = {r2:.3f}")
print(f"R² según scikit-learn: {r2_score(y_arr, pred):.3f}")

# Un modelo con una sola columna: bmi (la de índice 2)
solo_bmi = X_est[:, [2]]
modelo_bmi = LinearRegression().fit(solo_bmi, y_arr)
print(f"R² usando solo bmi: {modelo_bmi.score(solo_bmi, y_arr):.3f}")

plt.scatter(pred, residuos, s=16)
plt.axhline(0, color="C2")
plt.xlabel("Predicción del modelo")
plt.ylabel("Residuo (real menos predicho)")
plt.title("Residuos frente a predicciones")
plt.show()