On the E ﬃ cient Monitoring of Multivariate Processes with Unknown Parameters

: Control charts are commonly used tools that deal with monitoring of process parameters in an e ﬃ cient manner. Multivariate control charts are more practical and are of greater importance for timely detection of assignable causes in multiple quality characteristics. This study deals with multivariate memory control charts to address smaller shifts in process mean vector. By adopting a new homogeneous weighting scheme, we have designed an e ﬃ cient structure for multivariate process monitoring. We have also investigated the e ﬀ ect of an estimated variance covariance matrix on the proposed chart by considering di ﬀ erent numbers and sizes of subgroups. We have evaluated the performance of the newly proposed multivariate chart under di ﬀ erent numbers of quality characteristics and varying sample sizes. The performance measures used in this study include average run length, standard deviation run length, extra quadratic loss, and relative average run length. The performance analysis revealed that the proposed control chart outperforms the usual scheme under both known and estimated parameters. An application of the study proposal is also presented using a data set related to Olympic archery, for the monitoring of the location of arrows over the concentric rings on the archery board.


Introduction
Statistical process control (SPC) is a collection of useful tools that handle the variations occurring in process parameters. It helps us to monitor any irregularities happening in the process behavior with respect to parameters of concern such as location and dispersion. It is comprised of many useful tools including a cause-and-effect diagram, check sheet, control chart, histogram, pareto chart, scatter diagram, and flowchart (stratification). A combination of all such tools is formally known as an SPC took-kit. All of these SPC tools have their own specific objectives, for example, the cause-and-effect diagram identifies the possible causes for a problem and carries out sorting of ideas into useful categories, the check sheet is used to collect and analyze data, the pareto chart filters the more significant factors, the scatter diagram shows the relationships among different factors in the process, and the flow chart separates the data gathered from various sources. The control chart gets the highest rank in the SPC tool-kit because of its ability to study/analyze process abnormalities over time. These abnormalities cause variations in the process output that may lead to deterioration of the process quality. These variations may be natural (random) or un-natural (systematic) and control charts have the ability to differentiate between these two types of variations. The natural variations are acceptable to the process, but un-natural variations are of concern to the process engineers, as delay in their identification may result into huge loss to the process owners.
In real process scenarios, we come across different quality characteristics of interest that need our attention in a process. This leads us to simultaneous handling of such variables in a multivariate set up, named multivariate SPC. The multivariate control charts are very important tools in timely identification of any changes in multiple quality characteristics of interest in a process. In the recent decades, multivariate control charts received attention because of advancements in computational facilities used as main barriers in their implementation. The three popular categories of univariate control charts (Shewhart, exponential weighted moving average (EWMA), and cumulative sum (CUSUM)) are equally effective in multivariate setup, namely, multivariate Shewhart (the well know Hotelling's T2, [1,2]); multivariate EWMA, that is, MEWMA ( [3,4]); and multivariate CUSUM, that is, MCUSUM ( [5][6][7]).
The literature on multivariate SPC control charts spans multiple directions that might be of concern in process monitoring such as detections of shifts in mean vector and variance covariance matrix, interpretations of out of control signals, principal components, regression adjustments and so on. The authors of [8] discussed multivariate charts for the monitoring of individual observations. The study of [3] gave a review on multivariate control charts. Other studies [5,6,9] have extended the ideology of the study of [10] towards designing quality control charts with the scope of monitoring two or more related characteristics jointly. The work of [9] suggested the use of simultaneously classical CUSUM control charts for related observations. The authors of [5] recommended two multivariate CUSUM control charts for jointly monitoring of related characteristics. First, the multivariate CUSUM chart converts the related characteristics into scalar quantity, and then formulates the structure of the control chart, while the second CUSUM procedure forms a CUSUM vector directly from the observations. Later, the work of [6] also recommended two multivariate CUSUM control charts, name MC 1 and MC 2 . The input statistic of the MC 1 control chart is based on a quadratic form using a vector of accumulated observations. The MC 2 control chart uses the scalar quantity computed through a quadratic form based on the current observation vector.
The work of [11] extended the idea of [12] for related quality characteristics, and proposed a multivariate exponential weighted moving average (MEWMA) control chart. Some of the modification in multivariate classical control charts can be seen in [13][14][15][16][17][18][19][20][21][22][23] and the references therein. The authors of [24] extended the study of [25] known as a homogeneously weighted moving average (HWMA) control chart for the multivariate setup, and named it the MHWMA control chart. The HWMA and MHWMA control charts assign a specific weight to the current observation, and the remaining weight is equally assigned to the rest of the observations for effective process monitoring.
The remaining part of the article is arranged as follows: Section 2 describes the preliminaries and some existing multivariate control charts available; Section 3 provides the charring structure of proposed multivariate control chart; Section 4 provides performance measures and evaluation of the proposed chart and estimation effects, along with a comparative analysis with some existing charts; Section 5 provides an illustrative example that demonstrates the implementation of the proposal; and finally, Section 6 concludes the study and provides some future recommendations.

Preliminaries and Existing Multivariate Control Charts
Let X be a p dimensional vector X p×1 having multivariate normal distribution with µ (mean vector) and Σ (variance-covariance matrix), written as follows: X ∼ N p (µ, Σ). Here, µ is p × 1 mean vector µ p×1 and Σ is a p × p matrix Σ p×p . For our study purposes, we will use µ 0 for the known mean vector and Σ 0 for the known variance-covariance matrix. The X, µ, Σ are defined, respectively, as follows: · · · · · · . . . · · · σ 1p σ 21 . . .
. . x ipt be the t th sample matrix consisting of x ijt as the i th (i = 1, 2, · · · , n) observation of the j th (j = 1, 2, · · · , p) quality charectersitic on the t th (t = 1, 2, · · · , m) sample. The multivariate data structure can be presented as given in Table A1 (Appendix A).
Let X t and S t be the p dimensional t th sample mean vector and sample variance-covariance matrix (X p×1 and S p×p ), respectively. These are defined as follows: x ikt − X kt x ikt − X kt .
In practice, when µ 0 and Σ 0 are unknown, both are estimated from m in-control Phase I for each sample of sample size n as follows: ∀ j l. The diagonal elements of are means of sample variances associated to j th quality characteristics p and off-diagonal are means of sample covariances, obtained from m samples. It is to be missioned that X refers to the mean vector of a particular sub-group, while X refers to overall the mean vector of all the m sub-groups.
We provide brief details of some common multivariate chats and propose a new control chart (for known and unknown parameters) for the mean vector.

χ 2 and Hotelling's T 2 Charts
The Shewhart multivariate process control chart based on the weighted Mahalanobis distance of the sample mean from the process mean is (cf. [2]) as follows: As X ∼ N p (µ 0 , Σ 0 ), the control charting statistics (1) follows χ 2 t distribution (Chi-square probability distribution with t degrees of freedom), which gives an out-of-control signal, if α,p is the p-dimensional α th upper percentage point of the χ 2 distribution) is a specified control limit. The study of [1] proposed the special case of (1) to monitor sample point as By replacing µ 0 and Σ 0 with X and S, respectively, the expression (1) becomes the following: Under the assumption that the m samples are independent and the joint distributions of the p variables is the multivariate normal, X − X at n = 1, which is called Hotelling's T 2 chart for individuals (cf. [7]). The Hotelling's T 2 chart gives an out-of-control

