1.3 Números aleatorios en R

La generación de números pseudoaleatorios en R es una de las mejores disponibles en paquetes estadísticos. Entre las herramientas implementadas en el paquete base de R podemos destacar:

  • set.seed(entero): permite establecer la semilla (y el generador).

  • RNGkind(): selecciona el generador.

  • rdistribución(n,...): genera valores aleatorios de la correspondiente distribución. Por ejemplo, runif(n, min = 0, max = 1), generaría n valores de una uniforme. Se puede acceder al listado completo de las funciones disponibles en el paquete stats mediante el comando ?distributions.

  • sample(): genera muestras aleatorias de variables discretas y permutaciones (se tratará en el Capítulo 5).

  • simulate(): genera realizaciones de la respuesta de un modelo ajustado.

Además están disponibles otros paquetes que implementan distribuciones adicionales (ver CRAN Task View: Probability Distributions). Entre ellos podríamos destacar los paquetes distr (clases S4; con extensiones en otros paquetes) y distr6 (clases R6).

La semilla se almacena en .Random.seed:

  • Inicialmente no existe. La recomendación es establecerla con set.seed(), en caso contrario se generará a partir del reloj del sistema4 cuando se necesite.

  • Se almacena como un objeto oculto en el entorno de trabajo (o entorno global .GlobalEnv). Con las opciones por defecto de R, si al terminar una sesión almacenamos el entorno (en un fichero .RData), al iniciar una nueva sesión se restaurará también la semilla (y se podría continuar con las simulaciones).

  • Es un vector de enteros cuya estructura depende del tipo de generador (en la Sección 2.2 se dan algunos detalles sobre la configuración por defecto), por lo que no debería ser modificado manualmente.
    Puede ser recomendable almacenar (el objeto completo) antes de generar simulaciones, e.g. seed <- .Random.seed. Esto permite reproducir los resultados y facilita la depuración de posibles errores.

En la mayoría de los ejemplos de este libro se generan todos los valores de una vez, se guardan y se procesan vectorialmente (normalmente empleando la función apply). En problemas mas complejos, en los que no es necesario almacenar todas las simulaciones, puede ser preferible emplear un bucle para generar y procesar cada simulación iterativamente. Por ejemplo podríamos emplear el siguiente esquema:

# Fijar semilla
set.seed(1)
for (isim in 1:nsim) {
  seed <- .Random.seed
  # Si se produce un error, podremos depurarlo ejecutando:
  #  .Random.seed <- seed
  ...
  # Generar valores pseudoaleatorios
  ...
}

o alternativamente fijar la semilla en cada iteración, por ejemplo:

for (isim in 1:nsim) {
  set.seed(isim)
  ...
  # Generar valores pseudoaleatorios
  ...
}

1.3.1 Opciones

Normalmente no nos va a interesar cambiar las opciones por defecto de R para la generación de números pseudoaleatorios. Para establecer estas opciones podemos emplear los argumentos kind = NULL, normal.kind = NULL y sample.kind = NULL en las funciones RNGkind() o set.seed(). A continuación se muestran las distintas opciones (resaltando en negrita los valores por defecto):

  • kind especifica el generador pseudoaleatorio (uniforme):

    • “Wichmann-Hill”: Ciclo \(6.9536\times10^{12}\)

    • “Marsaglia-Multicarry”: Ciclo mayor de \(2^{60}\)

    • “Super-Duper”: Ciclo aprox. \(4.6\times10^{18}\) (S-PLUS)

    • “Mersenne-Twister”: Ciclo \(2^{19937}-1\) y equidistribution en 623 dimensiones.

    • “Knuth-TAOCP-2002”: Ciclo aprox. \(2^{129}\).

    • “Knuth-TAOCP”

    • “user-supplied”: permite emplear generadores adicionales.

  • normal.kind selecciona el método de generación de normales (se tratará más adelante): “Kinderman-Ramage”, “Buggy Kinderman-Ramage”, “Ahrens-Dieter”, “Box-Muller”, “Inversion” , o “user-supplied”.

  • sample.kind selecciona el método de generación de uniformes discretas (el empleado por la función sample(), que cambió ligeramente5 a partir de la versión 3.6.0 de R): “Rounding” (versión anterior a 3.6.0) o “Rejection”.

Estas opciones están codificadas (con índices comenzando en 0) en el primer componente de la semilla:

set.seed(1)
.Random.seed[1]
## [1] 10403

Los dos últimos dígitos se corresponden con el generador, las centenas con el método de generación de normales y las decenas de millar con el método uniforme discreto.

1.3.2 Paquetes de R

Otros paquetes de R que pueden ser de interés:

  • setRNG contiene herramientas que facilitan operar con la semilla (dentro de funciones,…).

  • random permite la descarga de números “true random” desde RANDOM.ORG.

  • randtoolbox implementa generadores más recientes (rngWELL) y generación de secuencias cuasi-aleatorias.

  • RDieHarder implementa diversos contrastes para el análisis de la calidad de un generador y varios generadores.

  • Runuran interfaz para la librería UNU.RAN para la generación (automática) de variables aleatorias no uniformes (ver Hörmann et al., 2004).

  • rsprng, rstream y rlecuyer implementan la generación de múltiples secuencias (para programación paralela).

  • gls, rngwell19937, randaes, SuppDists, lhs, mc2d, fOptions, …

