Métodos bayesianos variacionales

Ajustar Compartir Imprimir Citar

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. Por lo general, se usan en modelos estadísticos complejos que consisten en variables observadas (generalmente denominadas "datos"), así como parámetros desconocidos y variables latentes, con varios tipos de relaciones entre los tres tipos de variables aleatorias, como podría ser descrito por un modelo gráfico. Como es típico 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 para dos propósitos:

Proporcionar una aproximación analítica a la probabilidad posterior de las variables no observadas, con el fin de hacer inferencia estadística sobre estas variables.
Deducir un límite inferior para la probabilidad marginal (a veces denominada evidencia) de los datos observados (es decir, la probabilidad marginal de los datos dado el modelo, con la marginación realizada sobre las variables no observadas). Esto se usa 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 haya sido el que generó los datos. (Consulte también el artículo sobre el factor de Bayes).

En el primer propósito (el de aproximar una probabilidad posterior), el Bayes variacional es una alternativa a los métodos de muestreo de Monte Carlo, en particular, los métodos de Monte Carlo de la 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ícil de evaluar directamente o de muestra. En particular, mientras que las técnicas de Monte Carlo proporcionan una aproximación numérica al posterior exacto utilizando un conjunto de muestras, el Bayes variacional proporciona una solución analítica exacta localmente óptima para una aproximación del posterior.

El bayesiano variacional puede verse como una extensión del algoritmo EM (expectativa-maximización) desde la estimación máxima a posteriori (estimación MAP) del valor único más probable de cada parámetro hasta la estimación totalmente bayesiana que calcula (una aproximación a) toda la distribución posterior de los parámetros y variables latentes. Al igual que en EM, encuentra un conjunto de valores de parámetros óptimos y tiene la misma estructura alterna que EM, basada en un conjunto de ecuaciones entrelazadas (mutuamente dependientes) que no se pueden resolver analíticamente.

Para muchas aplicaciones, Bayes variacional produce soluciones de precisión comparable al muestreo de Gibbs a mayor velocidad. Sin embargo, derivar el conjunto de ecuaciones utilizado para actualizar los parámetros de forma iterativa a menudo requiere una gran cantidad de trabajo en comparación con derivar las 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 ${mathbf {Z}}={Z_{1}puntos Z_{n}}$ dados algunos datos $mathbf{X}$ se aproxima mediante la llamada distribución variacional, ${displaystyle Q(mathbf{Z}):}$ $P({mathbf {Z}}mid {mathbf {X}})aprox. Q({mathbf {Z}}).$

La distribución $Q({mathbf {Z}})$ está restringida a pertenecer a una familia de distribuciones de forma más simple que $P({mathbf {Z}}mid {mathbf {X}})$ (eg una familia de distribuciones Gaussianas), seleccionadas con la intención de hacerlas $Q({mathbf {Z}})$ similares a las verdaderas posteriores, $P({mathbf {Z}}mid {mathbf {X}})$ .

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

Divergencia KL

El tipo más común de Bayes variacional utiliza la divergencia de Kullback-Leibler (KL-divergencia) de Q de P como la elección de la función de disimilitud. Esta elección hace que esta minimización sea manejable. La divergencia KL se define como ${displaystyle D_{mathrm {KL} }(Qparallel P)triangleq sum _{mathbf {Z} }Q(mathbf {Z})log {frac {Q(mathbf {Z})}{P(mathbf {Z} mid mathbf {X})}}.}$

Tenga en cuenta que Q y P están invertidos de lo que cabría esperar. Este uso de la divergencia KL invertida es conceptualmente similar al algoritmo de maximización de expectativas. (Usar la divergencia KL de la otra manera produce el algoritmo de propagación de expectativas).

Dificultad

Las técnicas variacionales se utilizan típicamente 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} '),dmathbf {Z} '}}}$

La marginación ${matemáticas Z}$ de calcular ${displaystyle P(mathbf {X})}$ en el denominador es típicamente intratable porque, por ejemplo, el espacio de búsqueda de ${matemáticas Z}$ es combinatoriamente grande. Por lo tanto, buscamos una aproximación, usando ${displaystyle Q(mathbf {Z})aprox. P(mathbf {Z} mid mathbf {X})}$ .

Límite inferior de evidencia

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 D_{mathrm {KL} }(Qparallel 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})izquierda[log P(mathbf {X})derecha]}$

Como ${displaystyle P(mathbf {X})}$ es una constante con respecto a ${matemáticas Z}$ y ${displaystyle sum_{mathbf {Z} }Q(mathbf {Z})=1}$ como $Q({mathbf {Z}})$ es una distribución, tenemos ${displaystyle D_{mathrm {KL} }(Qparallel P)=sum_{mathbf {Z} }Q(mathbf {Z})left[log Q(mathbf {Z})- log P(mathbf {Z},mathbf {X})right]+log P(mathbf {X})}$

que, de acuerdo con la definición de valor esperado (para una variable aleatoria discreta), se puede escribir de la siguiente manera ${displaystyle D_{mathrm {KL} }(Qparallel 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} }(Qparallel P)-mathbb {E}_{mathbf {Q} }left[log Q( mathbf {Z})-log P(mathbf {Z},mathbf {X})right]=D_{mathrm {KL} }(Qparallel P)+{mathcal {L}}(Q)}$

Como la evidencia logarítmica $log P({mathbf {X}})$ se fija con respecto a $q$ , maximizar el término final ${ matemáticas {L}} (Q)$ minimiza la divergencia KL $q$ de $PAGS$ . Mediante la elección adecuada de $q$ , ${ matemáticas {L}} (Q)$ se vuelve manejable para calcular y maximizar. Por lo tanto, tenemos una aproximación analítica $q$ para el posterior $P({mathbf {Z}}mid {mathbf {X}})$ y un límite inferior ${ matemáticas {L}} (Q)$ para la evidencia logarítmica $log P({mathbf {X}})$ (dado que la divergencia KL no es negativa).

El límite inferior ${ matemáticas {L}} (Q)$ 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 $nombre del operador {E}_{{Q}}[log P({mathbf {Z}},{mathbf {X}})]$ más la entropía de $q$ . El término ${ matemáticas {L}} (Q)$ también se conoce como Evidencia límite inferior, abreviado como ELBO, para enfatizar que es un límite inferior en el registro de evidencia de los datos.

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: ${displaystyle D_{mathrm {KL} }(Qparallel P)geq D_{mathrm {KL} }(Qparallel Q^{*})+D_{mathrm {KL} }(Q^{ *}parallel P),forall Q^{*}in {mathcal {C}}}$

donde ${ matemáticas {C}}$ es un conjunto convexo y la igualdad se cumple si: ${displaystyle Q=Q^{*}triangleq arg min_{Qin {mathcal {C}}}D_{mathrm {KL} }(Qparallel 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: ${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 })}}Correcto),}$

donde la constante de normalización 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}) paralelo 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})}}derecha).}$

El término ${ estilo de visualización zeta ( mathbf {X})}$ a menudo se denomina límite inferior de evidencia (ELBO) en la práctica, ya que ${displaystyle P(mathbf {X})geq zeta (mathbf {X})=exp({mathcal {L}}(Q^{*}))}$ , como se muestra arriba.

Al intercambiar los roles de ${ estilo de visualización mathbf {Z} _ {1}}$ y ${ estilo de visualización mathbf {Z} _ {2},}$ podemos calcular iterativamente los márgenes aproximados ${displaystyle q^{*}(mathbf {Z} _{1})}$ y ${displaystyle q^{*}(mathbf {Z} _{2})}$ del modelo verdadero ${displaystyle P(mathbf {Z} _{1}mid mathbf {X})}$ y ${displaystyle P(mathbf {Z} _{2}mid mathbf {X}),}$ respectivamente. Aunque se garantiza que este esquema iterativo converge monótonamente, el convergido ${displaystyle Q^{*}}$ es solo un minimizador local de ${displaystyle D_{mathrm {KL} }(Qparallel P)}$ .

