Práctica de Montecarlo: La panadería de Pierre

por J.R

Desde hace tiempo buscaba un ejemplo práctico y ni demasiado sencillo o demasiado complicado del método de Montecarlo, y finalmente encontré un ejemplo financiero bastante singular en este PDF (todo el crédito de la elaboración o publicación del ejercicio es para su autor, recomiendo que lo descargues). Me pareció una buena oportunidad para escribir el primer post técnico de éste blog, así que he elaborado un pequeño tutorial del método basado en ese ejemplo y el lenguaje de programación Python. El método igual puede implementarse en cualquier otro lenguaje, no es mi intención que sea exclusivo para pythonistas. Para este fin explico primero el concepto de forma general, y luego un poco los atajos y peculiaridades que utilizo de python que lo diferencian de otros lenguajes. Espero resulte de provecho para el mayor número de personas.

El mejor resumen que se me ocurre del método de Montecarlo, como quisiera que lo explicaran en otro sitio, es que si conoces las probabilidades en porcentajes de que ocurra algo, entonces puedes obtener los números concretos si simulas que eso ocurre 10000 veces o más. Es como una historia artificial, en la cual simulas que la realidad ocurre de forma muy acelerada, de forma que días o años pueden pasar en milisegundos. Así puedes sacar conclusiones sin esperar tanto y sin necesariamente conocer a profundidad los detalles involucrados. Y de hecho, el método es precisamente utilizado en los casos en que esos detalles sobrepasan las técnicas analíticas conocidas, en campos tan diversos como las finanzas, la guerra y la ciencia. Es considerado, además, uno de los algoritmos más importantes del siglo XX.

La raíz del truco se encuentra en la generación de números aleatorios, y afortunadamente casi todos los  lenguajes de programación poseen uno. En Python, por ejemplo, para generar un número aleatorio basta con importar el módulo random, y hacer una llamada a random.random(). Esto genera de forma aleatoria un decimal entre 0 y 1. Es conveniente mencionar que los números generados de esa manera realmente son pseudoaleatorios y la distribución estadística es uniforme (todos los casos son igual de probables), pero lo que particularmente nos interesa en el tutorial es que ese intervalo entre 0 y 1 convenientemente también es la forma en que se expresan en formal decimal los porcentajes. Supongamos que en un salón de clases el 60% de los alumnos es de estatura media, un 20% son altos y otro 20% son bajos. Podríamos pensar entonces que los números son 0.2 (altura alta), 0.6 (altura media) y 0.2 (altura baja). El truco para cuadrar los números aleatorios dentro de esos decimales está en calcular las probabilidades acumuladas.  La altura media es la alta más ella misma, y la altura alta es es la suma de todas (podría hacerse al revés sin problemas, pero sigo el estilo del ejercicio de la panadería):

        Altura               Alta   Media  Baja 
Probabilidad porcentual:     0.2    0.6    0.2
Probabilidad acumulada :     0.2    0.8    1.0

Dado que la suma total siempre equivale a 1, al calcular un número aleatorio decimal entre 0 y 1 podemos encasillarlo dentro de las probabilidades comparando si es menor que alguna de las probabilidades acumuladas. Es allí, en la idea de utilizar los números aleatorios para eso y hacer el experimento artificial miles de veces en una PC, donde está el corazón del asunto. El resto de este tutorial es simplemente una implementación en Python aplicada a un par de casos particulares. En el primero, antes de ir a ayudar a Pierre con su panadería vamos a explorar un poco las alturas comentadas antes y crear en el proceso las mismas funciones que utilizaremos después.

Tenemos tres necesidades generales y, por lo tanto, elaboraremos tres funciones para :

  • Convertir una lista de porcentajes decimales (las alturas) en un intervalo de probabilidades acumuladas.
  • Encontrar la posición que ocupa un número aleatorio en dicho intervalo (nosotros recordaremos el significado de esa posición gracias a otra lista)
  • Generar de la manera más versátil posible posiciones aleatorias en base a una distribución de porcentajes (una combinación de ambas funciones anteriores)

