1 Introducción.

Este proyecto tiene como propósito documentar y organizar un conjunto de códigos diseñados para la generación de números pseudoaleatorios, así como sus aplicaciones dentro del contexto de la simulación estocástica. A través de distintas técnicas y algoritmos, se busca explorar el comportamiento de estos generadores, su estructura matemática y su utilidad práctica en la simulación de fenómenos aleatorios. Esta documentación está orientada tanto a registrar el desarrollo técnico del proyecto como a facilitar su comprensión, reproducción y análisis.

2 Generadores

2.1 Generador Congruencial Multiplicativo.

El Generador Congruencial Multiplicativo (GCM) es un tipo de generador de números pseudoaleatorios definido por la recurrencia:
\[x_{n+1} = (b * x_n ) \: mod \: m \]
donde:

  •  \(x_0\)  es la semilla inicial
  •  \(b\)  es el multiplicador
  •  \(m\)  es el módulo

Cada número pseudoaleatorio \(u_n\) se obtiene normalizando el valor generado:

\[ u_n = \frac{x_n}{m} \]

Se han implementado varios generadores con diferentes parámetros clásicos:

  • RANDU: \[ m = 2^{31}, \quad b = 65539, \quad x_0 = 1 \]

  • MINSTD (Minimal Standard Generator): \[ m = 2^{31}, \quad b = 16807, \quad x_0 = 1 \]

  • RANF: \[ m = 2^{48}, \quad b = 44485709377909, \quad x_0 = 1 \]

  • ZX81 (basado en el generador de la computadora ZX81): \[ m = 2^{16} + 1, \quad b = 75, \quad x_0 = 1 \]

class GCM:
  def __init__(self, modulo, multiplicador, semilla):
    self.m = modulo
    self.b = multiplicador
    self.x0 = semilla
  def muestra(self, n):
    L = []
    x = self.x0
    for i in range(n):
      x = (self.b * x) % self.m
      L.append(x / self.m)
    self.x0 = x
    return L

RANDU = GCM(2**31, 65539, 1)
MINSTD = GCM(2**31, 16807, 1)
RANF = GCM(2**48, 44485709377909, 1)
ZX81 = GCM((2**16)+1, 75, 1)
## RANDU (10 muestras):
## [3.051897510886192e-05, 0.00018310965970158577, 0.0008239871822297573, 0.003295936156064272, 0.012359732296317816, 0.04449496837332845, 0.15573221957311034, 0.533938602078706, 0.8020416363142431, 0.006802399177104235]
## MINSTD (10 muestras):
## [7.826369255781174e-06, 0.1315377880819142, 0.7556042927317321, 0.44134794222190976, 0.7348649236373603, 0.8747715731151402, 0.2858293461613357, 0.9338209335692227, 0.7284304979257286, 0.7313786377198994]
## RANF (10 muestras):
## [0.15804498821804103, 0.8251314258663776, 0.3368007872298229, 0.8651650404273745, 0.07383628059541891, 0.5474493739795072, 0.5753776055978399, 0.6525021829214701, 0.7365476028938538, 0.5057520856606494]
## ZX81 (10 muestras):
## [0.0011443917176556754, 0.08582937882417566, 0.43720341181317424, 0.7902558859880678, 0.2691914491050857, 0.18935868288142577, 0.20190121610693196, 0.14259120801989716, 0.6943406014922868, 0.07554511192150999]

2.2 Generador Congruencial Lineal (GCL)

El Generador Congruencial Lineal (GCL) es una extensión del generador multiplicativo, e incluye una constante aditiva. Se define por la siguiente recurrencia:

\[ x_{n+1} = (b \cdot x_n + a) \: \bmod \: m \]

donde:

  • \(x_0 \quad\) es la semilla inicial
  • \(b \quad\) es el multiplicador
  • \(a \quad\) es la constante aditiva
  • \(m \quad\) es el módulo

Cada número pseudoaleatorio se normaliza así:

\[ u_n = \frac{x_n}{m} \]

2.2.0.1 Parámetros usados en el ejemplo:

\[ m = 13, \quad b = 17, \quad a = 33, \quad x_0 = 77 \]

class GCL:
  def __init__(self, modulo, multiplicador, constante, semilla):
    self.m = modulo
    self.b = multiplicador
    self.a = constante
    self.x0 = semilla
  def muestra(self, n):
    L = []
    x = self.x0
    for i in range(n):
      x = (self.b * x + self.a) % self.m
      L.append(x / self.m)
    self.x0 = x
    return L

gen1 = GCL(13, 17, 33, 77)
## GCL (muestra):
## [0.23076923076923078, 0.46153846153846156, 0.38461538461538464, 0.07692307692307693, 0.8461538461538461]

2.3 Generador Congruencial Polinomial (GCP)

El Generador Congruencial Polinomial (GCP) es una generalización de los generadores congruenciales, en el que se utiliza un polinomio evaluado en la semilla para generar nuevos valores. Su fórmula se define como:

\[ x_{n+1} = \left( \sum_{i=0}^{k} b_i \cdot x_n^i \right) \bmod m \]

donde:

  • \(x_0 \quad\) es la semilla inicial
  • \(b_0, b_1, \dots, b_k \quad\) son los coeficientes del polinomio
  • \(m \quad\) es el módulo

Y el valor pseudoaleatorio normalizado se calcula como:

\[ u_n = \frac{x_n}{m} \]

2.3.0.1 Parámetros usados en el ejemplo:
  • Polinomio: \(P(x) = 5 - x^2 + x^3\)
  • Coeficientes: \(b = [5,\ 0,\ -1,\ 1]\)
  • Módulo: \(m = 17\)
  • Semilla: \(x_0 = 2\)
class GCP:
    def __init__(self, modulo, coeficientes, semilla):
        self.m = modulo 
        self.b = coeficientes
        self.x0 = semilla

    def muestra(self, n):
        L = []
        x = self.x0
        for _ in range(n):
            x = (sum(self.b[i] * (x ** i) for i in range(len(self.b)))) % self.m
            L.append(x / self.m)
            self.x0 = x
        return L

gen1 = GCP(17, [5, 0, -1, 1], 2)
## GCP (1 muestra):
## [0.5294117647058824]

2.4 Generador Congruencial Multiplicativo Combinado (GCMC)

El Generador Congruencial Multiplicativo Combinado (GCMC) utiliza dos generadores congruenciales multiplicativos independientes y combina sus salidas para mejorar la calidad estadística del generador. La fórmula combinada se define como:

\[ \begin{aligned} x_{n+1} &= (b_1 \cdot x_n) \bmod m_1 \\ y_{n+1} &= (b_2 \cdot y_n) \bmod m_2 \\ z_{n+1} &= (x_{n+1} - y_{n+1}) \bmod m_1 \end{aligned} \]

donde:

  • \(x_0, y_0\) son las semillas iniciales
  • \(b_1, b_2\) son los multiplicadores
  • \(m_1, m_2\) son los módulos

El valor pseudoaleatorio combinado se normaliza como:

\[ u_n = \frac{z_n}{m_1} \]

2.4.0.1 Parámetros usados en el ejemplo:

\[ m_1 = 123, \quad m_2 = 131, \quad b_1 = 45, \quad b_2 = 28, \quad x_0 = 21, \quad y_0 = 11 \]

class GCMC:
  def __init__(self, modulo1, modulo2, multiplicador1, multiplicador2, semilla1, semilla2):
    self.m1 = modulo1
    self.m2 = modulo2
    self.b1 = multiplicador1
    self.b2 = multiplicador2
    self.x0 = semilla1
    self.y0 = semilla2

  def muestra(self, n):
    L1 = []
    L2 = []
    L = []
    x = self.x0
    y = self.y0
    for i in range(n):
      x = (self.b1 * x) % self.m1
      y = (self.b2 * y) % self.m2
      z = (x - y) % self.m1
      L1.append(x / self.m1)
      L2.append(y / self.m2)
      L.append(z / self.m1)
    self.x0 = x
    self.y0 = y
    return L

gencomb1 = GCMC(123, 131, 45, 28, 21, 11)
## GCMC (5 muestras):
## [0.3089430894308943, 0.8455284552845529, 0.6097560975609756, 0.34959349593495936, 0.3983739837398374]

2.5 Generador Congruencial Combinado Multivariado (GCCM)

El Generador Congruencial Combinado Multivariado (GCCM) es una extensión combinada de generadores congruenciales multiplicativos que considera múltiples valores pasados (tipo ventana de memoria) para mejorar la calidad del número generado. Está definido por:

\[ \begin{aligned} x_{n+1} &= \left( \sum_{i=0}^{k-1} b_{1,i} \cdot x_{n-i} \right) \bmod m_1 \\ y_{n+1} &= \left( \sum_{i=0}^{k-1} b_{2,i} \cdot y_{n-i} \right) \bmod m_2 \\ z_{n+1} &= (x_{n+1} - y_{n+1}) \bmod m_1 \end{aligned} \]

donde:

  • \(\mathbf{x}_0, \dots, \mathbf{x}_{k-1} \quad\) son los valores iniciales (semillas) del generador X
  • \(\mathbf{y}_0, \dots, \mathbf{y}_{k-1} \quad\) son los valores iniciales (semillas) del generador Y
  • \(b_{1,i}, b_{2,i} \quad\) son los coeficientes multiplicadores de cada generador
  • \(m_1, m_2 \quad\) son los módulos

Y la salida pseudoaleatoria es:

\[ u_n = \frac{z_n}{m_1} \]

2.5.0.1 Parámetros usados en el ejemplo:
  • \(m_1 = 13,\quad m_2 = 94\)
  • \(\mathbf{b_1} = [56, 21],\quad \mathbf{b_2} = [23, 43]\)
  • \(\mathbf{x_0} = [4, 45],\quad \mathbf{y_0} = [65, 90]\)
class GCCM:
    def __init__(self, modulo1, modulo2, multiplicador1, multiplicador2, semillasx, semillasy):
        self.m1 = modulo1
        self.m2 = modulo2
        self.b1 = multiplicador1
        self.b2 = multiplicador2
        self.x0 = semillasx
        self.y0 = semillasy

    def muestra(self, n):
        listaz = []
        k = len(self.b1)
        x = self.x0
        y = self.y0

        for i in range(n):
            nuevo_x = sum(self.b1[i] * x[i] for i in range(k)) % self.m1
            nuevo_y = sum(self.b2[i] * y[i] for i in range(k)) % self.m2
            z = (nuevo_x - nuevo_y) % self.m1
            x.insert(0, nuevo_x)
            self.y0.insert(0, nuevo_y)
            x.pop()
            y.pop()
            listaz.append(z / self.m1)
        return listaz

