Métodos bayesianos variacionales

Métodos matemáticos utilizados en la inferencia bayesiana y el aprendizaje automático

Los métodos bayesianos variacionales son una familia de técnicas para aproximar integrales intratables que surgen en la inferencia bayesiana y el aprendizaje automático . Se utilizan normalmente en modelos estadísticos complejos que constan de variables observadas (normalmente denominadas "datos"), así como parámetros desconocidos y variables latentes , con diversos tipos de relaciones entre los tres tipos de variables aleatorias , como podría describirse mediante un modelo gráfico . Como es habitual en la inferencia bayesiana, los parámetros y las variables latentes se agrupan como "variables no observadas". Los métodos bayesianos variacionales se utilizan principalmente con dos fines:

1. Proporcionar una aproximación analítica a la probabilidad posterior de las variables no observadas, con el fin de realizar inferencia estadística sobre estas variables.
2. Para derivar un límite inferior para la probabilidad marginal (a veces llamada evidencia ) de los datos observados (es decir, la probabilidad marginal de los datos dado el modelo, con marginalización realizada sobre variables no observadas). Esto se utiliza típicamente para realizar la selección de modelos , la idea general es que una mayor probabilidad marginal para un modelo dado indica un mejor ajuste de los datos por ese modelo y, por lo tanto, una mayor probabilidad de que el modelo en cuestión sea el que generó los datos. (Véase también el artículo sobre el factor de Bayes ).

En el primer caso (el de aproximar una probabilidad posterior), el método Bayes variacional es una alternativa a los métodos de muestreo de Monte Carlo (en particular, los métodos de Monte Carlo de cadena de Markov , como el muestreo de Gibbs ) para adoptar un enfoque totalmente bayesiano para la inferencia estadística sobre distribuciones complejas que son difíciles de evaluar directamente o de muestrear . En particular, mientras que las técnicas de Monte Carlo proporcionan una aproximación numérica a la probabilidad posterior exacta utilizando un conjunto de muestras, el método Bayes variacional proporciona una solución analítica exacta y localmente óptima para una aproximación de la probabilidad posterior.

El Bayes variacional puede considerarse una extensión del algoritmo de expectativa-maximización (EM) desde la estimación de máxima verosimilitud (ML) o máxima a posteriori (MAP) del valor más probable de cada parámetro hasta la estimación completamente bayesiana, que calcula (una aproximación a) la distribución posterior completa de los parámetros y las variables latentes. Al igual que en el EM, encuentra un conjunto de valores óptimos de los parámetros y tiene la misma estructura alternante que el EM, basada en un conjunto de ecuaciones interconectadas (mutuamente dependientes) que no se pueden resolver analíticamente.

Para muchas aplicaciones, el método Bayes variacional produce soluciones de precisión comparable al muestreo de Gibbs a mayor velocidad. Sin embargo, derivar el conjunto de ecuaciones utilizadas para actualizar los parámetros de manera iterativa a menudo requiere una gran cantidad de trabajo en comparación con la derivación de ecuaciones de muestreo de Gibbs comparables. Este es el caso incluso para muchos modelos que son conceptualmente bastante simples, como se demuestra a continuación en el caso de un modelo básico no jerárquico con solo dos parámetros y sin variables latentes.

Derivación matemática

Problema

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

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

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

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

Divergencia KL

El tipo más común de Bayes variacional utiliza la divergencia de Kullback-Leibler (divergencia KL) de Q a partir de P como la opción de función de disimilitud. Esta opció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} )}}.}$

Obsérvese que Q y P son inversos a lo que cabría esperar. Este uso de la divergencia KL invertida es conceptualmente similar al algoritmo de maximización de expectativas . (El uso de la divergencia KL de la otra forma produce el algoritmo de propagación de expectativas ).

Dificultad

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 marginalización sobre el cálculo en el denominador suele ser intratable, porque, por ejemplo, el espacio de búsqueda de es combinatoriamente grande. Por lo tanto, buscamos una aproximación, utilizando . ${\displaystyle \mathbf {Z} }$ ${\displaystyle P(\mathbf {X} )}$ ${\displaystyle \mathbf {Z} }$ ${\displaystyle Q(\mathbf {Z} )\approx P(\mathbf {Z} \mid \mathbf {X} )}$

Límite inferior de evidencia

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

${\displaystyle 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]}$

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

${\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 ), se puede escribir 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 \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)}$

Como la evidencia logarítmica se fija con respecto a , maximizar el término final minimiza la divergencia KL de desde . Con la elección adecuada de , se vuelve factible de calcular y maximizar. Por lo tanto, tenemos tanto una aproximación analítica para el posterior como un límite inferior para la evidencia logarítmica (ya que la divergencia KL no es negativa). ${\displaystyle \log P(\mathbf {X} )}$ ${\displaystyle Q}$ ${\displaystyle {\mathcal {L}}(Q)}$ ${\displaystyle Q}$ ${\displaystyle P}$ ${\displaystyle Q}$ ${\displaystyle {\mathcal {L}}(Q)}$ ${\displaystyle Q}$ ${\displaystyle P(\mathbf {Z} \mid \mathbf {X} )}$ ${\displaystyle {\mathcal {L}}(Q)}$ ${\displaystyle \log P(\mathbf {X} )}$

El límite inferior se conoce como energía libre variacional (negativa) en analogía con la energía libre termodinámica porque también se puede expresar como una energía negativa más la entropía de . El término 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) en la evidencia logarítmica de los datos. ${\displaystyle {\mathcal {L}}(Q)}$ ${\displaystyle \operatorname {E} _{Q}[\log P(\mathbf {Z} ,\mathbf {X} )]}$ ${\displaystyle Q}$ ${\displaystyle {\mathcal {L}}(Q)}$

