Métodos bayesianos variacionales

En la inferencia variacional , la distribución posterior sobre un conjunto de variables no observadas${\displaystyle \mathbf {Z} =\{Z_{1}\dots Z_{n}\}}$dados algunos datos${\displaystyle \mathbf {X} }$se aproxima mediante una denominada distribución variacional ,${\displaystyle Q(\mathbf {Z} ):}$

${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )\approx Q(\mathbf {Z} ).}$

La distribución${\displaystyle Q(\mathbf {Z} )}$está restringido a pertenecer a una familia de distribuciones de forma más simple que${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$(por ejemplo, una familia de distribuciones gaussianas), seleccionadas con la intención de hacer${\displaystyle Q(\mathbf {Z} )}$similar a la verdadera posterior,${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$.

La similitud (o disimilitud) se mide en términos de una función de disimilitud.${\displaystyle d(Q;P)}$y por lo tanto la inferencia se realiza seleccionando la distribución${\displaystyle Q(\mathbf {Z} )}$que minimiza${\displaystyle d(Q;P)}$.

El tipo más común de Bayes variacional utiliza la divergencia de Kullback-Leibler (divergencia KL) de Q con respecto a P como función de disimilitud. Esta elección hace que esta minimización sea manejable. La divergencia KL se define como

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)\triangleq \sum _{\mathbf {Z} }Q(\mathbf {Z} )\log {\frac {Q(\mathbf {Z} )}{P(\mathbf {Z} \mid \mathbf {X} )}}.}$

Nótese que Q y P están invertidos respecto a lo que cabría esperar. Este uso de la divergencia KL invertida es conceptualmente similar al algoritmo de maximización de la expectativa . (Utilizar la divergencia KL en sentido contrario produce el algoritmo de propagación de la expectativa ).

Las técnicas variacionales se utilizan normalmente para formar una aproximación para:

${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{P(\mathbf {X} )}}={\frac {P(\mathbf {X} \mid \mathbf {Z} )P(\mathbf {Z} )}{\int _{\mathbf {Z} }P(\mathbf {X} ,\mathbf {Z} ')\,d\mathbf {Z} '}}}$

La marginación sobre${\displaystyle \mathbf {Z} }$calcular${\displaystyle P(\mathbf {X} )}$en el denominador suele ser intratable, porque, por ejemplo, el espacio de búsqueda de${\displaystyle \mathbf {Z} }$es combinatoriamente grande. Por lo tanto, buscamos una aproximación, utilizando${\displaystyle Q(\mathbf {Z} )\approx P(\mathbf {Z} \mid \mathbf {X} )}$.

Dado que${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} ,\mathbf {Z} )}{P(\mathbf {X} )}}}$, la divergencia KL anterior también se puede escribir como

${\displaystyle {\begin{array}{rl}D_{\mathrm {KL} }(Q\parallel P)&=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log {\frac {Q(\mathbf {Z} )}{P(\mathbf {Z} ,\mathbf {X} )}}+\log P(\mathbf {X} )\right]\\&=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log P(\mathbf {X} )\right]\end{array}}}$

Porque${\displaystyle P(\mathbf {X} )}$es una constante con respecto a${\displaystyle \mathbf {Z} }$y${\displaystyle \sum _{\mathbf {Z} }Q(\mathbf {Z} )=1}$porque${\displaystyle Q(\mathbf {Z} )}$es una distribución, tenemos

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\log P(\mathbf {X} )}$

que, según la definición de valor esperado (para una variable aleatoria discreta ), puede escribirse de la siguiente manera:

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)=\mathbb {E} _{\mathbf {Q} }\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\log P(\mathbf {X} )}$

que se puede reorganizar para convertirse en

${\displaystyle {\begin{array}{rl}\log P(\mathbf {X} )&=D_{\mathrm {KL} }(Q\parallel P)-\mathbb {E} _{\mathbf {Q} }\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]\\&=D_{\mathrm {KL} }(Q\parallel P)+{\mathcal {L}}(Q)\end{array}}}$

Como evidencia del registro${\displaystyle \log P(\mathbf {X} )}$está fijo con respecto a${\displaystyle Q}$, maximizando el término final${\displaystyle {\mathcal {L}}(Q)}$minimiza la divergencia KL de${\displaystyle Q}$de${\displaystyle P}$. Mediante la elección apropiada de${\displaystyle Q}$,${\displaystyle {\mathcal {L}}(Q)}$se vuelve manejable para calcular y maximizar. Por lo tanto, tenemos una aproximación analítica.${\displaystyle Q}$para la parte posterior${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$y un límite inferior${\displaystyle {\mathcal {L}}(Q)}$para la evidencia del registro${\displaystyle \log P(\mathbf {X} )}$(dado que la divergencia KL no es negativa).

El límite inferior${\displaystyle {\mathcal {L}}(Q)}$se conoce como energía libre variacional (negativa) por analogía con la energía libre termodinámica porque también puede expresarse como una energía negativa${\displaystyle \operatorname {E} _{Q}[\log P(\mathbf {Z} ,\mathbf {X} )]}$más la entropía de${\displaystyle Q}$. El término${\displaystyle {\mathcal {L}}(Q)}$También se conoce como Límite Inferior de Evidencia , abreviado como ELBO , para enfatizar que es un límite inferior (en el peor de los casos) de la evidencia logarítmica de los datos.

Mediante el teorema pitagórico generalizado de la divergencia de Bregman , del cual la divergencia KL es un caso especial, se puede demostrar que: ^{[ 1 ] [ 2 ]}

${\displaystyle D_{\mathrm {KL} }(Q\parallel P)\geq D_{\mathrm {KL} }(Q\parallel Q^{*})+D_{\mathrm {KL} }(Q^{*}\parallel P),\forall Q^{*}\in {\mathcal {C}}}$

dónde ${\displaystyle {\mathcal {C}}}$es un conjunto convexo y la igualdad se cumple si:

${\displaystyle Q=Q^{*}\triangleq \arg \min _{Q\in {\mathcal {C}}}D_{\mathrm {KL} }(Q\parallel P).}$

En este caso, el minimizador global${\displaystyle Q^{*}(\mathbf {Z} )=q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})q^{*}(\mathbf {Z} _{2})=q^{*}(\mathbf {Z} _{2}\mid \mathbf {Z} _{1})q^{*}(\mathbf {Z} _{1}),}$con${\displaystyle \mathbf {Z} =\{\mathbf {Z_{1}} ,\mathbf {Z_{2}} \},}$se puede encontrar de la siguiente manera: ^{[ 1 ]}

${\displaystyle {\begin{array}{rl}q^{*}(\mathbf {Z} _{2})&={\frac {P(\mathbf {X} )}{\zeta (\mathbf {X} )}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}\\&={\frac {1}{\zeta (\mathbf {X} )}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right),\end{array}}}$

en la que la constante de normalización es:

${\displaystyle {\begin{array}{rl}\zeta (\mathbf {X} )&=P(\mathbf {X} )\int _{\mathbf {Z} _{2}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}\\&=\int _{\mathbf {Z} _{2}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right).\end{array}}}$

El término${\displaystyle \zeta (\mathbf {X} )}$En la práctica, a menudo se le llama límite inferior de la evidencia ( ELBO ), ya que${\displaystyle P(\mathbf {X} )\geq \zeta (\mathbf {X} )=\exp({\mathcal {L}}(Q^{*}))}$, ^{[ 1 ]} como se muestra arriba.

Al intercambiar los roles de${\displaystyle \mathbf {Z} _{1}}$y${\displaystyle \mathbf {Z} _{2},}$podemos calcular iterativamente el aproximado${\displaystyle q^{*}(\mathbf {Z} _{1})}$y${\displaystyle q^{*}(\mathbf {Z} _{2})}$de los marginales del modelo verdadero${\displaystyle P(\mathbf {Z} _{1}\mid \mathbf {X} )}$y${\displaystyle P(\mathbf {Z} _{2}\mid \mathbf {X} ),}$respectivamente. Aunque este esquema iterativo tiene garantizada la convergencia monótona, ^{[ 1 ]} la convergencia${\displaystyle Q^{*}}$es solo un minimizador local de${\displaystyle D_{\mathrm {KL} }(Q\parallel P)}$.

Si el espacio restringido${\displaystyle {\mathcal {C}}}$está confinado dentro de un espacio independiente, es decir${\displaystyle q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})=q^{*}(\mathbf {Z_{1}} ),}$El esquema iterativo anterior se convertirá en la denominada aproximación de campo medio.${\displaystyle Q^{*}(\mathbf {Z} )=q^{*}(\mathbf {Z} _{1})q^{*}(\mathbf {Z} _{2}),}$como se muestra a continuación.

Colocamos distribuciones previas conjugadas sobre la media desconocida.${\displaystyle \mu }$y precisión${\displaystyle \tau }$, es decir, la media también sigue una distribución gaussiana, mientras que la precisión sigue una distribución gamma . En otras palabras:

${\displaystyle {\begin{aligned}\tau &\sim \operatorname {Gamma} (a_{0},b_{0})\\\mu |\tau &\sim {\mathcal {N}}(\mu _{0},(\lambda _{0}\tau )^{-1})\\\{x_{1},\dots ,x_{N}\}&\sim {\mathcal {N}}(\mu ,\tau ^{-1})\\N&={\text{number of data points}}\end{aligned}}}$

Los hiperparámetros${\displaystyle \mu _{0},\lambda _{0},a_{0}}$y${\displaystyle b_{0}}$en las distribuciones previas son valores fijos y dados. Se pueden establecer en números positivos pequeños para dar distribuciones previas amplias que indiquen ignorancia sobre las distribuciones previas de${\displaystyle \mu }$y${\displaystyle \tau }$.

Se nos da${\displaystyle N}$puntos de datos${\displaystyle \mathbf {X} =\{x_{1},\ldots ,x_{N}\}}$y nuestro objetivo es inferir la distribución posterior${\displaystyle q(\mu ,\tau )=p(\mu ,\tau \mid x_{1},\ldots ,x_{N})}$de los parámetros${\displaystyle \mu }$y${\displaystyle \tau .}$

La probabilidad conjunta de todas las variables se puede reescribir como

${\displaystyle p(\mathbf {X} ,\mu ,\tau )=p(\mathbf {X} \mid \mu ,\tau )p(\mu \mid \tau )p(\tau )}$

donde los factores individuales son

${\displaystyle {\begin{aligned}p(\mathbf {X} \mid \mu ,\tau )&=\prod _{n=1}^{N}{\mathcal {N}}(x_{n}\mid \mu ,\tau ^{-1})\\p(\mu \mid \tau )&={\mathcal {N}}\left(\mu \mid \mu _{0},(\lambda _{0}\tau )^{-1}\right)\\p(\tau )&=\operatorname {Gamma} (\tau \mid a_{0},b_{0})\end{aligned}}}$

dónde

${\displaystyle {\begin{aligned}{\mathcal {N}}(x\mid \mu ,\sigma ^{2})&={\frac {1}{\sqrt {2\pi \sigma ^{2}}}}e^{\frac {-(x-\mu )^{2}}{2\sigma ^{2}}}\\\operatorname {Gamma} (\tau \mid a,b)&={\frac {1}{\Gamma (a)}}b^{a}\tau ^{a-1}e^{-b\tau }\end{aligned}}}$

Supongamos que${\displaystyle q(\mu ,\tau )=q(\mu )q(\tau )}$, es decir, que la distribución posterior se factoriza en factores independientes para${\displaystyle \mu }$y${\displaystyle \tau }$Este tipo de suposición subyace al método bayesiano variacional. La verdadera distribución posterior no se comporta de esta manera (de hecho, en este caso sencillo, se sabe que es una distribución gaussiana-gamma ), por lo que el resultado que obtenemos será una aproximación.

Entonces

${\displaystyle {\begin{aligned}\ln q_{\mu }^{*}(\mu )&=\operatorname {E} _{\tau }\left[\ln p(\mathbf {X} \mid \mu ,\tau )+\ln p(\mu \mid \tau )+\ln p(\tau )\right]+C\\&=\operatorname {E} _{\tau }\left[\ln p(\mathbf {X} \mid \mu ,\tau )\right]+\operatorname {E} _{\tau }\left[\ln p(\mu \mid \tau )\right]+\operatorname {E} _{\tau }\left[\ln p(\tau )\right]+C\\&=\operatorname {E} _{\tau }\left[\ln \prod _{n=1}^{N}{\mathcal {N}}\left(x_{n}\mid \mu ,\tau ^{-1}\right)\right]+\operatorname {E} _{\tau }\left[\ln {\mathcal {N}}\left(\mu \mid \mu _{0},(\lambda _{0}\tau )^{-1}\right)\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\ln \prod _{n=1}^{N}{\sqrt {\frac {\tau }{2\pi }}}e^{-{\frac {(x_{n}-\mu )^{2}\tau }{2}}}\right]+\operatorname {E} _{\tau }\left[\ln {\sqrt {\frac {\lambda _{0}\tau }{2\pi }}}e^{-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}}\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}\left({\frac {1}{2}}(\ln \tau -\ln 2\pi )-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right)\right]+\operatorname {E} _{\tau }\left[{\frac {1}{2}}(\ln \lambda _{0}+\ln \tau -\ln 2\pi )-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}{\frac {1}{2}}(\ln \tau -\ln 2\pi )\right]+\operatorname {E} _{\tau }\left[{\frac {1}{2}}(\ln \lambda _{0}+\ln \tau -\ln 2\pi )\right]+C_{2}\\&=\operatorname {E} _{\tau }\left[\sum _{n=1}^{N}-{\frac {(x_{n}-\mu )^{2}\tau }{2}}\right]+\operatorname {E} _{\tau }\left[-{\frac {(\mu -\mu _{0})^{2}\lambda _{0}\tau }{2}}\right]+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right\}+C_{3}\end{aligned}}}$

En la derivación anterior,${\displaystyle C}$,${\displaystyle C_{2}}$y${\displaystyle C_{3}}$referirse a valores que son constantes con respecto a${\displaystyle \mu }$. Tenga en cuenta que el término${\displaystyle \operatorname {E} _{\tau }[\ln p(\tau )]}$no es una función de${\displaystyle \mu }$y tendrá el mismo valor independientemente del valor de${\displaystyle \mu }$Por lo tanto, en la línea 3 podemos absorberlo en el término constante al final. Hacemos lo mismo en la línea 7.

La última línea es simplemente un polinomio cuadrático en${\displaystyle \mu }$Dado que este es el logaritmo de${\displaystyle q_{\mu }^{*}(\mu )}$, podemos ver que${\displaystyle q_{\mu }^{*}(\mu )}$en sí misma es una distribución gaussiana .

Con una cierta cantidad de matemáticas tediosas (expandir los cuadrados dentro de las llaves, separar y agrupar los términos que involucran${\displaystyle \mu }$y${\displaystyle \mu ^{2}}$y completando el cuadrado sobre${\displaystyle \mu }$), podemos derivar los parámetros de la distribución gaussiana:

${\displaystyle {\begin{aligned}\ln q_{\mu }^{*}(\mu )&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\sum _{n=1}^{N}(x_{n}^{2}-2x_{n}\mu +\mu ^{2})+\lambda _{0}(\mu ^{2}-2\mu _{0}\mu +\mu _{0}^{2})\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{\left(\sum _{n=1}^{N}x_{n}^{2}\right)-2\left(\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}\mu ^{2}\right)+\lambda _{0}\mu ^{2}-2\lambda _{0}\mu _{0}\mu +\lambda _{0}\mu _{0}^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right\}+C_{3}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu \right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)(\lambda _{0}+N)\mu \right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu \right)\right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu +\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}-\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right)\right\}+C_{4}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu ^{2}-2\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)\mu +\left({\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right)\right\}+C_{5}\\&=-{\frac {\operatorname {E} _{\tau }[\tau ]}{2}}\left\{(\lambda _{0}+N)\left(\mu -{\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}\right\}+C_{5}\\&=-{\frac {1}{2}}(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\left(\mu -{\frac {\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}}{\lambda _{0}+N}}\right)^{2}+C_{5}\end{aligned}}}$

Tenga en cuenta que todos los pasos anteriores se pueden acortar utilizando la fórmula para la suma de dos ecuaciones cuadráticas .

En otras palabras:

${\displaystyle {\begin{aligned}q_{\mu }^{*}(\mu )&\sim {\mathcal {N}}(\mu \mid \mu _{N},\lambda _{N}^{-1})\\\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\end{aligned}}}$

La derivación de${\displaystyle q_{\tau }^{*}(\tau )}$Es similar a lo anterior, aunque omitimos algunos detalles en aras de la brevedad.

${\displaystyle {\begin{aligned}\ln q_{\tau }^{*}(\tau )&=\operatorname {E} _{\mu }[\ln p(\mathbf {X} \mid \mu ,\tau )+\ln p(\mu \mid \tau )]+\ln p(\tau )+{\text{constant}}\\&=(a_{0}-1)\ln \tau -b_{0}\tau +{\frac {1}{2}}\ln \tau +{\frac {N}{2}}\ln \tau -{\frac {\tau }{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]+{\text{constant}}\end{aligned}}}$

Al exponenciar ambos lados, podemos ver que${\displaystyle q_{\tau }^{*}(\tau )}$es una distribución gamma . Específicamente:

${\displaystyle {\begin{aligned}q_{\tau }^{*}(\tau )&\sim \operatorname {Gamma} (\tau \mid a_{N},b_{N})\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\end{aligned}}}$

Recapitulemos las conclusiones de las secciones anteriores:

${\displaystyle {\begin{aligned}q_{\mu }^{*}(\mu )&\sim {\mathcal {N}}(\mu \mid \mu _{N},\lambda _{N}^{-1})\\\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N)\operatorname {E} _{\tau }[\tau ]\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\end{aligned}}}$

y

${\displaystyle {\begin{aligned}q_{\tau }^{*}(\tau )&\sim \operatorname {Gamma} (\tau \mid a_{N},b_{N})\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\end{aligned}}}$

En cada caso, los parámetros de la distribución de una de las variables dependen de las esperanzas calculadas con respecto a la otra variable. Podemos desarrollar las esperanzas utilizando las fórmulas estándar para las esperanzas de los momentos de las distribuciones gaussiana y gamma:

${\displaystyle {\begin{aligned}\operatorname {E} [\tau \mid a_{N},b_{N}]&={\frac {a_{N}}{b_{N}}}\\\operatorname {E} \left[\mu \mid \mu _{N},\lambda _{N}^{-1}\right]&=\mu _{N}\\\operatorname {E} \left[X^{2}\right]&=\operatorname {Var} (X)+(\operatorname {E} [X])^{2}\\\operatorname {E} \left[\mu ^{2}\mid \mu _{N},\lambda _{N}^{-1}\right]&=\lambda _{N}^{-1}+\mu _{N}^{2}\end{aligned}}}$

Aplicar estas fórmulas a las ecuaciones anteriores es trivial en la mayoría de los casos, pero la ecuación para${\displaystyle b_{N}}$requiere más trabajo:

${\displaystyle {\begin{aligned}b_{N}&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[\sum _{n=1}^{N}(x_{n}-\mu )^{2}+\lambda _{0}(\mu -\mu _{0})^{2}\right]\\&=b_{0}+{\frac {1}{2}}\operatorname {E} _{\mu }\left[(\lambda _{0}+N)\mu ^{2}-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu +\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\operatorname {E} _{\mu }[\mu ^{2}]-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\operatorname {E} _{\mu }[\mu ]+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\left(\lambda _{N}^{-1}+\mu _{N}^{2}\right)-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu _{N}+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\\\end{aligned}}}$

Podemos entonces escribir las ecuaciones de los parámetros de la siguiente manera, sin ninguna expectativa:

${\displaystyle {\begin{aligned}\mu _{N}&={\frac {\lambda _{0}\mu _{0}+N{\bar {x}}}{\lambda _{0}+N}}\\\lambda _{N}&=(\lambda _{0}+N){\frac {a_{N}}{b_{N}}}\\{\bar {x}}&={\frac {1}{N}}\sum _{n=1}^{N}x_{n}\\a_{N}&=a_{0}+{\frac {N+1}{2}}\\b_{N}&=b_{0}+{\frac {1}{2}}\left[(\lambda _{0}+N)\left(\lambda _{N}^{-1}+\mu _{N}^{2}\right)-2\left(\lambda _{0}\mu _{0}+\sum _{n=1}^{N}x_{n}\right)\mu _{N}+\left(\sum _{n=1}^{N}x_{n}^{2}\right)+\lambda _{0}\mu _{0}^{2}\right]\end{aligned}}}$

Tenga en cuenta que existen dependencias circulares entre las fórmulas para${\displaystyle \lambda _{N}}$y${\displaystyle b_{N}}$Esto sugiere naturalmente un algoritmo similar al EM :

1. Calcular${\displaystyle \sum _{n=1}^{N}x_{n}}$y${\displaystyle \sum _{n=1}^{N}x_{n}^{2}.}$Utilice estos valores para calcular${\displaystyle \mu _{N}}$y${\displaystyle a_{N}.}$
2. Inicializar${\displaystyle \lambda _{N}}$a algún valor arbitrario.
3. Utilice el valor actual de${\displaystyle \lambda _{N},}$junto con los valores conocidos de los otros parámetros, para calcular${\displaystyle b_{N}}$.
4. Utilice el valor actual de${\displaystyle b_{N},}$junto con los valores conocidos de los otros parámetros, para calcular${\displaystyle \lambda _{N}}$.
5. Repita los dos últimos pasos hasta la convergencia (es decir, hasta que ninguno de los valores haya cambiado más que una pequeña cantidad).

Entonces tenemos valores para los hiperparámetros de las distribuciones aproximadas de los parámetros posteriores, que podemos usar para calcular cualquier propiedad que queramos de la distribución posterior, por ejemplo, su media y varianza, una región de máxima densidad del 95 % (el intervalo más pequeño que incluye el 95 % de la probabilidad total), etc.

Se puede demostrar que este algoritmo garantiza la convergencia a un máximo local.

Cabe destacar también que las distribuciones posteriores tienen la misma forma que las distribuciones previas correspondientes. No partimos de esta premisa; la única suposición que hicimos fue que las distribuciones se factorizaban, y la forma de las distribuciones se derivó de forma natural. Resulta (véase más adelante) que el hecho de que las distribuciones posteriores tengan la misma forma que las distribuciones previas no es una coincidencia, sino un resultado general siempre que las distribuciones previas pertenezcan a la familia exponencial , lo cual ocurre con la mayoría de las distribuciones estándar.

El ejemplo anterior muestra el método mediante el cual se deriva la aproximación variacional-bayesiana a una densidad de probabilidad posterior en una red bayesiana dada:

1. Describe la red con un modelo gráfico , identificando las variables observadas (datos).${\displaystyle \mathbf {X} }$y variables no observadas ( parámetros)${\displaystyle {\boldsymbol {\Theta }}}$y variables latentes${\displaystyle \mathbf {Z} }$) y sus distribuciones de probabilidad condicionales . El método bayesiano variacional construirá entonces una aproximación a la probabilidad posterior.${\displaystyle p(\mathbf {Z} ,{\boldsymbol {\Theta }}\mid \mathbf {X} )}$La aproximación tiene la propiedad básica de ser una distribución factorizada, es decir, un producto de dos o más distribuciones independientes sobre subconjuntos disjuntos de las variables no observadas.
2. Divida las variables no observadas en dos o más subconjuntos, sobre los cuales se derivarán los factores independientes. No existe un procedimiento universal para hacer esto; crear demasiados subconjuntos produce una aproximación deficiente, mientras que crear muy pocos hace que todo el procedimiento variacional bayesiano sea intratable. Normalmente, la primera división consiste en separar los parámetros y las variables latentes; a menudo, esto es suficiente por sí solo para producir un resultado manejable. Supongamos que las particiones se llaman${\displaystyle \mathbf {Z} _{1},\ldots ,\mathbf {Z} _{M}}$.
3. Para una partición dada${\displaystyle \mathbf {Z} _{j}}$Escribe la fórmula para la distribución que mejor se aproxime.${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )}$utilizando la ecuación básica${\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{i\neq j}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}}$.
4. Complete la fórmula para la distribución de probabilidad conjunta utilizando el modelo gráfico. Cualquier distribución condicional de componentes que no involucre ninguna de las variables en${\displaystyle \mathbf {Z} _{j}}$pueden ignorarse; se incorporarán al término constante.
5. Simplifica la fórmula y aplica el operador de expectativa, siguiendo el ejemplo anterior. Idealmente, esto debería simplificarse en expectativas de funciones básicas de variables que no están en${\displaystyle \mathbf {Z} _{j}}$(p. ej., primeros o segundos momentos brutos , esperanza de un logaritmo, etc.). Para que el procedimiento variacional bayesiano funcione correctamente, estas esperanzas deben poder expresarse analíticamente como funciones de los parámetros y/o hiperparámetros de las distribuciones de estas variables. En todos los casos, estos términos de esperanza son constantes con respecto a las variables en la partición actual.
6. La forma funcional de la fórmula con respecto a las variables en la partición actual indica el tipo de distribución. En particular, al exponenciar la fórmula se obtiene la función de densidad de probabilidad (FDP) de la distribución (o al menos, una proporcional a ella, con una constante de normalización desconocida ). Para que el método sea viable, debe ser posible reconocer la forma funcional como perteneciente a una distribución conocida. Puede requerirse una manipulación matemática significativa para convertir la fórmula en una forma que coincida con la FDP de una distribución conocida. Cuando esto se logra, la constante de normalización puede restablecerse por definición, y las ecuaciones para los parámetros de la distribución conocida pueden derivarse extrayendo las partes apropiadas de la fórmula.
7. Cuando todas las expectativas pueden reemplazarse analíticamente con funciones de variables que no están en la partición actual, y la función de densidad de probabilidad (PDF) se expresa en una forma que permite su identificación con una distribución conocida, el resultado es un conjunto de ecuaciones que expresan los valores de los parámetros óptimos como funciones de los parámetros de las variables en otras particiones.
8. Cuando este procedimiento se puede aplicar a todas las particiones, el resultado es un conjunto de ecuaciones interrelacionadas que especifican los valores óptimos de todos los parámetros.
9. A continuación, se aplica un procedimiento de maximización de la esperanza (EM), seleccionando un valor inicial para cada parámetro y recorriendo una serie de pasos. En cada paso, se actualizan las ecuaciones, modificando cada parámetro sucesivamente. Se garantiza la convergencia.

Debido a todas las manipulaciones matemáticas involucradas, es fácil perder de vista el panorama general. Lo importante es:

1. La idea del método bayesiano variacional es construir una aproximación analítica a la probabilidad posterior del conjunto de variables no observadas (parámetros y variables latentes), dados los datos. Esto significa que la forma de la solución es similar a otros métodos de inferencia bayesiana , como el muestreo de Gibbs , es decir, una distribución que busca describir todo lo que se sabe sobre las variables. Al igual que en otros métodos bayesianos, pero a diferencia, por ejemplo, del algoritmo de expectativa-maximización (EM) u otros métodos de máxima verosimilitud , ambos tipos de variables no observadas (parámetros y variables latentes) se tratan de la misma manera, es decir, como variables aleatorias . Las estimaciones para las variables se pueden obtener mediante los métodos bayesianos estándar, por ejemplo, calculando la media de la distribución para obtener una estimación puntual o derivando un intervalo creíble , una región de máxima densidad, etc.
2. La "aproximación analítica" significa que se puede escribir una fórmula para la distribución posterior. La fórmula generalmente consiste en un producto de distribuciones de probabilidad bien conocidas, cada una de las cuales se factoriza sobre un conjunto de variables no observadas (es decir, es condicionalmente independiente de las otras variables, dados los datos observados). Esta fórmula no es la verdadera distribución posterior, sino una aproximación a ella; en particular, generalmente coincidirá bastante bien en los momentos de menor orden de las variables no observadas, por ejemplo, la media y la varianza .
3. El resultado de todas las manipulaciones matemáticas es (1) la identidad de las distribuciones de probabilidad que componen los factores, y (2) fórmulas interdependientes para los parámetros de estas distribuciones. Los valores reales de estos parámetros se calculan numéricamente, mediante un procedimiento iterativo alterno muy similar al algoritmo EM.

El método bayesiano variacional (VB) se compara frecuentemente con el algoritmo de maximización de la esperanza (EM). El procedimiento numérico es bastante similar, ya que ambos son procedimientos iterativos alternados que convergen sucesivamente hacia valores óptimos de los parámetros. Los pasos iniciales para derivar los respectivos procedimientos también son vagamente similares, pues ambos parten de fórmulas para densidades de probabilidad e implican una cantidad considerable de manipulaciones matemáticas.

Sin embargo, existen varias diferencias. La más importante es qué es lo que se está calculando.

• El algoritmo EM calcula estimaciones puntuales de la distribución posterior de aquellas variables aleatorias que pueden clasificarse como "parámetros", pero solo estimaciones de las distribuciones posteriores reales de las variables latentes (al menos en el "EM suave", y a menudo solo cuando las variables latentes son discretas). Las estimaciones puntuales calculadas son las modas de estos parámetros; no se dispone de ninguna otra información.
• Por otro lado, VB calcula estimaciones de la distribución posterior real de todas las variables, tanto parámetros como variables latentes. Cuando se necesitan obtener estimaciones puntuales, generalmente se utiliza la media en lugar de la moda, como es habitual en la inferencia bayesiana. En consecuencia, los parámetros calculados en VB no tienen la misma significancia que los de EM. EM calcula los valores óptimos de los parámetros de la propia red bayesiana. VB calcula los valores óptimos de los parámetros de las distribuciones utilizadas para aproximar los parámetros y las variables latentes de la red bayesiana. Por ejemplo, un modelo de mezcla gaussiana típico tendrá parámetros para la media y la varianza de cada uno de los componentes de la mezcla. EM estimaría directamente los valores óptimos para estos parámetros. VB, sin embargo, primero ajustaría una distribución a estos parámetros —normalmente en forma de una distribución a priori , por ejemplo, una distribución gamma inversa escalada normal— y luego calcularía los valores para los parámetros de esta distribución a priori, es decir, esencialmente hiperparámetros . En este caso, VB calcularía estimaciones óptimas de los cuatro parámetros de la distribución gamma inversa escalada normal que describe la distribución conjunta de la media y la varianza del componente.