[贝叶斯必读]Bayesian Cognitive Modelling(using R)

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

• http://www.amazon.com/Bayesian-Cognitive-Modeling-Practical-Course/dp/1107018455/ref=sr_1_1?ie=UTF8&qid=1369156307&sr=8-1&keywords=bayesian+cognitive+modeling+a+practical+course
• http://bayesmodels.com/
  

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享2 收藏7 回帖

关键词：Modelling Cognitive Bayesian modelli modell

1. # R2WinBUGS Example;
   
2. # When you work through the code for the first time,
   
3. # execute each command one at a time to better understand
   
4. # what it does.
   
5. 
   
6. # clears workspace:  
   
7. rm(list=ls())
   
8. 
   
9. # sets working directories:
   
10. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
    
11. 
    
12. library(R2WinBUGS)
    
13. bugsdir <- "C:/Program Files/WinBUGS14"
    
14. 
    
15. k <- 5
    
16. n <- 10
    
17. 
    
18. data <- list("k", "n") # to be passed on to WinBUGS
    
19. 
    
20. myinits <- list(
    
21.   list(theta = 0.1), #chain 1 starting value
    
22.   list(theta = 0.9)) #chain 2 starting value
    
23. 
    
24. # parameters to be monitored:       
    
25. parameters <- c("theta")
    
26. 
    
27. # The following command calls WinBUGS with specific options.
    
28. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
29. samples <- bugs(data, inits=myinits, parameters,
    
30.                                  model.file ="Rate_1.txt",
    
31.                                  n.chains=2, n.iter=20000, n.burnin=1, n.thin=1,
    
32.                                  DIC=T, bugs.directory=bugsdir,
    
33.                                  codaPkg=F, debug=F)
    
34. # Now the values for the monitored parameters are in the "samples" object,
    
35. # ready for inspection.
    
36. 
    
37. # The commands below are useful for a quick overview:
    
38. print(samples)  # a rough summary
    
39. plot(samples)   # a visual representation
    
40. names(samples)  # summarizes the variables
    
41. samples$summary # more detailed summary
    
42. chain <- 1
    
43. samples$sims.array[1:15,chain,]# array: element, chain, column (theta/deviance)
    
44. 
    
45. # Collect posterior samples across all chains:
    
46. theta <- samples$sims.list$theta
    
47. 
    
48. # Now let's plot a histogram for theta.
    
49. # First, some options to make the plot look better:
    
50. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
51.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
52. Nbreaks <- 80
    
53. y       <- hist(theta, Nbreaks, plot=F)
    
54. plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
    
55.      xlim=c(0,1), ylim=c(0,10), xlab="Rate", ylab="Posterior Density")
    
56. # NB. ylim=c(0,10) defines the range of the y-axis. Adjust the upper value
    
57. # in case your posterior distribution falls partly outside this range.

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
   
6. 
   
7. library(R2WinBUGS)
   
8. bugsdir <- "C:/Program Files/WinBUGS14"
   
9. 
   
10. k1 <- 5
    
11. k2 <- 7
    
12. n1 <- 10
    
13. n2 <- 10
    
14. 
    
15. data <- list("k1", "k2", "n1", "n2") # to be passed on to WinBUGS
    
16. myinits <-        list(
    
17.   list(theta1 = 0.1, theta2 = 0.9))
    
18. 
    
19. # parameters to be monitored:       
    
20. parameters <- c("delta", "theta1", "theta2")
    
21. 
    
22. # The following command calls WinBUGS with specific options.
    
23. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
24. samples <- bugs(data, inits=myinits, parameters,
    
25.                                  model.file ="Rate_2.txt",
    
26.                                  n.chains=1, n.iter=10000, n.burnin=1, n.thin=1,
    
27.                                  DIC=T, bugs.directory=bugsdir,
    
28.                                  codaPkg=F, debug=F)
    
29. # Now the values for the monitored parameters are in the "samples" object,
    
30. # ready for inspection.
    
31. 
    
32. # Collect posterior samples:
    
33. delta <- samples$sims.list$delta
    
34. 
    
