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")
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)
m
qnorm(1 - 0.1)  # P(Z < 1.281552) = 0.9
qnorm(0.1) # P(Z < -1.281552) = 0.1
qnorm(0.1) # P(Z < -1.281552) = 0.1
qnorm(1 - 0.1)  # P(Z < 1.281552) = 0.9
x <- 40000
mu <- 36500
sigma <- 5000
x
z <- (x - mu)/sigma
z
sigma
mu
1 - pnorm(x, mean = mu, sd = sigma)  # P(x >= 40000)
z
pnorm(z)
1 - pnorm(z) # # P(z >= (40000 - 36500) / 5000)
1 - pnorm(z) # # P(z >= (40000 - 36500) / 5000)
1 - pnorm(x, mean = mu, sd = sigma)  # P(x >= 40000)
# 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
p
n <- 100
n
n * p
(n * p >= 5)
(n * (1 - p) >= 5)
mu <- n * p
mu
sigma <- sqrt(n * p * (1 - p))
sigma
dbinom(x, size = n, prob = 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)
z1 <- (12 - 0.5 - mu)/sigma
z2 <- (12 + 0.5 - mu)/sigma
pnorm(z2) - pnorm(z1)
pnorm((13 + 0.5 - mu)/sigma) # If P(X <= n) use P(X < n + 0.5)
# 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)
# 20 random numbers from U(0, 10)
runif(20, 0, 10)
0:10
# 5 random numbers from integer 0,...,10 without replacement
sample(0:10, 5)
# 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)
# 7 random numbers from four-digit integer 1000,...,9999 without replacement
sample(1000:9999, 7)
# 7 random numbers from four-digit integer 1000,...,9999 without replacement
sample(1000:9999, 7)
source("C:/Users/hanmi/Downloads/Stat_Rcode_ch02-09_new/Stat_Rcode_ch02-09_new/Chapter07.R", encoding = 'BIG5', echo=TRUE)
# 5 random numbers from integer 0,...,10 with replacement
sample(0:10, 5, replace = T)
# 5 random numbers from integer 0,...,10 with replacement
sample(0:10, 5, replace = T)
# 5 random numbers from integer 0,...,10 with replacement
sample(0:10, 5, replace = T)
c("H", "T")
# Flip a fair coin (H, T) 1000 times
coin <- sample(c("H", "T"), 1000, replace = T)
coin
table(coin)
flip/sum(flip)
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 a fair die (1,..., 6)
sample(1:6, 1)
# 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)
sample(1:6, 1) + sample(1:6, 1)
set.seed(12345)
sample(1:6, 1) + sample(1:6, 1)
set.seed(12345)
sample(1:6, 1) + sample(1:6, 1)
sample(1:6, 1) + sample(1:6, 1)
set.seed(12345)
sample(1:6, 1) + sample(1:6, 1)
sample(1:6, 1) + sample(1:6, 1)
sample(1:6, 1) + sample(1:6, 1)
sample(1:6, 1) + sample(1:6, 1)
# selection subset from a data set
head(iris)
tail(iris)
dim(iris)
nrow(iris)
1:nrow(iris)
id <- sampleA(1:nrow(iris), 10)
id <- sampleA(1:nrow(iris), 10)
id <- sample(1:nrow(iris), 10)
id
iris.subset <- iris[id, ]
iris.subset
# 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
# set.seed(12345)
n <- 30
n
1:N
id <- sample(1:N, n)
id
id <- sample(1:N, n)
EAI.subset <- EAI[id, ]
id <- sample(1:N, n)
id
id <- sample(1:N, n)
EAI.subset <- EAI[id, ]
EAI.subset <- EAI[id, ]
EAI.subset
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")
EAI.subset$"Annual Salary"
# xbar = Sample mean annual salary
mean(EAI.subset$"Annual Salary")
# s = Sample standard deviation for annual salary
sd(EAI.subset$"Annual Salary")
EAI$"Training Program"
EAI$"Training Program" == "Yes"
which(EAI$"Training Program" == "Yes")
length(which(EAI$"Training Program" == "Yes"))
nrow(EAI)
# 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
n <- 30
sample(1:nrow(EAI), n)
sample(1:nrow(EAI), n)
sample(1:nrow(EAI), n)
sample(1:nrow(EAI), n)
EAI.subset <- EAI[sample(1:nrow(EAI), n), ]
EAI.subset
xbar <- mean(EAI.subset$"Annual Salary")
pbar <- length(which(EAI.subset$"Training Program" == "Yes"))/n
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
myfun()
result <- replicate(500, myfun())
myfun()
myfun()
myfun()
myfun()
myfun()
result <- replicate(500, myfun())
result
sample.mean  <- result[1, ]
sample.mean
sample.proportion  <-  result[2, ]
sample.proportion
# 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")
sample.mean
# Expected Value of xbar
mean(sample.mean)
mean(EAI$"Annual Salary")
# Expected Value of xbar
mean(sample.mean)
mean(EAI$"Annual Salary")
# sd of xbar
sd(sample.mean)
sigma <- sd(EAI$"Annual Salary")
sigma <- sd(EAI$"Annual Salary")
sigma/sqrt(n)
sigma.xbar.fin <- sqrt((N - n)/(N - 1)) * (sigma/sqrt(n))
sigma.xbar.fin
sigma.xbar.inf <- sigma/sqrt(n)
sigma.xbar.inf
result[2, ]
# mean and sd of sample.proportion
sample.proportion  <-  result[2, ]
mean(sample.proportion)
sd(sample.proportion)
sample.proportion  <-  result[2, ]
mean(sample.proportion)
sd(sample.proportion)
sd(sample.proportion)
# CENTRAL LIMIT THEOREM: U(0, 1)
par(mfrow = c(1, 4))
x <- seq(-0.5, 1.5, 0.01)
x
plot(x, dunif(x), type = "s", main = "Population Distribution, U(0, 1)")
runif(2)
runif(5)
runif(30)
mean(runif(2))
runif(5)
mean(runif(5))
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
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")
n <- 100
xbar <- 82
sigma <- 20 # sigma Known
n <- 100
xbar <- 82
sigma <- 20 # sigma Known
sigma.xbar <- sigma/sqrt(n)
sigma.xbar
alpha <- 0.05  # 95% significance level
qnorm(alpha/2)
qnorm(alpha/2)
z.alpha <- qnorm(1 - alpha/2)
z.alpha
# Margin of Error
margin.error <- z.alpha * sigma.xbar
margin.error
c(-margin.error, margin.error)
# 100(1 - alpha)%  confidence coefficient for population mean
c(xbar - margin.error, xbar + margin.error)
# if you read excel
# library(readxl)
LloydS <- read_excel("data/chap08/LloydS.xlsx")
LloydS
str(LloydS)
head(LloydS)
tail(LloydS)
N <- length(LloydS$"Amount Spent")
N
n <- 30
id <- sample(1:N, n)
LloydS$"Amount Spent"[id]
xbar <- mean(LloydS$"Amount Spent"[id])
xbar
sigma.xbar <- 20/sqrt(n)
sigma.xbar
alpha <- 0.05  # 95% significance level
qnorm(alpha/2)
z.alpha <- qnorm(1 - alpha/2)
z.alpha
c(xbar - margin.error, xbar + margin.error)
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(10000)
count
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
z.alpha
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
mu
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(10000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:10000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
count <- numeric(10000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:10000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(10000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:10000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(15000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:10000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(15000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:10000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
# suppose true mean mu = mean(LloydS$"Amount Spent")
count <- numeric(15000)
mu <- mean(LloydS$"Amount Spent")
mu
n <- 30
alpha <- 0.05  # 95% significance level
z.alpha <- qnorm(1 - alpha/2)
for(i in 1:15000){
xbar <- mean(LloydS$"Amount Spent"[sample(1:N, n)])
sigma.xbar <- sigma/sqrt(n)
margin.error <- z.alpha * sigma.xbar
lower <- xbar - margin.error
upper <- xbar + margin.error
if(lower <= mu & mu <= upper) count[i] <- 1
}
mean(count)
dt(x = 1.812, df = 10)
pt(q = 1.812, df = 10)
qt(p = 0.05, df =10)
qt(p = 1- 0.05, df =10)
qt(p = 0.05, df =10, lower.tail = FALSE)
rt(5, df =10)
NewBalance <- read_excel("data/chap08/NewBalance.xlsx")
head(NewBalance)
str(NewBalance)
n <- length(NewBalance$NewBalance)
NewBalance
NewBalance
NewBalance$NewBalance
n <- length(NewBalance$NewBalance)
n
alpha <- 0.05
xbar <- mean(NewBalance$NewBalance)
xbar
s <- sd(NewBalance$NewBalance)
s
qt(1-alpha/2, df = n -1)
qt(alpha/2, df = n -1)
E <- qt(1-alpha/2, df = n -1) * s/sqrt(n)
E
round(c(xbar - E, xbar + E))
head(Scheer)
# Example
#Training Time In Days for a Sample of 20 Scheer Industries Employees
Scheer <- read_excel("data/chap08/Scheer.xlsx")
head(Scheer)
str(Scheer)
hist(Scheer$Days, main = "Histogram of Training Times for the Scheer Industries Sample")
# Example
#Training Time In Days for a Sample of 20 Scheer Industries Employees
Scheer <- read_excel("data/chap08/Scheer.xlsx")
hist(Scheer$Days, main = "Histogram of Training Times for the Scheer Industries Sample")
hist(Scheer$Days, main = "Histogram of Training Times for the Scheer Industries Sample")
n <- length(Scheer$Days)
n
alpha <- 0.05
alpha
xbar <- mean(Scheer$Days)
xbar
s <- sd(Scheer$Days)
s
E <- qt(1-alpha/2, df = n -1) * s/sqrt(n)
round(c(xbar - E, xbar + E), 1)
c(xbar - E, xbar + E)
round(c(xbar - E, xbar + E), 1)
margin.error <- 2
margin.error
alpha <- 0.05
z.alpha <- qnorm(1-alpha/2)
z.alpha
sigma <- 9.65
(z.alpha^2 * sigma^2)/ margin.error^2
ceiling((z.alpha^2 * sigma^2)/ margin.error^2)
n <- 900
n
x <- 396
p.bar <- x/n
p.bar
z.value <- qnorm(1-alpha/2)
z.value
E <- z.value * sqrt(p.bar * (1-p.bar)/n)
E
E
c(p.bar - E, p.bar + E)
p.bar <- 0.44
E <- 0.025
alpha <- 0.05
z <- qnorm(1-alpha/2)
z
n <- (z^2 * p.bar * (1-p.bar))/E^2
ceiling(n)
p.bar <- 0.5
n <- (z^2 * p.bar * (1-p.bar))/E^2
ceiling(n)