MCUSUM C Chart
The study of [5] proposed two multivariate CUSUM charts. One of these charts (the superior one) is based on the charting statistics: ∀ k > 0. The MCUSUM chart gives an out-of-control alarm if where h 2 is chosen to achieve a desired value of in-control average run length (ARL).

MCUSUM PR Chart
The work of [6] proposed two invariant multivariate CUSUM charts. Both charts are based on the quadratic forms of the mean vector and can be distinguished with respect to the accumulation (i.e., the sum).
MC 1 Chart: The MC 1 accumulates the X vectors and then generates the quadratic forms, while in MC 2 , initially, the quadratic form is produced for each X, and then these are accumulated at a later stage.
where n t the number of subgroups since the most recent renewal (i.e., zero value) of the CUSUM. As C 1 t /n t can be written as the vector C 1 t /n t represents the difference between the accumulated sample average and the target mean. Therefore, at time t, the multivariate process mean may be estimated as C 1 t /n t + µ 0 . The norm of C 1 t is a measure of the distance between the estimate and the target, mathematically defined as follows: A multivariate control chart can be constructed by defining MC 1 chart as and The multivariate process thorough this multivariate CUSUM chart will be declared as out-of-control if MC 1 t > h 3 , where h 3 is chosen to achieve a desired value of in-control ARL. MC 2 Chart: In the MC 1 chart, if accumulation function is applied after taking the square of the distance of each sample mean from µ 0 , it results in a new structure. As a result, we evaluate D 2 t of the t th sample mean as follows: having a χ 2 with p degrees of freedom (on-target) and a non-central χ 2 (off-target). A one-sided CUSUM structure is given as follows: where MC 2 0 = 0. An out-of-control signal is alarmed if MC 2 t > h 3 , where h 3 is set according to a desired in-control ARL.

MEWMA L Chart
The work of [11] developed the multivariate EWMA (MEWMA L ) chart to quickly detect shifts in the process mean vector µ from the target vector µ 0 , and its statistic is defined as follows: where, initially, Z 0 = µ 0 , λ = diag λ 1 , λ 2 , . . . , λ p of the smoothing constant with 0 < λ j ≤ 1, j = 1, 2, . . ., p (usually in practice λ 1 = λ 2 = . . . = λ p = λ), and I is the identity matrix (cf. [11]). The process gives an out-of-control signal in MEWMA L chart when where h 4 is selected to get the chosen in-control average run length. The term Σ Z t is an asymptotic variance-covariance matrix (as t → ∞ ) defined as 2−λ Σ 0 is the time varying variance-covariance matrix. Symbolically, for the MEWMA L chart for asymptotic and time varying cases, we may write as MEWMA L as and MEWMA L tv , respectively. It is to be noted that, when λ = 1, the MEWMA L chart is converted to Hotelling's T 2 chart. (cf. [11]).

MHWMA A Chart
The study of [24] proposed a multivariate homogeneously weighted moving average (MHWMA A ) chart for the mean vector. The plotting statistic of MHWMA A chart is given as follows: where X t−1 represents the sample average of the previous information up to and including the t − 1 observation, X 0 = µ 0 . For the case of equal λs, the λ will be λ = diag(λ, λ, . . . , λ) in Equation (12).
The chart indicates a signal if where h 5 and λ are set based on the in-control average run length (ARL), and Σ H t refers to the variance-covariance matrix defined as follows: If p = 1, the monitoring statistic in Equation (12) becomes H t = λX t + (1 − λ)X t−1 and the variance of H t in Equation (14) becomes . Further details can be seen in [24].

The Proposed MHWMA p Control Chart for Subgroups
The proposed multivariate homogeneous weighted moving average (MHWMA P ) assigns a certain weight to current sample and the remaining weight is uniformly distributed among the previous samples. The charting statistic of this chart, for subgroups of size n, is defined as follows: This can also be extended as follows: where X t denotes the sample mean vector for t th sample, and X t−1 is the mean vector based on the means of t − 1 preceding samples. Initially, the value of X 0 is set equal to the in-control sample mean vector, that is, X 0 = µ 0 . As mentioned earlier, the quality characteristic X t ∼ N p (µ 0 , Σ 0 ). So, it is a well-known fact that the distribution of the sample mean vector X t is also normal with µ 0 and Σ 0 n , that is, X t ∼ N p µ 0 , Σ 0 n (cf. [27]). The process signals in MHWMA P chart when where h 6 is set according to a pre-specified in-control ARL, and Σ MH t is the variance-covariance matrix at time t defined as follows: If p = 1, we have MH t = λX t + (1 − λ)X t−1 (cf. Equation (15)), and the variance of MH t in Equation (14) becomes . The mean (E(MH t )) and variance (V(MH t )) under in-control situation are derived in the Appendix A.

Estimation of Unknown Parameters for MHWMA P Chart for Subgroups
According to [4], the sampling distribution of charting statistics of a chart should be accounted for variation because it may seriously affect the in-control and out-of-control performance. The denial of this variation may cause a significant increase in the number of false alarms in the case of small sample size data in Phase I. The estimation error is inversely proportional to sample size of Phase I data ( [4]).
As in the previous Section 2.1, we usually assume normality and known parameters such as mean vector and variance-covariance matrix (cf. [24]). However, these parameters, that is, µ 0 and Σ 0 , are not always known in practice, and hence estimation is used, which ultimately affects charts' performance. In this section, we study the performance of the MHWMA chart under the estimation case. When the µ 0 and Σ 0 are unknown, they are usually estimated based on m subgroup (of a fixed size n) from an in-control N p (µ 0 , Σ 0 ), where m > 1 and m(n − 1) > p. The estimators of µ 0 and Σ 0 areμ 0 = X andΣ 0 = S. The study of [27] showed that the in-control estimator of the mean vector In this study, we discuss two cases, namely Case I and Case II: Case I: The values of µ 0 and Σ 0 are known and the control limits for Case I are defined in (19); Case II: The values of µ 0 and Σ 0 are estimated from Phase I asμ = X. andΣ = S (derived in Appendix A-III). On the basis ofμ andΣ, the process gives an out-of-control signal in MHWMA P chart when where h 6 is chosen to achieve a desired in-control ARL performance measure, and Σ MH t is the covariance matrix at time point t defined as follows: The estimated mean (E(MH t )) and estimated variance (V(MH t )) under the in-control situation are derived in the Appendix A.

Performance Evaluations and Comparisons
The current section comprises of a set of performance measures that are based on the run length observations. These performance measures, namely, average run length (ARL), standard deviation of the run length (SDRL), extra quadratic loss (EQL), relative average run length (RARL), and performance comparison index (PCI), are analyzed to examine and compare the detection capability of different charts. The ARL studies the performance of a chart at each shift separately, while EQL, RARL, and PCI deal with the overall shift range. The term δ can be defined as δ = n(µ 1 − µ 0 ) Σ −1 0 (µ 1 − µ 0 ), where µ 1 is out-of-control mean vector (cf. [24]). These measures are defined as follows.
Average Run Length (ARL) is the expected number of points prior to a point indicates an out-of-control signal by the chart. The ARL 0 and ARL 1 are two types of ARL mentioned as the in-control and out-of-control situation of a process, respectively. ARL 0 and ARL 1 are the average number of samples desired to get a false alarm in an in-control situation and a shift in an out-of-control situation, respectively. ARL 1 ≤ ARL 0 , a chart with a lower ARL 1 is considered as an efficient chart relative to another chart (cf. [28]). The best chart acquires a lower number of samples to detect the shift in the minimum ARL value. ARL has been used as a performance measure in many research articles such as [29,30].
Along with ARL, we have also reported standard deviation of run length (SDRL) to assess the behavior of the run length distribution (cf. [31]), which measure of the spread in the run length around its center.
The distribution of δ is uniform with density function (δ max − δ min ) −1 over the entire shift range [δ min , δ max ] (cf. [28,35]). The probability distribution of δ (other than normal distribution) has a partial impact on EQL-based performance of charts (cf. [36]). The best performing chart has a minimum EQL value as compared with any competing chart. A control chart exhibiting superior EQL performance under uniform distribution of shift also offers superior performance under non-uniform cases. Relative Average Run Length (RARL) inspects the overall performance of a chart relative a benchmark chart (cf. [33]). Mathematically, it is defined as RARL = (δ max − δ min ) −1 δ max δ min ARL(δ)/ ARL opt (δ)dδ, where ARL(δ) and ARL opt (δ) are the ARL of a specific chart and the standard charts, respectively, at δ. For a standard chart, the RARL is equal to 1.
Performance Comparison Index (PCI) evaluates the performance of control charts based on the EQL. It is defined as the ratio of EQL value of a specific chart to EQL value of optimum chart (cf. [37]). Analogous with RARL, PCI may also get two values: PCI = 1 and PCI > 1 for best and competing charts, respectively. Mathematically, PCI = EQL/EQL opt . As RARL and PCI have similar interpretations, they may be used instead of each other. One may see more details of these performance measures in [32,37] and the references therein.

