¿Por qué Monte Carlo? Te estarás preguntando qué tiene que ver un famoso casino de Mónaco con las matemáticas. La respuesta está en el azar. Durante la década de 1940, científicos como Stanislaw Ulam comenzaron a usar simulación aleatoria para resolver problemas difíciles, y bautizaron el método en honor al casino de Monte Carlo, donde la suerte y la probabilidad son reyes. De hecho, se cuenta que Ulam tenía un tío tan aficionado a las apuestas que siempre hablaba de ir a Monte Carlo – ¡y así quedó el nombre pegadizo para este método basado en juegos de azar!

Pero no necesitas viajar a un casino para entender el método de Monte Carlo. Podemos explicarlo con algo mucho más cotidiano. Imagina que tienes una bolsa llena de bolitas de colores (o caramelos). Algunas son rojas y otras de otros colores. Si sacas una bolita al azar, ¿cómo podrías averiguar la proporción de bolitas rojas en la bolsa sin mirar dentro? Una forma es repetir el experimento muchas veces: sacas una bolita, anotas si es roja o no, la devuelves y vuelves a mezclar, y así sucesivamente. A medida que haces más y más extracciones, la fracción de veces que obtienes una roja debería acercarse a la proporción real de rojas en la bolsa. ¡Eso, en esencia, es el método de Monte Carlo! Usas la simulación repetitiva para estimar un resultado que sería difícil de calcular directamente.

En resumen, el método de Monte Carlo se llama así por la conexión con el azar (como en un casino), y su idea fundamental se puede entender como hacer muchos experimentos aleatorios para aproximar respuestas. Es un enfoque poderoso, utilizado hoy en día en campos tan diversos como la física, la economía o la informática, cada vez que un problema es demasiado complejo para resolverlo con fórmulas exactas. Ahora sí, comencemos nuestro viaje didáctico por este método de forma cercana y motivadora, ¡verás que tú también puedes dominarlo paso a paso!

Simular, contar, promediar: la idea básica de Monte Carlo

El corazón del método de Monte Carlo se resume en tres acciones: simular experimentos aleatorios, contar resultados y promediar para obtener una estimación. Esto suena simple, ¡y lo es! Veamos qué significa cada paso:

• Simular: significa generar datos aleatorios que imitan el fenómeno que queremos estudiar. Por ejemplo, si quisiéramos simular sacar bolitas de una bolsa, podemos usar un generador de números aleatorios para decidir el color de cada bolita simulada.
• Contar: después de cada simulación, registramos lo que pasó. ¿Salió bolita roja o no? ¿El “experimento” fue éxito o fracaso? ¿El valor simulado cumple cierta condición? Vamos contando estos resultados.
• Promediar: finalmente, calculamos la proporción o promedio de interés a partir de nuestros datos simulados. Ese promedio es nuestra estimación Monte Carlo para el resultado deseado.

¿Por qué funciona esto? Gracias a un principio estadístico fundamental llamado ley de los grandes números. En términos sencillos: cuantos más experimentos aleatorios hagas, más se acercará el promedio de tus resultados simulados al valor verdadero (teórico) que buscas. Al principio, con pocos intentos, tu estimación puede estar algo desviada, pero si repites la simulación cientos o miles de veces, el promedio “se estabiliza” alrededor del valor correcto.

En la práctica, para simular utilizamos típicamente números aleatorios uniformes entre 0 y 1. “Uniforme” significa que cualquier número entre 0 y 1 es igualmente probable. Puedes pensar en ello como un generador de "decimales sorpresa". A partir de esos números uniformes se puede generar cualquier otro tipo de comportamiento aleatorio (veremos cómo más adelante). Lo genial es que lenguajes como R nos dan funciones listas para obtener estos aleatorios uniformes fácilmente.

Entonces, resumido en una línea: el método Monte Carlo consiste en usar la computadora (o cualquier fuente de aleatoriedad) para "lanzar dados virtuales" muchas veces, registrar qué ocurre en cada lanzamiento, y luego sacar el promedio de esos resultados para estimar la respuesta a nuestro problema. ¡Así de simple es la idea central!

