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

 

 

 




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




(#387). BÚSQUEDA DE SOLUCIONES (I): MÉTODO DE LA BISECCIÓN

[MONOTEMA] Vamos a comenzar una serie de posts sobre la búsqueda de soluciones  numéricas (también llamadas raíces o ceros)  para una función. El objetivo es que los alumnos de marketing se familiaricen con técnicas matemáticas y de computación que les ayuden a resolver problemas más complejos en el futuro. Para ello vamos a emplear el software Maxima, de libre acceso.

Función de partida

Imaginemos que conocemos la forma de una función que relaciona la inversión en publicidad con los beneficios. En entradas anteriores hemos explicado que esa relación es probablemente no lineal, y que puede adoptar diversas formas. Vamos a partir de la siguiente función:

 Lo primero que vamos a hacer es representar esa función en Maxima:

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)”]);

En esta hipotética situación, el analista de marketing puede computar el máximo y el cero de la función, de la siguiente forma:

derivada: diff(funcion, x);
maximo:solve(derivada,x), numer;
cero: solve(funcion,x),numer;
f_maximo: ev(funcion,x=1.5), numer;
f_cero:ev(funcion,x=3.03297097167559), numer;
plot2d([[discrete, [[1.5,f_maximo], [3.03297097167559, f_cero]]], funcion],[x,0,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]);

De este modo tenemos que el máximo es , y el cero es cuando (hay otra solución al ser una ecuación de segundo grado pero no nos interesa porque se refiere a una inversión en publicidad negativa, lo cual no es posible).

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

b387_3Es evidente que nos interesa conocer el valor de inversión que maximiza los beneficios, pero también nos puede resultar atractivo conocer el punto a partir del cual incrementar más la inversión produce pérdidas. Y este es el objetivo en el que nos vamos a centrar.

Para ello implementar el método de la bisección, que puede consultarse en Burden, Faires & Burden (2017).

Método de la bisección