Monte Carlo Simulation Procedure
The Monte Carlo simulation technique is used to get am estimated solution of certain mathematical and statistical complixities. The literature on this simulation technique can be seen in [38][39][40], and the references therein. In order to find the values ARLs and SDRLs, a Monte Carlo simulation study based on 50,000 replications is carried out in R language (cf. [41]) by setting ARL 0 ≈ 200; p = 2, 3, 4, 5, and 10; k 1 = k 2 = 0.50; and λ = 0.10. For a suitable number of Monte Carlo simulation iterations in quality control, one may see [40,42,43]. The aforementioned computational methods are quite popular for control charts (cf. [44,45]). The calculation technique is described as follows: • Generate the random samples from multivariable normal distribution and compute the values of the charting statistics Hotelling's T 2 , MCUSUM C , MC 1 , MEWMA L as , MEWMA L tv , and MHWMA P . In order to compute other performance measures, we applied famous numerical integration procedures (such as Simpson/Trapezoidal methods) through these steps. i.
Get EQL by integrating the result (obtained in step i) over δ range using the Simpson or Trapezoidal integration technique. iii.
Take the ratio of the ARLs of a specific to optimum/standard charts and multiply with the inverse of entire range of δ. iv.
Calculate RARL by integrating the output over δ range. v.
Get PCI by taking the ratio of EQLs of a particular chart to the optimum/standard chart.

Performance Analysis and Comparisons
This section provides a comparative analysis of understudy and proposed charts using two methods, namely, (i) point to point shifts approach and (ii) entire shift range approach. The first approach is used to calculate ARL values, while the second is used to obtain the values of EQL, RARL, and PCI measures. The ARL values are obtained for understudy charts and the log(ARL) values are plotted along with δ values for a better comparative view (cf. Figures 1-5). On the basis of the ARL 1 values, we calculated the cumulative EQL, RARL, and PCI values, namely, CEQL, CRARL, and CPCI, respectively, and the results are given in Tables 1-3. On the basis of Tables 1-3, we have also provided the superiority order understudy charts at varying choices of p and δ in Table 4.                 Table 4. Superiority order of multivariate charts at varying p and δ.
In order to study the estimation effect of parameters on ARL 0 and SDRL 0 , the related results are provided in Tables 5 and 6, respectively. For this, we produced the ARL 0 and SDRL 0 values at p (p = 2, 3, 5, 10) n (n = 3, 5, 10, 15) λ (λ = 0.10, 0.20, 0.30, 0.40, 0.50, 0.60); m (m = 20, 30, 50, 100, 200, 300, 500, and ∞). For the suitable choice of subgroup size K, the literature may be recommended in this direction, for example, see [46,47] and the references therein. It is to be noted that the estimated ARL 0 . and SDRL 0 results converge to the results for the known parameter if m = ∞.   The best performing chart will exhibit the lowest ARL curve, which means that a lower number of samples are required for a shift in the detection.The findings are summarized as follows.