Ejemplo: Caramelos de colores (probabilidad tipo Bernoulli)

Pasemos a un ejemplo concreto y dulce. Imagina que tu bolsa tiene caramelos y un 25% de ellos son de color rojo (es decir, 1 de cada 4 en promedio debería ser rojo). Si metes la mano y sacas un caramelo al azar, la probabilidad de que sea rojo es 0.25 (25%). Ahora bien, ¿cómo simularías ese proceso con la computadora usando R? Vamos a hacerlo paso a paso.

En R existe una función para generar números uniformes entre 0 y 1 llamada runif(). Podemos usar esos números para decidir el color de un caramelo simulado. La idea es: si el número aleatorio generado es menor que 0.25, consideramos que "sacamos" un caramelo rojo (éxito); si es mayor o igual a 0.25, entonces no es rojo. En términos de variable aleatoria, estamos definiendo una variable de tipo Bernoulli (que vale 1 con probabilidad 0.25 y 0 con probabilidad 0.75).

La sintaxis de R nos permite hacer esto de forma concisa. Generemos, por ejemplo, 10 extracciones simuladas de caramelos y veamos cuántos resultan rojos:

set.seed(123)
runif(5)
runif(4, min = -1, max = 1)

set.seed(1)
as.integer(runif(10) < 0.25)

set.seed(1)
rbinom(10, size = 1, prob = 0.25)

# Simular 10 extracciones de caramelos (1 = rojo, 0 = no rojo)
as.integer(runif(10) < 0.25)
## [1] 0 0 1 0 0 1 0 0 0 1

En el resultado anterior, 1 indica que ese intento obtuvo un caramelo rojo, y 0 que no. Podemos ver, por ejemplo, que en esas 10 simulaciones salieron 3 caramelos rojos (los valores "1" en la salida). ¿Es eso mucho o poco? Bueno, con tan solo 10 intentos, tuvimos 3 rojos, lo cual es un 30%. El verdadero porcentaje era 25%, pero con tan pocos intentos es normal que la estimación (30% en este caso) fluctúe un poco.

Repitamos la simulación con un número más grande de extracciones, digamos 1000, y calculemos la fracción de caramelos rojos simulados (que debería acercarse a 0.25):

# Repetir el experimento con 1000 extracciones
mean(as.integer(runif(1000) < 0.25))
## [1] 0.247

¡Ahora sí! Al simular 1000 extracciones, el promedio (proporción de veces que obtuvimos rojo) fue aproximadamente 0.247, muy cercano al 0.25 esperado. Si aumentamos a 10,000 o 100,000 simulaciones, esta media estaría aún más cerca de 0.25 (probablemente algo como 0.2503, 0.2498, etc.). Aquí estamos viendo en acción la idea de "simular, contar, promediar": generamos datos aleatorios según la probabilidad dada, contamos cuántos éxitos obtenemos, y el promedio (número de éxitos dividido entre el total de intentos) nos da una estimación de la probabilidad.

Este sencillo ejemplo de los caramelos demuestra cómo Monte Carlo nos permite aproximar probabilidades. No necesitábamos matemáticas avanzadas para llegar al 25% – lo obtuvimos simulando. Por supuesto, en este caso conocíamos de antemano la respuesta (0.25), pero el método se vuelve realmente útil cuando no sabemos la respuesta teórica o es muy difícil de calcular. Por ejemplo, ¿y si la probabilidad no fuera tan evidente? ¿O si el proceso físico fuera complicado? Monte Carlo nos ofrece una herramienta para abordarlo: basta con que sepamos cómo simular el proceso en pequeño, y luego dejar que la ley de los grandes números haga su magia mediante muchos intentos.

Ejemplo: Lanzando dardos a un tablero (estimación de áreas por Monte Carlo)

