(#454). NO LINEALIDAD Y MASCARILLAS

Uno de los objetivos de los primeros días de clase en mis asignaturas de toma de decisiones en marketing es que los alumnos entiendan la dificultad de los escenarios no lineales, tan característicos de la gran mayoría de fenómenos que se dan en la realidad. En general, los humanos tenemos muchos problemas para manejar nuestro pensamiento en este tipo de escenarios, por eso los estudiantes deben poner especial énfasis en comprenderlos.

En este artículo vamos a realizar varias simulaciones relacionadas con la desgracia que nos ha tocado vivir, el coronavirus que provoca la COVID-19, para mostrar con fines meramente didácticos la importancia de aplicar las herramientas del pensamiento y análisis no lineal. Es importante resaltar que el ejemplo que vamos a comentar es simplemente eso, un ejemplo, y que no tiene la entidad de un estudio científico, ya que es demasiado simple e incompleto. Sin embargo, ofrece unos resultados que, al menos, estimulan el pensamiento crítico.

Infectados en España

En la siguiente figura se muestran los infectados en España hasta el 26 de marzo, intervalo de tiempo donde el crecimiento era exponencial .

Como se puede ver, podemos ajustar la curva a través de un polinomio de orden 3 donde y representa la población de infectados y t el tiempo en días. Sin embargo, en aras de facilitar nuestro razonamiento y el análisis posterior, podemos realizar una aproximación más burda a la curva usando el siguiente razonamiento:

Lo que nos dice la expresión anterior es que el cambio temporal en el número infectados  es proporcional a la población que hay en el momento anterior al diferencial de tiempo, con un factor de 33.55%. Esto es algo que se aproxima bastantea la realidad, ya que en esos primeros 33 días el crecimiento del número de infectados estaba en torno al 42% en promedio.

No nos debemos preocupar demasiado por perder algo de exactitud porque los datos tampoco son perfectos, y hay retrasos en los reportes y otros factores que afectan. Por ejemplo, entre los días 12 y 13 el crecimiento fue solo del 6%, mientras que entre los días 14 y 15 fue del 91%.

Si discretizamos el tiempo, la expresión anterior es equivalente a esta (para cambios en t muy pequeños):

Es decir:

No obstante, seguiremos manteniendo el enfoque de continuidad para los análisis. De este modo, hay que resolverla ecuación diferencial para conocer la trayectoria del número de infectados. Resolviendo esa ecuación diferencial (ver este post como ayuda), obtenemos que:

El número de personas infectadas el día 33 era de 64059, mientras que la aproximación exponencial nos reporta una cifra de 64312.

Mascarillas

Durante muchas semanas, tanto la OMS como algunos responsables gubernamentales no consideraron pertinente recomendar el uso de mascarillas, desoyendo las indicaciones de científicos como, por ejemplo, Eric Feigl-Ding.

Las mascarillas no son perfectas, no proveen una protección total, pero incluso las más sencillas pueden tener porcentajes de efectividad en la contención. Si estipulamos un rango de efectividad entre el 0 ya el 100%, entonces si el 100% de la población no lleva habitualmente mascarilla (aquí se excluye obviamente el personal sanitario en el trato de enfermos), su protección es del 0%.

De este modo, podemos simular varios escenarios en los cuales se pasa de que nadie lleve mascarilla, a que diferentes porcentajes de la población las lleven (10%, 20% y 30%, respectivamente). Para aquellos que llevan mascarilla hemos estipulado dos niveles de protección (20% y 50%), en función del tipo de mascarilla y su comportamiento de riesgo. Los datos se muestran en la siguiente figura.

El siguiente paso ahora es modificar la ecuación diferencial para tener en cuenta una nueva variables que considere el uso de mascarillas.

donde es el valor promedio que describe el uso de mascarillas y el riesgo asociado a ellas.

donde es la frecuencia relativa de personas que llevan mascarilla normalizada en la escala [0,1], lo que quiere decir que si el 100% no lleva mascarilla le correspondería un 1. Por su parte, es el nivel de protección de cada mascarilla, también normalizado en una escala [0,1].

Hay que tener cuidado en este punto, porque ahora la ecuación de la trayectoria es una exponencial modificada con la nueva especificación:

De este modo si el 100% de personas no lleva mascarilla:

que es nuestra ecuación anterior de partida.

Y si el 100% de personas llevara mascarilla con un nivel de protección del 100% (algo imposible), entonces:

con lo que el número de infectados sería constante, es decir, no habría nuevos infectados.

Primer análisis de escenarios

Podemos comenzar un primer análisis con los 4 escenarios descritos en la figura que representa la distribución de las mascarillas. El primero de ellos, donde el 100% no las lleva (excepto el personal sanitario), ya lo tenemos. El objetivo es ahora estimar el número de infectados cuando una parte pequeña de la población lleva mascarilla. Los resultados se muestran a continuación. 

Como se puede apreciar, pequeños cambios en la distribución de uso de mascarillas producen enormes cambios en el número de infectados. Fijémonos en el escenario en que 90% no lleva mascarilla, un 5% lleva mascarillas que protegen un 20%, y un 5% lleva mascarillas que protegen un 50%. Tan sólo ese pequeño cambio produce un descenso del 37.8% en el número de infectados pasando a ser 40018. El descenso es del 75% en el último escenario, en el que el 70% de la gente sigue sin llevar mascarillas, pero el 15% lleva mascarillas de protección 20% y el 15% restante lleva mascarillas de protección 50%.

Riesgo simétrico

Hasta ahora hemos asumido que el riesgo de contagiarse o contagiar es constante para cada nivel de protección. Sin embargo, esta situación puede no ser correcta.

Imaginemos que el 100% de la población lleva mascarillas que protegen el 50%, entonces:

Pero ahora imaginemos que el 50% no lleva mascarilla y el otro 50% lleva mascarillas con la máxima protección (100%), entonces: 

El valor promedio es el mismo cuando la situación en cuanto a contagios podría ser muy diferente. Otras formas de riesgo simétrico podrían ser no lineales, por ejemplo, en modo de U-invertida. Así, esa función no lineal actuaría como una función de pesos para ponderar:

Como función de pesos podemos proponer la siguiente:

que da el máximo peso a los casos extremos, es decir, la máxima influencia sobre los contagios se produce cuando ninguno lleva mascarilla o cuando todos la llevan.

Dejamos como ejercicio para los estudiantes más curiosos el calcular ahora la trayectoria ante los 4 escenarios propuestos.

Riesgo asimétrico

El desafío ahora es concebir que el riesgo es asimétrico, una función convexa en forma de J que otorgue mayor peso a los casos en los que la protección es mayor a partir de cierto umbral. Podemos proponer la siguiente función:

Y ahora podemos a volver a realizar el análisis de escenarios, obteniendo que para el caso del escenario 2 (10% usan mascarilla), los infectados se habrían reducido un 23.4%; para el escenario 3 (20% usan mascarilla), los infectados se habrían reducido un 40.3%; y para el escenario 4 (30% usan mascarilla), lo infectados se habrían reducido un 53.5%.

Evidentemente, podemos simular un escenario que reflejara mucho más el compromiso de gobernantes y ciudadanos con el uso de mascarillas, haciendo que el 50% llevara mascarilla con una protección del 50%. Esto habría reducido los infectados un 83%.

Conclusión

Este artículo tiene que interpretarse como un mero ejercicio con fines de didácticos para estudiantes universitarios, y nunca como un estudio científico cuyas conclusiones sean robustas, ya que lo que se plantea es una gran simplificación de la situación, y se dejan muchas variables fuera.

La idea es hacer ver a los estudiantes que cuando se manejan ecuaciones no lineales nuestra mente tiene dificultades en inferir qué va a suceder, a no ser que nos ayudemos de herramientas matemáticas.

Lo que hemos visto es que incluso con pequeñas acciones como el incremento débil del uso de mascarillas justo cuando aparecieron los primeros casos, se podría haber reducido considerablemente el número de infectados, siempre dentro de las aproximaciones y asunciones que hemos realizado.

Sin embargo, insistimos en que lo importante de este artículo es motivar a los estudiantes a que profundicen en las dinámicas no lineales de los fenómenos que nos rodean, y cómo las asimetrías pueden afectar a los resultados.

Posts relacionados




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




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