Skip to content
Snippets Groups Projects
Commit 3fbc566a authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Lança aula de aceitação e rejeição.

parent 414d3e7e
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
tutoriais/04-rel-va.Rmd
+ 5
− 0
View file @ 3fbc566a
...@@ -379,6 +379,11 @@ lines(pbinom(xx, size = size, prob = prob) ~ xx, ...@@ -379,6 +379,11 @@ lines(pbinom(xx, size = size, prob = prob) ~ xx,
col = 2) col = 2)
``` ```
# Conteúdo complementar
* Univariate Distribution Relationships:
<http://www.math.wm.edu/~leemis/chart/UDR/UDR.html>.
# Próxima aula # Próxima aula
* Método Box-Muller para geração de números da Normal. * Método Box-Muller para geração de números da Normal.
......
tutoriais/06-met-ac-rej.Rmd 0 → 100644
+ 358
− 0
View file @ 3fbc566a
---
title: "Método da aceitação e rejeição"
author: Prof. Walmes M. Zeviani
date: 2017-08-25
bibliography: ../config/Refs.bib
csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
# Método da aceitação e rejeição
Quando se deseja gerar números aleatórios de distribuições de
probabilidade, nem sempre é possível aplicar o método da transformação
integral da probabilidade. Alguns situação são:
* Casos em que não se tem expressão para a função de densidade
acumulada, $F$, como é o caso da distribuição gaussiana e gama.
* Casos em que mesmo que se conheça $F(x)$ não possível chegar à uma
expressão para a função $F^{-1}$.
Se a função densidade de probabilidade for conhecida, $f$, então é
possível gerar números aleatórios dessa variável aleatória caracterizada
por $f$ pelo método da aceitação e rejeição. É necessário satisfazer
dois requisitos:
1. Ter um bom gerador de números uniformes.
2. Ter um bom gerador de números de uma variável aleatória
representada por uma distribuição $D$, escolhida de tal maneira que
existe uma constante $M$ tal que a densidade de $g$ que caracteriza
a distribuição $D$ satisfaz $f(x) \leq M g(x)$ para todo $x$ do
domínio de $f$.
O seguinte algoritmo permite gerar números aleatórios de uma
distribuição de probabilidade caracterizada pela função densidade $f$:
1. Gerar $y$ como sendo uma ocorrência da variável aleatória
representada por $D$.
2. Gerar $u$ como sendo uma ocorrência de uma uniforme padrão.
3. Se $u \leq f(y)/(M g(y))$ considerar que $x = y$ é um valor da
distribuição de probabilidade alvo cuja densidade é $f$, caso
contrário, descartar $y$.
4. Repetir até atingir o número de valores desejado $n$.
<!-- .
Por que o método gera números de acordo com a distribuição alvo?
$$
\begin{aligned}
\Pr(X\leq k)\, &= \Pr(Y\leq k| U\leq q)\newline
& \text{em que } q= \frac{f(y)}{M g(y)} \newline
&= \frac{\Pr(Y\leq k, U\leq q)}{\Pr(U\leq q)}\newline
& \text{passando para integrais} \newline
&= \frac{
\int_{-\infty}^{k} \left(\int_{0}^q h(u) \, \text{d}u \right) g(y) \, \text{d}y}{
\int_{-\infty}^{\infty} \left(\int_{0}^q h(u) \, \text{d}u \right) g(y) \, \text{d}y} \newline
& \text{tem-se que } h(u) =1 \text{ é a densidade da uniforme, então} \newline
&= \frac{
\int_{-\infty}^{k} q\, g(y) \, \text{d}y}{
\int_{-\infty}^{\infty} q\, g(y) \, \text{d}y} \newline
& \text{substituindo o valor de } q \newline
&= \frac{
\int_{-\infty}^{k} \frac{f(y)}{M g(y)}\, g(y) \, \text{d}y}{
\int_{-\infty}^{\infty} \frac{f(y)}{M g(y)}\, g(y) \, \text{d}y} \newline
&= \frac{
\frac{1}{M}\int_{-\infty}^{k} f(y) \, \text{d}y}{
\frac{1}{M}\int_{-\infty}^{\infty} f(y) \, \text{d}y} \newline
&= \frac{
\int_{-\infty}^{k} f(y) \, \text{d}y}{
1} \newline
\Pr(X\leq k)\, &= \int_{-\infty}^{k} f(y) \, \text{d}y.
\end{aligned}
$$
-->
```{r}
#-----------------------------------------------------------------------
# Seja X uma v.a. com f(x) = 1.5 * x^2, -1 < x < 1. Simular valores
# desta.
# Gráfico da f.d.p da v.a. X, f(x).
curve(1.5 * (x ^ 2), -1, 1)
# Tem integral 1?
integrate(function(x) 1.5 * (x^2), lower = -1, upper = 1)
#-----------------------------------------------------------------------
# Considere como g() a densidade de uma uniforme entre -1 e 1. Então se
# a base é 2, o altura deve ser 0.5 para ter produto 1, assim g(y) =
# 0.5, -1 < y < 1. Qual deve ser um valor de M para garantir que f(x) <=
# (M g(x)) para todo x dentro do [-1, 1]? O valor de M tem que ser 3,
# pois 3 * 0.5 <= sup_{x} f(x) = 1.5.
curve(1.5 * (x^2), -1, 1, col = 2)
curve(0.5 + 0 * x, add = TRUE, lty = 2)
curve(3 * 0.5 + 0 * x, add = TRUE, lty = 2, lwd = 2)
legend("bottomright",
legend = c("f(x)", "g(x)", "M g(x)"),
lty = c(1, 2, 2),
col = c(2, 1, 1),
lwd = c(1, 1, 2),
bty = "n")
# Criando os elementos necessários.
f <- function(x) 1.5 * x^2
g <- function(x) 0.5 + 0 * x
M <- 3
x <- NULL
# 1. Gerar y cuja densidade é g()
y <- runif(n = 1, -1, 1)
y
# 2. Gerar u de uma uniforme padrão.
u <- runif(n = 1)
u
# 3. Comparar e decidir.
r <- f(y)/(M * g(y))
r
if (u < r) {
x <- y
print("u < r então valor aceito.")
} else {
print("u >= r então valor descartado.")
}
```
# Simular números gaussianos usando a Cauchy
```{r}
#-----------------------------------------------------------------------
# Gerar da gaussiana padrão usando a cauchy.
par(mfrow = c(1, 2))
# Gaussiana e Cauchy.
curve(dnorm(x), -4, 4)
curve(dcauchy(x), add = TRUE, lty = 2)
# Razão entre elas.
curve(dnorm(x)/dcauchy(x), -4, 4)
abline(v = c(-1, 1), lty = 2)
layout(1)
dnorm(1)/dcauchy(1)
M <- sqrt(2 * pi/exp(1))
M
curve(M * dcauchy(x), -4, 4,
lty = 2,
ylim = c(0, M * dcauchy(0)),
ylab = "Densidade")
curve(dnorm(x), add = TRUE)
legend("topright",
legend = c("f(x)", "M g(x)"),
lty = c(1, 2),
bty = "n")
#-----------------------------------------------------------------------
# Aplicando o método.
M <- sqrt(2 * pi/exp(1))
f <- function(x) dnorm(x, 0, 1)
g <- function(x) dcauchy(x, 0, 1)
# A inversa da F para a Cauchy é: F^{-1} = tan(pi * (1 - u)). No
# entanto, vamos usar o gerador pronto rcauchy().
# 1. Gerar y cuja densidade é g, uma Cauchy.
y <- rcauchy(n = 1)
y
# 2. Gerar u de uma uniforme padrão.
u <- runif(n = 1)
u
# 3. Comparar e decidir.
r <- f(y)/(M * g(y))
r
if (u < r) {
x <- y
print("u < r, então valor aceito.")
} else {
print("u >= r, então valor descartado.")
}
#-----------------------------------------------------------------------
# Aplicar dentro de um laço condicional para obter N valores.
n <- 1 # Contador de valores aceitos.
l <- 1 # Contador de ciclos.
N <- 999 # Total de número à gerar.
# Vetor vazio.
# x <- vector(mode = "numeric", length = N)
x <- numeric(N)
while (n < N) {
y <- rcauchy(n = 1)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
if (u < w) {
x[n] <- y
n <- n + 1
}
l <- l + 1
}
plot(ecdf(x))
curve(pnorm(x), add = TRUE, col = 2)
# Taxa de aceitação.
n/l # Observada.
1/M # Teórica.
```
# Ilustrando o funcionamento do método
```{r, eval = FALSE}
#-----------------------------------------------------------------------
# Ilustrando.
n <- 1
N <- 25
curve(M * g(x), -6, 6, lty = 2, ylim = c(0, M * g(0)))
curve(f, add = TRUE)
while (n < N) {
y <- rcauchy(n = 1)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
if (u < w) {
n <- n + 1
rug(y)
points(y, u * M * g(y))
} else {
points(y, u * M * g(y), pch = 4)
}
Sys.sleep(1)
}
```
```{r, fig.width = 7, fig.height = 7, fig.show = "none", message = FALSE, warning = FALSE}
#-----------------------------------------------------------------------
# Gerar o gif.
library(animation)
loop <- function() {
M <- sqrt(2 * pi/exp(1))
f <- function(x) dnorm(x, 0, 1)
g <- function(x) dcauchy(x, 0, 1)
n <- 0L
l <- 1L
N <- 20L
q <- 0
set.seed(123)
x <- numeric(N)
while (n < N) {
curve(M * g(x), -6, 6,
lty = 2, ylim = c(0, M * g(0)),
ylab = NA,
xlab = NA,
n = 303)
curve(f, add = TRUE, n = 303)
legend("topleft",
legend = c("f(x)", "M g(x)"),
lty = c(1, 2),
bty = "n")
legend("topright",
title = "Decisão:",
legend = c("Aceitar", "Rejeitar"),
lty = 1,
col = c(3, 2),
bty = "n")
legend("right",
title = "Valor:",
legend = c("Aceito", "Rejeitado"),
pch = c(1, 4),
col = c(3, 2),
bty = "n")
y <- rcauchy(n = 1)
title(sub = substitute("Candidato: " * y == yy,
list(yy = round(y, 4))))
rug(y, lwd = 2, col = 1)
segments(y, 0, y, f(y), col = 3)
segments(y, f(y), y, M * g(y), col = 2)
segments(y - 0.2, 0, y + 0.2, 0, col = 1)
segments(y - 0.2, f(y), y + 0.2, f(y), col = 3)
segments(y - 0.2, M * g(y), y + 0.2, M * g(y), col = 2)
u <- runif(n = 1)
w <- f(y)/(M * g(y))
d <- u < w
if (d) {
x[n] <- y
n <- n + 1L
l <- l + 1L
q <- n/l
points(y, u * M * g(y), pch = 1, col = 3)
} else {
l <- l + 1L
q <- n/l
points(y, u * M * g(y), pch = 4, col = 2)
}
mtext(side = 3,
text = substitute(u * ' = ' * uu <= r * ' = ' * rr * '?',
list(uu = round(u, 4),
rr = round(w, 4))))
msg <- sprintf("iterações: %d \t aceitos: %d \t taxa ac.: %0.2f",
l, n, q)
title(main = list(msg, font = 1))
# Sys.sleep(0.5)
} # END while(n < N)
} # END loop().
# Verifica.
# loop()
saveHTML(loop(),
img.name = "met-ac-rej",
imgdir = "met-ac-rej",
htmlfile = "met-ac-rej.html",
autobrowse = FALSE,
verbose = FALSE,
title = "Método da aceitação e rejeição",
ani.width = 600,
ani.height = 600)
```
<!-- Ajusta automaticamente a altura e largura do iframe. -->
<script>
function resizeIframe(obj) {
obj.style.height = obj.contentWindow.document.body.scrollHeight + 'px';
obj.style.width = obj.contentWindow.document.body.scrollWidth + 'px';
}
</script>
<iframe
src="met-ac-rej.html"
style="display: block; margin: auto;"
frameborder="0"
scrolling="no"
onload="resizeIframe(this)">
</iframe>
# Conteúdo complementar
* Ler
[Generating Random Numbers and Random Variables](http://www.mathematik.uni-kl.de/~seifried/eng/teaching/teach_07ws/ozgur.pdf).
* Ler
[Playing with R: Rejection Sampling](http://playingwithr.blogspot.com.br/2011/06/rejection-sampling.html).
* Assitir
[Accept-Reject method](https://www.youtube.com/watch?v=KoNGH5PkXDQ).
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment