Geração de números aleatórios uniformes

class: center, middle, inverse, title-slide

# Geração de números aleatórios uniformes
### Fernando Mayer
### 2019-09-06 .footnotesize[(última atualização 2022-02-14)]

---

# Introdução

- Simulação computacional de processos estocásticos depende da
geração de números aleatórios (GNA).

- Inferência Bayesiana, jogos digitais, otimização estocástica...

Objetivos:

- Discutir a importância da geração de números aleatórios em Estatística.

- Descrever os principais algorítmos para GNA uniformes.

---
# A distribuição uniforme

Se uma variável aleatória `$X$` tem distribuição uniforme contínua entre
`$a$` e `$b$`, denotamos `$X \sim \text{U}(a, b)$`, e sua função de
densidade é

$$
f(x; a, b) = \frac{1}{b - a} \cdot I(a \leq x \leq b), \quad -\infty < a <
b < \infty.
$$

A função de probabilidade é

$$
F(x; a, b) = \Pr(X < x) = \int_{-\infty}^{x} f(x; a, b)\, \text{d}x =
`\begin{cases}
0, & x < a \\
\dfrac{x - a}{b - a}, & a \leq x < b\\
1, & x \geq b.
\end{cases}`
$$

A média e variância são funções dos parâmetros

$$
\text{E}(X) = \frac{a + b}{2} \quad\text{e}\quad \text{V}(X) =
\frac{(b - a)^2}{12}.
$$

---
# Importância da distribuição uniforme

- A distribuição tem pouca importância do ponto de vista de
modelagem.

- Porém, para simulação computacional, ela é central.

- Números aleatório da Uniforme são:

- O principal ingrediente para GNA de outras distribuições.

- Usados em várias aplicações de Monte Carlo
  (e.g. integração MC, inf. Bayesiana por MCMC, etc.)

- Utilizados em otimização estocástica.

- Empregados em jogos digitais (e.g. poker online).

---
# Geração de números aleatórios

## Números aleatórios reais

- São gerados por dispositivos via processos físicos.
  - Globo de sorteio da Mega-Sena e bingos.
  - Lançamento de um dado, moeda.
  - Roleta de cassino.

- Para uso em maior escala, são geralmente baseados em fenômenos
que geram um nível baixo de rúido.
  - Ruído termal.
  - Ruído fotoelétrico.
  - Ruído acústico.
  - Movimento browniano de partículas.

- Os sinais precisam ser traduzidos para números (intensidade).
- A distribuição da VA pode não ter forma conhecida.
- A VA pode ter correlação serial.
- São usados em criptografia.

Veja:

- https://www.random.org/
- https://www.fsf.org/blogs/gnu-press/neug-trng

---
# Geração de números aleatórios

## Números aleatórios reais

