Categoría: Unidad 8 Aproximación por Regresión

De una estación meteorológica se obtiene un archivo.csv con los datos de los sensores disponibles durante una semana.

2021OctubreEstMetorologica.csv

Para analizar el comportamiento de la variable de temperatura, se requiere disponer de un modelo polinomial que describa la temperatura a lo largo del día, para un solo día.

Como valores de variable independiente utilice un equivalente numérico decimal de la hora, minuto y segundo.

a. Realice la lectura del archivo, puede usar las instrucciones descritas en el  enlace: Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad (CCPG1001: Fundamentos de programación)

b. Seleccione los datos del día 1 del mes para realizar la gráfica, y convierta tiempo en equivalente decimal.

c. Mediante pruebas, determine el grado del polinomio más apropiado para minimizar los errores.

Desarrollo

literal a

Se empieza con las instrucciones del enlace dadas añadiendo los parpametros de entrada como undia = 0 y grado del polinomio como gradom = 2 como el ejercicio anterior.

literal b

En el modelo polinomial, el eje x es un número decimal, por lo que los valores de hora:minuto:segundo se convierte a un valor decimal. Para el valor decimal de xi se usa la unidad de horas y las fracciones proporcionales de cada minuto y segundo.

literal c

Se inicia con el valor de gradom = 2, observando el resultado se puede ir incrementando el grado del polinomio observando los parámetros de error y coeficiente de determinación hasta cumplir con la tolerancia esperada segun la aplicación del ejercicio.

Resultados obtenidos:

columnas en tabla: 
Index(['No', 'Date', 'Time', 'ColdJunc0', 'PowerVolt', 'PowerKind', 'WS(ave)',
       'WD(ave)', 'WS(max)', 'WD(most)', 'WS(inst_m)', 'WD(inst_m)',
       'Max_time', 'Solar_rad', 'TEMP', 'Humidity', 'Rainfall', 'Bar_press.',
       'Soil_temp', 'fecha', 'horadec'],
      dtype='object')
ymedia =  25.036805555555553
 f = -6.57404678141012e-6*x**6 + 0.00052869201494877*x**5 - 0.0152875582352169*x**4 + 0.184200388015364*x**3 - 0.761164009032398*x**2 + 0.393278389794015*x + 22.1142936414255
coef_determinacion r2 =  0.9860621424061621
98.61% de los datos
     se describe con el modelo

Instrucciones en Python

# lecturas archivo.csv de estación meteorológica
import numpy as np
import sympy as sym
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator

# INGRESO
narchivo = "2021OctubreEstMetorologica.csv"
sensor = 'TEMP'
undia  = 2 # dia a revisar
gradom = 6 # grado del polinomio

# PROCEDIMIENTO
tabla = pd.read_csv(narchivo, sep=';',decimal=',')
n_tabla = len(tabla)

# fechas concatenando columnas de texto
tabla['fecha'] = tabla['Date']+' '+tabla['Time']

# convierte a datetime
fechaformato = "%d/%m/%Y %H:%M:%S"
tabla['fecha'] = pd.to_datetime(tabla['fecha'],
                                format=fechaformato)
# serie por dia, busca indices
diaIndice = [0] # indice inicial
for i in range(1,n_tabla-1,1):
    i0_fecha = tabla['fecha'][i-1]
    i1_fecha = tabla['fecha'][i]
    if i0_fecha.day != i1_fecha.day:
        diaIndice.append(i)
diaIndice.append(len(tabla)-1) # indice final
m = len(diaIndice)

# horas decimales en un día
horadia = tabla['fecha'].dt.hour
horamin = tabla['fecha'].dt.minute
horaseg = tabla['fecha'].dt.second
tabla['horadec'] = horadia + horamin/60 + horaseg/3600

# PROCEDIMIENTO Regresión Polinomial grado m
m = gradom
# Selecciona dia
i0 = diaIndice[undia]
i1 = diaIndice[undia+1]
# valores a usar en regresión
xi = tabla['horadec'][i0:i1]
yi = tabla[sensor][i0:i1]
n  = len(xi)

# llenar matriz a y vector B
k = m + 1
A = np.zeros(shape=(k,k),dtype=float)
B = np.zeros(k,dtype=float)
for i in range(0,k,1):
    for j in range(0,i+1,1):
        coeficiente = np.sum(xi**(i+j))
        A[i,j] = coeficiente
        A[j,i] = coeficiente
    B[i] = np.sum(yi*(xi**i))

# coeficientes, resuelve sistema de ecuaciones
C = np.linalg.solve(A,B)

# polinomio
x = sym.Symbol('x')
f = 0
fetiq = 0
for i in range(0,k,1):
    f = f + C[i]*(x**i)
    fetiq = fetiq + np.around(C[i],4)*(x**i)

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