Si el espacio restringido ${ matemáticas {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.

Aproximación del campo medio

$Q({mathbf {Z}})$ 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 $mathbf{Z}$ en ${mathbf {Z}}_{1}puntos {mathbf {Z}}_{M}$ , $Q({mathbf {Z}})=prod_{{i=1}}^{M}q_{i}({mathbf {Z}}_{i}mid {mathbf {X}})$

Se puede demostrar usando el cálculo de variaciones (de ahí el nombre "Bayes variacional") que la "mejor" distribución $q_{j}^{{*}}$ para cada uno de los factores $q_{j}$ (en términos de la distribución que minimiza la divergencia KL, como se describe arriba) se puede expresar como: $q_{j}^{{*}}({mathbf {Z}}_{j}mid {mathbf {X}})={frac {e^{{operatorname {E}_{{i neq j}}[ln p({mathbf {Z}},{mathbf {X}})]}}}{int e^{{operatorname {E}_{{ineq j} }[ln p({mathbf {Z}},{mathbf {X}})]}},d{mathbf {Z}}_{j}}}$

donde $operatorname {E}_{{ineq j}}[ln p({mathbf {Z}},{mathbf {X}})]$ es la expectativa del logaritmo de la probabilidad conjunta de los datos y las variables latentes, tomada sobre todas las variables que no están en la partición: consulte para una derivación de la distribución $q_{j}^{{*}}({mathbf {Z}}_{j}mid {mathbf {X}})$ .

En la práctica, solemos trabajar en términos de logaritmos, es decir: $ln q_{j}^{{*}}({mathbf {Z}}_{j}mid {mathbf {X}})=operatorname {E}_{{ineq j}}[ ln p({mathbf {Z}},{mathbf {X}})]+{text{constante}}$

La constante en la expresión anterior está relacionada con la constante de normalización (el denominador en la expresión anterior para $q_{j}^{{*}}$ ) 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.).

Usando las propiedades de las expectativas, la expresión $operatorname {E}_{{ineq j}}[ln p({mathbf {Z}},{mathbf {X}})]$ generalmente se puede simplificar 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 más altos como la varianza) de las variables latentes que no están en la partición actual (es decir, variables latentes no incluidas en ${mathbf {Z}}_{j}$ ). 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 sugiere naturalmente un algoritmo iterativo, muy parecido a EM (algoritmo de maximización de expectativas), en el que las expectativas (y posiblemente los momentos más altos) de las variables latentes se inicializan de alguna manera (quizás al azar), y luego los parámetros de cada distribución son calculada a su vez usando 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á.

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 de las variables en cuestión, normalmente se puede determinar la familia de la distribución (que a su vez determina la valor de la constante). La fórmula de 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. Por lo general, estas expectativas se pueden simplificar en funciones de expectativas de las propias variables (es decir, los medios); a veces, expectativas de variables al cuadrado (que pueden estar relacionadas con la varianza de las variables), o expectativas de potencias superiores (es decir, momentos superiores) también aparecen. 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 se pueden expresar como una serie de ecuaciones con dependencias mutuas 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, que en la mayoría de los casos está garantizado para converger. Un ejemplo hará más claro este proceso. 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 se pueden expresar como una serie de ecuaciones con dependencias mutuas 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, que en la mayoría de los casos está garantizado para converger. Un ejemplo hará más claro este proceso. 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 se pueden expresar como una serie de ecuaciones con dependencias mutuas 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, que en la mayoría de los casos está garantizado para converger. Un ejemplo hará más claro este proceso. El resultado es que las fórmulas para los parámetros de las distribuciones de cada variable se pueden expresar como una serie de ecuaciones con dependencias mutuas 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, que en la mayoría de los casos está garantizado para converger. Un ejemplo hará más claro este proceso. El resultado es que las fórmulas para los parámetros de las distribuciones de cada variable se pueden expresar como una serie de ecuaciones con dependencias mutuas 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, que en la mayoría de los casos está garantizado para converger. Un ejemplo hará más claro este proceso.

Una fórmula de dualidad para la inferencia variacional

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

Teorema Considere dos espacios de probabilidad ${displaystyle (Theta,{mathcal {F}},P)}$ y ${displaystyle (Theta,{mathcal {F}},Q)}$ con ${ Displaystyle Q ll P}$ . Suponga que existe una medida de probabilidad dominante común $lambda$ tal que ${displaystyle Plllambda}$ y ${displaystyle Qlllambda}$ . Sea $h$ cualquier variable aleatoria de valor real ${displaystyle (Theta,{mathcal {F}},P)}$ que satisfaga ${displaystyle hen L_{1}(P)}$ . Entonces se cumple la siguiente igualdad ${displaystyle log E_{P}[exp h]={text{sup}}_{Qll P}{E_{Q}[h]-D_{text{KL}}(Q paralela P)}.}$

Además, el supremo en el lado derecho se alcanza si y solo 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 $q$ , donde ${displaystyle p(theta)=dP/dlambda }$ y ${displaystyle q(theta)=dQ/dlambda }$ denotan las derivadas de Radon-Nikodym de las medidas de probabilidad $PAGS$ y $q$ con respecto a $lambda$ , respectivamente.

Un ejemplo básico

Considere un modelo bayesiano simple no jerárquico que consta de un conjunto de observaciones iid de una distribución gaussiana, con media y varianza desconocidas. A continuación, analizamos este modelo con gran detalle para ilustrar el funcionamiento del método variacional de Bayes.

Por conveniencia matemática, en el siguiente ejemplo trabajamos en términos de precisión, es decir, el recíproco de la varianza (o en una Gaussiana multivariada, la inversa de la matriz de covarianza), en lugar de la varianza en sí. (Desde un punto de vista teórico, la precisión y la varianza son equivalentes ya que existe una correspondencia uno a uno entre los dos).

El modelo matemático

Colocamos distribuciones previas conjugadas en la media $mu$ y la precisión desconocidas $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{alineado}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{número de puntos de datos}}end{alineados}}}$

Los hiperparámetros $mu_0, lambda_0, a_0$ y $b_{0}$ en las distribuciones anteriores 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 $mu$ y $tau$ .

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

