################################################################################
# Stat-Tutorial-06-CentralLimitTheorem.r                                       #
#                                                                              #
# This is a tutorial on the central limit theorem.                             #
#                                                                              #
# R. Labouriau                                                                 #
# Last revised: Fall 2020                                                      #
################################################################################
# Copyright © 2019 by Rodrigo Labouriau

# Place to save figures
# OG <- "~/Intro2R/Temp"

# Uniform distribution
# Remark: Note that the expectation of a uniform distributed random variable 
#         is 1/2 = 0.5 and the variance is 1/12 (I did not make this calculation
#         in the lectures, but believe me, please)
n.rep <- 200
X <- numeric(n.rep)
n.observations <- 1
for(i in 1:n.rep){
 x <- runif(n.observations)
 X[i] <- (sqrt(n.observations)*(mean(x) - 0.5)) / sqrt(1/12)
}
# pdf(paste(OG, "Figure-Ch5-100.pdf", sep=""))
# hist(X)
qqnorm(X); qqline(X)
# dev.off()

n.observations <- 500
for(i in 1:n.rep){
 x <- runif(n.observations)
 X[i] <- (sqrt(n.observations)*(mean(x) - 0.5)) / sqrt(1/12)
}

# pdf(paste(OG, "Figure-Ch5-101.pdf", sep=""))
# hist(X)
qqnorm(X); qqline(X)
# dev.off()

# Poisson distribution
n.rep <- 100
X <- numeric(n.rep)
L <- 4              # This will be the intensity or lambda parameter.
n.observations <- 1
for(i in 1:n.rep){
 x <- rpois(n=n.observations, lambda=L)
 X[i] <- (sqrt(n.observations)*(mean(x) - L)) / sqrt(L)
}
# hist(X)
# pdf(paste(OG, "Figure-Ch5-102.pdf", sep=""))
qqnorm(X); qqline(X)
# dev.off()

n.observations <- 10000
for(i in 1:n.rep){
 x <- rpois(n=n.observations, lambda=9)
 X[i] <- (sqrt(n.observations)*(mean(x) - L)) / sqrt(L)
}

# hist(X)
# pdf(paste(OG, "Figure-Ch5-103.pdf", sep=""))
qqnorm(X); qqline(X)
# dev.off()


# Binomial distribution
n.rep <- 100
X <- numeric(n.rep)
n.observations <- 1
for(i in 1:n.rep){
 x <- rbinom(n=n.observations, size=5, prob=0.5)
 X[i] <- (sqrt(n.observations)*(mean(x) - 5/2)) / sqrt(0.5*(1-0.5)*5)
}
hist(X)
qqnorm(X); qqline(X)

n.observations <- 10000
for(i in 1:n.rep){
 x <- rbinom(n=n.observations, size=5, prob=0.5)
 X[i] <- (sqrt(n.observations)*(mean(x) - 5/2)) / sqrt(0.5*(1-0.5)*5)
}

hist(X)
qqnorm(X); qqline(X)

# Binomial with rare events

n.rep <- 100
X <- numeric(n.rep)
n.observations <- 1
for(i in 1:n.rep){
 x <- rbinom(n=n.observations, size=5, prob=0.05)
 X[i] <- (sqrt(n.observations)*(mean(x) - (5*0.05))) / sqrt(0.05*(1-0.05)*5)
}
hist(X)
qqnorm(X); qqline(X)

n.observations <- 100
for(i in 1:n.rep){
 x <- rbinom(n=n.observations, size=5, prob=0.05)
 X[i] <- (sqrt(n.observations)*(mean(x) - (5*0.05))) / sqrt(0.05*(1-0.05)*5)
}

hist(X)
qqnorm(X); qqline(X)

# The end