35. # Now let's plot a histogram for delta.
    
36. # First, some options to make the plot look better:
    
37. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
38.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
39. Nbreaks <- 80
    
40. y       <- hist(delta, Nbreaks, plot=F)
    
41. plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
    
42.      xlim=c(-1,1), ylim=c(0,10), xlab="Difference in Rates", ylab="Posterior Density")
    
43. 
    
44. # mean of delta:
    
45. mean(delta)
    
46. # median of delta:
    
47. median(delta)
    
48. # mode of delta, estimated from the "density" smoother:
    
49. density(delta)$x[which(density(delta)$y==max(density(delta)$y))]
    
50. # 95% credible interval for delta:
    
51. quantile(delta, c(.025,.975))

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
   
6. 
   
7. library(R2WinBUGS)
   
8. bugsdir <- "C:/Program Files/WinBUGS14"
   
9. 
   
10. k1 <- 5
    
11. k2 <- 7
    
12. n1 <- 10
    
13. n2 <- 10
    
14. 
    
15. data <- list("k1", "k2", "n1", "n2") # to be passed on to WinBUGS
    
16. myinits <-        list(
    
17.   list(theta = 0.5))
    
18. 
    
19. # parameters to be monitored:       
    
20. parameters <- c("theta")
    
21. 
    
22. # The following command calls WinBUGS with specific options.
    
23. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
24. samples <- bugs(data, inits=myinits, parameters,
    
25.                                  model.file ="Rate_3.txt",
    
26.                                  n.chains=1, n.iter=1000, n.burnin=1, n.thin=1,
    
27.                                  DIC=T, bugs.directory=bugsdir,
    
28.                                  codaPkg=F, debug=T)
    
29. # NB. This is the first time we use debug=T. It is often useful to study the
    
30. # WinBUGS output before returning to R. When you are do exploring the chains
    
31. # in WinBUGS, you have to close (not minimize) WinBUGS in order to return to R.
    
32. 
    
33. theta <- samples$sims.list$theta
    
34.                                 
    
35. # Now let's plot a histogram for theta.
    
36. # First, some options to make the plot look better:
    
37. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
38.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
39. Nbreaks <- 80
    
40. y       <- hist(theta, Nbreaks, plot=F)
    
41. plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
    
42.      xlim=c(0,1), ylim=c(0,10), xlab="Rate", ylab="Posterior Density")

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
   
6. 
   
7. library(R2WinBUGS)
   
8. bugsdir <- "C:/Program Files/WinBUGS14"
   
9. 
   
10. k <- 1
    
11. n <- 15
    
12. # Uncomment for Trompetter Data
    
13. # k <- 24
    
14. # n <- 121
    
15. 
    
16. data <- list("k", "n") # to be passed on to WinBUGS
    
17. myinits <-        list(
    
18.   list(theta = 0.5, thetaprior = 0.5))
    
19. 
    
20. # parameters to be monitored:       
    
21. parameters <- c("theta", "thetaprior", "postpredk", "priorpredk")
    
22. 
    
23. # The following command calls WinBUGS with specific options.
    
24. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
25. samples <- bugs(data, inits=myinits, parameters,
    
26.                                  model.file ="Rate_4.txt",
    
27.                                  n.chains=1, n.iter=5000, n.burnin=1, n.thin=1,
    
28.                                  DIC=T, bugs.directory=bugsdir,
    
29.                                  codaPkg=F, debug=T)
    
30. # Now the values for the monitored parameters are in the "samples" object,
    
31. # ready for inspection.
    
32. 
    
33. ######################Plots########################################################################
    
34. theta      <- samples$sims.list$theta
    
35. thetaprior <- samples$sims.list$thetaprior
    
36. postpredk  <- samples$sims.list$postpredk
    
37. priorpredk <- samples$sims.list$priorpredk
    
38. 
    
39. layout(matrix(c(1,2),2,1))
    
40. layout.show(2)
    
41. #Prior and posterior of theta
    
42. plot(density(theta, from=0, to=1), zero.line=F, axes=F, main="", xlab="", ylab="", xlim=c(0,1), ylim=c(0,6))
    
43. axis(1, at=c(0,0.2,0.4,0.6,0.8,1), lab=c("0","0.2","0.4","0.6","0.8","1"),cex.axis=0.8)
    
44. mtext("Rate", side=1, line=2.25, cex=1.2)
    
45. axis(2, at=c(0,2,4,6),cex.axis=0.8)
    
46. mtext("Density", side=2, line=2.25, cex=1.2)
    
47. lines(density(thetaprior, from=0, to=1), lty=3, col="gray")
    
48. legend(0.6,5.75, c("Prior", "Posterior"), lty=c(3,1), col=c ("grey", "black"))
    
49. 
    
50. #Prior and posterior predictive
    
51. mybreaks <- seq(from=-.5,to=n+1,by=1)
    
52. my.at    <- seq(from=0,to=n,by=1)
    
53. hist(postpredk,breaks=mybreaks,freq=F, right=F, ylab="", xlab="", ylim=c(0,0.3),main="", axes=F )
    
54. axis(1, at=my.at,lab=my.at,cex.axis=0.8)
    
55. mtext("Success Count", side=1, line=2.25, cex=1.2)
    
56. axis(2,at=c(0,0.1,0.2,0.3),lab=c("0","0.1","0.2","0.3"),cex.axis=0.8)
    
57. mtext("Mass", side=2, line=2.25, cex=1.2)
    
58. hist(priorpredk, breaks=mybreaks,freq=F,right=F,add=T, lty=3,border="grey")
    
59. legend(8,0.3, c("Prior", "Posterior"), lty=c(3,1),col=c("grey", "black"))

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
   
6. 
   
7. library(R2WinBUGS)
   
8. bugsdir <- "C:/Program Files/WinBUGS14"
   
9. 
   
10. k1 <- 0
    
11. k2 <- 10
    
12. n1 <- 10
    
13. n2 <- 10
    
14. 
    
15. data <- list("k1", "k2", "n1", "n2") # to be passed on to WinBUGS
    
16. myinits <-        list(
    
17.   list(theta = 0.5))
    
18. 
    
19. # parameters to be monitored:       
    
20. parameters <- c("theta", "postpredk1", "postpredk2")
    
21. 
    
22. # The following command calls WinBUGS with specific options.
    
23. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
24. samples <- bugs(data, inits=myinits, parameters,
    
25.                                  model.file ="Rate_5.txt",
    
26.                                  n.chains=1, n.iter=1000, n.burnin=1, n.thin=1,
    
27.                                  DIC=T, bugs.directory=bugsdir,
    
28.                                  codaPkg=F, debug=T)
    
29. 
    
30. theta      <- samples$sims.list$theta
    
31. postpredk1 <- samples$sims.list$postpredk1
    
32. postpredk2 <- samples$sims.list$postpredk2
    
33.                                                                   
    
34. # Two-panel plot.
    
35. layout(matrix(c(1,2),1,2))
    
36. layout.show(2)
    
37. # First, a histogram for theta.
    
38. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
39.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
40. Nbreaks <- 80
    
41. y       <- hist(theta, Nbreaks, plot=F)
    
42. plot(c(y$breaks, max(y$breaks)), c(0,y$density,0), type="S", lwd=2, lty=1,
    
43.      xlim=c(0,1), ylim=c(0,10), xlab="Theta", ylab="Density")
    
44. # let's plot a density estimate over this:
    
45. lines(density(theta), col="red", lwd=2)
    
46. 
    
47. # Second plot, the data space (predictives)
    
48. plot(k1,k2,type="p", pch=4, cex=2, lwd=2, xlab="Success Count 1", ylab="Success Count 2",
    
49.      xlim=c(-1, n1+1), ylim=c(-1,n2+1))        
    
50. nsamples <- length(theta)
    
51. sc <- 10
    
52. for (i in 0:n1)
    
53. {
    
54.   for (j in 0:n2)
    
55.   {
    
56.     match.preds <- sum(postpredk1==i & postpredk2==j)/nsamples
    
57.     if (match.preds > 0)
    
58.     {
    
59.       points(i,j, pch=0, cex=sc*sqrt(match.preds))
    
60.     }
    
61.   }
    
62. }

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/Binomial")
   
6. 
   
7. library(R2WinBUGS)
   
8. bugsdir <- "C:/Program Files/WinBUGS14"
   
9. 
   
10. nmax <- 500
    
11. k    <- c(16,18,22,25,27)
    
12. m    <- length(k)
    
13. 
    
14. data <- list("nmax", "k", "m") # to be passed on to WinBUGS
    
15. myinits <-        list(
    
16.   list(theta = 0.5, n = nmax/2))
    
17. 
    
18. # parameters to be monitored:       
    
19. parameters <- c("theta", "n")
    
20. 
    
21. # The following command calls WinBUGS with specific options.
    
22. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
23. samples <- bugs(data, inits=myinits, parameters,
    
24.                                  model.file ="Survey.txt",
    
25.                                  n.chains=1, n.iter=5000, n.burnin=1, n.thin=1,
    
26.                                  DIC=T, bugs.directory=bugsdir,
    
27.                                  codaPkg=F, debug=T)
    
28. # Now the values for the monitored parameters are in the "samples" object,
    
29. # ready for inspection.
    
30. 
    
31. theta      <- samples$sims.list$theta
    
32. n          <- samples$sims.list$n
    
33. 
    
34. ## First calculate MLE:
    
35. cc <- -Inf
    
36. ind <- 0
    
37. 
    
38. for (i in 1:length(n))
    
39. {
    
40.         logL <- 0
    
41.         for(j in 1:m)
    
42.         {               
    
43.                 logL <- logL+lgamma(n[i]+1)-lgamma(k[j]+1)-lgamma(n[i]-k[j]+1)
    
44.                 logL <- logL+k[j]*log(theta[i])+(n[i]-k[j])*log(1-theta[i])
    
45.         }
    
46.         if (logL>cc)
    
47.   {
    
48.                 ind <- i
    
49.                 cc <- logL
    
50.         }
    
51. }
    
52. # end MLE
    
53. 
    
54. ######################Plots########################################################################
    
55. layout(matrix(c(2,0,1,3),2,2,byrow=T), width=c(2/3, 1/3), heights=c(1/3,2/3))
    
56. xhist <- hist(n,plot=F)
    
57. yhist <- hist(theta,plot=F)
    
58. top <- max(c(xhist$counts, yhist$counts))
    
59. xrange <- c(0,nmax)
    
60. yrange <- c(0,1)
    
61. 
    
62. par(mar=c(5,5,1,1))
    
63. plot(n,theta,xlim=xrange, ylim=yrange,ylab="", xlab="")
    
64. axis(1)
    
65. mtext("Number of Surveys", side=1,line=2.25, cex=1.2)
    
66. axis(2, cex=1.2)
    
67. las=0
    
68. mtext("Rate of Return", side=2 ,line=2.25, cex=1.2)
    
69. las=1
    
70. points(mean(n),mean(theta), col="red", lwd=3, pch=4) #expectation
    
71. points(n[ind],theta[ind], col="green", lwd=3, pch=10) #Maximum Likelihood
    
72. 
    
73. par(mar=c(0,4,1,1))
    
74. barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0,col="lightblue")
    
75. 
    
76. par(mar=c(4,0,1,3))
    
77. barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE,col="lightblue")

复制代码

1. # clears workspace:  
   
2. rm(list=ls())
   
3. 
   
4. # sets working directories:
   
5. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/DataAnalysis")
   
6. library(R2WinBUGS)
   
7. bugsdir <- "C:/Program Files/WinBUGS14"
   
8. 
   
9. c <- scan("changepointdata.txt")
   
10. n <- length(c)
    
11. t <- 1:n
    
12. 
    
13. data <- list("c", "n", "t") # to be passed on to WinBUGS
    
14. myinits <- list(
    
15.   list(mu = c(1,1), lambda = 1, tau = n/2))
    
16. 
    
17. # parameters to be monitored:       
    
18. parameters <- c("mu","sigma","tau")
    
19. 
    
20. # The following command calls WinBUGS with specific options.
    
21. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
22. samples <- bugs(data, inits=myinits, parameters,
    
23.                                  model.file ="ChangeDetection.txt",
    
24.                                  n.chains=1, n.iter=1000, n.burnin=1, n.thin=1,
    
25.                                  DIC=T, bugs.directory=bugsdir,
    
26.                                  codaPkg=F, debug=T)
    
27. # Now the values for the monitored parameters are in the "samples" object,
    
28. # ready for inspection.
    
29. 
    
30. plot(samples)
    
31. mean.tau <- samples$mean$tau
    
32. mean.mu1 <- samples$mean$mu[1]
    
33. mean.mu2 <- samples$mean$mu[2]
    
34. 
    
35. #some plotting options to make things look better:
    
36. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
37.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
38. # the plotting:
    
39. plot(c, type="l", main="", ylab="Values", xlab="Samples")
    
40. lines(c(1, mean.tau), c(mean.mu1,mean.mu1), lwd=2, col="red")
    
41. lines(c(mean.tau+1,length(c)), c(mean.mu2,mean.mu2), lwd=2, col="red")

复制代码

1. #### After 949 failed attempts, Cha Sa-Soon from Korea finally passed   
   
2. #### her theoretical drivers exam. What can we infer about theta, the
   
3. #### probability of answering any one question correctly?
   
4.      
   
5. # clears workspace:  
   
6. rm(list=ls())
   
7. 
   
8. # sets working directories:
   
9. setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/DataAnalysis")
   
10. library(R2WinBUGS)
    
11. bugsdir = "C:/Program Files/WinBUGS14"
    
12. 
    
13. nattempts  <- 950   
    
14. nfails     <- 949   
    
15. n          <- 50    # Number of questions
    
16. z = c(rep(NA,nfails),30)
    
17. y = c(rep(1,nfails),0)
    
18. 
    
19. data <- list("nattempts","n","z","y") # to be passed on to WinBUGS
    
20. 
    
21. myinits <- list(
    
22.   list(theta = 0.5))
    
23. 
    
24. parameters <- c("theta")
    
25. 
    
26. # The following command calls WinBUGS with specific options.
    
27. # For a detailed description see Sturtz, Ligges, & Gelman (2005).
    
28. samples <- bugs(data, inits=myinits, parameters,
    
29.                                  model.file ="ChaSaSoon.txt",
    
30.                                  n.chains=1, n.iter=4100, n.burnin=100, n.thin=1,
    
31.                                  DIC=T, bugs.directory=bugsdir, #DIC=T or else error message
    
32.                                  codaPkg=F, debug=T)
    
33. # Now the values for theta are in the "samples" object, ready for inspection.
    
34. 
    
35. # Collect all samples in "theta":
    
36. theta <- samples$sims.list$theta
    
37. 
    
38. # Plot the posterior for theta:
    
39. par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    
40.     font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
    
41. plot(density(theta), ylim=c(0,200), xlim=c(0.2,0.5), lty=1, lwd=2,
    
42.      axes=F, main=" ", xlab=" ", ylab="Posterior Density")
    
43. axis(1, at = c(0.20, 0.30, 0.40, 0.50), lab=c("0.20", "0.30", "0.40", "0.50"))
    
44. axis(2)
    
45. mtext(expression(paste(theta," for Cha Sa-soon")), side=1, line = 2.8, cex=2)
    
46. 
    
47. # plot 95% confidence interval
    
48. x0 <- quantile(theta,p=c(.025,.975))[[1]]
    
49. x1 <- quantile(theta,p=c(.025,.975))[[2]]
    
50. arrows(x0, 150, x1, 150, length = 0.05, angle = 90, code = 3, lwd=2)
    
51. text((x0+x1)/2, 170, labels="95%", cex = 1.5)
    
52. text(x0-.03, 150, labels=as.character(round(x0,2)), cex = 1.5)
    
53. text(x0+.04, 150, labels=as.character(round(x1,2)), cex = 1.5)

复制代码

感谢楼主的分享，顶一个