Rao-Blackwellization de Gibbs Sampler

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.

mcmc monte-carlo gibbs point-estimation rao-blackwell mscnvrsy
fonte

Respostas:

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. $\theta = (\mu, \phi$ $E[\theta \mid \text{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(\phi \mid \mu, data)$ $\phi^*$ $\phi^*$ $\mu$

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 é $\phi$ $\mu$ $f(\mu \mid data)$ $f(\mu \mid data)$

Como você pode agora estimar a média posterior de usando apenas essas amostras da margem posterior de ? $\phi$ $\mu$

Isso é feito via Rao-Blackwellization.

\begin{aligned} E [ϕ ∣ d a t a] & = \int ϕ f (μ, ϕ ∣ d a t a) d μ d ϕ \\ = \int ϕ f (ϕ ∣ μ, d a t a) f (μ ∣ d a t a) d μ d ϕ \\ = \int ϕ^{*} f (μ ∣ d a t a) d μ . \end{aligned}

$\begin{align*} E[\phi \mid data]& = \int \phi \; f(\mu, \phi \mid data) d\mu \, d\phi\\ & = \int \phi \; f(\phi \mid \mu, data) f(\mu \mid data) d\mu \, d\phi\\ & = \int \phi^* f(\mu \mid data) d\mu. \end{align*}$

Assim, suponha que tenhamos obtido as amostras partir da margem posterior de . Então $X_1, X_2, \dots X_N$ $\mu$

\hat{ϕ} = \frac{1}{N} \sum_{i = 1}^{N} ϕ^{*} (X_{i}),

$\hat{\phi} = \dfrac{1}{N} \sum_{i=1}^{N} \phi^*(X_i),$

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

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 $\theta = (\mu, \phi)$ $y$

μ ∣ ϕ, y \sim N (ϕ^{2} + 2 y, y^{2})

$\mu \mid \phi, y \sim N(\phi^2 + 2y, y^2)$

ϕ ∣ μ, y \sim G a m m a (2 μ + y, y + 1)

$\phi \mid \mu, y \sim Gamma(2\mu + 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 . $f(\mu, \phi \mid y)$ $(\mu_1, \phi_1), (\mu_2, \phi_2), \dots, (\mu_N, \phi_N)$ $\phi$ $\phi$

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

E [ϕ | μ, y] = \frac{2 μ + y}{y + 1} = ϕ^{*} .

$E[\phi | \mu, y] = \dfrac{2 \mu + y}{y + 1} = \phi^*.$

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

\hat{ϕ} = \frac{1}{N} \sum_{i = 1}^{N} \frac{2 μ_{i} + y}{y + 1} .

$\hat{\phi} = \dfrac{1}{N} \sum_{i=1}^{N} \dfrac{2 \mu_i + y}{y + 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. $\phi$ $\phi$ $\mu$

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 (μ, ϕ ∣ d a t a)

$f(\mu, \phi \mid data)$

f (μ ∣ ϕ, d a t a)

$f(\mu \mid \phi, data)$

f (ϕ ∣ μ, d a t a)

$f(\phi \mid \mu, data)$

μ

$\mu$

ϕ^{*}

$\phi^*$

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

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 . $\pi_2(\theta_2|y)$

\begin{array}{rcl} π_{2} (θ_{2} | y) & = & \int π (θ_{1}, θ_{2} | y) d θ_{1} \\ = & \int π_{2 | 1} (θ_{2} | θ_{1}, y) π_{1} (θ_{1} | y) d θ_{1} \\ = & E (π_{2 | 1} (θ_{2} | θ_{1}, y)) \end{array}

$\begin{eqnarray*} \pi_2(\theta_2|y)&=&\int \pi(\theta_1,\theta_2|y)d\theta_1\\ &=&\int \pi_{2|1}(\theta_2|\theta_1,y)\pi_1(\theta_1|y)d\theta_1\\ &=&E(\pi_{2|1}(\theta_2|\theta_1,y)) \end{eqnarray*}$

θ_{2}

$\theta_2$

θ_{2}

$\theta_2$

θ_{2}

$\theta_2$

θ_{1}

$\theta_1$

θ_{2}

$\theta_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

V a r (X) = E [V a r (X | Y)] + V a r [E (X | Y)],

$Var(X)=E[Var(X|Y)]+Var[E(X|Y)],$

V a r (X | Y)

$Var(X|Y)$

E {(X - E (X | Y))^{2} | Y}

$E\left\{(X-E(X|Y))^2|Y\right\}$

V a r (E (X | Y)) = E [(E (X | Y) - E (X))^{2}]

$Var(E(X|Y))=E\left[(E(X|Y)-E(X))^2\right]$

V a r (X) \geq V a r [E (X | Y)] .

$Var(X)\geq Var[E(X|Y)].$

(θ_{1 i}, θ_{2 i})

$(\theta_{1i},\theta_{2i})$

π_{2} (θ_{2} | y)

$\pi_2(\theta_2|y)$

{\hat{π}}_{2} (θ_{2} | y) = \frac{1}{M} \sum_{i = 1}^{M} π_{2 | 1} (θ_{2} | θ_{1 i}, y)

$\hat{\pi}_2(\theta_2|y)=\frac{1}{M}\sum_{i=1}^M\pi_{2|1}(\theta_2|\theta_{1i},y)$

θ_{2 i}

$\theta_{2i}$

θ_{2}

$\theta_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 é . $X$ $Y$ $\rho$

π (x, y) \propto \exp {- \frac{1}{2 (1 - ρ^{2})} (x^{2} + y^{2} - 2 ρ x y)}

$\pi(x,y)\propto\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2+y^2-2\rho x y)\right\}$

Y \sim N (0, 1)

$Y\sim N(0,1)$

Y

$Y$

X = x

$X=x$

N (ρ x, 1 - ρ^{2})

$N(\rho x,1-\rho^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 RB $M$ $(X,Y)$ $Y$ $y$

{\hat{π}}_{Y} (y) = \frac{1}{M} \sum_{i = 1}^{M} \frac{1}{\sqrt{1 - ρ^{2}} \sqrt{2 π}} \exp {- \frac{1}{2 (1 - ρ^{2})} (y - ρ x_{i})^{2}}

$\hat\pi_Y(y)=\frac{1}{M}\sum_{i=1}^M\frac{1}{\sqrt{1-\rho^2}\sqrt{2\pi}}\exp\left\{-\frac{1}{2(1-\rho^2)}(y-\rho x_i)^2\right\}$

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):

Christoph Hanck
fonte