###############################################
## Author: Joshua M. Tebbs
## Date: 27 July 2016
## Update: 4 April 2024
## STAT 110 course notes: R Code Chapter 18
###############################################

# Figure 18.1
# Page 184
# Normal distribution
x = seq(2,11.2,0.01)
pdf = dnorm(x,6.6,1.3)
plot(x,pdf,type="l",xlab="Hours of sleep per night",ylab="",xaxp=c(2.7,10.5,6))
abline(h=0)
# Shading (right)
x = seq(8,11.2,0.01)
y = dnorm(x,6.6,1.3)
polygon(c(8,x,11.2),c(0,y,0),col="lightblue")
points(x=8,y=0,pch=19,cex=1.5)
text(10,0.05,0.1357,cex=1.25)


# Figure 18.2
# Page 186
# Simulate sample proportions from a SRS
n = 1000 # sample size 
B = 10000 # number of samples that we will simulate
p = 0.30 # population proportion 
# This code generates the samples and calculates the sample proportions
binomial.data = rbinom(B,n,p) 
sample.prop = binomial.data/n    
# Make histogram of 10000 sample proportions (one from each sample)
bins = seq(0.24,0.36,0.005)  
hist(sample.prop,breaks=bins,xlab="Sample proportion",ylab="Count",main = "",col="lightblue",xaxp=c(0.24,0.36,12))


# Figure 18.3
# Page 187
# Run code above for Figure 18.2
hist(sample.prop,prob=TRUE,breaks=bins,xlab="Sample proportion",ylab="Density",main = "",col="lightblue",xaxp=c(0.24,0.36,12),ylim=c(0,30))
# Add estimate of population density curve 
lines(sort(sample.prop),dnorm(sort(sample.prop),0.3,sqrt(0.3*0.7/1000)),col="red",lwd=2)


# Figure 18.4
# Page 188
y = seq(5,40,0.1)
pdf = dgamma(y,21.83,scale=0.92)
plot(y,pdf,type="l",xlab="BMI",ylab="Population density curve")
abline(h=0)
bmi = rgamma(25,21.83,scale=0.92) # SRS of size n=25 children
points(x=bmi,y=rep(0,25),pch=19,cex=1)
points(x=mean(bmi),y=0,pch=19,cex=2,col="red")
# Calculate sample mean and sample standard deviation
mean(bmi)
sd(bmi)


# Figure 18.5
# Page 189
# Simulate sample means from a SRS
B = 10000 # number of samples
n = 25 # sample size
# Generate the data
data = matrix(rgamma(n*B,21.83,scale=0.92),nrow=B,ncol=n)
# Calculate sample mean for each sample
sample.mean = apply(data,1,mean)
# Make histogram of 10,000 sample means (one from each sample)
bins = seq(16,24,0.25)  
hist(sample.mean,breaks=bins,xlab="Sample mean",ylab="Count",main="",col="lightblue",xaxp=c(16,24,8))


# Figure 18.6
# Page 190
# Run code above for Figure 18.5
hist(sample.mean,prob=TRUE,breaks=bins,xlab="Sample mean",ylab="Density",main = "",col="lightblue",xaxp=c(16,24,8))
# Add estimate of population density curve 
lines(sort(sample.mean),dnorm(sort(sample.mean),20.1,0.86),col="red",lwd=2)