{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3884a69a",
   "metadata": {},
   "source": [
    "# Semana 2 — Generación de Variables Aleatorias: Transformada Inversa y Aceptación-Rechazo\n",
    "\n",
    "**Objetivo.** Derivar e implementar el método de la **transformada inversa** para distribuciones continuas de uso común (Uniforme, Exponencial, Weibull, Rayleigh, Logística, Pareto y Triangular) y el método de **aceptación y rechazo** para una densidad polinómica y una distribución Beta, validando cada implementación mediante histogramas contra la densidad teórica.\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15be43d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "np.random.seed(2026)\n",
    "N = 50_000  # tamaño de cada muestra a generar\n",
    "\n",
    "resultados = {}  # aquí guardamos (x, pdf, rango) de cada distribución de la Parte A\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02c55718",
   "metadata": {},
   "source": [
    "## 0. Método de la transformada inversa\n",
    "\n",
    "Si $X$ es una variable aleatoria continua con función de distribución acumulada (CDF) $F$ estrictamente creciente e invertible, y $U \\sim \\text{Uniforme}(0,1)$, entonces\n",
    "\n",
    "$$X = F^{-1}(U)$$\n",
    "\n",
    "tiene exactamente la distribución de $X$ (teorema de la transformada inversa).\n",
    "\n",
    "**Receta general** para simular $X \\sim F$:\n",
    "1. Generar $U \\sim \\text{Uniforme}(0,1)$.\n",
    "2. Despejar $x$ de $u = F(x)$ para obtener $F^{-1}(u)$ en forma cerrada.\n",
    "3. Calcular $X = F^{-1}(U)$.\n",
    "\n",
    "Este método es exacto y muy eficiente **siempre que $F^{-1}$ tenga forma cerrada**. En este notebook vas a aplicarlo a 7 distribuciones continuas, en orden de dificultad creciente. Para cada una vamos a graficar dos histogramas lado a lado: el de los números pseudoaleatorios base $U(0,1)$ y el de la variable transformada $X = F^{-1}(U)$, cada uno con su densidad teórica superpuesta."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88eeca2b",
   "metadata": {},
   "source": [
    "## 1. Herramienta de visualización reutilizable \n",
    "\n",
    "Construye una función auxiliar que vas a reutilizar en todo el notebook: dibuja el histograma normalizado de una muestra junto con la curva de la densidad teórica correspondiente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab445580",
   "metadata": {},
   "outputs": [],
   "source": [
    "def histograma_vs_teorica(muestras, pdf_teorica, rango, nombre, bins=50, color=\"#4C72B0\", ax=None):\n",
    "    \"\"\"\n",
    "    Dibuja el histograma normalizado (density=True) de `muestras`, restringido\n",
    "    a `rango`, superponiendo la curva de la densidad teórica `pdf_teorica`\n",
    "    (función vectorizada de x).\n",
    "\n",
    "    Si `ax` es None crea una figura nueva; si se pasa un `ax` existente,\n",
    "    dibuja sobre él (útil para paneles con varias subgráficas).\n",
    "    \"\"\"\n",
    "    # TODO 1: si ax es None, crea (fig, ax) nuevos con plt.subplots(figsize=(5, 3.5))\n",
    "    # TODO 2: dibuja el histograma de `muestras` en `ax` (bins=bins, range=rango, density=True)\n",
    "    # TODO 3: crea xs = np.linspace(rango[0], rango[1], 400) y dibuja pdf_teorica(xs) sobre ax\n",
    "    # TODO 4: pon título=nombre, etiqueta los ejes (\"x\", \"densidad\") y una leyenda\n",
    "    # TODO 5: retorna ax\n",
    "    raise NotImplementedError(\"Implementa histograma_vs_teorica\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "u_demo = np.random.random(5000)\n",
    "histograma_vs_teorica(u_demo, pdf_teorica=lambda t: np.ones_like(t), rango=(0, 1),\n",
    "                       nombre=\"Prueba: U(0,1)\")\n",
    "plt.show()\n",
    "print(\"OK: si ves un histograma aproximadamente plano en 1.0, la función funciona.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "767525c0",
   "metadata": {},
   "source": [
    "## 2. Uniforme(a, b) \n",
    "\n",
    "$$f(x) = \\frac{1}{b-a}, \\ a \\le x \\le b \\qquad F(x) = \\frac{x-a}{b-a} \\qquad F^{-1}(u) = a + (b-a)\\,u$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a84d38",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_uniforme(u, a, b):\n",
    "    \"\"\"F^{-1}(u) = a + (b - a) * u\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "a_unif, b_unif = 2.0, 8.0\n",
    "u = np.random.random(N)\n",
    "x = F_inv_uniforme(u, a_unif, b_unif)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))\n",
    "histograma_vs_teorica(u, pdf_teorica=lambda t: np.ones_like(t), rango=(0, 1),\n",
    "                       nombre=\"U(0,1) — números base\", ax=axes[0])\n",
    "histograma_vs_teorica(x, pdf_teorica=lambda t: np.full_like(t, 1 / (b_unif - a_unif)),\n",
    "                       rango=(a_unif, b_unif), nombre=f\"Uniforme({a_unif}, {b_unif})\", ax=axes[1])\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "resultados[\"Uniforme(a,b)\"] = dict(\n",
    "    x=x, pdf=lambda t: np.full_like(t, 1 / (b_unif - a_unif)), rango=(a_unif, b_unif))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ee1dc38",
   "metadata": {},
   "source": [
    "## 3. Exponencial(λ) \n",
    "\n",
    "$$f(x) = \\lambda e^{-\\lambda x}, \\ x \\ge 0 \\qquad F(x) = 1 - e^{-\\lambda x} \\qquad F^{-1}(u) = -\\frac{\\ln(1-u)}{\\lambda}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58e54ca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_exponencial(u, lam):\n",
    "    \"\"\"F^{-1}(u) = -ln(1 - u) / lam\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "lam = 2.0\n",
    "u = np.random.random(N)\n",
    "x = F_inv_exponencial(u, lam)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))\n",
    "histograma_vs_teorica(u, pdf_teorica=lambda t: np.ones_like(t), rango=(0, 1),\n",
    "                       nombre=\"U(0,1) — números base\", ax=axes[0])\n",
    "histograma_vs_teorica(x, pdf_teorica=lambda t: lam * np.exp(-lam * t), rango=(0, 5 / lam),\n",
    "                       nombre=f\"Exponencial(λ={lam})\", ax=axes[1])\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "resultados[\"Exponencial(λ)\"] = dict(x=x, pdf=lambda t: lam * np.exp(-lam * t), rango=(0, 5 / lam))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e7f51eb",
   "metadata": {},
   "source": [
    "## 4. Weibull(α, β) \n",
    "\n",
    "$$f(x) = \\frac{\\alpha}{\\beta}\\left(\\frac{x}{\\beta}\\right)^{\\alpha-1} e^{-(x/\\beta)^\\alpha}, \\ x \\ge 0 \\qquad F(x) = 1 - e^{-(x/\\beta)^\\alpha} \\qquad F^{-1}(u) = \\beta\\big(-\\ln(1-u)\\big)^{1/\\alpha}$$\n",
    "\n",
    "Esta vez también te toca escribir tú las líneas que generan `u` y `x` (la función de densidad teórica y la figura ya están dadas)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8510406b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_weibull(u, alpha, beta):\n",
    "    \"\"\"F^{-1}(u) = beta * (-ln(1-u))**(1/alpha)\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "alpha_w, beta_w = 2.0, 3.0\n",
    "# TODO: genera u = np.random.random(N)\n",
    "# TODO: calcula x = F_inv_weibull(u, alpha_w, beta_w)\n",
    "\n",
    "\n",
    "def pdf_weibull(t, alpha=alpha_w, beta=beta_w):\n",
    "    return (alpha / beta) * (t / beta) ** (alpha - 1) * np.exp(-(t / beta) ** alpha)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))\n",
    "histograma_vs_teorica(u, pdf_teorica=lambda t: np.ones_like(t), rango=(0, 1),\n",
    "                       nombre=\"U(0,1) — números base\", ax=axes[0])\n",
    "histograma_vs_teorica(x, pdf_teorica=pdf_weibull, rango=(0, 8),\n",
    "                       nombre=f\"Weibull(α={alpha_w}, β={beta_w})\", ax=axes[1])\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "resultados[\"Weibull(α,β)\"] = dict(x=x, pdf=pdf_weibull, rango=(0, 8))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4348e61b",
   "metadata": {},
   "source": [
    "## 5. Rayleigh(σ) \n",
    "\n",
    "$$f(x) = \\frac{x}{\\sigma^2} e^{-x^2/(2\\sigma^2)}, \\ x \\ge 0 \\qquad F(x) = 1 - e^{-x^2/(2\\sigma^2)} \\qquad F^{-1}(u) = \\sigma\\sqrt{-2\\ln(1-u)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3012966b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_rayleigh(u, sigma):\n",
    "    \"\"\"F^{-1}(u) = sigma * sqrt(-2 ln(1-u))\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "sigma_r = 2.0\n",
    "# TODO: genera u = np.random.random(N)\n",
    "# TODO: calcula x = F_inv_rayleigh(u, sigma_r)\n",
    "\n",
    "\n",
    "def pdf_rayleigh(t, sigma=sigma_r):\n",
    "    return (t / sigma ** 2) * np.exp(-t ** 2 / (2 * sigma ** 2))\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))\n",
    "histograma_vs_teorica(u, pdf_teorica=lambda t: np.ones_like(t), rango=(0, 1),\n",
    "                       nombre=\"U(0,1) — números base\", ax=axes[0])\n",
    "histograma_vs_teorica(x, pdf_teorica=pdf_rayleigh, rango=(0, 8),\n",
    "                       nombre=f\"Rayleigh(σ={sigma_r})\", ax=axes[1])\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "resultados[\"Rayleigh(σ)\"] = dict(x=x, pdf=pdf_rayleigh, rango=(0, 8))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c6edeb9",
   "metadata": {},
   "source": [
    "## 6. Logística(μ, s) \n",
    "\n",
    "$$f(x) = \\frac{e^{-(x-\\mu)/s}}{s\\left(1+e^{-(x-\\mu)/s}\\right)^2} \\qquad F(x) = \\frac{1}{1+e^{-(x-\\mu)/s}} \\qquad F^{-1}(u) = \\mu + s\\ln\\!\\left(\\frac{u}{1-u}\\right)$$\n",
    "\n",
    "Ahora, además de `F_inv_logistica` y de generar `u`, `x`, te toca armar tú la figura 1x2 completa y llamar a `histograma_vs_teorica` para ambos histogramas (usa como referencia las secciones anteriores)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76fc6e4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_logistica(u, mu, s):\n",
    "    \"\"\"F^{-1}(u) = mu + s * ln(u / (1 - u))\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "mu_l, s_l = 0.0, 1.0\n",
    "# TODO: genera u y calcula x = F_inv_logistica(u, mu_l, s_l)\n",
    "\n",
    "\n",
    "def pdf_logistica(t, mu=mu_l, s=s_l):\n",
    "    z = np.exp(-(t - mu) / s)\n",
    "    return z / (s * (1 + z) ** 2)\n",
    "\n",
    "# TODO: arma una figura 1x2 (fig, axes = plt.subplots(1, 2, figsize=(10, 3.5)))\n",
    "# TODO: llama histograma_vs_teorica para `u` (rango (0,1), pdf uniforme) en axes[0]\n",
    "# TODO: llama histograma_vs_teorica para `x` (rango (-8,8), pdf_teorica=pdf_logistica) en axes[1]\n",
    "# TODO: fig.tight_layout() y plt.show()\n",
    "\n",
    "\n",
    "resultados[\"Logística(μ,s)\"] = dict(x=x, pdf=pdf_logistica, rango=(-8, 8))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9155ee3",
   "metadata": {},
   "source": [
    "## 7. Triangular(a, b, c) \n",
    "\n",
    "La triangular es la primera distribución **por partes** del notebook: su CDF cambia de fórmula en $x=c$ (la moda).\n",
    "\n",
    "$$\n",
    "F(x) = \\begin{cases}\n",
    "\\dfrac{(x-a)^2}{(b-a)(c-a)}, & a \\le x \\le c \\\\[6pt]\n",
    "1 - \\dfrac{(b-x)^2}{(b-a)(b-c)}, & c < x \\le b\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Invirtiendo cada tramo (evaluado en el punto de quiebre $u^* = F(c) = \\frac{c-a}{b-a}$):\n",
    "\n",
    "$$\n",
    "F^{-1}(u) = \\begin{cases}\n",
    "a + \\sqrt{u\\,(b-a)(c-a)}, & u \\le u^* \\\\[6pt]\n",
    "b - \\sqrt{(1-u)\\,(b-a)(b-c)}, & u > u^*\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Implementa esto de forma **vectorizada** (sin ciclos de Python) usando `np.where`. Esta sección está casi completamente abierta: solo tienes la densidad teórica de referencia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d21a057d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_inv_triangular(u, a, b, c):\n",
    "    \"\"\"Inversa por partes de la CDF triangular(a, b, c), vectorizada con np.where.\"\"\"\n",
    "    # TODO 1: calcula el punto de quiebre u_estrella = (c - a) / (b - a)\n",
    "    # TODO 2: retorna np.where(u <= u_estrella, <tramo izquierdo>, <tramo derecho>)\n",
    "    #         usando las fórmulas de la celda de teoría\n",
    "    raise NotImplementedError(\"Implementa F_inv_triangular\")\n",
    "\n",
    "a_t, b_t, c_t = 0.0, 10.0, 3.0\n",
    "# TODO: genera u y calcula x = F_inv_triangular(u, a_t, b_t, c_t)\n",
    "\n",
    "\n",
    "def pdf_triangular(t, a=a_t, b=b_t, c=c_t):\n",
    "    t = np.asarray(t, dtype=float)\n",
    "    izquierda = 2 * (t - a) / ((b - a) * (c - a))\n",
    "    derecha = 2 * (b - t) / ((b - a) * (b - c))\n",
    "    return np.where(t < c, izquierda, derecha)\n",
    "\n",
    "# TODO: arma una figura 1x2 y llama histograma_vs_teorica para `u` (rango (0,1))\n",
    "#       y para `x` (rango (a_t, b_t), pdf_teorica=pdf_triangular)\n",
    "\n",
    "\n",
    "resultados[\"Triangular(a,b,c)\"] = dict(x=x, pdf=pdf_triangular, rango=(a_t, b_t))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde56374",
   "metadata": {},
   "source": [
    "## 8. Bonus — Panel resumen (Nivel avanzado, opcional)\n",
    "\n",
    "Arma una única figura 2x4 que muestre, de un solo vistazo, el histograma final (con su densidad teórica) de las 7 distribuciones trabajadas. Vuelve a usar `histograma_vs_teorica`, esta vez pasándole siempre un `ax` de una grilla de subplots — no hace falta regenerar ninguna muestra, todo está guardado en `resultados`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d420703",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 4, figsize=(16, 7))\n",
    "# TODO: recorre resultados.items() en paralelo con axes.flat (zip) y llama\n",
    "#       histograma_vs_teorica(r[\"x\"], pdf_teorica=r[\"pdf\"], rango=r[\"rango\"], nombre=nombre, ax=ax)\n",
    "# TODO: apaga (ax.axis(\"off\")) los ejes sobrantes que no se usaron\n",
    "fig.suptitle(\"Resumen — Transformada inversa aplicada a 7 distribuciones\", y=1.02)\n",
    "fig.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "257a256a",
   "metadata": {},
   "source": [
    "## 10.Método de aceptación y rechazo\n",
    "\n",
    "Cuando $F^{-1}$ no tiene forma cerrada (o es muy costosa de evaluar), se puede simular $X \\sim f$ usando una densidad auxiliar $g$ (fácil de simular) y una constante $c$ tal que\n",
    "\n",
    "$$f(x) \\le c\\, g(x) \\quad \\text{para todo } x \\text{ en el soporte.}$$\n",
    "\n",
    "**Algoritmo:**\n",
    "1. Generar $Y \\sim g$ (la \"propuesta\").\n",
    "2. Generar $U \\sim \\text{Uniforme}(0,1)$, independiente de $Y$.\n",
    "3. Si $U \\le \\dfrac{f(Y)}{c\\,g(Y)}$: aceptar $X = Y$. Si no, volver al paso 1.\n",
    "\n",
    "La probabilidad de aceptar en cada intento es $1/c$: entre más ajustada sea la cota $c\\,g(x)$ alrededor de $f(x)$, más eficiente es el método. En los dos casos de esta parte usamos como propuesta $g(x) = \\text{Uniforme}(0,1)$, la más simple posible, porque ambas densidades objetivo tienen soporte $[0,1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da53904a",
   "metadata": {},
   "source": [
    "## 11. Implementación genérica de aceptación-rechazo \n",
    "\n",
    "Implementa una función reutilizable para los dos casos siguientes. Por eficiencia, genera los candidatos **por lotes** (arreglos de NumPy) en vez de uno a la vez con un ciclo de Python puro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78812707",
   "metadata": {},
   "outputs": [],
   "source": [
    "def aceptacion_rechazo(f, c, rango, n, tam_lote=5000):\n",
    "    \"\"\"\n",
    "    Genera n muestras de la densidad objetivo f(x) en `rango`, usando\n",
    "    aceptación-rechazo con propuesta g(x) = Uniforme(rango) y cota c tal\n",
    "    que f(x) <= c * g(x) para todo x en rango.\n",
    "\n",
    "    Retorna (muestras, tasa_aceptacion_empirica).\n",
    "    \"\"\"\n",
    "    # TODO 1: calcula g = 1 / (b - a), la densidad de la propuesta Uniforme(rango)\n",
    "    # TODO 2: inicializa una lista `aceptados`, y contadores total_aceptado=0, total_propuesto=0\n",
    "    # TODO 3: en un ciclo while (mientras total_aceptado < n):\n",
    "    #           - genera un lote de propuestas y = np.random.uniform(a, b, size=tam_lote)\n",
    "    #           - genera u = np.random.random(tam_lote)\n",
    "    #           - calcula el criterio de aceptación: u <= f(y) / (c * g)\n",
    "    #           - agrega y[criterio] a `aceptados` y actualiza los contadores\n",
    "    # TODO 4: concatena `aceptados` con np.concatenate y recorta a los primeros n\n",
    "    # TODO 5: calcula tasa = total_aceptado / total_propuesto\n",
    "    # TODO 6: retorna (muestras, tasa)\n",
    "    raise NotImplementedError(\"Implementa aceptacion_rechazo\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "np.random.seed(0)\n",
    "m_test, tasa_test = aceptacion_rechazo(f=lambda x: np.ones_like(x), c=1.0, rango=(0, 1), n=2000)\n",
    "assert len(m_test) == 2000\n",
    "assert 0.8 <= tasa_test <= 1.0  # con f = g, la tasa esperada es 1.0\n",
    "np.random.seed(2026)  # restaura la semilla del curso\n",
    "print(\"OK: aceptacion_rechazo pasa la prueba básica. tasa =\", round(tasa_test, 3))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ef175ba",
   "metadata": {},
   "source": [
    "## 12. Caso A — Densidad polinómica \n",
    "\n",
    "$$f(x) = 6x(1-x), \\quad 0 \\le x \\le 1$$\n",
    "\n",
    "Es una densidad válida ($\\int_0^1 6x(1-x)\\,dx = 1$) y, al ser un polinomio, no hace falta ninguna distribución \"con nombre\" para definirla. Como propuesta usamos $g(x) = \\text{Uniforme}(0,1)$; la cota $c$ se halla maximizando $f$ analíticamente:\n",
    "\n",
    "$$f'(x) = 6 - 12x = 0 \\;\\Rightarrow\\; x^* = 0.5, \\qquad c = f(0.5) = 1.5$$\n",
    "\n",
    "(Curiosidad: esta densidad coincide con una Beta(2,2) — pero aquí la tratamos como lo que es, un polinomio simple, sin usar `scipy.stats.beta`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fc29c4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_polinomica(x):\n",
    "    \"\"\"f(x) = 6x(1-x), 0 <= x <= 1\"\"\"\n",
    "    # TODO: implementa la fórmula anterior\n",
    "    raise NotImplementedError\n",
    "\n",
    "c_poli = None  # TODO: reemplaza con la cota c hallada analíticamente en la celda de teoría\n",
    "\n",
    "x_poli, tasa_poli = aceptacion_rechazo(f_polinomica, c_poli, rango=(0, 1), n=N)\n",
    "print(f\"Tasa de aceptación empírica: {tasa_poli:.4f}  (esperada: {1/c_poli:.4f})\")\n",
    "\n",
    "u_prop = np.random.uniform(0, 1, N)  # solo para ilustrar el histograma de propuestas\n",
    "\n",
    "# TODO: arma una figura 1x2 y llama histograma_vs_teorica para `u_prop` (rango (0,1))\n",
    "#       y para `x_poli` (rango (0,1), pdf_teorica=f_polinomica)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "018d50f8",
   "metadata": {},
   "source": [
    "## 13. Caso B — Beta(α, β) \n",
    "\n",
    "Ahora usamos una Beta(α=2.5, β=3.5) genuina (parámetros no enteros), cuya CDF involucra la función beta incompleta y **no tiene inversa en forma cerrada** — este es exactamente el tipo de caso donde el método de aceptación y rechazo es indispensable (a diferencia de todas las distribuciones de la Parte A).\n",
    "\n",
    "Usa `scipy.stats.beta.pdf(x, a, b)` como $f$, propuesta $g(x)=\\text{Uniforme}(0,1)$, y encuentra la cota $c = \\max_{x\\in[0,1]} f(x)$ **numéricamente**: evalúa $f$ sobre una grilla fina de puntos en $[0,1]$ y toma el máximo (no hace falta resolver ninguna derivada a mano esta vez). Esta sección está casi completamente abierta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5c0fb52",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_beta, beta_beta = 2.5, 3.5\n",
    "\n",
    "def f_beta(x):\n",
    "    # TODO: retorna stats.beta.pdf(x, alpha_beta, beta_beta)\n",
    "    raise NotImplementedError\n",
    "\n",
    "# TODO: crea una grilla fina de puntos en [0,1] (p.ej. np.linspace(0, 1, 20_001))\n",
    "# TODO: calcula c_beta = f_beta(grilla).max()\n",
    "# TODO: imprime la cota c_beta encontrada\n",
    "\n",
    "# TODO: llama x_beta, tasa_beta = aceptacion_rechazo(f_beta, c_beta, rango=(0, 1), n=N)\n",
    "# TODO: imprime la tasa de aceptación empírica y la esperada (1/c_beta)\n",
    "\n",
    "u_prop = np.random.uniform(0, 1, N)\n",
    "\n",
    "# TODO: arma una figura 1x2 y llama histograma_vs_teorica para `u_prop` (rango (0,1))\n",
    "#       y para `x_beta` (rango (0,1), pdf_teorica=f_beta)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec7e6c2",
   "metadata": {},
   "source": [
    "## 14. Preguntas de discusión\n",
    "\n",
    "Responde en la celda de texto siguiente, apoyándote en los resultados de las secciones anteriores.\n",
    "\n",
    "1. De las 7 distribuciones de la Parte A, ¿en cuál fue más difícil obtener $F^{-1}$ en forma cerrada y qué tuviste que hacer distinto para implementarla, comparado con las demás?\n",
    "2. ¿Por qué no se puede usar el método de la transformada inversa para la Beta(2.5, 3.5) de la Parte B?\n",
    "3. En el caso polinómico, ¿qué pasaría con la tasa de aceptación si usaras una cota $c$ más grande de la necesaria (por ejemplo $c=3$ en lugar de $1.5$)? ¿Seguiría siendo válido el método?\n",
    "4. Comparando el histograma de las propuestas $Y\\sim U(0,1)$ con el histograma de las muestras aceptadas en los dos casos de la Parte B, ¿qué le hace el filtro de aceptación-rechazo a la forma de la distribución?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f14ee0bc",
   "metadata": {},
   "source": [
    "_Escribe aquí tus respuestas a las 4 preguntas anteriores._"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
