# Example of a Normal quantile plot of data x to provide a visual
# assessment of its conformity with a normal (data is standardized first)
# 
# The ordered data values are posterior point estimates of the underlying
# quantile function.  So, if you plot the ordered data values (y-axis) against
# the exact theoretical quantiles (x-axis), you get a scatter that should be close
# to a straight line if the data look like a random sample from the theoretical
# distribution.  This function chooses the normal as the theory, to provide a
# graphical/visual assessment of how normal the data appear.
#

# To help with assessing the relevance of sampling variability on just
# "how close" to the normal the data appears, we add (very) approximate
# posterior 95% intervals for the uncertain quantile function at each
# point qqbayes below is an S-Plus funtion.  Download this material and
# save it # in a file: call the file qqbayes.ssc

# Then fire up S-Plus and *load* this function using the S-Plus command:
#                          source("qqbayes.ssc") 
# Then you can use the function with any data. e.g., to see how a random
# sample of size 250 from  a Beta(5,3) distribution compares with a normal
# (a normal with the same mean and variance as the Beta sample), just enter
#                                 qqbayes( rbeta(250,5,3) )
# Use 
#                                 qqbayes( rnorm(n,0,1) )
# several times for a specific sample size n -- to get a "feel" for how much
# sampling variability there is in a sample of size n that DOES come from a normal.
# Repeat with different values of n, e.g., n=30, n=100, n=500, etc.
#
# If you read in the Neural Noise data into a vector noise, simply use
#                                 qqbayes( noise )
# What'd you think? How about the neural signal data? How about the
# exchange rate returns?  etc. 
#

qqbayes<-function(mydata) {
   n<-length(mydata) 
   p<-(1:n)/(n+1)
   x<-(mydata-mean(mydata))/sqrt(var(mydata))
   x<-sort(x)
   z<-qnorm(p)
 
   plot(z,x,xlab="standard normal quantiles",ylab="ordered data")
   abline(0,1,col=2)
  
   s<-1.96*sqrt(p*(1-p)/n)
   pl<-p-s; i<-pl<1&pl>0
   lines(z[i],quantile(x,probs=pl[i]),col=3)
   pl<-p+s; i<-pl<1&pl>0
   lines(z[i],quantile(x,probs=pl[i]),col=3)
   title("Normal quantile plot with approx 95% posterior intervals",cex=0.75)
   }