(#398). AJUSTE DE FUNCIONES (II): CAPACIDAD EXPLICATIVA Y PARTICIÓN DE LA FUNCIÓN

[MONOTEMA] Continuamos con una introducción sencilla al análisis de datos. Tras comparar los resultados de las splines cúbicas y el método de mínimos cuadrados, debemos ahora plantearnos algunas cuestiones sobre la idoneidad de lo que hemos hecho hasta ahora, y las opciones que aparecen entonces. Lo vamos a hacer de forma muy simple, para estimular el interés de los estudiantes hacia cuestiones más complejas.

Datos de partida

Recordemos que tenemos los siguientes datos analizados con splines y mínimos cuadrados:

load (interpol);
p:[[1,1],[3,2],[4,4], [6,5],[7,7],[8,9],[9,2],[10,10],[11,10],[12,10]];
splines: cspline(p);

x_:[1,3,4,6,7,8,9,10,11,12];
fx:[1,2,4,5,7,9,2,10,10,10];
x_aproximar:2;
x_interpolar:x_aproximar;
fx_estimada:0.8212147134302823*z+0.1693755346449957;
fx_estimada_x_aproximar: ev(fx_estimada, z=x_aproximar), numer;
f_splines: ev(splines,x=x_interpolar), numer;
plot2d([[discrete, x_, fx], [discrete, [[x_aproximar,fx_estimada_x_aproximar],
[x_interpolar,f_splines]]],fx_estimada,splines],
[x,0,15],[y,0,10], [style, points, points, lines, lines], [color, green, red, magenta, orange, magenta],
[xlabel, “Número de empleados”],
[ylabel, “Satisfacción del cliente”], [legend, false]);

b397_4Capacidad explicativa

La capacidad explicativa del modelo ajustado por mínimos cuadrados se suele medir con el coeficiente de determinación (R-cuadrado), que es una medida de la varianza explicada por el modelo, es decir, cuantifica la importancia de la varianza de error, y que está entre 0 y 1 (siendo 1 el ajuste perfecto). Se obtiene por el cociente entre la varianza de los datos ajustados frente a la de las imágenes de los datos observados.

En Maxima podemos obtenerlo así:

load(descriptive);
varianza_fx:var1(fx_traspuesta), numer;
ajuste: transpose(ev(fx_estimada,z=x));

varianza_fx_estimada:var1(ajuste);

R_cuadrado:(varianza_fx_estimada/varianza_fx);

El resultado es:

Pero debemos plantearnos dos cuestiones fundamentales tras inspeccionar gráficamente los datos. En primer lugar, la hipótesis de relación lineal parece no cumplirse a medida que el número de empleados crece. Y en segundo lugar el dato [9,2] puede ser un dato atípico, un outlier que por diferentes razones (sean aleatorias o no), nos está condicionando la estimación, y que quizá habría que estudiar obviar.

Parte lineal

¿Cómo sería el ajuste si sólo nos centramos en el rango de datos [1,8]? Vamos a verlo:

x:matrix([1,3,4,6,7,8]);
fx:matrix([1,2,4,5,7,9]);
x_traspuesta: transpose(x);
fx_traspuesta:transpose(fx);
n: length(x_traspuesta);
datos: zeromatrix(n,1);
datos_xcuad:x_traspuesta.x;
datos_y:fx_traspuesta;
datos_xy:fx_traspuesta.x;
datos_x: x_traspuesta;
print(datos_xcuad,datos_y, datos_xy,datos_x);
suma_xcuad:sum(datos_xcuad[i,i],i,1,n);
suma_y:sum(datos_y[i,1],i,1,n);
suma_xy:sum(datos_xy[i,i],i,1,n);
suma_x:sum(datos_x[i,1],i,1,n);
cuad_suma_x:suma_x^2;
alpha:((suma_xcuad*suma_y)-(suma_xy*suma_x))/((n*suma_xcuad)-cuad_suma_x), numer;
beta:((n*suma_xy)-(suma_x*suma_y))/((n*suma_xcuad)-cuad_suma_x), numer;
fx_estimada:alpha+beta*z;

Lo que nos da la siguiente recta de ajuste (la estimación de la función).

Entonces podemos hacer la representación gráfica del ajuste de los datos:

b398_3

La recta azul es la nueva función de ajuste, cuya capacidad explicativa es:

El ajuste es mucho mejor, claramente. Lo que sucede es que sólo podemos hacer esta acción si podemos justificar muy bien esa estimación de mínimos cuadrados para sólo una parte de los datos empíricos.

Los outliers pueden surgir aleatoriamente, y es parte de la variabilidad global, pero en detemrinados casos podemos intentar eliminarlos si sospechamos que hay “algo raro” en ese dato (alguna variación sistemática que ha producido un valor extraño).

Viendo gráficamente los datos, y teniendo en cuenta el contexto teórico, está claro de que hay un punto de saturación donde por más que se añadan vendedores no se incrementa el nivel de satisfacción (es máximo ya con 10 vendedores). Por tanto, no tiene mucho sentido intentar modelar esa parte de los datos junto al resto.

Conclusión

Hemos visto que, a veces, y si la teoría y la perspicacia del investigador lo permiten, podemos descartar algunos datos para obtener ajustes mucho mejores, y de este modo, realizar predicciones más fiables en el rango de datos que nos interese. No siempre, desde luego, ello es recomendable, pero la idea es que a los datos hay que mirarlos muy bien antes de analizarlos.

Hay dos aspectos fundamentables que no hemos tocado aquí, pero que son inherentes a cualquier éxito en los análisis: (1) El modelo teórico tiene que estar causalmente bien especificado, no importa la R-cuadrado (aunque sea alto) si no es teóricamente plausible esa relación; (2) El planteamiento de un modelo estadístico es también el de sus asunciones, las cuales no hemos testado en su totalidad, como por ejemplo, si la varianza de error es homogénea.

En futuros posts iremos aclarando mejor estas y otras cuestiones.

Todos los posts relacionados

 

 

 




(#397). AJUSTE DE FUNCIONES (I): SPLINES VS MÍNIMOS CUADRADOS

[MONOTEMA] Tras explicar varios de los más empleados métodos de interpolación para buscar aproximarnos a la función que describe los datos empíricos, hemos visto que las splines cúbicas nos ofrecen mucha flexibiliad. Sin embargo, y aunque pueda parecer paradójico, el hecho de buscar una función de interpolación que pase por todos los nodos puede ser contraproducente, porque podemos estar dando demasiada importancia al ruido, cuando lo que buscamos es la señal. Es decir, cuando en los datos hay errores aleatorios (y también algún error sistemático puntual), buscar aproximaciones con funciones que no pasen necesariamente por esos nodos pero que tengan otras propiedades (como la minimización de una variable de error), puede ser más útil.

El método básico para realizar este ajuste es el de mínimos cuadrados, un potente método no paramétrico que minimiza la suma al cuadrado de los errores. Y es lo que vamos a ver en este post, aunque da manera muy simplificada ya que está sobradamente explicado en multitud de textos especializados de prácticamente todas las áreas científicas.

Así, compararemos gráficamente el ajuste por mínimos cuadrados con el método de las splines cúbicas, para estimular la reflexión de los alumnos.

Datos de partida

Usaremos los mismos datos que en el método de Lagrange, y en el de todos los ejemplos sobre interpolación:

x f(x)
1 1
3 2
4 4
6 5
7 7
8 9
9 2
10 10
11 10
12 10

Estos datos relacionan la cantidad de empleados utilizados en un gran supermercado (x) con la saltisfacción del cliente . La satisfacción del cliente se mide en una escala de 0 a 10, donde 0 es el valor mínimo y 10 el valor máximo.

x_:[1,3,4,6,7,8,9,10,11,12];
fx_:[1,2,4,5,7,9,2,10,10,10];
plot2d([discrete, x_, fx_],
[x,0,15],[y,0,10], [style, points], [color,green],
[xlabel, “Número de empleados”],
[ylabel, “Satisfacción del cliente”], [legend, false]);

Splines cúbicos

Una vez que ya sabemos programar las splines paso a paso, podemos actuar de forma más directa con Maxima, pidiendo que use todos los datos disponibles:

load (interpol);

p:[[1,1],[3,2],[4,4], [6,5],[7,7],[8,9],[9,2],[10,10],[11,10],[12,10]];

splines: cspline(p);

Y lo podemos graficar de la siguiente forma:

x_:[1,3,4,6,7,8,9,10,11,12];
fx:[1,2,4,5,7,9,2,10,10,10];
plot2d([[discrete, x_, fx], splines],
[x,0,15],[y,0,10], [style, points, lines], [color,green, orange],
[xlabel, “Número de empleados”],
[ylabel, “Satisfacción del cliente”], [legend, false]);

b397_2

Vemos que las splines hacen aparentemente un “buen trabajo” de ajuste, pero tenemos ese punto “problemático” [9,2] que hace que la curva de interpolación cambie drásticamente. Sin embargo, en este caso no nos afecta demasiado a nuestro objetivo de interpolar puntos que no están en los datos en el intervalo [1,12]. Tanto el nodo 2 como el nodo 5 se pueden evaluar con la curva de interpolación ya que las splines van de nodo a nodo, por lo que ese comportamiento “raro” en el punto [9,2] no afecta prácticamente a las evaluaciones entre los otros nodos.

Mínimos cuadrados lineales

Este método se basa en la minimización de la función de error:

donde  se corresponde con las imágenes de los datos brutos , y son parámetros a estimar.

Minimizar esa función requiere que:

Es importante señalar que se asume que el error tiene de media cero, es decir, es aleatorio (ruido blanco).

De este modo:

Así, podemos programarlo con Maxima usando el siguiente código:

x:matrix([1,3,4,6,7,8,9,10,11,12]);
fx:matrix([1,2,4,5,7,9,2,10,10,10]);
x_traspuesta: transpose(x);
fx_traspuesta:transpose(fx);
n: length(x_traspuesta);
datos: zeromatrix(n,1);
datos_xcuad:x_traspuesta.x;
datos_y:fx_traspuesta;
datos_xy:fx_traspuesta.x;
datos_x: x_traspuesta;
print(datos_xcuad,datos_y, datos_xy,datos_x);
suma_xcuad:sum(datos_xcuad[i,i],i,1,n);
suma_y:sum(datos_y[i,1],i,1,n);
suma_xy:sum(datos_xy[i,i],i,1,n);
suma_x:sum(datos_x[i,1],i,1,n);
cuad_suma_x:suma_x^2;
alpha:((suma_xcuad*suma_y)-(suma_xy*suma_x))/((n*suma_xcuad)-cuad_suma_x), numer;
beta:((n*suma_xy)-(suma_x*suma_y))/((n*suma_xcuad)-cuad_suma_x), numer;
fx_estimada:alpha+beta*x_;

Lo que nos da la siguiente recta de ajuste (la estimación de la función).

Entonces podemos hacer la representación gráfica del ajuste de los datos:

kill (all);
x_:[1,3,4,6,7,8,9,10,11,12];
fx:[1,2,4,5,7,9,2,10,10,10];
x_aproximar:2;
fx_estimada:0.8212147134302823*x+0.1693755346449957;
fx_estimada_x_aproximar: ev(fx_estimada, x=x_aproximar), numer;
plot2d([[discrete, x_, fx], [discrete, [[x_aproximar,fx_estimada_x_aproximar]]],fx_estimada],
 [x,0,15],[y,0,10], [style, points, points, lines], [color,green, red, orange],
[xlabel, “Número de empleados”],
    [ylabel, “Satisfacción del cliente”], [legend, false]);

b397_3Para 2 empleados la función de ajuste estimada nos dice que la satisfacción será de 1.81.

En próximos posts hablaremos con más detalle de la evaluación de este ajuste, pero ahora simplemente los estudiantes pueden hacer un ejercicio de reflexión ante el visionado de la superposición de las dos funciones computadas, la primera con splines cúbicas y la segunda con el método de los mínimos cuadrados.

load (interpol);
p:[[1,1],[3,2],[4,4], [6,5],[7,7],[8,9],[9,2],[10,10],[11,10],[12,10]];
splines: cspline(p);

x_:[1,3,4,6,7,8,9,10,11,12];
fx:[1,2,4,5,7,9,2,10,10,10];
x_aproximar:2;
x_interpolar:x_aproximar;
fx_estimada:0.8212147134302823*x+0.1693755346449957;
fx_estimada_x_aproximar: ev(fx_estimada, x=x_aproximar), numer;
f_splines: ev(splines,x=x_interpolar), numer;
plot2d([[discrete, x_, fx], [discrete, [[x_aproximar,fx_estimada_x_aproximar],
        [x_interpolar,f_splines]]],fx_estimada,splines],
 [x,0,15],[y,0,10], [style, points, points, lines, lines], [color, green, red, magenta, orange, magenta],
[xlabel, “Número de empleados”],
    [ylabel, “Satisfacción del cliente”], [legend, false]);

b397_4

Cuando el número de empleados es 2, la estimación del nivel de satisfacción difiere en función del método empleado. Además, hemos realizado el ajuste de mínimos cuadrados sin tener en cuenta que aparentemente hay una ruptura de la linealidad cuando el número de empleados se incrementa. En posteriores posts, analizaremos con más tranquilidad estos detalles.

Conclusión

Hemos comparado los resultados obtenidos con splines cúbicas con respecto al ajuste de mínimos cuadrados ordinarios, sin entrar en profundidades sobre la idoneidad de ambos análisis. Las predicciones que podemos hacer en base a los datos empíricos divergen, lo que nos debe advertir de lo prudentes que hemos de ser a la hora de analizar los datos.

Todos los posts relacionados

 

 

 




(#396). INTERPOLACIÓN (IV): MÉTODO DE SPLINE CÚBICO

[MONOTEMA] Los métodos de interpolación vistos hasta ahora no permitían aproximar a la función por tramos. El método de spline sí que lo puede hacer, y es muy potente para conseguir buenos resultados en cualquier tipo de análisis de datos. Seguiremos a Burden, Faires & Burden (2017), y lo enfocaremos, como siempre, de manera simplificada.

Datos de partida

Usaremos los mismos datos que en el método de Lagrange:

x f(x)
1 1
3 2
4 4
6 5
7 7
8 9
9 2
10 10
11 10
12 10

Estos datos relacionan la cantidad de empleados utilizados en un gran supermercado (x) con la saltisfacción del cliente . La satisfacción del cliente se mide en una escala de 0 a 10, donde 0 es el valor mínimo y 10 el valor máximo.

x_:[1,3,4,6,7,8,9,10,11,12];
fx_:[1,2,4,5,7,9,2,10,10,10];
plot2d([discrete, x_, fx_],
[x,0,15],[y,0,10], [style, points], [color,green],
[xlabel, “Número de empleados”],
[ylabel, “Satisfacción del cliente”], [legend, false]);

Método de spline cúbico

(1) Objetivo: Aproximarse numéricamente a la función , de la que conocemos sólo ciertos datos por medio de varios polinomios que pasen por los i puntos conocidos, es decir, calcularemos polinomios por tramos de datos. Emplearemos, además, splines cúbicos, es decir, polinomios de orden 3.

(2) Condiciones: La función debe ser continua en el intervalo considerado y además tiene que ser doblemente diferenciable.

(3) Descripción rápida: Partimos de la serie de datos , los cuales es conveniente ordenar de manera creciente por los nodos ,  y escogemos los j subintervalos para los cuales vamos a calcular splines .

La condiciones que deben satisfacerse son las siguientes:

a) son los j polinomios que se calculan para los i=j+1 nodos.

b) y para cada j. Es decir, el spline cúbico debe pasar por todas las imágenes de los datos, y esto hace que dos splines consecutivos pasen por un nodo común, por lo que al evaluar ambos splines en el mismo punto el resultado deba ser el mismo.

c) para cada j=1…n-2.  No sólo las evaluaciones de los splines, sino también sus derivadas primera y segunda deben coincidir en los nodos de conexión.

