8.2 Regresión por Mínimos cuadrados con Python

Referencia: Chapra 17.1.2 p 469. Burden 8.1 p499

Criterio de ajuste a una línea recta por mínimos cuadrados

Una aproximación de la relación entre los puntos xi, yi por medio de un polinomio de grado 1, busca minimizar la suma de los errores residuales de los datos.

y_{i,modelo} = p_1(x) = a_0 + a_1 x_i

Se busca el valor mínimo para los cuadrados de las diferencias entre los valores de yi con la línea recta.

S_r = \sum_{i=1}^{n} \Big( y_{i,medida} - y_{i,modelo} \Big)^2 S_r= \sum_{i=1}^{n} \Big( y_{i} - (a_0 + a_1 x_i) \Big)^2

para que el error acumulado sea mínimo, se deriva respecto a los coeficientes de la recta a0 y a1 y se iguala a cero,

\frac{\partial S_r}{\partial a_0}= (-1)2 \sum_{i=1}^{n} \Big( y_{i} - a_0 - a_1 x_i \Big) \frac{\partial S_r}{\partial a_1}= (-1)2 \sum_{i=1}^{n} \Big( y_{i} - a_0 - a_1 x_i \Big)x_i 0 = \sum_{i=1}^{n} y_i - \sum_{i=1}^{n} a_0 - \sum_{i=1}^{n} a_1 x_i 0= \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} a_0x_i - \sum_{i=1}^{n} a_1 x_i^2

simplificando,

\sum_{i=1}^{n} a_0 + a_1 \sum_{i=1}^{n} x_i = \sum_{i=1}^{n} y_{i} a_0 \sum_{i=1}^{n} x_i + a_1 \sum_{i=1}^{n} x_i^2 = \sum_{i=1}^{n} y_i x_i

que es un conjunto de dos ecuaciones lineales simultaneas con dos incógnitas a0 y a1, los coeficientes del sistema de ecuaciones son las sumatorias que se obtienen completando la siguiente tabla:

Tabla de datos y columnas para ∑ en la fórmula
xi yi xiyi xi2 yi2
x0 y0 x0y0 x02 y02
∑ xi ∑ yi ∑ xiyi ∑ xi2 ∑ yi2
\begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix}

cuya solución es:

a_1 = \frac{n \sum x_i y_i - \sum x_i \sum y_i}{n \sum x_i^2 - \big( \sum x_i \big) ^2 } a_0 = \overline{y} - a_1 \overline{x}

usando la media de los valores en cada eje para encontrar a0

Coeficiente de correlación

El coeficiente de correlación se puede obtener con las sumatorias anteriores usando la siguiente expresión:

r= \frac{n \sum x_i y_i - \big( \sum x_i \big) \big( \sum y_i\big)} {\sqrt{n \sum x_i^2 -\big(\sum x_i \big) ^2 }\sqrt{n \sum y_i^2 - \big( \sum y_i \big)^2}}

En un ajuste perfecto, Sr = 0 y r = r2 = 1, significa que la línea explica
el 100% de la variabilidad de los datos.

Aunque el coeficiente de correlación ofrece una manera fácil de medir la bondad del ajuste, se deberá tener cuidado de no darle más significado del que ya tiene.

El solo hecho de que r sea “cercana” a 1 no necesariamente significa que el ajuste sea “bueno”. Por ejemplo, es posible obtener un valor relativamente alto de r cuando la relación entre y y x no es lineal.

Los resultados indicarán que el modelo lineal explicó r2 % de la incertidumbre original


Algoritmo en Python

Siguiendo con los datos propuestos del ejemplo en Chapra 17.1 p470:

xi = [1, 2, 3, 4, 5, 6, 7] 
yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

Aplicando las ecuaciones para a0 y a1 se tiene el siguiente resultado para los datos de prueba:

 f =  0.839285714285714*x + 0.0714285714285712
coef_correlación   r  =  0.9318356132188194
coef_determinación r2 =  0.8683176100628931
86.83% de los datos está descrito en el modelo lineal
>>>

con las instrucciones:

# mínimos cuadrados, regresión con polinomio grado 1
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xi = [1,   2,   3,  4,  5,   6, 7]
yi = [0.5, 2.5, 2., 4., 3.5, 6, 5.5]

# PROCEDIMIENTO
xi = np.array(xi,dtype=float)
yi = np.array(yi,dtype=float)
n  = len(xi)

# sumatorias y medias
xm  = np.mean(xi)
ym  = np.mean(yi)
sx  = np.sum(xi)
sy  = np.sum(yi)
sxy = np.sum(xi*yi)
sx2 = np.sum(xi**2)
sy2 = np.sum(yi**2)

# coeficientes a0 y a1
a1 = (n*sxy-sx*sy)/(n*sx2-sx**2)
a0 = ym - a1*xm

# polinomio grado 1
x = sym.Symbol('x')
f = a0 + a1*x

fx = sym.lambdify(x,f)
fi = fx(xi)

# coeficiente de correlación
numerador = n*sxy - sx*sy
raiz1 = np.sqrt(n*sx2-sx**2)
raiz2 = np.sqrt(n*sy2-sy**2)
r = numerador/(raiz1*raiz2)

# coeficiente de determinacion
r2 = r**2
r2_porcentaje = np.around(r2*100,2)

# SALIDA
# print('ymedia =',ym)
print(' f = ',f)
print('coef_correlación   r  = ', r)
print('coef_determinación r2 = ', r2)
print(str(r2_porcentaje)+'% de los datos')
print('     está descrito en el modelo lineal')

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
# plt.stem(xi,yi,bottom=ym,linefmt ='--')
plt.plot(xi,fi, color='orange',  label=f)

# lineas de error
for i in range(0,n,1):
    y0 = np.min([yi[i],fi[i]])
    y1 = np.max([yi[i],fi[i]])
    plt.vlines(xi[i],y0,y1, color='red',
               linestyle ='dotted')
plt.legend()
plt.xlabel('xi')
plt.title('minimos cuadrados')
plt.show()

Coeficiente de correlación con Numpy

Tambien es posible usar la librería numpy para obtener el resultado anterior,

>>> coeficientes = np.corrcoef(xi,yi)
>>> coeficientes
array([[1.        , 0.93183561],
       [0.93183561, 1.        ]])
>>> r = coeficientes[0,1]
>>> r
0.9318356132188195