Estatística Computacional II

leg.ufpr.br/~walmes/ensino/EC2

1 Método da aceitação e rejeição

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:

  1. Ter um bom gerador de números uniformes.
  2. Ter um bom gerador de números de uma variável aleatória representada por uma distribuição \(D\), escolhida de tal maneira que existe uma constante \(M\) tal que a densidade de \(g\) que caracteriza a distribuição \(D\) satisfaz \(f(x) \leq M g(x)\) para todo \(x\) do domínio de \(f\).

O seguinte algoritmo permite gerar números aleatórios de uma distribuição de probabilidade caracterizada pela função densidade \(f\):

  1. Gerar \(y\) como sendo uma ocorrência da variável aleatória representada por \(D\).
  2. Gerar \(u\) como sendo uma ocorrência de uma uniforme padrão.
  3. Se \(u \leq f(y)/(M g(y))\) considerar que \(x = y\) é um valor da distribuição de probabilidade alvo cuja densidade é \(f\), caso contrário, descartar \(y\).
  4. Repetir até atingir o número de valores desejado \(n\).
#-----------------------------------------------------------------------
# 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."

2 Simular números gaussianos usando a Cauchy

#-----------------------------------------------------------------------
# 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

3 Ilustrando o funcionamento do método

#-----------------------------------------------------------------------
# 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)

4 Conteúdo complementar

Estatística Computacional II Prof. Walmes Zeviani