La probabilidad conjunta

La probabilidad conjunta de todas las variables se puede reescribir como $p({mathbf {X}},mu,tau)=p({mathbf {X}}mid mu,tau)p(mu mid tau)p(tau)$

donde están los factores individuales $begin{align} p(mathbf{X}mid mu,tau) & = prod_{n=1}^N mathcal{N}(x_nmid mu,tau^{-1}) \ p(mumid tau) & = mathcal{N} left (mumid mu_0, (lambda_0 tau)^{-1} right) \ p(tau) & = nombre del operador{Gamma}(taumid a_0, b_0) end{align}$

dónde ${begin{alineado}{mathcal {N}}(xmid mu,sigma ^{2})&={frac {1}{{sqrt {2pi sigma ^{2}} }}}e^{{{frac {-(x-mu)^{2}}{2sigma ^{2}}}}}\nombre del operador {Gamma}(tau mid a,b)&={frac {1}{Gamma (a)}}b^{a}tau ^{{a-1}}e^{{-btau }}end{alineado}}$

Aproximación factorizada

Suponga que $q(mu,tau)=q(mu)q(tau)$ , es decir, que la distribución posterior se factoriza en factores independientes para $mu$ y $tau$ . Este tipo de suposición subyace en el método bayesiano variacional. De hecho, la verdadera distribución posterior no 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.

Derivación de q (μ)

Después ${displaystyle {begin{alineado}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}\ &=nombre del operador {E} _{tau}left[ln prod _{n=1}^{N}{sqrt {frac {tau }{2pi }}}e^{-{ frac {(x_{n}-mu)^{2}tau }{2}}}right]+operatorname {E} _{tau }left[ln {sqrt {frac { lambda _{0}tau}{2pi}}}e^{-{frac {(mu -mu _{0})^{2}lambda _{0}tau}{2 }}}right]+C_{2}\&=nombre del operador {E}_{tau}left[sum_{n=1}^{N}left({frac {1}{ 2}}(ln tau -ln 2pi)-{frac {(x_{n}-mu)^{2}tau }{2}}right)right]+operatorname { E} _{tau }left[{frac {1}{2}}(ln lambda _{0}+ln tau -ln 2pi)-{frac {(mu - mu _{0})^{2}lambda _{0}tau}{2}}right]+C_{2}\&=nombre del operador {E}_{tau}left[sum_{n=1}^{N}-{frac {(x_{n}-mu)^{2}tau }{2} }right]+nombre del operador {E} _{tau }left[-{frac {(mu -mu _{0})^{2}lambda _{0}tau }{2} }right]+nombre del operador {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 2pi) 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} derecha}+C_{3}end{alineado}}}=nombre del operador {E}_{tau}left[sum_{n=1}^{N}-{frac {(x_{n}-mu)^{2}tau }{2} }right]+nombre del operador {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{ alineado}}}=nombre del operador {E}_{tau}left[sum_{n=1}^{N}-{frac {(x_{n}-mu)^{2}tau }{2} }right]+nombre del operador {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{ alineado}}}$

En la derivación anterior, $C$ , $C_{2}$ y $C_{3}$ se refieren a valores que son constantes con respecto a $mu$ . Tenga en cuenta que el término $nombre del operador {E}_{{tau }}[ln p(tau)]$ no es una función de $mu$ y tendrá el mismo valor independientemente del valor de $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 $mu$ . Dado que este es el logaritmo de $q_{mu}^{*}(mu)$ , podemos ver que $q_{mu}^{*}(mu)$ en sí mismo es una distribución gaussiana.

