(#420). LA PARADOJA DE CHEVALIER DE MERE Y LOS ESTUDIANTES DE MARKETING

[MONOTEMA] Uno de los objetivos fundamentales de mi labor como profesor de marketing es hacer ver a los estudiantes que necesitan una formación científica muy alta para entender mejor los fenómenos sociales a los que se enfrentan. Por ello, hacemos énfasis en la importancia de las matemáticas, la lógica y la heurística para la toma de decisiones, pero con especial interés en estimular el pensamiento y la formación en estadística.

La paradoja de Chevalier de Mere es una buena forma de mostrar cómo los razonamientos aparentemente lógicos y sencillos pueden ser erróneos, con lo que se necesita de un cierto dominio de herramientas matemáticas para llegar al resultado correcto.

La paradoja

Aris Spanos, uno de los grandes maestros en econometría de nuestros días, en su libro Probability Theory and Statistical Inference, la describe de la siguiente forma:

En el siglo XVII, Pascal envió una carta a Fermat en relación a un problema que le había planteado un noble (y experto jugador) llamado Chevalier de Mere.

De Mere observó la siguiente regularidad derivada de los datos empíricos:

  • La probabiliad de obtener al menos un 6 en 4 lanzamientos de un dado es mayor de 0.5
  • La probabilidad de obtener al menos un 6 doble (un 12) en 24 lanzamientos de dos datos es menor que 0.5

De Mere se preguntaba cómo eso era posible si, siguiendo un razonamiento aparentemente lógico, ambas probabilidades debían ser iguales, es decir, conseguir un 6 sobre 4 lanzamientos de un dado debería ser lo mismo que conseguir un doble 6 en 24 lanzamientos de 2 dados, ya que 4/6 es igual a 24/36.

La solución

Aris Spanos, sin embargo, argumenta que no existe tal paradoja, porque cuando uno analiza los datos con suficiente rigor, las probabilidades difieren.

  • La probabilidad de obtener un doble 6 es:

  • La probabilidad de obtener un doble 6 en n lanzamientos es:

 

  • La probabilidad de no obtener un doble 6 en n lanzamientos es:

 

  • La probabilidad de obtener al menos un doble 6 en n lanzamientos es: 
  • Para  

 

En el caso de un dado:

  • La probabilidad de obtener un  6 es:

 

  • La probabilidad de obtener un 6 en m lanzamientos es: 

  • La probabilidad de no obtener un  6 en m lanzamientos es:

 

  • La probabilidad de obtener al menos un 6 en m lanzamientos es: 
  • Para 

 

Por tanto, el resultado empírico de de Mere era correcto, pero su razonamiento lógico fallaba.

Una explicación más clara

Podemos obtener una explicación más clara sobre la solución del problema si consultamos el (recomendable) artículo de  Basulto y Camúñez (2007).

No obstante, podemos simplemente proponer un razonamiento que nos lleve a concluir que ambas probabilidades son diferentes.

Podemos reordenar las ecuaciones anteriores de la siguiente forma:

Si llamamos:

Entonces:

Tomando logaritmos neperianos a ambos lados:

Si ahora suponemos que , es decir que ambas probabilidades son iguales, entonces:

De este modo:

Así, la única manera en que las probabilidades pueden ser iguales es cuando se cumple la relación anterior, es decir, que la razón del número de ensayos sea igual a la inversa de la razón de los logaritmos neperianos de las probabiliades de no obtener el número deseado.

Para el caso del ejemplo de Spanos:

lo que es, obviamente, imposible, y por tanto ambas probabilidades no pueden ser iguales.

Para que ambas fueran iguales, entonces para el caso de dos dados la probabilidad de no obtener un doble 6 debería ser 34.926/36 en lugar de 35/36.

También se puede ver de otro modo, y es que si los ensayos son los mismos, por ejemplo 4, entonces:

que obviamente no es posible.

Un gráfico que puede ayudar a comprender lo que está sucediendo es el que relaciona la probabilidad de no obtener el número deseado cuando se lanzan k dados:

Como se aprecia, el incremento de probabilidad es no lineal (en realidad la figura debería ser dibujada con puntos discretos para k=1, k=2, etc, pero se ha dibujado contínua para que se aprecie mejor la no linealidad).

Conclusión

La paradoja de de Mere no es en realidad una paradoja, sino una muestra de que razonamientos aparentemente sencillos y lógicos no dan el resultado correcto. Para llegar a resolver este tipo de problemas, necesitamos herramientas matemáticas que nos ayuden a no equivocarnos (o a hacerlo lo menos posible).

Los estudiantes de marketing deben apresurarse a interesarse por este sistema de pensamiento para tomar decisiones, sobre todo cuando han de analizar datos empíricos.

Todos los posts relacionados




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




(#385). LAS IMITACIONES PERJUDICAN A LA MARCA ORIGINAL POR SOBRE EXPOSICIÓN MÁS QUE POR SUSTITUCIÓN

[REVISIÓN DE ARTÍCULO] En este estudio publicado en el International Journal of Research in Marketing, los autores analizan un interesante dilema sobre los efectos de las imitaciones de diseños originales en la industria de la moda.

Es una prácticamente común en esta industria que otras marcas y cadenas de distribuidores copien diseños originales de marcas más exclusivas, y luego los vendan a un precio significativamente inferior.

Un ejemplo es el que se muestra a continuación, donde Forever 21 imitó los diseños de Trovata, algo que incluso se llevó a juicio. Forever 21 ha recibido más de 50 demandas de diseñadores por copia de sus creaciones, las cuales, como comentan los autores, la mayoría no han prosperado.

b385_2

En Estados Unidos no existe el nivel de protección a los diseños originales, como sí ocurre en la Unión Europea. Uno de los motivos a esa desprotección es que algunos stakeholders tienen la creencia de que las imitaciones favorecen a la marca copiada, ya que se produce un proceso de aceleración en la difusión del diseño, que actúa como una plataforma publicitaria que incentiva la demanda también del original.

b385_3Por contra, está claro que también se produce un efecto de sustitución, en el que potenciales compradores del producto original se decantan por la imitación. El debate está en que para algunos el proceso de aceleración es más fuerte que el de sustitución (la marca original se ve beneficiada), mientras que para otros el de sutitución es más fuerte que el de aceleración (la marca original se ve perjudicada).

Sin embargo, los autores proponen que, además, existe un efecto de sobre exposición, por el cual cuando se produce la imitación, el crecimiento de las ventas hace que ese diseño sea ubícuo, y en un mercado como el de la moda, eso es un arma de doble filo, porque puede hacer que ese elemento dfierenciador se pierda, y por tanto a largo plazo bajen ostensiblemente las ventas de ese producto (al menos en el segmento de los que buscan exclusividad).

Aceleración, sustitución y sobre exposición son fuerzas que se relacionan de manera dinámica, y el objetivo de esta investigación es tratar de esclarecer cuál será su efecto en los resultados en el largo plazo, es decir, si la resultante beneficiará o no a la marca original.

Como los autores bien indican, las imitaciones son diferentes a las copias prácticamente exactas de productos (productos falsos). En las imitaciones, el diseño se vende con otra marca, mientras que en las falsificaciones el producto se vende como si fuera el original cuando no lo es realmente. En la Unión Europea existe protección no sólo contra las falsificaciones, sino también contra las imitaciones (por periodo de 3 años para los diseños no registrados y por periodo de 10 años para los registrados).

Las cifras sobre este tipo de fraude son demoledoras. Según datos de 2012-2013, el mercado de imitaciones sólo en Estados Unidos es de unos $11300 millones, mientras que el de las imitaciones en la Unión Europea es de $17500 millones, cifras que incluso los autores piensan que son demasiado conservadoras.

Las diferencias en precio entre los originales y las imitaciones puede variar en rango. Los autores, en la muestra que consideran encuentran un promedio de un ratio de 40.5 para los bolsos, de 5.1 para calzado y de 5.0 para la ropa. Es decir, aproximadamente un bolso de imitación costaría 40 veces menos que el original, mientras que calzado o ropa costarían 5 veces menos.

Metodología

Los autores emplean un modelo dinámico de simulación, donde ajustan varios parámetros tras recopilar datos sobre precios y ventas de 20 productos imitados (bolsos, calzado y ropa).

Resultados e implicaciones

Los autores encuentran que, a largo plazo, los efectos de la sobre exposición hacen que la marca original se vea perjudicada, lo que lleva a concluir a los autores que se hace necesaria una protección legal en Estados Unidos que defienda las creaciones originales.

Limitaciones/Comentarios

El problema de este artículo es el de la no consideración de variables y supuestos que afectan a los modelos de simulación estimados. Por ello, no he puesto especial énfasis en la metodología y resultados, que son dependientes de factores no tenidos en cuenta, como las campañas de publicidad, los acuerdos con celebridades, la estacionalidad de las venteas, etc. Los autores reconocen esas limitaciones, ciertamente.

Por tanto, para mí lo fundamental de este estudio es la exposición inicial sobre los posibles efectos de las imitaciones, ya que muestran la extrema complejidad de dilucidar los efectos que pueden producir sobre la marca copiada, algo que, si no se dispone de una gran diversdad de datos es tarea prácticamente imposible analizar empíricamente.

No obstante, esta investigación muestra lateralmente dos problemas quizá con mayor trasfondo: (1) la desvergüenza e impunidad de ciertas marcas que basan parte de su negocio en la explotación de imitaciones; (2) el papel del consumidor desde el punto de vista ético cuando compra imitaciones a sabiendas y, por supuesto, cuando adquiere falsificaciones.

LEE EL ARTÍCULO ORIGINAL AQUÍ:

Appel, G. et al. (2018). On the monetary impact of fashion design piracy. International Journal of Research in Marketing, doi: 10.1016/j.ijresmar.2018.08.003

Indicadores de calidad de la revista*

 

Impact Factor (2017)

Cuartil

Categoría

Thomson-Reuters (JCR)

2.593

Q1

BUSINESS

Scimago (SJR)

2.528

Q1 MARKETING

* Es simplemente un indicador aproximado para valorar la calidad de la publicación

Todos los posts relacionados




(#382). HORMESIS EN EL FACE VALUE DE CUPONES DE DESCUENTO

[REVISIÓN DE ARTÍCULO] En este estudio publicado en el Journal of Marketing, los autores analizan la relación existente entre el valor monetario de un cupón de descuento y el consiguiente gasto en uno de los posibles productos en los que puede emplearse.

El cupón de descuento es una clásica herramienta de promoción de ventas, y se puede diseñar de tal forma que se empleee con varios productos de una misma marca, independientemente del precio que tengan estos. Por ejemplo, un cupón de descuento de $50 para un producto Apple que permita emplearlo en cualquiera de sus dispositivos a la venta.

Lo que se plantean los autores es si el “face value” del cupón, es decir, si la magnitud del descuento tiene una relación lineal con el importe del producto que se compra y donde se hace efectivo ese descuento. Aunque una relación lineal sería intuitivamente esperable, los autores encuentran una asociación en forma de U invertida, en un claro patrón de hormesis.

Los autores argumentan que existe evidencia en la literatura de marketing de que pueden ocurrir casos contrapuestos en la mente del consumidor. En primer lugar, su presupuesto mental se puede incrementar, con lo que estaría dispuesto a gastar más dinero del inicialmente planeado. Sin embargo, en segundo lugar, el consumidor puede centrarse en la comparativa de ahorro entre productos de diverso precio dada la cantidad fija de descuento. Por tanto, en el primer caso existiría una fuerza que dirigiría al consumidor a gastar más, y en el segundo una fuerza que lo llevaría a gastar menos. Para explicar el resultado de esa tensión, los autores proponen que debe haber un valor umbral a partir del cual prevalece una fuerza sobre otra, es decir, por debajo de ese valor se tiende a incrementar el gasto, y por encima a reducirlo.

Metodología y resultados

Los autores realizaron 5 estudios diferentes para probar su hipótesis de asociación en forma de U invertida, así como para investigar los factores moderadores en esa relación.

El primero de ellos fue un estudio de campo, mientras que el resto fueron experimentos realizados en un entorno controlado.

Para el estudio de campo, los autores recogieron datos de 106 restaurantes en China entre mayo de 2005 y marzo de 2008, con un total de 26660 consumidores registrados. Los autores dividieron a los restaurantes en función del gasto medio por persona que los clientes habían reportado, es decir, una forma de diferenciar aquellos más caros de los más baratos. Los resultados se muestran en la siguiente figura:

b382_2De este modo, la hipótesis lineal se cumple en restaurantes más baratos, pero el efecto de hormesis se presenta en los más caros.

El segundo estudio trató de considear la propensión al ahorro de los consumidores. Como era de esperar, los que tenían más orientación al ahorro presentaban una asociación no lineal con la probabilidad de elegir la opción más cara (ver figura):

b382_3No obstante, puede haber una contradicción entre los resultados de estos dos primeros estudios, porque en el primero de ellos la curva es no lineal en situación de precios caros, mientras que en el segundo lo es para aquellos que más desean ahorrar (probablemente son más sensibles al precio).

En el tercer estudio los autores analizaron el rol de una sobrecarga de información (cuando se provee de mucha información sobre un producto hasta el punto de que es difícil procesarla en su totalidad), en el contexto de elección en un restaurante. Como se aprecia en la siguiente figura, la hormesis aparece en la situación en que hay menos sobrecarga de información, por lo que los consumidores tienen una mayor capacidad para no distraerse y realizar comparaciones entre el ahorro a la hora de elegir diferentes opciones de producto.

b382_4

En el cuarto estudio, los autores mostraron que la tendencia a realizar comparaciones modera la relación entre la magnitud del cupón y el precio del producto elegido, tal y como se aprecia a continuación:

b382_5En el quinto estudio, finalmente, los autores les pidieron a los participantes que explícitamente indicaran su preferencia por un producto específico de la gama ofrecida sobre la que gastar el cupón. Para aquellos con una débil preferencia inicial, la curva vuelve a ser no lineal:

b382_6

Limitaciones/Comentarios

La asociación entre la magnitud de un cupón y el gasto subsiguiente cuando se puede elegir entre diferentes productos de la misma marca es compleja, ya que diversos factores moderadores influyen en que se produzca un efecto de hormesis, una curva en forma de U invertida. Así, el nivel de precio de los productos, la orientación al ahorro de los consumidores, la sobrecarga de información, la tendencia a comparar y las preferencias iniciales cambian la forma de la relación.

Básicamente, la lección que podemos aprender de esta investigación es que no siempre el tamaño de los descuentos hace que se produzca una elección más cara. Así, si lo que pretende la marca o el distribuidor es incentivar la compra de los productos más caros (los que probablemente dejan mayor margen), en situaciones donde el consumidor tiene tendencia al ahorro, los precios son elevados, o el cliente es más maduro y medita diversas comparaciones, se puede producir una respuesta a la dosis en forma de U invertida.

¿Dónde está el punto que maximiza el efecto? Como siempre que hablamos de hormesis, es muy complejo saberlo. Los autores se refieren a ello como una de las limitaciones de su estudio.

LEE EL ARTÍCULO ORIGINAL AQUÍ:

Jia, H. M et al. (2018). Do Consumers Always Spend More When Coupon Face Value is Larger? The Inverted U-Shaped Effect of Coupon Face Value on Consumer Spending Level. Journal of Marketing, doi 10.1509/jm.14.0510

Indicadores de calidad de la revista*

Impact Factor (2017) Cuartil Categoría
Thomson-Reuters (JCR) 7.338 Q1 BUSINESS
Scimago (SJR) 8.616 Q1 MARKETING

* Es simplemente un indicador aproximado para valorar la calidad de la publicación

Todos los posts relacionados




(377). AVANZANDO EN LA COMPLEJIDAD; DOBLE DEPENDENCIA

[MONOTEMA] Una vez introducidos en los campos direccionales y en la filosofía de análisis de las ecuaciones diferenciales, vamos a complicar un poco más el caso a estudiar. Ahora, tenemos que el crecimiento de la población de individuos, es decir, la derivada de la población con respecto al tiempo no sólo depende de la propia población (y del tiempo, claro), sino que lo hace también de otra población que actúa como limitante.

Pensemos en una población de presas (conejos) y una de predadores (zorros). A la primera la llamaremos x, y a la segunda y. En este post vamos a mostrar qué sucede cuando existe esa dependencia, mostrando la propuesta de Lotka-Volterra.

Ecuaciones de Lotka-Volterra

Comenzamos especificando la ecuación:

Ahora tenemos dos parámetros más, alfa y beta. El primero representa la tasa de natalidad de las presas (conejos), y el segundo la de mortalidad cuando se encuentran con los predadores (zorros). Por tanto, el crecimiento dependerá de la diferencia entre los que nazcan y los que mueran.

Del mismo modo, podemos pensar que la población de zorros depende de la de conejos, si no hay otras presas que cazar. La ecuación de crecimiento para los zorros sería:

Aquí representa una especie de tasa de transferencia energética en cada interacción entre zorros y conejos, y la tasa de mortalidad de los zorros.

Se crea, por tanto, una interdependencia entre ambas poblaciones, siendo complejo a nivel intuitivo predecir cómo evolucionarán en el tiempo.

Solución de la ecuación

Podemos dividir ambas ecuaciones, y llegar a esta expresión:

De este modo, podemos intentar resolverla con Maxima. Para ello, hemos de especificar un valor para los 4 parámetros, y también unas condiciones iniciales para las dos variables: conejos y zorros.

En cuanto al valor de los parámetros, especificamos lo siguiente:

Y también exploramos dos condiciones iniciales. En la primera (y=6 , x=2), y en la segunda (y=2 , x=2).  Podemos ya, entonces, programar el siguiente código:

A:0.5;B:0.5;C:0.5;D:0.5;
ecuacionLotkaVolterra: ‘diff(y,x) -(((C*x*y)-(D*y))/((A*x)-(B*x*y)))=0;
solucionLotkaVolterra: ode2(ecuacionLotkaVolterra,y,x);
solucionfinal_1: ic1(solucionLotkaVolterra,y=6,x=2);
solucionfinal_2: ic1(solucionLotkaVolterra,y=2,x=2);
load(implicit_plot);
implicit_plot([solucionfinal_1,solucionfinal_2] ,[y,0,8],[x,0,8], [ylabel, “Zorros”], [xlabel, “Conejos”], [legend, false]);

Hay que observar que, aunque el tiempo se incremente indefinidamente la dinámica del sistema está regida por el conjunto límite que representa cada una de las curvas solución. Las poblaciones crecen y decrecen siempre de la misma forma.

Podemos representar los campos direccionales también:

load(drawdf)$
drawdf([x*(0.5-0.5*y),-y*(0.5-0.5*x)],[x,0,8], [y,0,8],[xlabel=”Conejos”], [ylabel=”Zorros”],[trajectory_at,2,6],[trajectory_at,2,2]);

b377_3

Conejos y zorros en el tiempo

Nos falta ahora conocer la evolución de los conejos y zorros en el tiempo. En realidad ya hemos visto que deber ser un recorrido periódico, aunque nos queda la duda de qué sucede realmente en el caso en que y=6 y x=2, porque parece que hay un momento en el que la trayectoria indica que la población de conejos y zorros es cero.

Para representar esa evolución temporal de las poblaciones, hemos de emplear el método de Runge-Kutta para obtener una solución numérica.  Vamos a ayudarnos del material creado por el profesor Renato Álvarez.

Para el caso de 2 conejos y 2 zorros (x=2, y=2), obtenemos lo siguiente:

kill (all);

sol:rk([x*(0.5-0.5*y),-y*(0.5-0.5*x)],[x,y],[2,2],[t,0,40,0.1]);

x:makelist([sol[k][1],sol[k][2]],k,1,length(sol))$
y:makelist([sol[k][1],sol[k][3]],k,1,length(sol))$
ciclo:makelist([sol[k][2],sol[k][3]],k,1,length(sol));

wxplot2d([[discrete,x],[discrete,y]],[legend, “Conejos”, “Zorros”],
[xlabel, “Tiempo”], [ylabel, “Población de conejos y zorros”] );

b377_5Pero para el caso de 2 conejos y 6 zorros (x=2, y=6), obtenemos lo siguiente:

kill (all);

sol:rk([x*(0.5-0.5*y),-y*(0.5-0.5*x)],[x,y],[2,6],[t,0,40,0.1]);

x:makelist([sol[k][1],sol[k][2]],k,1,length(sol))$
y:makelist([sol[k][1],sol[k][3]],k,1,length(sol))$
ciclo:makelist([sol[k][2],sol[k][3]],k,1,length(sol));

wxplot2d([[discrete,x],[discrete,y]],[legend, “Conejos”, “Zorros”],
[xlabel, “Tiempo”], [ylabel, “Población de conejos y zorros”] );

b377_4

Por tanto, los conejos se extinguen en esta segunda situación y, si no hay un repoblamiento o entran nuevos conejos en el hábitat, la población de zorros también se extinguirá.

Un recomendable texto para introducirse en el significado de la interdependencia puede encontrarse aquí.  Al mismo tiempo, los estudiantes pueden explorar de manera interactiva el comportamiento del sistema acoplado en las demostraciones de Wolfram.

Como puede verse en estos dos “pantallazos”, los resultados concuerdan con nuestra programación (los colores están cambiados)

b377_6b377_7

Conclusión

En esta simple introducción a la dinámica de la interacción para ilustrar la complejidad, los alumnos de marketing pueden comprobar que, cuando existe acoplamiento, el comportamiento de las poblaciones es extremadamente difícil de inferir sin ayuda de las matemáticas.

Los parámetros pueden variar (aquí los hemos mantenido todos igual a 0.5), y las condiciones iniciales también. Es más, el modelo se puede complicar con la introducción de competencia en el hábitat (los conejos compiten por recursos limitados de alimentación) y otras restricciones en las ecuaciones.

Todo esto nos sirve para entender que cuando se planea una intervención sobre un sistema acoplado que está en homeóstasis, las consecuencias sobre el sistema pueden ser contra intuitivas o indeseables  Y aquí el alumno puede imaginar sistemas biológicos, económicos, políticos, etc.

Todos los posts relacionados