MCMC_Gibbs/MH采样的简单案例R代码实现

练习：用程序实现正态分布均值、方差的后验分布抽样。

---

题目背景

记$X_i,i=1...,n$来自一维正态$N(\mu ,{\sigma }^2)$的独立同分布样本，$\mu$具有先验分布$N(a,b)$，${\sigma }^2$具有先验分布$invGamma(c,d)$.
似然函数： π ( X ∣ μ , σ 2 ) ∝ ( σ 2 ) − n / 2 exp ⁡ { − ∑ i ( X i − μ ) 2 2 σ 2 } \pi(X|\mu,\sigma^2) \propto (\sigma^2)^{-n/2} \exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2}\} π(X∣μ,σ2)∝(σ2)−n/2exp{−i∑​2σ2(Xi​−μ)2​}

补充$invGamma(\alpha ,\lambda )$的密度函数形式：
p ( x ; α , λ ) = λ α Γ ( α ) x − ( α + 1 ) exp ⁡ { − λ x } , x ≥ 0 p(x;\alpha,\lambda)=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{-(\alpha+1)}\exp \{ -\frac{\lambda}x \}, \quad x\ge 0 p(x;α,λ)=Γ(α)λα​x−(α+1)exp{−xλ​},x≥0

Gibbs抽样（详细公式推导）

$\mu$的条件（后验）分布：*
p ( μ ∣ X , σ 2 ) ∝ π ( X ∣ μ , σ 2 ) π ( μ ) ∝ exp ⁡ { − ∑ i ( X i − μ ) 2 2 σ 2 } exp ⁡ { − ( μ − a ) 2 2 b } ∝ exp ⁡ { − ( μ − μ n ) 2 2 σ n 2 } ∼ N ( μ n , σ n 2 ) p(\mu|X,\sigma^2) \propto \pi(X|\mu,\sigma^2)\pi(\mu) \propto \exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2} \} \exp\{ -\frac{(\mu-a)^2}{2b}\} \\ \propto \exp\{ -\frac{(\mu-\mu_n)^2}{2\sigma_n^2} \} \sim N(\mu_n,\sigma_n^2) p(μ∣X,σ2)∝π(X∣μ,σ2)π(μ)∝exp{−i∑​2σ2(Xi​−μ)2​}exp{−2b(μ−a)2​}∝exp{−2σn2​(μ−μn​)2​}∼N(μn​,σn2​)
${\mu }_n=(\frac{{\sum }_i{X}_i}{{\sigma }^2}+\frac{a}{b})/(\frac{n}{{\sigma }^2}+\frac{1}{b})$,${\sigma }_n^2=(\frac{n}{{\sigma }^2}+\frac{1}{b}{)}^{-1}$.

${\sigma }^2$的条件（后验）分布：
p ( σ 2 ∣ X , μ ) ∝ π ( X ∣ μ , σ 2 ) π ( σ 2 ) ∝ σ − n exp ⁡ { − ∑ i ( X i − μ ) 2 2 σ 2 } σ − 2 ( c + 1 ) exp ⁡ { − d σ 2 } ∝ σ − 2 ( n / 2 + c + 1 ) exp ⁡ { − 1 σ 2 [ ∑ i ( X i − μ ) 2 2 + d ] } ∼ i n v G a m m a ( c n , d n ) p(\sigma^2|X,\mu) \propto \pi(X|\mu,\sigma^2)\pi(\sigma^2) \\ \propto \sigma^{-n}\exp\{ - \sum_i \frac{(X_i-\mu)^2}{2\sigma^2} \} \sigma^{-2(c+1)}\exp \{ -\frac{d}{\sigma^2}\} \\ \propto \sigma^{-2(n/2+c+1)} \exp\{-\frac{1}{\sigma^2}[\sum_i \frac{(X_i-\mu)^2}{2}+d ] \} \sim invGamma(c_n,d_n) p(σ2∣X,μ)∝π(X∣μ,σ2)π(σ2)∝σ−nexp{−i∑​2σ2(Xi​−μ)2​}σ−2(c+1)exp{−σ2d​}∝σ−2(n/2+c+1)exp{−σ21​[i∑​2(Xi​−μ)2​+d]}∼invGamma(cn​,dn​)
$c_n=n/2+c$, $d_n={\sum }_i\frac{(X_i-\mu {)}^2}{2}+d$

Gibbs采样R代码实现

$\mu$的更新函数

#给定观测值和其他参数，更新参数mu
updt_mu_r = function(X,sigma2,prior){  #输入X，上一步迭代中的sigma，超参数a,b的先验值
  n = length(X)
  sumX = sum(X)   
  tmp = 1 / (n/sigma2 + 1/prior[2])  #mu的后验分布方差
  mu_n = ( sumX/sigma2 + prior[1]/prior[2] ) * tmp #mu的后验分布均值
  sigma_n = sqrt(tmp) #mu的后验分布标准差
  return( rnorm(1,mu_n,sigma_n) ) #返回一个从mu的后验分布中的抽样值即为mu的更新值
}

${\sigma }^2$的更新函数

updt_sigma2_r = function(X,mu,prior){ #输入X，上一步迭代中的mu，超参数c,d的先验值
  n = length(X)
  shape = 0.5*n + prior[1]  #sigma^2的后验分布形状参数
  scale = sum((X-mu)^2)    
  scale = scale/2 + prior[2] #sigma^2的后验分布尺度参数
  return( 1/rgamma(1,shape,scale) )  #返回一个从sigma^2的后验分布中的抽样值即为sigma^2的更新值
}

主函数

mcmc_Gaussian_r = function(X, n_sample, burn_in, thin, init, prior_mu, prior_sigma){
#输入X，采样的量，燃烧期，保存的间隔数，参数初始值，mu的先验分布中的超参数值，sigma2的先验分布中的超参数值。
  
  mu_sample = rep(0,n_sample)
  sigma2_sample = rep(0,n_sample)
  #赋予初值
  mu = init[1]  
  sigma2 = init[2]
  #燃烧期
  for(i in 1:burn_in){     #反复更新采样后验均值和方差
    mu = updt_mu_r(X, sigma2, prior_mu)
    sigma2 = updt_sigma2_r(X,mu,prior_sigma)
  }
  
  if(thin<=0) thin=1  #thin不能小于1，表示每隔thin-1次迭代记录一次数据。
  
  for(i in 1:n_sample){
    for(j in 1:thin){
      mu = updt_mu_r(X, sigma2, prior_mu)
      sigma2 = updt_sigma2_r(X,mu,prior_sigma)
    }
    mu_sample[i] = mu
    sigma2_sample[i] = sigma2
  }
  
  out = list(mu =  mu_sample, sigma2 = sigma2_sample)
  return(out)  #返回列表，out[[1]]是mu的模拟值，out[[2]]是sigma2的模拟值
}

案例求解

#模拟从一个均值为1，方差为4的正态分布中抽样
a = rnorm(5000,1,2)  
res2 = mcmc_Gaussian_r(a,n_sample = 1000, burn_in = 100, thin = 1,
                       init = c(10,10), prior_mu = c(0,100),prior_sigma = c(1,1))
#返回结果res为一个列表                       

部分结果展示：

> res2[[1]][1:10]
 [1] 0.9787088 0.9541495 1.0188802 0.9628411 1.0027464 0.9689442 0.9722291 0.9486590
 [9] 0.9947923 0.9592428
> res2[[2]][1:10]
 [1] 3.840518 3.901282 3.968341 3.875591 3.974695 3.921036 3.888012 4.023784 3.890451
[10] 3.919828

很接近真实的均值和方差。

