library(qcc)
## Package 'qcc' version 2.7
## Type 'citation("qcc")' for citing this R package in publications.
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(grid)
# Question 2 #
# Set parameter #
X1 <- c(44.5, 44.2, 41.6, 41.5, 45.4, 44.0, 44.2, 43.5, 42.8, 43.7, 45.4, 41.7, 44.9, 48.6)
X2 <- c(45.0, 43.8, 44.8, 43.6, 41.9, 48.6, 44.9, 43.5, 41.7, 44.0, 42.8, 43.6, 44.3, 43.8)
X <- c(X1, X2) # data #
mu0 <- round(44.00,digits = 2) # mean #
lam <- 0.3 # lambda #
L <- 3 # L #
sigma <- 1.8 # sigma #
# Calculate EWMA value #
EWMA <- matrix(0,length(X),1)
rownames(EWMA) <- EWMA
for(i in 1:length(X)){
if(i == 1){
EWMA[i,1] <- lam* X[i] +(1 - lam)*mu0}
else{
EWMA[i,1] <- lam* X[i] +(1 - lam)*EWMA[i-1,1]
}}
# Calculate UCL and LCL #
d <- sigma / sqrt(length(X))
temp <- lam/(1 - lam)
UCL <- round( mu0 + L * d * sqrt(temp) ,digits = 2)
LCL <- round( mu0 - L * d * sqrt(temp) ,digits = 2)
Result <- matrix(0,length(EWMA[,1]),1)
Result[c(which(EWMA > UCL),which(EWMA < LCL)),1] = 3
Result[-c(which(EWMA > UCL),which(EWMA < LCL)),1] = 4
EWMA1 <- cbind(EWMA,Result)
# Trail Chart #
EWMA1 <- EWMA1 %>% as.data.frame()
colnames(EWMA1) <- c('EWMA','Result')
EWMA1$Result <- EWMA1$Result %>% as.factor()
EWMA1$Result <- gsub('3','Out of control',EWMA1$Result)
EWMA1$Result <- gsub('4','In-control',EWMA1$Result)
EWMA1$Sample <- c(1:length(EWMA1[,1]))
Echart <- ggplot(EWMA1)
Echart +
geom_point(data = EWMA1 , aes(x=Sample, y=EWMA, colour = Result )) +
scale_color_manual(values = c('#252525', '#ef3b2c'))+
ggtitle("EWMA Chart \n Trail Chart")+
geom_line(data = EWMA1 , aes(x=Sample, y=EWMA),lty = 2)+
geom_hline(yintercept = c(UCL,LCL,mu0)) +
scale_y_continuous(breaks=c(UCL, mu0, LCL),
labels=c(paste0("UCL = ",UCL),paste0("CL = ",mu0),paste0("LCL = ",LCL)))
# Accroding trail chart, there are seven outlier in this operation #
# So, we need to remove the outlier and re-calculate #
# remove Outlier #
rem = 1
while(rem > 0){
d <- sigma / sqrt(length(EWMA))
temp <- lam/(1 - lam)
UCL <- round(mu0 + L * d * sqrt(temp), digits = 2)
LCL <- round(mu0 - L * d * sqrt(temp), digits = 2)
rem <- c(which(EWMA > UCL),which(EWMA < LCL))
if(length(rem) != 0){
EWMA <- EWMA[-rem,]}
rem <- length(rem)
}
# Trace Chart #
EWMA <- EWMA %>% as.data.frame()
colnames(EWMA) <- c('EWMA')
EWMA$Sample <- c(1:length(EWMA[,1]))
Echart <- ggplot(EWMA)
Echart +
geom_point(data = EWMA , aes(x=Sample, y=EWMA ),colour = 'Black') +
ggtitle("EWMA Chart \n Trace Chart")+
geom_line(data = EWMA , aes(x=Sample, y=EWMA),lty = 2)+
geom_hline(yintercept = c(UCL,LCL,mu0)) +
scale_y_continuous(breaks=c(UCL, mu0, LCL),
labels=c(paste0("UCL = ",UCL),paste0("CL = ",mu0),paste0("LCL = ",LCL)))
# Question 4 #
X1 <- c(68,64,63,61,64,62,68,62,65,61,67,68,70,65,62,64,67,69,69,63)
X1 <- matrix(X1,4,5)
X2 <- c(37,34,36,35,34,37,38,38,35,36,38,39,34,37,35,34,36,39,39,37,35,32,35,34,34)
X2 <- matrix(X2,5,5)
X3 <- c(82,84,84,83,83,86,86,81,87,87,82,86,81,82,87,87,85,83,89,85)
X3 <- matrix(X3,4,5)
X4 <- c(5,8,4,4,7,8,6,9,9,9,9,7,3,2,8,4,9,2,7,6,2,2,8,5,5)
X4 <- matrix(X4,5,5)
X5 <- c(22,26,28,26,24,23,25,24,28,28,28,27,26,23,25,21,21,29,25,29,24,23,27,27,26)
X5 <- matrix(X5,5,5)
range <- function(X){
R = max(X) - min(X)
return(R)
}
X <- rbind(X1,X2,X3,X4,X5) # data #
# Calculating #
X <- X %>% as.data.frame %>% mutate(X_bar = apply(X,1,mean),R = apply(X,1,range))
# Standardize #
Rs <- matrix(0,length(X[,1]),1)
Xs <- matrix(0,length(X[,1]),1)
for(i in 1:4){
Rs[i,1] <- X$R[i]/mean(X$R[1:4])
Xs[i,1] <-(X$X_bar[i] - mean(X$X_bar[1:4]))/mean(X$R[1:4])
}
for(i in 5:9){
Rs[i,1] <- X$R[i]/mean(X$R[5:9])
Xs[i,1] <-(X$X_bar[i] - mean(X$X_bar[5:9]))/mean(X$R[5:9])
}
for(i in 10:13){
Rs[i,1] <- X$R[i]/mean(X$R[10:13])
Xs[i,1] <-(X$X_bar[i] - mean(X$X_bar[10:13]))/mean(X$R[10:13])
}
for(i in 14:18){
Rs[i,1] <- X$R[i]/mean(X$R[14:18])
Xs[i,1] <-(X$X_bar[i] - mean(X$X_bar[14:18]))/mean(X$R[14:18])
}
for(i in 19:23){
Rs[i,1] <- X$R[i]/mean(X$R[19:23])
Xs[i,1] <-(X$X_bar[i] - mean(X$X_bar[19:23]))/mean(X$R[19:23])
}
len <- matrix(1:length(X[,1]),length(X[,1]),1)
X <- cbind(len,X,Rs,Xs)
colnames(X) <- c('num','x1','x2','x3','x4','x5','X_bar','R','Rs','Xs')
# Draw Standardize R Chart #
# n = 5, UCL = D4(2.114), LCL = D4(0), CL = 0 #
UCL <- 2.114
LCL <- 0
mu0 <- 1
Rplot <- ggplot()
Rplot <- Rplot +
theme(plot.background = element_rect(colour = "black", size = 3, linetype = 4,
fill = "white"),
plot.title = element_text(colour = "black", face = "bold", size = 20, vjust = 1),
plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), "inches"))+
geom_point(data = X,aes(x = num,y = Rs)) +
ggtitle("Standardized R Chart") +
xlab("Number") +
geom_line(data = X , aes(x=num, y=Rs),lty = 1) +
geom_hline(yintercept = c(UCL,LCL,mu0)) +
geom_vline(xintercept = c(4.2,9.2,13.2,18.2),linetype="dotted") +
scale_y_continuous(breaks = c(UCL, mu0, LCL),
labels = c(paste0("UCL = ",UCL),paste0("CL = ",mu0),paste0("LCL = ",LCL))) +
scale_x_continuous(breaks = c(1:length(X[,1])),
labels = 1:length(X[,1]))
grob1 <- grobTree(textGrob(c("Format 1","Format 2","Format 3","Format 4","Format 5"), x=c(0.019,0.22,0.41,0.585,0.82),
y=0.9, hjust=0,gp=gpar(col="black", fontsize=8)))
Rplot1 <- Rplot + annotation_custom(grob1)
Rplot1
# Accroding the standardize R chart, there is no outlier #
# Draw Standardize X Chart #
# n = 5, UCL = A2(0.577), LCL = -A2(0.577), CL = 0 #
UCL <- 0.577
LCL <- -0.577
mu0 <- 0
Rplot <- ggplot()
Rplot <- Rplot +
theme(plot.background = element_rect(colour = "black", size = 3, linetype = 4,
fill = "white"),
plot.title = element_text(colour = "black", face = "bold", size = 20, vjust = 1),
plot.margin = unit(c(0.2, 0.2, 0.2, 0.2), "inches"))+
geom_point(data = X,aes(x = num,y = Xs)) +
ggtitle("Standardized X Chart") +
xlab("Number") +
geom_line(data = X , aes(x=num, y=Xs),lty = 1) +
geom_hline(yintercept = c(UCL,LCL,mu0)) +
geom_vline(xintercept = c(4.2,9.2,13.2,18.2),linetype="dotted") +
scale_y_continuous(breaks = c(UCL, mu0, LCL),
labels = c(paste0("UCL = ",UCL),paste0("CL = ",mu0),paste0("LCL = ",LCL))) +
scale_x_continuous(breaks = c(1:length(X[,1])),
labels = 1:length(X[,1]))
grob1 <- grobTree(textGrob(c("Format 1","Format 2","Format 3","Format 4","Format 5"), x=c(0.019,0.22,0.41,0.585,0.82),
y=0.9, hjust=0,gp=gpar(col="black", fontsize=8)))
Rplot1 <- Rplot + annotation_custom(grob1)
Rplot1
# Accroding the standardize X_bar chart, there is no outlier #