gencombm1 = GCCM(13, 94, [56, 21], [23, 43], [4, 45], [65, 90])
## GCCM (5 muestras):
## [0.38461538461538464, 0.9230769230769231, 0.5384615384615384, 0.6153846153846154, 0.6923076923076923]

2.6 Generador Congruencial Lineal Multivariado (GCLM)

El Generador Congruencial Lineal Multivariado (GCLM) utiliza una combinación lineal de varios estados anteriores para generar el siguiente valor de la secuencia. Su fórmula general es:

\[ x_{n+1} = \left( \sum_{i=0}^{k-1} b_i \cdot x_{n-i} \right) \bmod m \]

donde:

  • \(x_0, x_1, \dots, x_{k-1} \quad\) son los valores iniciales del generador (semillas)
  • \(b_0, b_1, \dots, b_{k-1} \quad\) son los multiplicadores
  • \(m \quad\) es el módulo

El valor generado se normaliza en el intervalo \((0, 1)\) como:

\[ u_n = \frac{x_n}{m} \]

2.6.0.1 Parámetros usados en el ejemplo:
  • \(\mathbf{x_0} = [7, 11, 26]\)
  • \(\mathbf{b} = [3, 5, 2]\)
  • \(m = 69\)
class GCLM:
    def __init__(self, semilla, multiplicadores, modulo):
        self.x0 = semilla
        self.b = multiplicadores
        self.m = modulo
    def muestra(self, n):
        resultados = []
        estado = self.x0
        k = len(estado)

        for _ in range(n):
            x = sum(self.b[i] * self.x0[i] for i in range(k)) % self.m
            estado.insert(0, x)
            estado.pop()
            resultados.append(x / self.m)
        return resultados

genlinm1 = GCLM([7, 11, 26], [3, 5, 2], 69)
## GCLM (5 muestras):
## [0.855072463768116, 0.391304347826087, 0.6521739130434783, 0.6231884057971014, 0.9130434782608695]

2.7 Funciones auxiliares para generadores binarios.

Estas funciones permiten representar números enteros como listas de bits (0 y 1), manipular matrices binarias y convertir resultados de vuelta a enteros.

# Convierte un número entero a su representación binaria como lista de bits (LSB primero)
def repbin(n):
    L = []
    while n > 0:
        r = n % 2
        n = n // 2
        L.append(r)
    return L

# Convierte un número entero a binario con longitud fija, rellenando con ceros a la izquierda si es necesario
def repbinlen(n, length):
    L = []
    while n > 0:
        r = n % 2
        n = n // 2
        L.append(r)
    while len(L) < length:
        L.append(0)
    return L

# Convierte una lista de bits a su valor entero correspondiente
def get_num(list_bin):
    n = 0
    for i in range(len(list_bin)):
        if list_bin[i] == 1:
            n += 2 ** i
    return n

# Multiplica una matriz binaria por un vector binario (operaciones XOR y AND bit a bit)
def multiplicar(matriz, vector):
    n = len(matriz)
    resultado = [0 for _ in range(n)]
    for i in range(n):  # Iterar sobre filas
        for j in range(n):  # Iterar sobre columnas
            resultado[i] ^= matriz[i][j] & vector[j]
    return resultado

2.8 Generador LFSR de tamaño 4 (Linear Feedback Shift Register)

El Generador de Retroalimentación Lineal (LFSR) es un generador pseudoaleatorio basado en operaciones matriciales binarias, especialmente útil en criptografía y simulaciones donde se requiere eficiencia. Su estructura se basa en la siguiente relación matricial:

\[ x_{n+1} = A \cdot x_n \]

donde:

  • \(x_n\) es un vector binario de longitud 4 que representa el estado actual.
  • \(A \in \{0,1\}^{4 \times 4}\) es la matriz de retroalimentación (shift + XOR).
  • La salida \(u_n\) se normaliza como:
    \[ u_n = \frac{x_n}{2^4} = \frac{x_n}{16} \]
2.8.0.1 Parámetros usados:
  • Semilla inicial: \(x_0 = 78\)
  • Matriz de transición:

\[ A = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{bmatrix} \]

class linear_feedback_generator_size_4:
    def __init__(self, semilla):
        self.matriz = [[False, True, False, False],
                       [False, False, True, False],
                       [False, False, False, True],
                       [True, False, False, True]]
        self.x0 = semilla
    def muestra(self, n):
        lista = []
        semilla = self.x0
        for i in range(n):
            nueva_semilla = get_num(multiplicar(self.matriz, repbinlen(semilla, 4)))
            lista.append(nueva_semilla / 16)
            semilla = nueva_semilla
        self.x0 = nueva_semilla
        return lista

genlf4 = linear_feedback_generator_size_4(78)
## LFSR (muestras):
## [0.9375, 0.4375, 0.6875, 0.3125, 0.625]

2.9 Generador LFSR de tamaño 8 (Linear Feedback Shift Register)

El Generador de Retroalimentación Lineal (LFSR) de tamaño 8 utiliza una matriz binaria de transición para desplazar los bits de estado y aplicar una combinación lineal (XOR) que determina el nuevo valor. Este tipo de generador es eficiente y produce secuencias pseudoaleatorias con propiedades deseables cuando se diseñan correctamente.

El generador se basa en la siguiente operación matricial:

\[ x_{n+1} = A \cdot x_n \]

donde:

  • \(x_n \quad\) es un vector binario de longitud 8 que representa el estado actual
  • \(A \quad\) es una matriz de retroalimentación de tamaño \(8 \times 8\)
  • El valor decimal generado se normaliza al intervalo \([0,1)\):

\[ u_n = \frac{x_n}{256} \]

2.9.0.1 Parámetros usados:
  • \(x_0 :\) Semilla inicial
  • Matriz de retroalimentación:

\[ A = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 \end{bmatrix} \]

class LinearFeedbackGeneratorSize8:
    def __init__(self, semilla):
        self.matriz = [
            [0, 1, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0, 0],
            [0, 0, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, 0, 1, 0, 0, 0],
            [0, 0, 0, 0, 0, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 0, 0, 1],
            [1, 1, 0, 0, 1, 1, 1, 1]
        ]
        self.x0 = semilla

    def muestra(self, n):
        lista = []
        semilla = self.x0
        for _ in range(n):
            nueva_semilla = get_num(multiplicar(self.matriz, repbinlen(semilla, 8)))
            lista.append(nueva_semilla / 256)
            semilla = nueva_semilla
        self.x0 = semilla
        return lista

genlf8 = LinearFeedbackGeneratorSize8(365)
## LFSR-8 (10 muestras):
## [0.7109375, 0.35546875, 0.17578125, 0.0859375, 0.04296875, 0.01953125, 0.5078125, 0.25390625, 0.125, 0.5625]

2.10 Generador LFSR de tamaño 12 (Linear Feedback Shift Register)

El Generador de Retroalimentación Lineal (LFSR) de tamaño 12 utiliza una matriz binaria de transición que realiza desplazamientos de bits en un vector de estado y aplica retroalimentación mediante una combinación lineal (XOR) en la última fila.

La evolución del sistema se describe mediante la operación:

\[ x_{n+1} = A \cdot x_n \]

donde:

  • \(x_n \quad\) es un vector binario de longitud 12 que representa el estado actual
  • \(A \quad\) es una matriz de retroalimentación de tamaño \(12 \times 12\)
  • El valor generado se convierte a decimal y se normaliza al intervalo \([0, 1)\):

\[ u_n = \frac{x_n}{4096} \]

2.10.0.1 Parámetros usados:
  • Semilla inicial: \(x_0 = 365\)
  • Matriz de retroalimentación:

\[ A = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \]

class LinearFeedbackGeneratorSize12:
    def __init__(self, semilla):
        self.matriz = [
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0]
        ]
        self.x0 = semilla

    def muestra(self, n):
        lista = []
        semilla = self.x0
        for _ in range(n):
            nueva_semilla = get_num(multiplicar(self.matriz, repbinlen(semilla, 12)))
            lista.append(nueva_semilla / 4096)
            semilla = nueva_semilla
        self.x0 = semilla
        return lista

genlf12 = LinearFeedbackGeneratorSize12(365)
## LFSR-12 (10 muestras):
## [0.04443359375, 0.022216796875, 0.010986328125, 0.50537109375, 0.252685546875, 0.126220703125, 0.56298828125, 0.781494140625, 0.890625, 0.9453125]

2.11 Generador LFSR de tamaño 16 (Linear Feedback Shift Register)

El Generador de Retroalimentación Lineal (LFSR) de tamaño 16 opera con una matriz de transición binaria que aplica desplazamientos de bits y una retroalimentación lineal basada en XOR. Este tipo de generador es eficiente para producir secuencias pseudoaleatorias con buena dispersión cuando se elige adecuadamente la matriz.

El proceso de generación sigue la relación:

\[ x_{n+1} = A \cdot x_n \]

donde:

  • \(x_n \quad\) es un vector binario de longitud 16 que representa el estado actual
  • \(A \quad\) es una matriz de retroalimentación de tamaño \(16 \times 16\)
  • El valor generado se normaliza como:

\[ u_n = \frac{x_n}{65536} \]

2.11.0.1 Parámetros usados:
  • Semilla inicial: \(x_0 = 365\)
  • Matriz de retroalimentación:

\[ A = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & 0 & 0 & \cdots & 1 \end{bmatrix} \] Debido a su tamaño, la matriz se describe parcialmente. Es una matriz de desplazamiento (superdiagonal) con la siguiente fila de retroalimentación:

Última fila: \([1,\ 0,\ 0,\ \dots,\ 1,\ 0,\ 1,\ 1,\ 0]\)

class LinearFeedbackGeneratorSize16:
    def __init__(self, semilla):
        self.matriz = [
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0]
        ]
        self.x0 = semilla

    def muestra(self, n):
        lista = []
        semilla = self.x0
        for _ in range(n):
            nueva_semilla = get_num(multiplicar(self.matriz, repbinlen(semilla, 16)))
            lista.append(nueva_semilla / 65536)
            semilla = nueva_semilla
        self.x0 = semilla
        return lista

genlf16 = LinearFeedbackGeneratorSize16(365)
## LFSR-16 (10 muestras):
## [0.502777099609375, 0.2513885498046875, 0.1256866455078125, 0.062835693359375, 0.0314178466796875, 0.0157012939453125, 0.507843017578125, 0.2539215087890625, 0.126953125, 0.5634765625]

2.12 Generador LFSR de tamaño 20 (Linear Feedback Shift Register)

El Generador de Retroalimentación Lineal (LFSR) de tamaño 20 aplica una operación de desplazamiento y retroalimentación XOR definida por una matriz binaria. Su estructura permite generar secuencias pseudoaleatorias con una longitud máxima cercana a \(2^{20} - 1\), dependiendo de la matriz usada.

