A Bayesian Control Chart for Monitoring Process Variance

Featured Application: The results in this article are applicable to the signal detection in a service process and the monitoring of an intelligent automated process. Abstract: Automation in the service industry is emerging as a new wave of industrial revolution. Standardization and consistency of service quality is an important part of the automation process. The quality control methods widely used in the manufacturing industry can provide service quality measurement and service process monitoring. In particular, the control chart as an online monitoring technique can be used to quickly detect whether a service process is out of control. However, the control of the service process is more difﬁcult than that of the manufacturing process because the variability of the service process comes from widespread and complex factors. First of all, the distribution of the service process is usually non-normal or unknown. Moreover, the skewness of the process distribution can be time-varying, even if the process is in control. In this study, a Bayesian procedure is applied to construct a Phase II exponential weighted moving average (EWMA) control chart for monitoring the variance of a distribution-free process. We explore the sampling properties of the new monitoring statistic, which is suitable for monitoring the time-varying process distribution. The average run lengths (ARLs) of the proposed Bayesian EWMA variance chart are calculated, and they show that the chart performs well. The simulation studies for a normal process, exponential process, and the mixed process of normal and exponential distribution prove that our chart can quickly detect any shift of a process variance. Finally, a numerical example of bank service time is used to illustrate the application of the proposed Bayesian EWMA variance chart and conﬁrm the performance of the process control.


Introduction
Due to the rapid development of artificial intelligence, service automation is emerging as a new wave of industrial revolution. Standardization and consistency of service quality are important parts of the automation process. The quality control methods widely used in the manufacturing industry can provide service quality measurement and service process monitoring. In particular, online process monitoring of the control chart can be used to quickly detect signals that indicates when the service process is out of control. However, the control of a service process is more difficult than that of a manufacturing process because the variability of a service process comes from widespread and complex factors [1][2][3]. First of all, the distribution of the service process is usually non-normal or unknown. Moreover, the skewness of process distribution is time-varying, even if the process is in control.
Because Shewhart variables control charts depend on the normality assumption, several studies have searched for effective methods to measure the performance of control charts with non-normal or unknown distribution data. Among them, most studies deal the new monitoring statistic and find that it is suitable for monitoring the time-varying process. The average run lengths (ARLs) of the proposed Bayesian control chart are calculated and show that the chart performs well. The simulation studies for a normal process, exponential process, and the mixed process of normal and exponential distribution prove that our chart can quickly detect any shift of process variance. Finally, a numerical example of service time is used to illustrate the application of the proposed Bayesian EWMA variance chart and confirm the performance of process control.
The article is organized as follows. In Section 2, we describe the preliminary settings of our control scheme and the distribution of a new statistic derived from the Bayesian approach. Then, the construction of a newly proposed Bayesian variance chart and its performance are discussed in Section 3. In Section 4, we propose a Bayesian exponentially weighted moving average (EWMA) variance chart. In Section 5, a simulation study for a normal process, exponential process, and the mixed process of normal and exponential distributions is performed to test the detection performance of the proposed Bayesian EWMA variance chart. In Section 6, a numerical example of a service system in a bank branch is used to illustrate how to construct the proposed variance chart. Finally, we summarize the findings and provide recommendations in Section 7.

Preliminary Settings
Assuming that a service process is of concern and a critical quality characteristic to be monitored, X, has a mean µ and variance σ 2 . In order to monitor the process variance, a random sample of size 2n, X 1 , X 2 , . . . , X 2n is taken from the service process. For convenience in the calculation of Equation (1), the sample size is even; if not, one observation is deleted. Thus, we define that: It is easy to show that E(Y j ) = σ 2 for j = 1, 2, . . . , n. When M is the total number of instances in which Y j > σ 2 , then M follows from a binomial distribution with parameters n and p, which is defined as: The value of p depends on the population of X and is always between 0.2 and 0.5 [23]. If the population of X is normal distribution, then p = 0.3173, which is independent of variance parameter σ 2 ; if the population of X is exponential distribution, then p = e − √ 2 = 0.2431, which is also independent of rate parameter and hence σ 2 . When p is fixed and known, the binomial distribution can be used to determine the control limits in a process control scheme. However, p is usually unknown in practice because the process distribution may be non-normal and unknown, especially in service processes. Moreover, in some special situations, p may not be fixed even under in-control. The reason for this is that the value of p is not only associated with the shift of variance but also with the skewness of the distribution. Randomized skewness of the distribution will disturb or mislead the judgment about whether the process variance shifts. Thus, we will consider the random behavior of p from the Bayesian perspective. The Bayesian approach views parameter p as a random variable that describes the previous information about parameter p stated in a prior distribution from which we can obtain the posterior distribution or posterior predictive distribution.
Because the beta distribution is the conjugate prior probability distribution for the binomial distributions, it is a suitable model for the random behavior of parameter p. Assuming that the prior distribution of p is a beta distribution: 4 of 19 where the shape parameters α and β are chosen to reflect any existing belief or information.
The statistic M defined above follows a binomial distribution with parameters n and p: Therefore, the marginal distribution of M is the beta-binomial distribution with parameters n, α, and β: Given the prior distribution of p, the probability function of M is then written as: and the mean as well as variance are: and respectively. Notice that the variability of the beta-binomial distribution is larger than that of the binomial distribution with p = α/(α + β), and the difference will get smaller as (α + β) increases or n decreases. Figure 1 depicts several probability functions of the beta-binomial distributions for comparison. Panel A of Figure 1 shows that, compared with the binomial distribution, the beta-binomial distributions are more dispersive. Panel B of Figure 1 shows that the skewness of the beta-binomial distributions increases as α/(α + β) decreases from 0.5.
Because the beta distribution is the conjugate prior probability distribution for the binomial distributions, it is a suitable model for the random behavior of parameter p. Assuming that the prior distribution of p is a beta distribution: where the shape parameters α and β are chosen to reflect any existing belief or information. The statistic M defined above follows a binomial distribution with parameters n and p: Therefore, the marginal distribution of M is the beta-binomial distribution with parameters n, α, and β: Given the prior distribution of p, the probability function of M is then written as: and the mean as well as variance are: and respectively. Notice that the variability of the beta-binomial distribution is larger than that of the binomial distribution with p = α/(α + β), and the difference will get smaller as (α + β) increases or n decreases. Figure 1 depicts several probability functions of the beta-binomial distributions for comparison. Panel A of Figure 1 shows that, compared with the binomial distribution, the beta-binomial distributions are more dispersive. Panel B of Figure 1 shows that the skewness of the beta-binomial distributions increases as α/(α + β) decreases from 0.5.  On the other hand, the prior distribution of p will be updated to a posterior distribution: because statistic M has been observed. Now, we then propose to collect a new sample and obtain a new statistic, M new , which also follows a binomial distribution: Then, the posterior predictive distribution of M new is a beta-binomial distribution:

The Construction of the Shewhart-Type Variance Chart with the Bayesian Approach
Assuming that the prior distribution of p is the beta distribution with parameters α and β, the expected value of p is expressed as E(p) = α/(α + β). Once process variance σ 2 shifts upwards/downwards, the expected process proportion E(p) will change upwards/downwards. Thus, monitoring process variance shifts is equivalent to monitoring the changes in the expected process proportion.

The Parameters of the Prior Distribution of p and Process Variance Are Known
In this subsection, we first develop the variance chart under known process variance σ 0 2 and parameters of the prior distribution of p, α 0 , and β 0 . For the in-control process, the monitoring statistic, M t , defined as the number of Y j 's > σ 0 2 at time t, follows from the beta-binomial distribution with parameters n, α 0 , and β 0 , written as: The mean and standard deviation of M t are and respectively. Thus, we propose a Shewhart-type Bayesian variance chart with the following limits: and plot M t , t = 1, 2, . . . . If either M t ≥ UCL M or M t ≤ LCL M , the process is deemed to show some out-of-variance-control signals.

