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

题目背景

记Xi,i=1...,n{X_i,i=1...,n}Xi,i=1...,n来自一维正态N(μ,σ2)N(\mu,\sigma^2)N(μ,σ2)的独立同分布样本，μ\muμ具有先验分布N(a,b)N(a,b)N(a,b)，σ2\sigma^2σ2具有先验分布invGamma(c,d)invGamma(c,d)invGamma(c,d).
似然函数：π(X∣μ,σ2)∝(σ2)−n/2exp⁡{−∑i(Xi−μ)22σ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(α,λ)invGamma(\alpha,\lambda)invGamma(α,λ)的密度函数形式：
p(x;α,λ)=λαΓ(α)x−(α+1)exp⁡{−λx},x≥0p(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(Xi−μ)22σ2}exp⁡{−(μ−a)22b}∝exp⁡{−(μ−μn)22σn2}∼N(μn,σn2)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)
μn=(∑iXiσ2+ab)/(nσ2+1b)\mu_n=(\frac{\sum_i X_i}{\sigma^2}+\frac{a}{b})/(\frac{n}{\sigma^2} + \frac{1}{b})μn=(σ2∑iXi+ba)/(σ2n+b1),σn2=(nσ2+1b)−1\sigma_n^2=(\frac{n}{\sigma^2}+\frac{1}{b})^{-1}σn2=(σ2n+b1)−1.

σ2\sigma^2σ2的条件（后验）分布：
p(σ2∣X,μ)∝π(X∣μ,σ2)π(σ2)∝σ−nexp⁡{−∑i(Xi−μ)22σ2}σ−2(c+1)exp⁡{−dσ2}∝σ−2(n/2+c+1)exp⁡{−1σ2[∑i(Xi−μ)22+d]}∼invGamma(cn,dn)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)
cn=n/2+cc_n=n/2+ccn=n/2+c, dn=∑i(Xi−μ)22+dd_n=\sum_i \frac{(X_i-\mu)^2}{2}+ddn=∑i2(Xi−μ)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的更新值
}

σ2\sigma^2σ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] = musigma2_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=0for(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算法进行抽样。

MCMC_Gibbs/MH采样的简单案例R代码实现相关推荐

MCMC(二)：MCMC采样和M-H采样
在前面<MCMC(一):蒙特卡罗方法和马尔科夫链>中,我们讲到给定一个概率平稳分布 π \pi π, 很难直接找到对应的马尔科夫链状态转移矩阵 P P P.而只要解决这个问题,我们就可以找 ...
Android复习12【广播接收者-BroadcastReceiver（简单案例-发送广播、静态注册、动态注册、本地广播、代码示例(别处登陆踢用户下线)、常用系统广播总结、音乐播放器）】
2020-04-28[11周-周二] 音乐播放器Android代码下载:https://wws.lanzous.com/ifqzihaxvij 目录简单案例-发送广播 2)动态注册实例(监听网 ...
MCMC详解2——MCMC采样、M-H采样、Gibbs采样（附代码）
MCMC是一种随机采样方法,用来处理一些复杂运算的近似求解.在HMM.LDA等模型中都有重要应用. 上一篇 MCMC详解1--蒙特卡洛方法目录 1,马尔科夫链模型 1.1 马尔科夫链转移矩阵的性质 ...
独家 | 利用Auto ARIMA构建高性能时间序列模型（附Python和R代码）
作者:AISHWARYA SINGH 翻译:陈之炎校对:丁楠雅本文共3400字,建议阅读10+分钟. 本文介绍了ARIMA的概念,并带你用Python和R训练一个数据集实现它. 简介想象你现在有 ...
手把手教你用Prophet快速进行时间序列预测（附Prophet和R代码）
作者:ANKIT CHOUDHARY 翻译:王雨桐校对:丁楠雅本文约3000字,建议阅读12分钟. 本文将通过拆解Prophet的原理及代码实例来讲解如何运用Prophet进行时间序列预测. 简介 ...
机器学习算法清单！附Python和R代码
来源:数据与算法之美本文约6000字,建议阅读8分钟. 通过本文为大家介绍了3种机器学习算法方式以及10种机器学习算法的清单,学起来吧~ 前言谷歌董事长施密特曾说过:虽然谷歌的无人驾驶汽车和机器人 ...
10 种机器学习算法的要点（附 Python 和 R 代码）（转载）
10 种机器学习算法的要点(附 Python 和 R 代码)(转载) from:https://zhuanlan.zhihu.com/p/25273698 前言谷歌董事长施密特曾说过:虽然谷歌的无人 ...
ViewPager 实现页面左右滑动的简单案例1
ViewPager 实现页面左右滑动的简单案例主要Activity: <RelativeLayoutxmlns:android="http://schemas.android.com ...
10 种机器学习算法的要点（附 Python 和 R 代码）
前言谷歌董事长施密特曾说过:虽然谷歌的无人驾驶汽车和机器人受到了许多媒体关注,但是这家公司真正的未来在于机器学习,一种让计算机更聪明.更个性化的技术. 也许我们生活在人类历史上最关键的时期:从使用大 ...

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

文章目录