變分貝葉斯方法

變分貝葉斯方法是近似貝葉斯推理與機器學習中較難積分的一系列技術，通常用於由觀測變量及未知參數、潛變量組成的複雜統計模型，之間存在的各種關係可用圖模式描述。貝葉斯推理中，參數和潛變量一般歸為「未觀察變量」。變分貝葉斯方法主要用於兩個目的：

- 提供未觀測變量後驗概率的分析近似，以便對變量統計推斷
- 得出觀測數據的邊緣似然的下界（即給定模型的數據的邊緣機率，邊緣化未觀測變量）。通常用於模型選擇，一般認為給定模型的邊緣似然越高，擬合效果越好，產生數據的可能性越大（另見貝葉斯因子）。

就目的1（近似後驗概率），變分貝葉斯是蒙特卡洛法（特別是馬爾可夫鏈蒙特卡洛，如吉布斯採樣）的替代，可採用完全貝葉斯的方法，對難以取值或抽樣的複雜分布進行統計推斷。蒙特卡洛技術利用一組樣本對精確後驗值進行數值近似，而變分貝葉斯可求近似後驗的局部最優的精確解析解。

變分貝葉斯可視作最大期望算法（EM）的推廣，從估計參數的單一最似然值的最大後驗概率（MAP估計），推廣到計算參數及潛變量的整個近似後驗分布的完全貝葉斯估計。貝葉斯估計同EM一樣，也能找到一組最優參數值，且與EM有相同的交替結構——基於不能求解析解的相依方程組。

很多應用中，變分貝葉斯能更快得到與吉布斯採樣精度相當的解，但推導迭代方程組需要更多工作。即便是很多概念上非常簡單的模型也如此，下面以只有2個參數、沒有潛變量的基本不分層模型為例說明。

## 數學推導

### 問題

變分推斷中，數據$\mathbf{X}$的一組未觀測變量$\mathbf{Z} = \{Z_1 \dots Z_n\}$的後驗分布近似於所謂變分分布$Q(\mathbf{Z})$

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

分布$Q(\mathbf{Z})$屬於在形式上比$P(\mathbf{Z}\mid \mathbf{X})$簡單的分布族（如高斯分布族），這是為使$Q(\mathbf{Z})$與真實後驗$P(\mathbf{Z}\mid \mathbf{X})$更相似。

相似性用不相似函數$d(Q; P)$衡量，推斷（inference）通過求分布$Q(\mathbf{Z})$使$d(Q; P)$最小化進行。

### KL散度

最常見的變分貝葉斯法以Q與P的KL散度作為不相似函數，很適合最小化。KL散度的定義是

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

注意Q和P是相反的。反向KL散度在概念上類似於期望最大化演算法（KL散度則產生期望傳播演算法）。

### 難解性

變分技術通常用於近似以下參數：

$$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'}$$

對$\mathbf Z$邊緣化（以計算分母中的$P(\mathbf X)$）通常是難解的，因為$\mathbf Z$的搜索空間在組合上可能很大。因此，我們用$Q(\mathbf Z) \approx P(\mathbf Z \mid \mathbf X)$求近似。

### 證據下界

考慮到$P(\mathbf Z \mid \mathbf X) = \frac{P(\mathbf X, \mathbf Z)}{P(\mathbf X)}$，KL散度也可寫作

$$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]$$

因為$Q(\mathbf{Z})$是分布，所以$\sum_\mathbf{Z} Q(\mathbf{Z}) = 1$；又由於$P(\mathbf{X})$對於$\mathbf Z$是常數，所以

$$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})$$

根據（離散隨機變量）期望的定義，上式可以寫作

$$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})$$

重排為

$$\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)$$

由於對數證據$\log P(\mathbf{X})$對於Q是定值，最大化末項$\mathcal{L}(Q)$會使Q對P的KL散度最小化。適當選擇Q可以讓$\mathcal{L}(Q)$更容易計算，更容易最大化，因此既有後驗$P(\mathbf{Z}\mid \mathbf{X})$的解析近似Q，也有對數證據$\log P(\mathbf{X})$的下界$\mathcal{L}(Q)$（KL散度非負）。

下界$\mathcal{L}(Q)$與熱力學自由能類似，也稱作（負）變分自由能，因為它也可以表為負能量$\operatorname{E}_{Q}[\log P(\mathbf{Z},\mathbf{X})]$加Q的熵。$\mathcal{L}(Q)$項也稱作證據下界（ELBO），是數據的對數證據的下界。

### 證明

