一.算法推导¶

Gibbs算法其实就是单分量MH算法中，建议分布取满条件概率分布的特殊情况，即：

$$q(x_j^{(i-1)}\rightarrow x_j^{'(i)} \mid x_{-j}^{(i)})=p(x_j^{'(i)} \mid x_{-j}^{(i)})$$

$$\alpha(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})=min\{1,\frac{p(x_j^{'(i)}\mid x_{-j}^{(i)})q(x_j^{'(i)}\rightarrow x_j^{(i-1)}\mid x_{-j}^{(i)})}{p(x_j^{(i-1)}\mid x_{-j}^{(i)})q(x_j^{(i-1)}\rightarrow x_j^{'(i)}\mid x_{-j}^{(i)})}\}=min\{1,\frac{p(x_j^{'(i)}\mid x_{-j}^{(i)})p(x_j^{(i-1)}\mid x_{-j}^{(i)})}{p(x_j^{(i-1)}\mid x_{-j}^{(i)})p(x_j^{'(i)}\mid x_{-j}^{(i)})}\}=min\{1,1\}=1$$

二.Gibbs采样流程¶

（1）初始化样本$x^{(0)}=(x_1^{(0)},x_2^{(0)},...,x_k^{(0)})$

（2）对$i=1,2,...,m,m+1,...n$：

（3）返回样本集$\{x_{m+1},x_{m+2},...,x_n\}$

三.案例¶

In [1]:
import os
os.chdir('../')
from ml_models import utils
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
#定义均值，协方差
u=np.asarray([0,0])
sigma=np.asarray([[1,0.5],
[0.5,1]])

In [3]:
import copy
#采样的样本量
nums=1000
count=0
points=[]
#采样x0
point=[np.random.randn(),np.random.randn()]
points.append(point)
while count<nums:
new_point=copy.deepcopy(point)
for k in (0,1):
# 按照满条件概率分布采样
new_point[k]=(np.random.randn()+0.5*new_point[(k+1)%2])*0.75
point=new_point
points.append(point)
count+=1

In [4]:
utils.plot_contourf(np.asarray(points[-400:]),lambda x:utils.gaussian_nd(x,u,sigma),8)

D:\app\Anaconda3\lib\site-packages\matplotlib\contour.py:967: UserWarning: The following kwargs were not used by contour: 'linewidth'
s)


四.小结¶

In [ ]: