library(mvtnorm)
set.seed(123)
主要函數
LogPosterior (N(0,I) , X 2 ( 10 ) X^2(10) X2(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重合效果很好。
反思和改進
- 繪圖的時候應該去掉前面的樣本，噪聲
- 先驗和target越近越好，是以bodyfat這個資料集要修改先驗
- beta.step 和 sigma2.step 要多吃嘗試，過大會波動劇烈，過小會很慢很慢。
- 高階圖例這樣畫 pch=c(1,NA,NA),lty=c(NA,1,1) 巧妙利用NA可以散點圖和折線圖疊加
![go語言的并發原理（goroutine）[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)