Con una cierta cantidad de matemáticas tediosas (expandiendo los cuadrados dentro de las llaves, separando y agrupando los términos que involucran $mu$ y $mu^{2}$ y completando el cuadrado sobre $mu$ ), podemos derivar los parámetros de la distribución gaussiana: $begin{align} 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_nmu + mu^2) + lambda_0(mu^2-2mu_0mu + mu_0^ 2) right } + C_3 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} left{ left(sum_{n=1}^N x_n^2right)-2left(sum_{n=1}^N x_nright)mu + left (sum_{n=1}^N mu^2 right) + lambda_0 mu^2-2lambda_0mu_0mu + lambda_0mu_0^2 right} + C_3 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} left{ (lambda_0+N)mu^2 -2left(lambda_0mu_0 + sum_{n=1}^N x_nright)mu + left(sum_{n=1} ^N x_n^2right) + lambda_0mu_0^2 right} + C_3 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} left{ (lambda_0+N)mu^2 -2left(lambda_0mu_0 + sum_{ n=1}^N x_nright)mu right} + C_4 \ &= -frac{operatorname{E}_{tau}[tau]}{2} left{ ( lambda_0+N)mu^2 -2left(frac{lambda_0mu_0 + sum_{n=1}^N x_n}{lambda_0+N} right)(lambda_0+N) mu derecha} + C_4 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} izquierda{ (lambda_0+N)izquierda(mu^2 -2 left(frac{lambda_0mu_0 + sum_{n=1}^N x_n}{lambda_0+N}right) muright) right} + C_4 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} left{ (lambda_0+N)left(mu^2 -2left(frac{lambda_0 mu_0 + sum_{n=1}^N x_n}{lambda_0+N}right) mu + left(frac{lambda_0mu_0 + sum_{n=1}^N x_n}{ lambda_0+N}right)^2 - left(frac{lambda_0mu_0 + sum_{n=1}^N x_n}{lambda_0+N}right)^2right) right} + C_4 \ &= -frac{nombre del operador{E}_{tau}[tau]}{2} left{ (lambda_0+N)left(mu^2 -2left(frac{lambda_0mu_0 + sum_{n=1}^N x_n}{lambda_0+N}right) mu + left(frac{lambda_0mu_0 + sum_{n=1}^ N x_n}{lambda_0+N}right)^2 right) right} + C_5 \ &= -frac{operatorname{E}_{tau}[tau]}{2} izquierda{ (lambda_0+N)izquierda(mu-frac{lambda_0mu_0 + sum_{n=1}^N x_n}{lambda_0+N}derecha)^2 derecha} + C_5 \ &= -frac{1}{2} (lambda_0+N)nombre del operador{E}_{tau}[tau] left(mu-frac{lambda_0mu_0 + sum_{n=1 }^N x_n}{lambda_0+N}right)^2 + C_5 end{align}$

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

En otras palabras: ${displaystyle {begin{alineado}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)nombre del operador {E} _{tau}[tau]\{bar {x}}&={ fracción {1}{N}}sum _{n=1}^{N}x_{n}end{alineado}}}$

Derivación de q(τ)

La derivación de $q_{tau}^{*}(tau)$ es similar a la anterior, aunque omitimos algunos de los detalles en aras de la brevedad. $begin{align} ln q_tau^*(tau) &= operatorname{E}_{mu}[ln p(mathbf{X}mid mu,tau) + ln p (mumid tau)] + ln p(tau) + text{constante} \ &= (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{constante} end{align}$

Exponenciando ambos lados, podemos ver que $q_{tau}^{*}(tau)$ es una distribución gamma. Específicamente: ${begin{alineado}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}}nombre del operador {E}_{mu } left[sum _{{n=1}}^{N}(x_{n}-mu)^{2}+lambda _{0}(mu -mu _{0})^{ 2}right]end{alineado}}$

Algoritmo para calcular los parámetros

Recapitulemos las conclusiones de los apartados anteriores: ${displaystyle {begin{alineado}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)nombre del operador {E} _{tau}[tau]\{bar {x}}&={ fracción {1}{N}}sum _{n=1}^{N}x_{n}end{alineado}}}$

y ${begin{alineado}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}}nombre del operador {E}_{mu } left[sum _{{n=1}}^{N}(x_{n}-mu)^{2}+lambda _{0}(mu -mu _{0})^{ 2}right]end{alineado}}$

En cada caso, los parámetros de distribución sobre una de las variables dependen de las expectativas que se tomen con respecto a la otra variable. Podemos expandir las expectativas, usando las fórmulas estándar para las expectativas de momentos de las distribuciones gaussiana y gamma: $begin{align} operatorname{E}[taumid a_N, b_N] &= frac{a_N}{b_N} \ operatorname{E} left [mumidmu_N,lambda_N^{ -1} right ] &= mu_N \ operatorname{E}left[X^2 right] &= operatorname{Var}(X) + (operatorname{E}[X])^2 operatorname{E} left [mu^2midmu_N,lambda_N^{-1} right] &= lambda_N^{-1} + mu_N^2 end{align}$

Aplicar estas fórmulas a las ecuaciones anteriores es trivial en la mayoría de los casos, pero la ecuación $b_{N}$ requiere más trabajo: $begin{align} b_N &= b_0 + frac{1}{2} operatorname{E}_mu left[sum_{n=1}^N (x_n-mu)^2 + lambda_0(mu - mu_0)^2right] \ &= b_0 + frac{1}{2} operatorname{E}_mu left[ (lambda_0+N)mu^2 -2 left (lambda_0mu_0 + sum_{n=1}^N x_n right)mu + left(sum_{n=1}^N x_n^2 right) + lambda_0mu_0^2 right] \ &= b_0 + frac{1}{2} left[ (lambda_0+N)operatorname{E}_mu[mu^2] -2 left (lambda_0mu_0 + sum_{ n=1}^N x_n right)nombre del operador{E}_mu [mu] + left (sum_{n=1}^N x_n^2 right) + lambda_0mu_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_0mu_0^2 right] \ end {alinear}$

