########################################################################################
## This code is for learning normal distribution
## [EMS] Chapter 5
## created on Jan. 05, 2007 by Leena Choi
## modified on Jan. 13, 2007 by Leena Choi
########################################################################################

########################################################################################
## normal distributions and standard normal distributions
########################################################################################

#library(grDevices)

## normal distribution with mean=mu and sd=sigma
set.seed(191)
mu <- 171.5
sigma <- 6.5
height <- rnorm(500, mean=mu, sd=sigma)
dd <- density(height, adjust = 1.5,)
yrange <- range(dd$y)
hist(height, probability = T, nclass=20, axes=F, ylim=c(0.0025, 0.07), ylab=" ", xlab="Height", main="Distribution of Height")
lines(dd$x, dd$y, col=2, lwd=2)
abline(v=0, lty=1, col=4, lwd=2)
text(182, 0.053, "Normal density curve", col=2)
axis(1)

par(mfrow=c(1,1))
x <- seq(150, 195, by=0.5)
y <- dnorm(x, mean=mu, sd=sigma)
plot(x, y, lwd=2, xlim=c(148, 198), ylim=c(0.0025, 0.07), col=2, type="l", xlab="Height (cm)", ylab="Density", bty="n")
text(mu, 0.066, labels = paste("\m =",mu), pos=1, vfont=c("serif symbol", "plain"), cex=1.2, col=3)
abline(v=mu, col=3, lty=2, lwd=1)

## convert to z-score
z <- (x - mu)/sigma
plot(z, y, lwd=2, xlim=c(-4, 4), ylim=c(0.0025, 0.07), col=2, type="l", xlab="z-score", ylab="Density", bty="n")
text(0, 0.066, labels = paste("\m =",0), pos=1, vfont=c("serif symbol", "plain"), cex=1.2, col=3)
abline(v=0, col=3, lty=2, lwd=1)

## comparison
par(mfrow=c(2,1))
plot(x, y, lwd=2, xlim=c(148, 198), ylim=c(0.0025, 0.07), col=2, type="l", xlab="Height (cm)", ylab="Density", bty="n")
text(mu, 0.075, labels = paste("\m =",mu), pos=1, vfont=c("serif symbol", "plain"), cex=1.2, col=3)
abline(v=mu, col=3, lty=2, lwd=1)

plot(z, y, lwd=2, xlim=c(-4, 4), ylim=c(0.0025, 0.07), col=2, type="l", xlab="z-score", ylab="Density", bty="n")
text(0, 0.075, labels = paste("\m =",0), pos=1, vfont=c("serif symbol", "plain"), cex=1.2, col=3)
abline(v=0, col=3, lty=2, lwd=1)

## z-score ; Standard normal deviate (SND)
par(mfrow=c(1,1))
plot(z, y, lwd=2, xlim=c(-4, 4), ylim=c(0.0025, 0.07), col=2, type="l", xlab="z-score", ylab="Density", bty="n")
abline(v=0, col=3, lty=2, lwd=1)
text(0.1, 0.065, labels = paste("\m =",0), vfont=c("serif symbol", "plain"), cex=1.2, col=3)

u.1s <- 1
l.1s <- -1

segments(u.1s, -0.002, u.1s, 0.0365, col=6, lwd=1)
segments(l.1s, -0.002, l.1s, 0.0365, col=6, lwd=1)
text(1, 0.043, labels = paste("+1","\s", sep=""), pos=4, vfont=c("serif symbol", "plain"), cex=1, col=6)
text(-1, 0.043, labels = paste("-1","\s", sep=""), pos=2, vfont=c("serif symbol", "plain"), cex=1, col=6)

u.2s <- 2
l.2s <- -2

segments(u.2s, -0.002, u.2s, 0.0082, col=5, lwd=1)
segments(l.2s, -0.002, l.2s, 0.0082, col=5, lwd=1)
text(2, 0.009, labels = paste("+2","\s", sep=""), pos=4, vfont=c("serif symbol", "plain"), cex=1, col=5)
text(-2, 0.009, labels = paste("-2","\s", sep=""), pos=2, vfont=c("serif symbol", "plain"), cex=1, col=5)

u.3s <- 3
l.3s <- -3
segments(u.3s, -0.002, u.3s, 0.0006, col=4, lwd=2)
segments(l.3s, -0.002, l.3s, 0.0006, col=4, lwd=2)
text(3, 0.002, labels = paste("+3","\s", sep=""), pos=4, vfont=c("serif symbol", "plain"), cex=1, col=4)
text(-3, 0.002, labels = paste("-3","\s", sep=""), pos=2, vfont=c("serif symbol", "plain"), cex=1, col=4)

## Density Curve
par(mfrow=c(1,1))
plot(x, y, lwd=2, xlim=c(148, 198), ylim=c(0.0025, 0.07), col=2, type="l", xlab="Height (cm)", ylab="Density", bty="n")
abline(v=mu, col=3, lty=2, lwd=1)
text(mu+3, 0.065, labels = paste("\m =",mu), vfont=c("serif symbol", "plain"), cex=1.2, col=3)

u.1s <- mu + sigma
l.1s <- mu - sigma

segments(u.1s, -0.002, u.1s, 0.0365, col=6, lwd=1)
segments(l.1s, -0.002, l.1s, 0.0365, col=6, lwd=1)
text(180, 0.043, labels = paste("\\+-1s"), vfont=c("serif symbol", "plain"), cex=1, col=6)

u.2s <- mu + 2*sigma
l.2s <- mu - 2*sigma

segments(u.2s, -0.002, u.2s, 0.0082, col=5, lwd=1)
segments(l.2s, -0.002, l.2s, 0.0082, col=5, lwd=1)
text(187, 0.009, labels = paste("\\+-2s"), vfont=c("serif symbol", "plain"), cex=1, col=5)

u.3s <- mu + 3*sigma
l.3s <- mu - 3*sigma
segments(u.3s, -0.002, u.3s, 0.0006, col=4, lwd=2)
segments(l.3s, -0.002, l.3s, 0.0006, col=4, lwd=2)
text(194, 0.002, labels = paste("\\+-3s"), vfont=c("serif symbol", "plain"), cex=1, col=4)