The In-Control ARL of the Shewhart-Type Bayesian Variance Chart
The average run length (ARL) is used to measure the performance of a control chart. The in-control ARL, denoted by ARL 0 , of the Bayesian variance chart is associated with n, α 0 , and β 0 . Under the in-control process, the chance that the proposed chart alarms a false signal is calculated by: (17) So, ARL 0 = 1/ALARM 0 . Table 1 presents the ARL 0 s of the proposed Bayesian variance chart for various combinations of n, α 0 , and β 0 . (α 0 , β 0 ) E 0 (p) n = 2 n = 3 n = 5 n = 10 n = 15 n = 20 n = 25  There are many infinite ARL 0 s in Table 1, indicating that the proposed Bayesian variance chart does not release false alarm under an in-control process. However, there also appear many ARL 0 s smaller than 370.4, indicating that the proposed chart frequently releases false alarm under the in-control process. The inconsistent results may be attributed to the discontinuity of the monitoring statistic.

The Out-of-Control ARL of the Shewhart-Type Bayesian Variance Chart
When the process is out-of-variance-control and parameters α 0 and β 0 are changed to α 1 and β 1 , the chance that the proposed chart alarms a true signal is calculated by: (18) So, ARL 1 = 1/ALARM 1 . Table 2 presents the ARL 1 s of the proposed Bayesian variance chart for various combinations of n, α 1 , and β 1 under the in-control parameters (α 0 , β 0 ) = (5, 10), (3,9), (46,100), and (32, 100), respectively. There are many infinite ARL 1 s in Table 2, which indicates that the proposed Bayesian variance chart fails to correctly detect any shift of process variance under the out-of-control process. A larger ARL 1 shows that the proposed chart has poor detecting ability for a shift of process variance.

The Parameters of the Prior Distribution of p and Process Variance Are Unknown
When the in-control process variance, σ 0 2 , is unknown, the following preliminary sample data can be used to estimate it. Assuming that are obtained from T sampling periods, each with 2n observations, then On the other hand, when there is no information about the parameters of a prior distribution of p, a non-informative prior distribution, Beta(1, 1), is suggested initially. Then, the preliminary sample data can be used to obtain the posterior distribution of p: where m t represents the number of Y j 's >σ 2 at the tth sampling period. Because a new statistic, M T+t , follows from a binomial distribution with parameters n and p, the posterior predictive distribution of M T+t is the beta-binomial distribution with parameters n,α, and β, written as: Thus, the in-control expected value of p is estimated with: The upper and lower control limits of M T+t are obtained, respectively, as follows: The M t for the T preliminary samples are plotted on the resulting Bayesian variance chart. If no points fall outside the control limits, then the process is regarded as the incontrol state. Thus, the M T+t for the subsequent monitoring samples are plotted in order. If either M T+t ≥ UCL M or M T+t ≤ LCL M , the process is deemed to show some out-of-variancecontrol signals. In the following discussions, we assume that the number of preliminary samples, T, is large enough so thatσ 2 andα α+β are almost the same as the true values of σ 0 2 and E 0 (p), respectively.

The Construction of the EWMA Variance Chart with the Bayesian Approach
Because observed statistic M t is discrete and the beta-binomial distribution is asymmetric, the in-control ARLs of the proposed Shewhart-type Bayesian variance chart cannot always satisfy the required in-control ARLs, such as 370.4. The same reason results in the anomaly that the out-of-control ARLs of the proposed Shewhart-type Bayesian variance chart do not inversely change with n as they normally should.
The EWMA chart has been used to monitor small shifts of process variance quickly and effectively. Moreover, the characteristics of the moving average of observed statistics may ease the disadvantages of discreteness and asymmetry.
We define a new EWMA statistic as: where λ is the smoothing parameter. Because E(M t ) = nα α+β and M 0 is set to be the mean of M t , the mean of EWMA t is specified as: On the other hand, because Var(M t ) = nαβ(α+β+n) , the variance of EWMA t is specified as: Thus, central limit theorem would ensure that the EWMA t will asymptotically follow from a continuous symmetric distribution with [15] and