Ahora veremos un ejemplo visual y divertido: imaginar lanzar dardos a una diana pintada en un tablero cuadrado. Supongamos que tenemos un tablero cuadrado y dibujamos en su interior un círculo más pequeño. Queremos estimar el área (superficie) del círculo sin ponernos a medir con fórmulas geométricas, solo mediante lanzamientos aleatorios de dardos.

La idea es similar a la de los caramelos, pero en vez de contar "rojo o no rojo", contaremos "el dardo cayó dentro del círculo o fuera de él (pero dentro del cuadrado)". Cada lanzamiento de dardo al azar al tablero es un experimento. Si un dardo cae dentro del círculo, lo consideramos un "éxito". Si cae en la parte del cuadrado que está fuera del círculo, es un "fracaso" para nuestro conteo. ¿Qué probabilidad tiene un dardo de caer dentro del círculo? Eso será precisamente la proporción área_círculo/area_cuadrado, porque asumimos que cada punto del tablero tiene la misma probabilidad de ser alcanzado.

En otras palabras, fracción de dardos dentro ≈ (área del círculo) / (área del cuadrado). Despejando, podemos estimar el área del círculo como área_cuadrado * (número de dardos dentro / número total de dardos). Si conocemos el área del cuadrado (por ejemplo, 1 unidad cuadrada si tomamos el lado = 1), entonces básicamente la fracción de aciertos dentro es ya una estimación del área del círculo en esas mismas unidades.

Hagámoslo animado para que se entienda mejor. En la siguiente figura simularemos lanzar 1000 dardos de manera aleatoria sobre un cuadrado de lado 1 (área = 1). Dentro de ese cuadrado hemos dibujado un círculo de radio 0.4 (es decir, bastante más pequeño que el cuadrado). El área real de ese círculo, que podríamos calcular con la fórmula geométrica, sería π*(0.4)^2 ≈ 0.503 unidades cuadradas (aproximadamente el 50.3% del área del cuadrado). Veamos cómo Monte Carlo puede acercarse a ese valor sin usar la fórmula, solo con puntos al azar:

Observa en la animación cómo los dardos (puntos) van cubriendo el tablero de forma aleatoria. Algunos caen dentro del círculo negro y se colorean de verde, mientras otros quedan fuera y se marcan en rojo. Al final del experimento, podemos calcular la proporción de dardos verdes (dentro del círculo). Si todo va bien, esa proporción debería acercarse bastante a ~0.503, que es el área real del círculo respecto al cuadrado.

De hecho, mientras más dardos lancemos, más precisa será la estimación. Si solo lanzáramos 10 dardos, la proporción podría variar mucho (tal vez 0.4, tal vez 0.6, dependiendo del azar). Con 1000 dardos, esperamos algo mucho más cercano a ~0.5. ¡Con 100,000 dardos virtuales podríamos obtener quizá 0.5028 o 0.5041, súper cerca! Esa es la fortaleza de Monte Carlo: repetición masiva para afinar resultados.

Podemos replicar este experimento con código R también. Vamos a simular 1000 dardos numéricamente y ver qué obtenemos. En este caso, generaremos pares de coordenadas (x, y) al azar dentro del cuadrado [0,1]×[0,1]. Luego comprobaremos qué proporción de esos puntos cae dentro del círculo de radio 0.4 centrado en el centro del cuadrado (coordenada del centro (0.5, 0.5)). Un punto está dentro del círculo si su distancia al centro <= 0.4, es decir, si (x-0.5)^2 + (y-0.5)^2 <= 0.4^2.

# Simular 1000 lanzamientos de dardo en el cuadrado unitario
n <- 1000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)
dentro <- as.integer((x - 0.5)^2 + (y - 0.5)^2 <= 0.4^2)
mean(dentro)          # fracción de dardos dentro del círculo
## [1] 0.508
mean(dentro) * 1      # área estimada del círculo (lado del cuadrado = 1, área=1)
## [1] 0.508

dentro <- as.integer((x - 0.5)^2 + (y - 0.5)^2 <= 0.4^2)

