library(qcc)
## Package 'qcc' version 2.7
## Type 'citation("qcc")' for citing this R package in publications.
# Xbar - R Chart #
# Input data #
x1 = c(6,10,7,8,9,12,16,7,9,10,7,8,16,7,15,9,12,4,8,8,9,12,7,8,14,9,12,10,4,8)
x2 = c(10,15,14,5,6,9,7,12,10,4,8,9,7,10,5,12,12,16,8,4,7,9,8,9,9,5,5,8,10,7)
x3 = c(6,9,12,10,4,9,12,4,8,16,7,8,9,9,10,8,10,7,7,9,16,7,8,5,10,12,14,7,8,9)
x4 = c(8,7,12,12,16,7,8,16,7,7,9,10,7,8,9,14,4,8,10,10,8,9,9,7,8,11,7,8,16,8)
x5 = c(16,15,14,4,4,8,16,6,9,7,8,16,6,9,5,9,7,10,5,9,4,9,12,8,9,10,10,7,9,10)
X = cbind(x1,x2,x3,x4,x5)

# R Chart #
RR = function(a){
R = matrix(0,length(a[,1]),1)
for(i in 1 : length(x1)){
  R[i,1] = max(a[i,]) - min(a[i,])
}
R_bar = mean(R)
R_UCL = 2.114 * R_bar
R_CL = R_bar
R_LCL = 0 * R_bar
for( i in 1:length(R)){
  if(R[i,1] >= R_UCL || R[i,1] <= R_LCL){
    X = X[-i,]; RR(X)
  }}
return(X)}

X = RR(X)
R = matrix(0,length(X[,1]),1)
for(i in 1 : length(x1)){
  R[i,1] = max(X[i,]) - min(X[i,])
}
R_bar = mean(R)
R_UCL = 2.114 * R_bar
R_CL = R_bar
R_LCL = 0 * R_bar

qcc.options(bg.margin = 'white')
qcc.options(bg.figure = 'white')
q.xbar = qcc(data = X, type = 'R')

# Xbar - R Chart #
X_b = matrix(0,length(X[,1]),1)
for(i in 1:length(X[,1])){
  X_b[i,1] = mean(X[i,])
}

XX = function(a){
  X_b = matrix(0,length(a[,1]),1)
  for(i in 1 : length(x1)){
    X_b[i,1] = mean(a[i,])
  }
  X_bar = mean(X_b)
  X_UCL = X_bar + 0.577 * R_bar
  X_CL = X_bar
  X_LCL = X_bar - 0.577 * R_bar
  for( i in 1:length(R)){
    if(X_b[i,1] >= X_UCL || X_b[i,1] <= X_LCL){
      X = X[-i,]; XX(X)
    }}
  return(X)}
X = XX(X)
qcc.options(bg.margin = 'white')
qcc.options(bg.figure = 'white')
q.xbar = qcc(data = X, type = 'xbar',data.name = 'Xbar - R Chart')

# 計算在規格上下限下,不合格率rate為多少 #
# USL規格上限,LSL規格下限 #
# 假設製成為常態下 #
USL = 10+4; LSL = 10-4; s = R_bar/2.326
rate = 1 - pnorm(USL,9.08,s) + pnorm(LSL,9.08,s)
rate
## [1] 0.2274874
# Xbar - S Chart #
# Input data #
x1 = c(6,10,7,8,9,12,16,7,9,10,7,8,16,7,15,9,12,4,8,8,9,12,7,8,14,9,12,10,4,8)
x2 = c(10,15,14,5,6,9,7,12,10,4,8,9,7,10,5,12,12,16,8,4,7,9,8,9,9,5,5,8,10,7)
x3 = c(6,9,12,10,4,9,12,4,8,16,7,8,9,9,10,8,10,7,7,9,16,7,8,5,10,12,14,7,8,9)
x4 = c(8,7,12,12,16,7,8,16,7,7,9,10,7,8,9,14,4,8,10,10,8,9,9,7,8,11,7,8,16,8)
x5 = c(16,15,14,4,4,8,16,6,9,7,8,16,6,9,5,9,7,10,5,9,4,9,12,8,9,10,10,7,9,10)
X = cbind(x1,x2,x3,x4,x5)

# S Chart #
SS = function(a){
  S = matrix(0,length(a[,1]),1)
  for(i in 1 : length(x1)){
    S[i,1] =sd(a[i,])
  }
  S_bar = mean(S)
  S_UCL = 2.089 * S_bar
  S_CL = S_bar
  S_LCL = 0 * S_bar
  for( i in 1:length(S)){
    if(S[i,1] >= S_UCL || S[i,1] <= S_LCL){
      X = X[-i,]; SS(X)
    }}
  return(X)}

X = SS(X)
S = matrix(0,length(X[,1]),1)
for(i in 1 : length(x1)){
  S[i,1] = sd(X[i,])
}
S_bar = mean(S)
S_UCL = 2.089 * S_bar
S_CL = S_bar
S_LCL = 0 * S_bar

qcc.options(bg.margin = 'white')
qcc.options(bg.figure = 'white')
q.xbar = qcc(data = X, type = 'S')

X_b = matrix(0,length(X[,1]),1)
for(i in 1:length(X[,1])){
  X_b[i,1] = mean(X[i,])
}

# Xbar - S chart #
XX = function(a){
  X_b = matrix(0,length(a[,1]),1)
  for(i in 1 : length(x1)){
    X_b[i,1] = mean(a[i,])
  }
  X_bar = mean(X_b)
  X_UCL = X_bar + 1.427 * S_bar
  X_CL = X_bar
  X_LCL = X_bar - 1.427 * S_bar
  for( i in 1:length(R)){
    if(X_b[i,1] >= X_UCL || X_b[i,1] <= X_LCL){
      X = X[-i,]; XX(X)
    }}
  return(X)}

X = XX(X)
qcc.options(bg.margin = 'white')
qcc.options(bg.figure = 'white')
q.xbar = qcc(data = X, type = 'xbar',data.name = 'Xbar - S Chart')

# 假設在製程穩定下,且為常態其分配會為N(9.08,3.154^2) #