Rao-Blackwellization de Gibbs Sampler

9

Atualmente, estou estimando um modelo de volatilidade estocástica com os métodos Monte Carlo da Cadeia de Markov. Assim, estou implementando os métodos de amostragem de Gibbs e Metropolis.

Supondo que eu considere a média da distribuição posterior, e não uma amostra aleatória, é isso que é comumente chamado de Rao-Blackwellization ?

No geral, isso resultaria em tomar a média sobre as médias das distribuições posteriores como estimativa de parâmetros.

mscnvrsy
fonte

Respostas:

10

Supondo que eu considere a média da distribuição posterior, e não uma amostra aleatória, é isso que é comumente chamado de Rao-Blackwellization?

Não estou muito familiarizado com os modelos de volatilidade estocástica, mas sei que, na maioria das situações, a razão pela qual escolhemos os algoritmos Gibbs ou MH para desenhar a partir do posterior é porque não conhecemos o posterior. Muitas vezes, queremos estimar a média posterior e, como não sabemos a média posterior, extraímos amostras da posterior e a estimamos usando a média da amostra. Portanto, não tenho certeza de como você conseguirá calcular a média da distribuição posterior.

Em vez disso, o estimador Rao-Blackwellized depende do conhecimento da média do condicional completo; mas mesmo assim a amostragem ainda é necessária. Eu explico mais abaixo.