# errores
ym = np.mean(yi)
xm = np.mean(xi)
errado = fi - yi

sr = np.sum((yi-fi)**2)
syx = np.sqrt(sr/(n-(m+1)))
st = np.sum((yi-ym)**2)

# coeficiente de determinacion
r2 = (st-sr)/st
r2_porcentaje = np.around(r2*100,2)

# SALIDA
print(' columnas en tabla: ')
print(tabla.keys())
print('ymedia = ',ym)
print(' f =',f)
print('coef_determinacion r2 = ',r2)
print(str(r2_porcentaje)+'% de los datos')
print('     se describe con el modelo')

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xi,fi, color='orange', label=fetiq)

# lineas de error
for i in range(i0,i1,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.xlabel('xi - hora en decimal')
plt.ylabel('yi - '+ sensor)
plt.legend()
etiq_titulo = sensor+ ' dia '+str(undia)
plt.title(etiq_titulo+': Regresion polinomial grado '+str(m))
plt.show()

Tarea

Determinar el polinomio de regresión para los días 3 y 5, y repetir el proceso para el sensor de Humedad (‘Humidity’)

Referencia: Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad en Fundamentos de Programación

Referencia: Chapra 17.2 p482, Burden 8.1 p501

Una alternativa a usar transformaciones ente los ejes para ajustar las curvas es usar regresión polinomial extendiendo el procedimiento de los mínimos cuadrados.

Suponga que la curva se ajusta a un polinomio de segundo grado o cuadrático

y = a_0 + a_1 x + a_2 x^2 +e

usando nuevamente la suma de los cuadrados de los residuos, se tiene,

S_r = \sum_{i=1}^n (y_i- (a_0 + a_1 x_i + a_2 x_i^2))^2

para minimizar los errores se deriva respecto a cada coeficiente: a₀, a₁, a₂

\frac{\partial S_r}{\partial a_0} = -2\sum (y_i- a_0 - a_1 x_i - a_2 x_i^2) \frac{\partial S_r}{\partial a_1} = -2\sum x_i (y_i- a_0 - a_1 x_i - a_2 x_i^2) \frac{\partial S_r}{\partial a_2} = -2\sum x_i^2 (y_i- a_0 - a_1 x_i - a_2 x_i^2)

Se busca el mínimo de cada sumatoria, igualando a cero la derivada y reordenando, se tiene el conjunto de ecuaciones:

a_0 (n) + a_1 \sum x_i + a_2 \sum x_i^2 = \sum y_i a_0 \sum x_i + a_1 \sum x_i^2 + a_2 \sum x_i^3 =\sum x_i y_i a_0 \sum x_i^2 + a_1 \sum x_i^3 + a_2 \sum x_i^4 =\sum x_i^2 y_i

con incógnitas a₀, a₁ y a₂, cuyos coeficientes se pueden evaluar con los puntos observados. Se puede usar un médoto directo de la unidad 3 para resolver el sistema de ecuaciones Ax=B.

\begin{bmatrix} n & \sum x_i & \sum x_i^2 \\ \sum x_i & \sum x_i^2 & \sum x_i^3 \\ \sum x_i^2 & \sum x_i^3 & \sum x_i^4 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \\ \sum x_i^2 y_i \end{bmatrix}

C = np.linalg.solve(A,B)

y = C_0 + C_1 x + C_2 x^2

El error estándar se obtiene como:

S_{y/x} = \sqrt{\frac{S_r}{n-(m+1)}}

siendo m el grado del polinomio usado, para el caso presentado m = 2

S_r = \sum{(yi-fi)^2}

S_r es la suma de los cuadrados de los residuos alrededor de la línea de regresión.

xi	yi	(yi-ymedia)2	(yi-fi)2
…	…
…	…
	∑yi	St	Sr

r^2 = \frac{S_t-S_r}{S_t} S_t = \sum{(yi-ym)^2}

siendo S_t la suma total de los cuadrados alrededor de la media para la variable dependiente y. Este valor es la magnitud del error residual asociado con la variable dependiente antes de la regresión.

El algoritmo se puede desarrollar de forma semejante a la presentada en la sección anterior,

Ejercicio Chapra 17.5 p484

Los datos de ejemplo para la referencia son:

xi = [0,   1,    2,    3,    4,   5]
yi = [2.1, 7.7, 13.6, 27.2, 40.9, 61.1]

resultado de algoritmo:

fx  =  1.86071428571429*x**2 + 2.35928571428571*x + 2.47857142857143
Syx =  1.117522770621316
coef_determinacion r2 =  0.9985093572984048
99.85% de los datos se describe con el modelo
>>> 

---

Instrucciones en Python

# regresion con polinomio grado m=2
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xi = [0,   1,    2,    3,    4,   5]
yi = [2.1, 7.7, 13.6, 27.2, 40.9, 61.1]
m  = 2

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

# llenar matriz a y vector B
k = m + 1
A = np.zeros(shape=(k,k),dtype=float)
B = np.zeros(k,dtype=float)
for i in range(0,k,1):
    for j in range(0,i+1,1):
        coeficiente = np.sum(xi**(i+j))
        A[i,j] = coeficiente
        A[j,i] = coeficiente
    B[i] = np.sum(yi*(xi**i))

# coeficientes, resuelve sistema de ecuaciones
C = np.linalg.solve(A,B)

# polinomio
x = sym.Symbol('x')
f = 0
fetiq = 0
for i in range(0,k,1):
    f = f + C[i]*(x**i)
    fetiq = fetiq + np.around(C[i],4)*(x**i)

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

# errores
ym = np.mean(yi)
xm = np.mean(xi)
errado = fi - yi

sr = np.sum((yi-fi)**2)
syx = np.sqrt(sr/(n-(m+1)))
st = np.sum((yi-ym)**2)

# coeficiente de determinacion
r2 = (st-sr)/st
r2_porcentaje = np.around(r2*100,2)

# SALIDA
print('ymedia = ',ym)
print(' f =',f)
print('coef_determinacion r2 = ',r2)
print(str(r2_porcentaje)+'% de los datos')
print('     se describe con el modelo')

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

# 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.xlabel('xi')
plt.ylabel('yi')
plt.legend()
plt.title('Regresion polinomial grado '+str(m))
plt.show()

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 x_i, y_i 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 y_i 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 a₀ y a₁ 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 a₀ y a₁, 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 a₀

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, S_r = 0 y r = r² = 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ó r² % 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 a₀ y a₁ 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

Referencia: Chapra 17.1 p 466. Burden 8.1 p498

Cuando los datos de un experimento presentan variaciones o errores sustanciales respecto al modelo matemático, la interpolación polinomial presentada en la Unidad 4 es inapropiada para predecir valores intermedios.

En el ejemplo de Chapra 17.1 p470, se presentan los datos de un experimento mostados en la imagen y la siguiente tabla:

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

Un polinomio de interpolación, por ejemplo de Lagrange de grado 6 pasará por todos los puntos, pero oscilando.

Una función de aproximación que se ajuste a la tendencia general, que no necesariamente pasa por los puntos de muestra puede ser una mejor respuesta. Se busca una «curva» que minimice las diferencias entre los puntos y la curva, llamada regresión por mínimos cuadrados.

---

Descripción con la media yi

Considere una aproximación para la relación entre los puntos x_i, y_i como un polinomio, grado 0 que sería la media de y_i. Para este caso, los errores se presentan en la gráfica:

Otra forma sería aproximar el comportamiento de los datos es usar un polinomio de grado 1. En la gráfica se pueden observar que para los mismos puntos el error disminuye considerablemente.

El polinomio de grado 1 recibe el nombre de regresión por mínimos cuadrados, que se desarrolla en la siguiente sección.

Instrucciones Python

# representación con la media
import numpy as np
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)

xm = np.mean(xi)
ym = np.mean(yi)

# SALIDA
print('ymedia = ', ym)

# grafica
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.stem(xi,yi,bottom=ym, linefmt = '--')
plt.xlabel('xi')
plt.ylabel('yi')
plt.legend()
plt.show()

---

Instrucciones Python – compara interpolación y regresión.

Para ilustrar el asunto y para comparar los resultados se usa Python, tanto para interpolación y mínimos cuadrados usando las librerías disponibles. Luego se desarrolla el algoritmo paso a paso.

# mínimos cuadrados
import numpy as np
import scipy.interpolate as sci
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)
yi = np.array(yi)
n = len(xi)

# polinomio Lagrange
px = sci.lagrange(xi,yi)
xj = np.linspace(min(xi),max(xi),100)
pj = px(xj)

# mínimos cuadrados
A = np.vstack([xi, np.ones(n)]).T
[m0, b0] = np.linalg.lstsq(A, yi, rcond=None)[0]
fx = lambda x: m0*(x)+b0
fi = fx(xi)

# ajusta límites
ymin = np.min([np.min(pj),np.min(fi)])
ymax = np.max([np.max(pj),np.max(fi)])

# SALIDA
plt.subplot(121)
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xj,pj,label='P(x) Lagrange')
plt.ylim(ymin,ymax)
plt.xlabel('xi')
plt.ylabel('yi')
plt.title('Interpolación Lagrange')

plt.subplot(122)
plt.plot(xi,yi,'o',label='(xi,yi)')
plt.plot(xi,fi,label='f(x)')
plt.ylim(ymin,ymax)
plt.xlabel('xi')
plt.title('Mínimos cuadrados')

plt.show()