MCMC-Gibbs抽样(R语言实现)

在贝叶斯统计中，不同于频率学派将参数视作为未知常数，贝叶斯学派将参数视作为随机变量。

若参数 θ \theta θ的先验分布为 P ( θ ) P(\theta) P(θ)，则可以根据获得的样本数据X和贝叶斯定理去跟新对于先验分布 P ( θ ) P(\theta) P(θ)的认知，得到参数 θ \theta θ的后验分布：
P ( θ ∣ X ) = P ( θ ) P ( X ∣ θ ) P ( X ) P(\theta|X)=\frac{P(\theta)P(X|\theta)}{P(X)} P(θ∣X)=P(X)P(θ)P(X∣θ)​
由于等式右侧的分母部分相较于后验分布 P ( θ ∣ X ) P(\theta|X) P(θ∣X)为一常数，在实际求贝叶斯后验分布时可将其如下正比于去除，方便推导后验分布。
P ( θ ∣ X ) ∝ P ( θ ) P ( X ∣ θ ) P(\theta|X)\propto P(\theta)P(X|\theta) P(θ∣X)∝P(θ)P(X∣θ)
基于上述参数的后验分布，可以进行统计推断。

但是问题在于对于多参数的联合后验分布往往没有标准形式，很难直接采样。因此可考虑从每个参数的完全条件分布采样，比如Gibbs抽样方法。

现具体介绍Gibbs抽样：

Gibbs抽样是一种马尔科夫链蒙特卡洛方法，即MCMC方法，马尔科夫链进行采样，蒙特卡洛方法进行近似积分，是一种结合蒙特卡洛方法和马尔可夫性的迭代算法。

假设参数向量为 ϕ = ( ϕ 1 , . . . , ϕ d ) ′ \phi = (\phi_1,...,\phi_d)' ϕ=(ϕ1​,...,ϕd​)′，在第t步迭代Gibbs样本时，对于每个参数 ϕ j t     ( j = 1 , . . . , d ) \phi_j^t\ \ \ (j=1,...,d) ϕjt​   (j=1,...,d)来自于条件分布 p ( ϕ j ∣ ϕ − j t − 1 , y ) p(\phi_j|\phi_{-j}^{t-1},y) p(ϕj​∣ϕ−jt−1​,y)，其中 ϕ − j t − 1 = ( ϕ 1 t , ϕ 2 t , . . . , ϕ j − 1 t , ϕ j + 1 t − 1 , . . . , ϕ d t − 1 ) \phi_{-j}^{t-1}=(\phi_1^t,\phi_2^t,...,\phi_{j-1}^t,\phi_{j+1}^{t-1},...,\phi_{d}^{t-1}) ϕ−jt−1​=(ϕ1t​,ϕ2t​,...,ϕj−1t​,ϕj+1t−1​,...,ϕdt−1​)，即前j-1个参数已跟新到第t步，第j个及第j个参数之后的参数跟新到第t-1，还未跟新到第t步，此时待跟新到第t步的参数为第j个参数 ϕ j \phi_j ϕj​，这个迭代算法直至收敛至目标的联合后验分布 P ( ϕ 1 , . . . , ϕ d ∣ y ) P(\phi_1,...,\phi_d|y) P(ϕ1​,...,ϕd​∣y).

R语言实践

Gibbs抽样获取二元正态分布的样本：

( X 1 X 2 ) ∼ ( [ μ 1 μ 2 ] [ σ 11 2 ρ σ 11 σ 22 ρ σ 11 σ 22 σ 22 2 ] ) \begin{pmatrix} X_1\\X_2 \end{pmatrix}\sim\begin{pmatrix} \begin{bmatrix} \mu_1\\\mu_2 \end{bmatrix}&\begin{bmatrix} \sigma_{11}^2&\rho\sigma_{11}\sigma_{22} \\\rho\sigma_{11}\sigma_{22} &\sigma_{22}^2 \end{bmatrix} \end{pmatrix} (X1​X2​​)∼([μ1​μ2​​]​[σ112​ρσ11​σ22​​ρσ11​σ22​σ222​​]​)

由多元正态的条件分布得性质可知:

X 1 ∣ X 2 ∼ N ( μ 1 + ρ σ 11 σ 22 ( x 2 − μ 2 ) , ( 1 − ρ 2 ) σ 11 2 ) X_1|X_2\sim N(\mu_1+\frac{\rho\sigma_{11}}{\sigma_{22}}(x_2-\mu_2),(1-\rho^2)\sigma_{11}^2) X1​∣X2​∼N(μ1​+σ22​ρσ11​​(x2​−μ2​),(1−ρ2)σ112​)
X 2 ∣ X 1 ∼ N ( μ 2 + ρ σ 22 σ 11 ( x 1 − μ 1 ) , ( 1 − ρ 2 ) σ 22 2 ) X_2|X_1\sim N(\mu_2+\frac{\rho\sigma_{22}}{\sigma_{11}}(x_1-\mu_1),(1-\rho^2)\sigma_{22}^2) X2​∣X1​∼N(μ2​+σ11​ρσ22​​(x1​−μ1​),(1−ρ2)σ222​)

#抽样总次数
N <- 5000
#丢弃的非平稳值
burn <- 1000
#存放采样值
X <- matrix(0,N,2)
#初始值
rho <- -0.5
mu1 <- 0
mu2 <- 0
sigma1 <- 1
sigma2 <- 2
s1 <- sqrt(1-rho^2)*sigma1
s2 <- sqrt(1-rho^2)*sigma2

X[1, ] <- c(mu1, mu2)
#迭代
for (i in 2:N) {
x2 <- X[i-1, 2]
m1 <- mu1 + rho * (x2 - mu2) * sigma1/sigma2
X[i, 1] <- rnorm(1, m1, s1)
x1 <- X[i, 1]
m2 <- mu2 + rho * (x1 - mu1) * sigma2/sigma1
X[i, 2] <- rnorm(1, m2, s2)
}
#丢掉未收敛的部分

x <- X[1000:N, ]

基于这个Gibbs采样样本可以计算联合分布的均值,变量间的相关系数等。

colMeans(X)

cor(X)

例二

采用Gibbs抽样跟新后验参数

以正态分布为例，假定联合先验是相互独立的，即 P ( θ , σ 2 ) = p ( θ ) p ( σ 2 ) P(\theta,\sigma^2)=p(\theta)p(\sigma^2) P(θ,σ2)=p(θ)p(σ2)，并服从半共轭先验分布：

θ ∼ N ( μ 0 , τ 0 2 ) \theta\sim N(\mu_0,\tau_0^2) θ∼N(μ0​,τ02​)

σ 2 ∼ I n v G a m m a ( v 0 / 2 , v 0 σ 0 2 / 2 ) \sigma^2\sim InvGamma(v_0/2,v_0\sigma_0^2/2) σ2∼InvGamma(v0​/2,v0​σ02​/2)

此时 σ 2 \sigma^2 σ2的边缘后验分布并无标准形式，但是可以较为容易的推出其条件后验分布如下：

p ( σ 2 ∣ θ , y 1 , . . . , y n ) ∝ p ( σ 2 ) p ( y 1 , . . . , y n ∣ θ , σ 2 ) \begin{align} p(\sigma^2|\theta,y_1,...,y_n)&\propto p(\sigma^2)p(y_1,...,y_n|\theta,\sigma^2)\\ \end{align} p(σ2∣θ,y1​,...,yn​)​∝p(σ2)p(y1​,...,yn​∣θ,σ2)​​

σ 2 ∣ θ , y 1 , . . . , y n ∼ I n v G a m m a ( v 0 + n 2 , ( v 0 + n ) ( v 0 σ 0 2 + ( n − 1 ) s 2 + n ( y ‾ − θ ) 2 ) 2 ( v 0 + n ) ) \sigma^2|\theta,y_1,...,y_n\sim InvGamma\left(\frac{v_0+n}{2},\frac{(v_0+n)(v_0\sigma^2_0+(n-1)s^2+n(\overline{y}-\theta)^2)}{2(v_0+n)}\right) σ2∣θ,y1​,...,yn​∼InvGamma(2v0​+n​,2(v0​+n)(v0​+n)(v0​σ02​+(n−1)s2+n(y​−θ)2)​)

按照上述思路可以同样推导得到 θ \theta θ的完全条件分布。

现执行Gibbs抽样，假设当前的参数状态 ϕ ( s ) = { θ ( s ) , σ 2 ( s ) } \phi^{(s)}=\{\theta^{(s)},\sigma^{2(s)}\} ϕ(s)={θ(s),σ2(s)}，则新的状态如下:

1:  θ ( s + 1 ) ∼ p ( θ ∣ σ 2 ( s ) , y 1 , . . . , y n ) 2:  σ 2 ( s + 1 ) ∼ p ( σ 2 ∣ θ ( s + 1 ) , y 1 , . . . , y n ) 3: 最后让 ϕ ( s + 1 ) = { θ ( s + 1 ) , σ 2 ( s + 1 ) } \begin{align} \text{1: }\theta^{(s+1)}\sim p(\theta|\sigma^{2(s)},y_1,...,y_n)\\ \text{2: }\sigma^{2(s+1)}\sim p(\sigma^2|\theta^{(s+1)},y_1,...,y_n)\\ \text{3: }\text{最后让}\phi^{(s+1)}=\{\theta^{(s+1)},\sigma^{2(s+1)}\}\end{align} 1: θ(s+1)∼p(θ∣σ2(s),y1​,...,yn​)2: σ2(s+1)∼p(σ2∣θ(s+1),y1​,...,yn​)3: 最后让ϕ(s+1)={θ(s+1),σ2(s+1)}​​
重复上述三个步骤的Gibbs抽样，我们最终会获得一个相依的参数序列 { ϕ ( 1 ) , . . . , ϕ ( s ) } \{\phi^{(1)},...,\phi^{(s)}\} {ϕ(1),...,ϕ(s)}，其中：

ϕ ( 1 ) = { θ ( 1 ) , σ 2 ( 1 ) } ϕ ( 2 ) = { θ ( s ) , σ 2 ( 2 ) } … ϕ ( s ) = { θ ( s ) , σ 2 ( s ) } \begin{align} \phi^{(1)}=\{\theta^{(1)},\sigma^{2(1)}\}\\ \phi^{(2)}=\{\theta^{(s)},\sigma^{2(2)}\}\\ \dots\\ \phi^{(s)}=\{\theta^{(s)},\sigma^{2(s)}\}\end{align} ϕ(1)={θ(1),σ2(1)}ϕ(2)={θ(s),σ2(2)}…ϕ(s)={θ(s),σ2(s)}​​
ϕ ( s ) \phi^{(s)} ϕ(s)仅通过 ϕ ( s − 1 ) \phi^{(s-1)} ϕ(s−1)与 ϕ ( 0 ) , . . . , ϕ ( s − 1 ) \phi^{(0)},...,\phi^{(s-1)} ϕ(0),...,ϕ(s−1)关联，所以 { ϕ ( 1 ) , . . . , ϕ ( s ) } \{\phi^{(1)},...,\phi^{(s)}\} {ϕ(1),...,ϕ(s)}为马尔可夫序列。

因此有

lim ⁡ s → ∞ P r ( ϕ ( s ) ∈ A ) → ∫ A p ( ϕ ) d ϕ \lim_{s\to \infty} Pr(\phi^{(s)}\in A)\overset{}{\rightarrow} \int_A p(\phi)d\phi s→∞lim​Pr(ϕ(s)∈A)→∫A​p(ϕ)dϕ

进一步对感兴趣的函数 g ( ϕ ) g(\phi) g(ϕ)求其期望 E ( g ( ϕ ) ) E(g(\phi)) E(g(ϕ))可以如下进行：

lim ⁡ n → ∞ 1 n ∑ s = 1 n g ( ϕ ( s ) ) → p E ( g ( ϕ ) ) = ∫ g ( ϕ ) p ( ϕ ) d ϕ \lim_{n\to \infty}\frac{1}{n}\sum_{s=1}^ng(\phi^{(s)})\overset{p}\rightarrow E(g(\phi)) = \int g(\phi)p(\phi)d\phi n→∞lim​n1​s=1∑n​g(ϕ(s))→pE(g(ϕ))=∫g(ϕ)p(ϕ)dϕ

我们把 1 n ∑ s = 1 n g ( ϕ ( s ) ) \frac{1}{n}\sum_{s=1}^ng(\phi^{(s)}) n1​∑s=1n​g(ϕ(s))称作 E ( g ( ϕ ) ) E(g(\phi)) E(g(ϕ))的蒙特卡洛近似。

假设初始值为mu为110，方差为5^2，收集到的数据为98, 100, 107, 110, 112, 117, 117, 120, 125, 130。

y<-c(98, 100, 107, 110, 112, 117, 117, 120, 125, 130)
mean.y <- mean(y)
var.y <- var(y)
mu0 <- 110
t20 <- 5^2
n <- 10
nu0<-1
s20<-100
#迭代次数
N <- 10000
#存放Gibbs采样值
PHI <- matrix(0,N,2)
#初始值
PHI[1,] <- phi<- c(mean.y,1/var.y)
#Gibbs抽样
for(s in 2:N){
  #从完全条件后验中生成新的theta值
  mun <- (mu0/t20+n*mean.y*phi[2])/(1/t20+n*phi[2])
  t2n <- 1/(1/t20 + n*phi[2])
  phi[1] <- rnorm(1,mun,sqrt(t2n))
  #从完全条件后验中生成新的1/sigma^2值
  nun<- nu0+n
  s2n<- (nu0*s20 + (n-1)*var.y + n*(mean.y-phi[1])^2 ) /nun
  phi[2]<- rgamma(1, nun/2, nun*s2n/2)
  PHI[s,]<-phi
}
#丢掉未收敛值
PHI.P<- PHI[1000:N,]

基于上述后验抽样的样本，估计参数95%的后验区间：

CI<-apply(PHI.P,2, function(x) quantile(x,probs=c(0.025,0.975)))
CI

可见 μ \mu μ的95%的后验区间为(106.5122,117.9113)