d) Si el spline es natural (frontera libre), entonces: 

Los splines de frontera condicionada no los vamos a comentar ya que es necesario conocer los valores de la derivada de la función original en los extremos, algo de lo que carecemos, debido a que sólo tenemos un conjunto de datos brutos y sus imágenes.

Estas condiciones suponen unas restricciones en la resolución de un sistema de ecuaciones con 4xj incógnitas. Es decir, 4 incógnitas por cada spline calculado. Así, para n=3 puntos obtendríamos j=2 splines cúbicos, de la forma:

Ahora tenemos que determinar esas constantes usando las restricciones comentadas anteriormente.

b) , ,

c)  , por lo que 

y , por lo que  

d) , por lo que

y , por lo que 

Ya sólo queda resolver el sistema de ecuaciones anterior, es decir, y considerando los  nodos [1,3,4] e imágenes [1,2,4]:

Implementación en Maxima

Vamos a realizar nuestra estimación en Maxima. Supongamos, como en el caso del método de Lagrange, que queremos aproximarnos a un valor que desconocemos, cuando el número de empleados es 2. Para ello vamos a coger los 3 primeros puntos, con la programación siguiente:

x_:[1,3,4];
fx:[1,2,4];
S1: a1+b1*(x-x_[1])+c1*(x-x_[1])^2+d1*(x-x_[1])^3;
S2: a2+b2*(x-x_[2])+c2*(x-x_[2])^2+d2*(x-x_[2])^3;

fx1:ev(S1,x=x_[1]);

fx2:ev(S1,x=x_[2]);

fx2_:ev(S2,x=x_[2]);

fx3:ev(S2,x=x_[3]);

S1prima: diff(S1,x);

S1prima_ev:ev(S1prima,x=x_[2]);

S2prima: diff(S2,x);

S2prima_ev:ev(S2prima,x=x_[2]);

S1primaprima: diff(S1prima,x);

S1primaprima_ev:ev(S1primaprima,x=x_[2]);

S2primaprima: diff(S2prima,x);

S2primaprima_ev:ev(S2primaprima,x=x_[2]);

S1primaprimainicio_ev:ev(S1primaprima,x=x_[1]);

S2primaprimafin_ev:ev(S2primaprima,x=x_[3]);

e1: fx1-fx[1];
e2: fx2-fx[2];
e3: fx2_-fx[2];
e4:fx3-fx[3];
e5:S1prima_ev-S2prima_ev;
e6:S1primaprima_ev-S2primaprima_ev;
e7:S1primaprimainicio_ev-0;
e8: S2primaprimafin_ev-0;

algsys ([e1, e2, e3, e4,e5,e6,e7,e8], [a1,b1,c1,d1,a2,b2,c2,d2]);

Aunque a priori parece un código algo complejo, lo cierto es que sigue paso a paso el algoritmo que deviene de las restricciones planteadas. Para formular el sistema de ecuaciones a resolver, hay que recordar que debemos escribir esas ecuaciones en  Maxima de la forma x-a=0 y no de la forma x=a.

