Tolerance Limit-based Estimation of the Proportion of Non-conforming Parts in a Multiple Stream Process

The conventional way to characterize the proportion of non-conforming parts in a process is to calculate process capability indices and transform them into a ratio. These widely used indices are able to give digestible information about the ratio of non-conforming parts if some assumptions are fulfilled. A correct estimation method should be based on the output distribution of the process, and the uncertainty of the parameter estimates should be considered, as well. In this article, a special case of the output distribution is examined: a mixture of normal distributions is considered. In practice, this output distribution appears if a multiple stream process is investigated. The novelty of this study is to apply the tolerance interval-based estimation method for the proportion of non-conforming parts in a case study of a multiple stream process and to qualify the limitations of the proposed estimation method. A simulation study is performed to investigate the bias, mean square error, and root mean square error of the estimates from the two estimation methods (process performance index-based and tolerance interval-based) for different sample sizes for each stream (N ). It was found that, if it may be assumed that the speed of the streams is equal in the case of the sample sizes investigated (N = 25, 50, 100 per head), the proposed (tolerance interval-based) method overestimates the proportion of non-conforming parts while the conventional (process performance index-based) method underestimates it. The tolerance-limit based estimation method has asymptotically better properties than the process performance index-based estimation method.


Introduction 1.Capability indices
Process capability indices (PCIs) are widely used in the manufacturing industry to quantify the capability in different processes [1,2]. These indices measure primarily the relationship between the specification interval (USL-LSL) and the variability of the process (6σ) [3]. The actual aim of the PCIs is to give information about the proportion of non-conforming parts without the explicit use of statistical terms.
One of the most popular indices is the potential capability (C P ): where USL is the upper specification limit, LSL is the lower specification limit, the σ is the square root of the variance of the (according to the conventional interpretation: normal) distribution of the quality characteristic [4]. In this term, the parameters of the distribution are assumed to be known. In the practice, the parameters of the distribution are usually estimated based on a dataset, thus the confidence interval of the C P can be calculated [5]: Where Ĉ P is the estimated value of the C P , χ α To use the model of the C P , some assumptions should be fulfilled [6]: The stability of the process -it means that assignable cause is not present.
• The quality characteristic of interest is a normally distributed random variable with well-defined parameters (μ and σ).
The expected value of the quality characteristic of interests is equal to the midpoint of the specification interval (μ = T = (USL -LSL)/2).
If the last assumption is not fulfilled the C PK index should be used [7]: C . The C PK index gives the so-called demonstrated excellence. It is also possible that only one-sided specification is given. In this case, the C PK is equal to the C PL or the C PU value.
The fractions of non-conforming may be transformed to capability indices [3]. To do this, it is necessary to define the probability model behind the index. As we have seen before, in the case of C P values this default (conventional) probability model is a normal distribution with well-defined parameters. If the quality characteristic is normally distributed and the process is centered in the midpoint of the specifications (USL + LSL)/2, then the proportion of non-conforming parts equals 2Φ(−3C P ), where Φ(•) denotes the standard normal distribution function. If the probability model behind the investigated process is different from normal, then the proportion of non-conforming parts can still be calculated according to this other model; it will correspond to the sum of the areas of the probability density of x below the LSL and above the USL.
In the model of C PK there is a single source of variability, and the quality characteristic of interest is also a normally distributed random variable. But the expected value of the distribution is not equal to the midpoint of the specification interval (T) i.e., the process is shifted [8]. Fig. 1 describes the probability model of C PK where the horizontal axis represents time (t), and the joint distribution is not connected to this axis. The presented Gauss curves here have identical distributions with the same parameters i.e., N(μ,σ S

T
). The question of interest is still the rate of non-conforming parts, however. C PK does not give the rate of non-conforming parts unequivocally as different rates of non-conforming could correspond to the same C PK values [9].
The process capability is often analyzed based on consecutive (and not single) samples. In this case, the variance may be estimated either from the whole dataset or based on the within-sample variance. If the process is not in control, these two estimates are not equal. It should be noted that in Shewhart concept, the process capability calculation makes sense only if the process is in control. In spite of that, in present manufacturing practice, the indices are used if the process is not stable in the Shewhart sense [10]. The long-term variability (σ LT ) is used rather than short-term variability in the concept of process performance index (P P ) [5].
Where σ LT means the standard deviation, which is obtained if the process is operated long run. The usual standard deviation of the sample is used to estimate σ LT .
The σ LT 2 contains both the within-sample variance and between sample variance and these variances are inseparable from each other. The short-term variance -which is used in the denominator of the C P index -is at the most as high as the long-term variance. Thus, the C P index (i.e., the process capability) is always higher or equal to the P P value (i.e., the process performance) [11]. The probability model of the P P (Fig. 2) assumes that there is only one normally distributed random variable which has σ LT 2 variance, and this assumption does not conform to the explanation behind it. In Fig. 2 the narrower normal curves represent short term distributions and the wider one represents the joint one which is considered to be (μ, σ LT 2 ). The horizontal axis is time, but the joint distribution is not connected to this axis i.e., it is not just in a later time.
If the expected value of the quality characteristic is not equal to the target value, the P PK index may be used (similarly to the C PK ): The curve of the joint distribution displayed in Fig. 2 is not the valid output distribution. This model concerns only the P P calculation which is totally incorrect.
As Montgomery referred to the work of Kotz and Lovelace the use of P P is statistical terrorism since specified values of P P and P PK are forced to apply by the supplier [5]. If the process is not in control, there is no way to calculate a simple process capability (or performance) index which would be able to give reliable information about the behavior of the process.
To calculate the proportion of non-conforming items the only sound way is to calculate the area (proportion) of the distribution below the LSL and above the USL values using the correct/appropriate probability model.

