The average run length (ARL), representing the average number of sam-ples taken before an action signal is given, is the most popular technique in evaluating a control chart or comparison of alternative control charts. In this section, we make further comparison of ARL0s for control charts constructed by empirical quantile and symmetric quantile.
Let α = (α1, ..., αk)0 be the percentage vector considered for constructing a quantile control chart. Let Qn(α) be the quantile control chart for monitoring the distributional shift of random variable X. We assume that when the pro-cess is in control, the quantile control chart satisfies the following asymptotic property
n(Qn(α) − Q0(α))0Σ−10 (Qn(α) − Q0(α)) → χ2(k) in distribution.
Hence we actually consider the following hypothesis
H0 : Q(α) = Q0(α) vs Q(α) 6= Q0(α), (3.15)
where Q(α) is the true population quantile vector estimated by statistic Qn(α).
Suppose that the significance level for the control chart is α. Then, the quantile control chart indicates rejecting H0 when
n(Qn(α) − Q0(α))0Σ−10 (Qn(α) − Q0(α) ≥ χ2α(k).
We note here that control chart by symmetric quantiles and by empirical quan-tiles have the same vector Q0(α) but have different asymptotic covariance ma-trices Σ0.
To study the performance of the control charts in terms of ARL, we consider linear function aX + b for representation of distributional shift. Note that the population quantile vector and covariance matrix for aX + b are, respectively, QaX+b(α) = aQ0(α) + b and ΣaX+b = a2Σ0. Then the probability that we claim for an out of control under the linear function aX + b is
PQaX+b(α),ΣaX+b(n(Qn(α) − Q0(α))0Σ−10 (Qn(α) − Q0(α)) ≥ χ2α(k))
= PQaX+b(α),ΣaX+b(n(Qn(α) − QaX+b(α) + (−(1 − a)Q0(α) + b))0Σ−1aX+b (Qn(α) − QaX+b(α) + (−(1 − a)Q0(α) + b)) ≥ 1
a2χ2α(k))
= P (χ2(k, n
a2(−(1 − a)Q0(α) + b)0Σ−10 (−(1 − a)Q0(α) + b)) ≥ 1
a2χ2α(k)) (3.16) (3.17) since√
nΣ−1/2aX+b0(Qn(α)−QaX+b(α)+(−(1−a)Q0(α)+b)) app∼ Nk(√
nΣ−1/2aX+b0(−(1−
a)Q0(α) + b), Ik) = Nk(
√n
a Σ−1/20 0(−(1 − a)Q0(α) + b), Ik) where χ2(k, λ) has a
noncentral chi-square distribution with non-centrality parameter λ. The ARL for quantile vector Qn(α) is
ARL = 1
P (χ2k(λ = an2(−(1 − a)Q0(α) + b)0Σ−10 (−(1 − a)Q0(α) + b)) ≥ a12χ2α(k)) where a level α control chart is expected to have ARL = α1 when the process is in control.
Since the asymptotic covariance matrices ΣaX+bfor control charts constructed by symmetric quantiles and by empirical quantiles are different that leads to different non-centrality parameters resulting varied performances in ARL0s. We then compute the ARL0s for the symmetric quantile control chart and the em-pirical quantile control chart for comparison.
Let us fix significance level at α = 0.005 and n = 1000, for a non-parametric study, we consider the Laplace distribution Lap(2) as the in control distribution with percentage vectors given as
k = 4 : (α1, α2, ..., α4) = (0.02, 0.05, 0.95, 0.98),
k = 10 : (α1, α2, ..., α10) = (0.02, 0.05, 0.13, 0.25, 0.37, 0.63, 0.75, 0.87, 0.95, 0.98)
For easiness of expression, we denote the ARL0s for empirical quantile charts and symmetric quantile charts, respectively, by ARLe and ARLs. The com-puted ARL0s are displayed in Table 3.2 and Table 3.3.
Table 3.2. Comparison of symmetric and empirical quantile charts by ARL (k=4)
(a, b) ARLs ARLe (a, b) ARLs ARLe
(1, 0) 200 200 (1.2, 0.5) 17.00 21.73 (1, 0.2) 172.14 196.81 (1.2, 1) 8.32 18.54 (1, 0.5) 90.89 181.30 (1.2, 2) 2.20 11.91 (1, 1) 22.74 138.83 (1.5, 1) 3.02 3.91 (1, 2) 2.68 61.77 (2, 1) 1.44 1.54 (1, 5) 1.00 5.32 (2, 2) 1.25 1.47
Table 3.3. Comparison of symmetric and empirical quantile charts by ARL (k=10)
(a, b) ARLs ARLe (a, b) ARLs ARLe
(1, 0) 200 200 (1.2, 0.5) 2.27 4.15 (1, 0.2) 102.4 140.7 (1.2, 1) 1.08 1.64 (1, 0.5) 13.05 37.01 (1.2, 2) 1 1 (1, 1) 1.42 3.97 (1.5, 0.5) 1.09 1.29
(1, 2) 1.00 1.01 (2, 1) 1 1
(1, 5) 1 1 (2, 2) 1 1
The case (a, b) = (1, 0) represents the process being in-control and both ARL0s are the expected number 200 for setting α = 0.05. Surprisingly ARL0ss are all smaller than the corresponding ARL0es unless they are number 10s. This
indicates that the symmetric quantiles based control chart can detect the dis-tributional shift with smaller number of samples.
We see that in this setting of coverage interval the symmetric quantiles based control chart is still more efficient than the empirical quantiles based control chart in detection of distributional shift.
3.6 Concluding Remarks
In contrast with empirical quantile based control chart of Grimshaw and Alt (1997), a symmetric quantile based control chart is proposed in this dissertation to monitor the popoulation-quantile vector aiming for monitoring more detailed features of a population distribution. The asymptotic theorem is also derived.
The symmetric quantile based control chart totally dominates the empirical one across all αi of a population quantile vector Q(α1, ..., αk), when the un-derlying distribution is Laplace distribution, which is widely used in modelling spectral vector of speech signals in speech recognition.
Chapter 4
Monitoring Nonlinear Profiles with Random Effects by Nonparametric Regression
In this chapter, we study nonlinear profile monitoring schemes. Principal component analysis is conducted, and a T2 chart and a combined chart based on principal component scores are studied as well as individual Principal Com-ponent charts.
4.1 Proposed Monitoring Schemes
4.1.1 A Motivated Example
This study was motivated by the aspartame example given in Kang and Albin (2000). Since no data are available, a profile of the form Y = I + M eN (x−1)2 + is used to mimic an aspartame profile. Then the idea is to perturb the parameters I, M, N randomly to create allowable profile-to-profile variations for an in-control process.
Thus the following random-effect model was considered to generate aspar-tame profiles:
Yj = I + M eN (xj−1)2 + j, j = 1, · · · , p, (4.1) where I ∼ N (µI, σI2), M ∼ N (µM, σM2 ), N ∼ N (µN, σ2N), ∼ N (0, σ2), and all the random components are independent of each other. Unfortunately, the response profile YYY = (Y1, · · · , Yp)0 of model (4.1) has a complicated distribution with mean µµµ = (µ1, · · · , µp)0 and covariance matrix ΣΣΣ as follows. For i, j = 1, · · · , p,
µj = E(Yj) = µI + µMeµN(xj−1)2+
σ2N (xj −1)4
2 ,
Cov(Yi, Yj) = σI2+ (µ2M + σM2 )
eµN[(xi−1)2+(xj−1)2]+
σ2N [(xi−1)2+(xj −1)2]2 2
−µ2MeµN(xi−1)2+
σ2N (xi−1)4
2 +µN(xj−1)2+σ2N (xj −1)
4
2 + σ2εδij, (4.2) where δij = 1 if i = j; and 0 otherwise. Note that, by (4.2), the covariance ma-trix ΣΣΣ will be changed if the mean of M or N shifts, a situation too complicated
to analyze the performance of the control charts under study.
So, instead, we model the aspartame profiles as realizations of a Gaussian stochastic process with the mean function
µ(x) = µI + µMeµN(x−1)2 (4.3) and a covariance function G(s, t), where s, t are in the domain of x. To re-tain a similar profile-to-profile variation as it would be in the random-effect model (4.1), we let the in-control profiles follow M V N (µµµ0, ΣΣΣ), where µµµ0 = (µ01, · · · , µ0p)0 with
µ0j = µI + µMeµN(xj−1)2, j = 1, · · · , p, (4.4) and ΣΣΣ is the covariance matrix given by equation (4.2).
When the mean function (4.3) is shifted, say, µI to µI+ασI, µM to µM+βσM, and µN to µN + γσN, µ0j is shifted from µI + µMeµN(xj−1)2 to
˜
µj ≡ (µI + ασI) + (µM + βσM)e(µN+γσN)(xj−1)2, j = 1, · · · , p.
Let µµeµ = (eµ1, · · · ,µep)0. Then the shift on the mean of YYY is δδδ ≡µeµµ − µµµ0.
4.1.2 Data Smoothing
In order to extend nonlinear profiles of a fixed parametric form to smooth profiles of a flexible nonparametric form, a smoothing technique is needed for de-noising sample profiles. The idea of smoothing is to fit a smooth function whose final form is determined by the data and the chosen level of smoothness for the
curve. One popular approach is to fit noisy data by splines. Frequently, cubic splines (i.e., piecewise cubic polynomials with continuous second derivatives) are used for such approximations. Two commonly used spline smoothing techniques are smoothing splines and B-spline regression, both are available in popular statistical packages like R, Splus, and others. Other smoothing techniques such as local polynomial smoothing and wavelets can be used as well. We remark based on our experiences that, by filtering out noises, the actual signals can be better extracted from the data and PCA can explore the variation among profiles a lot better. In particular, smoothing tends to be more advantageous as the noise level (σ2) gets larger.
4.1.3 Phase I Monitoring
Assume that a set of n historical profiles is available for Phase I analysis. We first apply a smoothing technique to each of the n profiles to filter out the noise, and then apply PCA to the smoothed profiles as follows. Denote the (p × 1) data vector of the i-th profile by yi and the usual sample covariance matrix of {yi, i = 1, · · · , n} by SSS. Apply the eigenanalysis to SSS. The eigenvector vr
corresponding to the r-th largest eigenvalue λr is the r-th principal component and Sir ≡ v0ryi is called the score of the r-th principal component of the i-th profile, r = 1, · · · , p, i = 1, · · · , n.
We select the number of the “effective” principal components by considering the total variation explained by the chosen principal components along with
the principle of parsimoniousness that we often use in the variable selection problem. Denote this number by K and the (K × 1) score vector (Si1, · · · , SiK)0 by sssi.
For Phase I monitoring, due to the dependency of the K PC-scores, we adopt the usual Hotelling T2 statistic described below. For the i-th profile, i = 1, . . . , n, the T2 statistic is defined as
Ti2 = (sssi− ¯sss)0BBB−1(sssi− ¯sss), (4.5) where ¯sss = Pn
i=1sssi/n and BBB = Pn
i=1(sssi− ¯sss)(sssi− ¯sss)0/(n − 1), the usual sample mean and sample covariance matrix of the score vectors.
Since score vectors are distributed as multivariate normal asymptotically (Anderson, 2003), according to Tracy et al. (1992), also Sullivan and Woodall (1996), we have
n
(n − 1)2Ti2 ∼ Beta K
2 ,n − K − 1 2
approximately.
Thus, an approximate α-level upper control limit can be set at the 100(1 − α) percentile of the beta distribution with K/2 and (n − K − 1)/2 as parameters.
For Phase I analysis, perform control-charting with the T2 statistic of the score vectors in (4.5) to detect the out-of-control profiles in the historical data set. If there are any, remove them and redo PCA and control-charting with the remaining profiles. Repeat this procedure until all the remaining profiles are within the control limit. These remaining profiles are considered as “in-control”
profiles and can be used to characterize the in-control process. The resulting
principal components and eigenvalues can then be used to set up the control limit for Phase II on-line monitoring.
4.1.4 Phase II Monitoring
As in most of Phase II studies, we assume the in-control process distribution of the profiles after de-noising has been characterized as Np(µµµ0, ΣΣΣ0), either from prior experiences or estimated from the Phase I analysis.
Our Phase II monitoring schemes are also based on PCA. Apply PCA to ΣΣΣ0 to obtain eigenvalues, λ1 ≥ · · · ≥ λp ≥ 0, and the corresponding eigenvectors, v1, · · · , vp. Similar to that in Phase I analysis, choose the number of effective principal components K based on the parsimoniousness and the total variation that the first K PCs account for. More specifically, since the r-th PC accounts for λr/Pp
r=1λr of the total variation, we can simply choose the first K such that PK
r=1λr/Pp
r=1λr reaches a desired level.
Now for each of the incoming profiles in Phase II monitoring, first smooth and then project it onto the first K PCs to obtain K PC-scores. Denote these scores by S1, · · · , SK. Since these scores are independent and Sr follows a normal distribution with mean v0rµ0 and variance λr when the process is in control, it is easy to construct a control chart for each of the K PC-scores accordingly. Denote the desired in-control false-alarm rate by α. Then the control limits for the r-th PC-score chart, which monitors the statistic Sr, is v0rµ0± Zα/2√
λr, r = 1, · · · , K.
If a particular mode of process change can be caught by one of the first K principal components, then we can use that particular PC-score chart to monitor it. However, very often a process shift is reflected in more than one principal component. When this happens, we can consider a combined chart scheme by combining all K PC-score charts. A combined chart scheme signals out-of-control when any of the K individual charts signals. Thus, the proposed combined chart is equivalent to monitoring the statistic
1≤r≤Kmax |Sr− v0rµ0
√λr
|.
This chart signals out-of-control when max1≤r≤K|(Sr − v0rµ0)/√
λr| > Zα0/2, where the individual false-alarm rate α0 should be chosen at the level of 1 − (1 − α)1/K so that the overall false-alarm rate is at the desired level α.
We can also consider a T2 chart by monitoring the statistic T2 =
K
X
r=1
(Sr − v0rµ0)2
λr , (4.6)
which follows the chi-square distribution with K degrees of freedom (denoted by χ2K) when the process is in control. Thus, the upper control limit is the 100(1 − α) percentile of χ2K.
4.1.5 ARL of the Proposed Schemes
We evaluate the performances of the proposed Phase II monitoring schemes described above in terms of ARL, the average run length. The ARL values of the individual PC-score chart can be computed as follows. Assume that the
mean of the profile has been shifted from µ0 to µ0 + δ. The probability of detecting the shift by the r-th PC-score chart is
p = 1 − P (|Sr − v0rµ0
where Φ is the cumulative distribution function of the standard normal variate Z and Zα/2 is the 100(1 − α/2) percentile of Z. Then the value 1/p is the ARL of the r-th PC-score chart.
Since the PC-scores S1, · · · , SK are independent, the ARL of the combined chart also can be computed easily by the reciprocal of
p = 1 − P ( max
Since T2 statistic in (4.6) follows a noncentral chi-square distribution with K degrees of freedom and non-centrality ξ = PK
r=1(v0rδ)2/λr (denoted by χ2K(ξ)).
Then the detecting power of the T2 chart can be easily calculated by p = P (T2 > χ2K,α) = P (χ2K(ξ) > χ2K,α),
where χ2K,α denotes the 100(1 − α) percentile of the central chi-square distribu-tion χ2K.