A plot example

This example is taken from “Introduction to the R Project for Statistical Computing for Use at the ITC” by David Rossiter, pp.133, http://cran.r-project.org/doc/contrib/Rossiter-RIntro-ITC.pdf

n <- 30; mu <- 180; sd <- 20; bin.width <- sd/3v <- rnorm(n, mu, sd)
hist(v, breaks = seq(mu-5*sd, mu+5*sd, by=bin.width))
points(x, dnorm(x, mu, sd)*(n*bin.width))

### visualise the variability of small random samples
# sample size
n <- 30
# number of plot rows, columns
rows <- 2; cols <- 2;
reps <- rows*cols
# parameters of the normal distribution
mu <- 180; sd <- 20
# set up graphic display
par(mfrow=c(rows, cols))
# number of s.d.'s for histogram displays
dd <- 3.5
# compute bin width from s.d and
# the number of bars for each
bin.width=sd/3
# scale x-axis
x.min <- mu-(sdd*sd); x.max <- mu+(sdd*sd)
# scale y-axis
y.max <- n*0.5*bin.width/sd

# compute and display each graph
for (i in 1:reps)
{
v <- rnorm(n, mu, sd)
hist(v, xlim=c(x.min,x.max),
ylim=c(0, y.max),breaks = seq(mu-5*sd, mu+5*sd, by=bin.width),
main="", xlab=paste("Sample",i)) ;
x <- seq(x.min, x.max, length=120)
# true normal distribution
points(x,dnorm(x, mu, sd)*(n*bin.width),type="l", col="blue", lty=1, lwd=1.8)
# distribution estimated from sample
points(x,dnorm(x, mean(v), sd(v))*(n*bin.width),type="l", col="red", lty=2, lwd=1.8)
# print sample params.
# and Pr(Type I error)
text(x.min, 0.9*y.max, paste("mean:", round(mean(v),2)),pos=4)
text(x.min, 0.8*y.max, paste("sdev:", round(sd(v),2)),pos=4)
text(x.min, 0.7*y.max,paste("Pr(t):", round((t.test(v, mu=mu))$p.value,2)),pos=4)
}

# clean up
par(mfrow=c(1,1))
rm(n, rows, cols, reps, mu, sd, v, i, sdd, bin.width, x.min, x.max, y.max, x)