###############################################################
#  以下MH算法，是通过LearnBayes包中的函数实现的，所以首先加载相关包
###############################################################
library(LearnBayes)
library(MASS)
library(mvtnorm)

logpost <- function(x, par) { #定义目标分布函数的对数形式
  mean <- par$mean
  sigma <- par$sigma
  return(dmvnorm(x, mean, sigma, log = T))
}
mean <- c(1, 1)
sigma <- matrix(c(1, 0.9, 0.9, 1), nrow = 2,
   ncol = 2)
par <- list(mean = mean, sigma = sigma)  #建议分布相关参数

varcov <- diag(c(0.3, 0.3))
proposal <- list(var = varcov, scale = 2)
start <- c(-4, 8)
m = 5000
set.seed(123)
rwmetrop1 <- function (logpost, proposal, start, m, ...) #对rwmetrop函数稍微修改
{
    pb = length(start)
    MM <- array(0, c(m, 2 * pb))
    Mpar = array(0, c(m, pb))
    b = matrix(t(start))
    lb = logpost(start, ...)
    a = chol(proposal$var)
    scale = proposal$scale
    accept = 0
    for (i in 1:m) {
        bc = b + scale * t(a) %*% array(rnorm(pb), c(pb, 1))
        
        lbc = logpost(t(bc), ...)
        prob = exp(lbc - lb)
        if (is.na(prob) == FALSE) {
            if (runif(1) < prob) {
                lb = lbc
                b = bc
                accept = accept + 1
            } else {
              MM[i, ] <- c(b, bc)
            }
        }
        Mpar[i, ] = b
    }
    accept = accept/m
    stuff = list(par = Mpar, accept = accept, MM = MM)
    return(stuff)
}

s <- rwmetrop1(logpost, proposal, start, m, par)
s1 <- rwmetrop1(logpost, proposal, start=c(5, -10), m, par)  #改变初始位置

###############################################################
#  用animation包中的函数saveHTML将以上MH算法动态化
###############################################################

saveHTML({
     ani.options(interval = 0.3)
   for(i in 1:200) {
  if(i <= 20) {
  mycontour(logpost, c(-10,10,-10,10),par, 
    xlim = c(-10, 8), ylim = c(-10, 8), nlevels=10, col = "red")
  if(i<=2) {

  points(s$par[1:i, 1],s$par[1:i, 2] ,col = "blue", cex = 2) 
  } else {
      lines(s$par[1:i, ], col = "blue", lwd = 2)}
  }
  else {
    mycontour(logpost, c(-10,10,-10,10),par, 
    xlim = c(-5, 6), ylim = c(-5, 6), nlevels=10, col = "red")
    lines(s$par[1:i, ], col = "blue", lwd = 2)}
  }},
    img.name = "brownian_motion_c", htmlfile = "index5.html", single.opts = "'controls': ['first', 'previous', 'play', 'next', 'last', 'loop', 'speed'], 'delayMin': 0", 
    description = c("Random walk of 10 points on the 2D plane", "(without the navigation panel)")
)


###############################################################
#  不论初始值设置在哪里，经过有限的几步后，都迅速转移到密度较高的中间区域
###############################################################
mycontour(logpost, c(-10,10,-10,10),par, 
    xlim = c(-8, 10), ylim = c(-8, 10), nlevels=10)
lines(s[[1]], col = rgb(0, 0, 0, 0.5))
mycontour(logpost, c(-10,10,-10,10),par, 
    xlim = c(-8, 10), ylim = c(-8, 10), nlevels=10)
lines(s1[[1]], col = rgb(0, 0, 0, 0.5))



###############################################################
#将接受的建议与拒绝建议分别显示，，建议分布产生的转移建议就像脱缰的“野马”，四
#处乱窜，而接受概率就像一条绳索，时时牵拽着这头野马，让它在目标分布所形成的“水草丰美”
#的区域漫步。
###############################################################
mycontour1 <- function (logf, limits, data, ...) 
{
    LOGF = function(theta, data) {
        if (is.matrix(theta) == TRUE) {
            val = matrix(0, c(dim(theta)[1], 1))
            for (j in 1:dim(theta)[1]) val[j] = logf(theta[j, 
                ], data)
        }
        else val = logf(theta, data)
        return(val)
    }
    ng = 50
    x0 = seq(limits[1], limits[2], len = ng)
    y0 = seq(limits[3], limits[4], len = ng)
    X = outer(x0, rep(1, ng))
    Y = outer(rep(1, ng), y0)
    n2 = ng^2
    Z = LOGF(cbind(X[1:n2], Y[1:n2]), data)
    Z = Z - max(Z)
    Z = matrix(Z, c(ng, ng))
    contour(x0, y0, Z, levels = -6, lwd = 1, 
        ...)
}
plot(0, type = 'n', xlim = c(-4, 6), ylim = c(-4, 6), xlab = '', ylab = '', 
  axes = F)
box()
lines(s[[1]][100:250,], col = rgb(0, 0, 0, 0.5), lwd = 2)
legend("topleft", legend = c("接受建议", "拒绝建议"), lty = c(1, 3), lwd = c(2, 1), 
  col = rgb(0, 0, 0, 0.5), bty = 'n' )
mycontour1(logpost, c(-4,5,-4,5),par, 
    xlim = c(-4, 10), ylim = c(-4, 10), add = T)
MM <- (s$MM)[100:250, ]
segments(MM[, 1], MM[, 2], MM[, 3], MM[, 4], col = rgb(0, 0, 0, 0.5), lty = 3)