Entonces podemos escribir las ecuaciones de los parámetros de la siguiente manera, sin ninguna expectativa: $begin{align} 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_0mu_0 + sum_{n =1}^N x_n right)mu_N + left (sum_{n=1}^N x_n^2 right) + lambda_0mu_0^2 right] end{align}$

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

Calcule $sum_{{n=1}}^{N}x_{n}$ y $sum _{{n=1}}^{N}x_{n}^{2}.$ Use estos valores para calcular $mu _{N}$ y $un}.$
Inicializar $lambda _{N}$ a algún valor arbitrario.
Use el valor actual de $lambda _{N},$ junto con los valores conocidos de los otros parámetros para calcular $b_{N}$ .
Use el valor actual de $b_{N},$ junto con los valores conocidos de los otros parámetros para calcular $lambda _{N}$ .
Repita los dos últimos pasos hasta la convergencia (es decir, hasta que ningún valor 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 posterior, por ejemplo, su media y varianza, una región de densidad más alta del 95% (el intervalo más pequeño que incluye 95 % de la probabilidad total), etc.

Se puede demostrar que este algoritmo está garantizado para converger a un máximo local.

Tenga en cuenta también que las distribuciones posteriores tienen la misma forma que las distribuciones anteriores correspondientes. No asumimos esto; la única suposición que hicimos fue que las distribuciones se factorizan, y la forma de las distribuciones siguió naturalmente. Resulta (ver más abajo) que el hecho de que las distribuciones posteriores tengan la misma forma que las distribuciones anteriores no es una coincidencia, sino un resultado general siempre que las distribuciones anteriores sean miembros de la familia exponencial, que es el caso de la mayoría de las distribuciones anteriores. distribuciones estándar.

Más discusión

Receta paso a paso

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

Describir la red con un modelo gráfico, identificando las variables observadas (datos) $mathbf{X}$ y no observadas (parámetros ${ símbolo de negrita Theta }$ y variables latentes $mathbf{Z}$ ) y sus distribuciones de probabilidad condicional. El bayesiano variacional construirá entonces una aproximación a la probabilidad posterior $p({mathbf {Z}},{boldsymbol Theta }mid {mathbf {X}})$ . 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.
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; la creación de demasiados subconjuntos produce una mala aproximación, mientras que la creación de muy pocos hace intratable todo el procedimiento variacional de Bayes. 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 tratable. Suponga que las particiones se llaman ${mathbf {Z}}_{1},ldots,{mathbf {Z}}_{M}$ .
Para una partición dada ${mathbf {Z}}_{j}$ , escriba la fórmula para la mejor aproximación de distribución $q_{j}^{{*}}({mathbf {Z}}_{j}mid {mathbf {X}})$ utilizando la ecuación básica $ln q_{j}^{{*}}({mathbf {Z}}_{j}mid {mathbf {X}})=operatorname {E}_{{ineq j}}[ ln p({mathbf {Z}},{mathbf {X}})]+{text{constante}}$ .
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 ${mathbf {Z}}_{j}$ puede ignorarse; se doblarán en el término constante.
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 ${mathbf {Z}}_{j}$ (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 esperados son constantes con respecto a las variables en la partición actual.
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, exponenciar 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. Es posible que se requiera una manipulación matemática significativa para convertir la fórmula en una forma que coincida con el PDF de una distribución conocida. Cuando se puede hacer esto, la constante de normalización se puede restablecer por definición y las ecuaciones para los parámetros de la distribución conocida se pueden derivar extrayendo las partes apropiadas de la fórmula.
Cuando todas las expectativas se pueden reemplazar analíticamente con funciones de variables que no están en la partición actual, y la PDF se pone 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 la parámetros de variables en otras particiones.
Cuando este procedimiento se puede aplicar a todas las particiones, el resultado es un conjunto de ecuaciones vinculadas entre sí que especifican los valores óptimos de todos los parámetros.
Luego se aplica un procedimiento de tipo maximización de expectativas (EM), seleccionando un valor inicial para cada parámetro y la iteración a través de una serie de pasos, donde en cada paso recorremos las ecuaciones, actualizando cada parámetro a su vez. Esto está garantizado para converger.

Puntos más importantes

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

La idea de 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. Como en otros métodos bayesianos, pero a diferencia, por ejemplo, de la maximización de expectativas (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. Luego, las estimaciones de las variables se pueden derivar en las formas bayesianas estándar, por ejemplo, calculando la media de la distribución para obtener una estimación de un solo punto o derivando un intervalo creíble, la región de mayor densidad, etc.
"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, por lo general coincidirá bastante en los momentos más bajos de las variables no observadas, por ejemplo, la media y la varianza.
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 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 a EM.

Comparado con la maximización de expectativas (EM)

El bayesiano variacional (VB) a menudo se compara con la maximización de expectativas (EM). El procedimiento numérico real es bastante similar, ya que ambos son procedimientos iterativos alternos que convergen sucesivamente en valores de parámetros óptimos. Los pasos iniciales para derivar los procedimientos respectivos también son vagamente similares, ambos comienzan con fórmulas para densidades de probabilidad y ambos involucran cantidades significativas de manipulaciones matemáticas.

Sin embargo, hay una serie de diferencias. Lo más importante es lo que se está calculando.

EM calcula estimaciones puntuales de la distribución posterior de aquellas variables aleatorias que se pueden categorizar como "parámetros", pero solo estimaciones de las distribuciones posteriores reales de las variables latentes (al menos en "EM suave", y a menudo solo cuando las variables latentes son discretas). Las estimaciones puntuales calculadas son las modas de estos parámetros; No hay otra informacion disponible.
VB, por otro lado, calcula estimaciones de la distribución posterior real de todas las variables, tanto parámetros como variables latentes. Cuando es necesario derivar estimaciones puntuales, generalmente se usa la media en lugar de la moda, como es normal en la inferencia bayesiana. Concomitante con esto, los parámetros calculados en VB notienen el mismo significado que los de EM. EM calcula los valores óptimos de los parámetros de la propia red Bayes. 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 de Bayes. Por ejemplo, un modelo típico de mezcla gaussiana 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. Sin embargo, VB primero ajustaría una distribución a estos parámetros, generalmente en forma de una distribución previa, por ejemplo, una distribución gamma inversa de escala normal, y luego calcularía los valores para los parámetros de esta distribución previa, es decir, esencialmente hiperparámetros. En este caso,

Un ejemplo más complejo

Imagine un modelo de mezcla bayesiano gaussiano descrito de la siguiente manera: ${begin{alineado}{mathbf {pi }}&sim operatorname {SymDir}(K,alpha _{0})\{mathbf {Lambda }}_{{i=1dots K}}&sim {mathcal {W}}({mathbf {W}}_{0},nu _{0})\{mathbf {mu }}_{{i=1 puntos K}}&sim {mathcal {N}}({mathbf {mu }}_{0},(beta _{0}{mathbf {Lambda }}_{i})^{ {-1}})\{mathbf {z}}[i=1dots N]&sim operatorname {Mult}(1,{mathbf {pi }})\{mathbf {x }}_{{i=1dots N}}&sim {mathcal {N}}({mathbf {mu }}_{{z_{i}}},{{mathbf {Lambda } }_{{z_{i}}}}^{{-1}})\K&={text{número de componentes de mezcla}}\N&={text{número de puntos de datos}}end{ alineado}}$

Nota:

SymDir() es la distribución simétrica de Dirichlet de dimension $k$ , con el hiperparámetro para cada componente establecido en $alpha _{0}$ . La distribución de Dirichlet es el conjugado anterior de la distribución categórica o distribución multinomial.
${mathcal {W}}()$ es la distribución de Wishart, que es el conjugado anterior de la matriz de precisión (matriz de covarianza inversa) para una distribución gaussiana multivariante.
Mult() es una distribución multinomial sobre una sola observación (equivalente a una distribución categórica). El espacio de estado es una representación "uno de K", es decir, un $k$ 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.
$mathcal{N}()$ es la distribución gaussiana, en este caso concretamente la distribución gaussiana multivariante.

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

${mathbf {X}}={{mathbf {x}}_{1},puntos,{mathbf {x}}_{N}}$ es el conjunto de $norte$ puntos de datos, cada uno de los cuales es un $D$ vector dimensional distribuido de acuerdo con una distribución gaussiana multivariada.
${mathbf {Z}}={{mathbf {z}}_{1},puntos,{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 $z_{{nk}}$ para $k=1puntos K$ , como se describe anteriormente.
$mathbf {pi}$ es la proporción de $k$ mezcla de los componentes de la mezcla.
${mathbf {mu}}_{{i=1dots K}}$ y ${mathbf {Lambda }}_{{i=1dots K}}$ especificar los parámetros (media y precisión) asociados a cada componente de la mezcla.

La probabilidad conjunta de todas las variables se puede reescribir como $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 ${begin{alineado}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 (Kalpha _{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{alineado}}$

dónde ${begin{alineado}{mathcal {N}}({mathbf {x}}mid {mathbf {mu }},{mathbf {Sigma }})&={frac {1}{ (2pi)^{{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 izquierda (-{frac {1}{2}}operatorname {Tr}({mathbf {W}}^{{-1}}{mathbf {Lambda }})derecha)\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{dimensionalidad de cada punto de datos}}end{alineado}}$

Suponga $q({mathbf {Z}},{mathbf {pi }},{mathbf {mu }},{mathbf {Lambda }})=q({mathbf {Z}})q({mathbf {pi }},{mathbf {mu }},{mathbf {Lambda }})$ que

Después ${begin{alineado}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{constante}}\&=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{constante}}\&=sum_{n =1}}^{N}sum_{{k=1}}^{K}z_{{nk}}ln rho_{{nk}}+{text{constante}}end{alineado }}$

donde hemos definido $ln rho _{{nk}}=nombre de operador {E}[ln pi _{k}]+{frac {1}{2}}nombre de operador {E}[ln |{mathbf { Lambda }}_{k}|]-{frac {D}{2}}ln(2pi)-{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 $ln q^{*}({mathbf{Z}})$ rendimientos $q^{*}({mathbf {Z}})propto prod_{{n=1}}^{N}prod_{{k=1}}^{K}rho_{nk }}^{{z_{{nk}}}}$

Requerir que esto sea normalizado termina requiriendo que la $rho _{{nk}}$ suma sea 1 sobre todos los valores de $k$ , produciendo $q^{*}({mathbf {Z}})=prod_{{n=1}}^{N}prod_{{k=1}}^{K}r_{{nk}}^ {{z_{{nk}}}}$

dónde $r_{{nk}}={frac {rho_{{nk}}}{sum_{{j=1}}^{K}rho_{{nj}}}}$

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

Además, notamos que $nombre del operador {E}[z_{{nk}}]=r_{{nk}},$

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

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

Después, ${begin{alineado}ln q^{*}({mathbf {pi }})&=ln p({mathbf {pi }})+operatorname {E}_{{{mathbf {Z}}}}[ln p({mathbf {Z}}mid {mathbf {pi }})]+{text{constante}}\&=(alpha _{0}- 1)sum_{{k=1}}^{K}ln pi_{k}+sum_{{n=1}}^{N}sum_{{k=1}}^ {K}r_{{nk}}ln pi _{k}+{text{constante}}end{alineado}}$

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

dónde $alpha _{k}=alpha _{0}+N_{k},$

dónde $N_{k}=sum_{{n=1}}^{N}r_{{nk}},$

Finalmente $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{constante}}$

Agrupando y leyendo términos que involucran ${mathbf {mu}}_{k}$ y ${mathbf {Lambda}}_{k}$ , el resultado es una distribución Gaussiana-Wishart dada por $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{alineado}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}&={ fracción {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{alineado}}}$

Finalmente, observe que estas funciones requieren los valores de $r_{{nk}}$ , que hacen uso de $rho _{{nk}}$ , que se define a su vez en base a $nombre del operador {E}[ln pi _{k}]$ , $nombre del operador {E}[ln |{mathbf {Lambda }}_{k}|]$ y $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})]$ . Ahora que hemos determinado las distribuciones sobre las que se toman estas expectativas, podemos derivar fórmulas para ellas: ${displaystyle {begin{alineado}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})]&=Dbeta_{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)+Dln 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{alineado}}}$

