Embed Size (px)
Citation preview
Adaptative Rejection Sampling - ARS
Adriana Lima
Silvio Cabral
18 de setembro de 2018
1 Rejection Sampling
A amostragem por rejeição (Rejection Sampling), baseia-se na amostragem de candidatos a partir de uma
distribuição mais fácil (g) e, em seguida, corrige-se a probabilidade de amostragem através da rejeição aleatória
de alguns candidatos.
Seja g a densidade da qual sabemos como amostrar e para a qual podemos facilmente calcular g(x). Seja
e(·) uma função (envelope), tal que e(x) = g(x)α ≥ f(x) para todo x para o qual f(x) > 0, para uma dada
constante α ≤ 1. A amostragem por rejeição é feita segundo o algoritmo:
Seguindo a regra de rejeição U > f(Y )/e(Y ), determinar se o valor Y será aceito ou não, é equivalente
a amostrar de U |Y = y ∼ U [0, e(y)] e veri�car se U |Y = y > f(y). Portanto, a regra de rejeição descarta
Y = y com probabilidade proporcional ao tamanho da barra acima de f(y), em relação ao seu comprimento
total.
Figura 1: Ilustração do Método de Amostragem por Rejeição para a Função Alvo f Usando o Envelope e.
1
2 Squeezed Rejection Sampling
A amostragem por rejeição requer que f seja avaliada para cada Y gerado. Porém, o custo computacional
para avaliar f pode ser alto, neste caso o uso amostragem por rejeição comprimida pode ser atrativo.
Uma estratégia para impedir a avaliação de f em alguns casos, é empregar uma função de compressão s,
tal que 0 ≤ s(x) ≤ f(x) em todo o suporte de f . Além disso, um envelope, e(·), também é utilizado.
A amostragem por rejeição comprimida (Squeezed Rejection Sampling) é realizada segundo o algoritmo:
Quando Y = y, o valor é aceito com probabilidade f(y)/e(y) e rejeitado com probabilidade 1−f(y)/e(y),
as mesmas probabilidades da amostragem por rejeição comum. Porém, quando veri�camos a condição de
U ≤ s(Y )/e(Y ), permitimos que aceitação de Y seja realizada avaliando s ao invés de f . Como s(x) < f(x),
em todo o suporte de f , obtemos a diminuição no número de avaliações de f .
Figura 2: Ilustração do Método de Amostragem por Rejeição Comprimida para a Função Alvo, f , Usando oEnvelope e e a Função de Compressão s. Keep Fist e Keep Later correspondem a linha 4 e 8 do algoritmo,respectivamente.
3 Adaptive Rejection Sampling
Uma das principais di�culdades da amostragem por rejeição é a construção de um envelope adequado.
Gilks e Wild propuseram uma estratégia de geração automática de envelope para amostragem de rejeição
comprimida, considerando que densidade é log-côncava, contínua e diferenciável. Onde as funções de envelope
2
(e) e compressão (s) são re�nada de forma iterativa.
Seja `(x) = log f(x). Sobre f , assumidos que a função é estritamente positiva, log-concava, contínua e
diferenciável.
O algoritmo é inicializado avaliando ` e `′ em k pontos, x1 < · · · < xk. Se o suporte de f se estende para
−∞, x1 deve ser escolhido tal que `′(x1) > 0. Caso o suporte de estenda até +∞, xk deve ser escolhido tal
que `′(xk) < 0.
De�na o casco superior de ` (denotado por e∗k), formado pela tangente de ` em cada ponto em Tk. As
tangentes em xi e xi+1 se intersectam no ponto
zi =`(xi+1)− `(xi)− xi+1`
′(xi+1) + xi`′(xi)
`′(xi)− `′(xi+1)(1)
para i = 1, . . . , k − 1.
Então,
e∗k = `(xi) + (x− xi)`′(xi) x ∈ [zi−1, zi] (2)
para i = 1, . . . , k, com z0 e zk de�nidos, respectivamente, como mínimo e máximo do suporte de f .
De�na a função de rejeição em Tk, como sendo ek(x) = exp {e∗k(x)}De�na o casco inferior de ` (denotado por s∗k), formado pelas retas que ligam os pontos (xi, `(xi)) e
(xi+1, `(xi+1)). Esse casco inferior é dado por
s∗k(x) =(xi+1 − x)`(xi) + (x− xi)`(xi+1)
xi+1 − xi(3)
para x ∈ [xi, xi+1]. Quando x < x1 ou x > xk, faça s∗k = −∞
Figura 3: Casco Inferior e Superior para `′(x) = logf(x) usado no Método de Amostragem por RejeiçãoAdaptativa quando k = 5.
Então, a função squeezing é dada por sk(x) = exp {s∗k(x)}A amostragem por rejeição adaptativa (ou Adaptive Rejection Sampling) é realizada segundo o algoritmo:
3
Gilks propôs um método similar para a geração automática do envelope, que não requer o cálculo de `′.
Dado o conjunto de pontos Tk, de�na a função Li(·) como sendo a reta que conecta os pontos (xi, `(xi)) e
(xi+1, `(xi+1)), para i = 1, . . . , k − 1. De�na
e∗k(x) =
min {Li−1(x), Li(x)} se x ∈ [xi, xi+1],
L1(x) se x < x1,
Lk−1(x) se x > xk,
(4)
por conversão, L0(x) = Lk(x) = +∞ A função de rejeição em Tk é dada por e(x) = exp {e∗k(x)} Como a
função squeezing de�nida em 3 não depende de `′, ela se mantém a mesma.
Figura 4: Casco Inferior e Superior para `′(x) = logf(x) usado no Método de Amostragem por RejeiçãoAdaptativa Livre de Derivada (não requer o cálculo de `′) quando k = 5.
O envelope proposto não é tão e�ciente quanto o quando utilizamos `′
4
Figura 5: Envelope e Função de Compressão para o Método de Amostragem por Rejeição Adaptativo.A Densidade Alvo é a Curva em Forma de Sino. O Primeiro Método Apresentado no Texto, usando `′,Produz o Envelope Representado pela Região Cinza Claro. O Método Livre de Derivada Produz o EnvelopeRepresentado pela Região Cinza Escuro. A Função de Compressão para Ambos Abordagens é Dada pelaCurva Pontilhada.
3.1 Exemplo
Suponha que desejamos simular de uma distribuição Beta(α = 2, β = 2). Sabemos que a distribuição
Beta(α = 2, β = 2) tem f.d.p. dada por:
f(x) = 6x(1− x) x ∈ (0, 1) e α, β > 0
0,0 0,2 0,4 0,6 0,8 1,0
0,0
0,5
1,0
1,5
x
f(x)
Figura 6: Densidade da Distribuição Beta(α = 2, β = 2)
A função de log-densidade é dada por:
`(x) = log(6) + log(x) + log(1− x).
Pode-se perceber através da Figura 8, a distribuição Beta(α = 2, β = 2) é log-côncava.
5
0,0 0,2 0,4 0,6 0,8 1,0
−2,
5−
2,0
−1,
5−
1,0
−0,
50,
00,
5
xl(x
)
Figura 7: Logaritmo da Densidade da Distribuição Beta(α = 2, β = 2)
A derivada da log-densidade é:
`′(x) =1
x− 1
1− xOs chutes iniciais escolhidos serão T3 = (0, 2; 0, 4; 0, 7).
O envelope para a log-densidade (e∗k), segundo o primeiro método proposto (Equação 2 ), é dado pelas
equações:
e∗1(x) = `(x1) + (x− x1)`′(x1) = −0, 0408 + (x− 0, 2)3, 75 se x ∈ [z0, z1]
e∗2(x) = `(x2) + (x− x2)`′(x2) = 0, 3646 + (x− 0, 4)0, 8333 se x ∈ [z1, z2]
e∗3(x) = `(x3) + (x− x3)`′(x3) = 0, 2311 + (x− 0, 7)(−1, 9048) se x ∈ [z2, z3]
onde,
z0 = mínimo do suporte def
z1 =`(x2)− `(x1)− x2`′(x2) + x1`
′(x1)
`′(x1)− `′(x2)= 0, 2819
z2 =`(x3)− `(x2)− x3`′(x3) + x2`
′(x2)
`′(x2)− `′(x3)= 0, 5599
z3 = máximo do suporte def
6
0,0 0,2 0,4 0,6 0,8 1,0
−2,
5−
2,0
−1,
5−
1,0
−0,
50,
00,
5
xl(x
)
Figura 8: Envelope (e∗∗k) e Função de Compressão (s∗∗k) para o Logaritmo da Densidade da DistribuiçãoBeta(α = 2, β = 2), segundo o ARS. A Log-Densidade é a Curva em Preto. A Função de Compressão é Dadapela Curva Azul, e o Envelope pela Curva em Vermelho.
Voltando para a função alvo, temos:
0,0 0,2 0,4 0,6 0,8 1,0
0,0
0,5
1,0
1,5
2,0
x
f(x)
Figura 9: Envelope (e) e Função de Compressão (s) para a Densidade da Distribuição Beta(α = 2, β = 2),segundo o ARS. A Densidade é a Curva em Preto. A Função de Compressão é Dada pela Curva Azul, e oEnvelope pela Curva em Vermelho.
Suponha agora que ao gerar de g(.), tenhamos obtido o valor 0, 8. Este valor não será aceito na primeira
condição, referente a função de compressão (linha 5 do algoritmo). Porém, suponha que o valor gerado
de U atende a condição do segundo critério de aceitação (linha 9 do algoritmo), então temos que T4 =
(0, 2; 0, 4; 0, 7; 0, 8). O evenvelope e a função de compressão atualizados, são mostrados na Figura 10 a seguir:
7
0,0 0,2 0,4 0,6 0,8 1,0
−2,
5−
2,0
−1,
5−
1,0
−0,
50,
00,
5
x
l(x)
(a) Grá�co da Log-Densidade
0,0 0,2 0,4 0,6 0,8 1,0
0,0
0,5
1,0
1,5
2,0
x
f(x)
(b) Grá�co da Densidade
Figura 10: Envelope e Função de Compressão para a Distribuição Beta(α = 2, β = 2), segundo o ARS. ALog-Densidade ou Densidade é a Curva em Preto. A Função de Compressão é Dada pela Curva Azul, e oEnvelope pela Curva em Vermelho.
Para a geração de uma amostra aleatória o algoritmo deveria ser executado até que a amostra atingisse
o tamanho desejado. Este exemplo, porém, é somente para ilustrar como se dá a mudança do envelope e da
função de compressão ao longo das iterações.
4 Utilizando o Adaptive Rejection Sampling no R
Utilizaremos o pacote AdapSamp.
install.packages("AdapSamp")
require(AdapSamp)
Este pacote além de gerar a.a. utilizando o Adaptive Rejection Sampling, também é capaz de executar o
Modi�ed Adaptive Rejection Sampling Algorithm, Concave-convex Adaptive Rejection Sampling Algorithm
e Adaptive Slice Sampling Algorithm. Sabe-se que o ARS só pode ser utilizado quando a função alvo
respeita algumas características, então, estes outros métodos são propostos como modi�cações do ARS onde
as suposições são mais �exíveis.
Para amostrar utilizando o ARS utiliza-se a função rARS(.), onde a sintaxe da função é:
rARS(n, formula,min = −Inf,max = Inf, sp),
onde n é o tamanho da amostra que se quer gerar; formula é a função de densidade; min,max é o
domínio da função de densidade, e, sp são os chutes iniciais (Tk).
O chute inicial não pode ser dado de forma aleatória, ele tem que obedecer algumas condições. Viu-se na
Seção 3 que x1 deve ser escolhido tal que `′(x1) > 0 e xk deve ser escolhido tal que `′(xk) < 0. Isto é, sendo
xMV o ponto de máximo, temos que x1 < xMV e xK > xMV . Dentre os xi(i = 1, ..., k), também não pode
estar o ponto de máximo.
8
5 Exemplos Computacionais
Apresentaremos três exemplos, que consistem em gerar amostras aleatórias com distribuições de densidade
Beta (Seção 5.1 ), Gama (Seção 5.2 ) e Normal (Seção 5.3 ). Em todos os exemplos vamos gerar uma amostra
de tamanho de n = 100 utilizando os algoritmos RS e ARS, além da função pronta no R. Compararemos
o tempo computacional gasto utilizando o pacote microbenchmark, além do quanto as médias e variâncias
amostrais se aproximaram da teórica. Foi �xado a semente set.seed(666).
require(microbenchmark)
n = 100
set.seed(666)
5.1 Beta(2,2)
Suponha que desejamos simular de uma distribuição Beta(α = 2, β = 2). Sabemos que a distribuição
Beta(α = 2, β = 2) tem f.d.p. dada por:
f(x) = 6x(1− x) x ∈ (0, 1) e α, β > 0,
com E(X) ==α
α+ β= 0, 5 e Var(X) = αβ
(α+β+1)(α+β)2 = 0, 05.
5.1.1 Geração dos Dados
� Rejection Sampling
Utilizando g(x) com distribuição Unif(0, 1) e α = 2/3, então o critério de aceitação é U ≤ 4y(1− y).
RS = function( n ){
k = 0
j = 0
x = NULL
while (k < n) {
u = runif(1)
j = j + 1
y = runif(1)
if ( u <= 4*y*(1-y) ) {
k = k + 1
x[k] = y
}
}
return(x)
}
9
amostra = RS(n)
� Função implementada no R
amostra1 = rbeta( n, 2, 2 )
� ARS utilizando somente o núcleo da função alvo, q(.)
amostra2 = rARS( n, "x*(1-x)", 0, 1, c(0.4,0.6) )
� ARS utilizando a função alvo, f(.)
amostra3 = rARS( n, "6*x*(1-x)", 0, 1, c(0.4,0.6) )
5.1.2 Comparando os Métodos
� Diferença entre a Média e Variância amostral e teórica
media = abs( cbind(mean(amostra),mean(amostra1),mean(amostra2),mean(amostra3)) - 0.5 )
variancia = abs( cbind(var(amostra),var(amostra1),var(amostra2),var(amostra3)) - 0.05 )
erro = round( rbind(media,variancia) , 4 )
colnames(erro) = c("RS", "R" , "ARS q(.)" , "ARS f(.)")
row.names(erro) = c("media", "variancia")
erro
## RS R ARS q(.) ARS f(.)
## media 0,0311 0,0083 0,017 0,0111
## variancia 0,0065 0,0019 0,006 0,0023
� Tempo Computacional
require(microbenchmark)
microbenchmark(
RS(n),
rbeta( n, 2, 2 ),
rARS( n, "x*(1-x)", 0, 1, c(0.4,0.6) ),
rARS( n, "6*x*(1-x)", 0, 1, c(0.4,0.6) ),
times = 100)
## Unit: microseconds
## expr min lq mean median
## RS(n) 3002 3339 5513 3487
10
## rbeta(n, 2, 2) 90 106 108 109
## rARS(n, "x*(1-x)", 0, 1, c(0.4, 0.6)) 338020 393906 448784 431552
## rARS(n, "6*x*(1-x)", 0, 1, c(0.4, 0.6)) 343301 400114 459812 432281
## uq max neval cld
## 3824 135783 100 a
## 113 140 100 a
## 493084 749263 100 b
## 500298 897200 100 b
5.2 Gama(2,1)
Suponha que desejamos simular de uma distribuição Gama(α = 2, β = 1). Sabemos que a distribuição
Gama(α = 2, β = 1) tem f.d.p. dada por:
f(x) = xe−x x ∈ [0,∞) e α, β > 0,
com E(X) = αβ = 2 e Var(X) = α
β2 = 2.
5.2.1 Geração dos Dados
� Rejection Sampling
Utilizando g(x) com distribuição Exp(λ = 1/2) e α = e/4, então o critério de aceitação é U ≤x2 exp
{− 1
2x+ 1}.
RS = function( n ){
k = 0
j = 0
x = NULL
while (k < n) {
u = runif(1)
j = j + 1
y = rexp(1 , 1/2 )
if ( u <= (y/2) * exp( -0.5*y + 1 )) {
k = k + 1
x[k] = y
}
}
return(x)
}
11
amostra = RS(n)
� Função implementada no R
amostra1 = rgamma( n, 2, 1 )
� ARS utilizando f(.)
amostra2 = rARS(n, "x*exp(-x)", 0, Inf, c(0.001,5) )
5.2.2 Comparando os Métodos
� Diferença entre a Média e Variância amostral e teórica
media = abs( cbind(mean(amostra),mean(amostra1),mean(amostra2)) - 2 )
variancia = abs( cbind(var(amostra),var(amostra1),var(amostra2)) - 2 )
erro = round( rbind(media,variancia) , 4 )
colnames(erro) = c("RS", "R" , "ARS f(.)")
row.names(erro) = c("media", "variancia")
erro
## RS R ARS f(.)
## media 0,22 0,25 0,024
## variancia 0,55 0,74 0,249
� Tempo Computacional
microbenchmark(
RS(n),
rgamma( n, 2, 1 ),
rARS(n, "x*exp(-x)", 0, Inf, c(0.001,5) ),
times = 100)
## Unit: microseconds
## expr min lq mean median
## RS(n) 3477 3771 5639 3904
## rgamma(n, 2, 1) 55 86 110 93
## rARS(n, "x*exp(-x)", 0, Inf, c(0.001, 5)) 666428 733282 902222 807052
## uq max neval cld
## 4224 77220 100 a
## 101 1728 100 a
## 955474 1981511 100 b
12
5.3 Normal(0,1)
Suponha que desejamos simular de uma distribuição Normal(µ = 0, σ2 = 1). Sabemos que a distribui-
ção Normal(µ = 0, σ2 = 1) tem f.d.p. dada por:
f(x) =1√2π
exp
[−x
2
2
], x ∈ (−∞,∞).
com E(X) = µ = 0 e Var(X) = σ2 = 1.
5.3.1 Geração dos Dados
� Rejection Sampling
Utilizando g(x) com distribuição Exp(1), então o critério de aceitação é U ≤ exp(− (y−1)2
2
).
RS = function( n ){
k = 0
j = 0
x = NULL
while (k < n) {
u = runif(1)
j = j + 1
y = rexp(1)
if ( u <= exp(-(y-1)^2/2) ) {
k = k + 1
u2 = runif(1)
if(u2 <= 0.5){
x[k] = y
}else{
x[k] = - y
}
}
}
return(x)
}
amostra = RS(n)
� Função implementada no R
13
amostra1 = rnorm( n )
� ARS utilizando q(.)
amostra2 = rARS(n,"exp(-(x)^2/2)",-Inf,Inf,c(-4,4))
� ARS utilizando f(.)
amostra3 = rARS(n,"exp(-(x)^2/2)/sqrt(2*pi)",-Inf,Inf,c(-4,4))
5.3.2 Comparando os Métodos
� Diferença entre a Média e Variância amostral e teórica
media = abs( cbind(mean(amostra),mean(amostra1),mean(amostra2),mean(amostra3)) - 0)
variancia = abs( cbind(var(amostra),var(amostra1),var(amostra2),var(amostra3)) - 1 )
erro = round( rbind(media,variancia) , 4 )
colnames(erro) = c("RS", "R" , "ARS q(.)" , "ARS f(.)")
row.names(erro) = c("media", "variancia")
erro
## RS R ARS q(.) ARS f(.)
## media 0,15 0,125 0,0031 0,018
## variancia 0,24 0,026 0,1497 0,037
� Tempo Computacional
microbenchmark(
RS(n),
rnorm( n ),
rARS(n,"exp(-(x)^2/2)",-Inf,Inf,c(-4,4)),
rARS(n,"exp(-(x)^2/2)/sqrt(2*pi)",-Inf,Inf,c(-4,4)),
times = 100)
## Unit: microseconds
## expr min lq
## RS(n) 3973 4304
## rnorm(n) 42 53
## rARS(n, "exp(-(x)^2/2)", -Inf, Inf, c(-4, 4)) 1028357 1141087
## rARS(n, "exp(-(x)^2/2)/sqrt(2*pi)", -Inf, Inf, c(-4, 4)) 1036102 1176834
## mean median uq max neval cld
## 6494 4460 4806 59143 100 a
## 59 60 63 205 100 a
## 1415153 1229886 1493530 3493872 100 b
## 1390250 1256117 1421721 3144997 100 b
14
6 Conclusão
O ARS é tão bom quanto os demais métodos, porém o tempo computacional gasto para a geração de amostras
aleatórias utilizando este algoritmo é muito maior. Isto acontece pois no ARS um novo envelope e uma nova
função de compressão tem que ser calculados toda vez que um novo valor é aceito na amostra.
A preferência pelo método ARS deve se dar quando o cálculo analítico do envelope, para a utilização
do Reject Sampling, é muito custoso. Isto é, quando achar uma função que exceda a função alvo em todo
domínio for muito complexo, é preferível optar pela utilização do ARS.
Referências
[1] Geof H. Givens, Jennifer A. Hoeting (2a Edição, 2005). Computational Statistics.WileyInterscience
[2] Dong Zhang (2018). AdapSamp: Adaptive Sampling Algorithms. R package version 1.1.1. https://CRAN.R-
project.org/package=AdapSamp
[3] Olaf Mersmann (2015). microbenchmark: Accurate Timing Functions. R package version 1.4-2.1.
https://CRAN.R-project.org/package=microbenchmark
15
Exercício Proposto
1. Dada a densidade
f(x|α, β) =Γ(α+ β)
Γ(α)Γ(β)xα−1(1− x)β−1I[0,1](x) (5)
(a) Considerando os pontos iniciais: P=(0.21, 0.64, 0.81), faça o grá�co da log-densidade, considerando
α = 2 e β = 5, com o casco inferior (função squeezing) e superior (envelope), segundo o ARS que
utiliza a derivada e o que não utiliza.
(b) No item anterior, o que acontece com as funções squeezing e o envelope se o valor 0.5 for aceito?
(c) Use o ARS para gerar uma amostra aleatória de tamanho 10000, considerando α = 2 e β = 5.
Compare os quartis da amostra gerada com os quantis teóricos.
16