根據布雷格曼散度（KL散度是其特例）的廣義勾股定理，可以證明

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

其中$\mathcal{C}$是凸集，且

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

時等號成立。這時，全局最小值$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)$，其中$\mathbf{Z}=\{\mathbf{Z_1},\mathbf{Z_2}\}$可以如下求得：

$$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),$$

其中歸一化常數為

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

$\zeta(\mathbf{X})$項在實踐中常稱作證據下界（ELBO），因為$P(\mathbf{X})\geq\zeta(\mathbf{X})=\exp(\mathcal{L}(Q^{*}))$，如上所示。

交換$\mathbf{Z}_1,\ \mathbf{Z}_2$的角色，可以迭代計算真實模型邊際$P(\mathbf{Z}_1\mid\mathbf{X}),\ P(\mathbf{Z}_2\mid\mathbf{X})$的近似$q^{*}(\mathbf{Z}_1),\ q^{*}(\mathbf{Z}_2)$。這種迭代保證單調收斂，但收斂到的$Q^{*}$只是$D_{\mathrm{KL}}(Q\parallel P)$的局部極小值。

若約束空間$\mathcal{C}$被限制在獨立空間內，即$q^{*}(\mathbf{Z}_1\mid\mathbf{Z}_2) = q^{*}(\mathbf{Z_1})$，則上述迭代將成為所謂平均場近似（mean field approximation）$Q^{*}(\mathbf{Z}) = q^{*}(\mathbf{Z}_1)q^{*}(\mathbf{Z}_2),$。

## 平均場近似

變分分布$Q(\mathbf{Z})$常被假定為潛變量的某劃分的因式分解，即對潛變量$\mathbf{Z}$的劃分$\mathbf{Z}_1 \dots \mathbf{Z}_M$。

$$Q(\mathbf{Z}) = \prod_{i=1}^M q_i(\mathbf{Z}_i\mid \mathbf{X})$$

由變分法可以證明，對每個因子$q_j$的最佳分布$q_j^{*}$（就KL散度最小的分布而言，如上所述）滿足

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

其中$\operatorname{E}_{q^*_{-j}} [\ln p(\mathbf{Z}, \mathbf{X})]$是數據與潛變量的對數聯合機率的期望值，是相對於不屬於劃分的所有變量的$q^*$而取。關於$q_j^{*}(\mathbf{Z}_j\mid \mathbf{X})$的推導見的引理 4.1。

實踐中，通常用對數表示：

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

上式中的常數與歸一化常數（$q_j^{*}$上述表達式中的分母）有關，通常通過檢查來恢復，因為表達式的其餘部分通常是已知的分布（如正態分布、伽馬分布等）。

利用期望的性質，式$\operatorname{E}_{q^*_{-j}} [\ln p(\mathbf{Z}, \mathbf{X})]$通常可簡為潛變量先驗分布的固定超參數及不屬於當前劃分的潛變量（即不屬於$\mathbf{Z}_j$的潛變量）的期望值的函數（有時還包括方差等高階矩）。這就在某一划分的變量分布參數同其他劃分的變量的期望值之間產生了循環依賴，自然就需要類似期望最大化演算法的迭代法，以某種方式（也許隨機）初始化潛變量期望（可能還有高階矩），再利用當前期望依次計算每個分布的參數、適當設置新算出的分布的期望。這種算法保證收斂。

換一種說法，對每個變量劃分，簡化其中變量分布的表達式、研究分布與相關變量的函數依賴關係，通常就能確定分布的族（反過來確定了常數）。分布的參數公式可用先驗分布的超參數（已知常數）表示，也可用其他劃分中變量函數的期望表示。通常這些期望可簡為變量本身的期望值函數（即平均值）；有時也會出現變量平方（與變量方差有關）或更高次冪高階矩）的期望。其他變量的分布一般來自已知族，相關期望的公式可查到；但公式取決於分布的參數，參數又取決於其他變量的期望，所以每個變量的分布參數公式都可表為一群變量間非線性相互依賴的方程。通常不能直接求解，不過可以用簡單的迭代演算法，大多數時候都能保證收斂。

## 變分推理對偶公式

下面的定理稱作變分推理對偶公式，解釋了變分貝葉斯方法中變分分布的一些重要性質。

3 考慮兩機率空間$(\Theta,\mathcal{F},P),\ (\Theta,\mathcal{F},Q)$，其中$Q \ll P$。假設存在共同的主導機率測度$\lambda$使$P \ll \lambda,\ Q \ll \lambda$。令h是$(\Theta,\mathcal{F},P)$上的任意實值隨機變量，滿足$h \in L_1(P)$，則下列等式成立

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