Estos resultados conducen a ${displaystyle r_{nk}propto {widetilde {pi }}_{k}{widetilde {Lambda }}_{k}^{1/2}exp left{-{frac { D}{2beta_{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 valores absolutos mediante la normalización $k$ para que los valores correspondientes sumen 1.

Tenga en cuenta que:

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

Esto sugiere un procedimiento iterativo que alterna entre dos pasos:

Un paso E que calcula el valor de $r_{{nk}}$ usar los valores actuales de todos los demás parámetros.
Un paso M que usa el nuevo valor de $r_{{nk}}$ para calcular nuevos valores de todos los demás parámetros.

Tenga en cuenta 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 $r_{{nk}}$ en el paso E corresponden estrechamente a las probabilidades posteriores de las variables latentes dados los datos, es decir $p({mathbf {Z}}mid {mathbf {X}})$ ; el cálculo de las estadísticas $N_{k}$ , ${{bar {{mathbf {x}}}}}_{k}$ , y ${mathbf{S}}_{k}$ corresponde estrechamente al cálculo de las correspondientes estadísticas de "recuento suave" 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.

Distribuciones de familia exponencial

Tenga en cuenta 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 correspondiente. distribución previa sobre la variable. Este es un resultado general que se cumple para todas las distribuciones anteriores derivadas de la familia exponencial.

Te puede interesar
Demografía de Gambia
(leer más)
Te puede interesar
Condado de Nassau, Florida
(leer más)
Te puede interesar
Metanálisis
Un metanálisis es un análisis estadístico que combina los resultados de múltiples estudios científicos. Los metanálisis se pueden realizar cuando hay... (leer más)
Más resultados...