Curso de ML EN

Capítulo 11 de 37 · intermedio

Regresión lineal

Qué cubre este capítulo

La regresión lineal es el modelo al que recurres primero cuando lo que quieres predecir es un número y no una etiqueta. Ajusta una línea recta — un plano, en realidad, en cuanto tienes más de un feature — a través de tus datos y lee las predicciones directamente de ella. Suena demasiado simple para importar, y justo por eso importa: es el baseline de batalla para todo problema de regresión, el modelo cuyos coeficientes de verdad puedes leer y discutir, y lo que cualquier modelo más sofisticado tiene que superar antes de que valga la pena el esfuerzo.

Lo construimos dos veces. Primero la solución cerrada — la ecuación normal, una sola resolución de sistema lineal que te entrega la mejor línea de un solo golpe. Después descenso de gradiente, que se arrastra hacia la misma respuesta un paso cuesta abajo a la vez. La forma cerrada es la que llevarías a producción; la versión con gradiente es la que querrías observar, porque es el mismo loop de optimización que entrena todo lo que viene después en el curso, y aquí es lo bastante simple para verlo en acción. Corremos ambas contra scikit-learn y confirmamos que las tres aterrizan en la misma línea.

Los datos reales son las viviendas de California: 20,640 grupos de bloques censales, ocho mediciones cada uno — ingreso mediano, antigüedad de las casas, ubicación, algunas más — y como target el valor mediano de la vivienda. También mantendremos a un lado un pequeño conjunto sintético, cuarenta puntos en una dimensión, para que de verdad puedas ver cómo se mueve la línea.

Un poco de historia

Los mínimos cuadrados son una de las ideas más viejas de las matemáticas aplicadas, y vienen con una auténtica pelea por la prioridad. Adrien-Marie Legendre publicó el método en 1805, en un apéndice sobre las órbitas de los cometas, y le dio el nombre que seguimos usando — moindres carrés, mínimos cuadrados. Carl Friedrich Gauss publicó su propia versión en 1809 y afirmó que lo venía usando desde 1795, para predecir la posición del planeta enano Ceres después de que desapareció detrás del sol. Los astrónomos habían perdido a Ceres; Gauss ajustó una curva al puñado de observaciones que tenían, les dijo hacia dónde apuntar sus telescopios, y ahí estaba. Esa es la historia fundacional de todo el campo, y es una regresión.

Gauss fue más lejos que Legendre — conectó los mínimos cuadrados con la distribución normal y demostró que era el estimador óptimo bajo ruido gaussiano, el resultado que hoy llamamos el teorema de Gauss-Markov. Francis Galton le dio el nombre de "regresión" casi un siglo después, en la década de 1880, estudiando cómo las estaturas de los hijos regresaban hacia la media de la población respecto a sus padres. La palabra se quedó aunque describe una peculiaridad de un dataset y no del método. Doscientos años después, la ecuación normal por la que discutieron Gauss y Legendre es la misma que está en el código de abajo, sin cambios.

La intuición

Tienes puntos, y quieres una línea que pase por ellos. No cualquier línea — la que quede tan cerca de todos los puntos como pueda, de modo que cuando llegue una x nueva puedas leer una y sensata. "Tan cerca como pueda" es todo el juego, y los mínimos cuadrados lo hacen preciso: mide la cercanía por la brecha vertical entre cada punto y la línea, eleva esas brechas al cuadrado para que no puedan cancelarse, y elige la línea que hace el total más pequeño.

Elevar al cuadrado es la decisión que define el método, y vale la pena detenerse en ella. Significa que un punto que está al doble de distancia de la línea aporta cuatro veces la penalización, así que el ajuste se dobla con fuerza para evitar los errores grandes y se encoge de hombros ante los pequeños. Eso es lo que hace que la línea se asiente sobre el grueso de los datos — y también es la razón por la que un solo outlier salvaje puede arrastrar toda la línea hacia él. Cada propiedad de la regresión lineal, buena y mala, se remonta a esa brecha al cuadrado.

Aquí está la idea sobre cuarenta puntos sintéticos. La línea es el ajuste de mínimos cuadrados: la única línea, de entre todas las posibles, donde las brechas verticales al cuadrado suman el total más pequeño.

Las matemáticas

Apila tus datos en una matriz XX con nn filas y dd columnas — una fila por ejemplo, una columna por feature. El modelo le da a cada fila una predicción como una suma ponderada de sus features más un intercepto:

y^i=w⊤xi+b,y^=Xw+b\hat{y}_i = w^{\top} x_i + b, \qquad \hat{y} = Xw + b

Aquí w∈Rdw \in \mathbb{R}^{d} es el vector de pesos — un coeficiente por feature — y bb es el intercepto escalar, la predicción cuando todos los features valen cero. Ajustar significa elegir ww y bb para minimizar el error cuadrático medio sobre los nn ejemplos:

L(w,b)=1n∑i=1n(y^i−yi)2L(w, b) = \frac{1}{n} \sum_{i=1}^{n} \left( \hat{y}_i - y_i \right)^2

Esta superficie es un tazón — convexa, un solo mínimo — así que "la mejor línea" está bien definida y no hay mínimos locales donde atorarse. Iguala las derivadas a cero y puedes resolver el mínimo directamente. Mete el intercepto dentro de los pesos pegándole a XX una columna de unos, y la solución es la ecuación normal:

w=(X⊤X)−1X⊤yw = \left( X^{\top} X \right)^{-1} X^{\top} y

Esa es la forma cerrada: una inversa de matriz, una línea de código, exacta. Cuando X⊤XX^{\top}X resulta caro o está mal condicionado no quieres formar la inversa — resuelves el sistema lineal en su lugar, o desciendes por el tazón usando el gradiente. El gradiente de la pérdida es donde el residual r=Xw+b−yr = Xw + b - y sale al frente:

∇wL=2nX⊤r,∂L∂b=2n∑i=1nri\nabla_w L = \frac{2}{n} X^{\top} r, \qquad \frac{\partial L}{\partial b} = \frac{2}{n} \sum_{i=1}^{n} r_i

El descenso de gradiente simplemente camina en contra de él, escalado por una tasa de aprendizaje η\eta, hasta que los pasos dejan de mover algo:

w←w−η ∇wL,b←b−η ∂L∂bw \leftarrow w - \eta \, \nabla_w L, \qquad b \leftarrow b - \eta \, \frac{\partial L}{\partial b}

El mismo destino que la ecuación normal, alcanzado a pie. Haremos ambos y verificaremos que coinciden.

Para qué sirve, para qué no

El argumento a favor de la regresión lineal es la interpretabilidad y la velocidad. Cada coeficiente es una oración: mantén todo lo demás fijo, mueve este feature una unidad, y la predicción se mueve exactamente esta cantidad. Esa es una afirmación que un experto del dominio puede confirmar o rechazar, lo cual vale más que una fracción de punto de accuracy en la mayoría de los lugares donde he trabajado. Se ajusta en forma cerrada, así que entrenar es una sola resolución lineal — sin tuning, sin epochs, sin semilla aleatoria. Y con el modelo así de restringido, hay muy poco que pueda hacer overfitting; es el extremo de alto sesgo y baja varianza del dial, el default seguro.

El argumento en contra es que el mundo muchas veces no es lineal, y la regresión lineal solo puede dibujar una superficie plana. Si la relación verdadera se curva, o dos features solo importan en combinación, un plano no puede capturarlo y ninguna cantidad de datos lo va a arreglar — tendrías que agregar la curva o la interacción como feature tú mismo. También se apoya en supuestos que fallan en silencio: espera que el ruido tenga una dispersión más o menos constante a lo largo del rango (homocedasticidad) y que los errores sean independientes, y cuando eso se rompe los coeficientes igual salen, pero la confianza que les pondrías no debería. Los features correlacionados (multicolinealidad) vuelven inestables los coeficientes — verás los signos ponerse raros en un minuto — y un solo outlier, gracias a esa pérdida al cuadrado, puede hacer palanca sobre toda la línea. Nada de esto lo hace inútil. Lo hace un baseline que entiendes, que es lo más útil que un modelo puede ser.

Los datos

