2  Variáveis Aleatórias Uniformes

Neste capítulo vamos considerar a geração de números aleatórios para o modelo probabilístico uniforme. A partir do modelo uniforme podemos gerar realizações de variáveis aleatórias de qualquer outro modelo probabilístico. Para gerar realizações de uma distribuição uniforme, precisamos gerar números aleatórios. Isso não pode ser realizado por máquinas. Na verdade qualquer sequência produzida por uma máquina é na verdade uma sequência previsível. Dessa forma, a denominamos de sequência de números pseudo-aleatórios.

Uma sequência de números será considerada “aleatória” do ponto de vista computacional se o programa que a gerar for diferente e estatisticamente não correlacionado com o programa que a usará. Assim, dois geradores de números aleatórios deveriam produzir os mesmos resultados nas suas aplicações. Se isso não ocorrer, um deles não pode ser considerado um bom gerador de números aleatórios (Press et al. 1992).

Os conceitos de números uniformes e números aleatórios podem ser muitas vezes confundidos. Números uniformes são aqueles que variam aleatoriamente em uma faixa determinada de valores com probabilidade constante. No entanto, devemos diferenciar números aleatórios uniformes de outros tipos de números aleatórios, como por exemplo, números aleatórios normais ou gaussianos. Estes outros tipos são geralmente provenientes de transformações realizadas nos números aleatórios uniformes. Então, uma fonte confiável para gerar números aleatórios uniformes determina o sucesso de métodos estocásticos de inferência e de todo processo de simulação Monte Carlo.

2.1 Números Aleatórios Uniformes

Números uniformes aleatórios são aqueles que, a princípio, se situam dentro de um intervalo real, geralmente, entre \(0\) e \(1\), para os quais não podemos produzir uma sequência previsível de valores e cuja função densidade é constante. Em vários programas de computadores estes números são gerados utilizando o comando random ou comandos similares. Em Pascal, por exemplo, se este comando for utilizado com o argumento \(n\), random(n), números aleatórios inteiros \(U\) do intervalo \(0 \leq U \le n-1\) são gerados e se o argumento \(n\) não for usado, os números gerados são valores aleatórios reais do intervalo \([0, 1)\).

Em geral, os programas utilizam o método congruencial. Vamos considerar os números uniformes inteiros \(U_1, U_2, U_3, \dots\) entre \(0\) e \(m-1\), em que \(m\) representa um grande número inteiro. Podemos gerar estes números utilizando o método congruencial por meio da relação recursiva:

\[ U_{i+1}=\left( a U_{i} + c \right) \mod m \qquad(2.1)\]

em que \(m\) é chamado de módulo, \(a\) e \(c\) são inteiros positivos denominados de multiplicador e incremento, respectivamente. O operador mod retorna o resto da divisão do argumento \(\left( a U_{i} + c \right)\) por \(m\). A sequência recorrente (Equation 2.1) se repete em um período que não é maior que \(m\), por razões óbvias. Se \(a\), \(c\) e \(m\) são adequadamente escolhidos, a sequência tem tamanho máximo igual a \(m\). A escolha do valor inicial \(U_0\) é também muito importante. O valor do número uniforme correspondente no intervalo de \(0\) a \(1\) é dado por \(U_{i+1}/m\), que é sempre menor que \(1\), mas podendo ser igual a zero.

Vamos apresentar um exemplo didático para ilustrar um gerador de números aleatórios. Sejam \(U_0=a=c=7\) e \(m=10\), logo,

\[U_1=(7\times 7+7) \mod 10 = 56 \mod 10 = 6\] \[U_2=(7\times6+7) \mod 10 = 49 \mod 10 = 9\] e assim sucessivamente. Obtemos a sequência de números aleatórios: \[ \{ 7, 6, 9, 0, 7, 6, 9, 0, 7, 6, 9, \cdots \} \] e verificamos que o período é igual a \(4\), \(\{7,6, 9, 0, \cdots\}\), que é menor do que \(m=10\).