Pruebas

Por el teorema de Pitágoras generalizado de la divergencia de Bregman , del cual la divergencia KL es un caso especial, se puede demostrar que: ^{[1] [2]}

Teorema de Pitágoras generalizado para la divergencia de Bregman ^[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}}}$

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

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

En este caso, el minimizador global se puede encontrar de la siguiente manera: ^[1] ${\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}),}$ ${\displaystyle \mathbf {Z} =\{\mathbf {Z_{1}} ,\mathbf {Z_{2}} \},}$

${\displaystyle 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),}$

en la que la constante normalizadora es:

${\displaystyle \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).}$

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

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

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

Aproximación del campo medio

Generalmente se supone que la distribución variacional se factoriza sobre alguna partición de las variables latentes, es decir, para alguna partición de las variables latentes en , ${\displaystyle Q(\mathbf {Z} )}$ ${\displaystyle \mathbf {Z} }$ ${\displaystyle \mathbf {Z} _{1}\dots \mathbf {Z} _{M}}$

${\displaystyle Q(\mathbf {Z} )=\prod _{i=1}^{M}q_{i}(\mathbf {Z} _{i}\mid \mathbf {X} )}$

Se puede demostrar utilizando el cálculo de variaciones (de ahí el nombre "Bayes variacional") que la "mejor" distribución para cada uno de los factores (en términos de la distribución que minimiza la divergencia KL, como se describió anteriormente) satisface: ^[3] ${\displaystyle q_{j}^{*}}$ ${\displaystyle q_{j}}$

${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )={\frac {e^{\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}}{\int e^{\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}\,d\mathbf {Z} _{j}}}}$

donde es la expectativa del logaritmo de la probabilidad conjunta de los datos y las variables latentes, tomada con respecto a todas las variables que no están en la partición: consulte el Lema 4.1 de ^[4] para obtener una derivación de la distribución . ${\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}$ ${\displaystyle q^{*}}$ ${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )}$

En la práctica, normalmente trabajamos en términos de logaritmos, es decir:

${\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}}$

La constante en la expresión anterior está relacionada con la constante de normalización (el denominador en la expresión anterior para ) y generalmente se restablece mediante inspección, ya que el resto de la expresión generalmente se puede reconocer como un tipo conocido de distribución (por ejemplo , gaussiana , gamma , etc.). ${\displaystyle q_{j}^{*}}$

Utilizando las propiedades de las expectativas, la expresión puede simplificarse usualmente en una función de los hiperparámetros fijos de las distribuciones previas sobre las variables latentes y de las expectativas (y a veces momentos superiores como la varianza ) de las variables latentes que no están en la partición actual (es decir, variables latentes no incluidas en ). Esto crea dependencias circulares entre los parámetros de las distribuciones sobre las variables en una partición y las expectativas de las variables en las otras particiones. Esto naturalmente sugiere un algoritmo iterativo , muy parecido a EM (el algoritmo de expectativa-maximización ), en el que las expectativas (y posiblemente los momentos superiores) de las variables latentes se inicializan de alguna manera (quizás aleatoriamente), y luego los parámetros de cada distribución se calculan a su vez utilizando los valores actuales de las expectativas, después de lo cual la expectativa de la distribución recién calculada se establece apropiadamente de acuerdo con los parámetros calculados. Se garantiza que un algoritmo de este tipo convergerá . ^[5] ${\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]}$ ${\displaystyle \mathbf {Z} _{j}}$

En otras palabras, para cada una de las particiones de variables, simplificando la expresión de la distribución sobre las variables de la partición y examinando la dependencia funcional de la distribución con respecto a las variables en cuestión, normalmente se puede determinar la familia de la distribución (que a su vez determina el valor de la constante). La fórmula para los parámetros de la distribución se expresará en términos de los hiperparámetros de las distribuciones anteriores (que son constantes conocidas), pero también en términos de expectativas de funciones de variables en otras particiones. Normalmente, estas expectativas se pueden simplificar en funciones de expectativas de las propias variables (es decir, las medias ); a veces también aparecen expectativas de variables al cuadrado (que pueden estar relacionadas con la varianza de las variables) o expectativas de potencias superiores (es decir, momentos superiores ). En la mayoría de los casos, las distribuciones de las otras variables serán de familias conocidas, y se pueden buscar las fórmulas para las expectativas relevantes. Sin embargo, esas fórmulas dependen de los parámetros de esas distribuciones, que a su vez dependen de las expectativas sobre otras variables. El resultado es que las fórmulas para los parámetros de las distribuciones de cada variable pueden expresarse como una serie de ecuaciones con dependencias mutuas y no lineales entre las variables. Por lo general, no es posible resolver este sistema de ecuaciones directamente. Sin embargo, como se describió anteriormente, las dependencias sugieren un algoritmo iterativo simple, cuya convergencia está garantizada en la mayoría de los casos. Un ejemplo aclarará este proceso.

Una fórmula de dualidad para la inferencia variacional

Ilustración gráfica del algoritmo de inferencia variacional de ascenso de coordenadas mediante la fórmula de dualidad ^[4]

El siguiente teorema se denomina fórmula de dualidad para la inferencia variacional. ^[4] Explica algunas propiedades importantes de las distribuciones variacionales utilizadas en los métodos de Bayes variacional.

Teorema Consideremos dos espacios de probabilidad y con . Supongamos que existe una medida de probabilidad dominante común tal que y . Sea cualquier variable aleatoria de valor real en que satisfaga . Entonces se cumple la siguiente igualdad ${\displaystyle (\Theta ,{\mathcal {F}},P)}$ ${\displaystyle (\Theta ,{\mathcal {F}},Q)}$ ${\displaystyle Q\ll P}$ ${\displaystyle \lambda }$ ${\displaystyle P\ll \lambda }$ ${\displaystyle Q\ll \lambda }$ ${\displaystyle h}$ ${\displaystyle (\Theta ,{\mathcal {F}},P)}$ ${\displaystyle h\in L_{1}(P)}$

${\displaystyle \log E_{P}[\exp h]={\text{sup}}_{Q\ll P}\{E_{Q}[h]-D_{\text{KL}}(Q\parallel P)\}.}$

Además, el supremo del lado derecho se alcanza si y sólo si se cumple

${\displaystyle {\frac {q(\theta )}{p(\theta )}}={\frac {\exp h(\theta )}{E_{P}[\exp h]}},}$

casi con seguridad con respecto a la medida de probabilidad , donde y denotan las derivadas de Radon-Nikodym de las medidas de probabilidad y con respecto a , respectivamente. ${\displaystyle Q}$ ${\displaystyle p(\theta )=dP/d\lambda }$ ${\displaystyle q(\theta )=dQ/d\lambda }$ ${\displaystyle P}$ ${\displaystyle Q}$ ${\displaystyle \lambda }$

Un ejemplo básico

Consideremos un modelo bayesiano no jerárquico simple que consiste en un conjunto de observaciones iid de una distribución gaussiana , con media y varianza desconocidas . ^[6] A continuación, trabajamos este modelo en gran detalle para ilustrar el funcionamiento del método Bayes variacional.

Para facilitar las matemáticas, en el siguiente ejemplo trabajamos en términos de precisión (es decir, el recíproco de la varianza o, en una ecuación gaussiana multivariante, el inverso de la matriz de covarianza ) en lugar de la varianza en sí. (Desde un punto de vista teórico, precisión y varianza son equivalentes, ya que existe una correspondencia uno a uno entre las dos).

El modelo matemático

Colocamos distribuciones previas conjugadas sobre la media y precisión desconocidas , 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 \mu }$ ${\displaystyle \tau }$

${\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 y 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 y . ${\displaystyle \mu _{0},\lambda _{0},a_{0}}$ ${\displaystyle b_{0}}$ ${\displaystyle \mu }$ ${\displaystyle \tau }$

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

La probabilidad conjunta

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 están los factores individuales

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

Aproximación factorizada

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

Derivación deq ( μ )

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, , y se refieren a valores que son constantes con respecto a . Nótese que el término no es una función de y tendrá el mismo valor independientemente del valor de . 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. ${\displaystyle C}$ ${\displaystyle C_{2}}$ ${\displaystyle C_{3}}$ ${\displaystyle \mu }$ ${\displaystyle \operatorname {E} _{\tau }[\ln p(\tau )]}$ ${\displaystyle \mu }$ ${\displaystyle \mu }$

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

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

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

Derivación deq(τ)

La derivación de es similar a la anterior, aunque omitimos algunos detalles por razones de brevedad. ${\displaystyle q_{\tau }^{*}(\tau )}$

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

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

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

Algoritmo para calcular los parámetros

Resumamos 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 sobre una de las variables dependen de las expectativas tomadas con respecto a la otra variable. Podemos ampliar las expectativas, utilizando las fórmulas estándar para las expectativas de 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 requiere más trabajo: ${\displaystyle b_{N}}$

${\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 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 y . Esto naturalmente sugiere un algoritmo similar a EM : ${\displaystyle \lambda _{N}}$ ${\displaystyle b_{N}}$

1. Calcular y utilizar estos valores para calcular y ${\displaystyle \sum _{n=1}^{N}x_{n}}$ ${\displaystyle \sum _{n=1}^{N}x_{n}^{2}.}$ ${\displaystyle \mu _{N}}$ ${\displaystyle a_{N}.}$
2. Inicializar con algún valor arbitrario. ${\displaystyle \lambda _{N}}$
3. Utilice el valor actual de junto con los valores conocidos de los otros parámetros para calcular . ${\displaystyle \lambda _{N},}$ ${\displaystyle b_{N}}$
4. Utilice el valor actual de junto con los valores conocidos de los otros parámetros para calcular . ${\displaystyle b_{N},}$ ${\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).

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

Se puede demostrar que se garantiza que este algoritmo convergerá a un máximo local.

Nótese también que las distribuciones posteriores tienen la misma forma que las distribuciones previas correspondientes. No asumimos esto; la única suposición que hicimos fue que las distribuciones se factorizan, y la forma de las distribuciones se deduce naturalmente. Resulta (ver abajo) 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 sean miembros de la familia exponencial , que es el caso de la mayoría de las distribuciones estándar.

Más discusión

Receta paso a paso

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. Describa la red con un modelo gráfico , identificando las variables observadas (datos) y las variables no observadas ( parámetros y variables latentes ) y sus distribuciones de probabilidad condicional . El método Bayes variacional construirá entonces una aproximación a la probabilidad posterior . La aproximación tiene la propiedad básica de que es una distribución factorizada, es decir, un producto de dos o más distribuciones independientes sobre subconjuntos disjuntos de las variables no observadas. ${\displaystyle \mathbf {X} }$ ${\displaystyle {\boldsymbol {\Theta }}}$ ${\displaystyle \mathbf {Z} }$ ${\displaystyle p(\mathbf {Z} ,{\boldsymbol {\Theta }}\mid \mathbf {X} )}$
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 de Bayes variacional sea intratable. Por lo general, la primera división es para separar los parámetros y las variables latentes; a menudo, esto es suficiente por sí solo para producir un resultado manejable. Suponga que las particiones se denominan . ${\displaystyle \mathbf {Z} _{1},\ldots ,\mathbf {Z} _{M}}$
3. Para una partición dada , escriba la fórmula para la mejor distribución aproximada utilizando la ecuación básica . ${\displaystyle \mathbf {Z} _{j}}$ ${\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )}$ ${\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. Se pueden ignorar todas las distribuciones condicionales de componentes que no involucren ninguna de las variables ; se incorporarán al término constante. ${\displaystyle \mathbf {Z} _{j}}$
5. Simplifique la fórmula y aplique 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 (por ejemplo, primer o segundo momento bruto , expectativa de un logaritmo, etc.). Para que el procedimiento de Bayes variacional funcione bien, estas expectativas generalmente deberían 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 expectativa son constantes con respecto a las variables en la partición actual. ${\displaystyle \mathbf {Z} _{j}}$
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, la exponenciación de la fórmula genera la función de densidad de probabilidad (PDF) de la distribución (o al menos, algo proporcional a ella, con una constante de normalización desconocida ). Para que el método general sea manejable, debería 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 PDF de una distribución conocida. Cuando esto se puede hacer, la constante de normalización se puede restablecer por definición y se pueden derivar ecuaciones para los parámetros de la distribución conocida 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 PDF puede ponerse en una forma que permite la 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 mutuamente vinculadas que especifican los valores óptimos de todos los parámetros.
9. Luego se aplica un procedimiento de tipo expectativa-maximización (EM), eligiendo un valor inicial para cada parámetro y repitiendo una serie de pasos, donde en cada paso recorremos las ecuaciones, actualizando cada parámetro a su vez. Se garantiza que esto convergerá.

Puntos más importantes

Debido a todas las manipulaciones matemáticas involucradas, es fácil perder de vista el panorama general. Los aspectos importantes son:

1. La idea del Bayes 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 bayesianos , 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, de la expectativa-maximización (EM) u otros métodos de máxima verosimilitud , ambos tipos de variables no observadas (es decir, parámetros y variables latentes) se tratan de la misma manera, es decir, como variables aleatorias . Las estimaciones para las variables se pueden derivar de las formas bayesianas estándar, por ejemplo, calculando la media de la distribución para obtener una estimación puntual única o derivando un intervalo creíble , la región de mayor densidad, etc.
2. "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 distribución posterior verdadera, sino una aproximación a ella; en particular, generalmente concordará bastante en los momentos más bajos 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 forman los factores y (2) fórmulas mutuamente dependientes para los parámetros de estas distribuciones. Los valores reales de estos parámetros se calculan numéricamente, a través de un procedimiento iterativo alterno muy parecido al EM.

En comparación con la expectativa-maximización (EM)

El método Bayes variacional (VB) suele compararse con el de expectativa-maximización (EM). El procedimiento numérico real es bastante similar, ya que ambos son procedimientos iterativos alternativos que convergen sucesivamente en valores óptimos de parámetros. Los pasos iniciales para derivar los respectivos procedimientos también son vagamente similares, ya que ambos comienzan con fórmulas para densidades de probabilidad y ambos implican cantidades significativas de manipulaciones matemáticas.

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

• El método EM calcula estimaciones puntuales de la distribución posterior de aquellas variables aleatorias que pueden clasificarse como "parámetros", pero sólo estimaciones de las distribuciones posteriores reales de las variables latentes (al menos en el método EM "suave", y a menudo sólo cuando las variables latentes son discretas). Las estimaciones puntuales calculadas son las modas de estos parámetros; no hay otra información disponible.
• VB, por otro lado, calcula estimaciones de la distribución posterior real de todas las variables, tanto parámetros como variables latentes. Cuando se necesitan derivar estimaciones puntuales, generalmente se utiliza la media en lugar de la moda, como es normal en la inferencia bayesiana. Concomitante con esto, los parámetros calculados en VB no tienen el mismo significado que los de EM. EM calcula valores óptimos de los parámetros de la propia red de Bayes. VB calcula valores óptimos de los parámetros de las distribuciones utilizadas para aproximar los parámetros y las variables latentes de la red de Bayes. 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 previa , por ejemplo, una distribución gamma inversa de escala normal ) y luego calcularía valores para los parámetros de esta distribución previa, es decir, esencialmente hiperparámetros . En este caso, VB calcularía estimaciones óptimas de los cuatro parámetros de la distribución gamma inversa de escala normal que describe la distribución conjunta de la media y la varianza del componente.

Un ejemplo más complejo

Modelo de mezcla gaussiana bayesiana que utiliza notación de placas . Los cuadrados más pequeños indican parámetros fijos; los círculos más grandes indican variables aleatorias. Las formas rellenas indican valores conocidos. La indicación [K] significa un vector de tamaño K ; [ D , D ] significa una matriz de tamaño D × D ; K solo significa una variable categórica con K resultados. La línea ondulada que sale de z y termina en una barra transversal indica un interruptor : el valor de esta variable selecciona, para las otras variables entrantes, qué valor utilizar de la matriz de tamaño K de valores posibles.

Imaginemos un modelo de mezcla gaussiana bayesiana descrito de la siguiente manera: ^[3]

${\displaystyle {\begin{aligned}\mathbf {\pi } &\sim \operatorname {SymDir} (K,\alpha _{0})\\\mathbf {\Lambda } _{i=1\dots K}&\sim {\mathcal {W}}(\mathbf {W} _{0},\nu _{0})\\\mathbf {\mu } _{i=1\dots K}&\sim {\mathcal {N}}(\mathbf {\mu } _{0},(\beta _{0}\mathbf {\Lambda } _{i})^{-1})\\\mathbf {z} [i=1\dots N]&\sim \operatorname {Mult} (1,\mathbf {\pi } )\\\mathbf {x} _{i=1\dots N}&\sim {\mathcal {N}}(\mathbf {\mu } _{z_{i}},{\mathbf {\Lambda } _{z_{i}}}^{-1})\\K&={\text{number of mixing components}}\\N&={\text{number of data points}}\end{aligned}}}$

Nota:

• SymDir() es la distribución de Dirichlet simétrica de dimensión , con el hiperparámetro para cada componente establecido en . La distribución de Dirichlet es la distribución conjugada previa de la distribución categórica o distribución multinomial . ${\displaystyle K}$ ${\displaystyle \alpha _{0}}$
• ${\displaystyle {\mathcal {W}}()}$es la distribución Wishart , que es la anterior conjugada de la matriz de precisión ( matriz de covarianza inversa ) para una distribución gaussiana multivariada .
• Mult() es una distribución multinomial sobre una única observación (equivalente a una distribución categórica ). El espacio de estados es una representación "uno de K", es decir, un vector dimensional en el que uno de los elementos es 1 (que especifica la identidad de la observación) y todos los demás elementos son 0. ${\displaystyle K}$
• ${\displaystyle {\mathcal {N}}()}$es la distribución gaussiana , en este caso específicamente la distribución gaussiana multivariada .

La interpretación de las variables anteriores es la siguiente:

• ${\displaystyle \mathbf {X} =\{\mathbf {x} _{1},\dots ,\mathbf {x} _{N}\}}$es el conjunto de puntos de datos, cada uno de los cuales es un vector dimensional distribuido según una distribución gaussiana multivariada . ${\displaystyle N}$ ${\displaystyle D}$
• ${\displaystyle \mathbf {Z} =\{\mathbf {z} _{1},\dots ,\mathbf {z} _{N}\}}$es un conjunto de variables latentes, una por punto de datos, que especifica a qué componente de la mezcla pertenece el punto de datos correspondiente, utilizando una representación vectorial "uno de K" con componentes para , como se describió anteriormente. ${\displaystyle z_{nk}}$ ${\displaystyle k=1\dots K}$
• ${\displaystyle \mathbf {\pi } }$son las proporciones de mezcla de los componentes de la mezcla. ${\displaystyle K}$
• ${\displaystyle \mathbf {\mu } _{i=1\dots K}}$y especificar los parámetros ( media y precisión ) asociados a cada componente de la mezcla. ${\displaystyle \mathbf {\Lambda } _{i=1\dots K}}$

La probabilidad conjunta de todas las variables se puede reescribir como

${\displaystyle p(\mathbf {X} ,\mathbf {Z} ,\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )=p(\mathbf {X} \mid \mathbf {Z} ,\mathbf {\mu } ,\mathbf {\Lambda } )p(\mathbf {Z} \mid \mathbf {\pi } )p(\mathbf {\pi } )p(\mathbf {\mu } \mid \mathbf {\Lambda } )p(\mathbf {\Lambda } )}$

donde están los factores individuales

${\displaystyle {\begin{aligned}p(\mathbf {X} \mid \mathbf {Z} ,\mathbf {\mu } ,\mathbf {\Lambda } )&=\prod _{n=1}^{N}\prod _{k=1}^{K}{\mathcal {N}}(\mathbf {x} _{n}\mid \mathbf {\mu } _{k},\mathbf {\Lambda } _{k}^{-1})^{z_{nk}}\\p(\mathbf {Z} \mid \mathbf {\pi } )&=\prod _{n=1}^{N}\prod _{k=1}^{K}\pi _{k}^{z_{nk}}\\p(\mathbf {\pi } )&={\frac {\Gamma (K\alpha _{0})}{\Gamma (\alpha _{0})^{K}}}\prod _{k=1}^{K}\pi _{k}^{\alpha _{0}-1}\\p(\mathbf {\mu } \mid \mathbf {\Lambda } )&=\prod _{k=1}^{K}{\mathcal {N}}(\mathbf {\mu } _{k}\mid \mathbf {\mu } _{0},(\beta _{0}\mathbf {\Lambda } _{k})^{-1})\\p(\mathbf {\Lambda } )&=\prod _{k=1}^{K}{\mathcal {W}}(\mathbf {\Lambda } _{k}\mid \mathbf {W} _{0},\nu _{0})\end{aligned}}}$

dónde

${\displaystyle {\begin{aligned}{\mathcal {N}}(\mathbf {x} \mid \mathbf {\mu } ,\mathbf {\Sigma } )&={\frac {1}{(2\pi )^{D/2}}}{\frac {1}{|\mathbf {\Sigma } |^{1/2}}}\exp \left\{-{\frac {1}{2}}(\mathbf {x} -\mathbf {\mu } )^{\rm {T}}\mathbf {\Sigma } ^{-1}(\mathbf {x} -\mathbf {\mu } )\right\}\\{\mathcal {W}}(\mathbf {\Lambda } \mid \mathbf {W} ,\nu )&=B(\mathbf {W} ,\nu )|\mathbf {\Lambda } |^{(\nu -D-1)/2}\exp \left(-{\frac {1}{2}}\operatorname {Tr} (\mathbf {W} ^{-1}\mathbf {\Lambda } )\right)\\B(\mathbf {W} ,\nu )&=|\mathbf {W} |^{-\nu /2}\left\{2^{\nu D/2}\pi ^{D(D-1)/4}\prod _{i=1}^{D}\Gamma \left({\frac {\nu +1-i}{2}}\right)\right\}^{-1}\\D&={\text{dimensionality of each data point}}\end{aligned}}}$

Supongamos que . ${\displaystyle q(\mathbf {Z} ,\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )=q(\mathbf {Z} )q(\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )}$

Entonces ^[3]

${\displaystyle {\begin{aligned}\ln q^{*}(\mathbf {Z} )&=\operatorname {E} _{\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } }[\ln p(\mathbf {X} ,\mathbf {Z} ,\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )]+{\text{constant}}\\&=\operatorname {E} _{\mathbf {\pi } }[\ln p(\mathbf {Z} \mid \mathbf {\pi } )]+\operatorname {E} _{\mathbf {\mu } ,\mathbf {\Lambda } }[\ln p(\mathbf {X} \mid \mathbf {Z} ,\mathbf {\mu } ,\mathbf {\Lambda } )]+{\text{constant}}\\&=\sum _{n=1}^{N}\sum _{k=1}^{K}z_{nk}\ln \rho _{nk}+{\text{constant}}\end{aligned}}}$

