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
Ler e fazer
Ler
##-----------------------------------------------------------------------------
## Método da congruência para gerar números uniformes.
x <- NULL
## Escolha errada dos parâmetros.
n <- 100; a <- 2; m <- 10; x[1] <- 3
for(i in 1:n){
x[i+1] <- (a*x[i])%%m
## print(x[i+1])
}
## O conjunto resultado é regular e não aleatório.
plot(x, type="o")
## Escolha correta (ou mais adequada).
a <- 16807; m <- 2^31-1
x[1] <- 33475890
for(i in 1:n){
x[i+1] <- (a*x[i])%%m
## print(x[i+1])
}
## O conjunto agora é aleatório.
plot(x, type="o")
## Passando para escala [0,1).
round(x/m, 3)
## [1] 0.016 0.995 0.527 0.649 0.365 0.435 0.636 0.397 0.587 0.759 0.568 0.361 0.336 0.520
## [15] 0.169 0.814 0.405 0.955 0.635 0.149 0.391 0.557 0.774 0.309 0.749 0.748 0.796 0.567
## [29] 0.251 0.710 0.209 0.986 0.735 0.546 0.840 0.199 0.848 0.108 0.223 0.504 0.154 0.560
## [43] 0.301 0.084 0.063 0.141 0.622 0.343 0.518 0.849 0.210 0.870 0.015 0.787 0.146 0.928
## [57] 0.725 0.397 0.106 0.398 0.767 0.407 0.927 0.983 0.645 0.811 0.169 0.248 0.408 0.778
## [71] 0.929 0.627 0.292 0.867 0.168 0.479 0.347 0.427 0.930 0.674 0.685 0.883 0.523 0.422
## [85] 0.933 0.475 0.208 0.856 0.265 0.017 0.681 0.099 0.922 0.831 0.996 0.266 0.648 0.546
## [99] 0.600 0.381 0.458
##-----------------------------------------------------------------------------
## Função que gera números uniformes pelo método da congruência.
rand <- function(n, seed=NULL){
a <- 16807; m <- 2^31-1
x <- vector(length=n+1, mode="integer")
if(is.null(seed)){
x[1] <- as.integer(Sys.time())
} else {
x[1] <- as.integer(seed)
}
for(i in 1:n){
x[i+1] <- (a*x[i])%%m
}
u <- x[-1]/m
return(u)
}
u <- rand(500)
par(mfrow=c(1,2))
hist(u, prob=TRUE); abline(h=1, col=2)
plot(ecdf(u), pch=NA); abline(a=0, b=1, col=2); layout(1)
##-----------------------------------------------------------------------------
## Para gerar variáveis aleatórias discretas limitadas.
## Dicotomizando em cara (se < 0.5) e coroa (se >= 0.5).
as.integer(u<0.5)
## [1] 0 1 1 0 0 0 0 0 0 0 1 1 0 0 1 1 0 0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 1 1 1 0 0 0 1 0 1 0
## [43] 1 1 0 1 0 0 0 0 1 0 1 0 1 1 0 0 1 1 1 1 1 1 1 0 1 1 0 0 1 1 1 1 1 1 0 1 0 0 0 1 0 1
## [85] 1 1 0 0 1 1 0 0 0 1 0 1 0 1 1 1 0 1 1 0 0 1 1 1 1 1 1 0 1 0 1 1 1 1 0 0 0 0 1 0 0 0
## [127] 0 1 0 1 0 0 1 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1 1 1 0 1 0 0
## [169] 0 0 0 1 0 1 1 1 1 1 1 1 0 0 0 0 1 0 1 0 1 1 1 1 0 1 1 1 1 0 0 1 0 1 1 1 0 1 1 1 1 0
## [211] 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 1 0 1 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 1
## [253] 0 0 0 0 1 0 0 0 0 1 0 1 0 1 1 1 1 1 1 1 0 1 0 1 0 1 1 1 0 0 1 1 1 0 0 1 1 0 1 0 1 1
## [295] 0 0 1 1 0 1 1 1 0 0 1 0 0 1 1 0 1 0 1 1 1 0 0 0 1 1 0 1 1 0 1 1 1 0 1 0 1 1 1 1 1 0
## [337] 1 1 0 1 0 0 0 1 1 1 0 1 0 1 0 1 1 1 0 0 1 0 0 0 0 1 1 1 1 1 0 1 0 1 0 0 1 0 0 1 0 0
## [379] 0 0 1 1 1 0 0 1 0 1 1 0 0 1 0 1 1 0 1 0 1 0 0 1 1 1 0 1 0 1 0 1 0 0 0 1 1 1 1 0 0 1
## [421] 0 1 1 0 0 1 0 1 1 1 1 1 1 0 0 1 0 1 1 1 0 0 1 1 0 1 1 0 1 1 1 0 1 0 1 0 1 0 1 0 0 1
## [463] 1 1 0 0 0 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 1 0 1 0 1 1 0 0 0 1 1 0 1 1 1
## Classificando para representar o lançamento de um dado.
clas <- seq(0, 1, length=7)
u <- rand(10000)
x <- findInterval(u, vec=clas)
tb <- table(x)
prop.table(tb)
## x
## 1 2 3 4 5 6
## 0.1682 0.1657 0.1649 0.1691 0.1689 0.1632
##-----------------------------------------------------------------------------
## Gerando visualização.
u <- rand(20)
clas <- seq(0, 1, length=5)
x <- findInterval(u, vec=clas)
plot(ecdf(seq_len(length(clas)-1)), main=NULL)
for(i in seq_along(u)){
rug(u[i], side=2, col=x[i])
points(x[i], u[i], col=x[i])
arrows(0, u[i], x[i], u[i], length=0.05, col=x[i], lty=2)
## Sys.sleep(time=0.3)
}
##-----------------------------------------------------------------------------
## Gerando um gif.
png(filename="fig/randDiscret%03d.png", family="Ubuntu Light",
type="cairo", res=150, width=1050, height=800)
for(i in seq_along(u)){
plot(ecdf(seq_len(length(clas)-1)), main=NULL)
rug(u[1:i], side=2)
points(x[1:i], u[1:i], col=x[1:i])
arrows(0, u[1:i], x[1:i], u[1:i], length=0.05, col=x[1:i], lty=2)
}
dev.off()
f <- list.files(path="./fig", pattern="randDsicret.*\\.png")
system("convert -delay 20 ./fig/randDiscret*.png ./fig/randDiscret.gif")
file.remove(paste0("./fig/", f))
##-----------------------------------------------------------------------------
## Gerando número aleatórios de outras distribuições pelo método da
## transformação integral da probabilidade.
u <- rand(10000)
## Números com distribuição exponencial.
lambda <- 2
x <- -log(1-u)/lambda
plot(ecdf(x))
curve((1-exp(-lambda*x))*(x>0), add=TRUE, col=2)
legend("right", legend=c("Empírica", "Teórica"),
col=1:2, lty=1, bty="n")
##-----------------------------------------------------------------------------
## Gerando visualização.
u <- rand(20)
x <- -log(1-u)/lambda
curve((1-exp(-lambda*x))*(x>0), 0, 1.25*max(x))
for(i in seq_along(u)){
rug(u[i], side=2)
points(x[i], u[i])
arrows(0, u[i], x[i], u[i], length=0.05, lty=2)
Sys.sleep(time=0.25)
arrows(x[i], u[i], x[i], 0, length=0.05, lty=2)
rug(x[i])
## Sys.sleep(time=0.25)
}
##-----------------------------------------------------------------------------
## Gerando um gif.
png(filename="fig/randExp%03d.png", family="Ubuntu Light",
type="cairo", res=150, width=1050, height=800)
for(i in seq_along(u)){
curve((1-exp(-lambda*x))*(x>0), 0, 1.25*max(x))
rug(u[1:i], side=2)
points(x[1:i], u[1:i])
arrows(0, u[1:i], x[1:i], u[1:i], length=0.05, lty=2)
Sys.sleep(time=0.25)
arrows(x[1:i], u[1:i], x[1:i], 0, length=0.05, lty=2)
rug(x[1:i])
Sys.sleep(time=0.25)
}
dev.off()
f <- list.files(path="./fig", pattern="randExp.*\\.png")
system("convert -delay 20 ./fig/randExp*.png ./fig/randExp.gif")
file.remove(paste0("./fig/", f))
##-----------------------------------------------------------------------------
## Como gerar números da distribuição Cauchy?
## browseURL("http://en.wikipedia.org/wiki/Cauchy_distribution")
F <- function(x) (1/pi)*atan(x)+(1/2)
curve(F, -3, 3)
x <- tan(pi*(1-rand(1000)))
plot(ecdf(x), xlim=c(-4,4))
curve(F, add=TRUE, col=2)
x <- rnorm(1000)/rnorm(1000)
plot(ecdf(x), xlim=c(-4, 4))
curve(F, add=TRUE, col=2)
##-----------------------------------------------------------------------------
## Como gerar números de uma distribuição Poisson? Usar a relação que
## tem com a exponencial.
##-----------------------------------------------------------------------------
## Usando recursos númericos para simular da normal via método da
## transformação integral da probabilidade.
x <- seq(-3,3,l=6)
Fx <- pnorm(x)
plot(Fx~x, type="l")
curve(pnorm(x), add=TRUE, col=2)
myFx <- approxfun(x=x, y=Fx)
xi <- seq(-3,3, by=0.1)
Fxi <- myFx(xi)
## plot(Fxi~xi)
myiFx <- approxfun(x=Fx, y=x)
u <- rand(10000)
xrand <- myiFx(u)
plot(ecdf(xrand))
curve(pnorm(x), add=TRUE, col=2)
##-----------------------------------------------------------------------------
## Gerando números normais usando como apoio o Teorema do Limite
## Central. Aqui fazemos a média de 6 números aleatórios que têm
## distribuição aproximadamente normal.
U <- matrix(rand(60000), ncol=6)
x <- apply(U, 1, mean)
qqnorm(x)
##-----------------------------------------------------------------------------
## Mudando o nível de aproximação da função densidade acumulada.
png(filename="fig/aprox%03d.png", family="Ubuntu Light",
type="cairo", res=150, width=1050, height=800)
for(i in 3:15){
x <- seq(-3,3,l=i)
Fx <- pnorm(x)
curve(pnorm(x), -3, 3, col=2)
lines(Fx~x, type="o")
## Sys.sleep(time=0.5)
}
dev.off()
f <- list.files(path="./fig", pattern="aprox.*\\.png")
system("convert -delay 10 ./fig/aprox*.png ./fig/aprox.gif")
file.remove(paste0("./fig/", f))
