Rao-Blackwellization de Gibbs Sampler

9

Actualmente estoy estimando un modelo de volatilidad estocástica con los métodos de Markov Chain Monte Carlo. De este modo, estoy implementando los métodos de muestreo de Gibbs y Metropolis.

Suponiendo que tomo la media de la distribución posterior en lugar de una muestra aleatoria, ¿es esto lo que comúnmente se conoce como Rao-Blackwellization ?

En general, esto resultaría en tomar la media sobre las medias de las distribuciones posteriores como estimación de parámetros.

mscnvrsy
fuente

Respuestas:

10

Suponiendo que tomo la media de la distribución posterior en lugar de una muestra aleatoria, ¿es esto lo que comúnmente se conoce como Rao-Blackwellization?

No estoy muy familiarizado con los modelos de volatilidad estocástica, pero sé que en la mayoría de los entornos, la razón por la que elegimos los algoritmos de Gibbs o MH para dibujar desde la parte posterior, es porque no conocemos la parte posterior. A menudo queremos estimar la media posterior, y como no conocemos la media posterior, extraemos muestras de la posterior y la estimamos utilizando la media de la muestra. Entonces, no estoy seguro de cómo podrá tomar la media de la distribución posterior.

En cambio, el estimador Rao-Blackwellized depende del conocimiento de la media del condicional completo; pero aun así se requiere muestreo. Os explico más a continuación.

Suponga que la distribución posterior se define en dos variables, ), de modo que desee estimar la media posterior: . Ahora, si hubiera una muestra de Gibbs disponible, podría ejecutar eso o ejecutar un algoritmo MH para muestrear desde la parte posterior.θ=(μ,ϕE[θdata]

Si puede ejecutar una muestra de Gibbs, entonces conoce en forma cerrada y conoce la media de esta distribución. Deje que eso significa ser . Tenga en cuenta que es una función de y los datos.f(ϕμ,data)ϕϕμ

Esto también significa que puede integrar desde la parte posterior, por lo que la parte posterior marginal de es (esto no se conoce completamente, pero se sabe hasta una constante). Ahora desea ejecutar una cadena de Markov tal que sea ​​la distribución invariable, y obtenga muestras de este posterior marginal. La pregunta esϕμf(μdata)f(μdata)

¿Cómo puede estimar ahora la media posterior de utilizando solo estas muestras del marginal posterior de ?ϕμ

Esto se hace a través de Rao-Blackwellization.

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

Supongamos que hemos obtenido muestras del margen posterior de . Entonces X1,X2,XNμ

ϕ^=1norteyo=1norteϕ(Xyo),

se llama estimador Rao-Blackwellized para . Lo mismo se puede hacer simulando también desde los márgenes articulares.ϕ

Ejemplo (puramente para demostración).

Suponga que tiene una articulación posterior desconocida para de la que desea muestrear. Sus datos son algunos , y tiene los siguientes condicionales completos θ=(μ,ϕ)y

μϕ,ynorte(ϕ2+2y,y2)
ϕμ,ysolunametrometrouna(2μ+y,y+1)

Ejecutas la muestra de Gibbs con estos condicionales y obtienes muestras de la articulación posterior . Deje que estas muestras sean . Puede encontrar la media muestral de s, y ese sería el estimador habitual de Monte Carlo para la media posterior de ..F(μ,ϕy)(μ1,ϕ1),(μ2,ϕ2),...,(μnorte,ϕnorte)ϕϕ

O tenga en cuenta que por las propiedades de la distribución Gamma

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

Aquí son los datos que se le proporcionan y, por lo tanto, se conocen. El estimador Rao Blackwellized sería entoncesy

ϕ^=1norteyo=1norte2μyo+yy+1.

Observe cómo el estimador para la media posterior de ni siquiera usa las muestras , y solo usa las muestras . En cualquier caso, como puede ver, todavía está utilizando las muestras que obtuvo de una cadena de Markov. Este no es un proceso determinista.ϕϕμ

Greenparker
fuente
Entonces, suponiendo que se conozca la distribución posterior del parámetro (que, según mi leal saber y entender, es cierto cuando se aplica el muestreo de Gibbs), tomar la media de la distribución en lugar de una muestra aleatoria sería el estimador Rao-Blackwellized. Espero haber entendido tu respuesta correctamente. Muchas gracias ya!
mscnvrsy
Eso es incorrecto. En el muestreo de Gibbs, no conoce la distribución posterior del parámetro, pero conoce el posterior condicional completo para cada parámetro. Hay una gran diferencia entre los dos. Arriba, la parte posterior es que se desconoce, y para que funcione la muestra de Gibbs necesita saber tanto como . Y también eres incorrecto en tu segunda comprensión. Todavía necesita tomar una muestra del margen posterior de , y luego calcular la media muestral de usando esas muestras para encontrar el estimador RB. F(μ,ϕreunatuna)F(μϕ,reunatuna)F(ϕμ,reunatuna)μϕ
Greenparker
@mscnvrsy Agregué un ejemplo para ayudar
Greenparker
Wow, muchas gracias por aclararme esto. Entonces, suponiendo que conozco las distribuciones condicionales completas, ¿puedo trabajar con los medios teóricos de las distribuciones condicionales y el promedio sobre estos medios teóricos (como E [phi | mu, y]) para obtener el estimador RB? ¿Esto minimizaría la varianza de mis estimaciones de parámetros?
mscnvrsy
Si estaba obteniendo muestras independientes, sí, minimizaría la varianza de los estimadores, sin embargo, dado que se trata de cadenas de Markov, generalmente se sabe que RB no necesariamente reduce la varianza, y hay algunos casos en los que la varianza incluso aumenta. Este artículo de Charlie Geyer dio algunos ejemplos a este punto.
Greenparker
9

La muestra de Gibbs se puede usar para mejorar la eficiencia de (digamos) muestras de un posterior marginal, . Nota Por lo tanto, el la densidad marginal de en algún valor es el valor esperado de la densidad condicional de dada en el punto .π2(θ2El |y)

π2(θ2El |y)=π(θ1,θ2El |y)reθ1=π2El |1(θ2El |θ1,y)π1(θ1El |y)reθ1=mi(π2El |1(θ2El |θ1,y))
θ2θ2θ2θ1θ2

Esto es interesante debido a la lema de descomposición de varianza donde la varianza condicional es . Además, . En particular, Una muestra de Gibbs nos dará realizaciones . El resultado es que es mejor estimar por que por alguna estimación de densidad de kernel convencional usando para el punto

Vunar(X)=mi[Vunar(XEl |Y)]+Vunar[mi(XEl |Y)],
Vunar(XEl |Y)mi{(X-mi(XEl |Y))2El |Y}Vunar(mi(XEl |Y))=mi[(mi(XEl |Y)-mi(X))2]
Vunar(X)Vunar[mi(XEl |Y)].
(θ1yo,θ2yo)π2(θ2El |y)
π^2(θ2El |y)=1METROyo=1METROπ2El |1(θ2El |θ1yo,y)
θ2yoθ2 - siempre que conozcamos las distribuciones condicionales (que es, por supuesto, por qué usamos el muestreo de Gibbs en primer lugar).

Ejemplo

Supongamos que e son normales bivariadas con medias cero, varianzas 1 y correlación . Es decir, Claramente, marginalmente, , pero supongamos que no sabemos esto. Es bien sabido que la distribución condicional de dado es .XYρ

π(X,y)Exp{-12(1-ρ2)(X2+y2-2ρXy)}
Ynorte(0 0,1)YX=Xnorte(ρX,1-ρ2)

Dadas algunas realizaciones de la estimación "Rao-Blackwell" de la densidad de en entonces es Como ilustración, comparemos una estimación de densidad del núcleo con el enfoque RBMETRO(X,Y)Yy

π^Y(y)=1METROyo=1METRO11-ρ22πExp{-12(1-ρ2)(y-ρXyo)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 la estimación de RB funciona mucho mejor (ya que explota la información condicional):

ingrese la descripción de la imagen aquí

Christoph Hanck
fuente