（补充）RCPP代码

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
void updt_mu(double &mu, NumericVector X, double sigma2, NumericVector prior){
  int n = X.size();
  double sumX = 0;
  double mu_n = 0;
  double sigma_n = 0;
  double tmp = 0;
  for(int i = 0; i < n; i++){
    sumX += X[i];
  }
  
  tmp = 1 / (n/sigma2 + 1/prior[1]);
  mu_n = ( sumX/sigma2 + prior[0]/prior[1] ) * tmp;
  sigma_n = sqrt(tmp);
  mu = rnorm(1,mu_n,sigma_n)[0]; 
}

// [[Rcpp::export]]
void updt_sigma2(double &sigma2, NumericVector X, double mu, NumericVector prior){
  int n = X.size();
  double shape = 0; 
  double scale = 0;
  shape = 0.5*n + prior[0];
  
  for(int i = 0; i < n; i++){
    scale += (X[i] - mu) * (X[i] - mu);
  }
  
  scale = scale/2 + prior[1];
  sigma2 = 1/rgamma(1,shape,1/scale)[0] ; 
}

// [[Rcpp::export]]
double loglik(NumericVector X, double mu, double sigma2){
  
  int n = X.size();
  double out = 0;
  
  for(int i=0; i < n; i++){
    out -= (X[i] - mu)*(X[i] - mu);
  }
  
  out = out/2/sigma2;
  out -=-0.5*n * log(2*3.1415926*sigma2);
  return out;
}

// [[Rcpp::export]]
List mcmc_Gaussian(NumericVector X, int n_sample, int burn_in, int thin, 
                   NumericVector init, NumericVector prior_mu, NumericVector prior_sigma){
  
  NumericVector mu_sample(n_sample);
  NumericVector sigma2_sample(n_sample);
  NumericVector loglk(n_sample);
  
  double mu = init[0];
  double sigma2 = init[1];
  
  for(int i = 0; i < burn_in; i++){
    updt_mu(mu, X, sigma2, prior_mu);
    updt_sigma2(sigma2, X, mu, prior_sigma);
  }
  
  if(thin<=0) thin=1;
  
  for(int i = 0; i < n_sample; i++){
/*
    if(i%100==0){
      std::cout << "iteration: " << i << " mu: " << mu << std::endl;
    }
*/
    
    for(int j = 0; j < thin; j++){
      updt_mu(mu, X, sigma2, prior_mu);
      updt_sigma2(sigma2, X, mu, prior_sigma);
    }
    
    mu_sample[i] = mu;
    sigma2_sample[i] = sigma2;
    loglk[i] =  loglik(X, mu, sigma2);
  }
  
  return List::create(
    _["mu"] = mu_sample,
    _["sigma2"] = sigma2_sample,
    _["loglik"] = loglk
  );
}

MH采样代码

MH采样原理不细说。
相同的案例。
直接上代码，来对比一下两种采样方法。
不好意思，这里只上RCPP代码，RCPP可以找资料学，不想学也不影响从代码中观看算法的逻辑！因为我也是RCPP入门级别选手~

小tips:老师说凡是遇到求密度，尤其是联合密度值，在代码中一定一定要取对数运算（不是对数变换，就是取对数后按照对数的计算规则进行后续一切运算），原因是密度在0~1之间，连续相乘后会非常接近0，很可能导致数值溢出啥的。

// [[Rcpp::export]]
//求X样本的联合概率密度的对数的函数，输入X,mu,sigma2，输出
double loglik(NumericVector X, double mu, double sigma2){
  int n = X.size();  //初始化变量n为数据X的维度大小
  double out = 0;  //初始化变量out=0

  for(int i=0; i < n; i++){              //i从0到n-1的循环
    out -= (X[i] - mu)*(X[i] - mu);      //这里在求X的联合密度的exp()内的一部分
  }

  out = out/2/sigma2;       //这里输出X的联合密度的exp()内的完整部分
  out -= 0.5*n * log(2*3.1415926*sigma2); //这里右边的式子计算的是X密度的exp()外面的数取对数，整个式子输出完整的X对数密度值
  return out;
}


