Motivação

A qualidade de uma simulação depende da qualidade dos números pseudoaleatórios utilizados.

Um bloco central para a simulação destes valores é a habilidade de gerar valores aleatórios uniformemente distribuídos no intervalo unitário.

Robert e Casella (2004) definem um gerador de números pseudoaleatórios uniforme como um algoritmo que, a partir de um valor inicial \(x_0\) e de uma transformação \(D\), produz uma sequência \((u_i) = (D(u_{i-1}))\) de valores no intervalo \((0, 1)\). Para todo \(n\), os valores \((u_1, ... , u_n)\) reproduzem o comportamento de uma amostra i.i.d. \((V_1, ... , V_n)\) de variáveis aleatórias uniformes quando comparados por meio de um conjunto habitual de testes.

Nesta aula, vamos nos concentrar neste problema através de uma classe particular de algoritmos geradores de números pseudoaleatórios chamados geradores congruenciais.

Geradores Congruenciais Multiplicativos

Definição

Um Gerador Congruencial Multiplicativo calcula valores sucessivos \(x_n\), \(n \geq 1\), de forma recursiva

\[ x_{n+1} = (a \; x_n) \mod m, \] a partir de um valor inicial \(x_0\), fixando \(a\) e \(m\) inteiros positivos. Os valores fixados são chamados de - \(x_0\): semente; - \(a\): multiplicador; e - \(m\): módulo.

De acordo com esta definição:

\(x_n \in \{0,1,2,...,m-1\}\); e
os números gerados

\[ u_n = \frac{x_n}{m}, \] chamados de números pseudoaleatórios, são tomados como um valor aleatório gerado de uma distribuição uniforme no intervalo \((0,1)\).

Exemplo 1

Encontrar a sequência de números pseudoaleatórios \(x_1,x_2,...,x_{10}\) supondo \[x_0=5 \qquad \mbox{ e } \qquad x_n = 3x_{n-1} \mod 150.\]

Solução: \[\begin{align*} \begin{array}{cccccccc} x_0 &=& 5 \\ x_1 &=& 3\;.\;5 &\mod 150 &=& 15 &\mod 150 &=& 15 \\ x_2 &=& 3\;.\;15 &\mod 150 &=& 45 &\mod 150 &=& 45 \\ x_3 &=& 3\;.\;45 &\mod 150 &=& 135 &\mod 150 &=& 135 \\ x_4 &=& 3\;.\;135 &\mod 150 &=& 405 &\mod 150 &=& 105 \\ x_5 &=& 3\;.\;105 &\mod 150 &=& 315 &\mod 150 &=& 15 \\ x_6 &=& 3\;.\;15 &\mod 150 &=& 45 &\mod 150 &=& 45 \\ x_7 &=& 135 \\ x_8 &=& 105 \\ x_9 &=& 15 \\ x_{10} &=& 45 \end{array} \end{align*}\] Sequência de valores pseudoaleatórios gerada: \[ \dfrac{15}{150} = 0,1 \qquad \qquad \dfrac{45}{150} = 0,3 \qquad \qquad \dfrac{135}{150} = 0,9 \qquad \qquad \dfrac{105}{150} = 0,7 \qquad \qquad \cdots \]

Implementação

gcmulti = function(n,a,m,seed){

  x = numeric(n+1)
  x[1] = seed

  for(i in 2:(n+1))
    x[i] = (a*x[i-1]) %% m

  u = x/m
  u = u[-1]
}

Voltando ao Exemplo 1

u = gcmulti(10,3,150,5)
u

##  [1] 0.1 0.3 0.9 0.7 0.1 0.3 0.9 0.7 0.1 0.3

Note que o número de valores distintos gerados é pequeno. Este número decorre da escolha dos valores de \(x_0\), \(a\) e \(m\).

Observação

De acordo com Ross(2022), a escolha de \(a\) e \(m\) deve satisfazer três critérios:

Para qualquer escolha de semente \(x_0\), a sequência resultante tem a ``aparência” de ser uma sequência de valores independentes de ua v.a. com distribuição uniforme em \((0,1)\).
Para qualquer escolha de semente \(x_0\), o número de valores gerados antes das repetições começarem deve ser grande.
Os valores podem ser calculados de forma eficiente pelo computador.

Uma recomendação que parece satisfazer estes critérios é utilizar

\(m \quad \longrightarrow \qquad\) maior número primo suportado pelo computador.

Exemplo 2

