[貝葉斯專著]Doing Bayesian Data Analysis: A Tutorial with R and BUGS

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

Doing Bayesian Data Analysis: A Tutorial with R and BUGS

John Kruschke

Book Description
Publication Date: Oct. 27 2010 |
ISBN-10: 0123814855 | ISBN-13: 978-0123814852 |
Edition: 1
There is an explosion of interest in Bayesian statistics, primarily because recently created computational methods have finally made Bayesian analysis tractable and accessible to a wide audience. Doing Bayesian Data Analysis, A Tutorial Introduction with R and BUGS, is for first year graduate students or advanced undergraduates and provides an accessible approach, as all mathematics is explained intuitively and with concrete examples. It assumes only algebra and 'rusty' calculus. Unlike other textbooks, this book begins with the basics, including essential concepts of probability and random sampling. The book gradually climbs all the way to advanced hierarchical modeling methods for realistic data. The text provides complete examples with the R programming language and BUGS software (both freeware), and begins with basic programming examples, working up gradually to complete programs for complex analyses and presentation graphics. These templates can be easily adapted for a large variety of students and their own research needs.The textbook bridges the students from their undergraduate training into modern Bayesian methods.
-Accessible, including the basics of essential concepts of probability and random sampling
-Examples with R programming language and BUGS software
-Comprehensive coverage of all scenarios addressed by non-bayesian textbooks- t-tests, analysis of variance (ANOVA) and comparisons in ANOVA, multiple regression, and chi-square (contingency table analysis).
-Coverage of experiment planning
-R and BUGS computer programming code on website
-Exercises have explicit purposes and guidelines for accomplishment
Product Details

• Hardcover: 672 pages
• Publisher: Academic Press; 1 edition (Oct. 27 2010)
• Language: English
• ISBN-10: 0123814855
• ISBN-13: 978-0123814852
• Product Dimensions: 23.6 x 19.3 x 3.3 cm
• Shipping Weight: 1.2 Kg
  

Product DescriptionReview"I think it fills a gaping hole in what is currently available, and will serve to create its own market as researchers and their students transition towards the routine application of Bayesian statistical methods.” -Prof. Michael lee, University of California, Irvine, and president of the Society for Mathematical Psychology
"Kruschke's text covers a much broader range of traditional experimental designs.has the potential to change the way most cognitive scientists and experimental psychologists approach the planning and analysis of their experiments" -Prof. Geoffrey Iverson, University of California, Irvine, and past president of the Society for Mathematical Psychology
"John Kruschke has written a book on Statistics. It's better than others for reasons stylistic. It also is better because itis Bayesian. To find out why, buy it -- it's truly amazin'!”-James L. (Jay) McClelland, Lucie Stern Professor & Chair, Dept. Of Psychology, Standford University

From the Back CoverThere is an explosion of interest in Bayesian statistics, primarily because recently created computational methods have finally made Bayesian analysis obtainable to a wide audience. Doing Bayesian Data Analysis, A Tutorial Introduction with R and BUGS, provides an accessible approach to Bayesian Data Analysis, as material is explained clearly with concrete examples. The book begins with the basics, including essential concepts of probability and random sampling, and gradually progresses to advanced hierarchical modeling methods for realistic data. The text delivers comprehensive coverage of all scenarios addressed by non-Bayesian textbooks- t-tests, analysis of variance (ANOVA) and comparisons in ANOVA, correlation, multiple regression, and chi-square (contingency table analysis).
This book is intended for first year graduate students or advanced undergraduates. It provides a bridge between undergraduate training and modern Bayesian methods for data analysis, which is becoming the accepted research standard. Prerequisite is knowledge of algebra and basic calculus.

• https://bbs.pinggu.org/thread-1048913-1-1.html
• https://sites.google.com/site/doingbayesiandataanalysis/1st-ed-stuff
• https://github.com/boboppie/kruschke-doing_bayesian_data_analysis
  

扫码加我 拉你入群

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

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

分享0 收藏4 回帖

关键词：Bayesian Tutorial Analysis Analysi alysis accessible recently audience interest because

1. graphics.off()
   
2. rm(list=ls(all=TRUE))
   
3. fnroot = "ANOVAonewayBrugs"
   
4. library(BRugs)         # Kruschke, J. K. (2010). Doing Bayesian data analysis:
   
5.                        # A Tutorial with R and BUGS. Academic Press / Elsevier.
   
6. #------------------------------------------------------------------------------
   
7. # THE MODEL.
   
8. 
   
9. modelstring = "
   
10. # BUGS model specification begins here...
    
11. model {
    
12.   for ( i in 1:Ntotal ) {
    
13.     y[i] ~ dnorm( mu[i] , tau )
    
14.     mu[i] <- a0 + a[x[i]]
    
15.   }
    
16.   #
    
17.   tau <- pow( sigma , -2 )
    
18.   sigma ~ dunif(0,10) # y values are assumed to be standardized
    
19.   #
    
20.   a0 ~ dnorm(0,0.001) # y values are assumed to be standardized
    
21.   #
    
22.   for ( j in 1:NxLvl ) { a[j] ~ dnorm( 0.0 , atau ) }
    
23.   atau <- 1 / pow( aSD , 2 )
    
24.   aSD <- abs( aSDunabs ) + .1
    
25.   aSDunabs ~ dt( 0 , 0.001 , 2 )
    
26. }
    
27. # ... end BUGS model specification
    
28. " # close quote for modelstring
    
29. # Write model to a file, and send to BUGS:
    
30. writeLines(modelstring,con="model.txt")
    
31. modelCheck( "model.txt" )
    
32. 
    
33. #------------------------------------------------------------------------------
    
34. # THE DATA.
    
35. 
    
36. # Specify data source:
    
37. dataSource = c( "McDonaldSK1991" , "SolariLS2008" , "Random" )[1]
    
38. # Load the data:
    
39. 
    
40. if ( dataSource == "McDonaldSK1991" ) {
    
41.   fnroot = paste( fnroot , dataSource , sep="" )
    
42.   datarecord = read.table( "McDonaldSK1991data.txt", header=T ,
    
43.                            colClasses=c("factor","numeric") )
    
44.   y = as.numeric(datarecord$Size)
    
45.   Ntotal = length(datarecord$Size)
    
46.   x = as.numeric(datarecord$Group)
    
47.   xnames = levels(datarecord$Group)
    
48.   NxLvl = length(unique(datarecord$Group))
    
49.   contrastList = list( BIGvSMALL = c(-1/3,-1/3,1/2,-1/3,1/2) ,
    
50.                        ORE1vORE2 = c(1,-1,0,0,0) ,
    
51.                        ALAvORE = c(-1/2,-1/2,1,0,0) ,
    
52.                        NPACvORE = c(-1/2,-1/2,1/2,1/2,0) ,
    
53.                        USAvRUS = c(1/3,1/3,1/3,-1,0) ,
    
54.                        FINvPAC = c(-1/4,-1/4,-1/4,-1/4,1) ,
    
55.                        ENGvOTH = c(1/3,1/3,1/3,-1/2,-1/2) ,
    
56.                        FINvRUS = c(0,0,0,-1,1) )
    
57. }
    
58. 
    
59. if ( dataSource == "SolariLS2008" ) {
    
60.   fnroot = paste( fnroot , dataSource , sep="" )
    
61.   datarecord = read.table("SolariLS2008data.txt", header=T ,
    
62.                            colClasses=c("factor","numeric") )
    
63.   y = as.numeric(datarecord$Acid)
    
64.   Ntotal = length(datarecord$Acid)
    
65.   x = as.numeric(datarecord$Type)
    
66.   xnames = levels(datarecord$Type)
    
67.   NxLvl = length(unique(datarecord$Type))
    
68.   contrastList = list( G3vOTHER = c(-1/8,-1/8,1,-1/8,-1/8,-1/8,-1/8,-1/8,-1/8) )
    
69. }
    
70. 
    