CASE I: Analysis under known parameters
It can be observed from the ARL curve for p = 2 (cf. in Figure 1) that MHWMA P and MEWMA L as charts outperform MEWMA L tv , MCUSUM C , MC 1 , and Hotelling's T 2 charts. The MHWMA P chart has smaller ARL values at δ = 0.50, and larger values at δ (δ = 1, 1.5, 2, 2.5, and 3) than the MEWMA L as chart. The MEWMA L tv chart has a smaller ARL value than MCUSUM at δ = 0.5, and vice versa at δ = 1. The MC 1 chart has smaller values than the MCUSUM C and MEWMA L tv charts, respectively, at δ = 1.5 − 3, while the traditional Hotelling's T 2 chart has minimum ARL values at all values of δ.
Overall, it can be observed that the MHWMA P chart outperforms the other charts at small shifts, while the MEWMA L as chart does with the increasing shift. The Hotelling's T 2 chart exhibits the worst performance at the shift range of 0.5 − 3.0, while the other charts exhibit performance in between the best and worst charts. The above performance order based on the cumulative performance measures provided in Tables 1-3 for p = 2 can be observed in Table 4.
It can be observed from the ARL curve for p = 3 given in Figure 2 that the MHWMA P and MEWMA L as charts outperform MEWMA L tv , MCUSUM C , MC 1 , and Hotelling's T 2 charts. The MHWMA P chart has smaller ARL values at δ = 0.50, and larger values at δ = 1, 1.5, 2, 2.5, and 3 than the MEWMA L as chart. The MEWMA L tv chart has a smaller ARL value than MCUSUM C at δ = 0.5. The MC 1 chart has smaller values than the MEWMA L tv and MCUSUM C charts, respectively, at δ = 1 and 1.5 and the MCUSUM C and MEWMA L tv charts, respectively, at δ = 2 − 3. The usual Hotelling's T 2 chart has minimum ARL values at all values of δ.
Overall, it can be observed that the MHWMA P . chart exhibits the best performance compared with the other charts at small shifts, while the MEWMA L as chart does with the increasing shift. The Hotelling's T 2 chart exhibits the worst performance at the entire shift range, while the other charts exhibit performance in between the best and worst charts. The above performance order based on the cumulative performance measures provided in Tables 1-3 for p = 3 can be observed in Table 4.
It can be observed from the ARL curve for p = 4 given in Figure 3 that the MHWMA P and MEWMA L as charts outperform MEWMA L tv , MCUSUM C , MC 1 , and Hotelling's T 2 charts. The MHWMA P chart has smaller ARL values at δ = 0.50 and 1, and larger values at δ = 1.5, 2, 2.5, and 3 than the MEWMA L as chart. The MCUSUM C chart has smaller ARL values than the MEWMA L tv and MC 1 charts, respectively, at δ = 0.5. The MC 1 chart has smaller values than the MCUSUM C and MEWMA L tv charts, respectively, at δ = 1 and the MEWMA L tv and MCUSUM C charts, respectively, at δ = 1.5 − 3. The usual Hotelling's T 2 chart has minimum ARL values at all values of δ.
Overall, it can be observed that the MHWMA P chart exhibits the best performance compared with the other charts at small shifts, while the MEWMA L as chart does with the increasing shift. The Hotelling's T 2 chart exhibits the worst performance at the entire shift range, while the other charts exhibit performance in between the best and worst charts. The above performance order based on the cumulative performance measures provided in Tables 1-3 for p = 4 can be observed in Table 4.
It can be observed from the ARL curve for p = 5 given in Figure 4 that the MHWMA P and MEWMA L as charts outperform MEWMA L tv , MCUSUM C , MC 1 , and Hotelling's T 2 charts. The MHWMA P chart has smaller ARL values at δ = 0.50 and 1, and larger values at δ = 1.5, 2, 2.5, and 3 than the MEWMA L as chart. The MCUSUM C chart has smaller ARL values than the MEWMA L tv and MC 1 charts, respectively, at δ = 0.5. The MEWMA L tv chart has smaller values than the MCUSUM C and MC 1 charts, respectively, at δ = 1 and the MC 1 and MCUSUM C charts, respectively, at δ = 1.5, 3. The usual Hotelling's T 2 chart has minimum ARL values at all values of δ.
Overall, it can be observed that the MHWMA P chart exhibits the best performance compared with the other charts at small shifts, while the MEWMA L as chart does with the increasing shift. The Hotelling's T 2 chart exhibits the worst performance at the entire shift range, while the other charts exhibit performance in between the best and worst charts. The above performance order based on the cumulative performance measures provided in Tables 1-3 for p = 5 can be observed in Table 4.
It can be observed from the ARL curve for p = 10 given in Figure 5 that MHWMA P performs better than the MCUSUM C , MC 1 , and MEWMA L as charts at δ = 0.5, 1, and 1.5, respectively, but shows inferior performance o the MEWMA L as chart at δ = 2 − 3, respectively. The MEWMA L as chart shows better performance than the MEWMA L tv and MCUSUM C charts at δ = 1. The MC 1 chart shows better performance than the MEWMA L as . and MEWMA L tv charts, respectively, at δ = 0.5 and the MEWMA L tv and MCUSUM C charts, respectively, at δ = 1.5, 3. The usual Hotelling's T 2 chart has minimum ARL values at all values of δ.
Overall, it can be observed that the MHWMA P chart exhibits the best performance compared with the other charts at small shifts, while the MEWMA L as chart does with the increasing shift. The Hotelling's T 2 chart exhibits the worst performance at the entire shift range, while the other charts exhibit performance in between the best and worst charts. The above performance order based on the cumulative performance measures provided in Tables 1-3 for p = 10 can be observed in Table 4.
Some specific findings regarding Case I are listed below:

CASE II: Analysis under estimated parameters
In order to examine the estimation effect of unknown parameters, we have observed that the ARL 0 and SDRL 0 values are affected because of the control limits based on estimated values of parameters. The detailed results about the estimation effect of parameters on ARL 0 and SDRL 0 values (at varying choices of p, λ, n, and m), provided in Tables 5 and 6, respectively, are given as follows: • The ARL 0 values are inversely related with the value of p in the case of estimated parameters. This means that, as the number of variables increases in the multivariate set up, the ARL 0 values decrease. For example, for λ = 0.1, m = 20, and n = 3, the ARL 0 = 91.12, 65.51, 48.25, 36.36, and 8.87 at p = 2, 3, 4, 5, and 10, respectively (see Table 5). We may also oberserv behavior of the SDRL 0 = 115.34, 83.32, 64.08, 51.48, and 19.13 at p = 2, 3, 4, 5, and 10, respectively (see Table 6).  Table 5). Similarly, we may observe the SDRL 0 results at n = 3 and m = 50, where the SDRL 0 = 131. 47, 163.29, 198.09, 212.99, 222.16, and 229.72 at λ = 0.10, 0.20, 0.30, 0.40, 0.50, and 0.60, respectively (see Table 6).

•
As the sample size n and/or sub-group size m increase, the false alarms produced by MHWMA P chart decrease, which means that the estimated version of ARL 0 converges to 200 as m → ∞ . For instance, for λ = 0.1, as n increases from 3 to 15, the ARL 0 increases from 110.26 to 169.61.72 at m = 30; as m increases from 30 to 500, the ARL 0 increases from 110.26 to 187.72 at n = 3 (see Table 5). Similar behavior of the SDRL 0 results can be observed as λ = 0.1, because, if n increases from 3 to 15, the SDRL 0 increases from 122.86 to 145.34 at n = 3; if m increases from 30 to 500, the SDRL 0 increases from 122.86 to 155.36 at n = 3 (see Table 6).
The Corrected Control Limits: As we observed in Case II above, ARL 0 values are affected owing to the estimation effects in our control limits. Consequently, it is necessary to revise the control limits' coefficient in order to attain the desired ARL 0 . For our study purposes, we have worked out the corrected control limits coefficients for our proposed chart at various choices of p (p = 2, 3, 5, 10), n (n = 3, 5, 10, 15), λ (λ = 0.10, 0.20), and m (m = 20, 30, 50, 100, 200, 300, 500, and ∞). The resulting coefficients are provided in Table 7 at some useful combinations of p, n, λ and m. Moreover, we have also produced some selective graphs of the control limits' coefficient h 6 versus m at various choices of p, n, and λ. The resulting graphs are provided in Figures 6 and 7. It may observed from these tables and figures that the corrected limits get closer to the known limit (limit of the case of known parameter) with the increase in either n or m or both of these quantities.  Table 5). Similar behavior of the 0 results can be observed as = 0.1, because, if increases from 3 to 15, the 0 increases from 122.86 to 145.34 at = 3; if increases from 30 to 500, the 0 increases from 122.86 to 155.36 at = 3 (see Table 6). The Corrected Control Limits: As we observed in Case II above, 0 values are affected owing to the estimation effects in our control limits. Consequently, it is necessary to revise the control limits' coefficient in order to attain the desired 0 . For our study purposes, we have worked out the corrected control limits coefficients for our proposed chart at various choices of ( = 2, 3, 5, 10), ( = 3, 5, 10, 15) , ( = 0.10, 0.20 ), and ( = 20, 30, 50, 100, 200, 300, 500, and ∞ ). The resulting coefficients are provided in Table 7 at some useful combinations of , , . Moreover, we have also produced some selective graphs of the control limits' coefficient versus at various choices of , , and . The resulting graphs are provided in Figures 6 and 7. It may observed from these tables and figures that the corrected limits get closer to the known limit (limit of the case of known parameter) with the increase in either or or both of these quantities.   Case when the parent distribution is heavy tailed: Throughout the paper so far, we have assumed that the variables of interest follow a multivariate normal distribution. This subsection provides the effect of deviation from normality in terms of a heavy tailed distribution, that is, multivariate t-distribution. The in-control ARL values for the proposed chart under known parameters case are given in Table 8, where the data are generated from multivariate t-distribution with v degrees of freedom. The table clearly shows that the ARL values of the proposed chart get close to 200 as we increase the value of v, and finally v → ∞ implies ARL → 200 . This is in line with the multivariate distributional theory that the multivariate T-distribution tends to multivariate normal distribution as v → ∞ . For the smaller values of v, the ARLs are substantially smaller than 200, implying that the constants evaluated in Table 7 are not appropriate.

An Application
Illustrating the proposed control charts using a real dataset is common practice; see, for example, [48]. For illustration purpose in the current study, we used a data set related to Olympic archery, for the monitoring of location of arrows over the concentric rings on the archery board (cf. [7]). Target archery is the most standard form of archery directed by the World Archery Federation, in which participants shoot at stationary circular targets at varying distances. The size of the target varies and is usually measured as a 122 cm face for a distance of 70 m, all containing 10 concentric rings, representing different sectors to score. The outermost two rings (rings 1 and 2) are white, rings 3 and 4 are black, rings 5 and 6 are blue, rings 7 and 8 are red, and rings 9 and 10 (the innermost) are gold (cf. Figure 8a).

