Rao-Blackwellization de Gibbs Sampler

J'évalue actuellement un modèle de volatilité stochastique avec les méthodes de Markov Chain Monte Carlo. Ainsi, j'implémente les méthodes d'échantillonnage de Gibbs et Metropolis.

En supposant que je prenne la moyenne de la distribution postérieure plutôt qu'un échantillon aléatoire, est-ce ce que l'on appelle communément Rao-Blackwellization ?

Dans l'ensemble, cela se traduirait par la prise de la moyenne sur la moyenne des distributions postérieures comme estimation des paramètres.

En supposant que je prenne la moyenne de la distribution postérieure plutôt qu'un échantillon aléatoire, est-ce ce que l'on appelle communément Rao-Blackwellization?

Je ne connais pas très bien les modèles de volatilité stochastique, mais je sais que dans la plupart des contextes, la raison pour laquelle nous choisissons les algorithmes de Gibbs ou MH pour dessiner à partir du postérieur, c'est parce que nous ne connaissons pas le postérieur. Souvent, nous voulons estimer la moyenne postérieure, et comme nous ne connaissons pas la moyenne postérieure, nous prélevons des échantillons sur la partie postérieure et l'estimons en utilisant la moyenne de l'échantillon. Donc, je ne sais pas comment vous pourrez prendre la moyenne de la distribution postérieure.

L'estimateur Rao-Blackwellized dépend plutôt de la connaissance de la moyenne du conditionnel complet; mais même dans ce cas, l'échantillonnage est toujours nécessaire. J'explique plus ci-dessous.

Supposons que la distribution postérieure soit définie sur deux variables, $\theta = (\mu, \phi$), de sorte que vous souhaitez estimer la moyenne postérieure: $E[\theta \mid \text{data}]$. Maintenant, si un échantillonneur Gibbs était disponible, vous pouvez l'exécuter ou exécuter un algorithme MH pour échantillonner à partir de la partie postérieure.

Si vous pouvez exécuter un échantillonneur Gibbs, alors vous savez $f(\phi \mid \mu, data)$sous forme fermée et vous connaissez la moyenne de cette distribution. Que cela signifie$\phi^*$. Notez que$\phi^*$ est fonction de $\mu$ et les données.

Cela signifie également que vous pouvez intégrer $\phi$ du postérieur, donc le postérieur marginal de $\mu$ est $f(\mu \mid data)$(ce n'est pas connu complètement, mais connu jusqu'à une constante). Vous voulez maintenant exécuter une chaîne de Markov telle que$f(\mu \mid data)$est la distribution invariante, et vous obtenez des échantillons de cette marge marginale. La question est

Comment pouvez-vous maintenant estimer la moyenne postérieure de $\phi$ en utilisant uniquement ces échantillons de la partie postérieure marginale de $\mu$?

Cela se fait via Rao-Blackwellization.

$$\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*}$$

Supposons donc que nous ayons obtenu des échantillons $X_1, X_2, \dots X_N$ du postérieur marginal de $\mu$. alors

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

est appelé l'estimateur Rao-Blackwellized pour $\phi$. La même chose peut être faite en simulant également à partir des marginaux conjoints.

Exemple (purement pour démonstration).

Supposons que vous ayez une articulation postérieure inconnue $\theta = (\mu, \phi)$à partir de laquelle vous souhaitez échantillonner. Vos données en sont$y$et vous disposez des conditions complètes suivantes

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

Vous exécutez l'échantillonneur Gibbs à l'aide de ces conditions et vous avez obtenu des échantillons de l'articulation postérieure $f(\mu, \phi \mid y)$. Que ces échantillons soient$(\mu_1, \phi_1), (\mu_2, \phi_2), \dots, (\mu_N, \phi_N)$. Vous pouvez trouver la moyenne de l'échantillon de la$\phi$s, et ce serait l'estimateur de Monte Carlo habituel pour la moyenne postérieure pour $\phi$..

Ou, notez que par les propriétés de la distribution Gamma

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

Ici $y$sont les données qui vous sont données et sont donc connues. L'estimateur Rao Blackwellized serait alors

$$\hat{\phi} = \dfrac{1}{N} \sum_{i=1}^{N} \dfrac{2 \mu_i + y}{y + 1}.$$

Remarquez comment l'estimateur de la moyenne postérieure de $\phi$ n'utilise même pas le $\phi$ échantillons, et utilise uniquement le $\mu$échantillons. Dans tous les cas, comme vous pouvez le voir, vous utilisez toujours les échantillons que vous avez obtenus d'une chaîne de Markov. Ce n'est pas un processus déterministe.

L'échantillonneur Gibbs peut ensuite être utilisé pour améliorer l'efficacité (par exemple) des échantillons à partir d'un postérieur marginal, appelez-le $\pi_2(\theta_2|y)$. Remarque

$$\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*}$$ Ainsi, la densité marginale de $\theta_2$ à une certaine valeur $\theta_2$ est la valeur attendue de la densité conditionnelle de $\theta_2$ donné $\theta_1$ à ce point $\theta_2$.

Ceci est intéressant en raison du lemme de décomposition de la variance

$$Var(X)=E[Var(X|Y)]+Var[E(X|Y)],$$ où la variance conditionnelle $Var(X|Y)$ est $E\left\{(X-E(X|Y))^2|Y\right\}$. Aussi,$Var(E(X|Y))=E\left[(E(X|Y)-E(X))^2\right]$. En particulier, $$Var(X)\geq Var[E(X|Y)].$$ Un échantillonneur Gibbs nous donnera des réalisations $(\theta_{1i},\theta_{2i})$. Le résultat est qu'il vaut mieux estimer$\pi_2(\theta_2|y)$ par $$\hat{\pi}_2(\theta_2|y)=\frac{1}{M}\sum_{i=1}^M\pi_{2|1}(\theta_2|\theta_{1i},y)$$ que par une estimation conventionnelle de la densité du noyau en utilisant le $\theta_{2i}$ pour le point $\theta_2$ - à condition de connaître les distributions conditionnelles (ce qui est bien sûr la raison pour laquelle nous utilisons l'échantillonnage de Gibbs en premier lieu).

Exemple

Supposer $X$ et $Y$ sont normaux bivariés avec des moyennes nulles, des variances 1 et une corrélation $\rho$. C'est,

$$\pi(x,y)\propto\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2+y^2-2\rho x y)\right\}$$ Clairement, marginalement, $Y\sim N(0,1)$, mais supposons que nous ne le savons pas. Il est bien connu que la distribution conditionnelle de$Y$ donné $X=x$ est $N(\rho x,1-\rho^2)$.

Étant donné certains $M$ réalisations de $(X,Y)$ l'estimation "Rao-Blackwell" de la densité de $Y$ à $y$ alors c'est

$$\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\}$$ À titre d'illustration, comparons une estimation de la densité du noyau à l'approche RB

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) 

Nous observons que l'estimation RB fait beaucoup mieux (car elle exploite les informations conditionnelles):