updated BEST

kdhakal · Feb 8, 2013 · 054b367 · 054b367
1 parent 94628c4
commit 054b367
Show file tree

Hide file tree

Showing 5 changed files with 53 additions and 30 deletions.
diff --git a/Bayesian estimation supersedes the t test/Kruschke2012JEPG.pdf b/Bayesian estimation supersedes the t test/Kruschke2012JEPG.pdf
diff --git a/Bayesian estimation supersedes the t test/src/BEST.R b/Bayesian estimation supersedes the t test/src/BEST.R
@@ -1,4 +1,4 @@
-# Version of May 26, 2012.
+# Version of January 19, 2013.
 # John K. Kruschke  
 # johnkruschke@gmail.com
 # http://www.indiana.edu/~kruschke/BEST/
@@ -11,13 +11,7 @@
 ### ******** SEE FILE BESTexample.R FOR INSTRUCTIONS **************
 ### ***************************************************************
 
-if ( .Platform$OS.type != "windows" ) { 
-  cat( " ** YOU APPEAR TO BE RUNNING ON A NON-WINDOWS OPERATING SYSTEM. **\n")
-  cat( " ** THE GRAPHICAL DISPLAYS MIGHT BE DISTORTED. IF YOU HAVE      **\n")
-  cat( " ** PROBLEMS, PLEASE USE THE SEARCH TERM 'graphics' AT          **\n")
-  cat( " ** http://doingbayesiandataanalysis.blogspot.com/              **\n")
-  windows <- function( ... ) X11( ... ) 
-}
+source("openGraphSaveGraph.R") # graphics functions for Windows, Mac OS, Linux
 
 BESTmcmc = function( y1, y2, numSavedSteps=100000, thinSteps=1, showMCMC=FALSE) { 
   # This function generates an MCMC sample from the posterior distribution.
@@ -97,8 +91,11 @@ BESTmcmc = function( y1, y2, numSavedSteps=100000, thinSteps=1, showMCMC=FALSE)
   #------------------------------------------------------------------------------
   # EXAMINE THE RESULTS
   if ( showMCMC ) {
-    windows()
+    openGraph(width=7,height=7)
     autocorr.plot( codaSamples[[1]] , ask=FALSE )
+    show( gelman.diag( codaSamples ) )
+    effectiveChainLength = effectiveSize( codaSamples ) 
+    show( effectiveChainLength )
   }
 
   # Convert coda-object codaSamples to matrix object for easier handling.
@@ -184,7 +181,7 @@ BESTplot = function( y1 , y2 , mcmcChain , ROPEm=NULL , ROPEsd=NULL ,
   nu = mcmcChain[,"nu"]
   if ( pairsPlot ) {
     # Plot the parameters pairwise, to see correlations:
-    windows()
+    openGraph(width=7,height=7)
     nPtToPlot = 1000
     plotIdx = floor(seq(1,length(mu1),by=length(mu1)/nPtToPlot))
     panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
@@ -204,7 +201,7 @@ BESTplot = function( y1 , y2 , mcmcChain , ROPEm=NULL , ROPEsd=NULL ,
   }
   source("plotPost.R")
   # Set up window and layout:
-  windows(width=6.0,height=8.0)
+  openGraph(width=6.0,height=8.0)
   layout( matrix( c(4,5,7,8,3,1,2,6,9,10) , nrow=5, byrow=FALSE ) )
   par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
 
@@ -235,7 +232,7 @@ BESTplot = function( y1 , y2 , mcmcChain , ROPEm=NULL , ROPEsd=NULL ,
   histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
                              -histBinWd ),
                         seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
-                             histBinWd ) ) )
+                             histBinWd ) , xLim ) )
   histInfo = hist( y1 , plot=FALSE , breaks=histBreaks )
   yPlotVec = histInfo$density 
   yPlotVec[ yPlotVec==0.0 ] = NA
@@ -260,7 +257,7 @@ BESTplot = function( y1 , y2 , mcmcChain , ROPEm=NULL , ROPEsd=NULL ,
   histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
                              -histBinWd ),
                         seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
-                             histBinWd ) ) )
+                             histBinWd ) , xLim ) )
   histInfo = hist( y2 , plot=FALSE , breaks=histBreaks )
   yPlotVec = histInfo$density 
   yPlotVec[ yPlotVec==0.0 ] = NA
@@ -470,7 +467,7 @@ makeData = function( mu1 , sd1 , mu2 , sd2 , nPerGrp ,
   y2 = c(y2,y2Out)
   #
   # Set up window and layout:
-  windows() # Plot the data.
+  openGraph(width=7,height=7) # Plot the data.
   layout(matrix(1:2,nrow=2))
   histInfo = hist( y1 , main="Simulated Data" , col="pink2" , border="white" , 
                    xlim=range(c(y1,y2)) , breaks=30 , prob=TRUE )

diff --git a/Bayesian estimation supersedes the t test/src/BEST1G.R b/Bayesian estimation supersedes the t test/src/BEST1G.R
@@ -1,5 +1,5 @@
 # MODIFIED FROM BEST.R FOR ONE GROUP INSTEAD OF TWO.
-# Version of September 3, 2012.
+# Version of January 19, 2013.
 # John K. Kruschke  
 # johnkruschke@gmail.com
 # http://www.indiana.edu/~kruschke/BEST/
@@ -12,13 +12,7 @@
 ### ******** SEE FILE BEST1Gexample.R FOR INSTRUCTIONS ************
 ### ***************************************************************
 
-if ( .Platform$OS.type != "windows" ) { 
-  cat( " ** YOU APPEAR TO BE RUNNING ON A NON-WINDOWS OPERATING SYSTEM. **\n")
-  cat( " ** THE GRAPHICAL DISPLAYS MIGHT BE DISTORTED. IF YOU HAVE      **\n")
-  cat( " ** PROBLEMS, PLEASE USE THE SEARCH TERM 'graphics' AT          **\n")
-  cat( " ** http://doingbayesiandataanalysis.blogspot.com/              **\n")
-  windows <- function( ... ) X11( ... ) 
-}
+source("openGraphSaveGraph.R") # graphics functions for Windows, Mac OS, Linux
 
 BEST1Gmcmc = function( y , numSavedSteps=100000, thinSteps=1, showMCMC=FALSE) { 
   # This function generates an MCMC sample from the posterior distribution.
@@ -93,8 +87,11 @@ BEST1Gmcmc = function( y , numSavedSteps=100000, thinSteps=1, showMCMC=FALSE) {
   #------------------------------------------------------------------------------
   # EXAMINE THE RESULTS
   if ( showMCMC ) {
-    windows()
+    openGraph(width=7,height=7)
     autocorr.plot( codaSamples[[1]] , ask=FALSE )
+    show( gelman.diag( codaSamples ) )
+    effectiveChainLength = effectiveSize( codaSamples ) 
+    show( effectiveChainLength )
   }
 
   # Convert coda-object codaSamples to matrix object for easier handling.
@@ -131,7 +128,7 @@ BEST1Gplot = function( y , mcmcChain , compValm=0.0 , ROPEm=NULL ,
 
   if ( pairsPlot ) {
     # Plot the parameters pairwise, to see correlations:
-    windows()
+    openGraph(width=7,height=7)
     nPtToPlot = 1000
     plotIdx = floor(seq(1,length(mu),by=length(mu)/nPtToPlot))
     panel.cor = function(x, y, digits=2, prefix="", cex.cor, ...) {
@@ -155,7 +152,7 @@ BEST1Gplot = function( y , mcmcChain , compValm=0.0 , ROPEm=NULL ,
 
   source("plotPost.R")
   # Set up window and layout:
-  windows(width=6.0,height=5.0)
+  openGraph(width=6.0,height=5.0)
   layout( matrix( c(3,3,4,4,5,5, 1,1,1,1,2,2) , nrow=6, ncol=2 , byrow=FALSE ) )
   par( mar=c(3.5,3.5,2.5,0.5) , mgp=c(2.25,0.7,0) )
 
@@ -185,7 +182,7 @@ BEST1Gplot = function( y , mcmcChain , compValm=0.0 , ROPEm=NULL ,
   histBreaks = sort( c( seq( histCenter-histBinWd/2 , min(xVec)-histBinWd/2 ,
                              -histBinWd ),
                         seq( histCenter+histBinWd/2 , max(xVec)+histBinWd/2 ,
-                             histBinWd ) ) )
+                             histBinWd ) , xLim ) )
   histInfo = hist( y , plot=FALSE , breaks=histBreaks )
   yPlotVec = histInfo$density 
   yPlotVec[ yPlotVec==0.0 ] = NA

diff --git a/Bayesian estimation supersedes the t test/src/BESTexample.R b/Bayesian estimation supersedes the t test/src/BESTexample.R
@@ -1,4 +1,4 @@
-# Version of May 31, 2012.
+# Version of January 19, 2013.
 # John K. Kruschke  
 # johnkruschke@gmail.com
 # http://www.indiana.edu/~kruschke/BEST/
@@ -31,6 +31,7 @@
 ###      HDIofMCMC.R
 ###      HDIofICDF.R
 ###      BESTexamplePower.R 
+###      openGraphSaveGraph.R
 ### 6. Make sure that R's working directory is the folder in which those 
 ###    files reside. In RStudio, use menu tabs Tools -> Set Working Directory.
 ###    If working in R, use menu tabs File -> Change Dir.
@@ -66,9 +67,9 @@ postInfo = BESTplot( y1 , y2 , mcmcChain , pairsPlot=TRUE )
 # Show detailed summary info on console:
 show( postInfo ) 
 # You can save the plot(s) using the pull-down menu in the R graphics window,
-# or on many platforms you can use R's saveplot() command:
-# savePlot( file="BESTexample" , type="eps" )
-# savePlot( file="BESTexample" , type="jpeg" )
+# or by using the following:
+# saveGraph( file="BESTexample" , type="eps" )
+# saveGraph( file="BESTexample" , type="jpeg" )
 
 # Save the data and results for future use:
 save( y1, y2, mcmcChain, postInfo, file="BESTexampleMCMC.Rdata" )

diff --git a/Bayesian estimation supersedes the t test/src/openGraphSaveGraph.R b/Bayesian estimation supersedes the t test/src/openGraphSaveGraph.R
@@ -0,0 +1,28 @@
+# openGraphSaveGraph.R
+# John K. Kruschke, 2013
+
+openGraph = function( width=7 , height=7 , mag=1.0 , ... ) {
+  if ( .Platform$OS.type != "windows" ) { # Mac OS, Linux
+    X11( width=width*mag , height=height*mag , type="cairo" , ... )
+  } else { # Windows OS
+    windows( width=width*mag , height=height*mag , ... )
+  }
+}
+
+saveGraph = function( file="saveGraphOutput" , type="pdf" , ... ) {
+  if ( .Platform$OS.type != "windows" ) { # Mac OS, Linux
+    if ( any( type == c("png","jpeg","jpg","tiff","bmp")) ) {
+      sptype = type
+      if ( type == "jpg" ) { sptype = "jpeg" }
+      savePlot( file=paste(file,".",type,sep="") , type=sptype , ... )     
+    }
+    if ( type == "pdf" ) {
+      dev.copy2pdf(file=paste(file,".",type,sep="") , ... )
+    }
+    if ( type == "eps" ) {
+      dev.copy2eps(file=paste(file,".",type,sep="") , ... )
+    }
+  } else { # Windows OS
+    savePlot( file=file , type=type , ... )
+  }
+}