California housing es el primer dataset real estándar de regresión, y viene incluido en scikit-learn, así que no hay nada que descargar. Cada una de las 20,640 filas es un grupo de bloques censales — de unos cientos a unos miles de personas — con ocho features: ingreso mediano, antigüedad mediana de las casas, promedio de cuartos y de recámaras, población, ocupación promedio, y latitud/longitud. El target es el valor mediano de la vivienda en el bloque, en unidades de $100,000. Una peculiaridad que conviene saber de entrada: el target está topado en 5.0, así que un muro de bloques queda exactamente en $500,000, lo cual va a regresar a morder el ajuste en la parte alta.

La señal más fuerte por mucho es el ingreso mediano, lo cual tiene sentido — las zonas más ricas tienen casas más caras. Aquí está el ingreso mediano contra el valor mediano de la vivienda en una muestra de bloques. La relación es claramente linealoide y claramente ruidosa, y puedes ver el tope: ese techo plano de puntos en 5.0.

Constrúyelo, una función a la vez

Seis funciones cortas, NumPy puro, en el orden en que realmente las escribirías: primero la pieza más pequeña que se pueda probar, luego las dos maneras de ajustar.

El modelo es una línea. Las predicciones son los features por los pesos más el intercepto — un producto matriz-vector y una suma:

def predict(X, w, b):
    """The linear model itself: yhat = Xw + b.

    X is (n, d) — n rows, d features. w is (d,), b is a scalar. The result is
    (n,), one predicted number per row. This one line is the whole model;
    everything else in the file is about choosing w and b.
    """
    return X @ w + b

Como eso es un producto de matrices normal y corriente, la misma función predice una fila o un dataset completo, un feature u ocho. Después, lo que estamos tratando de hacer pequeño — el error cuadrático medio, el promedio de los residuales al cuadrado:

def mse_loss(X, y, w, b):
    """Mean squared error: the average squared gap between guess and truth.

    This is the thing we minimize. Squaring makes every miss positive and
    punishes big misses far more than small ones, which is what pins the fit
    to the bulk of the data (and what makes outliers hurt).
    """
    resid = predict(X, w, b) - y
    return float(np.mean(resid ** 2))

La pérdida te dice qué tan equivocada está una línea; el R² te dice si eso sirve de algo, comparando tu error contra el error de simplemente adivinar la media cada vez. Es el número que la gente de verdad reporta, así que lo implementamos directamente:

def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of the variance we explain.

    1 minus (our squared error / the squared error of always guessing the
    mean). 1.0 is perfect, 0.0 is no better than the mean, negative is worse.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot

Ahora el ajuste. La forma cerrada resuelve el mínimo de un solo golpe. Le pegamos una columna de unos a los features para que el intercepto entre como un peso ordinario, y luego resolvemos las ecuaciones normales. Nota que resolvemos el sistema en lugar de invertir literalmente X⊤XX^{\top}X — misma respuesta, mejor comportamiento numérico:

def normal_equation(X, y):
    """The closed-form OLS solution: w = (XtX)^-1 Xt y, intercept and all.

    We glue a column of ones onto X so the intercept rides along as the first
    weight, then solve the normal equations. `np.linalg.solve` factors the
    matrix instead of forming the inverse explicitly — same answer, better
    conditioned. Splits the fitted vector back into (w, b) at the end.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                 # (n, d+1)
    theta = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b

Eso es la regresión lineal, terminada. Todo lo que sigue es la segunda ruta a la misma respuesta, la que podemos observar. El descenso de gradiente necesita la pendiente de la pérdida en el punto actual, para saber hacia dónde es cuesta abajo. Deriva el MSE y el residual sale al frente:

def gradient(X, y, w, b):
    """The slope of the MSE at the current (w, b): which way is downhill.

    Differentiate the mean squared error. The residual falls out front, so
    the gradient is just the residual projected back onto each feature (and
    summed, for the intercept), scaled by 2/n.
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y              # (n,)
    grad_w = (2.0 / n) * (X.T @ resid)        # (d,)
    grad_b = (2.0 / n) * float(resid.sum())
    return grad_w, grad_b

Un paso es simplemente moverse en contra de ese gradiente por la tasa de aprendizaje:

def gd_step(w, b, grad_w, grad_b, lr):
    """One downhill step: move against the gradient by the learning rate."""
    return w - lr * grad_w, b - lr * grad_b

Y el ajuste es empezar desde una línea plana en cero y dar esos pasos hasta que dejan de cambiar algo. Como la pérdida es un tazón, esto no puede atorarse — rueda hasta el mismo fondo al que la ecuación normal salta:

def gradient_descent(X, y, lr, n_iters, record=False):
    """Fit by repeatedly stepping downhill from w = 0, b = 0.

    With a small enough learning rate this walks to the same (w, b) the normal
    equation hands you in one shot — the MSE surface is a bowl, so there's a
    single minimum to fall into. When `record` is on we log (w, b, loss) BEFORE
    each step so the chapter can animate the line converging.
    """
    w = np.zeros(X.shape[1])
    b = 0.0
    history = []
    for _ in range(n_iters):
        if record:
            history.append((w.copy(), b, mse_loss(X, y, w, b)))
        grad_w, grad_b = gradient(X, y, w, b)
        w, b = gd_step(w, b, grad_w, grad_b, lr)
    if record:
        history.append((w.copy(), b, mse_loss(X, y, w, b)))
    return w, b, history

Míralo funcionar

Esto es descenso de gradiente ajustando la línea sintética, un paso real por frame. El panel de arriba son los datos con la línea actual dibujada encima; la línea empieza plana en cero — pendiente 0, intercepto 0 — y rota y se eleva hasta acomodarse conforme corren los pasos. El panel de abajo traza el MSE en un eje logarítmico, para que veas la pérdida caer por un precipicio al principio y luego aplanarse conforme la línea se asienta. La leyenda indica la iteración, la pendiente y el intercepto actuales, y el MSE actual.

Presiona play. Cada frame es una llamada al loop de gradiente-y-paso de impl.py, con tasa de aprendizaje 0.12. Reinicia y córrelo cuantas veces quieras.

Dos cosas para notar. La línea se mueve más rápido al principio, cuando los residuales son más grandes y el gradiente es empinado, y apenas se mueve al final — eso es la curva de pérdida aplanándose, el modelo arrastrándose el último tramo hasta el mínimo. Y el lugar donde aterriza no es la línea con la que generamos los datos. Construimos los puntos con pendiente 2.6 e intercepto 5.0, pero el ajuste de mínimos cuadrados sale con pendiente 3.00 e intercepto 5.10, porque cuarenta puntos ruidosos no fijan la línea con exactitud. Después de 60 iteraciones el descenso de gradiente ha igualado la solución cerrada a dos decimales; ambos coinciden en que la mejor línea a través de estos puntos tiene pendiente 3.00, no el 2.6 del que partimos. Esa brecha entre la línea generadora verdadera y el mejor ajuste a una muestra finita no es un bug — es de lo que trata la estadística.

La implementación completa

El archivo entero, sin librería, de arriba a abajo — el modelo, las dos métricas, la ecuación normal y el descenso de gradiente. Este es el código que la animación de arriba realmente ejecutó:

"""Ordinary least squares linear regression, built from scratch.

A linear model predicts a number as a weighted sum of the features plus an
intercept. "Fitting" means choosing the weights that make the mean squared
error smallest. There are two ways to get there and this file has both: the
closed-form normal equation (one linear solve, exact) and gradient descent
(take a step downhill, repeat) so we can watch the fit converge.

Pure NumPy — no ML library anywhere in this file. Every function below shows
up in the chapter one step at a time; the `# region:` markers are what the
book's include directives pull in.
"""

import numpy as np
import pandas as pd


# region: predict
def predict(X, w, b):
    """The linear model itself: yhat = Xw + b.

    X is (n, d) — n rows, d features. w is (d,), b is a scalar. The result is
    (n,), one predicted number per row. This one line is the whole model;
    everything else in the file is about choosing w and b.
    """
    return X @ w + b
# endregion


# region: mse_loss
def mse_loss(X, y, w, b):
    """Mean squared error: the average squared gap between guess and truth.

    This is the thing we minimize. Squaring makes every miss positive and
    punishes big misses far more than small ones, which is what pins the fit
    to the bulk of the data (and what makes outliers hurt).
    """
    resid = predict(X, w, b) - y
    return float(np.mean(resid ** 2))
# endregion


# region: r2
def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of the variance we explain.

    1 minus (our squared error / the squared error of always guessing the
    mean). 1.0 is perfect, 0.0 is no better than the mean, negative is worse.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot
# endregion