para computadores com tamanho de palavra de 32 bits: \[m = 2^{31}-1 = 2.147.483.647 \qquad \qquad a = 7^5 = 16.807\] Este caso particular é conhecido como Gerador de Park-Miller.

u = gcmulti(10000,16807,2^31-1,12345)

hist(u)

para computadores com tamanho de palavra de 36 bits: \[m = 2^{35}-31 \qquad \qquad a = 5^5 = 3125\]

u = gcmulti(10000,3125,2^35-1,12345)

hist(u)

Geradores Congruenciais Lineares Mistos

Definição

Um Gerador Congruencial Misto calcula valores sucessivos \(x_n\), \(n \geq 1\), de forma recursiva

\[ x_{n+1} = (a \; x_n+c) \mod m, \] a partir de um valor inicial \(x_0\), fixando \(a\), \(c\) e \(m\) inteiros positivos. Os valores fixados são chamados de - \(x_0\): semente; - \(a\): multiplicador; - \(c\): incremento; e - \(m\): módulo.

Novamente de acordo com a definição:

\(x_n \in \{0,1,2,...,m-1\}\); e
os números gerados

\[ u_n = \frac{x_n}{m}, \] chamados de números pseudoaleatórios, são tomados como um valor aleatório gerado de uma distribuição uniforme no intervalo \((0,1)\).

Exemplo 3

Encontrar a sequência de números pseudoaleatórios \(x_1,x_2,...,x_{10}\) supondo \[x_0=3 \qquad \mbox{ e } \qquad x_n = (5x_{n-1}+7) \mod 200.\]

Solução: \[\begin{align*} \begin{array}{cccccccc} x_0 &=& 3 \\ x_1 &=& (5\;.\;3+7) &\mod 200 &=& 22 &\mod 200 &=& 22 \\ x_2 &=& 5\;.\;22+7) &\mod 200 &=& 117 &\mod 200 &=& 117 \\ x_3 &=& (5\;.\;117+7) &\mod 200 &=& 592 &\mod 200 &=& 192 \\ x_4 &=& (5\;.\;192+7) &\mod 200 &=& 967 &\mod 200 &=& 167 \\ x_5 &=& (5\;.\;167+7) &\mod 200 &=& 842 &\mod 200 &=& 42 \\ x_6 &=& (5\;.\;42+7) &\mod 200 &=& 217 &\mod 200 &=& 17 \\ x_7 &=& (5\;.\;17+7) &\mod 200 &=& 92 &\mod 200 &=& 92 \\ x_8 &=& (5\;.\;92+7) &\mod 200 &=& 467 &\mod 200 &=& 67 \\ x_9 &=& (5\;.\;67+7) &\mod 200 &=& 342 &\mod 200 &=& 142 \\ x_{10} &=& (5\;.\;142+7) &\mod 200 &=& 717 &\mod 200 &=& 117 \\ \end{array} \end{align*}\] Sequência de valores pseudoaleatórios gerada: \[ \dfrac{22}{200} = 0,110 \qquad \qquad \dfrac{117}{200} = 0,585 \qquad \qquad \dfrac{192}{200} = 0,960 \qquad \qquad \dfrac{167}{200} = 0,835 \qquad \qquad \dfrac{42}{200} = 0,210 \] \[ \dfrac{17}{200} = 0,085 \qquad \qquad \dfrac{92}{200} = 0,460 \qquad \qquad \dfrac{67}{200} = 0,335 \qquad \qquad \dfrac{142}{200} = 0,710 \qquad \qquad \dfrac{117}{200} = 0,585 \]

Implementação

gcmisto = function(n,a,c,m,seed){

  x = numeric(n+1)
  x[1] = seed

  for(i in 2:(n+1))
    x[i] = (a*x[i-1]+c) %% m

  u = x/m
  u = u[-1]
}

Voltando ao Exemplo 3

u = gcmisto(10,5,7,200,3)
u

##  [1] 0.110 0.585 0.960 0.835 0.210 0.085 0.460 0.335 0.710 0.585

Novamente o número de valores distintos gerados decorre da escolha dos valores de \(x_0\), \(a\), \(c\) e \(m\).

A recomentadação é escolher \(m\) igual ao comprimento de palavra suportada pelo computador.

Exemplo 4

u = gcmisto(10000,16807,2,2^31-1,12345)

hist(u)

Exemplo 5: Relação entre os geradores

u1 = gcmulti(20,3,150,7)
u2 = gcmisto(20,3,0,150,7)

u1