Multiple stream processes and capability indices
Multiple stream processes are present in different areas of the industry: filling machines in cosmetics, beverage (food filling machine with multiple heads) or pharmaceutical industry (container filling on a multiple-head machine), production of rubber hoses by extrusion etc. [5,12,13].
In a multiple stream process, the machine has several streams producing in most cases identical number of units. However, the distribution of quality characteristic of the parts is often different for the streams. One source of this variability is a random cause since the quality characteristic is a random variable. The parameters of the distribution of the quality characteristic may change from stream to stream, however. The question of interest is still the non-conforming rate. As the PCIs are used generally, it is obvious to apply them here, as well. Calculating process capability i.e., the non-conforming fraction for multiple stream process is not simple, however. The simplest (and totally incorrect) way would be to calculate process performance (P P or P PK ) indices without considering the probability model behind the process. Fig. 3 illustrates schematically a possible correct probability model of the multiple stream process (with 4 streams, the output rates of the streams are not equal) together with the incorrect joint distribution considered when using the P P index.
The use of average C PK value has been proposed for multiple stream processes by Bothe [14]. This index takes into consideration both the non-centrality of the process and multiple streams (mixture distribution). The computation of the proportion of non-conforming items is based on the percentage of the distribution which is out of the specifications by every stream (every normal distribution). These individual percentages should be averaged and transformed into C PL and C PU values. The minimum from these gives the average C PK index. During this calculation, the following assumptions are made: • Each stream has the same output rate in the production rate. • The probability model is Gaussian for each stream.
• The parameters of the probability models are known.
In addition to all that, the parameters of the distributions are not known (in practice), they are to be estimated from a dataset. As it was mentioned before, the PCIs do not contain explicitly the uncertainty of the parameters. Only a confidence interval can be given for a PCI as a parameter. To define the probability model of a multiple stream  process, imagine a situation that there is a machine with multiple heads (the number of the heads is finite value), the parts coming from these heads are not separated from each other, however. The quality characteristic of this population is a random variable (following normal distribution in the simplest case), the actual value of which depends on the head from which they are originated.
Suppose that x is the quality characteristic of the product which belongs to any head of the machine. The density function of the mixture distribution is [15] where k denotes the identification number of the head, ϕ k (x) is the density function of the quality characteristics (x) of the parts from the k-th head, w k is the output rate of the k-th head. The number of heads is finite and countable.
In this case, the distribution function is a weighted sum, as well: where Φ k (x) is the distribution function of the quality characteristics (x) of the parts from the k-th head, w k is the output rate of the k-th head. For example, if there is a machine with 8 heads and the heads are working with the same speed, the output rate of every head is 1/8. This is the probability that a randomly chosen part is coming from the k-th head of the machine. Let us assume that the Φ k (x) is a distribution function of a normally distributed random variable e.g., with the same variances and different expected values for each head. The parameters of the mixture of univariate normal distributions are [16]: and σ µ σ µ µ Where the X k is the quality characteristic of the part from the k-th head, μ k is the expected value of the normal distribution for the k-th head, σ k 2 is the variance of the normal distribution for the k-th head. From the process capability point of view, the question of interest is the proportion of non-conforming items (out of specification) to the whole mixture distribution (the quality characteristic of a sample/part from either head). Thus, the question concerns the quantile (q) of the population: Where the p denotes the cumulative probability, which belongs to the q quantile of the mixture distribution, N(μ k ,σ k 2 ) means the density function of the k-th component of the mixture distribution.
Other papers suggest the use of global and local PCIs [17], multivariate PCI (in the case mixture of multivariate normal distributions) [18] to evaluate the capability of a multiple stream process. However, the real question is about the tail areas of a supposed distribution which are out of the specification limits. Thus, the novelty of this work is to apply a new, tolerance interval-based estimation method to calculate the proportion of non-conforming parts in multiple stream processes instead of the use of different type of PCIs. After giving the theoretical background of the suggested calculation method (Section 2), the steps of the calculation will be demonstrated with a case study (Section 3) along which the difference between the tolerance interval-based and P P -based method will be discussed (Section 3.1). Finally, a simulation study will be performed in Section 4 to compare the properties of the estimates (bias, standard error, root mean square error) of the two calculation methods.