Este generador se basa en la operación:

\[ x_{n+1} = A \cdot x_n \]

donde:

  • \(x_n \quad\) es un vector binario de 20 bits representando el estado actual
  • \(A \quad\) es la matriz de retroalimentación de tamaño \(20 \times 20\)
  • El resultado es convertido a decimal y normalizado como:

\[ u_n = \frac{x_n}{1048576} \]

2.12.0.1 Parámetros usados:
  • Semilla inicial: \(x_0 = 365\)
  • Matriz de retroalimentación:

Debido a su tamaño, la matriz se describe parcialmente. Es una matriz de desplazamiento (superdiagonal) con la siguiente fila de retroalimentación:

Última fila:
\([1,\ 0,\ 0,\ 1,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0,\ 0]\)

class LinearFeedbackGeneratorSize20:
    def __init__(self, semilla):
        self.matriz = [
            [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
            [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
            [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
        ]
        self.x0 = semilla

    def muestra(self, n):
        lista = []
        semilla = self.x0
        for _ in range(n):
            nueva_semilla = get_num(multiplicar(self.matriz, repbinlen(semilla, 20)))
            lista.append(nueva_semilla / 1048576)
            semilla = nueva_semilla
        self.x0 = semilla
        return lista

genlf20 = LinearFeedbackGeneratorSize20(365)
## LFSR-20 (10 muestras):
## [0.0001735687255859375, 8.678436279296875e-05, 4.291534423828125e-05, 2.09808349609375e-05, 1.049041748046875e-05, 4.76837158203125e-06, 0.5000019073486328, 0.2500009536743164, 0.625, 0.3125]

2.13 Generador MRG32k3a

El generador MRG32k3a es un generador múltiple recursivo combinado (Combined Multiple Recursive Generator) propuesto por Pierre L’Ecuyer. Es ampliamente utilizado por su largo período y buenas propiedades estadísticas.

El MRG32k3a combina dos generadores recursivos de orden 3 mediante la siguiente fórmula:

Sea:

  • \(m_1 = 4294967087\)
  • \(m_2 = 4294944443\)
  • \(a_{12} = 1403580,\quad a_{13} = 810728\)
  • \(a_{21} = 527612,\quad a_{23} = 1370589\)

La recursión para cada secuencia es:

\[ s_{n}^{(1)} = (a_{12} s_{n-2}^{(1)} - a_{13} s_{n-3}^{(1)}) \mod m_1 \]

\[ s_{n}^{(2)} = (a_{21} s_{n-1}^{(2)} - a_{23} s_{n-3}^{(2)}) \mod m_2 \]

Y la combinación final es:

\[ z_n = (s_{n}^{(1)} - s_{n}^{(2)}) \mod m_1 \]

El valor pseudoaleatorio generado se normaliza como:

\[ u_n = \begin{cases} \frac{z_n}{m_1} & \text{si } z_n > 0 \\ \frac{m_1 - 1}{m_1} & \text{si } z_n = 0 \end{cases} \]

2.13.0.1 Estado inicial por defecto:
  • \(x_0 = 1, x_1 = 1, x_2 = 20250409\)
  • \(y_0 = 1, y_1 = 1, y_2 = 1\)
class MRG32k3a:
    def __init__(self, x0=1, x1=1, x2=20250409, y0=1, y1=1, y2=1):
        self.m1 = 4294967087
        self.m2 = 4294944443
        self.a12 = 1403580
        self.a13n = 810728
        self.a21 = 527612
        self.a23n = 1370589
        self.s1 = [x2 % self.m1, x1 % self.m1, x0 % self.m1]
        self.s2 = [y2 % self.m2, y1 % self.m2, y0 % self.m2]

    def muestra(self, n=1):
        resultados = []
        for _ in range(n):
            p1 = (self.a12 * self.s1[1] - self.a13n * self.s1[2]) % self.m1
            self.s1 = [p1] + self.s1[:2]

            p2 = (self.a21 * self.s2[0] - self.a23n * self.s2[2]) % self.m2
            self.s2 = [p2] + self.s2[:2]

            z = (self.s1[0] - self.s2[0]) % self.m1
            resultados.append(z / self.m1 if z > 0 else (self.m1 - 1) / self.m1)
        return resultados
## [0.0003395772238661628, 0.31734079130092296, 0.49998174223494346, 0.22045602301957756, 0.4334982178642246]

3 Generadores de Variables con algunas distribuciones.

3.1 Generador de variables Bernoulli

El Generador de Bernoulli produce variables aleatorias que toman el valor \(1\) con probabilidad \(p\) (éxito), y el valor \(0\) con probabilidad \(1 - p\) (fracaso). Es una de las distribuciones discretas más simples y se basa en la comparación de una variable aleatoria uniforme con el parámetro \(p\):

\[ X \sim \text{Bernoulli}(p) \]

\[ P(X = 1) = p, \quad P(X = 0) = 1 - p \]

La implementación se realiza generando una secuencia de valores uniformemente distribuidos \(U \sim \mathcal{U}(0,1)\) y aplicando:

\[ X = \begin{cases} 1 & \text{si } U < p \\ 0 & \text{si } U \geq p \end{cases} \]

3.1.0.1 Parámetro usado:
  • \(p = 0.25\)
class GeneradorBernoulli:
    def __init__(self, probabilidad):
        self.p = probabilidad

    def muestra(self, n):
        uniformes = generador.muestra(n)
        return [1 if u < self.p else 0 for u in uniformes]
muestras_bernoulli = GeneradorBernoulli(0.80).muestra(5)
## Bernoulli (5 muestras, p = 0.25):
## [1, 1, 0, 1, 1]

3.2 Generador de variables Binomiales

La distribución binomial representa el número de éxitos en \(n\) ensayos independientes, cada uno con probabilidad de éxito \(p\). Es una generalización de la distribución Bernoulli:

\[ X \sim \text{Binomial}(n, p) \]

\[ P(X = k) = \binom{n}{k} p^k (1 - p)^{n - k}, \quad \text{para } k = 0, 1, \dots, n \]

Para simular esta distribución, se generan \(n\) variables independientes \(U_i \sim \mathcal{U}(0,1)\) y se cuentan cuántas cumplen \(U_i < p\):

\[ X = \sum_{i=1}^{n} \mathbb{I}_{\{U_i < p\}} \]

3.2.0.1 Parámetros usados:
  • \(p = 0.25\) — probabilidad de éxito en cada ensayo
  • \(n = 10\) — número de ensayos
class GeneradorBinomial:
    def __init__(self, p, n):
        self.p = p
        self.n = n

    def muestra(self, m):
        resultados = []
        for _ in range(m):
            uniformes = generador.muestra(self.n)
            ensayos = [1 if u < self.p else 0 for u in uniformes]
            resultados.append(sum(ensayos))
        return resultados
muestras_binomial = GeneradorBinomial(p=0.25, n=10).muestra(5)
## Binomial (5 muestras, n=10, p=0.25):
## [1, 2, 2, 1, 1]

3.3 Generador de variables Geométricas

La distribución geométrica modela la cantidad de ensayos necesarios hasta obtener el primer éxito, en una secuencia de ensayos de Bernoulli independientes con probabilidad de éxito \(p\). Su función de probabilidad es:

\[ X \sim \text{Geom}(p) \]

\[ P(X = k) = (1 - p)^{k - 1} \cdot p, \quad \text{para } k = 1, 2, 3, \dots \]

Para simular esta distribución, se generan variables uniformes hasta obtener el primer valor menor a \(p\), y se cuenta cuántos intentos fueron necesarios:

3.3.0.1 Parámetro usado:
  • \(p = 0.1\)
class GeneradorGeometrico:
    def __init__(self, p):
        self.p = p

    def muestra(self, n):
        resultados = []
        for _ in range(n):
            conteo = 0
            while True:
                conteo += 1
                u = generador.muestra(1)[0]
                if u < self.p:
                    resultados.append(conteo)
                    break
        return resultados
muestras_geometricas = GeneradorGeometrico(p=0.1).muestra(5)
## Geométrica (5 muestras, p=0.1):
## [32, 5, 6, 6, 9]

3.4 Generador de variables Binomiales Negativas

La distribución binomial negativa modela el número de ensayos necesarios hasta obtener \(r\) éxitos, con una probabilidad de éxito \(p\) en cada ensayo independiente. Su función de probabilidad es:

\[ X \sim \text{NegBin}(r, p) \]

\[ P(X = k) = \binom{k-1}{r-1} p^r (1 - p)^{k - r}, \quad \text{para } k = r, r+1, \dots \]

Para simularla, se generan ensayos hasta acumular \(r\) éxitos, y se registra el número total de intentos realizados.

3.4.0.1 Parámetros usados:
  • \(p = 0.1\) — probabilidad de éxito
  • \(r = 3\) — número de éxitos requeridos
class GeneradorBinomialNegativo:
    def __init__(self, p, r):
        self.p = p
        self.r = r

    def muestra(self, n):
        resultados = []
        for _ in range(n):
            ensayos = 0
            exitos = 0
            while exitos < self.r:
                ensayos += 1
                u = generador.muestra(1)[0]
                if u < self.p:
                    exitos += 1
            resultados.append(ensayos)
        return resultados
muestras_neg_bin = GeneradorBinomialNegativo(p=0.1, r=3).muestra(5)
## Binomial Negativa (5 muestras, p=0.1, r=3):
## [42, 52, 72, 20, 38]

3.5 Generador de variables Uniformes Discretas

La distribución uniforme discreta genera valores enteros dentro de un intervalo \([a, b]\), todos con igual probabilidad. Su función de probabilidad es:

\[ X \sim \text{Unif}(a, b) \]

\[ P(X = x) = \frac{1}{b - a + 1}, \quad \text{para } x = a, a+1, \dots, b \]

Se simula generando un número uniforme \(u \in (0, 1)\) y aplicando la transformación:

\[ X = \lfloor u \cdot (b - a + 1) \rfloor + a \]

3.5.0.1 Parámetros usados:
  • \(a = 0\)
  • \(b = 100\)
class GeneradorUniformeDiscreto:
    def __init__(self, a, b):
        self.a = a
        self.b = b

    def muestra(self, n):
        resultados = []
        rango = self.b - self.a + 1
        for _ in range(n):
            u = generador.muestra(1)[0]
            valor = int(u * rango) + self.a
            resultados.append(valor)
        return resultados
muestras_uniforme_discreto = GeneradorUniformeDiscreto(a=0, b=100).muestra(5)
## Uniforme Discreta (5 muestras, a=0, b=100):
## [84, 98, 28, 69, 66]

3.6 Generador de variables Exponenciales

La distribución exponencial modela el tiempo entre eventos en un proceso de Poisson. Su función de densidad es:

\[ X \sim \text{Exp}(\lambda) \]

\[ f(x) = \lambda e^{-\lambda x}, \quad x \geq 0 \]

Para simularla, se usa el método de la inversa:

\[ X = -\frac{1}{\lambda} \ln(1 - U) \]

donde \(U\) es una variable aleatoria uniforme en \((0,1)\).

3.6.0.1 Parámetro usado:
  • \(\lambda = 0.5\)
import math
class GeneradorExponencial:
    def __init__(self, lamb):
        self.lamb = lamb

    def simulacion(self):
        u = generador.muestra(1)[0]
        return -math.log(1 - u) / self.lamb

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
muestras_exponenciales = GeneradorExponencial(lamb=0.5).muestra(5)
## Exponencial (5 muestras, lambda = 0.5):
## [1.3039628125141254, 1.0660299057264315, 5.799319886774467, 0.9132174873842833, 1.6542525069840919]

3.7 Generador de variables Uniformes Continuas

La distribución uniforme continua genera valores dentro del intervalo \([a, b]\), con densidad constante. Su función de densidad es:

\[ X \sim \text{Unif}(a, b) \]

\[ f(x) = \frac{1}{b - a}, \quad \text{para } a \leq x \leq b \]

Para simularla, se utiliza la transformación directa de una variable uniforme \(U \in (0,1)\):

\[ X = a + (b - a) \cdot U \]

3.7.0.1 Parámetros usados:
  • \(a = 0\)
  • \(b = 1\)
class GeneradorUniformeContinua:
    def __init__(self, a, b):
        self.min = min(a, b)
        self.max = max(a, b)

    def simulacion(self):
        u = generador.muestra(1)[0]
        return self.min + (self.max - self.min) * u

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
muestras_uniforme_continua = GeneradorUniformeContinua(a=0, b=1).muestra(5)
## Uniforme Continua (5 muestras, a = 0, b = 1):
## [0.16519537114674052, 0.576132003313766, 0.49449745480668916, 0.6066155216616215, 0.894411314495831]

3.8 Generador de variables Poisson

La distribución de Poisson modela el número de eventos en un intervalo fijo de tiempo o espacio, cuando los eventos ocurren con una tasa constante \(\lambda\) y de forma independiente entre sí. Se denota:

\[ X \sim \text{Poisson}(\lambda) \]

Su función de probabilidad es:

\[ P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \dots \]

Para simular esta distribución, se utiliza un método basado en tiempos entre llegadas (suma de variables exponenciales):

\[ \text{Generar } t_i \sim \text{Exp}(\lambda), \quad \text{contar cuántos } t_i \text{ caben en } [0,1] \]

3.8.0.1 Parámetro usado:
  • \(\lambda = 15\)
import math
class GeneradorPoisson:
    def __init__(self, lamb):
        self.lamb = lamb

    def simulacion(self):
        n = 0
        t = 0.0
        while True:
            u = generador.muestra(1)[0]
            delta = -math.log(1 - u) / self.lamb
            if t + delta > 1.0:
                break
            t += delta
            n += 1
        return n

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
muestras_poisson = GeneradorPoisson(lamb=15).muestra(5)
## Poisson (5 muestras, lambda = 15):
## [22, 15, 8, 15, 8]

3.9 Generador de variables Poisson (Método Alternativo)

Esta implementación alternativa para simular variables aleatorias con distribución Poisson se basa en la suma de tiempos entre llegadas generados por una distribución exponencial:

\[ t_i \sim \text{Exp}(\lambda) \]

El número de eventos \(N\) que ocurren en el intervalo \([0, 1]\) se calcula sumando iterativamente los \(t_i\) hasta que la suma exceda 1:

\[ N = \max\left\{n \in \mathbb{N} \mid \sum_{i=1}^{n} t_i \leq 1 \right\} \]

3.9.0.1 Parámetro usado:
  • \(\lambda = 15\)
class GeneradorPoissonAlternativo:
    def __init__(self, parametro_lamb):
        self.parametro_lamb = parametro_lamb
        self.exponencial = GeneradorExponencial(parametro_lamb)

    def simulacion(self):
        n = 0
        t = self.exponencial.simulacion()
        while t <= 1.0:
            n += 1
            t += self.exponencial.simulacion()
        return n

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
muestras_poisson = GeneradorPoissonAlternativo(parametro_lamb=15).muestra(5)
## Poisson Alternativo (5 muestras, lambda = 15):
## [17, 11, 19, 14, 13]

3.10 Generador de variables distribución Normal

Este generador produce variables aleatorias con distribución normal \(\mathcal{N}(\mu, \sigma^2)\) usando el método polar de Box-Muller, el cual transforma dos variables uniformes independientes \(U_1, U_2 \sim \text{Unif}(0,1)\) en una variable normal estándar:

  1. Se generan dos valores: \[ V_1 = 2U_1 - 1, \quad V_2 = 2U_2 - 1 \]

  2. Se calcula: \[ R = V_1^2 + V_2^2 \]

  3. Si \(0 < R \leq 1\), se acepta y se obtiene: \[ Z = V_1 \cdot \sqrt{\frac{-2 \ln R}{R}} \]

  4. Finalmente, la variable normal con media y varianza específicas es: \[ X = \mu + \sigma Z \]

3.10.0.1 Parámetros usados:
  • \(\mu = 0\)
  • \(\sigma^2 = 1\)
class GeneradorNormal:
    def __init__(self, media, varianza):
        self.media = media
        self.varianza = varianza

    def simulacion(self):
        while True:
            u1 = generador.muestra(1)[0]
            u2 = generador.muestra(1)[0]
            v1 = 2 * u1 - 1
            v2 = 2 * u2 - 1
            r = v1 * v1 + v2 * v2
            if 0 < r <= 1:
                break
        factor = math.sqrt(-2 * math.log(r) / r)
        z1 = v1 * factor
        return self.media + math.sqrt(self.varianza) * z1

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
muestras_normales = GeneradorNormal(media=0, varianza=1).muestra(10)
## Normal (10 muestras, media = 0, varianza = 1):
## [0.20466767400891167, -0.5555880474927649, 0.42941046677492706, -1.0927184609735743, -0.2131141365416757, -0.4833419316066803, 0.1665598132548493, -1.105867018213305, 1.6056721581789444, -0.9083520442161196]

3.11 Generador de variables normales rectificadas

Este generador produce variables aleatorias con distribución normal rectificada \(\mathcal{N}^{+}(\mu, \sigma^2)\), es decir, una distribución normal donde todos los valores negativos son truncados a cero. Se basa en el método anterior para simular una variable \(X \sim \mathcal{N}(\mu, \sigma^2)\), pero se impone la condición:

\[ X = \begin{cases} \mu + \sigma Z & \text{si } \mu + \sigma Z > 0 \\ 0 & \text{en otro caso} \end{cases} \]

3.11.0.1 Parámetros usados:
  • \(\mu = 0\)
  • \(\sigma^2 = 1\)
class GeneradorNormalRectificado:
    def __init__(self, media, varianza):
        self.media = media
        self.varianza = varianza

    def simulacion(self):
        while True:
            u1 = generador.muestra(1)[0]
            u2 = generador.muestra(1)[0]
            v1 = 2 * u1 - 1
            v2 = 2 * u2 - 1
            r = v1*v1 + v2*v2
            if 0 < r <= 1:
                break
        factor = math.sqrt(-2 * math.log(r) / r)
        z = v1 * factor
        x = self.media + math.sqrt(self.varianza) * z
        return x if x > 0 else 0.0

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]

muestras_nrect = GeneradorNormalRectificado(media=0, varianza=1).muestra(10)
## Muestras Normal Rectificada N(0,1):
## [1.0880483983117055, 0.5355938290489028, 0.0, 0.0, 1.0333716691083992, 0.0, 0.0, 0.0, 1.4488558084152579, 1.9141408081884506]

3.12 Generador de variables normales acotadas

Este generador simula variables aleatorias con distribución normal truncada, es decir, variables \(X \sim \mathcal{N}(\mu, \sigma^2)\) que se restringen al intervalo \([a, b]\). La generación se realiza con el método anterior y se aplica rechazo si el valor cae fuera de los límites.

3.12.0.1 Fórmula base:

\[ X = \mu + \sigma Z, \quad \text{con } Z \sim \mathcal{N}(0, 1) \] Se acepta \(X\) solo si: \[ a \leq X \leq b \]

3.12.0.2 Parámetros usados:
  • \(\mu = 0\)
  • \(\sigma^2 = 1\)
  • \(a = -2\)
  • \(b = 2\)
class GeneradorNormalAcotada:
    def __init__(self, media, varianza, a, b):
        self.media = media
        self.varianza = varianza
        self.a = a
        self.b = b

    def simulacion(self):
        while True:
            u1 = generador.muestra(1)[0]
            u2 = generador.muestra(1)[0]
            v1 = 2 * u1 - 1
            v2 = 2 * u2 - 1
            r = v1*v1 + v2*v2
            if 0 < r <= 1:
                break
        factor = math.sqrt(-2 * math.log(r) / r)
        z = v1 * factor
        x = self.media + math.sqrt(self.varianza) * z
        if self.a <= x <= self.b:
            return x
        return self.simulacion()

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]

a = -2
b = 2
muestras_norm_acot = GeneradorNormalAcotada(media=0, varianza=1, a=a, b=b).muestra(10)
## Muestras Normal Acotada en [-2,2]:
## [-0.12381913653855846, 1.6104847902959503, -1.4359182477062387, 1.4956869166023063, 0.6696148747804004, 0.3400698127312014, 0.11977166074151453, 0.9345804605176253, 0.3217251398915182, -0.4933379139896068]

3.13 Generador de variables Chi-cuadrada

Este generador simula variables aleatorias con distribución \(\chi^2_k\), que representa la suma de los cuadrados de \(k\) variables normales estándar independientes: \[ X = \sum_{i=1}^k Z_i^2, \quad Z_i \sim \mathcal{N}(0,1) \]

3.13.0.1 Parámetro usado:
  • \(k = 5\) grados de libertad
class GeneradorChiCuadrada:
    def __init__(self, grados_libertad):
        self.n = grados_libertad
        self.normal = GeneradorNormal(0,1)

    def simulacion(self):
        suma = 0.0
        for _ in range(self.n):
            z = self.normal.simulacion()
            suma += z*z
        return suma

    def muestra(self, m):
        return [self.simulacion() for _ in range(m)]

resultados_chi2 = GeneradorChiCuadrada(grados_libertad=5).muestra(10)
## Muestras Chi-cuadrada df=5:
## [3.48645258705134, 2.979149258311065, 4.474218806855241, 1.7893772241680626, 12.966314136000925, 2.5688324423980404, 11.9944310132163, 4.477538796423481, 3.483840630072731, 1.2831985286592444]

3.14 Generador de variables Gamma \(\Gamma(k, \theta)\)

Este generador simula variables aleatorias con distribución Gamma con forma entera positiva \(k\) y escala \(\theta\).
Se basa en la propiedad de que una variable Gamma con \(k \in \mathbb{Z}^+\) puede representarse como la suma de \(k\) variables exponenciales independientes: \[ X = \sum_{i=1}^{k} E_i, \quad E_i \sim \text{Exponencial}(\theta) \]

3.14.0.1 Parámetros usados:
  • \(k = 2\): parámetro de forma
  • \(\theta = 5\): parámetro de escala
class GeneradorGamma:
    """
    Simula muestras de la distribución Gamma(k, theta) con k entero positivo
    usando la suma de k exponentiales (Erlang) y el generador uniforme global.
    """
    def __init__(self, k, theta):
        self.k = int(k)
        self.theta = theta
        if self.k < 1:
            raise ValueError("El parámetro k debe ser entero positivo.")

    def simulacion(self):
        """Una muestra individual de Gamma(k, theta)."""
        suma = 0.0
        for _ in range(self.k):
            u = generador.muestra(1)[0]
            suma += -math.log(1 - u) * self.theta
        return suma

    def muestra(self, n):
        """Genera n muestras de la distribución Gamma(k, theta)."""
        return [self.simulacion() for _ in range(n)]

k = 2
theta = 5
muestras_gamma = GeneradorGamma(k, theta).muestra(10)
## Muestras Gamma(k=2, θ=5):
## [21.465229154128988, 18.75899859818425, 2.679820105737583, 4.6475763747124414, 14.329735434503124, 0.9143796875589865, 8.637781818685564, 2.357575822086353, 4.26166488168549, 6.1985534244500515]

3.15 Generador de variables Cauchy

La distribución de Cauchy es una distribución continua con colas pesadas, definida por su mediana \(x_0\) y parámetro de escala \(\gamma\).
No tiene media ni varianza definidas.

Su función de densidad es:

\[ f(x) = \frac{1}{\pi \gamma \left[1 + \left(\frac{x - x_0}{\gamma}\right)^2\right]}, \quad x \in \mathbb{R} \]

Para simular esta distribución se usa la transformación de una variable uniforme \(U \sim \mathcal{U}(0,1)\):

\[ X = x_0 + \gamma \cdot \tan\left[ \pi (U - 0.5) \right] \]

3.15.0.1 Parámetros usados:
  • \(x_0 = 0\) — mediana
  • \(\gamma = 1\) — escala
class GeneradorCauchy:
    def __init__(self, x0, gamma):
        self.x0 = x0
        self.gamma = gamma

    def simulacion(self):
        u = generador.muestra(1)[0]
        return self.x0 + self.gamma * math.tan(math.pi * (u - 0.5))

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]

