理解MCMC

在贝叶斯统计中，经常需要计算后验概率，概率计算就涉及到积分问题。一种解决方法是用解析式得到后验概率直接计算，另一种是利用统计模拟来计算近似值。 考虑一个简单问题，我们对一个硬币反复投掷，对于出现正面的概率theta先主观设定为一个均匀分布，然后实际投掷14次，得到11次正面，要根据这个信息data来更新后验概率。下面用统计模拟来计算。

	myData=c(1,1,1,1,1,1,1,1,1,1,1,0,0,0)
	# 尝试1万个不同的参数
	tryn = 1e4
	Theta = sort(runif(tryn))
	pTheta = 1/tryn
	z = sum( myData==1 )
	N = length( myData )
	# 似然函数
	pDataGivenTheta = Theta^z * (1-Theta)^(N-z)
	pData = sum( pDataGivenTheta * pTheta )
	# 后验概率
	pThetaGivenData = pDataGivenTheta * pTheta / pData
	plot(x=Theta,y=pThetaGivenData,type='l')

view raw simulation1.R hosted with ❤ by GitHub

因为实际的数据表现是正面出现次数较多，所以后验概率会做出相应调整。这个统计模拟的例子非常简单，容易实施。但如果涉及的参数个数很多时，计算量就会非常大。此时就可以尝试另一种统计模拟方法，即MCMC(Markov chain Monte Carlo)。我们先考虑一个思想实验：

一个岛国下面有7个岛，各岛人数不一样。旅行者已经在其中一个岛上，他要去各岛游历，准备在人口多的岛上游玩的时间长一些。但他并不知道具体的人口数，可以询问本岛和相邻岛的市长。他的旅行计划如下：
首先每天会扔一枚硬币，如果正面向上，就计划去左岛，如果反面向上，就计划去右岛。
然后比较当前所在岛人口数a和计划去的目的岛人口数b，若b小于a，则实际去的概率为b/a，若b大于a，则实际去的概率为1。所以实际去的概率归纳为min(1,b/a)。
就这样，旅行者每天不断的在不同岛间移动，最终他在各岛上呆的时间之比，就会收敛为各岛人口之比，这样就通过一个随机模拟得到了一个概率分布。

	# 假设七个岛人口数分布如下
	# 1:7/sum(1:7)
	# 尝试10万次旅行
	trajLength = 1e5
	trajectory = rep( 0 , trajLength )
	trajectory[1] = 3 # 第一次在第3个岛上
	target = function(currentPosition,proposedJump){
	tartetposition = currentPosition+proposedJump
	p = min(1,tartetposition/currentPosition)
	return(p)
	}
	for(t in 1:(trajLength-1) ) {
	currentPosition = trajectory[t]
	proposedJump = sample(c(-1,1),1)
	probAccept = target(currentPosition,proposedJump)
	if ( runif(1) < probAccept ) {
	trajectory[ t+1 ] = currentPosition + proposedJump
	if ( trajectory[ t+1 ]>7 | trajectory[ t+1 ]<1) {
	trajectory[ t+1 ] = currentPosition }
	} else {
	trajectory[ t+1 ] = currentPosition
	}
	}
	# 通过MCMC得到人口分布
	res = table(trajectory)
	res/sum(res)
	# 下面我们用MCMC来解决开始的硬币问题
	# data
	myData=c(1,1,1,1,1,1,1,1,1,1,1,0,0,0)
	# P(theta)
	prior = function( theta ) {
	prior = runif(length(theta))
	prior[theta>1|theta<0]=0
	return( prior )
	}
	# P(data|theta)
	likelihood = function( theta , data ) {
	z=sum(data==1)
	N = length( data )
	pDataGivenTheta = theta^z * (1-theta)^(N-z)
	pDataGivenTheta[ theta > 1 | theta < 0 ] = 0
	return( pDataGivenTheta )
	}
	# P(theata|data)*P(theta)
	targetRelProb = function( theta , data ) {
	targetRelProb = likelihood( theta , data ) * prior( theta )
	return( targetRelProb )
	}
	trajLength = 1e5
	trajectory = rep( 0 , trajLength )
	trajectory[1] = 0.50
	burnIn = ceiling( .1 * trajLength )
	nAccepted = 0
	nRejected = 0
	set.seed(47405)
	for(t in 1:(trajLength-1) ) {
	currentPosition = trajectory[t]
	proposedJump=rnorm(1,mean=0,sd=0.1)
	probAccept = min( 1,
	targetRelProb( currentPosition + proposedJump , myData )
	/ targetRelProb( currentPosition , myData ) )
	if ( runif(1) < probAccept ) {
	trajectory[ t+1 ] = currentPosition + proposedJump
	if ( t > burnIn ) { nAccepted = nAccepted + 1 }
	} else {
	trajectory[ t+1 ] = currentPosition
	if ( t > burnIn ) { nRejected = nRejected + 1 }
	}
	}
	acceptedTraj = trajectory[ (burnIn+1) : length(trajectory) ]
	hist(acceptedTraj)

view raw MCMC.R hosted with ❤ by GitHub

参考资料：Doing Bayesian Data Analysis with R