Multivariate control chart based on PCA mix for variable and attribute quality characteristics

ABSTRACT Two types of control charts exist based on different quality characteristics: variable and attribute. These characteristics are commonly monitored using separate procedures. Only a few studies focused on the utilization of control charts to monitor a process with mixed characteristics. This study develops a new concept of the control chart based on a Principal Component Analysis (PCA) Mix, that is a PCA method that can jointly handle continuous and categorical data. The Kernel Density Estimation (KDE) method is used to estimate the control limit. Through simulation studies, the performance of the proposed chart is evaluated using the Average Run Length (ARL). control limits obtained from KDE produce a stable ARL0 at ~ 370 for For the shifted process, the proposed chart demonstrates excellent performance for an appropriate number of principal components used. Applications of the simulated process and real cases show that the proposed chart is sensitive to monitoring the shifted process.


Introduction
One of the most powerful tools in Statistical Process Control (SPC) is the control chart, which has been widely used in industries and services. Based on the type of monitored quality characteristics, there are two kinds of control charts: variable (interval or ratio) and attribute (category) (Montgomery, 2009). The variable control chart is appropriate to use if the quality characteristics of the product are measured on a numerical scale such as weight and height. On the other hand, if the quality characteristics are not measured on a numerical scale (in categorical form) such as color and softness, an attribute control chart can be employed (Sorooshian, 2013). The control chart is not only developed for the univariate case but also for the multivariate case where each characteristic cannot be monitored separately. Therefore, multiple quality characteristics should be monitored together using a multivariate control chart.
The recent development of the multivariate control chart based on the Shewhart approach includes robust Hotelling's T 2 control chart (Abu-Shawiesh, Kibria, & George, 2014;Alfaro & Ortega, 2009;Ali, Syed Yahaya, & CONTACT Muhammad Mashuri m_mashuri@statistika.its.ac.id This article has been republished with minor changes. These changes do not impact the academic content of the article. Saracco (2014). Although PCA Mix can be employed for different types of quality characteristics, the distribution of Principal Component Scores (PCs) produced by PCA Mix does not follow a multivariate normal distribution or even unknown distribution. Consequently, a conventional T 2 control chart constructed by PCA Mix generates many false alarms because it assumes a multivariate normal distribution. To overcome this issue, a Kernel Density Estimation (KDE) can be employed to calculate the control limit (Ahsan, Mashuri, Kuswanto, Prastyo, & Khusna, 2018;Chou, Mason, & Young, 2001;Phaladiganon, Kim, Chen, Baek, & Park, 2011;Phaladiganon et al., 2013).
The key contribution of this paper lies in forming a new procedure for monitoring the mixed characteristics in a single chart. To date, the product in a manufacturing process is only monitored using one of the characteristic types measured (only considering variable or attribute characteristics). If both characteristics are monitored, this is commonly done using separate procedures. This actually becomes ineffective because it needs to see the results from the two different procedures before finally making a decision on the quality of the product. For example, when the length and weight of a product have been stated in an in-control condition, it is also necessary to analyse the defects and level of softness from the product using different procedures. If the monitoring results of the two procedures give different results, the decision-maker will have difficulty in determining the quality of the product. Maybe the product will be declared as in-control if the two procedures produce the same results, or the decision-maker will prefer the results from one procedure when the two procedures yield different outcomes.
Based on the aforementioned reasons, this paper focuses on developing a new approach to monitor a process with a mixture of variable and attribute characteristics. This new approach is expected to simplify the monitoring process without reducing the level of precision and accuracy. A multivariate T 2 control chart based on PCA Mix is proposed, while KDE is employed to estimate its control limit. The performance of the proposed chart with a control limit calculated using KDE is compared with a proposed chart that employs conventional F distribution to calculate the control limit using the Average Run Length (ARL). Moreover, the application of the proposed chart to monitor simulated and real case data is presented in this study. This paper is organized as follows. Section 2 describes the charting procedure of a T 2 chart based on PCA Mix. Section 3 presents some simulation studies. Next, the application of the proposed chart for simulated data and real cases is displayed in section 4. The managerial implications of the proposed chart are discussed in section 5. Finally, section 6 summarizes the obtained results.

Charting procedures 2.1. PCA mix
Multivariate data analysis can be described as statistical methods to analyse data consisting of two or more quality characteristics. These quality characteristics can be either variable (interval or ratio) or attribute (category). PCA is one of the statistical methods for the dimensional reduction of continuous data, which in SPC are referred to as variable characteristics. MCA is an extension of Correspondence Analysis (CA) and is used to analyse the relationship pattern of several correlated categorical variables, which in SPC are referred to as attribute characteristics. MCA can be viewed as a generalization of the PCA method when the observations are categorical (Jolliffe, 2002). Thus, the different types of quality characteristics can be handled together by using the PCA Mix method as a combination of PCA and MCA.
The PCA Mix procedure used in this paper follows the approach proposed by Chavent et al. (2014). Let n Â p matrix X 1 and n Â q matrix X 2 consist of variable and attribute characteristics, respectively, where n is the number of observations, p is the number of variable characteristics, and q is the number of attribute characteristics. An indicator matrix G with dimensions n Â m contains binary coding from each level of attribute characteristics, where m is the number of total levels of attribute characteristics. An n Â ðp þ mÞ matrix Z ¼ ½Z 1 ; Z 2 consists of a real number element, where Z 1 and Z 2 are centred matrices of X 1 and G . Z , is calculated as where N ¼ 1 n I n is the weights of the rows of Z, M ¼ diag 1; :::; 1; n n 1 ; :::; n n m is the weights of the columns of Z, the first p columns of Z are weighted by 1, and the last m columns are weighted by n n s ; for s ¼ 1; 2; . . . ; m: The next step is solving the eigenvalue problem of Z , using the Generalized Singular Value Decomposition (GSVD) in Chavent et al. (2014) as where λ 1 ; λ 2 ; . . . ; λ r are the eigenvalues of Z , ; and r denotes the rank of Z , : Matrix U, which has n Â r dimensions, is an eigenvector of Z , , and V is the ðp þ mÞ Â r matrix of the eigenvectors of Z , : Thus, the principal component of PCA mix can be computed as with the size of n Â r:

Kernel density estimation
The kernel is a weighting function commonly used in nonparametric techniques. The kernel estimator was first introduced by Murray Rosenblatt (1956) and Parzen (1962) and was called the Rosenblatt-Parzen kernel density estimator. Chou et al. (2001) introduced KDE to estimate the distribution of T 2 statistics. Given n values of T 2 i statistics obtained from an in-control condition, the empirical density of T 2 can be expressed using the following kernel function: where K andĥ are the kernel function and estimated smoothing parameter, respectively.
Furthermore, the control limit of T 2 based on KDE can be determined by the percentile of kernel distribution (Phaladiganon et al., 2011). The cumulative distribution function for KDE can be written as Thus, the control limit of T 2 based on KDE is the 100ð1 À αÞ th percentile of the empirical density off h ðtÞ, as follows: The integral in equation (5) is not simple to solve owing to the complication of the kernel function. Thus, the trapezoidal rule (Burden & Faires, 2011), one of the numerical integration methods to approximate the definite value of an integral equation, is used to solve the kernel control limit.

T 2 control chart based on PCA
The conventional T 2 control chart is constructed based on PCA uses the first k PCs. The statistic of the T 2 chart based on PCA can be computed by using the following formula: where the first k PCs are y v ; v ¼ 1; :::; k, and λ v is the eigenvalue corresponding to the v-th PC. Under the assumption that the data follow a multivariate normal distribution, the control limit of a PCA-based T 2 control chart can be obtained as follows: CL ¼ kðn þ 1Þðn À 1Þ n 2 À nk F ðα;k;nÀkÞ ; where n is the number of observations, k is the number of PCs retained, and α is the false alarm rate.

T 2 control chart based on PCA mix
In this research, a T 2 control chart is constructed based on PCA Mix with KDE using the PCs obtained from PCA Mix, which is denoted as Y mix in equation . The proposed multivariate control chart is constructed by T 2 statistics for the first k PCs of Y mix : The T 2 statistic for PC Mix, which is denoted byT 2 ; is calculated as follows: where v ¼ 1; 2; :::; k,λ v is an eigenvalue that corresponds to the vth PCs from the incontrol process; and μ is the mean vector of the in-control process. Furthermore, the control limit of the proposed chart is calculated using the KDE method because the distribution of theT 2 statistic is unknown. The Gaussian KDE of theT 2 statistic is estimated as follows:f where the Gaussian kernel is defined as Moreover, the distribution function off h ðtÞ that is denoted byF h ðtÞ is calculated bŷ F h ðtÞ is calculated using the trapezoid rule method as follows: where π min and π max are the minimum and maximum value off h ðT 2 Þ: Hence, the control limit ofT 2 based on KDE is the 100ð1 À αÞ th percentile and is calculated as follows: Furthermore, the established control limit is used to monitor new observations.

Simulation studies
In this section, several simulation studies are conducted to investigate the distribution of Y mix ; to determine the control limit of the proposed chart, and to evaluate its performance. The variable characteristics are assumed to follow the multivariate normal distribution X 1 ,N p ðμ; AEÞ: Meanwhile, the attribute characteristics follow the multinomial distribution X 2 ,Mðn; θ 1 ; θ 2 ; :::; θ m Þ: The number of variable characteristics p used in this study are 5, 10, and 30; the number of principal components is k; and the number of observations n = 1000. One attribute characteristic follows the multinomial distribution Mðθ 1 ; θ 2 ; θ 3 Þ used in this study. There are three scenarios for values of θ 1 ; θ 2 ; and θ 3 : (1) θ 1 ; θ 2 ¼ 0:3 and θ 3 ¼ 0:4; (2) θ 1 ; θ 2 ¼ 0:1 and θ 3 ¼ 0:8; (3) θ 1 ; θ 2 ¼ 0:05 and θ 3 ¼ 0:9: The combinations θ 1 ; θ 2 ¼ 0:3; and θ 3 ¼ 0:4 parameters represent the slightly balanced class in the attribute characteristics, while the others represent imbalanced as well as extreme imbalanced classes. Furthermore, the KDE is used to calculate the control limit, and the ARL is used to evaluate the performance of the proposed chart. The ARLs are calculated through a simulation process for 1000 replications.

Distribution of Y mix
The simulation study in this section investigates whether the distribution of PC Mix follows the multivariate normal. The variable characteristics are generated by a multivariate normal distribution with mean vector μ ¼ 0 and covariance matrix AE¼I . The attribute characteristics are generated by a multinomial distribution with three types of parameter, as mentioned previously. Figure 1 shows a Q-Q plot of PC Mix for p = 5 and k = 4 with μ ¼ 0 and AE¼I: Figure 2 illustrates the Q-Q plot of PC Mix for p = 10 and k = 5. Meanwhile, Figure 3 displays the Q-Q plot of for p = 30 and k = 20. From these figures, it can be seen that the distribution of PC Mix does not follow the multivariate normal distribution confirmed by the statistic and the p values of a Royston multivariate normal test. However, there is an indication that the distribution of PC Mix will follow the multivariate normal distribution when the number of quality characteristics increases while the parameter for the attribute characteristic is balanced at θ 1 ; θ 2 ¼ 0:3 and θ 3 ¼ 0:4 (see Figure 3). The distribution of PC Mix deviates from a multivariate normal distribution when the parameter for the attribute characteristic is imbalanced, for example, when θ 1 ; θ 2 ¼ 0:1 and θ 3 ¼ 0:8 as well as θ 1 ; θ 2 ¼ 0:05 and θ 3 ¼ 0:9 . These results are also   supported by increasing the p-value when the number of quality characteristics is increasing and the parameters of the attribute characteristics are balanced.
In order to prove these findings, four statistical methods for multivariate normal distribution testing are employed: Mardiah, Henze-Zirkler, Royston, and Generalized Shapiro-Wilk. This simulation is repeated 1,000 times for various combinations of p, k, and values of θ . Furthermore, the probability of failing to reject the null hypothesis for each combination is calculated. Table 1 shows the probability of PCs to follow a multivariate normal distribution. In this simulation, various p and k are combined with θ 1 ; θ 2 ; and θ 3 to calculate the probability. The simulation process reveals that all combinations produce a zero probability of following a multivariate normal distribution. However, there is still a possibility that PCs will follow a normal multivariate distribution when using very large p and k as well as the balanced parameter of the attribute characteristics.
From these investigations, it can be concluded that PC Mix does not follow a multivariate normal distribution, nor does its e T 2 statistic follow the F distribution that the conventional T 2 -based PCA chart does. Hence, it is not appropriate to use the assumption of multivariate normal for this proposed chart.

Control limit
The findings in section 3.1 proved that PC Mix does not follow a multivariate normal distribution and that its T 2 statistic does not follow an F distribution. As a result, KDE is employed in order to estimate an unknown distribution. The control limit ofT 2 needs to be estimated using both the KDE method and conventional F distribution as in equation in order to compare the performance of these control limits by using ARL 0 criteria. Table 2 presents a comparison ofT 2 control limits with its ARL 0 using the KDE method and an F distribution. Since the simulation study is conducted using a significance level α ¼ 0.00273 that corresponds to a theoretical ARL 0 of 370, theT 2 control limits using the KDE method are more reliable than those using an F distribution. For several combinations of p, k, and parameters of attribute characteristics, theT 2 control limits using the KDE method always produce an empirical ARL 0 at about 370. Meanwhile, the ARL 0 of theT 2 control limits using an F distribution is not stable. However, for p ¼ 30; k ¼ 20 and the parameter of attribute variables θ 1 ; θ 2 ¼ 0:1andθ 3 ¼ 0:8, both the ARL 0 and the control limit using the KDET 2 method are coincidentally similar to those using an F distribution.

Performance of proposed chart
This section provides a performance evaluation of the proposedT 2 control chart using the KDE method. The out-of-control ARL is calculated by adding a shift for each variable characteristic μ shift ¼ μ þ δ μ ; where δ μ ¼ 0:1 and the attribute characteristics are shifted by θ shift ¼ ½θ 1 À δ θ ; θ 2 À δ θ ; θ 3 þ 2δ θ ; where δ θ ¼ 0:0025. Three kinds of numbers of PCs (representing low, moderate, and high) are evaluated for each number of variable characteristics. Moreover, the KDE control limits for several combinations of p, k, and the parameter of the attribute characteristic are taken from Table 2. A simulation study is conducted using significance level α ¼ 0.00273 that corresponds to a theoretical ARL 0 of 370. The empirical values of ARL 0 obtained from the simulation studies are presented in Table 4-7. The ARLs of the proposedT 2 control chart using the KDE method with θ 1 ; θ 2 ¼ 0:3 and θ 3 ¼ 0:4 are listed in Table 3. For various combinations of the number of variable characteristics and the number of PCs, the proposed control chart yields a stable ARL 0 at about 370. Meanwhile, ARL 1 decreases as the variable and attribute parameter shifts increase. For the number of variable characteristics p equal to 5 and 10, the ARL 1 value increases as the number of PCs used increases. In the other words, the ability of theT 2 control chart using the KDE method in detecting actual shifts decreases for a large number of PCs. Table 4 presents the ARLs of theT 2 control chart using the KDE method with θ 1 ; θ 2 ¼ 0:1 and θ 3 ¼ 0:8: The larger shift of a process leads to a better performance Table 3. ARLs for θ 1 ; θ 2 ¼ 0:3 and θ 3 ¼ 0:4 with several combinations of p and k.  of theT 2 control chart to detect shifts, which are confirmed by the ARL 1 value. As the PC number increases, the ability of the proposed control chart to detect a specified shift also increases. However, for p = 5 and k = 2, the ability of the proposed chart to detect small shifts is poor. This may be a result of the total variance, which can be explained by a small number of PCs that is not enough to represent the process data. Thus, for an imbalanced class of the attribute characteristic, theT 2 control chart detects shifts more quickly when the process has a larger shift and number of PCs. The ARLs of theT 2 control chart using the KDE method with θ 1 ; θ 2 ¼ 0:05 and θ 3 ¼ 0:9 are listed in Table 5. Based on the simulation results, the value of ARL 1 generally decreases as the number of PCs increases. However, the ability of the proposed chart to detect shifts in the process is poorer when the number of variable characteristics p used is smaller. Hence, for the extreme imbalanced class of attribute characteristics, the proposed chart exhibits better performance for large numbers of variable characteristics and components used.
This study also compared the performance of theT 2 control chart for several numbers of quality characteristics with a certain total variance explained by PCs. For p = 5, 70% of the total variation can be explained using the first four principal components. Meanwhile, for p = 10, this can be explained by the first six components, and p = 30 by the first 15 components. Table 6 presents an ARL comparison of the proposed control chart for several quality characteristic numbers with the total variance explained by PCs of 70% by the first k components. For several combinations of attribute parameters, ARL 1 of theT 2 control chart decreases as the number of variable characteristics increases. This points out that for certain total variances explained by an appropriate number of PCs, a higher number of variable characteristics resulted in a higher ability of theT 2 control chart to detect a specified shift of the process. Figure 5 illustrates an ARL comparison for various θ 1 ; θ 2 ; and θ 3 at different numbers of variable characteristics as summarized in Table 6. Figure 5 shows that the proposed chart has better performance in monitoring the process with extreme imbalanced Table 5. ARLs for θ 1 ; θ 2 ¼ 0:05 and θ 3 ¼ 0:9 with several combinations of p and k. parameters for small shifts in the process. However, it exhibits poorer performance in detecting large shifts in the process. On the other hand, for a process with balanced parameters, the proposed chart has better performance in detecting a large shift.