En una ejecución típica, obtuvimos una fracción de 0.508 (o 50.8% de los dardos dentro). Esto equivale a un área estimada de 0.508 para el círculo. El valor real era ~0.503, así que nuestro error fue de tan solo unas milésimas. Nada mal para una aproximación basada en tiros al azar, ¿no? Si repites el código, obtendrás valores cercanos que fluctúan ligeramente. Esa variabilidad es normal en Monte Carlo, pero podemos hacerla tan pequeña como queramos aumentando el número de simulaciones (a costa de más tiempo de cómputo, claro).

Este experimento del tablero y el círculo es un ejemplo clásico del método de Monte Carlo aplicado a estimación de áreas y de π. De hecho, quizás notes que si el círculo hubiera sido de diámetro igual al lado del cuadrado (o sea, un círculo que toca los bordes, inscripto en el cuadrado), la fracción dentro sería aproximadamente π/4 (~0.785) y multiplicando por 4 obtendríamos una estimación del famoso número π. ¡Así es, se puede estimar π lanzando dardos! Aquí usamos un círculo más pequeño deliberadamente para variar el ejemplo, pero el principio es el mismo.

Hemos aprendido que Monte Carlo puede estimar áreas o probabilidades geométricas muy eficazmente. Es como si en lugar de calcular integrales o usar fórmulas, contáramos puntos. Y a veces, contar puntos es más fácil que resolver ecuaciones complicadas. A continuación, profundizaremos un poco más presentando otra técnica interesante dentro de Monte Carlo: la transformación inversa, que sirve para simular tiempos de espera y otras distribuciones no uniformes.

Transformada inversa: simulando tiempos de espera (el ejemplo del autobús)

Hasta ahora hemos simulado variables bastante simples: básicamente casos de "éxito o fracaso" (rojo/no rojo, dentro/fuera) o coordenadas uniformes en un área. ¿Qué pasa cuando queremos simular una variable aleatoria que siga alguna distribución específica diferente de la uniforme? Por ejemplo, los tiempos de espera de un autobús: es conocido que si los autobuses pasan de manera "completamente aleatoria" por una parada, los tiempos de espera siguen aproximadamente una distribución exponencial. Esto significa que hay muchos tiempos de espera cortos, y pocos muy largos, decayendo exponencialmente la probabilidad a medida que el tiempo crece.

La distribución exponencial tiene una propiedad llamada "sin memoria": básicamente siempre "se olvida" de cuánto llevas esperando, la probabilidad de que falte cierto tiempo adicional para que llegue sigue siendo la misma. Sin entrar en detalles matemáticos profundos, la función de distribución acumulada (CDF) de un tiempo de espera exponencial con parámetro λ (lambda) es:

\(F(t)=1-e^{-\lambda t},\; t\ge 0\)

Por ejemplo, si λ = 0.2 (en unidades de "por minuto"), la media de la distribución sería 5 minutos (ya que la media de un exponencial es \(1/\lambda\)). Con esa λ, tendríamos F(5) = 1 - e^-0.2*5 ≈ 0.632, es decir, un 63.2% de probabilidad de que el bus haya llegado en 5 minutos o menos.

Bueno, ¿y cómo simulamos un tiempo aleatorio que siga esta distribución exponencial? Aquí es donde entra el poderoso método de la transformada inversa. Según este método, si U es un número uniforme ~ (0,1), entonces podemos generar una variable aleatoria con cualquier distribución dada aplicando la función inversa de la CDF a ese U. En otras palabras, para obtener un valor X que siga la distribución con CDF F, tomamos X = F-1(U).

En el caso de la distribución exponencial, tenemos F(t) = 1 - e^-λt. Si ponemos u = F(t) y despejamos t en términos de u, obtenemos la inversa:

\(t = F^{-1}(u) = -\tfrac{1}{\lambda}\,\ln(1-u)\)

