Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Variational Bayesian methods
open-in-new
<h2 id="mathematical-derivation">Mathematical derivation</h2> <h3>Problem</h3> <p>In <a href="/facts/Calculus_of_variations/Ya5Lsma3">variational</a> inference, the posterior distribution over a set of unobserved variables Z = { Z 1 … Z n } {\displaystyle \mathbf {Z} =\{Z_{1}\dots Z_{n}\}} given some data X {\displaystyle \mathbf {X} } is approximated by a so-called variational distribution, Q ( Z ) : {\displaystyle Q(\mathbf {Z} ):} </p> P ( Z ∣ X ) ≈ Q ( Z ) . {\displaystyle P(\mathbf {Z} \mid \mathbf {X} )\approx Q(\mathbf {Z} ).} <p>The distribution Q ( Z ) {\displaystyle Q(\mathbf {Z} )} is restricted to belong to a family of distributions of simpler form than P ( Z ∣ X ) {\displaystyle P(\mathbf {Z} \mid \mathbf {X} )} (e.g. a family of Gaussian distributions), selected with the intention of making Q ( Z ) {\displaystyle Q(\mathbf {Z} )} similar to the true posterior, P ( Z ∣ X ) {\displaystyle P(\mathbf {Z} \mid \mathbf {X} )} . </p><p>The similarity (or dissimilarity) is measured in terms of a dissimilarity function d ( Q ; P ) {\displaystyle d(Q;P)} and hence inference is performed by selecting the distribution Q ( Z ) {\displaystyle Q(\mathbf {Z} )} that minimizes d ( Q ; P ) {\displaystyle d(Q;P)} . </p> <h3>KL divergence</h3> <p>The most common type of variational Bayes uses the <a href="/facts/Kullback%25E2%2580%2593Leibler_divergence/nh7SjlPE">Kullback–Leibler divergence</a> (KL-divergence) of <i>Q</i> from <i>P</i> as the choice of dissimilarity function. This choice makes this minimization tractable. The KL-divergence is defined as </p> D K L ( Q ∥ P ) ≜ ∑ Z Q ( Z ) log Q ( Z ) P ( Z ∣ X ) . {\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} )}}.} <p>Note that <i>Q</i> and <i>P</i> are reversed from what one might expect. This use of reversed KL-divergence is conceptually similar to the <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">expectation–maximization algorithm</a>. (Using the KL-divergence in the other way produces the <a href="/facts/Expectation_propagation/sSeGGBfY">expectation propagation</a> algorithm.) </p> <h3>Intractability</h3> <p>Variational techniques are typically used to form an approximation for: </p> P ( Z ∣ X ) = P ( X ∣ Z ) P ( Z ) P ( X ) = P ( X ∣ Z ) P ( Z ) ∫ Z P ( X , Z ′ ) d Z ′ {\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} '}}} <p>The marginalization over Z {\displaystyle \mathbf {Z} } to calculate P ( X ) {\displaystyle P(\mathbf {X} )} in the denominator is typically intractable, because, for example, the search space of Z {\displaystyle \mathbf {Z} } is combinatorially large. Therefore, we seek an approximation, using Q ( Z ) ≈ P ( Z ∣ X ) {\displaystyle Q(\mathbf {Z} )\approx P(\mathbf {Z} \mid \mathbf {X} )} . </p> <h3>Evidence lower bound</h3> <p class="note">Main article: <a href="/facts/Evidence_lower_bound/cANNQv7y">Evidence lower bound</a></p> <p>Given that P ( Z ∣ X ) = P ( X , Z ) P ( X ) {\displaystyle P(\mathbf {Z} \mid \mathbf {X} )={\frac {P(\mathbf {X} ,\mathbf {Z} )}{P(\mathbf {X} )}}} , the KL-divergence above can also be written as </p> D K L ( Q ∥ P ) = ∑ Z Q ( Z ) [ log Q ( Z ) P ( Z , X ) + log P ( X ) ] = ∑ Z Q ( Z ) [ log Q ( Z ) − log P ( Z , X ) ] + ∑ Z Q ( Z ) [ log P ( X ) ] {\displaystyle {\begin{array}{rl}D_{\mathrm {KL} }(Q\parallel P)&=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log {\frac {Q(\mathbf {Z} )}{P(\mathbf {Z} ,\mathbf {X} )}}+\log P(\mathbf {X} )\right]\\&=\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]+\sum _{\mathbf {Z} }Q(\mathbf {Z} )\left[\log P(\mathbf {X} )\right]\end{array}}} <p>Because P ( X ) {\displaystyle P(\mathbf {X} )} is a constant with respect to Z {\displaystyle \mathbf {Z} } and ∑ Z Q ( Z ) = 1 {\displaystyle \sum _{\mathbf {Z} }Q(\mathbf {Z} )=1} because Q ( Z ) {\displaystyle Q(\mathbf {Z} )} is a distribution, we have </p> D K L ( Q ∥ P ) = ∑ Z Q ( Z ) [ log Q ( Z ) − log P ( Z , X ) ] + log P ( X ) {\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} )} <p>which, according to the definition of <a href="/facts/Expected_value/1XV0JKL8">expected value</a> (for a discrete <a href="/facts/Random_variable/TwTBXnLT">random variable</a>), can be written as follows </p> D K L ( Q ∥ P ) = E Q [ log Q ( Z ) − log P ( Z , X ) ] + log P ( X ) {\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} )} <p>which can be rearranged to become </p> log P ( X ) = D K L ( Q ∥ P ) − E Q [ log Q ( Z ) − log P ( Z , X ) ] = D K L ( Q ∥ P ) + L ( Q ) {\displaystyle {\begin{array}{rl}\log P(\mathbf {X} )&=D_{\mathrm {KL} }(Q\parallel P)-\mathbb {E} _{\mathbf {Q} }\left[\log Q(\mathbf {Z} )-\log P(\mathbf {Z} ,\mathbf {X} )\right]\\&=D_{\mathrm {KL} }(Q\parallel P)+{\mathcal {L}}(Q)\end{array}}} <p>As the <i>log-<a href="/facts/Model_evidence/4CYPcBYb">evidence</a></i> log P ( X ) {\displaystyle \log P(\mathbf {X} )} is fixed with respect to Q {\displaystyle Q} , maximizing the final term L ( Q ) {\displaystyle {\mathcal {L}}(Q)} minimizes the KL divergence of Q {\displaystyle Q} from P {\displaystyle P} . By appropriate choice of Q {\displaystyle Q} , L ( Q ) {\displaystyle {\mathcal {L}}(Q)} becomes tractable to compute and to maximize. Hence we have both an analytical approximation Q {\displaystyle Q} for the posterior P ( Z ∣ X ) {\displaystyle P(\mathbf {Z} \mid \mathbf {X} )} , and a lower bound L ( Q ) {\displaystyle {\mathcal {L}}(Q)} for the log-evidence log P ( X ) {\displaystyle \log P(\mathbf {X} )} (since the KL-divergence is non-negative). </p><p>The lower bound L ( Q ) {\displaystyle {\mathcal {L}}(Q)} is known as the (negative) variational free energy in analogy with <a href="/facts/Thermodynamic_free_energy/GtjtQqnq">thermodynamic free energy</a> because it can also be expressed as a negative energy E Q [ log P ( Z , X ) ] {\displaystyle \operatorname {E} _{Q}[\log P(\mathbf {Z} ,\mathbf {X} )]} plus the <a href="/facts/Entropy_(information_theory)/NLg4NLvt">entropy</a> of Q {\displaystyle Q} . The term L ( Q ) {\displaystyle {\mathcal {L}}(Q)} is also known as Evidence Lower Bound, abbreviated as <a href="/facts/Evidence_lower_bound/cANNQv7y">ELBO</a>, to emphasize that it is a lower (worst-case) bound on the log-evidence of the data. </p> <h3>Proofs</h3> <p>By the generalized <a href="/facts/Pythagorean_theorem/Lh6x2Kr5">Pythagorean theorem</a> of <a href="/facts/Bregman_divergence/wREwF8qo">Bregman divergence</a>, of which KL-divergence is a special case, it can be shown that:<a class="footnote-ref" id="fnref:1" href="#fn:1"><sup>1</sup></a><a class="footnote-ref" id="fnref:2" href="#fn:2"><sup>2</sup></a> </p> D K L ( Q ∥ P ) ≥ D K L ( Q ∥ Q ∗ ) + D K L ( Q ∗ ∥ P ) , ∀ Q ∗ ∈ C {\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}}} <p>where C {\displaystyle {\mathcal {C}}} is a convex set and the equality holds if: </p> Q = Q ∗ ≜ arg min Q ∈ C D K L ( Q ∥ P ) . {\displaystyle Q=Q^{*}\triangleq \arg \min _{Q\in {\mathcal {C}}}D_{\mathrm {KL} }(Q\parallel P).} <p>In this case, the global minimizer Q ∗ ( Z ) = q ∗ ( Z 1 ∣ Z 2 ) q ∗ ( Z 2 ) = q ∗ ( Z 2 ∣ Z 1 ) q ∗ ( Z 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}),} with Z = { Z 1 , Z 2 } , {\displaystyle \mathbf {Z} =\{\mathbf {Z_{1}} ,\mathbf {Z_{2}} \},} can be found as follows:<a class="footnote-ref" id="fnref:3" href="#fn:3"><sup>3</sup></a> </p> q ∗ ( Z 2 ) = P ( X ) ζ ( X ) P ( Z 2 ∣ X ) exp ( D K L ( q ∗ ( Z 1 ∣ Z 2 ) ∥ P ( Z 1 ∣ Z 2 , X ) ) ) = 1 ζ ( X ) exp E q ∗ ( Z 1 ∣ Z 2 ) ( log P ( Z , X ) q ∗ ( Z 1 ∣ Z 2 ) ) , {\displaystyle {\begin{array}{rl}q^{*}(\mathbf {Z} _{2})&={\frac {P(\mathbf {X} )}{\zeta (\mathbf {X} )}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}\\&={\frac {1}{\zeta (\mathbf {X} )}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right),\end{array}}} <p>in which the normalizing constant is: </p> ζ ( X ) = P ( X ) ∫ Z 2 P ( Z 2 ∣ X ) exp ( D K L ( q ∗ ( Z 1 ∣ Z 2 ) ∥ P ( Z 1 ∣ Z 2 , X ) ) ) = ∫ Z 2 exp E q ∗ ( Z 1 ∣ Z 2 ) ( log P ( Z , X ) q ∗ ( Z 1 ∣ Z 2 ) ) . {\displaystyle {\begin{array}{rl}\zeta (\mathbf {X} )&=P(\mathbf {X} )\int _{\mathbf {Z} _{2}}{\frac {P(\mathbf {Z} _{2}\mid \mathbf {X} )}{\exp(D_{\mathrm {KL} }(q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})\parallel P(\mathbf {Z} _{1}\mid \mathbf {Z} _{2},\mathbf {X} )))}}\\&=\int _{\mathbf {Z} _{2}}\exp \mathbb {E} _{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}\left(\log {\frac {P(\mathbf {Z} ,\mathbf {X} )}{q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})}}\right).\end{array}}} <p>The term ζ ( X ) {\displaystyle \zeta (\mathbf {X} )} is often called the <a href="/facts/Model_evidence/4CYPcBYb">evidence</a> lower bound (ELBO) in practice, since P ( X ) ≥ ζ ( X ) = exp ( L ( Q ∗ ) ) {\displaystyle P(\mathbf {X} )\geq \zeta (\mathbf {X} )=\exp({\mathcal {L}}(Q^{*}))} ,<a class="footnote-ref" id="fnref:4" href="#fn:4"><sup>4</sup></a> as shown above. </p><p>By interchanging the roles of Z 1 {\displaystyle \mathbf {Z} _{1}} and Z 2 , {\displaystyle \mathbf {Z} _{2},} we can iteratively compute the approximated q ∗ ( Z 1 ) {\displaystyle q^{*}(\mathbf {Z} _{1})} and q ∗ ( Z 2 ) {\displaystyle q^{*}(\mathbf {Z} _{2})} of the true model's marginals P ( Z 1 ∣ X ) {\displaystyle P(\mathbf {Z} _{1}\mid \mathbf {X} )} and P ( Z 2 ∣ X ) , {\displaystyle P(\mathbf {Z} _{2}\mid \mathbf {X} ),} respectively. Although this iterative scheme is guaranteed to converge monotonically,<a class="footnote-ref" id="fnref:5" href="#fn:5"><sup>5</sup></a> the converged Q ∗ {\displaystyle Q^{*}} is only a local minimizer of D K L ( Q ∥ P ) {\displaystyle D_{\mathrm {KL} }(Q\parallel P)} . </p><p>If the constrained space C {\displaystyle {\mathcal {C}}} is confined within independent space, i.e. q ∗ ( Z 1 ∣ Z 2 ) = q ∗ ( Z 1 ) , {\displaystyle q^{*}(\mathbf {Z} _{1}\mid \mathbf {Z} _{2})=q^{*}(\mathbf {Z_{1}} ),} the above iterative scheme will become the so-called mean field approximation Q ∗ ( Z ) = q ∗ ( Z 1 ) q ∗ ( Z 2 ) , {\displaystyle Q^{*}(\mathbf {Z} )=q^{*}(\mathbf {Z} _{1})q^{*}(\mathbf {Z} _{2}),} as shown below. </p> <h2 id="mean-field-approximation">Mean field approximation</h2> <p>The variational distribution Q ( Z ) {\displaystyle Q(\mathbf {Z} )} is usually assumed to factorize over some <a href="/facts/Partition_of_a_set/VeP9g8CA">partition</a> of the latent variables, i.e. for some partition of the latent variables Z {\displaystyle \mathbf {Z} } into Z 1 … Z M {\displaystyle \mathbf {Z} _{1}\dots \mathbf {Z} _{M}} , </p> Q ( Z ) = ∏ i = 1 M q i ( Z i ∣ X ) {\displaystyle Q(\mathbf {Z} )=\prod _{i=1}^{M}q_{i}(\mathbf {Z} _{i}\mid \mathbf {X} )} <p>It can be shown using the <a href="/facts/Calculus_of_variations/Ya5Lsma3">calculus of variations</a> (hence the name "variational Bayes") that the "best" distribution q j ∗ {\displaystyle q_{j}^{*}} for each of the factors q j {\displaystyle q_{j}} (in terms of the distribution minimizing the KL divergence, as described above) satisfies:<a class="footnote-ref" id="fnref:6" href="#fn:6"><sup>6</sup></a> </p> q j ∗ ( Z j ∣ X ) = e E q − j ∗ [ ln p ( Z , X ) ] ∫ e E q − j ∗ [ ln p ( Z , X ) ] d Z 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}}}} <p>where E q − j ∗ [ ln p ( Z , X ) ] {\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]} is the <a href="/facts/Expected_value/1XV0JKL8">expectation</a> of the logarithm of the <a href="/facts/Joint_probability/klX2ksGY">joint probability</a> of the data and latent variables, taken with respect to q ∗ {\displaystyle q^{*}} over all variables not in the partition: refer to Lemma 4.1 of<a class="footnote-ref" id="fnref:7" href="#fn:7"><sup>7</sup></a> for a derivation of the distribution q j ∗ ( Z j ∣ X ) {\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )} . </p><p>In practice, we usually work in terms of logarithms, i.e.: </p> ln q j ∗ ( Z j ∣ X ) = E q − j ∗ [ ln p ( Z , X ) ] + constant {\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}} <p>The constant in the above expression is related to the <a href="/facts/Normalizing_constant/KKTfybqr">normalizing constant</a> (the denominator in the expression above for q j ∗ {\displaystyle q_{j}^{*}} ) and is usually reinstated by inspection, as the rest of the expression can usually be recognized as being a known type of distribution (e.g. <a href="/facts/Gaussian_distribution/UapjjPyQ">Gaussian</a>, <a href="/facts/Gamma_distribution/lczcdJmw">gamma</a>, etc.). </p><p>Using the properties of expectations, the expression E q − j ∗ [ ln p ( Z , X ) ] {\displaystyle \operatorname {E} _{q_{-j}^{*}}[\ln p(\mathbf {Z} ,\mathbf {X} )]} can usually be simplified into a function of the fixed <a href="/facts/Hyperparameter_(Bayesian_statistics)/zqcLPQcb">hyperparameters</a> of the <a href="/facts/Prior_distribution/JQKAD4o0">prior distributions</a> over the latent variables and of expectations (and sometimes higher <a href="/facts/Moment_(mathematics)/GL7vJ1nb">moments</a> such as the <a href="/facts/Variance/ULBJKXD1">variance</a>) of latent variables not in the current partition (i.e. latent variables not included in Z j {\displaystyle \mathbf {Z} _{j}} ). This creates <a href="/facts/Circular_dependency/l6mh3qVC">circular dependencies</a> between the parameters of the distributions over variables in one partition and the expectations of variables in the other partitions. This naturally suggests an <a href="/facts/Iterative/AMc1p7BL">iterative</a> algorithm, much like EM (the <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">expectation–maximization algorithm</a>), in which the expectations (and possibly higher moments) of the latent variables are initialized in some fashion (perhaps randomly), and then the parameters of each distribution are computed in turn using the current values of the expectations, after which the expectation of the newly computed distribution is set appropriately according to the computed parameters. An algorithm of this sort is guaranteed to <a href="/facts/Limit_of_a_sequence/Ktge2RwL">converge</a>.<a class="footnote-ref" id="fnref:8" href="#fn:8"><sup>8</sup></a> </p><p>In other words, for each of the partitions of variables, by simplifying the expression for the distribution over the partition's variables and examining the distribution's functional dependency on the variables in question, the family of the distribution can usually be determined (which in turn determines the value of the constant). The formula for the distribution's parameters will be expressed in terms of the prior distributions' hyperparameters (which are known constants), but also in terms of expectations of functions of variables in other partitions. Usually these expectations can be simplified into functions of expectations of the variables themselves (i.e. the <a href="/facts/Mean/swcEd4Pg">means</a>); sometimes expectations of squared variables (which can be related to the <a href="/facts/Variance/ULBJKXD1">variance</a> of the variables), or expectations of higher powers (i.e. higher <a href="/facts/Moment_(mathematics)/GL7vJ1nb">moments</a>) also appear. In most cases, the other variables' distributions will be from known families, and the formulas for the relevant expectations can be looked up. However, those formulas depend on those distributions' parameters, which depend in turn on the expectations about other variables. The result is that the formulas for the parameters of each variable's distributions can be expressed as a series of equations with mutual, <a href="/facts/Nonlinear/4aYItALC">nonlinear</a> dependencies among the variables. Usually, it is not possible to solve this system of equations directly. However, as described above, the dependencies suggest a simple iterative algorithm, which in most cases is guaranteed to converge. An example will make this process clearer. </p> <h2 id="a-duality-formula-for-variational-inference">A duality formula for variational inference</h2> <p>The following theorem is referred to as a duality formula for variational inference.<a class="footnote-ref" id="fnref:9" href="#fn:9"><sup>9</sup></a> It explains some important properties of the variational distributions used in variational Bayes methods. </p><p>Theorem Consider two <a href="/facts/Probability_spaces/dEcv2VrU">probability spaces</a> ( Θ , F , P ) {\displaystyle (\Theta ,{\mathcal {F}},P)} and ( Θ , F , Q ) {\displaystyle (\Theta ,{\mathcal {F}},Q)} with Q ≪ P {\displaystyle Q\ll P} . Assume that there is a common dominating <a href="/facts/Probability_measure/8BjhgoPR">probability measure</a> λ {\displaystyle \lambda } such that P ≪ λ {\displaystyle P\ll \lambda } and Q ≪ λ {\displaystyle Q\ll \lambda } . Let h {\displaystyle h} denote any real-valued <a href="/facts/Random_variable/TwTBXnLT">random variable</a> on ( Θ , F , P ) {\displaystyle (\Theta ,{\mathcal {F}},P)} that satisfies h ∈ L 1 ( P ) {\displaystyle h\in L_{1}(P)} . Then the following equality holds </p> log E P [ exp h ] = sup Q ≪ P { E Q [ h ] − D KL ( Q ∥ P ) } . {\displaystyle \log E_{P}[\exp h]={\text{sup}}_{Q\ll P}\{E_{Q}[h]-D_{\text{KL}}(Q\parallel P)\}.} <p>Further, the supremum on the right-hand side is attained <a href="/facts/If_and_only_if/bYSxGJ66">if and only if</a> it holds </p> q ( θ ) p ( θ ) = exp h ( θ ) E P [ exp h ] , {\displaystyle {\frac {q(\theta )}{p(\theta )}}={\frac {\exp h(\theta )}{E_{P}[\exp h]}},} <p>almost surely with respect to probability measure Q {\displaystyle Q} , where p ( θ ) = d P / d λ {\displaystyle p(\theta )=dP/d\lambda } and q ( θ ) = d Q / d λ {\displaystyle q(\theta )=dQ/d\lambda } denote the Radon–Nikodym derivatives of the probability measures P {\displaystyle P} and Q {\displaystyle Q} with respect to λ {\displaystyle \lambda } , respectively. </p> <h2 id="a-basic-example">A basic example</h2> <p>Consider a simple non-hierarchical Bayesian model consisting of a set of <a href="/facts/Independent_identically_distributed/othIRaWt">i.i.d.</a> observations from a <a href="/facts/Gaussian_distribution/UapjjPyQ">Gaussian distribution</a>, with unknown <a href="/facts/Mean/swcEd4Pg">mean</a> and <a href="/facts/Variance/ULBJKXD1">variance</a>.<a class="footnote-ref" id="fnref:10" href="#fn:10"><sup>10</sup></a> In the following, we work through this model in great detail to illustrate the workings of the variational Bayes method. </p><p>For mathematical convenience, in the following example we work in terms of the <a href="/facts/Precision_(statistics)/xsMYkkjA">precision</a> — i.e. the reciprocal of the variance (or in a multivariate Gaussian, the inverse of the <a href="/facts/Covariance_matrix/kLhgjHC6">covariance matrix</a>) — rather than the variance itself. (From a theoretical standpoint, precision and variance are equivalent since there is a <a href="/facts/One-to-one_correspondence/j4gTuTmW">one-to-one correspondence</a> between the two.) </p> <h3>The mathematical model</h3> <p>We place <a href="/facts/Conjugate_prior/ScCFcs8b">conjugate prior</a> distributions on the unknown mean μ {\displaystyle \mu } and precision τ {\displaystyle \tau } , i.e. the mean also follows a Gaussian distribution while the precision follows a <a href="/facts/Gamma_distribution/lczcdJmw">gamma distribution</a>. In other words: </p> τ ∼ Gamma ( a 0 , b 0 ) μ | τ ∼ N ( μ 0 , ( λ 0 τ ) − 1 ) { x 1 , … , x N } ∼ N ( μ , τ − 1 ) N = number of data points {\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}}} <p>The <a href="/facts/Hyperparameter_(Bayesian_statistics)/zqcLPQcb">hyperparameters</a> μ 0 , λ 0 , a 0 {\displaystyle \mu _{0},\lambda _{0},a_{0}} and b 0 {\displaystyle b_{0}} in the prior distributions are fixed, given values. They can be set to small positive numbers to give broad prior distributions indicating ignorance about the prior distributions of μ {\displaystyle \mu } and τ {\displaystyle \tau } . </p><p>We are given N {\displaystyle N} data points X = { x 1 , … , x N } {\displaystyle \mathbf {X} =\{x_{1},\ldots ,x_{N}\}} and our goal is to infer the <a href="/facts/Posterior_distribution/fojsMD0D">posterior distribution</a> q ( μ , τ ) = p ( μ , τ ∣ x 1 , … , x N ) {\displaystyle q(\mu ,\tau )=p(\mu ,\tau \mid x_{1},\ldots ,x_{N})} of the parameters μ {\displaystyle \mu } and τ . {\displaystyle \tau .} </p> <h3>The joint probability</h3> <p>The <a href="/facts/Joint_probability/klX2ksGY">joint probability</a> of all variables can be rewritten as </p> p ( X , μ , τ ) = p ( X ∣ μ , τ ) p ( μ ∣ τ ) p ( τ ) {\displaystyle p(\mathbf {X} ,\mu ,\tau )=p(\mathbf {X} \mid \mu ,\tau )p(\mu \mid \tau )p(\tau )} <p>where the individual factors are </p> p ( X ∣ μ , τ ) = ∏ n = 1 N N ( x n ∣ μ , τ − 1 ) p ( μ ∣ τ ) = N ( μ ∣ μ 0 , ( λ 0 τ ) − 1 ) p ( τ ) = Gamma ( τ ∣ a 0 , b 0 ) {\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}}} <p>where </p> N ( x ∣ μ , σ 2 ) = 1 2 π σ 2 e − ( x − μ ) 2 2 σ 2 Gamma ( τ ∣ a , b ) = 1 Γ ( a ) b a τ a − 1 e − b τ {\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}}} <h3>Factorized approximation</h3> <p>Assume that q ( μ , τ ) = q ( μ ) q ( τ ) {\displaystyle q(\mu ,\tau )=q(\mu )q(\tau )} , i.e. that the posterior distribution factorizes into independent factors for μ {\displaystyle \mu } and τ {\displaystyle \tau } . This type of assumption underlies the variational Bayesian method. The true posterior distribution does not in fact factor this way (in fact, in this simple case, it is known to be a <a href="/facts/Gaussian-gamma_distribution/cW9fI2D4">Gaussian-gamma distribution</a>), and hence the result we obtain will be an approximation. </p> <h3>Derivation of <i>q</i>(<i>μ</i>)</h3> <p>Then </p> ln q μ ∗ ( μ ) = E τ [ ln p ( X ∣ μ , τ ) + ln p ( μ ∣ τ ) + ln p ( τ ) ] + C = E τ [ ln p ( X ∣ μ , τ ) ] + E τ [ ln p ( μ ∣ τ ) ] + E τ [ ln p ( τ ) ] + C = E τ [ ln ∏ n = 1 N N ( x n ∣ μ , τ − 1 ) ] + E τ [ ln N ( μ ∣ μ 0 , ( λ 0 τ ) − 1 ) ] + C 2 = E τ [ ln ∏ n = 1 N τ 2 π e − ( x n − μ ) 2 τ 2 ] + E τ [ ln λ 0 τ 2 π e − ( μ − μ 0 ) 2 λ 0 τ 2 ] + C 2 = E τ [ ∑ n = 1 N ( 1 2 ( ln τ − ln 2 π ) − ( x n − μ ) 2 τ 2 ) ] + E τ [ 1 2 ( ln λ 0 + ln τ − ln 2 π ) − ( μ − μ 0 ) 2 λ 0 τ 2 ] + C 2 = E τ [ ∑ n = 1 N − ( x n − μ ) 2 τ 2 ] + E τ [ − ( μ − μ 0 ) 2 λ 0 τ 2 ] + E τ [ ∑ n = 1 N 1 2 ( ln τ − ln 2 π ) ] + E τ [ 1 2 ( ln λ 0 + ln τ − ln 2 π ) ] + C 2 = E τ [ ∑ n = 1 N − ( x n − μ ) 2 τ 2 ] + E τ [ − ( μ − μ 0 ) 2 λ 0 τ 2 ] + C 3 = − E τ [ τ ] 2 { ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 } + C 3 {\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}}} <p>In the above derivation, C {\displaystyle C} , C 2 {\displaystyle C_{2}} and C 3 {\displaystyle C_{3}} refer to values that are constant with respect to μ {\displaystyle \mu } . Note that the term E τ [ ln p ( τ ) ] {\displaystyle \operatorname {E} _{\tau }[\ln p(\tau )]} is not a function of μ {\displaystyle \mu } and will have the same value regardless of the value of μ {\displaystyle \mu } . Hence in line 3 we can absorb it into the <a href="/facts/Constant_term/KRFu0aum">constant term</a> at the end. We do the same thing in line 7. </p><p>The last line is simply a quadratic polynomial in μ {\displaystyle \mu } . Since this is the logarithm of q μ ∗ ( μ ) {\displaystyle q_{\mu }^{*}(\mu )} , we can see that q μ ∗ ( μ ) {\displaystyle q_{\mu }^{*}(\mu )} itself is a <a href="/facts/Gaussian_distribution/UapjjPyQ">Gaussian distribution</a>. </p><p>With a certain amount of tedious math (expanding the squares inside of the braces, separating out and grouping the terms involving μ {\displaystyle \mu } and μ 2 {\displaystyle \mu ^{2}} and <a href="/facts/Completing_the_square/SN6xA8sc">completing the square</a> over μ {\displaystyle \mu } ), we can derive the parameters of the Gaussian distribution: </p> ln q μ ∗ ( μ ) = − E τ [ τ ] 2 { ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 } + C 3 = − E τ [ τ ] 2 { ∑ n = 1 N ( x n 2 − 2 x n μ + μ 2 ) + λ 0 ( μ 2 − 2 μ 0 μ + μ 0 2 ) } + C 3 = − E τ [ τ ] 2 { ( ∑ n = 1 N x n 2 ) − 2 ( ∑ n = 1 N x n ) μ + ( ∑ n = 1 N μ 2 ) + λ 0 μ 2 − 2 λ 0 μ 0 μ + λ 0 μ 0 2 } + C 3 = − E τ [ τ ] 2 { ( λ 0 + N ) μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) μ + ( ∑ n = 1 N x n 2 ) + λ 0 μ 0 2 } + C 3 = − E τ [ τ ] 2 { ( λ 0 + N ) μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) μ } + C 4 = − E τ [ τ ] 2 { ( λ 0 + N ) μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) ( λ 0 + N ) μ } + C 4 = − E τ [ τ ] 2 { ( λ 0 + N ) ( μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) μ ) } + C 4 = − E τ [ τ ] 2 { ( λ 0 + N ) ( μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) μ + ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) 2 − ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) 2 ) } + C 4 = − E τ [ τ ] 2 { ( λ 0 + N ) ( μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) μ + ( λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) 2 ) } + C 5 = − E τ [ τ ] 2 { ( λ 0 + N ) ( μ − λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) 2 } + C 5 = − 1 2 ( λ 0 + N ) E τ [ τ ] ( μ − λ 0 μ 0 + ∑ n = 1 N x n λ 0 + N ) 2 + C 5 {\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}}} <p>Note that all of the above steps can be shortened by using the formula for the <a href="/facts/Normal_distribution/UapjjPyQ">sum of two quadratics</a>. </p><p>In other words: </p> q μ ∗ ( μ ) ∼ N ( μ ∣ μ N , λ N − 1 ) μ N = λ 0 μ 0 + N x ¯ λ 0 + N λ N = ( λ 0 + N ) E τ [ τ ] x ¯ = 1 N ∑ n = 1 N x n {\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}}} <h3>Derivation of q(τ)</h3> <p>The derivation of q τ ∗ ( τ ) {\displaystyle q_{\tau }^{*}(\tau )} is similar to above, although we omit some of the details for the sake of brevity. </p> ln q τ ∗ ( τ ) = E μ [ ln p ( X ∣ μ , τ ) + ln p ( μ ∣ τ ) ] + ln p ( τ ) + constant = ( a 0 − 1 ) ln τ − b 0 τ + 1 2 ln τ + N 2 ln τ − τ 2 E μ [ ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 ] + constant {\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}}} <p>Exponentiating both sides, we can see that q τ ∗ ( τ ) {\displaystyle q_{\tau }^{*}(\tau )} is a <a href="/facts/Gamma_distribution/lczcdJmw">gamma distribution</a>. Specifically: </p> q τ ∗ ( τ ) ∼ Gamma ( τ ∣ a N , b N ) a N = a 0 + N + 1 2 b N = b 0 + 1 2 E μ [ ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 ] {\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}}} <h3>Algorithm for computing the parameters</h3> <p>Let us recap the conclusions from the previous sections: </p> q μ ∗ ( μ ) ∼ N ( μ ∣ μ N , λ N − 1 ) μ N = λ 0 μ 0 + N x ¯ λ 0 + N λ N = ( λ 0 + N ) E τ [ τ ] x ¯ = 1 N ∑ n = 1 N x n {\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}}} <p>and </p> q τ ∗ ( τ ) ∼ Gamma ( τ ∣ a N , b N ) a N = a 0 + N + 1 2 b N = b 0 + 1 2 E μ [ ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 ] {\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}}} <p>In each case, the parameters for the distribution over one of the variables depend on expectations taken with respect to the other variable. We can expand the expectations, using the standard formulas for the expectations of moments of the Gaussian and gamma distributions: </p> E [ τ ∣ a N , b N ] = a N b N E [ μ ∣ μ N , λ N − 1 ] = μ N E [ X 2 ] = Var ( X ) + ( E [ X ] ) 2 E [ μ 2 ∣ μ N , λ N − 1 ] = λ N − 1 + μ N 2 {\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}}} <p>Applying these formulas to the above equations is trivial in most cases, but the equation for b N {\displaystyle b_{N}} takes more work: </p> b N = b 0 + 1 2 E μ [ ∑ n = 1 N ( x n − μ ) 2 + λ 0 ( μ − μ 0 ) 2 ] = b 0 + 1 2 E μ [ ( λ 0 + N ) μ 2 − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) μ + ( ∑ n = 1 N x n 2 ) + λ 0 μ 0 2 ] = b 0 + 1 2 [ ( λ 0 + N ) E μ [ μ 2 ] − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) E μ [ μ ] + ( ∑ n = 1 N x n 2 ) + λ 0 μ 0 2 ] = b 0 + 1 2 [ ( λ 0 + N ) ( λ N − 1 + μ N 2 ) − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) μ N + ( ∑ n = 1 N x n 2 ) + λ 0 μ 0 2 ] {\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}}} <p>We can then write the parameter equations as follows, without any expectations: </p> μ N = λ 0 μ 0 + N x ¯ λ 0 + N λ N = ( λ 0 + N ) a N b N x ¯ = 1 N ∑ n = 1 N x n a N = a 0 + N + 1 2 b N = b 0 + 1 2 [ ( λ 0 + N ) ( λ N − 1 + μ N 2 ) − 2 ( λ 0 μ 0 + ∑ n = 1 N x n ) μ N + ( ∑ n = 1 N x n 2 ) + λ 0 μ 0 2 ] {\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}}} <p>Note that there are circular dependencies among the formulas for λ N {\displaystyle \lambda _{N}} and b N {\displaystyle b_{N}} . This naturally suggests an <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">EM</a>-like algorithm: </p> <ol><li>Compute ∑ n = 1 N x n {\displaystyle \sum _{n=1}^{N}x_{n}} and ∑ n = 1 N x n 2 . {\displaystyle \sum _{n=1}^{N}x_{n}^{2}.} Use these values to compute μ N {\displaystyle \mu _{N}} and a N . {\displaystyle a_{N}.} </li> <li>Initialize λ N {\displaystyle \lambda _{N}} to some arbitrary value.</li> <li>Use the current value of λ N , {\displaystyle \lambda _{N},} along with the known values of the other parameters, to compute b N {\displaystyle b_{N}} .</li> <li>Use the current value of b N , {\displaystyle b_{N},} along with the known values of the other parameters, to compute λ N {\displaystyle \lambda _{N}} .</li> <li>Repeat the last two steps until convergence (i.e. until neither value has changed more than some small amount).</li></ol> <p>We then have values for the hyperparameters of the approximating distributions of the posterior parameters, which we can use to compute any properties we want of the posterior — e.g. its mean and variance, a 95% highest-density region (the smallest interval that includes 95% of the total probability), etc. </p><p>It can be shown that this algorithm is guaranteed to converge to a local maximum. </p><p>Note also that the posterior distributions have the same form as the corresponding prior distributions. We did <i>not</i> assume this; the only assumption we made was that the distributions factorize, and the form of the distributions followed naturally. It turns out (see below) that the fact that the posterior distributions have the same form as the prior distributions is not a coincidence, but a general result whenever the prior distributions are members of the <a href="/facts/Exponential_family/1LkkqEIf">exponential family</a>, which is the case for most of the standard distributions. </p> <h2 id="further-discussion">Further discussion</h2> <h3>Step-by-step recipe</h3> <p>The above example shows the method by which the variational-Bayesian approximation to a <a href="/facts/Posterior_probability/fojsMD0D">posterior probability</a> density in a given <a href="/facts/Bayesian_network/Ry31vq6T">Bayesian network</a> is derived: </p> <ol><li>Describe the network with a <a href="/facts/Graphical_model/XxfmKhmM">graphical model</a>, identifying the observed variables (data) X {\displaystyle \mathbf {X} } and unobserved variables (<a href="/facts/Parameter/zFvVFMv9">parameters</a> Θ {\displaystyle {\boldsymbol {\Theta }}} and <a href="/facts/Latent_variable/ohg99BnW">latent variables</a> Z {\displaystyle \mathbf {Z} } ) and their <a href="/facts/Conditional_probability_distribution/0eGm3P9W">conditional probability distributions</a>. Variational Bayes will then construct an approximation to the posterior probability p ( Z , Θ ∣ X ) {\displaystyle p(\mathbf {Z} ,{\boldsymbol {\Theta }}\mid \mathbf {X} )} . The approximation has the basic property that it is a factorized distribution, i.e. a product of two or more <a href="/facts/Statistical_independence/NUzQtnUL">independent</a> distributions over disjoint subsets of the unobserved variables.</li> <li>Partition the unobserved variables into two or more subsets, over which the independent factors will be derived. There is no universal procedure for doing this; creating too many subsets yields a poor approximation, while creating too few makes the entire variational Bayes procedure intractable. Typically, the first split is to separate the parameters and latent variables; often, this is enough by itself to produce a tractable result. Assume that the partitions are called Z 1 , … , Z M {\displaystyle \mathbf {Z} _{1},\ldots ,\mathbf {Z} _{M}} .</li> <li>For a given partition Z j {\displaystyle \mathbf {Z} _{j}} , write down the formula for the best approximating distribution q j ∗ ( Z j ∣ X ) {\displaystyle q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )} using the basic equation ln q j ∗ ( Z j ∣ X ) = E i ≠ j [ ln p ( Z , X ) ] + constant {\displaystyle \ln q_{j}^{*}(\mathbf {Z} _{j}\mid \mathbf {X} )=\operatorname {E} _{i\neq j}[\ln p(\mathbf {Z} ,\mathbf {X} )]+{\text{constant}}} .</li> <li>Fill in the formula for the <a href="/facts/Joint_probability_distribution/klX2ksGY">joint probability distribution</a> using the graphical model. Any component conditional distributions that don't involve any of the variables in Z j {\displaystyle \mathbf {Z} _{j}} can be ignored; they will be folded into the constant term.</li> <li>Simplify the formula and apply the expectation operator, following the above example. Ideally, this should simplify into expectations of basic functions of variables not in Z j {\displaystyle \mathbf {Z} _{j}} (e.g. first or second raw <a href="/facts/Moment_(mathematics)/GL7vJ1nb">moments</a>, expectation of a logarithm, etc.). In order for the variational Bayes procedure to work well, these expectations should generally be expressible analytically as functions of the parameters and/or <a href="/facts/Hyperparameter_(Bayesian_statistics)/zqcLPQcb">hyperparameters</a> of the distributions of these variables. In all cases, these expectation terms are constants with respect to the variables in the current partition.</li> <li>The functional form of the formula with respect to the variables in the current partition indicates the type of distribution. In particular, exponentiating the formula generates the <a href="/facts/Probability_density_function/zvfybna4">probability density function</a> (PDF) of the distribution (or at least, something proportional to it, with unknown <a href="/facts/Normalization_constant/KKTfybqr">normalization constant</a>). In order for the overall method to be tractable, it should be possible to recognize the functional form as belonging to a known distribution. Significant mathematical manipulation may be required to convert the formula into a form that matches the PDF of a known distribution. When this can be done, the normalization constant can be reinstated by definition, and equations for the parameters of the known distribution can be derived by extracting the appropriate parts of the formula.</li> <li>When all expectations can be replaced analytically with functions of variables not in the current partition, and the PDF put into a form that allows identification with a known distribution, the result is a set of equations expressing the values of the optimum parameters as functions of the parameters of variables in other partitions.</li> <li>When this procedure can be applied to all partitions, the result is a set of mutually linked equations specifying the optimum values of all parameters.</li> <li>An <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">expectation–maximization</a> (EM) type procedure is then applied, picking an initial value for each parameter and the iterating through a series of steps, where at each step we cycle through the equations, updating each parameter in turn. This is guaranteed to converge.</li></ol> <h3>Most important points</h3> <p>Due to all of the mathematical manipulations involved, it is easy to lose track of the big picture. The important things are: </p> <ol><li>The idea of variational Bayes is to construct an analytical approximation to the <a href="/facts/Posterior_probability/fojsMD0D">posterior probability</a> of the set of unobserved variables (parameters and latent variables), given the data. This means that the form of the solution is similar to other <a href="/facts/Bayesian_inference/tEpK3zLx">Bayesian inference</a> methods, such as <a href="/facts/Gibbs_sampling/HDYRhxEq">Gibbs sampling</a> — i.e. a distribution that seeks to describe everything that is known about the variables. As in other Bayesian methods — but unlike e.g. in <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">expectation–maximization</a> (EM) or other <a href="/facts/Maximum_likelihood/0Yq2dpQD">maximum likelihood</a> methods — both types of unobserved variables (i.e. parameters and latent variables) are treated the same, i.e. as <a href="/facts/Random_variable/TwTBXnLT">random variables</a>. Estimates for the variables can then be derived in the standard Bayesian ways, e.g. calculating the mean of the distribution to get a single point estimate or deriving a <a href="/facts/Credible_interval/BeVTjGrm">credible interval</a>, highest density region, etc.</li> <li>"Analytical approximation" means that a formula can be written down for the posterior distribution. The formula generally consists of a product of well-known probability distributions, each of which <i>factorizes</i> over a set of unobserved variables (i.e. it is <a href="/facts/Conditionally_independent/udUXFaMw">conditionally independent</a> of the other variables, given the observed data). This formula is not the true posterior distribution, but an approximation to it; in particular, it will generally agree fairly closely in the lowest <a href="/facts/Moment_(mathematics)/GL7vJ1nb">moments</a> of the unobserved variables, e.g. the <a href="/facts/Mean/swcEd4Pg">mean</a> and <a href="/facts/Variance/ULBJKXD1">variance</a>.</li> <li>The result of all of the mathematical manipulations is (1) the identity of the probability distributions making up the factors, and (2) mutually dependent formulas for the parameters of these distributions. The actual values of these parameters are computed numerically, through an alternating iterative procedure much like EM.</li></ol> <h3>Compared with expectation–maximization (EM)</h3> <p>Variational Bayes (VB) is often compared with <a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">expectation–maximization</a> (EM). The actual numerical procedure is quite similar, in that both are alternating iterative procedures that successively converge on optimum parameter values. The initial steps to derive the respective procedures are also vaguely similar, both starting out with formulas for probability densities and both involving significant amounts of mathematical manipulations. </p><p>However, there are a number of differences. Most important is <i>what</i> is being computed. </p> <ul><li>EM computes point estimates of posterior distribution of those random variables that can be categorized as "parameters", but only estimates of the actual posterior distributions of the latent variables (at least in "soft EM", and often only when the latent variables are discrete). The point estimates computed are the <a href="/facts/Mode_(statistics)/vq5tGp0G">modes</a> of these parameters; no other information is available.</li> <li>VB, on the other hand, computes estimates of the actual posterior distribution of all variables, both parameters and latent variables. When point estimates need to be derived, generally the <a href="/facts/Mean/swcEd4Pg">mean</a> is used rather than the mode, as is normal in Bayesian inference. Concomitant with this, the parameters computed in VB do <i>not</i> have the same significance as those in EM. EM computes optimum values of the parameters of the Bayes network itself. VB computes optimum values of the parameters of the distributions used to approximate the parameters and latent variables of the Bayes network. For example, a typical Gaussian <a href="/facts/Mixture_model/WWe6TJlg">mixture model</a> will have parameters for the mean and variance of each of the mixture components. EM would directly estimate optimum values for these parameters. VB, however, would first fit a distribution to these parameters — typically in the form of a <a href="/facts/Prior_distribution/JQKAD4o0">prior distribution</a>, e.g. a <a href="/facts/Normal-scaled_inverse_gamma_distribution/2koIxNwD">normal-scaled inverse gamma distribution</a> — and would then compute values for the parameters of this prior distribution, i.e. essentially <a href="/facts/Hyperparameter_(Bayesian_statistics)/zqcLPQcb">hyperparameters</a>. In this case, VB would compute optimum estimates of the four parameters of the normal-scaled inverse gamma distribution that describes the joint distribution of the mean and variance of the component.</li></ul> <h2 id="a-more-complex-example">A more complex example</h2> <p>Imagine a Bayesian <a href="/facts/Gaussian_mixture_model/WWe6TJlg">Gaussian mixture model</a> described as follows:<a class="footnote-ref" id="fnref:11" href="#fn:11"><sup>11</sup></a> </p> π ∼ SymDir ( K , α 0 ) Λ i = 1 … K ∼ W ( W 0 , ν 0 ) μ i = 1 … K ∼ N ( μ 0 , ( β 0 Λ i ) − 1 ) z [ i = 1 … N ] ∼ Mult ( 1 , π ) x i = 1 … N ∼ N ( μ z i , Λ z i − 1 ) K = number of mixing components N = number of data points {\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}}} <p>Note: </p> <ul><li>SymDir() is the symmetric <a href="/facts/Dirichlet_distribution/Qq13f5g2">Dirichlet distribution</a> of dimension K {\displaystyle K} , with the hyperparameter for each component set to α 0 {\displaystyle \alpha _{0}} . The Dirichlet distribution is the <a href="/facts/Conjugate_prior/ScCFcs8b">conjugate prior</a> of the <a href="/facts/Categorical_distribution/xpc6RWM2">categorical distribution</a> or <a href="/facts/Multinomial_distribution/y58v1p9J">multinomial distribution</a>.</li> <li> W ( ) {\displaystyle {\mathcal {W}}()} is the <a href="/facts/Wishart_distribution/MEcqTump">Wishart distribution</a>, which is the conjugate prior of the <a href="/facts/Precision_matrix/xsMYkkjA">precision matrix</a> (inverse <a href="/facts/Covariance_matrix/kLhgjHC6">covariance matrix</a>) for a <a href="/facts/Multivariate_Gaussian_distribution/2Xfegqz2">multivariate Gaussian distribution</a>.</li> <li>Mult() is a <a href="/facts/Multinomial_distribution/y58v1p9J">multinomial distribution</a> over a single observation (equivalent to a <a href="/facts/Categorical_distribution/xpc6RWM2">categorical distribution</a>). The state space is a "one-of-K" representation, i.e., a K {\displaystyle K} -dimensional vector in which one of the elements is 1 (specifying the identity of the observation) and all other elements are 0.</li> <li> N ( ) {\displaystyle {\mathcal {N}}()} is the <a href="/facts/Gaussian_distribution/UapjjPyQ">Gaussian distribution</a>, in this case specifically the <a href="/facts/Multivariate_Gaussian_distribution/2Xfegqz2">multivariate Gaussian distribution</a>.</li></ul> <p>The interpretation of the above variables is as follows: </p> <ul><li> X = { x 1 , … , x N } {\displaystyle \mathbf {X} =\{\mathbf {x} _{1},\dots ,\mathbf {x} _{N}\}} is the set of N {\displaystyle N} data points, each of which is a D {\displaystyle D} -dimensional vector distributed according to a <a href="/facts/Multivariate_Gaussian_distribution/2Xfegqz2">multivariate Gaussian distribution</a>.</li> <li> Z = { z 1 , … , z N } {\displaystyle \mathbf {Z} =\{\mathbf {z} _{1},\dots ,\mathbf {z} _{N}\}} is a set of latent variables, one per data point, specifying which mixture component the corresponding data point belongs to, using a "one-of-K" vector representation with components z n k {\displaystyle z_{nk}} for k = 1 … K {\displaystyle k=1\dots K} , as described above.</li> <li> π {\displaystyle \mathbf {\pi } } is the mixing proportions for the K {\displaystyle K} mixture components.</li> <li> μ i = 1 … K {\displaystyle \mathbf {\mu } _{i=1\dots K}} and Λ i = 1 … K {\displaystyle \mathbf {\Lambda } _{i=1\dots K}} specify the parameters (<a href="/facts/Mean/swcEd4Pg">mean</a> and <a href="/facts/Precision_(statistics)/xsMYkkjA">precision</a>) associated with each mixture component.</li></ul> <p>The joint probability of all variables can be rewritten as </p> p ( X , Z , π , μ , Λ ) = p ( X ∣ Z , μ , Λ ) p ( Z ∣ π ) p ( π ) p ( μ ∣ Λ ) p ( Λ ) {\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 } )} <p>where the individual factors are </p> p ( X ∣ Z , μ , Λ ) = ∏ n = 1 N ∏ k = 1 K N ( x n ∣ μ k , Λ k − 1 ) z n k p ( Z ∣ π ) = ∏ n = 1 N ∏ k = 1 K π k z n k p ( π ) = Γ ( K α 0 ) Γ ( α 0 ) K ∏ k = 1 K π k α 0 − 1 p ( μ ∣ Λ ) = ∏ k = 1 K N ( μ k ∣ μ 0 , ( β 0 Λ k ) − 1 ) p ( Λ ) = ∏ k = 1 K W ( Λ k ∣ W 0 , ν 0 ) {\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}}} <p>where </p> N ( x ∣ μ , Σ ) = 1 ( 2 π ) D / 2 1 | Σ | 1 / 2 exp { − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) } W ( Λ ∣ W , ν ) = B ( W , ν ) | Λ | ( ν − D − 1 ) / 2 exp ( − 1 2 Tr ( W − 1 Λ ) ) B ( W , ν ) = | W | − ν / 2 { 2 ν D / 2 π D ( D − 1 ) / 4 ∏ i = 1 D Γ ( ν + 1 − i 2 ) } − 1 D = dimensionality of each data point {\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}}} <p>Assume that q ( Z , π , μ , Λ ) = q ( Z ) q ( π , μ , Λ ) {\displaystyle q(\mathbf {Z} ,\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )=q(\mathbf {Z} )q(\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )} . </p><p>Then<a class="footnote-ref" id="fnref:12" href="#fn:12"><sup>12</sup></a> </p> ln q ∗ ( Z ) = E π , μ , Λ [ ln p ( X , Z , π , μ , Λ ) ] + constant = E π [ ln p ( Z ∣ π ) ] + E μ , Λ [ ln p ( X ∣ Z , μ , Λ ) ] + constant = ∑ n = 1 N ∑ k = 1 K z n k ln ρ n k + constant {\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}}} <p>where we have defined </p> ln ρ n k = E [ ln π k ] + 1 2 E [ ln | Λ k | ] − D 2 ln ( 2 π ) − 1 2 E μ k , Λ k [ ( x n − μ k ) T Λ k ( x n − μ k ) ] {\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})]} <p>Exponentiating both sides of the formula for ln q ∗ ( Z ) {\displaystyle \ln q^{*}(\mathbf {Z} )} yields </p> q ∗ ( Z ) ∝ ∏ n = 1 N ∏ k = 1 K ρ n k z n k {\displaystyle q^{*}(\mathbf {Z} )\propto \prod _{n=1}^{N}\prod _{k=1}^{K}\rho _{nk}^{z_{nk}}} <p>Requiring that this be normalized ends up requiring that the ρ n k {\displaystyle \rho _{nk}} sum to 1 over all values of k {\displaystyle k} , yielding </p> q ∗ ( Z ) = ∏ n = 1 N ∏ k = 1 K r n k z n k {\displaystyle q^{*}(\mathbf {Z} )=\prod _{n=1}^{N}\prod _{k=1}^{K}r_{nk}^{z_{nk}}} <p>where </p> r n k = ρ n k ∑ j = 1 K ρ n j {\displaystyle r_{nk}={\frac {\rho _{nk}}{\sum _{j=1}^{K}\rho _{nj}}}} <p>In other words, q ∗ ( Z ) {\displaystyle q^{*}(\mathbf {Z} )} is a product of single-observation <a href="/facts/Multinomial_distribution/y58v1p9J">multinomial distributions</a>, and factors over each individual z n {\displaystyle \mathbf {z} _{n}} , which is distributed as a single-observation multinomial distribution with parameters r n k {\displaystyle r_{nk}} for k = 1 … K {\displaystyle k=1\dots K} . </p><p>Furthermore, we note that </p> E [ z n k ] = r n k {\displaystyle \operatorname {E} [z_{nk}]=r_{nk}\,} <p>which is a standard result for categorical distributions. </p><p>Now, considering the factor q ( π , μ , Λ ) {\displaystyle q(\mathbf {\pi } ,\mathbf {\mu } ,\mathbf {\Lambda } )} , note that it automatically factors into q ( π ) ∏ k = 1 K q ( μ k , Λ k ) {\displaystyle q(\mathbf {\pi } )\prod _{k=1}^{K}q(\mathbf {\mu } _{k},\mathbf {\Lambda } _{k})} due to the structure of the graphical model defining our Gaussian mixture model, which is specified above. </p><p>Then, </p> ln q ∗ ( π ) = ln p ( π ) + E Z [ ln p ( Z ∣ π ) ] + constant = ( α 0 − 1 ) ∑ k = 1 K ln π k + ∑ n = 1 N ∑ k = 1 K r n k ln π k + constant {\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}}} <p>Taking the exponential of both sides, we recognize q ∗ ( π ) {\displaystyle q^{*}(\mathbf {\pi } )} as a <a href="/facts/Dirichlet_distribution/Qq13f5g2">Dirichlet distribution</a> </p> q ∗ ( π ) ∼ Dir ( α ) {\displaystyle q^{*}(\mathbf {\pi } )\sim \operatorname {Dir} (\mathbf {\alpha } )\,} <p>where </p> α k = α 0 + N k {\displaystyle \alpha _{k}=\alpha _{0}+N_{k}\,} <p>where </p> N k = ∑ n = 1 N r n k {\displaystyle N_{k}=\sum _{n=1}^{N}r_{nk}\,} <p>Finally </p> ln q ∗ ( μ k , Λ k ) = ln p ( μ k , Λ k ) + ∑ n = 1 N E [ z n k ] ln N ( x n ∣ μ k , Λ k − 1 ) + constant {\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}}} <p>Grouping and reading off terms involving μ k {\displaystyle \mathbf {\mu } _{k}} and Λ k {\displaystyle \mathbf {\Lambda } _{k}} , the result is a <a href="/facts/Gaussian-Wishart_distribution/aA5kTsR1">Gaussian-Wishart distribution</a> given by </p> q ∗ ( μ k , Λ k ) = N ( μ k ∣ m k , ( β k Λ k ) − 1 ) W ( Λ k ∣ W k , ν 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})} <p>given the definitions </p> β k = β 0 + N k m k = 1 β k ( β 0 μ 0 + N k x ¯ k ) W k − 1 = W 0 − 1 + N k S k + β 0 N k β 0 + N k ( x ¯ k − μ 0 ) ( x ¯ k − μ 0 ) T ν k = ν 0 + N k N k = ∑ n = 1 N r n k x ¯ k = 1 N k ∑ n = 1 N r n k x n S k = 1 N k ∑ n = 1 N r n k ( x n − x ¯ k ) ( x n − x ¯ k ) T {\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}}} <p>Finally, notice that these functions require the values of r n k {\displaystyle r_{nk}} , which make use of ρ n k {\displaystyle \rho _{nk}} , which is defined in turn based on E [ ln π k ] {\displaystyle \operatorname {E} [\ln \pi _{k}]} , E [ ln | Λ k | ] {\displaystyle \operatorname {E} [\ln |\mathbf {\Lambda } _{k}|]} , and E μ k , Λ k [ ( x n − μ k ) T Λ k ( x n − μ 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})]} . Now that we have determined the distributions over which these expectations are taken, we can derive formulas for them: </p> E μ k , Λ k [ ( x n − μ k ) T Λ k ( x n − μ k ) ] = D β k − 1 + ν k ( x n − m k ) T W k ( x n − m k ) ln Λ ~ k ≡ E [ ln | Λ k | ] = ∑ i = 1 D ψ ( ν k + 1 − i 2 ) + D ln 2 + ln | W k | ln π ~ k ≡ E [ ln | π k | ] = ψ ( α k ) − ψ ( ∑ i = 1 K α i ) {\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}}} <p>These results lead to </p> r n k ∝ π ~ k Λ ~ k 1 / 2 exp { − D 2 β k − ν k 2 ( x n − m k ) T W k ( x n − m k ) } {\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\}} <p>These can be converted from proportional to absolute values by normalizing over k {\displaystyle k} so that the corresponding values sum to 1. </p><p>Note that: </p> <ol><li>The update equations for the parameters β k {\displaystyle \beta _{k}} , m k {\displaystyle \mathbf {m} _{k}} , W k {\displaystyle \mathbf {W} _{k}} and ν k {\displaystyle \nu _{k}} of the variables μ k {\displaystyle \mathbf {\mu } _{k}} and Λ k {\displaystyle \mathbf {\Lambda } _{k}} depend on the statistics N k {\displaystyle N_{k}} , x ¯ k {\displaystyle {\bar {\mathbf {x} }}_{k}} , and S k {\displaystyle \mathbf {S} _{k}} , and these statistics in turn depend on r n k {\displaystyle r_{nk}} .</li> <li>The update equations for the parameters α 1 … K {\displaystyle \alpha _{1\dots K}} of the variable π {\displaystyle \mathbf {\pi } } depend on the statistic N k {\displaystyle N_{k}} , which depends in turn on r n k {\displaystyle r_{nk}} .</li> <li>The update equation for r n k {\displaystyle r_{nk}} has a direct circular dependence on β k {\displaystyle \beta _{k}} , m k {\displaystyle \mathbf {m} _{k}} , W k {\displaystyle \mathbf {W} _{k}} and ν k {\displaystyle \nu _{k}} as well as an indirect circular dependence on W k {\displaystyle \mathbf {W} _{k}} , ν k {\displaystyle \nu _{k}} and α 1 … K {\displaystyle \alpha _{1\dots K}} through π ~ k {\displaystyle {\widetilde {\pi }}_{k}} and Λ ~ k {\displaystyle {\widetilde {\Lambda }}_{k}} .</li></ol> <p>This suggests an iterative procedure that alternates between two steps: </p> <ol><li>An E-step that computes the value of r n k {\displaystyle r_{nk}} using the current values of all the other parameters.</li> <li>An M-step that uses the new value of r n k {\displaystyle r_{nk}} to compute new values of all the other parameters.</li></ol> <p>Note that these steps correspond closely with the standard EM algorithm to derive a <a href="/facts/Maximum_likelihood/0Yq2dpQD">maximum likelihood</a> or <a href="/facts/Maximum_a_posteriori/Y5zHyqfX">maximum a posteriori</a> (MAP) solution for the parameters of a <a href="/facts/Gaussian_mixture_model/WWe6TJlg">Gaussian mixture model</a>. The responsibilities r n k {\displaystyle r_{nk}} in the E step correspond closely to the <a href="/facts/Posterior_probability/fojsMD0D">posterior probabilities</a> of the latent variables given the data, i.e. p ( Z ∣ X ) {\displaystyle p(\mathbf {Z} \mid \mathbf {X} )} ; the computation of the statistics N k {\displaystyle N_{k}} , x ¯ k {\displaystyle {\bar {\mathbf {x} }}_{k}} , and S k {\displaystyle \mathbf {S} _{k}} corresponds closely to the computation of corresponding "soft-count" statistics over the data; and the use of those statistics to compute new values of the parameters corresponds closely to the use of soft counts to compute new parameter values in normal EM over a Gaussian mixture model. </p> <h2 id="exponential-family-distributions">Exponential-family distributions</h2> <p>Note that in the previous example, once the distribution over unobserved variables was assumed to factorize into distributions over the "parameters" and distributions over the "latent data", the derived "best" distribution for each variable was in the same family as the corresponding prior distribution over the variable. This is a general result that holds true for all prior distributions derived from the <a href="/facts/Exponential_family/1LkkqEIf">exponential family</a>. </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Variational_message_passing/vptPNdXb">Variational message passing</a>: a modular algorithm for variational Bayesian inference.</li> <li><a href="/facts/Variational_autoencoder/I5wBbYpR">Variational autoencoder</a>: an artificial neural network belonging to the families of probabilistic graphical models and Variational Bayesian methods.</li> <li><a href="/facts/Expectation%25E2%2580%2593maximization_algorithm/32P2bdBc">Expectation–maximization algorithm</a>: a related approach which corresponds to a special case of variational Bayesian inference.</li> <li><a href="/facts/Generalized_filtering/ld6tXByk">Generalized filtering</a>: a variational filtering scheme for nonlinear state space models.</li> <li><a href="/facts/Calculus_of_variations/Ya5Lsma3">Calculus of variations</a>: the field of mathematical analysis that deals with maximizing or minimizing functionals.</li> <li>Maximum entropy discrimination: This is a variational inference framework that allows for introducing and accounting for additional large-margin constraints<a class="footnote-ref" id="fnref:13" href="#fn:13"><sup>13</sup></a></li></ul> <h2 id="external-links">External links</h2> <ul><li><a href="https://www.inference.phy.cam.ac.uk/mackay/itila/">The on-line textbook: Information Theory, Inference, and Learning Algorithms</a>, by <a href="/facts/David_J.C._MacKay/NK3D8JME">David J.C. MacKay</a> provides an introduction to variational methods (p. 422).</li> <li><a href="https://www.robots.ox.ac.uk/~sjrob/Pubs/fox_vbtut.pdf">A Tutorial on Variational Bayes</a>. Fox, C. and Roberts, S. 2012. Artificial Intelligence Review, <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1007%2Fs10462-011-9236-8">10.1007/s10462-011-9236-8</a>.</li> <li><a href="https://www.gatsby.ucl.ac.uk/vbayes/">Variational-Bayes Repository</a> A repository of research papers, software, and links related to the use of variational methods for approximate Bayesian learning up to 2003.</li> <li><a href="https://www.cse.buffalo.edu/faculty/mbeal/thesis/index.html">Variational Algorithms for Approximate Bayesian Inference</a>, by M. J. Beal includes comparisons of EM to Variational Bayesian EM and derivations of several models including Variational Bayesian HMMs.</li> <li><a href="https://www.cs.jhu.edu/~jason/tutorials/variational.html">High-Level Explanation of Variational Inference</a> by Jason Eisner may be worth reading before a more mathematically detailed treatment.</li> <li><a href="https://arxiv.org/abs/1803.10998">Copula Variational Bayes inference via information geometry (pdf)</a> by Tran, V.H. 2018. This paper is primarily written for students. Via <a href="/facts/Bregman_divergence/wREwF8qo">Bregman divergence</a>, the paper shows that Variational Bayes is simply a generalized Pythagorean projection of true model onto an arbitrarily correlated (copula) distributional space, of which the independent space is merely a special case.</li> <li><a href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076">An in depth introduction to Variational Bayes note</a>. Nguyen, D. 2023</li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Tran, Viet Hung (2018). "Copula Variational Bayes inference via information geometry". arXiv:1803.10998 [cs.IT]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> <li id="fn:2"><p>Adamčík, Martin (2014). "The Information Geometry of Bregman Divergences and Some Applications in Multi-Expert Reasoning". Entropy. 16 (12): 6338–6381. Bibcode:2014Entrp..16.6338A. doi:10.3390/e16126338. <a href="https://doi.org/10.3390%2Fe16126338" target="_blank">https://doi.org/10.3390%2Fe16126338</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li> <li id="fn:3"><p>Tran, Viet Hung (2018). "Copula Variational Bayes inference via information geometry". arXiv:1803.10998 [cs.IT]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li> <li id="fn:4"><p>Tran, Viet Hung (2018). "Copula Variational Bayes inference via information geometry". arXiv:1803.10998 [cs.IT]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li> <li id="fn:5"><p>Tran, Viet Hung (2018). "Copula Variational Bayes inference via information geometry". arXiv:1803.10998 [cs.IT]. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:5" class="footnote-back-ref">↩</a></p></li> <li id="fn:6"><p>Nguyen, Duy (15 August 2023). "AN IN DEPTH INTRODUCTION TO VARIATIONAL BAYES NOTE". doi:10.2139/ssrn.4541076. SSRN 4541076. Retrieved 15 August 2023. <a href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076" target="_blank">https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076</a> <a href="#fnref:6" class="footnote-back-ref">↩</a></p></li> <li id="fn:7"><p>Lee, Se Yoon (2021). "Gibbs sampler and coordinate ascent variational inference: A set-theoretical review". Communications in Statistics - Theory and Methods. 51 (6): 1–21. arXiv:2008.01006. doi:10.1080/03610926.2021.1921214. S2CID 220935477. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:7" class="footnote-back-ref">↩</a></p></li> <li id="fn:8"><p>Boyd, Stephen P.; Vandenberghe, Lieven (2004). Convex Optimization (PDF). Cambridge University Press. ISBN 978-0-521-83378-3. Retrieved October 15, 2011. <a href="978-0-521-83378-3" target="_blank">978-0-521-83378-3</a> <a href="#fnref:8" class="footnote-back-ref">↩</a></p></li> <li id="fn:9"><p>Lee, Se Yoon (2021). "Gibbs sampler and coordinate ascent variational inference: A set-theoretical review". Communications in Statistics - Theory and Methods. 51 (6): 1–21. arXiv:2008.01006. doi:10.1080/03610926.2021.1921214. S2CID 220935477. <a href="/wiki/ArXiv_(identifier)" target="_blank">/wiki/ArXiv_(identifier)</a> <a href="#fnref:9" class="footnote-back-ref">↩</a></p></li> <li id="fn:10"><p>Bishop, Christopher M. (2006). "Chapter 10". Pattern Recognition and Machine Learning. Springer. ISBN 978-0-387-31073-2. <a href="978-0-387-31073-2" target="_blank">978-0-387-31073-2</a> <a href="#fnref:10" class="footnote-back-ref">↩</a></p></li> <li id="fn:11"><p>Nguyen, Duy (15 August 2023). "AN IN DEPTH INTRODUCTION TO VARIATIONAL BAYES NOTE". doi:10.2139/ssrn.4541076. SSRN 4541076. Retrieved 15 August 2023. <a href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076" target="_blank">https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076</a> <a href="#fnref:11" class="footnote-back-ref">↩</a></p></li> <li id="fn:12"><p>Nguyen, Duy (15 August 2023). "AN IN DEPTH INTRODUCTION TO VARIATIONAL BAYES NOTE". doi:10.2139/ssrn.4541076. SSRN 4541076. Retrieved 15 August 2023. <a href="https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076" target="_blank">https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4541076</a> <a href="#fnref:12" class="footnote-back-ref">↩</a></p></li> <li id="fn:13"><p>Sotirios P. Chatzis, “Infinite Markov-Switching Maximum Entropy Discrimination Machines,” Proc. 30th International Conference on Machine Learning (ICML). Journal of Machine Learning Research: Workshop and Conference Proceedings, vol. 28, no. 3, pp. 729–737, June 2013. <a href="http://proceedings.mlr.press/v28/chatzis13.pdf" target="_blank">http://proceedings.mlr.press/v28/chatzis13.pdf</a> <a href="#fnref:13" class="footnote-back-ref">↩</a></p></li> </ol>