# out.width='.6\\linewidth'
hist(c(300,2000), 
     prob = T, 
     xlim = c(-3.2,3.2),
     ylim = c(0,0.45),
     xlab = "",
     main = "Standard Normal")

X <- rnorm(3000000, mean = 0, sd = 1)
lines(density(X), lwd = 2, col = "blue") 
abline(v = 0, lwd = 2, lty = 2, col = "red")


# 
nhanes <- read.delim("http://myweb.uiowa.edu/pbreheny/data/lipids.txt")


# eval=F
## hist(nhanes$LDL)


# out.width='.6\\linewidth'
{
hist(nhanes$LDL, prob = T)
xbar<- mean(nhanes$LDL)
std.dev<- sd(nhanes$LDL)
X <- rnorm(3000000, mean = xbar, sd = std.dev)
lines(density(X), lwd = 2, col = "blue") 
}



# 
xbar <- mean(nhanes$LDL)
std.dev <- sd(nhanes$LDL)


# 
# Using the pnorm function directly
pnorm(160, mean = xbar, sd = std.dev, lower.tail = FALSE)

# Or simply find the z-score first and plug it into the pnorm function
z <- (160 - xbar) / std.dev
pnorm(z, lower.tail = FALSE) # Note that we did not insert a mean/sd argument. R will default these arguments to the standard normal curve, mean  = 0, sd = 1


# 
1 - pnorm(160, mean = xbar, sd = std.dev, lower.tail = TRUE)


# 
sum(nhanes$LDL >= 160) / length(nhanes$LDL)


# 
# Finding the area of the lower tails: P(X < 189) - P(X < 160)
pnorm(189, mean = xbar, sd=std.dev, lower.tail = TRUE) - 
pnorm(160, mean = xbar, sd=std.dev, lower.tail = TRUE)

#OR

# Finding the area of the upper tails: P(X > 160) - P(X > 189)
pnorm(160, mean = xbar, sd=std.dev, lower.tail = FALSE) -
pnorm(189, mean = xbar, sd=std.dev, lower.tail = FALSE)

#OR

# Finding the area outside the desired interval and taking the complement: 1 - (P(X < 160) + P(X > 189))
1 - (pnorm(160, mean = xbar, sd=std.dev, lower.tail = TRUE) +
pnorm(189, mean = xbar, sd=std.dev, lower.tail = FALSE))


# 
qnorm(p = 0.06872828, mean = xbar, sd=std.dev, lower.tail = F)


# out.width='.6\\linewidth'
hist(nhanes$TRG)
abline(v=mean(nhanes$TRG), lty=1, lwd = 2, col="blue")
abline(v=median(nhanes$TRG), lty=2, lwd = 2, col="red")
legend(x = "topright",
       legend = c("Mean","Median"),
       lty = c(1, 2),
       col = c("blue", "red"))



# 
sample10 <- sample(nhanes$TRG, 10)
sample30 <- sample(nhanes$TRG, 30)
sample300 <- sample(nhanes$TRG, 300)

mean(sample10)
mean(sample30)
mean(sample300)


# echo = FALSE
ten <- NULL
thirty <- NULL
threehundred <- NULL
for (i in 1:30000){
  sample10 <- sample(nhanes$TRG, 10)
  sample30 <- sample(nhanes$TRG, 30)
  sample300 <- sample(nhanes$TRG, 300)

  ten[i] <- mean(sample10)
  thirty[i] <- mean(sample30)
  threehundred[i] <- mean(sample300)
}
hist(ten, main = "Means of 10 subject study", breaks = 15, xlab = "")


# echo = F
hist(thirty, main = "Means of 30 subject study", breaks = 15, xlab = "")
hist(threehundred, main = "Means of 300 subject study", breaks = 15, xlab = "")


# echo = F
x <- matrix(c(mean(ten), sd(ten), mean(thirty), sd(thirty), mean(threehundred), sd(threehundred), mean(nhanes$TRG), sd(nhanes$TRG)), nrow = 4, byrow = T)
rownames(x) <- c("sample of ten", "sample of thirty", "sample of three-hundred", "population")
colnames(x) <- c("mean", "sd")
x


# 
mean_LDL <- mean(nhanes$LDL)
sd_LDL <- sd(nhanes$LDL)

pnorm(q = 123, mean = mean_LDL, sd = sd_LDL, lower.tail = F)


# 
mean_LDL <- mean(nhanes$LDL)
sd_LDL <- sd(nhanes$LDL)

pnorm(q = 126, mean = mean_LDL, sd = sd_LDL, lower.tail = T) - pnorm(q = 118, mean = mean_LDL, sd = sd_LDL, lower.tail = T)


# 
mean_LDL <- mean(nhanes$LDL)
se_LDL <- sd(nhanes$LDL) / sqrt(50)

pnorm(q = 123, mean = mean_LDL, sd = se_LDL, lower.tail = F)


# 
mean_LDL <- mean(nhanes$LDL)
se_LDL <- sd(nhanes$LDL) / sqrt(50)

pnorm(q = 126, mean = mean_LDL, sd = se_LDL, lower.tail = T) - pnorm(q = 118, mean = mean_LDL, sd = se_LDL, lower.tail = T)


# 
mean_LDL <- mean(nhanes$LDL)
se_LDL <- sd(nhanes$LDL) / sqrt(50)

qnorm(0.025, mean = mean_LDL, sd = se_LDL)
qnorm(0.975, mean = mean_LDL, sd = se_LDL)