Este método tem a desvantagem de ser correlacionado serialmente. Se \(m\), \(a\) ou \(c\) não forem cuidadosamente escolhidos a correlação pode comprometer a sequência gerada. Por outro lado, o método tem a vantagem de ser muito rápido. Podemos perceber que a cada chamada do método, somente alguns poucos cálculos são executados. Escolhemos, em geral, o valor de \(m\) pelo maior inteiro que pode ser representado pela máquina de \(32\) bits, qual seja, \(2^{32}\). Um exemplo que foi utilizado por muitos anos nos computadores IBM mainframe, que representam uma péssima escolha é \(a= 65.539\) e \(m = 2^{31}\).

A correlação serial não é o único problema desse método. Os bits de maior ordem são mais aleatórios do que os bits de menor ordem (mais significantes). Devemos gerar inteiros entre 1 e 20 por \(j=1+int (20 \times random(semente))\), ao invés de usar o método menos acurado \(j=1+\mod (int(1000000\times random(semente)),20)\), que usa bits de menor ordem. Existem fortes evidências, empíricas e teóricas, que o método congruencial

\[ U_{i+1}= a U_{i} \mod m \qquad(2.2)\]

é tão bom quanto o método congruencial com \(c \ne 0\), se o módulo \(m\) e o multiplicador \(a\) forem escolhidos com cuidado (Press et al. 1992). Park and Miller (1988) propuseram um gerador “padrão” mínimo baseado nas escolhas: \[ a = 7^5 = 16.807 \quad m = 2^{31} - 1 = 2.147.483.647 \qquad(2.3)\]

Este gerador de números aleatórios não é perfeito, mas passou por todos os testes a qual foi submetido e tem sido usado como padrão para comparar e julgar outros geradores. Um problema que surge e que devemos contornar é que não é possível implementarmos diretamente em uma linguagem de alto-nível a equação (Equation 2.2) com as constantes de (Equation 2.3), pois o produto de \(a\) e \(U_i\) excede, em geral, o limite máximo de \(32\) bits para inteiros. Podemos usar um truque, devido a Schrage (1979), para multiplicar inteiros de \(32\) bits e aplicar o operador de módulo, garantindo portabilidade para implementação em praticamente todas as linguagens e todas as máquinas. O algoritmo de Schrage baseia-se na fatoração de \(m\) dada por:

\[ m= aq + r;\qquad \textrm{i.e.}, \qquad q = \lfloor m/a \rfloor; \quad r =m \mod a \]