71. if ( dataSource == "Random" ) {
    
72.   fnroot = paste( fnroot , dataSource , sep="" )
    
73.   #set.seed(47405)
    
74.   ysdtrue = 4.0
    
75.   a0true = 100
    
76.   atrue = c( 2 , -2 ) # sum to zero
    
77.   npercell = 8
    
78.   datarecord = matrix( 0, ncol=2 , nrow=length(atrue)*npercell )
    
79.   colnames(datarecord) = c("y","x")
    
80.   rowidx = 0
    
81.   for ( xidx in 1:length(atrue) ) {
    
82.     for ( subjidx in 1:npercell ) {
    
83.       rowidx = rowidx + 1
    
84.       datarecord[rowidx,"x"] = xidx
    
85.       datarecord[rowidx,"y"] = ( a0true + atrue[xidx] + rnorm(1,0,ysdtrue) )
    
86.     }
    
87.   }
    
88.   datarecord = data.frame( y=datarecord[,"y"] , x=as.factor(datarecord[,"x"]) )
    
89.   y = as.numeric(datarecord$y)
    
90.   Ntotal = length(y)
    
91.   x = as.numeric(datarecord$x)
    
92.   xnames = levels(datarecord$x)
    
93.   NxLvl = length(unique(x))
    
94.   # Construct list of all pairwise comparisons, to compare with NHST TukeyHSD:
    
95.   contrastList = NULL
    
96.   for ( g1idx in 1:(NxLvl-1) ) {
    
97.     for ( g2idx in (g1idx+1):NxLvl ) {
    
98.       cmpVec = rep(0,NxLvl)
    
99.       cmpVec[g1idx] = -1
    
100.       cmpVec[g2idx] = 1
     
101.       contrastList = c( contrastList , list( cmpVec ) )
     
102.     }
     
103.   }
     
104. }
     
105. 
     
106. # Specify the data in a form that is compatible with BRugs model, as a list:
     
107. ySDorig = sd(y)
     
108. yMorig = mean(y)
     
109. z = ( y - yMorig ) / ySDorig
     
110. datalist = list(
     
111.   y = z ,
     
112.   x = x ,
     
113.   Ntotal = Ntotal ,
     
114.   NxLvl = NxLvl
     
115. )
     
116. # Get the data into BRugs:
     
117. modelData( bugsData( datalist ) )
     
118. 
     
119. #------------------------------------------------------------------------------
     
120. # INTIALIZE THE CHAINS.
     
121. 
     
122. # Autocorrelation within chains is large, so use several chains to reduce
     
123. # degree of thinning. But we still have to burn-in all the chains, which takes
     
124. # more time with more chains (on serial CPUs).
     
125. nchain = 5
     
126. modelCompile( numChains = nchain )
     
127. 
     
128. if ( F ) {
     
129.    modelGenInits() # often won't work for diffuse prior
     
130. } else {
     
131.   #  initialization based on data
     
132.   theData = data.frame( y=datalist$y , x=factor(x,labels=xnames) )
     
133.   a0 = mean( theData$y )
     
134.   a = aggregate( theData$y , list( theData$x ) , mean )[,2] - a0
     
135.   ssw = aggregate( theData$y , list( theData$x ) ,
     
136.                   function(x){var(x)*(length(x)-1)} )[,2]
     
137.   sp = sqrt( sum( ssw ) / length( theData$y ) )
     
138.   genInitList <- function() {
     
139.     return(
     
140.         list(
     
141.             a0 = a0 ,
     
142.             a = a ,
     
143.             sigma = sp ,
     
144.             aSDunabs = sd(a)
     
145.         )
     
146.     )
     
147.   }
     
148.   for ( chainIdx in 1 : nchain ) {
     
149.     modelInits( bugsInits( genInitList ) )
     
150.   }
     
151. }
     
152. 
     
153. #------------------------------------------------------------------------------
     
154. # RUN THE CHAINS
     
155. 
     
156. # burn in
     
157. BurnInSteps = 10000
     
158. modelUpdate( BurnInSteps )
     
159. # actual samples
     
160. samplesSet( c( "a0" ,  "a" , "sigma" , "aSD" ) )
     
161. stepsPerChain = ceiling(5000/nchain)
     
162. thinStep = 750
     
163. modelUpdate( stepsPerChain , thin=thinStep )
     
164. 
     
165. #------------------------------------------------------------------------------
     
166. # EXAMINE THE RESULTS
     
167. 
     
168. source("plotChains.R")
     
169. source("plotPost.R")
     
170. 
     
171. checkConvergence = T
     
172. if ( checkConvergence ) {
     
173.    sumInfo = plotChains( "a0" , saveplots=T , filenameroot=fnroot )
     
174.    sumInfo = plotChains( "a" , saveplots=T , filenameroot=fnroot )
     
175.    sumInfo = plotChains( "sigma" , saveplots=T , filenameroot=fnroot )
     
176.    sumInfo = plotChains( "aSD" , saveplots=T , filenameroot=fnroot )
     
177. }
     
178. 
     
179. # Extract and plot the SDs:
     
180. sigmaSample = samplesSample("sigma")
     
181. aSDSample = samplesSample("aSD")
     
182. windows()
     
183. layout( matrix(1:2,nrow=2) )
     
184. par( mar=c(3,1,2.5,0) , mgp=c(2,0.7,0) )
     
185. plotPost( sigmaSample , xlab="sigma" , main="Cell SD" , breaks=30 )
     
186. plotPost( aSDSample , xlab="aSD" , main="a SD" , breaks=30 )
     
187. dev.copy2eps(file=paste(fnroot,"SD.eps",sep=""))
     
188. 
     
189. # Extract a values:
     
190. a0Sample = samplesSample( "a0" )
     
191. chainLength = length(a0Sample)
     
192. aSample = array( 0 , dim=c( datalist$NxLvl , chainLength ) )
     
193. for ( xidx in 1:datalist$NxLvl ) {
     
194.    aSample[xidx,] = samplesSample( paste("a[",xidx,"]",sep="") )
     
195. }
     
196. 
     
197. # Convert to zero-centered b values:
     
198. mSample = array( 0, dim=c( datalist$NxLvl , chainLength ) )
     
199. for ( stepIdx in 1:chainLength ) {
     
200.     mSample[,stepIdx ] = ( a0Sample[stepIdx] + aSample[,stepIdx] )
     
201. }
     
202. b0Sample = apply( mSample , 2 , mean )
     
203. bSample = mSample - matrix(rep( b0Sample ,NxLvl),nrow=NxLvl,byrow=T)
     
204. # Convert from standardized b values to original scale b values:
     
205. b0Sample = b0Sample * ySDorig + yMorig
     
206. bSample = bSample * ySDorig
     
207. 
     
208. # Plot b values:
     
209. windows(datalist$NxLvl*2.75,2.5)
     
210. layout( matrix( 1:datalist$NxLvl , nrow=1 ) )
     
211. par( mar=c(3,1,2.5,0) , mgp=c(2,0.7,0) )
     
212. for ( xidx in 1:datalist$NxLvl ) {
     
213.     plotPost( bSample[xidx,] , breaks=30 ,
     
214.               xlab=bquote(beta*1[.(xidx)]) ,
     
215.               main=paste("x:",xnames[xidx])  )
     
216. }
     
217. dev.copy2eps(file=paste(fnroot,"b.eps",sep=""))
     
218. 
     
219. # Display contrast analyses
     
220. nContrasts = length( contrastList )
     
221. if ( nContrasts > 0 ) {
     
222.    nPlotPerRow = 5
     
223.    nPlotRow = ceiling(nContrasts/nPlotPerRow)
     
224.    nPlotCol = ceiling(nContrasts/nPlotRow)
     
225.    windows(3.75*nPlotCol,2.5*nPlotRow)
     
226.    layout( matrix(1:(nPlotRow*nPlotCol),nrow=nPlotRow,ncol=nPlotCol,byrow=T) )
     
227.    par( mar=c(4,0.5,2.5,0.5) , mgp=c(2,0.7,0) )
     
228.    for ( cIdx in 1:nContrasts ) {
     
229.        contrast = matrix( contrastList[[cIdx]],nrow=1) # make it a row matrix
     
230.        incIdx = contrast!=0
     
231.        histInfo = plotPost( contrast %*% bSample , compVal=0 , breaks=30 ,
     
232.                 xlab=paste( round(contrast[incIdx],2) , xnames[incIdx] ,
     
233.                             c(rep("+",sum(incIdx)-1),"") , collapse=" " ) ,
     
234.                 cex.lab = 1.0 ,
     
235.                 main=paste( "X Contrast:", names(contrastList)[cIdx] ) )
     
236.    }
     
237.    dev.copy2eps(file=paste(fnroot,"xContrasts.eps",sep=""))
     
238. }
     
