Estatística Computacional II
|
Quando se deseja gerar números aleatórios de distribuições de probabilidade, nem sempre é possível aplicar o método da transformação integral da probabilidade. Alguns situação são:
Se a função densidade de probabilidade for conhecida, \(f\), então é possível gerar números aleatórios dessa variável aleatória caracterizada por \(f\) pelo método da aceitação e rejeição. É necessário satisfazer dois requisitos:
O seguinte algoritmo permite gerar números aleatórios de uma distribuição de probabilidade caracterizada pela função densidade \(f\):
#-----------------------------------------------------------------------
# Seja X uma v.a. com f(x) = 1.5 * x^2, -1 < x < 1. Simular valores
# desta.
# Gráfico da f.d.p da v.a. X, f(x).
curve(1.5 * (x ^ 2), -1, 1)
# Tem integral 1?
integrate(function(x) 1.5 * (x^2), lower = -1, upper = 1)
## 1 with absolute error < 1.1e-14
#-----------------------------------------------------------------------
# Considere como g() a densidade de uma uniforme entre -1 e 1. Então se
# a base é 2, o altura deve ser 0.5 para ter produto 1, assim g(y) =
# 0.5, -1 < y < 1. Qual deve ser um valor de M para garantir que f(x) <=
# (M g(x)) para todo x dentro do [-1, 1]? O valor de M tem que ser 3,
# pois 3 * 0.5 <= sup_{x} f(x) = 1.5.
curve(1.5 * (x^2), -1, 1, col = 2)
curve(0.5 + 0 * x, add = TRUE, lty = 2)
curve(3 * 0.5 + 0 * x, add = TRUE, lty = 2, lwd = 2)
legend("bottomright",
legend = c("f(x)", "g(x)", "M g(x)"),
lty = c(1, 2, 2),
col = c(2, 1, 1),
lwd = c(1, 1, 2),
bty = "n")
# Criando os elementos necessários.
f <- function(x) 1.5 * x^2
g <- function(x) 0.5 + 0 * x
M <- 3
x <- NULL
# 1. Gerar y cuja densidade é g()
y <- runif(n = 1, -1, 1)
y
## [1] 0.9658441
# 2. Gerar u de uma uniforme padrão.
u <- runif(n = 1)
u
## [1] 0.6891094
# 3. Comparar e decidir.
r <- f(y)/(M * g(y))
r
## [1] 0.9328548
if (u < r) {
x <- y
print("u < r então valor aceito.")
} else {
print("u >= r então valor descartado.")
}
## [1] "u < r então valor aceito."
#-----------------------------------------------------------------------
# Gerar da gaussiana padrão usando a cauchy.
par(mfrow = c(1, 2))
# Gaussiana e Cauchy.
curve(dnorm(x), -4, 4)
curve(dcauchy(x), add = TRUE, lty = 2)
# Razão entre elas.
curve(dnorm(x)/dcauchy(x), -4, 4)
abline(v = c(-1, 1), lty = 2)
layout(1)
dnorm(1)/dcauchy(1)
## [1] 1.520347
M <- sqrt(2 * pi/exp(1))
M
## [1] 1.520347
curve(M * dcauchy(x), -4, 4,
lty = 2,
ylim = c(0, M * dcauchy(0)),
ylab = "Densidade")
curve(dnorm(x), add = TRUE)
legend("topright",
legend = c("f(x)", "M g(x)"),
lty = c(1, 2),
bty = "n")
#-----------------------------------------------------------------------
# Aplicando o método.
M <- sqrt(2 * pi/exp(1))
f <- function(x) dnorm(x, 0, 1)
g <- function(x) dcauchy(x, 0, 1)
# A inversa da F para a Cauchy é: F^{-1} = tan(pi * (1 - u)). No
# entanto, vamos usar o gerador pronto rcauchy().
# 1. Gerar y cuja densidade é g, uma Cauchy.
y <- rcauchy(n = 1)
y
## [1] -1.948804
# 2. Gerar u de uma uniforme padrão.
u <- runif(n = 1)
u
## [1] 0.8539964
# 3. Comparar e decidir.
r <- f(y)/(M * g(y))
r
## [1] 0.5922061
if (u < r) {
x <- y
print("u < r, então valor aceito.")
} else {
print("u >= r, então valor descartado.")
}
## [1] "u >= r, então valor descartado."
#-----------------------------------------------------------------------
# Aplicar dentro de um laço condicional para obter N valores.
n <- 1 # Contador de valores aceitos.
l <- 1 # Contador de ciclos.
N <- 999 # Total de número à gerar.
# Vetor vazio.
# x <- vector(mode = "numeric", length = N)
x <- numeric(N)
while (n < N) {
y <- rcauchy(n = 1)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
if (u < w) {
x[n] <- y
n <- n + 1
}
l <- l + 1
}
plot(ecdf(x))
curve(pnorm(x), add = TRUE, col = 2)
# Taxa de aceitação.
n/l # Observada.
## [1] 0.6537958
1/M # Teórica.
## [1] 0.6577446
#-----------------------------------------------------------------------
# Ilustrando.
n <- 1
N <- 25
curve(M * g(x), -6, 6, lty = 2, ylim = c(0, M * g(0)))
curve(f, add = TRUE)
while (n < N) {
y <- rcauchy(n = 1)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
if (u < w) {
n <- n + 1
rug(y)
points(y, u * M * g(y))
} else {
points(y, u * M * g(y), pch = 4)
}
Sys.sleep(1)
}
#-----------------------------------------------------------------------
# Gerar o gif.
library(animation)
loop <- function() {
M <- sqrt(2 * pi/exp(1))
f <- function(x) dnorm(x, 0, 1)
g <- function(x) dcauchy(x, 0, 1)
n <- 0L
l <- 1L
N <- 20L
q <- 0
set.seed(123)
x <- numeric(N)
while (n < N) {
curve(M * g(x), -6, 6,
lty = 2, ylim = c(0, M * g(0)),
ylab = NA,
xlab = NA,
n = 303)
curve(f, add = TRUE, n = 303)
legend("topleft",
legend = c("f(x)", "M g(x)"),
lty = c(1, 2),
bty = "n")
legend("topright",
title = "Decisão:",
legend = c("Aceitar", "Rejeitar"),
lty = 1,
col = c(3, 2),
bty = "n")
legend("right",
title = "Valor:",
legend = c("Aceito", "Rejeitado"),
pch = c(1, 4),
col = c(3, 2),
bty = "n")
y <- rcauchy(n = 1)
title(sub = substitute("Candidato: " * y == yy,
list(yy = round(y, 4))))
rug(y, lwd = 2, col = 1)
segments(y, 0, y, f(y), col = 3)
segments(y, f(y), y, M * g(y), col = 2)
segments(y - 0.2, 0, y + 0.2, 0, col = 1)
segments(y - 0.2, f(y), y + 0.2, f(y), col = 3)
segments(y - 0.2, M * g(y), y + 0.2, M * g(y), col = 2)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
d <- u < w
if (d) {
x[n] <- y
n <- n + 1L
l <- l + 1L
q <- n/l
points(y, u * M * g(y), pch = 1, col = 3)
} else {
l <- l + 1L
q <- n/l
points(y, u * M * g(y), pch = 4, col = 2)
}
mtext(side = 3,
text = substitute(u * ' = ' * uu <= r * ' = ' * rr * '?',
list(uu = round(u, 4),
rr = round(w, 4))))
msg <- sprintf("iterações: %d \t aceitos: %d \t taxa ac.: %0.2f",
l, n, q)
title(main = list(msg, font = 1))
# Sys.sleep(0.5)
} # END while(n < N)
} # END loop().
# Verifica.
# loop()
saveHTML(loop(),
img.name = "met-ac-rej",
imgdir = "met-ac-rej",
htmlfile = "met-ac-rej.html",
autobrowse = FALSE,
verbose = FALSE,
title = "Método da aceitação e rejeição",
ani.width = 600,
ani.height = 600)
Estatística Computacional II | Prof. Walmes Zeviani |