em que \(\lfloor z \rfloor\) denota a parte inteira do número \(z\) utilizado como argumento. Para um número \(U_{i}\) entre 1 e \(m-1\) e para r pequeno, especificamente para \(r<q\), Schrage (1979) mostrou que ambos \(a(U_i \mod q)\) e \(r\lfloor U_i/q \rfloor\) pertencem ao intervalo \(0 \cdots m-1\) e que \[ aU_{i} \mod m=\left\{% \begin{array}{ll} a(U_i \mod q) - r\lfloor U_i/q \rfloor & \textrm{se maior que 0}\\ a(U_i \mod q) - r\lfloor U_i/q\rfloor + m & \textrm{caso contrário.}\\ \end{array} \right.% \qquad(2.4)\]

Computacionalmente observamos que a relação: \[ a(U_i \mod q) = a(U_i-q\lfloor U_i/q\rfloor)\] se verifica. No Python pode-se optar por usar o operador % (mod), que retorna o resto da operação entre dois inteiros e o operador // (div), que retorna o resultado do dividendo para operações com inteiros. A quantidade \(U_i \mod q = U_i - (U_i \quad \textrm{div} \quad q)\times q\) pode ser obtida em Python simplesmente com \(U_i \% q\). Atribuímos o resultado a uma variável qualquer definida como inteiro. Para aplicarmos o algoritmo de Schrage às constantes de (Equation 2.3) devemos usar os seguintes valores: \(q = 127.773\) e \(r = 2.836\).

A seguir apresentamos o algoritmo do gerador padrão mínimo de números aleatórios: 

# gerador padrão mínimo de números aleatórios adaptado de Park and
# Miller. Retorna desvios aleatórios uniformes entre 0 e 1. Fazer
# "sem" igual a qualquer valor inteiro para iniciar a sequência;
# "sem" não pode ser alterado entre sucessivas chamadas da sequência
# se "sem" for zero ou negativo, um valor dependente do valor do relógio
# do sistema no momento da chamada é usado como semente. Constantes
# usadas a = 7^5 = 16.807; m = 2^31 - 1 = 2.147.483.647
# e c = 0

def  gnup(sem, q, a, r, m):
  k = sem // q # divisão por inteiros
  sem = a * (sem % q) - r * k
  if sem < 0:
    sem = sem + m
  u =  sem / m
  res = {'sem': sem, 'u': u}
  return res

from datetime import datetime # obter o data/horário do sistema

def gnap(n, sem = 0):
  a = 16807; m = 2147483647
  q = 127773; r = 2836
  mr = 1 / m
  if sem <= 0:
    t = datetime.now()
    sem = t.second+t.minute*60+t.hour*3600+t.day*86400
  u = []
  for i in range(n):
    x = gnup(sem, q, a, r, m)
    u.append(x['u'])
    sem = x['sem']
  return u  
# Exemplos de uso
n = 5 
x = gnap(n,0)
# Formatando a saída para 5 casas decimais
print(["%0.5f" % v for v in x])
#  especificando a semente
y = gnap(n, 1001)
# Formatando a saída para 5 casas decimais
print(["%0.5f" % v for v in y])
# tempo de execução para cada número aleatório
n = 100000
t1 = datetime.now()
x = gnap(n)
t2 = datetime.now()
t = t2-t1
print('tempo médio (micros): ',t.microseconds / n)
print('tempo total (micros): ',t.microseconds)
['0.80348', '0.04011', '0.17663', '0.58029', '0.89324']
['0.00783', '0.66933', '0.36093', '0.10878', '0.30000']
tempo médio (micros):  0.28273
tempo total (micros):  28273

Foi computado o tempo médio e o tempo total para rodar 100000 números aleatórios uniformes. Para capturar o tempo, foi usado o pacote datetime e também para obter a diferença de tempo entre o início e o fim do processamento de nossa função. O bloco para fazer isso deve ser executado de uma só vez, ou seja, marcando o bloco e teclando enter no Positron.

Algumas considerações a respeito desse algoritmo: a) a função gnap (gerador de números aleatórios mínima) retorna um número real entre \(0\) e \(1\). Este tipo de especificação determina que as variáveis possuam precisão dupla. A precisão dupla (double-precision float) possui números na faixa de \(\pm 1,7,0 \times 10^{-308}\) \(\cdot\) \(\pm 1,7 \times 10^{308}\), ocupa \(24\) bytes de memória e possui \(15-16\) dígitos significantes; b) o valor da semente é definido pelo usuário e é passado como parâmetro para a função. Isso significa que a variável do programa que chama a função e que é passada como semente deve ser atualizada com o novo valor modificado em gnup. Se o seu valor inicial for zero ou negativo, a função atribui um inteiro dependente da hora do sistema no momento da chamada; c) a função tem dependência do pacote datetime, que foi usado para capturar a hora e dia do sistema.

A rotina é iniciada com os valores de \(n\) e da semente fornecidos pelo usuário. Se a semente for nula ou negativa, atribuímos um novo valor dependente do relógio do sistema no momento da chamada. A partir deste ponto o programa deve chamar reiteradas vezes a função gnap, que retorna o valor do número aleatório \(0\) e \(1\) utilizando o algoritmo descrito anteriormente, até que a sequência requerida pelo usuário seja completada. Nas sucessivas chamadas desta função, o valor da semente é sempre igual ao valor do último passo.

