BSC贝叶斯-MCMC

library(mvtnorm)
set.seed(123)

主要函数

LogPosterior (N(0,I) , $X^2(10)$ )

LogPosterior <- function(beta, sigma2, y, X){
  p <- length(beta)
  if (sigma2<=0){
    return(NA)
  }
  else{
    ## The log likelihood function 取对数log,加和sum
    LogLikelihood <- sum(dnorm(x = y, mean = X%*%beta, 
                               sd = sqrt(sigma2), log = TRUE))
    
    ## The priors -beta-概率密度
    LogPrior4beta <- dmvnorm(x = t(beta), mean = matrix(0, 1, p),
                             sigma = diag(p), log = TRUE)
    # dmvnorm计算多元正态密度函数
    
    ## The priors -sigma2-概率密度
    LogPrior4sigma2 <- dchisq(x = sigma2, df = 10, log = TRUE)
    LogPrior <- LogPrior4beta + LogPrior4sigma2
    
    ## The log posterior
    LogPosterior <- LogLikelihood + LogPrior
    return(LogPosterior)
  }
}

这里中间加了一个判断，因为是随机抽出的sigma2（来自正太）有一定可能是负数，所以加了一个if和else判断。否则会持续报错NANANANANANA。

Gibbs函数

Gibbs<-function(X,y,beta.ini,sigma2.ini,beta.step,sigma2.step,nIter){
  ## 储存sample
  p <- ncol(X) #维度
  beta <- matrix(NA, nIter, p)
  sigma2 <- matrix(NA, nIter, 1)
  acc <- matrix(NA, nIter, 2)
  #接受概率，2列分别是beta和sigma2
  beta[1, ] <- beta.ini
  sigma2[1, ] <- sigma2.ini
  
  for(i in 2:nIter){
    beta.current <- beta[i-1, ]
    sigma2.current <- sigma2[i-1,]
    
    ## The Metropolis Algorithm For beta
    beta.prop <- rmvnorm(n = 1,
      mean = matrix(beta.current, 1),
      sigma = diag(p)*beta.step)# FV
    #beta.prop <- rmvnorm(n = 1, mean = matrix(beta.current, 1), sigma = diag(p)) 
    
    logPosterior.beta.prop <- LogPosterior(
      beta = t(beta.prop),
      sigma2 = sigma2.current,
      y = y, X = X)
    
    logPosterior.beta.current <- LogPosterior(
      beta = beta.current,
      sigma2 = sigma2.current,
      y = y, X = X)
    
    logratio <- (logPosterior.beta.prop -
                   logPosterior.beta.current)
    
    acc.prob.beta <- min(exp(logratio), 1) #做差取指数，就是概率之比
    
    if(!is.na(acc.prob.beta) & runif(1) < acc.prob.beta)
    {
      beta.current <- beta.prop
    }
    acc[i, 1] <- acc.prob.beta
    beta[i, ] <- beta.current
    
    
    ## The Metropolis Algorithm For sigma2
    # sigma2.prop <- runif(1,0,2*sigma2.current)
    sigma2.prop <- rnorm(1, sigma2.current,sigma2.step) 
    
    logPosterior.sigma2.prop <- LogPosterior(
      beta = matrix(beta.current),
      sigma2 = sigma2.prop,
      y = y, X = X)
    
    logPosterior.sigma2.current <- LogPosterior(
      beta = matrix(beta.current),
      sigma2 = sigma2.current,
      y = y, X = X)
    
    logratio <- logPosterior.sigma2.prop -
      logPosterior.sigma2.current
    
    acc.prob.sigma2 <- min(exp(logratio), 1)
    
    if(!is.na(acc.prob.sigma2) & runif(1, 0, 1) < acc.prob.sigma2) # accept the proposal
    {
      sigma2.current <- sigma2.prop
    }
    
    ## Update the matrix
    acc[i, 2] <- acc.prob.sigma2
    sigma2[i, ] <- sigma2.current
  }
  return(list(beta=beta,acc=acc,sigma2=sigma2))
}

有关Gibbs：其中修改过sigma2.prop <- runif(1,0,2*sigma2.current)，但是可能是因为没有体现在old附近产生新样本，多次都无法收敛，所以用正太随机更适宜。

数值模拟

# real: beta -2 3 -1 2  sigma2 1
beta <- matrix(c(-2, 3, -1, 2))
X <- cbind(1, matrix(runif(1000*3), 1000))
y <- X%*%beta + rnorm(1000,0,1)
ou<-Gibbs(X,y,beta.ini = c(0,2,2,1),sigma2.ini = 0.5,
          beta.step = 0.01,sigma2.step=0.1,nIter = 10000)
