Universidade Federal do Paraná
Curso de Estatística
CE 089 - Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani

Aula 05

Tabela de conteúdo

• Geração de números aleatórios Gaussianos pelo método Box-Muller

Geração de números aleatórios Gaussianos pelo método Box-Muller

• Ler Gaussian integrals.
• Assitir A integral gaussiana.
• Ler Método jacobiano da transformação, página 32.
• Ler Box-Muller method

##-----------------------------------------------------------------------------
## Gerando (pares de) número Gaussianos a partir de (pares de) números
## uniformes com o método Box-Muller.

x1 <- runif(1000)
x2 <- runif(1000)

y <- c(sqrt(-2*log(x1))*sin(2*pi*x2),
       sqrt(-2*log(x1))*cos(2*pi*x2))

plot(ecdf(y))
curve(pnorm(x, 0, 1), add=TRUE, col=2)

##-----------------------------------------------------------------------------
## Forma alternativa que evita as funções sin() e cos(). Por outro lado
## tem o custo de testar se o ponto simulado está dentro de um
## círculo. Os número uniformes usados devem ser de -1 à 1.

n <- 0; N <- 5000;
u <- matrix(vector(mode="numeric", length=2*N), ncol=2)
x1 <- vector(mode="numeric", length=N)

while(n<N){
    us <- runif(2, -1, 1) ## u1 e u2
    r2 <- sum(us^2)       ## u1^2+u2^2
    if(r2<=1){
        n <- n+1
        u[n,] <- us ## Matriz com u1 e u2 nas colunas.
        x1[n] <- r2 ## Vetor com os valores de x1=u1^2+u2^2.
    }
}

## Para verificar que u1^2+u2^2 tem distribuição uniforme.
plot(ecdf(x1)); curve(punif(x), add=TRUE, col=2)

## Gera valores da normal.
## y <- sqrt(-2*log(x1))*(u[,2]*sqrt(x1)
y <- sqrt(-2*log(x1)/x1)*u[,2]

## Verifica a distribuição dos números gerados.
plot(ecdf(y))
curve(pnorm(x, 0, 1), add=TRUE, col=2)

##-----------------------------------------------------------------------------
## Qual procedimento custa mais caro?

proc1 <- function(N){
    x <- matrix(runif(2*N), ncol=2)
    y <- c(sqrt(-2*log(x[,1]))*sin(2*pi*x[,2]),
           sqrt(-2*log(x[,1]))*cos(2*pi*x[,2]))
    return(y)
}

proc1(N=10)

##  [1]  0.07496  0.60162 -1.11450 -1.07215 -0.26772 -1.38736  0.79206  0.01287  1.50679
## [10]  0.28034  0.50684  0.09619  0.24968 -0.64311 -0.51765  0.01357 -0.12955  0.24965
## [19]  1.31446  1.40407

proc2 <- function(N){
    n <- 0
    u <- matrix(vector(mode="numeric", length=2*N), ncol=2)
    x1 <- vector(mode="numeric", length=N)
    while(n<N){
        us <- runif(2, -1, 1) ## u1 e u2
        r2 <- sum(us^2)       ## u1^2+u2^2
        if(r2<=1){
            n <- n+1
            u[n,] <- us ## Matriz com u1 e u2 nas colunas.
            x1[n] <- r2 ## Vetor com os valores de x1=u1^2+u2^2.
        }
    }
    y <- sqrt(-2*log(x1)/x1)*u[,2]
}

proc2(N=10)

system.time(proc1(N=1000000))

##    user  system elapsed 
##   0.626   0.019   0.646

system.time(proc2(N=1000000))

##    user  system elapsed 
##  13.181   0.055  13.247