[貝葉斯專著]Bayesian Biostatistics(using WinBUGS)

1. # computation posterior measures stroke study
   
2. 
   
3. alpha <- 19
   
4. beta <- 133
   
5. 
   
6. pmode <- (alpha-1)/(alpha+beta-2);pmode
   
7. 
   
8. pmean <- alpha/(alpha+beta);pmean
   
9. 
   
10. pmedian <- qbeta(0.5,alpha,beta);pmedian
    
11. 
    
12. pvar <- alpha*beta/((alpha+beta)^2*(alpha+beta+1));pvar
    
13. psd <- sqrt(pvar);psd

复制代码

1. # gamma prior density to STM study
   
2. # plot posterior density and normal approximating
   
3. # based on first 10 children
   
4. 
   
5. attach(dmft)
   
6. objects(2)
   
7. 
   
8. par(mfrow=c(1,1))
   
9. par(pty="m")
   
10. 
    
11. dmft <- V1[1:10]
    
12. mdmft <- mean(dmft)
    
13. 
    
14. n <- length(dmft)
    
15. sdmft <- sum(dmft)
    
16. 
    
17. n;sdmft
    
18. 
    
19. x <- seq(1,4.5,0.001)
    
20. 
    
21. alphap <- 3+sdmft
    
22. betap <- 1+n
    
23. 
    
24. y <- dgamma(x, shape=alphap, rate = betap)
    
25. 
    
26. plot(x,y,xlab=expression(theta),ylab="",type="n",bty="l",
    
27.      cex.lab=1.8,cex=1.2,cex.axis=1.3)
    
28. lines(x,y,lwd=2)
    
29. text(locator(1),"Posterior",adj=0,cex=1.5)
    
30. 
    
31. postmode <- (alphap-1)/betap
    
32. postmean <- alphap/betap
    
33. postvar <- alphap/betap^2
    
34. 
    
35. postmedian <- qgamma(0.5,alphap,betap)
    
36. 
    
37. postmode;postmean;postmedian;postvar;sqrt(postvar)
    
38. 
    
39. #Bayesian CLT
    
40. 
    
41. bmean <- mdmft
    
42. bvar <- mdmft/n
    
43. bsd <- sqrt(bvar)
    
44. w <- dnorm(x, bmean, bsd)
    
45. lines(x,w,lwd=3,lty=8,col="red")
    
46. text(locator(1),"BCLT",col= "red",cex=1.5,adj=1)
    
47. 
    
48. 
    
49. # based on all children
    
50. 
    
51. 
    
52. dmft <- V1
    
53. mdmft <- mean(dmft)
    
54. n <- length(dmft)
    
55. 
    
56. x <- seq(2.15,2.35,0.001)
    
57. 
    
58. alphap <- 9761
    
59. betap <- 4352
    
60. 
    
61. y <- dgamma(x, shape=alphap, rate = betap)
    
62. 
    
63. par(mfrow=c(1,1))
    
64. par(pty="m")
    
65. 
    
66. plot(x,y,xlab=expression(theta),ylab="",type="n",bty="l",
    
67.      cex.lab=1.8,cex=1.2,cex.axis=1.3)
    
68. lines(x,y,lwd=2)
    
69. text(locator(1),"Posterior",adj=0,cex=1.5)
    
70. 
    
71. postmode <- (alphap-1)/betap
    
72. postmean <- alphap/betap
    
73. postvar <- alphap/betap^2
    
74. 
    
75. postmedian <- qgamma(0.5,alphap,betap)
    
76. 
    
77. postmode;postmean;postmedian;postvar;sqrt(postvar)
    
78. 
    
79. #Bayesian CLT
    
80. 
    
81. bmean <- mdmft
    
82. bvar <- mdmft/n
    
83. bsd <- sqrt(bvar)
    
84. w <- dnorm(x, bmean, bsd)
    
85. 
    
86. lines(x,w,lwd=3,lty=8,col="red")
    
87. text(locator(1),"BCLT",col= "red",cex=1.5,adj=1)

复制代码

1. # midpoint approach for combination log-normal & Poisson
   
2. 
   
3. attach(dmft)
   
4. objects(2)
   
5. 
   
6. par(mfrow=c(1,3))
   
7. par(pty="s")
   
8. par(mar=c(5,5,4,2)+0.1)
   
9. 
   
10. gridsize <- 0.01
    
11. gtheta <- seq(0.01,10,gridsize)
    
12. ntheta <- length(gtheta)
    
13. 
    
14. # prior = log-normal
    
15. 
    
16. sigma0 <- 0.5
    
17. mu0 <- log(2)
    
18. 
    
19. fact1 <- 1/(gtheta*sigma0*sqrt(2*pi))
    
20. fact2 <- exp(-0.5*(log(gtheta)-mu0)^2/sigma0^2)
    
21. priortheta <- fact1*fact2
    
22. 
    
23. 
    
24. plot(gtheta,priortheta,xlab=expression(theta),
    
25.      ylab="Prior density",type="n",bty="l",cex.lab=1.8,cex=1.2,cex.axis=1.3)
    
26. lines(gtheta,priortheta,lwd=2,col="red")
    
27. text(4,0.2,"Log-normal prior",adj=0,cex=1.5,col="red")
    
28. 
    
29. 
    
30. dmft <- V1[1:10]
    
31. mdmft <- mean(dmft)
    
32. sdmft <- sum(dmft)
    
33. n <- length(dmft)
    
34. 
    
35. poissonlik <- exp(-n*gtheta)*gtheta^sdmft
    
36. proppost <- exp(-n*gtheta)*fact2*gtheta^(sdmft-1)
    
37. 
    
38. plot(gtheta,proppost,xlab=expression(theta),ylab="k x Posterior density",
    
39.      cex.lab=1.8,cex=1.2,cex.axis=1.3,
    
40.      type="n",xlim=c(1,5),bty="l")
    
41. lines(gtheta,proppost,lwd=2,col="red")
    
42. 
    
43. # calculation AUC using midpoint of interval +
    
44. # function evaluation at the midpoint
    
45. 
    
46. left <- gtheta[1:(ntheta-1)]
    
47. right <- gtheta[2:ntheta]
    
48. mid <- (left+right)/2
    
49. 
    
50. fact2mid <- exp(-0.5*(log(mid)-mu0)^2/sigma0^2)
    
51. proppostmid <- exp(-n*mid)*fact2mid*mid^(sdmft-1)
    
52. auc <- sum(proppostmid)*gridsize
    
53. auc
    
54. 
    
55. postmid <- proppost/auc
    
56. 
    
57. 
    
58. plot(gtheta,postmid,xlab=expression(theta),ylab="Posterior density",
    
59.      cex.lab=1.8,cex=1.2,cex.axis=1.3,
    
60.      type="n",xlim=c(1,5),bty="l")
    
61. lines(gtheta,postmid,lwd=2,col="red")
    
62. 
    
63. 
    
64. 
    
65. # calculation of mean using midpoint of interval
    
66. 
    
67. meangrid <- proppostmid%*%mid*gridsize/auc
    
68. meangrid
    
69. 
    
70. # calculation of variance using midpoint of interval
    
71. 
    
72. vargrid <- proppostmid%*%mid^2*gridsize/auc - meangrid^2
    
73. vargrid

复制代码