Capítulo 1

Introducción a la Simulación

La Simulación se entiende como la reproducción de un fenómeno real mediante otro más sencillo y más adecuado para ser estudiado. En el caso de los fenómenos estadísticos la consideramos como una técnica que consiste en realizar experimentos de muestreo sobre el modelo de un sistema bajo determinadas condiciones.

1.1 Conceptos básicos.

Al utilizar la experimentación directa nos enfrentamos con varias problemáticas, como: el aumento del costo, cuando depende de muchas variables el proceso se relentiza y en muchas ocasiones se vuelve incierto al depender del futuro o del pasado, o en otros casos llega ésta a no ser ética, como la experimentación en seres humanos; entre muchas otras que se pueden presentar.

Cuando la realidad es demasiado compleja el estudio se vuelve muy tedioso, por ello es preferible trabajar con un modelo del sistema real. Se entiende por modelo como un conjunto de variables junto con ecuaciones matemáticas que las relacionan y restricciones sobre dichas variables. Al tratar el grado de certeza con el cuál se conocen los parámetros de un modelo matemático, se reconocen dos tipos de modelos:

  • Modelos deterministas: Se dice determinista cuando se tiene certeza de los valores de los parámetros del modelo.

  • Modelos estocásticos(con componente aleatoria): Se entiende como estocástico o probabilístico cuando los parámetros usados para caracterizar el modelo son variables aleatorias que tienen unos comportamientos estimados pero no se conoce con certidumbre previamente cuál será el valor que tomen, debido a que no se dispone de toda la información sobre las variables que influyen en el modelo (esto puede ocurrir simplemente por errores de medida).

La modelización es una etapa presente en la mayor parte de los trabajos de investigación (especialmente en las ciencias experimentales). El modelo debería considerar las variables más relevantes para explicar el fenómeno en estudio y las principales relaciones entre ellas. La inferencia estadística proporciona herramientas para estimar los parámetros y contrastar la validez de un modelo estocástico a partir de los datos observados.

La idea es emplear el modelo (suponiendo que es válido) para resolver el problema de interés. Si se puede obtener la solución de forma analítica, esta suele ser exacta (aunque en ocasiones solo se dispone de soluciones aproximadas, basadas en resultados asintóticos, o que dependen de suposiciones que pueden ser cuestionables) y a menudo la resolución también es rápida. Cuando la solución no se puede obtener de modo analítico (o si la aproximación disponible no es adecuada) se puede recurrir a la simulación.

En este documento nos centraremos en el caso de la Simulación Estocástica, en ella: las conclusiones se obtienen generando repetidamente simulaciones del modelo aleatorio. Muchas veces se emplea la denominación de Método Monte-Carlo 1 Ref.1 como sinónimo de simulación estocástica, pero normalmente se trata de métodos especializados que emplean simulación para resolver problemas que pueden no estar relacionados con un modelo estocástico de un sistema real. Un ejemplo de estos métodos especializados son los métodos de integración y optimización Monte-Carlo.

1.1.1 Ventajas e inconvenientes de la simulación.

Cuando se trata de Ventajas, se presentan algunas como:

• Cuando la resolución analítica no puede llevarse a cabo.

• Cuando existen medios de resolver analíticamente el problema pero dicha resolución es complicada y costosa (o solo proporciona una solución aproximada).

• Si se desea experimentar antes de que exista el sistema (pruebas para la construcción de un sistema).

• Cuando es imposible experimentar sobre el sistema real por ser dicha experimentación destructiva.

• En ocasiones en las que la experimentación sobre el sistema es posible pero no ética.

• En sistemas que evolucionan muy lentamente en el tiempo.

El principal incoveniente puede ser el tiempo de computación necesario, aunque gracias a la gran potencia de cálculo de los computadores actuales, se puede obtener rápidamente una solución aproximada en la mayor parte de los problemas susceptibles de ser modelizados.

Además siempre están presentes los posibles problemas debidos a emplear un modelo:

• La construcción de un buen modelo puede ser una tarea muy costosa (compleja, laboriosa y requerir mucho tiempo; un ejemplo de ellos son los modelos climáticos).

• Frecuentemente el modelo omite variables o relaciones importantes entre ellas (los resultados pueden no ser válidos para el sistema real).

• Resulta difícil conocer la precisión del modelo formulado.

Otro problema de la simulación es que se obtienen resultados para unos valores concretos de los parámetros del modelo, por lo que en principio resultaría complicado extrapolar las conclusiones a otras situaciones.

1.2 Generación de números (pseudo)aleatorios.

Los números pseudo-aleatorios son un números generados en un proceso que parece producir números al azar, pero no lo hace realmente; pues secuencias de números pseudo-aleatorios son generadas por un algoritmo completamente determinista que con frecuencia se realiza desde por un programa de ordenador. Ref.2

Aplicaciones:

• Estadística:

– Muestreo, remuestreo, …

– Aproximación de distribuciones (de estadísticos, estimadores, …)

– Realización de contrastes, intervalos de confianza, …

– Comparación de estimadores, contrastes, …

– Validación teoría (distribución asintótica,…)

– Inferencia Bayesiana

• Optimización: Algoritmos genéticos, …

• Análisis numérico: Aproximación de integrales, resolución de ecuaciones, …

• Computación: Diseño, verificación y validación de algoritmos,…

• Criptografía: Protocolos de comunicación segura, …

• Física: Simulación de fenómenos naturales, …

• …

En el documento más adelante nos centraremos en algunas de las aplicaciones de utilidad en Estadística.

1.3 Números aleatorios puros.

Una sucesión de números aleatorios puros (true random), se caracteriza por que no existe ninguna regla o plan que nos permita conocer sus valores.

Normalmente son obtenidos por procesos físicos (loterías, ruletas, ruidos,…) y se almacenaban en tablas de dígitos aleatorios. Por ejemplo, en 1955 la Corporación RAND publicó el libro A Million Random Digits with 100,000 Normal Deviates que contenía números aleatorios generados mediante una ruleta electrónica conectada a una computadora (ver Figura 1.1).

\[ \text {Figura 1.1: Líneas 10580-10594, columnas 21-40, del libro "A Million Random Digits with 100,000 Normal Deviates".}\]

El procedimiento para generar de forma manual números aleatorios en un rango de 1 a m era el siguiente:

• Se selecciona al azar un punto de inicio en la tabla y la dirección que se seguirá.

• Se agrupan los dígitos de forma que “cubran” el valor de m.

• Se va avanzado en la dirección elegida, seleccionando los valores menores o iguales que m y descartando el resto.

Hoy en día están disponibles generadores de números aleatorios “online”, por ejemplo:

• RANDOM.ORG: ruido atmosférico (ver paquete en R).

• HotBits: desintegración radiactiva.

Aunque para un uso profesional puede ser recomendable emplear generadores implementados mediante hardware:

• Intel Digital Random Number Generator

• An Overview of Hardware based True Random Number Generators

Sus principales aplicaciones hoy en día son en criptografía y juegos de azar, donde resulta especialmente importante su impredecibilidad.

El principal inconveniente para su aplicación en el campo de la Estadística (y en otros casos) es que los valores generados deberían ser independientes e idénticamente distribuidos con distribución conocida, algo que resulta difícil (o imposible) de garantizar. Siempre está presente la posible aparición de sesgos, principalmente debidos a fallos del sistema o interferencias. Por ejemplo, en el caso de la máquina RAND, fallos mecánicos en el sistema de grabación de los datos causaron problemas de aleatoriedad.

Otro inconveniente estaría relacionado con su reproducibilidad, por lo que habría que almacenarlos en tablas si se quieren volver a reproducir los resultados.

1.3.1 Alternativas:

A partir de la década de 1960, al disponer de computadoras de mayor velocidad, empezó a resultar más eficiente generar valores mediante software en lugar de leerlos de las tablas. Se distingue entre dos tipos de secuencias:

• números Pseudo-aleatorios: simulan realizaciones de una variable aleatoria (uniforme).

• números cuasi-aleatorios: secuencias determinísticas con una distribución más uniforme en el rango considerado (se podría pensar que son una única generación de una variable aleatoria).

Algunos problemas, como la integración numérica (más adelante se tratarán métodos de integración Monte Carlo), no dependen realmente de la aleatoriedad de la secuencia. Para evitar generaciones poco probables, se puede recurrir a secuencias cuasi-aleatorias, también denominadas sucesiones de baja discrepancia (hablaríamos entonces de métodos cuasi-Monte Carlo). La idea sería que la proporción de valores en una región cualquiera sea siempre aproximadamente proporcional a la medida de la región (como sucedería en media con la distribución uniforme, aunque no necesariamente para una realización concreta).

