MCMC马尔可夫蒙特卡洛模拟R语言实例

##MCMC实例
##由两个正态分布根据相应权重构成的新分布
target<-function(x) 0.5*dnorm(x,-1,0.7)+0.5*dnorm(x,1,0.8)
curve(target,col='blue',lty=2,lwd=2,xlim=c(-4,4))
##设计一个标准差为0.5，变量x作为均值的正态分布抽样公式q
##为后续从该分布中抽取随机数做准备
##标准差可以尝试不同值来比较最终的抽样效果
q=function(x) rnorm(1,mean=x,sd=0.5)
##设置一个包含了10000个值得数据集合x
x<-rep(NA,10000)
##覆盖初始变量值x[1]为1
##初始值可以尝试不同值来比较最终的抽样效果
x[1]<-1
##循环语句：将当前得x[i]值作为变量x录入建议正态分布抽样公式q中
##得到新值xj套入Metropolis公式，计算接受概率A
##引入一个服从[0,1]均匀分布的随机值，与A比较，当A大时，用xj值覆盖x[i+1],否则x[i]覆盖x[i+1]
for(i in 1:9999){
  xj<-q(x[i])
  A<-min(target(xj)/target(x[i]),1)
  if (runif(1)<A)
    x[i+1]<-xj else x[i+1]<-x[i]
}
##绘制更新后的数据集合直方图
hist(x,breaks=100)

 

老师让学习了解MCMC，看了一个讲解的比较好的视频分享给大家：【数之道】马尔可夫链蒙特卡洛方法是什么？十五分钟理解这个数据科学难点_哔哩哔哩_bilibili

简而言之：MCMC中蒙特卡洛发挥接受采样作用，马尔科夫链是找到蒙特卡洛中最为稳态的建议函数。

 后续还会继续学习更新。