An Application
Illustrating the proposed control charts using a real dataset is common practice; see, for example, [48]. For illustration purpose in the current study, we used a data set related to Olympic archery, for the monitoring of location of arrows over the concentric rings on the archery board (cf. [7]). Target archery is the most standard form of archery directed by the World Archery Federation, in which participants shoot at stationary circular targets at varying distances. The size of the target varies and is usually measured as a 122 cm face for a distance of 70 m, all containing 10 concentric rings, representing different sectors to score. The outermost two rings (rings 1 and 2) are white, rings 3 and 4 are black, rings 5 and 6 are blue, rings 7 and 8 are red, and rings 9 and 10 (the innermost) are gold (cf. Figure 8a). The scoring criterion is as follows: we add up all the points based on the arrows hitting different rings as mentioned above. The innermost ring gets a score of 10 and the outermost ring gets a score of 1, and for the inner rings, the score gradually decrease from 10 to 1 (cf. Figure 8b). The first 72 shots in finishes of three arrows of the ranking round of a particular archer with 24 subgroups were collected when the process was working under in-control state (i.e., Phase I) in order to estimate the unknown parameters. We used these estimated values of parameters and computed trail control limits. By plotting the Phase I data on these trail control limits, we found the first six points out of the control limits (which means that the process was not in-control when these samples were collected). Hence, we discarded first six samples from our Phase I data and revised our estimates to compute revised control limits. Now, we have 18 subgroups containing 54 shoots of the elimination round in Phase I, and another 18 additional subgroups consisting of 54 shoots of the elimination round with the same subgroup size, which were collected in Phase II. Setting in-control 0 = 200 under Hotelling's 2 , the MEWMA (MEWMA ), and MHWMA charts are made with ℎ 0 = 12.20, ℎ 4 = 11.21, ℎ 6 = 11.60, and = 0.25, respectively. In order to demonstrate the application of study, we plotted the values of charting statistics (calculated from Phase I and Phase II data) of the Hotelling's 2 , MEWMA , and MHWMA charts against their respective control limits in Figure 9. The scoring criterion is as follows: we add up all the points based on the arrows hitting different rings as mentioned above. The innermost ring gets a score of 10 and the outermost ring gets a score of 1, and for the inner rings, the score gradually decrease from 10 to 1 (cf. Figure 8b). The first 72 shots in finishes of three arrows of the ranking round of a particular archer with 24 subgroups were collected when the process was working under in-control state (i.e., Phase I) in order to estimate the unknown parameters. We used these estimated values of parameters and computed trail control limits. By plotting the Phase I data on these trail control limits, we found the first six points out of the control limits (which means that the process was not in-control when these samples were collected). Hence, we discarded first six samples from our Phase I data and revised our estimates to compute revised control limits. Now, we have 18 subgroups containing 54 shoots of the elimination round in Phase I, and another 18 additional subgroups consisting of 54 shoots of the elimination round with the same subgroup size, which were collected in Phase II. Setting in-control ARL 0 = 200 under Hotelling's T 2 , the MEWMA L as (MEWMA as ), and MHWMA P charts are made with h 0 = 12.20, h 4 = 11.21, h 6 = 11.60, and λ = 0.25, respectively. In order to demonstrate the application of study, we plotted the values of charting statistics (calculated from Phase I and Phase II data) of the Hotelling's T 2 , MEWMA L as , and MHWMA P charts against their respective control limits in Figure 9. From Figure 9, it can be observed that the Hotelling's 2 , MEWMA , and MHWMA charts are detecting an upward and downward shift in the process. The Hotelling's 2 chart indicates the first out-of-control signal at the last sample, while the MEWMA chart triggers the first out-of-control signal at the 24th sample. The proposed MHWMA chart detects the first out-of-control signal at the 24th sample (cf. Figure 9). Moreover, the Hotelling's 2 , MEWMA , and proposed MHWMA charts detect 1, 3, and 7 out-of-control signals (cf. Figure 9), respectively, which means that the proposed MHWMA chart has better detection ability than the other charts. These results display the supremacy of the proposed chart against the competing charts, and these results are also consistent with the findings of the study.
The property of better detection ability of MHWMAp chart claimed in Section 5 has a strong connection with the findings of Section 4, where the superiority of the proposal is established based on RL properties. These properties are based on the repeated iterations from a theoretical model by fixing ARL0 at a specific level for all of the charts under investigation, in order to maintain uniformity in performance evaluation. The RL analysis is also carried out at varying shifts, and the resulting outcomes are presented in tabular and graphical forms (cf. Tables 1-4 and Figures 1-5), where the study proposal exhibits a clear dominance over others. Therefore, the feature of better detection in our application leads to justified conclusions, which will be acceptable to practitioners and quality control engineers.