muestras_cauchy = GeneradorCauchy(x0=0, gamma=9).muestra(10)
## Muestras Cauchy:
## [21.457306344150044, 139.22382924419946, 6.609392917407011, -11.251811589457848, 10.741688560177748, 5.063569251370891, -29.904770196536756, 1.5980605084223518, -4.077301885170503, 12.100971098062239]

3.16 Promedios acumulados de la distribución Cauchy

Aunque muchas distribuciones tienen una media bien definida a la cual los promedios de muestras convergen por la Ley de los Grandes Números, esto no ocurre con la distribución Cauchy.
Esto se debe a que la distribución Cauchy no tiene esperanza matemática finita.

En este experimento se observa cómo los promedios acumulados de \(n\) simulaciones de Cauchy fluctúan intensamente y no se estabilizan con el crecimiento de \(n\):

import matplotlib.pyplot as plt

x0, gamma, n_max = 0, 1, 10000
cauchy_gen = GeneradorCauchy(x0, gamma)
promedios_cauchy = []
suma = 0.0

for i in range(1, n_max + 1):
    x = cauchy_gen.simulacion()
    suma += x
    promedios_cauchy.append(suma / i)

plt.plot(range(1, n_max + 1), promedios_cauchy)
plt.title("Promedios Acumulados de la Distribución Cauchy")
plt.xlabel("n")
plt.ylabel("Promedio acumulado")
plt.grid(True)
plt.show()

