Yield-Related Process Capability Indices for
Processes of Multiple Quality Characteristics
Jyh-Jen Horng Shiau,
a*
†Chia-Ling Yen,
aW. L. Pearn
band Wan-Tsz Lee
aProcess capability indices (PCIs) have been widely used in industries for assessing the capability of manufacturing processes. Castagliola and Castellanos (Quality Technology and Quantitative Management 2005, 2(2):201–220), viewing that there were no clear links between the definition of the existing multivariate PCIs and theoretical proportion of nonconforming product items, defined a bivariate Cpkand Cp(denoted by BCpkand BCp, respectively) based on the proportions of nonconforming product items over four convex polygons for bivariate normal processes with a rectangular specification region. In this paper, we extend their definitions to MCpkand MCpfor multivariate normal processes withflexible specification regions. To link the index to the yield, we establish a‘reachable’ lower bound for the process yield as a function of MCpk. An algorithm suitable for such processes is developed to compute the natural estimate of MCpkfrom process data. Furthermore, we construct via the bootstrap approach the lower confidence bound of MCpk, a measure often used by producers for quality assurance to consumers. As for BCp, wefirst modify the original definition with a simple preprocessing step to make BCpscale-invariant. A very efficient algorithm is developed for computing a natural estimator ^BCpof BCp. This new approach of BCpcan be easily extended to MCpfor multivariate processes. For BCp, we further derive an approximate normal distribution for ^BCp, which enables us to construct procedures for making statistical inferences about process capability based on data, including the hypothesis testing, confidence interval, and lower confidence bound. Finally, the proposed procedures are demonstrated with three real data sets. Copyright © 2012 John Wiley & Sons, Ltd.
Keywords: multivariate process capability indices; yield assurance index; normal approximation; lower confidence bound; bootstrap
1. Introduction
P
rocess capability indices (PCIs) are some measures developed for engineering convenience to quantify process performances. Facing stronger-than-ever competitions nowadays, manufacturing companies must assure product quality to their customers to stay competitive. By revealing how well the actual process is in conformance with the manufacturing specifications, PCIs have been widely used in industries as a metric of quality assurance in recent years, and the role is becoming more and more important. Univariate PCIs have been extensively studied in the literature. Some indices such as Cp, Cpu, Cpl, Cpk, Cpm, and Cpmk have shown good values in evaluating univariate processes. See, for example, Kane,1Chan et al.,2Pearn et al.,3Kotz and Johnson,4Kotz and Lovelace,5and Pearn and Kotz.6In many applications, especially in high-technology industries, processes are so complex that the product quality often is affected by multiple characteristics simultaneously. As a result, appropriate multivariate PCIs for assessing processes/products of more than one quality characteristic are desirable. Nonetheless, multivariate PCIs have received comparatively a lot less attention than univariate PCIs in the literature.
In recent years, more and more studies have been devoted to multivariate PCIs. Chan et al.7extended their univariate index Cpmin Chan et al.2to a multivariate version by measuring how far away from the target vector the process mean is in the Mahalanobis distance. Pearn et al.3proposed a multivariate version of Cpand Cpmwith an approach they claimed to be more natural than that of Chan et al.7Hubele et al.8proposed a process capability vector for bivariate normal processes, and later, Shahriari et al.9extended it to the multivariate case. Taam et al.10proposed a multivariate Cpas the ratio of two areas, the area of a modified specification (also called the modified engineering tolerance region by some researchers), defined as the largest ellipsoid centering at the target value and completely within the original specification over the area of the elliptical process region that covers 99.73% of the multivariate normal process. Considering the possible shift of the process mean from the target vector, Taam et al.10further modified this index by taking into account an adjustment factor that measures the closeness between the process mean and the target vector to define a
a
Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan
bDepartment of Industrial Engineering and Management, National Chiao Tung University, Hsinchu, Taiwan *Correspondence to: Jyh-Jen Horng Shiau, Institute of Statistics, National Chiao Tung University, Hsinchu, Taiwan. †E-mail: jyhjen@stat.nctu.edu.tw
Research Article
(wileyonlinelibrary.com) DOI: 10.1002/qre.1397 Published online 17 April 2012 in Wiley Online Library
multivariate Cpm index exactly the same way as Pearn et al.3Recently, Pan and Lee11revised Taam et al.’s modified engineering tolerance region by taking into account the correlation between multiple quality characteristics and proposed two new multivariate PCIs for Cpand Cpm, respectively, which could reflect more correctly the process precision and accuracy. Chen12proposed a multivari-ate PCI using the concept of a tolerance zone, which allowsflexible specifications and no assumptions on the process distribution. Pal13proposed a bivariate PCI as the ratio of the area of the specification rectangle and the 99.73% area of the process region, similar to the index proposed by Taam et al.10Bothe14proposed a multivariate Cpkindex defined as ZP/3, where P is the conforming proportion and ZPis the Pth quantile of the standard normal distribution. Wang et al.15compared three multiple PCIs proposed by Taam et al.,10 Chen,12and Shahriari et al.,9respectively, via graphical and computational examples. Wang and Chen16 and Wang and Du17 proposed multivariate PCIs using principal component analysis (PCA). Recently, Shinde and Khadse18 pointed out that the specification region corresponding to the principal components used in Wang and Chen’s PCI definition was not correct and suggested an alternative method for assessing multivariate process capability based on the empirical probability distribution of principal components. In an earlier work, Shinde and Khadse19reviewed and compared some multivariate PCIs based on fraction conforming interpretation. Gonzalez and Sanchez,20 by relating the actual variability of the process with the prespecified nonconforming proportion, proposed a unitary PCI that can be applied to non-centered univariate processes as well as to general multivariate processes.
Among univariate PCIs, Cpkcould be the most popular one, not only because it accounts for both process mean and variance when assessing the process capability but probably also because it links directly to the process yield by the following inequality given in Boyles:21 2Φ 3Cpk 1⩽%yield⩽Φ 3Cpk ; (1)
where %yield stands for the process yield, andΦ() is the cumulative distribution function (c.d.f.) of the standard normal distribution. Therefore, Cpkis sometimes referred to as a yield assurance index. It is well known that the lower bound, 2Φ(3Cpk) 1, is not a trivial bound; instead, it is a‘reachable’ lower bound in the sense that it can be reached by some processes.
However, in the literature, the link between multivariate PCIs and the product yield was seldom emphasized. One exception was the work of Castagliola and Castellanos.22 Castagliola and Castellanos,22 viewing that there were no clear links between the definition of the existing multivariate PCIs and theoretical proportion of nonconforming product items, defined two indices, BCpk and BCp, based on the proportions of nonconforming product items over four convex polygons for bivariate normal processes with a rectangular specification region. This definition of BCpk is rather interesting because it accounts for the relative position and the orientation of the process distribution with respect to the specification region. Specifically, as an extension of Cpk, BCpkquantifies the process capability based on the smallest conforming proportion of the four convex polygons formed by dividing the rectangular specification region with the two main axes (i.e., the principal components) of the process distribution.
As defined by Castagliola and Castellanos,22the relationship between BC
pkand BCpis analogous to that between the univariate Cpk and Cp. Specifically, because the notion of Cponly concerns with the process variability, BCpwas defined as the maximum BCpkvalue of all bivariate normal distributions with the same covariance matrix as the process covariance matrixΣ, which is analogous to the univariate case. Although this definition works for the univariate case, unfortunately it fails the scale-invariance property for multivariate processes. As a result, the value of BCpis not consistent when quality characteristics are measured in different units or scales.
The main purpose of this paper is to extend the notion of BCpk and BCp to multivariate PCIs (denoted by MCpk and MCp, respectively) for multivariate processes that may have more than two quality characteristics. We first define MCpk and then establish the same‘reachable’ lower bound for the process yield in terms of MCpkas in (1). Because the computation method of Castagliola and Castellanos22 is only for bivariate processes with rectangular specification regions, we further develop a computation method that can be implemented for multivariate processes withflexible specification regions. As for BCp, we modify the definition of Castagliola and Castellanos22 with a simple preprocessing step; then BC
pbecomes scale-invariant. Also, because the computation method provided in Castagliola and Castellanos22for ^BC
p, a natural estimate of BCp, is very time-consuming, a very efficient algorithm is developed. Moreover, we derive an approximate distribution for ^BCp, which enables us to provide statistical procedures for making inferences about process capability based on data, including hypothesis testing, confidence interval (CI), and lower confidence bound. The results of statistical inferences are very useful in decision making. In particular, the lower confidence bound is a measure of high practical value because it directly links to the quality assurance. Finally, we also extend BCpto MCpfor multivariate processes.
The rest of the paper is organized as follows. In Section 2, wefirst review the index BCpkproposed by Castagliola and Castellanos22 for bivariate processes. Then, we establish the link between this index and the product yield. We further extend this index to multivariate processes of more than two characteristics. After that, we give an algorithm for the estimation of MCpkand propose obtaining lower confidence bounds by bootstrap methods. For demonstration, we apply the methods to simulated examples in the bivariate case. We also study the distribution of a natural estimator of BCpkby simulation. In Section 3, for BCp, we show how to obtain a scale-invariant BCpand how and why ^BCp can be efficiently calculated. We further derive an approximate normal distribution for ^BCp and then use it to construct statistical procedures for inferences about process capability. In Section 4, as illustrative examples, we apply the proposed indices and inference procedures to the two sets of real (bivariate) data presented in Castagliola and Castellanos22and a trivariate dataset obtained from a stencil printing process described in Pan and Lee.11Finally, we conclude the paper with a brief summary and remark in Section 5.
2. Multivariate C
pkindex–a yield measuring process capability index
2.1. Alternative definition for Cpk
Assume that the quality characteristic X of a product item follows a normal distribution with meanm and variance s2, denoted by N(m, s2). Let [LSL, USL] be the corresponding lower and upper specification limits. The usual definition of Cpkas defined by Kane1is
Cpk¼ min USL m 3s ; m LSL 3s ; (2)
which accounts for not only the spread of the process but also the location of the process mean relative to the specification limits. Equivalent to (2), an alternative definition for Cpkwas proposed by Castagliola and Castellanos.22This definition is based on the lower and upper proportions of nonconforming product items, pL= P(X⩽ LSL) and pU= P(X⩾ USL), as follows. Because X follows N(m, s2), p L¼ Φ mLSLs and pU¼ Φ USLms . Thus, Cpkis equivalent to 1 3min Φ 1 p L ð Þ; Φ1 p U ð Þ ; (3)
because of the fact that the c.d.f.Φ() is a strictly increasing function. Similarly, the usual definition of Cp= (USL LSL)/6s proposed by Kane1is equivalent to 1 6 Φ 1 p U ð Þ Φ1ð ÞpL : 2.2. Castagliola and Castellanos’ definition of BCpk
Let X1and X2be the quality characteristics of interest with the specification limits [LSL1, USL1] for X1and [LSL2, USL2] for X2. These limits define a rectangular specification area. Assume that X = (X1, X2)Tfollows a bivariate normal distribution with meanm = (m1,m2)Tand variance-covariance matrixΣ. Applying eigenvalue–eigenvector decomposition to Σ, one can obtain two eigenvalues, l21⩾l22> 0, and the associated eigenvectors, v1and v2. Let R = [v1, v2]. Then RTR = I andΣ can be expressed as Σ = RV RT, where V is the diagonal matrix with diagonal elements l21 andl22. In fact, the matrix R represents the rotation matrix that rotates the original axes to the main axes of the bivariate normal distribution (see Figure 1); the directions of the eigenvectors v1and v2correspond to the main axes; and l21 and l22 are the variances of X projected onto the two main axes, respectively. More specifically, if we let Si¼ viTX , then Si N viTm; l
2 i
for i = 1, 2; moreover, S1and S2are independent. In the PCA context, v1and v2are called principal components, and S1and S2are called principal component scores.
Suppose we move the origin to the process meanm and rotate the two original axes to the directions of v1and v2. Then, the two new axes divide the plane into four regions, A1, A2, A3, and A4. Obviously, P(X2 Ai) = 1/4 for i = 1, 2, 3, 4. Denote the specification region by A. Let Qi= Ai∩ A and qi= P(X2 Qi) for i = 1, 2, 3, 4. Then, the probability that X is in Aibut not in the specification region A is pi= 1/4 qi. See Figure 1. In other words, pis are the proportions of nonconforming product items in the four quadrants of the new coordinate system.
By analogy to the alternative definition of Cpkgiven in (3), Castagliola and Castellanos 22 defined a bivariate Cpkas BCpk¼ 1 3min Φ 1 2p 1 ð Þ; Φ1 2p 2 ð Þ; Φ1 2p 3 ð Þ; Φ1 2p 4 ð Þ : (4)
This definition is similar to the univariate case as in (3), except that 0 ⩽ pi⩽ 1/4 for i = 1, 2, 3, 4 in the bivariate case, whereas 0 ⩽ pU, pL⩽ 1/2 in the univariate case. We will extend this definition to higher dimensions later.
1 LSL USL1 2 LSL 2 USL 2 X 1 X 4 Q 1 Q 3 Q 1 v 2 v 2 Q
Figure 1. Explaining diagram of BCpk
2.3. Nonconforming rate based on BCpk
According to the definition of BCpkin the last subsection, we can establish a link between the nonconforming rate (denoted by %NC) and BCpkas follows. First note that, by (4),
BCpk¼ 1 3max Φ 1 2p 1 ð Þ; Φ1 2p 2 ð Þ; Φ1 2p 3 ð Þ; Φ1 2p 4 ð Þ : BecauseΦ 1() is a strictly increasing function, we have
BCpk¼ 1 3Φ 1 2p max ð Þ; (5)
where pmax= max{p1, p2, p3, p4}. Then, by (5),
pmax¼ 1
2Φ 3BCpk
: (6)
Note that pmax⩽ %NC ⩽ 4pmax. Plugging (6) into this inequality, we obtain 1 2Φ 3BCpk ⩽%NC⩽2Φ 3BCpk : (7)
The upper bound in (7) is very useful and is not a loose bound, meaning that it is reachable for some processes. Usually, producers can take this upper bound as a metric of quality assurance to customers. For example, if the process is with BCpk= 1.00, one can guarantee that there will be at most 2700 non-conformities in 1,000,000 product items. On the other hand, the lower bound in (7) is quite conservative; nevertheless, it is a convenient bound, meaning when a practitioner obtains a BCpkvalue, a lower bound of the nonconforming rate as such is immediately available to him/her.
Table I gives the upper and lower bounds of the nonconforming rate %NC for various values of BCpk. Figure 2 plots the bounds. We can see the bounds drop sharply as BCpkincreases and soon carried out to near zero level when BCpk≥ 1.5.
The second inequality of (7) is equivalent to
2Φ 3BCpk
1⩽% yield; (8)
which provides a‘reachable’ lower bound for the yield. Note that this lower bound is the same as that in (1) for the univariate case. When the BCpkvalue of the process is available, producers can assure the yield level with this lower bound to their customers. 2.4. Extending BCpkto higher dimensions
We now generalize the alternative definition of BCpkto multivariate processes of k characteristics where k> 2. By the same notion for the bivariate case, dividing the Euclidean space Rkinto 2khyperquadrants by the k main axes (i.e., principal components) of the k-variate distribution, we can define a multivariate Cpkindex as
MCpk¼ 1 3min Φ 1 2k1p 1 ; Φ1 2k1p 2 ; . . . ; Φ1 2k1p 2k ¼ 1 3max Φ 1 2k1p 1 ; Φ1 2k1p 2 ; . . . ; Φ1 2k1p 2k ¼ 13Φ1 2k1pmax ; (9)
where piis the probability of a randomly selected sample that is in the ith hyperquadrant but not in the specification region for i = 1, 2, , 2kand p
max¼ max p1; p2; . . . ; p2k
. Equivalently,
Table I. Bounds of non-conformities based on BCpk
BCpk
Non-conformities (in ppm)
Lower bound Upper bound
0.60 17965.15956 71860.63823 0.80 4098.76796 16395.07185 1.00 674.94902 2699.79606 1.33 16.51832 66.07330 1.50 1.69884 6.79535 1.60 0.39666 1.58666 1.67 0.13608 0.54430 2.00 0.00049 0.00197
490
pmax¼ 1
2k1Φ 3MCpk
:
Because pmax⩽ %NC ⩽ 2kpmax, we can also obtain bounds for the nonconforming rate in the general multivariate case as 1 2k1Φ 3MCpk ⩽%NC⩽2Φ 3MCpk ; (10) which is equivalent to 2Φ 3MCpk 1⩽%yield⩽1 1 2k1Φ 3MCpk : (11) 2.5. Estimation of MCpk
With the definition of MCpkgiven in (9), a natural estimator of MCpkis ^ MCpk¼ 1 3Φ 1 2k1^p max ; where^pmax¼ max ^p1; ^p2; . . . ; ^p2k
with^pibeing the sample nonconforming proportion in the ith hyperquadrant for i = 1, 2,. . ., 2k. The method to calculate ^pi and^qi proposed by Castagliola and Castellanos22is a fairly complicated integration method that turns two-dimensional integrations over convex polygons into line integrals based on Green’s formula, which cannot be directly extended to higher dimensions. Next, we propose an algorithm to calculate MC^ pk from a sample of n k-dimensional quality characteristic vectors, x1, x2,. . ., xn.
Algorithm for Calculating ^MCpk:
1. Estimatem and Σ by the usual sample mean and sample covariance matrix as
^m ¼1 n Xn j¼1 xj and ^Σ ¼ 1 n 1 Xn j¼1 xj ^m xj ^m T :
2. Compute eigenvalues ^l21; ^l22; . . . ; ^l2kand eigenvectors^v1; ^v2; . . . ; ^vkof ^Σ. (Denote by A1; A2; . . . ; A2kthe 2khyperquadrants of the Euclidean space Rkusing the center^m as the origin and the k orthogonal directions ^v1; ^v2; . . . ; ^vkas the coordinate axes.) 3. Compute an estimate ^pi of pi for i = 1, 2,. . ., 2k by Monte Carlo integration as follows. Generate {X1, X2,. . ., XN}, a very
large number of data from Nk
ð
^m; ^ΣÞ
. For each Xj, j = 1, 2,. . ., N, determine which hyperquadrant it belongs to by the signs of Xj ^mT
^vl; l ¼ 1; 2; . . . ; k
n o
. Compute the proportion of Xjs that are in the intersection of the specification region A and the ith hyperquadrant Aiby
0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 1 2 3 4 5 6 7 x 104 pk NC (ppm) upb lwb
Figure 2. Bounds of non-conformities based on BCpk
^qi¼ ♯ Xj2 Ai∩A N (12) and compute ^pi¼ 1 2k ^qi for i¼ 1; 2; . . . ; 2k:
4. Compute^pmax¼ max ^p1; ^p2; . . . ; ^p2k
and then the estimate for MCpkby
^ MCpk¼ 1 3Φ 1 2k1^p max :
We remark that despite requiring intensive computation for higher dimensions because of the curse of dimensionality, this Monte Carlo method works for all dimensions and any shapes of the specification region.
However, when studying the distribution of ^MCpkempirically, we need to repeat the estimation procedure for each realization of ^
MCpk. Then the algorithm described earlier becomes computationally infeasible because each replication would need the generation of a huge amount of data from its own Nk
ð
^m; ^ΣÞ
, say, N = 1,000,000 or even 10,000,000 in the bivariate case, for Monte Carlo integration; and a large number of replications are needed to obtain a good approximation. To overcome this computation difficulty, we develop a procedure that requires generating N data from the standard multivariate normal distribution Nk(0, I) only once—which in turn can be carried out by simply generating N k data from the (univariate) standard normal distribution—as described in the following.BecauseΣ is assumed symmetric positive definite, there exists a unique symmetric positive definite matrix Σ1/2such thatΣ = (Σ1/2) (Σ1/2) (Golub and Van Loan,23 p. 395). To simplify the notation, (Σ1/2) 1 is denoted by Σ 1/2. It is well known that the affine transformation of Z =Σ 1/2(X m) transforms a random vector X following Nk(m, Σ) to a vector variate Z following the standard multivariate normal distribution Nk(0, I). With this, we can just generate N vector variates, {Zj, j = 1, 2,. . ., N} only once from Nk(0, I) and reuse them for all replications. Specifically, for a replication with sample mean ^m and sample covariance matrix ^Σ, we transform the specification region by the transformation T ð Þ ¼ Σ1=2ð ^mÞ . When the specification region is a rectangle (or cube), we only need to transform the vertices and then reconstruct the specification region in the transformed space. Then we can compute ^qis in (12) by ^qi¼ ♯ Zj2 A′i∩A′ N for i¼ 1; 2; . . . ; 2 k;
where A′i is the ith hyperquadrant of Rkin the standard coordinate system, and A′ is the transformed specification region T(A). For simplicity, we illustrate our method by examples in the bivariate case with rectangular specifications. The rectangular specification is the most widely used shape in real-life applications. Table II lists the distribution parameters and the specifications of four examples. For illustration, Figure 3 plots, for each case, a set of sample data with size n = 100, the rectangular specification region, and the two orthogonal lines crossed at the sample mean with the eigenvectors of ^Σ as their directions.
To evaluate how well ^BCpkestimates BCpk, we generate 1000 sets of data with size n = 100 for each case and compute ^BCpkfor each set of data with the aforementioned algorithm using N = 1, 000, 000. Tables III–VI present, respectively, for Cases 1–4, the true values of q1, q2, q3, q4, and BCpkas well as the sample mean and sample standard deviation of 1000 values of^q1,^q2,^q3,^q4, and ^BCpk. The bias defined as the difference between the sample mean and the true value is also included. The results indicate that ^BCpkis a reasonable estimator.
Table II. Parameters and specifications of four bivariate normal examples
Distribution parameters X1spec X2spec
m1 s21 m2 s22 r LSL1 USL1 LSL2 USL2 Case1 6.0 0.8 7.0 1.0 0.0 2.0 10.0 3.0 10.0 Case2 5.0 0.5 6.0 0.45 0.5 2.0 8.0 3.0 8.0 Case3 3.0 1.0 6.0 1.0 0.2 0.5 6.5 1.0 7.0 Case4 1.0 1.0 1.0 1.0 0.2 1.0 5.0 1.0 3.0
492
0 2 4 6 8 10 12 0 2 4 6 8 10 12 X1 X2 Case 1 0 2 4 6 8 10 12 0 2 4 6 8 10 12 X1 X2 Case 2 −2 0 2 4 6 8 0 2 4 6 8 10 X1 X2 Case 3 −2 0 2 4 6 8 −2 0 2 4 6 8 X1 X2 Case 4
Figure 3. Illustration of the four cases in Table II, displaying 100 data points, the corresponding principal components, and the specification region
Table III. Estimation results of Case 1
True value Sample mean Sample SD Bias
q1 0.2493360 0.2491962 0.0006310 0.0001398
q2 0.2493280 0.2491763 0.0006431 0.0001517
q3 0.2499820 0.2499745 0.0000264 0.0000075
q4 0.2499860 0.2499727 0.0000302 0.0000133
BCpk 1.0004445 0.9969328 0.07807013 0.0035117
Table IV. Estimation results of Case 2
True value Sample mean Sample SD Bias
q1 0.2499920 0.2492297 0.0136783 0.0007623
q2 0.2499970 0.2492439 0.0136790 0.0007531
q3 0.2499870 0.2492214 0.0136778 0.0007656
q4 0.2499990 0.2492467 0.0136792 0.0007523
BCpk 1.3488170 1.3100525 0.1081182 0.0387645
Table V. Estimation results of Case 3
True value Sample mean Sample SD Bias
q1 0.1613160 0.1606882 0.0177669 0.0006278 q2 0.1784460 0.1778262 0.0143900 0.0006198 q3 0.2462340 0.2459017 0.0023246 0.0003323 q4 0.2499080 0.2498783 0.0001174 0.0000297 BCpk 0.3084807 0.3094362 0.0459877 0.0009555
493
2.6. Statistical inference about MCpkvia bootstrap approach
For statistical inference about any PCI, practitioners often emphasize the lower confidence bound of the PCI rather than the usual CI, viewing from the quality assurance aspect. Among PCIs, lower confidence bounds of indices in the Cpk-family are of particular interests for practitioners because the yield corresponding to the lower confidence bound represents the worst yield at a certain confidence level. So far as we know, no lower confidence bound has been developed for the Cpk index for processes of multiple characteristics in the literature.
Note that, with a set of process data, we can only obtain an ^MCpk. To infer anything about MCpk, for example, its lower confidence bound, usually, we would need the distribution of ^MCpkor have many ^MCpks. Unfortunately, the distribution of ^MCpkis analytically intractable, and repeating experiments to obtain a number of estimates is not possible or economical for most of the applications. Under such circumstances, the bootstrap approach introduced by Efron24 is commonly used for statistical inferences. With the bootstrap method, one can repeat the resampling procedure many times to obtain‘bootstrap’ estimates of the parameter of interest without specific model assumptions to infer about the population parameter.
We emphasize here that the purpose of the bootstrap is not to obtain a better parameter estimate because the bootstrap distribution is always centered around the statistic calculated from the data (here the value of ^MCpk), not the unknown population value (MCpk). Rather, the bootstrap is useful and convenient for quantifying the behavior of a parameter estimate ( ^MCpk), for example, obtaining its standard error/bias, or calculating CIs.
In this study, for making inferences on MCpk, we also choose the bootstrap rather than other approaches to deal with this problem because it is simple, convenient to perform, and having sound statistical justifications for statistical inferences.
The bootstrap procedure can be briefly described as follows. Suppose that we have a random sample {X1, X2,. . ., Xn} of size n from a population with the c.d.f. Fθ, in whichθ is the parameter of interest. Let ^θ be an estimator of θ. Denote by X 1; X2; . . . ; Xn a resampled data set of size n obtained by resampling with replacement from {X1, X2,. . ., Xn}. Then, Xjs are independent and identically distributed following the empirical distribution with the c.d.f. Fnð Þ ¼x 1nPni¼11 Xð i⩽xÞ. Calculate ^θ with this resampled data set and denote it by ^θ. Repeat this for B times, and obtain ^nθj; j ¼ 1; 2; . . . ; Bo. Then, we can have inferences aboutθ based on the bootstrap estimates ^θ1;^θ2; . . . ;^θB.
In this section, taking BCpkas an illustrative example, we describe how to obtain lower confidence bounds for MCpkvia various bootstrap methods, including the basic bootstrap method, percentile bootstrap method, standard bootstrap method, and bias-corrected percentile bootstrap method. For more details about bootstrap methods, see, for example, Davison and Hinkly,25Efron,24 Efron and Tibshirani,26and Carpenter and Bithell.27
1. Basic bootstrap method
Following Davison and Hinkly,24one can obtain an approximate 100(1 a)% CI of θ by the basic bootstrap method as 2^θ ^θB 1a 2 ð Þ ½ ð Þ; 2^θ ^θð½Bð Þa2Þ ;
where [x] stands for the nearest integer of x, ^θð Þi is the ith ordered estimate from the bootstrap procedure, and ^θ is the estimate from the original sample. Analogously, an approximate 100(1 a)% lower confidence bound of the basic bootstrap method is 2^θ ^θð½B 1að ÞÞ.
2. Standard bootstrap method
The average and standard deviation of B bootstrap estimates ^θ1;^θ2; . . . ;^θBare θ ¼1BX B i¼1^θ i and S¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 B 1 Xn i¼1 ^θ i θ 2 s ;
respectively. One can use normal approximation to obtain an approximate 100(1a)% CI of θ based on the standard bootstrap method as θ Za 2S ; θ þ Za 2S
and the corresponding approximate 100(1a)% lower confidence bound as θ ZaS. Table VI. Estimation results of Case 4
True value Sample mean Sample SD Bias
q1 0.2364030 0.2290143 0.0098276 0.0073887 q2 0.0139380 0.0179137 0.0150866 0.0039757 q3 0.0000060 0.0198540 0.0008475 0.0003036 q4 0.0159930 0.0157176 0.0165061 0.0038610 BCpk 0.0000100 0.0005120 0.0014184 0.0005020
494
3. Percentile bootstrap method
The percentile bootstrap method simply takes the sample 100(a/2) and the 100(1 a/2) percentage points as the confidence bounds to construct an approximate 100(1 a)% CI as
^θ Bð Þa2 ½ ð Þ; ^θB 1a 2 ð Þ ½ ð Þ
and the corresponding approximate 100(1 a)% lower confidence bound as ^θð½Bð ÞaÞ. 4. Bias-corrected percentile bootstrap method
It is possible that the bootstrap distribution obtained using only a sample of the complete bootstrap distribution may be shifted higher or lower than would be expected; thus, the bias-corrected percentile (BCP) bootstrap method was suggested by Efron and Tibshirani25to correct this bias. First, using the distribution of ^θ, calculate the probability
p0¼ P ^θ
⩽^θ
by the proportion of ^θis satisfying ^θi⩽^θ. Second, calculate
z0¼ Φ1ð Þ;p0 PL;a=2¼ Φ 2z 0 za=2; PU;a=2¼ Φ 2z 0þ za=2: Finally, an approximate 100(1 a)% CI obtained by the BCP bootstrap method is
^θ BPL;a=2
½
ð Þ; ^θð½BPU;a=2Þ
and the corresponding approximate 100(1 a)% lower confidence bound is ^θ BPL;a
½
ð Þ.
These four bootstrap methods have their own advantages and disadvantages, see Carpenter and Bithell.26Briefly speaking, the percentile method is simple to calculate and often works well, especially when the sampling distribution is symmetrical; however, it may not have the correct coverage when the sampling distribution is skewed. Both basic bootstrap and standard bootstrap CIs require statistics with small bias and sampling distributions close to normal. When bias or skewness is present in the bootstrap distribution, one can use BCP bootstrap method; the basic, standard, and percentile intervals are inaccurate under these circumstances unless the sample sizes are very large. Carpenter and Bithell26provided a practical guide on the use of bootstrap CIs, including guides on when they should be used, which method should be chosen, and how they should be implemented. Which bootstrap methods are proper for a particular application depends on the sampling distribution, sample size, bootstrap number, or even the original data.
For demonstration, we apply these four bootstrap methods to Case 1 example given in Table II, which has a BCpkvery close to 1.00. Consider various sample sizes n = 30(10)100, 125(25)200, 250, 300. For each generated sample, we use the algorithm described earlier to obtain a ^BCpk. Here, for the Monte Carlo integration, we generate N = 1, 000, 000 data from N2(0, I). We then perform the bootstrap resampling B = 3000 times to obtain 3000 bootstrap estimates of BCpk. With these 3000 estimates, we obtain a lower confidence bound (LCB) for each of the four bootstrap methods. To see the performance of the proposed LCB, we repeat these steps for 300 times to obtain 300 LCBs.
To compare the four bootstrap methods, we consider three criteria: (i) the mean of LCB (the closer to the nominal value the better); (ii) the standard deviation of LCB (the smaller the better); and (iii) the coverage probability (the closer to the confidence level the better). Each of these three criteria has its own merit. As their sample version, Table VII lists the sample mean, sample standard deviation, and coverage rate of the simulated 300 LCBs obtained at 90% confidence level under various sample sizes n for each of the four bootstrap methods. To help us read the simulation results, for each sample size n, we highlight the mean LCB closest to 1.00, the smallest standard deviation, and the coverage rate closest to 90% among the four methods; we also‘teletype’ the worst ones. From Table VII, we observe the following:
• It is clearly seen that the mean LCB becomes closer to the nominal value 1.00 and the standard deviation becomes smaller as data size n becomes larger for all bootstrap methods.
• The percentile method performs the best, and the basic method performs the worst under the criteria (i) and (ii); more specifically, the performance orderings of the four methods are percentile > standard > BCP > basic under both criteria. • The coverage rates of the four methods are reasonable in general, but the ordering is mixed, no apparent patterns. The basic
bootstrap method seems performing the worst (especially when n is large), whereas the other three methods seem not too much different.
We remark that the ordering based on empirical LCB values is not necessarily the actual ordering because of the estimation error. In summary, all four methods are reasonable methods with the percentile method being the best and the basic method being the worst. Therefore, we recommend the percentile bootstrap method for computing the LCB of MCpk.
2.7. Empirical distribution function of ^BCpk
To acquire some idea about the sampling distribution of ^BCpk, we simulate 100, 000 ^BCpks for Case 1 (BCpk 1.00) and Case 2 (BCpk 1.33) examples in Table II by the algorithm described in Section 2.5 with sample data of size n = 500 and using N = 10, 000, 000 simulated N2(0, I) data for Monte Carlo integration.
Figure 4 displays the distributions of 100, 000 ^BCpks in the form of histograms for these two cases, which look fairly normal-like. To see if ^BCpkbehaves similarly to a normal distribution, we further calculate the empirical cumulative distribution function of the simulated ^BCpks and compare it with a normal distribution for each case. See Figure 5. In addition, Figure 6 presents the Q-Q plots of the 100,000 simulated ^BCpks for both cases. These plots suggest that the sampling distribution of ^BCpkis fairly close to a normal distribution.
3.
Multivariate C
pindex
—a variation measuring process capability index
3.1. Bivariate Cpindex: BCp
In the univariate case, the index Cpis defined as
Cp¼
USL LSL
6s ;
which only accounts for the process variation and totally ignores the location of the process meanm. So Cpis sometimes referred to as a variation measuring index. By their definitions, Cpk⩽ Cpand the equality holds only whenm = (LSL + USL)/2, that is, when the process is well centered in the specification region, the most desirable location.
Table VII. The sample mean, sample standard deviation, and coverage rate of 300 lower confidence bounds at 90% confidence level for an example with BCpk 1
Data size Sample mean Sample SD Coverage rate (%) Sample mean Sample SD Coverage rate (%)
Basic Standard 30 0.812937 0.16718 86.33 0.856004 0.13526 86.33 40 0.828009 0.13559 89.33 0.863823 0.11560 89.00 50 0.846850 0.12064 87.00 0.879382 0.10116 86.00 60 0.855107 0.10840 88.67 0.880893 0.09870 89.33 70 0.862006 0.09320 91.67 0.886200 0.08901 87.67 80 0.869648 0.08790 90.33 0.890566 0.08171 88.00 90 0.884494 0.08394 91.67 0.903152 0.07898 90.33 100 0.886191 0.07873 93.00 0.904637 0.07520 90.67 125 0.899050 0.06851 94.33 0.914252 0.06745 89.67 150 0.906501 0.06354 92.33 0.918361 0.06206 89.33 175 0.916951 0.05776 92.00 0.925228 0.05672 90.33 200 0.926768 0.05584 92.00 0.935849 0.05548 90.33 250 0.928409 0.04986 91.33 0.937738 0.04836 90.00 300 0.934398 0.04385 93.00 0.938921 0.04360 93.33
Percentile Bias-corrected percentile
30 0.862382 0.13089 86.33 0.853172 0.14994 85.00 40 0.869565 0.11272 89.33 0.860984 0.12254 88.33 50 0.883955 0.09887 85.33 0.874425 0.10969 86.00 60 0.885301 0.09680 89.33 0.878285 0.10110 88.67 70 0.890589 0.08901 87.33 0.882519 0.09169 89.00 80 0.894095 0.08039 87.67 0.888119 0.08327 89.33 90 0.906300 0.07792 90.33 0.901452 0.08044 90.67 100 0.907816 0.07405 90.33 0.902246 0.07539 91.67 125 0.916704 0.06675 89.67 0.912058 0.06693 92.33 150 0.920422 0.06145 89.00 0.917068 0.06243 89.67 175 0.926494 0.05673 90.67 0.923867 0.05735 91.33 200 0.937575 0.05518 88.67 0.935379 0.05571 89.00 250 0.938241 0.04828 89.67 0.937381 0.04846 89.67 300 0.940157 0.04361 92.00 0.940005 0.04439 92.33
Bold emphasis indicates the (best) method with the lower confidence bound closest to 1.00, the least standard deviation, and the coverage rate closest to 90% under the three criteria, respectively; on the other hand, the teletyped numbers indicates the worst method among the four methods.
Because of this, Castagliola and Castellanos22defined a new bivariate Cpindex, BCp, as the maximum value of BCpkover the process meanm and the angle θ of the rotation matrix R that rotates the original axes to the two principal components (i.e., eigenvectors of Σ) as described earlier in Section 2.2. Specifically,
BCp max
m;θ BCpk: (13)
Because BCpkitself is defined as the minimum of the four values, Φ 1(2p1), Φ 1(2p2), Φ 1(2p3), and Φ 1(2p4), the maximum value of BCpk is necessarily reached when Φ 1(2p1) = Φ 1(2p2) = Φ 1(2p3) = Φ 1(2p4), that is, when p1= p2= p3= p4 and p p1+ p2+ p3+ p4is minimum; or equivalently when q1= q2= q3= q4and q q1+ q2+ q3+ q4is maximum. As a result,
0.85 0.9 0.95 1 1.05 1.1 1.15 0 2 4 6 8 10 12 pk pk (a) Density Histogram Normal Plot 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 0 2 4 6 8 10 12 (b) Density Histogram Normal Plot
Figure 4. Histograms of 100,000 ^BCpks with a normal curve for (a) BCpk 1.00 (b) BCpk 1.33
0.85 0.9 0.95 1 1.05 1.1 1.15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) ECDF Empirical CDF ECDF Normal 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (b) ECDF Empirical CDF ECDF Normal pk pk
Figure 5. Comparing the empirical cumulative distribution function (c.d.f.) of 100,000 ^BCpks with a normal distribution for (a) BCpk 1.00 and (b) BCpk 1.33
BCp¼ 1 3Φ 1 2p 4 ¼ 1 3Φ 1 p 2 ¼ 1 3Φ 1 1 q 2 : (14)
Then, it remains tofind the m and θ such that the corresponding process has the previously mentioned property. Castagliola and Castellanos22proposed a solution (in fact, a procedure) by lettingm1= (LSL1+ USL1)/2,m2= (LSL2+ USL2)/2 and varying the rotation angleθ of the rotation matrix R to search for the optimal p. The definition and the procedure seem reasonable intuitively.
Unfortunately, the BCpdefined as in (13) is not scale-invariant, in the sense that the value of BCpvaries if we rescale the process variables X1and X2, or simply use another unit for the variables. Let us take the four examples in Table II to demonstrate this scaling problem. We rescale each case by X1′ ¼ 2X1and X′2¼ 3X2. Table VIII lists the parameters and specifications of each case after scaling. Table IX presents the values of BCpkand BCpalong with the corresponding values ofθ and p obtained by Castagliola and Castellanos’s
−5 0 5 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25
Standard Normal Quantile (a)
Quantile
QQ Plot of Sample Data versus Standard Normal
−5 0 5 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6
Standard Normal Quantile (b)
QQ Plot of Sample Data versus Standard Normal
pk
Quantile
pk
Figure 6. Q–Q plots of 100,000 ^BCpks for (a) BCpk 1.00 and (b) BCpk 1.33
Table VIII. Parameters and specifications of four examples after scaling
Distribution parameters X1′ spec X2′ spec
m1 s21 m2 s22 r LSL1 USL1 LSL2 USL2
Case1′ 12.0 3.2 21.0 9.0 0.0 4.0 20.0 9.0 30.0
Case2′ 10.0 2.0 18.0 4.05 0.5 4.0 16.0 9.0 27.0
Case3′ 6.0 4.0 18.0 9.0 0.2 1.0 13.0 3.0 21.0
Case4′ 2.0 4.0 3.0 9.0 0.2 2.0 10.0 3.0 9.0
Table IX. BCpand BCpkvalues of four examples before and after scaling
BCpk BCp Case 1 1.0004445 1.2528016 θ = 90∘, p = 0.0000428 Case 1′ 1.0004445 1.2098245 θ = 25∘, p = 0.0000710 Case 2 1.3488170 1.3862485 θ = 0∘(or 90∘,180∘, 270∘), p = 0.0000080 Case 2′ 1.3488170 1.393727 θ = 170∘, p = 0.000007 Case 3 0.3084807 0.9257124 θ = 0∘(or 90∘,180∘, 270∘), p = 0.0013710 Case 3′ 0.3084807 0.9334442 θ = 170∘, p = 0.0012763 Case 4 0.0000100 0.3353143 θ = 45∘, p = 0.0786108 Case 4′ 0.0000100 0.3550567 θ = 75∘, p = 0.0716998
498
procedure for all four cases. From Table IX, it is clearly seen thatθ, p, and BCpchange their values after rescaling whereas BCpkstays the same. The scale-invariance property of BCpk is apparent because rescaling will not change the probability of the quality characteristic vector X being in each specification polygon Qi.
In real practice, it is very common for people to use different units for quality characteristics. Any process assessment scheme definitely should not be affected by the unit used, which means, mathematically, a well-defined capability index should be invariant of scaling. Here, we propose a simple solution tofix this problem: rescale the data and the specifications such that the specification rectangle becomes a square centered at the origin (0,0), which in some sense is to‘equalize’ the importance of the variables. Let the specification region be the rectangle [LSL1, USL1] [LSL2, USL2]. As an example, we can transform the quality characteristic vector X = (X1, X2)Tinto X′ ¼ X1′; X2′ T by X1′ ¼ 1 USL1LSL1 X1 USL1þLSL1 2 and X2′ ¼ 1 USL2LSL2 X2 USL2þLSL2 2
. Then, the specification rectangle is transformed into the unit square 1
2; 1 2 1 2; 1 2
. In fact, the requirement of the unit length is not necessary, a square centered at the origin is sufficient. Then BCpbecomes scale-invariant because the distribution of X′ and the specification region will be the same no matter which scale was used originally.
Suppose X is a bivariate normal random vector following Nm1; m2; s2 1; s22; r . Then, X′ follows N 1 USL1LSL1 m1 USL1þLSL1 2 ; 1 USL2LSL2 m2 USL2þLSL2 2 ; s21 USL1LSL1 ð Þ2; s2 2 USL2LSL2 ð Þ2; rÞ. 3.2. Estimation of BCp
In this subsection, wefirst develop an algorithm to calculate ^BCp, and then derive an approximate normal distribution for ^BCp by Taylor expansion. Based on this normal approximation, we will develop procedures for statistical inference about BCp, including the following: (i) testing whether the process is capable or not by hypothesis testing; (ii) constructing a CI of BCpto obtain the precision of the estimate; and (iii) providing a lower confidence bound for practical usage in quality assurance.
As mentioned earlier, in the univariate case, if we move the process mean m to the middle of the two specification limits, then Cp= Cpk, that is, Cp¼ maxmCpk. For the bivariate case, it is slightly more complicated because we need to center and scale the variables first and then find maxθBCpk. By setting the transformed specification region to be 12;12
1 2;12
, it is obvious that we should move the process mean to the origin. For computing efficiency, instead of rotating the process distribution for each θ to find/calculate BCp by the Monte Carlo integration (i.e., varying the transformed process distribution against the fixed transformed specification region), we keep the distribution fixed and rotate the specification region so that only one set of N-simulated transformed process data is needed for all θ to perform Monte Carlo integration. Figure 7(a) shows the relative position of the square and the process distribution for θ = 0∘ (or 90∘, 180∘, 270∘) whereas Figure 7(b) shows that for θ = 45∘ (or 135∘, 225∘, 315∘).
Because BCpis the maximum value of BCpk, it is natural to repeatedly use the algorithm for calculating BCpkto calculate ^BCp as follows.
Algorithm for calculating ^BCp: 1. Transform process data by
2 X specification process (a) 1 X 2 X process (b) 1 X specification
Figure 7. The two relative positions of the specification square w.r.t. the bivariate distribution when (a) θ = 0∘and (b)θ = 45∘(corresponding to BC p)
X1′ ¼ 1 USL1 LSL1 X1 USL1þ LSL1 2 and X2′ ¼ 1 USL2 LSL2 X2 USL2þ LSL2 2
so that the specification region becomes the unit square 1 2;12 1 2;12 . 2. Compute the sample covariance matrix ^Σ from the transformed data. 3. Calculate the eigenvalues ^l21and ^l22of ^Σ.
4. Loopθ over 0 ⩽ θ < 360∘tofind the optimal angle: for each candidate value of θ, • rotate the square 1
2;12 1 2;12 by an angleθ;
• use the same approach for Monte Carlo integration as that in calculating ^BCpto compute^q1ð Þ, ^qθ 2ð Þ, ^qθ 3ð Þ, and ^qθ 4ð Þ, theθ probabilities that a bivariate random vector following N 0 ; 0; ^l21; ^l22; 0 falls in the rotated specification region intersecting the four quadrants, respectively;
• then, compute ^q θð Þ ¼ ^q1ð Þ þ ^qθ 2ð Þ þ ^qθ 3ð Þ þ ^qθ 4ð Þ, the probability that the bivariate random vector falls in the rotatedθ square (see Figure 7 for (a)θ = 0∘and (b)θ = 45∘).
Find an angleθ such that ^q θð Þ is maximized over 0 ⩽ θ < 360∘. Denote the optimalθ by θ*and^q θð Þ by ^q . 5. Compute ^BCp¼ 13Φ
1 1^q
2
.
In an earlier study, we applied this algorithm to various bivariate normal processes and found that we always gotθ*= 45∘for the optimal^q. Then, we started wondering: is it true that BCpequals BCpkwhen and only when the process mean is at the center of the specification square and the two axes of the process distribution are exactly the two crossed lines connecting the vertices of the square as depicted in Figure 7(b)? The answer is yes. To show this,first, it is fairly obvious to see that only two particular positions of the square can have q1= q2= q3= q4as depicted in Figure 7; that is, whenθ is (i) 0∘(also, 90∘, 180∘, 270∘) as in Figure 7(a), or (ii) 45∘ (also, 135∘, 225∘, 315∘) as in Figure 7(b). This observation was further confirmed by computer computation. Next, we show that q (45∘)> q(0∘).
Because BCpkis a function of the yield q, wefirst derive the formula for q. As before, by considering the relative position of the transformed process and specification region, we can assume, without loss of generality, that the process follows N 0; 0; l 21; l22; 0, and the specification square for θ = 0∘has verticesp2ffiffiffi=2; 0, 0 ;pffiffiffi2=2,pffiffiffi2=2; 0, and 0 ; pffiffiffi2=2; see Figure 7(a). Denote the yields for the cases ofθ = 0∘andθ = 45∘as a functionl1andl2by q0(l1,l2) and q*(l1,l2), respectively. Then, givenl1andl2, for the case ofθ = 0∘, we have
q0ðl1; l2Þ ¼ 2Φ 1 2l1 1 2Φ 1 2l2 1 ; because X1and X2are independent. For the case ofθ = 45∘, we have
qðl1; l2Þ ¼ 4 R ffiffip2=2 0 Rx1þ ffiffi2 p =2 0 1 2pl1l2 e x 2 1 2l21 x2 2 2l22dx 2dx1 ¼ 4R ffiffip2=2 0 1 ffiffiffiffiffiffi 2p p l1 e x 2 1 2l21 Φ x1þ ffiffiffi 2 p =2 l2 1=2 dx1 ¼ 4R ffiffip2=2 0 f x1 l1 Φ x1þ ffiffiffi 2 p =2 l2 dx1 2Φ ffiffiffi 2 p 2l1 þ 1:
Because it is difficult to show that q0(l1,l2)< q*(l1,l2) for alll1andl2analytically, we verify the claim by computation as follows. We calculate these two values by numerical integration for various values ofl1andl2. Figure 8 presents the ratio q*/q0as a function of l1/l2 for l1=.25,.375,.5, 1, 2, 4, 8 and l1/l2= 1, 2,. . ., 100. The solid line presents q*/q0 for l1= 1, the three dashed lines are forl1=.25,.375,.5, and the three dotted lines are forl1= 2, 4, 8, respectively. It is clearly seen that the values of q*/q0are all greater than 1.
Therefore, instead of looping over various values ofθ to calculate the optimal ^q, we can simply set θ*= 45∘and calculate ^BCp directly. This would save tremendous amount of computing time. Thus, we can simplify the computing algorithm by replacing the original Step 4 with the following Step 4*:
4*. Generate N bivariate normal variates from N 0 ; 0; ^l21; ^l22; 0 . Compute the proportion of the N bivariate data that fall in the unit square as depicted in Figure 7(b) to obtain^q.
3.3. Statistical inference about BCpvia normal approximation
The exact distribution of ^BCpis mathematically intractable. Fortunately, we can obtain a normal approximation to the distribution of ^
BCpby taking itsfirst-order Taylor expansion as follows.
The partial derivatives of q*(l1,l2) with respective tol1andl2are, respectively, @qðl1; l2Þ @l1 ¼ 4 Z ffiffip2=2 0 f x1 l1 x2 1 l2 1 l1 1 ! Φ x1þ ffiffiffi 2 p =2 l2 dx1þ f ffiffiffi 2 p 2l1 ffiffiffi 2 p l2 1 ! ; @qðl1; l2Þ @l2 ¼ 4 Z ffiffip2=2 0 f x1 l1 f x1þ ffiffiffi 2 p =2 l2 x1 ffiffiffi 2 p =2 l2 2 ! dx1;
which can be evaluated numerically when givenl1andl2.
Denote Q1ðl1; l2Þ @q@lðl11;l2Þand Q2ðl1; l2Þ @q@lðl12;l2Þ. Then, an approximate distribution of ^BCp can be obtained as
N BCp; Q2 1ðl1; l2Þl21þ Q22ðl1; l2Þl22 2 36n f 3BCp 2 ! (15)
by Taylor expansion. The derivation of (15) is given in the Appendix.
Plugging ^liforli, i = 1, 2, we can construct from the approximate normal distribution of ^BCpan approximate 100(1-a)% CI as
^ BCp Za=2 Q2 1 ^l1; ^l2 ^l2 1þ Q22 ^l1; ^l2 ^l2 2 1=2 ffiffiffiffiffiffiffiffi 72n p f 3 ^BCp (16)
and an approximate 100(1a)% lower confidence bound as
^ BCp Za Q2 1 ^l1; ^l2 ^l2 1þ Q 2 2 ^l1; ^l2 ^l2 2 1=2 ffiffiffiffiffiffiffiffi 72n p f 3 ^BCp ; (17)
where Zais theath upper quantile of the standard normal distribution.
To investigate whether the process capability meets customers’ demands or not, practitioners can perform hypothesis testing on BCp. Consider the hypotheses
H0: BCp⩽C process is not capableð Þ H1: BCp> C process is capableð Þ
where C> 0 is the preset acceptable capability level. A naive test can be conducted with the following test statistic:
0 10 20 30 40 50 60 70 80 90 100 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 λ1=1 λ 1=2 λ1=4 λ1=8 λ1=0.5 λ1=0.375 λ1=0.25 λ1/λ 2
Figure 8. Ratio q*(l1,l2)/q0(l1,l2) forl1=.25,.375,.5, 1, 2, 4, 8 andl1/l2= 1, 2,. . ., 100
Z¼ ffiffiffiffiffiffiffiffi 72n p f 3 ^BCp ^ BCp C Q2 1 ^l1; ^l2 ^l2 1þ Q22 ^l1; ^l2 ^l2 2 1=2
and the decision making rule at the significance level a is
reject H0; if Z> Za; do not reject H0; if Z⩽ Za:
3.4. The simulated examples
We apply the proposed methods for BCpto the four examples in Table II. To evaluate how well the estimation method performs, we generate 1000 sets of 100 bivariate normal data to obtain 1000 ^BCps for each case. Table X presents the sample mean and sample standard deviation of 1000 ^BCps as well as the true value of BCpfor each of Cases 1–4. With the sample size only 100, the estimation result is quite satisfactory. Table XI gives the 90% approximate LCB and 90% approximate CI of one simulated data set of size n = 100 for each case. The results indicate that the estimation method is satisfactory.
3.5. Multivariate Cpindex: MCp
The BCpcan be easily extended to the general case of k quality characteristics. After rescaling the data and the specifications as described earlier, the definition (13) becomes
MCp max m;R MCpk;
where R is any rotation matrix in the Euclidean space Rk. Analogous to (14), MCpsatisfies MCp¼ 1 3Φ 1 2k1 p 2k ¼ 1 3Φ 1 p 2 ¼ 1 3Φ 1 1 q 2 ;
where q is the probability of a random normal vector with mean 0 and covariance matrixΣ = diag(l1,l2,. . ., lk) falling in the origin-centered unit cube with the 2kvertices all at axes, analogous to Figure 7(b).
4.
Three real-life application examples
For demonstration, we employ our estimating methods given in the last two sections to three real-life industrial examples. Thefirst two examples were described in Chen12involving bivariate processes, and the third one involves a trivariate process described in Pan et al.28
The first example was originally presented by Sultan29 regarding an industrial process in which the Brinell hardness (X1) and the tensile strength (X2) are the quality characteristics. The Chen–Sultan’s data consist of 25 samples taken from a process with the specifications for hardness and tensile strength being [112.7, 241.3] and [32.7, 73.3], respectively. Figure 9 depicts the data and the
Table X. Summary of 1000 ^BCps
Case True value Sample mean Sample SD Bias
1 1.245320 1.2453743 0.0692768 0.000054
2 1.404933 1.404937 0.0882426 0.000004
3 0.927604 0.9277373 0.0532411 0.000133
4 0.326505 0.3280435 0.0230090 0.001539
Table XI. 90% approximate lower confidence bound and confidence interval of BCpfrom one data set
Case LCB CI 1 1.183757 [1.166305, 1.324336] 2 1.350439 [1.334990, 1.474876] 3 0.870863 [0.854778, 1.000430] 4 0.298699 [0.290817, 0.362191]
502
specification region. Table XII presents the estimate ^BCpkalong with the corresponding 90% CIs and 90% lower confidence bounds obtained by employing the four bootstrap methods with B = 3000 for thefirst example. Based on the ^BCpkvalue of 1.050281, we can obtain by inequality (7) the estimated upper bound of the nonconforming rate to be 1628 ppm. Also, taking the LCB of the percentile bootstrap method as an example, with LCB = 0.7977719, by inequality (8), we can say that, with 90% confidence, the product yield is at least 98.3303%.
The second example was originally presented by Pal13regarding a manufacturing process of a bobbin with the height (X
1) and the weight (X2) as quality characteristics. In this example, one hundred samples were taken from a process with the specification [40, 42] for height and [44, 46.5] for weight. The data and the specification region are displayed in Figure 10. The results of this example are presented in Table XIII. For this example, the ^BCpk¼ 0:9610820, which leads to an estimated upper bound of the nonconforming rate of 3936 ppm. Also, with the LCB of the percentile bootstrap method being 0.8825650, one can conclude that the yield of the product is at least 99.1896% with 90% confidence.
As to BCp, by our efficient algorithm given in Subsection 3.2 (with Step 4*), for Chen–Sultan’s example, we obtain ^BCp¼ 1:1228071 with the optimal^p ¼ 0:0001890. Based on the normal approximation, (16) gives a 90% CI (0.8577761, 1.3878381) and (17) gives a 90% lower confidence bound (0.9163140). Similarly, for Pal’s data, we obtain ^BCp¼ 1:1890221 with the optimal ^p ¼ 0:0000903, an approx-imate 90% CI (1.0748297, 1.3032146) and an approxapprox-imate 90% lower confidence bound (1.1000516).
The third example is related to the example presented in Pan and Lee11regarding the solder paste stencil printing process, a cost-effective process that has been widely used in traditional high-volume surface mount assembly. For more descriptions about the pro-cess, see Pan et al.28In this process, solder-deposited volume (X
1), area (X2), and height (X3) are the three quality characteristics with high correlation. There are two kinds of stencils (of the same pattern) with thickness, 0.1 mm (4mil) and 0.15 mm (6mil) (1mil = 0.0254
mm). And there arefive different aperture sizes, 30, 25, 20, 16, and 12 (mil). Pan and Lee11considered a QFP
4mil, 30process to demon-strate the effectiveness of their new multivariate PCIs, where QFP4mil, 30represents that the process is for quadflat package (QFP), one of the advanced packages, with the stencil thickness as 4mil and aperture size as 30. The specifications and target values for this QFP4mil, 30process are given in Table XIV. The authors kindly provided us with 150 measurements (not the same set of data as in Pan and Lee11) from the same process. Figure 11 displays the 150 data points and the specification region. We first check if the process follows a multivariate normal distribution by applying the Shapiro–Wilk normality test to these measurements, which
100 120 140 160 180 200 220 240 260 20 30 40 50 60 70 80 USL2=73.3 LSL2=32.7 LSL1=112.7 USL1=241.3 X1 (Hardness) X2 (Strength)
Figure 9. Data and specification region of Chen–Sultan’s example
Table XII. BCpkestimate, 90% confidence intervals, and 90% bootstrap lower confidence bounds of Chen–Sultan’s example ^
BCpk¼ 1:050281
^p1 ^p2 ^p3 ^p4
0.000312 0.000043 0.000407 0.00005
90% confidence interval
Basic Standard Percentile Bias-corrected percentile
(0.7543558, 1.3818883) (0.7196855, 1.3147242) (0.7186737,1.3462062) (0.8036442,1.4171286)
90% bootstrap lower confidence bound
Basic Standard Percentile Bias-corrected percentile
0.8451777 0.7853992 0.7977719 0.8544017
gives a p-value of 0.8489; so that the normality holds. The sample mean vector and the sample covariance matrix are, respectively, XT ¼ 0:075859; 0:817971; 0:097080ð Þ and S¼ 0:0000250 0:0002601 0:0000012 0:0002601 0:0028808 0:0000079 0:0000012 0:0000079 0:0000151 0 @ 1 A:
Applying the proposed procedures to this dataset, we obtain the estimation results (presented in Table XV) for MCpk, including the estimate ^MCpkand the 90% CI/LCB obtained via bootstrap with B = 3000. With ^MCpk¼ 0:9355062, the estimated upper bounds of the nonconforming rates is 5008 ppm by inequality (10). Again, taking the LCB of the percentile bootstrap method as an example, with LCB = 0.8620695, by thefirst inequality in (11), we can say that, with 90% confidence, the yield of the product is at least 99.4992%. As to MCp, by extending our efficient algorithm given in Subsection 3.2 to three dimensions (see Subsection 3.5 for the extension of Step 4*), we obtain ^MCp¼ 1:1615466 with the optimal ^p ¼ 0:0000616.
For practitioners’ convenience, MATLAB (MathWorks, Natick, MA, USA) programs for computing the CIs and lower confidence bounds using the bootstrap approach are provided in http://www.stat.nctu.edu.tw/~jjhs/MPCI.zip or upon request.
40 40.5 41 41.5 42 43.5 44 44.5 45 45.5 46 46.5 47 USL2=46.5 LSL2=44 USL1=42 LSL1=40 X1 (Height) X2 (Weight)
Figure 10. Data and specification region of Pal’s example
Table XIII. BCpkestimate, 90% confidence intervals, and 90% bootstrap lower confidence bounds of Pal’s example ^
BCpk¼ 0:9610820
^p1 ^p2 ^p3 ^p4
0.000984 0.000083 0.000011 0.000081
90% confidence interval
Basic Standard Percentile Bias-corrected percentile
(0.8523437,1.0597296) (0.8581923,1.0663559) (0.8624344,1.0698203 ) (0.8641023,1.0719583)
90% bootstrap lower confidence bound
Basic Standard Percentile Bias-corrected percentile
0.8778255 0.8811811 0.8825650 0.8842324
Table XIV. The target values and specifications for QFP4mil, 30process
Quality characteristic Target USL LSL
Deposited volume (X1) 0.0787 0.10250 0.0549
Deposited area (X2) 0.7870 0.96870 0.6052
Deposited height (X3) 0.1000 0.12765 0.07235
5. Conclusion
In this paper, we studied the bivariate PCIs, BCpkand BCp, proposed by Castagliola and Castellanos 22
and extended their notion to processes of more than two quality characteristics.
We established a link between BCpkand the process yield by showing that the same inequality as in the univariate case holds, that is, 2Φ(3BCpk) 1 ⩽ %yield. This lower bound provides a measure for quality assurance. With the same notion of BCpk, we defined an index MCpkfor processes of more than two characteristics and proved that the lower bound inequality 2Φ(3MCpk) 1 ⩽ %yield also holds. We also provided a new algorithm for computing the natural estimate ^MCpkof MCpk, which is more efficient than the algorithm provided for the bivariate case by Castagliola and Castellanos.22Moreover, the new algorithm can be used for processes with more general specification regions. For statistical inference, we utilized the bootstrap approach to obtain a lower confidence bound of BCpk. Among the four popular bootstrap methods under study, we recommend the percentile bootstrap method.
For BCp, we found that the original definition given in Castagliola and Castellanos22is not scale-invariant. We proposed a simple preprocessing step tofix the problem. By finding the exact situation when BCpk= BCp(also for MCpk= MCp), we developed an efficient algorithm for computing the natural estimate of BCp, which is a lot faster than the method given in Castagliola and Castellanos.22We further derived an approximate normal distribution for ^BCp by taking itsfirst-order Taylor expansion. This enabled us to derive statistical procedures for making inferences about process capability based on data, including hypothesis testing, confidence interval, and lower confidence bound.
Our simulation studies indicated that the sampling distribution of ^BCp is fairly close to a normal distribution. If one couldfind a suitable normal approximation for it, then statistical inferences about BCpk based on the normal approximation would be more computationally efficient and perhaps more statistically efficient than that obtained by the bootstrap approach.
Acknowledgements
The authors would like to express their gratitude to the editor and an anonymous referee for the careful review and constructive suggestions. The work was supported in part by the National Research Council of Taiwan, Grant No. NSC97-2118-M-009-002-MY2 and NSC99-2118-M-009-003-MY2. 0.050.0549 0.07 0.08 0.09 0.1025 0.11 0.55 0.6052 0.7 0.75 0.8 0.85 0.9 0.9687 1.05 0.06 0.07235 0.1 0.12765 USL 1 X1 (Deposited Volume) LSL1 LSL2 X2 (Deposited Area) USL2 LSL3 USL3 X3 (Deposited Height)
Figure 11. Data and specification region of stencil printing example
Table XV. MCpkestimate, 90% confidence intervals, and 90% bootstrap lower confidence bounds of stencil printing example ^
MCpk¼ 0:9355062
^p1 ^p2 ^p3 ^p4 ^p5 ^p6 ^p7 ^p8
0.000005 0.000011 0.000015 0.000014 0.000597 0.000602 0.000611 0.000626
90% confidence interval
Basic Standard Percentile Bias-Corrected Percentile
(0.8241617,1.0277553) (0.8362759,1.0406581) (0.8432572,1.0468508) (0.8447718,1.0490951)
90% bootstrap lower confidence bound
Basic Standard Percentile Bias-corrected percentile
0.8513668 0.8588471 0.8620695 0.8636374
References
1. Kane VE. Process capability indices. Journal of Quality Technology 1986; 18(1):41–45.
2. Chan LK, Cheng SW, Spring FA. A new measure of process capability: Cpm. Journal of Quality Technology 1988; 20(3):162–175.
3. Pearn WL, Kotz S, Johnson NL. Distributional and inferential properties of process capability indices. Journal of Quality Technology 1992; 24(4):216–231.
4. Kotz S, Johnson NL. Process Capability Indices. Chapman and Hall: London, 1993.
5. Kotz S, Lovelace CR. Process Capability Indices in Theory and Practice. Arnold: London, 1998.
6. Pearn WL, Kotz S. Encyclopedia and Handbook of Process Capability Indices: A Comprehensive Exposition of Quality Control Measures. World Scientific, Hackensack: New Jersey, 2006.
7. Chan LK, Cheng SW, Spring FA. Multivariate measure of process capability. Journal of Modeling and Simulation 1991; 11(1):1–6.
8. Hubele HF, Shahriari H, Cheng CS. A bivariate process capability vector. In Statistical Process Control in Manufacturing, Keats JB, Montgomery DC
(eds). Marcel Dekker: New York, 1991; 299–310.
9. Shahriari H, Hubele NF, Lawrence FP. A multivariate process capability vector. Proceedings of the 4th Industrial Engineering Research Conference, Institute of Industrial Engineers 1995; 304–309.
10. Taam W, Subbaiah P, Liddy JW. A note on multivariate capability indices. Journal of Applied Statistics 1993; 20(3):339–351.
11. Pan JN, Lee CY. New capability indices for evaluating the performance of multivariate manufacturing processes. Quality and Reliability Engineering International 2010; 26(1):3–15.
12. Chen H. A multivariate process capability index over a rectangular solid tolerance zone. Statistica Sinica 1994; 4(2):740–758. 13. Pal S. Performance evaluation of a bivariate normal process. Quality Engineering 1999; 11(3):379–386.
14. Bothe DR. Composite capability index for multiple product characteristics. Quality Engineering 1999; 12(2):253–258.
15. Wang FK, Hubele NF, Lawrence FP, Miskulin JD, Shahriari H. Comparison of three multivariate process capability indices. Journal of Quality
Technology 2000; 32(3):263–275.
16. Wang FK, Chen JC. Capability index using principal component analysis. Quality Engineering 1998–1999; 11(1):21–27.
17. Wang FK, Du T. Using principal component analysis in process performance for multivariate data. Omega: The International Journal of Management Science 2000; 28(2):185–194.
18. Shinde RL, Khadse KG. Multivariate process capability using principal componet analysis. Quality and Reliability Engineering International 2009; 25(1):69–77.
19. Shinde RL, Khadse KG. A review and comparison of some multivariate process capability indices based on fraction conforming interpretation. Statistical Methods 2005; 7:95–115.
20. Gonzalez I, Sanchez I. Capability indices and nonconforming proportion in univariate and multivariate processes. International Journal of Advanced
Manufacturing Technology 2009; 44(9):1036–1050.
21. Boyles RA. The Taguchi capability index. Journal of Quality Technology 1991; 23(1):17–26.
22. Castagliola P, Castellanos JG. Capability indices dedicated to the two quality characteristics case. Quality Technology & Quantitative Management 2005; 2(2):201–220.
23. Golub GH, Van Loan CF. Matrix Computations (2nd edn). The Johns Hopkins University Press, Baltimore: Maryland, 1989. 24. Efron B. Bootstrap methods: another look at the jackknife. The Annals of Statistics 1979; 7(1):1–26.
25. Davison AC, Hinkly DV. Bootstrap Methods and their Application. Cambridge University Press: United States of America, 1997.
26. Efron B, Tibshirani RJ. Bootstrap methods for standard errors, confidence interval, and other measures of statistical accuracy. Statistical Science 1986; 1(1):54–77.
27. Carpenter J, Bithell J. Bootstrap confidence intervals: when, which, what? A practical guide for medical statisticians. Statistics in Medicine 2000; 19(9):1141–1164.
28. Pan J, Tonkay GL, Storer RH, Sallade RM, Leandri DJ. Critical variables of solder paste stencil printing for Micro-BGA andfine-pitch QFP. IEEE
Transactions on Electronics Packaging Manufacturing 2004; 27:125–132.
29. Sultan TI. An acceptance chart for raw materials of two correlated properties. Quality Assurance 1986; 12(3):70–72. 30. Anderson TW. An Introduction to Multivariate Statistical Analysis (3rd edn). Wiley and Sons, Hoboken: New Jersey, 2003.
Appendix A
Derivation of the normal approximation given in (15)
Let f qð Þ ¼ BCp¼ 13Φ 1 1q 2 . Denote^q q ^l1; ^l2
. Then, fð Þ ¼ ^^q BCp. Expanding fð Þ at q by Taylor expansion, we have^q ^ BCp 1 3Φ 1 1 q 2 þ 1 6f Φ1 1 q 2 ^q qð Þ ¼ BCpþ 1 6f 3BCp ^q qð Þ: (A:1)
Because ^lis are eigenvalues of the sample covariance matrix ^Σ, by Anderson30(pages 473–474), we have ffiffiffi n p ^l1 ^l2 l1 l2 !d N 0; 1 2l 2 1 0 0 1 2l 2 2 2 6 4 3 7 5 0 B @ 1 C A as n!1 Then, by Theorem 4.2.3 of Anderson,30
ffiffiffi n p q ^l1; ^l2 q lð 1; l2Þ !d N 0;1 2 Q 2 1ðl1; l2Þl21þ Q 2 2ðl1; l2Þl22 as n! 1. Now denoting Zq pffiffiffin q ^l1; ^l2 q lð 1; l2Þ , we have, by (A.1),
506
^ BCp BCpþ 1 6pffiffiffinf 3BCp Zq and (15) is derived. Authors' biographies
Jyh-Jen Horng Shiau is currently a full professor in the Institute of Statistics at National Chiao Tung University, Taiwan. She holds a BS in Mathematics from National Taiwan University, Taipei, Taiwan, an MS in Applied Mathematics from the University of Maryland Baltimore County, an MS in Computer Science and PhD in Statistics from the University of Wisconsin-Madison. Formerly, she taught at Southern Methodist University, the University of Missouri at Columbia, and National Tsing Hua University, and worked at AT&T Bell Laboratories before she moved to Taiwan. She is a former managing editor of an international journal Quality Technology & Quantitative Management (2004-2006). Her primary research interests include industrial statistics, nonparametric and semiparametric regression, and functional data analysis. She is a lifetime member of the International Chinese Statistical Association.
Chia-Ling Yen received her PhD degree in Statistics from National Chao-Tung University in 2008. Her research interests include multivariate analysis and statistical process control.
W. L. Pearn received the PhD degree in Operations Research from the University of Maryland, College Park. He is a professor of Operations Research and Quality Assurance at National Chiao-Tung University (NCTU), Hsinchu, Taiwan, China. He was with AT&T Bell Laboratories as a member of the quality research staff before joining NCTU. His research interests include process capability, network optimization, and production management.
Wan-Tse Lee received her master’s degree in Statistics from National Chao-Tung University in 2009. She is currently working as a quality engineer in AU Optronics Corporation (AUO) in the Hsinchu Science and Industrial Park, Hsinchu, Taiwan. Her research interests include process capability and statistical process control.