O período de gnap é da ordem de \(2^{31}\approx2,15\times 10^{9}\), ou seja, a sequência completa é um pouco superior a 2 bilhões de números aleatórios. Assim, podemos utilizar gnap para alguns poucos propósitos práticos. Como já salientamos o gerador padrão mínimo de números aleatórios possui duas limitações básicas, quais sejam, sequência curta e correlação serial. Assim, como existem métodos para eliminar a correlação serial e que aumentam o período da sequência, recomendamos que sejam adotados. Claro que a função apresentada teve por objetivo ilustrar como podemos programar nossas próprias funções para gerarmos números aleatórios uniformes. O Python, no entanto, possui seu próprio gerador de números uniformes, que veremos na sequência. Um dos melhores e mais interessantes geradores de números aleatórios é o Mersenne Twister (MT). Mersenne Twister é um gerador de números pseudo-aleatórios desenvolvido por Makoto Matsumoto e Takuji Nishimura nos anos de 1996 e 1997 (Matsumoto and Nishimura 1998). O MT possui os seguintes méritos segundo seus desenvolvedores:

  • foi desenvolvido para eliminar as falhas dos diferentes geradores existentes;

  • possui a vantagem de apresentar o maior período e maior ordem de equidistribuição do que de qualquer outro método implementado. Ele fornece um período que é da ordem de \(2^{19.937}-1 \approx 4,3154 \times 10^{6001}\), e uma equidistribuição 623-dimensional;

  • é um dos mais rápido geradores existentes, embora complexo;

  • faz uso de forma muito eficiente da memória.

Existem muitas versões implementadas deste algoritmo, inclusive em Fortran e C e que estão disponíveis na internet. Felizmente, o Python já possui este algoritmo implementado. Por se tratar de um tópico mais avançado, que vai além do que pretendemos apresentar nestas notas de aulas, não descreveremos este tipo de procedimento para incorporações de funções escritas em outras linguagens.

2.2 Números Aleatórios Uniformes no Python

No Python podemos gerar números aleatórios uniformes contínuos utilizando uma função pré-programada. Os números aleatórios uniformes são gerados pelo comando random.uniform(low=0.0, high=1.0, size=None) da biblioteca numpy, em que None é para gerar apenas um valor, podendo ser substituído por \(n\), entre low e high (excluído), que são argumentos que delimitam o valor mínimo e máximo da sequência a ser gerada. O controle da semente para se gerar uma sequência reproduzível de números uniformes é dada pelo comando random.seed(seed=None) do numpy, em que o argumentoseed deve ser um número inteiro. O Python automaticamente determina a cada chamada uma nova semente. Conseguimos gerar diferentes sequências em cada chamada do comando gerador, sem nos preocuparmos com a semente aleatória. O gerador de números aleatórios uniformes usa o algoritmo Mersenne-Twister por padrão.

No programa apresentado a seguir ilustramos como podemos gerar \(n\) números aleatórios uniformes entre \(0\) e \(1\) de forma compacta, simples e eficiente:

import numpy as np
n = 5
x = np.random.uniform(low=0.0, high=1.0, size=n)
print(x)
# Fixando a semente
np.random.seed(seed=1000)
np.random.uniform(low=0.0, high=1.0, size=n)
np.random.seed(seed=1000)
np.random.uniform(low=0.0, high=1.0, size=n) # igual a anterior (mesma semente)
[0.34992956 0.09726348 0.17952792 0.2662993  0.91158536]
array([0.65358959, 0.11500694, 0.95028286, 0.4821914 , 0.87247454])
array([0.65358959, 0.11500694, 0.95028286, 0.4821914 , 0.87247454])

Fizemos um programa para comparar o tempo de execução das funções gnap e random.uniform() e retornamos o tempo médio para cada realização da variável aleatória. Desta forma verificamos que o algoritmo random.uniform é mais rápido de todos, conforme valores relativos apresentados a seguir. Obviamente temos que considerar que o algoritmo do numpy é uma função compilada. O algoritmognap por sua vez foram implementados em Python, que usa uma linguagem interpretada. As comparações de tempo podem ser vistas no programa a seguir. O programa utilizado foi:

# tempo de execução para cada número aleatório
# comparativo entre nosso gerador e o do np
n = 100000
t1 = datetime.now()
x = gnap(n)
t2 = datetime.now()
tgp = t2-t1
print('tempo médio gnap: ',tgp.microseconds / n)
t1 = datetime.now()
x = np.random.uniform(low=0.0, high=1.0, size=n)
t2 = datetime.now()
tnp = t2-t1
print('tempo médio numpy: ',tnp.microseconds / n)
# tempo relativo
tgp / tnp
tempo médio gnap:  0.26367
tempo médio numpy:  0.01228
21.471498371335503

Podemos fazer um histograma usando o programa a seguir. Para isso usamos a biblioteca matplotlib e o resultado foi:

import matplotlib.pyplot  as plt
n = 100000
x = np.random.uniform(low=0.0, high=1.0, size=n) # gnap(n) troque para testar
graf = plt.hist(x, bins='auto', color='#14e8f3',rwidth=0.95,alpha=0.7)

2.3 Exercícios

  1. Utilizar o gerador gnap para gerar \(n\) realizações de uma distribuição exponencial \(f(x)= \lambda e^{-\lambda x}\). Sabemos do teorema da transformação de probabilidades, que se \(U\) tem distribuição uniforme, \(X=F^{-1}(U)\) tem distribuição de probabilidade com densidade \(f(x)=F'(x)\); em que \(F(x)=\int_{-\infty}^{x}f(t)dt\) é a função de distribuição de \(X\) e \(F^{-1}(y)\) é a sua função inversa para o valor \(y\). Para a exponencial a função de distribuição de probabilidade é: \(F(x)=\int_{0}^{x}\lambda e^{-\lambda t}dt=1-e^{-\lambda x}\). Para obtermos a função inversa temos que igualar \(u\) a \(F(x)\) e resolver para x. Assim, \(u=1-e^{-\lambda x}\) e resolvendo para \(x\) temos: \(x=-\ln{(1-u)}/\lambda\). Devido à simetria da distribuição uniforme \(1-u\) pode ser trocado por \(u\). O resultado final é: \(x=-\ln{(u)}/\lambda\). Para gerar números da exponencial basta gerar números uniformes e aplicar a relação \(x=-\ln{(u)}/\lambda\). Fazer isso para construir uma função que gera \(n\) realizações exponenciais. Aplicar a função para obter amostras aleatórias da exponencial de tamanho \(n=100\) e obter o histograma da amostra simulada. Calcule a média e a variância e confronte com os valores teóricos da distribuição exponencial.

  2. Para gerar números de uma distribuição normal, cuja densidade é dada por \(f(x) = 1/(\sqrt{2\pi \sigma^2})\exp\{-(x-\mu)^2/(2\sigma^2)\}\), qual seria a dificuldade para podermos utilizar o teorema anunciado no exercício proposto anterior?

  3. Como poderíamos adaptar o algoritmo apresentados nesse capítulo para gerar números aleatórios uniformes utilizando os valores propostos por Park e Miller, ou seja, \(a=48.271\) e \(m=2^{31}-1\)? Implementar o algoritmo, tomando cuidado em relação aos novos multiplicador \(q\) e resto \(r\) da fatoração de \(m\)?

  4. Como você poderia propor um teste estatístico simples para avaliar a aleatoriedade da sequência de números uniformes gerados por esses algoritmos apresentados no capítulo? Implementar sua ideia.