Application to simulated process
In this section, the application of the proposed control chart is illustrated using the simulation data for several combinations of p, k, and parameters of attribute Table 6. ARLs ofT 2 control chart with k that can explain 70% of total variance for several p.   Table 7. Scenario of setting for simulation data. characteristics as summarized in Table 7. For this study, an in-control process with 70 observations for variable characteristics is assumed to follow a multivariate normal distribution N p ð0; IÞ, and the attribute characteristic follows a multinomial distribution Mðθ 1 ; θ 2 ; θ 3 Þ: The next 30 observations are the shifted process generated by multivariate normal distribution N p ð2; IÞ for variable characteristics, with the same process for attribute characteristics. The control limits used in this section are taken from Table 2. Figure 6a shows the application of the proposed chart for scenarios i, ii, and iii. For all parameters of the attribute characteristics, the chart quickly detects a shift at the 71st observation. When the parameter of the attribute characteristic is extremely imbalanced for θ 1 ; θ 2 ¼ 0:05 and θ 3 ¼ 0:9 (scenario iii), the shifted process is not seen as clearly as with the other settings. Figure 6b shows the application of the proposed chart for scenarios iv, v, and vi. Similar to the previous scenario, the chart quickly detects the shift at the 71st observation for all attribute characteristic parameters. Meanwhile, for an extreme imbalanced parameter of the attribute characteristic (scenario vi), it can be seen that the process is shifted, which is not clearly seen in the previous figure (scenario iii). Better results are found for scenarios vii, viii, and ix, as shown in Figure 6c. For all parameters of attribute characteristics, the proposed chart can clearly identify the shifts in the process. This is indicated by the fact that almost all of the shifted observations fall outside the control limit.

Application to a real case
In this section, the proposed chart is applied to a real case. The dataset used to evaluate the performance of the proposed chart is the Machine Failure dataset. This dataset can b e d o w n l o a d e d a t h t t p s : / / b i g m l . c o m / u s e r / c z u r i a g a / g a l l e r y / d a t a s e t / 587d062d49c4a16936000810. The dataset is collected by monitoring a one-year basis machine, and the information about failures of the process is recorded. This dataset consists of 16 variable characteristics and three attribute characteristics.
The number of variable characteristics used in this study is eight, while the number of attribute characteristics is two. The number of principal components and the level of significance are five and 0.00273, respectively. The first 400 samples are treated as incontrol observations (training dataset) and are used to estimate the control limit using KDE. By contrast, the 401st to 650th observations are treated as testing datasets and are monitored using the estimated control limit obtained from in-control observations. Figure 7 shows the utilization of theT 2 control chart to monitor the process for the testing dataset. The control limit obtained from KDE is 14.872. The original dataset has information about in-control and out-of-control labels for each observation so that the performance of the proposed control chart can be easily evaluated. Table 8 shows a Figure 7. Application ofT 2 control chart to monitor the 401st to 650th observation of machine failure dataset using KDE control limit. performance comparison between the proposed chart and conventional chart in monitoring the real case. Given the null hypothesis that the observations in the testing dataset are in control, the proposed control chart produces zero false alarms and has an accuracy detection of about 99.6%. The proposed control chart successfully detects two of three actual out-of-control observations. Compared to the proposed chart, the application of the conventional chart to monitor the process for the testing dataset is presented in Figure 8. Using the conventional F control limit (25.222), the conventional chart successfully detects the actual incontrol and out-of-control observations with an accuracy of 99.2% and a detection rate of about 33.3% (detecting 1 out of 3 actual out-of-control observations). As a note, the number of actual out-of-control observations in the testing dataset is only three. These findings show that the proposed chart has better performance than the conventional chart when monitoring actual out-of-control observations.
There are two reasons that the proposed chart has better performance. First, the proportion of the attribute characteristics of the testing dataset is balanced. The first attribute characteristic has eight categories with slightly balanced proportion (0.104, 0.223, 0.096, 0.100, 0.096, 0.127, 0.127, and 0.127), while the second attribute characteristic (four categories) has a balanced proportion (0.267, 0.259, 0.239, and 0.235). According to the simulation studies, for a small number of variable characteristics, the proposed chart is at the peak of its performance when it is used to monitor the balanced attribute characteristics (see Table 3 and Figure 6a). These facts reveal that a more accurate monitoring process will be achieved if these attribute characteristics are included in the analysis. Second, the KDE control limit captures the actual out-of-control observations better than the conventional F distribution control limit, which was developed under a multivariate normal assumption. However, the eight variable characteristics in the Machine Failure testing dataset do not follow the multivariate normal distribution confirmed by the small p-value of the Mardiah test < 0:0001 ð Þand the Q-Q plot, as shown in Figure 9. These facts indicate that the conventional control chart only has the best performance when an assumption is fulfilled. By using the KDE method, the empirical density of the e T 2 statistic can be calculated, and the control limit can be estimated by taking its 100ð1 À αÞ th percentile as defined in equation . The estimated KDE control limit is superior in detecting the actual out-of-control observations, whereas the conventional control limit is not.