Summary, Conclusions, and Future Recommendations
In this study, we proposed a new weighting scheme to design an efficient structure for proposed multivariate homogeneous weighted moving average MHWMA chart. We developed the control charting structure of the proposed chart along with some existing Hotelling's 2 , MCUSUM , 1 , MEWMA , and MEWMA charts. We investigated the performance of all charts under different numbers of quality characteristics and varying sample sizes and perceived that the performance of charts improves as the number of variables increases in the monitoring of the multivariate process. Overall, the proposed chart exhibits the best performance and Hotelling's 2 chart shows the worst performance. The performance of the proposed chart is seriously affected by the usage of a small number of Phase I samples at the estimation stage. Moreover, the in-control performance of the proposed chart is inversely related to the increment in the number of variables and smoothing parameter.
The current study can be extended by using different sampling strategies including ranked set sampling, two phases sampling, and so on in multivariate control charting. Moreover, the current From Figure 9, it can be observed that the Hotelling's T 2 , MEWMA L as , and MHWMA P charts are detecting an upward and downward shift in the process. The Hotelling's T 2 chart indicates the first out-of-control signal at the last sample, while the MEWMA L as chart triggers the first out-of-control signal at the 24th sample. The proposed MHWMA P chart detects the first out-of-control signal at the 24th sample (cf. Figure 9). Moreover, the Hotelling's T 2 , MEWMA L as , and proposed MHWMA P charts detect 1, 3, and 7 out-of-control signals (cf. Figure 9), respectively, which means that the proposed MHWMA P chart has better detection ability than the other charts. These results display the supremacy of the proposed chart against the competing charts, and these results are also consistent with the findings of the study.
The property of better detection ability of MHWMAp chart claimed in Section 5 has a strong connection with the findings of Section 4, where the superiority of the proposal is established based on RL properties. These properties are based on the repeated iterations from a theoretical model by fixing ARL0 at a specific level for all of the charts under investigation, in order to maintain uniformity in performance evaluation. The RL analysis is also carried out at varying shifts, and the resulting outcomes are presented in tabular and graphical forms (cf. Tables 1-4 and Figures 1-5), where the study proposal exhibits a clear dominance over others. Therefore, the feature of better detection in our application leads to justified conclusions, which will be acceptable to practitioners and quality control engineers.

Summary, Conclusions, and Future Recommendations
In this study, we proposed a new weighting scheme to design an efficient structure for proposed multivariate homogeneous weighted moving average MHWMA P chart. We developed the control charting structure of the proposed MHWMA P chart along with some existing Hotelling's T 2 , MCUSUM C , MC 1 , MEWMA L as , and MEWMA L tv charts. We investigated the performance of all charts under different numbers of quality characteristics and varying sample sizes and perceived that the performance of charts improves as the number of variables increases in the monitoring of the multivariate process. Overall, the proposed chart exhibits the best performance and Hotelling's T 2 chart shows the worst performance. The performance of the proposed chart is seriously affected by the usage of a small number of Phase I samples at the estimation stage. Moreover, the in-control performance of the proposed chart is inversely related to the increment in the number of variables and smoothing parameter.
The current study can be extended by using different sampling strategies including ranked set sampling, two phases sampling, and so on in multivariate control charting. Moreover, the current study may also be extended for different runs rules schemes of quality control. Last, but not least, the current study can be of potential use in manufacturing and services processes. Acknowledgments: The authors are grateful to the referees for their valuable comments that helped to improve the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature X
Vector of multivariate variables µ Mean vector Σ Variance-covariance matrix p Number of variables N p (µ, Σ) Multivariate normal distribution As the i th (i = 1, 2, · · · , n) observation of the j th (j = 1, 2, · · · , p) quality Charectersitic on the t th (t = 1, 2, · · · , m) sample X t t th sample mean vector S t t th sample variance-covariance matrix µ 0 Known mean vector Σ 0 Known variance-covariance matrix X Estimator of µ 0 S Estimator of Σ 0 χ 2 t Chi-qquare probability distribution with t degrees of freedom h 0 Control limit of Hotelling's T 2 charts when µ 0 and Σ 0 are used Plotting statistics of Hotelling's T 2 chart when X and S are estimated h 1 Control limit Hotelling's T 2 chart when X and S are used F v 1 ,v 2 Distribution with degree of freedom v 1 = p and v 2 = mn − m − p + Inverse variance covariance matric for the MHWMA p chart h 6 Control limit of the MHWMA p chart µ 0 Estimator of µ 0 Σ 0 Estimator of Σ 0 N p µ 0 , Σ 0 nm Multivariate normal distribution of X δ shift µ 1 Out-of-control mean vector δ max Maximum value of δ δ min Minimum value of δ Appendix A Table A1. Multivariate data structure.
Let X 1 , X 2 , . . . , X j be the j data vectors of size n, which are independently sampled from a multivariate normal distribution at time t with mean vector µ 0 and variance-covariance matrix Σ 0 .

1.
For an in-control process E X j = E X t = µ 0 E X j = E 1 n n i=1 X i j = E 1 n X 1j + X 2j + . . . + X n j E X j = 1 n E X 1j + E X 2 j + . . . + E X n j = 1 n [µ 0 + µ 0 + . . . + µ 0 ] = µ 0 (A1) Now, By substituting µ 0 from (A1), E X t becomes For an in-control process V X j = Σ 0 n and V X t = Σ 0 nt with V X 0 = 0 As X j−1 and X j are vectors consisting of independent random sample values, therefore, Covar X j−1 , X j = 0, so Now, As X t−1 and X t are independent sample means vector, therefore, Covar X t−1 , X t = 0, so By substituting Σ 0 n from (A3), V X t becomes As the variance of a constant is zero, hence V X 0 = V(µ 0 ) = 0