Por ejemplo, el paquete implementa métodos para la generación de secuencias cuasialeatorias (ver Figura 1.2).

library(randtoolbox)
n <- 2000
par.old <- par( mfrow=c(1,3))
plot(halton(n, dim = 2), xlab = 'x1', ylab = 'x2')
plot(sobol(n, dim = 2), xlab = 'x1', ylab = 'x2')
plot(torus(n, dim = 2), xlab = 'x1', ylab = 'x2')

par(par.old)

Figura 1.2: Secuencias cuasi-aleatorias bidimensionales obtenidas con los métodos de Halton (izquierda), Sobol (centro) y Torus (derecha).

En este documento sólo consideraremos los números pseudoaleatorios y por comodidad se eliminará el prefijo “pseudo” en algunos casos.

1.4 Números pseudoaleatorios.

La mayoría de los métodos de simulación se basan en la posibilidad de generar números pseudoaleatorios que imiten las propiedades de generaciones independientes de una distribución \(\mathcal{U}(0, 1)\).

El procedimiento habitual para obtiener estas secuencias es emplear un algoritmo recursivo denominado \(generador\):

\[x_i=f(x_{i-1},x_{i-2},\cdots,x_{i-k})\]

donde:

  • \(k\) es el orden del generador
  • \((x_{0},x_1,\cdots,x_{k-1})\) es la semilla (estado inicial).

El \(periodo\) o \(logitud\) del ciclo es la longitud de la secuencia antes de que vuelva a repetirse. Lo denotaremos por \(p\).

Los números de la sucesión serán predecibles, conociendo el algoritmo y la semilla.

  • Sin embargo, si no se conociesen, no se debería poder distinguir una serie de números pseudoaleatorios de una sucesión de números verdaderamente aleatoria (utilizando recursos computacionales razonables).

  • En caso contrario esta predecibilidad puede dar lugar a serios problemas. Ref.3

Como regla general, por lo menos mientras se está desarrollando un programa, interesa fijar la semilla de aleatorización.

  • Permite la reproductibilidad de los resultados.
  • Facilita la depuración del código.

Todo generador de números pseudoaleatorios mínimamente aceptable debe comportarse como si se tratase de una muestra genuina de datos independientes de una \(\mathcal{U}(0, 1)\). Otras propiedades de interés serían:

  • Reproducibilidad a partir de la semilla.

  • Periodo suficientemente largo.

  • Eficiencia (rapidez y requerimientos de memoria).

  • Portabilidad.

  • Generación de sub-secuencias (computación en paralelo).

  • Parsimonia, …

    “… random numbers should not be generated with a method chosen at random.”

    • Knuth, D.E. (TAOCP, 2002)

Gran cantidad de algoritmos:

  • Cuadrados medios, Lehmer,…
  • Congruenciales
  • Registros desfasados
  • Combinaciones
  • …

Código fuente disponible en múltiples librerias:

Nos centraremos en los generadores congruenciales, descritos en la Sección 3.1.

Capítulo 2

Números aleatorios en R

La generación de números pseudoaleatorios en R es uno de los mejores paquetes estadísticos disponibles. Entre las herramientas que están en el paquete R se encuentran:

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

  • RNGkind(): selecciona el generador.

  • rdistribución(n,…): genera valores aleatorios de la correspondiente distrubución. Por ejemplo: , generaría \(n\) valores distribuidos uniformemente.

  • sample(): genera muestras aleatorias de variables discretas y permutaciones.(ver cap. 6)

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

Además estám disponibles otros paquetes, como el paquete distr, que implementa distribuciones adicionales.

La semilla se almacena (en globalenv) en .Random.seed; es un vector de enteros cuya dimensión depende del tipo de generador:

  • No debe ser modificado manualmente; se guarda con el entorno de trabajo.

  • Si no se especifica con (o no existe) se genera a partir del reloj del sistema.

Nota: Puede ser recomendable (para depurar) almacenar antes de generar simulaciones, e.g. semilla <- .Random.seed.

En la mayoría de los ejemplos de este documento 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
  ...
}

Hay que tener en cuenta que .Random.seed se almacena en el entorno de trabajo.

2.1 Opciones.

La función permite seleccionar el tipo de generador (en negrita los valores por defecto):

  • especifica el generador aleatorio:

    – “Wichmann-Hill”: Ciclo \(6.9536 × 10^{12}\)

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

    – “Super-Duper”: Ciclo aprox. \(4.6 × 10^{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”

  • 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”.

2.2 Paquetes de R.

Otros paquetes de R que pueden ser de interes:

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

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

  • implementa generadores más recientes () y generación de secuencias cuasialeatorias.

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

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

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

  • , , , , , , ,..

2.3 Ejercicios

Ejercicio 2.1

Sea \((X,𝑌)\) es un vector aleatorio con distribución uniforme en el cuadrado \([−1, 1] \times [−1, 1]\) de área 4.

a) Aproximar mediante simulación \(𝑃(𝑋 + 𝑌 \leq 0)\) y compararla con la probabilidad teórica (obtenida aplicando la regla de Laplace \(\frac{area \ favorable}{area \ posible}\) ).

Generamos valores del proceso bidimensional:

set.seed(1)
nsim <- 10000
x <- runif(nsim, -1, 1)
y <- runif(nsim, -1, 1)

La probabilidad teórica es \(1/2\) y la aproximación por simulación es la frecuencia relativa del suceso en los valores generados:

indice <- (x+y < 0)
sum(indice)/nsim
## [1] 0.4996

Alternativamente (la frecuencia relativa es un caso particular de la media) se puede obtener de forma más simple como:

mean(indice)
## [1] 0.4996

Nota: R maneja internamente los valores lógicos como y .

b) Aproximar el valor de \(𝜋\) mediante simulación a partir de \(𝑃 (X^2 + Y^2 \leq 1)\).

set.seed(1)
n <- 10000
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
indice <- (x^2+y^2 < 1)
mean(indice)
## [1] 0.7806
pi/4
## [1] 0.7853982
pi_aprox <- 4*mean(indice)
pi_aprox
## [1] 3.1224

Generamos el correspondiente gráfico (ver Figura 2.1) (los puntos con color negro tienen distribución uniforme en el círculo unidad).

# Colores y símbolos dependiendo de si el índice correspondiente es verdadero:
color <- ifelse(indice, "black", "red")
simbolo <- ifelse(indice, 1, 4)
plot(x, y, pch = simbolo, col = color, xlim = c(-1, 1), ylim = c(-1, 1), xlab="X", ylab="Y", asp = 1)
# asp = 1 para dibujar circulo
symbols(0, 0, circles = 1, inches = FALSE, add = TRUE)
symbols(0, 0, squares = 2, inches = FALSE, add = TRUE)

Figura 2.1: Valores generados con distribución uniforme bidimensional, con colores y símbolos indicando si están dentro del círculo unidad.

Ejercicio 2.2 (Experimento de Bernoulli).

Consideramos el experimento de Bernoulli consistente en el lanzamiento de una moneda.

a) Empleando la función , obtener \(1000\) simulaciones del lanzamiento de una moneda , suponiendo que no está trucada. Aproximar la probabilidad de cara a partir de las simulaciones.

set.seed(1)
nsim <- 10000
x <- sample(c(cara = 1, cruz = 0), nsim, replace = TRUE, prob = c(0.5,0.5))
mean(x)
## [1] 0.4953
barplot(100*table(x)/nsim, ylab = "Porcentaje") # Representar porcentajes

b) En R pueden generarse valores de la distribución de Bernoulli mediante la función . Generar un gráfico de lineas considerando en el eje \(𝑋\) el número de lanzamientos (de 1 a 10000) y en el eje \(𝑌\) la frecuencia relativa del suceso cara (puede ser recomendable emplear la función ).

set.seed(1)
nsim <- 1000
p <- 0.4
x <- rbinom(nsim, size = 1, prob = p) # Simulamos una Bernouilli
n <- 1:nsim
# Alternativa programación: x <- runif(nsim) < p
mean(x)
## [1] 0.394
plot(n, cumsum(x)/n, type="l", ylab="Proporción de caras",
xlab="Número de lanzamientos", ylim=c(0,1))
abline(h=p, lty=2, col="red")

Figura 2.2: Frecuencias relativas de los valores generados con distribución Bernoulli (aproximaciones por simulación de las probabilidades teóricas).

\[ \text {Figura 2.3: Gráfico de convergencia de la aproximación por simulación a la probabilidad teórica.}\]

Ejercicio 2.3 (Simulación de un circuito).

Simular el paso de corriente a través del siguiente circuito, donde figuran las probabilidades de que pase corriente por cada uno de los interruptores:

Considerar que cada interruptor es una v.a. de Bernoulli independiente para simular 1000 valores de cada una de ellas.

Nota: R maneja internamente los valores lógicos como y . Recíprocamente, cualquier nº puede ser tratado como lógico (al estilo de C). El entero \(0\) es equivalente a y cualquier entero distinto de \(0\) a .

set.seed(1)
nsim <- 10000
x1 <- rbinom(nsim, size=1, prob=0.9)
x2 <- rbinom(nsim, size=1, prob=0.8)
z1 <- x1 | x2 # Operador lógico "O"
x3 <- rbinom(nsim, size=1, prob=0.6)
x4 <- rbinom(nsim, size=1, prob=0.5)
z2 <- x3 | x4
z3 <- z1 | z2
x5 <- rbinom(nsim, size=1, prob=0.7)
fin <- z3 & x5 # Operador lógico "Y"
mean(fin)
## [1] 0.6918

2.4 Tiempo de CPU.

La velocidad del generador suele ser una característica importante. 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.

• Rprof(fichero): permite evaluar el rendimiento muestreando la pila en intervalos para determinar en que funciones se emplea el tiempo de computación (de utilidad para optimizar la velocidad).

• Rprof(NULL): desactiva el muestreo.

• summaryRprof(fichero): muestra los resultados.

2.4.1 Utilidades tiempo de CPU

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("\nTiempo ú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) median(runif(n)) # Función de prueba... 
CPUtimeini() 
funtest(1000000) 
## [1] 0.5003826
CPUtimeprint() 
## 
## Tiempo última operación:
##    user  system elapsed 
##   0.050   0.004   0.055 
## Tiempo total operación:
##    user  system elapsed 
##   0.050   0.004   0.055
funtest(1000) 
## [1] 0.4809216
CPUtimeprint() 
## 
## Tiempo última operación:
##    user  system elapsed 
##   0.022   0.000   0.022 
## Tiempo total operación:
##    user  system elapsed 
##   0.072   0.004   0.077

2.4.2 Paquetes de R

Por ejemplo, se puede emplear el paquete npsp (fichero \(cpu.time.R\)):

• Se usa Call cpu.time(restart = TRUE) donde quieres empezar a contar.

• Se usa Call cpu.time() para imprimir / obtener tiempos totales y / o parciales (desde la última llamada a esta función) reales y de CPU.

# CPU utilidades de tiempo 
#-------------------------
.cpu.time.ini <- function() { 
time.ini <- structure(rep(0, 5), .Names = c("user.self", "sys.self", "elapsed", 
"user.child", "sys.child"), class = "proc_time")# proc.time() 
time.last <- time.ini 
function(..., reset = FALSE, total = TRUE, last = TRUE, flush = FALSE) { 
res <- list(time = proc.time()) 
if (reset) { 
time.ini <<- res$time 
time.last <<- time.ini 
res$last <- res$total <- 0 
if (total | last) cat("CPU time has been initialized.\n") 
} else { 
res$last <- res$time - time.last 
res$total <- res$time - time.ini 
if (last) { 
cat("Time of last operation:", ..., "\n") 
print(res$last) 
} 
if (total) { 
cat("Total time:\n") 
print(res$total) 
} 
if (flush) flush.console() 
time.last <<- res$time 
} 
return(invisible(res)) 
} 
} 

cpu.time <- .cpu.time.ini()

Ejemplo

cpu.time(reset = TRUE) 
## CPU time has been initialized.
res <- funtest(1000000) 
cpu.time('\nSample median of', 1000000, 'values =', res, total = FALSE) 
## Time of last operation: 
## Sample median of 1e+06 values = 0.4993599 
##    user  system elapsed 
##   0.047   0.008   0.055
res <- funtest(1000) 
cpu.time('\nSample median of', 1000, 'values =', res) 
## Time of last operation: 
## Sample median of 1000 values = 0.4632395 
##    user  system elapsed 
##   0.007   0.000   0.007 
## Total time:
##    user  system elapsed 
##   0.054   0.008   0.062

Otro paquete que puede ser de utilidad es microbenchmark (si se quieren estudiar con más detalle los tiempos de computación; aunque en este libro no será el caso…). 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.

Capítulo 3

Generación de números pseudoaleatorios

Como ya se comentó, los distintos métodos de simulación requieren disponer de secuencias de números pseudoaleatorios que imiten las propiedades de generaciones independientes de una distribución \(\mathcal{U}(0, 1)\). En primer lugar nos centraremos en el caso de los generadores congruenciales. A pesar de su simplicidad, podrían ser adecuados en muchos casos y constituyen la base de los generadores avanzados habitualmente considerados. Posteriormente se dará una visión de las diferentes herramientas para estudiar la calidad de un generador de números pseudoaleatorios.

3.1 Generadores congruenciales

En los generadores congruenciales lineales se considera una combinación lineal de los últimos \(k\) enteros generados y se calcula su resto al dividir por un entero fijo \(m\). En el método congruencial simple (de orden \(k = 1\)), partiendo de una semilla inicial \(x_0\), el algoritmo secuencial es el siguiente:

\[\begin{aligned} x_i & = (ax_{i-1} + c) \quad\text{mod} \ m\\ u_i & = \dfrac{x_i}{m}\\ i & = 1, 2, … \end{aligned}\]

donde \(a\) (multiplicador), \(c\) (incremento) y \(m\) (módulo) son parámetros enteros del generador fijados de antemano. Si \(c = 0\) el generador se denomina congruencial multiplicativo (Lehmer, 1951) y en caso contrario se dice que es mixto (Rotenburg, 1960).

Este método está implementado en el fichero RANDC.R, tratando de imitar el funcionamiento de un generador de R (aunque de un forma no muy eficiente):

# Generador congruencial de números pseudoaleatorios 
# ================================================== 
# -------------------------------------------------- 
# initRANDC(semilla,a,c,m) 
# Selecciona el generador congruencial 
# Por defecto RANDU de IBM con semilla del reloj 
# OJO: No se hace ninguna verificación de los parámetros 

initRANDC <- function(semilla=as.numeric(Sys.time()), a=2^16+3, c=0, m=2^31) {  
.semilla <<- as.double(semilla) %% m #Cálculos en doble precisión 
 .a <<- a 
 .c <<- c 
 .m <<- m 

return(invisible(list(semilla=.semilla,a=.a,c=.c,m=.m))) #print(initRANDC()) 
} 
# -------------------------------------------------- 
# RANDC() 
# Genera un valor pseudoaleatorio con el generador congruencial 
# Actualiza la semilla (si no existe llama a initRANDC) 
RANDC <- function() { 
if (!exists(".semilla", envir=globalenv())) initRANDC() 
.semilla <<- (.a * .semilla + .c) %% .m 
return(.semilla/.m) 
} 
# -------------------------------------------------- 
# RANDCN(n) 
# Genera un vector de valores pseudoaleatorios con el generador congruencial 
# (por defecto de dimensión 1000) 
# Actualiza la semilla (si no existe llama a initRANDC) 
RANDCN <- function(n=1000) { 
x <- numeric(n)
for(i in 1:n) x[i]<-RANDC() 
return(x) 
# return(replicate(n,RANDC())) # Alternativa más rápida
} 
initRANDC(543210) # Fijar semilla 543210 para reproductibilidad

Nota: Para evitar problemas computacionales, se recomienda emplear el método de Schrage (ver Bratley et al., 1987; L’Ecuyer, 1988)

Ejemplos:

• \(c = 0\), \(a = 2^{16} + 3 = 65539\) y \(m = 2^{31}\), generador RANDU de IBM (no recomendable).

• \(c = 0\), \(a = 7^5 = 16807\) y \(m = 2^{31} - 1\) (primo de Mersenne), Park y Miller (1988) “minimal standar”, empleado por las librerías IMSL y NAG.

• \(c = 0\), \(a = 48271\) y \(m = 2^{31} - 1\) actualización del “minimal standar” propuesta por Park, Miller y Stockmeyer (1993).

Los parámetros y la semilla determinan los valores generados:

\[x_i = (a_i x_0 + c \frac{ a^i - 1}{a - 1}) \text{ mod } m\]

A pesar de su simplicidad, una adecuada elección de los parámetros permite obtener de manera eficiente secuencias de números “aparentemente” i.i.d. \(\mathcal{U}(0, 1).\)

El procedimiento habitual solía ser escoger \(m\) de forma que la operación del módulo se pudiese realizar de forma muy eficiente, para posteriormente seleccionar \(c\) y \(a\) de forma que el período fuese lo más largo posible (o suficientemente largo).

Teorema 3.1 (Hull y Dobell, 1962)

Un generador congruencial tiene período máximo \((p = m)\) si y solo si:

  1. \(c\) y \(m\) son primos relativos (i.e. \(m.c.d.(c, m) = 1\)).

  2. \(a - 1\) es múltiplo de todos los factores primos de \(m\) (i.e. \(a \equiv 1\) mod \(q\), para todo \(q\) factor primo de \(m\)).

  3. Si \(m\) es múltiplo de 4, entonces \(a - 1\) también lo ha de ser (i.e. \(m \equiv 0\) mod 4 \(\Rightarrow a \equiv 1\) mod 4).

Algunas consecuencias:

• Si \(m\) primo, \(p = m\) si y solo si \(a=1\).

• Un generador multiplicativo no cumple la condición 1 (\(m.c.d\)(0, \(m\)) = \(m\)).

Teorema 3.2 Un generador multiplicativo tiene período máximo (\(p = m - 1\)) si:

  1. \(m\) es primo.
  2. \(a\) es una raiz primitiva de \(m\) (i.e. el menor entero \(q\) tal que \(a^q = 1\) mod \(m\) es \(q = m - 1\)).

Además de preocuparse de la longitud del ciclo, las secuencias generadas deben aparentar muestras i.i.d. \(\mathcal{U}(0, 1)\).

Uno de los principales problemas es que los valores generados pueden mostrar una clara estructura reticular. Este es el caso por ejemplo del generador RANDU de IBM muy empleado en la década de los 70 (ver Figura 3.1) 2. Por ejemplo, el conjunto de datos randu contiene 400 tripletas de números sucesivos obtenidos con la implementación de VAX/VMS 1.5 (1977).

system.time(u <- RANDCN(9999)) # Generar 
##    user  system elapsed 
##   0.024   0.000   0.024
xyz <- matrix(u, ncol = 3, byrow = TRUE) 
library(plot3D) 
points3D(xyz[,1], xyz[,2], xyz[,3], colvar = NULL, phi = 60, 
theta = -50, pch = 21, cex = 0.2)

\[ \text {Figura 3.1: Grafico de dispersión de tripletas del generador RANDU de IBM (contenidas en 15 planos).}\]

En general todos los generadores de este tipo van a presentar estructuras reticulares. Marsaglia (1968) demostró que las \(k\)-uplas de un generadores multiplicativo están contenidas en a lo sumo \((k!m)^{1/k}\) hiperplanos paralelos (para más detalles sobre la estructura reticular, ver por ejemplo Ripley, 1987, sección 2.7). Por tanto habría que seleccionar adecuadamente \(m\) y \(c\) (\(a\) solo influiría en la pendiente) de forma que la estructura reticular sea impreceptible teniendo en cuenta el número de datos que se pretende generar (por ejemplo de forma que la distancia mínima entre los puntos sea próxima a la esperada en teoría).

Se han propuesto diversas pruebas (ver Sección 3.2) para determinar si un generador tiene problemas de este tipo y se han realizado numerosos estudios para determinadas familias (e.g. Park y Miller, 1988, estudiaron que parámetros son adecuados para \(m = 2^{31} - 1\)).

• En cualquier caso, se recomienda considerar un “periodo de seguridad” \(\approx \sqrt{p}\) para evitar este tipo de problemas.

• Aunque estos generadores tiene limitaciones en su capacidad para producir secuencias muy largas de números i.i.d. \(\mathcal{U}(0, 1)\), es un elemento básico en generadores más avanzados.

3.1.1 Otros generadores

Se han considerado diversas extensiones del generador congruencial lineal simple:

• Lineal múltiple: \(x_i = a_0 + a_1x_{i-1} + a_2c_{i-2} + \cdots + a_k x _{i-k}\) mod \(m\), con periodo \(p \leq m^k -1\).

• No lineal: \(x_i = f (x_{i-1}, x_{i-2}, \cdots , x_{i-k})\) mod m. Por ejemplo \(x_i = a_0 + a_1 x_{i-1} + a_2 x^2_{i-1}\) mod \(m\).

• Matricial: \(x_i = A_0 + A_1 x_{i-1} + A_2 x_{i-2} + \dots + A_k x_{i-k}\) mod \(m\).

Un ejemplo de generador congruencia lineal múltiple es el denominado generador de Fibonacci retardado (Fibonacci-lagged generator; Knuth, 1969):

\[x_n = (x_{n-37} + x_{n-100}) \ \text{mod} \ 2^{30},\]

con un período aproximado de \(2^{129}\) y que puede ser empleado en R (lo cual no sería en principio recomendable; ver Knuth Recent News 2002) estableciendo kind a "Knuth-TAOCP-2002" o Knuth-TAOCP" en la llamada a set.seed() o RNGkind().

Un caso particular del generador lineal múltiple son los denominados generadores de registros desfasados (más relacionados con la Criptografía). Se generan bits de forma secuencial considerando \(m = 2\) y \(a_i \in {0, 1}\) y se van combinando \(l\) bits para obtener valores en el intervalo (0, 1), por ejemplo \(u_i = 0.\) \(x_{it+1}x_{it+2} … x_{it+l}\), siendo \(t\) un parámetro denominado aniquilación (Tausworthe, 1965). Los cálculos se pueden realizar rápidamente mediante operaciones lógicas (los sumandos de la combinación lineal se traducen en un “o” exclusivo XOR), empleando directamente los registros del procesador (ver por ejemplo, Ripley, 1987, Algoritmo 2.1).

Otras alternativas consisten en la combinanción de varios generadores, las más empleadas son:

• Combinar las salidas: por ejemplo \(u_i=\sum_{l=1}^L u_i^{(l)}\) mod 1, donde \(u_i^{(l)}\) es el \(i\)-ésimo valor obtenido con el generador \(l\).

• Barajar las salidas: por ejemplo se crea una tabla empleando un generador y se utiliza otro para seleccionar el índice del valor que se va a devolver y posteriormente actualizar.

El generador L’Ecuyer-CMRG (L’Ecuyer, 1999), empleado como base para la generación de múltiples secuencias en el paquete parallel, combina dos generadores concruenciales lineales múltiples de orden \(k = 3\) (el periodo aproximado es \(2^{191}\)).

Ejercicio 3.1 (Análisis de un generador congruencial).

Considera el generador congruencial definido por:

\[\begin{aligned} x_{n+1} &= (5x_n +1) \text{ mod } 512, \\ u_{n+1} &= \frac{x_{n+1}}{512}, n=0,1,\dots \end{aligned}\]

(de ciclo máximo).

a) Generar 500 valores de este generador, obtener el tiempo de CPU, representar su distribución mediante un histograma (en escala de densidades) y compararla con la densidad teórica.

