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

class: center, middle, inverse, title-slide

.title[
# Geração de números aleatórios não uniformes
]
.subtitle[
## Método da transformação integral de probabilidade
]
.author[
### Fernando Mayer
]
.date[
### 2019-09-10 .footnotesize[(última atualização 2024-05-21)]
]

---

# Introdução

- A GNA Uniformes são o ponto de partida para GNA de outras
distribuições.

- Do ponto de vista de simulação computacional, é importante
gerar números das distribuições de probabilidade mais comuns.

- Também é importante saber gerar números de distribuições desconhecidas
  ou da qual se sabe apenas seu núcleo.

### Objetivos

- Mostrar a GNA de VAs contínuas e discretas.

- Apresentar o método da transformação integral de probabilidades.

---
# Distribuição uniforme: propriedades

Se `$X \sim \text{U}(a, b)$`, então sua função de densidade é

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

A função de distribuição é

$$
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}`
$$

### Caso particular: `$X \sim \text{U}(0, 1)$`

$$
f(x) = 1, \quad 0 < x < 1
$$

$$
F(x) =
`\begin{cases}
0, & x < 0 \\
x, & 0 \leq x < 1\\
1, & x \geq b.
\end{cases}`
$$

---
# Distribuição uniforme: propriedades

```r
x <- seq(0, 1, 0.01)
fx <- dunif(x)
par(mfrow = c(1, 2))
plot(x, fx, type = "l")
plot(ecdf(x))
par(mfrow = c(1, 1))
```

---
# Transformação integral de probabilidade

### Definição (Transformação Integral de Probabilidade)

.content-box-yellow[
(a) Se `$X$` é uma variável aleatória com função de distribuição acumulada
`$F_X(x)$`, então `$U = F_X(X)$` é uniformemente distribuída no intervalo
`$(0,1)$`, isto é, `$P[U \leq u] = u$`, `$0 < u < 1$`.
]

.content-box-purple[
(b) Se `$U$` é uniformemente distribuída no intervalo `$(0,1)$`, então
`$X = F_{X}^{-1}(U)$` possui função de distribuição acumulada
`$F_X(\cdot)$`.
]

### Demonstração

.pull-left[
Parte (a)
$$
`\begin{align*}
P[U \leq u] &= P[F_X(X) \leq u] \\
 &= P[F_{X}^{-1}(F_X(X)) \leq F_{X}^{-1}(u)] \\
 &= P[X \leq F_{X}^{-1}(u)] \\
 &= F_X(F_{X}^{-1}(u)) \\
 &= u
\end{align*}`
$$
]

.pull-right[
Parte (b)
$$
`\begin{align*}
P[X \leq x] &= P[F_X^{-1}(U) \leq x] \\
 &= P[F_{X}(F_X^{-1}(U)) \leq F_{X}(x)] \\
 &= P[U \leq F_{X}(x)] \\
 &= F_X(x)
