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

 




(#391). BÚSQUEDA DE SOLUCIONES (V): MÉTODO DE LA POSICIÓN FALSA

[MONOTEMA] Una variante del método de la secante es el método de la posición falsa (Regula Falsi), que básicamente ejecuta el procedimiento de la secante pero cambiando levemente el algoritmo para que entre dos puntos de dos iteraciones sucesivas siempre se encuentra la solución, de forma análoga a como ocurre con el método de la bisección. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de la posición falsa

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

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

(3) Descripción rápida: Partimos del método de la secante:

Para calcular p2 empleamos la fórmula anterior (tras proveer las semillas p0 y p1). Sin embargo, para calcular p3, usaremos los valores (p2,p1,f(p2),f(p1)) sólo si f(p2) y f(p1) tienen signo opuesto. Si ambos tiene el mismo signo, entonces en lugar del par (p1,f(p1)), escogemos (p0,f(po)). Y se sigue ese mismo criterio para el resto de iteraciones.

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. Para emplear este criterio, hemos de añadir el requerimiento de que la función sea dos veces derivable.

Normalmente, el método de la posición falsa converge de forma más lenta que el de Newton y el de la secante.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de los dos valores con los que inicíabamos el método de la bisección: [2.5,3.5], que actúan como las primeras semillas en el método de la posición falsa, y vamos a adaptar en Maxima la propuesta de Ramírez (2008):

falsapos(arg):=block([a:arg[1],b:arg[2],f:arg[3],xm], xm:(a*ev(f,x=b)-b*ev(f,x=a))/(ev(f,x=b)-ev(f,x=a)),
if (ev(f,x=xm)*ev(f,x=a)<0) then [a,xm,f] else [xm,b,f] );
ini:[2.5,3.5,-(x^2)+(3*x)+0.1]; for i:1 thru 10 do (print(ini[1]),ini:falsapos(ini));

Y nos genera estos resultados:

2.5
2.95
3.021739130434783
3.031481481481481
3.03277399056109
3.03294493097818
3.032967529289182
3.032970516620634
3.032970911521146
3.032970963723678

Para trabajar un poco estos procedimientos, los estudiantes pueden intentar programar una rutina con Maxima donde se codifiquen los procedimientos de Newton, secante y posición falsa, y comparar las aproximaciones.

Conclusión

En este post hemos explicado brevemente en qué consiste el método de la posición falsa. Burden, Faires & Burden (2017) no recomiendan por lo general este método, pero no deja de ser interesante para trabajar las rutinas de programación y para entender los diferentes procedimientos numéricos de búsquedas de raíces.

Todos los posts relacionados




(#390). BÚSQUEDA DE SOLUCIONES (IV): MÉTODO DE LA SECANTE

[MONOTEMA] Una modificación del método de Newton permite una aproximación a este sin computar la derivada. Se trata del método de la secante, que explicaremos a continuación. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de la secante

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

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

(3) Descripción rápida: Partimos del método de Newton, que recordemos permite construir la sucesión donde:

Pero esa derivada puede ser compleja de calcular, por lo que podemos aproximarla de la siguiente forma:

 De este modo, obtenemos:

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. Para emplear este criterio, hemos de añadir el requerimiento de que la función sea dos veces derivable.

Normalmente, el método de la secante converge de forma más lenta que el de Newton, pero más rápido que otros, como el de la bisección.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de los dos valores con los que inicíabamos el método de la bisección: [2.5,3.5], que actúan como las primeras semillas en el método de la secante:

f(i):=-x[i]^2+3*x[i]+0.1;
x[0]:2.5;
x[1]:3.5;
for i:2 thru 8 do
    (
    x[i]:x[i-1]-((x[i-1]-x[i-2])*f(i-1)/(f(i-1)-f(i-2))),
    print(i-2,x[i-2]));
datos: makelist([[i-2], x[i-2]], i, 2, 8);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

b390_2

Como se puede observar, el método converge rápido, y también podemos establecer criterios de parada con programaciones similares a las realizadas con los otros métodos.

Otra forma de hacerlo

El método de la secante, como su propio nombre indica, traza rectas secantes a la función y que pasan por el eje X. La idea es ir aproximando secantes entre puntos de la función hasta converger en la solución. Para ello, es útil considerar la ecuación de la recta que une los puntos (a, f(a)) y (b, f(b)):

Podemos verlo gráficamente si estipulamos que a y b sean 2.5 y 3.5, respectivamente, que son los valores empleados como intervalos en el método de la bisección:

funcion:-(x^2)+(3*x)+0.1;
a:2.5;
fa:(ev(funcion,x=a));
b:3.5;
fb:(ev(funcion,x=b));
secante:((fb-fa)*(x-a)/(b-a))+fa, ratsimp;
f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion, secante],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines, lines], [color, red, green, orange], [legend, false]);