239. 
     
240. #==============================================================================
     
241. # Do NHST ANOVA and t tests:
     
242. 
     
243. theData = data.frame( y=y , x=factor(x,labels=xnames) )
     
244. aovresult = aov( y ~ x , data = theData ) # NHST ANOVA
     
245. cat("\n------------------------------------------------------------------\n\n")
     
246. print( summary( aovresult ) )
     
247. cat("\n------------------------------------------------------------------\n\n")
     
248. print( model.tables( aovresult , "means" ) , digits=4 )
     
249. windows()
     
250. boxplot( y ~ x , data = theData )
     
251. cat("\n------------------------------------------------------------------\n\n")
     
252. print( TukeyHSD( aovresult , "x" , ordered = FALSE ) )
     
253. windows()
     
254. plot( TukeyHSD( aovresult , "x" ) )
     
255. if ( T ) {
     
256.   for ( xIdx1 in 1:(NxLvl-1) ) {
     
257.     for ( xIdx2 in (xIdx1+1):NxLvl ) {
     
258.       cat("\n----------------------------------------------------------\n\n")
     
259.       cat( "xIdx1 = " , xIdx1 , ", xIdx2 = " , xIdx2 ,
     
260.            ", M2-M1 = " , mean(y[x==xIdx2])-mean(y[x==xIdx1]) , "\n" )
     
261.       print( t.test( y[x==xIdx2] , y[x==xIdx1] , var.equal=T ) ) # t test
     
262.     }
     
263.   }
     
264. }
     
265. cat("\n------------------------------------------------------------------\n\n")
     
266. 
     
267. #==============================================================================

复制代码

1. graphics.off()
   
2. rm(list=ls(all=TRUE))
   
3. fnroot = "ANOVAonewayBrugsSTZ"
   
4. library(BRugs)         # Kruschke, J. K. (2011). Doing Bayesian data analysis:
   
5.                        # A Tutorial with R and BUGS. Academic Press / Elsevier.
   
6. #------------------------------------------------------------------------------
   
7. # THE MODEL.
   
8. 
   
9. modelstring = "
   
10. # BUGS model specification begins here...
    
11. model {
    
12.   for ( i in 1:Ntotal ) {
    
13.     y[i] ~ dnorm( mu[i] , tau )
    
14.     mu[i] <- a0 + a[x[i]]
    
15.   }
    
16.   #
    
17.   tau <- pow( sigma , -2 )
    
18.   sigma ~ dunif(0,10) # y values are assumed to be standardized
    
19.   #
    
20.   a0 ~ dnorm(0,0.001) # y values are assumed to be standardized
    
21.   #
    
22.   for ( j in 1:NxLvl ) { a[j] ~ dnorm( 0.0 , atau ) }
    
23.   atau <- 1 / pow( aSD , 2 )
    
24.   aSD <- abs( aSDunabs ) + .1
    
25.   aSDunabs ~ dt( 0 , 0.001 , 2 )
    
26.   # Convert a0,a[] to sum-to-zero b0,b[] :
    
27.   for ( j in 1:NxLvl ) { m[j] <- a0 + a[j] }
    
28.   b0 <- mean( m[1:NxLvl] )
    
29.   for ( j in 1:NxLvl ) { b[j] <- m[j] - b0 }
    
30. }
    
31. # ... end BUGS model specification
    
32. " # close quote for modelstring
    
33. # Write model to a file, and send to BUGS:
    
34. writeLines(modelstring,con="model.txt")
    
35. modelCheck( "model.txt" )
    
36. 
    
37. #------------------------------------------------------------------------------
    
38. # THE DATA.
    
39. 
    
40. # Specify data source:
    
41. dataSource = c( "McDonaldSK1991" , "SolariLS2008" , "Random" )[1]
    
42. # Load the data:
    
43. 
    
44. if ( dataSource == "McDonaldSK1991" ) {
    
45.   fnroot = paste( fnroot , dataSource , sep="" )
    
46.   datarecord = read.table( "McDonaldSK1991data.txt", header=T ,
    
47.                            colClasses=c("factor","numeric") )
    
48.   y = as.numeric(datarecord$Size)
    
49.   Ntotal = length(datarecord$Size)
    
50.   x = as.numeric(datarecord$Group)
    
51.   xnames = levels(datarecord$Group)
    
52.   NxLvl = length(unique(datarecord$Group))
    
53.   contrastList = list( BIGvSMALL = c(-1/3,-1/3,1/2,-1/3,1/2) ,
    
54.                        ORE1vORE2 = c(1,-1,0,0,0) ,
    
55.                        ALAvORE = c(-1/2,-1/2,1,0,0) ,
    
56.                        NPACvORE = c(-1/2,-1/2,1/2,1/2,0) ,
    
57.                        USAvRUS = c(1/3,1/3,1/3,-1,0) ,
    
58.                        FINvPAC = c(-1/4,-1/4,-1/4,-1/4,1) ,
    
59.                        ENGvOTH = c(1/3,1/3,1/3,-1/2,-1/2) ,
    
60.                        FINvRUS = c(0,0,0,-1,1) )
    
61. }
    
62. 
    
63. if ( dataSource == "SolariLS2008" ) {
    
64.   fnroot = paste( fnroot , dataSource , sep="" )
    
65.   datarecord = read.table("SolariLS2008data.txt", header=T ,
    
66.                            colClasses=c("factor","numeric") )
    
67.   y = as.numeric(datarecord$Acid)
    
68.   Ntotal = length(datarecord$Acid)
    
69.   x = as.numeric(datarecord$Type)
    
70.   xnames = levels(datarecord$Type)
    
71.   NxLvl = length(unique(datarecord$Type))
    
72.   contrastList = list( G3vOTHER = c(-1/8,-1/8,1,-1/8,-1/8,-1/8,-1/8,-1/8,-1/8) )
    
73. }
    
74. 
    
75. if ( dataSource == "Random" ) {
    
76.   fnroot = paste( fnroot , dataSource , sep="" )
    
77.   #set.seed(47405)
    
78.   ysdtrue = 4.0
    
79.   a0true = 100
    
80.   atrue = c( 2 , -2 ) # sum to zero
    
81.   npercell = 8
    
82.   datarecord = matrix( 0, ncol=2 , nrow=length(atrue)*npercell )
    
83.   colnames(datarecord) = c("y","x")
    
84.   rowidx = 0
    
85.   for ( xidx in 1:length(atrue) ) {
    
86.     for ( subjidx in 1:npercell ) {
    
87.       rowidx = rowidx + 1
    
88.       datarecord[rowidx,"x"] = xidx
    
89.       datarecord[rowidx,"y"] = ( a0true + atrue[xidx] + rnorm(1,0,ysdtrue) )
    
90.     }
    
91.   }
    
92.   datarecord = data.frame( y=datarecord[,"y"] , x=as.factor(datarecord[,"x"]) )
    
93.   y = as.numeric(datarecord$y)
    
94.   Ntotal = length(y)
    
95.   x = as.numeric(datarecord$x)
    
96.   xnames = levels(datarecord$x)
    
97.   NxLvl = length(unique(x))
    
98.   # Construct list of all pairwise comparisons, to compare with NHST TukeyHSD:
    
99.   contrastList = NULL
    
100.   for ( g1idx in 1:(NxLvl-1) ) {
     
101.     for ( g2idx in (g1idx+1):NxLvl ) {
     
102.       cmpVec = rep(0,NxLvl)
     
103.       cmpVec[g1idx] = -1
     
104.       cmpVec[g2idx] = 1
     
105.       contrastList = c( contrastList , list( cmpVec ) )
     
106.     }
     
107.   }
     
108. }
     
109. 
     
110. # Specify the data in a form that is compatible with BRugs model, as a list:
     
111. ySDorig = sd(y)
     
112. yMorig = mean(y)
     
113. z = ( y - yMorig ) / ySDorig
     
114. datalist = list(
     
115.   y = z ,
     
116.   x = x ,
     
117.   Ntotal = Ntotal ,
     
118.   NxLvl = NxLvl
     
119. )
     
120. # Get the data into BRugs:
     
