张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

MCMC中的Metropolis Hastings抽样法

已有 31796 次阅读 2012-2-29 07:54 |个人分类:统计分析|系统分类:科研笔记| MCMC, Hastings

Metropolis Hastings抽样法示例
jinlongzhang01@gmail.com
Metropolis Hasting(下面简称MH)是蒙特卡罗马尔科夫链中一种重要的抽样方法。本文简要介绍MH算法,并给出一个实例。
MH算法在参数空间随机取值,作为起始点。按照参数的概率分布生成随机的参数,按照这一系列参数的组合,计算当前点的概率密度。依据当前点和起始点概率密度比值是否大于(0,1)之间的随机数来判断是否保留当前点。若当前点的概率密度大于该随机数,就称这种状态为接受状态,此时,在满足参数概率分布的前提下,继续随机抽取参数的组合,作为下一点,计算下一点的概率密度,并计算下一点概率密度和概率密度的比值,并继续循环。若当前点不能被接受,则继续在满足参数概率分布的前提下,继续生成随机数,作为新的参数组合,直到参数组合能够被接受为止。
 
下面是用MCMC MH抽样法生成满足一定条件二元正态分布的例子,供感兴趣的同仁参考。
问题:## 对于一个二元正态分布,给定输入数据点向量x (x包含两个元素,x1,x2), 平均值参数向量mu(x1,x2),和sigma(方差矩阵),写出二元正态分布的概率密度函数。

该问题引自: http://users.aims.ac.za/~ioana/
由于本人对MCMC理解不深,对一些概念的理解难免出现错误,望读者能够批评得阅读,若发现问题,请及时告知本人。

## 解答:

mvdnorm <- function(x, mu, sigma){

    #x减去mu

    x.minus.mu <- x - mu 

    exp.arg    <- -0.5 * sum(x.minus.mu * solve(sigma, x.minus.mu))  

    # det(sigma) sigma 的行列式

    return( 1 / (2 * pi * sqrt(det(sigma))) * exp(exp.arg) )

}

 

## 问题二

## 假设二元正态分布的参数如下: 

## 两个维度的平均值分别为 2 3

# 协方差矩阵为

# 4 1

# 1 4

# 尝试用蒙特卡洛马尔科夫链 Metropolis Hastings 抽样法生成后验分布,进行10000次随机抽样,并计算随机点的接受率。

 

# 答: 按照题意,有

mu <- c(2 ,3)

sigma <- matrix(c(4, 1, 1, 4), nrow = 2)

 

# 限制sampler在空间的移动速率,数值越大,变化越快,该数值的设定待进一步讨论。

 

sd.proposal <- 2

 

## 设定模拟的次数

n <- 10000

## 生成NA组成的矩阵,用于保存模拟的结果

x <- matrix(nrow = n, ncol = 2)

 

# 设定sampler的初始值,假定数据点从 0, 0开始 (实际上该sampler可以从任意点开始移动)

cur.x <- c(0, 0)

 

# 计算给定初始值时的概率密度

cur.f <- mvdnorm(cur.x, mu, sigma)

 

### 蒙特卡洛马尔科夫链

n.accepted <- 0

for(i in 1:n){

    new.x <- cur.x + sd.proposal * rnorm(2)  ## 随机生成x

    new.f <- mvdnorm(new.x, mu, sigma)       ## 计算概率密度

    if(runif(1) < new.f/cur.f){  

        ## new.f/cur.f 概率密度的比率 (0,1)之间的随机数相比

        ## 若该比率小于随机数,则接受该点

       n.accepted <- n.accepted + 1

       cur.x <- new.x

       cur.f <- new.f

       }

    x[i,] <- cur.x                            ## cur.x存到第i

}

 

#查看接受率

n.accepted/n

 #查看每个变量的随机变化情况

par(mfrow=c(2,1))

plot(x[,1], type="l", xlab="t", ylab="X_1", main="Sample path of X_1")

plot(x[,2], type="l", xlab="t", ylab="X_2", main="Sample path of X_2")

 

## 绘制椭圆概率密度图

library(MASS)

proline.density <- kde2d(x[,1], x[,2], h = 5)

par(mfrow = c(1, 1))

plot(x, col = "gray")

contour(proline.density, add = TRUE)

points(2,3, pch = 3)


抽样后的概率密度
Trace plot


https://blog.sciencenet.cn/blog-255662-542389.html

上一篇:将比对好的fasta序列转换成relaxed phylip格式
下一篇:于默奥的蓝天白云
收藏 IP: 213.21.84.*| 热度|

5 陈学雷 黄继红 徐姗 符海洋 余晓美

该博文允许注册用户评论 请点击登录 评论 (5 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 00:55

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部