R语言使用Metropolis- Hasting抽样算法进行逻辑回归

原文链接：http://tecdat.cn/?p=6761

在逻辑回归中，我们将二元因变量Y_i回归到协变量X_i上。下面的代码使用Metropolis采样来探索 beta_1和beta_2 的后验Yi到协变量Xi。

定义expit和分对数链接函数

logit<-function(x){log(x/(1-x))} 此函数计算beta_1，beta_2的联合后验。它返回后验的对数以获得数值稳定性。(β1,β2)(β1,β2)。它返回后验的对数获得数值稳定性。
log_post<-function(Y,X,beta){
prob1 <- expit(beta[1] + beta[2]*X)
like+prior}

这是MCMC的主要功能.can.sd是候选标准偏差。

Bayes.logistic<-function(y,X,
n.samples=10000,
can.sd=0.1){

keep.beta <- matrix(0,n.samples,2)
keep.beta[1,] <- beta

acc <- att <- rep(0,2)

for(i in 2:n.samples){

for(j in 1:2){

att[j] <- att[j] + 1

# 抽取候选:

canbeta <- beta
canbeta[j] <- rnorm(1,beta[j],can.sd)
canlp <- log_post(Y,X,canbeta)

# 计算接受率:

R <- exp(canlp-curlp)
U <- runif(1)
if(U<R){
acc[j] <- acc[j]+1
}
}
keep.beta[i,]<-beta

}
# 返回beta的后验样本和Metropolis的接受率

list(beta=keep.beta,acc.rate=acc/att)}

生成模拟数据

set.seed(2008)
n <- 100
X <- rnorm(n)
true.p <- expit(true.beta[1]+true.beta[2]*X)
Y <- rbinom(n,1,true.p)

拟合模型

burn <- 10000
n.samples <- 50000
fit <- Bayes.logistic(Y,X,n.samples=n.samples,can.sd=0.5)
tock <- proc.time()[3]

tock-tick## 3.72

结果

abline(true.beta[1],0,lwd=2,col=2)

abline(true.beta[2],0,lwd=2,col=2)

hist(fit$beta[,1],main="Intercept",xlab=expression(beta[1]),breaks=50)

abline(v=true.beta[2],lwd=2,col=2)

print("Posterior mean/sd")
## [1] "Posterior mean/sd"
print(round(apply(fit$beta[burn:n.samples,],2,mean),3))
## [1] -0.076 0.798
print(round(apply(fit$beta[burn:n.samples,],2,sd),3))
## [1] 0.214 0.268

##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6990 -1.1039 -0.6138 1.0955 1.8275
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07393 0.21034 -0.352 0.72521
## X 0.76807 0.26370 2.913 0.00358 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 138.47 on 99 degrees of freedom
## Residual deviance: 128.57 on 98 degrees of freedom
## AIC: 132.57
##
## Number of Fisher Scoring iterations: 4

非常感谢您阅读本文，有任何问题请在下面留言！