Motivação

Muitas vezes, queremos simular valores de uma variável aleatória \(X\) cuja distribuição é conhecida, mas a função de distribuição acumulada (FDA) inversa, \(F_X^{-1}(u)\), não possui forma analítica fechada (nem uma aproximação adequada que permita seu cálculo), inviabilizando o uso do método da transformação inversa.

Exemplos:

distribuição Beta;
distribuições multimodais;
inúmeras distribuições a posteriori em Inferência Bayesiana.

O Método da Aceitação-Rejeição surge como uma alternativa extremamente poderosa. A ideia central é simples: se não conseguimos gerar diretamente da distribuição alvo, geramos amostras de uma distribuição auxiliar e aceitamos essas amostras com uma probabilidade proporcional à sua importância sob a distribuição alvo.

Método da Aceitação-Rejeição para V.A. Contínuas

Seja \(X\) uma variável aleatória contínua com função de densidade alvo \(f_X\) que desejamos simular. Suponha que dispomos de uma densidade candidata \(q\) da qual sabemos simular eficientemente.

O método da aceitação-rejeição é um método de reamostragem, baseado em dois passos:

gerar valores da distribuição auxiliar conhecida \(q\); e
usar um mecanismo para ``corrigir’’ os valores, fazendo com que representem (ao menos de forma aproximada) a distribuição de interesse \(f_X\)

Em particular, o método para corrigir os valores exige que exista uma constante \(k \ge 1\) tal que: \[f_X(x) \le k \cdot q(x), \quad \forall x \text{ tal que } f_X(x) > 0.\]

Neste caso, o mecanismo de geração consiste em gerar um valor \(y\) de \(q\)e aceitar o valor gerado como sendo de \(f_X\) com probabilidade \(\displaystyle \frac{f_X(y)}{k \cdot q(y)}\) que representa a proporção de pontos gerados abaixo de \(f_X\).

Intuição

Algoritmo

Gere um valor \(y\) da v.a. \(Y\) com função de densidade de probabilidade \(q\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{f_X(y)}{k\cdot q(y)}\), tome \(x = y\); caso contrário, retorne ao passo 1.

Repita os passos 1, 2 e 3 quantas vezes forem necessárias até obter a amostra do tamanho desejado.

Proposição

Dado que um ponto foi gerado e aceito de acordo com o Método da Aceitação-Rejeiçao (conforme descrito acima), este é um representante do modelo com função de densidade de probabilidade \(f_X\).

Demonstração: (Casella & Robert, 2005)

Queremos provar que a densidade condicional de \(Y\), dado que foi aceito, é exatamente a densidade alvo \(f_X\).

Seja \(A\) o evento de aceitação, isto é, \(A = \left\{ U \le \frac{f_X(Y)}{k \cdot q(Y)} \right\}\).

De fato, a probabilidade de aceitação marginal, \(P(A)\), é dada por: \[\begin{align} P(A) & = \int_{-\infty}^{\infty} P\left( Y=y, A \right) dy \\ & = \int_{-\infty}^{\infty} P(Y=y)P(A|Y=y) dy \\ & = \int_{-\infty}^{\infty} P(Y=y)P\left(\left.U \le \frac{f_X(Y)}{k \cdot q(Y)}\right|Y=y\right) dy \\ & \underbrace{=}_{\star} \int_{-\infty}^{\infty} q(y) \dfrac{f_X(y)}{k \cdot q(y)} dy = \int_{-\infty}^{\infty} \dfrac{f_X(y)}{k} dy = \dfrac{1}{k} \end{align}\] \(^\star\) uma vez que \(U \sim \text{U}(0,1) \Rightarrow F_U(u) = u, \;0<u<1\).

Agora, calculamos a função de distribuição acumulada condicionada à aceitação: \[\begin{align} P(Y \le x \mid A) &= \dfrac{P(Y \le x, A)}{P(A)} = k \cdot P\left( Y \le x, A \right) = k \cdot P\left( Y \le x, U \le \frac{f_X(Y)}{k \cdot q(Y)} \right) \\ & = k \int_{-\infty}^x q(y) \left[\int_0^{\frac{f_X(y)}{k\cdot q(y)}}1du\right] dy = k \int_{-\infty}^x q(y) \left[\frac{f_X(y)}{k\cdot q(y)} \right] dy \\ & = k \int_{-\infty}^x \left[\frac{f_X(y)}{k} \right] dy = \int_{-\infty}^x f_X(y) dy = F_X(y); \end{align}\] ou seja, como a acumulada condicional é \(F_X(y)\), a variável aceita distribui-se exatamente como \(f_X\). \(\quad\quad \square\).