(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, donde  y son de distinto signo, es decir, .

(3) Descripción rápida: Escoger un intervalo  a priori donde e ir evaluando soluciones   en la mitad de ese intervalo y de los siguientes   hasta que se llegue a un cierto número de iteraciones  o a un valor prefijado de error . De este modo, se trata de ir evaluando la función en los valores de medios de cada intervalo sucesivo (donde la anchura de cada intervalo es la mitad del anterior) hasta aproximarse a la solución. Por el teorema de Bolzano, sabemos que vamos a encontrar esa raíz.

(4) Convergencia: Siempre converge a una solución, aunque su velocidad puede ser lenta.

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

Implementación en Maxima

Vamos a programar con Maxima el algoritmo de la bisección. Para ello, empleamos un bucle con 15 iteraciones. La inspección gráfica nos permite estipular un intervalo inicial [2.5,3.5]:

a[0]:2.5;
b[0]:3.5;
for i:1 thru 15 do(
p[i]:a[i-1]+((b[i-1]-a[i-1])/2),
f[i]: -(p[i])^2+3*(p[i])+0.1,
g[i]: -(a[i-1])^2+3*(a[i-1])+0.1,
if g[i]*f[i]<0
then (b[i]:p[i],
        a[i]:a[i-1])
else (a[i]:p[i],
                b[i]:b[i-1]),
Error_abs[i]:abs(p[i+1]-p[i]),
Error_rel[i]: Error_abs[i]/abs(p[i+1]));

datos: makelist([i,”Iteración”,a[i],b[i],p[i],f[i], Error_abs[i-1], Error_rel[i-1]], i, 1,15);
matriz_resultados:apply(matrix,datos);

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

b387_4_1Como puede apreciarse en las dos primeras columnas de resultados tenemos los extremos de cada intervalo del algoritmo. En la tercera está el valor de p, que siempre es la mitad del valor de esos extremos. la cuarta columna evalúa la función en cada p considerado, mientras que las columnas quinta y sexta corresponden, respectivamente, al error absoluto y al relativo.

Podemos ver aquí como las iteraciones nos acercan a la solución:

matriz_t: transpose(matriz_resultados);

p_list: matriz_t[5];

fp_list: matriz_t[6];

funcion: -x^2+3*x+0.1;
plot2d([[discrete, p_list,fp_list],funcion],[x,0,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]);

b387_5 Y también es interesante comprobar cómo disminuye el error relativo con el número de iteraciones:

Error_rel_list: matriz_t[8];
Iteraciones_list:matriz_t[1];

plot2d([discrete, Iteraciones_list,Error_rel_list],[x,1,15],[y,0,0.09],
     [xlabel, “Iteraciones”],
     [ylabel, “Error relativo”],
[style, lines], [color, red], [legend, false]);

b387_6 Criterio de parada

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.001:

a[0]:2.5;
b[0]:3.5;
for i:1 while Error_rel[i]>0.001
    do(
p[i]:a[i-1]+((b[i-1]-a[i-1])/2),
f[i]: -(p[i])^2+3*(p[i])+0.1,
g[i]: -(a[i-1])^2+3*(a[i-1])+0.1,
if g[i]*f[i]<0
then (b[i]:p[i],
        a[i]:a[i-1])
else (a[i]:p[i],
                b[i]:b[i-1]),
Error_abs[i]:abs(p[i+1]-p[i]),
Error_rel[i]: Error_abs[i]/abs(p[i+1]),
        print(i,p[i], f[i],Error_rel[i]));

Le hemos dicho al programa que fije un error relativo mayor que 0.001, 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.

Alternativamente podríamos conocer de forma conservadora las iteraciones necesarias para una precisión determinada empleado la siguiente fórmula:

Conclusión

En este post hemos explicado con carácter eminentemente práctico el método de la bisección para buscar una aproximación numérica a la solución de una ecuación. Es importante señalar que los errores absoluto y relativo que hemos definido en las programaciones con Maxima se refieren siempre a las estimaciones de la solución en sus respectivas iteraciones, pero no al error absoluto y relativo tomado cuando se conoce valor real de la solución.

En próximos posts, comentaremos otros métodos diferentes.

Todos los posts relacionados




(#373). PREDICCIONES CON BUENAS APROXIMACIONES

[MONOTEMA] En el contexto de gran incertidumbre en las ciencias sociales, los estudiantes de Administración de Empresas se introducen (sólo un poco) en diferentes técnicas de predicción. Sin embargo, es clave entender que, aunque se consigan buenas aproximaciones a las ecuaciones que rigen un modelo de predicción, el resultado final de ésta puede distar de la realidad de manera importante.

En este post, vamos a hacer unas simples simulaciones, con el fin de ir incorporando los conocimientos que se van adquiriendo en la carrera, para evaluar en qué medida una buena aproximación puede ser útil o no.

Las ventas de un producto

Imaginemos que las ventas de un producto a lo largo del próximo año van a ser las siguientes:

Enero 0
Febrero 2
Marzo 6
Abril 12
Mayo 20
Junio 30
Julio 42
Agosto 56
Septiembre 72
Octubre 90
Noviembre 110
Diciembre 132
Total año 572

Podemos representar esas ventas con el siguiente código en una sesión de wxMaxima:

draw2d(fill_color=orange,
yrange=[0,150],
    fill_density=1,
    xlabel=”Meses”,
    ylabel=”Unidades vendidas”,
bars([1,0,1],[2,2,1],[3,6,1],[4,12,1],[5,20,1],[6,30,1],
    [7,42,1],[8,56,1],[9,72,1],[10,90,1],
        [11,110,1],[12,132,1]));

Esas ventas, en realidad, no las sabemos, pero hemos tenido la destreza y el acierto de plantear un modelo matemático que las predice muy bien. Para ello, disponemos de la siguiente función:

donde y son las ventas y x el número de mes (desde 1 hasta 12).

Incorporamos la curva al gráfico anterior:

funcion:x^2-x;
draw2d(fill_color=grey,
yrange=[0,150],
    fill_density=0,
    xlabel=”Meses”,
    ylabel=”Unidades vendidas”,
bars([1,0,1],[2,2,1],[3,6,1],[4,12,1],[5,20,1],[6,30,1],
    [7,42,1],[8,56,1],[9,72,1],[10,90,1],
        [11,110,1],[12,132,1]), color=orange,line_width=4,explicit(funcion,x,0,12));

Y vemos como el área entre los puntos donde está definida la función (0 y 12), es una buena aproximación a las ventas totales:

draw2d(fill_color=grey,
yrange=[0,150],
    fill_density=0,
    xlabel=”Meses”,
    ylabel=”Unidades vendidas”,
bars([1,0,1],[2,2,1],[3,6,1],[4,12,1],[5,20,1],[6,30,1],
    [7,42,1],[8,56,1],[9,72,1],[10,90,1],
        [11,110,1],[12,132,1]), color=orange,fill_color=orange,
filled_func=true,line_width=4,
    explicit(funcion,x,0,12));

El cálculo del área bajo la curva puede hacerse fácilmente en Maxima con este código:

funcion:x^2-x;

integrate(funcion,x,0,12);

Y el valor es de 504. De este modo, hemos conseguido una estimación de las 572 ventas reales, con un error relativo del 11.9%.

Es sencillo ver que el área de cada barra (base x altura) representa las unidades vendidas (la base siempre es 1), y la suma del área de todas las barras es 572 unidades. Nuestra función hace un buen trabajo, pero está definida de manera continua en el intervalo (0,12), y no de manera discreta. He ahí la diferencia, y por eso la integral reporta una aproximación, que sería cada vez mejor en la medida en que el número de puntos discretos se incrementara (los rectángulos fueran cada vez más finos).

Es decir, aunque tenemos una muy buena función que describe los datos empíricos, cometemos errores de predicción. Esa función continua es una aproximación que hemos considerado necesaria, porque en algún momento podríamos tener la intención de predecir las ventas de manera quincenal, por ejemplo, a 15 de octubre. De este modo, el mes sería 10.5, y necesitamos una función continua definida en (0,12) para poder realizar la predicción.

Si queremos predecir lo que se venderá en un mes determinado, podemos emplear la integral definida entre el mes en cuestión y el mes anterior. Por ejemplo, si queremos saber lo que venderemos el mes 6, integraremos entre 5 y 6, y se queremos saber lo que venderemos hasta el mes 6, entonces integraremos entre 0 y 6.

funcion:x^2-x;

integrate(funcion,x,5,6);

integrate(funcion,x,0,6);

Así, para el mes de junio la previsión será (redondeando) de 25 unidades vendidas (30 serían las reales), mientras que en el primer semestre la predicción sería de 54, cuando la realidad sería 70.

Una buena aproximación (de la aproximación) empeora la predicción

Pese a que tenemos una función que se ajusta muy bien a la realidad, las predicciones no son exactas, y como hemos visto tienen errores. ¿Pero qué sucede si no hemos sido tan afortunados de encontrar una función tan buena como la anterior, y sólo es una buena aproximación a ella?.

Podemos simular de nuevo lo que pasaría, empleando el desarrollo de Taylor de la función descrita, en este caso, sólo hasta el primer término. Para ello, elegimos como entorno el punto medio del rango de datos, es decir, el mes 6 (junio). Usaremos el siguiente código:

funcion:x^2-x;

funcion_aproximada:taylor(funcion,x,6,1);

plot2d([funcion,funcion_aproximada],[x,0,12],[y,0,150]);

373_5Como se aprecia en el gráfico anterior, hemos realizado una aproximación lineal a la función cuadrática original. Visualmente se puede asegurar que la predicción será peor, y que habrá una infra estimación severa con esta nueva función lineal:

funcion_aproximada:taylor(funcion,x,6,1);

integrate(funcion_aproximada,x,0,12);

El resultado predicho para el total anual es de 360, muy lejos delos 572 reales, también muy lejos de los 504 de la función no lineal. De este modo, el error relativo sería del 37.1%, ciertamente inasumible.

Conclusión

Tomar decisiones en marketing es siempre complejo, debido a la alta incertidumbre asociada con los fenómenos sociales. Las herramientas matemáticas pueden ayudar a encontrar funciones que nos permitan acercarnos a una predicción aceptable.

Pero, como acabamos de ver, incluso aproximaciones excelentes nos proporcionan errores de cierta entidad. Si, además, no conseguimos llegar a encontrar el modelo adecuado y empleamos aproximaciones lineales puede que los errores sean inaceptables.

Obviamente, lo descrito aquí es sólo una forma de ilustrar este problema. Hay múltiples opciones para intentar encontrar buenos modelos. Sin embargo, la idea subyacente sigue siendo la misma; la complejidad lo “complica” todo.

Todos los posts relacionados




(#362). ¿HORMESIS ENTRE EL VOLUMEN DE ENTRENAMIENTO Y LA HIPERTROFIA?

[MONOTEMA] Emplear ejemplos de otras disciplinas científicas es una herramienta muy interesante para conseguir estimular el aprendizaje de los alumnos que cursan una determinada rama del conocimiento. En mis clases de marketing suelo utilizar, sobre todo, analogías con el ámbito de la salud y el deporte, donde también tengo experiencia como investigador.

Es cierto que, a veces, podemos abusar de ejemplos sacados del mundo del deporte, pensando que son familiares para todos, cuando en realidad siempre hay una parte del alumnado que les cuesta seguirlo, tal y como comenta Andrew Gelman en su blog.

Sin embargo, sigo pensando que las sinergias obtenidas de la unión de varias áreas de la ciencia, junto con el ejercicio de trasladar la explicación de un contexto a otro, son muy estimulantes y motivadoras, y producen un entendimiento global mucho más completo del problema.

En este post, vamos de nuevo a pararnos en un fenómeno no lineal: la hormesis, algo de lo que ya hemos hablado en post anteriores, pero que sigue siendo un concepto de difícil digestión para algunos estudiantes de marketing en el Grado de Administración y Dirección de Empresas. A través de la analogía con lo que sucede en el entrenamiento de hipertrofia muscular, vamos a intentar mejorar la comprensión de una de las más relevantes características de los sistemas complejos.

La explicación de Bayesian Bodybuilding

Menno Henselmans publicó recientemente este interesante artículo en su blog, donde aludía a resultados de sus propias investigaciones y también a otras publicaciones para realizar un dibujo estimulante sobre la relación entre el volumen de entrenamiento con pesas y la hipertrofia muscular.

Henselmans analiza la relación entre ambas variables tras consultar decenas de estudios publicados hasta la fecha. En esos estudios, se reporta el tamaño de efecto (magnitud de la diferencia sustantiva) y el volumen referido al número de series por grupo muscular ejecutadas en una semana de entreno.

Recordemos que el tamaño de efecto complementa al resultado de los clásicos test estadísticos. Por ejemplo, con muestras muy grandes las diferencias entre el grupo de experimental y de control pueden ser significativas estadísticamente, pero el tamaño de efecto (magnitud real de la diferencia) puede ser de escasa relevancia. Por contra, resultados estadísticamente no significativos por falta de potencia (como puede ocurrir con muestras pequeñas) pueden enmascarar tamaños de efecto importantes.

Por ejemplo, con muestras suficientemente grandes se podría decir que la estatura de una población de niños es mayor que otra cuando la diferencia es de 150 cm a 149.5 cm. Sería una diferencia estadísticamente significativa pero sustantivamente poco relevante; la magnitud del efecto es muy pequeña. Sin embargo, otro estudio con muestras pequeñas de otras poblaciones de niños podría reportar que ambas poblaciones tienen la misma estatura aunque la diferencia real fuera de 150 cm frente a 140 cm. Aquí la magnitud del efecto sería considerable, pero la baja potencia lo haría estadísticamente no significativa. Lo ideal, por tanto, es diseñar estudios que detecten estadísticamente magnitud de efectos, y los valoren según su importancia sustantiva.

La conclusión a las que llega Henselmans es que un incremento del volumen de entrenamiento se asocia linealmente con el incremento de la hipertrofia muscular. Sin embargo, en subpoblaciones distintas, como en deportistas inexpertos frente a deportistas experimentados, esa relación difiere, y muestra un carácter no lineal con rendimientos decrecientes. Gráficamente se entiende mejor:

Henselmans-meta-analysis-training-volume-untrainedHenselmans-meta-analysis-training-volume-trainedEs más, un entrenamiento de alto volumen puede llevar al sobre entrenamiento, como bien indica citando este estudio. De este modo, hay un patrón de respuesta a la dosis por el cual las máximas ganancias de hipertrofia se consiguen con un volumen medio de entrenamiento, pero no con los extremos testados. Esto puede ocurrir también en el entrenamiento de fuerza (no siempre hipertrofia y fuerza van de la mano).

Henselmans concluye el artículo con la identificación de diferentes variables que pueden afectar a la relación entre el volumen y la hipertrofia (ver siguiente figura).

Training-volume-vs-recovery-capacity

Esas variables son moderadoras en el sentido de que cambian el tamaño de efecto de la asociación entre ambas. De este modo, se complica mucho más el entendimiento de esa relación, porque en función de una o varias de esas variables, la curva de hormesis puede variar. En suma, cada individuo necesita un estudio pormenorizado de esas variables para diseñar una rutina de entrenamiento óptima. Sin embargo, para dar reglas generales hemos de tomar la distribución de tamaños de efecto, y dar una curva de hormesis promedio, pero que no tiene por qué funcionar con la situación de un individuo particular. Así, volvemos a darnos cuenta de cómo se pueden producir cambios importantes cuando se baja del “universo promedio” a los casos particulares. Y esto es el “abc” de la interpretación de los análisis estadísticos en ciencias sociales.

Matizaciones al análisis de Bayesian Bodybuilding

Pero lo más atractivo de todo esto es que quizá podríamos interpretar los resultados de los análisis de Henselmans de manera diferente.

Para ello, hemos pasado los datos de las gráficas a una hoja de cálculo. Hemos perdido precisión, sí, y no conocemos las varianzas de cada uno de los tamaños de efecto. Pero, aunque sea de manera aproximada (y no podemos ir más allá por esta limitación), podemos ofrecer un dibujo ligeramente divergente.

En cuanto a los deportistas inexpertos, podemos hacer un primer análisis del modelo lineal, y obtenemos estos resultados:

hiper1

Claramente se aprecia que no hay evidencias para un efecto lineal.  La regresión lineal simple por mínimos cuadrados ordinarios (con homocedasticidad) nos dice perfectamente que no existe asociación entre el volumen de entrenamiento (volumen) y la hipertrofia (efecto). Por tanto, contradecimos la conclusión de Henselmans.

Si añadimos un término cuadrático buscando la no linealidad, tampoco obtenemos significatividad:

hiper2

bayesian1

Y lo mismo ocurre con un análisis similar (no igual) realizado con polinomios fraccionales:

hiper3 bayesian2

No obstante, tanto el ajuste cuadrático como el de polinomios fraccionales nos muestran una posible relación de hormesis, donde la dosis justa de entrenamiento que maximiza la hipertrofia se da en los valores centrales.

Pero, aunque la relación sea no significativa, ¿se puede hablar de hormesis? No del todo, hay que ser prudente, porque tenemos poca muestra (39 observaciones). Quizá con una muestra más grande (más estudios publicados) se podría llegar a la significatividad, pero no lo sabemos. Siendo conservadores, diríamos que para deportistas inexpertos no hay una relación clara entre volumen de entrenamiento e hipertrofia, y que, por ejemplo, entrenar 10 series por músculo y semana ofrece las mismas ventajas que cualquier otro entrenamiento más exigente.

Pero, de nuevo, esto son recomendaciones generales; la heterogeneidad de los estudios (como puede apreciarse en los diagramas de dispersión mostrados por Henselmans) es ciertamente importante.

¿Y qué sucede para deportistas con experiencia?. Pues podríamos matizar las conclusiones de Henselmans en función de los datos (de nuevo insisto que no es una crítica  Henselmans, ya que no dispongo de datos originales), sino un ejercicio sencillo de interpretación, con fines docentes, nada más.

En este caso el efecto sobre la hipertrofia lo hemos llamado “e2” y el volumen de entrenamiento “v2”. Vemos que estamos cerca de una relación lineal. Como hay heterocedasticidad, al computar los errores estándar robustos la significatividad al 95% se pierde ligeramente.

hiper4

Si intentamos hacer un ajuste cuadrático y con polinomios fraccionales, de nuevo los resultados son no significativos. Sin embargo, al visionar las gráficas de ajuste vemos que existe una tendencia creciente, y no esos rendimientos decrecientes de los que hablaba Henselmans.

bayesian3 bayesian4

De este modo, conseguir ganancias de hipertrofia es más complicado para deportistas entrenados (algo lógico, porque se parte de un nivel inicial de hipertrofia superior), pero aquí la tendencia es lineal, y probablemente con más estudios podríamos conseguir un nivel de asociación significativo, es decir, a mayor volumen mayores ganancias de músculo. No obstante, y de nuevo siendo prudentes, con esos 44 estudios que tenemos no podemos concluir más que existe una evidencia limitada de que así ocurre.

Conclusión

Menno Henselmans, de Bayesian Bodybuilding, ha realizado un interesante análisis sobre la relación entre el volumen de entrenamiento y la ganancia muscular. Sus interpretaciones sobre linealidad, hormesis y rendimientos decrecientes pueden, quizá, ser matizadas, aunque como hemos indicado nuestro análisis es sólo aproximado, ya que no disponemos de los datos originales y ha habido contaminación por pasar los datos de los gráficos de dispersión a una hoja de cálculo. No obstante, y teniendo en cuenta esta limitación, hemos mostrado que las relaciones mostradas en gráficos no tienen por qué ser significativas a nivel estadístico, y que esas tendencias pueden ser reinterpretadas en función del tipo de ajuste elegido, es decir, el tipo de modelo empleado.

Sin embargo, esas tendencias son en sí valiosas, en el sentido en que es probable que muestren un efecto escondido por un tamaño pequeño de muestra. Si aceptamos esta hipótesis, podemos decir que existe una relación de hormesis entre el volumen y la hipertrofia en deportistas no entrenados, pero esa relación desaparece en individuos más expertos, donde la asociación es creciente y lineal. 

Como puede apreciarse, esto cambia de forma importante las recomendaciones generales sobre entrenamiento de hipertrofia en función de las características de cada individuo, y es ahí donde de nuevo Henselmans acierta en indicarnos los múltiples factores que pueden entrar en juego.

El traslado de este ejemplo al ámbito del marketing o de la administración de empresas, es sencillo, como vimos en posts anteriores. Es un buen ejercicio para los estudiantes el pensar en relaciones de variables en el ámbito de las ciencias sociales donde se pueden dar relaciones análogas a las descritas en este post, así como reflexionar sobre la complejidad de todos los fenómenos que nos rodean y que pretendemos entender, y la diferencia que existe entre encontrar patrones generales y la realidad a nivel de caso individual.

Todos los posts relacionados




(#346). MODELO NO LINEAL EN LAS VARIABLES PERO LINEAL EN LOS PARÁMETROS CON STATA Y MAXIMA

[MONOTEMA] Seguimos creando material auxiliar para comprender mejor las herramientas de toma de decisiones en marketing. En este post vamos a comentar uno de los aspectos que quizá cuesta un poco más de trabajo entender en las primeras clases de marketing, y es el de emplear el mismo modelo de regresión lineal cuando el modelo teórico es no lineal.

Esta aparente contradicción llama la atención de los estudiantes, y más cuando estamos incidiendo en el hecho de que la no linealidad está presente en la mayoría de los fenómenos que queremos estudiar. Sin embargo, la clave está en entender que la estimación de mínimos cuadrados ordinarios que solemos emplear, también se puede aplicar a modelos no lineales. Y es así, porque la no linealidad está en las variables pero no en los parámetros estimados. Este hecho es fundamental para comprender que con la misma herramienta con la que intentamos hacer predicciones en modelos lineales, también podemos emplearla en modelos no lineales.

Vamos a ilustrarlo con un ejemplo, usando códigos de Stata y MAXIMA, para facilitar que cualquier alumno pueda reproducirlos y hacer probaturas a su antojo.

Generación de la curva no lineal

Elegimos una curva normal de media=5 y desviación típica=1. Por tanto, transformamos la función de densidad normal:

en la siguiente función:

que tiene la característica forma de campana:

curva1Para generar el gráfico anterior, hemos empleado el siguiente código en wxMaxima:

funcion1:(1/(2*%pi)^0.5)*%e^(-0.5*(x-5)^2)$
plot2d(funcion1, [x,0,10],[style, lines],
[color, green], [legend, false],
[xlabel,”x”], [ylabel, “y”]);

Es evidente que se trata de una curva no lineal, pero podemos tratar de realizar predicciones muy ajustadas si planteamos el modelo correcto y lo estimamos por mínimos cuadrados ordinarios (con linealidad en los parámetros).

Ajuste del modelo con Stata

Vamos a emplear ahora uno de los programas más potentes de análisis de datos, y de los más indicados para los profesionales del marketing: Stata.

Para ello, abrimos un archivo en blanco de Stata, y utilizamos el interfaz de comandos que el programa proporciona, llamado “Do-file Editor”, que es simplemente un editor de texto donde se puede escribir un código (secuencia de instrucciones) para ejecutar.

set obs 1000
generate t = _n
drawnorm x, n(1000) means(5) sds(1)
gen y=(1/(2*3.1415)^0.5)*exp(-0.5*((x-5)^2))
gen lnx=ln(x)
gen lny=ln(y)
gen x_cuad=x*x
regress lny x x_cuad
gen pred=-13.4+(5*x)-(0.5*x_cuad)
gen transf_pred=exp(pred)
gen dif= y- transf_pred
tsset t
twoway (tsline dif, lcolor(red) lwidth(vvvthin) xtitle(t) ytitle(Diferencia entre el valor real y el predicho))

Explicamos de manera muy sencilla lo que hemos programado. En primer lugar le hemos dicho al programa que genere una muestra de 1000 observaciones, y una variable llamada “t” que es el identificador de cada una de ellas.

Después, le pedimos a Stata que genere para esos 1000 casos una variable Normal llamada “x” con media 5 y desviación típica 1, tal y como hemos mencionado en el apartado anterior.

Ahora debemos asignar una probabilidad a cada uno de esos valores de la variable x, y lo hacemos simplemente espeficiando una variable y que es la función de densidad Normal con media 5 y desviación típica 1.

Para ilustrar que a partir de los valores de x podemos explicar y/o predecir y, tenemos que especificar un modelo que relacione ambas variables. Ese modelo, es obvio que es el regido por la función de densidad Normal, pero ahí no tenemos nada claro cómo encajar nuestro marco de regresión lineal.

Lo que hacemos entonces es tomar logaritmos neperianos a ambos lados de la igualdad:

De este modo:

Si lo reordenamos:

Con esto ya tenemos una disposición muy interesante para realizar la ilustración que pretendíamos.

Si imaginamos ahora que no conocíamos esta especificación exacta de la realidad, pero tenemos el suficiente conocimiento teórico para postular un modelo que trate de describir esos datos, podemos entonces plantear el siguiente modelo de regresión no lineal en las variables, pero lineal en los parámetros:

siendo u el error aleatorio, ruido blanco que desconocemos y que no covaría con x. Entonces planteamos el modelo de regresión con Stata (regress lny x x_cuad), y obtenemos los parámetros estimados correspondientes:

salidastata

Claro, conseguimos un ajuste prácticamente perfecto, porque hemos hecho un poco de “trampa” a la hora de plantear el modelo (sabíamos el modelo subyacente antes de plantearlo).

Para ilustrar mejor cómo predice el modelo (y que se vean unas mínimas diferencias con las observaciones reales) cogemos sólo el primer decimal de la constante, es decir:

Y entonces generamos el valor predicho por el modelo (gen pred=-13.4+(5*x)-(0.5*x_cuad)). Por último, transformamos el valor predicho con el antilogaritmo neperiano, para compararlo con nuestras observaciones reales. Y generamos el gráfico con Stata.

Graph1

Como puede observarse, los errores son extremadamente pequeños, del orden de milésimas de unidad. Eso quiere decir, que nuestras predicciones son realmente buenas.

Ahora podemos ir a wxMaxima y comparar la curva original de datos observados (funcion1), con la curva de valores predichos (función2)

funcion1:(1/(2*%pi)^0.5)*%e^(-0.5*(x-5)^2)$
funcion2: %e^(-13.4+5*x-0.5*x^2)$
plot2d([funcion1,funcion2], [x, 0, 10],[style, lines, lines],
[color, green, orange], [legend, false],
[xlabel,”x”], [ylabel, “y”]);

Y entonces obtenemos la curva original (en verde) frente a la predicha (en naranja), ambas casi idénticas:

curva2

Conclusión

El ajuste por mínimos cuadrados ordinarios es una herramienta muy potente para intentar conocer la realidad que subyace a los datos, independientemente de que la relación entre las variables sea lineal o no lineal.

Sin embargo, y como veremos en posteriores posts, hay que conocer adecuadamente las asunciones que se vinculan a cada modelo propuesto, y testarlas vis a vis con los datos empíricos. Recordemos que las asunciones están dentro de la propia especificación del modelo, por lo que son inherentes a él.

En cualquier caso, y a modo meramente cualitativo, los estudiantes de marketing deben darse cuenta de la importancia que tiene dominar la estadística, para ayudar en la toma de decisiones.

Todos los posts relacionados




(#335). EXPLORANDO EL COMPORTAMIENTO CAÓTICO CON MAXIMA

[MONOTEMA] Continuamos analizando las relaciones no lineales entre variables en esta serie de posts para facilitar la comprensión acerca de la complejidad de los fenómenos sociales.

Hoy vamos comentar (sin entrar en detalles profundos) el comportamiento caótico de un sistema. La idea es ver que, incluso en sistemas muy sencillos, con ecuaciones simples, hacer predicciones puede ser una tarea muy complicada.

Vamos a servirnos parcialmente del artículo de Morante & Vallejo (2013), que de manera bastante amigable explican cómo poder simular este tipo de comportamientos en Maxima.

La ecuación lógistica

Emplearemos la ecuación logística (o aplicación logística) para ilustrar ese comportamiento:

donde x es una variable que está entre 0 y 1 y r está entre 0 y 4. En el contexto de los estudios poblacionales, x representa el porcentaje de población sobre el máximo posible, y r una tasa entre la reproducción y la mortandad. Así, si r<1, la mortandad es mayor que la reproducción, y el sistema acabará evolucionando hacia un valor de x=0 en un determinado tiempo t.

Evolución del sistema en función de r

En primer lugar, vamos a realizar varias simulaciones sobre la evolución temporal del sistema en función de varios valores de r: 0.8, 1.8, 2.8, 3.8. Para ello tomamos como valor fijo x=0.9, es decir, el valor de partida es que la población está casi al máximo de lo posible.

Para ello, escribimos el siguiente código en una sesión con wxMaxima:

F(r,x):=r*x*(1-x)$    /* Ecuación logística*/
load(dynamics)$    /* Carga del paquete para la simulación*/
x_0:0.9$    /* Condición inicial del sistema*/
dominio_x:[y,0,1]$   /* Rango del eje y*/
tiempo:30$    /* Intervalos de tiempo discreto */
r:0.8$   /* Valor asignado a r */
evolution(F(r,x),x_0,tiempo,dominio_x,[color, red],[xlabel, “tiempo”],[ylabel, “x_t_+_1”], [legend, “r=0.8”]);

r_08

Como era previsible, la población decae hasta que desaparece. Ahora, podemos probar el mismo código de Maxima pero con los valores r=1.8 y r=2.8:

r_18
r_28

En ambos casos se tiende a un punto de equilibiro en (r-1)/r. Sin embargo, cuando r=2.8, la respuesta del sistema comienza a oscilar entre ese punto, suavizándose progresivamente.

Pero cuando r=3.8, el resultado es el siguiente:

r_38

Como puede apreciarse, el sistema se vuelve mucho más oscilante, y no se adivina ninguna tendencia a la estabilidad.

Programación “a mano” de la evolución del sistema

Hemos empleado la función “evolution” de Maxima, que facilita mucho el trabajo. Sin embargo, podemos programar “a mano” el comportamiento del sistema, y así tratar de entenderlo mejor.

Para ello, escribimos el siguiente código en una sesión con wxMaxima:

x[t]:=r*(x[t-1])*(1-(x[t-1]))$        /* ecuación logistica*/
r:2.8$      /* fijamos el valor de r*/
x[1]:0.9$      /* fijamos la condición inicial para t=1*/
datos: makelist([t,x[t]],t,1,30)$     /* hacemos una lista con el resultado de 30 pasos de tiempo*/
matriz:apply(matrix,datos)$    /* convertimos esa lista en una matriz para poder emplearla luego*/
matriz_t: transpose(matriz)$  /* transponemos la matriz anterior para poder usarla en las gráficas*/
t_list: matriz_t[1]$  /* esa matriz traspuesta tiene 2 filas y 30 columnas. La primera fila es la correspondiente a los pasos de tiempo, y es la extraemos aquí */
x_list: matriz_t[2]$  /* extraemos la segunda fila de la matriz, que es la correspondiente al resultado de la función en cada paso de tiempo */
g[t]:=(r-1)/r$  /* creamos una nueva función que representa el punto de equilibrio del sistema, que obviamente depende del valor de r que hayamos fijado. Este valor es una constante, ya que no depende de t */
punto_eq:makelist([t,g[t]],t,1,30)$  /* seguimos el mismo proceso anterior y hacemos la lista de 30 pasos de tiempo */
matriz2:apply(matrix,punto_eq)$  /* convertimos esa lista en una matriz para poder emplearla luego*/
matriz2_t:transpose(matriz2)$  /* transponemos la matriz anterior para poder usarla en las gráficas*/
r_list:matriz2_t[2]$ /* extraemos la segunda fila de la matriz, que es la correspondiente al resultado de la función en cada paso de tiempo */
plot2d([[discrete, t_list,x_list],
[discrete,t_list, r_list]], [style,points,lines],
[xlabel, “tiempo”],[ylabel, “x_t_+_1”],
[x,0,30],[y,0,1],[color,red,blue],
[legend, “Evolución del sistema para r=2.8”, “Punto de equilibrio”]);

Lo que hemos hecho es simplemente reescribir la ecuación logística de esta forma:

Y además hemos añadido la línea que representa el punto de equilibrio. Cualquier lector puede hacer probaturas con otros valores de r y de condición inicial para ver gráficamente cómo evoluciona el sistema.

Dependencia sensible

Una de las características de los procesos caóticos (probablemente la más conocida, y quizá la más inquietante) es la extrema dependencia de las condiciones iniciales. Así, cambios ínfimos en el valor de la variable que inicia el ciclo de evolución hace evolucionar al sistema de forma totalmente diferente.

Para realizar una prueba visual, hemos tomado dos valores de inicio: 0.900 y 0.899, es decir, una milésima de diferencia entre ambos valores.

Para ello, escribimos el siguiente código en una sesión con wxMaxima:

x[t]:=r*(x[t-1])*(1-(x[t-1]))$
r:3.8$
x[1]:0.900$
datos: makelist([t,x[t]],t,1,30)$
matriz:apply(matrix,datos)$
matriz_t: transpose(matriz)$
t_list: matriz_t[1]$
x_list: matriz_t[2]$
g[t]:=r*(g[t-1])*(1-(g[t-1]))$
g[1]:0.899$
datos2:makelist([t,g[t]],t,1,30)$
matriz2:apply(matrix,datos2)$
matriz2_t:transpose(matriz2)$
g_list:matriz2_t[2]$
plot2d([[discrete, t_list,x_list],
[discrete,t_list, g_list]], [style,lines,lines],
[xlabel, “tiempo”],[ylabel, “x_t_+_1”],
[x,0,30],[y,0,1],[color,red,blue],
[legend, “Condición inicial=0.900”, “Condición inicial=0.899”]);

El código simplemente refleja la programación de dos funcions “x[t]” y “g[t]”, ambas idénticas, con la única diferencia de esa milésima en las condiciones iniciales. Los resultados para 30 pasos temporales se muestran en el siguiente gráfico:

caosComo se puede apreciar, ambos sistemas comienzan de manera muy parecida su evolución, pero a partir de 10 pasos temporales empiezan a diverger hasta que el comportamiento de uno hace imposible la predicción del comportamiento del otro.

Conclusión

Los sistemas caóticos pueden tener ecuaciones muy simples, pero pese a esa simplicidad los hacen puntualmente impredecibles, aunque para determinados sistemas las órbitas pueden encontrarse en un espacio delimitado, llamado atractor.

Para las ciencias sociales, y especialmente para el marketing, lo mostrado en este post debe hacer reflexionar a los estudiantes sobre qué es la impredecibilidad y cómo entenderla y afrontarla.

Todos los posts relacionados