##########################################################################
# 統計學 (一)                                                            #
# Anderson et al., Statistics for Business & Economics (14th Edition)    #
# Chapter 7: Sampling and Sampling Distributions                         #
# 吳漢銘 國立政治大學統計學系                                            #
# http://www.hmwu.idv.tw/                                                #
##########################################################################



######################################
# 7.2 Selecting a Sample             #
######################################

# 20 random numbers from U(0, 10)
runif(20, 0, 10)

# 5 random numbers from integer 0,...,10 without replacement
sample(0:10, 5)

# 7 random numbers from four-digit integer 1000,...,9999 without replacement
sample(1000:9999, 7)

# 5 random numbers from integer 0,...,10 with replacement
sample(0:10, 5, replace = T)

# Flip a fair coin (H, T) 1000 times
coin <- sample(c("H", "T"), 1000, replace = T)
coin
table(coin)

flip <- table(coin)
flip/sum(flip)


# Flip a unfair coin (H, T) 10 times with probs 0.7 and 0.3
sample(c("H", "T"), 10, replace = T, prob = c(0.7, 0.3))

# Toss a fair die (1,..., 6) 
sample(1:6, 1)

# Toss two fair dice and compute its sum
set.seed(12345)
sample(1:6, 1) + sample(1:6, 1)

set.seed(12345)
sum(sample(1:6, 2, replace = T))


# selection subset from a data set
head(iris)
tail(iris)
dim(iris)

id <- sample(1:nrow(iris), 10)
iris.subset <- iris[id, ]
dim(iris.subset)
iris.subset



######################################
# 7.3 Point Estimation               #
######################################
# Annual Salary and Training Program Status for a Simple Random Sample of 30 EAI Managers

# install.packages("readxl")
library(readxl)

## if you read excel file
EAI <- read_excel("data/chap07/EAI.xlsx")
str(EAI)
head(EAI)
tail(EAI)
N <- nrow(EAI)
N

## if you read csv file, run the following: 
# EAI <- read.csv("data/chap07/EAI.csv")
# EAI$"Annual Salary" <- EAI$Annual.Salary
# EAI$"Training Program" <- EAI$Training.Program

# set.seed(12345)
n <- 30
id <- sample(1:N, n)
EAI.subset <- EAI[id, ]
dim(EAI.subset)
head(EAI.subset)
names(EAI.subset)

# m = Population mean annual salary
mean(EAI$"Annual Salary")

# s = Population standard deviation for annual salary
sd(EAI$"Annual Salary")

# xbar = Sample mean annual salary
mean(EAI.subset$"Annual Salary")

# s = Sample standard deviation for annual salary
sd(EAI.subset$"Annual Salary")


# Population proportion having completed the management training program
p <- length(which(EAI$"Training Program" == "Yes"))/nrow(EAI)
p

# x: the number of managers in the sample who completed the management training program.
# pbar = Sample proportion having completed the management training program
x <- length(which(EAI.subset$"Training Program" == "Yes"))
x
pbar <- x / nrow(EAI.subset)
pbar 



############################################################
# 7.4 introduction to Sampling Distributions               #
############################################################

n <- 30
myfun <- function(){
  EAI.subset <- EAI[sample(1:nrow(EAI), n), ]
  xbar <- mean(EAI.subset$"Annual Salary") 
  pbar <- length(which(EAI.subset$"Training Program" == "Yes"))/n
  c(xbar, pbar)
}
myfun()
result <- replicate(500, myfun())
sample.mean  <- result[1, ]
sample.proportion  <-  result[2, ]


# FIGURE 7.1: Relative Frequency Histogram of xbar Values 
# from 500 Simple Random Samples of Size 30 Each
hist(sample.mean, xlab = "Values of xbar", ylab = "Relative Frequency")


# FIGURE 7.2: Relative Frequency Histogram of pbar Values 
# from 500 Simple Random Samples of Size 30 Each
hist(sample.proportion, xlab = "Values of pbar", ylab = "Relative Frequency")



############################################################
# 7.5 Sampling Distribution of xbar                        #
############################################################

# Expected Value of xbar
mean(sample.mean)
mean(EAI$"Annual Salary")


# sd of xbar
sd(sample.mean)
sigma <- sd(EAI$"Annual Salary")
sigma/sqrt(n)

# Finite Population
# N <- nrow(EAI)
# n <- 30
sigma.xbar.fin <- sqrt((N - n)/(N - 1)) * (sigma/sqrt(n))
sigma.xbar.fin

# Infinite Population
sigma.xbar.inf <- sigma/sqrt(n)
sigma.xbar.inf
  
  
# mean and sd of sample.proportion
sample.proportion  <-  result[2, ]
mean(sample.proportion)
sd(sample.proportion)     


# CENTRAL LIMIT THEOREM: U(0, 1)
par(mfrow = c(1, 4))
x <- seq(-0.5, 1.5, 0.01)
plot(x, dunif(x), type = "s", main = "Population Distribution, U(0, 1)")
no.rep <- 500
x.num2 <- replicate(no.rep, mean(runif(2)))
x.num5 <- replicate(no.rep, mean(runif(5)))
x.num30 <- replicate(no.rep, mean(runif(30)))
hist(x.num2, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 2")
points(density(x.num2), type = "l", col = "blue")
hist(x.num5, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 5")
points(density(x.num5), type = "l", col = "blue")
hist(x.num30, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 30")
points(density(x.num30), type = "l", col = "blue")


# CENTRAL LIMIT THEOREM: Exp(0, 1)
par(mfrow = c(1, 4))
x <- seq(0, 10, 0.1)
plot(x, dexp(x, rate = 1/2), type = "l", main = "Population Distribution, Exp(2)")
no.rep <- 500
x.num2 <- replicate(no.rep, mean(rexp(2, rate = 1/2)))
x.num5 <- replicate(no.rep, mean(rexp(5, rate = 1/2)))
x.num30 <- replicate(no.rep, mean(rexp(30, rate = 1/2)))
hist(x.num2, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 2")
points(density(x.num2), type = "l", col = "blue")
hist(x.num5, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 5")
points(density(x.num5), type = "l", col = "blue")
hist(x.num30, freq = F, xlab = "xbar", main = "Sampling Distribution of xbar, n = 30")
points(density(x.num30), type = "l", col = "blue")


