MCMC.R

# 假设七个岛人口数分布如下
# 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)