Universidade Federal do Paraná
Curso de Estatística
CE 089 -
Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
##-----------------------------------------------------------------------------
## 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