3.17 Descripción del método de Aceptación-Rechazo

El método de aceptación-rechazo es una técnica para generar muestras de una distribución cuya función de densidad \(f_X(x)\) es complicada de muestrear directamente. Se basa en usar una densidad auxiliar \(f_Y(x)\), de la que sí es sencillo extraer muestras, y una constante \(c\) tal que

\[ f_X(x) \;\le\; c\,f_Y(x), \quad \forall\,x\in\text{soporte}(f_X). \]

El procedimiento general es:

  1. Elegir una función auxiliar \(f_Y(x)\) que cumpla:

    • Su soporte contenga al de \(f_X(x)\).
    • Sea fácil generar \(y \sim f_Y(x)\) (por ejemplo, \(f_Y(x)\sim\mathcal{U}(a,b)\)).

  2. Calcular la constante \[ c \;=\; \sup_{x} \left\{ \frac{f_X(x)}{f_Y(x)} \right\}. \] En la práctica, a menudo se evalúa \(f_X(x)/f_Y(x)\) en puntos estratégicos o se aproxima por búsqueda en una rejilla.

  3. Repetir hasta obtener \(n\) muestras:

    1. Generar \(y \sim f_Y(x)\).
    2. Generar \(u \sim \mathcal{U}(0,1)\).
    3. Calcular la razón \[ r \;=\; \frac{f_X(y)}{c\,f_Y(y)}. \]
    4. Si \(u \le r\), aceptar \(y\) como muestra de \(f_X\); de lo contrario, rechazar \(y\) y volver al paso 3.1.

Cada muestra aceptada \(y\) tiene probabilidad de ser tomada proporcional a \(f_X(y)\), de modo que la colección final de valores sigue la distribución objetivo \(f_X(x)\).

3.17.0.1 Ventajas y desventajas
  • Ventaja: Permite simular densidades complicadas sin necesidad de la transformada inversa ni de métodos especializados.
  • Desventaja: La eficiencia depende de la constante \(c\). Si \(c\) es muy grande (es decir, \(f_X\) está “muy por debajo” de \(c\,f_Y\)), la proporción de rechazos será alta y se desperdiciarán muchos candidatos.

3.18 Simulación de una distribución personalizada por aceptación-rechazo

Queremos simular la siguiente función de densidad definida en el intervalo \([-1, 1]\):

\[ f_X(x) = \frac{\sqrt{3}}{\pi\,\bigl(x^2 + x + 1\bigr)}, \quad x \in [-1, 1]. \]

Utilizamos como densidad auxiliar

\[ f_Y(x) = \begin{cases} 0.5, & -1 \le x \le 1,\\ 0, & \text{en otro caso}, \end{cases} \]

es decir, \(f_Y(x)\sim\mathcal{U}(-1,1)\).

Utilizamos una densidad auxiliar uniforme \(f_Y(x) \sim \mathcal{U}(-1,1)\) y aplicamos el método de aceptación-rechazo usando el generador MRG32k3a.

import math
def f_x(x):
    return (math.sqrt(3)) / (math.pi * (x**2 + x + 1))
def f_y(x):
    return 0.5 if -1 <= x <= 1 else 0.0
def integral_fx(a, b, pasos=10000):
    h = (b - a) / pasos
    suma = 0.0
    for i in range(pasos):
        xi = a + (i + 0.5) * h
        suma += f_x(xi)
    return h * suma

integral = integral_fx(-1, 1)
print(f"Integral aproximada de f_X(x) en [-1, 1] = {integral:.6f}")
## Integral aproximada de f_X(x) en [-1, 1] = 1.000000
x_max = -0.5
fx_max = f_x(x_max)
fy_val = f_y(x_max)
c = fx_max / fy_val
print(f"f_X({x_max}) = {fx_max:.6f}, f_Y({x_max}) = {fy_val}")
## f_X(-0.5) = 0.735105, f_Y(-0.5) = 0.5
print(f"Constante c = f_X/f_Y = {c:.6f}")
## Constante c = f_X/f_Y = 1.470210
class MRG32k3a:
    def __init__(self, x0=1, x1=1, x2=20250409, y0=1, y1=1, y2=1):
        self.m1 = 4294967087
        self.m2 = 4294944443
        self.a12 = 1403580
        self.a13n = 810728
        self.a21 = 527612
        self.a23n = 1370589
        self.s1 = [x2 % self.m1, x1 % self.m1, x0 % self.m1]
        self.s2 = [y2 % self.m2, y1 % self.m2, y0 % self.m2]

    def muestra(self, n=1):
        resultados = []
        for _ in range(n):
            p1 = (self.a12 * self.s1[1] - self.a13n * self.s1[2]) % self.m1
            self.s1 = [p1] + self.s1[:2]
            p2 = (self.a21 * self.s2[0] - self.a23n * self.s2[2]) % self.m2
            self.s2 = [p2] + self.s2[:2]
            z = (self.s1[0] - self.s2[0]) % self.m1
            resultados.append(z / self.m1 if z > 0 else (self.m1 - 1) / self.m1)
        return resultados
class GeneradorUniformeContinua:
    def __init__(self, a, b):
        self.min = min(a, b)
        self.max = max(a, b)

    def simulacion(self):
        u = generador_principal.muestra(1)[0]
        return self.min + (self.max - self.min) * u

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]

generador_principal = MRG32k3a()
uniforme = GeneradorUniformeContinua(-1, 1)

def aceptar_rechazar(n, c):
    muestras = []
    while len(muestras) < n:
        y = uniforme.muestra(1)[0]               
        u = generador_principal.muestra(1)[0]
        if u <= f_x(y) / (c * f_y(y)):
            muestras.append(y)
    return muestras
muestras = aceptar_rechazar(10000, c)
import numpy as np
import matplotlib.pyplot as plt

plt.hist(muestras, bins=20, range=(-1,1), density=True, alpha=0.6, color='k', edgecolor='w')
x_vals = np.linspace(-1, 1, 1000)
fx_vals = [f_x(x) for x in x_vals]
plt.plot(x_vals, fx_vals, color='lime', lw=2)
plt.title('Simulación por Aceptación-Rechazo de $f_X(x)$ con MRG32k3a')
plt.xlabel('x')
plt.ylabel('Densidad')
plt.grid(True)
plt.show()

3.19 Simulación de Gamma \(\Gamma(\alpha,1)\) por Aceptación–Rechazo

La distribución Gamma con parámetro de forma \(\alpha>0\) y escala \(\theta=1\) tiene función de densidad:

\[ f_{\Gamma}(x;\,\alpha) \;=\; \begin{cases} \dfrac{x^{\alpha - 1}\,e^{-x}}{\Gamma(\alpha)}, & x>0,\\ 0, & \text{en otro caso}. \end{cases} \]

Para simular \(X \sim \Gamma(\alpha,1)\) utilizamos un generador auxiliar \(Y\) definido por una mezcla: - Con probabilidad \(p = \dfrac{e\,\alpha}{e + \alpha}\), \(Y \sim \mathrm{Beta}(\alpha,1)\) en \((0,1)\). - Con probabilidad \(1 - p\), \(Y \sim \mathrm{Exponencial}(1)\) en \((1,\infty)\).

La densidad auxiliar correspondiente es:

\[ f_{Y}(y; \alpha) \;=\; \begin{cases} p \,\dfrac{y^{\alpha - 1}}{B(\alpha,1)} = p\,y^{\alpha - 1}, & 0<y\le 1,\\ (1-p)\,e^{-y}, & y > 1,\\ 0, & \text{en otro caso}, \end{cases} \] donde \(B(\alpha,1)=1\) y \[ p = \frac{e\,\alpha}{e + \alpha}. \]

El método de aceptación–rechazo acepta \(Y\) como muestra de \(\Gamma(\alpha,1)\) si:

\[ U \;\le\; \frac{f_{\Gamma}(Y;\,\alpha)}{f_{Y}(Y;\,\alpha)}, \quad U \sim \mathcal{U}(0,1). \]

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma as gamma_func
from scipy.stats import beta, expon
import math

def f_gamma(x, alpha):
    return (x**(alpha - 1) * math.exp(-x)) / gamma_func(alpha) if x > 0 else 0.0

def f_y_gamma(y, alpha):
    p = (math.e * alpha) / (math.e + alpha)
    if 0 < y <= 1:
        return p * (y**(alpha - 1))
    elif y > 1:
        return (1 - p) * math.exp(-y)
    else:
        return 0.0
def simular_Y(alpha, generador_uniforme):
    p = (math.e * alpha) / (math.e + alpha)
    u = generador_uniforme.muestra(1)[0]
    if u <= p:
        return beta.rvs(alpha, 1)
    else:
        return expon.rvs()
def aceptar_rechazar_gamma(n, alpha, generador_uniforme):
    muestras = []
    while len(muestras) < n:
        y = simular_Y(alpha, generador_uniforme)
        u = generador_uniforme.muestra(1)[0]
        if u <= f_gamma(y, alpha) / f_y_gamma(y, alpha):
            muestras.append(y)
    return np.array(muestras)
def graficar_gamma(muestras, alpha, bins=100):
    x = np.linspace(0, np.max(muestras), 1000)
    fx = [f_gamma(val, alpha) for val in x]
    plt.figure(figsize=(10, 5))
    plt.hist(muestras, bins=bins, density=True, alpha=0.6, color='k', label='Simulado')
    plt.plot(x, fx, 'lime', label=f'Teórica Gamma($\\alpha$={alpha},1)', lw=2)
    plt.title(f"Simulación Gamma($\\alpha$={alpha},\\theta=1) por Aceptación–Rechazo")
    plt.xlabel("x")
    plt.ylabel("Densidad")
    plt.grid(True)
    plt.legend()
    plt.xlim(0, np.percentile(muestras, 99))
    plt.show()
alpha = 0.6
n = 10000
muestras_gamma = aceptar_rechazar_gamma(n, alpha, uniforme)
graficar_gamma(muestras_gamma, alpha, bins=100)

print(f"Muestras Gamma(α={alpha},1): media aproximada = {np.mean(muestras_gamma):.4f}, varianza aproximada = {np.var(muestras_gamma):.4f}")
## Muestras Gamma(α=0.6,1): media aproximada = 0.4965, varianza aproximada = 0.3695

3.20 Generador de variables Beta \(\mathrm{Beta}(\alpha,\beta)\)

La distribución Beta con parámetros de forma enteros \(\alpha>0\) y \(\beta>0\) tiene densidad en \((0,1)\):

\[ f_{\,\mathrm{Beta}}(x;\,\alpha,\beta) \;=\; \frac{x^{\alpha - 1}\,(1 - x)^{\beta - 1}}{B(\alpha,\beta)}, \quad B(\alpha,\beta) = \frac{\Gamma(\alpha)\,\Gamma(\beta)}{\Gamma(\alpha + \beta)}. \]

Una manera de simular \(X \sim \mathrm{Beta}(\alpha,\beta)\) es usar la relación con la distribución Gamma: - Si \(U \sim \Gamma(\alpha,1)\) y \(V \sim \Gamma(\beta,1)\) son independientes, entonces \[ X = \frac{U}{\,U + V\,} \;\sim\; \mathrm{Beta}(\alpha,\beta). \] - Para \(\alpha,\beta\in\mathbb{Z}^+\), cada \(\Gamma(k,1)\) puede obtenerse como suma de \(k\) variables exponenciales iid de parámetro 1: \[ U = \sum_{i=1}^{\alpha} E_i, \quad V = \sum_{j=1}^{\beta} F_j, \quad E_i, F_j \sim \mathrm{Exp}(1). \]

3.20.0.1 Parámetros específicos
  • \(\alpha = 2\)
  • \(\beta = 5\)
import math
import numpy as np

class GeneradorBeta:
    def __init__(self, alpha, beta):
        if not (isinstance(alpha, int) and alpha > 0 and isinstance(beta, int) and beta > 0):
            raise ValueError("alpha y beta deben ser enteros positivos")
        self.alpha = alpha
        self.beta  = beta
        self.expo  = GeneradorExponencial(1)  # Exponencial con parámetro λ = 1

    def simulacion(self):
        # Generar U = suma de alpha variables Exp(1)
        sum_x = 0.0
        for _ in range(self.alpha):
            sum_x += self.expo.simulacion()
        # Generar V = suma de beta  variables Exp(1)
        sum_y = 0.0
        for _ in range(self.beta):
            sum_y += self.expo.simulacion()
        # X = U / (U + V) ~ Beta(alpha, beta)
        return sum_x / (sum_x + sum_y)

    def muestra(self, n):
        return [self.simulacion() for _ in range(n)]
gen_beta = GeneradorBeta(alpha=2, beta=5)
muestras_beta = gen_beta.muestra(10)
print("Muestras Beta(2,5):", muestras_beta)
## Muestras Beta(2,5): [0.14305765581915905, 0.45397318654392466, 0.14382747984538258, 0.14135660285305562, 0.4964655165493996, 0.23837913465831212, 0.3088257338841022, 0.13944713869064626, 0.5291457224114471, 0.4755434895169242]

4 Procesos Estocásticos

4.1 Proceso de Poisson Homogéneo

Un Proceso de Poisson Homogéneo con tasa \(\lambda > 0\) es un proceso estocástico \(\{N(t)\}_{t \ge 0}\) tal que:

  1. \(N(0) = 0\).
  2. Incrementos en intervalos disjuntos son independientes.
  3. Para \(s < t\), la diferencia \(N(t) - N(s)\) sigue una distribución Poisson con parámetro \(\lambda (t - s)\): \[ P\bigl(N(t) - N(s) = k\bigr) = \frac{e^{-\lambda (t - s)} \bigl[\lambda (t - s)\bigr]^k}{k!}, \quad k = 0,1,2,\dots \]
  4. Los tiempos entre llegadas (inter-arrivals) son iid \(\mathrm{Exp}(\lambda)\).

Para simular eventos en \([0, T]\), se genera una sucesión de intervalos exponenciales \[ \Delta_i = -\frac{1}{\lambda} \ln\bigl(1 - U_i\bigr), \quad U_i \sim \mathcal{U}(0,1), \] y se acumulan hasta superar \(T\).

import math

class ProcesoPoissonHomogeneo:
    def __init__(self, lam):
        self.lam = lam  

    def simulacion(self, T):
        """
        Simula los tiempos de llegada de eventos en [0, T]
        en un Proceso de Poisson Homogéneo con tasa lambda.
        """
        tiempos = []
        t = 0.0
        while True:
            u = generador.muestra(1)[0]
            delta = -math.log(1 - u) / self.lam   
            t += delta
            if t > T:
                break
            tiempos.append(t)
        return tiempos

    def muestra(self, m, T):
        """
        Genera m trayectorias completas en [0, T],
        devolviendo una lista de listas de tiempos de llegada.
        """
        return [self.simulacion(T) for _ in range(m)]
import matplotlib.pyplot as plt
proceso = ProcesoPoissonHomogeneo(lam=3)
T = 15                                   

tiempos_llegada = proceso.simulacion(T)
tiempos = [0.0] + tiempos_llegada + [T]
cuentas = list(range(0, len(tiempos_llegada) + 1)) + [len(tiempos_llegada)]
plt.figure(figsize=(8, 4))
plt.step(tiempos, cuentas, where='post')
plt.xlabel('Tiempo')
plt.ylabel('Número de eventos $N(t)$')
plt.title(f'Proceso de Poisson Homogéneo (λ={proceso.lam}, T={T})')
plt.ylim(0, max(cuentas) + 1)
## (0.0, 48.0)
plt.grid(True)
plt.show()

4.2 Proceso de Wiener (Movimiento Browniano)

El Proceso de Wiener (o Movimiento Browniano) \(\{W(t)\}_{t \ge 0}\) es un proceso estocástico continuo que satisface:

  1. \(W(0) = 0\).
  2. Incrementos independientes: para \(0 \le s < t\), \(W(t) - W(s)\) es independiente de la historia previa.
  3. Incrementos distribuidos normalmente:
    \[ W(t) - W(s) \;\sim\; \mathcal{N}\bigl(0,\,t - s\bigr). \]

Para simularlo en un intervalo \([0, T]\) con \(n_{\text{pasos}}\) pasos, dividimos el tiempo en subintervalos de tamaño \[ \Delta t \;=\; \frac{T}{n_{\text{pasos}}}. \] Entonces, definimos \[ W(t_{j}) \;=\; W(t_{j-1}) \;+\; \sqrt{\Delta t}\,Z_{j}, \quad Z_{j} \sim \mathcal{N}(0,1), \] donde \(t_{j} = j\,\Delta t\).

En este ejemplo utilizamos dos generadores pseudoaleatorios en paralelo para aplicar el método polar de Box–Muller y obtener las variables \(Z_{j}\):
- MRG32k3a como generador primario,
- GeneradorUniformeContinua (o un generador uniforme global) como generador secundario.

import numpy as np
import math

