(#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




(#395). INTERPOLACIÓN (III): MÉTODO DE DIFERENCIAS DIVIDIDAS

[MONOTEMA] El polinomio de Lagrange se puede obtener también por medio de un método sencillo que emplea diferencias divididas de las imágenes con respecto a x. Este método también se conoce como “diferencias divididas de Newton”. 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 diferencias divididas

(1) Objetivo: Aproximarse numéricamente a la función , de la que conocemos sólo ciertos datos por medio de un polinomio que “pase” por los i puntos conocidos.

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Partimos de la serie de datos y definimos la primera diferencia dividida como

Si implementamos ese cálculo de manera recursiva llegamos a la fórmula general:

Así, construimos una tabla de n-1 diferencias divididas a partir de las cuales podemos hallar el polinimio de Lagrange:

(4) Estimación del error: Podemos emplear el mismo sistema que para el método de Lagrange

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 4 primeros puntos, con la programación siguiente:

x_:[1,3,4,6];
fx:[1,2,4,5];
valores: transpose(matrix(x_));
imagenes: transpose(matrix(fx));
[n,m]:matrix_size (valores);
funcion__1: zeromatrix(n,m);
for k:1 thru n-1 do
    (funcion__1[k,1]: (imagenes[k+1,1]-imagenes[k,1])/(valores[k+1,1]-valores[k,1]))$
print(funcion__1);
funcion__2: zeromatrix(n,m);
for k:1 thru n-2 do
    (funcion__2[k,1]: (funcion__1[k+1,1]-funcion__1[k,1])/((valores[k+2,1]-valores[k,1])))$
print(funcion__2);
funcion__3: zeromatrix(n,m);
for k:1 thru n-3 do
    (funcion__3[k,1]: (funcion__2[k+1,1]-funcion__2[k,1])/((valores[k+3,1]-valores[k,1])))$
print(funcion__3);
matriztotal: addcol(valores,imagenes,funcion__1,funcion__2, funcion__3);
p3:matriztotal[1,2]+((x-matriztotal[1,1])*matriztotal[1,3])+((x-matriztotal[1,1])*(x-matriztotal[2,1])*matriztotal[1,4])
+((x-matriztotal[1,1])*(x-matriztotal[2,1])*(x-matriztotal[3,1])*matriztotal[1,5]);

La tabla resultante (en forma de matriz) es:

b395_2donde en la primera columna tenemos los valores de   en la segunda los valores de y en las 3 restantes los valores computados por el procedimiento de Newton.

Por tanto, podemos construir el polinomio de Lagrange a partir de estos resultados, que nos lo da la parte de código ya implementado:

p3:matriztotal[1,2]+((x-matriztotal[1,1])*matriztotal[1,3])+((x-matriztotal[1,1])*(x-matriztotal[2,1])*matriztotal[1,4])
+((x-matriztotal[1,1])*(x-matriztotal[2,1])*(x-matriztotal[3,1])*matriztotal[1,5]);

El polinomio resultante es:

Que como vemos es un polinomio de orden 3, ya que hemos cogido 4 puntos para la interpolación. Si ahora evaluamos el polinomio en el punto de interés x=2:

solucion:ev(p3,x=2);

Obtenemos el resultado de:

que es exactamente lo mismo que lográbamos con la implementación directa de la fórmula del polinomio de Lagrange.

Diferencias hacia atrás

La obtención del polinomio de Lagrange puede lograrse cambiando la forma de recorrer la matriz de diferencias devididas. En lugar de hacerlo desde lo hacemos desde . Gráficamente quizá se entiende mejor: en verde el método hacia adelante y en rojo el método hacia atrás.

b395_3

b395_4

Los resultados son idénticos, y lo podemos comprobar con el siguiente código:

p3_atras:matriztotal[4,2]+((x-matriztotal[4,1])*matriztotal[3,3])+((x-matriztotal[4,1])*(x-matriztotal[3,1])*matriztotal[2,4])
+((x-matriztotal[4,1])*(x-matriztotal[3,1])*(x-matriztotal[2,1])*matriztotal[1,5]), ratsimp;

Evidentemente, tenemos que tener cuidado de comenzar a contar las “x” desde el final en lugar que desde el principio.

Nodos con igual espaciado

Cuando los puntos (nodos) se pueden ordenar con igual espaciado, podemos también obtener el polinomio de Lagrange a través de la siguiente fórmula:

donde  n es el orden del polinomio, y son las diferencias divididas entre las imágenes consecutivas. Además, donde es el dato para el que se quiere estimar su imagen y es el dato de partida elegido para la interpolación. Finalmente es la holgura del espaciado entre los datos.

En Maxima, podemos implementar el procedimiento anterior intentando inferir el valor de x=5 cuando se cogen los nodos equiespaciados [4,6,8,10]. El valor de la fórmula anterior quedaría así

x_:[4,6,8,10];
fx:[4,5,9,10];
valores: transpose(matrix(x_));
imagenes: transpose(matrix(fx));
[n,m]:matrix_size (valores);
funcion__1: zeromatrix(n,m);
for k:1 thru n-1 do
    (funcion__1[k,1]: imagenes[k+1,1]-imagenes[k,1])$
print(funcion__1);
funcion__2: zeromatrix(n,m);
for k:1 thru n-2 do
    (funcion__2[k,1]: funcion__1[k+1,1]-funcion__1[k,1])$
print(funcion__2);
funcion__3: zeromatrix(n,m);
for k:1 thru n-3 do
    (funcion__3[k,1]: funcion__2[k+1,1]-funcion__2[k,1])$
print(funcion__3);
matriztotal: addcol(valores,imagenes,funcion__1,funcion__2, funcion__3);
xs: 5;
xp:4;
h:2;
s:(xs-xp)/h;
p3_equiespaciado:matriztotal[1,2]+(s*matriztotal[1,3])+((s*(s-1)/2)*matriztotal[1,4])
+(s*(s-1)*(s-2)/6)*matriztotal[1,5];

Que nos da un valor:

Conclusión

El método de diferencias divididas proporciona una forma simplificada de obtener el polinomio de Lagrange. Sin embargo, está sujeto a las mismas limitaciones, es decir, hemos de ser muy cuidadosos a la hora de eleguir los nodos de interpolación.

Todos los posts relacionados




(#394). INTERPOLACIÓN (II): MÉTODO DE NEVILLE

[MONOTEMA] Una manera de simplificar los cálculos globales del método de Lagrange es emplear la propuesta de Neville, que genera los polonimios de interpolación de forma recursiva en varias etapas de cálculo.

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 Neville

(1) Objetivo: Aproximarse numéricamente a la función , de la que conocemos sólo ciertos datos  por medio de un polinomio que “pase” por los i puntos conocidos.

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Partimos de la serie de datos  y definimos  los elementos de la tabla de interpolación:

Básicamente se trata de ir interpolando entre pares de valoremos consecutivos en la primera columna, luego entre pares de valores que se llevan dos filas en la segunda, tres filas en la tercera, etc.

Así, construimos una tabla que es una matriz de n filas por n columnas que es triangular y cuyo último elemento de la diagonal coincide con el valor del polinomio de Lagrange entre todos los puntos considerados.

(4) Estimación del error: Podemos emplear el mismo sistema que para el método de Lagrange

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_: matrix([1,3,4]);
X_t:transpose(X_);
FX_:matrix([1,2,4]);
FX_t:transpose(FX_);
[n,m] : matrix_size (FX_t) ;
Neville1: zeromatrix(n,m);
Neville2: zeromatrix(n,m);
Neville3 : zeromatrix(n,m);
x:2;
for k:1 thru n do
    (
    Neville1[k,1]:FX_t[k,1],
 for k:2 thru n do
    (
    Neville2[k,1]:(((x-X_t[k-1,1])*FX_t[k,1])-(x-X_t[k,1])*FX_t[k-1,1])/(X_t[k,1]-X_t[k-1]),
 for k:3 thru n do
    (
     Neville3[k,1]:((((x-X_t[k-2,1])*Neville2[k,1]))-((x-X_t[k,1])*Neville2[k-1,1]))/(X_t[k,1]-X_t[k-2,1]))));
ResultadosNeville:addcol(X_t,Neville1,Neville2,Neville3);

La tabla resultante (en forma de matriz) es:

b394_2

donde en la primera columna tenemos los valores de y en las 3 restantes los valores computados por el procedimiento de Neville. Es decir los valores de los puntos más la tabla de Neville. Nótese que la columna dos son las imágenes  de la primera columna.

Por tanto, tenemos que, la estimación es 1 (el último elemento de la diagonal de las columnas de Neville), que coincide con el resultado calculado por el método de Lagrange.

En Maxima, lo que hemos hecho es crear una matriz de 3×3 (ya que usamos 3 puntos para la interporlación de x=2) con ceros a la que vamos a ir añadiendo en cada columna los valores calculados por el procedimiento de Neville.

Conclusión

El método de Neville realiza una interpolación iterada de manera secuencial, usando las estimaciones previas del polinomio de Lagrange, simplificando el proceso de cálculo.

Todos los posts relacionados

 




(#393). INTERPOLACIÓN (I): POLINOMIO DE LAGRANGE

[MONOTEMA] Para los alumnos de marketing es tarea fundamental conocer cómo se pueden hacer predicciones a partir de datos empíricos. Vamos a explicar en una serie de posts diversas formas de hacerlo, comenzando con una de las más conocidas, el polinomio de Lagrange. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017). El objetivo es, asimismo, que los estudiantes sean autónomos y sepan programar este tipo de métodos, y para ello proveeremos de códigos de Maxima.

Datos de partida

Imaginemos que partimos de los siguientes datos empíricos:

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.

Lo más indicado es proceder primeramente con una representación gráfica de los datos. Recordemos que el problema aquí es que conocemos los datos pero no conocemos cómo se relacionan, es decir, sabemos que existe una función , que depende obviamente de x, pero no sabemos cómo es.

En Maxima podemos inspeccionar esos datos 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_],
    [x,0,15],[y,0,10], [style, points], [color,green],
[xlabel, “Número de empleados”],
    [ylabel, “Satisfacción del cliente”], [legend, false]);

Vemos que en principio existe una asociación positiva entre ambas variables, pero nos interesaría saber cuál será el nivel de satisfacción si contamos con empleados que no estaban en los datos originales. Además, esa relación no está muy claro que sea lineal, al menos no en todo el dominio de la función, por lo que necesitamos recurrir a herramientas matemáticas para tratar de ser más precisos en el análisis y tener menos riesgo en la toma de decisiones.

Método de Lagrange

(1) Objetivo: Aproximarse numéricamente a la función , de la que conocemos sólo ciertos datos  por medio de un polinomio que “pase” por los i puntos conocidos.

(2) Condiciones: La función debe ser continua en el intervalo considerado.

(3) Descripción rápida: Partimos de dos puntos y . Definimos las funciones:

El polinomio de interpolación de Lagrange entre esos dos puntos es:

Por lo que tendremos un polinomio de grado m=1 que pasa por .Es decir, para n puntos vamos a considerar un polinomio de grado m=n-1

De forma más compacta y general:

Por ejemplo, para n=3 puntos, entonces el polinomio sería de orden m=2, y deberíamos definir n=3 sumandos en  , donde cada será la múltiplicación de m=2 términos. Es decir:

(4) Estimación del error: Para estimar el error necesitamos conocer la forma de la función real:

donde  es un número que está dentro  del intervalo de datos considerado. De este modo, únicamente si conocemos podremos derivarla n veces. Por tanto, si tenemos n=3 puntos habrá que hacer n=3 derivadas para obtener el error de un polinomio de orden m=2. En la mayoría de los casos no conoceremos esa función, ya que ese desconocimiento precisamente es la causa por la cual necesitamos interpolar. No obstante, si conocemos la función, podemos calcular una cota de error máximo para cualquier punto contenido dentro del intervalo considerado.

Implementación en Maxima

Vamos a realizar nuestra primera estimación en Maxima. Supongamos que queremos aproximarnos a un valor que desconocemos, cuando el número de empleados es 2. Para ello tenemos diversas opciones. La primera es considerar sólo un número de puntos cercano al objetivo, y la segunda considerar toda la información disponible. Ambas tienen ventajas e inconvenientes. Comencemos con la primera opción, calculando un polinomio de Lagrange de orden 2 que pase por los 3 primeros puntos de nuestros datos.

p(x):=(x-x1)*(x-x2)*fx0/((x0-x1)*(x0-x2))+
(x-x0)*(x-x2)*fx1/((x1-x0)*(x1-x2))+
(x-x0)*(x-x1)*fx2/((x2-x0)*(x2-x1));
expand(p(x));
x0:1;
x1:3;
x2:4;
fx0:1;
fx1:2;
fx2:4;
solucion1: expand(p(x));
f_solucion1: ev(solucion1,x=2);

x_:[1,3,4,6,7,8,9,10,11,12];
fx_:[1,2,4,5,7,9,2,10,10,10];
x_solucion1:2;
plot2d([[discrete, x_, fx_], [discrete, [[x_solucion1,f_solucion1]]],solucion1],
 [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]);

El polinomio resultante es:

que nos indica que la satisfación del cliente será:

Y su representación gráfica:

b393_3Vamos a ir comprobando qué ocurriría con la estimación a medida que incorporamos nueva información en el polinomio de Lagrange.

Si ahora usamos n=4 puntos:

kill (all);

p(x):=(x-x1)*(x-x2)*(x-x3)*fx0/((x0-x1)*(x0-x2)*(x0-x3))+
(x-x0)*(x-x2)*(x-x3)*fx1/((x1-x0)*(x1-x2)*(x1-x3))+
(x-x0)*(x-x1)*(x-x3)*fx2/((x2-x0)*(x2-x1)*(x2-x3))+
(x-x0)*(x-x1)*(x-x2)*fx3/((x3-x0)*(x3-x1)*(x3-x2));
expand(p(x));
x0:1;
x1:3;
x2:4;
x3:6;
fx0:1;
fx1:2;
fx2:4;
fx3:5;
solucion1: expand(p(x));
f_solucion1: ev(solucion1,x=2);

x_:[1,3,4,6,7,8,9,10,11,12];
fx_:[1,2,4,5,7,9,2,10,10,10];
x_solucion1:2;
plot2d([[discrete, x_, fx_], [discrete, [[x_solucion1,f_solucion1]]],solucion1],
 [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]);

El polinomio resultante es:

que nos indica que la satisfación del cliente será:

Y su representación gráfica:

Vemos como ahora la estimación cambia considerablemente al cambiar también la forma de la función de interpolación.

Pero si empleamos todos los puntos disponibles:

p(x):=(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx0/
((x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)*(x0-x8)*(x0-x9))+
(x-x0)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx1/
((x1-x0)*(x1-x2)*(x1-x3)*(x1-x4)*(x1-x5)*(x1-x6)*(x1-x7)*(x1-x8)*(x1-x9))+
(x-x0)*(x-x1)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx2/
((x2-x0)*(x2-x1)*(x2-x3)*(x2-x4)*(x2-x5)*(x2-x6)*(x2-x7)*(x2-x8)*(x2-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx3/
((x3-x0)*(x3-x1)*(x3-x2)*(x3-x4)*(x3-x5)*(x3-x6)*(x3-x7)*(x3-x8)*(x3-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx4/
((x4-x0)*(x4-x1)*(x4-x2)*(x4-x3)*(x4-x5)*(x4-x6)*(x4-x7)*(x4-x8)*(x4-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x6)*(x-x7)*(x-x8)*(x-x9)*fx5/
((x5-x0)*(x5-x1)*(x5-x2)*(x5-x3)*(x5-x4)*(x5-x6)*(x5-x7)*(x5-x8)*(x5-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x7)*(x-x8)*(x-x9)*fx6/
((x6-x0)*(x6-x1)*(x6-x2)*(x6-x3)*(x6-x4)*(x6-x5)*(x6-x7)*(x6-x8)*(x6-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x8)*(x-x9)*fx7/
((x7-x0)*(x7-x1)*(x7-x2)*(x7-x3)*(x7-x4)*(x7-x5)*(x7-x6)*(x7-x8)*(x7-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x9)*fx8/
((x8-x0)*(x8-x1)*(x8-x2)*(x8-x3)*(x8-x4)*(x8-x5)*(x8-x6)*(x8-x7)*(x8-x9))+
(x-x0)*(x-x1)*(x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6)*(x-x7)*(x-x8)*fx9/
((x9-x0)*(x9-x1)*(x9-x2)*(x9-x3)*(x9-x4)*(x9-x5)*(x9-x6)*(x9-x7)*(x9-x8)) ;
x0:1;
x1:3;
x2:4;
x3:6;
x4:7;
x5:8;
x6:9;
x7:10;
x8:11;
x9:12;
fx0:1;
fx1:2;
fx2:4;
fx3:5;
fx4:7;
fx5:9;
fx6:2;
fx7:10;
fx8:10;
fx9:10;
solucion1: expand(p(x));
f_solucion1: ev(solucion1,x=2);

El polinomio resultante es de noveno grado, y cuando lo evaluamos en x=2 nos indica que la satisfación del cliente es:

Y su representación gráfica:

b393_5

Rápidamente se puede ver el “disparate” de esta interpolación que nos lleva un resultado absolutamente sin sentido en el contexto de los datos.

Aunque en posteriores posts de esta serie comentaremos con mayor profundidad el concepto de sobreajuste, lo que ha sucedido es precisamente un ajuste innecesario de todos los puntos, que indica que hemos modelado no sólo la “señal” subyacente, sino el “ruido”. Ambos conceptos, señal y ruido, son fundamentales para entender la naturaleza de las estimaciones, inferencias e interpolaciones.

Una vez que ya hemos practicado con la programación “a mano”, podemos pedirle a Maxima que nos “ayude” con la interpolación usando una de sus funciones propias del programa, lo que resulta mucho más sencillo:

kill (all);

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

Que nos da la solución anterior con mucho menos esfuerzo de programación.

Conclusión

En este post hemos explicado brevemente en qué consiste el método de los polinomios de Lagrange para interpolar una función. Antes de su aplicación, hay que estudiar previamente la naturaleza de los datos, hacer inspecciones gráficas, y valorar en qué medida son necesarios escoger todos los datos, o sólo aquellos que estén en la vecindad del punto que se quiere interpolar.

Aprenderemos a evaluar en qué situaciones conviene su implementación o es mejor realizar ajustes estadísticos a los datos (como por mínimos cuadrados ordinarios), que nos permitan minimizar y analizar mejor la señal, pese a que podamos perder precisión para algún punto específico que queramos predecir.

Todos los posts relacionados