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

?

在邏輯回歸中，我們將二元響應(yīng)\（Y_i \）回歸到協(xié)變量\（X_i \）上。下面的代碼使用Metropolis采樣來(lái)探索\（\ beta_1 \）和\（\ beta_2 \）的后驗(yàn)YiYi到協(xié)變量XiXi。

?

定義expit和分對(duì)數(shù)鏈接函數(shù)

1. logit<-function(x){log(x/(1-x))} 此函數(shù)計(jì)算\（（\ beta_1，\ beta_2）\）的聯(lián)合后驗(yàn)。它返回后驗(yàn)的對(duì)數(shù)以獲得數(shù)值穩(wěn)定性。(β1,β2)(β1,β2)。它返回后驗(yàn)的對(duì)數(shù)以獲得數(shù)值穩(wěn)定性。

2. log_post<-function(Y,X,beta){

3. prob1 ?<- expit(beta[1] + beta[2]*X)

4. like+prior}

這是MCMC的主要功能.can.sd是候選標(biāo)準(zhǔn)偏差。

1. Bayes.logistic<-function(y,X,

2. n.samples=10000,

3. can.sd=0.1){

4. 


5. keep.beta ? ? <- matrix(0,n.samples,2)

6. keep.beta[1,] <- beta

7. 


8. acc ? <- att <- rep(0,2)

9. 


10. for(i in 2:n.samples){

11. 


12. for(j in 1:2){

13. 


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

15. 


16. # Draw candidate:

17. 


18. canbeta ? ?<- beta

19. canbeta[j] <- rnorm(1,beta[j],can.sd)

20. canlp ? ? ?<- log_post(Y,X,canbeta)

21. 


22. # Compute acceptance ratio:

23. 


24. R <- exp(canlp-curlp)

25. U <- runif(1)

26. if(U<R){

27. acc[j] <- acc[j]+1

28. }

29. }

30. keep.beta[i,]<-beta

31. 


32. }

33. # Return the posterior samples of beta and

34. # the Metropolis acceptance rates

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

生成一些模擬數(shù)據(jù)

1. set.seed(2008)

2. n ? ? ? ? <- 100

3. X ? ? ? ? <- rnorm(n)

4. true.p ? ?<- expit(true.beta[1]+true.beta[2]*X)

5. Y ? ? ? ? <- rbinom(n,1,true.p)

擬合模型

1. burn ? ? ?<- 10000

2. n.samples <- 50000

3. fit ?<- Bayes.logistic(Y,X,n.samples=n.samples,can.sd=0.5)

4. tock <- proc.time()[3]

5. 


6. tock-tick

結(jié)果

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

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

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

2. 


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

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

1. print("Posterior mean/sd")

2. ## [1] "Posterior mean/sd"

3. print(round(apply(fit$beta[burn:n.samples,],2,mean),3))

4. ## [1] -0.076 ?0.798

5. print(round(apply(fit$beta[burn:n.samples,],2,sd),3))

6. ## [1] 0.214 0.268

7. print(summary(glm(Y~X,family="binomial")))

8. ##

9. ## Call:

10. ## glm(formula = Y ~ X, family = "binomial")

11. ##

12. ## Deviance Residuals:

13. ## ? ? Min ? ? ? 1Q ? Median ? ? ? 3Q ? ? ?Max

14. ## -1.6990 ?-1.1039 ?-0.6138 ? 1.0955 ? 1.8275

15. ##

16. ## Coefficients:

17. ## ? ? ? ? ? ? Estimate Std. Error z value Pr(>|z|)

18. ## (Intercept) -0.07393 ? ?0.21034 ?-0.352 ?0.72521

19. ## X ? ? ? ? ? ?0.76807 ? ?0.26370 ? 2.913 ?0.00358 **

20. ## ---

21. ## Signif. codes: ?0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

22. ##

23. ## (Dispersion parameter for binomial family taken to be 1)

24. ##

25. ## ? ? Null deviance: 138.47 ?on 99 ?degrees of freedom

26. ## Residual deviance: 128.57 ?on 98 ?degrees of freedom

27. ## AIC: 132.57

28. ##

29. ## Number of Fisher Scoring iterations: 4

?

非常感謝您閱讀本文，有任何問(wèn)題請(qǐng)?jiān)谙旅媪粞裕?/h1>

標(biāo)簽：