Tolerance interval-based calculation of the nonconforming fraction for multiple stream processes
As we stated in our former paper [19], the theoretically sound way to calculate the proportion of non-conforming items is connected to the calculation of the tolerance limits [20]. A tolerance interval is a statistical interval that contains at least a specified proportion of the population with defined confidence [21].
For a normally distributed random variable (X) of which the parameters (μ and σ) are known, the one-sided upper tolerance limit can be given exactly with the equation.
x U = μ + Z 1-γ σ. The Z 1-γ denotes the 1 -γ probability quantile of the standard normal distribution. If the parameters of the distribution are not known, the form x U =x ̅ + k s may be used. Here the x ̅ is the sample mean and the s is the standard deviation of the random sample, k is a parameter. The task is to find the k value based on the following equation: It means that the upper x U =x ̅ + k s limit should be given for which at least the 1 -γ proportion of the population is lower with 1 -α confidence. As described by Owen [22], this problem goes back to the non-central t-distribution i.e., Thus, the parameter k can be calculated based on the non-central t-distribution if the number of samples (n), the proportion of the distribution (1 -γ) and the confidence limit (1 -α) are given. The choice of the degree of confidence (1 -α) should be based on the acceptable degree of confidence [21].
As we have seen before, with the calculation of capability indices the aim is to estimate the proportion of non-conforming parts in the process i.e., to determine the proportion of the population that is out of the specification limits. The proportion of non-conforming parts itself is the probability of quantile belonging to the specification limit. However, the parameters of the (mixture) distribution are not known (only estimates are available based on the sample), therefore the quantiles are uncertain. To handle this, during the estimation of the non-conformity rate the one-sided upper (lower) tolerance limit is supposed to be known which is equal to the USL (LSL) value and the proportion of the distribution belonging to this tolerance limits are to be calculated. This calculation method is a reversed situation compared to the usual tolerance interval calculation problems. To use the mentioned estimation method the first step is the calculation of the k parameter: If the parameter k is known, the 1 -γ probability quantile of the non-central t-distribution can be calculated from Eq. (12). The question is the value of the non-centrality parameter for degrees of freedom equal n -1. The calculated non-centrality parameter is divided by n to give the value of z 1-γ . The γ value is found from the Z-table, as the proportion of the z variable larger than the USL value with a certain confidence. The calculation method to define the ratio below the LSL value is the same, but the x ̅ -k's = LSL equation should be used to calculate the k' parameter. Table 1 shows the mass of samples (in g) from the 8 heads of an automatic food filling machine. The dataset is originated from a real industrial process, and it should be noted that we had no access to the real process to investigate it. We found it interesting enough for demonstration, however.

Case study based on real industrial process
To check the normality of the dataset normal Q-Q plots were created. Based on that, one outlier was detected and removed from the original dataset of HEAD1. Thus, in this case, the output rate of the heads differs from each other.
There are no estimates for the parameters of the distribution i.e., the situation is a Phase I study.
The X and MR-charts constructed (Fig. 4) show that there is a drift in the 5th head, thus the process is not in control. Also, other assignable causes are on the third head (samples 2 and 5). In view of we had no access to the real process, we could not investigate it and take corrective action.
For the sake of illustration of the methods, we considered only streams that were in control. Thus, the dataset of HEAD5 and two points from HEAD3 were ignored during the calculation of non-conformity rates. Without the mentioned points it makes sense to analyze the process capability (i.e., the process is considered in control).
The quality characteristic in every head varies according to a normal distribution with well-defined parameters. According to the standard, the maximum allowed negative The standard does not give an upper limit, but it is obviously not acceptable from the producer point of view to exceed a certain limit. Using this consideration, the maximum deviation is chosen as 3% in both directions. The nominal value is 375 g. This analysis aims to calculate the ratio of non-conforming for this process. The lower specification limit (LSL) is 363.75 g, the upper specification limit (USL) is 386.25 g.
On account of we would like to demonstrate the calculations with knowing distributional parameters, as well, these parameters were taken heuristically based on the dataset and are summarized in Table 2.

