Markov Chain Monte Carlo (MCMC) 介紹 Markov Chain Monte Carlo (MCMC)，一種適用於複雜分佈上採樣的模型。
2024/09/15 介紹 若要對分佈採樣或計算期望值問題，可以透過 accept-reject method 等 Monte Carlo 方法；但對於複雜、多變量、變數不獨立的分佈，可以透過 Markov Chain Monte Carlo (MCMC) 方法。
假設一隨機變量 x ∈ X x \in X x ∈ X p ( x ) p (x) p ( x ) p ( x ) p (x) p ( x ) X = { X 0 , X 1 , ⋯ , X t , ⋯ } X = \{ X_0, X_1, \cdots, X_t, \cdots \} X = { X 0 , X 1 , ⋯ , X t , ⋯ } S S S X X X p p p
因此核心問題就是如何建構此 Markov chain，一種方案是找出 transition function 並建構 reversible Markov chain 使其滿足 ergodic theorem，因此起始點就會與結果無關，最終都會收斂至唯一的穩定分佈，一般常用的方案有 Metropolis–Hastings algorithm 或 Gibbs sampling。
MCMC 在相鄰兩次的採樣是相關的，若需要獨立樣本可以考慮每隔數個週期後再接受結果，並且相對於 accept-reject method，並不需要找到已知且容易採樣的分佈進行拒絕，因為在高維度不容易尋找，一般而言 MCMC 方法的效率會高於 accept-reject method。
步驟 建構滿足 ergodic theorem 的 Markov chain，並使其穩定分佈與採樣目標 p ( x ) p (x) p ( x ) 選定狀態空間中任意起始點 x 0 x_0 x 0 { x 0 , x 1 , ⋯ , x t , ⋯ } \{ x_0, x_1, \cdots, x_t, \cdots \} { x 0 , x 1 , ⋯ , x t , ⋯ } 選定兩整數 m < n m < n m < n { x m + 1 , x m + 2 , ⋯ , x n } \{ x_{m + 1}, x_{m + 2}, \cdots, x_n \} { x m + 1 , x m + 2 , ⋯ , x n } 然而有幾個點需要注意
如何定義 Markov chain 使其滿足 ergodic theorem 如何選取收斂步數 m m m 如何選取迭代步數 n n n MCMC 與機器學習 MCMC 能用於機器學習中的模型學習與推理，使其在機器學習中扮演重要角色，尤其在 Bayesian learning；考慮觀測到的資料為隨機變量 y ∈ Y y \in Y y ∈ Y x ∈ X x \in X x ∈ X
p ( x ∣ y ) = p ( x ) p ( y ∣ x ) ∫ X P ( y ∣ x ′ ) p ( x ′ ) d x ′ \begin{align*}
p (x | y)
= \frac{p (x) p (y | x)}{\int_X P (y | x') p (x') dx'}
\end{align*} p ( x ∣ y ) = ∫ X P ( y ∣ x ′ ) p ( x ′ ) d x ′ p ( x ) p ( y ∣ x ) 在其中有 3 種能應用 MCMC 的場景
Normalization:
∫ X p ( y ∣ x ′ ) p ( x ′ ) d x ′ \begin{align*}
\int_X p (y | x') p (x') dx'
\end{align*} ∫ X p ( y ∣ x ′ ) p ( x ′ ) d x ′ Marginalization: 若存在隱藏變數 z ∈ Z z \in Z z ∈ Z
p ( x ∣ y ) = ∫ Z p ( x , z ∣ y ) d z \begin{align*}
p (x | y)
= \int_Z p (x, z | y) dz
\end{align*} p ( x ∣ y ) = ∫ Z p ( x , z ∣ y ) d z Expectation: 若需計算 f ( x ) f (x) f ( x )
E p [ f ( X ) ] = ∫ X f ( x ) p ( x ∣ y ) d x \begin{align*}
E_p [f (X)]
= \int_X f (x) p (x | y) dx
\end{align*} E p [ f ( X )] = ∫ X f ( x ) p ( x ∣ y ) d x 上述任一情況如果模型複雜，則計算將會非常困難，但能用 MCMC 做一個通用且有效率的選項
Metropolis–Hasting Algorithm 這邊介紹 MCMC 方法中常用的一個方案 Metropolis–Hasting algorithm。給定 pdf p ( x ) p (x) p ( x ) p ( x , x ′ ) p (x, x') p ( x , x ′ )
p ( x , x ′ ) = q ( x , x ′ ) α ( x , x ′ ) \begin{align*}
p (x, x')
= q (x, x') \alpha (x, x')
\end{align*} p ( x , x ′ ) = q ( x , x ′ ) α ( x , x ′ ) 其中 q ( x , x ′ ) q (x, x') q ( x , x ′ ) α ( x , x ′ ) \alpha (x, x') α ( x , x ′ ) q ( x , x ′ ) q (x, x') q ( x , x ′ ) α ( x , x ′ ) \alpha (x, x') α ( x , x ′ )
α ( x , x ′ ) = min { 1 , p ( x ′ ) q ( x ′ , x ) p ( x ) q ( x , x ′ ) } \begin{align*}
\alpha (x, x')
= \min \left\lbrace 1, \frac{p (x') q (x', x)}{p (x) q (x, x')} \right\rbrace
\end{align*} α ( x , x ′ ) = min { 1 , p ( x ) q ( x , x ′ ) p ( x ′ ) q ( x ′ , x ) } 使得
p ( x , x ′ ) = { q ( x , x ′ ) if p ( x ′ ) q ( x ′ , x ) ≥ p ( x ) q ( x , x ′ ) q ( x , x ′ ) p ( x ′ ) p ( x ) if p ( x ′ ) q ( x ′ , x ) > p ( x ) q ( x , x ′ ) \begin{align*}
p (x, x')
= \begin{cases}
q (x, x') & \text{if } p (x') q(x', x) \geq p (x) q (x, x') \\
q (x, x') \frac{p (x')}{p (x)} & \text{if } p (x') q(x', x) > p (x) q (x, x')
\end{cases}
\end{align*} p ( x , x ′ ) = { q ( x , x ′ ) q ( x , x ′ ) p ( x ) p ( x ′ ) if p ( x ′ ) q ( x ′ , x ) ≥ p ( x ) q ( x , x ′ ) if p ( x ′ ) q ( x ′ , x ) > p ( x ) q ( x , x ′ ) 而在 Markov chain p ( x , x ′ ) p (x, x') p ( x , x ′ ) t − 1 t - 1 t − 1 x x x x t − 1 = x x_{t - 1} = x x t − 1 = x q ( x , x ′ ) q (x, x') q ( x , x ′ ) x ′ x' x ′ α ( x , x ′ ) \alpha (x, x') α ( x , x ′ ) x t = x ′ x_t = x' x t = x ′ x t = x x_t = x x t = x
透過此方式建構的 p ( x , x ′ ) p (x, x') p ( x , x ′ ) p ( x , x ′ ) p (x, x') p ( x , x ′ ) p ( x ) p (x) p ( x )
演算法 在 Markov chain p ( x , x ′ ) p (x, x') p ( x , x ′ )
選定起始點 x 0 x_0 x 0 對 i = 1 , 2 , ⋯ , n i = \conti{n} i = 1 , 2 , ⋯ , n 設定狀態 x i − 1 = x x_{i - 1} = x x i − 1 = x q ( x , x ′ ) q (x, x') q ( x , x ′ ) x ′ x' x ′ 計算接受機率 α ( x , x ′ ) \alpha (x, x') α ( x , x ′ ) 採樣 u ∼ Unif ( 0 , 1 ) u \sim \unif (0, 1) u ∼ Unif ( 0 , 1 ) u ≤ α ( x , x ′ ) u \leq \alpha (x, x') u ≤ α ( x , x ′ ) x i = x ′ x_i = x' x i = x ′ x i = x x_i = x x i = x 獲取採樣集合 { x m + 1 , x m + 2 , ⋯ , x n } \{ x_{m + 1}, x_{m + 2}, \cdots, x_n \} { x m + 1 , x m + 2 , ⋯ , x n } Proposal Distribution 選取 Proposal Distribution q ( x , x ′ ) q (x, x') q ( x , x ′ )
對稱 (symmetic) 的，即 q ( x , x ′ ) = q ( x ′ , x ) q (x, x') = q (x', x) q ( x , x ′ ) = q ( x ′ , x ) 獨立 (independent) 的，即 q ( x , x ′ ) = q ( x ′ ) q (x, x') = q (x') q ( x , x ′ ) = q ( x ′ ) 若選取對稱分佈，即 q ( x , x ′ ) = q ( x ′ , x ) q (x, x') = q (x', x) q ( x , x ′ ) = q ( x ′ , x )
α ( x , x ′ ) = { 1 , p ( x ′ ) p ( x ) } \begin{align*}
\alpha (x, x')
= \left\lbrace 1, \frac{p (x')}{p (x)} \right\rbrace
\end{align*} α ( x , x ′ ) = { 1 , p ( x ) p ( x ′ ) } 此選取稱為 Metropolis selection；在這中有個特例，設定 q ( x , x ′ ) = q ( ∣ x − x ′ ∣ ) q (x, x') = q ( |x - x'|) q ( x , x ′ ) = q ( ∣ x − x ′ ∣ ) x x x x ′ x' x ′
q ( x , x ′ ) ∝ exp [ − ( x − x ′ ) 2 2 ] \begin{align*}
q (x, x')
\propto \exp \left[ - \frac{(x - x')^2}{2} \right]
\end{align*} q ( x , x ′ ) ∝ exp [ − 2 ( x − x ′ ) 2 ] 若選擇獨立分佈，即 q ( x , x ′ ) = q ( x ′ ) q (x, x') = q (x') q ( x , x ′ ) = q ( x ′ )
α ( x , x ′ ) = { 1 , p ( x ′ ) / q ( x ′ ) p ( x ) / q ( x ) } \begin{align*}
\alpha (x, x')
= \left\lbrace 1, \frac{p (x') / q(x')}{p (x) / q (x)} \right\rbrace
\end{align*} α ( x , x ′ ) = { 1 , p ( x ) / q ( x ) p ( x ′ ) / q ( x ′ ) } 這種抽樣方式簡單易實現，但收斂速度較慢，且一般 q ( x ) q (x) q ( x ) p ( x ) p (x) p ( x )
Full Conditional Distribution 若採樣目標為 k k k x = ( x 1 , x 2 , ⋯ , x k ) T x = (\contia{x}{k})^T x = ( x 1 , x 2 , ⋯ , x k ) T I ⊂ K = { 1 , 2 , ⋯ , k } I \subset K = \{ \conti{k} \} I ⊂ K = { 1 , 2 , ⋯ , k } p ( x I ∣ x − I ) p (x_{I} | x_{-I}) p ( x I ∣ x − I ) x I = { x i ∣ i ∈ I } x_{I} = \{ x_i | i \in I \} x I = { x i ∣ i ∈ I } x − I = { x i ∣ i ∉ I } x_{-I} = \{ x_i | i \not\in I \} x − I = { x i ∣ i ∈ I }
Full conditional distribution 擁有以下性質
p ( x I ∣ x − I ) = p ( x ) ∫ p ( x ) d x I ∝ p ( x ) \begin{align*}
p (x_I | x_{-I})
= \frac{p (x)}{\int p (x) dx_I}
\propto p (x)
\end{align*} p ( x I ∣ x − I ) = ∫ p ( x ) d x I p ( x ) ∝ p ( x ) 而任意兩 x , x ′ x, x' x , x ′
p ( x I ′ ∣ x − I ′ ) p ( x I ∣ x − I ) = p ( x ′ ) p ( x ) \begin{align*}
\frac{p (x_I' | x_{-I}') }{p (x_I | x_{-I}) }
= \frac{p (x')}{p (x)}
\end{align*} p ( x I ∣ x − I ) p ( x I ′ ∣ x − I ′ ) = p ( x ) p ( x ′ ) 這項性質被用於改進 Metropolis–Hasting Algorithm，其中計算 p ( x ′ ) p ( x ) \frac{p (x')}{p (x)} p ( x ) p ( x ′ ) p ( x I ′ ∣ x − I ′ ) p ( x I ∣ x − I ) \frac{p (x_I' | x_{-I}') }{p (x_I | x_{-I}) } p ( x I ∣ x − I ) p ( x I ′ ∣ x − I ′ )
Single-Component Metropolis–Hastings Algorithm 在 Metropolis–Hastings Algorithm 需對多變量分佈採樣，有時這種採樣非常困難，因此轉而對單一變數的條件分佈採樣，並遍歷所有變數以此簡化問題，這種方式被稱為 single-component Metropolis–Hastings algorithm
設採樣目標為 k k k x = ( x 1 , x 2 , ⋯ , x k ) T x = (\contia{x}{k})^T x = ( x 1 , x 2 , ⋯ , x k ) T x j x_j x j x x x j j j i i i
x ( i ) = ( x 1 ( i ) , x 2 ( i ) , ⋯ , x k ( i ) ) T \begin{align*}
x^{(i)}
= (x_1^{(i)}, x_2^{(i)}, \cdots, x_k^{(i)})^T
\end{align*} x ( i ) = ( x 1 ( i ) , x 2 ( i ) , ⋯ , x k ( i ) ) T 其中 x j ( i ) x_{j}^{(i)} x j ( i ) i i i x x x j j j
為了獲取採樣結果 { x ( 1 ) , x ( 2 ) , ⋯ , x ( n ) } \{ x^{(1)}, x^{(2)}, \cdots, x^{(n)} \} { x ( 1 ) , x ( 2 ) , ⋯ , x ( n ) } i − 1 i - 1 i − 1 j j j x j ( i − 1 ) x_j^{(i - 1)} x j ( i − 1 ) i i i j j j x j ′ ( i ) x_j'^{(i)} x j ′ ( i ) q ( x j ( i − 1 ) , x j ′ ( i ) ∣ x − j ( i ) ) q (x_j^{(i - 1)}, x_j'^{(i)} | x_{-j}^{(i)}) q ( x j ( i − 1 ) , x j ′ ( i ) ∣ x − j ( i ) ) x − j ( i ) x_{-j}^{(i)} x − j ( i ) j j j x ( i ) x^{(i)} x ( i ) j − 1 j - 1 j − 1
x − j ( i ) = ( x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) T \begin{align*}
x_{-j}^{(i)}
= (x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_{k}^{(i - 1)})^T
\end{align*} x − j ( i ) = ( x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) T 接著決定是否接受此次更新，即
α ( x j ( i − 1 ) , x j ′ ( i ) ∣ x − j ( i ) ) = min { 1 , p ( x j ′ ( i ) ∣ x − j ( i ) ) q ( x j ′ ( i ) , x j ( i − 1 ) ∣ x − j ( i ) ) p ( x j ( i − 1 ) ∣ x − j ( i ) ) q ( x j ( i − 1 ) , x j ( i ) ∣ x − j ( i ) ) } \begin{align*}
\alpha (x_{j}^{(i - 1)}, x_{j}'^{(i)} | x_{-j}^{(i)} )
= \min \left\lbrace 1, \frac{p (x_j'^{(i)} | x_{-j}^{(i)}) q (x_j'^{(i)}, x_j^{(i - 1)} | x_{-j}^{(i)})}{p (x_j^{(i - 1)} | x_{-j}^{(i)}) q (x_j^{(i - 1)}, x_j^{(i)} | x_{-j}^{(i)})} \right\rbrace
\end{align*} α ( x j ( i − 1 ) , x j ′ ( i ) ∣ x − j ( i ) ) = min { 1 , p ( x j ( i − 1 ) ∣ x − j ( i ) ) q ( x j ( i − 1 ) , x j ( i ) ∣ x − j ( i ) ) p ( x j ′ ( i ) ∣ x − j ( i ) ) q ( x j ′ ( i ) , x j ( i − 1 ) ∣ x − j ( i ) ) } 若接受此採樣，則 x j ( i ) = x j ′ ( i ) x_j^{(i)} = x_j'^{(i)} x j ( i ) = x j ′ ( i ) x j ( i ) = x j ( i − 1 ) x_j^{(i)} = x_{j}^{(i - 1)} x j ( i ) = x j ( i − 1 )
Gibbs Sampling Gibbs sampling 是一種 MCMC 常用的方式，是 single-component Metropolis–Hastings algorithm 的一種特例，但更容易實現，常用於估計多變量參數。其概念是從 joint probability distribution 定義 full conditional distribution，並在其中採樣。
給定 k k k p ( x 1 , x 2 , ⋯ , x k ) p (\contia{x}{k}) p ( x 1 , x 2 , ⋯ , x k ) x ( 0 ) = ( x 1 ( 0 ) , x 2 ( 0 ) , ⋯ , x k ( 0 ) ) T x^{(0)} = (x_1^{(0)}, x_2^{(0)}, \cdots, x_k^{(0)})^T x ( 0 ) = ( x 1 ( 0 ) , x 2 ( 0 ) , ⋯ , x k ( 0 ) ) T i i i j j j x j ( i ) x_j^{(i)} x j ( i ) p ( x j ∣ x − j ( i ) ) p (x_j | x_{-j}^{(i)}) p ( x j ∣ x − j ( i ) )
p ( x j ∣ x − j ( i ) ) = p ( x j ∣ x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) \begin{align*}
p (x_j | x_{- j}^{(i)})
= p (x_j | x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_{k}^{(i - 1)})
\end{align*} p ( x j ∣ x − j ( i ) ) = p ( x j ∣ x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) 其中定義 proposed distribution 是 full conditional distribution，即
q ( x , x ′ ) = p ( x j ′ ∣ x − j ) \begin{align*}
q (x, x')
= p (x_j' | x_{-j})
\end{align*} q ( x , x ′ ) = p ( x j ′ ∣ x − j ) 透過這個轉換使得 α ( x , x ′ ) = 1 \alpha (x, x') = 1 α ( x , x ′ ) = 1
Gibbs sampling 適用於 full conditional distribution 容易採樣，因此能提高採樣效率；single-component Metropolis–Hastings algorithm 即便 full conditional distribution 不易採樣也適用，但效率較低。
演算法 目標：給定 pdf p ( x ) p (x) p ( x ) f ( x ) f (x) f ( x ) x m + 1 , x m + 2 , ⋯ , x n x_{m + 1}, x_{m + 2}, \cdots, x_n x m + 1 , x m + 2 , ⋯ , x n f m n f_{mn} f mn 輸入：pdf p ( x ) p (x) p ( x ) f ( x ) f (x) f ( x ) m m m n n n 輸出：採樣結果 x m + 1 , x m + 2 , ⋯ , x n x_{m + 1}, x_{m + 2}, \cdots, x_n x m + 1 , x m + 2 , ⋯ , x n f m n f_{mn} f mn 設定起始點 x ( 0 ) = ( x 1 ( 0 ) , , x 2 ( 0 ) , ⋯ , x k ( 0 ) ) T x^{(0)} = (x_1^{(0)}, ,x_2^{(0)}, \cdots, x_k^{(0)})^T x ( 0 ) = ( x 1 ( 0 ) ,, x 2 ( 0 ) , ⋯ , x k ( 0 ) ) T
設 i = 1 , 2 , ⋯ , n i = \conti{n} i = 1 , 2 , ⋯ , n
採樣 x 1 ( i ) x_1^{(i)} x 1 ( i ) p ( x 1 ∣ x 2 ( i − 1 ) , x 3 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) p (x_1 | x_2^{(i - 1)}, x_3^{(i - 1)}, \cdots, x_k^{(i - 1)}) p ( x 1 ∣ x 2 ( i − 1 ) , x 3 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) … 採樣 x j ( i ) x_j^{(i)} x j ( i ) p ( x j ∣ x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) p (x_j | x_1^{(i)}, \cdots, x_{j - 1}^{(i)}, x_{j + 1}^{(i - 1)}, \cdots, x_k^{(i - 1)}) p ( x j ∣ x 1 ( i ) , ⋯ , x j − 1 ( i ) , x j + 1 ( i − 1 ) , ⋯ , x k ( i − 1 ) ) … 採樣 x k ( i ) x_k^{(i)} x k ( i ) p ( x k ∣ x 1 ( i ) , x 2 ( i ) , ⋯ , x k − 1 ( i ) ) p (x_k | x_1^{(i)}, x_2^{(i)}, \cdots, x_{k - 1}^{(i)}) p ( x k ∣ x 1 ( i ) , x 2 ( i ) , ⋯ , x k − 1 ( i ) ) 獲取採樣集合 { x ( m + 1 ) , x ( m + 2 ) , ⋯ , x ( n ) } \{ x^{(m + 1)}, x^{(m + 2)}, \cdots, x^{(n)} \} { x ( m + 1 ) , x ( m + 2 ) , ⋯ , x ( n ) }
計算樣本平均
f m n = 1 n − m ∑ i = m + 1 n f ( x ( i ) ) \begin{align*}
f_{mn}
= \frac{1}{n - m} \sum_{i = m + 1}^{n} f (x^{(i)})
\end{align*} f mn = n − m 1 i = m + 1 ∑ n f ( x ( i ) ) 範例 用 Gibbs sampling 對二元常態分佈取樣，設
X = ( X 1 , X 2 ) T ∼ N ( ( 0 0 ) , ( 1 ρ ρ 1 ) ) \begin{align*}
X
= (X_1, X_2)^T
\sim N \left( \begin{pmatrix}
0 \\
0
\end{pmatrix}, \begin{pmatrix}
1 & \rho \\
\rho & 1
\end{pmatrix} \right)
\end{align*} X = ( X 1 , X 2 ) T ∼ N ( ( 0 0 ) , ( 1 ρ ρ 1 ) ) 則 full conditional distribution 為
X 1 ∣ X 2 ∼ N ( ρ X 1 , ( 1 − ρ 2 ) ) X 2 ∣ X 1 ∼ N ( ρ X 2 , ( 1 − ρ 2 ) ) \begin{align*}
X_1 | X_2
& \sim N (\rho X_1, (1 - \rho^2)) \\
X_2 | X_1
& \sim N (\rho X_2, (1 - \rho^2))
\end{align*} X 1 ∣ X 2 X 2 ∣ X 1 ∼ N ( ρ X 1 , ( 1 − ρ 2 )) ∼ N ( ρ X 2 , ( 1 − ρ 2 )) 接著可透過 Gibbs sampling 採樣，設起始狀態 x ( 0 ) = ( x 1 ( 0 ) , x 2 ( 0 ) ) T x^{(0)} = (x_1^{(0)}, x_2^{(0)})^T x ( 0 ) = ( x 1 ( 0 ) , x 2 ( 0 ) ) T
第 1 次迭代:x 1 ( 1 ) ∼ N ( ρ x 2 ( 0 ) , ( 1 − ρ 2 ) ) x_1^{(1)} \sim N (\rho x_2^{(0)}, (1 - \rho^2)) x 1 ( 1 ) ∼ N ( ρ x 2 ( 0 ) , ( 1 − ρ 2 )) x 2 ( 1 ) ∼ N ( ρ x 1 ( 1 ) , ( 1 − ρ 2 ) ) x_2^{(1)} \sim N (\rho x_1^{(1)}, (1 - \rho^2)) x 2 ( 1 ) ∼ N ( ρ x 1 ( 1 ) , ( 1 − ρ 2 )) … 第 i i i x 1 ( i ) ∼ N ( ρ x 2 ( i − 1 ) , ( 1 − ρ 2 ) ) x_1^{(i)} \sim N (\rho x_2^{(i - 1)}, (1 - \rho^2)) x 1 ( i ) ∼ N ( ρ x 2 ( i − 1 ) , ( 1 − ρ 2 )) x 2 ( i ) ∼ N ( ρ x 1 ( i ) , ( 1 − ρ 2 ) ) x_2^{(i)} \sim N (\rho x_1^{(i)}, (1 - \rho^2)) x 2 ( i ) ∼ N ( ρ x 1 ( i ) , ( 1 − ρ 2 )) … 第 i i i x 1 ( n ) ∼ N ( ρ x 2 ( n − 1 ) , ( 1 − ρ 2 ) ) x_1^{(n)} \sim N (\rho x_2^{(n - 1)}, (1 - \rho^2)) x 1 ( n ) ∼ N ( ρ x 2 ( n − 1 ) , ( 1 − ρ 2 )) x 2 ( n ) ∼ N ( ρ x 1 ( n ) , ( 1 − ρ 2 ) ) x_2^{(n)} \sim N (\rho x_1^{(n)}, (1 - \rho^2)) x 2 ( n ) ∼ N ( ρ x 1 ( n ) , ( 1 − ρ 2 )) 選定收斂次數 m m m n n n { x ( m + 1 ) , x ( m + 2 ) , ⋯ , x ( n ) } \{ x^{(m + 1)}, x^{(m + 2)}, \cdots, x^{(n)} \} { x ( m + 1 ) , x ( m + 2 ) , ⋯ , x ( n ) }
採樣計算 在 Gibbs sampling 須對 full conditional distribution 進行多次採樣，因此能透過分佈的特徵來改進或加速採樣效率。以 Bayesian learning 作為範例，設 y y y x = ( α , θ , z ) x = (\alpha, \theta, z) x = ( α , θ , z ) α \alpha α θ \theta θ z z z
p ( x ∣ y ) = p ( α , θ , z ∣ y ) ∝ p ( z , y ∣ θ ) p ( θ ∣ α ) p ( α ) \begin{align*}
p (x | y)
= p (\alpha, \theta, z | y)
\propto p (z, y | \theta) p (\theta | \alpha) p (\alpha)
\end{align*} p ( x ∣ y ) = p ( α , θ , z ∣ y ) ∝ p ( z , y ∣ θ ) p ( θ ∣ α ) p ( α ) 其中 p ( z , y ∣ θ ) p (z, y | \theta) p ( z , y ∣ θ ) p ( θ ∣ α ) p (\theta | \alpha) p ( θ ∣ α ) p ( α ) p (\alpha) p ( α ) p ( x ∣ y ) p (x | y) p ( x ∣ y )
p ( α i ∣ α − i , θ , z , y ) ∝ p ( θ ∣ α ) p ( α ) p ( θ j ∣ θ − j , α , θ , z , y ) ∝ p ( z , y ∣ θ ) p ( θ ∣ α ) p ( z k ∣ z − k , α , θ , y ) ∝ p ( z , y ∣ θ ) \begin{align*}
p (\alpha_i | \alpha_{-i}, \theta, z, y)
& \propto p (\theta | \alpha) p (\alpha) \\
p (\theta_j | \theta_{-j}, \alpha, \theta, z, y)
& \propto p (z, y | \theta) p (\theta | \alpha) \\
p (z_k | z_{-k}, \alpha, \theta, y)
& \propto p (z, y | \theta)
\end{align*} p ( α i ∣ α − i , θ , z , y ) p ( θ j ∣ θ − j , α , θ , z , y ) p ( z k ∣ z − k , α , θ , y ) ∝ p ( θ ∣ α ) p ( α ) ∝ p ( z , y ∣ θ ) p ( θ ∣ α ) ∝ p ( z , y ∣ θ ) full conditional probability 是多個 conditional probability 相乘，每個機率分佈僅與少數變數相關，因此對 full conditional probability 可以僅透過數個 conditional probability 的採樣來執行，以此降低運算複雜度
參考資料 LI, Hang. Machine Learning Methods. Springer Nature, 2023. Lange, Kenneth, J. Chambers, and W. Eddy. Numerical analysis for statisticians. Vol. 1. New York: springer, 2010.