Con el comando “algsys” resolvemos el sistema de ecuaciones, y obtenemos la siguiente solución:

Por tanto, podemos construir las dos splines:

Y ahora hacemos la comprobación con Maxima. Esas splines evaluadas en los nodos tienen que dar sus imágenes correspondientes:

S1solucion: 1+(1/8)*(x-x_[1])^3;
S2solucion: 2+(3/2)*(x-x_[2])+(3/4)*(x-x_[2])^2-(1/4)*(x-x_[2])^3;

Comprobacionnodo1: ev(S1solucion,x=x_[1]);

Comprobacionnodo2: ev(S1solucion,x=x_[2]);

Comprobacionnodo2_: ev(S2solucion,x=x_[2]);

Comprobacionnodo3: ev(S2solucion,x=x_[3]);

Como vemos, los resultados coinciden,  lo que indica que hemos hecho bien la programación y ya tenemos las 2 splines cúbicas.

Podemos dibujar las splines y estimar el valor para el cual queríamos hacer la interpolación, es decir, ¿cuál sería el valor de satisfacción del cliente con 2 empleados en el establecimiento?.

Para ello empleamos la primea spline que es la que pasa por los dos nodos que contienen al valor en cuestión, y su resultado es 1.125. Con Maxima lo podemos hacer fácilmente, y además representarlo en un gráfico:

x_:[1,3,4,6,7,8,9,10,11,12];
fx:[1,2,4,5,7,9,2,10,10,10];
x_interpolar:2;
f_S1solucion: ev(S1solucion,x=x_interpolar), numer;
plot2d([[discrete, x_, fx], [discrete, [[x_interpolar,f_S1solucion]]],S1solucion,S2solucion],
[x,0,15],[y,0,10], [style, points, points, lines,lines], [color,green, red, orange,magenta],
[xlabel, “Número de empleados”],
[ylabel, “Satisfacción del cliente”], [legend, false]);

b396_2La curva naranja representa la primera spline y la curva rosa la segunda.

Podemos hacer una comprobación con el comando de Maxima “cspline” y veremos que nos da el mismo resultado:

S1solucion: 1+(1/8)*(x-x_[1])^3, expand;
S2solucion: 2+(3/2)*(x-x_[2])+(3/4)*(x-x_[2])^2-(1/4)*(x-x_[2])^3,expand;

load (interpol);

p:[[1,1],[3,2],[4,4]];

cspline(p);

Conclusión

Las splines cúbicas son una muy interesante forma de realizar interpolaciones, probablemente una de las que más ventajas conlleva y que permite realizar aproximaciones más fiables. Su procedimiento en Maxima es un poco laborioso, pero si seguirmos ordenadamente los pasos que surgen de las restricciones del método, no debemos tener problema.

Todos los posts relacionados