理解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') 

因为实际的数据表现是正面出现次数较多，所以后验概率会做出相应调整。这个统计模拟的例子非常简单，容易实施。但如果涉及的参数个数很多时，计算量就会非常大。此时就可以尝试另一种统计模拟方法，即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 ) { nAcceptednAccepted = nAccepted + 1 } 
  } else { 
    trajectory[ t+1 ] = currentPosition 
    if ( t > burnIn ) { nRejectednRejected = nRejected + 1 } 
  } 
} 
  
acceptedTraj = trajectory[ (burnIn+1) : length(trajectory) ] 
hist(acceptedTraj)