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



######################################
# Section 5.2, 5.3                   #
######################################
DiCarlo <- data.frame(NoSold = 0:5, 
                      Days = c(54, 117, 72, 42, 12, 3))
DiCarlo
N <- sum(DiCarlo$Days)
DiCarlo$Prob.Sold <- DiCarlo$Days/N
DiCarlo   # Table 5.4
attach(DiCarlo)

# Expected Value
mu <- sum(NoSold * Prob.Sold)
mu

# Variance
sigma2 <- sum((NoSold - mu)^2 * Prob.Sold)
sigma2
sqrt(sigma2)



######################################
# Section 5.4                        #
######################################
Automobiles <- matrix(c(21, 30, 24, 9, 2, 0,
                        21, 36, 33, 18, 2, 1,
                        9, 42, 9, 12, 3, 2,
                        3, 9, 6, 3, 5, 0), ncol = 6, byrow = T)
dimnames(Automobiles) <- list(Geneva = 0:3, Saratoga = 0:5)
Automobiles.tb <- as.table(Automobiles)
Automobiles.tb/sum(Automobiles.tb)

margin.table(Automobiles.tb, 1)
margin.table(Automobiles.tb, 2)
N <- sum(Automobiles.tb)
N

Automobiles.df <- as.data.frame(Automobiles.tb)

DailySales <- as.integer(as.character(Automobiles.df$Geneva)) + as.integer(as.character(Automobiles.df$Saratoga))
fs <- tapply(Automobiles.df$Freq, DailySales, sum)/N
s <- as.integer(names(fs))
Table5.9 <- data.frame(s, fs)
Table5.9

# expected value
m <- sum(s * fs)
m

# variance
s2 <- sum( (s - m) ^ 2 * fs)
s2

cal_m_var <- function(x, fx){
  m <- sum(x * fx)
  s2 <-  sum((x - m)^2 * fx)
  list(expected.value = m, variance = s2)
}

Geneva <- margin.table(Automobiles.tb, 1)/N  # Table 5.10
Geneva.stat <- cal_m_var(0:3, Geneva)
Geneva.stat

Saratoga <- margin.table(Automobiles.tb, 2)/N  
Saratoga.stat <-  cal_m_var(0:5, Saratoga)
Saratoga.stat

# covariance = (Var(X+Y) - Var(X) - Var(Y))/2
sigma.xy <- (s2 - Geneva.stat$variance - Saratoga.stat$variance)/2
sigma.xy


# correlation
rho.xy <- sigma.xy/sqrt(Geneva.stat$variance * Saratoga.stat$variance)
rho.xy


######################################
# Section 5.5                        #
######################################

# Martin Clothing Store Problem
# X ~ B(n , p)
Purchase <- c("S", "F")
n <- 3
p <- 0.3

dbinom.tb <- data.frame(x = 0:n, fx = dbinom(0:n, n, p))
dbinom.tb
plot(dbinom.tb, type = "h", lwd = 3)

# X ~ B(10, 0.4), f(3)
dbinom(x = 3, size = 10, prob = 0.4)

# X ~ B(10, 0.4), P( X <= 3)
pbinom(q = 3, size = 10, prob = 0.4)
pbinom(3, 10, 0.4)


# X ~ B(10, 0.4), F(q) = 0.2149908 => q = 3
qbinom(p = 0.2149908, size = 10, prob = 0.4)

# X ~ B(10, 0.4), x1, x2, ..., xn iid from X
rbinom(n = 5, size = 10, prob = 0.4)



######################################
# Section 5.6                        #
######################################

# X ~ Poi(10), f(x), x = 0, 1, 2, ...
mu <- 10
dpois.tb <- data.frame(x = 0:5, dpois = round(dpois(x = 0:5, lambda = mu), 4))
dpois.tb
plot(dpois.tb, type = "h", lwd = 3)

# X ~ Poi(10), F(5) = P(X <= 5)
ppois(q = 5, lambda = 10)

# X ~ Poi(10), f(q) = 0.0378
qpois(p = 0.0378, lambda = 10)

# X ~ Poi(10), x1, x2,..., x7 iid X
rpois(7, lambda = 10)



######################################
# Section 5.7                        #
######################################

# textbook:
# X ~ Hyper(N = 12, r = 5, n = 3), r: white balls

# R:
# X ~ Hyper(m = 5, n = 7, k = 3), m: white balls

m <- 5 # r
n <- 12 - m # (m + n = N)
k <- 3 # n
N <- m + n

# f(0), f(1)     
dhyper(0, m, n, k)
dhyper(1, m, n, k)

# F(x) = P( X <= 1)
phyper(1, m, n, k)

# mean
p <- m/N
k * p

x <- 0:k
fx <- dhyper(x, m, n, k)
mu <- sum(x * fx)
mu

# variance
k * p * (1 - p) * (N - k)/(N - 1)

sum( (x - mu)^2 * fx)

