Método da Aceitação-Rejeição
Disciplina: GET00204 - Estatística Computacional 2
Motivação
Suponha que queremos gerar valores de uma v.a. discreta \(X\), cuja função de probabilidade é dada por \[p_X(x_j) = p_j, \qquad j = 0,1,\ldots.\]
Na impossibilidade de gerar os valores diretamente de \(p_X\), o método da aceitação-rejeição é um método de reamostragem que propõe gerar um valor \(y\) de uma função de probabilidade auxiliar \(q\) e aceitar (ou não) o valor gerado com probabilidade proporcional à razão \(\dfrac{p_X(y)}{q(y)}\).
Método da Aceitação-Rejeição para V.A. Discretas
Suponha \(X\) uma v.a. discreta com f.p.
\[ P(X = x_j) = p_X(x_j) = p_j, \qquad j = 0,1,\ldots, \qquad \mbox{com} \qquad \sum_{j} p_j = 1. \]
Neste caso, se \(Y\) é uma v.a. discreta com f.p. \(P(Y = y) = q(y)\) da qual sabemos simular e, se existe uma constante \(k\) tal que
\[ \dfrac{p_X(x_j)}{q(x_j)} \leq k, \qquad \forall x_j \in Im(X), \]
então, aceitando \(x_j\) com probabilidade \(\dfrac{p_X(x_j)}{k\cdot q(x_j)}\), obtemos um valor representativo da função de probabilidade \(p_X\).
Algoritmo
Gere um valor \(y\) da v.a. \(Y\) com função de probabilidade \(q\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{p_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 probabilidade \(p_X\).
Demonstração:
Seja \(X\) uma v.a. discreta com função de probabilidade \(p_X\) e \(y\) um valor simulado de outra v.a. discreta \(Y\), com \(Im(X) \subset Im(Y)\), através da função de probabilidade \(q\). Neste caso, e considere o evento \(A=\{\mbox{ponto }y\mbox{ aceito}\}\) \[\begin{align} P(Y=y,A) = P(Y=y)P(A|Y=y) = q(y) \dfrac{p_X(y)}{k \cdot q(y)} = \dfrac{p_X(y)}{k}. \end{align}\]
Além disso, \[\begin{align} P(A) & = \sum_yP(Y=y,A) \\ & = \sum_y P(Y=y)P(A|Y=y) \\ & = \sum_y q(y) \dfrac{p_X(y)}{k \cdot q(y)} = \sum_y \dfrac{p_X(y)}{k} = \dfrac{1}{k}; \end{align}\] portanto, \[\begin{align} P(X=y|A) = \dfrac{P(Y=y,A)}{P(A)} = \dfrac{\dfrac{p_X(y)}{k}}{\dfrac{1}{k}} = p_X(y), \end{align}\] ou seja, os valores aceitos seguem o modelo com f.p. dada por \(p_X \quad\quad \square\).
Exemplo 1
Simular valores para a v.a. \(X\) com f.p. dada por
| \(x\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(p_X(x)\) | 0,11 | 0,12 | 0,09 | 0,08 | 0,12 | 0,10 | 0,09 | 0,09 | 0,10 | 0,10 |
Solução:
Inicialmente, note que, tomando \[ \displaystyle k = max_j\dfrac{p_X(j)}{q(j)} = max_j\dfrac{p_X(j)}{1/10} = 10 \cdot max_j p_X(j) = 10\cdot 0,12 = 1,2, \quad \mbox{ temos } \quad \dfrac{p_X(j)}{q(j)} \leq k, \forall j \in Im(X). \]
Algoritmo 1.1
Gere um valor \(y\) de \(Y \sim \text{U}\{1,2,...,10\}\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{p_X(y)}{1,2 \cdot 1/10} = \dfrac{p_X(y)}{0,12}\), tome \(x = y\); caso contrário, retorne ao passo 1.
pX = c(0.11, 0.12, 0.09, 0.08, 0.12, 0.10, 0.09, 0.09, 0.10, 0.10)
n = 10000
y = sample(1:10,n,replace = TRUE)
x = NULL
for (i in 1:n){
u = runif(1)
x[i] = ifelse(u < pX[y[i]]/0.12,y[i],NA)
}
# número de pontos aceitos
n.gerados = length(x[is.na(x)==FALSE])
n.gerados## [1] 8328table(x)/n.gerados## x
## 1 2 3 4 5 6 7
## 0.11107109 0.12584054 0.09834294 0.07396734 0.11551393 0.10218540 0.08909702
## 8 9 10
## 0.08729587 0.10302594 0.09365994Algoritmo 1.2 (supondo que não conhecemos a função sample)
Gere um valor \(u_1\) de \(U \sim \text{U}(0,1)\) e tome \(y = \left\lfloor(10\cdot u_1)\right\rfloor + 1\).
Gere um segundo valor \(u_2\) de \(U \sim \text{U}(0,1)\).
Se \(u_2 < \dfrac{p_X(y)}{1,2 \cdot 1/10} = \dfrac{p_X(y)}{0,12}\), tome \(x = y\); caso contrário, retorne ao passo 1.
pX = c(0.11, 0.12, 0.09, 0.08, 0.12, 0.10, 0.09, 0.09, 0.10, 0.10)
n = 10000
u1 = runif(n)
y = floor(10*u1)+1
x = NULL
for (i in 1:n){
u2 = runif(1)
x[i] = ifelse(u2 < pX[y[i]]/0.12,y[i],NA)
}
# número de pontos aceitos
n.gerados = length(x[is.na(x)==FALSE])
n.gerados## [1] 8314table(x)/n.gerados## x
## 1 2 3 4 5 6 7
## 0.11005533 0.12448881 0.08816454 0.07962473 0.12244407 0.09718547 0.09309598
## 8 9 10
## 0.09153236 0.09947077 0.09393794Observação
Como os valores da variável auxiliar \(Y\) são gerados de forma independente, tendo, cada um, probabilidade de aceitação igual a \(1/k\), a v.a. que conta o número de valores gerados até que o primeiro \(y\) seja aceito como representante de \(p_X\) segue uma distribuição Geométrica com parâmetro \(1/k\), cuja média é igual a \(k\).
pX = c(0.11, 0.12, 0.09, 0.08, 0.12, 0.10, 0.09, 0.09, 0.10, 0.10)
k = 5
n = 10000
x = NULL
N = NULL
for (j in 1:n){
cont = 0
aceitou = FALSE
while(!aceitou){
# gera uma nova proposta Y ~ Uniforme{1,...,10}
y = sample(1:10, 1)
# gera U ~ Uniforme(0,1)
u = runif(1)
# conta esta proposta
cont = cont + 1
# aceita ou rejeita a proposta
if (u <= pX[y] / (k/10)){
x[j] = y
aceitou = TRUE
}
}
# número de propostas necessárias para obter x[j]
N[j] = cont
}
# número de pontos aceitos
n.gerados = length(x[is.na(x)==FALSE])
n.gerados## [1] 10000table(x)/n.gerados## x
## 1 2 3 4 5 6 7 8 9 10
## 0.1127 0.1208 0.0920 0.0828 0.1160 0.1014 0.0855 0.0864 0.1009 0.1015k## [1] 5mean(N)## [1] 5.0217Visitando algumas Distribuições Discretas
Exemplo 2: Ditribuição Bernoulli
Simular valores de uma v.a. \(X \sim \text{Ber}(p)\), para \(p\) fixo.
Solução:
Sejam \(X\) e \(U\) duas v.a. discretas tais que \(X \sim \text{Ber}(p)\) e \(Y \sim \text{U}\{0,1\}\).
Neste caso, podemos tomar \[ \displaystyle k = max_j\dfrac{p_X(j)}{q(j)} = max_j\dfrac{p^j(1-p)^{1-j}}{1/2} = 2 \cdot max \{p,1-p\}. \]
Algoritmo 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{p_X(y)}{max(p,1-p)}\), tome \(x = y\); caso contrário, retorne ao passo 1.
n = 10000
p = 0.3
rbern.rej = function(n,p){
k = max(2*p, 2*(1-p))
x = rep(NA,n)
cont = 0
while(cont < n){
# proposta uniforme
y = sample(c(0,1),1,replace = TRUE)
u = runif(1)
px = (p^y)*((1-p)^(1-y))
gy = 0.5
if(u < px/(k*gy)){
x[cont] = y
cont = cont+1
}
}
x
}
x = rbern.rej(n,p)
table(x)/n## x
## 0 1
## 0.7015 0.2984Exemplo 3: Distribuição Geométrica
Seja \(X\) a v.a. que conta o número de realizações independentes de um ensaio de Bernoulli com probabilidade de sucesso \(p\) necessárias até que ocorra o primeiro sucesso.
Simular valores da v.a. \(X\), para \(p\) fixo.
Solução:
Neste caso, vamos denotar \(X \sim \text{Geo}(p)\), cuja função de probabilidade é dada por \[p_X(x) = p(1-p)^{x-1} = p q^{x-1}, \quad x = 1, 2, 3, \dots,\] onde denotamos \(q=1-p\).
Note que a imagem da v.a. de interesse é infinita, portanto utilizar a distribuição uniforme discreta como distribuição auxiliar gera um problema. Por simplicidade, faremos esta aproximação supondo um conjunto discreto \(\{1,...,M\}\) para \(M\) grande apenas para ilustrar o comportamento do algoritmo.
A função de probabilidade do modelo Geométrico é decrescente em \(x\); logo, utilizando a distribuição uniforme discreta no conjunto \(\{1,...,M\}\) como distribuição auxiliar podemos tomar \[ \displaystyle k = max_j\dfrac{p_X(j)}{q(j)} = \dfrac{p(1-p)^{1-1}}{1/M} = M\cdot p. \]
Algoritmo 3
Gere um valor \(y\) de \(Y \sim \text{U}\{1,...,M\}\).
Gere um valor \(u\) de \(U \sim \text{U}(0,1)\).
Se \(u < \dfrac{p_X(y)}{M\cdot p\cdot (1/M)} = (1-p)^{y-1}\), tome \(x = y\); caso contrário, retorne ao passo 1.
n = 10000
p = 0.3
rgeom.rej = function(n,p){
M = 100
k = M*p
x = rep(NA,n)
cont = 0
while(cont < n){
# proposta uniforme
y = sample(1:M,1,replace = TRUE)
u = runif(1)
px = p*((1-p)^(y-1))
gy = 1/M
if(u < px/(k*gy)){
x[cont] = y
cont = cont+1
}
}
x
}
x = rgeom.rej(n,p)
freq.observadas = prop.table(table(x))
valores.gerados = as.numeric(names(freq.observadas))
freq.teoricas = dgeom(valores.gerados-1,p)
bp = barplot(freq.observadas,
main = "Amostragem vs. Distribuição Teórica",
xlab = "Valores (x)",
ylab = "Probabilidade",
col = "lightblue",
ylim = c(0, max(freq.observadas) * 1.15))
points(bp,freq.teoricas,type="p",col=2,pch=19)
Exercício 1
Crie uma função que simule valores de uma v.a. discreta \(X \sim \text{Poi}(\lambda)\) através do método da aceitação-rejeição.
Exercício 2
Crie uma função que simule valores de uma v.a. discreta \(X \sim \text{Bin}(n,p)\) através do método da aceitação-rejeição.
Exercício 3
Compare os geradores implementados com as funções nativas 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.