# 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)
}