Gibbs Sampling

Gibbs Sampling

1 Gibbs算法

1.1 二阶段Gibbs

已知二维随机向量$(X_1,X_2)$服从联合概率分布$f(x_1,x_2)$,对应的边缘概率密度函数为
$$\begin{gathered} f(x_1)=\int f(x_1,x_2)d{x}_2 \\ f(x_2)=\int f(x_1,x_2)d{x}_1; \end{gathered}$$
则可以计算出随机变量$X_1$和$X_2$的条件概率分布：
$$\begin{array}{r}f(x_1|x_2)=\frac{f(x_1,x_2)}{f(x_2)}=\frac{f(x_1,x_2)}{\int f(x_1,x_2)d{x}_1}\\ f(x_2|x_1)=\frac{f(x_1,x_2)}{f(x_1)}=\frac{f(x_1,x_2)}{\int f(x_1,x_2)d{x}_2}\end{array}$$
二阶段Gibbs算法如下：

（1）给定初始值$(x_1^0,x_2^0)$，

（2）对应$t=2,3\dots N$,

• 从条件密度$f(x_1|x_2^{t-1})$采样得到$x_1^t$
• 从条件密度$f(x_2|x_1^t)$采样得到$x_2^t$

根据上面的算法得到样本 { x t = ( x 1 t , x 2 t ) ∣ 0 ≤ t ≤ N } \{\boldsymbol{x^t} = (x_1^t,x_2^t)|0\le t\le N\} {xt=(x1t​,x2t​)∣0≤t≤N}。不难发现，下一期$x^{t+1}$仅仅取决于当前$x^t$的取值，而与过去时刻取值无关，这种性质具有马尔可夫链的性质。因此，二阶段Gibbs采样也属于马尔科夫链蒙特卡洛模拟(MCMC)采样。为保证采样结果收敛到目标分布，前一段时间的采样应该舍去。

---

1.2 多阶段Gibbs

类似地可以得到多阶段Gibbs算法：

（1）给定初始值$(x_1^0,x_2^0\dots x_n^0)$，

（2）对应$t=2,3\dots T$,

• 从条件密度$f(x_1|x_2^{t-1},x_3^{t-1}\dots x_n^{t-1})$采样得到$x_1^t$;
• 从条件密度$f(x_2|x_1^t,x_3^{t-1}\dots x_n^{t-1})$采样得到$x_2^t$;
• ……
• 从条件密度$f(x_n|x_1^t,x_3^t\dots x_{n-1}^t)$采样得到$x_n^t$;

为保证采样结果收敛到目标分布，前一段时间的采样应该舍去。

---

2 代码实现

以二元正态分布为例。已知二元随机变量$(x_1,x_2)$满足二元联合正态分布，均值分别为${\mu }_1=2,{\mu }_2=1.5$，方差为${\sigma }_1=3,{\sigma }_2=5$，相关系数$\rho =-0.5$。使用Gibbs方法进行采样。对于一般二维正态分布
$$f(x,y)=\frac{1}{2\pi {\sigma }_x{\sigma }_y\sqrt{1-{\rho }^2}}\exp \left(-\frac{1}{2\left(1-{\rho }^2\right)}\left[\frac{{\left(x-{\mu }_x\right)}^2}{{\sigma }_x^2}-2\rho \frac{\left(x-{\mu }_x\right)\left(y-{\mu }_y\right)}{{\sigma }_x{\sigma }_y}+\frac{{\left(y-{\mu }_y\right)}^2}{{\sigma }_y^2}\right]\right)$$
在$x$条件下。$y$的条件分布也为正态分布(一维)：
$$f(y\mid x)=\frac{1}{\sqrt{2\pi }{\sigma }_{y\mid x}}\exp \left(-\frac{{\left(y-{\mu }_{y\mid x}\right)}^2}{2{\sigma }_{y\mid x}^2}\right)$$
其中条件均值：
$${\mu }_{y\mid x}={\mu }_y+\rho \frac{{\sigma }_y}{{\sigma }_x}\left(x-{\mu }_x\right)$$
条件方差：
$${\sigma }_{y\mid x}={\sigma }_y\sqrt{1-{\rho }^2}$$
即
$$\begin{array}{r}y|x\sim N({\mu }_y+\rho \frac{{\sigma }_y}{{\sigma }_x}\left(x-{\mu }_x\right),{\sigma }_y\sqrt{1-{\rho }^2})\\ x|y\sim N({\mu }_x+\rho \frac{{\sigma }_x}{{\sigma }_y}\left(y-{\mu }_x\right),{\sigma }_x\sqrt{1-{\rho }^2})\end{array}$$

#二阶段Gibbs抽样
x1  <- numeric()
x2 <- numeric()
rho = 0.6;mu1 = 1.1;sigma1 = 3;mu2 = 1.8;sigma2 = 4
x1[1] = 0
x2[1] = 0
N = 5000
for(t in 2:N){
  print(t)
  x1[t] = rnorm(1,mean = mu1+rho*sigma1/sigma2*(x2[t-1]-mu2),
                sd = sqrt(1-rho^2)*sigma1)
  x2[t] = rnorm(1,mean = mu2+rho*sigma2/sigma1*(x1[t]-mu1),
                sd = sqrt(1-rho^2)*sigma2)
}
cor(x1,x2)
plot(x1,x2)
grid(col = "black")
# 样本期望、方差
mean(x1[1001:2000])
mean(x2[1001:2000])
sd(x1[1001:2000])
sd(x2[1001:2000])

---

参考书籍：

黄长全，贝叶斯统计及其R实现[M]. 北京: 清华大学出版社, 2017.