Funciones de medición de error de predicción — Parte II —

En una entrada anterior hablé sobre las funciones de error para medir la capacidad de cualquier modelo mecánico de material para estimar datos experimentales. En esta entrada daré continuación a esto último, pero de manera práctica, es decir, al presentar la implementación en el lenguaje de programación Python 3.8.

Se ha mencionado que las ecuaciones más comunes para medir el error son las siguientes:

• Error Cuadrático Medio (MSE),
• $MSE = \frac{1}{n}\sum_{i=1}^{n}(\sigma_{m_i} - \sigma_{e_i})^{2}$

• Raíz del Error Cuadrático Medio (RMSE), 
• $RMSE = \sqrt{(\frac{1}{n})\sum_{i=1}^{n}(\sigma_{m_i} - \sigma_{e_i})^{2}}$
• Coeficiente de Determinación, 
• $R^2=1-{\sum_{i=1}^{n}(\sigma_{m_i}-\sigma_{e_i})^2}/{\sum_{i=1}^{n}(\sigma_{e_i}-\bar{\sigma_{m}})^2}$.

donde $\sigma_{m}$ y $\sigma_{e}$ refieren a los esfuerzos dados por el modelo y los medidos experimentalmente de manera respectiva. Asimismo,  $\bar{\sigma_{m}}$ representa la media de los valores contenidos en las predicciones dadas por un modelo.

Algo que no se mencionó, sin embargo, es qué miden las ecuaciones anteriores. Por esto, la siguiente lista de ítem enumera su significado, y propósito de manera sucinta:

1. El $MSE$ mide el promedio de la desviación cuadrada entre los valores ajustados con la observación de los datos reales. Los posibles valores que puede adoptar esta medida de error son $ 0 \leq MSE < \infty $. Los valores cercanos a cero indican un buen modelo (para el caso de modelos estacionarios).
2. Estadísticamente, por otra parte,  la $RMSE$ es la raíz cuadrada de la varianza de los residuos o, elementalmente, la raíz cuadrada del $MSE$. Los valores de esta medida pueden ser $ 0 \leq RMSE < \infty $. En este caso, también se desea que los valores sean cercanos a cero.
3. El valor del coeficiente de determinación R2 mide la cantidad de variación que representa el modelo ajustado. Se utiliza con frecuencia para comparar modelos y evaluar qué modelo proporciona el mejor ajuste a los datos. Los valores del coeficiente de determinación pueden ser $ -\infty < R^2 \leq 1$. En el caso de esta medida de error, se espera que mientras mas cercano sea a la unidad, se tiene un mejor modelo que estima datos experimentales, y los valores negativos indican modelos pobres en capacidad de estimación.

Antes de seguir con la implementación, debo hacer énfasis en que estas medidas de error no son las únicas, y en la estádistica existen muchas otras, como el $AIC$, $\text{Adj } r^2$, $BIC$, o $PIC$.

Implementación de las medidas de error en Python 3.8.

Para la implementación de las medidas de error que ya se han mencionado, usaré Python 3.8, Pandas, y scikit learn. Este último sólo para cotejar las funciones que aquí se presentan para obtener valores del error de un caso específico de material.

Para dar estructura a esto, dividiré las siguientes secciones en: (1) las funciones en Python, (2) los datos para validarlas, (3) el experimento con pandas, y (4) la validación con scikit learn.

(1) Las funciones en Python

Al observar nuestras ecuaciones, nos percataremos que existen sumas en cada una de ellas. Para llevar al código de una función, lo que podemos hacer es iterar en un ciclo for cada suma para cada una de las ecuaciones y luego operar sobre ellas. 

Comencemos con la función del Error Cuadrático Medio:

1 2 3 4 5 6

# MSE def mse(y_m, y_e): sum1 = 0.0 for i in range(len(y_m)): sum1 += ((y_m[i] - y_e[i]) ** 2) return sum1 / len(y_m)

En la primera línea lo que se hace es crear la función mse, que tiene por parámetros a dos variables. Para nuestro caso, estas vendrían a ser estructuras de datos tipo vector (aunque en el código no se especifica), y_m y y_e, que refieren a datos del modelo y observados en experimentación.

Una vez que se crea la función, incializamos la varible local sum, que contendrá el valor al hacer la suma de la ecuación del $MSE$. Consecuentemente, se hace el ciclo para llenarla con el cuadrado de la resta de los datos individuales para cada iteración i tal como se muestra en el código anterior.