\end{align*}`
$$
]

---
# Transformação integral de probabilidade

### Observação

Se `$F_X$` for estritamente crescente, então `$F_X^{-1}$` é bem definida
por
$$
F_X^{-1}(u) = x \quad \Longleftrightarrow \quad F_X(x) = u.
$$
No entanto, por definição, `$F_X$` pode ser não decrescente (constante em
algum intervalo), então `$F_X^{-1}$` não é bem definida.

Considere que, se `$x_1$` e `$x_2$`, `$x_1 < x_2$`, são os limites (inferior e
superior, respectivamente) da parte não descrescente de uma função,
então qualquer `$x$` que satisfaça `$x_1 \leq x \leq x_2$` satisfaz
`$F_X(x) = u$`. Esse problema é evitado definindo `$F_X^{-1}(u)$` para
`$0 < u < 1$` por
$$
F_X^{-1}(u) = \inf \\{ x \, : \, F_X(x) \geq u \\}
$$
ou seja, o menor valor (ínfimo) de `$x$` no subconjunto de valores
formados por `$F_X(x) \geq u$`. Dessa forma, por essa definição, se houver
um `$x$` que satisfaça `$x_1 \leq x \leq x_2$`, então `$F_X^{-1}(u) = x_1$`.

---
# Transformação integral de probabilidade

- Como `$F_X^{-1}(u) = x$`, então `$F_X^{-1}(u)$` possui a mesma
  distribuição de `$X$`.

- Portanto, para gerar um valor aleatório de `$X$`, determine
`$u \sim \text{U}(0,1)$` e calcule a inversa `$F_X^{-1}(u)$`.

- O método é de fácil implementação, desde que `$F_X^{-1}(u)$` seja fácil
  de ser determinada.

- O método serve para gerar variáveis contínuas ou discretas.

### Resumindo:

1. Determine a função inversa `$F_X^{-1}(u)$`
2. Escreva um comando ou função para calcular `$F_X^{-1}(u)$`
3. Para cada valor aleatório a ser gerado:
  - (a) Obtenha um valor aleatório `$u$` de `$\text{U}(0,1)$`
  - (b) Calcule `$x = F_X^{-1}(u)$`

---
class: inverse, middle
# Variáveis contínuas

---
# Variáveis contínuas

Considere a seguinte distribuição de probabilidade:

$$
f(x) = 3x^2, \quad 0 < x < 1
$$

A função de distribuição é

$$
F(x) = \int_0^x 3v^2\, dv = x^3
$$

A inversa é

$$
u = x^3 \quad \Rightarrow \quad x = u^{1/3} = F_X^{-1}(u)
$$

Uma implementação simples seria

```r
## Gera 1000 valores da uniforme
n <- 1000
u <- runif(n)
## Calcula a inversa
x <- u^(1/3)
```

---
# Variáveis contínuas

.pull-left[

```r
## Histograma dos valores gerados
hist(x, prob = TRUE)
## Modelo teórico
curve(3 * x^2, from = 0, to = 1, add = TRUE, col = 2)
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-5-1.png" width="90%" style="display: block; margin: auto;" />
]
.pull-right[

```r
## Acumulada empírica
plot(ecdf(x))
## Acumulada teórica
curve(x^3, from = 0, to = 1, add = TRUE, col = 2)
legend("left", legend = c("Empírica", "Teórica"),
       lty = 1, col = 1:2, bty = "n")
```

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

---
# Variáveis contínuas

### Exponencial

.pull-left-70[
Uma v.a. `$X$` segue o **modelo exponencial** com parâmetro `$\lambda > 0$`
se sua densidade é dada por
$$
f(x) = \lambda e^{-\lambda x} \cdot I(x > 0),
$$

sendo `$\lambda > 0$`. Denota-se por `$X \sim \text{Exp}(\lambda)$`.

A função de distribuição da Exponencial é
$$
F(x) = \int_{0}^{x} f(v)\, \text{d}v = 1 - \exp\{(-\lambda x)\}.
$$

A inversa de `$F(x)$` é
$$
`\begin{align*}
u &= 1 - \exp\{-\lambda x\}\\
1 - u &= \exp\{-\lambda x\}\\
\log(1 - u) &= -\lambda x\\
-\frac{\log(1 - u)}{\lambda} &= x\\
\therefore \quad x &= F^{-1}(u) = -\frac{\log(1 - u)}{\lambda}.
\end{align*}`
$$
]
.pull-right-30[
Implementação em R

```r
randexp <- function(n, lambda) {
    u <- runif(n)
    x <- -log(1 - u)/lambda
    return(x)
}
```
]

---
# Variáveis contínuas

.pull-left.code80[

```r
*x <- randexp(n = 1000, lambda = 0.5)
hist(x, freq = FALSE)
curve(dexp(x, rate = 0.5), add = TRUE, col = 2, from = 0)
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-8-1.png" width="90%" style="display: block; margin: auto;" />
]
.pull-right.code80[

```r
plot(ecdf(x))
Fx <- function(x, lambda) 1 - exp(-lambda * x)
curve(Fx(x, lambda = 0.5), add = TRUE, col = 2, from = 0)
## curve(pexp(x, rate = 0.5), add = TRUE, col = 2, from = 0)
plot(ecdf(rexp(x, rate = 0.5)), add = TRUE, col = 3)
legend("right", legend = c("Empírica", "Teórica", "rexp"),
       lty = 1, col = 1:3, bty = "n")
```

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

---
# Variáveis contínuas

### Uniforme no intervalo `$[a,b]$`

.pull-left-70[
Se `$X \sim U[a,b]$` então
$$
f(x) = \frac{1}{b - a}, \quad  a \leq x \leq b
$$
Sua função de distribuição é
$$
F(x) = \int_{a}^{x} f(v)\, \text{d}v = \frac{x-a}{b-a}, \quad
a \leq x \leq b
$$
Sua inversa é
$$
u = \frac{x-a}{b-a} \quad \Rightarrow \quad
x = a + (b-1)u
$$
]
.pull-right-30[

```r
randunif <- function(n, a, b) {
    u <- runif(n)
    x <- a + (b - 1) * u
    return(x)
}
```
]

---
# Variáveis contínuas

### Uniforme no intervalo `$[a,b]$`

.pull-left.code80[

```r
*x <- randunif(n = 1000, a = 1, b = 10)
hist(x, freq = FALSE)
curve(dunif(x, 1, 10), add = TRUE, col = 2, from = 1)
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-11-1.png" width="90%" style="display: block; margin: auto;" />
]
.pull-right.code80[

```r
plot(ecdf(x))
Fx <- function(x, a, b) (x - a)/(b - a)
curve(Fx(x, a = 1, b = 10), add = TRUE, col = 2,
      from = 1, to = 10)
## curve(punif(x, 1, 10), add = TRUE, col = 2, from = 1, to = 10)
plot(ecdf(runif(x, 1, 10)), add = TRUE, col = 3)
legend("right", legend = c("Empírica", "Teórica", "runif"),
       lty = 1, col = 1:3, bty = "n")
```

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

---
exclude: true
# Variáveis contínuas

### Gama

.pull-left[
Se `$X \sim G[r,\lambda]$` então
$$
f(x) = \frac{1}{\Gamma(r)\lambda^r} x^{r-1} e^{\frac{x}{\lambda}},
\quad  x > 0
$$
Sua função de distribuição é
$$
F(x) = \int_{0}^{x} f(v)\, \text{d}v =
\frac{1}{\Gamma(r)} \gamma \left(r, \frac{x}{\lambda} \right)
$$
Sua inversa é
$$
u = \frac{1}{\Gamma(r)} \gamma \left(r, \frac{x}{\lambda} \right)
\quad \Rightarrow \quad
x = \frac{u^{-(r+1)}}{\Gamma(r)\lambda^r}e^{-u/\lambda}
$$
Onde `$\gamma(\cdot)$` é a [função gama
incompleta](https://en.wikipedia.org/wiki/Incomplete_gamma_function).
]
.pull-right.code80[

```r
randgama <- function(n, r, lambda) {
    u <- runif(n)
    x <- (u^(-(r + 1))/(gamma(r) * lambda^r)) * exp(-u/lambda)
    return(x)
}
```
]

---
exclude: true
# Variáveis contínuas

### Gama

.pull-left.code80[

```r
xx <- rgamma(1000, 1, 2)
hist(xx)
```

```r
*x <- randgama(n = 1000, r = 1/200, lambda = 1/200)
hist(x, freq = FALSE)
curve(dunif(x, 1, 10), add = TRUE, col = 2, from = 1)
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-14-2.png" width="90%" style="display: block; margin: auto;" />
]
.pull-right.code80[

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

---
class: inverse, middle
# Variáveis discretas

---
# Variáveis discretas

A GNA de VAs discretas é feito por busca direta no intervalo ao qual o
valor uniforme pertence na distribuição acumulada.

Se `$X$` é uma VA discreta e

$$
`\begin{align*}
\ldots < x_{i-1} < x_i < x_{i+1} < \ldots
\end{align*}`
$$

são os **pontos de descontinuidade** de `$F_X(x)$`, então a transformação
inversa é

$$
`\begin{align*}
F_X^{-1}(u) = x_i, \quad \text{ onde } \quad
F_X(x_{i-1}) < u \leq F_X(x_{i}).
\end{align*}`
$$

Portanto, para cada valor:

1. Gere `$u$` de `$U(0,1)$`
2. Define `$x_i$` onde `$F_X(x_{i-1}) < u \leq F_X(x_{i})$`

**Obs.:** o passo 2 pode ser difícil para algumas distribuições.

---
# Variáveis discretas

### Bernoulli

.pull-left[
Se `$X \sim Ber(p)$` então
$$
P(X=x) = p^x (1-p)^{1-x}, \quad  x=0,1
$$
Sua função de distribuição é
$$
F_X(x) =
`\begin{cases}
1-p, & x = 0 \\
1,   & x = 1
\end{cases}`
$$
Portanto, se `$p=0.4$`, por exemplo,
$$
F_X^{-1}(u) =
`\begin{cases}
0,   & u \leq 0.6 \\
1, & u > 0.6
\end{cases}`
$$
]
.pull-right[

```r
p <- 0.4; q <- 1 - p
x <- rbinom(1000, size = 1, prob = p)
plot(ecdf(x), verticals = TRUE, col.vert = "red")
```

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

---
# Variáveis discretas

### Bernoulli

.pull-left[

```r
## Uma implementação elegante
n <- 1000
u <- runif(n)
*x <- as.integer(u > 0.6)
```

```r
## Uma forma mais longa (e menos eficiente)
x2 <- integer(n)
p <- 0.4
q <- 1 - p
for(i in 1:length(x2)) {
    if(u[i] < q) {
        x2[i] <- 0L
    } else {
        x2[i] <- 1L
    }
}
identical(x, x2)
```

```
# [1] TRUE
```
]

.pull-right[

```r
## Usando ifelse vetorizado
x3 <- ifelse(u < q, 0L, 1L)
identical(x2, x3)
```

```
# [1] TRUE
```

```r
## Usando seleção condicional
x4 <- integer(n)
x4[u >= q] <- 1L
identical(x3, x4)
```

```
# [1] TRUE
```
Outras formas de gerar da Bernoulli:

```r
rbinom(n, size = 1, prob = p)
sample(c(0, 1), size = n, replace = TRUE,
       prob = c(0.6, 0.4))
```
]

---
# Variáveis discretas

### Bernoulli

.pull-left.code80[

```r
addmargins(prop.table(table(x)))
```

```
# x
#     0     1   Sum 
# 0.578 0.422 1.000
```

```r
plot(prop.table(table(x)), ylim = c(0, 0.7), type = "h")
lines(x = c(0, 1) + .01, y = c(.6, .4),
      col = "red", type = "h", lwd = 2)
legend("topright", legend = c("Empírica", "Teórica"),
       lty = 1, col = 1:2, bty = "n")
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-21-1.png" width="75%" style="display: block; margin: auto;" />
]
.pull-right.code80[

```r
plot(ecdf(x))
curve(pbinom(x, size = 1, p = p), add = TRUE, type = "s", col = 2)
plot(ecdf(rbinom(x, size = 1, prob = p)), add = TRUE, col = 3)
legend("right", legend = c("Empírica", "Teórica", "rbinom"),
       lty = 1, col = 1:3, bty = "n")
```

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

---
# Variáveis discretas

### Uniforme discreta

.pull-left[
Se `$X$` tem distribuição uniforme discreta com `$k$` valores, o seu
suporte é o conjunto `$x \in \{1, 2, \ldots, k\}$`, `$k \geq 2$`. A função
de probabilidade é
$$
p(x) = \frac{1}{k} \cdot I(x \in \{1, 2, \ldots, k\}).
$$

A função de probabilidade acumulada é
$$
F(x) = \frac{x}{k}.
$$

Gera-se números da VA uniforme `$U \sim \text{U}(0, 1)$` e obtém-se
números da VA `$X$` de distribuição uniforme discreta usando

$$
X =
`\begin{cases}
1, & \text{se } U < 1/k \\
2, & \text{se } 1/k \leq U < 2/k \\
\vdots & \vdots \\
k, & \text{se } (k-1)/k \leq U < 1. \\
\end{cases}`
$$
]
.pull-right[
Um exemplo com `$X \sim \text{U}(k=10)$`

```r
x <- sample(1:10, size = 1000, replace = TRUE)
plot(ecdf(x), verticals = TRUE, col.vert = "red")
```

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

---
# Variáveis discretas

### Uniforme discreta

Uma implementação "ingênua" seria

```r
n <- 1000
*u <- runif(n)
x <- integer(n)
k <- 10
x[u < 1/k] <- 1
x[u >= 1/k & u < 2/k] <- 2
## Seria muito trabalhoso
```

Uma implementação usando conceitos de objetos do R

```r
xc <- cut(u, breaks = c(0, seq(1/k, (k-1)/k, 0.1), k/k),
          include.lowest = TRUE)
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] 
#        85       112       112       102       113        83        96        92 
# (0.8,0.9]   (0.9,1] 
#       101       104
```

```r
*x <- as.integer(xc)
```

---
# Variáveis discretas

### Uniforme discreta

.pull-left.code80[

```r
addmargins(prop.table(table(x)))
```

```
# x
#     1     2     3     4     5     6     7     8     9    10   Sum 
# 0.085 0.112 0.112 0.102 0.113 0.083 0.096 0.092 0.101 0.104 1.000
```

```r
plot(prop.table(table(x)), ylim = c(0, 0.15), type = "h")
lines(x = (1:10) + .1, y = rep(.1, 10),
      col = "red", type = "h", lwd = 2)
legend("topright", legend = c("Empírica", "Teórica"),
       lty = 1, col = 1:2, bty = "n")
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-26-1.png" width="75%" style="display: block; margin: auto;" />
]
.pull-right.code80[

```r
plot(ecdf(x))
Fx <- function(x, k) x/k
curve(Fx(x, k = 10), add = TRUE, col = 2, from = 1, to = 10, n = 10, type = "s")
plot(ecdf(sample(1:10, size = 1000, replace = TRUE)), add = TRUE, col = 3)
legend("right", legend = c("Empírica", "Teórica", "sample"),
       lty = 1, col = 1:3, bty = "n")
```

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

---
# Variáveis discretas

### Algoritmo de busca sequêncial

Uma forma mais geral de gerar valores aleatórios de **qualquer**
distribuição discreta.

Defina o vetor `$p$` de probabilidades (implicita ou explicitamente) e
gere `$u \sim U(0,1)$`. Com isso, faça:

1. Atribua `$x_0 = 1$`, `$p_0 = P[X = x_0]$` e `$F(x) = p_0$`
2. Enquanto `$u > F(x)$` faça:
  - `$x = x+1$`
  - `$F(x) = F(x) + p_x$`
3. Retorne `$x$`

---
# Variáveis discretas

### Algoritmo de busca sequêncial

.pull-left-60[
No caso da distribuição Uniforme discreta acima, com `$k=10$`, o algoritmo
seria então

1. Atribua `$x_0 = 1$`, `$p_0 = P[X = x_0] = 1/10 = 0.1$` e `$F(x) = p_0$`
2. Enquanto `$u > F(x)$` faça:
  - `$x = x+1$`
  - `$F(x) = F(x) + p_x$`
3. Retorne `$x$`
]
.pull-right-40[
Uma implementação básica seria

```r
randd <- function(x, px) {
    u <- runif(1)
    x <- 1
    Fx <- px[x]
    while(u > Fx) {
        x <- x + 1
        Fx <- Fx + px[x]
    }
    return(x)
}
```
]

---
# Variáveis discretas

### Algoritmo de busca sequêncial

Note que esse algoritmo possui três "restrições":
- o vetor de probabilidades `$p$` deve ser passado explicitamente,
- gera apenas valores de VAs que começam em 1,
- gera apenas **uma realização** da variável aleatória.

```r
k <- 10
x <- 1:k
px <- rep(1/k, length(x))
randd(x = x, px = px)
```

```
# [1] 4
```

Para gerar mais valores podemos usar a função `replicate()`.

```r
replicate(n = 10, randd(x = x, px = px))
```

```
#  [1] 9 5 8 3 3 1 3 2 1 2
```

---
# Variáveis discretas

Considere a seguinte distribuição de probabilidade

.pull-left-30[

| `$X$`| `$P(X)$`|
|---:|------:|
|   1|   0.20|
|   2|   0.15|
|   3|   0.25|
|   4|   0.40|
]
.pull-right-70[

```r
x <- 1:4
px <- c(.2, .15, .25, .4)
plot(x, cumsum(px), type = "s", ylim = c(0, 1), axes = FALSE)
axis(1, at = 1:4, labels = 1:4); axis(2); box()
points(x, cumsum(px), pch = 16)
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-32-1.png" width="70%" style="display: block; margin: auto;" />
]

---
# Variáveis discretas

Para gerar valores aleatórios desta distribuição, podemos usar a mesma
função definida acima, apenas passando os vetores de `$X$` e `$p_x$`
adequados.

.pull-left-40[

```r
randd(x = x, px = px)
```

```
# [1] 4
```

```r
replicate(10, randd(x = x, px = px))
```

```
#  [1] 1 4 4 1 1 2 3 3 4 1
```
]
.pull-right-60[

```r
xr <- replicate(n = 10000, randd(x = x, px = px))
plot(ecdf(xr), col = 2); lines(x, cumsum(px), type = "s")
points(x, cumsum(px), pch = 16)
legend("left", legend = c("Empírica", "Teórica"),
       lty = 1, col = 1:2, bty = "n")
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-34-1.png" width="70%" style="display: block; margin: auto;" />
]

---
# Variáveis discretas

### Algoritmo de busca sequencial - Poisson

.pull-left[
- No caso de uma distribuição como a Poisson, a aplicação direta deste
  algoritmo não é possível pois:
  - a Poisson assume um número infinito de valores discretos
  ( `$0, 1, 2, \ldots$` ),
  - as probabilidades dependem do parâmetro da distribuição.

- Assim, não é possível definir previamente os valores de `$X$` e `$p_x$`.

- Uma solução seria então calcular as probabilidades implicitamente
  através de alguma forma de **recursão**.
]
.pull-right[
Sabendo que

$$
P[X=x] = \frac{e^{-\lambda}\lambda^x}{x!}
$$

Então podemos determinar a probabilidade de `$x+1$` como sendo

$$
`\begin{align*}
P[X = x+1] &= \frac{e^{-\lambda}\lambda^{x+1}}{(x+1)!} \\
  &= \frac{e^{-\lambda}\lambda^{x}\lambda}{(x+1)x!} \\
  &= \frac{e^{-\lambda}\lambda^{x}}{x!} \frac{\lambda}{(x+1)} \\
  &= P[X=x] \frac{\lambda}{(x+1)} \\
\end{align*}`
$$

Dessa forma, podemos usar essa **atualização de probabilidade
recursiva** no algoritmo de busca sequencial.
]

---
# Variáveis discretas

### Algoritmo de busca sequêncial - Poisson

Uma adaptação do algoritmo geral para a Poisson seria então

.pull-left[
1. Gera `$u \sim U(0,1)$`
2. Atribua `$x_0 = 0$`, `$p_0 = P[X = x_0]$`, `$F(x) = p_0$` e `$p=p_0$`
2. Enquanto `$u > F(x)$` faça:
  - `$x = x+1$`
  - `$p = \frac{\lambda}{x} p$`
  - `$F(x) = F(x) + p$`
3. Retorne `$x$`
]
.pull-right[

```r
randpois <- function(lambda) {
    u <- runif(1)
    x <- 0
    p0 <- exp(-lambda)
    Fx <- p0
    p <- p0
    while(u > Fx) {
        x <- x + 1
        p <- (lambda/x) * p
        Fx <- Fx + p
    }
    return(x)
}
```
]

---
# Variáveis discretas

.pull-left-40[

```r
randpois(lambda = 5)
```

```
# [1] 5
```

```r
replicate(10, randpois(lambda = 5))
```

```
#  [1] 2 1 1 7 6 3 3 6 1 3
```
]
.pull-right-60[

```r
x <- replicate(10000, randpois(lambda = 5))
plot(ecdf(x))
curve(ppois(x, lambda = 5), add = TRUE, type = "s", col = 2)
plot(ecdf(rpois(x, lambda = 5)), add = TRUE, col = 3)
legend("right", legend = c("Empírica", "Teórica", "rpois"),
       lty = 1, col = 1:3, bty = "n")
```

<img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-37-1.png" width="80%" style="display: block; margin: auto;" />
]

---
class: inverse, middle
# Exercícios

---
# Exercícios

- Gerar valores aleatórios das distribuições **contínuas** especificadas
abaixo
  - `$F(x) = 1-(x-1)^{2}, \quad 0 < x < 1.$`
  - `$F(x) = x^{\theta}, \quad 0 < x < 1, \theta > 1.$`
  - `$F(x) = 1-\exp\{-x^2/2\tau^2\}, \quad x > 0, \tau > 0.$`
  - `$f(x) = \frac{1}{2} \sin(x), \quad 0 < x < \pi.$`
  - `$f(x) = \frac{\sec^2(x)}{\sqrt{3}} , \quad 0 < x < \pi/3.$`

> Dica: use programas de matemática simbólica para verificar os
resultados. Wolfram Alpha: <http://www.wolframalpha.com/>.

- Implemente o algoritmo de busca sequencial para gerar valores de uma
distribuição binomial.