##  [1] 0.14 0.42 0.26 0.78 0.34 0.02 0.06 0.18 0.54 0.62 0.86 0.58 0.74 0.22 0.66
## [16] 0.98 0.94 0.82 0.46 0.38

u2

##  [1] 0.14 0.42 0.26 0.78 0.34 0.02 0.06 0.18 0.54 0.62 0.86 0.58 0.74 0.22 0.66
## [16] 0.98 0.94 0.82 0.46 0.38

plot(u1,u2)

Propriedade do Período Máximo

O Gerador Congruencial Misto possui período máximo \(m\) se:

\(c\) e \(m\) são primos entre si;
\(a-1\) é divisível pelos fatores primos de \(m\);
se \(m\) é múltiplo de 4, então \(a-1\) também é.

Demonstração

O resultado decorre do Teorema de Hull-Dobell.

Exemplo 6: Verificando o período

Suponha um gerador congruencial misto com semente \(x_0=7\), multiplicador \(a=5\), incremento \(c=3\) e módulo \(m=16\).

Note que

\(c=3\) e \(m=16\) são primos entre si;
\(a-1=4\) é divisível por 2 (único fator primo de \(m=2^4\));
\(m=16\) é múltiplo de 4 e \(a-1=4\) também é.

Neste caso, o período máximo de \(m=16\) é atingido.

u = gcmisto(50,5,3,16,7)
u

##  [1] 0.3750 0.0625 0.5000 0.6875 0.6250 0.3125 0.7500 0.9375 0.8750 0.5625
## [11] 0.0000 0.1875 0.1250 0.8125 0.2500 0.4375 0.3750 0.0625 0.5000 0.6875
## [21] 0.6250 0.3125 0.7500 0.9375 0.8750 0.5625 0.0000 0.1875 0.1250 0.8125
## [31] 0.2500 0.4375 0.3750 0.0625 0.5000 0.6875 0.6250 0.3125 0.7500 0.9375
## [41] 0.8750 0.5625 0.0000 0.1875 0.1250 0.8125 0.2500 0.4375 0.3750 0.0625

Teste Visual de Independência

Suponha a sequência de v.a. \[U_1,U_2,...,U_n.\] Se esta sequência for composta por observações independentes e identicamente distribuídas no intervalo \((0,1)\), então o vetor bidimensional \[(U_i,U_{i+1})\] deve ter distribuição uniforme sobre o quadrado unitário \([0,1]×[0,1]\).

Neste caso, um gráfico de pontos \[(u_1,u_2),\;(u_2,u_3),\;...,\;(u_{n-1},u_n),\] conhecido como gráfico de dispersão de pares sucessivos (serial scatter plot ou lag plot), constitui uma das ferramentas gráficas mais simples para avaliar a independência serial de uma sequência pseudoaleatória.

Caso haja independência, os pontos do gráfico devem preencher o quadrado unitário de maneira aproximadamente homogênea, sem qualquer padrão visual. Idealmente, se o gerador produz valores aproximadamente independentes, espera-se observar uma “nuvem” uniforme de pontos, não existindo padrões como

alinhamentos;
faixas;
curvas;
regiões vazias;
aglomerados;

etc.

u = gcmisto(5000,1103515245,12345,2^31,123)

# Lembre-se que se U ~ U(0,1), então E[U] = 1/2 = 0,5 e Var(U) = 1/12 = 0,08333. 

mean(u)

## [1] 0.4951754

var(u)

## [1] 0.08556708

plot(u[-5000],u[-1],
     pch=19,
     cex=.2)

# Também esperamos correlações próximas de zero.
cor(u[-1],u[-length(u)])

## [1] -0.01912896

Exercício 1

Implemente um Gerador Congruencial Misto com:

\(a=7\)
\(c=5\)
\(m=32\)

e determine o período.

Compare histogramas produzidos por diferentes escolhas de parâmetros.
Verifique numericamente:

\[ E(U) = 0,5. \]

Faça um gráfico de dispersão de pares sucessivos e conclua sobre a independência dos valores gerados.
Pesquise os parâmetros utilizados pelo R para geração de números pseudoaleatórios.

Teste Qui-quadrado para verificar a Uniformidade dos números pseudoaleatórios

Gostaríamos de testar a hipótese \[ H_0: U \sim U(0,1), \] contra a hipótese alternativa \[ H_1: U \not\sim U(0,1). \]

A ideia do teste é comparar a distribuição observada dos valores gerados com a distribuição esperada caso a hipótese de uniformidade seja verdadeira.

Para isso, vamos realizar um teste de aderência qui-quadrado para verificar se os números pseudoaleatórios gerados possuem distribuição uniforme no intervalo \((0,1)\).

Para realizar o teste, o primeiro passo é dividir os dados em classes. A divisão em classes permite transformar os valores contínuos gerados pelo algoritmo em frequências observadas, que podem ser comparadas com as frequências esperadas sob a hipótese nula, através da estatística:

\[ \chi^2= \sum_{i=1}^{k} \frac{(O_i-E_i)^2}{E_i}, \] onde:

\(O_i\) representa a frequência observada na classe \(i\);
\(E_i\) representa a frequência esperada;
\(k=10\) é o número de classes.

Quando os números gerados seguem uma distribuição uniforme, espera-se que as diferenças entre valores observados e esperados sejam pequenas, produzindo um valor pequeno para \(\chi^2\).

Exemplo 7

Geração dos número pseudo-aleatórios:

n = 10000
u = gcmisto(n,16807,2,2^31-1,12345)

Divisão dos valores em classes e obtenção dos valores observados:

# cria 11 pontos igualmente espaçados entre 0 e 1, que serão utilizados para determinar 0 classes:
breaks = seq(0,1,length=11)

# faz as contagens em cada classe gerada
obs = table(cut(u,breaks=breaks))
obs

## 
##   (0,0.1] (0.1,0.2] (0.2,0.3] (0.3,0.4] (0.4,0.5] (0.5,0.6] (0.6,0.7] (0.7,0.8] 
##      1011       950      1036       989       969       968       967       998 
## (0.8,0.9]   (0.9,1] 
##      1058      1054

Obtenção da estatística de teste:

De fato, sob \(H_0\), todos os intervalos possuem a mesma probabilidade \[ P(U\in[0\;,\;0,1))= P(U\in[0,1\;,\;0,2))= \cdots= 0,1. \] Logo, a a frequência esperada em cada classes seria \[ E_i=n\times0,1. \]

Assim, esperamos aproximadamente \(10.000*0,1 = 1000\) observações em cada intervalo.

chisq.test(obs)

## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 13.396, df = 9, p-value = 0.1455

Interpretação do resultado

O resultado produzido pelo R apresenta o valor da estatística qui-quadrado, os graus de liberdade e o p-valor.

Para um dado nível de significância \(\alpha\), se

\[ p\text{-valor}>\alpha, \]

não rejeitamos a hipótese nula; ou seja não existe evidência estatística suficiente para afirmar que os valores gerados deixam de seguir uma distribuição uniforme em \((0,1)\).

\[ p\text{-valor}\leq \alpha, \]

rejeitamos a hipótese nula; portanto existem evidências estatísticas de que os valores gerados não seguem uma distribuição uniforme em \((0,1)\).

Aplicação: Simulação de uma variável Bernoulli

u = gcmulti(10000,16807,2^31-1,100)

x = ifelse(u < 0.3,1,0)

mean(x)

## [1] 0.3014

Resumindo

Nenhum gerador é verdadeiramente aleatório.

Deseja-se:

grande período;
independência aproximada;
boa distribuição.

Exercício 2

Implemente o gerador de Park-Miller.
Compare os histogramas de um GCL e de um GCM.
Calcule autocorrelações de ordens:

Aplique o teste qui-quadrado para diferentes tamanhos amostrais.
Pesquise outros testes de aleatoriedade disponíveis.

Referências

ROSS, Sheldon M. Simulation. 6. ed. London: Academic Press, 2022.

ROBERT, Christian P.; CASELLA, George. Monte Carlo Statistical Methods. 2. ed. New York: Springer, 2004.

Geração de Números Pseudoaleatórios

Disciplina: GET00204 - Estatística Computacional 2

Motivação

Geradores Congruenciais Multiplicativos

Definição

Exemplo 1

Implementação

Voltando ao Exemplo 1

Observação

Exemplo 2

Geradores Congruenciais Lineares Mistos

Definição

Exemplo 3

Implementação

Voltando ao Exemplo 3

Exemplo 4

Exemplo 5: Relação entre os geradores

Propriedade do Período Máximo

Demonstração

Exemplo 6: Verificando o período

Teste Visual de Independência

Exercício 1

Teste Qui-quadrado para verificar a Uniformidade dos números pseudoaleatórios

Exemplo 7

Aplicação: Simulação de uma variável Bernoulli

Resumindo

Exercício 2

Referências