Finalmente, la función regresa el valor de la suma dividida entre el número total de datos del vector y_m. Esto se logra mediante la función len() en Python, que, dada una estructura como los vectores, regresa a su vez el tamaño de ésta.

Dado que la $RMSE$ es simplemente la raíz del $MSE$, sería tentador aprovecharnos de la función de este y elevarlo a una potencia de 0.5. Aunque esto puede ser eficiente, me gustaría agregar la función para Python.

Función de la Raíz del Error Cuadrático Medio:

1 2 3 4 5 6

# RMSE def rmse(y_m, y_e): sum1 = 0.0 for i in range(len(y_m)): sum1 += ((y_m[i] - y_e[i]) ** 2) return (sum1 / len(y_m)) ** 0.5

Como se observa, el único cambio fue regresar el valor, tal como se dijo, elevado a 0.5. Finalmente, la función del Coeficiente de Determinación usa dos sumas sum1, y sum2:

1 2 3 4 5 6 7 8 9

# R^2 def r2coef(y_m, y_e): sum1 = 0.0 sum2 = 0.0 mean = sum(y_m) / len(y_m) for i in range(len(y_m)): sum1 += ((y_m[i] - y_e[i]) ** 2) sum2 += ((y_e[i] - mean) ** 2) return 1.0 - (sum1 / sum2)

Se notará que aquí que simplemente se aprovecha una función sum() para encontrar la media del vector y_m. Aunque hay muchas formas de escribir las funciones precedentes en Python (incluso en una línea). La idea del desglose de línea a línea y el uso del ciclo for es con el propósito de que estas sean más inteligibles y susceptibles a adaptarse en otros lenguajes de programación.

(2) Los datos

La Figura 1 muestra la gráfica de los datos a tensión uniaxial de un caucho vulcanizado. Para el propósito de esta entrada, los datos corresponden a una matriz de que tiene tres columnas: (a) una correspondiente a los datos experimentales de la deformación, (b) una que contiene los datos del esfuerzo experimental, y (c) una que corresponde a datos obtenidos por un modelo hiperelástico (línea oscura). Los datos son un archivo CSV y pueden ser descargados de esta liga.

Figura 1. Gráfica de esfuerzo/deformación de un caucho natural (Treloar, 1944).

(3) El experimento con pandas

Toda vez que se han definido las funciones y los datos con los que se trabajará, lo siguiente corresponde a hacer uso de la libreria pandas para cargar estos últimos y crear una estructura con la cual llevar a cabo el cálculo del error.

1 2 3 4 5 6 7 8 9 10 11 12 13

import pandas as pd data = pd.read_csv('treloar_uniaxial.csv') s_e = data.iloc[:, 1] s_m = data.iloc[:, 2] error_MSE = mse(s_m, s_e) r2 = r2coef(s_m, s_e) error_RMSE = rmse(s_m, s_e) print(error_MSE) print(error_RMSE) print(r2)

La primera línea importa pandas; la tercera línea carga los datos a través de la función read_csv() de dicha librería; en la variables s_e y s_m de las líneas 4 y 5 almacenamos los datos del esfuerzo experimental y del modelo respectivamente; en las líneas 7, 8 y 9 mandamos a ejecutar las funciones y almacenamos sus valores de retorno en las variables que se muestran; finalmente, las líneas 11, 12 y 13 mandan a imprimir en pantalla los errores calculados:

0.02960504181818179 #MSE
0.17206115720342516 #RMSE
0.9920784500779709 #R^2

Estos resultados indican que se tiene un $MSE = 0.0296$, un $RMSE = 0.1720$, y un $R^2 = 0.9920$. Interpretando las cifras, se puede llegar a la conclusión de que es un modelo muy preciso, pues las estimaciones sugieren errores que se ubican menores al 2% respecto a los datos experimentales.

(4) La validación con scikit learn

Scikit learn es una librería muy usada para llevar a cabo cómputos de diversa naturaleza en la ciencia de datos. Como tal, contiene medidas de error ya implementadas para ser de uso fácil. En este caso, ya que es una librería bien construida y mantenida por equipos de ingeniería, quise validar mis funciones con las que ya tiene integradas scikit learn. En este caso, sólo se validará el $MSE$ y el $R^2$:

1 2 3 4 5 6 7 8 9 10 11 12 13

