# script to shade area under normal distribution curve and plot
# results; data is a vector that contains the sample's data,
# lowlim and uplim are the lower and upper limits for the areas
# to shade, color is the color to use for shading, and overlay =
# TRUE will add the shaded area to the existig plot

shadeArea = function(data, lowlim, uplim, color = "blue", 
                     overlay = FALSE, ...) {
  
  # determine means and standard deviation
  avg = mean(data)
  stddev = sd(data)
  
  # set lower and upper boundaries for plotting along the x-axis
  lower = avg - 4 * stddev
  upper = avg + 4 * stddev
  
  # create 101 points along the x-axis distributed equally
  # between the lower and the upper boundaries
  x = seq(lower, upper, (upper - lower)/100)
  
  # plot the normal distribution curve
  if (overlay == FALSE) {
    plot(x,dnorm(x, mean = avg,sd = stddev), 
       type = "l", bty = "n", xaxs = "i", 
       yaxs = "i", yaxt = "n", ...)
  }
  
  # set size of polygons and plot polygons for shading
  dx = seq(lowlim, uplim, 0.01)
  polygon(c(lowlim, dx, uplim), 
          c(0, dnorm(dx, avg, stddev), 0), 
          border = NA, col = color)
}