Más comúnmente se escribe \(X = -\tfrac{1}{\lambda}\,\ln(1-U)\). Y como 1-U es también uniforme (restarle a 1 no cambia la distribución), a veces lo verás simplemente como \(-\tfrac{1}{\lambda}\,\ln U\).

¿Qué significa esto en la práctica? Que podemos tomar un número uniforme aleatorio U y "transformarlo" con esa fórmula para obtener un X que se distribuye exponencialmente. ¡Voilà! La magia de Monte Carlo nos permite simular tiempos de espera (y muchas otras cosas) a partir de uniformes.

Hagamos un ejemplo concreto: supongamos que el tiempo de espera de nuestro autobús sigue Exp(λ = 0.2 1/minuto), es decir, en promedio 5 minutos de espera. Generemos algunos tiempos de espera simulados usando la fórmula anterior. También usaremos la función incorporada rexp() de R que genera exponenciales, para comparar:

# Simular 5 tiempos de espera de autobús con λ = 0.2 (media = 5 min)
set.seed(123)  # fijamos una semilla para reproducibilidad
U <- runif(5)
tiempos <- -log(1 - U) / 0.2
tiempos
## [1]  2.251  2.948 13.168  6.041  1.634

# Verificar con la función rexp de R (debería dar resultados de similar distribución)
rexp(5, rate = 0.2)
## [1]  2.572  5.899 13.494  2.233  2.481

En la salida vemos dos conjuntos de 5 números. Los primeros (tiempos) los obtuvimos transformando uniformes con la fórmula -log(1-U)/0.2. Los segundos los generó directamente rexp. Ambos "lucen" como posibles tiempos en minutos: la mayoría son algunos minutos (2.2, 2.9, 6.0, etc.) pero de vez en cuando aparece uno más grande (13 minutos en el ejemplo). No son iguales en ambas listas (cada simulación aleatoria produce valores diferentes), pero estadísticamente provienen de la misma distribución. Lo importante es que hemos aprendido a simular una variable continua no trivial (tiempo de espera exponencial) usando Monte Carlo.

La técnica de la transformada inversa se puede aplicar siempre que conozcamos la CDF de la distribución y que podamos invertirla matemáticamente. Para algunas distribuciones hay que usar métodos más avanzados cuando no existe una inversa cerrada, pero en muchísimos casos (como exponencial, uniforme, algunas gammas, etc.) sí podemos obtenerla y es una herramienta excelente para generar valores aleatorios con la distribución deseada.

¿Por qué nos detuvimos a ver esto? Porque muchas aplicaciones de Monte Carlo requieren generar datos aleatorios con ciertas distribuciones (no todo es 0/1 o uniformes). Por ejemplo, simulación de precios en economía (podrían seguir distribuciones normales), simulación de fallos en máquinas (podrían seguir distribuciones Weibull), etc. Con el método de Monte Carlo y la transformada inversa en tu arsenal, ¡puedes crear prácticamente cualquier tipo de variable aleatoria que necesites para un experimento simulado!

Integración Monte Carlo: calculando el área bajo una curva

Hemos visto Monte Carlo aplicado a probabilidades y a áreas geométricas. Vamos ahora a atacar un problema clásico de integración usando Monte Carlo. La integración, desde el punto de vista geométrico, no es más que calcular áreas bajo curvas. ¿Y qué acabamos de hacer con los dardos? ¡Calcular áreas! Solo que en lugar de una curva teníamos un círculo. Pero podemos usar la misma idea para aproximar integrales difíciles.

Supongamos que queremos calcular la siguiente integral definida:

I = ∫₀^π/2 cos(x) dx.

Si recuerdas un poco de cálculo, puede que identifiques que la integral de cos(x) es sin(x), y evaluar sin(x) de 0 a π/2 nos da sin(π/2) - sin(0) = 1 - 0 = 1. Así que en este caso sabemos que el área bajo la curva g(x) = cos(x) desde 0 hasta π/2 es exactamente 1 (en unidades cuadradas). Pero pretendamos por un momento que no supiéramos resolver esa integral analíticamente, o que g(x) fuera una función más complicada cuyo área no se puede obtener con una fórmula simple. Monte Carlo viene al rescate.