此外，若且唯若（supremum）

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

關於機率測度Q幾乎是確定的，其中$p(\theta) = dP/d\lambda,\ q(\theta) = dQ/d\lambda$分別表示機率測度P、Q對$\lambda$的Radon–Nikodym導數，右式才取上界。

## 基本例子

考慮簡單的不分層貝葉斯模型，由來自正態分布（均值、方差未知）的獨立同分布觀測量組成。下面我們將詳細介紹這模型，以說明變分貝葉斯方法的工作原理。

為了數學上的方便，下例中我們用精度（即方差倒數；多元正態分布時是協方差矩陣的逆。理論上精度和方差等價，因為兩者間是雙射）。

### 數學模型

將共軛先驗分布置於未知均值$\mu$、精度$\tau$，即均值也遵循高斯分布，而精度遵循伽馬分布。

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

先驗分布超參數$\mu_0, \lambda_0, a_0$、$b_0$是已知定值，可設為小正數，以得到較寬的先驗分布，表明我們對$\mu$、$\tau$的先驗分布一無所知。

已知N個數據點$\mathbf{X} = \{x_1, \ldots, x_N\}$，而我們的目標是推斷參數$\mu,\ \tau$的後驗分布$q(\mu, \tau)=p(\mu,\tau\mid x_1, \ldots, x_N)$。

### 聯合機率

所有變量的聯合機率可重寫為

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

個體因子是

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

其中

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

### 因式分解近似

假設$q(\mu,\tau) = q(\mu)q(\tau)$，即後驗分布因式分解為$\mu,\ \tau$的獨立因子。這種假設是變分貝葉斯方法的基礎。事實上，真正的後驗分布並非如此（這種簡單情形我們知道它是正態-伽馬分布），因此所得結果將是近似值。

### 導出q(μ)

則

$$\begin{align} \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 \\ <pre> &= \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 </pre> \end{align}$$

上面的推導中，$C,\ C_2,\ C_3$指對$\mu$為常數的值。注意$\operatorname{E}_{\tau}[\ln p(\tau)]$項不是$\mu$的函數，無論$\mu$的值是多少，它都不變。因此，第三行中可以把它吸收到末尾的常數項，第七行也如此。

最後一行是$\mu$的二次多項式。由於這是$q_\mu^*(\mu)$的對數，我們可以看到$q_\mu^*(\mu)$本身是正態分布。

通過一些繁瑣的計算（展開括號里的平方、分離出涉及$\mu$、$\mu^2$的項、對$\mu$應用配方法），可得到正態分布的參數：

$$\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 \\ <pre> &= -\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 </pre> \end{align}$$

注意上面所有步驟用二次和公式都能化簡。

即

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

### 推導q(τ)

$q_\tau^*(\tau)$的推導大致相同，簡潔起見略去細節。

$$\begin{align} \ln q_\tau^*(\tau) &= \operatorname{E}_{\mu}[\ln p(\mathbf{X}\mid \mu,\tau) + \ln p(\mu\mid \tau)] + \ln p(\tau) + \text{constant} \\ <pre> &= (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} </pre> \end{align}$$

兩側取指數，可見$q_\tau^*(\tau)$服從伽馬分布。

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

### 計算參數

回憶前幾節的結論：

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

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

每種情形下，某變量的分布參數取決於對另一變量的期望。可用正態分布和伽馬分布矩期望的標準公式推廣期望：

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

大多數時候將這些公式應用到上述方程不困難，但$b_N$需要更多工作：

$$\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)^2\right] \\ <pre> &= 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] \\ </pre> \end{align}$$

這樣就可以寫出參數方程如下，不需要期望：

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

注意$\lambda_N$和$b_N$存在相互依賴關係，自然產生了類似最大期望算法的演算法：

- 計算$\sum_{n=1}^N x_n,\ \sum_{n=1}^N x_n^2.$用所得值計算$\mu_N,\ a_N.$
- 初始化$\lambda_N$為任意值。
- 用$\lambda_N$的現值和其它參數的已知值計算$b_N$。
- 用$b_N$的現值和其它參數的已知值計算$\lambda_N$。
- 重複最後兩部，直到收斂（即直到兩值的更新量不超過閾值）。

然後就有了後驗參數近似分布的超參數值，可用於計算後驗的任何屬性，如均值、方差、95%最高密度區域（包含總機率95%的最小區間）等等。