Suponha que a distribuição posterior seja definida em duas variáveis, ), de modo que você queira estimar a média posterior: . Agora, se um amostrador Gibbs estivesse disponível, você poderia executá-lo ou executar um algoritmo MH para amostrar a partir do posterior.E [ θ dados ]θ=(μ,ϕE[θdata]

Se você pode executar um amostrador Gibbs, conhece em formato fechado e conhece a média dessa distribuição. Que isso signifique ser . Observe que é uma função de e dos dados.ϕ ϕ μf(ϕμ,data)ϕϕμ

Isso também significa que você pode integrar out partir do posterior, de modo que o marginal posterior de seja (isso não é conhecido completamente, mas é conhecido como constante). Agora você deseja executar uma cadeia de Markov de forma que seja a distribuição invariante e obtenha amostras desse marginal posterior. A questão éu f ( u | d um t um ) f ( u | d um t um )ϕμf(μdata)f(μdata)

Como você pode agora estimar a média posterior de usando apenas essas amostras da margem posterior de ?μϕμ

Isso é feito via Rao-Blackwellization.

E[ϕdata]=ϕf(μ,ϕdata)dμdϕ=ϕf(ϕμ,data)f(μdata)dμdϕ=ϕf(μdata)dμ.

Assim, suponha que tenhamos obtido as amostras partir da margem posterior de . Então u & Phi; = 1X1,X2,XNμ

ϕ^=1Ni=1Nϕ(Xi),

é chamado de estimador Rao-Blackwellized para . O mesmo pode ser feito simulando também os marginais comuns.ϕ

Exemplo (Puramente para demonstração).

Suponha que você tenha uma articulação desconhecida posterior para da qual você deseja provar. Seus dados são alguns e você tem os seguintes condicionais completos y μ | φ , y ~ N ( φ 2 + 2 y , y 2 ) φ | μ , y ~ L um m m um ( 2 μ + y , y + 1 )θ=(μ,ϕ)y

μϕ,yN(ϕ2+2y,y2)
ϕμ,yGamma(2μ+y,y+1)

Você executa o amostrador de Gibbs usando esses condicionais e obtém amostras da articulação posterior . Que essas amostras sejam . Você pode encontrar a média amostral dos s, e esse seria o estimador de Monte Carlo usual para a média posterior para .( μ 1 , & Phi 1 ) , ( μ 2 , φ 2 ) , ... , ( μ N , φ N ) φ φf(μ,ϕy)(μ1,ϕ1),(μ2,ϕ2),,(μN,ϕN)ϕϕ

Ou, observe que pelas propriedades da distribuição Gamma

E[ϕ|μ,y]=2μ+yy+1=ϕ.

Aqui são os dados fornecidos a você e, portanto, são conhecidos. O estimador Rao Blackwellized seria entãoy

ϕ^=1Ni=1N2μi+yy+1.

Observe como o estimador para a média posterior de nem usa as amostras e apenas usa as amostras . De qualquer forma, como você pode ver, ainda está usando as amostras obtidas de uma cadeia de Markov. Este não é um processo determinístico.ϕϕμ

Greenparker
fonte
Então, assumindo que a distribuição posterior do parâmetro seja conhecida (o que, de acordo com o meu conhecimento, é verdadeiro ao aplicar a amostragem de Gibbs), calcular a média da distribuição em vez de uma amostra aleatória seria o estimador Rao-Blackwellized? Espero ter entendido sua resposta corretamente. Muito obrigado já!
22416 mscnvrsy
Isso está incorreto. Na amostragem de Gibbs, você não conhece a distribuição posterior do parâmetro, mas conhece o posterior condicional completo para cada parâmetro. Existe uma grande diferença entre os dois. Acima, o posterior é que é desconhecido, e para o amostrador de Gibbs para o trabalho que você precisa saber tanto e . E você também está incorreto em seu segundo entendimento. Você ainda precisa colher uma amostra da margem posterior de e, em seguida, calcular a média da amostra de usando essas amostras para encontrar o estimador RB. f(μ,ϕdata)f(μϕ,data)f(ϕμ,data)μϕ
Greenparker
@mscnvrsy eu adicionei um exemplo para ajuda
Greenparker
Uau, muito obrigado por esclarecer isso para mim. Então, assumindo que conheço todas as distribuições condicionais, posso trabalhar com os meios teóricos das distribuições condicionais e calcular a média desses meios teóricos (como E [phi | mu, y]) para obter o estimador RB? Isso minimizaria a variação das minhas estimativas de parâmetros?
mscnvrsy
Se você estivesse obtendo amostras independentes, sim, minimizaria a variação dos estimadores, no entanto, como você está lidando com cadeias de Markov, é geralmente conhecido que o RB não necessariamente reduz a variação e há alguns casos em que a variação aumenta. Este artigo de Charlie Geyer deu alguns exemplos a este ponto.
Greenparker
9

O amostrador Gibbs pode então ser usado para melhorar a eficiência de (digamos) amostras de um posterior marginal, chame-o . Nota Assim, o a densidade marginal de com algum valor é o valor esperado da densidade condicional de fornecida no ponto .π2(θ2|y)

π2(θ2|y)=π(θ1,θ2|y)dθ1=π2|1(θ2|θ1,y)π1(θ1|y)dθ1=E(π2|1(θ2|θ1,y))
θ2θ2θ2θ1θ2

Isso é interessante devido ao lema da decomposição da variância onde a variância condicional é . Além disso, . Em particular, Um amostrador de Gibbs nos dará realizações . O resultado é que é melhor estimar por que em algumas estimativas convencionais de densidade de kernel usando o para o ponto

Var(X)=E[Var(X|Y)]+Var[E(X|Y)],
Var(X|Y)E{(XE(X|Y))2|Y}Var(E(X|Y))=E[(E(X|Y)E(X))2]
Var(X)Var[E(X|Y)].
(θ1i,θ2i)π2(θ2|y)
π^2(θ2|y)=1Mi=1Mπ2|1(θ2|θ1i,y)
θ2iθ2 - desde que conheçamos as distribuições condicionais (razão pela qual usamos a amostragem de Gibbs em primeiro lugar).

Exemplo

Suponha que e sejam bivariados normais com médias zero, variâncias 1 e correlação . Ou seja, Claramente, marginalmente, , mas vamos fingir que não sabemos disso. É sabido que a distribuição condicional de dada é .XYρ

π(x,y)exp{12(1ρ2)(x2+y22ρxy)}
YN(0,1)YX=xN(ρx,1ρ2)

Dadas algumas realizações de a estimativa "Rao-Blackwell" da densidade de em é Como ilustração, vamos comparar uma estimativa de densidade do kernel com a abordagem RBM(X,Y)Yy

π^Y(y)=1Mi=1M11ρ22πexp{12(1ρ2)(yρxi)2}
library(mvtnorm)

rho <- 0.5
R <- 50
xy <- rmvnorm(n=R, mean=c(0,0), sigma= matrix(c(1,rho,rho,1), ncol=2))
x <- xy[,1]
y <- xy[,2]

kernel_density <- density(y, kernel = "gaussian")
plot(kernel_density,col = "blue",lty=2,main="Rao-Blackwell estimates from conditional normals",ylim=c(0,0.4))
legend(1.5,.37,c("Kernel","N(0,1)","Rao-Blackwell"),lty=c(2,1,3),col=c("blue","black","red"))
g <- seq(-3.5,3.5,length=100)
lines(g,dnorm(g),lty=1) # here's what we pretend not to know

density_RB <- rep(0,100)
for(i in 1:100) {density_RB[i] <- mean(dnorm(g[i], rho*x, sd = sqrt(1-rho^2)))}
lines(g,density_RB,col = "red",lty=3) 

Observamos que a estimativa RB se sai muito melhor (pois explora as informações condicionais):

insira a descrição da imagem aqui

Christoph Hanck
fonte