Calculation of the proportion of non-conforming parts 3.1.1 The case of distribution parameters known Tolerance limit-based calculation method
Assume that the parameters of the distributions are known. As the amount of data is not sufficiently large to assume the parameters as known, we did this just for the sake of illustration. Distributional tests confirmed that the quality characteristic (mass) is normally distributed for each head. In view of the expected values and the variances of the distributions are assumed to be known the calculation is based on the standard normal distribution.
Thus, the z value is calculated for each head with the following formula: If x is substituted with the LSL (or with the USL) value, the result is the z value from which the proportion below the LSL (or above the USL) are found. For example, for the first head: . (15) According to the z-table, the proportion below the LSL value is 0.0262 while the proportion above the USL is 0.0004. The rate of non-conforming parts for the first head is the sum of these probabilities, thus 0.0262 + 0.0004 = 0.0266. As during this calculation, the parameters of the distributions are assumed to be known, this value is considered to be the true rate of non-conforming parts of HEAD1.
Extending this calculation to the other heads, the results i.e., the true rates of non-conforming parts for the other heads are shown in Table 3. The total non-conforming rate is their average, thus 0.0223.

Illustration of the drawbacks of the P P -based estimation method
To illustrate the drawbacks of the P P index in the case of multiple stream processes, we apply it to the dataset even if it is not sound since it does not take into account the mixture distribution of the quality characteristic.
In this case, it is assumed that is the data come from a single distribution (what is not true) with arbitrarily assumed expected value (374.00) and variance (15.63). The expression "arbitrarily assumed" means that the expected value and variance are calculated based on the average of the expected values and variances of single heads. It should be highlighted, that these values of the parameters concern P PL is multiplied by 3 gives the z value, which is 2.592. Based on the Z-table the probability of being below 2.592 is 0.0048. In the same way, the z value belonging to P PU is 3.099, the probability of finding z above this value is 9.71 • 10 -4 . The non-conforming rate is the sum of 0.0048 and 9.71 • 10 -4 , thus equal to 0.0057.
Comparing the two results it is well seen that the result from the P P -based method (0.0057) is significantly lower than the true rate of non-conforming parts (0.0223). Thus, the P P -based method underestimates the non-conforming rate in this case.

