Universidade Federal do Paraná
Curso de Gradução em Estatística

CE 089 - Estatística Computacional II

http://www.leg.ufpr.br/ce089

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

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.


Uma proposta de gerador

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:

  • Temos um bom gerador de números aleatórios?
  • O que se pode dizer pela análise do gráfico 1?
  • O gráfico 2 une os pontos do gráfico 1 serialmente. Isso acrescenta alguma informação?
  • O gráfico 3 dá detalhes da série para os primeiros 50 valores. Existe algo curioso nessa série que seja indesejável? O que seria e o que a causa?
  • O gráfico 4 é o partial auto correlation function usado em séries temporais para verificar o grau de dependência serial que existe. O que esse gráfico está dizendo?
  • Considerando todos os aspectos, temos um bom gerador de números aleatórios?

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

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"