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ϕjt1,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}) ϕjt1=(ϕ1t,ϕ2t,...,ϕj1t,ϕj+1t1,...,ϕdt1),即前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,...,ϕdy).

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} (X1X2)([μ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) X1X2N(μ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) X2X1N(μ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) σ2InvGamma(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,...,ynInvGamma(2v0+n,2(v0+n)(v0+n)(v0σ02+(n1)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)} ϕ(s1) ϕ ( 0 ) , . . . , ϕ ( s − 1 ) \phi^{(0)},...,\phi^{(s-1)} ϕ(0),...,ϕ(s1)关联,所以 { ϕ ( 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 slimPr(ϕ(s)A)Ap(ϕ)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 nlimn1s=1ng(ϕ(s))pE(g(ϕ))=g(ϕ)p(ϕ)dϕ

我们把 1 n ∑ s = 1 n g ( ϕ ( s ) ) \frac{1}{n}\sum_{s=1}^ng(\phi^{(s)}) n1s=1ng(ϕ(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)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

统计包

谢谢您的赞赏,祝您安康

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值