En este Notebook revisaremos brevemente un par de resultados que nos ayudarán a generar números aleatorios de diversas distribuciones. Comenzamos definiendo qué es una sucesión de números aleatorios uniformemente distribuidos en el intervalo $[0,1]$.
<aside> 💡
Definición (Sucesión de números aleatorios)
Una sucesión $\{ x_n \}_{n\in \mathbb{N}^{+}}$ de números reales en el intervalo $[0,1]\subset \mathbb{R}$ es una sucesión de números aleatorios si existen
Un espacio de probabilidad $(\Omega, \mathcal{F}, P)$,
Una sucesión $\{ U_n \}_{n\in \mathbb{N}^{+}}$ de variables aleatorias independientes e idénticamente distribuidas con distribución $\mathbf{U}([0,1])$, y
$\omega \in \Omega$ tal que $\forall n \in \mathbb{N}^{+}$, $x_n = U_n(\omega)$.
</aside>
Intuitivamente, podemos pensar en números aleatorios uniformemente distribuidos como las realizaciones de variables aleatorias uniformes.
Más adelante veremos que una consecuencia del Principio fundamental de la simulación es que cualquier variable aleatoria puede ser simulada usando únicamente la distribución uniforme en $[0,1]$. Por lo anterior, el problema de simular variables aleatorias se reduce a saber cómo simular variables aleatorias uniformes (es decir, cómo generar números aleatorios uniformemente distribuidos en $[0,1]$).
Un procedimiento usual para generar números (pseudo)-aleatorios es generar números enteros $y_n$ a través de la relación congruencial de recurrencia,
$$ \begin{equation} y_{n+1} \equiv a y_n + b \hspace{0.5 cm} \text{mod } N \end{equation} $$
y luego construir la sucesión de números pseudo-aleatorios como sigue,
$$ x_{n} = \frac{y_n}{N}, \hspace{0.5 cm} y_n \in \{0, \dots, N-1 \} $$
donde $mcd(a, N) = 1$ (es decir, $a$ y $N$ son primos relativos).
Note que, por construcción, $x_n \in [0,1)$ para toda $n$, por lo que son buenos candidatos para ser números (pseudo)-aleatorios uniformemente distribuidos en $(0,1)$.
La parte fina del generador lineal congruencial es saber elegir valores de $a, b$ y $N$ adecuados. Para ver los detalles teóricos véase [1].
Por simplicidad, optaremos por hacer un generador lineal congruencial siguiendo algunas ideas del generador Mersenne Twister y del generador de Lehmer.
El generador Lehmer supone que $b=0$, y el generador Mersenne Twister que $N$ es un número primo de Mersenne. Además necesitamos que $mcd(y_0, N) = 1$ y que el orden multiplicativo módulo N de $a$ sea grande (aunque esto es importante para obtener el mayor periodo, no es obligatorio).
La relación congruencial de recurrencia se reduce a,
$$ \begin{equation} y_{n+1} \equiv a y_n \hspace{0.5 cm} \text{mod } N \end{equation} $$
con $N = 2 ^{p} - 1$ con $p$ un número primo. Al valor inicial $y_0$ se le conoce como semilla (seed, en inglés).
Implementación del generador lineal congruencial
from math import gcd
import matplotlib.pyplot as plt
import numpy as np
from typing import List
class RandomGenerator:
def __init__(self, seed: int, N: int = (2**43)-1, a: int = 321456,
b: int = 0):
"""
:param N: Módulo. Un número primo de Mersenne. Por defecto, p = 43, N = 2^{p} - 1
:param a: Magnitud.
:param b: Salto. Por defecto su valor es 0 (Lehmer generator).
:param seed: Semilla para generar el valor inicial.
"""
# Verificamos que a y X_0 sean primos relativos con N.
assert gcd(N, a) == 1, print(f"mcd({a},{N}) != 1")
assert gcd(N, seed) == 1, print(f"mcd({seed},{N}) != 1")
self.y_0: int = seed
self.__N: int = N
self.__a: int = a
self.__b: int = b
def generate(self, size: int) -> np.ndarray[float]:
""" Genera un arreglo con self.size números pseudo-aleatorios con
distribución uniforme.
:param size: Tamaño de la muestra.
"""
# Arreglo de tamaño self.size para guardar los números generados.
U: np.ndarray[float] = np.zeros(size)
# Primer elemento y = y_1 = a * y_0 + b mod N.
y: int = (self.__a * self.y_0 + self.__b) % self.__N
# Dividimos por N para asegurar que el resultado esté entre 0 y 1
# y guardamos en el arreglo.
U[0] = y / self.__N
# Construimos los demás elementos recursivamente
for i in range(1, size):
y: float = (self.__a * y + self.__b) % self.__N
U[i] = y / self.__N
return U
Probamos el generador