Estatística Computacional II

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

1 O método de Box-Müller

Assista ao vídeo A integral Gaussiana no perfil do Luiz Chamon para compreender a passagem de coordenadas cartesianas para coordenadas polares.

De acordo com o vídeo, a integral em \(\text{d}x \text{d}y\) pode ser escrita em coordenadas polares, \(\text{d}r \text{d}\theta\), \[ \begin{align*} I^2 = 2\pi &= \int_\mathbb{R} \int_\mathbb{R} \exp\{-(x^2 + y^2)/2\}\, \text{d}x \text{d}y\\ &= \int_{0}^{2\pi} \int_{0}^{\infty} r\exp\{-r^2/2\}\, \text{d}r \text{d}\theta \end{align*} \]

Para gerar números aleatórios da distribuição normal padrão, precisamos de valores para o raio \(r\) e o ângulo \(\theta\) de tal forma a poder convertê-los em valores \(x\) e \(y\) usando as expressões \[ \begin{cases} x = r \cos(\theta) \\ y = r \sin(\theta). \end{cases} \]

Valores para o ângulo podem ser obtidos pelo produto de uma \(\text{Uniforme}(0, 1)\) por \(2\pi\), ou seja \(\theta \sim 2\pi U(0, 1)\).

Para o raio, precisamos verificar qual a distribuição adequada. Note que a função de distribuição acumulada da variável \(R\) pode ser determinada a partir da integral \[ \begin{align*} \Pr(R \leq r) &= \int_{0}^{2\pi} \int_{0}^{r} (2\pi)^{-1} s\exp\{-s^2/2\}\, \text{d}s \text{d}\theta \\ &= (2\pi)^{-1} \int_{0}^{2\pi} \text{d}\theta \int_{0}^{r} s\exp\{-s^2/2\}\, \text{d}s \\ &= (2\pi)^{-1} (2\pi) \int_{0}^{r} s\exp\{-s^2/2\}\, \text{d}s\\ &= \int_{0}^{r} s\exp\{-s^2/2\}\, \text{d}s. \end{align*} \]

Para resolver a integral, considere \(u = s^2\) e com isso \(\text{d}s = \text{d}u/2s\). Os limites de integração são alterados: quando \(s = 0\) tem-se que \(u = s^2 = 0\); quando \(s = r\) tem-se que \(u = s^2 = r^2\). Dessa forma, a integral fica \[ \begin{align*} \Pr(R \leq r) &= \int_{0}^{r} s\exp\{-s^2/2\}\, \text{d}s\\ &= \int_{0}^{r^2} s\exp\{-u/2\}\, (2s)^{-1}\text{d}u\\ &= \frac{1}{2} \int_{0}^{r^2} \exp\{-u/2\}\, \text{d}u\\ &= \frac{1}{2} (-2) (\exp\{-u/2\}) \biggr\rvert_{0}^{r^2}\\ &= -(\exp\{-u/2\}) \biggr\rvert_{0}^{r^2}\\ &= -(\exp\{-r^2/2\} - \exp\{-0^2/2\})\\ &= 1 - \exp\{-r^2/2\}. \end{align*} \]

O resultado obtido foi a função de distribuição acumulada da variável aleatória \(R\), ou seja, \(F(r)\). Se formos capazes de inverter essa função, poderemos gerar números aleatórios de \(R\) usando números uniformes.

A inversa da função \(F(r) = 1 - \exp\{-r^2/2\}\) é \[ \begin{align*} u &= 1 - \exp\{-r^2/2\}\\ \log(1 - u) &= -r^2/2\\ r &= \pm \sqrt{-2 \log(1 - u)} \end{align*} \]

Como os valores de \(u\) são realizações de uma uniforme padrão, podemos simplificar esse resultado para \[ r = \sqrt{-2 \log(u)} \]

Assim, para gerar números da normal padrão, usamos \[ \begin{cases} x = \sqrt{-2 \log(u_1)}\, \cos(2\pi u_2)\\ y = \sqrt{-2 \log(u_1)}\, \sin(2\pi u_2), \end{cases} \] em que \(u_1\) e \(u_2\) são números da Uniforme padrão.

2 Gerando números da Normal

Utilizando os resultados obtidos, geramos números aleatórios com a implementação imediata dos mesmos.

# Números uniformes.
n <- 1000
u1 <- runif(n)
u2 <- runif(n)

# Funções dos números uniformes.
a <- sqrt(-2 * log(u1))
b <- 2 * pi * u2

# Números da normal padrão.
x1 <- a * cos(b)
x2 <- a * sin(b)

# Verificação.
qqnorm(c(x1, x2))

FERREIRA (2013) aponta que pode-se evitar o uso das funções trigonométricas de seno e cosseno usando números uniformes dentro do círculo trigonométrico de raio unitário. No entanto, essa solução deixa de ser vetorizável e pode não mais eficiente do que a versão vetorizável que usa as funções trigonométricas.

#--------------------------------------------
# Evitando as funções trigonométricas (mas deixando de ser vetorial).

# Cria vetores para serem prenchidos.
n <- 1000
w1 <- numeric(n)
w2 <- numeric(n)
v2 <- numeric(n)

i <- 1
while (i <= n) {
    # Gera os números e obtém o raio.
    u1 <- runif(1, -1, 1)
    u2 <- runif(1, -1, 1)
    r2 <- u1^2 + u2^2
    # Testa se os valores são admissíveis e incrementa contador.
    if (r2 < 1 & r2 != 0) {
        w1[i] <- u1
        w2[i] <- u2
        v2[i] <- r2
        i <- i + 1
    }
}

a <- sqrt(-2 * log(v2)/v2)
x1 <- a * w1
x2 <- a * w2

# Verificação.
qqnorm(c(x1, x2))

3 Comparando as implementações

rnorm_vec <- function(n) {
    m <- ceiling(n/2)
    u1 <- runif(m)
    u2 <- runif(m)
    a <- sqrt(-2 * log(u1))
    b <- 2 * pi * u2
    x <- c(a * cos(b), a * sin(b))
    if (n %% 2) x else x[-1]
}

rnorm_loop <- function(n) {
    m <- ceiling(n/2)
    w1 <- numeric(m)
    w2 <- numeric(m)
    v2 <- numeric(m)
    i <- 1
    while (i <= m) {
        u1 <- runif(1, -1, 1)
        u2 <- runif(1, -1, 1)
        r2 <- u1^2 + u2^2
        if (r2 < 1 & r2 != 0) {
            w1[i] <- u1
            w2[i] <- u2
            v2[i] <- r2
            i <- i + 1
        }
    }
    a <- sqrt(-2 * log(v2)/v2)
    x <- c(a * w1, a * w2)
    if (n %% 2) x else x[-1]
}

library(rbenchmark)

n <- 1000
benchmark(rnorm_vec(n = n),
          rnorm_loop(n = n),
          replications = 500)
##                test replications elapsed relative user.self sys.self
## 2 rnorm_loop(n = n)          500   1.079    22.02     1.076        0
## 1  rnorm_vec(n = n)          500   0.049     1.00     0.048        0
##   user.child sys.child
## 2          0         0
## 1          0         0
n <- 2
benchmark(rnorm_vec(n = n),
          rnorm_loop(n = n),
          replications = 5000)
##                test replications elapsed relative user.self sys.self
## 2 rnorm_loop(n = n)         5000   0.042      1.4     0.040        0
## 1  rnorm_vec(n = n)         5000   0.030      1.0     0.032        0
##   user.child sys.child
## 2          0         0
## 1          0         0

Esse resultado mostra que, embora a primeira implementação use funções trigonométricas e isso seja pouco eficiente, o fato de executar vetorialmente compensa tal custo até mesmo para obter apenas um par de números. Quando maior é a quantidade de números, maior vantagem da implementação vetorial.

4 Próxima aula

Referências

FERREIRA, D. F. Estatística computacional em java. Editora UFLA, 2013.

Estatística Computacional II Prof. Walmes Zeviani