Simulacion
Simulacion
Simulacion
2020-04-21
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 . . . . . . . . . . . . . . . . . . . . . . 10
1.4 Números pseudoaleatorios . . . . . . . . . . . . . . . . . . . . . . 12
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 . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3
4 ÍNDICE GENERAL
Referencias 207
Bibliografía básica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
Bibliografía complementaria . . . . . . . . . . . . . . . . . . . . . . . . 207
ÍNDICE GENERAL 5
A Enlaces 211
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
• Cuando es imposible experimentar sobre el sistema real por ser dicha ex-
perimentación destructiva.
Aplicaciones:
• Estadística:
– Muestreo, remuestreo, …
– Aproximación de distribuciones (de estadísticos, estimadores, …)
– Realización de contrastes, intervalos de confianza, …
– Comparación de estimadores, contrastes, …
– Validación teoría (distribución asintótica,…)
– Inferencia Bayesiana
• Optimización: Algoritmos genéticos, …
• Análisis numérico: Aproximación de integrales, resolución de ecuaciones,
…
• Computación: Diseño, verificación y validación de algoritmos,…
• Criptografía: Protocolos de comunicación segura, …
• Física: Simulación de fenómenos naturales, …
• …
En los capítulos 8 y 9 nos centraremos en algunas de las aplicaciones de utilidad
en Estadística.
12 CAPÍTULO 1. INTRODUCCIÓN A LA SIMULACIÓN
Una sucesión de números aleatorios puros (true random), se caracteriza por que
no existe ninguna regla o plan que nos permita conocer sus valores.
Figura 1.1: Líneas 10580-10594, columnas 21-40, del libro *A Million Random
Digits with 100,000 Normal Deviates*.
• http://www.fourmilab.ch/hotbits
• http://software.intel.com.
• http://spectrum.ieee.org
1.3.1 Inconvenientes:
1.3.2 Alternativas:
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
par(par.old)
donde:
• k es el orden del generador.
• (x0 , x1 , · · · , xk−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 p.
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).
1.4. NÚMEROS PSEUDOALEATORIOS 15
Números aleatorios en R
17
18 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
2.1 Opciones
– “Knuth-TAOCP”
– “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).
( )
b) Aproximar el valor de π mediante simulación a partir de P X 2 + Y 2 ≤ 1 .
set.seed(1)
n <- 10000
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
indice <- (x^2+y^2 < 1)
mean(indice)
## [1] 0.7806
pi/4
## [1] 0.7853982
pi_aprox <- 4*mean(indice)
pi_aprox
## [1] 3.1224
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)
2.3. EJERCICIOS 21
1.0
0.5
0.0
Y
−0.5
−1.0
Ejercicio 2.2.
## [1] 0.4953
barplot(100*table(x)/nsim) # Representar porcentajes
22 CAPÍTULO 2. NÚMEROS ALEATORIOS EN R
50
40
30
20
10
0
0 1
## [1] 0.394
plot(n, cumsum(x)/n, type="l", ylab="Proporción de caras",
xlab="Número de lanzamientos", ylim=c(0,1))
abline(h=p, lty=2, col="red")
2.3. EJERCICIOS 23
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:
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 nsim = 10000 jugadas de este
juego y calcular la proporción de veces que se gana la apuesta (obtener al
menos un 6 en n lanzamientos), usando n = 4. Comparar el resultado con
la probabilidad teórica 1 − (5/6)n .
set.seed(1)
n <- 4
nsim <- 10000
mean(replicate(nsim, deMere(n)))
## [1] 0.5148
1-(5/6)^n
2.4. TIEMPO DE CPU 25
## [1] 0.5177469
La velocidad del generador suele ser una característica importante. Para evaluar
el rendimiento están disponibles en R distintas herramientas:
Por ejemplo, podríamos emplear las siguientes funciones para ir midiendo los
tiempos de CPU durante una simulación:
CPUtimeini <- function() {
.tiempo.ini <<- proc.time()
.tiempo.last <<- .tiempo.ini
}
## [1] 0.5003313
CPUtimeprint()
##
## Tiempo última operación:
## user system elapsed
## 0.06 0.00 0.06
## Tiempo total operación:
## user system elapsed
## 0.06 0.00 0.06
funtest(1000)
## [1] 0.5050682
CPUtimeprint()
##
## Tiempo última operación:
## user system elapsed
## 0.02 0.00 0.02
## Tiempo total operación:
## user system elapsed
## 0.08 0.00 0.08
2.4.2 Paquetes de R
#' @export
cpu.time <- .cpu.time.ini()
Ejemplo:
cpu.time(reset = TRUE)
Generación de números
pseudoaleatorios con
distribución uniforme
xi = (axi−1 + c) mod m
xi
ui =
m
i = 1, 2, . . .
29
30CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNI
# --------------------------------------------------
# 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
}
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 k-uplas de un generadores multipli-
3.1. GENERADORES CONGRUENCIALES (LINEALES) 33
1/k
cativo están contenidas en a lo sumo (k!m) hiperplanos paralelos (para más
detalles sobre la estructura reticular, ver por ejemplo Ripley, 1987, sección 2.7).
Por tanto habría que seleccionar adecuadamente m y c (a solo influiría en la
pendiente) de forma que la estructura reticular sea impreceptible teniendo en
cuenta el número de datos que se pretende generar (por ejemplo de forma que
la distancia mínima entre los puntos sea próxima a la esperada en teoría).
Se han propuesto diversas pruebas (ver Sección 3.2) para determinar si un 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 m = 231 − 1).
• En cualquier caso, se recomienda considerar un “periodo de seguridad”
√
≈ p para evitar este tipo de problemas.
• Aunque estos generadores tiene limitaciones en su capacidad para producir
secuencias muy largas de números i.i.d. U(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 m = 2 y ai ∈ {0, 1} y se van com-
binando l bits para obtener valores en el intervalo (0, 1), por ejemplo ui =
34CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNI
Ejercicio 3.1.
Histogram of u
1.0
0.8
0.6
Density
0.4
0.2
0.0
## [1] 0.4999609
## [1] 0.402
mean((0.4 < u) & (u < 0.8)) # Alternativa
## [1] 0.402
36CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNI
#-------------------------------------------------------------------------------
chisq.test.cont <- function(x, distribution = "norm", nclasses = floor(length(x)/5),
output = TRUE, nestpar = 0, ...) {
# Funciones distribución
q.distrib <- eval(parse(text = paste("q", distribution, sep = "")))
d.distrib <- eval(parse(text = paste("d", distribution, sep = "")))
# Puntos de corte
q <- q.distrib((1:(nclasses - 1))/nclasses, ...)
tol <- sqrt(.Machine$double.eps)
xbreaks <- c(min(x) - tol, q, max(x) + tol)
# Gráficos y frecuencias
if (output) {
xhist <- hist(x, breaks = xbreaks, freq = FALSE, lty = 2, border = "grey50")
curve(d.distrib(x, ...), add = TRUE)
} else {
xhist <- hist(x, breaks = xbreaks, plot = FALSE)
}
# Cálculo estadístico y p-valor
O <- xhist$counts # Equivalente a table(cut(x, xbreaks)) pero más eficiente
E <- length(x)/nclasses
DNAME <- deparse(substitute(x))
METHOD <- "Pearson's Chi-squared test"
STATISTIC <- sum((O - E)^2/E)
names(STATISTIC) <- "X-squared"
PARAMETER <- nclasses - nestpar - 1
names(PARAMETER) <- "df"
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
# Preparar resultados
classes <- format(xbreaks)
classes <- paste("(", classes[-(nclasses + 1)], ",", classes[-1], "]",
sep = "")
RESULTS <- list(classes = classes, observed = O, expected = E, residuals = (O -
E)/sqrt(E))
if (output) {
cat("\nPearson's Chi-squared test table\n")
print(as.data.frame(RESULTS))
}
if (any(E < 5))
warning("Chi-squared approximation may be incorrect")
structure(c(list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL,
method = METHOD, data.name = DNAME), RESULTS), class = "htest")
}
Histogram of x
1.0
0.8
0.6
Density
0.4
0.2
0.0
##
## Pearson's Chi-squared test table
## classes observed expected residuals
## 1 (-1.490116e-08, 1.000000e-01] 51 50 0.1414214
## 2 ( 1.000000e-01, 2.000000e-01] 49 50 -0.1414214
## 3 ( 2.000000e-01, 3.000000e-01] 49 50 -0.1414214
## 4 ( 3.000000e-01, 4.000000e-01] 50 50 0.0000000
## 5 ( 4.000000e-01, 5.000000e-01] 51 50 0.1414214
## 6 ( 5.000000e-01, 6.000000e-01] 51 50 0.1414214
## 7 ( 6.000000e-01, 7.000000e-01] 49 50 -0.1414214
## 8 ( 7.000000e-01, 8.000000e-01] 50 50 0.0000000
## 9 ( 8.000000e-01, 9.000000e-01] 50 50 0.0000000
## 10 ( 9.000000e-01, 9.980469e-01] 50 50 0.0000000
##
## Pearson's Chi-squared test
##
## data: u
## X-squared = 0.12, df = 9, p-value = 1
Como se muestra en la Figura 3.3 el histograma de la secuencia generada es muy
plano (comparado con lo que cabría esperar de una muestra de tamaño 500 de
3.2. ANÁLISIS DE LA CALIDAD DE UN GENERADOR 39
Ejercicio 3.2.
0.4
0.2
0.0
# Test de Kolmogorov-Smirnov
ks.test(u, "punif", 0, 1)
##
## One-sample Kolmogorov-Smirnov test
##
## data: u
## D = 0.0033281, p-value = 1
## alternative hypothesis: two-sided
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 (ui , ui+k ), con k = 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 UNI
## 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:
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
##
44CAPÍTULO 3. GENERACIÓN DE NÚMEROS PSEUDOALEATORIOS CON DISTRIBUCIÓN UNI
## data: estadistico
## D = 0.023499, p-value = 0.6388
## alternative hypothesis: two-sided
Análisis de los p-valores:
# Histograma
hist(pvalor, freq=FALSE)
abline(h=1) # curve(dunif(x,0,1), add=TRUE)
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
Adicionalmente, si queremos estudiar la proporción de rechazos (el tamaño del
contraste) para los posibles valores de α, podemos emplear la distribución em-
pírica del p-valor (proporción de veces que resultó menor que un determinado
valor):
# Distribución empírica
curve(ecdf(pvalor)(x), type = "s", lwd = 2,
xlab = 'Nivel de significación', ylab = 'Proporción de rechazos')
3.2. ANÁLISIS DE LA CALIDAD DE UN GENERADOR 45
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.
10(2k− 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)))
3.3. EJERCICIOS 47
# -------------------------------------------------
# RANDVN()
# Genera un valor pseudoaleatorio con el generador de Von Neumann
# Actualiza la semilla (si no existe llama a initRANDVN)
RANDVN <- function() {
if (!exists(".semilla", envir=globalenv())) initRANDVN()
z <- .semilla^2
.semilla <<- trunc((z-trunc(z/.aux)*.aux)/.aux2)
return(.semilla/10^.n)
}
# -------------------------------------------------
# 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
49
50 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
## [1] 0.5047
Podemos generar un gráfico con la evolución de la aproximación con el siguiente
código:
plot(cumsum(rx)/1:nsim, type="l", lwd=2, xlab="Número de generaciones",
ylab="Proporción muestral", ylim=c(0,1))
abline(h = mean(rx), lty = 2)
# valor teórico
abline(h = p)
1.0
0.8
Proporción muestral
0.6
0.4
0.2
0.0
Número de generaciones
Una suposición crucial es que las variables Xi deben tener varianza finita (real-
mente esta suposición puede relajarse: E (|Xi |) < ∞). 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
( ) Sb2
Vd
ar X n =
n
con:
1 ∑( )2
n
Sbn2 = Xi − X .
n − 1 i=1
p̂n (1 − p̂n )
Vd
ar (p̂n ) = ,
n−1
Xn − µ d
Zn = → N (0, 1)
√σ
n
i.e. lim FZn (z) = Φ(z). Por tanto, un intervalo de confianza asintótico para µ
n→∞
es: ( )
Sbn Sbn
IC1−α (µ) = X n − z1−α/2 √ , X n + z1−α/2 √ .
n n
Sbn
Podemos considerar que z1−α/2 √ es la precisión obtenida (con nivel de con-
n
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
forma que:
Sbn
z1−α/2 √ < ε
n
3.1. j = j + 1.
⌈( / )2 ⌉
3.2. nj = z1−α/2 Sbnj−1 ε .
y calcular X nj y Sbnj .
nj
3.3. Generar {Xi }i=nj−1 +1
/
√
4. Devolver X nj y z1−α/2 Sbnj nj .
Sbn
z1−α/2 √ < ε X n .
n
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"
56 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
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:
4.5. EL PROBLEMA DE LA DEPENDENCIA 57
esterr[nsim]
## [1] 0.004599203
En este caso al haber dependencia positiva se produce una subestimación del
verdadero error estandar.
acf(as.numeric(rx))
Series as.numeric(rx)
1.0
0.8
0.6
ACF
0.4
0.2
0.0
0 10 20 30 40
Lag
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
Esta forma de proceder podría ser adecuada para tratar de aproximar la preci-
4.5. EL PROBLEMA DE LA DEPENDENCIA 59
sión:
esterr[nrxi]
## [1] 0.02277402
pero no sería eficiente para aproximar la media. Siempre será preferible emplear
todas las observaciones.
Por ejemplo, se podría pensar en considerar las medias de grupos de 24 valores
consecutivos y suponer que hay independencia entre ellas:
rxm <- rowMeans(matrix(rx, ncol = lag, byrow = TRUE))
nrxm <- length(rxm)
n <- 1:nrxm
est <- cumsum(rxm)/n
esterr <- sqrt(cumsum((rxm-est)^2))/n # Error estándar
plot(est, type="l", lwd=2, ylab="Probabilidad",
xlab=paste("Número de simulaciones /", lag + 1), ylim=c(0,0.6))
abline(h = est[length(rxm)], lty=2)
lines(est + 2*esterr, lty=2) # OJO! Supone independencia
lines(est - 2*esterr, lty=2)
abline(h = 1/3, col="darkgray") # Prob. teor. cadenas Markov
0.6
0.5
0.4
Probabilidad
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:
60 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓ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 Xi y Xi+k es aproximadamente nula, y consideramos las
subsecuencias (lotes) (Xt+1 , Xt+2 , . . . , Xt+k ) con t = (j − 1)k, j = 1, . . . , m y
n = mk. Entonces:
( )
( ) 1∑ ∑ ∑
n m ik
1 1
V ar X̄ = V ar Xi = V ar Xt
n i=1 m j=1 k
t=(i−1)k+1
1 ∑ ∑ ( )
m ik
1 1
≈ 2 V ar Xt ≈ V ar X̄k
m j=1 k m
t=(i−1)k+1
## [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
}
4.5. EL PROBLEMA DE LA DEPENDENCIA 61
La idea sería considerar las medias de las series como una muestra independiente
de una nueva variable y estimar su varianza de la forma habitual:
# Media de cada secuencia
n <- 1:nsim
est <- apply(rxm, 1, function(x) cumsum(x)/n)
matplot(n, est, type = 'l', lty = 3, col = "lightgray",
ylab="Probabilidad", xlab="Número de simulaciones")
# Aproximación
mest <- apply(est, 1, mean)
lines(mest, lwd = 2)
abline(h = mest[nsim], lty = 2)
# Precisión
mesterr <- apply(est, 1, sd)/sqrt(nsec)
lines(mest + 2*mesterr, lty = 2)
lines(mest - 2*mesterr, lty = 2)
# Prob. teor. cadenas Markov
abline(h = 1/3, col="darkgray")
1.0
0.8
0.6
Probabilidad
0.4
0.2
0.0
Número de simulaciones
# Aproximación final
mest[nsim] # mean(rxm)
## [1] 0.3089
62 CAPÍTULO 4. ANÁLISIS DE RESULTADOS DE SIMULACIÓN
# 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.
Xt = µ + ρ ∗ (Xt−1 − µ) + εt
Podemos tener en cuenta que en este caso la varianza es:
σε2
var(Xt ) = E(Xt2 ) − µ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 X 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
U = F (X) ∼ U (0, 1) ,
ya que:
( ) ( )
G (u) = P (Y ≤ u) = P (F (X) ≤ u) = P X ≤ F −1 (u) = F F −1 (u) = u.
F −1 (U ) ∼ X
67
68 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
2. Devolver X = F −1 (U ).
ln (1 − u)
1 − e−λx = u ⇔ x = − ,
λ
el algoritmo para simular esta variable mediante el método de inversión es:
1. Generar U ∼ U (0, 1).
ln (1 − U )
2. Devolver X = − .
λ
En el último paso podemos emplear directamente U en lugar de 1 − U , ya que
1 − U ∼ U (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 F (x) F −1 (U ) simplificada
ln (1 − U ) ln U
exp (λ) λe−λx , si 1 − e−λx − −
λ λ
(λ > 0) x≥0 ( ( ))
1 1 arctan x 1
Cauchy + tan π U − tan πU
π (1 + x2 ) 2( π ) 2
2( x) 2 x2 ( √ ) ( √ )
Triangular 1− , x− a 1− 1−U a 1− U
a a a 2a
en (0, a) si 0 ≤ x ≤ a ( )a
aba b b b
Pareto , si 1−
xa+1 x (1 − U )
1/a U 1/a
(a, b > 0) x≥b
1/α 1/α
(− ln (1 − U )) (− ln U )
αλα xα−1 e−(λx) 1, − e−(λx)
α α
Weibull
λ λ
(λ, α > 0) si x ≥ 0 strut
Ejercicio 5.1.
Ventajas:
• Aplicable, en principio, a cualquier distribución continua.
Inconvenientes:
• Puede no ser posible encontrar una expresión explícita para F −1 (u) .
• Aún disponiendo de una expresión explícita para F −1 (u), su evaluación
directa puede requerir mucho tiempo de computación.
Alternativas:
5.1. MÉTODO DE INVERSIÓN 71
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.
∑4 ∑4
siendo A (x) = i=0 ai xi y B (x) = i=0 bi x
i
con:
72 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
a0 = −0.322232431088 b0 = 0.0993484626060
a1 = −1 b1 = 0.588581570495
a2 = −0.342242088547 b2 = 0.531103462366
a3 = −0.0204231210245 b3 = 0.103537752850
a4 = −0.0000453642210148 b4 = 0.0038560700634
4. Devolver X.
{ } ∫ b
Area de (x, y) ∈ R2 : a < x < b; 0 ≤ y ≤ f (x)
P (a < X < b) = = f (x) dx
Area de Af a
Por tanto, si (T, Y ) sigue una distribución uniforme en Acg , aceptando los va-
lores de (T, Y ) que pertenezcan a Af (o a Af ∗ ) se obtendrán generaciones con
distribución uniforme sobre Af (o Af ∗ ) y la densidad de la primera componente
T será f .
74 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
5.2.1 Algoritmo
f (x) ≤ c · g (x) , ∀x ∈ R.
Ejercicio 5.2.
1 x2
f (x) = √ e− 2 , x ∈ R,
2π
5.2. MÉTODO DE ACEPTACIÓN RECHAZO 75
√
2e
copt = ≃ 1. 3155.
π
# 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 (Af ) 1
p= = .
area (Acg ) c
Por tanto:
1
E (N ) = =c
p
es el número medio de iteraciones del algoritmo (el número medio de pares de
variables (T, U ) 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 c más eficiente
será el algoritmo (el caso de c = 1 se correspondería con g = f y no tendría senti-
do emplear este algoritmo). El principal problema con este método es encontrar
una densidad auxiliar g de forma que:
f (x)
copt = max .
{x:g(x)>0} g (x)
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
}
# 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
cat("Número de generaciones = ", ngen)
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
s(x) ≤ f (x).
Algoritmo:
1. Generar U ∼ U (0, 1) y T ∼ g.
2. Si c · U · g (T ) ≤ s (T ) devolver X = T ,
en caso contrario
2.a. si c · U · g (T ) ≤ f (T ) devolver X = T ,
2.b. en caso contrario volver al paso 1.
Cuanto mayor sea el área bajo s (x) (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 g y
s 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.
Se puede ver como una modificación del método de aceptación rechazo, de espe-
cial interés cuando el soporte no es acotado.
84 CAPÍTULO 5. SIMULACIÓN DE VARIABLES CONTINUAS
λ −λ|x|
f (x) = e , ∀x ∈ R,
2
se deduce que:
1 1
f (x) = f1 (x) + f2 (x)
2 2
siendo:
{ {
λe−λx si x ≥ 0 λeλx si x < 0
f1 (x) = f2 (x) =
0 six < 0 0 six ≥ 0
Algoritmo:
1. Generar U, V ∼ U (0, 1).
5.5. MÉTODOS ESPECÍFICOS PARA LA GENERACIÓN DE ALGUNAS DISTRIBUCIONES NOTABLES85
Observaciones:
Algoritmo:
1. Generar Y ∼ h.
2. Generar X ∼ g(·|Y ).
Simulación de variables
discretas
xi x1 x2 ··· xn ···
P (X = xi ) p1 p2 ··· pn ···
Considerando como partida una U (0, 1), la idea general consiste en asociar a
cada posible valor xi de X 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
Q (u) ≤ x ⇐⇒ u ≤ F (x).
i <- 1
while (Fx[i] < U[j]) i <- i + 1
X[j] <- x[i]
ncomp <<- ncomp + i
}
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.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))
92 CAPÍTULO 6. SIMULACIÓN DE VARIABLES DISCRETAS
## [1] 0.00234375
max(abs(res$psim - res$pteor) / res$pteor)
## [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
i 1 2 ··· n ···
P (I = i) p1 p2 ··· pn ···
∑
n−1
ipi + (n − 1) pn .
i=1
Ejercicio 6.2.
ncomp <- 0
ind <- order(fmp, decreasing=TRUE)
rx <- rfmp(x[ind], fmp[ind], nsim)
## [1] 3.08969
sum((1:length(x))*fmp[ind]) # Valor teórico
## [1] 3.083984
Como ya se comentó, en R se recomienda emplear la función sample (implementa
eficientemente el método de Alias):
system.time( rx <- sample(x, nsim, replace = TRUE, prob = fmp) )
j0 = ⌊mU ⌋ + 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:
n
E (N ) ≤ 1 +
m
Algoritmo 6.3 (de simulación mediante una tabla guía)
Inicialización:
1. Hacer F1 = p1 .
2. Desde i = 2 hasta n hacer Fi = Fi−1 + pi .
Cálculo de la tabla guía:
1. Hacer g1 = 1 e i = 1.
2. Desde j = 2 hasta m hacer
2.a Mientras (j − 1)/m > Fi hacer i = i + 1.
2.b gj = i
Simulación mediante tabla guía:
1. Generar U ∼ U (0, 1).
2. Hacer j = ⌊mU ⌋ + 1.
3. Hacer i = gj .
4. Mientras U > Fi hacer i = i + 1.
5. Devolver X = xi .
Ejercicio 6.3.
Diseñar una rutina que permita generar nsim valores de una distribución dis-
creta usando una tabla guía. Repetir el ejercicio anterior empleando esta rutina
con m = n − 1.
rfmp.tabla <- function(x, prob = 1/length(x), m, nsim = 1000) {
# Simulación v.a. discreta a partir de función de masa de probabilidad
# por tabla guia de tamaño m
6.3. MÉTODO DE ALIAS 95
# Inicializar tabla y FD
Fx <- cumsum(prob)
g <- rep(1,m)
i <- 1
for(j in 2:m) {
while (Fx[i] < (j-1)/m) i <- i+1
g[j] <- i
}
# Generar valores
X <- numeric(nsim)
U <- runif(nsim)
for(j in 1:nsim) {
i <- g[floor(U[j]*m)+1]
while (Fx[i] < U[j]) i <- i + 1
X[j] <- x[i]
}
return(X)
}
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
Diseñar una rutina que permita generar nsim 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
Los métodos anteriores están pensados para variables que toman un número
finito de valores. Si la variable discreta tiene dominio infinito no se podrían
almacenar las probabilidades acumuladas, aunque en algunos casos podrían cal-
cularse de forma recursiva.
e−λ λj−1
pj = P (X = xj ) = P (X = j − 1) = , j = 1, 2, . . .
(j − 1)!
6.5. CÁLCULO DIRECTO DE LA FUNCIÓN CUANTIL 99
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) ) )
}
6.6. OTROS MÉTODOS 101
# Empleando ifelse se complica un poco más pero el resultado es una función vectorial.
curve(fdistr(x), from = -0.1, to = 1.1, type = 's',
main = 'Función de distribución')
# Discontinuidades en 0 y 1/5
abline(h = c(1/10, 2/10, 3/10), lty = 2)
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 U ∼ U (0, 1)
2. Devolver X = Q (U )
En este caso concreto:
1. Generar U ∼ U (0, 1)
1
2. Si U < 10 devolver X = 0
3. Si U < 2
10 devolver X = 2(U − 1
10 )
3 2
4. Si U < 10 devolver X = 10
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
t
Si µ = (µ1 , µ2 , . . . , µd ) es un vector (de medias) y Σ es una matriz d×d definida
positiva (de varianzas-covarianzas), el vector aleatorio X sigue una distribución
normal multivariante con esos parámetros, X ∼ Nd (µ, Σ), si su función de
densidad es de la forma:
( )
1 1 t −1
f (x) = exp − (x − µ) Σ (x − µ) ,
(2π)n/2 |Σ|1/2 2
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
2d 1 1
U 1 d (T) ≤ 1C (T) ,
Vd (1) 2d [−1,1] Vd (1) d
7.3. FACTORIZACIÓN DE LA MATRIZ DE COVARIANZAS 109
o, lo que es lo mismo, U 1[−1,1]d (T) ≤ 1Cd (T). Dado que el número aleatorio
U está en el intervalo (0, 1) y que las funciones indicadoras valen 0 ó 1, esta
condición equivale a que 1[−1,1]d (T) = 1Cd (T), es decir, a que T ∈ Cd , es decir,
que se verifique:
T12 + T22 + · · · + Td2 ≤ 1.
( )
d
Por otra parte, la simulación de T ∼ U [−1, 1] puede hacerse trivialmente
mediante Ti ∼ U (−1, 1) para cada i = 1, 2, . . . , d, ya que las componentes son
independientes. Como el valor de U es superfluo en este caso, el algoritmo queda:
1. Simular V1 , V2 , . . . , Vd ∼ U (0, 1) independientes.
2. Para i = 1, 2, . . . , d hacer Ti = 2Vi − 1.
3. Si T12 + T22 + · · · + Td2 > 1 entonces volver al paso 1.
t
4. Devolver X = (T1 , T2 , . . . , Td ) .
Ver el Ejercicio 2.1 para el caso de d = 2.
Usando las fórmulas del “volumen” de una “esfera” d-dimensional:
d/2 d
π r
si d es par
Vd (r) = (d/2)!
⌊ 2 ⌋+1 π ⌊ 2 ⌋ rd
d d
2 si d es impar
1 · 3 · 5···d
puede verse que el número medio de iteraciones del algoritmo, dado por la
d
constante c = Vd2(1) , puede llegar a ser enórmemente grande. Así, si d = 2
se tiene c = 1.27, si d = 3 se tiene c = 1.91, si d = 4 entonces c = 3.24 y
para d = 10 resulta c = 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.
Cov(AX) = AAt .
Proposición 7.1
Si X ∼ Nd (µ, Σ) y A es una matriz p × d, de rango máximo, con p ≤ d,
entonces:
( )
AX ∼ Np Aµ, AΣAt .
Y = µ + HΛ1/2 Z ∼ Nd (µ, Σ) .
Y = µ + LZ ∼ Nd (µ, Σ) .
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: 0x000000002e2a61e8>
## <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.
es decir: { √
2
π 1 − x21 si x1 ∈ [−1, 1]
f1 (x1 ) =
0 si x1 ∈
/ [−1, 1]
Además:
1 [ √ √ ]
f (x1 , x2 ) 1
f2 (x2 |x1 ) = = √ 2 = √
π
, si x2 ∈ − 1 − x1 , 1 − x1
2 2
f1 (x1 ) 2 1−x1 2 1 − x21
π
tendríamos que:
1
f (x1 , x2 ) = √
2πσ1 σ2 1 − ρ2
( ( ))
1 (x1 − µ1 )2 (x2 − µ2 )2 2ρ(x1 − µ1 )(x2 − µ2 )
exp − + −
2(1 − ρ2 ) σ12 σ22 σ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 (x1 − µ1 )2
f1 (x1 ) = √ exp −
σ1 2π 2σ12
f (x1 , x2 )
f2 (x2 |x1 ) =
f1 (x1 )
( )2
x2 − µ 2 − − µ1 ) σ2
1 σ1 ρ(x1
= √ exp −
σ2 2π(1 − ρ2 ) 2σ2 (1 − ρ )
2 2
( ) ( )
Por tanto X1 ∼ N µ1 , σ12 y X2 |X1 ∼ N µ2 + σσ21 ρ(X1 − µ1 ), σ22 (1 − ρ2 ) , y
el algoritmo sería el siguiente:
Algoritmo 7.4 (de simulación de una normal bidimensional) 1. Simular
Z1 , Z2 ∼ N (0, 1) independientes.
2. Hacer X1 = µ1 + σ1 Z1 .
√
3. Hacer X2 = µ2 + σ2 ρZ1 + σ2 1 − ρ2 Z 2 .
Este algoritmo es el mismo que obtendríamos con la factorización de la matrix
de covarianzas ya que Σ = LLt con:
( 2 )
σ1 √ 0
L=
σ2 ρ σ2 1 − ρ2
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:
3
2
1
0
y
−1
−2
−3
Figura 7.4: Realizaciones condicionales del proceso funcional del Ejemplo 7.7.
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
0.4
0.2
0.0
# Bucle simulación
nsim <- 1 # 1000
for (i in 1:nsim) {
y <- L %*% rnorm(n)
# calcular estadísticos, errores,...
}
120CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
## [,1]
## [1,] -0.62645381
## [2,] -0.05969442
## [3,] -0.98014198
## [4,] 1.09215113
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
## krige.conv: model with constant mean
## krige.conv: sampling from the predictive distribution (conditional simulations)
7.5. SIMULACIÓN CONDICIONAL E INCONDICIONAL 121
simul. cond. 1
1.0
0.8
0.6
Y Coord
0.4 0.2
0.0
simul. cond. 2
1.0
0.8
0.6
Y Coord
0.4 0.2
0.0
simul. cond. 3
1.0
0.8
0.6
Y Coord
0.4 0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
X Coord
simul. cond. 3
1.0
0.8
0.6
Y Coord
0.4 0.2
0.0
par(par.old)
Los valores en las posiciones {(0, 0), (0, 1), (1, 0), (1, 1)} coinciden con los gene-
rados anteriormente.
7.6. SIMULACIÓN BASADA EN CÓPULAS 123
365
360
355
350
Time
Además de las cópulas Gausianas, es una de las familias de cópulas más utiliza-
das. Son de la forma:
( d )
∑
−1
C(x1 , x2 , . . . , xd ) = Ψ Ψ (Fi (xi )) ,
i=1
Ejemplos:
• Cópula producto o independiente: Ψ(x) = − ln(x),
α
• Cópula de Gumbel: Ψ(x) = (− ln(x)) ; α ≥ 1
7.6.2 Simulación
2. Obtener V = CU−1 (W )
( )
3. Devolver F1−1 (U ), F2−1 (V )
Ejercicio 7.1.
( ( ) )− α1
Cu−1 (w) ≡ u−α w− α+1 − 1 + 1
α
,
diseñar una rutina que permita generar una muestra de tamaño n 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)
}
1.0
0.8
0.6
v
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)
7.6. SIMULACIÓN BASADA EN CÓPULAS 127
1.0
1.0
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)
0.4
0.2
0.0
y[,1]
128CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
1.0
0.8
0.6
y[,3]
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]
2
1
0
0 2 4 6 8
exp1
7.7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES DISCRETAS129
# 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]
par(par.old)
z <- 1:10
xy <- matrix(z, ncol = 2)
xy
## [,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
7.7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES DISCRETAS131
## buena 320 14
## mala 80 36
## 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
## buena 0.71111111 0.03111111
## mala 0.17777778 0.08000000
## [1] 1 2 3 4
Con probabilidades:
pz <- df$p
pz
as.vector(pij)
## 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
## 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
7.7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES DISCRETAS133
El código anterior puede ser empleado para simular tablas de contingencia. Aun-
que en estos casos se suele fijar el total de la tabla (o incluso las frecuencias
marginales). En este caso, sólo habría que fijar el nº de simulaciones al total de
la tabla:
nsim <- sum(n)
set.seed(1)
rz <- sample(z, nsim, replace = TRUE, prob = pz)
rtable <- table(rz) # Tabla de frecuencias unidimensional
matrix(rtable, ncol = 2) # Tabla de frecuencias bidimensional
## [,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
Consideraríamos como probabilidades:
pind <- (rowSums(tabla) %o% colSums(tabla))/(sum(tabla)^2)
matrix(pind, nrow = nrow(tabla))
## [,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
134CAPÍTULO 7. SIMULACIÓN DE DISTRIBUCIONES MULTIVARIANTES
##
## 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
7.8. EJERCICIOS PROPUESTOS 135
Ejercicio 7.4.
Aplicaciones de la
simulación en Inferencia
Estadística
137
138CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
1∑
n
µ̂ = X = Xi
n i=1
es: ( )
σ
X∼N µ, √
n
Confirmar este resultado mediante simulación, para ello:
a) Crear un conjunto de datos muestras con 500 muestras de tamaño n = 10
de una N (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:
8.1. DISTRIBUCIÓN EN EL MUESTREO 139
mux <- 1
sdx <- 2
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
es: ( )
σ σ
IC1−α (µ) = X − z1−α/2 √ , X + z1−α/2 √ .
n n
## [1] 480
8.2. INTERVALOS DE CONFIANZA 143
## [1] 96
100*(1 - alfa) # Proporción teórica bajo normalidad
## [1] 95
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
Como ejemplo ilustrativo, generamos el gráfico de los primeros 100 inter-
valos:
m <- 100
tmp <- muestras2[1:m,]
attach(tmp)
color <- ifelse(cob,"blue","red")
plot(1:m, mean, col = color, ylim = c(min(ici),max(ics)),
xlab = "Muestra", ylab = "IC")
arrows(1:m, ici, 1:m, ics, angle = 90, length = 0.05, code = 3, col = color)
abline(h = muexp, lty = 3)
detach(tmp)
2.5
2.0
1.5
IC
1.0
0.5
0.0
0 20 40 60 80 100
Muestra
1.4
1.2
1.0
0.8
Densidad
0.6
0.4
0.2
0.0
Medias
Ejercicio 8.4.
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)
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)
}
# Gráfico coberturas:
plot(p.teor, p.cov, type = "l", ylim = c(1 - 4 * alpha, 1))
abline(h = 1 - alpha, lty = 2)
8.2. INTERVALOS DE CONFIANZA 149
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:
150CAPÍ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
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))
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
Histogram of pvalor
4
3
Density
2
1
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
Valores iniciales:
8.3. CONTRASTES DE HIPÓTESIS 153
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
}
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)
0.6
0.4
0.2
0.0
Nivel de significación
Valores iniciales:
set.seed(54321)
nx <- 30
ratex <- 1
8.3. CONTRASTES DE HIPÓTESIS 155
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
## Proporción de rechazos al 10% = 0.008
Histogram of pvalor
2.0
1.5
Density
1.0
0.5
0.0
pvalor
156CAPÍ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
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
}
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)
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 )
}
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)
0.6
0.4
0.2
0.0
shape
0.2
0.1
0.0
0 5 10
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)
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")
8.5. REMUESTREO BOOTSTRAP 163
0.3
Error cuadrático
0.2
0.1
0.0
Media Mediana
8.5.1 Idea:
# Distribución bootstrap
curve(ecdf(data)(x), ylab = "FD", type = "s", lwd = 2)
# Distribución teórica
abline(a = 0, b = 1, lty = 2)
1.0
0.8
0.6
FD
0.4
0.2
0.0
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
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)
## [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)
8.5. REMUESTREO BOOTSTRAP 167
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)
2
1
0
stat.boot − stat.dat
168CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
## [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)
}
set.seed(1)
stat.dat <- fstatistic(dat)
stat.boot <- boot.strap(dat, nboot, fstatistic)
library(boot)
# ?boot
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)
170CAPÍTULO 8. APLICACIONES DE LA SIMULACIÓN EN INFERENCIA ESTADÍSTICA
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)
8.5. REMUESTREO BOOTSTRAP 171
0.2
* * * * * *
*** * **** *** * *** * *********** **** ** ******************** *** * * ** *
* ** * ** * * *
5, 10, 16, 50, 84, 90, 95 %−iles of (T*−t)
* * ** **** ** *
* * **** ********** *
********************* **** *
* ** ** * ** ** *
*
*** * **** *** * * *** * * ** ** * ** * * ** * * * * * * *
** * * ****** **** * *** * ****************** *********************** ** **
0.1
* ***** * ** * * * ** * *
0.0
* ** * *
** * * ****** *** * *** ** ************** ** ** * * * * * ** ** * *
***** ******* **** ** *** * **** * *** ** *
*** * ****** *** ** ** **************** **
************** ***** *** * * ***** * ** ** * ** ** *
* *
* * * * * * * *
** * * ***** *** * *** ****************** *********************** ** * ****** * *** ** * ** ** *
−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
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
Integración y Optimización
Montecarlo
173
174 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
1∑
n
I≈ h (xi ) (b − a)
n i=1
Ejercicio 9.1.
Crear una función que implemente la integración Monte Carlo clásica para apro-
ximar integrales del tipo:
∫ b
I= h(x)dx.
a
curve(fun, 0, 1)
abline(h = 0, lty = 2)
abline(v = c(0, 1), lty = 2)
set.seed(1)
mc.integral0(fun, 0, 1, 20)
## [1] 0.977663
9.1. INTEGRACIÓN MONTE CARLO (CLÁSICA) 175
3
fun(x)
2
1
0
mc.integral0(fun, 0, 1, 100)
## [1] 0.7311169
mc.integral0(fun, 0, 1, 100)
## [1] 0.8304299
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)))
}
176 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
set.seed(1)
mc.integral(fun, 0, 1, 5000)
## $valor
## [1] 0.8142206
##
## $error
## [1] 0.03087305
abline(h = 4/5, lty = 2)
1.1
1.0
Media y rango de error
0.9
0.8
0.7
0.6
0.5
Iteraciones
## $valor
## [1] 0.8142206
##
## $error
## [1] 0.0309005
1∑
n
θ≈ h (xi )
n i=1
## [1] 0.7967756
# error <- 2*sd(h(x))/sqrt(nsim)
error <- 2*sd(x)/sqrt(nsim)
error
## [1] 0.004728174
Ejercicio 9.2.
Aproximar: ∫ ∞
1 x2
ϕ(t) = √ e− 2 dx,
t 2π
para t = 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(1)
178 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
## [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
Como ya se comentó anteriormente, sería preferible generar más valores donde
hay mayor “área”, pero en este caso f concentra la densidad en una región
que no resulta de utilidad. Por ese motivo puede ser preferible recurrir a una
densidad auxiliar que solvente este problema.
puede ser preferible generar observaciones de una densidad g que tenga una
forma similar al producto hf .
Si Y ∼ g:
∫ ∫
h (x) f (x)
θ= h (x) f (x)dx = g(x)dx = E (q (Y )) .
g(x)
h(x)f (x)
siendo q (x) = g(x) .
Si y1 , y2 , . . . , yn i.i.d. Y ∼ g:
1∑ 1∑
n n
θ≈ q (yi ) = w(yi )h (yi ) = θ̂g
n i=1 n i=1
f (x)
con w (x) = g(x) .
La idea básica es que si la densidad g tiene colas más pesadas que la densidad f
con mayor facilidad puede dar lugar a “simulaciones”
( ) con varianza finita (podría
emplearse en casos en los que no existe E h2 (X) ; ver Sección 4.1 en el Tema
4 de Análisis de resultados).
La distribución de los pesos w(yi ) debería ser homogénea para evitar datos
influyentes (inestabilidad).
Ejercicio 9.3.
1.0e−05
5.0e−06
0.0e+00
x
180 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
## [1] 3.144811e-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)
5e−06
4e−06
Aproximación
3e−06
2e−06
1e−06
Iteraciones
## [1] 1.371154e-07
Empleando la aproximación tradicional:
est <- mean(rnorm(nsim) > 4.5)
est
## [1] 0
9.2. MUESTREO POR IMPORTANCIA 181
sqrt(est * (1 - est)/nsim)
## [1] 0
Ejercicio 9.4.
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 = "Iteraciones
abline(h = pcauchy(6) - pcauchy(2), lty = 2)
182 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
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))
1200
1000
800
600
400
200
0
9.2. MUESTREO POR IMPORTANCIA 183
∑
n
w(yi )h (yi )
i=1
θ≈ ∑
n .
w(yi )
i=1
Ejercicio 9.5.
Histogram of rx
0.3
Density
0.2
0.1
0.0
−4 −2 0 2
rx
Ejercicio 9.6.
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.30
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(...) {
input <- c(...)
output <- f(...)
## Hacer algo ...
points(input[1], input[2], col = "lightgray")
if(best.f > output) {
lines(rbind(best.par, input), lwd = 2, col = "blue")
best.f <<- output
best.par <<- input
# points(best.par[1], best.par[2], col = "blue", pch = 20)
# cat("par = ", best.par, "value = ", best.f, "\n")
}
## ... con output e input
return(output)
}
}
0 0 −550
−8
0 0
−70 −65 0
−1
−1350
−60
20
−1250
0
−1
−500 1 00 −1150
−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
50 00
−16
50 −8
50 00
−6 −7
−2
−550
−2 −1 0 1 2 3 4 5
µ1
menor probabilidad las regiones donde la función objetivo es mayor. Por ejemplo,
se podría pensar en generar valores de acuerdo a una densidad (tranformación
Boltzman-Gibbs):
g(x) ∝ exp (−f (x)/T ) ,
donde T > 0 es un parámetro (denominado temperatura) seleccionado de forma
que se garantize la integrabilidad.
Entre los métodos de optimización Monte Carlo podríamos destacar:
• Métodos con gradiente aleatorio.
• Temple simulado.
• Algoritmos genéticos.
• Montecarlo EM.
• …
9.4.1 Algoritmo:
Ejercicio 9.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 parameters.
# For "SANN" maxit gives the total number of function evaluations: there is no other stopping c
# Defaults to 10000.
res <- optim(starts[j, ], tee.optim2d(like), method = "SANN", control = list(temp = 100, maxit
points(res$par[1],res$par[2], pch = 19)
190 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
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
3 −650
−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 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
40
−165 0 00
0 00 50 −8
−2
−550 −6 −6
−2 −1 0 1 2 3 4 5
µ1
192 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
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
194 CAPÍTULO 9. INTEGRACIÓN Y OPTIMIZACIÓN MONTECARLO
Capítulo 10
Técnicas de reducción de la
varianza
• Variables antitéticas.
• Muestreo estratificado.
• Variables de control.
• Métodos de remuestreo.
• Condicionamiento.
• …
195
196 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
θ = E (Z)
Para aproximar: ∫ 1
I= h (x) dx,
0
a partir de x1 , x2 , . . . , xn i.i.d. U (0, 1). Podemos emplear:
( )
h (U ) + h (1 − U )
I=E
2
1 ∑n
≈ (h (xi ) + h (1 − xi )) .
2n i=1
Y = 2µ − X
10.2. VARIABLES ANTITÉTICAS 197
Crear una función que implemente la técnica de variables antitéticas para apro-
ximar integrales del tipo:
∫ b
I= h (x) dx.
a
Emplearla para aproximar:
( ) ∫ 2
U (0,2) 1 x
E e = e dx ≈ 3.194,
0 2
y representar gráficamente la aproximación en función de n. Función objetivo:
a <- 0; b <- 2
ftn <- function(x) return(exp(x)/(b-a))
curve(ftn, a, b, ylim=c(0,4))
abline(h=0,lty=2)
abline(v=c(a,b),lty=2)
4
3
ftn(x)
2
1
0
## [1] 3.194528
Para la aproximación por integración Monte Carlo podemos emplear la función
del capítulo anterior:
198 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
200 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 201
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 n de una U (0, 1), se
pueden generar l = nk observaciones (1 ≤ k ≤ n) de la forma:
( )
(j − 1) j
Uj1 , . . . , Ujl ∼ U , para j = 1, ..., k.
k k
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.
Ejemplo 10.1. Muestreo estratificado de una exponencial (libro Ricardo)
1. Para i = 1, 2, . . . , 10:
2. Generar Ui :
3. Devolver Xi = − ln Ui .
• V ar (Xi ) = 0.0214644 si i = 1, 2, 3, 4,
• V ar (Xi ) = 0.229504 si i = 5, 6, 7, 8, 9 y
• V ar (X10 ) = 1.
Como consecuencia:
( ) 1 ∑
10
V ar X = 2 V ar (Xi ) = 0.022338
10 i=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 203
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
Ejercicio 10.3.
204 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
Cov(X, Y )
α∗ = − ,
V ar(Y )
( )
con V ar(X ∗ ) = V ar(X) 1 − ρ2 (X, Y ) (lo que supone una reducción del
100ρ2 (X, Y ) %).
En la práctica normalmente α∗ no es conocida. Para estimarlo se puede realizar
ajuste lineal de X sobre Y (a partir de los datos simulados Xi e Yi , 1 ≤ i ≤ n):
SXY
• Si x̂ = β̂0 + β̂1 y es la recta ajustada, con β̂1 = y β̂0 = X − β̂1 Y , la
SY2
estimación sería:
α̂∗ = −β̂1
∗
• Si µY = 0 ⇒ θ̂ = X = β̂0 .
## [1] 3.194528
Aproximación clásica por simulación:
set.seed(54321)
nsim <- 1000
u <- runif(nsim, a, b)
expu <- exp(u)
mean(expu)
## [1] 3.182118
Con variable control:
plot(u, expu)
reg <- lm(expu ~ u)$coef
abline(reg, col='blue')
7
6
5
expu
4
3
2
1
## (Intercept)
## 3.204933
Lo siguiente ya no sería necesario:
expuc <- expu - reg[2]*(u-1)
mean(expuc)
## [1] 3.204933
Estimación del porcentaje de reducción en la varianza:
206 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
100*(var(expu)-var(expuc))/var(expu)
## [1] 93.91555
E (X) − E (Y ) = E (X − Y ) .
1∑
n
X −Y = (Xi − Yi )
n i=1
( ) 1
V ar X − Y = (V ar (X) + V ar (Y ))
n
y tendríamos que:
( ) 1 ∑
N
1
V ar X − Y = 2 V ar (Xi − Yi ) = V ar (Xi − Yi )
n i=1 n
1
= (V ar (Xi ) + V ar (Yi ) − 2Cov (Xi , Yi ))
n
1
≤ (V ar (Xi ) + V ar (Yi ))
n
## [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
Bibliografía complementaria
Libros
209
210 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
Gentle, J.E. (1998). Random number generation and Monte Carlo methods.
Springer-Verlag.
Hörmann, W. et al. (2004). Automatic Nonuniform Random Variate Generation.
Springer.
Knuth, D.E. (1969). The Art of Computer Programming. Volume 2. Addison-
Wesley.
Knuth, D.E. (2002). The Art of Computer Programming. Volume 2, third edition,
ninth printing. Addison-Wesley.
Law, A.M. y Kelton, W.D. (1991). Simulation, modeling and analysis. McGraw-
Hill.
Moeschlin, O., Grycko, E., Pohl, C. y Steinert, F. (1998). Experimental stochas-
tics. Springer-Verlag.
Nelson, R. (1995). Probability, stochastic processes, and queueing theory: the
mathematics of computer performance modelling. Springer-Verlag.
Pardo, L. y Valdés, T. (1987). Simulación. Aplicaciones prácticas a la empresa.
Díaz de Santos.
Robert, C.P. y G. Casella (2004). Monte Carlo statistical methods. Springer.
Artículos
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.
212 CAPÍTULO 10. TÉCNICAS DE REDUCCIÓN DE LA VARIANZA
Apéndice A
Enlaces
213
214 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 X1 , . . . , Xn (o los errores de un modelo) consituyen una muestra
aleatoria simple de una variable aleatoria X. 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 Fθ (x), 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.
217
218 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
H0 simple H0 compuesta
{ {
H0 : F = N (0, 1) H0 : F = N (µ, σ 2 )
H1 : F ̸= N (0, 1) H1 : F ̸= N (µ, σ 2 )
B.1.1 Histograma
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
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 221
1.0
0.8
distribution function
0.6
0.4
0.2
0.0
10 15 20 25 30
siendo x(i) los cuantiles observados y qi = F0−1 (pi ) los esperados2 bajo H0 .
Ejemplo:
qqnorm(datos)
qqline(datos)
2 Típicamente
{ (i−0.5)
}
pi = n
: i = 1, · · · , n .
222 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
B.1. MÉTODOS DE BONDAD DE AJUSTE 223
∑
k
(ni − ei )2
χ2 = ∼ χ2k−r−1 , si H0 cierta
i=1
ei aprox.
ni
(ni −ei )2
Clases observadas pi bajo H0 ei bajo H0 ei
(n1 −e1 )2
C1 n1 p1 e1 e1
.. .. .. .. ..
. . . . .
224 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
ni
(ni −ei )2
Clases observadas pi bajo H0 ei bajo H0 ei
(nk −ek )2
Ck nk pk ek
∑ ∑ ∑ ek
∑k
χ2 = i=1 (ni −e
2
i)
Total i ni = n i pi = 1 i ei =n ei
i=1
ei
##
## 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. n > 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 ei = n · pi 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).
B.1. MÉTODOS DE BONDAD DE AJUSTE 225
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 ei .
• Si la variable es ordinal (o numérica) debe juntarse la que causó el proble-
ma con una de las adyacentes.
Si la variable de interés es continua, una forma de garantizar que ei ≥ 5 consiste
en tomar un número de intervalos k ≤ ⌊n/5⌋ y de forma que sean equiprobables
pi = 1/k, considerando los puntos críticos xi/k de la distribución bajo H0 .
Por ejemplo, se podría emplear la siguiente función (que imita a las incluídas
en R):
#-------------------------------------------------------------------------------
# chisq.test.cont(x, distribution, nclasses, output, nestpar,...)
#-------------------------------------------------------------------------------
# Realiza el test ji-cuadrado de bondad de ajuste para una distribución continua
# discretizando en intervalos equiprobables.
# Parámetros:
# distribution = "norm","unif",etc
# nclasses = floor(length(x)/5)
# output = TRUE
# nestpar = 0= nº de parámetros estimados
# ... = parámetros distribución
# Ejemplo:
# chisq.test.cont(x, distribution="norm", nestpar=2, mean=mean(x), sd=sqrt((nx-1)/nx)*sd(x))
#-------------------------------------------------------------------------------
chisq.test.cont <- function(x, distribution = "norm", nclasses = floor(length(x)/5),
output = TRUE, nestpar = 0, ...) {
# Funciones distribución
q.distrib <- eval(parse(text = paste("q", distribution, sep = "")))
d.distrib <- eval(parse(text = paste("d", distribution, sep = "")))
# Puntos de corte
q <- q.distrib((1:(nclasses - 1))/nclasses, ...)
tol <- sqrt(.Machine$double.eps)
xbreaks <- c(min(x) - tol, q, max(x) + tol)
# Gráficos y frecuencias
if (output) {
xhist <- hist(x, breaks = xbreaks, freq = FALSE, lty = 2, border = "grey50")
curve(d.distrib(x, ...), add = TRUE)
} else {
xhist <- hist(x, breaks = xbreaks, plot = FALSE)
}
# Cálculo estadístico y p-valor
226 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
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
## 2 (14.49908,16.94725] 3 5.125 -0.9386680
## 3 (16.94725,18.77800] 4 5.125 -0.4969419
## 4 (18.77800,20.41732] 6 5.125 0.3865103
## 5 (20.41732,22.05663] 4 5.125 -0.4969419
## 6 (22.05663,23.88739] 8 5.125 1.2699625
## 7 (23.88739,26.33556] 4 5.125 -0.4969419
## 8 (26.33556,30.77000] 6 5.125 0.3865103
##
## Pearson's Chi-squared test
##
## data: datos
## X-squared = 3.6829, df = 5, p-value = 0.5959
tribución empírica Fn :
( ) i
Teniendo en cuenta que Fn X(i) = n:
{ }
i i−1
Dn = max − F0 (X(i) ), F0 (X(i) ) −
1≤i≤n n n
{ + −
}
= max Dn,i , Dn,i
1≤i≤n
p = P (Dn ≥ d) ≤ α.
##
## One-sample Kolmogorov-Smirnov test
##
## data: datos
## D = 0.13239, p-value = 0.4688
## alternative hypothesis: two-sided
Si H0 es compuesta, el procedimiento habitual es estimar los parámetros desco-
nocidos por máxima verosimilitud y emplear F̂0 en lugar de F0 . Sin embargo,
al proceder de esta forma es de esperar que F̂0 se aproxime más que F0 a la
distribución empírica, por lo que los cuantiles de la distribución de Dn pueden
ser demasiado conservativos (los p-valores tenderán a ser mayores de lo que de-
berían) y se tenderá a aceptar la hipótesis nula (puede ser preferible aproximar
B.2. DIAGNOSIS DE LA INDEPENDENCIA 229
##
## 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
Típicamente Cov(X1 , X2 ) > 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 H0 en contrastes).
# Matriz de covarianzas
t.dist <- as.matrix(dist(t))
t.cov <- exp(-t.dist)
# 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
3
2
1
0
x
−1
−2
Datos independientes
Datos dependientes
−3
## [1] 0.8067621
var(x2)
## [1] 0.1108155
En el caso de datos dependientes se produce una clara subestimación de la
varianza
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
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.
B.2. DIAGNOSIS DE LA INDEPENDENCIA 233
1.0
2
1.0
0.5
1
0.5
0.0
0
−0.5
−1
0.0
−1.0
−2
par(old.par)
0.5
1
X_t+1
X_t+1
X_t+1
0.0
0.5
−0.5
−1
0.0
−1.0
−2
par(old.par)
Ejemplo
234 APÉNDICE B. BONDAD DE AJUSTE Y ALEATORIEDAD
30
25
X_t+1
20
15
10
10 15 20 25 30
X_t
## [1] 0.01344127
B.2.4 El correlograma
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
• 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).
old.par <- par(mfrow = c(1, 3))
acf(x1, main = 'Independencia')
acf(x2, main = 'Dependencia positiva')
acf(x3, main = 'Dependencia negativa')
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)
2
|r(k)| < √
n
Ejemplo
acf(datos) # correlaciones
Series datos
1.0
0.8
0.6
0.4
ACF
0.2
0.0
−0.2
0 5 10 15
Lag
3 En algunos campos, como en estadística espacial, en lugar del covariograma se suele em-
Series x2
0.10
0.08
0.06
ACF (cov)
0.04
0.02
0.00
0 5 10 15 20
Lag
+ + + + − − − + + + − − + + + + + + − − −−
+ + ++ − − − + + + −− + + + + ++ − − −−
| {z } | {z } | {z } |{z} | {z } | {z }
1 2 3 4 5 6
##
## 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.
• 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
##
## 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.
donde ⌊u⌋ denota la parte entera de u. De esta forma se consigue una sucesión
de enteros aleatorios supuestamente independientes con distribución uniforme
en {1, . . . , K}.
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
B.3. CONTRASTES ESPECÍFICOS PARA GENERADORES ALEATORIOS241
K!
P (C = c) = S(d, c),
(K − c)!K d
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).
B.3. CONTRASTES ESPECÍFICOS PARA GENERADORES ALEATORIOS243
P (Q = 5) = 0.03840000, P (Q = 6) = 0.07680000,
P (Q = 7) = 0.09984000, P (Q = 8) = 0.10752000,
P (Q = 9) = 0.10450944, P (Q = 10) = 0.09547776,
P (Q = 11) = 0.08381645, P (Q = 12) = 0.07163904,
P (Q = 13) = 0.06011299, P (Q = 14) = 0.04979157,
P (Q = 15) = 0.04086200, P (Q = 16) = 0.03331007,
P (Q = 17) = 0.02702163, P (Q = 18) = 0.02184196,
P (Q = 19) = 0.01760857, P (Q ≥ 20) = 0.07144851.
Integración numérica
245
246 APÉNDICE C. INTEGRACIÓN NUMÉRICA
4
3
fun(x)
2
1
0
return(trapezoid.vec(f.vec, h))
}
trapezoid(fun, 0, 1, 20)
## [1] 0.8033325
El error en esta aproximación se corresponde con:
(b − a)3 ′′
f (ξ),
12n2
para algún a ≤ ξ ≤ b (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 máximo
3
′′
absoluto es (b−a)
12n2 maxa≤ξ≤b |f (ξ)|. En el caso general multidimensional sería
O(n− d ).
2
simpson(fun, 0, 1, 20)
248 APÉNDICE C. INTEGRACIÓN NUMÉRICA
## [1] 0.8000033
El máximo error (en el caso unidimensional) viene dado por la expresión:
(b − a)5
(4)
max f (ξ).
180n4 a≤ξ≤b
level = 1
S.old <- (b-a) * (fun(a) + fun(b))/2
S.new <- quadrature_internal(S.old, fun,
C.2. INTEGRACIÓN NUMÉRICA BIDIMENSIONAL 249
quadrature(fun, 0, 1)
## [1] 0.8
Fuente: r-blogger Guangchuang Yu
C.1.4 Comandos de R
require(MASS)
area(fun, 0, 1)
## [1] 0.8000001
.
Consideraremos como ejemplo:
∫ ∫
1 1 ( )
x2 − y 2 dxdy = 0
−1 −1
.
f2d <- function(x,y) x^2 - y^2
nx = 21
ny = 21
x <- seq(ax, bx, length = nx)
y <- seq(ay, by, length = ny)
z <- outer(x, y, f2d)
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
Error O(n− d ).
2
## [1] -8.881784e-18
252 APÉNDICE C. INTEGRACIÓN NUMÉRICA
C.2.3 Comandos de R
Fuente: tolstoy.newcastle.edu.au.
Alternativamente se podría emplear la función adaptIntegrate del paquete
cubature.