Antes de seguir quisiera confesar que encuentro más cómodo programar usando nombres de variables en inglés (a pesar que el PDF anterior esté en español). Me parece que si ya las palabras claves del lenguaje están en inglés, escribir el resto en el mismo idioma hace que sea más fácil de leer. En todo caso sólo es algo personal,  éstas son las funciones, a continuación comento sus detalles.

distribution_interval, aunque breve, utiliza algunas características únicas de Python. enumerate es la forma peculiar que tiene Python para hacer un ciclo for tradicional recordando posiciones (en este caso es “n”, “d” es cada elemento del intervalo), y dists[0:n+1] es el list slicing que obtiene el elemento n y todos los anteriores de la lista.

interval_position simplemente recorre el intervalo dado (generado por distribution_interval) y revisa cual elemento es menor que el número aleatorio rn. Si retorna 0 entonces el número es menor que el primer elemento, y así sucesivamente. Como comenté antes, seremos nosotros los que recordaremos el significado de cada posición.

position_generator utiliza las funciones anteriores, creando el intervalo a partir de la distribución porcentual y enviándole un número pseudoaleatorio a interval_position, encapsulando todo dentro de un generador. Éste objeto característico de Python funciona como un generador (… perdón) de valores, y tiene dos ventajas: No genera una lista inmensa de items, sino que los emite uno a uno, lo cual resulta bastante útil si el generador es infinito, ya que la memoria de la PC no lo es. Realmente es un artificio bastante útil de los desarrolladores del lenguaje, aunque confieso que realmente lo uso aquí porque me permitirá ser vanidoso y hacer la simulación en una sola línea más adelante.

El código de la simulación sería el siguiente:

La simulación ocurre principalmente en la línea simulation = heights[generator.next()]. Allí el generador emite un valor posicional (0, 1 o 2) dentro del intervalo de distribuciones acumuladas, y en el mismo paso lo usamos para obtener la cadena con esa misma posición en la lista de heights. Ya sabemos entonces si en ese experimento resultó una altura alta, media o baja, por lo que queda es anotar que caso fué, y al final dividir los tres resultados por el número de veces que realizamos el método. La salida del programa debería ser más o menos similar a esto:

Alta    0.1986
Media   0.5967
Baja    0.2047

Los resultados no tienen por qué ser iguales, ya que estamos trabajando con números aleatorios. Si bien parece que no logramos realmente nada y obtuvimos lo que ya sabíamos, ésto resulta mucho más útil cuando tratamos de predecir el resultado de varias probabilidades que dependen unas de otras. Y ahora si, vamos con el ejemplo. El PDF que vinculé al inicio del artículo nos presenta el siguiente enunciado:

La panadería de Pierre hace y vende pan francés. Cada mañana, la panadería satisface la demanda del día con pan recién horneado. Pierre puede hacer el pan únicamente en lotes de una docena de panes. Cada pan tiene un costo de fabricación de 25c de peso. Supondremos, por simplicidad, que la demanda diaria total de pan también se presenta en múltiplos de 12. Los datos demuestran que esta demanda varía de 36 a 96 panes diarios. Un pan se vende a 40c de peso y si sobra al final del día se vende a una cocina de beneficencia a un precio de recuperación de 10c de peso por pan. Si la demanda es mayor que la oferta, suponemos que hay un costo de ganancia perdida de 15c de peso/pan, debido a la pérdida de clientes que van con los competidores, etc. Los registros de la panadería muestran que la demanda diaria se puede clasificar en tres tipos: alta, media y baja. Estas demandas se presentan con probabilidades de .30, .45 y .25, respectivamente. La distribución de la demanda por categorías aparece en la Tabla. Pierre quisiera determinar el número óptimo de panes que debe hacer cada día para maximizar la ganancia (ingresos + ingresos de recuperación – costo de fabricación – costo de ingresos perdidos).

Y la tabla de cada tipo de demanda, algo preformateada para que sea fácil usarla en Python:

Demanda: [  36,   48,   60,   72,   84,   96]
Alta     (0.05, 0.10, 0.25, 0.30, 0.20, 0.10)
Media    (0.10, 0.20, 0.30, 0.25, 0.10, 0.05)
Baja     (0.15, 0.25, 0.35, 0.15, 0.05, 0.05)

Una de las cosas bellas de los números es que simplifican mucho las cosas. Si eres como yo, quizás al principio luego de leer el enunciado te cueste ver cómo lo que hicimos antes puede aplicarse a eso. Realmente este caso es idéntico al anterior, sólo que con un paso adicional: Primero determinamos con los números .30, .45 y .25 si la demanda es alta, media o baja, y después determinamos cuántos panes fueron vendidos ese día de acuerdo con la distribución de cada tipo de demanda. Esto cambia las posibilidades, por ejemplo, es más posible que en un día de demanda baja se compren 36 panes (0.15) que en un día de demanda alta (0.05).

Aquello corresponde a la realidad externa de la panadería, la demanda de los consumidores. La simulación principalmente se encarga de eso, pero también tenemos que recordar la otra realidad, la interna de la panadería de Pierre. En primer lugar, que tendremos que repetir la simulación con cada tipo de demanda por docena, como si Pierre decidiera arbitrariamente un día producir 36, 48 o 96 panes, y que los panes son producidos en docenas por la forma y el tamaño de su horno (mejor pensarlo así que suponer que es para simplificar el ejercicio). Y, por último, que como toda pequeña empresa su panadería tiene una forma específica de calcular sus ingresos, por lo que tenemos que hacer una función de finanzas para saber como le fue en un determinado día. Esta la podemos hacer en base a la oferta y la demanda del siguiente modo:

Realmente las finanzas no son mi fuerte, la única razón por la que sé que eso es correcto es porque al final mis resultados coinciden con los del PDF. En un caso real, por supuesto, habría que ser mucho más cuidadoso, y de ser posible siempre probar nuestras simulaciones con datos preexistentes para estar seguros de lo que estamos haciendo. Con esa función y las anteriores generales, sin embargo, ya estamos listos para realizar la simulación definitiva.

Aquí creamos 4 generadores, uno para saber el tipo de demanda, y otros tres para la distribución de cada uno. Los ubicamos en una lista con la misma posición en que ordenamos los tipos de demanda, y luego hacemos algo muy similar al caso de las alturas. En este caso, toda la magia es realizada en la línea

demand = positions[demand_amount[demand_type.next()].next()]

Es exactamente lo mismo que realizamos con las alturas, sólo que utilizando dos generadores y dos listas en vez de una. Para simulaciones más complejas sería conveniente crear un objeto que las ordenara automáticamente, pero aún con dos resulta manejable en una sola línea.

Posteriormente, con total += profit(supply, demand) vamos acumulando las ganancias, y repetimos la simulación 10000 veces, como si pasaran todos esos días mientras Pierre produce una cantidad de panes a capricho sin pensar en si perderá o ganará dinero ese día. No tiene porqué hacerlo e ir a la quiebra, claro, porque todo esto ocurre en una simulación en nuestro CPU. El resultado de correr el código anterior debería ser similar (no igual) a esta salida:

36  1.33704
48  4.35042
60  6.45066
72  6.831
84  6.0399
96  4.60692

La producción de panes diaria que le proporciona mayores beneficios es 72, y resulta que es igual que cuando resuelves el problema de forma analítica, según el PDF:

ingresos pan

Ya sea por curiosidad o porque nos vayamos a dedicar a eso, sería conveniente también saber como obtener estos valores exactos según cálculos de probabilidad. Pero eso ya sería un artículo distinto, aquí termina este tutorial, espero que haya sido entretenido. Para mayores detalles sobre el método de Montecarlo el PDF de donde obtuve el ejercicio es bastante extenso, y también puedes leer calculando cosas mientras jugamos a la ruleta de Pybonacci, donde el mismo método es utilizado para calcular areas de polígonos, prueba de gran versatilidad que tiene. ¡Cualquier detalle o sugerencia por favor comentar!

Anuncios