Universidade Federal do Paraná
Curso de Estatística
CE 089 -
Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
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:
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:
O seguinte algoritmo permite gerar números aleatórios de uma distribuição de probabilidade caracterizada pela função densidade \(f\):
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} \]
##-----------------------------------------------------------------------------
## Seja X uma v.a. com f(x) = 1.5*x^2, -1<x<1. Simular valores desta.
## Gráfico de f(x).
curve(1.5*(x^2), -1, 1)
## Tem integral 1?
integrate(function(x) 1.5*(x^2), -1, 1)
## 1 with absolute error < 1.1e-14
##-----------------------------------------------------------------------------
## Considere g como 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(x) = 0.5,
## -1<x<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
## [1] -0.6343
## 2. Gerar u de uma uniforme padrão.
u <- runif(n=1); u
## [1] 0.7738
## 3. Comparar e decidir.
w <- f(y)/(M*g(y)); w
## [1] 0.4023
if(u<w){
x <- y
print("u<w, então valor aceito")
} else {
print("u>=w, então valor descartado")
}
## [1] "u>=w, então valor descartado"
##-----------------------------------------------------------------------------
## Gerar da gaussiana padrão usando a cauchy.
par(mfrow=c(2,1))
## 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)))
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 entando,
## 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.
w <- f(y)/(M*g(y)); w
if(u<w){
x <- y
print("u<w, então valor aceito")
} else {
print("u>=w, 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)
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.
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)
}