class ProcesoWiener:
    def __init__(self, T, n_pasos, n_trayectorias=1):
        self.T = T
        self.n_pasos = n_pasos  
        self.dt = T / n_pasos                    
        self.n_trayectorias = n_trayectorias     

    def simulacion(self):
        tiempos = np.linspace(0, self.T, self.n_pasos + 1)
        trayectorias = np.zeros((self.n_trayectorias, self.n_pasos + 1))
        generador_principal = MRG32k3a()
        generador_secundario = GeneradorUniformeContinua(0, 1)  

        for i in range(self.n_trayectorias):
            j = 1
            while j <= self.n_pasos:
                u1 = generador_principal.muestra(1)[0]  
                u2 = generador_secundario.muestra(1)[0]
                v1 = 2 * u1 - 1
                v2 = 2 * u2 - 1
                r2 = v1 * v1 + v2 * v2
                if r2 == 0 or r2 > 1:
                    continue
                factor = math.sqrt(-2 * math.log(r2) / r2)
                z1 = v1 * factor
                z2 = v2 * factor
                trayectorias[i, j] = trayectorias[i, j-1] + math.sqrt(self.dt) * z1
                j += 1
                if j <= self.n_pasos:
                    trayectorias[i, j] = trayectorias[i, j-1] + math.sqrt(self.dt) * z2
                    j += 1
        return tiempos, trayectorias
import matplotlib.pyplot as plt

T = 1.0
n_pasos = 50
n_trayectorias = 20
pw = ProcesoWiener(T=T, n_pasos=n_pasos, n_trayectorias=n_trayectorias)
t, W = pw.simulacion()
plt.figure(figsize=(8, 4))
for traj in W:
    plt.plot(t, traj)
plt.title('Trayectorias de Proceso de Wiener')
plt.xlabel('t')
plt.ylabel('W(t)')
plt.grid(True)
plt.show()

5 Algunas aplicaciones

5.1 Estimador de \(\pi\) por Monte Carlo con Intervalo de Confianza t-Student

Para estimar \(\pi\), usamos el método de Monte Carlo generando puntos uniformes en el cuadrado \([-1,1]^2\). Si \((X,Y)\) es uniforme en \([-1,1]^2\), entonces

\[ P\bigl(X^2 + Y^2 \le 1\bigr) \;=\; \frac{\text{área del círculo de radio 1}}{\text{área del cuadrado}} \;=\; \frac{\pi}{4}. \]

Por lo tanto, si \(N\) puntos caen dentro del círculo, una estimación de \(\pi\) es