# region: normal_equation
def normal_equation(X, y):
    """The closed-form OLS solution: w = (XtX)^-1 Xt y, intercept and all.

    We glue a column of ones onto X so the intercept rides along as the first
    weight, then solve the normal equations. `np.linalg.solve` factors the
    matrix instead of forming the inverse explicitly — same answer, better
    conditioned. Splits the fitted vector back into (w, b) at the end.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                 # (n, d+1)
    theta = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b
# endregion


# region: gradient
def gradient(X, y, w, b):
    """The slope of the MSE at the current (w, b): which way is downhill.

    Differentiate the mean squared error. The residual falls out front, so
    the gradient is just the residual projected back onto each feature (and
    summed, for the intercept), scaled by 2/n.
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y              # (n,)
    grad_w = (2.0 / n) * (X.T @ resid)        # (d,)
    grad_b = (2.0 / n) * float(resid.sum())
    return grad_w, grad_b
# endregion


# region: gd_step
def gd_step(w, b, grad_w, grad_b, lr):
    """One downhill step: move against the gradient by the learning rate."""
    return w - lr * grad_w, b - lr * grad_b
# endregion


# region: gradient_descent
def gradient_descent(X, y, lr, n_iters, record=False):
    """Fit by repeatedly stepping downhill from w = 0, b = 0.

    With a small enough learning rate this walks to the same (w, b) the normal
    equation hands you in one shot — the MSE surface is a bowl, so there's a
    single minimum to fall into. When `record` is on we log (w, b, loss) BEFORE
    each step so the chapter can animate the line converging.
    """
    w = np.zeros(X.shape[1])
    b = 0.0
    history = []
    for _ in range(n_iters):
        if record:
            history.append((w.copy(), b, mse_loss(X, y, w, b)))
        grad_w, grad_b = gradient(X, y, w, b)
        w, b = gd_step(w, b, grad_w, grad_b, lr)
    if record:
        history.append((w.copy(), b, mse_loss(X, y, w, b)))
    return w, b, history
# endregion


def load_data(path="../data/california_housing.csv"):
    """California housing: 20,640 block groups, 8 features, median house value.

    Returns (X, y, feature_names). y is the median house value in units of
    $100,000, the column the model learns to predict.
    """
    df = pd.read_csv(path)
    target = "MedHouseVal"
    features = [c for c in df.columns if c != target]
    X = df[features].to_numpy(float)
    y = df[target].to_numpy(float)
    return X, y, features

La versión de librería

Nadie escribe a mano la ecuación normal en producción, y una vez que la entiendes tú tampoco deberías. sklearn.linear_model.LinearRegression es el mismo modelo con el mismo objetivo — mínimos cuadrados ordinarios — y ajusta el intercepto por default:

def sklearn_fit(X_train, y_train):
    """Fit OLS with sklearn. Returns (weights, intercept) — the twin of our
    normal_equation, so the coefficients line up one for one."""
    model = LinearRegression()            # fit_intercept=True by default
    model.fit(X_train, y_train)
    return model.coef_, float(model.intercept_)

Por debajo, sklearn no forma (X⊤X)−1(X^{\top}X)^{-1} más de lo que lo hacemos nosotros. Resuelve el sistema de mínimos cuadrados con una SVD, que es más robusta que la ecuación normal cuando los features son casi colineales — y California housing tiene algo de eso, con el promedio de cuartos y el promedio de recámaras moviéndose juntos y la latitud y longitud trazando la costa. Con datos bien condicionados las dos rutas dan coeficientes idénticos; con datos feos la SVD se degrada con más gracia. Para la evaluación lo dejamos ajustar y reportamos R² y MSE de prueba:

def sklearn_score(X_train, y_train, X_test, y_test):
    """Fit on train, report test R^2 and test MSE — the face-off numbers."""
    model = LinearRegression()
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    r2 = float(model.score(X_test, y_test))
    mse = float(((pred - y_test) ** 2).mean())
    return r2, mse

La otra cosa que vale la pena leer del modelo ajustado son los coeficientes, porque aquí es donde la regresión lineal se gana el pan. El ingreso mediano sale alrededor de 0.44 — una unidad extra de ingreso mediano (eso es $10,000) predice unos $43,900 más de valor mediano de vivienda, con todo lo demás fijo. La latitud y la longitud salen ambas fuertemente negativas, alrededor de -0.42 y -0.44, que es el modelo descubriendo el mapa: los precios caen conforme te mueves al norte y tierra adentro desde la costosa costa sur. Y el promedio de recámaras y el promedio de cuartos salen con signos opuestos (+0.65 y -0.11) aunque ambos son "más casa" — eso es multicolinealidad, los dos features correlacionados repartiéndose el crédito de una manera que hace difícil leer cada coeficiente por su cuenta. La interpretabilidad es el punto de venta, pero viene con letras chiquitas.

Desde cero contra la librería

Divide los datos 80/20 — 16,512 bloques para entrenar, 4,128 apartados — ajusta sobre la mitad de entrenamiento de tres maneras, y evalúa sobre la mitad apartada. Nuestra ecuación normal, nuestro descenso de gradiente (sobre features estandarizados, ya que las escalas crudas difieren por órdenes de magnitud), y sklearn:

Las tres barras miden lo mismo: R² de prueba de 0.585, MSE de prueba de 0.528. No es que estén cerca, son idénticas — nuestra ecuación normal iguala los coeficientes de sklearn a doce decimales, y el descenso de gradiente, con suficientes pasos, aterriza en las mismas predicciones. Solo hay una línea de mínimos cuadrados, y todo método correcto la encuentra. Las cifras exactas viven en results.json, que se regenera cada vez que cambia el código, para que la prosa y la gráfica no puedan desviarse de lo que el código produjo.

Un R² de 0.585 significa que el modelo explica alrededor del 59% de la varianza en los valores de las viviendas a partir de estos ocho features — decente para un ajuste lineal sacado de la caja, y lejísimos del techo. Para ver por dónde se le fuga, grafica los residuales contra los valores ajustados. En un ajuste limpio esto es una nube sin forma centrada en cero; que haya estructura aquí significa que al modelo se le está escapando algo:

La nube no es informe, y esa es la lección. Hay una raya diagonal de puntos bajando por el lado derecho — esos son los bloques topados, los que en realidad valen $500,000 o más y que el modelo, sin tener idea del tope, predice tan alto como $717,000 para luego comerse un residual negativo enorme. Y la dispersión de los residuales se ensancha conforme crece el valor ajustado, una forma de abanico: eso es heterocedasticidad, el supuesto de varianza constante fallando a plena vista. El modelo está más equivocado, y equivocado de forma menos predecible, en las casas caras que en las baratas. La regresión lineal te entregará coeficientes con gusto sobre datos como estos; la gráfica de residuales es donde te dice lo que no pudo ajustar.

Conclusiones

Ajusta la línea primero. El mismo consejo que con el umbral, un peldaño más arriba: antes del gradient boosting y la red neuronal, corre mínimos cuadrados ordinarios, lee el R², y mira los coeficientes. Ese número es la vara que tu modelo más sofisticado tiene que superar, y esos coeficientes son una verificación de cordura gratis sobre si los features siquiera apuntan en la dirección correcta. Si el ingreso sale negativo, tienes un bug o un problema de datos, y quieres enterarte por un modelo que puedes leer, no por una caja negra.

Las dos construcciones de este capítulo son las dos mitades del resto del curso. La ecuación normal es la última vez que obtendrás la respuesta en forma cerrada — es un regalo que la pérdida de la regresión lineal sea un tazón convexo con una fórmula para el fondo. Todo lo más difícil renuncia a eso, y cuando lo hace, el descenso de gradiente es la herramienta que sigue funcionando: el loop exacto que viste ajustar esta línea es el que entrena la regresión logística, y el boosting, y toda red neuronal más adelante, solo que sobre una superficie de pérdida sin fórmula y, eventualmente, sin garantía de un mínimo único. Apréndelo aquí, en el único problema donde puedes verificar su trabajo contra una respuesta exacta.

Y quédate con la gráfica de residuales. Los coeficientes te dicen lo que el modelo aprendió; los residuales te dicen lo que no pudo. El abanico y la raya diagonal de esa última gráfica son el modelo siendo honesto sobre sus propios supuestos — linealidad y varianza constante — justo donde se rompen. Recurre a algo más pesado cuando los residuales muestren estructura que la línea no puede absorber, lo cual con datos reales es la mayoría de las veces. Pero empieza aquí, porque un ajuste lineal que entiendes vale más que uno complicado que no.