donde hemos definido

${\displaystyle \ln \rho _{nk}=\operatorname {E} [\ln \pi _{k}]+{\frac {1}{2}}\operatorname {E} [\ln |\mathbf {\Lambda } _{k}|]-{\frac {D}{2}}\ln(2\pi )-{\frac {1}{2}}\operatorname {E} _{\mathbf {\mu } _{k},\mathbf {\Lambda } _{k}}[(\mathbf {x} _{n}-\mathbf {\mu } _{k})^{\rm {T}}\mathbf {\Lambda } _{k}(\mathbf {x} _{n}-\mathbf {\mu } _{k})]}$

Exponenciando ambos lados de la fórmula para obtener rendimientos ${\displaystyle \ln q^{*}(\mathbf {Z} )}$

${\displaystyle q^{*}(\mathbf {Z} )\propto \prod _{n=1}^{N}\prod _{k=1}^{K}\rho _{nk}^{z_{nk}}}$

Requerir que esto se normalice termina requiriendo que la suma sea 1 para todos los valores de , lo que da como resultado ${\displaystyle \rho _{nk}}$ ${\displaystyle k}$

${\displaystyle q^{*}(\mathbf {Z} )=\prod _{n=1}^{N}\prod _{k=1}^{K}r_{nk}^{z_{nk}}}$

