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

# Figure 17.1
# Page 170
# Simulate 1000 coin flips
flip = rbinom(1000,1,0.5)
# Create plot of the proportion of tails
# Initialise
proportion.tails = flip*0
for (i in 1:1000){
               proportion.tails[i] = sum(flip[1:i])/i
               }
plot(proportion.tails,ylab="Proportion of tails",xlab="Number of tosses",type="o",pch=16)
abline(h=0.5,lty=2)


# Figure 17.2
# Page 171
# Females (left)
female = rbinom(26630,1,0.0804)
# Create plot of the proportion of tails
# Initialise
proportion.positive = female*0
for (i in 1:26630){
               proportion.positive[i] = sum(female[1:i])/i
               }
plot(proportion.positive,ylab="Proportion of positive female cases",xlab="Number of females tested",type="o",pch=16)
# Males (right)
male= rbinom(7481,1,0.1730)
# Create plot of the proportion of tails
# Initialise
proportion.positive = male*0
for (i in 1:7481){
               proportion.positive[i] = sum(male[1:i])/i
               }
plot(proportion.positive,ylab="Proportion of positive male cases",xlab="Number of males tested",type="o",pch=16)

# Figure 17.3 
# Page 174
# Random numbers between 1 and 100
integers = seq(1:100)
data = sample(integers,100000,replace=T)
count <- function(data, n){ 
               length((which(data == n))) 
               }
counts = integers*0
for(i in 1:100) counts[i] <- count(data,i)
prob = counts/100000
plot(integers,prob,type="h",xlab="",ylab="Proportion of times selected",xaxp=c(0,100,10),ylim=c(0,0.02))
abline(h=0)
abline(h=0.01,lty=2)


# Figure 17.4
# Page 178
# Normal density curve (left)
x = seq(70,180,0.1)
pdf = dnorm(x,125,15)
plot(x,pdf,type="l",xlab="Systolic blood pressure (mm Hg)",ylab="",xaxp=c(80,170,6))
abline(h=0)
# Illustrate LLN using SBP example (right)
# Simulate data
sbp = rnorm(10000,125,16)
# Initialise
sample.mean = sbp*0
for (i in 1:10000){
               sample.mean[i] = mean(sbp[1:i])
               }
plot(sample.mean,ylab="Sample mean",xlab="Number of males sampled",type="o",pch=16)
abline(h=125,lty=2)