Construction of the EWMA Variance Control Chart
Assuming that the in-control process variance is known to be σ 0 2 and the prior distribution of p is also known to be Beta(α 0 , β 0 ), we can further construct a Bayesian EWMA variance chart with the following control limits: where (k 1 , k 2 ) are coefficients of the control limits. The adjustment allowing unequal coefficients of the control limits conforms to the asymmetry of beta-binomial distribution.
Next, we will discuss how to choose suitable coefficients (k 1 , k 2 ) to satisfy a required in-control average run length. The ARLs are again used to measure the performance of the proposed Bayesian EWMA variance chart. We evaluate the ARLs of the proposed Bayesian EWMA variance chart by a Monte Carlo simulation approach. In order to ensure that the in-control ARLs are at least 370.4, the coefficients of the control limits (k 1 , k 2 ) are chosen with a procedure as follows.
Determine the least k 1 to make sure the ARL(k 1 ) is larger than 740.8; 6.
Determine the least k 2 to make sure the ARL(k 2 ) is larger than 370.4. Table 3 shows the coefficients of the control limits (k 1 , k 2 ) of the proposed Bayesian EWMA variance chart with λ = 0.05, where the ARL 0 s are slightly beyond 370.4, for various different n, α 0 , and β 0 . The results in Table 3 show that k 1 and k 2 are close, while E 0 (p) = 0.5; however, k 1 changes inversely with E 0 (p) and k 2 changes positively with E 0 (p). The in-control area gets bigger due to the increasing skewness of the beta-binomial distribution.

Evaluation of the EWMA Variance Control Chart
To realize the performance of the proposed Bayesian EWMA variance chart, the outof-control ARL will be calculated. First of all, the coefficients of the control limits, k 1 and k 2 , are determined from the procedure described in the previous subsection. When the process is out-of-variance-control and the parameters of the prior distribution of p are changed to α 1 and β 1 , the ARL 1 s of the proposed Bayesian EWMA variance chart are evaluated by the Monte Carlo simulation approach. The procedure of evaluation is presented as follows: 1.
Repeat step 3 1000 times and obtain the average run length, ARL 1 . Table 4 shows the ARL 1 s of the proposed Bayesian EWMA variance chart with the combination of k 1 and k 2 for various different n, α 1 , and β 1 under the in-control parameters (α 0 , β 0 ) = (5, 10), (3,9), (46,100), (32, 100), and smoothing parameter λ = 0.05.   (9,9) 0.5000 The results in Table 4 show that the ARL 1 decreases quickly as E 1 (p) shifts away from E 0 (p) in each set of sample size, documenting that the proposed Bayesian EWMA variance chart possesses excellent detecting power for any shift of the process variance.

Performance Comparison of the EWMA Variance Control Chart
In general, the EWMA type control chart performs better than the Shewhart-type chart for small shifts. The single sampling EWMA variance (SS EWMA-V) chart has proved better performance than the NLE, CEW, NP-M, and IRC charts from previous studies [23,34,35]. To demonstrate that the proposed Bayesian EWMA variance chart outperforms the existing charts, the SS EWMA-V chart is considered for comparison. Table 5 presents the out-of-control average run lengths of the proposed Bayesian EWMA variance chart and the SS EWMA-V chart with a sample size of 10 (n = 5) and smoothing parameter 0.05 (λ = 0.05). The results of Table 5 show that, while E 0 (p) and E 1 (p) are both equal to p 0 and p 1 , respectively, the ARL 1 s of the Bayesian EWMA variance chart are all smaller than those of the SS EWMA-V chart, showing that our proposed chart could detect a shift of variance more quickly than the SS EWMA-V chart.   Notes: sample size is 10 (n = 5); smoothing parameter is 0.05 (λ = 0.05); a Denotes an in-control process.
Hence, we conclude that the out-of-control detection performance of the proposed Bayesian EWMA variance chart is always better than those of the SS EWMA-V chart, no matter whether the shift from p 0 /E 0 (p) to p 1 /E 1 (p) is small, medium, or large.