The case of distribution parameters not known Tolerance limit-based estimation method
If the parameters of the distributions are not known, as the sample size is not large enough, the uncertainty of the parameters is to be considered, as well. The theoretically sound calculation method is based on the tolerance limit calculation in this case. The question is (still) the probability of the distribution below the LSL and above the USL value. As shown in Section 2 the first task is to calculate the k parameter from the x ̅ + ks = LSL equation. Table 4 contains the averages and the standard deviations for the heads. The out-of-control points are discarded, as in Table 2, but instead of the assumed true parameter values here the parameters of the distributions are estimated, using the number of data points shown in the table. .
. The non-centrality parameter was calculated with qt() function of the R software (version 4.0.3.) giving 7.513 in this case. This value multiplied by n gives the z 1-γ value which is 1.533. From the Z table, the 1-γ is equal to 0.9374, thus the γ is 0.0626. This is the tail area of the distribution which is above the USL value can be found based on the 95% confidence tolerance limit. To give the proportion of the opposite side of the distribution which is under the LSL value (γ') the calculation method is the same. The result is γ' = 0.0874. Thus, for the first head, the estimated proportion of non-conforming items is γ +γ' =0.0626 + 0.0874 = 0.1501.
The estimated rates of non-conforming from the tolerance limit-based method for every head are detailed in Table 5.
The weighted average of these values gives the estimated proportion of non-conforming items for this multiple stream process: it is 0.0746. The P P -based estimation method For the sake of illustration of the incorrect use of the P P index, calculate the proportion of non-conforming items without care of the real structure of the dataset (process) and the uncertainty of the estimated parameters of the distribution. Using this method, the P PK is calculated and for the estimation of the ratio of non-conforming P PU and P PL is used. Averaging for the whole dataset the mean is 374.05 and the standard deviation is 4.09.
Following the steps of the calculation (in Section 4.1.2) presented earlier the estimated value of the proportion of non-conforming parts is 0.0059 + 0.0014 = 0.0073.
Comparing the results of the two estimation methods the capability index-based estimation method gives a much lower result for the proportion of non-conforming  parts. Earlier the true ratio of non-conforming parts was calculated as well (0.0223), thus the bias of the two estimation methods can be obtained. The bias of the P P -based method is lower than the bias of the tolerance limit-based method, but we note that with this method the proportion of non-conforming parts is underestimated. It should be mentioned that these results could be compared only because it was supposed the parameters are known. In practice, the parameters are not known, only better or worse estimates are available depending on the sample.

Simulation study
To investigate the bias, standard error and root mean square error of the estimates by the two calculation methods a simulation study was performed using R (version 4.0.3). First, 7 (from HEAD1-4 and HEAD6-8) normally distributed populations (with one million repetition each) were simulated with the parameters given in Table 2.
Based on the dataset of these simulated populations the true proportion of non-conforming parts was calculated (USL = 386.25, LSL = 363.75). The following steps were repeated 10000 times: random samples of N elements were chosen from the 7 populations, and the proportion of non-conforming parts were calculated with the • tolerance interval-based and • P PK -based methods.
The average of the estimates, the bias (the deviation of the average from the known value of the rate of non-conforming parts), the standard error and the root mean square error of the estimates for N = 25, 50, 100 and 5000 are shown in Table 6.
In the case of the sample sizes investigated the tolerance interval-based method overestimates the rate of non-conforming parts while the P PK based method underestimates it. The absolute value of the bias of the former method is one order of magnitude higher than that of the estimate obtained by the P PK based method. This bias is decreased with increasing the sample size in the case of the tolerance limit-based estimation method and it goes to the true value of the proportion of non-conforming parts. Therefore, the tolerance limit-based estimation method gives asymptotically better estimates than the P PK -based estimation method.
The standard error of the tolerance limit-based method is higher than that of the standard error obtained with the P PK methods in the case of lower sample sizes. The tolerance interval-based method recognizes the uncertainty of the parameters of the distribution contrary to the P PK -based method, where this uncertainty is not considered. Because of this, the variability of the estimate of the tolerance limit-based method is higher.
To give a more reasonable estimate of the proportion of non-conforming parts instead of the conventional P P -based calculation the suggested tolerance limit-based estimation method may be improved. The tolerance-limit based estimation method requires setting the confidence level, but it should be harmonized with the sample size. In this study, the confidence level was fixed at 95%, thus with the harmonization of the sample size and the confidence level, the bias of the proposed method may be reduced.

Conclusions
In this study, two ways to estimate the proportion of non-conforming parts in a multiple stream process were investigated. The estimation methods in the literature apply capability indices to estimate the rate of non-conforming parts in multiple stream processes. Despite that, the method proposed here uses the tolerance limit calculation as the basis of the estimation. According to this, our proposed calculation method considers the uncertainty of the estimated distribution parameters, as well.
A case study with a multiple head's food filler machine was presented for the comparison of the conventional (process performance index-based) and tolerance limit-based estimation methods. In this situation, the process performance index -based method underestimated (significantly) the proportion of non-conforming parts while the tolerance limit-based method overestimated it: the estimated value was almost two times higher than the true rate of non-conforming parts. It means that using the tolerance interval method we are on the safe side. To see the general properties (bias, standard error, root mean square error) of the two estimation methods depending on the sample sizes simulation study was performed. The most important conclusion is that the tolerance limit-based method is able to give better estimates (the bias is lower while the standard error has a similar magnitude) if the sample size is high enough. It seems, the proposed method gives asymptotically better estimates for the proportion of non-conforming parts than the conventional process performance index-based estimation method.
By lower sample sizes the bias of the tolerance limit-based method has a positive sign while the bias of the process performance index-based estimation method is negative. If the acceptance of the conformity of an industrial process during a validation procedure is based on the estimated value of the proportion of non-conforming parts, the overestimation is less dangerous. Thus, the above-mentioned fact is still valid i.e. if the tolerance limit-based estimation method is used for the estimation of the proportion of non-conforming items, we are on the safe side.