原文链接：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

## elapsed ##    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)

 hist(fit$beta\[,2\],main="Slope",xlab=expression(beta\[2\]),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

十分感谢您浏览本文，有任何问题请在上面留言！

最受欢迎的见解

1.应用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行

2.R语言中的Stan概率编程MCMC采样的贝叶斯模型

3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

4.R语言BUGS JAGS贝叶斯剖析 马尔科夫链蒙特卡洛办法（MCMC）采样

5.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

6.R语言Gibbs抽样的贝叶斯简略线性回归仿真剖析

7.R语言用Rcpp减速Metropolis-Hastings抽样预计贝叶斯逻辑回归模型的参数

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

9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长