Process Simulations
The out-of-control variance detection performance of our proposed Bayesian EWMA variance chart is further evaluated under the normal process, exponential process, and the mixed process of normal and exponential distribution. The process with mixed distribution demonstrates a situation where the process sometimes is normally distributed and sometimes is exponentially distributed, both with the same variance but with different skewness.
For any normal distribution with standard deviation σ, while σ shifts to dσ, the probability that statistic Y j is larger than σ 2 can be computed as: where Φ(.) represents the cumulative probability function of a standard normal variable and d represents a shift factor. Panel A of Table 6 presents the relationship between d and p in the normal process to look at how p varies with the shifts of process variance. For any exponential distribution with standard deviation σ, while σ shifts to dσ, the probability that statistic Y j is larger than σ 2 can be computed as: where d represents a shift factor. Panel B of Table 6 presents the relationship between d and p in the exponential process to look at how p varies with the shifts of process variance. Both processes in Table 6 indicate that p is positively associated with d.
Suppose that the mixed process is an equally-weighted mixture of a normal distribution and exponential distribution with same variance σ 2 . Thus, for the mixed distribution with standard deviation σ, while σ shifts to dσ, the probability that statistic Yj is larger than σ 2 is expressed as: Because probability p N is the same for any normal distribution, the in-control normal process is set a standard normal distribution in the following simulations, without loss of generality. Similarly, probability p E is the same for any exponential distribution, so the in-control exponential process is set an exponential distribution with standard deviation 1, without loss of generality.
The following simulation procedure will be applied to three processes to evaluate the performance of the proposed Bayesian EWMA variance chart:
α 0 and β 0 are chosen such that E 0 (p) is as close to 0.3173 as possible in the normal process (0.2431 in the exponential process or 0.2802 in the mixed process); 3. k 1 and k 2 are determined by the procedure described in Section 4.1; 4.
Repeat step 6 1000 times and obtain the average run length, ARL 1 .
While simulating random samples of size 2n from the equally-weighted mixed process, we assume that observations within one sample come from the same distribution; that is, the skewness of distribution varies only between samples. Table 7 shows the simulation results of the ARL 1 s of the proposed Bayesian EWMA variance chart under a normal process, exponential process, and the equally-weighted mixed process of normal and exponential distributions. a denotes the in-control process; b denotes the coefficients of the control limits (k 1 , k 2 ); c denotes the preset ARL 0 .
Panel A of Table 7 indicates that the simulated ARL 0 s are around the preset ARL 0 s in a normal process. Panel B of Table 7 indicates that the simulated ARL 0 s, except for n = 2, are around the preset ARL 0 s in an exponential process. However, Panel C of Table 7 indicates that the simulated ARL 0 s are significantly different from the ARL 0 s in the mixed process. The less robust ARL 0 s in the mixed process result from the possibility that the prior distribution of p may be bimodal. Nevertheless, the ARL 1 s in Table 7 reversely (positively) change with d when d > 1 (<1) and quickly reduce to 1 as d departs from 1, demonstrating that the proposed Bayesian EWMA variance chart possesses excellent detecting power for any upward or downward shift of the process variance.

An Example for Demonstration
A banking example from Yang et al. [21] is applied to illustrate the proposed Bayesian EWMA variance chart. The banking industry regards service time as an important quality characteristic. Under the demand for standardization and uniformity of service time, a control chart will be constructed to monitor the service time and detect whether the variance is offset.
Service time usually comes from a complex service process, whose distribution may not be normal. Sometimes, the distribution of the service time is time-varying even under the in-control condition. Because we know little about the nature of the service process, the proposed Bayesian EWMA variance is applicable to the banking case.
In order to measure the efficiency in the service system of a bank branch, the incontrol sampling service time was measured from 10 counters every day for 15 days. So, 15 samples of size 2n = 10 were collected. A histogram shows that these samples are right-skewed. Table 8 (21) as well asÊ 0 (p) = 23/77 = 0.2987 by Equation (22). Past experience has confirmed that 30.0969 is very close to the true value of in-control process variance, σ 0 2 , and 0.2987 is also very close to the true value of the in-control mean of p, E 0 (p). When the smoothing parameter λ is 0.05, referring to the procedure described in Section 4.1, we obtain k 1 = 3.07 and k 2 = 2.86. Thus, the control limits are calculated as UCL EWMA = 2.0094 and LCL EWMA = 1.0129 through Equations (31) and (32). From Equation (30), we preset EWMA 0 = 1.4935, then EWMA t is computed in order by Equation (25) and presented in Table 8. Figure 2 shows the proposed Bayesian EWMA variance chart. The observed statistics EWMA t s are plotted in order. Because all the EWMA t s fall into the region between control limits, the chart reveals that the process is in control.  Figure 2 shows the proposed Bayesian EWMA variance chart. The observed statistics EWMA t s are plotted in order. Because all the EWMA t s fall into the region between control limits, the chart reveals that the process is in control. In order to illustrate the detecting power for the variance change, ten new samples of size 10 were collected from a new automatic service system of the bank branch. Installing an automatic service system makes more consistent service time. The process variance will be reduced significantly. Table 9 lists the service time of the ten new out-of-variance-control samples, M t , and EWMA t . Note that the out-of-control process variance is estimated with = 2.9437, showing a significant downward shift in the process variance. Hence, the observed M t s are all 0 and the observed EWMA t s decrease with t. In order to illustrate the detecting power for the variance change, ten new samples of size 10 were collected from a new automatic service system of the bank branch. Installing an automatic service system makes more consistent service time. The process variance will be reduced significantly. Table 9 lists the service time of the ten new out-of-variance-control samples, M t , and EWMA t . Note that the out-of-control process variance is estimated withσ 2 = 2.9437, showing a significant downward shift in the process variance. Hence, the observed M t s are all 0 and the observed EWMA t s decrease with t.
We then plotted the ten new observed statistics EWMA t s in order on the proposed control chart, see Figure 3. The proposed Bayesian EWMA variance chart detected out-ofcontrol signals from the first sample onward (samples 1-10). That is, our chart can quickly detect an out-of-control signal when the variability of the new service time is significantly reduced because of the improved new automatic service system.
We then plotted the ten new observed statistics EWMA t s in order on the proposed control chart, see Figure 3. The proposed Bayesian EWMA variance chart detected out-of-control signals from the first sample onward (samples 1-10). That is, our chart can quickly detect an out-of-control signal when the variability of the new service time is significantly reduced because of the improved new automatic service system.

Conclusions
In this paper, we propose a new control chart, the Bayesian EWMA variance chart, to monitor changes in the process variance when the distribution of a critical quality characteristic is non-normal, unknown, or skewness-time-varying. The average run lengths (ARLs) of the proposed chart are calculated and show that the chart performs well. When compared with the well-performing SS EWMA-V chart, the out-of-control detection performance of the proposed Bayesian EWMA variance chart is always better than that of the SS EWMA-V chart, no matter whether the shift is small, medium, or large. The simulation studies for a normal process, exponential process, and the mixed process of normal and exponential distribution prove that our chart can quickly detect a small, medium, or large shift of process variance. Moreover, a numerical example of service time is used to illustrate the application of the proposed Bayesian EWMA variance chart and confirm the performance of process control.
However, when monitoring the mixed process, the in-control ARLs seem to be biased due to the occurrence of the bimodal random behavior of p, which will limit the use of the proposed chart, especially in situations that have a large dispersion of skewness. A future study could consider multimodal prior distribution to catch the random behavior of p resulting from the skewness-time-varying process. Moreover, this study does not test or compare the effects of different weights under a mixed distribution. This may be a good issue for further studies.

Conclusions
In this paper, we propose a new control chart, the Bayesian EWMA variance chart, to monitor changes in the process variance when the distribution of a critical quality characteristic is non-normal, unknown, or skewness-time-varying. The average run lengths (ARLs) of the proposed chart are calculated and show that the chart performs well. When compared with the well-performing SS EWMA-V chart, the out-of-control detection performance of the proposed Bayesian EWMA variance chart is always better than that of the SS EWMA-V chart, no matter whether the shift is small, medium, or large. The simulation studies for a normal process, exponential process, and the mixed process of normal and exponential distribution prove that our chart can quickly detect a small, medium, or large shift of process variance. Moreover, a numerical example of service time is used to illustrate the application of the proposed Bayesian EWMA variance chart and confirm the performance of process control.
However, when monitoring the mixed process, the in-control ARLs seem to be biased due to the occurrence of the bimodal random behavior of p, which will limit the use of the proposed chart, especially in situations that have a large dispersion of skewness. A future study could consider multimodal prior distribution to catch the random behavior of p resulting from the skewness-time-varying process. Moreover, this study does not test or compare the effects of different weights under a mixed distribution. This may be a good issue for further studies.

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