{ "cells": [ { "cell_type": "markdown", "id": "b1f26c1e", "metadata": {}, "source": [ "# Simulación de Monte Carlo paso a paso: estimar π lanzando dardos\n", "\n", "**Objetivo.** Explicar, con el ejemplo más simple posible, las **6 etapas generales** de una simulación de Monte Carlo y sus **4 componentes fundamentales**, calcular **varias métricas de análisis de resultados** (no solo la media) y entender, con una demostración empírica, **por qué el Teorema del Límite Central (TLC)** es lo que hace válidos los intervalos de confianza que reportamos.\n", "\n", "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3b677926", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib.patches import Circle, Rectangle, FancyBboxPatch, FancyArrowPatch\n", "from matplotlib import animation\n", "from IPython.display import HTML\n", "\n", "np.random.seed(2026)\n" ] }, { "cell_type": "markdown", "id": "f272eef9", "metadata": {}, "source": [ "## 1. El mapa general: las 6 etapas de una simulación de Monte Carlo\n", "\n", "**Objetivo general (el mismo de la infografía de clase):** transformar el problema de interés en la estimación de un valor esperado de una variable de interés.\n", "\n", "El siguiente diagrama resume el flujo que vamos a seguir. Cada \"Paso\" de este notebook corresponde a una de estas 6 etapas — y al final vamos a conectar todo con los **4 componentes** fundamentales de cualquier simulación (generador de números pseudoaleatorios, distribuciones de entrada, generadores de variables aleatorias, y estimador de respuesta)." ] }, { "cell_type": "code", "execution_count": null, "id": "71fab7e5", "metadata": {}, "outputs": [], "source": [ "def dibujar_paso(ax, x, y, w, h, numero, texto, color=\"#F2994A\"):\n", " box = FancyBboxPatch((x, y), w, h, boxstyle=\"round,pad=0.02,rounding_size=0.08\",\n", " facecolor=\"white\", edgecolor=color, lw=2, zorder=2)\n", " ax.add_patch(box)\n", " ax.add_patch(Circle((x + 0.14, y + h - 0.14), 0.14, facecolor=color, edgecolor=\"none\", zorder=3))\n", " ax.text(x + 0.14, y + h - 0.14, str(numero), ha=\"center\", va=\"center\", color=\"white\",\n", " fontsize=11, fontweight=\"bold\", zorder=4)\n", " ax.text(x + w / 2, y + h / 2 - 0.05, texto, ha=\"center\", va=\"center\", fontsize=8.3, zorder=4)\n", "\n", "def dibujar_pipeline():\n", " fig, ax = plt.subplots(figsize=(11, 6.5))\n", " ax.set_xlim(0, 11); ax.set_ylim(0, 7); ax.axis(\"off\")\n", " pasos = [\n", " (0.3, 4.6, \"Definir variables\\nde entrada y\\ndistribuciones\"),\n", " (3.9, 4.6, \"Crear modelo\\nmatemático /\\nfórmula\"),\n", " (7.5, 4.6, \"Generar\\nmuestras\\naleatorias\"),\n", " (7.5, 1.0, \"Calcular el\\nmodelo múltiples\\nveces (iterar)\"),\n", " (3.9, 1.0, \"Recopilar y\\nanalizar\\nresultados\"),\n", " (0.3, 1.0, \"Visualizar y\\ntomar\\ndecisiones\"),\n", " ]\n", " w, h = 3.0, 1.7\n", " for i, (x, y, t) in enumerate(pasos, 1):\n", " dibujar_paso(ax, x, y, w, h, i, t)\n", "\n", " def flecha(x1, y1, x2, y2):\n", " ax.add_patch(FancyArrowPatch((x1, y1), (x2, y2), arrowstyle=\"-|>\", mutation_scale=18,\n", " color=\"#B35A1F\", lw=2, zorder=1))\n", " flecha(3.3, 5.45, 3.9, 5.45)\n", " flecha(6.9, 5.45, 7.5, 5.45)\n", " flecha(9.0, 4.6, 9.0, 2.7)\n", " flecha(7.5, 1.85, 6.9, 1.85)\n", " flecha(3.9, 1.85, 3.3, 1.85)\n", "\n", " fig.suptitle(\"Las 6 etapas de una simulación de Monte Carlo\", fontsize=13, fontweight=\"bold\")\n", " fig.tight_layout()\n", " plt.show()\n", "\n", "dibujar_pipeline()\n" ] }, { "cell_type": "markdown", "id": "468e52d4", "metadata": {}, "source": [ "## 2. El ejemplo: estimar π lanzando dardos\n", "\n", "Imagina un cuadrado de lado 2 (área $=4$) con un círculo de radio 1 inscrito (área $=\\pi \\cdot 1^2 = \\pi$). Si lanzas dardos **uniformemente al azar** sobre el cuadrado, la probabilidad de que un dardo caiga dentro del círculo es\n", "\n", "$$P(\\text{dentro del círculo}) = \\frac{\\text{área del círculo}}{\\text{área del cuadrado}} = \\frac{\\pi}{4}$$\n", "\n", "Si lanzamos $n$ dardos y contamos cuántos caen dentro, podemos **despejar** $\\pi$:\n", "\n", "$$\\hat\\pi = 4 \\cdot \\frac{\\text{dardos dentro}}{n}$$\n", "\n", "Este es, literalmente, el objetivo general de la infografía: convertimos el problema \"¿cuánto vale $\\pi$?\" en el problema de estimar un valor esperado — aquí, $E[I] = P(\\text{dentro}) = \\pi/4$, donde $I$ es un indicador binario." ] }, { "cell_type": "markdown", "id": "d64cb0d8", "metadata": {}, "source": [ "## Paso 1 — Variables de entrada y distribuciones\n", "\n", "La posición de cada dardo se describe con dos variables aleatorias independientes:\n", "\n", "$$X \\sim \\text{Uniforme}(-1, 1), \\qquad Y \\sim \\text{Uniforme}(-1, 1)$$\n", "\n", "(el cuadrado va de $-1$ a $1$ en ambos ejes). No hay código en este paso — es puramente la definición conceptual del problema, antes de programar nada." ] }, { "cell_type": "markdown", "id": "b8274725", "metadata": {}, "source": [ "## Paso 2 — Crear el modelo matemático / fórmula \n", "\n", "El \"modelo\" aquí es solo una pregunta geométrica: ¿el punto $(x,y)$ cae dentro del círculo unitario? Un punto está dentro si $x^2+y^2 \\le 1$." ] }, { "cell_type": "code", "execution_count": null, "id": "18e35106", "metadata": {}, "outputs": [], "source": [ "def dentro_del_circulo(x, y):\n", " \"\"\"Retorna True (o un arreglo de booleanos) si (x, y) cae dentro del círculo unitario.\"\"\"\n", " # TODO: retorna la condición x**2 + y**2 <= 1\n", " raise NotImplementedError(\"Implementa dentro_del_circulo\")\n", "\n", "# --- prueba rápida (no la modifiques) ---\n", "assert dentro_del_circulo(0, 0) == True # el centro siempre está adentro\n", "assert dentro_del_circulo(1, 1) == False # esquina del cuadrado: está afuera\n", "assert dentro_del_circulo(1, 0) == True # justo sobre el borde del círculo\n", "print(\"OK: dentro_del_circulo pasa la prueba básica.\")\n" ] }, { "cell_type": "markdown", "id": "e1105fd5", "metadata": {}, "source": [ "## Paso 3 — Generar muestras aleatorias \n", "\n", "Igual que en la Semana 2, generamos $U \\sim U(0,1)$ y usamos la transformada inversa de la Uniforme para llevarla al rango $[-1,1]$: $X = -1 + 2U$." ] }, { "cell_type": "code", "execution_count": null, "id": "004ba343", "metadata": {}, "outputs": [], "source": [ "n = 10_000 # número de dardos (iteraciones) — el mismo valor de la infografía\n", "\n", "# TODO: genera U1 = np.random.random(n) y U2 = np.random.random(n)\n", "# TODO: calcula X = -1 + 2*U1 y Y = -1 + 2*U2 (transformada inversa de la Uniforme)\n", "\n", "\n", "print(\"Primeros 5 dardos (x, y):\")\n", "for i in range(5):\n", " print(f\" ({X[i]:.3f}, {Y[i]:.3f})\")\n", "\n", "assert X.shape == (n,) and Y.shape == (n,)\n", "assert np.all((X >= -1) & (X <= 1)) and np.all((Y >= -1) & (Y <= 1))\n", "print(\"OK: las coordenadas están dentro de [-1, 1].\")\n" ] }, { "cell_type": "markdown", "id": "d983eabe", "metadata": {}, "source": [ "## Paso 4 — Calcular el modelo múltiples veces (iterar) \n", "\n", "Aplica `dentro_del_circulo` a los `n` dardos **de una sola vez**, de forma vectorizada (sin ciclos de Python)." ] }, { "cell_type": "code", "execution_count": null, "id": "2a565b65", "metadata": {}, "outputs": [], "source": [ "# TODO: calcula I = dentro_del_circulo(X, Y).astype(float)\n", "\n", "\n", "print(f\"De {n:,} dardos, {int(I.sum()):,} cayeron dentro del círculo.\")\n", "assert I.shape == (n,)\n", "assert set(np.unique(I)) <= {0.0, 1.0}\n", "print(\"OK: el indicador solo toma valores 0 y 1.\")\n" ] }, { "cell_type": "markdown", "id": "c0ab4867", "metadata": {}, "source": [ "## Paso 5a — Recopilar y analizar resultados: más allá de la media \n", "\n", "La media (o proporción) es solo el punto de partida. Calcula también la dispersión, la precisión del estimador, el error contra el valor real (¡aquí lo conocemos!) y qué tan grande es esa incertidumbre en términos relativos." ] }, { "cell_type": "code", "execution_count": null, "id": "727ec795", "metadata": {}, "outputs": [], "source": [ "# TODO: p_hat = I.mean()\n", "# TODO: pi_hat = 4 * p_hat\n", "\n", "# TODO: s_I = I.std(ddof=1) # desviación estándar de los dardos individuales\n", "# TODO: SE_pi = 4 * s_I / np.sqrt(n) # error estándar del estimador pi_hat\n", "\n", "# TODO: ic_bajo = pi_hat - 1.96 * SE_pi\n", "# TODO: ic_alto = pi_hat + 1.96 * SE_pi\n", "\n", "# TODO: error_abs = abs(pi_hat - np.pi)\n", "# TODO: error_rel = error_abs / np.pi\n", "\n", "# TODO: CV = SE_pi / pi_hat\n", "\n", "# TODO: imprime las 8 métricas anteriores con una etiqueta clara para cada una\n" ] }, { "cell_type": "markdown", "id": "f66dd8a1", "metadata": {}, "source": [ "## Paso 5b — Cuando la mediana y los percentiles sí importan \n", "\n", "El indicador $I$ (0/1) es demasiado simple para que la mediana o los percentiles aporten algo nuevo. Pero cada dardo también tiene una **distancia al centro** $R=\\sqrt{X^2+Y^2}$, una variable **continua** — aquí sí tiene sentido comparar media vs. mediana, y mirar percentiles." ] }, { "cell_type": "code", "execution_count": null, "id": "26dae7c3", "metadata": {}, "outputs": [], "source": [ "# TODO: R = np.sqrt(X**2 + Y**2)\n", "\n", "# TODO: media_R = R.mean()\n", "# TODO: mediana_R = np.median(R)\n", "# TODO: std_R = R.std(ddof=1)\n", "# TODO: CV_R = std_R / media_R\n", "# TODO: percentiles_R = np.percentile(R, [5, 25, 50, 75, 95])\n", "\n", "# TODO: imprime las 5 métricas anteriores con una etiqueta clara para cada una\n" ] }, { "cell_type": "markdown", "id": "40e4bc0d", "metadata": {}, "source": [ "## Paso 5c — Tabla de referencia: ¿qué métrica uso y cuándo?\n", "\n", "| Métrica | Qué mide | ¿Cuándo usarla? |\n", "|---|---|---|\n", "| **Media / proporción** ($\\bar Y$, $\\hat p$) | El valor \"central\" esperado del resultado | Siempre es el punto de partida: resume el resultado típico |\n", "| **Desviación estándar** ($s$) | Qué tan dispersos son los resultados **individuales** | Para entender el riesgo/variabilidad intrínseca del sistema (¿qué tan impredecible es UNA corrida?) |\n", "| **Error estándar** ($SE = s/\\sqrt n$) | Qué tan preciso es el **estimador** de la media (no los datos crudos) | Para decidir cuántas réplicas $n$ hacen falta, y para reportar la incertidumbre de la estimación |\n", "| **Intervalo de confianza** ($\\bar Y \\pm 1.96\\, SE$) | Rango donde se espera esté el valor real, con 95% de confianza | Para comunicar resultados a quien toma decisiones, con un nivel de certeza explícito |\n", "| **Error relativo vs. valor real** | Qué tan lejos está la estimación del valor correcto | Solo cuando se conoce la respuesta real — para *validar* el método antes de aplicarlo a un caso donde no se conoce |\n", "| **Coeficiente de variación** ($CV = SE/\\bar Y$) | Dispersión **relativa**, sin unidades | Para comparar variabilidad entre escenarios distintos, o entre variables con escalas diferentes |\n", "| **Mediana / percentiles** | El centro \"robusto\" y los casos extremos de la distribución completa | Cuando la distribución es asimétrica o tiene colas pesadas (tiempos de espera, pérdidas), o interesa el peor/mejor caso — **no** aportan nada nuevo para variables binarias como $I$ |" ] }, { "cell_type": "markdown", "id": "f6896c49", "metadata": {}, "source": [ "## 10. El Teorema del Límite Central — ¿por qué funciona el intervalo de confianza?\n", "\n", "El intervalo de confianza del Paso 5a ($\\hat\\pi \\pm 1.96\\cdot SE$) asume que $\\hat\\pi$ se distribuye, aproximadamente, como una Normal. Pero cada dardo individual es apenas una variable binaria (0/1) — ¡nada parecido a una campana de Gauss! ¿Por qué entonces podemos usar la fórmula de la Normal?\n", "\n", "La respuesta es el **Teorema del Límite Central (TLC)**: si $I_1,\\dots,I_n$ son variables aleatorias independientes e idénticamente distribuidas, con media $\\mu$ y varianza $\\sigma^2$ finita, entonces para $n$ suficientemente grande\n", "\n", "$$\\bar I_n = \\frac{1}{n}\\sum_{i=1}^n I_i \\ \\approx\\ \\text{Normal}\\left(\\mu,\\ \\frac{\\sigma^2}{n}\\right)$$\n", "\n", "**sin importar la forma de la distribución original de cada $I_i$.** Como $\\hat\\pi = 4\\bar I_n$ es una transformación lineal de $\\bar I_n$, también es aproximadamente Normal. Esto es lo que hace posible la simulación de Monte Carlo en general: no importa qué tan \"fea\" o discreta sea la salida de un solo replicado — el **promedio de muchos replicados** se comporta de forma predecible (Normal), y eso es lo que permite construir intervalos de confianza válidos.\n", "\n", "Vamos a comprobarlo empíricamente: en lugar de correr el experimento una sola vez con $n=10\\,000$ dardos, lo corremos $M$ veces de forma **independiente** (con menos dardos por corrida, para que la variabilidad entre corridas sea visible) y graficamos el histograma de las $M$ estimaciones de $\\hat\\pi$ obtenidas, junto con la Normal teórica que predice el TLC." ] }, { "cell_type": "code", "execution_count": null, "id": "847ef069", "metadata": {}, "outputs": [], "source": [ "M = 3000 # número de réplicas independientes del experimento completo\n", "n_rep = 1000 # dardos por réplica (menos que los 10,000 de antes, para ver la variabilidad)\n", "\n", "# TODO: genera U1_r, U2_r de forma (M, n_rep) con np.random.random((M, n_rep))\n", "# TODO: calcula Xr = -1 + 2*U1_r, Yr = -1 + 2*U2_r\n", "# TODO: calcula Ir = dentro_del_circulo(Xr, Yr).astype(float)\n", "# TODO: calcula pi_hats = 4 * Ir.mean(axis=1) -> un pi_hat por cada una de las M réplicas\n", "\n", "\n", "# TODO: calcula media_pi_hats y std_pi_hats (media y std de pi_hats, ddof=1)\n", "\n", "# TODO: calcula la normal teórica predicha por el TLC:\n", "# p_teorico = np.pi / 4\n", "# sigma_teorico = 4 * np.sqrt(p_teorico * (1 - p_teorico) / n_rep)\n", "\n", "# TODO: imprime media_pi_hats, std_pi_hats y sigma_teorico, comparándolos\n", "\n", "\n", "# TODO: arma un histograma (density=True) de pi_hats, superpón la curva Normal teórica\n", "# (usa la fórmula de la densidad Normal con media=np.pi y desviación=sigma_teorico)\n", "# y una línea vertical en np.pi\n" ] }, { "cell_type": "markdown", "id": "c56e6c6e", "metadata": {}, "source": [ "## Paso 6 — Visualizar y tomar decisiones: los dardos en vivo (dado)\n", "\n", "Con todo lo anterior implementado, animamos los dardos cayendo sobre el cuadrado (azules = dentro del círculo, rojos = fuera) mientras, en el panel derecho, $\\hat\\pi(n)$ converge al valor real de $\\pi$." ] }, { "cell_type": "code", "execution_count": null, "id": "675f6bcb", "metadata": {}, "outputs": [], "source": [ "N_ANIM = 2500\n", "frames_n = list(range(10, 101, 10)) + list(range(120, 501, 20)) + list(range(550, N_ANIM + 1, 50))\n", "\n", "pi_hat_n = 4 * np.cumsum(I[:N_ANIM]) / np.arange(1, N_ANIM + 1)\n", "\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), dpi=70)\n", "ax1.add_patch(Rectangle((-1, -1), 2, 2, fill=False, edgecolor=\"black\", lw=1.2))\n", "ax1.add_patch(Circle((0, 0), 1, fill=False, edgecolor=\"black\", lw=1.2))\n", "scat_in = ax1.scatter([], [], s=5, color=\"#2D9CDB\", label=\"dentro\")\n", "scat_out = ax1.scatter([], [], s=5, color=\"#EB5757\", label=\"fuera\")\n", "ax1.set_xlim(-1.05, 1.05); ax1.set_ylim(-1.05, 1.05)\n", "ax1.set_aspect(\"equal\"); ax1.axis(\"off\")\n", "ax1.legend(loc=\"lower center\", ncol=2, fontsize=7)\n", "\n", "line2, = ax2.plot([], [], color=\"#4C72B0\", lw=1.3)\n", "ax2.axhline(np.pi, color=\"black\", ls=\"--\", lw=1, label=f\"π real = {np.pi:.4f}\")\n", "ax2.set_xlim(1, N_ANIM); ax2.set_ylim(2.6, 3.6)\n", "ax2.set_xlabel(\"n (dardos)\"); ax2.set_ylabel(\"π̂(n)\")\n", "ax2.legend(fontsize=7)\n", "\n", "titulo = fig.suptitle(\"\")\n", "\n", "def actualizar(i):\n", " k = frames_n[i]\n", " mask_in = I[:k] == 1\n", " scat_in.set_offsets(np.column_stack([X[:k][mask_in], Y[:k][mask_in]]))\n", " scat_out.set_offsets(np.column_stack([X[:k][~mask_in], Y[:k][~mask_in]]))\n", " line2.set_data(np.arange(1, k + 1), pi_hat_n[:k])\n", " titulo.set_text(f\"n = {k} π̂(n) = {pi_hat_n[k - 1]:.4f}\")\n", " return []\n", "\n", "ani = animation.FuncAnimation(fig, actualizar, frames=len(frames_n), blit=False)\n", "plt.close(fig)\n", "HTML(ani.to_jshtml(fps=10))\n" ] }, { "cell_type": "markdown", "id": "54dd897f", "metadata": {}, "source": [ "## Los 4 componentes de la simulación, aplicados a este ejemplo\n", "\n", "| Componente | En este ejemplo |\n", "|---|---|\n", "| **1. Generador de números pseudoaleatorios (LCG)** | `np.random.random()` — en la Semana 2 implementaste tu propio generador congruencial; aquí usamos el de NumPy, que sigue la misma idea |\n", "| **2. Definición de distribuciones de entrada** | $X, Y \\sim \\text{Uniforme}(-1,1)$ (Paso 1) |\n", "| **3. Generadores de variables aleatorias** (transformada inversa / aceptación-rechazo / Box-Muller) | Transformada inversa de la Uniforme: $X=-1+2U$ (Paso 3) — el caso más simple de lo que viste en la Semana 2 |\n", "| **4. Estimador de respuesta** | $\\hat\\pi = 4\\cdot\\bar I_n$, con su error estándar e intervalo de confianza justificados por el TLC (Pasos 5 y 10) |" ] }, { "cell_type": "markdown", "id": "cba684e6", "metadata": {}, "source": [ "## Preguntas de discusión\n", "\n", "1. En el Paso 5b, la media y la mediana de $R$ (distancia al centro) no coinciden exactamente. ¿Qué te dice eso sobre la forma de la distribución de $R$? ¿Por qué esta pregunta no tendría sentido si la hiciéramos sobre el indicador $I$ (0/1) del Paso 5a?\n", "2. Si quisieras reportar \"qué tan riesgoso/variable es un solo lanzamiento de dardo\", ¿usarías la desviación estándar de $I$ o el error estándar $SE(\\hat\\pi)$? Justifica la diferencia entre ambas.\n", "3. En la sección 10 (TLC), ¿qué le pasaría al ancho del histograma de las $M$ estimaciones si aumentaras `n_rep` (dardos por réplica) de 1000 a 10 000? ¿Y si mantienes `n_rep` fijo pero aumentas $M$?\n", "4. ¿Por qué es válido usar la fórmula del intervalo de confianza ($\\pm 1.96\\cdot SE$) aunque cada dardo individual sea una variable binaria y no una Normal?\n", "5. Compara este ejemplo con el sistema de 5 componentes de la semana pasada: ¿qué métrica de las presentadas aquí sería más útil para reportar el tiempo de falla de ese sistema, y por qué?" ] }, { "cell_type": "markdown", "id": "96d23ec8", "metadata": {}, "source": [ "_Escribe aquí tus respuestas a las 5 preguntas anteriores._" ] }, { "cell_type": "markdown", "id": "fb04b950", "metadata": {}, "source": [ "## Conclusiones\n", "\n", "_Escribe 2-3 conclusiones sobre lo observado en este notebook._" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }