library(MASS) rm(list=ls()) ############## First DSfQM Multivariate SPC ##############

mean vector

mean_v <- c(0, 0, 0) mean_v2 <- c(1, 2, 1)

covariance matrix

cov_v <- c(1, 0.9, 0.9,

•        0.9, 1, 0.9,
  

•        0.9, 0.9, 1)
  

cov_m <- matrix(cov_v, nrow=3, byrow=TRUE)

phase1 dataset

set.seed(221013) data_phase1 <- mvrnorm(n = 50, mean_v, cov_m) alpha_CL <- 0.005

phase2 dataset

data_phase2_1 <- mvrnorm(n = 20, mean_v, cov_m) #phase2 first mean vector (0,0,0) data_phase2_2 <- mvrnorm(n = 20, mean_v2, cov_m) #phase2 first mean vector (1,2,1)

math

T2_all <- rep(0, nrow(data_phase1)) for(i in 1:nrow(data_phase1)){

• T2_all[i] = t(data_phase1[i, ]-mean_v) %*% solve(cov_m) %*% (data_phase1[i, ]-mean_v)
  

• }

UCL

alpha_CL <- 0.005 m <- nrow(data_phase1) p <- ncol(data_phase1) UCL_beta <- ((m-1)^2/m) * qbeta((1-alpha_CL), p/2, (m-p-1)/2) UCL_beta <- round(UCL_beta, 3)

phase1 plot

par=(mai= c(1,1,1,1)) plot <- plot(T2_all, ylab=expression(T^2), pch=19, type="o", ylim = c(0, 13)) abline(h=UCL_beta, lty=3)

phase1 + phase2 dataset

data_all <- rbind(data_phase1, data_phase2_1, data_phase2_2) # rbind T2_all <- rep(0, nrow(data_all))

math

for(i in 1:nrow(data_all)){

• T2_all[i] = t(data_all[i, ]-mean_v) %*% solve(cov_m) %*% (data_all[i, ]-mean_v)
  

• }

T2_all [1] 1.7004121 4.9174964 1.8412709 7.4735680 0.9182962 0.9218589 0.5539093 2.8362162 [9] 5.0544994 0.8462186 0.6163710 0.7825867 3.4256396 1.4827785 0.7058302 2.2598263 [17] 1.0455508 1.0944138 3.2912419 7.4165530 0.3804331 4.1048206 3.9767988 0.7068873 [25] 3.2431984 2.5882892 0.5371815 2.4052199 0.3558945 5.1932590 6.1760059 0.9805274 [33] 6.4068551 1.7620251 0.9644944 0.8427350 0.4827468 4.5114542 0.5590087 2.4626906 [41] 1.2253159 0.5823863 6.4274547 5.6381058 0.5222864 0.1675730 9.8272575 5.3796506 [49] 3.2861222 2.1295909 2.1337302 2.6828213 4.9186597 6.3304750 0.5826896 3.0631043 [57] 1.9575197 0.7384016 4.0993319 3.2519713 0.2551689 0.9739668 3.1720690 2.0688032 [65] 2.4367682 0.3599036 0.9055491 1.0974007 1.5174754 1.8486463 16.0563151 4.3709364 [73] 17.9637827 7.8418643 7.2986299 13.1158701 14.9838862 6.7810319 10.0576957 18.6068708 [81] 16.0341166 6.9071602 10.1408299 4.9735880 9.2427243 18.1967559 15.6227545 10.7633004 [89] 1.5567915 14.2399518

phase 1+2 plot

par=(mai= c(1,1,1,1)) plot <- plot(T2_all, ylab=expression(T^2), pch=19, type="o", ylim = c(0, 20)) abline(h=UCL_beta, lty=3)

phase 1+2 plot : length of alarm

length(which(T2_all>UCL_beta)) [1] 9 data_phase2_1 <- vector(mode = "list", length = 200) data_phase2_2 <- vector(mode = "list", length = 200) data_all <- vector(mode = "list", length = 200) T2_all <- vector(mode = "list", length = 200) RL <-vector(mode = "list", length = 200) RL [[1]] NULL

[[2]] NULL

[[3]] NULL

[[4]] NULL

[[5]] NULL

[[6]] NULL

[[7]] NULL

[[8]] NULL

[[9]] NULL

[[10]] NULL

[[11]] NULL

[[12]] NULL

[[13]] NULL

[[14]] NULL

[[15]] NULL

[[16]] NULL

[[17]] NULL

[[18]] NULL

[[19]] NULL

[[20]] NULL

[[21]] NULL

[[22]] NULL

[[23]] NULL

[[24]] NULL

[[25]] NULL

[[26]] NULL

[[27]] NULL

[[28]] NULL

[[29]] NULL

[[30]] NULL

[[31]] NULL

[[32]] NULL

[[33]] NULL

[[34]] NULL

[[35]] NULL

[[36]] NULL

[[37]] NULL

[[38]] NULL

[[39]] NULL

[[40]] NULL

[[41]] NULL

[[42]] NULL

[[43]] NULL

[[44]] NULL

[[45]] NULL

[[46]] NULL

[[47]] NULL

[[48]] NULL

[[49]] NULL

[[50]] NULL

[[51]] NULL

[[52]] NULL

[[53]] NULL

[[54]] NULL

[[55]] NULL

[[56]] NULL

[[57]] NULL

[[58]] NULL

[[59]] NULL

[[60]] NULL

[[61]] NULL

[[62]] NULL

[[63]] NULL

[[64]] NULL

[[65]] NULL

[[66]] NULL

[[67]] NULL

[[68]] NULL

[[69]] NULL

[[70]] NULL

[[71]] NULL

[[72]] NULL

[[73]] NULL

[[74]] NULL

[[75]] NULL

[[76]] NULL

[[77]] NULL

[[78]] NULL

[[79]] NULL

[[80]] NULL

[[81]] NULL

[[82]] NULL

[[83]] NULL

[[84]] NULL

[[85]] NULL

[[86]] NULL

[[87]] NULL

[[88]] NULL

[[89]] NULL

[[90]] NULL

[[91]] NULL

[[92]] NULL

[[93]] NULL

[[94]] NULL

[[95]] NULL

[[96]] NULL

[[97]] NULL

[[98]] NULL

[[99]] NULL

[[100]] NULL

[[101]] NULL

[[102]] NULL

[[103]] NULL

[[104]] NULL

[[105]] NULL

[[106]] NULL

[[107]] NULL

[[108]] NULL

[[109]] NULL

[[110]] NULL

[[111]] NULL

[[112]] NULL

[[113]] NULL

[[114]] NULL

[[115]] NULL

[[116]] NULL

[[117]] NULL

[[118]] NULL

[[119]] NULL

[[120]] NULL

[[121]] NULL

[[122]] NULL

[[123]] NULL

[[124]] NULL

[[125]] NULL

[[126]] NULL

[[127]] NULL

[[128]] NULL

[[129]] NULL

[[130]] NULL

[[131]] NULL

[[132]] NULL

[[133]] NULL

[[134]] NULL

[[135]] NULL

[[136]] NULL

[[137]] NULL

[[138]] NULL

[[139]] NULL

[[140]] NULL

[[141]] NULL

[[142]] NULL

[[143]] NULL

[[144]] NULL

[[145]] NULL

[[146]] NULL

[[147]] NULL

[[148]] NULL

[[149]] NULL

[[150]] NULL

[[151]] NULL

[[152]] NULL

[[153]] NULL

[[154]] NULL

[[155]] NULL

[[156]] NULL

[[157]] NULL

[[158]] NULL

[[159]] NULL

[[160]] NULL

[[161]] NULL

[[162]] NULL

[[163]] NULL

[[164]] NULL

[[165]] NULL

[[166]] NULL

[[167]] NULL

[[168]] NULL

[[169]] NULL

[[170]] NULL

[[171]] NULL

[[172]] NULL

[[173]] NULL

[[174]] NULL

[[175]] NULL

[[176]] NULL

[[177]] NULL

[[178]] NULL

[[179]] NULL

[[180]] NULL

[[181]] NULL

[[182]] NULL

[[183]] NULL

[[184]] NULL

[[185]] NULL

[[186]] NULL

[[187]] NULL

[[188]] NULL

[[189]] NULL

[[190]] NULL

[[191]] NULL

[[192]] NULL

[[193]] NULL

[[194]] NULL

[[195]] NULL

[[196]] NULL

[[197]] NULL

[[198]] NULL

[[199]] NULL

[[200]] NULL

for(i in 1:200){

• data_phase2_1[[i]] <- mvrnorm(n = 20, mean_v, cov_m)
  

• data_phase2_2[[i]] <- mvrnorm(n = 20, mean_v2, cov_m)
  

• data_all[[i]] <- rbind(data_phase1, data_phase2_1[[i]], data_phase2_2[[i]])# rbind
  

• T2_all[[i]] <- rep(0, nrow(data_all[[i]]))
  

• for(j in 1:nrow(data_all[[i]])){
  

•     T2_all[[i]][j] = t(data_all[[i]][j, ]-mean_v) %*% solve(cov_m) %*% (data_all[[i]][j, ]-mean_v)
  

• RL[[i]] <- which(T2_all[[i]]>UCL_beta)
  

• }
  

• }

which(T2_all[[3]] > UCL_beta) [1] 71 76 77 78 79 83 84 86 89 RL <- unlist(RL) RL [1] 74 76 77 78 79 82 89 74 79 80 84 85 87 88 89 71 76 77 78 79 83 84 86 89 71 72 74 77 79 84 [31] 60 71 72 73 77 78 79 80 82 85 88 89 52 74 75 76 80 83 87 90 73 74 75 76 78 79 82 84 85 86 [61] 87 88 90 75 78 79 81 87 88 89 71 74 75 76 80 82 86 87 88 89 90 72 76 80 81 83 86 90 53 71 [91] 72 74 75 77 79 81 82 83 86 87 77 78 82 85 86 71 73 74 76 78 83 86 87 89 72 73 76 77 79 80 [121] 82 88 73 75 76 90 72 75 76 77 80 85 86 87 88 89 90 71 72 73 75 77 79 80 82 83 85 86 88 74 [151] 77 79 81 82 83 87 88 89 76 79 80 81 84 86 88 89 90 72 74 76 78 80 82 83 84 87 88 73 76 78 [181] 80 87 90 71 72 76 77 80 81 88 90 73 77 79 80 83 88 89 73 74 77 78 79 81 83 85 89 72 73 75 [211] 76 79 80 82 83 85 88 89 52 71 76 77 78 80 81 82 83 85 88 90 74 78 83 88 89 72 75 78 80 88 [241] 90 73 75 78 79 82 85 87 72 73 75 76 81 82 84 85 90 71 73 78 81 85 89 72 74 78 79 80 89 90 [271] 73 75 77 78 79 84 86 89 90 71 76 78 80 82 83 90 71 72 75 77 79 80 83 84 86 88 90 72 73 76 [301] 77 79 80 83 85 86 88 89 75 76 77 79 83 84 85 87 88 90 72 74 75 76 77 82 84 86 88 89 54 71 [331] 72 73 74 77 80 83 84 71 75 78 81 82 84 85 90 71 72 74 75 76 79 82 84 85 86 87 88 71 74 75 [361] 76 82 87 89 71 72 75 76 77 82 83 84 89 90 72 76 79 82 83 85 89 72 74 76 78 79 80 84 85 86 [391] 88 90 71 73 75 77 78 79 83 84 86 71 73 77 79 80 81 83 84 85 87 88 89 64 71 73 75 78 81 85 [421] 89 90 75 76 82 85 88 89 74 80 82 83 84 85 86 89 69 73 77 78 81 82 85 86 71 72 73 76 77 78 [451] 80 82 84 85 87 88 90 71 74 75 79 86 71 72 76 77 80 81 83 86 87 90 78 81 82 84 85 87 89 90 [481] 74 76 78 79 83 84 86 87 89 75 78 79 80 81 83 84 87 88 89 90 72 73 77 79 82 83 84 87 88 89 [511] 75 76 79 82 86 89 90 60 67 72 73 74 81 82 84 85 87 88 75 76 77 78 79 80 82 83 71 72 74 77 [541] 79 82 84 87 77 78 79 80 85 86 87 88 73 75 77 81 84 88 89 90 71 74 76 78 81 83 84 86 90 75 [571] 77 78 80 83 84 85 88 90 71 72 73 75 77 78 79 80 81 82 84 76 78 82 83 85 86 87 89 90 71 72 [601] 73 74 75 78 79 80 81 83 89 90 74 75 76 77 78 80 82 83 85 86 89 71 77 78 81 82 84 86 88 74 [631] 75 77 80 83 85 86 87 88 89 69 71 77 84 87 88 90 71 74 76 77 82 85 87 88 90 71 75 77 79 82 [661] 86 88 89 73 74 75 76 77 78 79 87 88 89 72 73 74 78 83 86 90 71 72 75 78 82 83 84 85 90 71 [691] 74 75 76 78 83 84 85 86 87 90 71 74 75 77 80 85 86 87 76 79 80 81 86 88 71 73 76 77 78 79 [721] 81 83 84 85 86 87 88 89 73 78 79 81 85 88 89 90 71 79 80 81 82 83 89 53 71 74 75 79 82 83 [751] 85 86 90 72 74 77 78 86 88 74 83 85 88 89 71 72 73 75 76 79 82 83 87 88 72 73 74 75 77 79 [781] 80 82 83 85 87 88 90 71 72 76 81 82 83 85 58 73 74 75 76 78 82 83 87 89 71 75 76 77 79 81 [811] 84 88 90 73 76 77 82 84 87 90 71 72 74 76 77 80 81 82 83 84 88 90 71 73 75 81 82 85 86 87 [841] 88 89 90 72 73 77 81 82 84 87 89 90 54 72 74 75 80 81 82 71 72 73 75 77 78 80 82 86 87 88 [871] 72 75 76 78 81 82 85 86 90 72 73 81 83 84 86 75 78 80 81 82 84 86 87 64 73 74 81 82 85 86 [901] 89 71 73 74 75 76 77 82 85 86 88 90 71 75 77 78 81 83 87 88 90 75 77 78 79 81 84 85 87 89 [931] 69 71 72 74 75 78 79 81 84 85 88 90 72 73 76 77 79 81 84 89 90 72 73 74 75 76 78 80 82 83 [961] 84 86 87 71 73 74 80 83 86 89 72 74 75 76 80 81 83 84 88 71 73 74 76 77 80 84 85 88 73 75 [991] 77 79 80 81 84 89 72 73 77 78 [ reached getOption("max.print") -- omitted 756 entries ] ARL <- mean(RL) print(ARL) [1] 80.13781