# View(ou$beta)
plot(ou$beta[, 1], type = "l", lwd = 2, col = 6, ylim = c(-4, 10))
lines(ou$beta[, 2], type = "l", lwd = 2, col = 2)
lines(ou$beta[, 3], type = "l", lwd = 2, col = 3)
lines(ou$beta[, 4], type = "l", lwd = 2, col = 5)
abline(h = -2, lwd = 2, col = 1, lty = 3)
abline(h = 3, lwd = 2, col = 1, lty = 3)
abline(h = -1, lwd = 2, col = 1, lty = 3)
abline(h = 2, lwd = 2, col = 1, lty = 3)
lines(ou$sigma2, type = "l", lwd = 2, col ="yellow2", ylim = c(-4, 10))
abline(h=1,lty = 3)
legend("topright",c("beta1","beta2","beta3","beta4","sigma2","real")
       ,col = c(6,2,3,5,"yellow2",1),lty = c(1,1,1,1,1,3)
       ,ncol=2)

数值模拟的真值为beta = -2 3 -1 2 sigma2=1，可以看出结果不错。千步以后就收敛了。

中间2个关键的参数beta.step = 0.01,sigma2.step=0.1，是通过反复试验得到的，最开始设置过大，一直不接受，始终在初始值。

Bodyfat

OLS结果

# bodyfat
data<-read.csv("Bodyfat.txt")
OLS.bodyfat<-lm(bodyfat~Density,data=data)
beta.ols=OLS.bodyfat$coefficients
sigma2.ols=summary(OLS.bodyfat)$sigma**2

先是做了OLS，发现Density的t检验通过，所以选择了二者进行回归。

beta.ols和sigma2.ols在后续调整先验的时候用到。

修改先验

## 修改先验函数
LogPosterior <- function(beta, sigma2, y, X){
  p <- length(beta)
  if (sigma2<=0){
    return(NA)
  }
  else{
    ## The log likelihood function 取对数log,加和sum
    LogLikelihood <- sum(dnorm(x = y, mean = X%*%beta, 
                               sd = sqrt(sigma2), log = TRUE))
    
    ## The priors -beta-概率密度
    LogPrior4beta <- dmvnorm(x = t(beta), mean =matrix(beta.ols,1,2),
                             sigma = diag(p), log = TRUE)
    # dmvnorm计算多元正态密度函数
    
    ## The priors -sigma2-概率密度
    LogPrior4sigma2 <- dchisq(x = sigma2, df = sigma2.ols, log = TRUE)
    LogPrior <- LogPrior4beta + LogPrior4sigma2
    
    ## The log posterior
    LogPosterior <- LogLikelihood + LogPrior
    return(LogPosterior)
  }
}

LogPrior4beta和LogPrior4sigma2进行了修改。均值改为OLS的beta，sigma2的自由度也改为了sigma2.ols

Bodyfat模拟

X <- cbind(1, matrix(data$Density))
y <- as.matrix(data$bodyfat)
ou<-Gibbs(X,y,beta.ini = c(0,2),sigma2.ini = 0.5,beta.step = 100,sigma2.step =100,nIter = 10000)

# View(ou$beta)
plot(ou$beta[, 1], type = "l", lwd = 2, col = 6,ylim = c(-500,500))
lines(ou$beta[, 2], type = "l", lwd = 2, col = 2)
abline(h = 477.650, lwd = 2, col = 1, lty = 3)
abline(h =-434.360, lwd = 2, col = 1, lty = 3)
lines(ou$sigma2, type = "l", lwd = 2, col ="yellow2", ylim = c(-4, 10))
abline(h=1.707688,lty = 3)
legend("topright",c("beta1","beta2","sigma2","real")
       ,col = c(6,2,"yellow2",1),lty = c(1,1,1,3)
       ,ncol=2)

模型的$\beta$均收敛至真实值（虚线）。sigma2有一段大幅跳跃。

plot(data$Density,data$bodyfat)
for (i in 500:10000){
    abline(ou$beta[i,1],ou$beta[i,2],col=3,lwd=0.1)
}
abline(OLS.bodyfat$coefficients[1],OLS.bodyfat$coefficients[2],col=2)
legend("topright",c("样本点","Gibbs","OLS"),col=c("black",3,2),pch=c(1,NA,NA),lty=c(NA,1,1))

值得注意的是我去除了前1000个Gibbs的sample，上图可以看到OLS和Gibbs重合效果很好。

反思和改进

1. 绘图的时候应该去掉前面的样本，噪声

2. 先验和target越近越好，所以bodyfat这个数据集要修改先验

3. beta.step 和 sigma2.step 要多吃尝试，过大会波动剧烈，过小会很慢很慢。

4. 高阶图例这样画 pch=c(1,NA,NA),lty=c(NA,1,1) 巧妙利用NA可以散点图和折线图叠加