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



#################################################
# 6.1: Uniform Probability Distribution         #
#################################################
# more details: https://r-coder.com/uniform-distribution-r/

# X ~ U(0, 1), 0 <= x <= 1, f(x) = 1/(b-a)=1
dunif(0.5) # uniform density f(0.5)
punif(0.5) # cumulative probability
qunif(0.6) # quantile
runif(10) # random numbers


# plot U(0, 1) density function 
x <- seq(-1, 2, length = 100)
y <- dunif(x)
plot(x, y, xlab = "x", ylab = "f(x)", main = "Uniform (0, 1)",
     type = "s", col = "blue")


# X ~ U(120, 140), 120 <= x <= 140, f(x) = 1/20
dunif(130, min = 120, max = 140) 
punif(130, min = 120, max = 140) 
qunif(0.6, 120, 140) 
runif(10, 120, 140) 


# plot U(120, 140) density function 
x <- seq(110, 150, length = 100)
y <- dunif(x, 120, 140)
plot(x, y, type = "s", lwd = 3, ylim = c(0, .2), col='blue',
     xlab = 'x', ylab = 'Probability', main = 'Uniform Distribution Plot')


# X ~ U(120, 140), P(120 <= x <= 130)
punif(130, 120, 140) - punif(120, 120, 140)


# X ~ U(120, 140), E(X) = 130, V(X) = 33.33
x <- runif(10000, 120, 140) 
mean(x)
var(x)
hist(x, main = "Random numbers from U(120, 140)")




#################################################
# 6.2: Normal Probability Distribution          #
#################################################

# standard normal, Z ~ N(0, 1)
x <- seq(-5, 5, 0.1)
plot(x, dnorm(x), type = "l", col = "red", lwd = 3)
curve(dnorm, -5, 5)

par(mfrow = c(2,2))
plot(x, dnorm(x), type = "l", col = "blue", 
     main = "Standard Normal Probability Density Function (pdf)")
plot(x, pnorm(x), type = "l", col = "blue", 
     main = "Cumulative Distribution Function (cdf)")
p <- seq(0, 1, 0.1)
plot(p, qnorm(p), type = "l", col = "blue", 
     main = "Quantile Function (qf)")
hist(rnorm(1000), main = "Random Numbers from N(0, 1)")