La técnica es muy similar a la de los dardos. Imagina el área bajo cos(x) entre 0 y π/2 como la silueta de una colina o montaña. Está encima del eje x (entre y=0 y y=1, ya que cos(x) va de 1 a 0 en ese intervalo) y debajo de la curva cos(x). Dibujemos ese cuadro: un rectángulo que enmarque la región (de 0 a π/2 en el eje x, y de 0 a 1 en el eje y), y la curva cos(x) trazando un arco descendente. Monte Carlo integrará "tirando puntos" dentro de ese rectángulo y contando cuántos quedan debajo de la curva.

En la siguiente animación, cada punto verde representa una muestra aleatoria que cayó por debajo de la curva cos(x) (es decir, y <= cos(x) en esa coordenada), mientras que los puntos rojos quedan por encima de la curva (y > cos(x), o sea fuera del área bajo la función). La fracción de puntos que quedan bajo la curva, multiplicada por el área total del rectángulo, nos dará la estimación del área buscada. El rectángulo aquí tiene base π/2 y altura 1, así que su área es π/2 ≈ 1.5708. Esperamos que ~63.66% de los puntos caigan bajo la curva (porque 1 es ~63.66% de 1.5708). Veámoslo en acción:

Tal como se esperaba, a medida que se lanzan muchos puntos, la nube verde delimita la forma de la curva cos(x). La estimación de área se va acercando a 1.0 (que es el valor real). En Monte Carlo, no es raro tener un pequeño error del orden de centésimas con unos pocos miles de puntos, pero ese error se reduce si seguimos aumentando simulaciones.

Podemos hacer el cálculo Monte Carlo también por código, para verificar numéricamente el resultado. Utilizaremos el enfoque formal de integración Monte Carlo: generar n puntos x uniformes en [0, π/2], evaluar g(x) = cos(x) en ellos, y promediar esos valores multiplicando por el largo del intervalo (π/2). En fórmula: \( \hat{I}_n = \frac{\pi/2 - 0}{n} \sum_{i=1}^n cos(x_i) \). El término (π/2 - 0) es la longitud del intervalo, aquí π/2.

# Estimar I = ∫[0, π/2] cos(x) dx mediante Monte Carlo
n <- 10000
x <- runif(n, min = 0, max = pi/2)
y <- cos(x)
est <- (pi/2) * mean(y)
est
## [1] 0.9987

Con 10,000 muestras uniformes, la estimación resultó ser ~0.9987, lo cual está a apenas ~0.13% de error respecto al valor verdadero 1. Si aumentáramos n, podríamos afinar aún más. Notemos cómo este método no calculó ninguna integral simbólica ni hizo divisiones en rectángulos al estilo tradicional; simplemente lanzó "muestreos aleatorios" bajo la curva.

Una forma intuitiva de entender la integración Monte Carlo es pensar que estamos contando puntos debajo de la montaña. Imagina que llovieran gotas uniformemente sobre todo el rectángulo que enmarca la región bajo la curva. La proporción de gotas que llegan al suelo (que tocan la zona bajo la curva y no se "quedan en el aire" por encima de la curva) es exactamente el área de la región bajo la curva dividida por el área total del rectángulo. Esto es análogo a lo que hicimos con los dardos y con los caramelos. Monte Carlo aprovecha esa idea para resolver integrales: "muestreamos" puntos y vemos qué fracción cae en la región de interés.

Aunque aquí integramos cos(x) que era fácil, en problemas reales a veces g(x) puede ser una función complicada o alta dimensional donde los métodos tradicionales fallan. Monte Carlo no se inmuta: mientras podamos generar puntos aleatorios y comprobar si caen debajo de la superficie (en un sentido general), podemos estimar el valor de integrales, volúmenes, expectativas, etc. En la sección final, te propondré algunos ejercicios para que practiques estas ideas con otros ejemplos interesantes.