\[ \hat{\pi} \;=\; 4 \,\frac{\#\{(X_i,Y_i) : X_i^2 + Y_i^2 \le 1\}}{N}. \]

Para obtener un intervalo de confianza al 95% sobre la media de \(\hat{\pi}\), repetimos el experimento \(n\) veces obteniendo estimaciones \(\hat{\pi}_1, \dots, \hat{\pi}_n\). Sea

\[ Y_n \;=\; \frac{1}{n} \sum_{i=1}^n \hat{\pi}_i, \qquad S_n^2 \;=\; \frac{1}{n-1} \sum_{i=1}^n (\hat{\pi}_i - Y_n)^2. \]

Con \(t_{\alpha/2,\,n-1}\) el cuantil de la distribución t de Student con \(n-1\) grados de libertad, el margen de error para un nivel de confianza del \(100(1-\alpha)\%\) es

\[ \text{margen} \;=\; t_{\,0.975,\,n-1} \,\frac{S_n}{\sqrt{n}}. \]

El intervalo de confianza al 95% para \(\pi\) es entonces

\[ \bigl(Y_n - \text{margen},\, Y_n + \text{margen}\bigr). \]

5.1.0.1 Parámetros específicos
  • \(N\): número de puntos por repetición (cada \(\hat{\pi}_i\) se basa en \(N\) puntos).
  • \(n\): número de repeticiones para calcular el IC.
  • Generador de números pseudoaleatorios: \(\texttt{MRG32k3a}\).
import math
import numpy as np
from scipy.stats import t

class EstimadorPi:
    """
    Estima π usando Monte Carlo en [-1,1]^2,
    calculando un intervalo de confianza al 95% vía t-Student.
    """
    def __init__(self, generador):
        self.generador = generador

    def estimar_pi(self, N):
        dentro = 0
        for _ in range(N):
            x = 2 * self.generador.muestra(1)[0] - 1  
            y = 2 * self.generador.muestra(1)[0] - 1
            if x**2 + y**2 <= 1:
                dentro += 1
        return 4 * dentro / N

    def experimento(self, N_puntos, n_repeticiones):
        estimaciones = [self.estimar_pi(N_puntos) for _ in range(n_repeticiones)]
        Yn = np.mean(estimaciones)
        Sn2 = np.var(estimaciones, ddof=1)
        Sn = math.sqrt(Sn2)
        t_critico = t.ppf(0.975, df=n_repeticiones - 1)
        margen = t_critico * Sn / math.sqrt(n_repeticiones)
        IC = (Yn - margen, Yn + margen)
        return {
            "media": Yn,
            "varianza_muestral": Sn2,
            "intervalo_confianza": IC,
            "estimaciones": estimaciones
        }

    def analizar_varianza_conforme_n(self, N_puntos, max_n):
        varianzas = []
        ns = range(2, max_n + 1)
        for n in ns:
            resultados = self.experimento(N_puntos, n)
            varianzas.append(resultados["varianza_muestral"])
        return ns, varianzas
generador_pi = MRG32k3a()
estimador = EstimadorPi(generador_pi)
resultado = estimador.experimento(N_puntos=1000, n_repeticiones=50)
print(f"Estimación media de π: {resultado['media']:.6f}")
## Estimación media de π: 3.145120
print(f"Varianza muestral: {resultado['varianza_muestral']:.6f}")
## Varianza muestral: 0.002043
print(f"Intervalo de confianza 95%: {resultado['intervalo_confianza']}")
## Intervalo de confianza 95%: (np.float64(3.132272998495353), np.float64(3.157967001504647))

5.2 Simulación de un Sistema de Colas M/M/1

Un sistema M/M/1 representa un proceso de colas con:

  1. Llegadas según un proceso de Poisson con tasa \(\lambda > 0\).
  2. Tiempos de servicio exponenciales con tasa \(\mu > 0\).
  3. Servidor único.
  4. Disciplina FIFO (primero en entrar, primero en salir).

Los tiempos entre llegadas \(\{\tau_i\}\) son iid con distribución \(\mathrm{Exp}(\lambda)\):

\[ \tau_i = -\frac{1}{\lambda} \ln(1 - U_i), \quad U_i \sim \mathcal{U}(0,1). \]

Los tiempos de servicio \(\{s_i\}\) son iid con distribución \(\mathrm{Exp}(\mu)\):

\[ s_i = -\frac{1}{\mu} \ln(1 - V_i), \quad V_i \sim \mathcal{U}(0,1). \]

Para simular hasta un tiempo máximo \(T\), se alternan eventos de llegada y salida:

  • Si la próxima llegada ocurre antes que la próxima salida, se genera un evento de llegada.
  • En caso contrario, se genera un evento de salida.

Se registran los tiempos de espera en la cola para cada cliente \(i\),
donde el tiempo de espera \(W_i = t_{\text{inicio servicio},i} - t_{\text{llegada},i}\).

import math
import numpy as np
generador_principal = MRG32k3a()

class Exponencial:
    """
    Genera un tiempo exponencial con parámetro lambda
    usando el generador uniforme MRG32k3a.
    """
    def __init__(self, parametro_lamb):
        self.parametro_lamb = parametro_lamb

    def simulacion(self):
        u = generador_principal.muestra(1)[0]  
        return -math.log(1 - u) / self.parametro_lamb


class SistemaColasMM1:
    """
    Simula un sistema de colas M/M/1 hasta tiempo_max:
    - tasa_llegada = λ    (Procesos de Poisson)
    - tasa_servicio = μ   (Exp(μ) para servicio)
    """
    def __init__(self, tasa_llegada, tasa_servicio, generador):
        self.lambda_   = tasa_llegada
        self.mu        = tasa_servicio
        self.generador = generador
        self.llegadas  = Exponencial(self.lambda_)
        self.servicios = Exponencial(self.mu)

    def simular(self, tiempo_max):
        tiempo           = 0.0
        cliente          = 0
        servidor_ocupado = False
        fila             = []
        tiempos_espera   = []
        siguiente_llegada = self.llegadas.simulacion()
        siguiente_salida  = float('inf')

        while tiempo < tiempo_max:
            if siguiente_llegada < siguiente_salida:
                tiempo = siguiente_llegada
                cliente += 1
                if not servidor_ocupado:
                    servidor_ocupado = True
                    siguiente_salida = tiempo + self.servicios.simulacion()
                    tiempos_espera.append(0.0)
                else:
                    fila.append(tiempo)
                siguiente_llegada = tiempo + self.llegadas.simulacion()
            else:
                tiempo = siguiente_salida
                if fila:
                    inicio   = fila.pop(0)
                    espera   = tiempo - inicio
                    tiempos_espera.append(espera)
                    siguiente_salida = tiempo + self.servicios.simulacion()
                else:
                    servidor_ocupado = False
                    siguiente_salida  = float('inf')

        return {
            "clientes_atendidos": cliente,
            "tiempos_espera": tiempos_espera,
            "prom_espera": np.mean(tiempos_espera) if tiempos_espera else 0.0,
            "max_espera": max(tiempos_espera) if tiempos_espera else 0.0
        }

generador = generador_principal
sistema = SistemaColasMM1(tasa_llegada=0.5, tasa_servicio=0.8, generador=generador)
resultado = sistema.simular(tiempo_max=120)
print("Clientes atendidos:    ", resultado["clientes_atendidos"])
## Clientes atendidos:     48
print("Promedio de espera:    ", resultado["prom_espera"])
## Promedio de espera:     1.8298604435737016
print("Máximo tiempo de espera:", resultado["max_espera"])
## Máximo tiempo de espera: 9.326055689964058

5.3 Comparación de Estimadores de \(p\) (Frecuentista vs Bayesiano)

Dado un parámetro de probabilidad \(p \in (0,1)\) para una distribución de Bernoulli, queremos comparar dos estimadores:

  1. Estimador frecuentista: \[ \hat{p}_1 = \frac{k}{n}, \quad k = \sum_{i=1}^n X_i,\quad X_i \sim \text{Bernoulli}(p). \]

  2. Estimador bayesiano con prior \(\mathrm{Beta}(1,1)\) (uniforme): \[ \hat{p}_2 = \frac{k + 1}{\,n + 2\,}. \]

En cada repetición, generamos \(n\) muestras de una Bernoulli con parámetro \(p\), calculamos \(k\), y obtenemos \(\hat{p}_1\) y \(\hat{p}_2\).
Luego repetimos el experimento \(R\) veces para observar la variabilidad de ambos estimadores.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def comparar_estimadores_con_zoom(p, n, repeticiones):
    p1_list = []
    p2_list = []
    
    for _ in range(repeticiones):
        gen_bernoulli = GeneradorBernoulli(p)
        muestras = gen_bernoulli.muestra(n)
        k = sum(muestras)
        # Estimador frecuentista
        p1 = k / n
        # Estimador bayesiano con prior Beta(1,1)
        p2 = (k + 1) / (n + 2)
        p1_list.append(p1)
        p2_list.append(p2)
    
    
    indices = np.arange(1, repeticiones + 1)
    y_min = min(min(p1_list), min(p2_list))
    y_max = max(max(p1_list), max(p2_list))
    margen = 0.02
    y_min_zoom = max(0, y_min - margen)
    y_max_zoom = min(1, y_max + margen)
    
    plt.figure(figsize=(10, 5))
    plt.scatter(indices, p1_list, color='blue', s=10, alpha=0.6, label='p1 (frecuentista)')
    plt.scatter(indices, p2_list, color='red',  s=10, alpha=0.6, label='p2 (bayesiano)')
    plt.axhline(y=p, color='green', linestyle='--', linewidth=2, label=f'p real = {p}')
    plt.ylim(y_min_zoom, y_max_zoom)
    plt.title(f"Estimaciones de $\\hat{{p}}_1$ y $\\hat{{p}}_2$ (p={p}, n={n}, repeticiones={repeticiones})")
    plt.xlabel("Repetición")
    plt.ylabel("Estimador de $p$")
    plt.legend()
    plt.grid(True)
    plt.show()
    
    return p1_list, p2_list
p_true = 0.3
n_muestras = 50
n_reps = 200
p1_list, p2_list = comparar_estimadores_con_zoom(p_true, n_muestras, n_reps)

indices = np.arange(1, n_reps + 1)
tabla_estimaciones = pd.DataFrame({
    "Iteración": indices,
    "p1 (frecuentista)": p1_list,
    "p2 (bayesiano)"   : p2_list
})

print(tabla_estimaciones)
##      Iteración  p1 (frecuentista)  p2 (bayesiano)
## 0            1               0.32        0.326923
## 1            2               0.32        0.326923
## 2            3               0.28        0.288462
## 3            4               0.30        0.307692
## 4            5               0.36        0.365385
## ..         ...                ...             ...
## 195        196               0.32        0.326923
## 196        197               0.22        0.230769
## 197        198               0.36        0.365385
## 198        199               0.36        0.365385
## 199        200               0.30        0.307692
## 
## [200 rows x 3 columns]

5.3.1 Conclusión sobre los estimadores de \(p\)

En la gráfica obtenida para \(p=0.3\), \(n=50\) y \(R=200\) repeticiones, observamos las estimaciones:

  • \(\hat{p}_1 = \dfrac{k}{n}\) (frecuentista), marcado en azul.
  • \(\hat{p}_2 = \dfrac{k + 1}{\,n + 2\,}\) (bayesiano con prior \(\mathrm{Beta}(1,1)\)), marcado en rojo.
  • La línea horizontal discontinua verde indica el valor real \(p = 0.3\).

A partir de la dispersión de puntos, podemos extraer las siguientes conclusiones:

  1. Sesgo de los estimadores
    • El estimador frecuentista \(\hat{p}_1\) fluctúa alrededor de 0.3 y, en promedio, muestra un sesgo muy pequeño.
    • El estimador bayesiano \(\hat{p}_2 = \frac{k+1}{n+2}\) tiende a ser ligeramente “empujado” hacia el valor 0.5 cuando \(k\) es pequeño, y hacia 0 cuando \(k\) es cercano a \(n\). Sin embargo, para \(n=50\) este efecto es leve: las nubes rojas (puntos bayesianos) están casi superpuestas a las azules (puntos frecuentistas), indicando que ambos estimadores son prácticamente equivalentes en promedio.
  2. Variabilidad (dispersión)
    • Tanto \(\hat{p}_1\) como \(\hat{p}_2\) muestran dispersión en torno a 0.3; la mayoría de las estimaciones se encuentran en el intervalo aproximado \([0.20,\;0.40]\).
    • A simple vista, las nubes azul y roja tienen amplitud similar, lo que sugiere que la varianza muestral de ambos estimadores es comparable cuando \(n=50\).
  3. Ventajas del estimador bayesiano para \(n\) pequeño
    • En las pocas repeticiones donde \(k=0\) (ningún “éxito” en las 50 pruebas), \(\hat{p}_1 = 0\) mientras que \(\hat{p}_2 = \tfrac{1}{52} \approx 0.019\). Este empuje hacia valores intermedios evita estimaciones exactas en los extremos (\(0\) o \(1\)), lo cual puede ser útil si \(n\) fuera considerablemente más pequeño.
    • De igual forma, cuando \(k = 50\) (todos “éxitos”), \(\hat{p}_1 = 1\) pero \(\hat{p}_2 = \tfrac{51}{52} \approx 0.981\). Nuevamente, el estimador bayesiano “suaviza” la estimación en los casos extremos.
  4. Comportamiento alrededor del valor real
    • Ambas estimaciones se agrupan alrededor de la línea \(p = 0.3\).
    • No se aprecia un sesgo sistemático fuerte: la nube de puntos se distribuye simétricamente por encima y por debajo de 0.3.
    • La línea verde (valor real) cruza la distribución de puntos en lugares donde la densidad de puntos es mayor, confirmando que los dos métodos son consistentes.
  5. Interpretación práctica
    • Para \(n = 50\) y \(p = 0.3\), la diferencia entre \(\hat{p}_1\) y \(\hat{p}_2\) es mínima.
    • Si el tamaño de muestra fuera mucho menor, el método bayesiano aportaría estimaciones menos extremas (evitando 0 o 1).
    • En escenarios de muestras moderadamente grandes (\(n \ge 50\)), el estimador frecuentista es adecuado y más sencillo, pero el bayesiano sigue siendo una buena alternativa si se desea incorporar un sesgo “suave” inicial.

6 Ecuación Diferencial Estocástica

6.1 Simulación de una Ecuación Diferencial Estocástica (SDE)

Queremos simular la SDE de la forma

\[ dX_t \;=\; f(t,\,X_t)\,dt \;+\; g(t,\,X_t)\,dW_t, \]

donde \(W_t\) es un proceso de Wiener estándar. Usaremos el esquema de Euler–Maruyama y el método de Box–Muller para generar los incrementos \(dW\).

  • Dividimos el intervalo \([0, T]\) en \(N\) pasos de tamaño \(\Delta t = T/N\).
  • El esquema de Euler–Maruyama aproxima: \[ X_{t_{i}} \;=\; X_{t_{i-1}} \;+\; f\bigl(t_{i-1},\,X_{t_{i-1}}\bigr)\,\Delta t \;+\; g\bigl(t_{i-1},\,X_{t_{i-1}}\bigr)\,\Delta W_{i}, \] donde \(\Delta W_{i} = \sqrt{\Delta t}\,Z_{i}\) y \(Z_{i} \sim \mathcal{N}(0,\,1)\).

Para obtener \(Z_{i}\), usamos dos generadores globales: - MRG32k3a como generador principal (uniformes \(U_1\)).
- LFSR20 como generador secundario (uniformes \(U_2\)).

Ambos producen números en \((0,1)\). Aplicamos el método polar de Box–Muller:

  1. Generamos \(U_1,\,U_2 \sim \mathcal{U}(0,1)\).
  2. Definimos \(V_1 = 2\,U_1 - 1\), \(V_2 = 2\,U_2 - 1\).
  3. Calculamos \(R^2 = V_1^2 + V_2^2\); si \(R^2 \in (0,1]\),
    \[ Z = V_1 \;\sqrt{\frac{-2 \ln R^2}{R^2}} \;\sim\; \mathcal{N}(0,1). \]
import numpy as np
import math
import matplotlib.pyplot as plt
generador_principal   = MRG32k3a()
generador_secundario  = LinearFeedbackGeneratorSize20(87685)

class SimuladorSDE:
    """
    Simula la SDE dX = f(t, X)·dt + g(t, X)·dW usando
    Euler–Maruyama y Box–Muller para generar los incrementos dW.

    Parámetros:
      - f: función determinista f(t, x)
      - g: función de difusión g(t, x)
      - x0: valor inicial X(0)
      - T: horizonte de tiempo
      - N: número de pasos discretos
      - semilla: (opcional) semilla entera para reproducibilidad
    """

    def __init__(self, f, g, x0, T, N, semilla=None):
        self.f       = f
        self.g       = g
        self.x0      = x0
        self.T       = T
        self.N       = N
        self.dt      = T / N
        self.tiempo  = np.linspace(0, T, N + 1)
        self.semilla = semilla

    def simular(self):
        if self.semilla is not None:
            generador_principal.__init__(
                x0=self.semilla,
                x1=self.semilla + 1,
                x2=self.semilla + 2,
                y0=self.semilla + 3,
                y1=self.semilla + 4,
                y2=self.semilla + 5
            )
            generador_secundario.__init__(semilla=self.semilla + 12345)

        X = np.zeros(self.N + 1)
        X[0] = self.x0

        for i in range(1, self.N + 1):
            t = self.tiempo[i - 1]
            while True:
                u1 = generador_principal.muestra(1)[0]
                u2 = generador_secundario.muestra(1)[0]
                v1 = 2 * u1 - 1
                v2 = 2 * u2 - 1
                r2 = v1 * v1 + v2 * v2
                if 0 < r2 <= 1:
                    break
            z  = v1 * math.sqrt(-2 * math.log(r2) / r2) 
            dW = z * math.sqrt(self.dt)

            X[i] = (
                X[i - 1]
                + self.f(t, X[i - 1]) * self.dt
                + self.g(t, X[i - 1]) * dW
            )

        return self.tiempo, X

    def graficar(self, tiempo, X, titulo="Simulación de SDE"):
        plt.figure(figsize=(10, 5))
        plt.plot(tiempo, X, label="Trayectoria simulada")
        plt.xlabel("Tiempo")
        plt.ylabel("X(t)")
        plt.title(titulo)
        plt.grid(True)
        plt.legend()
        plt.show()

def f_gbm(t, x):
    mu = 0.1
    return mu * x

def g_gbm(t, x):
    sigma = 0.3
    return sigma * x
sim_gbm = SimuladorSDE(
    f=f_gbm,
    g=g_gbm,
    x0=1.0,
    T=1.0,
    N=500,
    semilla=2025
)

tiempo_gbm, X_gbm = sim_gbm.simular()
sim_gbm.graficar(tiempo_gbm, X_gbm, titulo="Geometric Brownian Motion (Euler–Maruyama)")