Gerando Números Aleatórios Fernando Henrique Ferraz Pereira da Rosa Vagner Aparecido Pedro Junior MAP-131 Laboratório de Matemática Aplicada Prof. Dr. Eduardo Colli 21 de novembro de 2002 1 Introdução A disponibilidade de números aleatórios é extremamente útil em uma variedade de situações, como na simulações de fenômenos fı́sicos (indo desde a fı́sica nuclear até engenharia de sistemas), em amostragem de uma população, na programação de computadores, na tomada de decisões ou até mesmo em entretenimento (bingos, loterias ou jogos). Na área de simulação, onsideremos por exemplo, a modelagem do tempo de acesso de um disco rı́gido, num computador pessoal. Podemos determinar que a duração desse evento irá ficar numa faixa conhecida, digamos de 0 a 200ms, de acordo com caracteristicas fı́sicas inerentes ao próprio disco rigido. Entretanto, o valor real desse evento vai depender de vários fatores, como a posição da cabeça de leitura quando a requisição é feita pelo sistema operacional, detalhes espúrios da implementação do suporte à controladora no sistema e até mesmo da temperatura e outras condições ambientais. Podemos considerar então que esse tempo de acesso é uma variável aleatória seguindo uma distribuição conveniente, e simulamos com o objetivo de ver se nosso modelo adere a realidade. Para fazermos essa simulação precisamos de números aleatórios que sigam essa dada distribuição, e para isso precisamos saber primeiro como gerar esses números aleatórios. Algumas fontes de números aleatórios são o lançamento de dados, a retirada de bolas numeradas de uma urna (com reposição), o uso de uma roleta ou ainda ruı́do eletrônico cuja saida é quantizada periodicamente [1]. Entretanto na esmagadora maioria das vezes usa-se o que foi convecionado chamar de números pseudo-aleatórios. Isso se dá pois o uso de números realmente aleatórios causa um impecı́lio nas simulações assim como no desenvolvimento de programas para computador: não se tem como repetir uma dada seqüência de números, com a intenção de verificar a simulação ou tentar corrigir um problema no programa. Além da possibilidade se repetir uma dada seqüência, de [1] temos que existem algumas outras propriedades importantes num bom gerador de números (pseudo)aleatórios, a saber: 1. Os números gerados devem seguir uma distribuição uniforme, pois números aleatórios de verdade seguem essa distribuição. 2. Os números devem ser estatisticamente independentes entre si. O valor de um número na seqüência não deve afetar o valor do próximo (na prática a maioria dos geradores usa seqüências recursivas, então há essa dependência dos valores anteriores, mas isso não é estatisticamente significafivo, dai o destaque para independência estatı́stica). 3. A seqüência não deve se repetir nunca. Isso é teoricamente impossı́vel mas na prática um perı́odo de repetição suficientemente grande é o bastante. 4. A geração desses números deve ser rápida, de modo a poupar recursos computacionais para as simulações em si. O uso de um algoritmo para gerar um número aleatório parece violar o princı́pio básico da aleatoridade, por isso é que se convenciona chamar esses números de sintéticos ou pseudo-aleatórios. Esses números correspondem às propriedades acima mas começam sempre de um valor inicial chamado semente 1 (seed ) e continuam numa maneira totalmente determinı́stica, por isso deve ser tomado muito cuidado para que a aleatoridade esteja presente. No desenvolvimento do trabalho vamos dar uma olhada nos principais algoritmos para geração de números aleatórios e em seguida vamos estudar também alguns testes de aleatoridade, que nos fornecem um meio objetivo de verificar a aleatoridade de uma seqüência. Vamos também implementar alguns desses métodos, e fazer testes para verificar como eles se comportam. 2 2 2.1 Algoritmos para geração de números pseudoaleatórios Método do quadrado do meio Esse método foi inventado por John von Neumann. Começa-se com uma semente (seed), esse número é então elevado ao quadrado, os dı́gitos do centro são usados como próximo elemento da seqüência. Caso o número de dı́gitos que fique a esquerda seja maior que os que fiquem a direita não há problema, simplesmente fixamos para qual lado vamos fazer o corte. Exemplo: começando a partir de x0 = 44214 vamos gerar uma seqüência de números aleatórios de 5 dı́gitos: = = = = = = x0 (44214)2 (54877)2 (11485)2 (19052)2 (29787)2 44214 1954877796 3011485129 0131905225 0362978704 0887265369 ⇒ ⇒ ⇒ ⇒ ⇒ = = = = = x1 x2 x3 x4 x5 48777 14851 19052 29787 72653 Veja 4.1 para a implementação em c do algoritmo do quadrado do meio (midsquare). Esse algoritmo tem diversas desvantagens. Em primeiro lugar geralmente as seqüências geradas por esse algoritmo se repetem em pouco tempo. Em segundo, sempre que um 0 for gerado todos os números seguintes da seqüência serão 0 também. Exemplo: tomando x0 = 121: x0 (121)2 (464)2 (529)2 (984)2 (825)2 (062)2 (384)2 (745)2 (502)2 (200)2 (000)2 2.2 = = = = = = = = = = = = 121 014641 215296 279841 968256 680625 003844 147456 555025 252004 040000 000000 ⇒ x1 ⇒ x2 ⇒ x3 ⇒ x4 ⇒ x5 ⇒ x6 ⇒ x7 ⇒ x8 ⇒ x9 ⇒ x10 ⇒ x11 = = = = = = = = = = = 464 529 984 825 062 384 745 502 200 000 000 Método linear congruente A grande maioria dos geradores de números pseudo-aleatórios hoje em dia são modificações do método linear congruente. Ele é dado através duma relação de recorrência, a saber: xn+1 = (axn + c) mod m, n ≥ 0 (1) O valor inicial x0 é chamado de semente (seed ), a é o multiplicador, c o incremento e m é o módulo. A escolha desses valores afeta drasticamente o 3 comprimento do perı́odo da seqüência de números aleatórios. É fácil notar que o perı́odo será sempre menor que m, e que portanto, para aumentarmos o perı́odo devemos usar m o maior possı́vel. Esse algoritmo apesar de extremamente simples é poderossı́mo, dada uma escolha adequada dos parâmetros a, c e m. Em [2] temos o seguinte Teorema que nos informa qual a escolha ótima desses parâmetros: Teorema A O método linear congruente tem um perı́odo m se e somente se: i) c e m são primos entre si. ii) b = a − 1 é um múltiplo de p, para todo p primo divisor de m iii) b é um múltiplo de 4, se m é um múltiplo de 4. Uma prova desse Teorema para um caso especial foi dada há pelo menos 100 anos atrás, por M. Greenberger e o caso geral foi dado por Hull e Dobel. Não iremos entrar em detalhes na sua demonstração matemática, que envolve diversos conceitos avançados de Álgebra. Em 4.2 implementamos em c esse algoritmo para geração de números aleatórios. Fizemos algumas modificações e escolhas particulares dos parâmetros a, m e c de modo que o usuário precisa entrar apenas com o número de números aleatórios que ele quer gerar (entre 0 e 1) e opcionalmente uma semente. 2.2.1 Escolhendo os parâmetros Em primeiro lugar, levamos em consideração o valor de m. Como já foi observado, o perı́odo da seqüência na melhor das hipóteses vai ser igual a m (ou seja, a cada m elementos a seqüência vai se repetir). Baseados no número de elementos n que o usuário quer gerar, escolhemos então um valor de m suficientemente maior que n mas não necessariamente sempre grande, de forma a poupar recursos computacionais, tendo m como função de n. Chegamos na seguinte relação: √ m = (menor primo( 10n))2 (2) Onde a função menor primo(x) retorna o maior número primo d, tal que d<x. Com (2) garantimos que m será suficientemente maior que n e ao mesmo tempo que m é da forma m = p2 com p primo, e logo os únicos divisores de m serão 1, m e p. Convenientemente então, escolhemos a de acordo com a seguinte relação: √ a = menor primo( 10n) + 1. (3) Seja: √ q = menor primo( 10n) Por (2) sabemos que q é o único divisor primo de m, e da condição (ii) do Teorema A sabemos que a−1 deve ser múltiplo de p para todo p primo divisor de m. Ora, temos então que a deve ser da forma: a − 1 = nq ⇒ a = nq + 1 (n ∈ N) Fazendo n = 1 chegamos na relação (3) acima. 4 m Escolhemos c aleatório entre 0 e 10 , c diferente de q. Dessa forma completamos a condição (i) do Teorema A. A condição (iii) não precisa ser levada em conta, já que escolhemos m múltiplo apenas de q, q primo. E por fim, caso o usuário não tenha escolhido ele próprio uma semente, escolhemos x0 aleatório entre 0 e m. Obtivemos resultados excelentes com essa implementação, conseguindo seqüências melhores que de alguns pacotes estatı́sticos disponı́veis. Mais a frente, quando entrarmos na parte de Testes de Aleatoridade vamos falar disso com mais detalhes. 2.3 Outros métodos Claramente, o método linear congruente não é o único meio de se conseguir boas seqüências de números aleatórios. Nessa sessão vamos dar uma olhada em alguns algoritmos outros algoritmos disponı́veis. Uma versão quadrática do método linear congruente é dada por: xn+1 = (dx2n + axn + c) mod m. (4) Com a, d e m escolhidos convenientemente, de forma similar que o método linear, também pode-se garantir que a seqüência vai ter perı́odo m. Outro método quadrático interessante proposto, foi: xn+1 = xn (xn + 1) mod 2e , n ≥ 0. (5) Com m uma potência de 2 e x0 mod 4 = 2, também fornece um perı́odo longo. Outras generalizações podem ser feitas, obtendo xn+1 a partir de um polinômio que depende de xn , xn−1 , ..., xk . Aplicam-se aı́ também regras às constantes e pode-se garantir perı́odos suficientemente longos. 5 3 Testes de Aleatoriedade Nós acabamos de estudar alguns métodos para geração de seqüências que parecem aleatórias. Vimos que seguindo alguns critérios nas escolhas dos algoritmos e parâmetros podemos aumentar o perı́odo da seqüência o quanto quisermos, mas somente isso não garante que nossa seqüência será útil para propósitos práticos. Considere a seqüência: 0, 1, 2, 3, 4, 5, 6, ..., n − 1, n O perı́odo pode ficar tão grande quanto quisermos, mas essa seqüência não é aleatória de forma alguma. Qual deve ser o critério de decisão para que possamos concluir com confiança se uma seqüência é aleatória ou não? Se for dada uma caneta e um papel para uma pessoa e pedir-se para essa pessoa escrever uma seqüência de números aleatórios, provavelmente a pessoa não conseguirá escrever uma seqüência satisfatoriamente aleatória. Nós tendemos a evitar coisas que não parecem aleatórias, como dois números repetidos, enquando que estatisticamente se espera que isso ocorra com freqüência. Ao mesmo tempo que se mostrarmos uma tabela de números aleatórios para algúem essa pessoa vai dizer que os números não são aleatórios. Um outro exemplo da nossa errônea intuição é com relação à expansão decimal do número π. Um estatı́stico diria que é apenas uma seqüência de números aleatórios, enquanto um numerólogo vai encontrar provavelmente toda a história da humanidade ali. Em [3] podemos procurar qualquer seqüência de dı́gitos em π. Pegando as 4 primeiras notas da primeira prova de MAP-131 e colocando em seqüência: 4.85.32.13.8 = 48532138 Obtemos o resultado: The string 48532138 was found at position 3135428 counting from the first digit after the decimal point. O que significa que o resultado da prova já estava predestinado por π ou... precisamos de um conceito mais rigoroso de aleatoridade. A seguir veremos os principais testes utilizados para se verificar a aleatoridade de uma seqüência. Em geral supõe-se que a seqüência não é aleatória, e se ela passar satisfatoriamente por uma meia dúzia de testes, consideramos ela satisfatoriamente aleatória. 3.1 Testes de Qui-Quadrado O teste de Qui-Quadrado é provavelmente o teste estatı́stico mais conhecido, e serve para diversas situações, como teste de homogeneidade de subpopulações, aderência de modelos ou independência de varı́aveis. Sua essência se baseia na soma dos resı́duos em relação a um modelo esperado, e pode ser esquematizado, de acordo com [4] da seguinte maneira: Dada uma variável X para a qual temos uma amostra de valores, desejamos verificar se essa varı́avel se adequa ou não a uma dada distribuição. Dividimos esses valores em k categorias, montando a tabela de freqüência: Categoria 1 Freq. Observada o1 6 2 o2 3 o3 .... ... k ok A partir do modelo que supomos ser adequado, construı́mos a tabela de freqüências esperadas de acordo com as categorias acima: Categoria 1 Freq. Esperada e1 2 e2 3 e3 .... ... k ek Calculamos então uma medida de distância entre os dados observados e o que se esperaria sob a hipótese de que ele segue o modelo que supomos, através da seguinte expressão: k X (oi − ei )2 Q2 = ei i=1 Sabe-se que para uma amostra razoavelmente grande, essa quantidade segue uma distribuição de Qui-Quadrado, com k − 1 graus de liberdade. Calculamos então: P (χ2k−1 > Q2obs ) (6) Caso esse valor seja muito pequeno (≤ 0.01) rejeitamos a hipótese de que nossos dados seguem o modelo sugerido. Para testar seqüências de números aleatórios, temos duas aplicações dos testes de Qui-Quadrado. Em primeiro lugar podemos fazer um teste de homogeneidade de subpopulações, dividindo os números aleatórios que geramos em seqüências de m valores e verificando se a distribuição dos números entre esses grupos é homogênea ou não. Outro teste possı́vel é pegar toda a seqüência de números aleatórios, dividı́la em k classes, e verificar as freqüências observadas. Supomos então que essa varı́avel segue distribuição uniforme contı́nua, calculamos as freqüências esperadas sob essa hipótese, e verificamos então através de (6) se nossa suposição está correta. 3.2 Testes usando Monte Carlo No trabalho anterior, usamos métodos Monte Carlo para estimar o valor de π. Um dos experimentos consistia em jogar pontos aleatórios numa circunferência circunscrita num quadrado, e fazendo uma razão incluindo as áreas das duas formas geométricas, estimavamos π. Quanto mais aleatoriamente conseguimos jogar os pontos na circunferência, ou em outras palavras, quanto mais aleatórios são os números que conseguimos gerar, melhor vai ser nossa estimativa de π. Para testar a aleatoridade de uma seqüência podemos então fazer o teste do quadrado e ver quão boa é a estimativa de π. 3.3 Teste de Kolmogorov-Smirnov Um outro teste disponı́vel é o teste de Kolmogorov-Smirnov, que pode ajudar a determinar se um dado conjunto de valores provêm de uma dada distribuição ou não. Geramos os números aleatórios, normalmente entre 0 e 1, e verificamos se esse conjunto de dados adere a uma distribuição uniforme no intervalo [0,1]. 3.4 Testes seriais Os testes que vimos até agora verificam a uniformidade de uma seqüencia e suas subseqüências, mas ainda assim, não são capazes de detectar o nı́vel de 7 dependência entre um número e seu antecessor, para isso temos alguns procedimentos especı́ficos. Um dos testes mais comuns desse gênero é o chamado teste do Poker. Neste teste considera-se n grupos de 5 inteiros sucessivos, e observamos em quais classisficações se encaixam as seqüências: 5 4 3 2 1 diferentes diferentes diferentes diferentes diferentes = = = = = todos os números diferentes entre si. um par. dois pares ou três iguais. quatro iguais. todos os 5 iguais. Fazemos então um teste de Qui-Quadrado, supondo que o número de casos vai ser equivalente para cada quı́ntupla. 3.5 Teste de Fernando e Vagner Analisando as seqüências, tivemos a idéia de um outro teste para verificar a dependência entre os números numa seqüência. Dado um vetor num com n elementos, tiramos uma cópia de num em num2, acrescentamos um 0 no final de num e outro 0 no final de num2. Fazemos então a diferença num − num2, armazenando isso num vetor residuo. Verificamos então como se comporta a distribuição desses resı́duos. Caso não haja uma dependência entre os números, esses resı́duos devem se distribuir de forma normal, sem nenhuma tendência especı́fica. Caso haja uma relação forte de dependência, espera-se uma distribuição assimétrica acentuada para um dos lados, assemelhando-se a uma exponencial ou similar. No caso de uma progressão aritmética, como a seqüência proposta no inı́cio dessa sessão (0, 1, 2, ..n) esse teste descobre na hora que a seqüência não é aleatória, pois todos os resı́duos serão iguais, de forma que quando plotarmos o histograma de resı́duos todos eles vão estar concentrados numa região só. 3.6 Testando uma seqüência de números aleatórios Por fim, vamos gerar números aleatórios usando nossa implementação do algoritmo linear congruente, e realizar testes para verificar como essa seqüência se comporta. Começamos gerando 50000 números entre 0 e 1, com uma semente não especificada. $ numeros 50000 | tail -50000 > num.txt Importamos esses dados no pacote estatı́stico R: > num <- scan("num.txt") E para fins de controle, geramos mais 50000 números aleatorios com distribuição uniforme em [0,1] mas não pela nossa implementação do algoritmo linear congruente e sim pela função equivalente em R: > rnum <- runif(50000,0,1) 8 Verificando as estatı́sticas descritivas observamos que os dados se assemelham: > summary(num) Min. 1st Qu. Median Mean 3rd Qu. Max. 6.105e-06 2.498e-01 4.993e-01 4.997e-01 7.494e-01 9.999e-01 > summary(rnum) Min. 1st Qu. Median Mean 3rd Qu. Max. 6.294e-05 2.492e-01 4.996e-01 4.984e-01 7.475e-01 1.000e+00 Verifiquemos os histogramas para essas duas seqüências: 9 Realizando um teste de Kolmorov-Smirnov nos dois conjuntos de dados, obtemos: One-sample Kolmogorov-Smirnov test data: num D = 0.0014, p-value = 1 One-sample Kolmogorov-Smirnov test data: rnum D = 0.0044, p-value = 0.3001 Indicando que a nossa implementação do método linear congruentes gerou números mais próximos da distribuição uniforme que o próprio gerador do R. Por fim, verifiquemos através do teste que nós implementamos o nı́vel de dependência entre os números subseqüentes. Geramos os vetores auxiliares, calculamos os resı́duos, e gerando os histogramas, obtemos: 10 O que mostra novamente que a nossa implementação do Algoritmo Linear Congruente foi bem sucedida, gerando a distribuição esperada de resı́duos. 11 4 4.1 Anexos Implementação do algoritmo do quadrado do meio /*********************************************************/ /* midsquare.c */ /* Este programa implementa o algoritmo de midsquare */ /* de von Neumann, para geracao de numeros aleatorios */ /* Fernando Henrique Ferraz Pereira da Rosa */ /* Vagner Aparecido Pedro Junior */ /* Data: 12/11/2002 */ /* */ /*********************************************************/ #include <stdio.h> #include <stdlib.h> int conta_digitos(int num); int midcut(int num,int nd); int main(int argc,char **argv) { int seed,nd,next,qtd,i; if (argc != 3) { printf("Uso: \n" " %s semente num, onde semente e’ um inteiro e \n" " num o numero de numeros aleatorios que voce quer gerar.\n", argv[0]); return(0); } sscanf(argv[1],"%d",&seed); sscanf(argv[2],"%d",&qtd); nd = conta_digitos(seed); printf("x[0] = %d\n",seed); for(i=1;i<=qtd;i++) { next = midcut(seed * seed,nd); printf("(%4d)^2 = %10d => x[%3d] = %d\n",seed,seed*seed,i,next); seed = next; } return(0); } /* conta_digitos: conta o numero de digitos * da representacao decimal de num */ int conta_digitos(int num) { 12 int cont; for(cont=1;num%10 != num;cont++) num /= 10; return(cont); } /* midcut: retorna o numero que esta no centro * de num, com nd digitos. * o cutpoint e’ para a esquerda, quando a divisao * nao e’ exata */ int midcut(int num,int nd) { int i,res,volta; char *string,*midpoint; nd *= 2; res = nd - conta_digitos(num); string = malloc((nd+1) * sizeof(char)); if (res > 0) { for(i=0;i<res;i++) { string[i] = ’0’; midpoint = &string[i+1]; } sprintf(midpoint,"%d",num); } else sprintf(string,"%d",num); nd /= 2; if (nd % 2 == 0) { /*divisao exata, primeiro tiramos nd/2 elementos do comeco *e depois nd/2 do final*/ for(i=0;i<nd/2;i++) string[i] = ’ ’; for(i=3*nd/2;i<2*nd+1;i++) string[i] = ’ ’; } else { /*como o corte e’ a esquerda, tiramos nd-nd/2 elementos da esquerda *e nd/2 elementos da direita.*/ for(i=0;i<nd-nd/2;i++) string[i] = ’ ’; for(i=2*nd-nd/2;i<2*nd;i++) string[i] = ’ ’; } sscanf(string,"%d",&volta); return(volta); } 13 4.2 Implementação do algoritmo linear congruente /*********************************************************/ /* linearcong.c */ /* Este programa implementa o algoritmo linear */ /* congruente. */ /* Fernando Henrique Ferraz Pereira da Rosa */ /* Vagner Aparecido Pedro Junior */ /* Data: 21/11/2002 */ /* */ /*********************************************************/ #include #include #include #include <stdio.h> <math.h> <stdlib.h> <time.h> int verifica_primo(int n); int sorteia_c(unsigned int aux,unsigned int m); int main (int argc,char **argv) { unsigned int n,aux,x0,c,next,m,i,a; srand(time(0)); if (argc != 2 && argc != 3) { printf("Uso: \n %s n [semente]\n" " com numero o n o numero de numeros aleatorios que voce quer\n" " gerar e opcionalmente uma semente.\n\n",argv[0]); return(0); } sscanf(argv[1],"%d",&n); aux = sqrt(n*10); while (verifica_primo(aux) == 0) aux--; m = aux*aux; a = aux+1; c = sorteia_c(aux,m); if (argc == 3) sscanf(argv[2],"%d",&x0); else x0 = ((double)rand()/RAND_MAX)*m; printf("x0 = %d\n" "a = %d\n" "c = %d\n" "aux = %d\n" "m = %d\n",x0,a,c,aux,m); 14 for (i=0;i<n;i++) { next = (a*x0 + c) % m; x0 = next; printf("%.20f\n",(double)next/m); } return(0); } /* verifica_primo: verifica se n e’ um numero e’ primo * Pos: retorna 1 se for primo. * retorna 0 caso contrario. */ int verifica_primo(int n) { int i,pausa; pausa = n/2; if (n == 0) return(0); for(i=2;i < pausa;i++) if (n % i == 0) return(0); return(1); } /* sorteia_c: acha um valor para c, tal que * c e m sao primos entre si. */ int sorteia_c(unsigned int aux,unsigned int m) { int num; num = ((double)rand()/RAND_MAX)*m/10; while (num == aux) num = ((double)rand()/RAND_MAX)*m/10; return(num); } 15 Referências [1] GRAYBEAL J.W. e POOCH W. U. Simulation: Principle and Methods. Winthrop Publishers, Inc. 1980. [2] KNUTH D.E. The Art of Computer Programming vol. 2 : Seminumerical Algorithms. Reading,Mass.: Addilson-Wesley, 1991. [3] ANDERSEN G. D. The π Search Page http://www.angio.net/pi/ piquery#award. [4] MAGALHÃES. M. N. e LIMA A.C.P. Noções de Probabilidade e Estatı́stica. Edusp, 2002. Sobre A versão eletrônica desse arquivo pode ser obtida em http://www.feferraz. net Copyright (c) 1999-2005 Fernando Henrique Ferraz Pereira da Rosa. É dada permiss~ ao para copiar, distribuir e/ou modificar este documento sob os termos da Licença de Documentaç~ ao Livre GNU (GFDL), vers~ ao 1.2, publicada pela Free Software Foundation; Uma cópia da licença em está inclusa na seç~ ao intitulada "Sobre / Licença de Uso". 16