initRANDC(321, 5, 1, 512) # Establecer semilla y parámetros
nsim <- 500
system.time(u <- RANDCN(nsim)) # Generar
##    user  system elapsed 
##   0.000   0.001   0.001
hist(u, freq = FALSE)
abline(h = 1) # Densidad uniforme

\[ \text {Figura 3.2: Histograma de los valores generados}\]

En este caso concreto la distribución de los valores generados es aparentemente más uniforme de lo que cabría esperar, lo que induciría a sospechar de la calidad de este gerenador.

b) Calcular la media de las simulaciones (mean) y compararla con la teórica.

La aproximación por simulación de la media teórica es:

mean(u)
## [1] 0.4999609

La media teórica es \(0.5\). Error absoluto \(3.90625 \times 10^{-5}\).

c) Aproximar (mediante simulación) la probabilidad del intervalo (0.4; 0.8) y compararla con la teórica.

La probabilidad teórica es \(0.8 - 0.4 = 0.4\) La aproximación mediante simulación:

sum((0.4 < u) & (u < 0.8))/nsim
## [1] 0.402
mean((0.4 < u) & (u < 0.8)) # Alternativa
## [1] 0.402

3.2 Análisis de la calidad de un generador

Para verificar si un generador tiene las propiedades estadísticas deseadas hay disponibles una gran cantidad de test de hipótesis y métodos gráficos, incluyendo métodos genéricos (de bondad de ajuste y aleatoriedad) y contrastes específicos para generadores aleatorios. Se trata principalmente de contrastar si las muestras generadas son i.i.d. \(\mathcal{U}(0, 1)\) (análisis univariante). Aunque los métodos más avanzados tratan normalmente de contrastar si las \(d\)-uplas:

\[(U_{t+1}, U_{t+2}, \dots, U_{t+d});\text{ } t=(i-1)d, i=1, \dots ,m\]

son i.i.d. \(\mathcal{U}(0, 1)^d\) (uniformes independientes en el hipercubo; análisis multivariante).