121. modelData( bugsData( datalist ) )
     
122. 
     
123. #------------------------------------------------------------------------------
     
124. # INTIALIZE THE CHAINS.
     
125. 
     
126. # Autocorrelation within chains is large, so use several chains to reduce
     
127. # degree of thinning. But we still have to burn-in all the chains, which takes
     
128. # more time with more chains (on serial CPUs).
     
129. nchain = 5
     
130. modelCompile( numChains = nchain )
     
131. 
     
132. if ( F ) {
     
133.    modelGenInits() # often won't work for diffuse prior
     
134. } else {
     
135.   #  initialization based on data
     
136.   theData = data.frame( y=datalist$y , x=factor(x,labels=xnames) )
     
137.   a0 = mean( theData$y )
     
138.   a = aggregate( theData$y , list( theData$x ) , mean )[,2] - a0
     
139.   ssw = aggregate( theData$y , list( theData$x ) ,
     
140.                   function(x){var(x)*(length(x)-1)} )[,2]
     
141.   sp = sqrt( sum( ssw ) / length( theData$y ) )
     
142.   genInitList <- function() {
     
143.     return(
     
144.         list(
     
145.             a0 = a0 ,
     
146.             a = a ,
     
147.             sigma = sp ,
     
148.             aSDunabs = sd(a)
     
149.         )
     
150.     )
     
151.   }
     
152.   for ( chainIdx in 1 : nchain ) {
     
153.     modelInits( bugsInits( genInitList ) )
     
154.   }
     
155. }
     
156. 
     
157. #------------------------------------------------------------------------------
     
158. # RUN THE CHAINS
     
159. 
     
160. # burn in
     
161. BurnInSteps = 10000
     
162. modelUpdate( BurnInSteps )
     
163. # actual samples
     
164. samplesSet( c( "a0" ,  "a" , "b0" , "b" , "sigma" , "aSD" ) )
     
165. stepsPerChain = ceiling(5000/nchain)
     
166. thinStep = 40
     
167. modelUpdate( stepsPerChain , thin=thinStep )
     
168. 
     
169. #------------------------------------------------------------------------------
     
170. # EXAMINE THE RESULTS
     
171. 
     
172. source("plotChains.R")
     
173. source("plotPost.R")
     
174. 
     
175. checkConvergence = T
     
176. if ( checkConvergence ) {
     
177.   # Merely for insight, take a look at the a0,a[] values:
     
178.   sumInfo = plotChains( "a0" , saveplots=T , filenameroot=fnroot )
     
179.   sumInfo = plotChains( "a" , saveplots=T , filenameroot=fnroot )
     
180.    sumInfo = plotChains( "b0" , saveplots=T , filenameroot=fnroot )
     
181.    sumInfo = plotChains( "b" , saveplots=T , filenameroot=fnroot )
     
182.    readline("Press any key to clear graphics and continue")
     
183.    graphics.off()
     
184.    sumInfo = plotChains( "sigma" , saveplots=T , filenameroot=fnroot )
     
185.    sumInfo = plotChains( "aSD" , saveplots=T , filenameroot=fnroot )
     
186. }
     
187. 
     
188. # Extract and plot the SDs:
     
189. sigmaSample = samplesSample("sigma")
     
190. aSDSample = samplesSample("aSD")
     
191. windows()
     
192. layout( matrix(1:2,nrow=2) )
     
193. par( mar=c(3,1,2.5,0) , mgp=c(2,0.7,0) )
     
194. histInfo = plotPost( sigmaSample , xlab="sigma" , main="Cell SD" , breaks=30 , col="skyblue" )
     
195. histInfo = plotPost( aSDSample , xlab="aSD" , main="a SD" , breaks=30 , col="skyblue" )
     
196. savePlot( file=paste(fnroot,"SD.eps",sep="") , type="eps" )
     
197. 
     
198. readline("Press any key to clear graphics and continue")
     
199. graphics.off()
     
200. 
     
201. # Extract b values:
     
202. b0Sample = samplesSample( "b0" )
     
203. chainLength = length(b0Sample)
     
204. bSample = array( 0 , dim=c( datalist$NxLvl , chainLength ) )
     
