Generación de números aleatorios

Los números aleatorios son un ingrediente básico en la simulación de todo sistema discreto. La mayoría de lenguajes de programación poseen una subrutina, objeto o función que generan números aleatorios. De la misma forma los softwares de simulación generan números aleatorios que son usados para generar variables aleatorias.

Propiedades de los números aleatorios.

Una secuencia de números aleatorios \(R_1, R_2,..,\) tiene dos importantes propiedades estadísticas:

1. Uniformidad

2. Independencia.

Cada número aleatorio \(R_i\) debe ser una muestra independiente extraída de una distribución uniforme continua entre cero y uno. La función de densidad de probabilidad está dada por la Ecuación 1:

\[\begin{align} f(x) = \left\{ \begin{array}{lr} 1, & 0 \leq x \leq 1 \\ 0, & de~lo~contrario \end{array} \right. \end{align} \tag{1}\]

La función de densidad de probabilidad se muestra en Figura 1:

Figura 1: Función de densidad de probabilidad para números aleatorios

El valor esperado para cada \(R_i\), \(E(R)\) está dado por la Ecuación 2:

\[\begin{align} E(R) &= \int_0^1 x~\mathrm{d}x \\ E(R) &= \frac{x^2}{2} \Big|_0^1 \\ E(R) &= \frac{1}{2} \end{align} \tag{2}\]

La varianza \(V(R)\) está dada por la Ecuación 3:

\[\begin{align} V(R) &= \int_0^1 x^2~\mathrm{d}x - [E(R)]^2 \\ V(R) &= \frac{x^3}{3} \Big|_0^1 - \left[ \frac{1}{2} \right]^2\\ V(R) &= \frac{1}{3} - \frac{1}{4} = \frac{1}{12} \end{align} \tag{3}\]

Algunas consecuencias de las propiedades de uniformidad e independencia son las siguientes:

1. Si el intervalo \([0, 1]\) es dividido en \(n\) clases o subintervalos de igual longitud, el número esperado de observacioenes en cada intervalo será igual a \(\frac{N}{n}\) donde \(N\) es el número total de observaciones.

2.La probabilidad de observar un valor en un intervalo determinado es independiente de los valores anteriores extraídos.

Generación de números

En la práctica se utilizan número pseudoaleatorios en los procesos de simulación. La palabra pseudo se utiliza para dar a entender que generar números aleatorios mediante algún método elimina la posibilidad de una verdadera aleatoriedad. Si se conoce el algoritmo y la semilla, el conjunto de números aleatorios \(R_i\) se podrá reproducir.

El propósito, igualmente, es reprodcir una secuencia de números entre \(0\) y \(1\) que simule o imite las propiedades de distribución uniforme e independencia tanto como sea posible.

En la generación de números pesudoaleatorios, cierto tipo de problemas pueden ocurrir. Estos errores, o desviaciones a la aleatoriedad, están todos relacionados con las propiedades previamente expuestas, algunos ejemplos de los errores son los siguientes:

1. Los números aleatorios generados podrían no estar uniformemente distribuídos.

2. Los números aleatorios generados podrían poseer una distribución discreta y no continua.

3. Los números aleatorios generados podrían tener media muy alta o muy baja.

4. La varianza de los números aleatorios generados podría ser muy alta o muy baja.

5. Podría haber dependencia.

Las desviaciones con respecto a la uniformidad y la independencia de un determinado esquema o algoritmo de generación a menudo pueden detectarse mediante pruebas estadísticas. Si se detectan tales desviaciones, el algoritmo de generación debe abandonarse en favor de un generador aceptable.

Un buen algoritmo de generación de números aleatorios debería cumplir:

1. Tener un ciclo lo suficientemente amplio. La longitud del ciclo o período de un algoritmo de generación de número pseudoaleatorios representa la longitud de la secuencia de números aleatorios antes de que los números empiecen a repetirse. DDe esta manera se se van a generar \(10.000\) eventos, el ciclo o período debería ser varias veces mayor.

2. Los números aleatorios deben poder reproducirse. Dado el punto de partida (o las condiciones), debe ser posible generar el mismo conjunto de números aleatorios, con total independencia del sistema que se esté simulando. Esto es útil a efectos de depuración y es un medio para facilitar las comparaciones entre sistemas.

3. Los números aleatorios generados deben aproximarse mucho a las propiedades estadísticas ideales de uniformidad e independencia.

4. La generación no debería representar mayor costo compuacional y debería ser rápida.

Métodos para generar números aleatorios.

Para poder generar entradas probabilísticas en un modelo de simulación se debe contar con un generador de números pseudoaleatorios. Los métodos para generar números pseudoaleatorios se clasifican de acuerdo con la forma en cómo son obtenidos.

Métodos físicos, manuales y tablas.

Han sido los métodos precursores para generar números pseudoaleatorios. Los números generados através de estos métodos tienen la característica de que es imposible reproducir su secuencia. Dentro de los métodos físicos se encuentran:

1. Discos aleaotorios.

2. Lanzamiento de monedas, dados, ruletas, cartas.

3. Tablas de números aleatorios.

Métodos para generar números pseudoaleatorios.

Métodos aritméticos

Método de los cuadrados medios

Método propuesto por Jhon Von Neumann que requiere una semilla y genera una secuencia de números \(R_i\) de la siguiente manera.

1. Establecer la semilla \(x_0\). \(x_0\) debe contener una cantidad de dígitos mayor a \(3\).

2. Hallar \(z_0\)

\[\begin{align} z_0 = x_0^2 \end{align}\]

3. Se toman los dígitos centrales de \(z_0\), generalmente \(4\) dígitos. Corresponde a la semilla para el siguiente número pseudoaleatorio.

\[\begin{align} x_1=digitos(z_{n-1}) \end{align}\]

4. Hallar \(R_i\)

\[\begin{align} R_i = \frac{x_i}{10000} \end{align}\]

El valor del denominador de la división puede cambiar en múltiplos de \(10\)

4. Repetir literales a a c, teniendo en cuenta el último \(R_i\) generado

Ejemplo de generación de pseudoaleatorios usando cuadrados medios.

Generar los primeros \(3\) números aleatorios a partir de \(x_0=3450\)

Se realiza el procedimiento para hallar \(R_1\).

1. \(x_0=3450\)

2. \(z_0 = 3450^2 = 11902500\)

3. Se escogen los dígitos: \(x_1=9025\)

4. \(R_1 = \frac{9025}{10000} = 0.9025\)

Se realiza el procedimiento para hallar \(R_2\).

1. \(x_1=9025\)

2. \(z_1 = 9025^2 = 81450625\)

3. Se escogen los dígitos: \(x_2=4506\)

4. \(R_2 = \frac{4506}{10000} = 0.4506\)

Se realiza el procedimiento para hallar \(R_3\).

1. \(x_2=4506\)

2. \(z_2 = 4506^2 = 20304036\)

3. Se escogen los dígitos: \(x_3=3040\)

4. \(R_3 = \frac{3040}{10000} = 0,3040\)

Desventaja del método no es capaz de generar una secuencia amplia de números pseudoaleatorios, se degenera de manera, relativamente rápida.

Algoritmo en R

#Definimos una función llamada cuadrados_medios.
#La función recibe: semilla: valor inicial para comenzar la generación, n: cantidad de números pseudoaleatorios que queremos generar.

cuadrados_medios <- function(semilla, n) {
  #Convertimos la semilla a texto mediante as.character()
  #y contamos cuántos dígitos tiene con nchar().
  #Por ejemplo, si semilla = 5731, entonces d = 4.
  d <- nchar(as.character(semilla))
  
  #Creamos un vector llamado semillas para almacenar
  #todas las semillas que se obtendrán durante el proceso.
  #Se necesitan n + 1 posiciones porque almacenaremos:
  #la semilla inicial + las n nuevas semillas.
  semillas <- numeric(n + 1)
  #Guardamos la semilla inicial en la primera posición del vector semillas.
  semillas[1] <- semilla
  #Creamos un vector llamado números para almacenar los números pseudoaleatorios generados.
  #Estos serán los valores R_i que estarán entre 0 y 1.
  numeros <- numeric(n)
  #Iniciamos un ciclo for que se repetirá n veces. Cada repetición genera una nueva semilla y un nuevo número pseudoaleatorio.
  for (i in 1:n) {
    #Elevamos al cuadrado la semilla actual. Por ejemplo: # 5731^2 = 32844361
    cuadrado <- semillas[i]^2
    # Convertimos el resultado del cuadrado a texto. Además, sprintf() permite agregar ceros a la izquierda para garantizar que el número tenga exactamente 2*d dígitos. Si d = 4, necesitamos que el cuadrado tenga 8 dígitos.
    cuadrado <- sprintf(paste0("%0", 2 * d, "d"), cuadrado)
    #Calculamos la posición donde comienzan los d dígitos centrales. Para una semilla de 4 dígitos: d/2 + 1 = 3. Por tanto, comenzaremos a extraer desde el tercer dígito.
    inicio <- d / 2 + 1
    # Calculamos la posición donde terminan los d dígitos centrales. Si d = 4:  inicio = 3, fin = 3 + 4 - 1 = 6, por tanto, extraemos los dígitos de la posición 3 a la 6.
    fin <- inicio + d - 1
    # Extraemos los d dígitos centrales del cuadrado. substr() permite extraer una parte de una cadena de texto.Por ejemplo: 32844361, Los cuatro dígitos centrales son: 8443
    nueva_semilla <- substr(cuadrado, inicio, fin)
    # Convertimos la nueva semilla de texto a número. Por ejemplo: "8443" -> 8443, Esto permite utilizarla en la siguiente iteración.
    semillas[i + 1] <- as.numeric(nueva_semilla) 
    #Transformamos la nueva semilla en un número pseudoaleatorio entre 0 y 1. Para una semilla de 4 dígitos, U_i = X_i / 10^4, Por ejemplo: 8443 / 10000 = 0.8443
    numeros[i] <- semillas[i + 1] / 10^d
  }
  
  #Creamos un data.frame para presentar los resultados de una manera organizada.
  resultados <- data.frame(
    #Número de la iteración. Va desde 1 hasta n.
    Iteracion = 1:n,
    #Semilla utilizada al comienzo de cada iteración.
    Semilla = semillas[1:n],
    #Cuadrado de la semilla.
    Cuadrado = semillas[1:n]^2,
    #Nueva semilla obtenida a partir de los dígitos centrales del cuadrado.
    Nueva_semilla = semillas[2:(n + 1)],
    #Número pseudoaleatorio generado entre 0 y 1.
    Numero = numeros
  )
  
  #La función devuelve el data.frame con todos los resultados.
  return(resultados)
}

#Usamos la función:
cuadrados_medios(semilla=3450, n=3)

  Iteracion Semilla Cuadrado Nueva_semilla Numero
1         1    3450 11902500          9025 0.9025
2         2    9025 81450625          4506 0.4506
3         3    4506 20304036          3040 0.3040

Método de los productos medios

Corresponde a un método bastante similar al de los cuadrados medios con la diferencia que se inicia con dos semillas. El procedimiento de generación de números pseudoaleatorios mediante el método se muestra a continuación:

Para iteraciones \(i=0,1,2,3... m\) y semillas \(j=1,2\):

1. Se definen los valores \(x_{0,1}\) y \(x_{0,2}\) como semillas iniciales, con cantidad de dígitos mayor a \(3\).

2. Se halla producto \(z_0\)

\[\begin{align} z_0 = x_{0,1}*x_{0,1} \end{align}\]

En general el producto se puede hallar como:

\[\begin{align} z_i = x_{i,1}*x_{i,2} \end{align}\]

3. Se halla \(x_1\) tomando los dígitos centrales de \(z_0\), generalmente \(4\) dígitos.

\[\begin{align} x_1 = digitos(z_0) \end{align}\]

En general:

\[\begin{align} x_{i+1}= digitos(z_i) \end{align}\]

4. Hallar \(R_{i+1}\)

\[\begin{align} R_{i+1} = \frac{x_{i+1}}{10000} \end{align}\]

El valor del denominador de la división puede cambiar en múltiplos de \(10\)

5. Nuevas semillas.

Se establecen las nuevas semillas como:

\[\begin{align} x_{i+1,~1} &= x_{i,~2}\\ x_{i+1,~2} &= x_{i+1} \end{align}\]

Iterar para el siguiente \(i\)

Ejemplo de generación de pseudoaleatorios usando productos medios.

Generar los primeros \(3\) números aleatorios a partir de \(x_{0,0}=4725\) y \(x_{0,1}=5250\)

1. Semillas iniciales

Para \(i=0\)

\[\begin{align} x_{0,1}=4725\\ x_{0,2}=5250 \end{align}\]

2. Se halla producto \(z_0\)

\[\begin{align} z_0 &= x_{0,1}*x_{0,2}\\ z_0 &= 4725*5250 = 24806250 \end{align}\]

3. Se halla \(x_1\) tomando los dígitos centrales de \(z_0\), generalmente \(4\) dígitos.

\[\begin{align} x_1 = 8062 \end{align}\]

4. Hallar \(R_1\)

\[\begin{align} R_1 = \frac{8062}{10000}=0.8062 \end{align}\]

5. Nuevas semillas.

\[\begin{align} x_{1,1} &= 5250\\ x_{1,2} &= 8062 \end{align}\]

Para \(i = 1\)

1. Semillas

\[\begin{align} x_{1,1} &= 5250\\ x_{1,2} &= 8062 \end{align}\]

2. \(z_1\)

\[\begin{align} z_1 &= x_{1,1}*x_{1,2}\\ z_1 &= 5250*8062 = 42325500 \end{align}\]

3. \(x_2\)

\[\begin{align} x_{2}= 3255 \end{align}\]

4. Hallar \(R_{2}\)

\[\begin{align} R_{2} &= \frac{x_{2}}{10000}\\ R_{2} &= \frac{3255}{10000} = 0,3255 \end{align}\]

5. Nuevas semillas.

\[\begin{align} x_{2,~1} &= x_{1,~2} = 8062\\ x_{2,~2} &= x_{2} = 3255 \end{align}\]

Para \(i = 2\)

1. Semillas

\[\begin{align} x_{2,~1} &= 8062\\ x_{2,~2} &= 3255 \end{align}\]

2. \(z_2\)

\[\begin{align} z_2 &= x_{2,1}*x_{2,2}\\ z_2 &= 8062*3255 = 26241810 \end{align}\]

3. \(x_3\)

\[\begin{align} x_{3}= 2418 \end{align}\]

4. Hallar \(R_{3}\)

\[\begin{align} R_{3} &= \frac{x_{3}}{10000}\\ R_{3} &= \frac{2418}{10000} = 0,2418 \end{align}\]

Algoritmo en R

#Definimos una función llamada productos_medios.
#La función recibe:
#semilla1: primera semilla inicial.
#semilla2: segunda semilla inicial.
#n: cantidad de números pseudoaleatorios que queremos generar.

productos_medios <- function(semilla1, semilla2, n) {
  
  #Convertimos ambas semillas a texto mediante as.character()
  #y contamos cuántos dígitos tiene cada una con nchar().
  #El método requiere que ambas semillas tengan la misma cantidad de dígitos.
  d1 <- nchar(as.character(semilla1))
  d2 <- nchar(as.character(semilla2))
  
  #Verificamos que las dos semillas tengan la misma cantidad de dígitos.
  #Si no tienen la misma cantidad, detenemos la función y mostramos un mensaje.
  if (d1 != d2) {
    stop("Las dos semillas deben tener la misma cantidad de dígitos.")
  }
  
  #Guardamos la cantidad de dígitos de las semillas.
  #Por ejemplo, si semilla1 = 5731 y semilla2 = 8452, entonces d = 4.
  d <- d1
  
  #Verificamos que la cantidad de dígitos sea par.
  #Esto facilita la extracción de los dígitos centrales.
  if (d %% 2 != 0) {
    stop("Las semillas deben tener una cantidad par de dígitos.")
  }
  
  #Creamos un vector llamado semillas para almacenar
  #las dos semillas iniciales y las n nuevas semillas generadas.
  #Por eso necesitamos n + 2 posiciones.
  semillas <- numeric(n + 2)
  
  #Guardamos la primera semilla inicial.
  semillas[1] <- semilla1
  
  #Guardamos la segunda semilla inicial.
  semillas[2] <- semilla2
  
  #Creamos un vector llamado numeros para almacenar
  #los números pseudoaleatorios generados.
  #Estos serán los valores U_i entre 0 y 1.
  numeros <- numeric(n)
  
  #Iniciamos un ciclo for que se repetirá n veces.
  #En cada repetición multiplicamos dos semillas consecutivas
  #para generar una nueva semilla.
  for (i in 1:n) {
    
    #Multiplicamos las dos semillas utilizadas en la iteración.
    #En la primera iteración se multiplica semilla1 por semilla2.
    producto <- semillas[i] * semillas[i + 1]
    
    #Convertimos el producto a texto.
    #sprintf() agrega ceros a la izquierda para garantizar
    #que el resultado tenga exactamente 2*d dígitos.
    #Si d = 4, el producto se representará con 8 dígitos.
    producto_texto <- sprintf(paste0("%0", 2 * d, "d"), producto)
    
    #Calculamos la posición donde comienzan los d dígitos centrales.
    #Por ejemplo, si d = 4:
    #inicio = 4/2 + 1 = 3.
    inicio <- d / 2 + 1
    
    #Calculamos la posición donde terminan los d dígitos centrales.
    #Si d = 4:
    #fin = 3 + 4 - 1 = 6.
    fin <- inicio + d - 1
    
    #Extraemos los d dígitos centrales del producto.
    #Por ejemplo, si el producto es 32844361,
    #los cuatro dígitos centrales son 8443.
    nueva_semilla <- substr(producto_texto, inicio, fin)
    
    #Convertimos la nueva semilla de texto a número.
    #Por ejemplo: "8443" se convierte en 8443.
    semillas[i + 2] <- as.numeric(nueva_semilla)
    
    #Transformamos la nueva semilla en un número pseudoaleatorio
    #entre 0 y 1.
    #Para una semilla de 4 dígitos:
    #U_i = X_i / 10^4.
    numeros[i] <- semillas[i + 2] / 10^d
  }
  
  #Creamos un data.frame para presentar los resultados
  #de una manera organizada.
  resultados <- data.frame(
    
    #Número de la iteración. Va desde 1 hasta n.
    Iteracion = 1:n,
    
    #Primera semilla utilizada en cada multiplicación.
    Semilla_1 = semillas[1:n],
    
    #Segunda semilla utilizada en cada multiplicación.
    Semilla_2 = semillas[2:(n + 1)],
    
    #Producto de las dos semillas.
    Producto = semillas[1:n] * semillas[2:(n + 1)],
    
    #Nueva semilla obtenida a partir de los dígitos centrales.
    Nueva_semilla = semillas[3:(n + 2)],
    
    #Número pseudoaleatorio generado entre 0 y 1.
    Numero = numeros
  )
  
  #La función devuelve el data.frame con todos los resultados.
  return(resultados)
}


#Usamos la función.
#Generamos 3 números pseudoaleatorios utilizando
#las semillas iniciales 3450 y 6789.
productos_medios(semilla1 = 4725, semilla2 = 5250, n = 3)

  Iteracion Semilla_1 Semilla_2 Producto Nueva_semilla Numero
1         1      4725      5250 24806250          8062 0.8062
2         2      5250      8062 42325500          3255 0.3255
3         3      8062      3255 26241810          2418 0.2418

Método de multiplicador constante

Corresponde a un método similar al generador de números mediante productos medios. Este método usa una constante multiplicativa denotada \(k\) que reemplaza una de las semillas y una semilla inicial \(x_0\). El procedimiento para hallar números aleatorios \(R_i\) mediante el presente método es el siguiente:

1. Establecer constante \(k\) y semilla \(x_0\). Tanto \(k\) como \(x_0\) deberían ser números con más de \(3\) dígitos. \(i =0,1,2...m\) corresponde a las iteraciones del método.

2. Halla el producto \(z_0\)

\[\begin{align} z_0 = k*x_0 \end{align}\]

En general, para una iteración \(i\):

\[\begin{align} z_i = k*x_i \end{align}\]

3. Hallar \(x_1\) seleccionando los dígitos centrales de \(z_0\)

\[\begin{align} x_1 = digitos(z_0) \end{align}\]

En general, para una iteración \(i\):

\[\begin{align} x_{i+1} = digitos(z_i) \end{align}\]

4. Hallar el número pseudoaleatorio \(R_{i+1}\)

\[\begin{align} R_{i+1} = \frac{x_{i+1}}{10000} \end{align}\]

El valor del denominador de la división puede cambiar en múltiplos de \(10\)

Ejemplo de generación de pseudoaleatorios usando multiplicador constante.

Generar los primeros \(3\) números pseudoaleatorios usando el método de multiplicador constante con \(k=5323\) y \(x_0 = 8506\).

Para \(i= 0\)

\[\begin{align} k&=5323\\ x_0 &= 8506 \end{align}\]

b.\(z_0\)

\[\begin{align} z_0 &= 5323*8506\\ z_0 &= 45277438 \end{align}\]

3. \(x_1\)

\[\begin{align} x_1 = 2774 \end{align}\]

4. Hallar el número pseudoaleatorio \(R_{i+1}\)

\[\begin{align} R_{1} &= \frac{2774}{10000}\\ R_1 &= 0,2774 \end{align}\]

Para \(i= 1\)

\[\begin{align} k&=5323\\ x_1 &= 2774 \end{align}\]

b.\(z_1\)

\[\begin{align} z_1 &= 5323*2774\\ z_1 &= 14766002 \end{align}\]

3. \(x_2\)

\[\begin{align} x_2 = 7660 \end{align}\]

4. Hallar el número pseudoaleatorio

\[\begin{align} R_{2} &= \frac{7660}{10000}\\ R_2 &= 0,7660 \end{align}\]

Para \(i= 2\)

\[\begin{align} k&=5323\\ x_2 &= 7660 \end{align}\]

2. \(z_2\)

\[\begin{align} z_2 &= 5323*7660\\ z_2 &= 40774180 \end{align}\]

3. \(x_3\)

\[\begin{align} x_3 = 7741 \end{align}\]

4. Hallar el número pseudoaleatorio

\[\begin{align} R_{3} &= \frac{7741}{10000}\\ R_3 &= 0,7741 \end{align}\]

Algoritmo en R

#Definimos una función llamada multiplicador_constante.
#La función recibe:
#semilla: valor inicial para comenzar la generación.
#multiplicador: valor constante que se multiplicará por cada semilla.
#n: cantidad de números pseudoaleatorios que queremos generar.

multiplicador_constante <- function(semilla, multiplicador, n) {
  
  #Convertimos la semilla a texto mediante as.character()
  #y contamos cuántos dígitos tiene con nchar().
  #Por ejemplo, si semilla = 5731, entonces d = 4.
  d <- nchar(as.character(semilla))
  
  #Verificamos que la cantidad de dígitos de la semilla sea par.
  #Esto facilita la extracción de los dígitos centrales.
  if (d %% 2 != 0) {
    stop("La semilla debe tener una cantidad par de dígitos.")
  }
  
  #Creamos un vector llamado semillas para almacenar
  #todas las semillas que se obtendrán durante el proceso.
  #Se necesitan n + 1 posiciones porque almacenaremos:
  #la semilla inicial + las n nuevas semillas.
  semillas <- numeric(n + 1)
  
  #Guardamos la semilla inicial en la primera posición
  #del vector semillas.
  semillas[1] <- semilla
  
  #Creamos un vector llamado numeros para almacenar
  #los números pseudoaleatorios generados.
  #Estos serán los valores U_i que estarán entre 0 y 1.
  numeros <- numeric(n)
  
  #Iniciamos un ciclo for que se repetirá n veces.
  #Cada repetición multiplica la semilla actual por
  #el mismo multiplicador constante.
  for (i in 1:n) {
    
    #Multiplicamos la semilla actual por el multiplicador constante.
    #Por ejemplo, si la semilla es 5731 y el multiplicador es 5123:
    #5731 * 5123 = 29359913.
    producto <- semillas[i] * multiplicador
    
    #Convertimos el resultado del producto a texto.
    #sprintf() permite agregar ceros a la izquierda para garantizar
    #que el número tenga una cantidad suficiente de dígitos.
    producto_texto <- sprintf(
      paste0("%0", 2 * d, "d"),
      producto
    )
    
    #Calculamos la posición donde comienzan los d dígitos centrales.
    #Para una semilla de 4 dígitos:
    #d/2 + 1 = 3.
    #Por tanto, comenzaremos a extraer desde el tercer dígito.
    inicio <- d / 2 + 1
    
    #Calculamos la posición donde terminan los d dígitos centrales.
    #Si d = 4:
    #inicio = 3.
    #fin = 3 + 4 - 1 = 6.
    #Por tanto, extraemos los dígitos desde la posición 3 hasta la 6.
    fin <- inicio + d - 1
    
    #Extraemos los d dígitos centrales del producto.
    #Por ejemplo, si el producto es 29359913,
    #los cuatro dígitos centrales son 3599.
    nueva_semilla <- substr(producto_texto, inicio, fin)
    
    #Convertimos la nueva semilla de texto a número.
    #Por ejemplo: "3599" se convierte en 3599.
    #Esto permite utilizarla en la siguiente iteración.
    semillas[i + 1] <- as.numeric(nueva_semilla)
    
    #Transformamos la nueva semilla en un número pseudoaleatorio
    #entre 0 y 1.
    #Para una semilla de 4 dígitos:
    #U_i = X_i / 10^4.
    #Por ejemplo:
    #3599 / 10000 = 0.3599.
    numeros[i] <- semillas[i + 1] / 10^d
  }
  
  #Creamos un data.frame para presentar los resultados
  #de una manera organizada.
  resultados <- data.frame(
    
    #Número de la iteración. Va desde 1 hasta n.
    Iteracion = 1:n,
    
    #Semilla utilizada al comienzo de cada iteración.
    Semilla = semillas[1:n],
    
    #Multiplicador constante utilizado en todas las iteraciones.
    Multiplicador = rep(multiplicador, n),
    
    #Producto de la semilla por el multiplicador constante.
    Producto = semillas[1:n] * multiplicador,
    
    #Nueva semilla obtenida a partir de los dígitos centrales.
    Nueva_semilla = semillas[2:(n + 1)],
    
    #Número pseudoaleatorio generado entre 0 y 1.
    Numero = numeros
  )
  
  #La función devuelve el data.frame con todos los resultados.
  return(resultados)
}


#Usamos la función.
#Generamos 3 números pseudoaleatorios utilizando
#la semilla inicial 3450 y el multiplicador constante 6789.
multiplicador_constante(
  semilla = 8506,
  multiplicador = 5323,
  n = 3
)

  Iteracion Semilla Multiplicador Producto Nueva_semilla Numero
1         1    8506          5323 45277438          2774 0.2774
2         2    2774          5323 14766002          7660 0.7660
3         3    7660          5323 40774180          7741 0.7741

Métodos congruenciales

Los métodos congruenciales provienen del empleo de la relación fundamental de congruencia. El propósito del uso de métodos congruenciales es la generación, en tiempos cortos, de sucesiones de números peusoaleatorios con períodos máximos.

Método congruencial lineal.

La ecuación recurrente para generar números pseudoaleatorios es la que se muestra en la Ecuación 4

\[\begin{align} x_{i+1} &= (ax_i + c)~ \text{mod}~ m\\ \\ R_i &= \frac{x_i}{m-1}\\ \\ i&= 1,2,3,...,n\\ \\ a, c, m:~~enteros \end{align} \tag{4}\]

Donde:

• \(x_0\): semilla

• \(a\): constante multiplicativa

• \(c\): constante aditiva.

• \(m\): modulador

Ejemplo de generación método congruencial lineal

Generar \(3\) números pseudoaleatorios utilizando el método congruencial lineal con:

\[\begin{align} x_0 &= 17\\ a &= 11\\ c &= 15\\ m &=100 \end{align}\]

Generando \(R_1\)

\[\begin{align} x_1 &= (ax_0 + c) \text{mod}~m\\ \\ x_1 &= (11*17 + 15) \text{mod}~100\\ \\ x_1 &= (202) \text{mod}~100\\ \\ x_1 &= 2 \\ R_1 &= \frac{2}{99} = 0.02020202 \end{align}\]

Generando \(R_2\)

\[\begin{align} x_1 &= 2\\ \\ x_2 &= (11*2 +15) \text{mod}~100\\ \\ x_2 &= (37) \text{mod}~100\\ \\ x_2 &= 37 \\ R_2 &= \frac{37}{99} = 0.3737374 \end{align}\]

Generando \(R_3\)

\[\begin{align} x_2 &= 37\\ \\ x_3 &= (11*37 +15) \text{mod}~100\\ \\ x_3 &= (422) \text{mod}~100\\ \\ x_3 &= 22 \\ R_3 &= \frac{22}{99} = 0.2222222 \end{align}\]

Algoritmo en R

#Indicamos que R muestre hasta 12 dígitos significativos
#en los resultados numéricos.
options(digits = 4)


#Definimos una función llamada congruencial_lineal.
#La función recibe:
#semilla: valor inicial para comenzar la generación.
#multiplicador: valor a utilizado en el método.
#incremento: valor c utilizado en el método.
#modulo: valor m utilizado en el método.
#n: cantidad de números pseudoaleatorios que queremos generar.

congruencial_lineal <- function(
  semilla,
  multiplicador,
  incremento,
  modulo,
  n
) {
  
  #Creamos un vector llamado semillas para almacenar
  #todas las semillas que se obtendrán durante el proceso.
  #Se necesitan n + 1 posiciones porque almacenaremos:
  #la semilla inicial + las n nuevas semillas.
  semillas <- numeric(n + 1)
  
  #Guardamos la semilla inicial en la primera posición
  #del vector semillas.
  semillas[1] <- semilla
  
  #Creamos un vector llamado numeros para almacenar
  #los números pseudoaleatorios generados.
  #Estos serán los valores U_i que estarán entre 0 y 1.
  numeros <- numeric(n)
  
  #Iniciamos un ciclo for que se repetirá n veces.
  #Cada repetición genera una nueva semilla mediante
  #la fórmula del método congruencial lineal.
  for (i in 1:n) {
    
    #Calculamos el producto entre el multiplicador y
    #la semilla actual.
    producto <- multiplicador * semillas[i]
    
    #Sumamos el incremento al producto.
    #De esta manera obtenemos:
    #a * X_i + c.
    suma <- producto + incremento
    
    #Aplicamos el operador módulo.
    #El resultado es el residuo de dividir la suma
    #entre el módulo m.
    #La nueva semilla queda entre 0 y m - 1.
    semillas[i + 1] <- suma %% modulo
    
    #Transformamos la nueva semilla en un número
    #pseudoaleatorio entre 0 y 1.
    #La fórmula es:
    #U_i = X_(i+1) / m.
    numeros[i] <- semillas[i + 1] / modulo
  }
  
  #Creamos un data.frame para presentar los resultados
  #de una manera organizada.
  resultados <- data.frame(
    
    #Número de la iteración. Va desde 1 hasta n.
    Iteracion = 1:n,
    
    #Semilla utilizada al comienzo de cada iteración.
    Semilla = semillas[1:n],
    
    #Multiplicador utilizado en el método.
    Multiplicador = rep(multiplicador, n),
    
    #Incremento utilizado en el método.
    Incremento = rep(incremento, n),
    
    #Módulo utilizado en el método.
    Modulo = rep(modulo, n),
    
    #Producto entre el multiplicador y la semilla.
    Producto = multiplicador * semillas[1:n],
    
    #Resultado de sumar el incremento.
    Suma = multiplicador * semillas[1:n] + incremento,
    
    #Nueva semilla obtenida después de aplicar el módulo.
    Nueva_semilla = semillas[2:(n + 1)],
    
    #Número pseudoaleatorio generado entre 0 y 1.
    Numero = numeros
  )
  
  #La función devuelve el data.frame con todos los resultados.
  return(resultados)
}


#Usamos la función.
#Generamos 10 números pseudoaleatorios utilizando:
#semilla inicial = 7
#multiplicador = 5
#incremento = 3
#módulo = 16
#n = 10

congruencial_lineal(
  semilla = 7,
  multiplicador = 5,
  incremento = 3,
  modulo = 16,
  n = 10
)

   Iteracion Semilla Multiplicador Incremento Modulo Producto Suma
1          1       7             5          3     16       35   38
2          2       6             5          3     16       30   33
3          3       1             5          3     16        5    8
4          4       8             5          3     16       40   43
5          5      11             5          3     16       55   58
6          6      10             5          3     16       50   53
7          7       5             5          3     16       25   28
8          8      12             5          3     16       60   63
9          9      15             5          3     16       75   78
10        10      14             5          3     16       70   73
   Nueva_semilla Numero
1              6 0.3750
2              1 0.0625
3              8 0.5000
4             11 0.6875
5             10 0.6250
6              5 0.3125
7             12 0.7500
8             15 0.9375
9             14 0.8750
10             9 0.5625

Condiciones para períodos largos con método congruencial lineal

Para obtener un período grande en la generación de números aleatorios con el método congruencial linea, se debe elegir los parámetros \(a\), \(c\) y \(m\) de manera adecuada, cumpliendo ciertas condiciones teóricas:

1. \(c\) y \(m\) son primos entre sí.

Lo anterior significa que su máximo común divisor es \(1\). Recuerde que para el cálculo de máximo común divisor entre dos números se hace:

1. Descomponer los dos números en sus factores primos.

2. Señalar factores comunes.

3. Entre los comunes escoger factores con menor exponente.

4. Multiplicar factores escogidos.

Por ejemplo, si se tienen \(c=4\) y \(m=8\) el mínimo común divisor será:

c <- 4
m <- 8

library(gmp)

Adjuntando el paquete: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

mcd<-gcd(c,m) 
mcd

[1] 4

Como el mínimo común divisor entre \(c=4\) y \(m=8\) es \(4\), esto es, \(mcd(4,8) = 4 \neq 1\), entonces el período de generación será corto.

Ahora, si se tienen \(c=3\) y \(m=8\) el mínimo común divisor será:

c <- 3
m <- 8

library(gmp)
mcd<-gcd(c,m) 
mcd

[1] 1

Como el mínimo común divisor entre \(c=3\) y \(m=8\) es \(1\), esto es, \(mcd(3,8) = 1\), entonces el generador tiene más posibilidades de producir una secuencia larga antes de repetirse.

2. \(a-1\) es múltiplo de todos los factores primos de \(m\)

Sí se descompone \(m\) en sus factores primos, entonces \(a-1\) debe ser divisible por cada uno de esos factores.

Suponga que \(m = 8\), si se descompone \(m\) es sus factores primos se obtiene

\[\begin{align} 8 &= 2*2*2 \\ 8 &= 2^3 \end{align}\]

Por lo que el único factor primo de \(8\) es \(2\), esto implica que \(a-1\) es múltiplo de \(2\). Por lo tanto \(a\) debe ser impar

Ejemplo de generación con período grande

• \(m=2^32\)

• \(a=1664525\)

• \(c=101390423\)

Condición 1: para el máximo común divisor:

c <- 101390423
m <- 2^32

library(gmp)
mcd<-gcd(c,m) 
mcd

[1] 1

Entonces \(c\) y \(m\) son primos entre sí.

Condición 2: para los factores primos de \(m\)

m <- 2^32

library(numbers)
factores_primos<-primeFactors(m)
factores_primos

 [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

El factor primo de \(m\) es \(2\), entonces \(a\) debe ser impar.

Generando números con los parámetros anteriores:

# Cantidad de números a generar.
n<- 100

#Parámetros
a <- 1664525
c <- 101390423
m <- 2^32
semilla <- 17

#Vector de almacenamiento de números
x <- numeric(n)
# Asignar semilla
x[1] <- semilla

# Ciclo de generación
for (i in 2:n){
  x[i] <- (a*x[i-1] + c) %% m
  }
R <- x/(m-1)
R

  [1] 3.958e-09 3.020e-02 6.701e-01 8.858e-01 2.391e-01 9.373e-01 1.088e-01
  [8] 8.209e-01 9.922e-02 5.009e-01 3.928e-01 1.913e-01 5.314e-01 6.324e-02
 [15] 7.484e-01 8.275e-01 4.802e-01 1.058e-01 1.282e-01 9.636e-01 5.980e-01
 [22] 3.800e-01 8.850e-01 5.907e-01 5.697e-01 4.731e-01 5.956e-01 2.545e-01
 [29] 6.886e-01 1.368e-01 2.986e-01 7.466e-01 5.737e-01 9.563e-01 6.636e-01
 [36] 3.962e-01 1.692e-01 1.885e-01 8.567e-01 9.121e-01 4.170e-01 8.473e-01
 [43] 1.173e-01 9.803e-01 2.427e-01 3.642e-01 3.131e-01 6.957e-02 7.454e-01
 [50] 4.016e-01 9.883e-01 1.481e-01 7.325e-01 6.489e-01 6.694e-01 2.643e-01
 [57] 9.897e-01 6.763e-01 4.475e-01 1.949e-01 9.273e-01 1.781e-01 9.840e-02
 [64] 6.141e-01 4.928e-01 1.804e-01 2.602e-01 4.855e-01 4.905e-01 4.392e-02
 [71] 4.308e-01 6.316e-01 6.995e-01 3.379e-01 8.034e-01 6.827e-01 4.299e-01
 [78] 9.147e-01 2.038e-01 8.590e-01 4.080e-01 7.628e-01 6.462e-01 2.888e-01
 [85] 1.313e-01 4.661e-01 2.727e-01 8.156e-01 8.439e-02 3.476e-01 6.922e-01
 [92] 4.984e-01 4.545e-01 6.259e-01 9.845e-01 2.569e-01 2.401e-01 1.630e-01
 [99] 8.846e-01 3.651e-01

Recomendaciones:

Para que el generador provea un período amplio:

• \(m=2^h\)
• \(a = 4k+1\)

Con \(k\) entero, \(c\) y \(m\) primos entre sí.

Método congruencial multiplicativo

La ecuación recurrente para generar números pseudoaleatorios es la que se muestra en la Ecuación 5

\[\begin{align} x_{i+1} &= (ax_i)~ \text{mod}~ m\\ \\ R_i &= \frac{x_i}{m-1}\\ \\ i &= 1,2,3,...,n\\ \\ a &= 3+8k\\ \\ m &= 2^h\\ \\ x_0 &~ impar \end{align} \tag{5}\]

Ejemplo de generación método congruencial multiplicativo

Generar \(3\) números pseudoaleatorios utilizando el método congruencial multiplicativo con:

\[\begin{align} x_0 &= 3\\ k &= 2\\ g &= 5 \end{align}\]

Con lo anterior:

\[\begin{align} a &= 3+ 8*2 = 19 \end{align}\]

\[\begin{align} m = 2^5 = 32 \end{align}\]

Generando \(R_1\)

\[\begin{align} x_1 &= (ax_0) \text{mod}~m\\ \\ x_1 &= (19*3) \text{mod}~32\\ \\ x_1 &= (57) \text{mod}~32\\ \\ x_1 &= 25 \\ R_1 &= \frac{25}{31} = 0.8064516 \end{align}\]

Generando \(R_2\)

\[\begin{align} x_1 &= 25\\ \\ x_2 &= (19*25) \text{mod}~32\\ \\ x_2 &= (475) \text{mod}~32\\ \\ x_2 &= 27 \\ R_2 &= \frac{27}{31} = 0.8709677 \end{align}\]

Generando \(R_3\)

\[\begin{align} x_2 &= 27\\ \\ x_3 &= (19*27) \text{mod}~32\\ \\ x_3 &= (513) \text{mod}~32\\ \\ x_3 &= 1 \\ R_3 &= \frac{1}{31} = 0.003225806 \end{align}\]

Algoritmo en R

options(digits = 12)

#Definimos una función llamada congruencial_multiplicativo.
#La función recibe:
#semilla: valor inicial para comenzar la generación.
#multiplicador: valor a utilizado en el método.
#modulo: valor m utilizado en el método.
#n: cantidad de números pseudoaleatorios que queremos generar.

congruencial_multiplicativo <- function(
  semilla,
  multiplicador,
  modulo,
  n
) {
  
  #Creamos un vector llamado semillas para almacenar
  #todas las semillas que se obtendrán durante el proceso.
  #Se necesitan n + 1 posiciones porque almacenaremos:
  #la semilla inicial + las n nuevas semillas.
  semillas <- numeric(n + 1)
  
  #Guardamos la semilla inicial en la primera posición
  #del vector semillas.
  semillas[1] <- semilla
  
  #Creamos un vector llamado numeros para almacenar
  #los números pseudoaleatorios generados.
  #Estos serán los valores U_i que estarán entre 0 y 1.
  numeros <- numeric(n)
  
  #Iniciamos un ciclo for que se repetirá n veces.
  #Cada repetición genera una nueva semilla mediante
  #la fórmula del método congruencial multiplicativo.
  for (i in 1:n) {
    
    #Multiplicamos la semilla actual por el multiplicador.
    #De esta manera obtenemos:
    #a * X_i.
    producto <- multiplicador * semillas[i]
    
    #Aplicamos el operador módulo.
    #El resultado es el residuo de dividir el producto
    #entre el módulo m.
    #La nueva semilla queda entre 0 y m - 1.
    semillas[i + 1] <- producto %% modulo
    
    #Transformamos la nueva semilla en un número
    #pseudoaleatorio entre 0 y 1.
    #La fórmula es:
    #U_i = X_(i+1) / m.
    numeros[i] <- semillas[i + 1] / (modulo-1)
  }
  
  #Creamos un data.frame para presentar los resultados
  #de una manera organizada.
  resultados <- data.frame(
    
    #Número de la iteración. Va desde 1 hasta n.
    Iteracion = 1:n,
    
    #Semilla utilizada al comienzo de cada iteración.
    Semilla = semillas[1:n],
    
    #Multiplicador utilizado en el método.
    Multiplicador = rep(multiplicador, n),
    
    #Módulo utilizado en el método.
    Modulo = rep(modulo, n),
    
    #Producto entre el multiplicador y la semilla.
    Producto = multiplicador * semillas[1:n],
    
    #Nueva semilla obtenida después de aplicar el módulo.
    Nueva_semilla = semillas[2:(n + 1)],
    
    #Número pseudoaleatorio generado entre 0 y 1.
    Numero = numeros
  )
  
  #La función devuelve el data.frame con todos los resultados.
  return(resultados)
}


#Usamos la función.
#Generamos 10 números pseudoaleatorios utilizando:
#semilla inicial = 12345
#multiplicador = 1103515245
#módulo = 2^31
#n = 10

congruencial_multiplicativo(
  semilla = 3,
  multiplicador = 19,
  modulo = 32,
  n = 3
)

  Iteracion Semilla Multiplicador Modulo Producto Nueva_semilla          Numero
1         1       3            19     32       57            25 0.8064516129032
2         2      25            19     32      475            27 0.8709677419355
3         3      27            19     32      513             1 0.0322580645161

A partir de número aleatorio \(R_8\) el ciclo de repite, el algoritmo cumplió su ciclo de vida. En este método el ciclo de vida está dado por:

\[\begin{align} 2^{g-2} \end{align}\]