1.3.3 Tiempo de CPU

La velocidad del generador suele ser una característica importante (también medir los tiempos, de cada iteración y de cada procedimento, en estudios de simulación). Para evaluar el rendimiento están disponibles en R distintas herramientas:

  • proc.time(): permite obtener tiempo de computación real y de CPU.

    tini <- proc.time()
# Código a evaluar
tiempo <- proc.time() - tini

  • system.time(expresión): muestra el tiempo de computación (real y de CPU) de expresión.

Por ejemplo, podríamos emplear las siguientes funciones para ir midiendo los tiempos de CPU durante una simulación:

CPUtimeini <- function() {
  .tiempo.ini <<- proc.time()
  .tiempo.last <<- .tiempo.ini
}

CPUtimeprint <- function() {
  tmp <- proc.time()
  cat("Tiempo última operación:\n")
  print(tmp-.tiempo.last)
  cat("Tiempo total operación:\n")
  print(tmp-.tiempo.ini)
  .tiempo.last <<- tmp
}

Llamando a CPUtimeini() donde se quiere empezar a contar, y a CPUtimeprint() para imprimir el tiempo total y el tiempo desde la última llamada a una de estas funciones. Ejemplo:

funtest <- function(n) mad(runif(n)) 
CPUtimeini()
result1 <- funtest(10^6)
CPUtimeprint()
## Tiempo última operación:
##    user  system elapsed 
##    0.16    0.00    0.16 
## Tiempo total operación:
##    user  system elapsed 
##    0.16    0.00    0.16
result2 <- funtest(10^3)
CPUtimeprint()
## Tiempo última operación:
##    user  system elapsed 
##    0.01    0.00    0.01 
## Tiempo total operación:
##    user  system elapsed 
##    0.17    0.00    0.17

Hay diversos paquetes que implementan herramientas similares, por ejemplo:

  • El paquete tictoc:

    • tic("mensaje"): inicia el temporizador y almacena el tiempo de inicio junto con el mensaje en una pila.

    • toc(): calcula el tiempo transcurrido desde la llamada correspondiente a tic().

      library(tictoc)
## Timing nested code
tic("outer")
   result1 <- funtest(10^6)
   tic("middle")
      result2 <- funtest(10^3)
      tic("inner")
         result3 <- funtest(10^2)
      toc() # inner
      ## inner: 0 sec elapsed
         toc() # middle
      ## middle: 0 sec elapsed
      toc() # outer
      ## outer: 0.14 sec elapsed
      ## Timing in a loop and analyzing the results later using tic.log().
tic.clearlog()
for (i in 1:10)
{
   tic(i)
   result <- funtest(10^4)
   toc(log = TRUE, quiet = TRUE)
}
# log.txt <- tic.log(format = TRUE)
# log.lst <- tic.log(format = FALSE)
log.times <- do.call(rbind.data.frame, tic.log(format = FALSE))
str(log.times)
      ## 'data.frame':    10 obs. of  4 variables:
##  $ tic         : num  4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26
##  $ toc         : num  4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26 4.26
##  $ msg         : chr  "1" "2" "3" "4" ...
##  $ callback_msg: chr  "1: 0 sec elapsed" "2: 0 sec elapsed" "3: 0 sec elapsed"..
      tic.clearlog()

# timings <- unlist(lapply(log.lst, function(x) x$toc - x$tic))
log.times$timings <- with(log.times, toc - tic)
summary(log.times$timings)
      ##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

  • La función cpu.time() del paquete simres:

    • cpu.time(restart = TRUE): inicia el temporizador y almacena el tiempo de inicio.

    • cpu.time(): calcula el tiempo (real y de CPU) total (desde tiempo de inicio) y parcial (desde la última llamada a esta función).

Hay que tener en cuenta que, por construcción, aunque se realicen en la mismas condiciones (en el mismo equipo), los tiempos de CPU en R pueden variar “ligeramente” entre ejecuciones. Si se quieren estudiar tiempos de computación de forma más precisa, se recomendaría promediar los tiempos de varias ejecuciones. Para ello se pueden emplear las herramientas del paquete microbenchmark. No obstante, para los fines de este libro no será necesaria tanta precisión.

Finalmente, si los tiempos de computación no fuesen asumibles, para identificar los cuellos de botella y mejorar el código para optimizar la velocidad, podríamos emplear la función Rprof(fichero). Esta función permite evaluar el rendimiento muestreando la pila en intervalos para determinar en que funciones se emplea el tiempo de computación. Después de ejecutar el código, llamando a Rprof(NULL) se desactiva el muestreo y con summaryRprof(fichero) se muestran los resultados (para analizarlos puede resultar de utilidad el paquete proftools).

  1. y del identificador asignado por el sistema operativo al proceso.↩︎

  2. Para evitar problemas de redondeo con tamaños extremadamente grandes; ver bug PR#17494.↩︎