// [[Rcpp::export]]
//MH更新参数的函数，注意到MH一次性更新所有参数，而Gibbs是逐个参数进行更新。
//该函数输入上一步迭代中的mu,sigma2，数据X，学习率lr，实际更新次数count。

void updt_para_MH(double &mu, double &sigma2, NumericVector X, double lr, double &count){
  //同时更新mu和sigma2，lr为学习率。更新步长采样自正态随机数（乘上学习率）
  NumericVector z0 = rnorm(2) * lr; 
  double mu_new = mu + z0[0];
  double sigma2_new = sigma2 + z0[1];

  double acc = 0;

  if(sigma2_new > 0){  //必须保证采样空间是正确的：sigma2为正数
    NumericVector test = runif(1);  //用来进行接受拒绝采样的随机数
    double logL_new = loglik(X,mu_new,sigma2_new); //计算新的对数似然
    double logL_old = loglik(X,mu,sigma2); //上一步参数对应的对数似然

    if( log(test[0]) < (logL_new - logL_old) ){ //如果均匀分布随机数小于新旧密度之比，则接受采样值
      acc = 1;   //确定更新
      count++;  //count记录真实的更新次数
    }
  }

  if(acc==1){     //确定更新
    mu = mu_new;
    sigma2 = sigma2_new;
  }
}


// [[Rcpp::export]]
//主函数：输入X，采样的次数，燃烧期，采样间隔记录的步数，参数初始值，学习率；输出参数MH采样值
List mcmc_Gaussian_MH(NumericVector X, int n_sample, int burn_in, int thin,
                      NumericVector init, double lr){
  //定义以及初始化变量               
  NumericVector mu_sample(n_sample);
  NumericVector sigma2_sample(n_sample);
  NumericVector loglk(n_sample);
  double mu = init[0];
  double sigma2 = init[1];
  double count = 0;              //这里初始化count=0
 //利用上一个参数更新函数进行参数的更新
  for(int i = 0; i < burn_in; i++){
    //updt_para_MH(mu,sigma2,X,lr,);   //偷懒省去燃烧期，更新过程与下面类似~
  }

  if(thin<=0) thin=1;       //thin必须为正整数

  for(int i = 0; i < n_sample; i++){   //循环n_sample次
    for(int j = 0; j < thin; j++){
    updt_para_MH(mu,sigma2,X,lr,count);   //更新参数！
    if(i%100==0){             #i每100进行下述操作
      if(count > 30){          #如果实际更新的次数大于30次进行下述操作
        lr = lr * 1.1;             #学习率增大0.1倍
      }else if(count < 20){  #否则学习率缩小0.1倍
        lr = lr * 0.9;
      }
      //上面的操作实际上是随着采样过程的进行，每100次采样中，当参数实际更新次数<20次时，学习率就衰减为原来的0.9倍；实际更新次数在20~30时，不改变学习率；实际更新次数>30时，学习率增大为原来的1.1倍.原因可能是怕刚开始走太快偏离方向于是慢慢走，后来趋于稳定变化很慢于是加大步伐？不懂这些实操技术，有些玄妙~
      std::cout << count/100 << " lr: " << lr << std::endl;  //不用管，用来在控制台输出学习率的变化情况~
      count = 0;    //count重置为0（该变量的时候加了"&"有关，像是内置参数）
     }
    }
    mu_sample[i] = mu;                 //记录  
    sigma2_sample[i] = sigma2;
    loglk[i] =  loglik(X, mu, sigma2);
  }
  //定义输出列表
  return List::create(
    _["mu"] = mu_sample,
    _["sigma2"] = sigma2_sample,
    _["loglik"] = loglk
  );
}

总结：

当参数的后验分布有显式的表达时，MCMC抽样不是必要的，在进行推断的时候有可能根据分布直接进行计算(如Gibbs)；当后验分布没有显式表达时（归一化参数在很多时候是算不出来的），可以通过MCMC算法进行抽样。