Managerial implication
In this globalization era, monitoring the quality of products using a control chart is essential for continuous process improvement. The main objective of a control chart is to monitor the process and reduce variations in the process. This monitoring process is done by identifying and eliminating assignable causes. An out-of-control observation issued by a control chart indicates the presence of assignable causes. To solve this problem, a root-cause analysis is conducted to reveal the reasons for the assignable causes. However, the conventional monitoring process using a control chart only considers one type of quality characteristic. For example, a variable control chart will be employed to monitor only numerical measurements such as weight and height, while an attribute control chart monitors categorical data such as the type of defect, color, and taste.
The findings in this study align with the concept of continuous improvement and the adaptive monitoring process with a mixed monitoring scheme using the proposed chart. This mixed scheme involves not only one type of quality characteristic but also both types of quality characteristic. Through simulation studies and application to a real case, this scheme was proven effective in monitoring shifts in the process. The decisionmaker in a company can take fast corrective actions for any assignable causes that occur owing to the sensitivity of the proposed scheme. Furthermore, the in-control observations obtained from past monitoring processes need to be stored in a database in order to improve the production process in the future.
Along with the development of the company, the production capability and level of accuracy of the production machinery will increase. As an implication of this development, the monitoring control limit needs to be adjusted. Using the KDE method, the historical incontrol data production is used to generate new control limits by fitting its statistical distributions. This new adjusted control limit will help the company to improve its production quality with stricter criteria for in-control products. Thus, this monitoring scheme guarantees precision and accuracy in monitoring the manufacturing process, and ensures the adaptability of the monitoring system to evolve continuously with developments in science and technology.

Conclusions
In this paper, a new control chart to monitor a process with variable and attribute characteristics is presented using the PCA Mix concept. Simulation studies showed that for several combinations of number of quality characteristics, number of PCs, and values of multinomial parameters, theT 2 control limit using the KDE method always results in an ARL 0 at about 370 for α ¼ 0:00273 . Meanwhile, the ARL 0 of theT 2 control chart using a conventional control limit is not stable. For the shifted process, the ARL 1 of the proposed chart is quickly decreased as the shift from both variable and attribute characteristics increases. TheT 2 control chart has better performance when it uses an appropriate number of principal components. Moreover, the applications for simulation data and real cases also show that the proposed chart is sensitive to monitoring the shifted process. The present work can be extended to learn how the mixed process with both variable and attribute characteristics affects the monitoring system not only in industries and services but also in other fields such as healthcare and network security systems. This study will also trigger the development of multivariate mixed charts for non-normal distribution. In addition, the bootstrap resampling method can be considered to replace the KDE for control limit estimation.