from sklearn.metrics import mean_squared_error #carga de funcion MSE from sklearn.metrics import r2_score #carga de función r2 error_MSE = mse(s_m, s_e) errorScikit = mean_squared_error(s_m, s_e) r2 = r2coef(s_m, s_e) r2scikit = r2_score(s_m, s_e) print('MSE de nuestra funcion es: {}'.format(error_MSE)) print('MSE scikit es: {}'.format(errorScikit)) print('R2 de nuestra funcion es: {}'.format(r2)) print('R2 scikit es: {}'.format(r2scikit))

Y el resultado es:

MSE de nuestra funcion es: 0.02960504181818179
MSE scikit es: 0.02960504181818179
R2 de nuestra funcion es: 0.9920784500779709
R2 scikit es: 0.9929096818928284

Tal como se observa, los resultados son muy similares. Por lo que la implementación de las funciones de error, es ciertamente transparente y se obtienen resultados correctos. Aunque las funciones pueden refinarse para incluir estados de excepción para contemplar errores en el código, lo elemental se ha cubierto.

Conclusión

En esta entrada se han presentado los códigos para implementar las funciones de error usando los datos experimentales del esfuerzo de un material y las estimaciones de un modelo. Se ha visto que los resultados de las mismas son muy similares a los que se obtienen con librerías completas y bien mantenidas por una comunidad profesional. 

Por lo anterior, debe tomarse esta entrada como una forma de ilustrar el funcionamiento de las medidas del error que se usan en el modelado constitutivo de materiales. Se recomienda, sin embargo, hacer uso de las que ya existen como las que tiene integradas scikit learn en su código.

Finalmente, el código completo se encuentra al final de esta entrada. 

Saludos, 

Alejandro 

Referencias

L. R. G. Treloar (1944). Stress-Strain Data for Vulcanized Rubber under Various Types of Deformation. Rubber Chemistry and Technology; 17 (4): 813–825. doi: https://doi.org/10.5254/1.3546701

# © 2020 Alejandro E. Rodriguez.  All rights reserved. For personal use only.
# It is not permitted to distribute verbatim copies of this document,
# and changing it is not allowed. Contact the owner of the code at blogmateria.com blog
# for further permissions.

# © 2020 Alejandro E. Rodríguez.  Todos los derechos reservados. Sólo para uso personal. 
# No está permitido distribuir copias de este documento, y no está permitido cambiarlo. 
# Contacta con el propietario del código en el blog de blogmateria.com para obtener más permisos.

import pandas as pd
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from matplotlib import pyplot as plt

data = pd.read_csv('treloar_uniaxial.csv')
s_e = data.iloc[:, 1]
s_m = data.iloc[:, 2]


def mse(y_m, y_e):
    sum1 = 0.0
    for i in range(len(y_m)):
        sum1 += ((y_m[i] - y_e[i]) ** 2)
    return sum1 / len(y_m)


def rmse(y_m, y_e):
    sum1 = 0.0
    for i in range(len(y_m)):
        sum1 += ((y_m[i] - y_e[i]) ** 2)
    return (sum1 / len(y_m)) ** 0.5


def r2coef(y_m, y_e):
    sum1 = 0.0
    sum2 = 0.0
    mean = sum(y_m) / len(y_m)
    for i in range(len(y_m)):
        sum1 += ((y_m[i] - y_e[i]) ** 2)
        sum2 += ((y_e[i] - mean) ** 2)
    return 1.0 - (sum1 / sum2)


error_MSE = mse(s_m, s_e)
errorScikit = mean_squared_error(s_m, s_e)

error_RMSE = rmse(s_m, s_e)

r2 = r2coef(s_m, s_e)
r2scikit = r2_score(s_m, s_e)

print('MSE de nuestra funcion es: {}'.format(error_MSE))
print('MSE scikit es: {}'.format(errorScikit))
print('R2 de nuestra funcion es: {}'.format(r2))
print('R2 scikit es: {}'.format(r2scikit))

plt.plot(data.iloc[:, 0], data.iloc[:, 2], 'k', label='Modelo', linewidth=2)
plt.plot(data.iloc[:, 0], data.iloc[:, 1], 'ro', label='Experimental')
plt.legend(loc="upper left")
plt.ylabel('Esfuerzo [MPa]')
plt.xlabel('Deformación [%]')
plt.show()