Observação

Assim como no caso discreto, o número de tentativas até obter uma aceitação segue uma distribuição Geométrica com probabilidade de sucesso \(p = 1/k\). Logo, o número esperado de iterações por amostra aceita é exatamente \(k\). Para otimizar o algoritmo, devemos escolher \(q(y)\) de tal forma que o valor de \(k\) seja o menor possível (\(k \to 1\)).

Exemplo 1

Deseja-se gerar valores de uma v.a. com f.d.p. dada por

\[ f_X(x)=2x, \qquad 0<x<1. \] Supondo que sabemos gerar apenas da distribuição uniforme padrão, simule valores de \(X\) aravés do Método da Aceitação-Rejeição.

Solução:

Como \(2x\le 2\) \(0<x<1\), escolhendo \(q(x)=1\), \(0<x<1\), temos \[ \frac{f_X(x)}{q(x)} = \frac{2x}{1} \leq 2, \forall x; \] logo, tomaremos \(k=2\).

Algoritmo 1.1

Gere um valor \(y\) de \(Y \sim \text{U}(0,1)\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{2y}{2 \cdot 1} = y\), tome \(x = y\); caso contrário, retorne ao passo 1.

n = 10000
y = runif(n)

x = NULL
for (i in 1:n){
  u = runif(1)
  x[i] = ifelse(u < y[i],y[i],NA)  
}

# número de pontos aceitos
n.gerados = length(x[is.na(x)==FALSE]) 
n.gerados

## [1] 5008

x = x[is.na(x)==FALSE]
hist(x,prob=TRUE,main="",ylab=expression(f[X](x)))
curve(2*x,add=T,col=2,lwd=2,lty=2)

Vamos refazer o exemplo, com uma escolha mais arbitrária de \(k\) para verificar o efeito na taxa de aceitação dos pontos propostos.

Solução:

Algoritmo 1.2

Gere um valor \(y\) de \(Y \sim \text{U}(0,1)\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{2y}{k \cdot 1} = y\), tome \(x = y\); caso contrário, retorne ao passo 1.

n = 10000
y = runif(n)

k = 5

x = NULL
for (i in 1:n){
  u = runif(1)
  x[i] = ifelse(u < 2*y[i]/k,y[i],NA)  
}

# número de pontos aceitos
n.gerados = length(x[is.na(x)==FALSE]) 
n.gerados

## [1] 1952

x = x[is.na(x)==FALSE]
hist(x,prob=TRUE,main="",ylab=expression(f[X](x)))
curve(2*x,add=T,col=2,lwd=2,lty=2)

Alguns erros comuns na implementação do Método de Aceitação-Rejeição:

Escolher \(q(x)\) com suporte menor que o de \(f_X(x)\).
Esquecer de calcular corretamente a constante \(k\).
Escolher um \(k\) excessivamente grande.
Confundir a densidade alvo com a distribuição proposta.

Visitando algumas Distribuições Contínuas

Exemplo 2: Distribuição Beta

Simular valores de uma v.a. \(X \sim \text{Beta}(2,4)\), cuja função de densidade de probabilidade é dada por \[f_X(x) = 20x(1-x)^3, \; x \in (0,1).\]

Solução:

Considere a densidade auxiliar \(U(0,1)\), ou seja, \(q(x)=1\), \(x \in (0,1)\).

Qual o menor valor da constante \(k\) tal que \(f_X(x) \leq kq(x)\), ou, equivalentemente, \(\dfrac{f_X(x)}{q(x)} \leq k\)?

\[\dfrac{d}{dx}\left(\dfrac{f_X(x)}{q(x)}\right) = 0 \quad \Leftrightarrow \quad x=\dfrac{1}{4} \quad \Rightarrow \quad k = \dfrac{135}{64};\]

o que leva a uma probabilidade de aceitação \[\dfrac{f_X(y)}{k\cdot q(y)} = \dfrac{20y(1-y)^3}{\dfrac{135}{64}}\]

Algoritmo 2.1

Gere um valor \(y\) de \(Y \sim \text{U}(0,1)\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{f(y)}{k\cdot q(y)} = \dfrac{20y(1-y)^3}{\dfrac{135}{64}}\), tome \(x = y\); caso contrário, retorne ao passo 1.

f = function(x){
  20*x*(1-x)^3  
}  

# Conferindo o valor de k calculado  
optimize(f,lower=0,upper=1,maximum=TRUE)

## $maximum
## [1] 0.2499993
## 
## $objective
## [1] 2.109375

# Note que f e a f.densidade da Beta(2,4)
curve(f,0,1,lwd=2)
curve(dbeta(x,shape1=2,shape2=4),0,1,col=2,lwd=2,lty=2,add=TRUE)

# Vamos usar a proposta uniforme conforme descrito no exemplo:
x.gerado = NULL
x.aceito = NULL
  
M = 1000
j = 0
cont = 0  # apenas para sabermos quantos pontos precisaram ser gerados
while(j<M){
  cont = cont+1
  x.gerado[cont] = runif(1,0,1)
  u = runif(1,0,1)
  p.aceit = 20*x.gerado[cont]*(1-x.gerado[cont])^3/(135/64)
  if(u<p.aceit){
    j=j+1
    x.aceito[j] = x.gerado[cont]}     
}

taxa.aceitacao = length(x.aceito)/length(x.gerado)
taxa.aceitacao

## [1] 0.490918

hist(x.aceito,freq=FALSE,main="",ylab=expression(f[X](x))) 
curve(dbeta(x,shape1=2,shape2=4),0,1,col=2,lwd=2,lty=2,add=TRUE)

Vamos refazer o exemplo fingindo não perceber que a densidade de interesse era beta. Neste caso, vamos utilizar a função ``dbeta” para simular valores de uma distribuição beta como densidade auxiliar, digamos, \(Y \sim \text{Beta}(2,2)\).

Solução:

Considere agora a densidade auxiliar \(\text{Beta}(2,2)\), ou seja, \(q(x)=6x(1-x)\), \(x \in (0,1)\).

Algoritmo 2.2

Gere um valor \(y\) de \(Y \sim \text{Beta}(2,2)\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{f(y)}{k\cdot q(y)} = \dfrac{20y(1-y)^3}{k \cdot 6y(1-y)} = \dfrac{10}{3k}(1-y)^2\), tome \(x = y\); caso contrário, retorne ao passo 1.

f = function(x){
  20*x*(1-x)^3  
}

q = function(x,a,b){
  dbeta(x,a,b)  
}

a = 2
b = 2
razao = function(x,a=a,b=b){
  f(x)/q(x,a,b)
}

# Obtendo k
k = optimize(razao,a=a,b=b,lower=0,upper=1,maximum=TRUE)$objective
k

## [1] 3.332893

# Comparando f e kq
curve(k*dbeta(x,shape1=a,shape2=b),0,1,col=3,lwd=2,lty=3,ylab=expression(f[X](x)))
curve(f,0,1,lwd=2,add=TRUE)

# Vamos usar a proposta uniforme conforme descrito no exemplo:
x.gerado = NULL
x.aceito = NULL

M = 1000
j = 0
cont = 0  # apenas para sabermos quantos pontos precisaram ser gerados
while(j<M){
  cont = cont+1
  x.gerado[cont] = rbeta(1,a,b)
  u = runif(1,0,1)
  p.aceit = (10/(3*k))*(1-x.gerado[cont])^2
  if(u<p.aceit){
    j=j+1
    x.aceito[j] = x.gerado[cont]}     
}

taxa.aceitacao = length(x.aceito)/length(x.gerado)
taxa.aceitacao

## [1] 0.30012

hist(x.aceito,freq=FALSE,main="",ylab=expression(f[X](x))) 
curve(dbeta(x,shape1=2,shape2=4),0,1,col=2,lwd=2,lty=2,add=TRUE)

Exemplo 3: Distribuição Cauchy

Gerar valores da distribuição Cauchy, cuja função de densidade de probabilidade é dada por \[\displaystyle f_X(x) = \frac{1}{\pi\sigma}\left[1+\left(\frac{x-\mu}{\sigma}\right)^2\right]^{-1}.\]

Vamos considerar duas densidades auxiliares: - uma distribuição uniforme \(\text{U}(\mu-c\sigma\;,\;\mu+c\sigma)\)}}, onde \(c\) será tal que a maior parte das caudas da distribuição sejam contempladas; e - uma distribuição t-Student {\(t_{\nu}(\mu,\sigma)\).

Solução:

Usando a proposta Uniforme:

f = function(x,mu,sigma){
  dcauchy(x,mu,sigma)  
}

# Usando uma q Uniforme U( mu-c*sigma , mu+c*sigma )
q = function(x,mu,sigma,c){
  ( 1/(2*c*sigma) )*I(x>(mu-c*sigma))*I(x<(mu+c*sigma))
}

mu = 0
sigma = 1
c = 8

razao = function(x,mu,sigma,c){
  f(x,mu,sigma)/q(x,mu,sigma,c)
}

# Obtendo k
k = optimize(razao,mu=mu,sigma=sigma,c=c,lower=-10,upper=10,maximum=TRUE)$objective
k

## [1] 5.092958

curve(f(x,mu=mu,sigma=sigma),-10,10,lwd=2,xlab=expression(x),ylab=expression(f[X](x)))
curve(k*q(x,mu=mu,sigma=sigma,c=c),col=2,lwd=2,lty=2,add=TRUE)

# Vamos usar a proposta Uniforme conforme descrito no exemplo:
x.gerado = NULL
x.aceito = NULL

M = 10000
j = 0
cont = 0  # apenas para sabermos quantos pontos precisaram ser gerados
while(j<M){
  cont = cont+1
  x.gerado[cont] = runif(1, mu-c*sigma , mu+c*sigma)
  u = runif(1,0,1)
  p.aceit = f(x.gerado[cont],mu=mu,sigma=sigma)/(k*q(x.gerado[cont],mu=mu,sigma=sigma,c=c))
  if(u<p.aceit){
    j=j+1
    x.aceito[j] = x.gerado[cont]}     
}

taxa.aceitacao = length(x.aceito)/length(x.gerado)
taxa.aceitacao

## [1] 0.1833651

hist(x.aceito,freq=FALSE,main="",xlab=expression(x),ylab=expression(f[X](x)),ylim=c(0,0.3)) 
curve(dcauchy(x,location=mu,scale=sigma),col=2,lwd=2,lty=2,add=TRUE)

Usando a proposta t-Student:

f = function(x,mu,sigma){
  dcauchy(x,mu,sigma)  
}

# Usando uma q t-Student t_nu(mu)
q = function(x,nu,mu){
  dt(x,nu,mu)
}

nu = 1/2
mu = 0
sigma = 1

razao = function(x,nu,mu,sigma){
  f(x,mu,sigma)/q(x,nu,mu)
}

# Obtendo k
k = optimize(razao,nu=nu,mu=mu,sigma=sigma,lower=-10,upper=10,maximum=TRUE)$objective
k

## [1] 1.345297

curve(f(x,mu=mu,sigma=sigma),-10,10,lwd=2,xlab=expression(x),ylab=expression(f[X](x)))
curve(k*q(x,nu=nu,mu=mu),col=2,lwd=2,lty=2,add=TRUE)

# Vamos usar a proposta Uniforme conforme descrito no exemplo:
x.gerado = NULL
x.aceito = NULL

M = 10000
j = 0
cont = 0  # apenas para sabermos quantos pontos precisaram ser gerados
while(j<M){
  cont = cont+1
  x.gerado[cont] = rt(1,df=nu,ncp=mu)
  u = runif(1,0,1)
  p.aceit = f(x.gerado[cont],mu=mu,sigma=sigma)/(k*q(x.gerado[cont],nu=nu,mu=mu))
  if(u<p.aceit){
    j=j+1
    x.aceito[j] = x.gerado[cont]}     
}

taxa.aceitacao = length(x.aceito)/length(x.gerado)
taxa.aceitacao

## [1] 0.7401925

hist(x.aceito[x.aceito>-10&x.aceito<10],freq=FALSE,main="",xlab=expression(x),ylab=expression(f[X](x))) 
curve(dcauchy(x,location=mu,scale=sigma),col=2,lwd=2,lty=2,add=TRUE)

Exercício 1

Implemente um gerador para:

\[ f_X(x)=3x^2, \qquad 0<x<1. \]

Estude empiricamente a taxa de aceitação para diferentes escolhas de distribuições auxiliares (prospostas).

Exercício 2

Crie uma função que simule valores de uma v.a. contínua \(X \sim \text{Beta}(a,b)\) através do Método da Aceitação-Rejeição, utilizando uma distribuição proposta uniforme padrão.
Compare o gerador implementado com a função nativa do R . Para isto utilize gráficos e medidas resumo sobre as amostras obtidas através das diferentes funções.

Referências

ROSS, Sheldon M. Simulation. 6. ed. London: Academic Press, 2022.

Método da Aceitação-Rejeição

Disciplina: GET00204 - Estatística Computacional 2

Motivação

Método da Aceitação-Rejeição para V.A. Contínuas

Intuição

Algoritmo

Proposição

Observação

Exemplo 1

Algoritmo 1.1

Algoritmo 1.2

Alguns erros comuns na implementação do Método de Aceitação-Rejeição:

Visitando algumas Distribuições Contínuas

Exemplo 2: Distribuição Beta

Algoritmo 2.1

Algoritmo 2.2

Exemplo 3: Distribuição Cauchy

Exercício 1

Exercício 2

Referências