b390_3Esta primera recta secante, corta a la función en f(a) y f(b), y ya nos quedamos bastante cerca de la solución. Simplemente hay que seguir trazando secantes iterando entre puntos consecutivos, es decir, el siguiente paso sería realizar el mismo cálculo pero tomando f(b) en lugar de f(a), y f(m) en lugar de f(b), siendo f(m) el valor de la función para el punto m (el punto de corte de la primera secante con el eje X). En Maxima, una forma de computarlo se puede encontrar en Ramírez (2008). Nosotros vamos a adaptar ese código de la siguiente manera:

secante(arg):=block([a:arg[1],b:arg[2],f:arg[3],xsec], xsec:(a*ev(f,x=b)-b*ev(f,x=a))/(ev(f,x=b)-ev(f,x=a)),
[b,xsec,f] );

ini:[2.5,3.5,-(x^2)+(3*x)+0.1];

for i:1 thru 10 do (print(float(ini[1])),ini:secante(ini));

El resultado es:

2.5
3.5
2.95
3.021739130434783
3.033284564740307
3.0329698187459
3.032970971557676
3.032970971675589

Conclusión

En este post hemos explicado brevemente en qué consiste el método de la secante, una de las variantes del método de Newton, que provee una convergencia más lenta, pero que tiene la ventaja de no tener que calcular derivadas en cada iteración.

Todos los posts relacionados




(#389). BÚSQUEDA DE SOLUCIONES (III): METODO DE NEWTON

[MONOTEMA] Seguimos en la búsqueda de los ceros de una función con aproximaciones numéricas. Tras ver el método de la bisección y el método del punto fijo,  vamos a continuar con el método de Newton, también conocido como Newton-Raphson. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en los ejemplos anteriores, donde se relaciona la inversión en publicidad con los beneficios.

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
[3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”],
[ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método de Newton

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

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

(3) Descripción rápida: Partimos del primer polinomio de Taylor para  expandido alrededor de  y evaluado en  . Suponemos que es una buena aproximación a , es decir, que la semilla inicial está relativamente cerca de la solución. Cuando , entonces, obviando el término de error de la serie de Taylor:

A partir de esa fórmula podemos construir la sucesión donde:

(4) Convergencia: Un criterio suficiente (pero no necesario) para que lo haga es que:

Si se cumple la fórmula anterior cuando está en un intervalo alrededor de la raíz entonces la solución converge en cualquier punto de ese intervalo. La convergencia es no lineal (cuadrática), y esta es una de las razones por las que es uno de los procedimientos más empleados.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error relativo:

Implementación en Maxima

Vamos a implementar el algoritmo partiendo de 3 semillas diferentes, para ver que en función de la elección de ellas nos puede llevar a una de las dos soluciones posibles.  Para ello, empleamos un bucle con 10 iteraciones, tomando como semillas: 3, 0.1 y 8. Esto hace que evaluemos 3 funciones diferentes (f1, f2, f3) a la hora de analizar su convergencia.

f1(i):= -x1[i]^2+3*x1[i]+0.1;
derivada_f1(i):=-(2*x1[i])+3;
x1[0]:3;
for i:1 thru 10 do(
x1[i]:x1[i-1]-(f1(i-1)/derivada_f1(i-1)));
f2(i):= -x2[i]^2+3*x2[i]+0.1;
derivada_f2(i):=-(2*x2[i])+3;
x2[0]:0.2;
for i:1 thru 10 do(
x2[i]:x2[i-1]-(f2(i-1)/derivada_f2(i-1)));
f3(i):= -x3[i]^2+3*x3[i]+0.1;
derivada_f3(i):=-(2*x3[i])+3;
x3[0]:8;
for i:1 thru 10 do(
x3[i]:x3[i-1]-(f3(i-1)/derivada_f3(i-1)));
datos: makelist([i,”Iteración”,x1[i],f1(i-1), x2[i],f2(i-1),
    x3[i], f3(i-1)], i, 1,10);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

 b389_2

El método converge bastante rápido, incluso cuando la semilla se aleja de la raíz, y como puede observarse en función de la elección de la semilla converge hacia una de las dos soluciones.

Criterio de parada

Al igual que hicimos en el método de la bisección, podemos estipular que nos baste con una tolerancia de error, es decir, pedirle al programa que pare las iteraciones cuando se cumpla una condición determinada, por ejemplo, que el error relativo sea menor que una cota que fijemos. De este modo, si fijamos que el error relativo sea como máximo 0.0001:

f1(i):= -p[i]^2+3*p[i]+0.1;
derivada_f1(i):=-(2*p[i])+3;
p[0]:8;
for i:1 while Error_rel[i]>0.0001 do(
p[i]:p[i-1]-(f1(i-1)/derivada_f1(i-1)),
Error_abs[i]:abs(p[i]-p[i-1]),
Error_rel[i]: Error_abs[i]/abs(p[i]),
    print(i,p[i],Error_rel[i]));

Le hemos dicho al programa que fije un error relativo mayor que 0.0001, por lo que parará en esa iteración, y de esta manera sabremos que en la siguiente iteración el error relativo estará por debajo de esa cota. Si ejecutamos esas instrucciones en Maxima vemos que para en la iteración número 5, por lo que incluso con semillas alejadas de la raíz la convergencia puede ser muy rápida.

Conclusión

En este post hemos explicado brevemente en qué consiste el método de Newton, un algoritmo muy empleado en diversas aplicaciones estadísticas, como para obtener los estimadores que minimizan una función de error (hacen que su derivada sea cero, por lo que es realmente obtener la aproximación numérica a la raíz de una función). El método de Newton tiene, además, diversas variantes que comentaremos en futuras entradas.

Todos los posts relacionados




(#388). BÚSQUEDA DE SOLUCIONES (II): METODO DEL PUNTO FIJO

[MONOTEMA] Continuamos buscando los ceros de una función con aproximaciones numéricas. Tras ver el método de la bisección vamos a continuar con el método del punto fijo. Lo haremos, como siempre, desde el punto de vista práctico, y siguiendo a Burden, Faires & Burden (2017).

Función de partida

Emplearemos la misma función que en el método de la bisección, donde se relaciona la inversión en publicidad con los beneficios.

Lo primero que vamos a hacer es representar esa función en Maxima, tal y como hicimos en el post anterior:

funcion: -x^2+3*x+0.1;
plot2d(funcion,[x,0,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”], [color, green]);

Pero ahora vamos a expandir el eje de las X, porque nos interesa por motivos ilustrativos representar el corte en el tramo negativo:

funcion: -x^2+3*x+0.1;
plot2d(funcion,[x,-1,4],[y,-3,3],
[xlabel, “Inversión en publicidad (miles de euros)”], [ylabel, “Beneficios (miles de euros)”],
    [color, green]);

b388_2Y vemos que tenemos dos ceros (las dos veces que corta al eje X). Dijimos que la parte negativa no estaba dentro del dominio de la función, pero nos interesa representarla para ver lo que ocurre con la iteración del método del punto fijo.

De este modo tenemos dos ceros de la función  y .

Podemos representar esos puntos en el gráfico anterior de la siguiente manera:

f_ceronegativo:ev(funcion,x=-0.03297097167558977), numer;
f_ceropositivo:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[-0.0329709716755897,f_ceronegativo],
                [3.03297097167559, f_ceropositivo]]], funcion],[x,-1,4],[y,-3,3],
     [xlabel, “Inversión en publicidad (miles de euros)”],
     [ylabel, “Beneficios (miles de euros)”],
[style, points, lines], [color, red, green], [legend, false]);

Método del punto fijo

(1) Objetivo: Aproximarse numéricamente a la solución p de una función , tal que

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

(3) Descripción rápida: Si la función , entonces podemos construir la función . De este modo, hemos de escribir nuestra función original como . Esto significa para que haya un cero , es decir, que haya una x despejada en uno de los lados de la igualdad.

(4) Convergencia: No siempre converge. Un criterio suficiente (pero no necesario) para que lo haga es que  , donde es una constante positiva. Si el método converge cuadráticamente.

(5) Estimación del error: Hay diversas opciones, aunque una de las más recomendables es el error realtivo:

Implementación en Maxima

El primer paso, por tanto es crear la expresión a partir de la función original. Suele haber varias formas de hacerlo, en este caso:

Por tanto, tenemos que:

Dibujemos con Maxima todas las funciones:

f(x):= -x^2+3*x+0.1;
g1(x):=((x^2)/3)-(0.1/3);
g2(x):=(-0.1/(-x+3));
plot2d([f(x),g1(x),g2(x), x],[x,-1,4],[y,-3,5]
, [xlabel, “Inversión en publicidad (miles de euros)”],
     [ylabel, “Beneficios (miles de euros)”],
[style, lines, lines, lines], [color, green, red, blue, black],
     [legend, “f(x)”, “g1(x)”, “g2(x)”, “y=x”]);

b388_4Como se puede apreciar, la recta corta a  y en el valor de x que corresponde a las soluciones de . No vamos a escoger porque tenemos una singularidad en x=3, y por tanto la función no es continua en un entorno cercano a la solución positiva.

Vamos a programar con Maxima el algoritmo, que es bastante sencillo. Para ello, empleamos un bucle con 15 iteraciones, comenzando en 3 como semilla inicial.

g1(i):=(((x[i])^2)/3)-(0.1/3);
x[0]:3;
for i:1 thru 15 do(
x[i]:g1(i-1),print(i,”_Iteración”, x[i-1],g1(i-1)));
datos: makelist([i,”Iteración”,x[i-1],g1(i-1)], i, 1,15);
matriz_resultados:apply(matrix,datos);

Le hemos dicho al programa que nos retorne en el output siguiente:

Como puede observarse, cada imagen de x es el valor de prueba de la siguiente iteración, hasta aproximarnos a la solución. Pero, atención, la solución es negativa, la que no nos interesaba por no tener sentido en nuestro ejemplo.

¿Por qué no converge a la solución positiva pese a que hemos escogido una semilla muy cercana a la solución real?

Si cogemos un intervalo en las proximidades de la solución y lo comparamos con , veremos que no se cumple esa condición. Por ejemplo, en el intervalo [2.5, 3.5] que es el que usamos en el método de la bisección, podemos coger x=3 y comprobar lo que sucede:

En cualquier otro punto de ese intervalo no se cumple el criterio suficiente para la convergencia. No obstante, podría haber convergencia (ya que el criterio no es necesario), pero hemos comprobado que no es así.

El método del punto fijo se va hacia la solución en el eje negativo de las X, donde sí que en un intervalo cercano se cumple esa condición. Por ejemplo, en el intervalo [-1,0].

Aceleración de Aitken

La convergencia se puede acerlerar imlpementando el método de Aitken, donde se construye una nueva sucesión a partir de los 2 primeros términos de la original (la calculada por el método numérico, en este caso del punto fijo):

Podemos programar en Maxima el método de Aitken en conjunción con el de punto fijo, de la siguiente forma:

g1(i):=(((p[i])^2)/3)-(0.1/3);
p[0]:3;
for i:1 thru 15 do(
p[i]:g1(i-1),
p_aitken[i+2]:p[i-1]-(((p[i]-p[i-1])^2)/(p[i+1]-2*p[i]+p[i-1])));
datos: makelist([i,”Iteración”,p[i-1],g1(i-1), p_aitken[i+2]], i, 1,15);
matriz_resultados:apply(matrix,datos);

Al implementar ese código en Maxima, se puede ver que la convergencia es más rápida.

Conclusión

En este post hemos dado unas pinceladas sobre la implementación del método de punto fijo, encontrando que pueden existir problemas de convergencia para determinadas funciones. Asimismo, la aceleraciónde Aitken puede emplearse para intentar aproximarse a la solución con menor coste computacional.

Todos los posts relacionados