O pacote
[random](https://cran.r-project.org/web/packages/random/index.html) do R
faz uma conexão com o site https://www.random.org/ para trazer números
aleatórios "reais".

```r
library(random)
randomNumbers(n = 10, min = 1, max = 10, col = 10, base = 10) # com duplicação
```

```
#      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
# [1,]  4  3  2  2  4 10  6  8  7   3
```

```r
randomSequence(min = 1, max = 10, col = 10) # sem duplicação
```

```
#      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
# [1,]  5  2  4  1 10  3  7  8  6   9
```

```r
randomStrings(n = 4, len = 5, digits = TRUE, upperalpha = TRUE,
              loweralpha = TRUE, unique = TRUE)
```

```
#      V1     
# [1,] "Vc5n5"
# [2,] "8EqrC"
# [3,] "e9Ttb"
# [4,] "t5YiA"
```

Na prática, para aplicações em estatística, os números aleatório reais
não tem muita utilidade, já que suas propriedades são desconhecidas.

---
# Geração de números aleatórios

## Números pseudo aleatórios

.pull-left[
- São gerados por programas de computador, ou seja, algorítmos.

- Os número não são realmente aleatórios, pois, dado o algorítmo
gerador, pode-se prever os números (tem reprodutibilidade).

- Do ponto de vista Estatístico, os números são imprevisíveis.

- Bons GNA terão propriedades interessantes para aplicações
Estatística.
]

.pull-right[
- Uma sequência de números gerados por um algoritmo, completamente
  determinado pelo valor inicial do algoritmo: **semente** (*seed*) ou
  **chave** (*key*).

- Embora o determinismo possa não ser desejado, ele possui algumas
  vantagens:
  - Em simulação, permite repetir o mesmo processo usando a mesma
    semente
  - São mais eficientes, pois podem produzir grandes sequências em pouco
    tempo

- Por construção, os algoritmos são **periódicos**: a sequência se
  repete a longo prazo
  - No entanto, os algoritmos modernos possuem períodos tão longos, que
    a periodicidade pode ser desprezada
]

---
# Geração de números pseudo aleatórios

.pull-left[
Um gerador de números pseudo aleatórios é uma estrutura do tipo
$$
\Xi = (S, s_0, f, U, g)
$$
onde:

- `$S$` é um conjunto **finito** de estados
- `$s_0$` é o estado incial ou **semente**
- `$f: S \rightarrow S$` é uma função de transformação
- `$U$` é um conjunto **finito** de símbolos de saída
- `$g: S \rightarrow U$` é uma função de saída
]

.pull-right[
De maneira geral, a geração de números aleatórios é da seguinte forma:

1. Defina o estado inicial (semente) `$s_0$`
2. Para `$i = 1, 2, 3, \ldots$` faça `$s_i = f(s_{i-1})$`
3. Para cada `$i$`, calcule `$u_i = g(s_i)$`

O **período** é definido pelo menor inteiro `$p$`, de forma que
$$
s_{p+i} = s_i
$$

Por exemplo, a sequência

```r
6 9 0 7 6 9 0 7 6 9 0 7
```

possui período `$p=4$`, pois `$s_{4+1} = s_1$`, `$s_{4+2} = s_2$`, e assim
sucessivamente.

]

---
# Geração de números pseudo aleatórios

Dessa forma, **aparentemente** qualquer função `$f(\cdot)$` recursiva
poderia ser usada como um gerador. Por exemplo:

$$
`\begin{align*}
f(x_i) = a x_{i-1}
\end{align*}`
$$

Uma implementação computacional seria:

```r
naive <- function(n, x0, a) {
 x <- integer(n + 1)
 x[1] <- x0
 for(i in 2:length(x)) {
 x[i] <- (a * x[i - 1])
 }
 return(x[-1])
}
```

Gerando 10 números dessa sequência com `$s_0=2$` e `$a=345$`

```r
naive(n = 10, x0 = 2, a = 345)
```

```
#  [1] 6.900000e+02 2.380500e+05 8.212725e+07 2.833390e+10 9.775196e+12
#  [6] 3.372443e+15 1.163493e+18 4.014050e+20 1.384847e+23 4.777723e+25
```

- É um gerador **válido**, mas não necessariamente **bom**.
- Como saber se é um bom gerador?

---
# Geração de números pseudo aleatórios

O que faz um gerador ser bom?

- A sequência deve ser intuitivamente aleatória

- Essa aleatoriedade deve ser estabelecida teoricamente, ou passar por
  testes de aleatoriedade

- Deve-se conhecer alguma coisa sobre as propriedades teóricas do
  gerador

- A distribuição da VA é matematicamente conhecida
(e.g. distribuição Uniforme).

- O gerador produz VAs que não apresentam correlação serial, ou
seja, as realizações são independentes umas das outras.

> Isso indica que geradores *ad hoc* (como o exemplo acima) devem ser
> evitados.

---
# Geração de números pseudo aleatórios: métodos

- Existem vários algorítmos para GNA. No R, existem 7 opções para
o algorítmo de GNA.

- Veja `help(Random)` para consultar a respectiva documentação.

- Veja a seção *Random number generators (RNG)* na
[CRAN Task View: Probability
Distributions](http://cran-r.c3sl.ufpr.br/web/views/Distributions.html)

Algorítmos disponíveis no R:

- Wichmann-Hill.

- Marsaglia-Multicarry.

- Super-Duper.

- **Mersenne-Twister** (*default*).

- Knuth-TAOCP-2002.

- Knuth-TAOCP.

- L'Ecuyer-CMRG.

---
# Geração de números pseudo aleatórios: métodos

O método gerador padrão no R é o Mersenne-Twister. Por exemplo

```r
runif(10)
```

```
#  [1] 0.4976992 0.7176185 0.9919061 0.3800352 0.7774452 0.9347052 0.2121425
#  [8] 0.6516738 0.1255551 0.2672207
```
gera uma sequência de números uniformes baseada nesse algoritmo.

Através da função `set.seed(<k>)`, podemos especificar um número inteiro
`$k$` como semente ( `$s_0$` )

```r
set.seed(1)
runif(10)
```

```
#  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
#  [7] 0.94467527 0.66079779 0.62911404 0.06178627
```
Sempre que essa semente for usada, a mesma sequência será obtida.

```r
set.seed(1)
runif(10)
```

```
#  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
#  [7] 0.94467527 0.66079779 0.62911404 0.06178627
```

---
# Geração de números pseudo aleatórios: métodos

A função `set.seed()` também serve para especificar outro método
gerador. Por exemplo

```r
set.seed(1, kind = "Knuth-TAOCP-2002")
runif(10)
```

```
#  [1] 0.4701631 0.1585074 0.4828353 0.2083405 0.5682372 0.4294435 0.2085656
#  [8] 0.4578193 0.6270014 0.8554341
```
É diferente do padrão

```r
set.seed(1, kind = "Mersenne-Twister")
runif(10)
```

```
#  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968
#  [7] 0.94467527 0.66079779 0.62911404 0.06178627
```

- **Atenção**: alterar o método através da função `set.seed()` fará com
  que o último método seja utilizado durante toda a sessão
- Para mudar o método gerador permanentemente em uma sessão, use a função
`RNGkind()` (mais seguro, veja os detalhes)
- É aconselhável manter o padrão (Mersenne-Twister) por ser o que possui
  as melhores propriedades.

---
# Gerador congruencial linear (GCL)

O **método congruencial linear** foi proposto por Lehmer (1951), e
serviu como base para o desenvolvimento dos demais métodos.

É baseado na seguinte expressão recursiva

$$
`\begin{align*}
f(x_{i}) = (a x_{i-1} + c) \text{ mod } m
\end{align*}`
$$

onde:

- `$a \in \mathbb{N}^{+}$` é chamado de multiplicador,
- `$c \in \mathbb{N}^{+}$` é chamado de incremento,
- `$m \in \mathbb{N}^{+}$` é chamado de módulo.

Um caso particular é quando `$c=0$`:

$$
`\begin{align*}
f(x_{i}) = a x_{i-1} \text{ mod } m
\end{align*}`
$$

Nesse caso é chamado de **gerador congruencial multiplicativo**.

Para obter valores uniformes no intervalo `$[0,1]$` ainda é necessário
dividir os valores `$x_i$` por `$m$`:

$$
g(x_i) = \frac{x_i}{m}
$$

---
# Gerador congruencial linear (GCL)

.pull-left-40[
Implementação básica
.code90[

```r
rcl <- function(n, x0, m, a, c, unit = TRUE) {
 x <- integer(n + 1)
 x[1] <- x0
 for(i in 2:length(x)) {
 x[i] <- (a * x[i - 1] + c) %% m
 }
 if(unit) x <- x/m
 return(x[-1])
}
```
]
]
.pull-right-60[
Alguns resultados

.code90[

```r
rcl(n = 10, x0 = 1, m = 1e6, a = 1, c = 1)
```

```
#  [1] 2.0e-06 3.0e-06 4.0e-06 5.0e-06 6.0e-06 7.0e-06 8.0e-06
#  [8] 9.0e-06 1.0e-05 1.1e-05
```

```r
rcl(n = 10, x0 = 1, m = 1e6, a = 143, c = 1)
```

```
#  [1] 0.000144 0.020593 0.944800 0.106401 0.215344 0.794193
#  [7] 0.569600 0.452801 0.750544 0.327793
```

```r
## Congruencial multiplicativo
rcl(n = 10, x0 = 1, m = 1e6, a = 1, c = 0)
```

```
#  [1] 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06 1e-06
# [10] 1e-06
```

```r
rcl(n = 10, x0 = 1, m = 1e6, a = 143, c = 0)
```

```
#  [1] 0.000143 0.020449 0.924207 0.161601 0.108943 0.578849
#  [7] 0.775407 0.883201 0.297743 0.577249
```
]
]

---
# Gerador congruencial linear (GCL)

.pull-left-40[
- Naturalmente, a escolha de `$a$`, `$c$` e `$m$` deve ser feita com cuidado
- Não são quaisquer valores que gerarão uma boa sequência
- Em geral:
  - `$m$` deve ser um número inteiro grande
  - `$c$` e `$m$` devem ser relativamente primos
  - `$a$` de tal forma que `$ax \text{ mod } m \neq 0$`
- **O período será no máximo `$m$`**, por isso deve ser um número grande
  - No entanto, dependendo da combinação de valores das constantes, o
    período pode ser bem menor
]
.pull-right-60[

```r
## Periodo = 4 < m = 10
rcl(n = 12, x0 = 7, m = 10, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 6 9 0 7 6 9 0 7 6 9 0 7
```

```r
## Outros períodos
rcl(n = 12, x0 = 8, m = 10, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 3 8 3 8 3 8 3 8 3 8 3 8
```

```r
rcl(n = 12, x0 = 1, m = 10, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 4 5 2 1 4 5 2 1 4 5 2 1
```

```r
rcl(n = 12, x0 = 0, m = 10, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 7 6 9 0 7 6 9 0 7 6 9 0
```

```r
rcl(n = 12, x0 = 2, m = 10, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 1 4 5 2 1 4 5 2 1 4 5 2
```

```r
rcl(n = 12, x0 = 2, m = 100, a = 7, c = 7, unit = FALSE)
```

```
#  [1] 21 54 85  2 21 54 85  2 21 54 85  2
```
]

---
# Gerador congruencial linear (GCL)

Algumas combinações destes valores já foram estabelecidas na literatura.

Por exemplo, Park e Miller (1988) definiram `$m = 2^{31} - 1$`, `$a=7^5$` e
`$c=0$` como uma boa combinação.

```r
rcl(n = 12, x0 = 7, m = 2^31 - 1, a = 7^5, c = 0, unit = FALSE)
```

```
#  [1]     117649 1977326743  621132276  452154665 1566311569 1143995257
#  [7]  707192808 1615021558 1621510873 1165762081 1469983786 1365616214
```

```r
rcl(n = 12, x0 = 7, m = 2^31 - 1, a = 7^5, c = 0, unit = TRUE)
```

```
#  [1] 5.478458e-05 9.207645e-01 2.892373e-01 2.105509e-01 7.293707e-01
#  [6] 5.327143e-01 3.293123e-01 7.520530e-01 7.550748e-01 5.428503e-01
# [11] 6.845145e-01 6.359146e-01
```

Alguns outros métodos disponíveis.
.small[
| GNA             | m            | a           | c          |
|-----------------|--------------|-------------|------------|
| Knuth-Lewis     | `$2^{32}$`     | `$1664525$`   | `$1.01e9^1$` |
| Lavaux-Jenssens | `$2^{48}$`     | `$31167285$`  | `$1$`        |
| Haynes          | `$2^{64}$`     | `$6.36e17^2$` | `$1$`        |
| Marsaglia       | `$2^{32}$`     | `$69069$`     | `$0$`        |
| Park-Miller     | `$2^{31} - 1$` | `$16807$`     | `$0$`        |
|                 |              |             |            |
]

---
# Gerador Mersenne Twister

- Proposto por Matsumoto e Nishimura (1998)

- É o método padrão implementado no R
  - Criado para gerar números uniformes: `runif()`
  - Como os métodos para gerar valores de outras distribuições também
    depende do gerador da uniforme, funções como `rnorm()` também são
    afetadas (caso o gerador seja alterado). **Detalhes no próximo
    slide**.

- Atualmente é o que possui as melhores propriedades de um bom gerador

- O nome vem do fato de que seu período é um número [Primo de
  Mersenne](https://pt.wikipedia.org/wiki/N%C3%BAmero_de_Mersenne),
  i.e., um primo da seuinte forma: `$M_p = 2^p - 1$`

- O gerador trabalha internamente com **números binários**, por isso é
  altamente eficiente computacionalmente (por isso também, o método não
  é simples)

- O período da sequência é de
`$2^{19927} - 1 \approx 4.3 \times 10^{6001}$`

---
# Gerador Mersenne Twister

Como a mudança de método gerador interfere nas funções:
.pull-left.code80[

```r
RNGkind()
```

```
# [1] "Mersenne-Twister" "Inversion"        "Rejection"
```

```r
set.seed(1, kind = "Mersenne-Twister")
RNGkind()
```

```
# [1] "Mersenne-Twister" "Inversion"        "Rejection"
```

```r
runif(5)
```

```
# [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
```

```r
rnorm(5)
```

```
# [1]  1.2724293  0.4146414 -1.5399500 -0.9285670 -0.2947204
```

```r
set.seed(1, kind = "Knuth-TAOCP-2002")
RNGkind()
```

```
# [1] "Knuth-TAOCP-2002" "Inversion"        "Rejection"
```

```r
runif(5)
```

```
# [1] 0.4701631 0.1585074 0.4828353 0.2083405 0.5682372
```
]
.pull-right.code80[

```r
rnorm(5)
```

```
# [1] -0.1777910 -0.1059292  1.0600282  0.7387614 -0.2825928
```

```r
set.seed(1, kind = "Mersenne-Twister")
RNGkind()
```

```
# [1] "Mersenne-Twister" "Inversion"        "Rejection"
```

```r
runif(5)
```

```
# [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819
```

```r
rnorm(5)
```

```
# [1]  1.2724293  0.4146414 -1.5399500 -0.9285670 -0.2947204
```
]

---
# Testes de aleatoriedade

---
# Testes de aleatoriedade

Existem diversos tipos de teste para a aleatoriedade de uma sequência.

- [Dieharder: A Random Number Test
  Suite](http://webhome.phy.duke.edu/~rgb/General/dieharder.php): uma
  série de testes, adotado como "padrão" atualmente

- Testes visuais
  - Gráfico da sequência
  - Histograma
  - ACF (autocorrelação)

- Testes estatísticos
  - Teste de qui-quadrado
  - Teste não paramétrico de Kolmogorov-Smirnov

Gera uma sequência como exemplo

```r
x <- rcl(n = 1000, x0 = 1, m = 2^12, a = 125, c = 1, unit = TRUE)
```

---
# Testes visuais

```r
plot(x)
hist(x)
```

---
# Testes visuais

ACF para correlação serial

```r
acf(x)
```

---
# Teste de qui-quadrado

- Compara frequência observada por classes, com o que seria esperado
- Testa a hipótese nula de que as frequências observadas e esperadas são
  iguais

```r
## Divide os dados em 10 classes de igual tamanho
xc <- cut(x, breaks = seq(0, 1, 0.1), include.lowest = TRUE)
## Com 1000 dados, deveriam haver 100 observações em cada classe.
## Estas são as frequências observadas
table(xc)
```

```
# xc
#   [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] 
#       100        96        98        85       105        93        97       125 
# (0.8,0.9]   (0.9,1] 
#       107        94
```

```r
## A função chisq.test() faz o teste usando esta tabela
chisq.test(x = table(xc))
```

```
# 
# 	Chi-squared test for given probabilities
# 
# data:  table(xc)
# X-squared = 10.38, df = 9, p-value = 0.3206
```

- Não rejeita a hipótese nula de que as frequências observadas e
  esperadas são iguais

---
# Teste de Kolmogorov-Smirnov

.pull-left[
- Teste "de aderência" não paramétrico

- Compara a distribuição acumulada empírica dos dados com a acumulada de
  alguma distribuição de referência

- Calcula a maior distância entre estas duas distribuições

- Testa a hipótese nula de que a acumulada da distribuição empírica (dos
  dados) "adere" (é igual) à distribuição teórica

```r
ks.test(x, "punif")
```

```
# 
# 	One-sample Kolmogorov-Smirnov test
# 
# data:  x
# D = 0.031059, p-value = 0.2896
# alternative hypothesis: two-sided
```

- Não rejeita a hipótese nula de que a distribuição acumulada empírica
  segue uma uniforme
]
.pull-right[
Visualização da distribuição acumulada empírica da sequência gerada
(linha preta) com uma distribuição acumulada de uma uniforme gerada pela
função `punif()` (linha vermelha - usada como referência).

```r
plot(ecdf(x))
plot(ecdf(punif(seq(0, 1, length.out = 1000))),
     add = TRUE, col = 2)
```

<img src="data:image/png;base64,#figures/02_RNG_uniforme/unnamed-chunk-23-1.png" width="90%" style="display: block; margin: auto;" />
]

---
class: center, middle, inverse
# Desafio

---
# GNA de von Neumann

### Artigo

von Neumann, J. (1951). Various Techniques Used in Connection with
Random Digits. In the Monte Carlo Method (ed. A.S. Householder et al.,
36-38). **Nat. Bur. Standards Appl.** Math Ser. no. 12.

### Algorítmo

- Definir um número `$u_0$` de quatro digitos decimais e atribuir
`$i = 0$`.
- Calcular `$u_i^2$`. Agregar zeros à esquerda, quando necessário, para
manter representação com 4 digitos: `$u_i^2 = d_7 d_6... d_0$`, onde cada
`$d_j$` é um inteiro entre 0 e 9.
- Fazer `$u_{i+1} = d_5 d_4 d_3 d_2$`.
- Incrementar `$i$` fazendo `$i = i+1$`.
- Repetir 2-4 até obter a sequência de tamanho desejado.

Implemente e inspecione as propriedades deste GNA Uniformes.