# Figure 6.4: three normal distributions
x <- seq(-30, 35, 0.1)
plot(x, dnorm(x, mean = -10, sd = 5), type = "l", col = "blue", 
     ylab = "Density", xlim = c(-35, 40),
     main = "Normal Distribution with Same Standard Deviations and 
     Different Means")
points(x, dnorm(x, mean = 0, sd = 5), type = "l")
points(x, dnorm(x, mean = 20, sd = 5), type = "l")
abline(v = c(-10, 0, 20), col = "gray")


# Figure 6.5: two normal distributions
x <- seq(-25, 27, 0.1)
plot(x, dnorm(x, mean = 2, sd = 5), type = "l", col = "blue", 
     ylab = "Density",
     main = "Normal Distribution with Same Mean (mu = 2) and 
     Different Standard Deviations, sigma = 5 and 10.")
points(x, dnorm(x, mean = 2, sd = 10), type = "l")
abline(v = 2, col = "gray")


# Figure 6.6: Area Under the Curve for Any Normal Distribution
x <- seq(-4, 4, 0.1)
y <- dnorm(x)
plot(x, y, type = "l", xlab = "", ylab = "Density", col = "blue", lwd = 2, 
     ylim = c(0, 0.8), xaxt = "n", main = "Area Under the Curve for Any Normal Distribution")
polygon(c(x, -3), c(y, 0), col = "lightblue")

arrows(c(-0.5, -0.5, -0.5), c(0.7, 0.6, 0.5), c(-3, -2, -1), c(0.7, 0.6, 0.5), 
       angle = 10)
arrows(c(0.5, 0.5, 0.5), c(0.7, 0.6, 0.5), c(3, 2, 1), c(0.7, 0.6, 0.5), 
       angle = 10)
text(c(0, 0, 0), c(0.7, 0.6, 0.5), c("99.74%", "95.44%", "68.26%"))
abline(v = c(-3, -2, -1, 1, 2, 3), lty = 1, col = "gray")

lab <- c(expression(mu - 3 * sigma), expression(mu - 2 * sigma), 
         expression(mu - sigma), expression(mu), 
         expression(mu + sigma), expression(mu + 2 * sigma),
         expression(mu + 3  * sigma))
axis(1, at = c(-3:3), labels = lab)


                             
# Figure 6.8: Cumulative Probability for Normal Distribution Corresponding
#             to P(z <= 1.00)
# Z ~ N(0, 1), P(Z <= 1.00)
pnorm(1)

# Figure 6.10: Z ~ N(0, 1), P(-1.00 <= Z <= 1.00)
pnorm(1) - pnorm(-1)

# Figure 6.11: Z ~ N(0, 1), P(Z > 1.58)
1 - pnorm(1.58)

# Figure 6.12: Z ~ N(0, 1), P(Z <= q) = 0.1
qnorm(0.1) # P(Z < -1.281552) = 0.1
qnorm(1 - 0.1)  # P(Z < 1.281552) = 0.9


# page 294
# CONVERTING TO THE STANDARD NORMAL RANDOM VARIABLE
# Grear Tire Company Problem
# x = number of miles the tires will last ~ N(mu, sigma^2)
# mean: mu = 36,500 miles and 
# standard deviation: sigma = 5000
# P(x >= 40000) = ?
x <- 40000
mu <- 36500
sigma <- 5000
z <- (x - mu)/sigma
1 - pnorm(z) # # P(z >= (40000 - 36500) / 5000) 
1 - pnorm(x, mean = mu, sd = sigma)  # P(x >= 40000) 



###############################################################
# 6.3 Normal Approximation of Binomial Probabilities          #
###############################################################

# X ~ B(n , p), making errors
# p = history of making errors in 10% of its invoices. 
# n = 100, A sample of 100 invoices has been taken, 
# compute the probability that 12 invoices contain errors. 
p <- 0.1
n <- 100

(n * p >= 5)
(n * (1 - p) >= 5)

mu <- n * p
sigma <- sqrt(n * p * (1 - p))
x <- 12
dbinom(x, size = n, prob = p)


# P(x = 12) = P(11.5 <= x <= 12.5) = P(0.50 <= z <= 0.83)
z1 <- (12 - 0.5 - mu)/sigma
z2 <- (12 + 0.5 - mu)/sigma
pnorm(z2) - pnorm(z1)

# compute the probability of 13 or fewer errors in the sample of 100 invoices
pbinom(13, size = n, prob = p) # P(X <= 13)
pnorm((13 + 0.5 - mu)/sigma) # If P(X <= n) use P(X < n + 0.5)
pnorm(13 + 0.5, mean = mu, sd = sigma)


#################################################
# 6.4 Exponential Probability Distribution      #
#################################################

# Figure 6.17: Exponential Distribution for the Schips Loading Dock Example
# x: loading time for a truck at the Schips loading dock ~ Exponential Distribution
# x ~ Exp(mu)
# mu = 15, average loading time is 15 minutes
mu <- 15

# note that in R, rate = 1/mean  or rate = 1/mu
# dexp(x, rate = 1, log = FALSE)
# pexp(q, rate = 1, lower.tail = TRUE, log.p = FALSE)
# qexp(p, rate = 1, lower.tail = TRUE, log.p = FALSE)
# rexp(n, rate = 1)

x <- seq(0, 30, 0.1)
plot(x, dexp(x, rate = 1/15), type = "l",  col = "blue", lwd= 2, 
     xlab = "Loading Time", ylab = "f(x)", 
     main = "Exponential Distribution for the Schips Loading Dock Example")
segments(c(6, 18), c(0,0),
         c(6, 18), dexp(c(6, 18), rate = 1/15))
text(10, 0.02, expression(P(paste(6 <= x) <= 18)))


# X ~ Exp(15), P(X <= 6), note that rate = 1/mean
pexp(6, rate = 1/15)

# X ~ Exp(15), P(6 <= X <= 18)
pexp(18, rate = 1/15) - pexp(6, rate = 1/15)