En esta sección emplearemos únicamente métodos genéricos, ya que también pueden ser de utilidad para evaluar generadores de variables no uniformes y para la construcción de modelos del sistema real (e.g. para modelar variables que se tratarán como entradas del modelo general). Sin embargo, los métodos clásicos pueden no ser muy adecuados para evaluar generadores de números pseudoaleatorios (e.g. L’Ecuyer y Simard, 2007).

Hay que destacar algunas diferencias entre el uso de este tipo de métodos en inferencia y en simulación. Por ejemplo, si empleamos un constrate de hipótesis del modo habitual, desconfiamos del generador si la muestra (secuencia) no se ajusta a la distribución teórica \((p-\text{valor} \leq \alpha)\). En este caso además, también se sospecha si se ajusta demasiado bien a la distribución teórica \((p-\text{valor}\geq 1 - \alpha)\), lo que indicaría que no reproduce adecuadamente la variabilidad.

Uno de los contrastes más conocidos es el test ji-cuadrado de bondad de ajuste (chisq.test para el caso discreto). Aunque si la variable de interés es continua, habría que discretizarla (con la correspondiente perdida de información). Por ejemplo, se podría emplear la siguiente función (que imita a las incluídas en R):

#-------------------------------------------------------------------------------
# chisq.test.cont(x, distribution, nclasses, output, nestpar,...)
#-----------------------------------------------------------------------
# Realiza el test ji-cuadrado de bondad de ajuste para una distribución
# continua discretizando en intervalos equiprobables.
# Parámetros:
# distribution = "norm","unif", etc
# nclasses = floor(length(x)/5)
# output = TRUE
# nestpar = 0 = nº de parámetros estimados
# ... = parámetros distribución
# Ejemplo:
# chisq.test.cont(x, distribution = "norm", nestpar = 2,
# mean = mean(x), sd = sqrt((nx - 1) / nx) * sd(x))
#-----------------------------------------------------------------------
chisq.test.cont <- function(x, distribution = "norm",
nclasses = floor(length(x)/5), output = TRUE, nestpar = 0, ...) {
# Funciones distribución
q.distrib <- eval(parse(text = paste("q", distribution, sep = "")))
d.distrib <- eval(parse(text = paste("d", distribution, sep = "")))
# Puntos de corte
q <- q.distrib((1:(nclasses - 1))/nclasses, ...)
tol <- sqrt(.Machine$double.eps)
xbreaks <- c(min(x) - tol, q, max(x) + tol)
# Gráficos y frecuencias
if (output) {
xhist <- hist(x, breaks = xbreaks, freq = FALSE,
lty = 2, border = "grey50")
curve(d.distrib(x, ...), add = TRUE)
} else {
xhist <- hist(x, breaks = xbreaks, plot = FALSE)
}
#Cálculo estadístico y p-valor
O <- xhist$counts # Equivalente a table(cut(x, xbreaks)) pero más eficiente
E <- length(x)/nclasses
DNAME <- deparse(substitute(x))
METHOD <- "Pearson's Chi-squared test"
STATISTIC <- sum((O - E)^2/E)
names(STATISTIC) <- "X-squared"
PARAMETER <- nclasses - nestpar - 1
names(PARAMETER) <- "df"
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
# Preparar resultados
classes <- format(xbreaks)
classes <- paste("(", classes[-(nclasses + 1)], ",", classes[-1], "]",
sep = "")
RESULTS <- list(classes = classes, observed = O, expected = E, residuals = (O -
E)/sqrt(E))
if (output) {
cat("\nPearson's Chi-squared test table\n")
print(as.data.frame(RESULTS))
}
if (any(E < 5))
warning("Chi-squared approximation may be incorrect")
structure(c(list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL,
method = METHOD, data.name = DNAME), RESULTS), class = "htest")
}

Continuando con el generador congruencial anterior, obtendríamos:

chisq.test.cont(u, distribution = "unif",
nclasses = 10, nestpar = 0, min = 0, max = 1)

## 
## Pearson's Chi-squared test table
##                          classes observed expected  residuals
## 1  (-1.490116e-08, 1.000000e-01]       51       50  0.1414214
## 2  ( 1.000000e-01, 2.000000e-01]       49       50 -0.1414214
## 3  ( 2.000000e-01, 3.000000e-01]       49       50 -0.1414214
## 4  ( 3.000000e-01, 4.000000e-01]       50       50  0.0000000
## 5  ( 4.000000e-01, 5.000000e-01]       51       50  0.1414214
## 6  ( 5.000000e-01, 6.000000e-01]       51       50  0.1414214
## 7  ( 6.000000e-01, 7.000000e-01]       49       50 -0.1414214
## 8  ( 7.000000e-01, 8.000000e-01]       50       50  0.0000000
## 9  ( 8.000000e-01, 9.000000e-01]       50       50  0.0000000
## 10 ( 9.000000e-01, 9.980469e-01]       50       50  0.0000000
## 
##  Pearson's Chi-squared test
## 
## data:  u
## X-squared = 0.12, df = 9, p-value = 1

Figura 3.3: Gráfico resultante de aplicar la función ‘chisq.test.cont()‘ comparando el histograma de los valores generados con la densidad uniforme.

Como se muestra en la Figura 3.3 el histograma de la secuencia generada es muy plano (comparado con lo que cabría esperar de una muestra de tamaño 500 de una uniforme), y consecuentemente el \(p\)-valor del contraste ji-cuadrado es prácticamente 1, lo que indicaría que este generador no reproduce adecuadamente la variabilidad de una distribución uniforme.

Otro contraste de bondad de ajuste muy conocido es el test de Kolmogorov-Smirnov, implementado en ks.test.

Ejercicio 3.2 (Análisis de un generador congruencial, continuación).

Continuando con el generador congruencial del Ejercicio anterior:

initRANDC(321, 5, 1, 512)
nsim <- 500
system.time(u <- RANDCN(nsim))

a) Realizar el contraste de Kolmogorov-Smirnov para estudiar el ajuste a una \(\mathcal{U}(0, 1)\).

Este contraste de hipótesis compara la función de distribución bajo la hipótesis nula con la función de distribución empírica, representadas en la Figura 3.4:

# Distribución empírica
curve(ecdf(u)(x), type = "s", lwd = 2)
curve(punif(x, 0, 1), add = TRUE)

\[ \text {Figura 3.4: Comparación de la distribución empírica de la secuencia generada con la función de distribución uniforme.}\]

Podemos realizar el contraste con el siguiente código:

# Test de Kolmogorov-Smirnov
ks.test(u, "punif", 0, 1)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  u
## D = 0.0033281, p-value = 1
## alternative hypothesis: two-sided

b) Obtener el gráfico secuencial y el de dispersión retardado, ¿se observa algún problema?

  • Gráfico secuencial:
plot(as.ts(u))

\[ \text {Figura 3.5: Gráfico secuencial de los valores generados}\]

  • Gráfico de dispersión retardado:
plot(u[-nsim],u[-1])

\[ \text {Figura 3.6: Gráfico de dispersión retardado de los valores generados}\]

c) Estudiar las correlaciones del vector \((u_i, u_{i+k}),\) con \(k = 1, \dots , 10\). Contrastar si son nulas.

  • Correlaciones:
acf(u)

\[ \text{Figura 3.7: Autocorrelaciones de los valores generados.}\]

Test de Ljung-Box:

Box.test(u, lag = 10, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  u
## X-squared = 22.533, df = 10, p-value = 0.01261

3.2.1 Repetición de contrastes

Los contrastes se plantean habitualmente desde el punto de vista de la inferencia estadística en la práctica: se realiza una prueba sobre la única muestra disponible. Si se realiza una única prueba, en las condiciones de \(H_0\) hay una probabilidad \(\alpha\) de rechazarla. En simulación tiene mucho más sentido realizar un gran número de pruebas:

• La proporción de rechazos debería aproximarse al valor de \(\alpha\)(se puede comprobar para distintos valores de \(\alpha\)).

• La distribución del estadístico debería ajustarse a la teórica bajo \(H_0\) (se podría realizar un nuevo contraste de bondad de ajuste).

• Los \(p\)-valores obtenidos deberían ajustarse a una \(\mathcal{U}(0, 1)\) (se podría realizar también un contraste de bondad de ajuste).

Este procedimiento es también el habitual para validar un método de contraste de hipótesis por simulación.

Ejemplo 3.1

Consideramos el generador congruencial RANDU:

# Valores iniciales
initRANDC(543210) # Fijar semilla para reproductibilidad
# set.seed(543210)
n <- 500
nsim <- 1000
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)
# Realizar contrastes
for(isim in 1:nsim) {
u <- RANDCN(n) # Generar
# u <- runif(n)
tmp <- chisq.test.cont(u, distribution="unif",
nclasses=100, output=FALSE, nestpar=0, min=0, max=1)
estadistico[isim] <- tmp$statistic
pvalor[isim] <- tmp$p.value
}