dónde

${\displaystyle r_{nk}={\frac {\rho _{nk}}{\sum _{j=1}^{K}\rho _{nj}}}}$

En otras palabras, es un producto de distribuciones multinomiales de observación única y factores sobre cada individuo , que se distribuye como una distribución multinomial de observación única con parámetros para . ${\displaystyle q^{*}(\mathbf {Z} )}$ ${\displaystyle \mathbf {z} _{n}}$ ${\displaystyle r_{nk}}$ ${\displaystyle k=1\dots K}$

Además, observamos que

${\displaystyle \operatorname {E} [z_{nk}]=r_{nk}\,}$

que es un resultado estándar para distribuciones categóricas.

Ahora, considerando el factor , note que se incluye automáticamente en debido a la estructura del modelo gráfico que define nuestro modelo de mezcla gaussiana, que se especifica arriba. ${\displaystyle q(\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )}$ ${\displaystyle q(\mathbf {\pi } )\prod _{k=1}^{K}q(\mathbf {\mu } _{k},\mathbf {\Lambda } _{k})}$

Entonces,

${\displaystyle {\begin{aligned}\ln q^{*}(\mathbf {\pi } )&=\ln p(\mathbf {\pi } )+\operatorname {E} _{\mathbf {Z} }[\ln p(\mathbf {Z} \mid \mathbf {\pi } )]+{\text{constant}}\\&=(\alpha _{0}-1)\sum _{k=1}^{K}\ln \pi _{k}+\sum _{n=1}^{N}\sum _{k=1}^{K}r_{nk}\ln \pi _{k}+{\text{constant}}\end{aligned}}}$

Tomando la exponencial de ambos lados, la reconocemos como una distribución de Dirichlet ${\displaystyle q^{*}(\mathbf {\pi } )}$

${\displaystyle q^{*}(\mathbf {\pi } )\sim \operatorname {Dir} (\mathbf {\alpha } )\,}$

dónde

${\displaystyle \alpha _{k}=\alpha _{0}+N_{k}\,}$

dónde

${\displaystyle N_{k}=\sum _{n=1}^{N}r_{nk}\,}$

Finalmente

${\displaystyle \ln q^{*}(\mathbf {\mu } _{k},\mathbf {\Lambda } _{k})=\ln p(\mathbf {\mu } _{k},\mathbf {\Lambda } _{k})+\sum _{n=1}^{N}\operatorname {E} [z_{nk}]\ln {\mathcal {N}}(\mathbf {x} _{n}\mid \mathbf {\mu } _{k},\mathbf {\Lambda } _{k}^{-1})+{\text{constant}}}$

Agrupando y leyendo términos que involucran y , el resultado es una distribución Gaussiana-Wishart dada por ${\displaystyle \mathbf {\mu } _{k}}$ ${\displaystyle \mathbf {\Lambda } _{k}}$

${\displaystyle q^{*}(\mathbf {\mu } _{k},\mathbf {\Lambda } _{k})={\mathcal {N}}(\mathbf {\mu } _{k}\mid \mathbf {m} _{k},(\beta _{k}\mathbf {\Lambda } _{k})^{-1}){\mathcal {W}}(\mathbf {\Lambda } _{k}\mid \mathbf {W} _{k},\nu _{k})}$

Dadas las definiciones