可以證明這種算法保證收斂到局部極大值。

還要注意，後驗分布與相應的先驗分布有同樣的形式。只需假設分布可以分解，分布形式便隨之而來。事實證明，後驗、先驗形式相同並非巧合，而是先驗分布為指數族成員時的一般結果，大多數標準分布都如此。

## 進一步討論

### 分步法

上面例子展示了給定貝葉斯網路中推導後驗概率密度的變分貝葉斯近似的方法：

- 用圖模式描述網路，確定觀察變量（數據）$\mathbf{X}$和未觀察變量（參數$\boldsymbol\Theta$、潛變量$\mathbf{Z}$）及其條件概率分布。之後，變分貝葉斯構建後驗概率$p(\mathbf{Z},\boldsymbol\Theta\mid\mathbf{X})$，這近似具有可分解分布的基本特性，即是多個獨立分布在不相交的未觀測變量子集上的積。
- 將未觀測變量劃分為多個子集，在其上推導出獨立因子。這種方法沒有通用的程序，子集太多會使近似結果不佳，子集過少會使整個變分貝葉斯方法變得難以實現。第一種分割方法是將參數和潛變量分開，一般足以產生可操作的結果。劃分記作$\mathbf{Z}_1,\ldots,\mathbf{Z}_M$。
- 對給定的劃分$\mathbf{Z}_j$，用基本方程$\ln q_j^{*}(\mathbf{Z}_j\mid \mathbf{X}) = \operatorname{E}_{i \neq j} [\ln p(\mathbf{Z}, \mathbf{X})] + \text{constant}$寫出最佳近似分布$q_j^{*}(\mathbf{Z}_j\mid \mathbf{X})$。
- 用圖模式填寫聯合機率分布公式。任何不涉及$\mathbf{Z}_j$中變量的分量條件分布都可忽略，它們將被摺疊到常數項中。
- 按上面例子，簡化公式並應用期望算子。理想情況下這應該簡化為不屬於$\mathbf{Z}_j$的變量的基本函數的期望（如第一或第二原始矩、對數期望等等）。為讓變分貝葉斯方法順利運行，這些期望值通常應可以解析地表為變量分布的參數和/或超參數的函數。這些期望項對當前劃分的變量都是常數。
- 式對當前劃分中變量的函數形式表明了分布的類型。特別地，取公式指數，便得到分布的機率密度函數（PDF）（或至少是與之成正比的函數，帶未知歸一化常數）。要使方法可操作，應能識別出函數形式所屬的已知分布；要將公式轉為匹配已知分布的PDF的形式，可能需要大量計算。若能做到，就可根據定義恢復歸一化常數，並提取公式中的適當部分得到已知分布的參數方程。
- 期望都可用非當前劃分變量的函數進行解析替換、並將PDF轉為可與已知分布相對應的形式，此時結果應是一組方程，將最優參數值表為其他劃分變量參數的函數。
- 若這程序適用於所有劃分，就會產生相互依賴的方程組，指明所有參數的最優值。
- 最大期望算法型程序可為每個參數選取初值，再迭代。每一步中，我們都會在方程中循環，依次更新每個參數。可以確證收斂。

## 另見

- 變分消息傳遞：半分貝葉斯推理的模式化演算法
- 變分自編碼器：屬於機率圖模式和變分貝葉斯法的ANN
- 最大期望算法：與變分貝葉斯推理特例相對應的一種相關方法
- 廣義濾波：非線性狀態空間模型的變分濾波框架
- 變分法：數學分析中處理最值化函數的領域

## 外部連結

- The on-line textbook: Information Theory, Inference, and Learning Algorithms , by David J.C. MacKay provides an introduction to variational methods (p. 422).
- A Tutorial on Variational Bayes . Fox, C. and Roberts, S. 2012. Artificial Intelligence Review, 10.1007/s10462-011-9236-8.
- Variational-Bayes Repository A repository of research papers, software, and links related to the use of variational methods for approximate Bayesian learning up to 2003.
- Variational Algorithms for Approximate Bayesian Inference , by M. J. Beal includes comparisons of EM to Variational Bayesian EM and derivations of several models including Variational Bayesian HMMs.
- High-Level Explanation of Variational Inference by Jason Eisner may be worth reading before a more mathematically detailed treatment.
- Copula Variational Bayes inference via information geometry (pdf) by Tran, V.H. 2018. This paper is primarily written for students. Via Bregman divergence, 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.
- An in depth introduction to Variational Bayes note . Nguyen, D. 2023