Proporción de rechazos:

# cat("\nProporción de rechazos al 1% =", sum(pvalor < 0.01)/nsim, "\n")
cat("\nProporción de rechazos al 1% =", mean(pvalor < 0.01), "\n")
## 
## Proporción de rechazos al 1% = 0.014
# cat("Proporción de rechazos al 5% =", sum(pvalor < 0.05)/nsim, "\n")
cat("Proporción de rechazos al 5% =", mean(pvalor < 0.05), "\n")
## Proporción de rechazos al 5% = 0.051
# cat("Proporción de rechazos al 10% =", sum(pvalor < 0.1)/nsim, "\n")
cat("Proporción de rechazos al 10% =", mean(pvalor < 0.1), "\n")
## Proporción de rechazos al 10% = 0.112

Análisis del estadístico contraste:

# Histograma
hist(estadistico, breaks = "FD", freq=FALSE)
curve(dchisq(x,99), add=TRUE)

# Test ji-cuadrado
# chisq.test.cont(estadistico, distribution="chisq", nclasses=20, nestpar=0, df=99)
# Test de Kolmogorov-Smirnov
ks.test(estadistico, "pchisq", df=99)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  estadistico
## D = 0.023499, p-value = 0.6388
## alternative hypothesis: two-sided

Análisis de los \(p\)-valores:

# Histograma
hist(pvalor, freq=FALSE)
abline(h=1) # curve(dunif(x,0,1), add=TRUE)

# Test ji-cuadrado
# chisq.test.cont(pvalor, distribution="unif", nclasses=20, nestpar=0, min=0, max=1)
# Test de Kolmogorov-Smirnov
ks.test(pvalor, "punif", min=0, max=1)
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pvalor
## D = 0.023499, p-value = 0.6388
## alternative hypothesis: two-sided

Adicionalmente, si queremos estudiar la proporción de rechazos (el tamaño del contraste) para los posibles valores de \(\alpha\), podemos emplear la distribución empírica del \(p\)-valor (proporción de veces que resultó menor que un determinado valor):

# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2,
xlab = 'Nivel de significación', ylab = 'Proporción de rechazos')
abline(a = 0, b = 1, lty = 2) 
curve(punif(x, 0, 1), add = TRUE)

3.2.2 Baterías de contrastes

Hay numerosos ejemplos de generadores que pasaron diferentes test de uniformidad y aleatoriedad pero que fallaron estrepitosamente al considerar nuevos contrastes diseñados específicamente para generadores aleatorios (ver Marsaglia et al., 1990). Por este motivo, el procedimiento habitual en la práctica es aplicar un número más o menos elevado de contrastes (de distinto tipo y difíciles de pasar, e.g. Marsaglia y Tsang, 2002), de forma que si el generador los pasa tendremos mayor confianza en que sus propiedades son las adecuadas. Este conjunto de pruebas es lo que se denomina batería de contrastes. Una de las primeras se introdujo en Knuth (1969) y de las más recientes podríamos destacar:

• Diehard tests (The Marsaglia Random Number CDROM, 1995): http://www.stat.fsu.edu/pub/diehard (versión archivada el 2016-01-25).

• Dieharder (Brown y Bauer, 2003): Dieharder Page, paquete RDieHarder.

• TestU01 (L’Ecuyer y Simard, 2007): http://simul.iro.umontreal.ca/testu01/tu01.html.

• NIST test suite (National Institute of Standards and Technology, USA, 2010): http://csrc.nist. gov/groups/ST/toolkit/rng.

Para más detalles, ver por ejemplo cita3:

• Marsaglia, G. y Tsang, W.W. (2002). Some difficult-to-pass tests of randomness. Journal of Statistical Software, 7(3), 1-9.

• Demirhan, H. y Bitirim, N. (2016). CryptRndTest: an R package for testing the cryptographic randomness. The R Journal, 8(1), 233-247.

Capítulo 4

Análisis de resultados de simulación

En este capítulo nos centraremos en la aproximación mediante simulación de la media teórica de un estadístico a partir de la media muestral de una secuencia de simulaciones de dicho estadístico. La aproximación de una probabilidad sería un caso particular considerando una variable de Bernouilli.

En primer lugar se tratará el análisis de la convergencia y la precisión de la aproximación por simulación. Al final del capítulo se incluye una breve introducción a los problemas de estabilización y dependencia (con los que nos solemos encontrar en simulación dinámica y MCMC).

4.1 Convergencia

Supongamos que estamos interesados en aproximar la media teórica \(\mu = E(x)\) a partir de una secuencia i.i.d. \(X_1, X_2, \dots, X_n\) mediante la media muestral \(\overline{X_n}\). Una justificación teórica de la validez de la aproximación obtenida mediante simulación es la ley (débil) de los grandes números:

• Si \(X_1, X_2, \cdots\) es una secuencia de v.a.’s independientes con:

\[ E(X_i)= \mu \text{ y } Var(X_i)=\sigma ^2 < \infty , \] entonces \(\overline{X_n} = (X_1+ \cdots + X_n)/n\) converge en probabilidad a \(\mu\). i.e. para cualquier \(\varepsilon > 0\):

\[ \lim_{n \to \infty} P( \vert \overline{X_n} - \mu \vert < \varepsilon) = 1 \]

• La ley fuerte establece la convergencia casi segura.

Ejemplo 4.1. Aproximación de una probabilidad

Simulamos una distribución de Bernoulli de parámetro \(p = 0.5\):

p <- 0.5
set.seed(1)
nsim <- 10000
# nsim <- 100
rx <- runif(nsim) <= p

La aproximación por simulación de \(p\) será:

mean(rx)
## [1] 0.5047

Podemos generar un gráfico con la evolución de la aproximación con el siguiente código:

plot(cumsum(rx)/1:nsim, type="l", lwd=2, xlab="Número de generaciones",
ylab="Proporción muestral", ylim=c(0,1))
abline(h = mean(rx), lty = 2)
# valor teórico
abline(h = p)

\[ \text{Figura 4.1: Aproximación de la proporción en función del número de generaciones.}\]

4.1.1 Detección de problemas de convergencia

Una suposición crucial es que las variables \(X_i\) deben tener varianza finita (realmente esta suposición puede relajarse: \(E (\vert X_i\vert) < \infty\)). En caso contrario la media muestral puede no converger a una constante. Un ejemplo conocido es la distribución de Cauchy:

set.seed(1)
nsim <- 10000
rx <- rcauchy(nsim)
plot(cumsum(rx)/1:nsim, type="l", lwd=2,
xlab="Número de generaciones", ylab="Media muestral")

Para detectar problemas de convergencia es recomendable representar la evolución de la aproximación de la característica de interés (sobre el número de generaciones), además de realizar otros análisis descriptivos de las simulaciones. Por ejemplo, en este caso podemos observar los valores que producen estos saltos mediante un gráfico de cajas:

boxplot(rx)

4.2 Estimación de la precisión

En el caso de la media muestral \(\overline{𝑋}_ 𝑛\), un estimador insesgado de \(𝑉 𝑎𝑟 (𝑋_𝑛) = 𝜎^2/𝑛\) es la cuasi-varianza muestral:

\[\widehat{Var} (𝑋_𝑛) =\frac{\widehat{𝑆}^2}{ 𝑛}\]

Figura 4.2: Evolución de la media muestral de una distribución de Cauchy en función del número de generaciones.

\[\text{Figura 4.3: Gráfico de cajas de 10000 generaciones de una distribución de Cauchy.}\]

con:

\[\widehat{S}^2_n =\dfrac{1}{n-1}\sum_{i=1 }^{n}(X_{i}-\overline{X})^2\]

En el caso de una proporción \[\hat{𝑝}_n:\]

\[\widehat{𝑉𝑎𝑟} (\hat{𝑝}_n ) =\dfrac{\hat{𝑝}_n(1-\hat{p}_n)}{n-1}\]

aunque se suele emplear la varianza muestral.

Los valores obtenidos servirían como medidas básicas de la precisión de la aproximación, aunque su principal aplicación es la construcción de intervalos de confianza.

4.3 Teorema Central del Limite

Si \(X_1, X_2, \cdots\) es una secuencia de v.a.’s indeoendientes con \(E(X_i)=\mu\) y \(Var(X_i)=\sigma^2<\infty\), entonces:

\[Z_n=\frac{\bar{X_n}-\mu}{\frac{\sigma}{\sqrt{n}}}\xrightarrow{d}N(0,1)\] i.e. \[\lim_{n \to \infty}F_{Z_n}(z)=\varPhi(z)\] por tanto, un intervalo de confianza asintotico para \(\mu\) es:

\[IC_{1-\alpha}(\mu)= \left( \bar{X_n}-z_{1-\alpha/2}\frac{\hat{S_n}}{\sqrt{n}} , \bar{X_n}+z_{1-\alpha/2}\frac{\hat{S_n}}{\sqrt{n}} \right)\].

Podemos considerar que \(z_{1-\alpha/2}\frac{\hat{S_n}}{\sqrt{n}}\) es la precisión obtenida (con nivel de confianza \(1 − \alpha\)). La convergencia de la aproximación, además de ser aleatoria, se podría considerar lenta. La idea es que para doblar la precisión (disminuir el error a la mitad), necesitaríamos un número de generaciones cuatro veces mayor. Pero una ventaja, es que este error no depende del número de dimensiones (en el caso multidimensional puede ser mucho más rápida que otras alternativas numéricas).

Ejemplo 4.2. Aproximación de la media de una distribución normal

xsd <- 1
xmed <- 0
set.seed(1)
nsim <- 1000
rx <- rnorm(nsim, xmed, xsd)

La aproximación por simulación de la media será:

mean(rx)
## [1] -0.01164814

Como medida de la precisión de la aproximación podemos considerar (se suele denominar error de la aproximación):

2*sd(rx)/sqrt(nsim)
## [1] 0.06545382

(es habitual emplear 2 en lugar de 1.96, lo que se correspondería con \(1 − \alpha = 0.9545\) en el caso de normalidad). Podemos añadir también los correspondientes intervalos de confianza al gráfico de convergencia:

n <- 1:nsim
est <- cumsum(rx)/n
# Errores estándar calculados con la varianza muestral por comodidad:
esterr <- sqrt(cumsum((rx-est)^2))/n
plot(est, type = "l", lwd = 2, xlab = "Número de generaciones",
ylab = "Media y rango de error", ylim = c(-1, 1))
abline(h = est[nsim], lty=2)
lines(est + 2*esterr, lty=3)
lines(est - 2*esterr, lty=3)
abline(h = xmed)

\[\text{Figura 4.4: Gráfico de convergencia incluyendo el error de la aproximación.}\]

4.4 Determinación del número de generaciones

Normalmente el valor de 𝑛 se toma del orden de varias centenas o millares. En los casos en los que la simulación se utiliza para aproximar una característica central de la distribución (como una media) puede bastar un número de simulaciones del orden de \(n = 100, 200, 500\). Sin embargo, en otros casos pueden ser necesarios valores del tipo \(n = 1000, 2000, 5000, 10000\). En muchas ocasiones puede interesar obtener una aproximación con un nivel de precisión fijado. Para una precisión absoluta \(\varepsilon\), se trata de determinar \(n\) de forma que:

\[z_{1-\alpha/2}\frac{\hat{S_n}}{\sqrt{n}}<\varepsilon\] Un algoritmo podría ser el siguiente:

  1. Hacer \(j=0\) y fijar un tamaño inicial \(𝑛_0\) (e.g. \(30\) ó \(60\)).

  2. Generar \(\{X_i\}_{i=1}^{n_0}\) y calcular \(\overline{X_{n_0}}\) y \(\widehat{S}_{n_0}\) .

  3. Mientras \(z_{1−\alpha/2}\widehat{S}_{n_j}/\sqrt{n_j}>\varepsilon\) hacer:

. \(𝑗 = 𝑗 + 1\).

. \[n_j =\lceil (z_{1−\alpha/2}\widehat{S}_{n_{j-1}}/\varepsilon)^2\rceil\].

. Generar \(\{X_i\}_{i=n_{j-1}+1}^{n_0}\) y calcular \(\overline{X_{n_j}}\) y \(\widehat{S_{n_j}}\).

  1. devolver \(\overline{X_{n_j}}\) y \(z_{1−\alpha/2}\widehat{S}_{n_j}/\sqrt{n_j}\).

Para una precisión relativa \(\varepsilon |\mu|\) se procede análogamente de forma que:

\[z_{1−\alpha /2}\widehat{S}_{n_j}<\varepsilon |\overline{X_n}|\].

4.5 El problema de la dependencia

En el caso de dependencia, la estimación de la precisión se complica:

\[Var(\overline{X}) = \frac{1}{n^2}\left ( \sum_{i=1}^{n}Var(X_{i})+2 \sum_{i<j}Cov(X_{i},X_{j})\right ).\]

Ejemplo 4.3. Aproximación de una proporción bajo dependencia (cadena de Markov)

Supongamos que en A Coruña llueve de media \(1/3\) días al año, y que la probabilidad de que un día llueva solo depende de lo que ocurrió el día anterior, siendo \(0.94\) si el día anterior llovió y \(0.03\) si no.

Podemos generar valores de la variable indicadora de día lluvioso con el siguiente código:

# Variable dicotómica 0/1 (FALSE/TRUE)
set.seed(1)
nsim <- 10000
alpha <- 0.03 # prob. de cambio si seco
beta <- 0.06 # prob. de cambio si lluvia
rx <- logical(nsim) # x == "llueve"
rx[1] <- FALSE # El primer día no llueve
for (i in 2:nsim)
  rx[i] <- if (rx[i-1]) runif(1) > beta else runif(1) < alpha

Se podría pensar en emplear las expresiones anteriores:

n <- 1:nsim
est <- cumsum(rx)/n
esterr <- sqrt(est*(1-est)/(n-1))  # OJO! Supone independencia
plot(est, type="l", lwd=2, ylab="Probabilidad",
     xlab="Número de simulaciones", ylim=c(0,0.6))
abline(h = est[nsim], lty=2)
lines(est + 2*esterr, lty=2)
lines(est - 2*esterr, lty=2)
abline(h = 1/3, col="darkgray")  # Prob. teor. cadenas Markov

Figura 4.5: Gráfico de convergencia incluyendo el error de la aproximación (calculado asumiendo independencia).

El gráfico de autocorrelaciones sugiere que si tomamos 1 de cada 25 podemos suponer independencia.

La aproximación de la proporción sería correcta (es consistente):

est[nsim]
## [1] 0.3038

Sin embargo, al ser datos dependientes esta aproximación del error estandar no es adecuada:

esterr[nsim]
## [1] 0.004599203

En este caso al haber dependencia positiva se produce una subestimación del verdadero error estandar.

acf(as.numeric(rx))

\[\text{Figura 4.6: Correlograma de la secuencia indicadora de días de lluvia.}\] 

lag <- 24
xlag <- c(rep(FALSE, lag), TRUE)
rxi <- rx[xlag]
acf(as.numeric(rxi))

Figura 4.7: Correlograma de la subsecuencia de días de lluvia obtenida seleccionando uno de cada 25 valores. 

nrxi <- length(rxi)
n <- 1:nrxi
est <- cumsum(rxi)/n
esterr <- sqrt(est*(1-est)/(n-1))
plot(est, type="l", lwd=2, ylab="Probabilidad",
xlab=paste("Número de simulaciones /", lag + 1), ylim=c(0,0.6))
abline(h = est[length(rxi)], lty=2)
lines(est + 2*esterr, lty=2) # Supone independencia
lines(est - 2*esterr, lty=2)
abline(h = 1/3, col="darkgray") # Prob. teor. cadenas Markov

Figura 4.8: Gráfico de convergencia de la aproximación de la probabilidad a partir de la subsecuencia de días de lluvia (calculando el error de aproximación asumiendo independencia).

Esta forma de proceder podría ser adecuada para tratar de aproximar la precisión:

esterr[nrxi]
## [1] 0.02277402

pero no sería eficiente para aproximar la media. Siempre será preferible emplear todas las observaciones.

Por ejemplo, se podría pensar en considerar las medias de grupos de 24 valores consecutivos y suponer que hay independencia entre ellas:

rxm <- rowMeans(matrix(rx, ncol = lag, byrow = TRUE))
nrxm <- length(rxm)
n <- 1:nrxm
est <- cumsum(rxm)/n
esterr <- sqrt(cumsum((rxm-est)^2))/n # Error estándar
plot(est, type="l", lwd=2, ylab="Probabilidad",
xlab=paste("Número de simulaciones /", lag + 1), ylim=c(0,0.6))
abline(h = est[length(rxm)], lty=2)
lines(est + 2*esterr, lty=2) # OJO! Supone independencia
lines(est - 2*esterr, lty=2)
abline(h = 1/3, col="darkgray") # Prob. teor. cadenas Markov

\[\text{Figura 4.9: Gráfico de convergencia de las medias por lotes.}\]

Esta es la idea del método de medias por lotes (batch means; macro-micro replicaciones) para la estimación de la varianza. En el ejemplo anterior se calcula el error estándar de la aproximación por simulación de la proporción:

esterr[nrxm]
## [1] 0.01569017

pero si el objetivo es la aproximación de la varianza (de la variable y no de las medias por lotes), habrá que reescalarlo adecuadamente. Supongamos que la correlación entre \(x_i\) y \(x_{i+k}\) es aproximadamente nula, y consideramos las subsecuencias (lotes) \((x_{t+1}, x_{t+2},..., x_{t+k})\) con \(t = (j-1)k\) y \(n=mk\). Entonces:

\[Var(\overline{X})=Var\left (\frac{1}{n}\sum_{i=1}^{n}X_{i} \right )=Var\left (\frac{1}{m}\sum_{j=1}^{m}\left (\frac{1}{k}\sum_{t=(i-1)k+1}^{ik}X_{t} \right ) \right )\] \[\approx \frac{1}{m^2}\sum_{j=1}^{m}Var\left (\frac{1}{k}\sum_{t=(i-1)k+1}^{ik}X_{t}\right )\approx \frac{1}{m} Var(\overline{X_k})\] donde \(\overline{X_k}\) es la media de una subsecuencia de longitud \(k\).

var.aprox <- nsim * esterr[length(rxm)]^2
var.aprox
## [1] 2.461814

Obtenida asumiendo independencia entre las medias por lotes, y que será una mejor aproximación que asumir independencia entre las generaciones de la variable:

var(rx)
## [1] 0.2115267

Alternativamente se podría recurrir a la generación de múltiples secuencias independientes entre sí:

# Variable dicotómica 0/1 (FALSE/TRUE)
set.seed(1)
nsim <- 1000
nsec <- 10
alpha <- 0.03 # prob de cambio si seco
beta <- 0.06 # prob de cambio si lluvia
rxm <- matrix(FALSE, nrow = nsec, ncol= nsim)
for (i in 1:nsec) {
# rxm[i, 1] <- FALSE # El primer día no llueve
# rxm[i, 1] <- runif(1) < 1/2 # El primer día llueve con probabilidad 1/2
rxm[i, 1] <- runif(1) < 1/3 # El primer día llueve con probabilidad 1/3 (ideal)
for (j in 2:nsim)
rxm[i, j] <- if (rxm[i, j-1]) runif(1) > beta else runif(1) < alpha
}

La idea sería considerar las medias de las series como una muestra independiente de una nueva variable y estimar su varianza de la forma habitual:

# Media de cada secuencia
n <- 1:nsim
est <- apply(rxm, 1, function(x) cumsum(x)/n)
matplot(n, est, type = 'l', lty = 3, col = "lightgray",
ylab="Probabilidad", xlab="Número de simulaciones")

# Aproximación
mest <- apply(est, 1, mean)
lines(mest, lwd = 2)
abline(h = mest[nsim], lty = 2)

# Precisión
mesterr <- apply(est, 1, sd)/sqrt(nsec)
lines(mest + 2*mesterr, lty = 2)
lines(mest - 2*mesterr, lty = 2)

# Prob. teor. cadenas Markov
abline(h = 1/3, col="darkgray")

\[\text{Figura 4.10: Gráfico de convergencia de la media de 10 secuencias generadas de forma independiente.}\] 

#Aproximación final
mest[nsim] # mean(rxm)
## [1] 0.3089
#Error estándar
mesterr[nsim]
## [1] 0.02403491

Este tipo de problemas se tratarán en la diagnosis de algoritmos de simulación Monte Carlo de Cadenas de Markov (MCMC). Aparecen también en la simulación dinámica (por eventos o cuantos).

4.5.1 Periodo de calentamiento

En el caso de simulación de datos dependientes (simulación dinámica) pueden aparecer problemas de estabilización. Puede ocurrir que el sistema evolucione lentamente en el tiempo hasta alcanzar su distribución estacionaria, siendo muy sensible a las condiciones iniciales con las que se comienzó la simulación. En tal caso resulta conveniente ignorar los resultados obtenidos durante un cierto período inicial de tiempo (denominado período de calentamiento o estabilización), cuyo único objeto es conseguir que se estabilice la distribución de probabilidad.

Como ejemplo comparamos la simulación del Ejemplo 4.3 con la obtenida considerando como punto de partida un día lluvioso (con una semilla distinta para evitar dependencia).

set.seed(2)
nsim <- 10000
rx2 <- logical(nsim)
rx2[1] <- TRUE # El primer día llueve
for (i in 2:nsim)
  rx2[i] <- if (rx2[i-1]) runif(1) > beta else runif(1) < alpha
n <- 1:nsim
est <- cumsum(rx)/n
est2 <- cumsum(rx2)/n
plot(est, type="l", ylab="Probabilidad", xlab="Número de simulaciones", ylim=c(0,0.6))
lines(est2, lty = 2)
# Ejemplo periodo calentamiento nburn = 2000
abline(v = 2000, lty = 3)
# Prob. teor. cadenas Markov
abline(h = 1/3, col="darkgray")

En estos casos puede ser recomendable ignorar los primeros valores generados (por ejemplo los primeros 2000) y recalcular los estadísticos deseados.

También trataremos este tipo de problemas en la diagnosis de algoritmos MCMC.

Ejemplo 4.4. Simulación de un proceso autorregresivo (serie de tiempo)

\[X_{t}=\mu+\rho *(X_{t-1}-\mu)+\varepsilon_{t}\] Podemos tener en cuenta que en este caso la varianza es:

\[Var(X_t)=E(X^2_{t})-\mu^2=\frac{\sigma ^2_{\varepsilon}}{1-\rho^2}\]

Establecemos los parámetros:

nsim <- 200 # Numero de simulaciones
xmed <- 0 # Media
rho <- 0.5 # Coeficiente AR
nburn <- 10 # Periodo de calentamiento (burn-in)

Se podría fijar la varianza del error:

evar <- 1
# Varianza de la respuesta
xvar <- evar / (1 - rho^2)

pero la recomendación sería fijar la varianza de la respuesta:

xvar <- 1
# Varianza del error
evar <- xvar*(1 - rho^2)

Para simular la serie, al ser un \(AR(1)\), normalmente simularíamos el primer valor

rx[1] <- rnorm(1, mean = xmed, sd = sqrt(xvar))

o lo fijamos a la media (en este caso nos alejamos un poco de la distribución estacionaria, para que el “periodo de calentamiento” sea mayor). Después generamos los siguientes valores de forma recursiva:

set.seed(1)
x <- numeric(nsim + nburn)
# Establecer el primer valor
x[1] <- -10
# Simular el resto de la secuencia
for (i in 2:length(x))
 x[i] <- xmed + rho*(x[i-1] - xmed) + rnorm(1, sd=sqrt(evar))
x <- as.ts(x)
plot(x)
abline(v = nburn, lty = 2)

\[\text{Figura 4.11: Ejemplo de una simulación de una serie de tiempo autorregresiva.}\]

y eliminamos el periodo de calentamiento:

rx <- x[-seq_len(nburn)]

Para simular una serie de tiempo en R se puede emplear la función del paquete base

En este caso el periodo de calentamiento se establece mediante el parámetro (que se fija automáticamente a un valor adecuado).

Por ejemplo, podemos generar este serie autoregressiva con:

rx2 <- arima.sim(list(order = c(1,0,0), ar = rho), n = nsim, n.start = nburn, sd = sqrt(evar))

La recomendación es fijar la varianza de las series simuladas si se quieren comparar resultados considerando distintos parámetros de dependencia.

4.6 Observaciones

• En el caso de que la característica de interés de la distribución de \(𝑋\) no sea la media, los resultados anteriores no serían en principio aplicables.

• Incluso en el caso de la media, las “bandas de confianza” obtenidas con el TCL son puntuales (si generamos nuevas secuencias de simulación es muy probable que no estén contenidas).

• En muchos casos (p.e. la generación de múltiples secuencias de simulación puede suponer un coste computacional importante), puede ser preferible emplear un método de remuestreo.

#–#

  1. El método de Montecarlo es un método de simulación que permite calcular estadísticamente el valor final de una secuencia de sucesos no deterministas (sujetos a variabilidad), como es el caso del plazo o el coste de un proyecto. En la práctica este análisis consiste en ejecutar varias veces los diferentes sucesos variando aleatoriamente su valor en función de la función estadística que los define, dando como resultado un conjunto de valores finales. Este conjunto permite calcular el valor medio y la variabilidad.↩︎

  2. Alternativamente se podría utilizar la función del paquete , y rotar la figura (pulsando con el ratón) para ver los hiperplanos: ↩︎

  3. También puede ser de interés el enlace Randomness Tests: A Literature Survey y la entidad certificadora (gratuita) en línea CAcert.↩︎