Universidade Federal do Paraná
Curso de Gradução em Estatística
Prof. Dr. Walmes M. Zeviani
Lab. de Estatística e Geoinformação - LEG
Departamento de Estatística - UFPR
Geração de números aleatórios é ramo importante da Estatística Computacional. A geração de números aleatórios permitem fazer de sorteios de brindes até produzir jogos de azar mas não é somente isso, é claro.
A geração de números aleatórios permite simulação variáveis aleatórias e, assim, de modelos estocásticos. Permite a simulação de dados com propriedades conhecidas, por exemplo, provenientes de determinadas distribuições de probabilidade para avaliar o desempenho de um teste para detecção de outliers ou o grau com o que o afastamento da normalidade influencia o desempenho do teste t. Com a simulação, pode-se investigar tamanhos de amostra para sitauções complexas que não dispõem de abordagem teórica. Dá para se avaliar níveis níveis cobertura de intervalos de confiança ou grau de vício de estimadores, bem como acessar a qualidade de um método de classificação e a velocidade de processamento de certa tarefa.
A simulação de variaáveis aleatórias pode ser feita de várias formas. A mais conhecida e geral é por meio da transformação integral da probabilidade que estabelece que variáveis aleatórias podem ser geradas a partir de variáveis aleatórias uniformes (uniforme contínua no intervalo [0, 1]) usando a inversa da função de distribuição acumulada (\(F^{-1}\)).
Antes de simularmos v.a. de distribuições de probabilidade, vamos primeiro simular números aleatórios e definir o que é uma sequência aleatória de valores.
Vamos começar com a tentativas obter um gerador de números aleatórios pelo emprego recursivo da função \(sin(x)\). Tomaremos cuidado de fazê-la ter domínio e contra domínio no intervalo unitário para que seja usada recursivamente sem problemas.
##----------------------------------------------------------------------
## Uma função para geração de números aleatórios no intervalo (0, 1).
## Divide para uma disposição 2x2.
par(mfrow=c(2,2))
## Domínio: [0, 2*pi].
## Imagem: [-1, 1].
curve(sin(x), from=0, to=2*pi)
## Domínio: [0, 1].
## Imagem: [-1, 1].
curve(sin(2*pi*x), from=0, to=1)
## Domínio: [0, 1].
## Imagem: [0, 1].
curve(0.5+sin(2*pi*x)/2, from=0, to=1,
ylim=c(0, 1), asp=1)
abline(v=0:1, h=0:1, lty=2, col=2)
layout(1)
A função a ser usada é a \(0.5+0.5\sin(2\pi x)\). Para gerar números aleatórios, vamos fazer uso recursivo dessa função \[ x_{i+1} = 0.5+0.5\sin(2\pi x_{i}). \]
Um valor inicial (x0
) é necessário para começar o processo. Sempre que o mesmo valor for usado, a mesma sequência será gerada. Daí que muitos preferem usar o termo pseudo aleatório para geradores assim. Embora a sequência gerada tenha propriedades aleatórias, na realidade é possível prever o próximo valor se conhecer o mecanísmo usado.
##----------------------------------------------------------------------
rand0 <- function(n, x0){
##-------------------------------------------
## Interrompe se x0 estiver fora do intervalo.
##
stopifnot(x0 > 0 & x0 < 1)
n <- n+1
##-------------------------------------------
## Cria um vetor vazio para receber os números.
##
x <- vector(mode="numeric", length=n)
x[1] <- x0
for (i in 2:n){
x[i] <- 0.5+sin(2*pi*x[i-1])/2
}
##-------------------------------------------
## Retorna a série, exceto o número inicial.
return(x[-1])
}
## Usando o gerador.
x <- rand0(n=1000, x0=0.3) ## com 5*pi
## Gráfico da série de valores.
par(mfrow=c(2,2))
plot(x, main=1)
plot(x, type="l", main=2)
plot(x[1:50], type="l", main=3)
pacf(x, main=4)
layout(1)
Interprete os gráficos e responda:
A seguir, como usar esses números para gerar os correspondentes ao lançamento de uma moeda justa, cara e coroa. E como gerar os correspondentes ao lançamento de um dado justo.
##----------------------------------------------------------------------
hist(x)
abline(v=0.5, col=2)
## Frequência de cara/coroa.
table(x<0.5)
FALSE TRUE
500 500
hist(x)
abline(v=(1:5)/6, col=2)
## Frequência das faces do dado.
table(cut(x, breaks=(0:6)/6))
(0,0.167] (0.167,0.333] (0.333,0.5] (0.5,0.667] (0.667,0.833]
180 104 216 216 102
(0.833,1]
182
##----------------------------------------------------------------------
## A melhor forma de verificar se a distribuição é uniforme, é pelo
## gráfico de distribuição de probabilidades acumulada.
plot(ecdf(x))
abline(v=c(0.1, 0.4), lty=2)
## Acima tem-se a delimitação de uma região aproximadamente distribuição
## uniforme.
u <- x[x>0.1 & x<0.4]
plot(ecdf(u))
## Só que quantos números restaram?
length(u)
[1] 203
## Poderia usar o outro lado também. Lembrar de subtrair.
plot(ecdf(x))
abline(v=c(0.6, 0.9), lty=2)
v <- x[x>0.6 & x<0.9]
u <- c(u, v-0.5)
plot(ecdf(u))
## Só que quantos números restaram?
length(u)
[1] 402
## Isso indica que para ter o número desejado de valores aleatório, deve
## se considerar essa fração de descarte. Isso implica em mais custo.
Considerando que o proveitamento é de 40%, 60% é descarte. Quantos valores simular para ter 1000 números válidos? O que fazer para que esses números entre 0.1 e 0.4 se tornem números entre 0 e 1?
##----------------------------------------------------------------------
## Visualizando.
x <- x[1:40]
tim <- 0.05
for (i in 2:length(x)){
plot(x[2:i]~x[1:(i-1)], xlim=c(0, 1), ylim=c(0,1), asp=1)
curve(0.5+sin(2*pi*x)/2, 0, 1, col="gray", add=TRUE)
arrows(x[i-1], 0, x[i-1], x[i], length=0.1); Sys.sleep(tim)
arrows(x[i-1], x[i], 0, x[i], length=0.1); Sys.sleep(tim)
arrows(0, x[i], x[i], 0, length=0.1); Sys.sleep(tim)
Sys.sleep(0.2)
}
plot(x)
## Função que quando chamada reproduz o for.
fun <- function(){
for (i in 2:length(x)){
plot(x[2:i]~x[1:(i-1)], xlim=c(0, 1), ylim=c(0,1), asp=1)
curve(0.5+sin(2*pi*x)/2, 0, 1, col="gray", add=TRUE)
arrows(x[i-1], 0, x[i-1], x[i], length=0.1)
arrows(x[i-1], x[i], 0, x[i], length=0.1)
arrows(0, x[i], x[i], 0, length=0.1)
}
}
##----------------------------------------------------------------------
## Animação.
library(animation)
saveGIF(fun(),
movie.name="gna.gif",
ani.width=400,
ani.height=400)
saveHTML(fun(),
img.name="gna", imgdir="gna",
interval=0.05,
htmlfile="gna.html",
verbose=FALSE)
O método da congruência tambpem é um recursivo. Esse método é reconhecido por ser um bm gerador de números uniformes. A recusividade é \[ x_{i+1} = (a x_{i}+c) \text{ mod } m \] em que todos os termos são números inteiros. Acontece que o método da conguência depende da escolha das constantes \(a\) e \(m\). O valor de \(c\) geralmente é zero. Essa fórmula recursiva vai produzir números aleatórios (inteiros) entre 0 e \(m\). Para que estejam no intervalo unitário basta dividir por \(m\) ao final.
Do que foi dito, você já percebe que quanto maior o valor de \(m\), mais números serão permitidos. Se queremos uma v.a. contínua, quanto mais valores melhor. Vamos ver como funciona com algumas escolhar arbitrárias de \(a\) e \(m\).
##----------------------------------------------------------------------
## Método da congruência.
rand1 <- function(n, x0, A=5, M=50, divide=TRUE){
stopifnot(x0 < M)
n <- n+1
x <- vector(mode="numeric", length=n)
x[1] <- x0
for (i in 2:n){
x[i] <- (A*x[i-1])%%M
}
if (divide){
return(x[-1]/M)
} else {
return(x[-1])
}
}
x <- rand1(n=100, x0=11, A=6, M=50, divide=FALSE)
## plot(x, type="o", ylim=c(0, 50))
Vemos aí uma má escolha de \(a\) e \(m\), certo? Explique porque é uma escolha ruim. Que tipos de números devem ser usados para \(a\) e \(m\). Estude a função e dê um palmite. Abaixo uma escolha mais apropriada.
##----------------------------------------------------------------------
## Uma escolha melhor para 'a' e 'm'.
x <- rand1(n=100, x0=73, A=16507, M=2^32-1)
plot(x)
plot(ecdf(x))
Faça gráfico e verifique se este método de fato gera uma sequência aleatória de números uniformes. Tente encontrar melhores valores para \(a\) e \(m\).
options(width=90)
devtools::session_info()
Sys.time()
setting value
version R version 3.2.2 (2015-08-14)
system x86_64, linux-gnu
ui X11
language en_US
collate en_US.UTF-8
tz Brazil/East
date 2015-10-15
package * version date source
devtools 1.9.1 2015-09-11 CRAN (R 3.2.2)
digest 0.6.8 2014-12-31 CRAN (R 3.2.2)
evaluate 0.8 2015-09-18 CRAN (R 3.2.2)
formatR 1.2.1 2015-09-18 CRAN (R 3.2.2)
htmltools 0.2.6 2014-09-08 CRAN (R 3.2.0)
knitr * 1.11 2015-08-14 CRAN (R 3.2.2)
lattice * 0.20-33 2015-07-14 CRAN (R 3.2.1)
latticeExtra * 0.6-26 2013-08-15 CRAN (R 3.1.0)
magrittr 1.5 2014-11-22 CRAN (R 3.2.2)
memoise 0.2.1 2014-04-22 CRAN (R 3.1.0)
RColorBrewer * 1.1-2 2014-12-07 CRAN (R 3.2.2)
rmarkdown * 0.7 2015-06-13 CRAN (R 3.2.1)
stringi 0.5-5 2015-06-29 CRAN (R 3.2.2)
stringr 1.0.0 2015-04-30 CRAN (R 3.2.2)
yaml 2.1.13 2014-06-12 CRAN (R 3.1.0)
[1] "2015-10-15 18:58:09 BRT"