{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9f1071c5",
   "metadata": {},
   "source": [
    "# Semana 2 — Generadores de números pseudoaleatorios U(0,1): LCG\n",
    "\n",
    "**Objetivo.** Implementar y comparar generadores congruenciales lineales — método **multiplicativo** y método **mixto** — cada uno con un conjunto de parámetros \"bueno\" y uno \"malo\", generando una secuencia de $N = 100\\,000$ números pseudoaleatorios $U(0,1)$ por caso. Evaluar la calidad de cada generador mediante su **histograma** y su **autocorrelograma**.\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a30423de",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "np.random.seed(2026)  # fija la semilla del curso para reproducibilidad general del notebook\n",
    "\n",
    "N = 1_000    # tamaño de cada secuencia pseudoaleatoria a generar\n",
    "MAXLAG = 3000  # rezago máximo a evaluar en el autocorrelograma\n",
    "BINS = 50      # número de bins para el histograma \"estándar\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d541eb63",
   "metadata": {},
   "source": [
    "## 0. Marco teórico — Generadores congruenciales lineales\n",
    "\n",
    "Un **generador congruencial lineal (LCG)** produce una secuencia de enteros $X_0, X_1, X_2, \\dots$ mediante la recurrencia\n",
    "\n",
    "$$X_{i+1} = (a X_i + c) \\bmod m, \\qquad U_i = \\frac{X_i}{m} \\in [0, 1)$$\n",
    "\n",
    "donde $a$ es el **multiplicador**, $c$ es el **incremento**, $m$ es el **módulo** y $X_0$ es la **semilla**. Según el valor de $c$ se distinguen dos casos particulares que vamos a implementar hoy:\n",
    "\n",
    "- **Método congruencial multiplicativo** ($c = 0$): $X_{i+1} = (a X_i) \\bmod m$. El periodo máximo posible es $m - 1$ (se alcanza si $m$ es primo y $a$ es raíz primitiva módulo $m$), o a lo sumo $m/4$ si $m$ es una potencia de 2.\n",
    "- **Método congruencial mixto** ($c \\neq 0$): $X_{i+1} = (a X_i + c) \\bmod m$. Por el **teorema de Hull–Dobell**, este generador alcanza periodo completo $m$ si y solo si: (1) $\\gcd(c, m) = 1$; (2) $a - 1$ es divisible por cada factor primo de $m$; (3) $a - 1$ es divisible por 4 si $m$ es divisible por 4.\n",
    "\n",
    "Un generador \"bueno\" **no** es simplemente el que produce un histograma que luce uniforme: también necesita **periodo suficientemente largo** (mucho mayor que la cantidad de números que se van a usar en la simulación) y **baja correlación serial** entre valores separados por distintos rezagos. Un generador \"malo\" puede tener un histograma marginal que luce perfectamente uniforme y aun así ser inservible porque su secuencia se repite cada pocos miles de valores — ese defecto solo se detecta con el **autocorrelograma**. Ese es justamente el punto que vamos a comprobar en este notebook.\n",
    "\n",
    "Vamos a construir y comparar 4 generadores, cada uno con $N = 100\\,000$ números:\n",
    "\n",
    "| Generador | Tipo | $a$ | $c$ | $m$ |\n",
    "|---|---|---:|---:|---:|\n",
    "| Multiplicativo bueno | multiplicativo | 16807 | 0 | $2^{31}-1$ |\n",
    "| Multiplicativo malo | multiplicativo | 1001 | 0 | $2^{14}$ |\n",
    "| Mixto bueno | mixto | 1664525 | 1013904223 | $2^{32}$ |\n",
    "| Mixto malo | mixto | 1001 | 1 | 2048 |\n",
    "\n",
    "El \"multiplicativo bueno\" corresponde al generador de Lehmer / Park–Miller (\"minimal standard\"), una referencia clásica en la literatura de simulación (p. ej. Banks et al.). El \"mixto bueno\" usa los parámetros sugeridos por *Numerical Recipes* (basados en Knuth). Los dos generadores \"malos\" usan módulos pequeños con parámetros que **no** garantizan periodo completo, para que puedas observar sus defectos de primera mano.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10a55f1d",
   "metadata": {},
   "source": [
    "## 1. Implementación del generador congruencial general \n",
    "\n",
    "Implementa una función genérica que sirva para los 4 generadores (multiplicativo y mixto son casos particulares, según `c`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd204d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generador_congruencial(semilla, a, c, m, n):\n",
    "    \"\"\"\n",
    "    Genera una secuencia de n números pseudoaleatorios U(0,1) usando el\n",
    "    generador congruencial lineal:\n",
    "        X_{i+1} = (a * X_i + c) mod m\n",
    "        U_i     = X_i / m\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    semilla : int -- estado inicial X0 (0 <= semilla < m)\n",
    "    a : int       -- multiplicador\n",
    "    c : int       -- incremento (c = 0 => método multiplicativo; c != 0 => método mixto)\n",
    "    m : int       -- módulo\n",
    "    n : int       -- cantidad de números a generar\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    u : np.ndarray, forma (n,)            -- secuencia en [0, 1)\n",
    "    x : np.ndarray, forma (n,), int64     -- estados enteros X_i (útiles para el periodo)\n",
    "    \"\"\"\n",
    "    # TODO 1: crea un arreglo `x` de enteros (dtype=np.int64) de tamaño n\n",
    "    # TODO 2: inicializa `estado = semilla`\n",
    "    # TODO 3: itera n veces: estado = (a * estado + c) % m; guarda el estado en x[i]\n",
    "    # TODO 4: calcula u = x / m\n",
    "    # TODO 5: retorna (u, x)\n",
    "    raise NotImplementedError(\"Implementa generador_congruencial\")\n",
    "\n",
    "# --- prueba rápida de tu implementación (no la modifiques) ---\n",
    "u_test, x_test = generador_congruencial(semilla=1, a=5, c=3, m=16, n=8)\n",
    "print(\"Estados: \", x_test)\n",
    "print(\"U(0,1):  \", np.round(u_test, 4))\n",
    "assert u_test.shape == (8,)\n",
    "assert np.all((u_test >= 0) & (u_test < 1))\n",
    "assert len(np.unique(x_test)) == 8  # con estos parámetros el periodo es 16 (completo)\n",
    "print(\"OK: la implementación pasa la prueba básica.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc3530a8",
   "metadata": {},
   "source": [
    "## 2. Los cuatro generadores a comparar \n",
    "\n",
    "Usa `generador_congruencial` para producir las 4 secuencias de tamaño `N` descritas en la tabla de la sección 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f938adf",
   "metadata": {},
   "outputs": [],
   "source": [
    "PARAMETROS = {\n",
    "    \"Multiplicativo bueno\": dict(semilla=1, a=16807,   c=0,          m=2**31 - 1),\n",
    "    \"Multiplicativo malo\":  dict(semilla=1, a=1001,    c=0,          m=2**14),\n",
    "    \"Mixto bueno\":          dict(semilla=1, a=1664525, c=1013904223, m=2**32),\n",
    "    \"Mixto malo\":           dict(semilla=1, a=1001,    c=1,          m=2048),\n",
    "}\n",
    "\n",
    "secuencias = {}\n",
    "# TODO: para cada (nombre, p) en PARAMETROS.items():\n",
    "#   - llama u, x = generador_congruencial(n=N, **p)\n",
    "#   - guarda en secuencias[nombre] un diccionario con u, x y los parámetros de p\n",
    "#     (pista: dict(u=u, x=x, **p))\n",
    "\n",
    "\n",
    "# TODO: imprime, para cada generador: a, c, m, min(u), max(u) y media(u)\n",
    "\n",
    "\n",
    "assert set(secuencias.keys()) == set(PARAMETROS.keys())\n",
    "assert all(len(d[\"u\"]) == N for d in secuencias.values())\n",
    "print(\"OK: las 4 secuencias fueron generadas.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76f2b4b",
   "metadata": {},
   "source": [
    "## 3. Detección empírica del periodo \n",
    "\n",
    "Antes de graficar, detecta el periodo de cada secuencia: el número de pasos que tarda el generador en repetir un estado ya visitado. Esto te va a ayudar a interpretar el autocorrelograma más adelante."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d34772f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detectar_periodo(x, tope=200_000):\n",
    "    \"\"\"\n",
    "    Detecta el periodo de una secuencia de estados enteros `x` buscando la\n",
    "    primera repetición de un estado ya visto (dentro de los primeros `tope`\n",
    "    estados).\n",
    "\n",
    "    Retorna el periodo (int) si se encontró una repetición, o None si no se\n",
    "    observó ninguna repetición dentro de `tope` estados (es decir, el\n",
    "    periodo real es >= tope).\n",
    "    \"\"\"\n",
    "    # TODO 1: crea un diccionario vacío `vistos` para guardar {valor: índice}\n",
    "    # TODO 2: recorre x[:tope] con su índice i\n",
    "    # TODO 3: si el valor ya está en `vistos`, retorna (i - vistos[valor])\n",
    "    # TODO 4: si no, guárdalo: vistos[valor] = i\n",
    "    # TODO 5: si terminas el recorrido sin repeticiones, retorna None\n",
    "    raise NotImplementedError(\"Implementa detectar_periodo\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert detectar_periodo(np.array([3, 7, 2, 3, 7, 2, 3])) == 3\n",
    "assert detectar_periodo(np.array([1, 2, 3, 4, 5])) is None\n",
    "print(\"OK: detectar_periodo pasa la prueba básica.\")\n",
    "\n",
    "periodos = {}\n",
    "# TODO: calcula periodos[nombre] = detectar_periodo(d[\"x\"]) para cada (nombre, d) en secuencias.items()\n",
    "\n",
    "\n",
    "# TODO: imprime, para cada generador, su módulo m y el periodo detectado\n",
    "#       (si periodos[nombre] es None, indica que el periodo es >= al tope usado)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67895c7b",
   "metadata": {},
   "source": [
    "## 4. Scatter \n",
    "\n",
    "Grafica, en una cuadrícula 2x2, el histograma de cada una de las 4 secuencias (normalizado como densidad) junto con la densidad teórica de $U(0,1)$, que es constante e igual a 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d79fcf60",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(11, 8))\n",
    "for ax, (nombre, d) in zip(axes.flat, secuencias.items()):\n",
    "    # TODO: dibuja el scatter de d[\"u\"] en `ax`\n",
    "    # \n",
    "    # TODO: pon como título `nombre` y etiqueta los ejes (x: \"observación\", y: \"u\")\n",
    "    pass\n",
    "fig.suptitle(f\"Scatter plots de las secuencias U(0,1)  (N = {N:,}, {BINS} bins)\", y=1.02)\n",
    "fig.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aa0cd07",
   "metadata": {},
   "source": [
    "**Para pensar (no hay que responder aquí todavía):** con `BINS = 50` bins, ¿los cuatro histogramas lucen razonablemente uniformes? Guarda esta observación — la vas a necesitar en la sección de preguntas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a26fe88f",
   "metadata": {},
   "source": [
    "## 5. Autocorrelograma \n",
    "\n",
    "Implementa el coeficiente de autocorrelación muestral y grafícalo para cada generador hasta `MAXLAG` rezagos:\n",
    "\n",
    "$$r_k = \\frac{\\displaystyle\\sum_{i=0}^{n-k-1} (u_i - \\bar u)(u_{i+k} - \\bar u)}{\\displaystyle\\sum_{i=0}^{n-1} (u_i - \\bar u)^2}, \\qquad k = 0, 1, \\dots, \\text{MAXLAG}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd841c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorrelacion(u, max_lag):\n",
    "    \"\"\"\n",
    "    Calcula el coeficiente de autocorrelación muestral r_k para\n",
    "    k = 0, 1, ..., max_lag. Retorna un arreglo `r` de longitud max_lag + 1\n",
    "    (r[0] siempre es 1.0).\n",
    "    \"\"\"\n",
    "    # TODO 1: centra la serie restando la media: centrado = u - u.mean()\n",
    "    # TODO 2: calcula la varianza total (denominador, es el mismo para todo k):\n",
    "    #         varianza_total = dot(centrado, centrado)\n",
    "    # TODO 3: para cada k en 0..max_lag calcula:\n",
    "    #         r[k] = dot(centrado[:n-k], centrado[k:]) / varianza_total\n",
    "    # TODO 4: retorna el arreglo r\n",
    "    raise NotImplementedError(\"Implementa autocorrelacion\")\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "r_test = autocorrelacion(np.array([1.0, 2.0, 3.0, 4.0, 5.0]), max_lag=2)\n",
    "assert np.isclose(r_test[0], 1.0)\n",
    "print(\"OK: autocorrelacion pasa la prueba básica. r_test =\", np.round(r_test, 4))\n",
    "\n",
    "autocorrelaciones = {}\n",
    "# TODO: calcula autocorrelaciones[nombre] = autocorrelacion(d[\"u\"], MAXLAG) para cada generador\n",
    "\n",
    "\n",
    "banda = 1.96 / np.sqrt(N)  # banda de referencia aproximada para ruido blanco\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(11, 8))\n",
    "for ax, (nombre, r) in zip(axes.flat, autocorrelaciones.items()):\n",
    "    # TODO: grafica r vs. el rezago k (lags = np.arange(len(r)))\n",
    "    # TODO: dibuja las bandas de referencia horizontales en +banda y -banda, y una línea en 0\n",
    "    # TODO: si periodos[nombre] no es None y es <= MAXLAG, dibuja una línea vertical\n",
    "    #       en x = periodos[nombre] para marcar el periodo detectado\n",
    "    # TODO: título = nombre, etiquetas de ejes (\"rezago k\", \"r_k\")\n",
    "    pass\n",
    "fig.suptitle(f\"Autocorrelograma  (N = {N:,}, rezago máximo = {MAXLAG})\", y=1.02)\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# TODO: imprime, para cada generador, el máximo |r_k| para k>=1 y en qué rezago k ocurre\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "467fdc3d",
   "metadata": {},
   "source": [
    "## 6. Prueba de bondad de ajuste chi-cuadrado (Nivel avanzado)\n",
    "\n",
    "Complementa el análisis visual con una prueba formal de uniformidad sobre los conteos por bin del histograma: $H_0$: los datos provienen de una $U(0,1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdd3482a",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "print(f\"{'Generador':22s} | {'chi2':>10s} | {'p-value':>10s} | decisión (alpha=0.05)\")\n",
    "for nombre, d in secuencias.items():\n",
    "    # TODO: calcula los conteos por bin con np.histogram (bins=BINS, range=(0,1))\n",
    "    # TODO: calcula el conteo esperado bajo H0 (uniforme): esperado = N / BINS\n",
    "    # TODO: usa stats.chisquare(conteos, f_exp=[esperado]*BINS) para obtener (chi2, p)\n",
    "    # TODO: define `decision` como \"no se rechaza H0 (uniforme)\" si p >= alpha,\n",
    "    #       o \"se rechaza H0\" en caso contrario\n",
    "    # TODO: imprime la fila: nombre, chi2, p, decision\n",
    "    pass\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}