205. for ( xidx in 1:datalist$NxLvl ) {
     
206.    bSample[xidx,] = samplesSample( paste("b[",xidx,"]",sep="") )
     
207. }
     
208. # Convert from standardized b values to original scale b values:
     
209. b0Sample = b0Sample * ySDorig + yMorig
     
210. bSample = bSample * ySDorig
     
211. 
     
212. # Plot b values:
     
213. windows(datalist$NxLvl*2.75,2.5)
     
214. layout( matrix( 1:datalist$NxLvl , nrow=1 ) )
     
215. par( mar=c(3,1,2.5,0) , mgp=c(2,0.7,0) )
     
216. for ( xidx in 1:datalist$NxLvl ) {
     
217.     histInfo = plotPost( bSample[xidx,] , breaks=30 , col="skyblue" ,
     
218.               xlab=bquote(beta*1[.(xidx)]) ,
     
219.               main=paste("x:",xnames[xidx])  )
     
220. }
     
221. savePlot( file=paste(fnroot,"b.eps",sep="") , type="eps" )
     
222. 
     
223. # Display contrast analyses
     
224. nContrasts = length( contrastList )
     
225. if ( nContrasts > 0 ) {
     
226.    nPlotPerRow = 5
     
227.    nPlotRow = ceiling(nContrasts/nPlotPerRow)
     
228.    nPlotCol = ceiling(nContrasts/nPlotRow)
     
229.    windows(3.75*nPlotCol,2.5*nPlotRow)
     
230.    layout( matrix(1:(nPlotRow*nPlotCol),nrow=nPlotRow,ncol=nPlotCol,byrow=T) )
     
231.    par( mar=c(4,0.5,2.5,0.5) , mgp=c(2,0.7,0) )
     
232.    for ( cIdx in 1:nContrasts ) {
     
233.        contrast = matrix( contrastList[[cIdx]],nrow=1) # make it a row matrix
     
234.        incIdx = contrast!=0
     
235.        histInfo = plotPost( contrast %*% bSample , compVal=0 , breaks=30 ,
     
236.                 xlab=paste( round(contrast[incIdx],2) , xnames[incIdx] ,
     
237.                             c(rep("+",sum(incIdx)-1),"") , collapse=" " ) ,
     
238.                 cex.lab = 1.0 ,
     
239.                 main=paste( "X Contrast:", names(contrastList)[cIdx] ),
     
240.                 col="skyblue" )
     
241.    }
     
242.    savePlot( file=paste(fnroot,"xContrasts.eps",sep="") , type="eps" )
     
243. }
     
244. 
     
245. #==============================================================================
     
246. # Do NHST ANOVA and t tests:
     
247. 
     
248. theData = data.frame( y=y , x=factor(x,labels=xnames) )
     
249. aovresult = aov( y ~ x , data = theData ) # NHST ANOVA
     
250. cat("\n------------------------------------------------------------------\n\n")
     
251. print( summary( aovresult ) )
     
252. cat("\n------------------------------------------------------------------\n\n")
     
253. print( model.tables( aovresult , "means" ) , digits=4 )
     
254. windows()
     
255. boxplot( y ~ x , data = theData )
     
256. cat("\n------------------------------------------------------------------\n\n")
     
257. print( TukeyHSD( aovresult , "x" , ordered = FALSE ) )
     
258. windows()
     
259. plot( TukeyHSD( aovresult , "x" ) )
     
260. if ( T ) {
     
261.   for ( xIdx1 in 1:(NxLvl-1) ) {
     
262.     for ( xIdx2 in (xIdx1+1):NxLvl ) {
     
263.       cat("\n----------------------------------------------------------\n\n")
     
264.       cat( "xIdx1 = " , xIdx1 , ", xIdx2 = " , xIdx2 ,
     
265.            ", M2-M1 = " , mean(y[x==xIdx2])-mean(y[x==xIdx1]) , "\n" )
     
266.       print( t.test( y[x==xIdx2] , y[x==xIdx1] , var.equal=T ) ) # t test
     
267.     }
     
268.   }
     
269. }
     
270. cat("\n------------------------------------------------------------------\n\n")
     
271. 
     
272. #==============================================================================

复制代码

第二版有吗？

谢谢分享

谢谢分享