${\displaystyle {\begin{aligned}\beta _{k}&=\beta _{0}+N_{k}\\\mathbf {m} _{k}&={\frac {1}{\beta _{k}}}(\beta _{0}\mathbf {\mu } _{0}+N_{k}{\bar {\mathbf {x} }}_{k})\\\mathbf {W} _{k}^{-1}&=\mathbf {W} _{0}^{-1}+N_{k}\mathbf {S} _{k}+{\frac {\beta _{0}N_{k}}{\beta _{0}+N_{k}}}({\bar {\mathbf {x} }}_{k}-\mathbf {\mu } _{0})({\bar {\mathbf {x} }}_{k}-\mathbf {\mu } _{0})^{\rm {T}}\\\nu _{k}&=\nu _{0}+N_{k}\\N_{k}&=\sum _{n=1}^{N}r_{nk}\\{\bar {\mathbf {x} }}_{k}&={\frac {1}{N_{k}}}\sum _{n=1}^{N}r_{nk}\mathbf {x} _{n}\\\mathbf {S} _{k}&={\frac {1}{N_{k}}}\sum _{n=1}^{N}r_{nk}(\mathbf {x} _{n}-{\bar {\mathbf {x} }}_{k})(\mathbf {x} _{n}-{\bar {\mathbf {x} }}_{k})^{\rm {T}}\end{aligned}}}$

Por último, observe que estas funciones requieren los valores de , que hacen uso de , que se define a su vez en función de , , y . Ahora que hemos determinado las distribuciones sobre las que se toman estas expectativas, podemos derivar fórmulas para ellas: ${\displaystyle r_{nk}}$ ${\displaystyle \rho _{nk}}$ ${\displaystyle \operatorname {E} [\ln \pi _{k}]}$ ${\displaystyle \operatorname {E} [\ln |\mathbf {\Lambda } _{k}|]}$ ${\displaystyle \operatorname {E} _{\mathbf {\mu } _{k},\mathbf {\Lambda } _{k}}[(\mathbf {x} _{n}-\mathbf {\mu } _{k})^{\rm {T}}\mathbf {\Lambda } _{k}(\mathbf {x} _{n}-\mathbf {\mu } _{k})]}$

${\displaystyle {\begin{aligned}\operatorname {E} _{\mathbf {\mu } _{k},\mathbf {\Lambda } _{k}}[(\mathbf {x} _{n}-\mathbf {\mu } _{k})^{\rm {T}}\mathbf {\Lambda } _{k}(\mathbf {x} _{n}-\mathbf {\mu } _{k})]&=D\beta _{k}^{-1}+\nu _{k}(\mathbf {x} _{n}-\mathbf {m} _{k})^{\rm {T}}\mathbf {W} _{k}(\mathbf {x} _{n}-\mathbf {m} _{k})\\\ln {\widetilde {\Lambda }}_{k}&\equiv \operatorname {E} [\ln |\mathbf {\Lambda } _{k}|]=\sum _{i=1}^{D}\psi \left({\frac {\nu _{k}+1-i}{2}}\right)+D\ln 2+\ln |\mathbf {W} _{k}|\\\ln {\widetilde {\pi }}_{k}&\equiv \operatorname {E} \left[\ln |\pi _{k}|\right]=\psi (\alpha _{k})-\psi \left(\sum _{i=1}^{K}\alpha _{i}\right)\end{aligned}}}$

Estos resultados conducen a

${\displaystyle r_{nk}\propto {\widetilde {\pi }}_{k}{\widetilde {\Lambda }}_{k}^{1/2}\exp \left\{-{\frac {D}{2\beta _{k}}}-{\frac {\nu _{k}}{2}}(\mathbf {x} _{n}-\mathbf {m} _{k})^{\rm {T}}\mathbf {W} _{k}(\mathbf {x} _{n}-\mathbf {m} _{k})\right\}}$

Estos se pueden convertir de valores proporcionales a absolutos normalizándolos de modo que los valores correspondientes sumen 1. ${\displaystyle k}$

Tenga en cuenta que:

1. Las ecuaciones de actualización de los parámetros , , y de las variables y dependen de las estadísticas , , y , y estas estadísticas a su vez dependen de . ${\displaystyle \beta _{k}}$ ${\displaystyle \mathbf {m} _{k}}$ ${\displaystyle \mathbf {W} _{k}}$ ${\displaystyle \nu _{k}}$ ${\displaystyle \mathbf {\mu } _{k}}$ ${\displaystyle \mathbf {\Lambda } _{k}}$ ${\displaystyle N_{k}}$ ${\displaystyle {\bar {\mathbf {x} }}_{k}}$ ${\displaystyle \mathbf {S} _{k}}$ ${\displaystyle r_{nk}}$
2. Las ecuaciones de actualización de los parámetros de la variable dependen del estadístico , que depende a su vez de . ${\displaystyle \alpha _{1\dots K}}$ ${\displaystyle \mathbf {\pi } }$ ${\displaystyle N_{k}}$ ${\displaystyle r_{nk}}$
3. La ecuación de actualización para tiene una dependencia circular directa de , , y así como una dependencia circular indirecta de , y a través de y . ${\displaystyle r_{nk}}$ ${\displaystyle \beta _{k}}$ ${\displaystyle \mathbf {m} _{k}}$ ${\displaystyle \mathbf {W} _{k}}$ ${\displaystyle \nu _{k}}$ ${\displaystyle \mathbf {W} _{k}}$ ${\displaystyle \nu _{k}}$ ${\displaystyle \alpha _{1\dots K}}$ ${\displaystyle {\widetilde {\pi }}_{k}}$ ${\displaystyle {\widetilde {\Lambda }}_{k}}$

Esto sugiere un procedimiento iterativo que alterna entre dos pasos:

1. Un paso E que calcula el valor utilizando los valores actuales de todos los demás parámetros. ${\displaystyle r_{nk}}$
2. Un paso M que utiliza el nuevo valor de para calcular nuevos valores de todos los demás parámetros. ${\displaystyle r_{nk}}$

Obsérvese que estos pasos se corresponden estrechamente con el algoritmo EM estándar para derivar una solución de máxima verosimilitud o máxima a posteriori (MAP) para los parámetros de un modelo de mezcla gaussiana . Las responsabilidades en el paso E se corresponden estrechamente con las probabilidades posteriores de las variables latentes dados los datos, es decir ; el cálculo de las estadísticas , , y se corresponde estrechamente con el cálculo de las estadísticas de "conteo suave" correspondientes sobre los datos; y el uso de esas estadísticas para calcular nuevos valores de los parámetros se corresponde estrechamente con el uso de conteos suaves para calcular nuevos valores de parámetros en EM normal sobre un modelo de mezcla gaussiana. ${\displaystyle r_{nk}}$ ${\displaystyle p(\mathbf {Z} \mid \mathbf {X} )}$ ${\displaystyle N_{k}}$ ${\displaystyle {\bar {\mathbf {x} }}_{k}}$ ${\displaystyle \mathbf {S} _{k}}$

Distribuciones de familia exponencial

Obsérvese que en el ejemplo anterior, una vez que se supuso que la distribución sobre las variables no observadas se factorizaba en distribuciones sobre los "parámetros" y distribuciones sobre los "datos latentes", la "mejor" distribución derivada para cada variable estaba en la misma familia que la distribución previa correspondiente sobre la variable. Este es un resultado general que se cumple para todas las distribuciones previas derivadas de la familia exponencial .

Véase también

• Paso de mensajes variacional : un algoritmo modular para la inferencia bayesiana variacional.
• Autocodificador variacional : una red neuronal artificial perteneciente a las familias de modelos gráficos probabilísticos y métodos bayesianos variacionales.
• Algoritmo de expectativa-maximización : un enfoque relacionado que corresponde a un caso especial de inferencia bayesiana variacional.
• Filtrado generalizado : un esquema de filtrado variacional para modelos de espacio de estados no lineales.
• Cálculo de variaciones : el campo del análisis matemático que trata de maximizar o minimizar funcionales.
• Discriminación de máxima entropía: este es un marco de inferencia variacional que permite introducir y tener en cuenta restricciones adicionales de gran margen ^[7]

Referencias

1. Tran, Viet Hung (2018). "Inferencia bayesiana variacional de cópula mediante geometría de la información". arXiv : 1803.10998 [cs.IT].
2. Adamčík, Martin (2014). "La geometría de la información de las divergencias de Bregman y algunas aplicaciones en el razonamiento multiexperto". Entropy . 16 (12): 6338–6381. Bibcode :2014Entrp..16.6338A. doi : 10.3390/e16126338 .
3. Nguyen, Duy (15 de agosto de 2023). "UNA INTRODUCCIÓN EN PROFUNDIDAD A LA NOTA VARIACIONAL DE BAYES". doi :10.2139/ssrn.4541076. SSRN  4541076 . Consultado el 15 de agosto de 2023 .
4. Lee, Se Yoon (2021). "Inferencia variacional de ascenso de coordenadas y muestreador de Gibbs: una revisión de la teoría de conjuntos". Comunicaciones en estadística: teoría y métodos . 51 (6): 1–21. arXiv : 2008.01006 . doi :10.1080/03610926.2021.1921214. S2CID  220935477.
5. Boyd, Stephen P.; Vandenberghe, Lieven (2004). Optimización convexa (PDF) . Cambridge University Press. ISBN 978-0-521-83378-3. Recuperado el 15 de octubre de 2011 .
6. Bishop, Christopher M. (2006). "Capítulo 10". Reconocimiento de patrones y aprendizaje automático . Springer. ISBN 978-0-387-31073-2.
7. Sotirios P. Chatzis, “Máquinas de discriminación de máxima entropía con conmutación infinita de Markov”, Proc. 30.ª Conferencia Internacional sobre Aprendizaje Automático (ICML). Journal of Machine Learning Research: Workshop and Conference Proceedings, vol. 28, n.º 3, págs. 729–737, junio de 2013.

Enlaces externos

• El libro de texto en línea: Teoría de la información, inferencia y algoritmos de aprendizaje, de David JC MacKay, ofrece una introducción a los métodos variacionales (p. 422).
• Un tutorial sobre Bayes variacional. Fox, C. y Roberts, S. 2012. Artificial Intelligence Review, doi :10.1007/s10462-011-9236-8.
• Repositorio Variacional-Bayes Un repositorio de artículos de investigación, software y enlaces relacionados con el uso de métodos variacionales para el aprendizaje bayesiano aproximado hasta 2003.
• Algoritmos variacionales para inferencia bayesiana aproximada, de MJ Beal, incluye comparaciones de EM con EM bayesiano variacional y derivaciones de varios modelos, incluidos HMM bayesianos variacionales.
• Puede que valga la pena leer Explicación de alto nivel de la inferencia variacional de Jason Eisner antes de un tratamiento matemáticamente más detallado.
• Inferencia bayesiana variacional por cópula mediante geometría de la información (pdf) por Tran, VH 2018. Este artículo está escrito principalmente para estudiantes. A través de la divergencia de Bregman , el artículo muestra que el bayesiano variacional es simplemente una proyección pitagórica generalizada del modelo verdadero sobre un espacio distribucional correlacionado arbitrariamente (cópula), del cual el espacio independiente es simplemente un caso especial.
• Nota sobre una introducción detallada al método Bayesiano variacional. Nguyen, D. 2023