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)) ``` <img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-2-1.png" width="60%" style="display: block; margin: auto;" /> --- # Transformação integral de probabilidade ### Definição (Transformação Integral de Probabilidade) <!-- Mood, pg. 202; Casella e Berger, pg. 49 --> .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\)`. <!-- Tem que ter os graficos aqui!!! --> --- # 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) ``` <img src="data:image/png;base64,#figures/03_RNG_inversa/unnamed-chunk-14-1.png" width="90%" style="display: block; margin: auto;" /> ```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[ ```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-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.