這篇回答節選自我的專欄 【機器學習中的數學：機率圖與隨機過程】 ，我們來仔細介紹和分析一下Gibbs采樣。
也歡迎關註我的知乎帳號 @石溪 ，將持續釋出機器學習數學基礎及演算法套用等方面的精彩內容。
1.Gibbs方法核心流程
在 M-H 演算法的基礎上，我們再介紹一下 Gibbs 采樣。
吉布士采樣是一種針對高維分布的采樣方法，假設我們待采樣的高維目標分布為 p(z_1,z_2,z_3,...,z_m) ，吉布士采樣的目標就是從這個高維分布中采樣出一組樣本 z^{(1)},z^{(2)},z^{(3)},...,z^{(N)} ，使得這一組樣本能夠服從 m 維分布 p(z) 。
我們先約定一個記法：按照上述假設，從分布中采得的樣本 z 是一個擁有 m 維內容值的樣本值， z_i 表示某個樣本 z 的第 i 維內容， z_{-i} 表示除開第 i 維樣本外，剩余的 m-1 維內容： z_1,z_2,...,z_{i-1},z_{i+1},...,z_m
具體如何采樣，核心在七個字：一維一維的采樣：
那麽何謂一維一維的采樣呢？我們以三維隨機變量 z 為例： p(z)=p(z_1,z_2,z_3)
首先在0時刻，我們給定一個初值，也就是 z^{(0)} 的三個內容的值： z^{(0)}_1,z^{(0)}_2,z^{(0)}_3
好，下面我們來看下一個時刻，也就是時刻1的采樣值如何生成：
z_1^{(1)}\sim p(z_1|z_2^{(0)},z_3^{(0)})
z_2^{(1)}\sim p(z_2|z_1^{(1)},z_3^{(0)})
z_3^{(1)}\sim p(z_2|z_1^{(1)},z_2^{(1)})
這樣就獲得了時刻1的包含三維內容值的完整采樣值。
那麽接著往下遞迴，概況起來，透過 t 時刻的采樣值遞推到 t+1 時刻，也是如上述方法那樣一維一維的生成采樣值的各個內容：固定其他維的內容值，透過條件機率，依機率生成 t+1 時刻的某一維度的內容值。
z_1^{(t+1)}\sim p(z_1|z_2^{(t)},z_3^{(t)})
z_2^{(t+1)}\sim p(z_2|z_1^{(t+1)},z_3^{(t)})
z_3^{(t+1)}\sim p(z_2|z_1^{(t+1)},z_2^{(t+1)})
這樣就由$t$時刻的采樣值 z^{(t)} 得到了下一時刻 t+1 的采樣值 z^{(t+1)}
按照該方法迴圈下去，就可以采樣出 N 個樣本值： z^{(1)},z^{(2)},z^{(3)},...,z^{(N)} ，丟棄掉前面一部份燃燒期的樣本，剩下的樣本就服從我們的目標高維分布。
2.Gibbs采樣的合理性
這就是Gibbs采樣的核心流程，有兩個大的問題。
第一個問題是：為什麽透過這種采樣方法采得的樣本能夠服從目標分布 p(z) ？
我們拿Gibbs采樣和M-H方法進行類比，我們會發現其實Gibbs采樣就是M-H采樣的一種特殊情況，這個結果就證明了Gibbs的可行性，註意邏輯是這樣的，M-H我們已經證明是符合細致平穩條件的，而Gibbs是M-H的一種特殊情況，因此Gibbs也符合細致平穩條件，因此可行。
首先，我們采樣的目標分布是 p(z) ，而我們實際使用的提議轉移過程 Q 就是機率 p 本身： Q(z^*,z)=p(z_i|z_{-i}^*) ，你發現了沒，在這個采樣過程，好像沒有接受率 \alpha 的出現。
我們重新看看接受率 \alpha(z^*,z) 的定義式：\alpha(z^*,z)=min(1,\frac{p(z^*)Q(z^*,z)}{p(z)Q(z,z^*)})
重點就在 \frac{p(z^*)Q(z^*,z)}{p(z)Q(z,z^*)} 這個式子： p(z)=p(z_i,z_{-i})=p(z_i|z_{-i})p(z_{-i})
\frac{p(z^*)Q(z^*,z)}{p(z)Q(z,z^*)}=\frac{[p(z_i^*|z_{-i}^*)p(z_{-i}^*)]p(z_i|z_{-i}^*)}{[p(z_i|z_{-i})p(z_{-i})]p(z_i^*|z_{-i})}
這個式子如何化簡？我們發現，在上面的轉移過程 z\rightarrow z^* 當中，我們是透過固定采樣值 z 的 -i 維內容，單采 z^* 的第 i 維內容值，那麽在一次內容采樣前後 z_{-i} 和 z_{-i}^* 是一碼事兒，那我們把上面式子中所有的 z_{-i}^* 都替換成 z_{-i} ，就有：
\frac{p(z^*)Q(z^*,z)}{p(z)Q(z,z^*)}=\frac{[p(z_i^*|z_{-i}^*)p(z_{-i}^*)]p(z_i|z_{-i}^*)}{[p(z_i|z_{-i})p(z_{-i})]p(z_i^*|z_{-i})}\\=\frac{[p(z_i^*|z_{-i})p(z_{-i})]p(z_i|z_{-i})}{[p(z_i|z_{-i})p(z_{-i})]p(z_i^*|z_{-i})}=1
也就是說Gibbs采樣中的接受率 \alpha 為：
\alpha(z^*,z)=min(1,\frac{p(z^*)Q(z^*,z)}{p(z)Q(z,z^*)})=min(1,1)=1
因此滿足細致平穩條件的原始等式：
p(z)Q(z,z^*)\alpha(z,z^*)=p(z^*)Q(z^*,z)\alpha(z^*,z) 就等效為：
p(z)Q(z,z^*)*1=p(z^*)Q(z^*,z)*1\Rightarrow p(z)Q(z,z^*)=p(z^*)Q(z^*,z) 其中： Q(z^*,z)=p(z_i|z_{-i}^*)
因此 p(z) 和 Q(z,z^*) 滿足細致平穩條件，按照 Q 描述的狀態轉移過程所得到的采樣樣本集，在數量足夠多的情況下，就近似為目標分布 p(z)
但是我們看出來了，我們把提議轉移機率矩陣 Q 就定位了 p 的滿條件機率，那麽很顯然就有一個限制條件，那就是高維分布 p 的所有滿條件機率都是易於采樣的。
關註 @石溪 知乎帳號，分享更多機器學習數學基礎精彩內容。
3.Gibbs采樣程式碼試驗
說了這麽多，我們舉個具體的例子，比如目標分布是一個二維的高斯分布，他的均值向量： \mu=\begin{bmatrix}3\\3\end{bmatrix} ，共變異數矩陣 \Sigma=\begin{bmatrix}1&0.6\\0.6&1\end{bmatrix} ，即$z$包含兩維內容 (z_1,z_2) ， z \sim N(\mu, \Sigma)
我們透過吉布士采樣的方法，來獲取服從上述分布的一組樣本值：
很顯然，這個二維高斯分布的目標分布是可以進行Gibbs采樣的，因為高斯分布有一個好的性質，我們之前反復說過，那就是高斯分布的條件分布也是高斯分布，因此保證了一維一維采樣的過程是順暢的。
我們先給出二維高斯分布的條件機率分布公式：
z 包含兩維內容 (z_1,z_2) ， z \sim N(\mu, \Sigma) ，其中 \mu=\begin{bmatrix}\mu_1\\\mu_2\end{bmatrix} ， \Sigma=\begin{bmatrix}\delta_{11}&\delta_{12}\\\delta_{21}&\delta_{22}\end{bmatrix}
那麽條件分布 z_1|z_2 和 z_2|z_1 分別服從下面的形式：
z_1|z_2\sim N(\mu_1+\frac{\delta_{12}}{\delta_{22}}(z_2-\mu_2),\delta_{11}-\frac{\delta^2_{12}}{\delta_{22}})
z_2|z_1\sim N(\mu_2+\frac{\delta_{12}}{\delta_{11}}(z_1-\mu_2),\delta_{22}-\frac{\delta^2_{12}}{\delta_{11}})
這就是二維高斯分布采樣過程中，逐維采樣所依從的滿條件分布，他是高斯分布，因此操作尤其簡便，下面我們來看程式碼：
程式碼片段：
import
numpy
as
np
import
matplotlib.pyplot
as
plt
import
seaborn
seaborn
.
set
()
#依照z1|z2的條件高斯分布公式，給定z2的條件情況下采樣出z1
def
p_z1_given_z2
(
z2
,
mu
,
sigma
):
mu
=
mu
[
0
]
+
sigma
[
1
][
0
]
/
sigma
[
0
][
0
]
*
(
z2
-
mu
[
1
])
sigma
=
sigma
[
0
][
0
]
-
sigma
[
1
][
0
]
**
2
/
sigma
[
1
][
1
]
return
np
.
random
.
normal
(
mu
,
sigma
)
#依照z2|z1的條件高斯分布公式，給定z1的條件情況下采樣出z2
def
p_z2_given_z1
(
z1
,
mu
,
sigma
):
mu
=
mu
[
1
]
+
sigma
[
0
][
1
]
/
sigma
[
1
][
1
]
*
(
z1
-
mu
[
0
])
sigma
=
sigma
[
1
][
1
]
-
sigma
[
0
][
1
]
**
2
/
sigma
[
0
][
0
]
return
np
.
random
.
normal
(
mu
,
sigma
)
#Gibbs采樣過程
def
gibbs_sampling
(
mu
,
sigma
,
samples_period
):
samples
=
np
.
zeros
((
samples_period
,
2
))
z2
=
np
.
random
.
rand
()
*
10
for
i
in
range
(
samples_period
):
z1
=
p_z1_given_z2
(
z2
,
mu
,
sigma
)
z2
=
p_z2_given_z1
(
z1
,
mu
,
sigma
)
samples
[
i
,
:]
=
[
z1
,
z2
]
return
samples
#目標分布p(z)
mus
=
np
.
array
([
3
,
3
])
sigmas
=
np
.
array
([[
1
,
.
7
],
[
.
7
,
1
]])
#確定總的采樣期和燃燒期
burn_period
=
int
(
1e4
)
samples_period
=
int
(
1e5
)
#得到采樣樣本，舍棄燃燒期樣本點
samples
=
gibbs_sampling
(
mus
,
sigmas
,
samples_period
)[
burn_period
:]
plt
.
plot
(
samples
[:,
0
],
samples
[:,
1
],
'ro'
,
alpha
=
0.05
,
markersize
=
1
)
plt
.
show
()
執行結果：
這個采樣結果的分布，看上去是滿意的，以上就是Gibbs采樣的原理和演示過程。
此內容節選自我的專欄【機器學習中的數學：機率圖與隨機過程】,前三節可免費試讀，歡迎訂閱！
當然還有【機器學習中的數學（全集）】系列專欄，歡迎大家閱讀，配合食用，效果更佳~
有訂閱的問題可咨詢微信：zhangyumeng0422
