### Filename: Normal2_2008.R
### Paul Johnson March 28, 2008
### This code should be available somewhere in http://pj.freefaculty.org. If it is not
### email me <pauljohn@ku.edu>
set.seed(234234)
sampleSize <- 500
mymean <- 0
mystd <- 1.5
myx <- seq( mymean - 3*mystd, mymean+ 3*mystd, length.out=sampleSize)
## Creates an empty container to hold samples:
mySamples <- list()
## for loop creates 100 samples and stores them
for (i in 1:100) {
mySamples[[i]] <- rnorm(sampleSize, mean = mymean, sd = mystd)
}
par(mfcol=c(5,2))
for (i in 1:10) {
hist(mySamples[[i]],main=paste("Sample", i) )
}
### force ranges to match
for (i in 1:10) {
mytitle <- paste("Sample", i, "mean", round(mean(mySamples[[i]]),2),"std.dev.", round(sd(mySamples[[i]]),2) )
myhist <- hist(mySamples[[i]],main= mytitle, freq=F, xlim=c(-6,6), xlab="Normally Distributed Variable" )
}
par(mfcol=c(1,1))
### sapply applies the mean function to each sample. The "s" in
### sapply stands for "simplify".
myMeans <- sapply(mySamples, mean)
h1 <- hist(myMeans, freq=FALSE, plot=TRUE, main=paste("Sampling Distribution of ",100," Samples"))
mtext(text=paste("True Mean = ", mymean, "True Standard Deviation = ", mystd) , side=3)
### Note the range of the means:
summary(myMeans)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.14600 -0.05190 -0.00112 -0.00535 0.04040 0.12100
### Compare against the range in one sample, say the 5th
summary(mySamples[[5]])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.640 -1.060 0.014 -0.067 0.874 5.460
text(h1$breaks[1], 5, pos = 4, paste("Mean of Means =", round(mean(myMeans),3),"\n Std. of Means ", round(sd(myMeans),3)), xpd = TRUE)
par(mfcol=c(2,1))
myhist <- hist(mySamples[[1]],main= mytitle, freq=F, xlim=c(-6,6), xlab="Normally Distributed Variable" )
hist(myMeans,freq=F, main=paste("Sampling Distribution of ",100," Samples"), xlim=c(-6,6))
### Another way. write a function to create the data and draw it.
### Use lapply to call it over and over
sampleSize <- 500
par(mfrow=c(5,2))
createDist <-function(i){
z <- rnorm(sampleSize, mean = mymean, sd = mystd)
title <- paste("Run", i, "Normal mean=",mymean," std.dev=",mystd)
hist(z,breaks=20,main=title)
}
createDist(1)
createDist(2)
createDist(3)
createDist(4)
createDist(5)
createDist(6)
createDist(7)
createDist(8)
createDist(9)
createDist(10)
lapply(1:10,createDist);
## [[1]]
## $breaks
## [1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
## [15] 2.0 2.5 3.0 3.5 4.0
##
## $counts
## [1] 2 2 4 5 14 20 28 54 70 53 70 58 37 24 27 14 10 8
##
## $density
## [1] 0.008 0.008 0.016 0.020 0.056 0.080 0.112 0.216 0.280 0.212 0.280
## [12] 0.232 0.148 0.096 0.108 0.056 0.040 0.032
##
## $mids
## [1] -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25
## [12] 0.75 1.25 1.75 2.25 2.75 3.25 3.75
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[2]]
## $breaks
## [1] -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
## [15] 2.5 3.0 3.5 4.0 4.5 5.0
##
## $counts
## [1] 1 5 4 13 25 30 49 68 76 74 45 45 30 12 11 5 4 2 1
##
## $density
## [1] 0.004 0.020 0.016 0.052 0.100 0.120 0.196 0.272 0.304 0.296 0.180
## [12] 0.180 0.120 0.048 0.044 0.020 0.016 0.008 0.004
##
## $mids
## [1] -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75
## [12] 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[3]]
## $breaks
## [1] -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
## [15] 2.5 3.0 3.5 4.0 4.5
##
## $counts
## [1] 1 3 8 19 13 44 34 47 71 81 61 49 33 13 11 10 1 1
##
## $density
## [1] 0.004 0.012 0.032 0.076 0.052 0.176 0.136 0.188 0.284 0.324 0.244
## [12] 0.196 0.132 0.052 0.044 0.040 0.004 0.004
##
## $mids
## [1] -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75
## [12] 1.25 1.75 2.25 2.75 3.25 3.75 4.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[4]]
## $breaks
## [1] -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
## [15] 2.5 3.0 3.5 4.0 4.5
##
## $counts
## [1] 1 3 4 11 21 37 44 65 78 61 59 44 28 26 12 2 2 2
##
## $density
## [1] 0.004 0.012 0.016 0.044 0.084 0.148 0.176 0.260 0.312 0.244 0.236
## [12] 0.176 0.112 0.104 0.048 0.008 0.008 0.008
##
## $mids
## [1] -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75
## [12] 1.25 1.75 2.25 2.75 3.25 3.75 4.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[5]]
## $breaks
## [1] -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5
## [15] 3.0 3.5 4.0 4.5 5.0 5.5
##
## $counts
## [1] 4 5 15 21 31 48 47 71 67 56 62 31 16 14 7 4 0 0 1
##
## $density
## [1] 0.016 0.020 0.060 0.084 0.124 0.192 0.188 0.284 0.268 0.224 0.248
## [12] 0.124 0.064 0.056 0.028 0.016 0.000 0.000 0.004
##
## $mids
## [1] -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75 1.25
## [12] 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[6]]
## $breaks
## [1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
## [15] 2.0 2.5 3.0 3.5 4.0 4.5
##
## $counts
## [1] 1 2 2 5 17 20 29 37 60 71 62 59 53 38 19 13 8 3 1
##
## $density
## [1] 0.004 0.008 0.008 0.020 0.068 0.080 0.116 0.148 0.240 0.284 0.248
## [12] 0.236 0.212 0.152 0.076 0.052 0.032 0.012 0.004
##
## $mids
## [1] -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25
## [12] 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[7]]
## $breaks
## [1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
## [15] 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5
##
## $counts
## [1] 1 2 1 5 10 19 36 38 51 79 67 53 47 38 32 13 5 2 0 0 1
##
## $density
## [1] 0.004 0.008 0.004 0.020 0.040 0.076 0.144 0.152 0.204 0.316 0.268
## [12] 0.212 0.188 0.152 0.128 0.052 0.020 0.008 0.000 0.000 0.004
##
## $mids
## [1] -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25
## [12] 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[8]]
## $breaks
## [1] -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5
## [15] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5
##
## $counts
## [1] 1 0 0 1 7 8 12 18 37 43 57 72 72 60 41 27 24 9 7 1 3
##
## $density
## [1] 0.004 0.000 0.000 0.004 0.028 0.032 0.048 0.072 0.148 0.172 0.228
## [12] 0.288 0.288 0.240 0.164 0.108 0.096 0.036 0.028 0.004 0.012
##
## $mids
## [1] -5.75 -5.25 -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75
## [12] -0.25 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[9]]
## $breaks
## [1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
## [15] 2.0 2.5 3.0 3.5 4.0 4.5 5.0
##
## $counts
## [1] 2 2 3 1 14 29 30 45 57 71 67 56 57 26 19 8 8 3 1 1
##
## $density
## [1] 0.008 0.008 0.012 0.004 0.056 0.116 0.120 0.180 0.228 0.284 0.268
## [12] 0.224 0.228 0.104 0.076 0.032 0.032 0.012 0.004 0.004
##
## $mids
## [1] -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25
## [12] 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
##
## [[10]]
## $breaks
## [1] -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
## [15] 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
##
## $counts
## [1] 2 0 1 5 3 13 29 35 53 52 58 56 63 40 47 24 12 4 2 0 1
##
## $density
## [1] 0.008 0.000 0.004 0.020 0.012 0.052 0.116 0.140 0.212 0.208 0.232
## [12] 0.224 0.252 0.160 0.188 0.096 0.048 0.016 0.008 0.000 0.004
##
## $mids
## [1] -5.25 -4.75 -4.25 -3.75 -3.25 -2.75 -2.25 -1.75 -1.25 -0.75 -0.25
## [12] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75
##
## $xname
## [1] "z"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"
par(mfrow=c(1,1))
### Another way to show a sampling distribution. This uses R's built in "Replicate"
sampleSize <- 10000
repls <- 500
mytitle <- paste("Replications=",repls,"SampleSize=",sampleSize,"Mean=",mymean,"Std.Dev=",mystd)
means <- replicate(repls,mean(rnorm(sampleSize,mean=mymean,sd=mystd)))
myhist <- hist(means, main=mytitle, breaks=15)
# Display some results in text inside the histo
# I know there should be a better way, but don't know it
# Note that in text, \n has a "next line" effect
mylabel <- paste("E(x)=",mymean,"\nobserved
mean=",round(mean(means),2),"\nobserved std.dev \nof means=",round(sd(means),2))
# This puts text into the histogram, guessing
# the positions from the 6'th bin's position
text(myhist$breaks[6],.8*max(myhist$counts),mylabel ,pos=2)
sampleSize <- 10000
repls <- 500
mytitle <- paste("uniforms")
means <- replicate(repls,mean(runif(sampleSize,min=0,max=222)))
par(mfcol=c(5,2))
for (i in 1:10)
hist(runif(sampleSize, min=0, max=222))
##dev.new()
myhist <- hist(means, main=mytitle, breaks=15)
# Display some results in text inside the histo
# I know there should be a better way, but don't know it
# Note that in text, \n has a "next line" effect
mylabel <- paste("E(x)=",mymean,"\nobserved mean=",round(mean(means),2),"\nobserved std.dev \nof means=",round(sd(means),2))
# This puts text into the histogram, guessing
# the positions from the 6'th bin's position
text(myhist$breaks[4], .8*max(myhist$counts), mylabel ,pos=2)
## TODO: superimpose Normal() on that histogram
