Tipos de líneas artísticas en simulación
Tipos de líneas artísticas en simulación
2020-03-12
2
Índice general
Prólogo 5
1 Introducción a la simulación 7
1.1 Conceptos básicos . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Generación de números (pseudo)aleatorios . . . . . . . . . . . . . 9
1.3 Números aleatorios puros . . . . . . . . . . . . . . . . . . . . . . 9
1.4 Números pseudoaleatorios . . . . . . . . . . . . . . . . . . . . . . 11
2 Números aleatorios en R 15
2.1 Opciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2 Paquetes de R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Tiempo de CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3
4 ÍNDICE GENERAL
Referencias 205
Bibliografía básica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
Bibliografía complementaria . . . . . . . . . . . . . . . . . . . . . . . . 205
ÍNDICE GENERAL 5
A Enlaces 207
7
8 ÍNDICE GENERAL
Capítulo 1
Introducción a la simulación
Podríamos definir la Simulación como una técnica que consiste en realizar expe-
rimentos de muestreo sobre el modelo de un sistema, con el objetivo de recopilar
información bajo determinadas condiciones.
9
10 CAPÍTULO 1. INTRODUCCIÓN A LA SIMULACIÓN
RAND publicó el libro A Million Random Digits with 100,000 Normal Devia-
tes que contenía números aleatorios generados mediante una ruleta electrónica
conectada a una computadora (ver Figura 1.1).
Figura 1.1: Líneas 10580-10594, columnas 21-40, del libro *A Million Random
Digits with 100,000 Normal Deviates*.
1.3.1 Inconvenientes:
Para su aplicación en el campo de la Estadística:
• Es necesario conocer su distribución.
• Los datos generados deberían ser i.i.d.
(sería también muy recomendable para otros casos) y en general:
• Reproductivilidad.
1.4. NÚMEROS PSEUDOALEATORIOS 13
1.3.2 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 (en el Capítulo 9 se tratarán
métodos de integración Monte Carlo), no dependen realmente de la aleatorie-
dad de la secuencia. Para evitar generaciones poco probables, se puede recurrir
a secuencias cuasi-aleatorias, también denominadas sucesiones de baja discre-
pancia (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 distribu-
ción uniforme, aunque no necesariamente para una realización concreta).
Por ejemplo, el paquete randtoolbox implementa métodos para la generación
de secuencias cuasi-aleatorias (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)
1.0
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
x2
x2
x2
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
x1 x1 x1
donde:
• 𝑘 es el orden del generador.
• (𝑥0 , 𝑥1 , ⋯ , 𝑥𝑘−1 ) es la semilla (estado inicial).
El periodo o longitud del ciclo es la longitud de la secuencia antes de que vuelva
a repetirse. Lo denotaremos por 𝑝.
Los números de la sucesión serán predecibles, conociendo el algoritmo y la semi-
lla.
• 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
(e.g. http://eprint.iacr.org/2007/419).
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 com-
portarse como si se tratase de una muestra genuina de datos independientes de
una 𝒰(0, 1). Otras propiedades de interés serían:
• Reproducibilidad a partir de la semilla.
• Periodo suficientemente largo.
1.4. NÚMEROS PSEUDOALEATORIOS 15
Números aleatorios en R
17
18 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
2.1 Opciones
La función RNGkind(kind = NULL, normal.kind = NULL) permite seleccionar
el tipo de generador (en negrita los valores por defecto):
• kindespecifica el generador aleatorio:
– “Wichmann-Hill”: Ciclo 6.9536 × 1012
– “Marsaglia-Multicarry”: Ciclo mayor de 260
– “Super-Duper”: Ciclo aprox. 4.6 × 1018 (S-PLUS)
– “Mersenne-Twister”: Ciclo 219937 − 1 y equidistribution en 623
dimensiones.
– “Knuth-TAOCP-2002”: Ciclo aprox. 2129 .
– “Knuth-TAOCP”
– “user-supplied”
• normal.kindselecciona el método de generación de normales (se trata-
rá más adelante). “Kinderman-Ramage”, “Buggy Kinderman-Ramage”,
“Ahrens-Dieter”, “Box-Muller”, “Inversion” , o “user-supplied”
2.2. PAQUETES DE R 19
2.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 numeros “true random” desde RAN-
DOM.ORG.
• randtoolbox implementa generadores más recientes (rngWELL) y genera-
ció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 se-
cuencias (para programación paralela).
• gls, rngwell19937, randaes, SuppDists, lhs, mc2d, fOptions, …
2.3 Ejercicios
Ejercicio 2.1.
## [1] 0.4996
mean(indice) # Alternativa
## [1] 0.4996
Nota: . R maneja internamente los valores lógicos como 1 (TRUE) y 0
(FALSE).
20 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
## [1] 0.7806
pi/4
## [1] 0.7853982
pi_aprox <- 4*mean(indice)
pi_aprox
## [1] 3.1224
Gráfico
# Colores y símbolos depediendo 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)
1.0
0.5
0.0
Y
−0.5
−1.0
Ejercicio 2.2.
2.3. EJERCICIOS 21
## [1] 0.4953
barplot(100*table(x)/nsim) # Representar porcentajes
50
40
30
20
10
0
0 1
## [1] 0.394
22 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
1.0
0.8
Proporción de caras
0.6
0.4
0.2
0.0
Número de lanzamientos
Ejercicio 2.3.
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:
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
Ejercicio 2.4.
n <- 4
lanz <- sample(1:6, replace=TRUE, size=n)
lanz
## [1] 3 5 1 6
6 %in% lanz
## [1] TRUE
b) Utilizar la función anterior para simular 𝑛𝑠𝑖𝑚 = 10000 jugadas de este
juego y calcular la proporción de veces que se gana la apuesta (obtener al
menos un 6 en 𝑛 lanzamientos), usando 𝑛 = 4. Comparar el resultado con
la probabilidad teórica 1 − (5/6)𝑛 .
24 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
set.seed(1)
n <- 4
nsim <- 10000
mean(replicate(nsim, deMere(n)))
## [1] 0.5148
1-(5/6)^n
## [1] 0.5177469
## [1] 0.5003313
CPUtimeprint()
##
## Tiempo última operación:
## user system elapsed
## 0.08 0.00 0.08
## Tiempo total operación:
## user system elapsed
## 0.08 0.00 0.08
funtest(1000)
## [1] 0.5050682
CPUtimeprint()
##
## Tiempo última operación:
## user system elapsed
## 0 0 0
## Tiempo total operación:
## user system elapsed
## 0.08 0.00 0.08
2.4.2 Paquetes de R
Por ejemplo, se puede emplear el paquete npsp (fichero cpu.time.R):
• Call cpu.time(restart = TRUE) where you want to start counting.
• Call cpu.time() to print/get total and/or partial (since the last call to
this function) real and CPU times.
# CPU time utilities
# ------------------
#' @seealso
#' \code{\link{proc.time}}, \code{\link{system.time}}, \code{\link{flush.console}}.
#' @export
cpu.time <- .cpu.time.ini()
Ejemplo:
cpu.time(reset = TRUE)
Generación de números
pseudoaleatorios con
distribución uniforme
29
30CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIF
# --------------------------------------------------
# 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
}
𝑎𝑖 − 1
𝑥𝑖 = (𝑎𝑖 𝑥0 + 𝑐 ) mod 𝑚
𝑎−1
library(plot3D)
points3D(xyz[,1], xyz[,2], xyz[,3], colvar = NULL, phi = 60, theta = -50, pch = 21, cex
z
En general todos los generadores de este tipo van a presentar estructuras reti-
culares. Marsaglia (1968) demostró que las 𝑘-uplas de un generadores multipli-
3.1. GENERADORES CONGRUENCIALES (LINEALES) 33
1/𝑘
cativo están contenidas en a lo sumo (𝑘!𝑚) 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 𝑚 y 𝑐 (𝑎 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 gene-
rador 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 𝑚 = 231 − 1).
• En cualquier caso, se recomienda considerar un “periodo de seguridad”
√
≈ 𝑝 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. 𝒰(0, 1), es un elemento básico en
generadores más avanzados.
con un período aproximado de 2129 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().
El generador Mersenne-Twister (Matsumoto y Nishimura, 1998), empleado por
defecto en R, de periodo 219937 − 1 y equidistribution en 623 dimensiones, se
puede expresar como un generador congruencial matricial lineal.
Un caso particular del generador lineal múltiple son los denominados genera-
dores de registros desfasados (más relacionados con la Criptografía). Se gene-
ran bits de forma secuencial considerando 𝑚 = 2 y 𝑎𝑖 ∈ {0, 1} y se van com-
binando 𝑙 bits para obtener valores en el intervalo (0, 1), por ejemplo 𝑢𝑖 =
34CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIF
Ejercicio 3.1.
Histogram of u
1.0
0.8
0.6
Density
0.4
0.2
0.0
## [1] 0.4999609
La media teórica es 0.5. Error absoluto 3.90625 × 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
𝑑
son i.i.d. 𝒰 (0, 1) (uniformes independientes en el hipercubo; análisis multiva-
riante). En el Apéndice B se describen algunos de estos métodos.
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 pa-
ra 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 pseudoalea-
torios (e.g. L’Ecuyer y Simard, 2007). La recomendación sería emplear baterías
de contrastes recientes, como las descritas en la Subsección 3.2.2.
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 (𝑝-valor ≤ 𝛼). En este caso además, también se
sospecha si se ajusta demasiado bien a la distribución teórica (𝑝-valor ≥ 1 − 𝛼),
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)
#-------------------------------------------------------------------------------
chisq.test.cont <- function(x, distribution = "norm", nclasses = floor(length(x)/5),
output = TRUE, nestpar = 0, ...) {
# Funciones distribución
3.2. ANÁLISIS DE LA CALIDAD DE UN GENERADOR 37
##
## Pearson's Chi-squared test table
38CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIF
Histogram of x
1.0
0.8
0.6
Density
0.4
0.2
0.0
Ejercicio 3.2.
0.4
0.2
0.0
##
40CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIF
Gráfico secuencial:
plot(as.ts(u))
1.0
0.8
0.6
as.ts(u)
0.4
0.2
0.0
Time
1.0
0.8
0.6
u[−1]
0.4
0.2
0.0
u[−nsim]
c) Estudiar las correlaciones del vector (𝑢𝑖 , 𝑢𝑖+𝑘 ), con 𝑘 = 1, … , 10. Contras-
tar si son nulas.
Correlaciones:
acf(u)
Series u
1.0
0.8
0.6
ACF
0.4
0.2
0.0
0 5 10 15 20 25
Lag
Test de Ljung-Box:
Box.test(u, lag = 10, type = "Ljung")
##
## Box-Ljung test
##
42CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIFORME
## data: u
## X-squared = 22.533, df = 10, p-value = 0.01261
Ejemplo 3.1.
# 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")
3.2. ANÁLISIS DE LA CALIDAD DE UN GENERADOR 43
##
## 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")
Histogram of estadistico
0.030
0.025
0.020
Density
0.015
0.010
0.005
0.000
estadistico
# 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
44CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNIF
Histogram of pvalor
1.2
1.0
0.8
Density
0.6
0.4
0.2
0.0
pvalor
# 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
1.0
0.8
Proporción de rechazos
0.6
0.4
0.2
0.0
Nivel de significación
3.3 Ejercicios
Ejercicio 3.3.
𝑥2𝑖 𝑘 𝑘
𝑥𝑖+1 = ⌊(𝑥2𝑖 − ⌊ 𝑘 ⌋ 10(2𝑘− 2 ) ) /10 2 ⌋
10(2𝑘− 2 )
# -------------------------------------------------
# initRANDVN(semilla,n)
# Inicia el generador
# n número de digitos centrales, 4 por defecto (debe ser un nº par)
# Por defecto semilla del reloj
# OJO: No se hace ninguna verificación de los parámetros
initRANDVN <- function(semilla = as.numeric(Sys.time()), n = 4) {
.semilla <<- as.double(semilla) %% 10^n # Cálculos en doble precisión
.n <<- n
.aux <<- 10^(2*n-n/2)
.aux2 <<- 10^(n/2)
return(invisible(list(semilla=.semilla,n=.n)))
}
# -------------------------------------------------
# RANDVN()
3.3. EJERCICIOS 47
# -------------------------------------------------
# RANDVNN(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 initRANDVN)
RANDVNN <- function(n = 1000) {
x <- numeric(n)
for(i in 1:n) x[i] <- RANDVN()
return(x)
# return(replicate(n,RANDVN())) # Alternativa más rápida
}
Ejercicio 3.4.
Análisis de resultados de
simulación
Work in progress…
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 intro-
ducció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 𝜇 = 𝐸 (𝑋)
a partir de una secuencia i.i.d. 𝑋1 , 𝑋2 , ⋯, 𝑋𝑛 mediante la media muestral 𝑋̄ 𝑛 .
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 𝑋1 , 𝑋2 , ⋯ es una secuencia de v.a.’s independientes con:
49
50 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
## [1] 0.5047
0.6
0.4
0.2
0.0
Número de generaciones
Una suposición crucial es que las variables 𝑋𝑖 deben tener varianza finita (real-
mente esta suposición puede relajarse: 𝐸 (|𝑋𝑖 |) < ∞). 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")
8
6
Media muestral
4
2
0
−2
Número de generaciones
15000
10000
5000
0
−5000
𝑆̂2
𝑉̂
𝑎𝑟 (𝑋 𝑛 ) =
𝑛
con:
𝑛
1 2
𝑆̂𝑛2 = ∑ (𝑋𝑖 − 𝑋) .
𝑛 − 1 𝑖=1
𝑝̂ (1 − 𝑝𝑛̂ )
𝑉̂
𝑎𝑟 (𝑝𝑛̂ ) = 𝑛 ,
𝑛−1
𝑋𝑛 − 𝜇 𝑑
𝑍𝑛 = → 𝑁 (0, 1)
√𝜎
𝑛
i.e. lim 𝐹𝑍𝑛 (𝑧) = Φ(𝑧). Por tanto, un intervalo de confianza asintótico para 𝜇
𝑛→∞
es:
𝑆̂ 𝑆̂
𝐼𝐶1−𝛼 (𝜇) = (𝑋 𝑛 − 𝑧1−𝛼/2 √𝑛 , 𝑋 𝑛 + 𝑧1−𝛼/2 √𝑛 ) .
𝑛 𝑛
𝑆̂
Podemos considerar que 𝑧1−𝛼/2 √𝑛 es la precisión obtenida (con nivel de con-
𝑛
fianza 1 − 𝛼).
La convergencia de la aproximación, además de ser aleatoria, se podría conside-
rar 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 multidimen-
sional puede ser mucho más rápida que otras alternativas numéricas).
xsd <- 1
xmed <- 0
set.seed(1)
nsim <- 1000
rx <- rnorm(nsim, xmed, xsd)
## [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 − 𝛼 =
0.9545 en el caso de normalidad). Podemos añadir también los correspondientes
intervalos de confianza al gráfico de convergencia:
54 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
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)
1.0
0.5
Media y rango de error
0.0
−0.5
−1.0
Número de generaciones
𝑛
3.3. Generar {𝑋𝑖 }𝑖=𝑛
𝑗
y calcular 𝑋 𝑛𝑗 y 𝑆̂𝑛𝑗 .
𝑗−1 +1
√
4. Devolver 𝑋 𝑛𝑗 y 𝑧1−𝛼/2 𝑆̂𝑛𝑗 / 𝑛𝑗 .
𝑆̂
𝑧1−𝛼/2 √𝑛 < 𝜀 ∣𝑋 𝑛 ∣ .
𝑛
Supongamos que en A Coruña llueve de media 1/3 días al año, y que la probabi-
lidad 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
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")
0.6
0.5
0.4 # Prob. teor. cadenas Markov
Probabilidad
0.3
0.2
0.1
0.0
Número de simulaciones
## [1] 0.3038
Sin embargo, al ser datos dependientes esta aproximación del error estandar no
es adecuada:
esterr[nsim]
## [1] 0.004599203
acf(as.numeric(rx))
Series as.numeric(rx)
1.0
0.8
0.6
ACF
0.4
0.2
0.0
0 10 20 30 40
Lag
Esta forma de proceder podría ser adecuada para tratar de aproximar la preci-
sión:
esterr[nrxi]
## [1] 0.02277402
58 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
Series as.numeric(rxi)
1.0
0.8
0.6
ACF
0.4
0.2
0.0
0 5 10 15 20 25
Lag
0.3
0.2
0.1
0.0
Número de simulaciones / 25
pero no sería eficiente para aproximar la media. Siempre será preferible emplear
4.5. EL PROBLEMA DE LA DEPENDENCIA 59
0.3
0.2
0.1
0.0
Número de simulaciones / 25
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 𝑋𝑖 y 𝑋𝑖+𝑘 es aproximadamente nula, y consideramos las
60 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
1 𝑛 1 𝑚 1 𝑖𝑘
̄ = 𝑉 𝑎𝑟 (
𝑉 𝑎𝑟 (𝑋) ∑ 𝑋𝑖 ) = 𝑉 𝑎𝑟 ⎛
⎜ ∑⎛ ⎜ ∑ 𝑋𝑡 ⎞
⎟⎞
⎟
𝑛 𝑖=1 𝑚 𝑗=1 𝑘
⎝ ⎝ 𝑡=(𝑖−1)𝑘+1 ⎠⎠
1 𝑚 1 𝑖𝑘
1
≈ ∑ 𝑉 𝑎𝑟 ⎛
⎜ ∑ 𝑋𝑡 ⎞
⎟ ≈ 𝑉 𝑎𝑟 (𝑋̄ 𝑘 )
𝑚2 𝑗=1 𝑘 𝑚
⎝ 𝑡=(𝑖−1)𝑘+1 ⎠
donde 𝑋̄ 𝑘 es la media de una subsecuencia de longitud 𝑘.
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 va-
riable:
var(rx)
## [1] 0.2115267
Alternativamente se podría recurrir a la generación de múltiples secuencias in-
dependientes 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",
4.5. EL PROBLEMA DE LA DEPENDENCIA 61
0.4
0.2
0.0
Número de simulaciones
# Aproximación final
mest[nsim] # mean(rxm)
## [1] 0.3089
# Error estándar
mesterr[nsim]
## [1] 0.02403491
Como ejemplo comparamos la simulación del Ejemplo 4.3 con la obtenida con-
siderando 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")
4.5. EL PROBLEMA DE LA DEPENDENCIA 63
0.6
0.5
0.4
Probabilidad
0.3
0.2
0.1
0.0
Número de simulaciones
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.
𝑋𝑡 = 𝜇 + 𝜌 ∗ (𝑋𝑡−1 − 𝜇) + 𝜀𝑡
Podemos tener en cuenta que en este caso la varianza es:
𝜎𝜀2
var(𝑋𝑡 ) = E(𝑋𝑡2 ) − 𝜇2 = .
1 − 𝜌2
xvar <- 1
# Varianza del error
evar <- xvar*(1 - rho^2)
−6
−8
−10
Time
rx <- x[-seq_len(nburn)]
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.
66 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
Capítulo 5
Simulación de variables
continuas
𝑈 = 𝐹 (𝑋) ∼ 𝒰 (0, 1) ,
ya que:
𝐹 −1 (𝑈 ) ∼ 𝑋
67
68 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
2. Devolver 𝑋 = 𝐹 −1 (𝑈 ).
1 − 𝑒−𝜆𝑥 si 𝑥 ≥ 0
𝐹 (𝑥) = {
0 si 𝑥 < 0
ln (1 − 𝑢)
1 − 𝑒−𝜆𝑥 = 𝑢 ⇔ 𝑥 = − ,
𝜆
el algoritmo para simular esta variable mediante el método de inversión es:
1. Generar 𝑈 ∼ 𝒰 (0, 1).
ln (1 − 𝑈 )
2. Devolver 𝑋 = − .
𝜆
En el último paso podemos emplear directamente 𝑈 en lugar de 1 − 𝑈 , ya que
1 − 𝑈 ∼ 𝒰 (0, 1). Esta última expresión para acelerar los cálculos es la que
denominaremos forma simplificada.
El código para implementar este algoritmo en R podría ser el siguiente:
# Simular vector exp(lambda)
tini <- proc.time()
lambda <- 2
nsim <- 10^5
set.seed(1)
U <- runif(nsim)
X <- -log(U)/lambda # -log(1-U)/lambda
Como se observa en la Figura 5.1 se trata de un método exacto (si está bien
implementado) y la distribución de los valores generados se aproxima a la dis-
tribución teórica como cabría esperar con una muestra de ese tamaño.
5.1. MÉTODO DE INVERSIÓN 69
2.5
2.0
1.5
Density
1.0
0.5
0.0
0 1 2 3 4 5
Forma
Nombre Densidad 𝐹 (𝑥) 𝐹 −1 (𝑈 ) simplificada
ln (1 − 𝑈 ) ln 𝑈
exp (𝜆) 𝜆𝑒−𝜆𝑥 , si 1 − 𝑒−𝜆𝑥 − −
𝜆 𝜆
(𝜆 > 0) 𝑥≥0
1 1 arctan 𝑥 1
Cauchy + tan (𝜋 (𝑈 − )) tan 𝜋𝑈
𝜋 (1 + 𝑥2 ) 2 𝜋 2
2 𝑥 2 𝑥2 √ √
Triangular (1 − ), (𝑥 − ) 𝑎 (1 − 1 − 𝑈 ) 𝑎 (1 − 𝑈 )
𝑎 𝑎 𝑎 2𝑎
en (0, 𝑎) si 0 ≤ 𝑥 ≤ 𝑎
𝑎
𝑎𝑏𝑎 𝑏 𝑏 𝑏
Pareto , si 1−( )
𝑥𝑎+1 𝑥 (1 − 𝑈 )
1/𝑎 𝑈 1/𝑎
(𝑎, 𝑏 > 0) 𝑥≥𝑏
1/𝛼 1/𝛼
𝛼 𝛼 (− ln (1 − 𝑈 )) (− ln 𝑈 )
Weibull 𝛼𝜆𝛼 𝑥𝛼−1 𝑒−(𝜆𝑥) ,1 − 𝑒−(𝜆𝑥) \
𝜆 𝜆
(𝜆, 𝛼 > 0) si 𝑥 ≥ 0
Ejercicio 5.1.
1.0
0.8
0.6
Density
0.4
0.2
0.0
−4 −2 0 2 4
Figura 5.2: Distribución de los valores generados de una doble exponencial me-
diante el método de inversión.
generar).
• Utilizar una aproximación a 𝐹 −1 (𝑢) (inversión aproximada).
4 4
siendo 𝐴 (𝑥) = ∑𝑖=0 𝑎𝑖 𝑥𝑖 y 𝐵 (𝑥) = ∑𝑖=0 𝑏𝑖 𝑥𝑖 con:
𝑎0 = −0.322232431088 𝑏0 = 0.0993484626060
𝑎1 = −1 𝑏1 = 0.588581570495
𝑎2 = −0.342242088547 𝑏2 = 0.531103462366
𝑎3 = −0.0204231210245 𝑏3 = 0.103537752850
𝑎4 = −0.0000453642210148 𝑏4 = 0.0038560700634
72 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
4. Devolver 𝑋.
𝐴𝑓 = {(𝑥, 𝑦) ∈ ℝ2 ∶ 0 ≤ 𝑦 ≤ 𝑓 (𝑥)} .
𝑏
Area de {(𝑥, 𝑦) ∈ ℝ2 ∶ 𝑎 < 𝑥 < 𝑏; 0 ≤ 𝑦 ≤ 𝑓 (𝑥)}
𝑃 (𝑎 < 𝑋 < 𝑏) = = ∫ 𝑓 (𝑥) 𝑑𝑥
Area de 𝐴𝑓 𝑎
Por tanto, si (𝑇 , 𝑌 ) sigue una distribución uniforme en 𝐴𝑐𝑔 , aceptando los va-
lores de (𝑇 , 𝑌 ) que pertenezcan a 𝐴𝑓 (o a 𝐴𝑓 ∗ ) se obtendrán generaciones con
distribución uniforme sobre 𝐴𝑓 (o 𝐴𝑓 ∗ ) y la densidad de la primera componente
𝑇 será 𝑓.
74 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
5.2.1 Algoritmo
Supongamos que 𝑓 es la densidad objetivo y 𝑔 es una densidad auxiliar (fácil de
simular y similar a 𝑓), de forma que existe una constante 𝑐 > 0 tal que:
𝑓 (𝑥) ≤ 𝑐 ⋅ 𝑔 (𝑥) , ∀𝑥 ∈ ℝ.
Ejercicio 5.2.
2𝑒
𝑐opt = √ ≃ 1. 3155.
𝜋
𝑔 (𝑇 ) 2𝑒 𝜋 𝑇2 𝑇2 1
𝑐⋅𝑈 ⋅ = √ 𝑈 √ exp ( − |𝑇 |) = 𝑈 ⋅ exp ( − |𝑇 | + ) ,
𝑓 (𝑇 ) 𝜋 2 2 2 2
# Grafico
curve(c.opt * ddexp(x), xlim = c(-4, 4), lty = 2)
curve(dnorm(x), add = TRUE)
76 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
0.6
0.5
c.opt * ddexp(x)
0.4
0.3
0.2
0.1
0.0
−4 −2 0 2 4
ngen <- 0
system.time(x <- rnormARn(nsim))
##
## Nº de generaciones = 13163
## Nº medio de generaciones = 1.3163
## Proporción de rechazos = 0.2402948
Histogram of x
0.4
0.3
Density
0.2
0.1
0.0
−2 0 2 4
area (𝐴𝑓 ) 1
𝑝= = .
area (𝐴𝑐𝑔 ) 𝑐
Por tanto:
1
𝐸 (𝑁 ) = =𝑐
𝑝
es el número medio de iteraciones del algoritmo (el número medio de pares de
variables (𝑇 , 𝑈 ) que se necesitan generar, y de comparaciones, para obtener una
simulación de la densidad objetivo).
Es obvio, por tanto, que cuanto más cercano a 1 sea el valor de 𝑐 más eficiente
será el algoritmo (el caso de 𝑐 = 1 se correspondería con 𝑔 = 𝑓 y no tendría senti-
do emplear este algoritmo). El principal problema con este método es encontrar
una densidad auxiliar 𝑔 de forma que:
𝑓 (𝑥)
𝑐opt = max .
{𝑥∶𝑔(𝑥)>0} 𝑔 (𝑥)
Ejercicio 5.3.
## $maximum
## [1] -0.999959
##
## $objective
## [1] 1.315489
# NOTA: Cuidado con los límites
# optimize(f=function(x){dnorm(x)/ddexp(x)}, maximum=TRUE, interval=c(0,2))
## [1] 1.315489
d) Aproximar el parámetro óptimo de la densidad auxiliar numéricamente
(normalmente comenzaríamos por este paso).
# Obtención de valores c y lambda óptimos aproximados
fopt <- function(lambda) {
# Obtiene c fijado lambda
optimize(f = function(x){dnorm(x)/ddexp(x,lambda)},
maximum=TRUE, interval=c(0,2))$objective
}
𝐿(x|𝜃)𝜋(𝜃)
𝜋(𝜃|x) =
∫ 𝐿(x|𝜃)𝜋(𝜃)𝑑𝜃
𝑛
siendo 𝐿(x|𝜃) la función de verosimilitud (𝐿(x|𝜃) = ∏ 𝑓(𝑥𝑖 |𝜃) suponiendo
𝑖=1
i.i.d.). Es decir:
𝜋(𝜃|x) ∝ 𝐿(x|𝜃)𝜋(𝜃).
Ejercicio 5.4.
set.seed(54321)
x <- rnorm(n, mean = mu0)
# Función de verosimilitud
lik <- function(mu){prod(dnorm(x, mean = mu))}
# Cota óptima
# Estimación por máxima verosimilitud
emv <- optimize(f = lik, int = range(x), maximum = TRUE)
emv
## $maximum
## [1] 0.7353805
##
## $objective
## [1] 3.303574e-08
c <- emv$objective
## [1] 0.7353958
y por tanto:
c <- lik(mean(x))
c
## [1] 3.303574e-08
Finalmente podríamos emplear el siguiente código para generar simulacio-
nes de la distribución a posteriori mediante aceptación-rechazo a partir de
la distribución de Cauchy:
ngen <- nsim
Y <- rcauchy(nsim)
ind <- (c*runif(nsim) > sapply(Y, lik)) # TRUE si no verifica condición
# Volver a generar si no verifica condición
while (sum(ind)>0){
le <- sum(ind)
ngen <- ngen + le
Y[ind] <- rcauchy(le)
ind[ind] <- (c*runif(le) > sapply(Y[ind], lik)) # TRUE si no verifica condición
}
{ # Número generaciones
5.3. MODIFICACIONES DEL MÉTODO DE ACEPTACIÓN RECHAZO 81
Distribución a posteriori
1.2
1.0
0.8
Density
0.6
0.4
0.2
0.0
𝑠(𝑥) ≤ 𝑓(𝑥).
Algoritmo:
1. Generar 𝑈 ∼ 𝒰 (0, 1) y 𝑇 ∼ 𝑔.
2. Si 𝑐 ⋅ 𝑈 ⋅ 𝑔 (𝑇 ) ≤ 𝑠 (𝑇 ) devolver 𝑋 = 𝑇 ,
en caso contrario
2.a. si 𝑐 ⋅ 𝑈 ⋅ 𝑔 (𝑇 ) ≤ 𝑓 (𝑇 ) devolver 𝑋 = 𝑇 ,
2.b. en caso contrario volver al paso 1.
Cuanto mayor sea el área bajo 𝑠 (𝑥) (más próxima a 1) más efectivo será el
algoritmo.
Se han desarrollado métodos generales para la construcción de las funciones 𝑔 y
𝑠 de forma automática (cada vez que se evalúa la densidad se mejoran las apro-
ximaciones). Estos métodos se basan principalmente en que una transformación
de la densidad objetivo es cóncava o convexa.
𝑘
con ∑𝑗=1𝑝𝑗 = 1, 𝑝𝑗 ≥ 0 y 𝑓𝑗 densidades (sería también válido para funciones
de distribución o variables aleatorias, incluyendo el caso discreto).
Algoritmo:
1. Generar 𝐽 con distribución 𝑃 (𝐽 = 𝑗) = 𝑝𝑗 .
2. Generar 𝑋 ∼ 𝑓𝐽 .
𝜆 −𝜆|𝑥|
𝑓 (𝑥) = 𝑒 , ∀𝑥 ∈ ℝ,
2
se deduce que:
1 1
𝑓 (𝑥) = 𝑓1 (𝑥) + 𝑓2 (𝑥)
2 2
siendo:
𝜆𝑒−𝜆𝑥 si 𝑥 ≥ 0 𝜆𝑒𝜆𝑥 si 𝑥 < 0
𝑓1 (𝑥) = { 𝑓2 (𝑥) = {
0 si𝑥 < 0 0 si𝑥 ≥ 0
Algoritmo:
1. Generar 𝑈 , 𝑉 ∼ 𝒰 (0, 1).
2. Si 𝑈 < 0.5 devolver 𝑋 = − ln (1 − 𝑉 ).
3. En caso contrario devolver 𝑋 = ln 𝑉 .
Observaciones:
5.5. MÉTODOS ESPECÍFICOS PARA LA GENERACIÓN DE ALGUNAS DISTRIBUCIONES NOTABLES85
𝑓(𝑥) = ∫ 𝑔(𝑥|𝑦)ℎ(𝑦)𝑑𝑦
Algoritmo:
1. Generar 𝑌 ∼ ℎ.
2. Generar 𝑋 ∼ 𝑔(⋅|𝑌 ).
Este algoritmo es muy empleado en Inferencia Bayesiana y en la simulación de
algunas variables discretas (como la Binomial Negativa, denominada también
distribución Gamma–Poisson, o la distribución Beta-Binomial), ya que el resul-
tado sería válido cambiando las funciones de densidad 𝑓 y 𝑔 por funciones de
masa de probabilidad.
Simulación de variables
discretas
Considerando como partida una 𝒰 (0, 1), la idea general consiste en asociar a
cada posible valor 𝑥𝑖 de 𝑋 un subintervalo de (0, 1) de logitud igual a la corres-
pondiente probabilidad. Por ejemplo, como ya se mostró en capítulos anteriores,
es habitual emplear código de la forma:
x <- runif(nsim) < p
Esta función del paquete base implementa eficientemente el método “alias” que
describiremos más adelante en la Sección 6.3.
87
88 CAPÍTULO 6. SIMULACIÓN DE VARIABLES DISCRETAS
𝑄 (𝑢) ≤ 𝑥 ⟺ 𝑢 ≤ 𝐹 (𝑥).
4. Devolver 𝑋 = 𝑥𝐼 .
return(X)
}
Ejercicio 6.1.
x <- 0:n
fmp <- dbinom(x, n, p)
ncomp <- 0
system.time( rx <- rfmp(x, fmp, nsim) )
## [1] 5.00322
El valor teórico es n*p = 5.
Número medio de comparaciones:
ncomp/nsim
## [1] 6.00322
# Se verá más adelante que el valor teórico es sum((1:length(x))*fmp)
0.25
0.20
frecuencia relativa
0.15
0.10
0.05
0.00
0 1 2 3 4 5 6 7 8 9 10
valores
## x psim pteor
## 1 0 0.00107 0.0009765625
## 2 1 0.00990 0.0097656250
## 3 2 0.04432 0.0439453125
## 4 3 0.11778 0.1171875000
## 5 4 0.20425 0.2050781250
## 6 5 0.24375 0.2460937500
## 7 6 0.20454 0.2050781250
## 8 7 0.11898 0.1171875000
## 9 8 0.04419 0.0439453125
## 10 9 0.01023 0.0097656250
## 11 10 0.00099 0.0009765625
# Errores
max(abs(res$psim - res$pteor))
## [1] 0.00234375
max(abs(res$psim - res$pteor) / res$pteor)
92 CAPÍTULO 6. SIMULACIÓN DE VARIABLES DISCRETAS
## [1] 0.09568
Nota: . Puede ocurrir que no todos los valores sean generados en la simulación.
En el código anterior si length(x) > length(psim), la sentencia res$pteor
<- fmp gererará un error. Alternativamente se podría emplear por ejemplo:
res <- data.frame(x = x, pteor = fmp, psim = 0)
res.sim <- table(rx)/nsim
index <- match(names(res.sim), x)
res$psim[index] <- res.sim
𝑖 1 2 ⋯ 𝑛 ⋯
𝑃 (ℐ = 𝑖) 𝑝1 𝑝2 ⋯ 𝑝𝑛 ⋯
𝐸 (ℐ) = ∑ 𝑖𝑝𝑖 .
𝑖
Ejercicio 6.2.
ncomp <- 0
6.2. MÉTODO DE LA TABLA GUÍA 93
## [1] 3.08969
sum((1:length(x))*fmp[ind]) # Valor teórico
## [1] 3.083984
𝑗−1 𝑗
𝐼𝑗 = [𝑢𝑗 , 𝑢𝑗+1 ) = [ , ) para 𝑗 = 1, 2, … , 𝑚
𝑚 𝑚
𝑗−1
𝑔𝑗 = 𝑄ℐ (𝑢𝑗 ) = inf {𝑖 ∶ 𝐹𝑖 ≥ 𝑢𝑗 = }
𝑚
𝑗0 = ⌊𝑚𝑈 ⌋ + 1
94 CAPÍTULO 6. SIMULACIÓN DE VARIABLES DISCRETAS
En este caso, puede verse que una cota del número medio de comparaciones es:
𝑛
𝐸 (𝑁 ) ≤ 1 +
𝑚
Ejercicio 6.3.
Diseñar una rutina que permita generar 𝑛𝑠𝑖𝑚 valores de una distribución dis-
creta usando una tabla guía. Repetir el ejercicio anterior empleando esta rutina
con 𝑚 = 𝑛 − 1.
6.3. MÉTODO DE ALIAS 95
𝑥𝑖 con prob. 𝑞𝑖
𝑄(𝑖) = {
𝑥𝑎𝑖 con prob. 1 − 𝑞𝑖
0.25
0.20
frecuencia relativa
0.15
0.10
0.05
0.00
0 1 2 3 4 5 6 7 8 9 10
valores
Ejercicio 6.4.
Diseñar una rutina que permita generar 𝑛𝑠𝑖𝑚 valores de una distribución dis-
creta usando el método de Alias. Repetir el ejercicio anterior empleando esta
rutina.
rfmp.alias <- function(x, prob = 1/length(x), nsim = 1000) {
# Inicializar tablas
a <- numeric(length(x))
q <- prob*length(x)
low <- q < 1
high <- which(!low)
low <- which(low)
while (length(high) && length(low)) {
l <- low[1]
h <- high[1]
a[l] <- h
q[h] <- q[h] - (1 - q[l])
if (q[h] < 1) {
high <- high[-1]
low[1] <- h
} else low <- low[-1]
} # while
# Generar valores
V <- runif(nsim)
i <- floor(runif(nsim)*length(x)) + 1
return( x[ ifelse( V < q[i], i, a[i]) ] )
}
0.25
0.20
frecuencia relativa
0.15
0.10
0.05
0.00
0 1 2 3 4 5 6 7 8 9 10
valores
𝑒−𝜆 𝜆𝑗−1
𝑝𝑗 = 𝑃 (𝑋 = 𝑥𝑗 ) = 𝑃 (𝑋 = 𝑗 − 1) = , 𝑗 = 1, 2, …
(𝑗 − 1)!
6.5. CÁLCULO DIRECTO DE LA FUNCIÓN CUANTIL 99
Ejercicio 6.5.
⎧ 0 si 𝑥 < 0
{ 𝑥+ 1 si 𝑥 ∈ [0, 15 )
𝐹 (𝑥) = ⎨ 2 101
𝑥 + 10 si 𝑥 ∈ [ 15 , 10
9
]
{
⎩ 1 en otro caso
Función de distribución:
fdistr <- function(x) {
ifelse(x < 0, 0,
ifelse(x < 1/5, x/2 + 1/10,
ifelse(x <= 9/10, x + 1/10, 1) ) )
}
# Empleando ifelse se complica un poco más pero el resultado es una función vectorial.
6.6. OTROS MÉTODOS 101
Función de distribución
1.0
0.8
0.6
fdistr(x)
0.4
0.2
0.0
Nota: Esta variable toma los valores 0 y 1/5 con probabilidad 1/10.
a) Diseña un algoritmo basándote en el método de inversión generalizado
para generar observaciones de esta variable.
El algoritmo general es siempre el mismo. Empleando la función cuantil:
el algoritmo sería:
1. Generar 𝑈 ∼ 𝒰 (0, 1)
2. Devolver 𝑋 = 𝑄 (𝑈 )
En este caso concreto:
1. Generar 𝑈 ∼ 𝒰 (0, 1)
1
2. Si 𝑈 < 10 devolver 𝑋 = 0
2 1
3. Si 𝑈 < 10 devolver 𝑋 = 2(𝑈 − 10 )
3 2
4. Si 𝑈 < 10 devolver 𝑋 = 10
1
5. En caso contrario devolver 𝑋 = 𝑈 − 10
102 CAPÍTULO 6. SIMULACIÓN DE VARIABLES DISCRETAS
Ejemplo:
set.seed(1)
nsim <- 10^4
system.time(simx <- rx(nsim))
Histogram of simx
2.5
2.0
1.5
Density
1.0
0.5
0.0
simx
1.0
0.8
0.6
ecdf(simx)(x)
0.4
0.2
0.0
Simulación de
distribuciones
multivariantes
105
106CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
𝑡
Si 𝜇 = (𝜇1 , 𝜇2 , … , 𝜇𝑑 ) es un vector (de medias) y Σ es una matriz 𝑑 × 𝑑 definida
positiva (de varianzas-covarianzas), el vector aleatorio X sigue una distribución
normal multivariante con esos parámetros, X ∼ 𝒩𝑑 (𝜇, Σ), si su función de
densidad es de la forma:
1 1
𝑓(x) = exp (− (x − 𝜇)𝑡 Σ−1 (x − 𝜇)) ,
(2𝜋)𝑛/2 |Σ|1/2 2
donde |Σ| es el determinante de Σ.
Si la matriz de covarianzas es diagonal Σ = 𝑑𝑖𝑎𝑔 (𝜎12 , 𝜎22 , … , 𝜎𝑑2 ), entonces las
componentes 𝑋𝑖 ∼ 𝒩 (𝜇𝑖 , 𝜎𝑖2 ) son independientes y podemos simular el vector
aleatorio de forma trivial, por ejemplo mediante el siguiente algoritmo:
Algoritmo 7.1 (de simulación de normales independientes) 1. Simular
𝑍1 , 𝑍2 , … , 𝑍𝑑 ∼ 𝒩 (0, 1) independientes.
2. Para 𝑖 = 1, 2, … , 𝑑 hacer 𝑋𝑖 = 𝜇𝑖 + 𝜎𝑖 𝑍𝑖 .
Las funciones implementadas en el paquete base de R permiten simular fácil-
mente en el caso independiente ya que admiten vectores como parámetros. Por
ejemplo en el caso bidimensional con 𝑋1 ∼ 𝒩 (0, 1) y 𝑋2 ∼ 𝒩 (−1, 0.52 ):
f1 <- function(x) dnorm(x)
f2 <- function(x) dnorm(x, -1, 0.5)
curve(f1, -3, 3, ylim = c(0, f2(-1)), ylab = "fdp")
curve(f2, add = TRUE, lty = 2)
0.8
0.6
0.4
fdp
0.2
0.0
−3 −2 −1 0 1 2 3
set.seed(1)
rnorm(2, c(0, -1), c(1, 0.5))
## X1 X2
## [1,] -0.6264538 -0.9081783
## [2,] -0.8356286 -0.2023596
## [3,] 0.3295078 -1.4102342
## [4,] 0.4874291 -0.6308376
## [5,] 0.5757814 -1.1526942
1
4 si 𝑥, 𝑦 ∈ [−1, 1]
𝑔 (𝑥, 𝑦) = {
0 en otro caso
1
2𝑑 si 𝑥𝑖 ∈ [−1, 1] , para todo 𝑖 = 1, 2, … , 𝑑
𝑔 (𝑥1 , 𝑥2 , … , 𝑥𝑑 ) = {
0 en otro caso
2𝑑 1 1
𝑈 1 𝑑 (T) ≤ 1 (T) ,
𝑉𝑑 (1) 2𝑑 [−1,1] 𝑉𝑑 (1) 𝐶𝑑
7.3. FACTORIZACIÓN DE LA MATRIZ DE COVARIANZAS 109
o, lo que es lo mismo, 𝑈 1[−1,1]𝑑 (T) ≤ 1𝐶𝑑 (T). Dado que el número aleatorio
𝑈 está en el intervalo (0, 1) y que las funciones indicadoras valen 0 ó 1, esta
condición equivale a que 1[−1,1]𝑑 (T) = 1𝐶𝑑 (T), es decir, a que T ∈ 𝐶𝑑 , es decir,
que se verifique:
𝑇12 + 𝑇22 + ⋯ + 𝑇𝑑2 ≤ 1.
𝑑
Por otra parte, la simulación de 𝑇 ∼ 𝒰 ([−1, 1] ) puede hacerse trivialmente
mediante 𝑇𝑖 ∼ 𝒰 (−1, 1) para cada 𝑖 = 1, 2, … , 𝑑, ya que las componentes son
independientes. Como el valor de 𝑈 es superfluo en este caso, el algoritmo queda:
1. Simular 𝑉1 , 𝑉2 , … , 𝑉𝑑 ∼ 𝒰 (0, 1) independientes.
2. Para 𝑖 = 1, 2, … , 𝑑 hacer 𝑇𝑖 = 2𝑉𝑖 − 1.
3. Si 𝑇12 + 𝑇22 + ⋯ + 𝑇𝑑2 > 1 entonces volver al paso 1.
𝑡
4. Devolver X = (𝑇1 , 𝑇2 , … , 𝑇𝑑 ) .
Ver el Ejercicio 2.1 para el caso de 𝑑 = 2.
Usando las fórmulas del “volumen” de una “esfera” 𝑑-dimensional:
⎧ 𝜋𝑑/2 𝑟𝑑
{ si 𝑑 es par
𝑉𝑑 (𝑟) = ⎨ (𝑑/2)!
⌊𝑑 𝑑
2 ⌋+1 𝜋 ⌊ 2 ⌋ 𝑟𝑑
{ 2 si 𝑑 es impar
⎩ 1 ⋅ 3 ⋅ 5⋯𝑑
puede verse que el número medio de iteraciones del algoritmo, dado por la
𝑑
constante 𝑐 = 𝑉 2(1) , puede llegar a ser enórmemente grande. Así, si 𝑑 = 2
𝑑
se tiene 𝑐 = 1.27, si 𝑑 = 3 se tiene 𝑐 = 1.91, si 𝑑 = 4 entonces 𝑐 = 3.24 y
para 𝑑 = 10 resulta 𝑐 = 401.5 que es un valor que hace que el algoritmo sea
tremendamente lento en dimensión 10. Esto está relacionado con la maldición de
la dimensionalidad (curse of dimensionality), a medida que aumenta el número
de dimensiones el volumen de la “frontera” crece exponencialmente.
𝐶𝑜𝑣(𝐴X) = 𝐴𝐴𝑡 .
Proposición 7.1
Si X ∼ 𝒩𝑑 (𝜇, Σ) y 𝐴 es una matriz 𝑝 × 𝑑, de rango máximo, con 𝑝 ≤ 𝑑,
entonces:
𝐴X ∼ 𝒩𝑝 (𝐴𝜇, 𝐴Σ𝐴𝑡 ) .
Y = 𝜇 + 𝐻Λ1/2 Z ∼ 𝒩𝑑 (𝜇, Σ) .
Y = 𝜇 + 𝐿Z ∼ 𝒩𝑑 (𝜇, Σ) .
nsim <- 20
n <- 100
x <- seq(0, 1, length = n)
# Media
mu <- sin(2*pi*x)
# Covarianzas
x.dist <- as.matrix(dist(x))
x.cov <- exp(-x.dist)
## [,1]
## 1 -0.6264538
## 2 -0.5307633
## 3 -0.5797968
## 4 -0.2844357
## 5 -0.1711797
## 6 -0.2220796
y para simular nsim:
z <- matrix(rnorm(nsim * n), nrow = n)
y <- mu + L %*% z
3
2
1
0
y
−1
−2
−3
Figura 7.2: Realizaciones del proceso funcional del Ejemplo 7.4, obtenidas a
partir de la factorización de Cholesky.
## if (EISPACK)
## stop("'EISPACK' is no longer supported by R", domain = NA)
## eS <- eigen(Sigma, symmetric = TRUE)
## ev <- eS$values
## if (!all(ev >= -tol * abs(ev[1L])))
## stop("'Sigma' is not positive definite")
## X <- matrix(rnorm(p * n), n)
## if (empirical) {
## X <- scale(X, TRUE, FALSE)
## X <- X %*% svd(X, nu = 0)$v
## X <- scale(X, FALSE, TRUE)
## }
## X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*%
## t(X)
## nm <- names(mu)
## if (is.null(nm) && !is.null(dn <- dimnames(Sigma)))
## nm <- dn[[1L]]
## dimnames(X) <- list(nm, NULL)
## if (n == 1)
## drop(X)
## else t(X)
## }
## <bytecode: 0x0000000031a8e548>
## <environment: namespace:MASS>
7.4. MÉTODO DE LAS DISTRIBUCIONES CONDICIONADAS 113
−1
−2
−3
−4
Figura 7.3: Realizaciones del proceso funcional del Ejemplo 7.4, obtenidas em-
pleando la función MASS::mvrnorm.
𝑓1,…,𝑖 (𝑥1 , 𝑥2 , … , 𝑥𝑖 )
𝑓𝑖 (𝑥𝑖 |𝑥1 , 𝑥2 , … , 𝑥𝑖−1 ) = ,
𝑓1,…,𝑖−1 (𝑥1 , 𝑥2 , … , 𝑥𝑖−1 )
+√1−𝑥21
1 2√1 − 𝑥21
𝑓1 (𝑥1 ) = ∫ 𝑑𝑥2 = si 𝑥1 ∈ [−1, 1] ,
−√1−𝑥21 𝜋 𝜋
es decir:
2√
1 − 𝑥21 si 𝑥1 ∈ [−1, 1]
𝑓1 (𝑥1 ) = { 𝜋
0 si 𝑥1 ∉ [−1, 1]
Además:
1
𝑓 (𝑥1 , 𝑥2 ) 𝜋 1
𝑓2 (𝑥2 |𝑥1 ) = = = , si 𝑥2 ∈ [−√1 − 𝑥21 , √1 − 𝑥21 ]
𝑓1 (𝑥1 ) 2√1−𝑥21 2√1 − 𝑥21
𝜋
tendríamos que:
1
𝑓(𝑥1 , 𝑥2 ) =
2𝜋𝜎1 𝜎2 √1 − 𝜌2
1 (𝑥 − 𝜇 )2 (𝑥 − 𝜇 )2 2𝜌(𝑥1 − 𝜇1 )(𝑥2 − 𝜇2 )
exp (− 2
( 1 2 1 + 2 2 2 − ))
2(1 − 𝜌 ) 𝜎1 𝜎2 𝜎1 𝜎2
de donde se deduce (ver e.g. Cao, 2002, p.88; o ecuaciones (7.1) y (7.2) en la
Sección 7.5.1) que:
1 (𝑥1 − 𝜇1 )2
𝑓1 (𝑥1 ) = √ exp (− )
𝜎1 2𝜋 2𝜎12
𝑓 (𝑥1 , 𝑥2 )
𝑓2 (𝑥2 |𝑥1 ) =
𝑓1 (𝑥1 )
2
1 ⎛ (𝑥2 − 𝜇2 − 𝜎𝜎2 𝜌(𝑥1 − 𝜇1 )) ⎞
= ⎜
exp ⎜− 1
⎟
⎟
𝜎2 √2𝜋(1 − 𝜌2 ) 2𝜎22 (1 − 𝜌2 )
⎝ ⎠
𝜎2
Por tanto 𝑋1 ∼ 𝒩 (𝜇1 , 𝜎12 ) y 𝑋2 |𝑋1 ∼ 𝒩 (𝜇2 + 𝜎1 𝜌(𝑋1 − 𝜇1 ), 𝜎22 (1 − 𝜌2 )), y el
algoritmo sería el siguiente:
Algoritmo 7.4 (de simulación de una normal bidimensional) 1. Simular
𝑍1 , 𝑍2 ∼ 𝒩 (0, 1) independientes.
2. Hacer 𝑋1 = 𝜇1 + 𝜎1 𝑍1 .
3. Hacer 𝑋2 = 𝜇2 + 𝜎2 𝜌𝑍1 + 𝜎2 √1 − 𝜌2 𝑍2 .
Este algoritmo es el mismo que obtendríamos con la factorización de la matrix
de covarianzas ya que Σ = 𝐿𝐿𝑡 con:
𝜎12 0
𝐿=( )
𝜎2 𝜌 𝜎2 √1 − 𝜌2
X1 𝜇 Σ Σ12
X=( ) , 𝜇 = ( 1 ) , Σ = ( 11 ),
X2 𝜇2 Σ21 Σ22
suponiendo que X1 se corresponde con los valores observados y X2 con los que
se pretende simular, entonces puede verse (e.g. Ripley, 1987) que la distribución
de X2 |X1 es normal con:
−1
𝐸 (X2 |X1 ) = 𝜇2 + Σ21 Σ11 (X1 − 𝜇1 ) , (7.1)
−1
𝐶𝑜𝑣 (X2 |X1 ) = Σ22 − Σ21 Σ11 Σ12 . (7.2)
3
2
1
0
y
−1
−2
−3
Figura 7.4: Realizaciones condicionales del proceso funcional del Ejemplo 7.7.
Empleando las herramientas del paquete geoR, resulta muy fácil obtener una si-
mulación incondicional en una rejilla en el cuadrado unidad mediante la función
grf:
library(geoR)
n <- 4
set.seed(1)
z <- grf(n, grid="reg", cov.pars=c(1,1))
## x y
## [1,] 0 0
## [2,] 1 0
## [3,] 0 1
## [4,] 1 1
z$data
7.5. SIMULACIÓN CONDICIONAL E INCONDICIONAL 119
0.4
0.2
0.0
# Bucle simulación
nsim <- 1 # 1000
for (i in 1:nsim) {
y <- L %*% rnorm(n)
120 CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
## [,1]
## [1,] -0.62645381
## [2,] -0.05969442
## [3,] -0.98014198
## [4,] 1.09215113
Para generar simulaciones condicionales podemos emplear la función
krige.conv. Por ejemplo, para generar 4 simulaciones en la rejilla regu-
lar 10 × 10 en el cuadrado unidad [0, 1] × [0, 1] condicionadas a los valores
generados en el apartado anterior podríamos emplear el siguiente código:
# Posiciones simulación condicional
nnx <- c(20, 20)
nn <- prod(nnx)
ndata.s <- expand.grid(x = seq(0, 1, len = nnx[1]), y = seq(0, 1, len = nnx[2]))
plot(data.s, type = "p", pch = 20, asp = 1)
points(ndata.s)
1.0
0.8
0.6
y
0.4
0.2
0.0
# Simulación condicional
set.seed(1)
s.out <- output.control(n.predictive = 4)
kc <- krige.conv(z, loc = ndata.s, output = s.out,
krige = krige.control(type.krige="SK", beta = 0, cov.pars = c(1, 1)))
## krige.conv: results will be returned only for prediction locations inside the border
7.5. SIMULACIÓN CONDICIONAL E INCONDICIONAL 121
1.0
0.8
0.8
0.6
0.6
Y Coord
Y Coord
0.4
0.4
0.2
0.2
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
X Coord X Coord
1.0
0.8
0.8
0.6
0.6
Y Coord
Y Coord
0.4
0.4
0.2
0.2
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
X Coord X Coord
par(par.old)
Los valores en las posiciones {(0, 0), (0, 1), (1, 0), (1, 1)} coinciden con los gene-
rados anteriormente.
375
370
Atmospheric concentration of CO2
365
360
355
350
Time
Ejemplos:
• Cópula producto o independiente: Ψ(𝑥) = − ln(𝑥),
1
• Cópula de Clayton: Ψ(𝑥) = 𝛼 (𝑥−𝛼 − 1) ; 𝛼 > 0,
𝛼
• Cópula de Gumbel: Ψ(𝑥) = (− ln(𝑥)) ; 𝛼 ≥ 1
7.6.2 Simulación
Las cópulas pueden facilitar notablemente la simulación de la distribución con-
junta. Si (𝑈 , 𝑉 ) ∼ 𝐶(⋅, ⋅) (marginales uniformes):
Ejercicio 7.1.
1
𝛼 −𝛼
𝐶𝑢−1 (𝑤) ≡ (𝑢−𝛼 (𝑤− 𝛼+1 − 1) + 1) ,
diseñar una rutina que permita generar una muestra de tamaño 𝑛 de esta
distribución.
rcclayton <- function(alpha, n) {
val <- cbind(runif(n), runif(n))
val[, 2] <- (val[, 1]^(-alpha) *
(val[, 2]^(-alpha/(alpha + 1)) - 1) + 1)^(-1/alpha)
return(val)
}
0.4
0.2
0.0
2.0
Density fu
1.5
1.0
nction
0.5
0.0
0.8
0.6
rcu
0.8
nif
0.4
0.6
[2]
# Distribuciones marginales
par.old <- par(mfrow = c(1, 2))
hist(rcunif[,1], freq = FALSE)
abline(h = 1)
hist(rcunif[,2], freq = FALSE)
abline(h = 1)
0.8
0.8
0.6
0.6
Density
Density
0.4
0.4
0.2
0.2
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
rcunif[, 1] rcunif[, 2]
par(par.old)
1.0
0.8
0.6
y[,2]
0.4
0.2
0.0
y[,1]
y[,2]
1.0
0.4
0.8
0.6
0.2
0.4
0.2
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0
y[,1]
4
3
exp2
2
1
0
0 2 4 6 8
exp1
# Distribuciones marginales
par.old <- par(mfrow = c(1, 2))
hist(rcexp[,1], freq = FALSE)
curve(dexp(x,1), add = TRUE)
hist(rcexp[,2], freq = FALSE)
curve(dexp(x,2), add = TRUE)
1.2
1.0
0.6
0.8
Density
Density
0.4
0.6
0.4
0.2
0.2
0.0
0.0
0 2 4 6 8 0 1 2 3 4 5
rcexp[, 1] rcexp[, 2]
128CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
par(par.old)
## [,1] [,2]
## [1,] 1 6
## [2,] 2 7
## [3,] 3 8
## [4,] 4 9
## [5,] 5 10
as.vector(xy)
## [1] 1 2 3 4 5 6 7 8 9 10
Si la variable discreta multidimensional no tiene soporte finito (no se podría
guardar la f.m.p. en una tabla), se podrían emplear métodos de codificación
más avanzados (ver sección 6.3 del libro de R. Cao).
## n particulas calidad
## 1 320 no buena
## 2 14 si buena
## 3 80 no mala
## 4 36 si mala
En lugar de estar en el formato de un conjunto de datos (data.frame) puede
que los datos estén en formato de tabla (table, matrix):
tabla <- xtabs(n ~ calidad + particulas)
tabla
## particulas
## calidad no si
## buena 320 14
## mala 80 36
Lo podemos convertir directamente a data.frame:
as.data.frame(tabla)
## n particulas calidad p
## 1 320 no buena 0.71111111
## 2 14 si buena 0.03111111
## 3 80 no mala 0.17777778
## 4 36 si mala 0.08000000
En formato tabla:
pij <- tabla/sum(tabla)
pij
## particulas
## calidad no si
130 CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
## [1] 1 2 3 4
Con probabilidades:
pz <- df$p
pz
## particulas calidad
## [1,] "no" "buena"
## [2,] "no" "buena"
## [3,] "no" "buena"
## [4,] "si" "mala"
## [5,] "no" "buena"
## [6,] "si" "mala"
Alternativamente, si queremos trabajar con data.frames:
etiquetas <- df[c('particulas', 'calidad')]
rxy <- etiquetas[rz, ]
head(rxy)
## particulas calidad
## 1 no buena
7.7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES DISCRETAS131
## 1.1 no buena
## 1.2 no buena
## 4 si mala
## 1.3 no buena
## 4.1 si mala
# Si se quieren eliminar las etiquetas de las filas:
row.names(rxy) <- NULL
head(rxy)
## particulas calidad
## 1 no buena
## 2 no buena
## 3 no buena
## 4 si mala
## 5 no buena
## 6 si mala
## [,1] [,2]
## [1,] 321 78
## [2,] 15 36
Aunque puede ser preferible emplear directamente rmultinom si se van a generar
muchas:
ntsim <- 1000
rtablas <- rmultinom(ntsim, sum(n), pz)
rtablas[ , 1:5] # Las cinco primeras simulaciones
## particulas
## calidad no si
## buena 320 14
## mala 80 36
## [,1] [,2]
## [1,] 0.6597531 0.08246914
## [2,] 0.2291358 0.02864198
rtablas <- rmultinom(ntsim, sum(n), pind)
rtablas[ , 1:5] # Las cinco primeras simulaciones
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabla
## X-squared = 60.124, df = 1, p-value = 8.907e-15
Histogram of simstat
2.0
1.5
Density
1.0
0.5
0.0
0 2 4 6 8 10 12
simstat
Ejercicio 7.4.
Aplicaciones de la
simulación en Inferencia
Estadística
135
136CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
1 𝑛
𝜇̂ = 𝑋 = ∑𝑋
𝑛 𝑖=1 𝑖
es:
𝜎
𝑋 ∼ 𝑁 (𝜇, √ )
𝑛
Confirmar este resultado mediante simulación, para ello:
a) Crear un conjunto de datos muestras con 500 muestras de tamaño 𝑛 = 10
de una 𝑁 (1, 2). Añadir al conjunto de datos las estimaciones de la media
y desviación típica obtenidas con cada una de las muestras.
Valores iniciales:
set.seed(54321) # Fijar semilla para reproductibilidad
nsim <- 500
nx <- 10
Valores teóricos:
mux <- 1
sdx <- 2
8.1. DISTRIBUCIÓN EN EL MUESTREO 137
Histogram of muestras$mean
0.6
0.5
0.4
Densidad
0.3
0.2
0.1
0.0
−1 0 1 2
Medias
Ejercicio 8.2.
Valores teóricos:
lambda <- 1
muexp <- 1/lambda
sdexp <- muexp
Estimaciones:
muestras2$mean <- rowMeans(muestras2[,1:nx])
muestras2$sd <- apply(muestras2[,1:nx], 1, sd)
Histogram of muestras2$mean
1.4
1.2
1.0
0.8
Densidad
0.6
0.4
0.2
0.0
Medias
## [1] 480
100*ncob/nsim # Proporción de intervalos
## [1] 96
100*(1 - alfa) # Proporción teórica bajo normalidad
## [1] 95
8.2. INTERVALOS DE CONFIANZA 141
1
0
−1
0 10 20 30 40 50
Muestra
detach(tmp)
## [1] 469
100*ncob/nsim # Proporción de intervalos
## [1] 93.8
100*(1 - alfa) # Proporción teórica bajo normalidad
## [1] 95
1.0
0.5
0.0
0 20 40 60 80 100
Muestra
detach(tmp)
1.4
1.2
1.0
0.8
Densidad
0.6
0.4
0.2
0.0
Medias
Ejercicio 8.4.
𝑝(1
̃ − 𝑝)̃ 𝑝(1
̃ − 𝑝)̃
𝑎
𝐼𝐶1−𝛼 (𝑝) = (𝑝̃ − 𝑧1−𝛼/2 √ , 𝑝̃ + 𝑧1−𝛼/2 √ ),
𝑛̃ 𝑛̃
𝑛𝑝 + 𝑎
siendo 𝑛̃ = 𝑛 + 2𝑎, 𝑝̃ = .
𝑛̃
En el caso de 𝑎 = 2 se denomina IC Agresti-Coull.
a) Teniendo en cuenta que la v.a. 𝑋 = 𝑛𝑝̂ ∼ ℬ(𝑛, 𝑝), obtener y represen-
tar gráficamente la cobertura teórica del intervalo de confianza estándar
(𝑎 = 0) de una proporción para una muestra de tamaño 𝑛 = 30, 𝛼 =
0.05 y distintos valores de 𝑝 (p.teor <- seq(1/n, 1 - 1/n, length =
1000)).
Parámetros:
8.2. INTERVALOS DE CONFIANZA 145
n <- 30
alpha <- 0.05
adj <- 0 # (adj <- 2 para Agresti-Coull)
Probabilidades teóricas:
m <- 1000
p.teor <- seq(1/n, 1 - 1/n, length = m)
Posibles resultados:
x <- 0:n
p.est <- (x + adj)/(n + 2 * adj)
ic.err <- qnorm(1 - alpha/2) * sqrt(p.est * (1 - p.est)/(n + 2 * adj))
lcl <- p.est - ic.err
ucl <- p.est + ic.err
Gráfico coberturas:
plot(p.teor, p.cov, type = "l", ylim = c(1 - 4 * alpha, 1))
abline(h = 1 - alpha, lty = 2)
146CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
1.00
0.95
p.cov
0.90
0.85
0.80
p.teor
Parámetros:
n <- 30
alpha <- 0.05
adj <- 2 # Agresti-Coull
# Probabilidades teóricas:
m <- 1000
p.teor <- seq(1/n, 1 - 1/n, length = m)
# Posibles resultados:
x <- 0:n
p.est <- (x + adj)/(n + 2 * adj)
ic.err <- qnorm(1 - alpha/2) * sqrt(p.est * (1 - p.est)/(n + 2 * adj))
lcl <- p.est - ic.err
ucl <- p.est + ic.err
# Recorrer prob. teóricas:
p.cov <- numeric(m)
for (i in 1:m) {
# cobertura de los posibles intervalos
cover <- (p.teor[i] >= lcl) & (p.teor[i] <= ucl)
# prob. de los posibles intervalos
p.rel <- dbinom(x[cover], n, p.teor[i])
# prob. total de cobertura
p.cov[i] <- sum(p.rel)
}
8.2. INTERVALOS DE CONFIANZA 147
# Gráfico coberturas:
plot(p.teor, p.cov, type = "l", ylim = c(1 - 4 * alpha, 1))
abline(h = 1 - alpha, lty = 2)
1.00
0.95
p.cov
0.90
0.85
0.80
p.teor
Parámetros:
n <- 30
alpha <- 0.05
adj <- 2 #' (2 para Agresti-Coull)
set.seed(54321)
nsim <- 500
# Probabilidades teóricas:
m <- 1000
p.teor <- seq(1/n, 1 - 1/n, length = m)
Representar:
plot(p.teor, p.cov, type = "l", ylim = c(1 - 4 * alpha, 1))
abline(h = 1 - alpha, lty = 2)
1.00
0.95
p.cov
0.90
0.85
0.80
p.teor
sx <- 1
nsim <- 1000
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)
Realizar contrastes
for(isim in 1:nsim) {
rx <- rnorm(nx, mx, sx)
tmp <- ks.test(rx, "pnorm", mean(rx), sd(rx))
estadistico[isim] <- tmp$statistic
pvalor[isim] <- tmp$p.value
}
Proporción de rechazos:
{
cat("\nProporción de rechazos al 1% =", mean(pvalor < 0.01), "\n")
cat("Proporción de rechazos al 5% =", mean(pvalor < 0.05), "\n")
cat("Proporción de rechazos al 10% =", mean(pvalor < 0.1), "\n")
}
##
## Proporción de rechazos al 1% = 0
## Proporción de rechazos al 5% = 0
## Proporción de rechazos al 10% = 0.001
Análisis de los p-valores:
hist(pvalor, freq=FALSE)
abline(h=1, lty=2) # curve(dunif(x,0,1), add=TRUE)
Histogram of pvalor
4
3
Density
2
1
0
pvalor
150CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2,
main = 'Tamaño del contraste', ylab = 'Proporción de rechazos',
xlab = 'Nivel de significación')
abline(a=0, b=1, lty=2) # curve(punif(x, 0, 1), add = TRUE)
1.0
0.8
Proporción de rechazos
0.6
0.4
0.2
0.0
Nivel de significación
Valores iniciales:
set.seed(54321)
nx <- 30
mx <- 0
sx <- 1
nsim <- 1000
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)
Realizar contrastes
for(isim in 1:nsim) {
rx <- rnorm(nx, mx, sx)
# tmp <- ks.test(rx, "pnorm", mean(rx), sd(rx))
tmp <- lillie.test(rx)
estadistico[isim] <- tmp$statistic
pvalor[isim] <- tmp$p.value
}
8.3. CONTRASTES DE HIPÓTESIS 151
Proporción de rechazos:
{
cat("\nProporción de rechazos al 1% =", mean(pvalor < 0.01), "\n")
cat("Proporción de rechazos al 5% =", mean(pvalor < 0.05), "\n")
cat("Proporción de rechazos al 10% =", mean(pvalor < 0.1), "\n")
}
##
## Proporción de rechazos al 1% = 0.01
## Proporción de rechazos al 5% = 0.044
## Proporción de rechazos al 10% = 0.089
Histogram of pvalor
1.2
1.0
0.8
Density
0.6
0.4
0.2
0.0
pvalor
# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2, main = 'Tamaño del contraste',
ylab = 'Proporción de rechazos', xlab = 'Nivel de significación')
abline(a=0, b=1, lty=2) # curve(punif(x, 0, 1), add = TRUE)
152CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
1.0
0.8
Proporción de rechazos
0.6
0.4
0.2
0.0
Nivel de significación
Realizar contrastes
for(isim in 1:nsim) {
rx <- rexp(nx, ratex)
tmp <- ks.test(rx, "pexp", 1/mean(rx))
estadistico[isim] <- tmp$statistic
pvalor[isim] <- tmp$p.value
}
Proporción de rechazos:
{
cat("\nProporción de rechazos al 1% =", mean(pvalor < 0.01), "\n")
cat("Proporción de rechazos al 5% =", mean(pvalor < 0.05), "\n")
cat("Proporción de rechazos al 10% =", mean(pvalor < 0.1), "\n")
}
##
## Proporción de rechazos al 1% = 0
## Proporción de rechazos al 5% = 0.004
8.3. CONTRASTES DE HIPÓTESIS 153
Histogram of pvalor
2.0
1.5
Density
1.0
0.5
0.0
pvalor
# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2,
main = 'Tamaño del contraste', ylab = 'Proporción de rechazos',
xlab = 'Nivel de significación')
abline(a=0, b=1, lty=2) # curve(punif(x, 0, 1), add = TRUE)
0.6
0.4
0.2
0.0
Nivel de significación
Simulación:
set.seed(54321)
nx <- 30
ratex <- 1
nsim <- 500
estadistico <- numeric(nsim)
pvalor <- numeric(nsim)
Realizar contrastes
for(isim in 1:nsim) {
rx <- rexp(nx, ratex)
# tmp <- ks.test(rx, "pexp", 1/mean(rx))
tmp <- ks.exp.sim(rx, nsim = 200)
estadistico[isim] <- tmp$statistic
pvalor[isim] <- tmp$p.value
}
8.3. CONTRASTES DE HIPÓTESIS 155
Proporción de rechazos:
{
cat("\nProporción de rechazos al 1% =", mean(pvalor < 0.01), "\n")
cat("Proporción de rechazos al 5% =", mean(pvalor < 0.05), "\n")
cat("Proporción de rechazos al 10% =", mean(pvalor < 0.1), "\n")
}
##
## Proporción de rechazos al 1% = 0.008
## Proporción de rechazos al 5% = 0.058
## Proporción de rechazos al 10% = 0.106
Histogram of pvalor
1.2
1.0
0.8
Density
0.6
0.4
0.2
0.0
pvalor
# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2,
main = 'Tamaño del contraste', ylab = 'Proporción de rechazos',
xlab = 'Nivel de significación')
abline(a=0, b=1, lty=2) # curve(punif(x, 0, 1), add = TRUE)
156CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
1.0
0.8
Proporción de rechazos
0.6
0.4
0.2
0.0
Nivel de significación
for (i in seq_along(shapex)) {
rx <- matrix(rweibull(nx*nsim, shape = shapex[i], scale = 1/ratex), ncol=nx)
preject[i] <- mean( apply(rx, 1, ks.test.p) <= alfa )
preject2[i] <- mean( apply(rx, 1, ks.exp.sim.p) <= alfa )
}
8.4. COMPARACIÓN DE ESTIMADORES 157
plot(shapex, preject, type="l", main = paste("Potencia del contraste ( alfa =", alfa, ")"),
xlab = "shape", ylab = "Proporción de rechazos")
lines(shapex, preject2, lty = 2)
abline(h = alfa, v = 1, lty = 3)
1.0
0.8
Potencia del contraste ( alfa = 0.1 )
Proporción de rechazos
0.6
0.4
0.2
0.0
shape
Función de densidad:
curve(dnorm(x, 0, 1), -3, 12, ylab = 'densidad', lty = 3)
curve(0.95*dnorm(x, 0, 1) + 0.05*dnorm(x, 3, 3), add = TRUE)
0.4
0.3
densidad
0.2
0.1
0.0
0 5 10
str(dat.sim[,1])
Histogram of dat.sim[, 1]
40
30
Frequency
20
10
0
−2 0 2 4 6
dat.sim[, 1]
## [1] 0.1459986
sd(mean.sim)
## [1] 0.1349537
mean(median.sim) # Coincide con el sesgo (media teórica es 0)
## [1] 0.04453509
sd(median.sim)
## [1] 0.1300611
Sesgo:
boxplot(mean.sim-xmed, median.sim-xmed,
names=c("Media","Mediana"), ylab="Sesgo")
abline(h = 0, lty = 2)
160CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
0.6
0.4
0.2
Sesgo
0.0
−0.2
−0.4
Media Mediana
Error cuadrático:
boxplot((mean.sim-xmed)^2, (median.sim-xmed)^2,
names=c("Media","Mediana"), ylab="Error cuadrático")
0.3
Error cuadrático
0.2
0.1
0.0
Media Mediana
# SE mediana
summary((median.sim-xmed)^2)
0.4
0.2
0.0
# Muestra bootstrap
xboot <- sample(data, replace=TRUE)
𝑡
Si x = (𝑥1 , 𝑥2 , ⋯ , 𝑥𝑛 ) es una muestra i.i.d. 𝐹𝜃 y 𝜃 ̂ = 𝑇 (x) es un estimador de
𝜃. Para 𝑏 = 1, … , 𝐵 ∶
𝑡
• x∗𝑏 = (𝑥∗𝑏1 , 𝑥∗𝑏2 , ⋯ , 𝑥∗𝑏𝑛 ) muestra bootstrap obtenida mediante muestreo
con reemplazamiento de {𝑥1 , 𝑥2 , ⋯ , 𝑥𝑛 }.
Para información adicional sobre bootstrap ver p.e.: Davison, A.C. and Hinkley,
D.V. (1997). Bootstrap Methods and Their Application. Cambridge University
Press
8.5.3 Ejemplo
Como ejemplo consideraremos una muestra de tamaño 𝑛 = 100 de una normal
estándar para tratar de aproximar el sesgo y error estándar de la media y la
mediana mediante bootstrap.
# dat <- dat.sim[, 1]
set.seed(54321)
ndat <- 100
datmed <- 0
datsd <- 1
dat <- rnorm(ndat, mean=datmed, sd=datsd)
Es habitual tomar nboot + 1 entero múltiplo de 100 e.g. nboot = 999 ó 1999
(facilita el cálculo de percentiles, orden (nboot + 1)*p)
El valor del estadístico en la muestra es:
stat.dat <- mean(dat)
## [1] -0.09274131
Bootstrap percentil:
hist(stat.boot, freq=FALSE, ylim = c(0,4))
abline(v=mean.boot, lwd=2)
# abline(v=stat.dat)
# Distribución poblacional
curve(dnorm(x, datmed, datsd/sqrt(ndat)), lty=2, add=TRUE)
abline(v=datmed, lwd=2, lty=2)
Histogram of stat.boot
4
3
Density
2
1
0
stat.boot
Bootstrap natural/básico:
hist(stat.boot-stat.dat, freq=FALSE, ylim = c(0,4))
abline(v=mean.boot-stat.dat, lwd=2)
# Distribución poblacional
# Distribución teórica de stat.dat - stat.teor
curve(dnorm(x, 0, datsd/sqrt(ndat)), lty=2, add=TRUE)
abline(v=0, lwd=2, lty=2)
8.5. REMUESTREO BOOTSTRAP 165
4
3
Density
2
1
0
stat.boot − stat.dat
## [1] 0.0006390906
# error estándar
sd(stat.boot)
## [1] 0.1044501
# error estándar teórico
datsd/sqrt(ndat)
## [1] 0.1
Versión “optimizada” para R:
boot.strap <- function(dat, nboot=1000, statistic=mean)
{
ndat <- length(dat)
dat.boot <- sample(dat, ndat*nboot, replace=T)
dat.boot <- matrix(dat.boot, ncol=nboot, nrow=ndat)
stat.boot <- apply(dat.boot, 2, statistic)
}
# max(dat)
}
set.seed(1)
stat.dat <- fstatistic(dat)
stat.boot <- boot.strap(dat, nboot, fstatistic)
Función estadístico:
boot.f <- function(data, indices){
# data[indices] será la muestra bootstrap
mean(data[indices])
}
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = dat, statistic = boot.f, R = nboot)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.0933804 0.0006390906 0.1049474
names(stat.boot)
8.5.5 Gráficos
hist(stat.boot$t, freq=FALSE)
Histogram of stat.boot$t
4
3
Density
2
1
0
stat.boot$t
plot(stat.boot)
Histogram of t
0.2
4
0.1
3
0.0
Density
−0.1
t*
2
−0.2
1
−0.3
−0.4
0
jack.after.boot(stat.boot)
168CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
0.2
* * * * * *
*** * **** *** * *** * *********** **** ** ******************** *** * * ** *
* ** * ** * * *
0.1
* ***** * ** * * * ** * *
0.0
** * * ** ** ***** * *** * ***************** ********************** *** ** ******** *** ** * ** ** * * *
*
−0.1
* ** * *
** * * ****** *** * *** ** ************** ** ** * * * * * ** ** * *
***** ******* **** ** *** * **** * *** ** *
*** * ****** *** ** ** **************** **
************** ***** *** * * ***** * ** ** * ** ** *
* *
* * * * * * * *
** * * ***** *** * *** ****************** *********************** ** * ****** * *** ** * ** ** *
−0.2
** *
95 58 96 80 91 29 55 34 84 83 3697 7022 100 38 50 42 75 21
33 30 63 89 47 27 31 11 46 88 37 1 28 51 16 9 45 67 73 62
15 66 99 4 69 65 68 78 24 92 821498 81 48 18 20 86 90 52
−0.3
94 25 60 79 40 5735 6 59 76 43 7787 44 74 26 23 13 17 72
7 85 3 53 41 10 2 32 39 5
64 61 54 93 71 49 12 56 19 8
−2 −1 0 1 2
Ejercicio 8.7.
set.seed(1)
boot(dat, boot.f, nboot, trim=0.2)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = dat, statistic = boot.f, R = nboot, trim = 0.2)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.1448019 0.006095294 0.119021
Integracion y Optimizacion
Montecarlo
1 𝑛
𝐼 = 𝐸 (ℎ (𝒰 (0, 1))) ≈ ∑ ℎ (𝑥𝑖 )
𝑛 𝑖=1
171
172 CAPÍTULO 9. INTEGRACION Y OPTIMIZACION MONTECARLO
𝑏 𝑏
1
𝐼 = ∫ ℎ (𝑥) 𝑑𝑥 = (𝑏 − 𝑎) ∫ ℎ (𝑥) 𝑑𝑥 = (𝑏 − 𝑎)𝐸 (ℎ (𝒰 (𝑎, 𝑏))) .
𝑎 𝑎 (𝑏 − 𝑎)
1 𝑛
𝐼≈ ∑ ℎ (𝑥𝑖 ) (𝑏 − 𝑎)
𝑛 𝑖=1
9.1.1 Ejercicio 1
Crear una función que implemente la integración Monte Carlo clásica para apro-
ximar integrales del tipo:
𝑏
𝐼 = ∫ ℎ (𝑥) 𝑑𝑥.
𝑎
Función a integrar:
fun <- function(x) ifelse((x > 0) & (x < 1), 4 * x^3, 0)
# return(4 * x^3)
curve(fun, 0, 1)
abline(h=0,lty=2)
abline(v=c(0,1),lty=2)
9.1. INTEGRACIÓN MONTE CARLO (CLÁSICA) 173
4
3
fun(x)
2
1
0
set.seed(54321)
mc.integral0(fun, 0, 1, 20)
## [1] 0.8051635
mc.integral0(fun, 0, 1, 100)
## [1] 1.134583
mc.integral0(fun, 0, 1, 100)
## [1] 0.935731
La función mc.integral0 no es adecuada para analizar la convergencia de la
aproximación por simulación. Una alternativa más eficiente para representar
gráficamente la convergencia:
mc.integral <- function(fun, a, b, n, plot = TRUE) {
fx <- sapply(runif(n, a, b), fun) * (b - a)
if (plot) {
estint <- cumsum(fx)/(1:n)
esterr <- sqrt(cumsum((fx - estint)^2))/(1:n)
plot(estint, ylab = "Media y rango de error", type = "l", lwd = 2, ylim = mean(fx) +
2 * c(-esterr[1], esterr[1]), xlab = "Iteraciones")
lines(estint + 2 * esterr, col = "darkgray", lwd = 2)
lines(estint - 2 * esterr, col = "darkgray", lwd = 2)
valor <- estint[n]
abline(h = valor)
return(list(valor = valor, error = 2 * esterr[n]))
} else return(list(valor = mean(fx), error = 2 * sd(fx)/sqrt(n)))
174 CAPÍTULO 9. INTEGRACION Y OPTIMIZACION MONTECARLO
set.seed(54321)
mc.integral(fun, 0, 1, 5000)
## $valor
## [1] 0.9939051
##
## $error
## [1] 0.03206515
abline(h = 1, lty = 2)
1.4
1.2
Media y rango de error
1.0
0.8
0.6
Iteraciones
## $valor
## [1] 0.9939051
##
## $error
## [1] 0.03209283
1 𝑛
𝜃≈ ∑ ℎ (𝑥𝑖 )
𝑛 𝑖=1
2
1
0
Los pasos serían simular x con densidad fun y aproximar la integral por
mean(h(x)):
x <- rfun(nsim) # Habría que implementarlo (por ejemplo por el método de inversión...)
res <- mean(h(x)) # Aproximación por simulación
error <- 2*sd(h(x))/sqrt(nsim)
9.1.3 Ejercicio 2
Aproximar:
∞
1 𝑥2
𝜙(𝑡) = ∫ √ 𝑒− 2 𝑑𝑥,
𝑡 2𝜋
para 𝑡 = 4.5, empleando integración Monte Carlo (aproximación tradicional con
dominio infinito).
# h <- function(x) x > 4.5
# f <- function(x) dnorm(x)
set.seed(54321)
nsim <- 10^3
x <- rnorm(nsim)
mean(x > 4.5) # mean(h(x))
## [1] 0
pnorm(-4.5) # valor teórico P(X > 4.5)
## [1] 3.397673e-06
De esta forma es dificil que se generen valores (en este caso ninguno) en la región
que interesaría para la aproximación de la integral:
any(x > 4.5)
## [1] FALSE
Sería preferible generar más valores donde hay mayor “área”…
puede ser preferible generar observaciones de una densidad 𝑔 que tenga una
forma similar al producto ℎ𝑓.
Si 𝑌 ∼ 𝑔:
ℎ (𝑥) 𝑓(𝑥)
𝜃 = ∫ ℎ (𝑥) 𝑓(𝑥)𝑑𝑥 = ∫ 𝑔(𝑥)𝑑𝑥 = 𝐸 (𝑞 (𝑌 )) .
𝑔(𝑥)
ℎ(𝑥)𝑓(𝑥)
siendo 𝑞 (𝑥) = 𝑔(𝑥) .
Si 𝑦1 , 𝑦2 , … , 𝑦𝑛 i.i.d. 𝑌 ∼ 𝑔:
1 𝑛 1 𝑛
𝜃≈ ∑ 𝑞 (𝑦𝑖 ) = ∑ 𝑤(𝑦𝑖 )ℎ (𝑦𝑖 ) = 𝜃𝑔̂
𝑛 𝑖=1 𝑛 𝑖=1
9.2. MUESTREO POR IMPORTANCIA 177
𝑓(𝑥)
con 𝑤 (𝑥) = 𝑔(𝑥) .
La idea básica es que si la densidad 𝑔 tiene colas más pesadas que la densidad 𝑓
con mayor facilidad puede dar lugar a “simulaciones” con varianza finita (podría
emplearse en casos en los que no existe 𝐸 (ℎ2 (𝑋)); ver Sección Convergencia en
el Tema 4 de Analisis de resultados).
La distribución de los pesos 𝑤(𝑦𝑖 ) debería ser homogénea para evitar datos
influyentes (inestabilidad).
9.2.1 Ejercicio 3
Aproximar la integral del ejercicio anterior empleando muestreo por importan-
cia considerando como densidad auxiliar una exponencial de parámetro 𝜆 = 1
truncada en 𝑡:
1.5e−05
dnorm(x) y dexp(x−4.5)*k
1.0e−05
5.0e−06
0.0e+00
## [1] 3.243321e-06
pnorm(-4.5) # valor teórico
## [1] 3.397673e-06
plot(cumsum(w)/1:nsim, type = "l", ylab = "Aproximación", xlab = "Iteraciones")
abline(h = pnorm(-4.5), lty = 2)
9.2. MUESTREO POR IMPORTANCIA 179
4e−06
3e−06
Aproximación
2e−06
1e−06
Iteraciones
## [1] 1.386216e-07
## [1] 0
sqrt(est * (1 - est)/nsim)
## [1] 0
9.2.2 Ejercicio 4
Aproximar 𝑃 (2 < 𝑋 < 6) siendo 𝑋 ∼ 𝐶𝑎𝑢𝑐ℎ𝑦(0, 1) empleando muestreo por
importancia y considerando como densidad auxiliar la normal estandar 𝑌 ∼
𝑁 (0, 1). Representar gráficamente la aproximación y estudiar los pesos 𝑤(𝑦𝑖 ).
Nota: En este caso van a aparecer problemas (la densidad auxiliar debería tener
colas más pesadas que la densidad objetivo; sería adecuado si intercambiaramos
las distribuciones objetivo y auxiliar, como en el ejercicio siguiente).
## [1] 0.09929348
pcauchy(6) - pcauchy(2) # Valor teórico
## [1] 0.09501516
Si se estudia la convergencia:
plot(cumsum(w * (y > 2) * (y < 6))/1:nsim, type = "l", ylab = "Aproximación", xlab = "I
abline(h = pcauchy(6) - pcauchy(2), lty = 2)
0.30
0.25
0.20
Aproximación
0.15
0.10
0.05
0.00
Iteraciones
La distribución de los pesos debería ser homogénea. Por ejemplo, si los reesca-
lamos para que su suma sea el número de valores generados, deberían tomar
valores en torno a uno:
boxplot(nsim * w/sum(w))
9.2. MUESTREO POR IMPORTANCIA 181
1200
1000
800
600
400
200
0
9.2.4 Ejercicio 5
Generar 1000 simulaciones de una distribución (aprox.) 𝑁 (0, 1) mediante re-
muestreo del muestreo por importancia de 105 valores de una 𝐶𝑎𝑢𝑐ℎ𝑦(0, 1).
Se trata de simular una normal a partir de una Cauchy (Sampling Importance
Resampling). En este caso f(y) = dnorm(y) y g(y) = dcauchy(y), al revés
del ejercicio anterior…
# Densidad objetivo
# f <- dnorm # f <- function(x) ....
Histogram of rx
0.3
Density
0.2
0.1
0.0
−4 −2 0 2
rx
La idea original consiste en buscar los ceros de su primera derivada (o del gra-
diente) empleando una aproximación iterativa:
9.3.1 Ejercicio 6
La mixtura de distribuciones normales:
1 3
𝑁 (𝜇1 , 1) + 𝑁 (𝜇2 , 1),
4 4
set.seed(12345)
p.sim <- rbinom(nsim, 1, 0.25)
data <- rnorm(nsim, mu1*p.sim + mu2*(1-p.sim), sd1*p.sim + sd2*(1-p.sim))
Histogram of data
0.25
0.20
Density
0.15
0.10
0.05
0.00
−2 0 2 4 6
data
Si queremos capturar los valores en los que se evalúa esta función, podemos
proceder de forma similar a como se describe en el capítulo Function operators
del libro “Advanced R” de Hadley Wickham: Behavioural FOs leave the inputs
and outputs of a function unchanged, but add some extra behaviour.
tee <- function(f) {
function(...) {
input <- if(nargs() == 1) c(...) else list(...)
output <- f(...)
# Hacer algo ...
# ... con output e input
return(output)
}
}
En este caso queremos representar los puntos en los que el algoritmo de op-
timización evalúa la función objetivo (especialmente como evoluciona el valor
óptimo)
tee.optim2d <- function(f) {
best.f <- Inf # Suponemos que se va a minimizar (opción por defecto)
best.par <- c(NA, NA)
function(...) {
9.3. OPTIMIZACIÓN MONTE CARLO 185
5
0 0 −550
−8
00 −70 −65 0
−1
−1350
−60
20
−1250
0
−1
−500 10 −1150
0
−1000 −1050
4
−450
−950
−900
−850
−800
−750
−700
−650
3
−600
2
−400 −400
µ2
−450
1 −500
−550
−600
−650
−700
−750
0
−800
−850
−900
−950
−1000
−1050
−1150 −1100 50
−6
−1
−1200
00
−1350 −1300 −7
−1
25
−1500 −14
0
5 0 00
−16
50 −8
50 00
−6 −7
−2
−550
−2 −1 0 1 2 3 4 5
µ1
• Temple simulado.
• Algoritmos genéticos.
• Montecarlo EM.
• …
9.4. TEMPLE SIMULADO 187
9.4.1 Algoritmo:
temp <- TEMP.INIT
place <- INIT.PLACEMENT()
cost.place <- COST(place)
while(temp < TEMP.FINAL) {
while(LOOP.CRITERION()) {
place.new <- PERTURB(place, temp)
cost.new <- COST(place.new)
cost.inc <- cost.new - cost.place
temp <- SCHEDULE(temp)
if ((cost.inc < 0) || (runif(1) > exp(-(cost.inc/temp)))) break
}
place <- place.new
cost.place <- cost.new
# temp <- SCHEDULE(temp)
}
COST <- function(place, ...) {...}
188 CAPÍTULO 9. INTEGRACION Y OPTIMIZACION MONTECARLO
Este algoritmo se puede ver como una adaptación del método de Metropolis-
Hastings que veremos más adelante (Tema 12 Introducción a los métodos de
cadenas de Markov Monte Carlo).
9.4.2 Ejercicio 7
set.seed(1)
for (j in 1:nstarts){
# Normalmente llamaríamos a optim(start, like, method = "SANN")
# Note that the "SANN" method depends critically on the settings of the control param
# For "SANN" maxit gives the total number of function evaluations: there is no other
res <- optim(starts[j, ], tee.optim2d(like), method = "SANN", control = list(temp = 1
points(res$par[1],res$par[2], pch = 19)
cat('par = ', res$par, 'value =', res$value, '\n')
}
9.4. TEMPLE SIMULADO 189
5
0 0 0 −550
−80 −70 −65 0 −13
−60 50
−11
−500 50 −1200
−9
50 −1000 −1050
4
−450
−900
−800 −850
−750
−700
−650
3
−600
2
−400 −400
µ2
−450
1
−500
−550
−600
−650
−700
0
−750
−800
−850
−900
−950
−1000
−1100 −10
50
−1
−1150
−1300 −1250 00
−7
−1
20
−1500 −1
0
40
−165 0 00
0 00 50 −8
−2
−550 −6 −6
−2 −1 0 1 2 3 4 5
µ1
set.seed(1)
for (j in 1:nstarts) {
sar <- SA(like, starts[j, ])
lines(sar$the[, 1], sar$the[, 2], lwd = 2, col = "blue")
points(sar$the[sar$it, 1], sar$the[sar$it, 2], pch = 19)
}
5
0 0 0 −550
−80 −70 −65 0 −13
−60 50
−11
−500 50 −1200
−9
50 −1000 −1050
4
−450
−900
−800 −850
−750
−700
−650
3
−600
2
−400 −400
µ2
−450
1
−500
−550
−600
−650
−700
0
−750
−800
−850
−900
−950
−1000
−1100 −10
50
−1
−1150
−1300 −1250 00
−7
−1
20
−1500 −1
0
4 00
−165 00
0 60
0 50 −8
−2
−550 − −6
−2 −1 0 1 2 3 4 5
µ1
9.5.1 Ejercicio 8
Repetir el ejercicio anterior empleando la función DEOptim.
Optimización con algoritmo genético implementado en DEoptim:
require(DEoptim)
5
0 0 0 −550
−80 −70 −65 0 −13
−60 50
−11
−500 50 −1200
−9
50 −1000 −1050
4
−450
−900
−800 −850
−750
−700
−650
3
−600
2
−400 −400
µ2
1 −450
−500
−550
−600
−650
−700
0
−750
−800
−850
−900
−950
−1000
−1100 −10
50
−1
−1150
−1300 −1250 00
−7
−1
20
−1500 −1
0
4 00
−165 00
0 00 50 −8
−2
−550 −6 −6
−2 −1 0 1 2 3 4 5
µ1
Capítulo 10
Técnicas de reducción de la
varianza
• Variables antitéticas.
• Muestreo estratificado.
• Variables de control.
• Métodos de remuestreo.
• Condicionamiento.
• …
193
194 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
𝜃 = 𝐸 (𝑍)
𝑋+𝑌 1
𝑉 𝑎𝑟 ( ) = (𝑉 𝑎𝑟 (𝑋) + 𝑉 𝑎𝑟 (𝑌 ) + 2𝐶𝑜𝑣 (𝑋, 𝑌 ))
2 4
𝜎2 1
= + 𝐶𝑜𝑣 (𝑋, 𝑌 )
2𝑛 2𝑛
2
𝜎
= (1 + 𝜌 (𝑋, 𝑌 )) ,
2𝑛
que el equivalente a una muestra unidimensional independiente con el mismo
número de observaciones 2𝑛 (con una reducción del −100𝜌 (𝑋, 𝑌 ) %).
𝑌 = 2𝜇 − 𝑋
Crear una función que implemente la técnica de variables antitéticas para apro-
ximar integrales del tipo:
𝑏
𝐼 = ∫ ℎ (𝑥) 𝑑𝑥.
𝑎
2
1
0
## [1] 3.194528
Para la aproximación por integración Monte Carlo podemos emplear la función
del capítulo anterior:
196 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
set.seed(54321)
res <- mc.integral(ftn, a, b, 500)
abline(h = teor)
4.5
4.0
Media y rango de error
3.5
3.0
2.5
2.0
Iteraciones
res
## $valor
## [1] 3.184612
##
## $error
## [1] 0.1619886
set.seed(54321)
res <- mc.integrala(ftn, a, b, 500)
4.5
4.0
Media y rango de error
3.5
3.0
2.5
2.0
Iteraciones
res
## $valor
## [1] 3.222366
##
## $error
198 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
## [1] 0.1641059
## [1] -1.307067
set.seed(54321)
res <- mc.integrala2(ftn, a, b, 500)
res
## $valor
## [1] 3.222366
##
## $error
## [1] 0.05700069
## [1] 64.81191
En este caso puede verse que la reducción teórica de la varianza es del 96.7%
10.3. ESTRATIFICACIÓN 199
10.3 Estratificación
Si se divide la población en estratos y se genera en cada uno un número de
observaciones proporcional a su tamaño (a la probabilidad de cada uno) nos
aseguramos de que se cubre el dominio de interés y se puede acelerar la conver-
gencia.
• Por ejemplo, para generar una muestra de tamaño 𝑛 de una 𝒰 (0, 1), se
pueden generar 𝑙 = 𝑛𝑘 observaciones (1 ≤ 𝑘 ≤ 𝑛) de la forma:
(𝑗 − 1) 𝑗
𝑈𝑗1 , … , 𝑈𝑗𝑙 ∼ 𝒰 ( , ) para 𝑗 = 1, ..., 𝑘.
𝑘 𝑘
Si en el número de obsevaciones se tiene en cuenta la variabilidad en el estrato
se puede obtener una reducción significativa de la varianza.
1. Para 𝑖 = 1, 2, … , 10:
2. Generar 𝑈𝑖 :
2a. Generar 𝑈 ∼ 𝑈 (0, 1).
2b. Si 𝑖 ≤ 4 hacer 𝑈𝑖 = 0.4 ⋅ 𝑈 + 0.6.
2c. Si 4 < 𝑖 ≤ 9 hacer 𝑈𝑖 = 0.5 ⋅ 𝑈 + 0.1.
2d. Si 𝑖 = 10 hacer 𝑈𝑖 = 0.1 ⋅ 𝑈 .
3. Devolver 𝑋𝑖 = − ln 𝑈𝑖 .
No es difícil probar que:
• 𝑉 𝑎𝑟 (𝑋𝑖 ) = 0.0214644 si 𝑖 = 1, 2, 3, 4,
• 𝑉 𝑎𝑟 (𝑋𝑖 ) = 0.229504 si 𝑖 = 5, 6, 7, 8, 9 y
• 𝑉 𝑎𝑟 (𝑋10 ) = 1.
Como consecuencia:
1 10
𝑉 𝑎𝑟 (𝑋) = ∑ 𝑉 𝑎𝑟 (𝑋𝑖 ) = 0.022338
102 𝑖=1
que es bastante menor que 1 (la varianza en el caso de muestreo aleatorio simple
no estratificado).
set.seed(54321)
res <- mc.integral(ftn, a, b, 500)
abline(h = teor)
10.3. ESTRATIFICACIÓN 201
4.5
4.0
Media y rango de error
3.5
3.0
2.5
2.0
Iteraciones
res
## $valor
## [1] 3.184612
##
## $error
## [1] 0.1619886
set.seed(54321)
mc.integrale(ftn, a, b, 500, 50)
## $valor
## [1] 3.193338
##
## $error
## [1] 0.1597952
set.seed(54321)
mc.integrale(ftn, a, b, 500, 100)
## $valor
## [1] 3.193927
##
## $error
## [1] 0.1599089
De esta forma no se tiene en cuenta la variabilidad en el estrato. El tamaño de
las submuestras debería incrementarse hacia el extremo superior.
Ejercicio 10.3.
202 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
𝐶𝑜𝑣(𝑋, 𝑌 )
𝛼∗ = − ,
𝑉 𝑎𝑟(𝑌 )
con 𝑉 𝑎𝑟(𝑋 ∗ ) = 𝑉 𝑎𝑟(𝑋) (1 − 𝜌2 (𝑋, 𝑌 )) (lo que supone una reducción del
100𝜌2 (𝑋, 𝑌 ) %).
En la práctica normalmente 𝛼∗ no es conocida. Para estimarlo se puede realizar
ajuste lineal de 𝑋 sobre 𝑌 (a partir de los datos simulados 𝑋𝑖 e 𝑌𝑖 , 1 ≤ 𝑖 ≤ 𝑛):
𝑆
• Si 𝑥̂ = 𝛽0̂ + 𝛽1̂ 𝑦 es la recta ajustada, con 𝛽1̂ = 𝑋𝑌 y 𝛽0̂ = 𝑋 − 𝛽1̂ 𝑌 , la
𝑆𝑌2
estimación sería:
𝛼∗̂ = −𝛽 ̂
1
∗
• Si 𝜇𝑌 = 0 ⇒ 𝜃 ̂ = 𝑋 = 𝛽0̂ .
## [1] 3.194528
10.4. VARIABLES DE CONTROL 203
## [1] 3.182118
4
3
2
1
## (Intercept)
## 3.204933
## [1] 3.204933
100*(var(expu)-var(expuc))/var(expu)
## [1] 93.91555
𝐸 (𝑋) − 𝐸 (𝑌 ) = 𝐸 (𝑋 − 𝑌 ) .
1 𝑛
𝑋−𝑌 = ∑ (𝑋 − 𝑌𝑖 )
𝑛 𝑖=1 𝑖
1 𝑁 1
𝑉 𝑎𝑟 (𝑋 − 𝑌 ) = ∑ 𝑉 𝑎𝑟 (𝑋𝑖 − 𝑌𝑖 ) = 𝑉 𝑎𝑟 (𝑋𝑖 − 𝑌𝑖 )
𝑛2 𝑖=1 𝑛
1
= (𝑉 𝑎𝑟 (𝑋𝑖 ) + 𝑉 𝑎𝑟 (𝑌𝑖 ) − 2𝐶𝑜𝑣 (𝑋𝑖 , 𝑌𝑖 ))
𝑛
1
≤ (𝑉 𝑎𝑟 (𝑋𝑖 ) + 𝑉 𝑎𝑟 (𝑌𝑖 ))
𝑛
## [1] 1.97439
# Aprox da precisión
var(x)
## [1] 3.669456
MC con variables antitéticas:
# xa <-
# mean(xa) # Aprox por MC da media (valor teor 1/lambda = 2)
# var(xa) # Aprox da precisión supoñendo independencia
# corr <- cor(x1,x2)
# var(xa)*(1 + corr) # Estimación varianza supoñendo dependencia
Bibliografía básica
Cao, R. (2002). Introducción a la simulación y a la teoría de colas. NetBiblo.
Gentle, J.E. (2003). Random number generation and Monte Carlo methods.
Springer‐Verlag.
Jones, O. et al. (2009). Introduction to Scientific Programming and Simulation
Using R. CRC.
Ripley, B.D. (1987). Stochastic Simulation. John Wiley & Sons.
Robert, C.P. y G. Casella (2010). Introducing Monte Carlo Methods with R.
Springer.
Ross, S.M. (1999).Simulación. Prentice Hall.
Suess, E.A. y Trumbo, B.E. (2010). Introduction to probability simulation and
Gibbs sampling with R. Springer.
Bibliografía complementaria
Libros
Azarang, M. R. y García Dunna, E. (1996). Simulación y análisis de modelos
estocásticos. McGraw-Hill.
Bratley, P., Fox, B.L. y Schrage L.E. (1987). A guide to simulation. Springer-
Verlag.
Devroye, L. (1986). Non-uniform random variate generation. Springer-Verlag.
Evans, M. y Swartz, T. (2000). Approximating integrals via Monte Carlo and
determinstic methods. Oxford University Press.
Gentle, J.E. (1998). Random number generation and Monte Carlo methods.
Springer-Verlag.
207
208 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
Artículos
Demirhan, H. y Bitirim, N. (2016). CryptRndTest: an R package for testing the
cryptographic randomness. The R Journal, 8(1), 233-247.
Downham, D.Y. (1970). Algorithm AS 29: The runs up and down test. Journal
of the Royal Statistical Society. Series C (Applied Statistics), 19(2), 190-192.
L’Ecuyer, P. (1999). Good parameters and implementations for combined mul-
tiple recursive random number generators. Operations Research, 47, 159–164.
L’Ecuyer, P. y Simard, R. (2007). TestU01: A C library for empirical testing
of random number generators. ACM Transactions on Mathematical Software
(TOMS), 33(4), 1-40.
Marsaglia, G. y Tsang, W.W. (2002). Some difficult-to-pass tests of randomness.
Journal of Statistical Software, 7(3), 1-9.
Marsaglia, G., Zaman, A. y Tsang, W.W. (1990). Toward a universal random
number generator. Stat. Prob. Lett., 9(1), 35-39.
Matsumoto, M. y Nishimura, T. (1998). Mersenne Twister: A 623-dimensionally
equidistributed uniform pseudo-random number generator, ACM Transactions
on Modeling and Computer Simulation, 8, 3–30.
Park, S.K. y Miller , K.W. (1988). Random number generators: good ones are
hard to find. Communications of the ACM, 31(10), 1192-1201.
Park, S.K., Miller, K.W. y Stockmeyer, P.K. (1993). Technical correspondence.
Communications of the ACM, 36(7), 108-110.
Apéndice A
Enlaces
209
210 APÉNDICE A. ENLACES
Bondad de Ajuste y
Aleatoriedad
En los métodos clásicos de inferencia estadística es habitual asumir que los valo-
res observados 𝑋1 , … , 𝑋𝑛 (o los errores de un modelo) consituyen una muestra
aleatoria simple de una variable aleatoria 𝑋. Se están asumiendo por tanto
dos hipótesis estructurales: la independencia (aleatoriedad) y la homogeneidad
(misma distribución) de las observaciones (o de los errores). Adicionalmente,
en inferencia paramétrica se supone que la distribución se ajusta a un modelo
paramétrico específico 𝐹𝜃 (𝑥), siendo 𝜃 un parámetro que normalmente es desco-
nocido.
Uno de los objetivos de la inferencia no paramétrica es desarrollar herramientas
que permitan verificar el grado de cumplimiento de las hipótesis anteriores1 .
Los procedimientos habituales incluyen métodos descriptivos (principalmente
gráficos), contrastes de bondad de ajuste (también de homogeneidad o de datos
atípicos) y contrastes de aleatoriedad.
En este apéndice se describen brevemente algunos de los métodos clásicos, prin-
cipalmente con la idea de que pueden ser de utilidad para evaluar resultados de
simulación y para la construcción de modelos del sistema real (e.g. para mode-
lar variables que se tratarán como entradas del modelo general). Se empleará
principalmente el enfoque de la estadística no paramétrica, aunque también se
mostrarán algunas pequeñas diferencias entre su uso en inferencia y en simula-
ción.
Los métodos genéricos no son muy adecuados para evaluar generadores aleato-
rios (e.g. L’Ecuyer y Simard, 2007). La recomendación sería emplear baterías de
1 El otro objetivo de la inferencia estadística no paramétrica es desarrollar procedimientos
alternativos (métodos de distribución libre) que sean válidos cuando no se verifica alguna de
las hipótesis estructurales.
213
214 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
𝐻0 simple 𝐻0 compuesta
𝐻0 ∶ 𝐹 = 𝒩(0, 1) 𝐻0 ∶ 𝐹 = 𝒩(𝜇, 𝜎2 )
{ {
𝐻1 ∶ 𝐹 ≠ 𝒩(0, 1) 𝐻1 ∶ 𝐹 ≠ 𝒩(𝜇, 𝜎2 )
B.1.1 Histograma
Se agrupan los datos en intervalos 𝐼𝑘 = [𝐿𝑘−1 , 𝐿𝑘 ) con 𝑘 = 1, … , 𝐾 y a cada
intervalo se le asocia un valor (altura de la barra) igual a la frecuencia absoluta
𝑛
de ese intervalo 𝑛𝑘 = ∑𝑖=1 1 (𝑋𝑖 ∈ [𝐿𝑘−1 , 𝐿𝑘 )), si la longitud de los intervalos
es constante, o proporcional a dicha frecuencia (de forma que el área coincida
con la frecuencia relativa y pueda ser comparado con una función de densidad):
𝑛𝑖
𝑓𝑛̂ (𝑥) =
𝑛 (𝐿𝑘 − 𝐿𝑘−1 )
Histogram of datos
0.08
0.06
Density
0.04
0.02
0.00
5 10 15 20 25 30 35
datos
n = 50 n = 500
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
n = 1000 n = 10000
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
⎧ 0 si 𝑥 < 𝑋(1)
{
𝐹𝑛 (𝑥) = ⎨ 𝑛𝑖 si 𝑋(𝑖) ≤ 𝑥 < 𝑋(𝑖+1)
{
⎩ 1 si 𝑋(𝑛) ≤ 𝑥
Ejemplo:
fn <- ecdf(datos)
curve(ecdf(datos)(x), xlim = extendrange(datos), type = 's',
ylab = 'distribution function', lwd = 2)
curve(pnorm(x, mean(datos), sd(datos)), add = TRUE)
B.1. MÉTODOS DE BONDAD DE AJUSTE 217
1.0
0.8
distribution function
0.6
0.4
0.2
0.0
10 15 20 25 30
siendo 𝑥(𝑖) los cuantiles observados y 𝑞𝑖 = 𝐹0−1 (𝑝𝑖 ) los esperados2 bajo 𝐻0 .
Ejemplo:
qqnorm(datos)
qqline(datos)
2 Típicamente (𝑖−0.5)
{𝑝𝑖 = 𝑛 ∶ 𝑖 = 1, ⋯ , 𝑛}.
218 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
30
25
Sample Quantiles
20
15
10
−2 −1 0 1 2
Theoretical Quantiles
require(car)
qqPlot(datos, "norm")
30
25
datos
20
15
10
38
10
−2 −1 0 1 2
norm quantiles
## [1] 10 38
𝐻0 ∶ 𝐹 = 𝐹0
{
𝐻1 ∶ 𝐹 ≠ 𝐹 0
B.1. MÉTODOS DE BONDAD DE AJUSTE 219
𝑛𝑖
(𝑛𝑖 −𝑒𝑖 )2
Clases observadas 𝑝𝑖 bajo 𝐻0 𝑒𝑖 bajo 𝐻0 𝑒𝑖
(𝑛1 −𝑒1 )2
𝐶1 𝑛1 𝑝1 𝑒1 𝑒1
⋮ ⋮ ⋮ ⋮ ⋮
(𝑛𝑘 −𝑒𝑘 )2
𝐶𝑘 𝑛𝑘 𝑝𝑘 𝑒𝑘 𝑒𝑘
𝑘 2
Total ∑𝑖 𝑛𝑖 = 𝑛 ∑ 𝑖 𝑝𝑖 = 1 ∑𝑖 𝑒𝑖 = 𝑛 𝜒2 = ∑𝑖=1 (𝑛𝑖 −𝑒
𝑒𝑖
𝑖)
𝑘
(𝑛𝑖 − 𝑒𝑖 )2
∑ ≥ 𝜒2𝑘−𝑟−1,1−𝛼
𝑖=1
𝑒𝑖
##
## Chi-squared test for given probabilities
##
## data: table(x)
## X-squared = 9.2, df = 4, p-value = 0.05629
La distribución exacta del estadístico del contraste es discreta (se podría aproxi-
mar por simulación, por ejemplo empleando los parámetros simulate.p.value
= TRUE y B = 2000 de la función chisq.test(); ver también el Ejercicio 7.2
de la Sección 7.7.3 para el caso del contraste chi-cuadrado de independencia).
Para que la aproximación continua 𝜒2 sea válida:
• El tamaño muestral debe ser suficientemente grande (p.e. 𝑛 > 30).
• La muestra debe ser una muestra aleatoria simple.
• Los parámetros deben estimarse (si es necesario) por máxima verosimili-
tud.
• Las frecuencias esperadas 𝑒𝑖 = 𝑛 ⋅ 𝑝𝑖 deberían ser todas ≥ 5 (realmente
esta es una restricción conservadora, la aproximación puede ser adecuada
si no hay frecuencias esperadas inferiores a 1 y menos de un 20% inferiores
a 5).
Si la frecuencia esperada de alguna clase es < 5, se suele agrupar con otra clase
(o con varias si no fuese suficiente con una) para obtener una frecuencia esperada
≥ 5:
• Cuando la variable es nominal (no hay una ordenación lógica) se suele
agrupar con la(s) que tiene(n) menor valor de 𝑒𝑖 .
B.1. MÉTODOS DE BONDAD DE AJUSTE 221
Histogram of x
0.10
0.08
0.06
Density
0.04
0.02
0.00
10 15 20 25 30
##
## Pearson's Chi-squared test table
## classes observed expected residuals
## 1 ( 9.06000,14.49908] 6 5.125 0.3865103
B.1. MÉTODOS DE BONDAD DE AJUSTE 223
𝑖
Teniendo en cuenta que 𝐹𝑛 (𝑋(𝑖) ) = 𝑛:
𝑖 𝑖−1
𝐷𝑛 = max { − 𝐹0 (𝑋(𝑖) ), 𝐹0 (𝑋(𝑖) ) − }
1≤𝑖≤𝑛𝑛 𝑛
+ −
= max {𝐷𝑛,𝑖 , 𝐷𝑛,𝑖 }
1≤𝑖≤𝑛
𝑝 = 𝑃 (𝐷𝑛 ≥ 𝑑) ≤ 𝛼.
##
## One-sample Kolmogorov-Smirnov test
##
## data: datos
## D = 0.13239, p-value = 0.4688
## alternative hypothesis: two-sided
Por ejemplo:
ks.test(datos, pnorm, mean(datos), sd(datos)) # One-sample Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: datos
## D = 0.097809, p-value = 0.8277
## alternative hypothesis: two-sided
library(nortest)
lillie.test(datos)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: datos
## D = 0.097809, p-value = 0.4162
B.2. DIAGNOSIS DE LA INDEPENDENCIA 225
Típicamente 𝐶𝑜𝑣(𝑋1 , 𝑋2 ) > 0 por lo que con los métodos “clásicos” (basados
en independencia) se suelen producir subestimaciones de las varianzas (IC más
estrechos y tendencia a rechazar 𝐻0 en contrastes).
Ejemplo: datos simulados
Consideramos un proceso temporal estacionario con dependencia exponencial (la
dependencia entre las observaciones depende del “salto” entre ellas; ver Ejemplo
7.4 en la Sección 7.3).
n <- 100 # Nº de observaciones
t <- seq(0, 1, length = n)
mu <- rep(0, n) # Media
# mu <- 0.25 + 0.5*t
# mu <- sin(2*pi*t)
# Matriz de covarianzas
t.dist <- as.matrix(dist(t))
t.cov <- exp(-t.dist)
226 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
# str(t.cov)
# num [1:100, 1:100] 1 0.99 0.98 0.97 0.96 ...
z <- rnorm(n)
x1 <- mu + z # Datos independientes
x2 <- mvrnorm(1, mu, t.cov) # Datos dependientes
−1
−2
Datos independientes
Datos dependientes
−3
## [1] 0.8067621
var(x2)
## [1] 0.1108155
par(old.par)
Es habitual que este tipo de análisis se realice sobre los residuos de un modelo
de regresión (e.g. datos <- residuals(modelo))
Este gráfico también podría servir para detectar dependencia temporal:
• Valores próximos muy parecidos (valores grandes seguidos de grandes y
viceversa) indicarían una posible dependencia positiva.
• Valores próximos dispares (valores grandes seguidos de pequeños y vice-
versa) indicarían una posible dependencia negativa.
228 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
30
30
25
25
as.ts(datos)
datos
20
20
15
15
10
10
0 10 20 30 40 0 10 20 30 40
Index Time
0.5
1
0.5
0.0
0
−0.5
−1
0.0
−1.0
−2
par(old.par)
1.0
2
1.0
0.5
1
X_t+1
X_t+1
X_t+1
0.0
0.5
−0.5
−1
0.0
−1.0
−2
0.0 0.5 1.0 −2 −1 0 1 2 −1.0 −0.5 0.0 0.5 1.0
par(old.par)
Ejemplo
# Gráfico de dispersion retardado
plot(datos[-length(datos)], datos[-1], xlab = "X_t", ylab = "X_t+1")
30
25
X_t+1
20
15
10
10 15 20 25 30
X_t
## [1] 0.01344127
1.0
1.0
0.8
0.8
0.6
0.6
U_t+1
U_t+1
0.4
0.4
0.2
0.2
0.0
0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
U_t U_t
B.2.4 El correlograma
Para estudiar si el grado de relación (lineal) entre 𝑋𝑖 e 𝑋𝑖+𝑘 podemos utilizar
el coeficiente de correlación:
• Comando R:acf(x)
En caso de independencia es de esperar que las autocorrelaciones muestrales sean
próximas a cero (valores “grandes” indicarían dependencia positiva o negativa
según el signo).
B.2. DIAGNOSIS DE LA INDEPENDENCIA 231
1.0
1.0
0.8
0.8
0.5
0.6
0.6
ACF
ACF
ACF
0.4
0.4
0.0
0.2
0.2
−0.5
−0.2
−0.2
−1.0
0 5 10 15 20 0 5 10 15 20 0 5 10 15 20
par(old.par)
1
𝑟(𝑘) ∼ 𝑁 (𝜌(𝑘), )
𝑎𝑝𝑟𝑜𝑥. 𝑛
2
|𝑟(𝑘)| < √
𝑛
Ejemplo
acf(datos) # correlaciones
232 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
Series datos
1.0
0.8
0.6
0.4
ACF
0.2
0.0
−0.2
0 5 10 15
Lag
Series x2
0.10
0.08
0.06
ACF (cov)
0.04
0.02
0.00
0 5 10 15 20
Lag
𝐻0 ∶ La muestra es aleatoria
{
𝐻1 ∶ La muestra no es aleatoria
𝑅 − 𝐸(𝑅)
𝑝 ≃ 2𝑃 (𝑍 ≥ ∣ ∣)
√𝑉 𝑎𝑟(𝑅)
library(tseries)
runs.test(as.factor(datos > median(datos)))
##
## Runs Test
##
## data: as.factor(datos > median(datos))
## Standard Normal = -0.4422, p-value = 0.6583
## alternative hypothesis: two.sided
Alternativamente, para evitar el cálculo del punto de corte (la mediana), reque-
rido para dicotomizar la variable continua, se podría emplear una modificación
de este contraste, el denominado test de rachas ascendentes y descendentes, en
el que se generan los valores + y − dependiendo de si el valor de la secuencia
es mayor o menor que el anterior (ver e.g. Downham, 1970). Este contraste es
más adecuado para generadores aleatorios.
𝑝 = 𝑃 (𝜒2𝑚 ≥ 𝑄)
• Comandos R:
Box.test(x, type=Ljung)
Box.test(x, lag, type=Ljung)
Ejemplo
Box.test(datos, type="Ljung") # Contrasta si la primera autocorrelación es nula
##
## Box-Ljung test
B.3. CONTRASTES ESPECÍFICOS PARA GENERADORES ALEATORIOS235
##
## data: datos
## X-squared = 0.0078317, df = 1, p-value = 0.9295
Box.test(datos, lag=5, type="Ljung") # Contrasta si las 5 primeras autocorrelaciones son nulas
##
## Box-Ljung test
##
## data: datos
## X-squared = 1.2556, df = 5, p-value = 0.9394
NOTA: Cuando se trabaja con residuos de un modelo lineal, para contrastar
que la primera autocorrelación es cero, es preferible emplear el test de Durbin-
Watson implementado en la función dwtest() del paquete lmtest.
set.seed(1)
u <- runif(1000)
k <- 10
x <- floor(k*u) + 1
# Test chi-cuadrado
f <- table(factor(x, levels = seq_len(k)))
chisq.test(f)
##
## Chi-squared test for given probabilities
##
## data: f
## X-squared = 10.26, df = 9, p-value = 0.3298
# Equivalente a
# source("Test Chi-cuadrado continua.R")
# chisq.test.cont(u, distribution = "unif", nclasses = k, output = FALSE, min = 0, max
𝐾!
𝑃 (𝐶 = 𝑐) = 𝑆(𝑑, 𝑐),
(𝐾 − 𝑐)!𝐾 𝑑
para obtener las frecuencias esperadas de cada clase y confrontarlas con las ob-
servadas vía el estadístico chi-cuadrado (e.g. Knuth, 2002, Sección 3.3.2, p. 65).
Existen varias elecciones comunes de 𝐾 y 𝑀 . Tomando 𝐾 = 5 con clases 𝑄 = 5,
𝑄 = 6, …, 𝑄 = 19, 𝑄 ≥ 20, las probabilidades vendrían dadas por:
𝑃 (𝑄 = 5) = 0.03840000, 𝑃 (𝑄 = 6) = 0.07680000,
𝑃 (𝑄 = 7) = 0.09984000, 𝑃 (𝑄 = 8) = 0.10752000,
𝑃 (𝑄 = 9) = 0.10450944, 𝑃 (𝑄 = 10) = 0.09547776,
𝑃 (𝑄 = 11) = 0.08381645, 𝑃 (𝑄 = 12) = 0.07163904,
𝑃 (𝑄 = 13) = 0.06011299, 𝑃 (𝑄 = 14) = 0.04979157,
𝑃 (𝑄 = 15) = 0.04086200, 𝑃 (𝑄 = 16) = 0.03331007,
𝑃 (𝑄 = 17) = 0.02702163, 𝑃 (𝑄 = 18) = 0.02184196,
𝑃 (𝑄 = 19) = 0.01760857, 𝑃 (𝑄 ≥ 20) = 0.07144851.
Integración numérica
239
240 APÉNDICE C. INTEGRACIÓN NUMÉRICA
4
3
fun(x)
2
1
0
trapezoid(fun, 0, 1, 20)
## [1] 1.0025
El error en esta aproximación se corresponde con:
(𝑏 − 𝑎)3 ″
𝑓 (𝜉),
12𝑛2
para algún 𝑎 ≤ 𝜉 ≤ 𝑏 (dependiendo del signo de la segunda derivada, i.e. de
si la función es cóncava o convexa, el error será negativo ó positivo). El error
3
máximo absoluto es (𝑏−𝑎) ″
12𝑛2 max𝑎≤𝜉≤𝑏 |𝑓 (𝜉)|. En el caso general multidimensional
2
sería 𝑂(𝑛− 𝑑 ).
𝑏 (𝑛/2)−1 𝑛/2
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 ≈ [𝑓(𝑥0 ) + 2 ∑ 𝑓(𝑥2𝑗 ) + 4 ∑ 𝑓(𝑥2𝑗−1 ) + 𝑓(𝑥𝑛 )],
𝑎 3 𝑗=1 𝑗=1
simpson(fun, 0, 1, 20)
## [1] 1
242 APÉNDICE C. INTEGRACIÓN NUMÉRICA
(𝑏 − 𝑎)5
max ∣𝑓 (4) (𝜉)∣ .
180𝑛4 𝑎≤𝜉≤𝑏
4
En el caso general multidimensional sería 𝑂(𝑛− 𝑑 ).
level = 1
S.old <- (b-a) * (fun(a) + fun(b))/2
S.new <- quadrature_internal(S.old, fun,
a, (a+b)/2, b, tol, level+1)
return(S.new)
}
C.2. INTEGRACIÓN NUMÉRICA BIDIMENSIONAL 243
quadrature(fun, 0, 1)
## [1] 1
Fuente r-blogger Guangchuang Yu?
C.1.4 Comandos de R
integrate(fun, 0, 1) # Permite límites infinitos
require(MASS)
area(fun, 0, 1)
## [1] 1
.
Consideraremos como ejemplo:
1 1
∫ ∫ (𝑥2 − 𝑦2 ) 𝑑𝑥𝑑𝑦 = 0
−1 −1
.
f2d <- function(x,y) x^2 - y^2
hx <- x[2]-x[1]
hy <- y[2]-y[1]
# persp3D(z = z, x = x, y = y)
persp3D.f2d <- function(f2d, ax=-1, bx=1, ay=-1, by=1, nx=21, ny=21, ...) {
x <- seq(ax, bx, length = nx)
y <- seq(ay, by, length = ny)
hx <- x[2]-x[1]
hy <- y[2]-y[1]
z <- outer(x, y, f2d)
persp3D(x, y, z, ...)
}
0.5
0.5
0.0
z
0.0
−0.5
1.0
−1.0
0.5
−0.5 −0.5
0.0
x0.0
y
−0.5
0.5
1.0−1.0
## [1] -8.881784e-18
C.2.3 Comandos de R
Suponiendo que la función es vectorial, podemos emplear:
integrate( function(y) {
sapply(y, function(y) {
integrate(function(x) f2d(x,y), ax, bx)$value }) },
ay, by)
Fuente tolstoy.newcastle.edu.au.
Alternativamente se podría emplear la función adaptIntegrate del paquete
cubature.