包含區間, 多變量管制圖及構面性資料統計品質管制

(1)

國立交通大學

統計學研究所

博

博士論文

士論文

包含區間

包含區間,

,

, 多變量管制圖

多變量管制圖

及

構面性資料統計品質管制

Coverage Interval, Multivariate Control Chart

and Profile Monitoring

研

究

生

生 :

生

:

: 林碩慧

林碩慧

指導教授

指導教授 :

:

: 陳鄰安

陳鄰安

教授

指導教授

指導教授 :

:

: 洪志真

洪志真

教授

中華民國一百

中華民國一百零一

零一

零一年

年

年一

一

一月

月

(2)

包含區間

包含區間,

,

, 多變量管制圖

多變量管制圖

及

構面性資料統計品質管制

Coverage Interval, Multivariate Control Chart

and Profile Monitoring

研研研

研究究究究生生生生：：：：林碩慧林碩慧林碩慧林碩慧 Student：：：Shuo-Huei Lin ：指

指指

指導導導導教教教教授授：授授：：：陳鄰安陳鄰安陳鄰安陳鄰安 Advisor：：：：Dr. Lin-An Chen 指

指指

指導導導導教教教教授授：授授：：：洪志真洪志真洪志真洪志真 Advisor：：：：Dr. Jyh-Jen Horng Shiau

國

立

交

通

大

學

統計學研究所

博

士

論

文

A Thesis

Submitted to Institute of Statistics College of Science

National Chiao Tung University in Partial Fulfillment of the Requirements

for the Degree of Ph.D in

Statistics Jan 2012

Hsinchu, Taiwan, Republic of China

中華民國一百

中華民國一百零一年一

零一年一

零一年一月

月

(3)

包含區間

包含區間,

,,

,多變量管制圖及

多變量管制圖及

構面性資料統計品質管制

研究生：林碩慧指導教授：陳鄰安博士指導教授：洪志真博士國立交通大學統計學研究所摘要摘要摘要摘要包含區間一般以經驗分位數來估計，本論文提出以對稱型分位數來估計，並且在論文中可看出對稱型分位數包含區間在非對稱分佈及離群值資料上有比經驗分位數包含區間有較短的長度及較好的穩健度。在對稱分佈尤其是厚尾型對稱分佈, 對稱型分位數包含區間有較小的變異。本論文亦提出以對稱型分位數建構多分位點管制圖，並且探討其大樣本理論。本論文亦探討非線性混合性構面型資料統計品質管制。我们以主成份分析來建構統計品

質管制Phase-I及Phase-II的監控統計量。在Phase-I，我們採用主成份計分建構的 T2

監控統計量。在Phase-II，各別主成份計分圖、主成份計分建構的 T2

管制圖及聯合型主成份計分管制圖被提出及比較，在應用面亦有所建議。

(4)

Abstract

Classically the non-parametric coverage interval is estimated by empirical quan-tiles. We introduce an alternative way for estimating the coverage interval by symmetric quantiles of Chen and Chiang (1996). We further show that this alternative estimator has a better precision in the sense that its asymptotic variances are smaller than the classical one.

In an attempt to develop a scheme for monitoring a vector of distributional quantiles, we propose a symmetric-quantiles-based control chart. Compara-tive studies in terms of the asymptotic covariance matrix and the average run length show that the proposed control chart is more efficient than the classical empirical-quantiles-based control chart.

The monitoring of process/product profiles is presently a growing and promis-ing area of research in statistical process control. We focus on developpromis-ing mon-itoring schemes for nonlinear profiles with random effects in this study. We utilize the technique of principal components analysis to analyze the covariance structure of the profiles and propose monitoring schemes based on principal component (PC) scores. In the Phase I analysis of historical data, due to the

dependency of the PC-scores, we adopt the usual Hotelling T2 chart to check

the stability. For Phase II monitoring, we study individual PC-score control charts, a combined chart scheme that combines all the PC-score charts, and a

(5)

T2 chart. Although an individual PC-score chart may be perfect for monitoring a particular mode of variation, a chart that can detect general shifts, such as

the T2 chart and the combined chart scheme, is more feasible in practice. The

performances of the schemes under study are evaluated in terms of the average run length.

(6)

Chapter 1 Introduction

Statistical control charts have been widely used in many real-life applications and various control chart techniques have been developed to handle different process scenarios and characteristics of data in this literature. In this disserta-tion, two different types of control charts, a symmetric-quantile based control chart and a nonparametric profile control chart, are proposed and studied. In addition to control charts, other applications of coverage interval are also in-troduced.

1.1 Coverage Interval

In statistics, we generally have a random sample, X1, ..., Xn, from a

distri-bution with the probability density function (pdf) f (x, θ), x ∈ <, where < is the sample space and θ is an unknown parameter (scalar or vector) with the parameter space Θ. The statistical inferences presented in the literature

(10)

gen-erally involve estimation and hypothesis about the unknown parameter θ. In many practical problems, we face to deal with statistical problems for inferences about a probabilistic unknown interval or multivariate region, not this single point of parameter θ.

Let X be a random variable or a statistic which is formulated from the

random sample X1, ..., Xn. A coverage interval is an interval with two confidence

limits which covers values in the population of random quantity W in some probabilistic sense. For example, C(θ) = (a(θ), b(θ)) is a γ coverage interval

if it satisfies γ = Pθ(X ∈ C(θ)) for θ ∈ Θ. Treating γ coverage interval

C(θ) as the unknown parameter interval of interest for inferences has provided many very important applications in various areas of sciences but still lacking communications in theory and application (see Huang, Chen and Welsh (2010) for a brief review).

This coverage interval has been referred with various terminologies in differ-ent areas of applications, for example, the coverage interval C(θ) is called the reference interval in laboratory chemistry to refer to population-based reference values obtained from a well-defined group of reference individuals. Laboratory test results are commonly compared to a reference interval before care-givers make physiological assessments, medical diagnoses, or management decisions. An individual who is being screened for a disorder based on a measurement is suspected to be abnormal if his/her measurement value lies outside the reference interval.

(11)

An experiment for measuring a measurand, a particular quantity subject to measurement, is a method through a process that tries to gain or discover

knowledge of the measurand. Measurements always have errors and

there-fore uncertainties. General rules for evaluating and reporting uncertainty in measurement have been published by the most important and internationally widespread metrological publication – ISO (the International Standards Orga-nization) Guide to the Expression of Uncertainty in Measurement. The mea-surement result should be reported with a specified confidence as an uncertainty interval defining the range of values that could reasonably be attributed to the measurand. Hence, a coverage interval in physical science is called the uncer-tainty interval. This topic is interesting but not studied in this dissertation.

The control chart is one of the basic quality improvement tools in statistical process control (SPC) for a process to identify special causes of variation and signal the need to take necessary corrective actions. If the variation is due to common causes alone, the process is said to be in statistical control. When special causes are present, the process is said to be out of control. Suppose that W is the statistic to be monitored and (a(θ), b(θ)) is the γ coverage interval. By letting LCL = a(θ) and U CL = b(θ) be the lower control limit and upper control limit, respectively, we may say that special causes may be present when an observation of W falls outside the interval (LCL, U CL).

A control chart is a statistical scheme devised for the purpose of checking and then monitoring the statistical stability of a process. Various control chart

(12)

techniques have been developed and widely applied for process control. Widely used Shewhart control chart, EWMA, and CUSUM chart are under the para-metric (normal distribution) statistical thinking to perform process control by monitoring the parameters. The design of control chart is based on the concept of the coverage interval of a certain statistic such as sample mean or sample range under the assumption of a specified distribution, mostly the normal dis-tribution.

1.2 Coverage Interval by Symmetric Quantile

The coverage interval can be estimated either parametrically or non-parametrically. The parametric method classically assumes that the underlying distribution

of the measurement variable is normal whereas, recently, Chen, Huang and Chen (2007) has proposed a technique for constructing coverage intervals for asymmetric distributions. On the other hand, the non-parametric approach estimates the quantiles (percentiles) directly; the most popular technique for estimating the unknown quantiles is through the empirical quantiles. Most authorities now recommend the nonparametric method because it makes no as-sumptions concerning the type of the reference variable. It is easy and reliable whether the reference variable follows either a normal or non-normal distribu-tion (see these points in Reed, Henry and Mason (1971) and Solberg (1986)).

(13)

di-agnose the disease with precision. Some factors that may increase the precision have been considered. Number one hundred and twenty or more of healthy subjects required for the determination of coverage interval has been recom-mended by International Federation of Clinical Chemistry. The determination of the confidence interval of the quantile, that is, the limits within which true quantile is located with a specified confidence, is strongly recommended. How-ever, Friedberg et al. (2007) has observed that analytic imprecision is a very important factor for the quality of the established coverage interval. Hence, searching for an alternative technique in developing coverage interval to in-crease the analytic precision of the computed coverage interval is an interesting and important topic.

We consider the non-parametric approach in constructing coverage interval. Is there an alternative technique producing a coverage interval of better preci-sion than that constructed by empirical quantile? For improving the efficiency of the trimmed mean for estimating the location parameter, Kim (1992) and Chen and Chiang (1996) introduced the symmetric quantile to construct an alternative trimmed mean. They observed that this trimmed mean can have asymptotic variances very close to the Cramer Rao lower bounds for several distributions, including heavy tail ones. One aim in this research is to con-struct an alternative coverage interval by symmetric quantiles and show that it does gain better precision than the classical version constructed by empirical quantiles. This surprising result has been published in Metrologia (Lin, Chan

(14)

and Chen (2008)).

1.3 Multivariate Control Chart by Symmetric

Quan-tiles

It is known that the performance of a normal-based control chart is seriously

degraded when the underlying distribution departs from normal. However,

manufacturing processes with non-normal quality characteristics variable are very common (see, for example, Cheng and Thaga (2006), Shiling and Nelson (1976) and Kanji and Arif (2000)). Nonparametric control charts are then suggested because it makes no assumption concerning the distribution of the monitoring variable. For example, Janacek and Meikle (1997) considered the median chart, Liu and Tang (1996) considered the boostrap control chart, and Grimshaw and Alt (1997) considered using quantile function to construct the control chart.

Among the existing nonparametric control charts, the quantile- based control chart by Grimshaw and Alt (1997) is a nonparametric control chart that can simultaneously monitor distribution parameters. Unlike the confidence interval that estimates the range in which a population parameter falls, the control chart is a coverage interval that estimates the range which covers a certain percentage of the population of a specific statistic. Because coverage intervals are based upon only a sample of the entire population, we cannot be 100% confident that

(15)

the interval will contain the specified proportion of the statistic’s population. It is interesting to study if there is an alternative quantile technique that can be used to construct a quantile control chart of better efficiency in some sense than that constructed by empirical quantiles. Again, one of our aims in this research is to construct an alternative control chart by symmetric quantiles and show that it does gain better efficiency than the classical version constructed by empirical quantiles.

1.4 Nonparametric Profile Control Chart

Statistical process control has been widely applied in many areas, especially in industries. Another topic in this research is about profile monitoring. Most statistical process control applications deal with the quality of a process or product can be adequately represented by the distribution of a univariate qual-ity characteristic or a vector of correlated qualqual-ity characteristics. But in many practical situations, the focus point is the relationship between a response vari-able and some explanatory varivari-ables. Thus, at each sampling stage, a collection of data points that can represent the relationship is a data curve or data pro-file. And thus a control chart designed for profile data is called a profile control chart.

The monitoring of process/product profiles is presently a growing and promis-ing area of research in statistical process control. One of the aims in this

(16)

re-search is to develop monitoring schemes for nonlinear profiles with random effects. We utilize the technique of principal components analysis to analyze the covariance structure of the profiles and propose monitoring schemes based on principal component (PC) scores.

Profile monitoring is a relatively new research area in quality control. Kang and Albin (2000) studied the problem of linear profile monitoring and proposed two control schemes by modelling the profiles with the simple linear regression

model, Y = A0+A1x+, where Y is the response variable and x is the

indepen-dent variable; A0 and A1 are the parameters to be estimated; the noise variables

’s are independent and normally distributed with mean zero and common

vari-ance σ2. By centering the x-values to make the least squares estimators of the

Y -intercept (A0) and slope (A1) independent of each other, Kim et al. (2003)

proposed a combined-chart scheme in which three EWMA charts designed re-spectively for detecting shifts in intercept, slope, and standard deviation (σ) are used simultaneously. Mahmound and Woodall (2004) presented and compared several control charts for Phase I analysis of linear profiles and applied some of the charts to a calibration application. For more discussions on linear profile monitoring, see the review paper by Woodall et al. (2004).

Shiau and Weng (2004) extended the above linear profile monitoring schemes to a scheme suitable for profiles of more general forms via nonparametric re-gression. No assumptions are made for the form of the profiles except the smoothness. The nonparametric regression model considered is Y = g(x) + ,

(17)

where g(x) is a smooth function and is the random error as before. Spline regression was adopted as the curve fitting/smoothing technique for its sim-plicity. They proposed an EWMA chart for detecting mean shifts, an R chart for variation changes, and an EWMSD (standard deviation) chart for variation increases.

Note that the models described above all consist of a deterministic line/curve plus random noises. It does not account for some allowable profile-to-profile variations that we often observe in many profile data, e.g., the aspartame ex-ample and VDP exex-ample, where these profile-to-profile variations should be con-sidered as caused by common causes. A monitoring scheme constructed based on the afore-mentioned “fixed-effect” model may interpret these common-cause variations as caused by some special causes and signal many false alarms. Thus, we need a suitable model that can cope with these common-cause variations and construct a monitoring scheme accordingly.

For this, Shiau, Lin, and Chen (2006) considered a random-effect linear model to develop monitoring schemes for linear profiles. Similarly, Jensen et al. (2006) proposed a linear mixed (effects) model for linear profiles. Williams et al. (2003) fitted nonlinear profiles by nonlinear parametric regression and

then monitored profiles with some T2 statistics of the estimated parameters.

Later, Williams et al. (2007) extended this methodology to nonlinear profiles with a nonconstant variance at set points to analyze a set of heteroscedastic dose-response profiles. Adopting a random-effect parametric nonlinear

(18)

regres-sion model for profiles, Shiau, Yen, and Feng (2006) proposed a robust nonlinear profile monitoring scheme. Jensen et al. (2009) proposed using nonlinear mixed models to model nonlinear profiles. Note that the parametric approaches men-tioned above all need to pre-specify a parametric functional form for profiles, a task often not so easy for practitioners. Qiu, Zou and Wang (2010) proposed a novel control chart, which dealt with mixed effect profile data without assump-tions on a parametric functional form. Their control chart is based on local linear kernel smoothing of profile data and on the EWMA weighting scheme at different time points and the heteroscedasticity of observations within each profile.

We extend the nonparametric fixed-effect model of Shiau and Weng (2004) to a random-effect model in order to incorporate some profile-to-profile variabil-ity as caused by common causes. With the random-effect model, we focus on the covariance structure and use the principal components analysis (PCA) to analyze it. Ding et al. (2006) also considered modelling profiles nonparametri-cally for a Phase I analysis, but proposed using ICA (independent components analysis) instead of PCA for monitoring profiles that are in clusters, a situation PCA may fail to preserve the clustering feature of the original data.

PCA is very useful in summarizing and interpreting a set of profile data with the same equally spaced x-values for each profile. We remark that the smoothing step described above can relax this requirement for profile data since the equally-spaced data can be obtained from the smoothed profiles easily.

(19)

Some pioneer works on analyzing curves with PCA include Castro et al. (1986), Rice and Silverman (1991), Jones and Rice (1992), and others. For applications, Shiau and Lin (1999) analyzed a set of accelerated LED degradation profiles to estimate the mean lifetime of the product with the techniques of nonparametric regression and PCA.

For Phase I profile monitoring, we propose using the usual Hotelling T2

chart, a commonly used control chart designed for multivariate process data, by treating the principal component (PC) scores of a profile obtained from PCA as the multivariate data. For Phase II process monitoring, we propose and study three monitoring schemes constructed by utilizing the eigenvalues and eigenvectors obtained from PCA to compute the PC-scores of each incoming profile, including individual PC-score charts, a combined chart that combines

all of the PC-score charts and a T2 chart (different from the T2 chart of Phase

I). The performances of these monitoring schemes are evaluated in terms of the average run length (ARL).

(20)

Chapter 2 A Nonparametric Coverage Interval

In this chapter, we propose a coverage interval estimation based on sym-metric quantiles. The coverage precision of the proposed symsym-metric coverage interval is studied and comparisons with empirical quantile coverage interval is also conducted to demonstrate that the symmetric quantile coverage interval is superior to its empirical counterpart in practical usage when the underlying distribution is symmetric.

2.1 Symmetric Coverage Interval

For random variable Y with cumulative distribution function F , the λth quantile is defined as

(21)

The classical 1 − α coverage interval is defined as C(1 − α) = (F−1(α 2), F −1 (1 − α 2)).

Suppose that we now have a random sample y1, ..., yn from distribution F . The

corresponding empirical 1 − α coverage interval is

Cn(1 − α) = (Fn−1( α 2), F −1 n (1 − α 2)) (2.1)

where we let F_n−1 be the empirical quantile, a quantile function with distribution

function of the sample type as Fn(y) = 1_nPn_i=1I(yi ≤ y).

Unlike the way in which the empirical quantile is constructed based on the cumulative distribution function, the so-called symmetric quantile of Chen and Chiang (1996) is formulated based on a folded distribution function. Let us consider the folded cumulative function about µ, known or unknown, as

Fs(a) = P (|y − µ| ≤ a), a ≥ 0.

Extending from Chen and Chiang (1996), we define the 1 − α symmetric cov-erage interval as

Cs(1 − α) = (µ − Fs−1(1 − α), µ + F

−1

s (1 − α))

where F_s−1(λ) = inf{a : Fs(a) ≥ λ}. If F is continuous, the 1 − α symmetric

coverage interval satisfies 1 − α = P (µ − F_s−1(1 − α) ≤ y ≤ µ + F_s−1(1 − α)). If

we further assume that F is symmetric at µ, it can be seen that

(22)

the classical one and the symmetric one are identical in the sense of containing the same set of reference individuals.

We interpret the folded cumulative function and the symmetric coverage interval through Figure2.1.

0 2 4 6 8 10 0.0 0.1 0.2 0.3 0.4 µµ − Fs−−1₍₍_0.8₎₎ µµ + Fs−−1((0.8)) µµ µµ+a µµ−a

Folded distribution function and symmetric coverage interval for Gamma distribution ΓΓ(2,2)

Figure 2.1: The symmetric coverage interval.

Considering the Gamma distribution Γ(2, 2) which has probability density function (pdf) as the curve in the figure, we consider the folded distribution about the median. With this distribution, the median µ is 3.36. For a given a > 0, the value of this folded distribution at a represents the probability of a region as the part of shadow. Suppose that the interest is 80% coverage

interval. For this continuous distribution, we search F_s−1(0.8) = a∗ such that

(23)

Hence the 80% symmetric coverage interval is

Cs(0.8) = (µ − Fs−1(0.8), µ + F

−1

s (0.8))

= (1.15, 5.57)

Let ˆµ be an estimate of µ. We may define the sample type 1 − α symmetric

coverage interval as Csn(1 − α) = (ˆµ − Fsn−1(1 − α), ˆµ + F −1 sn (1 − α)) (2.3) where Fsn(a) = _n1 Pn

i=1I(|yi − ˆµ| ≤ a) is the sample type folded cumulative

distribution function and F_sn−1(1 − α) = inf{a : Fsn(a) ≥ 1 − α}.

Let’s give a simple example to describe the construction of sample symmetric coverage interval. Suppose that we have a set of observations that are ordered as

−5, −3, −2, −1, −0.5, 0.5, 1, 3, 50, 100.

We want to construct 80% empirical and symmetric coverage intervals. With

F_n−1(0.1) = −5 and F_n−1(0.9) = 50, the 80% empirical coverage interval is

Cn(0.8) = (−5, 50). (2.4)

For construction of symmetric coverage interval, we choose sample median as the estimate of µ. That is,

ˆ µ = F_n−1(0.5) = inf{a : 1 10 10 X i=1 I(yi ≤ a) ≥ 0.5} = −0.5.

Let’s denote residuals ei = yi− ˆµ, i = 1, ..., 10. The residuals are

(24)

The sample type folded cumulative distribution function is Fsn(a) = 1 10 10 X i=1 I(|ei| ≤ a).

For examples, Fsn(0) = ₁₀1, Fsn(1) = ₁₀1[I(| − 0.5| ≤ 1) + I(|0| ≤ 1) + I(|1| ≤

1)] = ₁₀3. Then we have F_sn−1(0.8) = inf{a : 1 10 10 X i=1 I(|ei| ≤ a) ≥ 0.8} = 4.5.

This indicates that the 80% symmetric coverage interval is

Csn(0.8) = (ˆµ − Fsn−1(0.8), ˆµ + F −1

sn (0.8))

= (−0.5 − 4.5, −0.5 + 4.5)

= (−5, 4). (2.5)

Comparing the resulted sample empirical and symmetric coverage intervals in (2.4) and (2.5), it is seen the benefit for using the latter one for that it is shorter than the former one. This would happen very often when the observa-tions are drawn from asymmetric distribuobserva-tions.

(25)

2.2 Precision Study of Symmetric Coverage Interval

The equality of (2.2) does not hold when the underlying distribution F is not symmetric so that there is no fair criterion to compare their corresponding sample coverage intervals. Hence, we may set the case that F is symmetric to compare the precision of these two coverage intervals through the asymptotic variances of their sample type coverage intervals.

We consider that µ is the median parameter and let ˆµ be the sample median

as ˆ µ = arginf_µ∈R n X i=1 |yi − µ|.

Suppose that we assume that F is continuous and symmetric at µ. From

Ruppert and Carroll (1980), we have a Bahadur representation for this sample median as n1/2(ˆµ − µ) = n−1/2 1 f (µ) n X i=1 (0.5 − I(yi ≤ µ)) + op(1). (2.6)

On the other hand, a Barhadur representation for F_sn−1(1 − α) developed by

Chen and Chiang (1996) is

n1/2(F_sn−1(1 − α) − (F−1(1 − α 2) − µ)) = 1 2f (F−1_{(1 −} α 2)) n−1/2 n X i=1 {1 − α − I(F−1(α 2) ≤ yi ≤ F −1_{(1 −} α 2)} + op(1). (2.7)

The assumption of symmetric distribution indicates that ˆµ − F_sn−1(1 − α) and

ˆ

(26)

we have a Bahadur representation for ˆµ − F_sn−1(1 − α) as n1/2((ˆµ − F_sn−1(1 − α)) − F−1(α 2)) = n −1/2 n X i=1 {[− 1 2f (µ) − 1 − α 2f (F−1_{(1 −} α 2)) ] I(yi ≤ F−1( α 2)) + [− 1 2f (µ) + α 2f (F−1_{(1 −} α 2)) ]I(F−1(α 2) ≤ yi ≤ µ) + [ 1 2f (µ) + α 2f (F−1_{(1 −} α 2)) ]I(µ < yi ≤ F−1(1 − α 2)) + [ 1 2f (µ) − 1 − α 2f (F−1_{(1 −} α 2)) ]I(yi ≥ F−1(1 − α 2))} + op(1). (2.8)

Since y1, ..., yn is a random sample from distribution F , we may see that the

asymptotic variance of n1/2(ˆµ − F_sn−1(1 − α) − F−1(α₂)) is σ_s2 = α 2[( 1 2f (µ) + 1 − α 2f (F−1(1 − α₂))) 2_{+ (} 1 2f (µ) − 1 − α 2f (F−1(1 − α₂))) 2_{] + (1 −} α 2)[(− 1 2f (µ) + α 2f (F−1(1 − α₂))) 2 + ( 1 2f (µ) + α 2f (F−1_{(1 −} α 2)) )2]. (2.9)

On the other hand, in this situation that y has a continuous and symmetric

distribution, we may see that n1/2(F_n−1(α₂) − F−1(α₂)) and n1/2(F_n−1(1 − α₂) −

F−1(1 −α₂)) also have the same asymptotic variance (see, for example, Sen and

Singer (1993, p168)) as σ2_e = α 2(1 − α 2)f −2_(F−1_{(1 −} α 2)). (2.10)

Since these two sample coverage intervals estimate the same population cov-erage interval, it is fair that we evaluate the efficiency of the symmetric type

(27)

coverage interval defined as the following Ef f = σ 2 e σ2 s . (2.11)

Let’s consider several distributions for computation of asymptotic variances of (2.9) and (2.10) to compare their corresponding efficiencies of (2.11) where distributions include standard normal distribution N (0, 1), t-distribution t(r) where r is the degrees of freedom, Cauchy distribution (Cauchy(s), s > 0) with pdf

f (y) = 1

π s

y2_{+ s}2, y ∈ R

and the Laplace distribution (Lap(b)) with pdf

f (y) = 1

2be

−|y|_b _{, y ∈ R.}

(28)

Table 2.1. The efficiencies, Ef f , of symmetric coverage interval dist/1 − α 0.6 0.8 0.9 0.95 0.98 N (0, 1) 0.87 0.87 1.02 1.21 1.48 t(r) r = 1 0.98 1.78 2.13 2.1 2.04 r = 5 0.89 1.01 1.31 1.61 1.85 r = 10 0.88 0.94 1.16 1.42 1.7 Cauchy(s) s = 1, 5, 10 0.98 1.78 2.13 2.1 2.04 Lap(b) b = 1, 5, 10 1.2 1.6 1.8 1.9 1.96

It is relatively efficient to use the empirical quantile to construct coverage inter-val when the quantile percentage is close to 0.5 in either direction. This means that when we want a 1 − α coverage interval with coverage probability 1 − α as the value of 0.6 or even smaller. The one estimated by empirical quantiles is the right choice. On the other hand, we see that it gains more precision to use symmetric quantile to construct coverage interval when 1 − α is with the value of 0.8 or larger. This alternative coverage interval is then attractive since it is very common that we apply coverage interval only for large 1 − α, for example, the reference interval in medical diagnosis chooses the value of 0.95. In fact, the case that when the underlying distribution is the Laplace one the coverage interval constructed by symmetric quantiles totally dominate the one

(29)

by empirical quantiles.

This interesting result is not surprising. The surprising fact is that, unlike estimation of location and scale parameters that have been received very much attention in statistical literature for proposing techniques and developing the-ories in gaining better precisions, the attention for developing alternative ways in constructing coverage intervals for gaining better precision than the classical one has not been paid in statistical and metrological literatures.

2.3 Concluding Remarks

As for constructing a coverage interval, a symmetric quantile coverage in-terval performs better than the commonly used empirical quantile coverage interval both empirically and theoretically. The robust estimation on the dis-tribution mass and the symmetric folding feature of the symmetric quantile coverage interval prevent the estimation of coverage interval heavily impacted by outliers. When an underlying distribution is not symmetric, with center on the distribution median, using symmetric quantiles to construct a coverage interval can cover more percentage of the higher density part compared with the other side with lower distribution density, and thus a symmetric quantile coverage interval gives a shorter coverage interval than the empirical one. Even when the distribution is symmetric, we still can see the symmetric quantile coverage interval is with smaller asymptotical variance than the empirical one

(30)

(31)

Chapter 3 Multivariate Control Chart by

Symmetric Quantiles

In this chapter, we apply symmetric quantile techniques in constructing con-trol chart schemes. A multivariate concon-trol chart by symmetric quantiles is proposed and its asymptotic characteristic is studied. Comparisons between empirical quantile based control chart and symmetric quantile control chart are also studied through efficiency in terms of asymptotic variances and ARL. Symmetric quantile based multivariate control chart is more efficient than the empirical quantile based multivariate control chart in detecting small distribu-tional shifts even when underlying distribution is not symmetric.

(32)

3.1 Empirical Quantiles Based Control Chart for

Mon-itoring Quantile Vector

With the interest of control charts comparison, let us first introduce, in de-tail, the empirical quantiles based control chart of Grimshaw and Alt (1997).

For percentages α1, ..., αk with α1 < α2 < ... < αk, let us consider the

popula-tion quantile vector

Q(α1, ..., αk) =           F−1(α1) F−1(α2) ... F−1(αk)           (3.1)

for monitoring that can be estimated by the corresponding empirical quantile vector Qe(α1, ..., αk) =           F_n−1(α1) F_n−1(α2) ... F_n−1(αk)           .

Grimshaw and Alt (1997) proposed to apply Qe(α1, ..., αk) to monitor the

pop-ulation quantile vector Q(α1, ..., αk). The asymptotic property of Qe(α1, ..., αk)

relies on the empirical quantiles F_n−1(αj)’s. A Bahhadur representation for

F_n−1(α) (see, for example, Ruppert and Carroll (1980)) is

√ n(F_n−1(α) − F−1(α)) = 1 f (F−1(α))n −1/2 n X i=1 (α − I(yi ≤ F−1(α))) + op(1)

(33)

where f is the probability density function of variable y and I is the indicator

function. This leads to that n1/2(Qe(α1, ..., αk) − Q(α1, ..., αk)) converges to a

normal distribution Nk(0k, Σ) with Σ = (σjke ) and where σije =

αi(1−αj) f (F−1_(α_i_{))f (F}−1_(α_j₎₎ for i ≤ j. This further implies that the following

n(Qe(α1, ..., αk) − Q(α1, ..., αk))0Σ−1(Qe(α1, ..., αk) − Q(α1, ..., αk)) → χ2(k) holds asymptotically in distribution.

Suppose that we have a training sample yij, i = 1, ..., n, j = 1, ..., m that

rep-resents an in-control data set of m samples of size n from distribution F so that

estimate of Q(α1, ..., αk) and Σ are available. Generally we let Qej(α1, ..., αk)

and ˆΣj be estimates, respectively, based on sample yij, i = 1, ..., n and define

Q0(α1, ..., αk) = _m1 P m j=1Qej(α1, ..., αk) and Σ0 = 1 m Pm j=1Σˆj. Treated

esti-mates Q0 and Σ0 as true values of Q(α1, ..., αk) and Σ, the control statistic and

upper control limit proposed by Grimshaw and Alt (1997) are

Control statistic Te = n(Qe(α1, ..., αk) − Q0(α1, ..., αk))0Σ−10 (Qe(α1, ..., αk) − Q0(α1, ..., αk))

U CLe = χ2α (3.2)

(3.3)

where χ2_α satisfies 1 − α = P (χ2(k) ≤ χ2_α). With this proposal, if a sample

point y1, ..., yn has a value of Te lying below the upper limit U CLe and a set

of T_e0s does not exhibit any systematic pattern, we say that the process is in

(34)

With Q(α1, ..., αk) of (3.1) as the target to be monitored, the idea behind this study is why should we estimate it by the empirical quantiles.

3.2 Symmetric Quantile Control Chart

From the study in chapter 2, we anticipate that using the symmetric

quan-tiles to estimate the distributional characteristic vector Q(α1, ..., αk) is more

efficient than that of empirical quantiles. Considering a number ` decreasing

percentages γ1 > γ2 > ... > γ`, we define its corresponding 2` symmetric

quan-tile vector and population symmetric quanquan-tiles, respectively, as

Qsn(γ1, ..., γ`) =                         Fsn(−)(γ1) Fsn(−)(γ2) ... Fsn(−)(γ`) Fsn(+)(γ`) Fsn(+)(γ`−1) ... Fsn(+)(γ1)                         and Qs(γ1, ..., γ`) =                         Fs(−)(γ1) Fs(−)(γ2) ... Fs(−)(γ`) Fs(+)(γ`) Fs(+)(γ`−1) ... Fs(+)(γ1)                         .

where, the γ symmetric quantile pair is defined as

{F_s(−)(γ), F_s(+)(γ)} = {µ − F_s−1(γ), µ + F_s−1(γ)} (3.4)

(35)

µ. We may define the sample type γ symmetric quantile pair as

{F_sn(−)(γ), F_sn(+)(γ)} = {ˆµ − F_sn−1(γ), ˆµ + F_sn−1(γ)} (3.5)

with sample folded distribution function Fsn(a) = _n1Pn_i=1I(|yi − ˆµ| ≤ a)

and

F_sn−1(γ) = inf{a : Fsn(a) ≥ γ}. (3.6)

From Chen and Chiang (1996), we may see that n1/2(Qsn(γ1, ..., γ`)−Qs(γ1, ..., γ`))

is asymptotically normal N2`(02`, Σs) for some matrix Σs that will be given

ex-plicitly latter. This further implies that the following

n(Qsn(γ1, ..., γ`) − Qs(γ1, ..., γ`))0Σ−1s (Qsn(γ1, ..., γ`) − Qs(γ1, ..., γ`)) → χ2(2`)

holds asymptotically in distribution.

Again, from the training sample yij, i = 1, ..., n, j = 1, ..., m that

repre-sents an in-control data set of m samples of size n, we let Qs0(γ1, ..., γ`) =

1 m Pm j=1Qsn,j(γ1, ..., γ`) and Σs0 = 1 m Pm j=1Σˆs,j where Qsn,j(γ1, ..., γ`) and ˆΣs,j

are estimates of, respectively, Qs(γ1, ..., γ`) and Σs. Let us denote these two

(36)

and upper control limit as

Control statistic Ts = n(Qsn(γ1, ..., γ`) − Qs0(γ1, ..., γ`))0Σ−1s0 (Qsn(γ1, ..., γ`) − Qs0(γ1, ..., γ`))

U CL = χ2_α(2`) (3.7)

The asymptotic covariance matrix Σs varies with the predetermined

estima-tor of µ. From a comparison in Chen and Chiang (1996), although various estimators of µ lead to different asymptotic distributions for their correspond-ing symmetric quantiles, however, their performances in constructcorrespond-ing statistical procedures such as trimmed means are very competitive. Hence from here

af-ter, for simplicity, we consider that µ is the median parameter and let ˆµ be the

`1-norm estimate as ˆ µ = arginf_µ∈R n X i=1 |yi − µ|.

A description of the asymptotic covariance matrix Σs under this setting of

predetermined estimator will be developed in next subsection.

3.3 Derivation of Asymptotic Covariance Matrix of

Sym-metric Quantile Vector

We assume that F is continuous and symmetric at median µ. From Ruppert

(37)

as n1/2(ˆµ − µ) = n−1/2 1 f (µ) n X i=1 (0.5 − I(yi ≤ µ)) + op(1). (3.8)

Furthermore, a Barhadur representation for F_sn−1(γ) of (3.6) that has been

de-veloped by Chen and Chiang (1996) is

n1/2(F_sn−1(γ) − F−1(γ)) = 1 2f (F−1(1+γ₂ ))n −1/2 n X i=1 (γ − I(−F−1(1 + γ 2 ) ≤ yi ≤ F −1₍1 + γ 2 )) + op(1). (3.9)

From (3.8) and (3.9), Bahadur representations, respectively, for Fsn(−)(γ) and

Fsn(+)(γ) are n1/2(F_sn(−)(γ) − F−1(1 − γ 2 )) = n −1/2 n X i=1 {[− 1 2f (µ) − γ 2f (F−1₍1+γ 2 )) ] I(yi ≤ F−1( 1 − γ 2 )) + [− 1 2f (µ) + 1 − γ 2f (F−1₍1+γ 2 )) ]I(F−1(1 − γ 2 ) ≤ yi ≤ µ) + [ 1 2f (µ) + 1 − γ 2f (F−1₍1+γ 2 )) ]I(µ < yi ≤ F−1( 1 + γ 2 )) + [ 1 2f (µ) − γ 2f (F−1₍1+γ 2 )) ]I(yi ≥ F−1( 1 + γ 2 ))} + op(1). (3.10) and n1/2(F_sn(+)(γ) − F−1(1 + γ 2 )) = n −1/2 n X i=1 {[− 1 2f (µ) + γ 2f (F−1(1+γ₂ ))] I(yi ≤ F−1( 1 − γ 2 )) + [− 1 2f (µ) − 1 − γ 2f (F−1₍1+γ 2 )) ]I(F−1(1 − γ 2 ) ≤ yi ≤ µ) + [ 1 2f (µ) − 1 − γ 2f (F−1₍1+γ 2 )) ]I(µ < yi ≤ F−1( 1 + γ 2 )) + [ 1 2f (µ) + γ 2f (F−1₍1+γ 2 )) ]I(yi ≥ F−1( 1 + γ 2 ))} + op(1). (3.11)

(38)

By letting ` = k₂ and α1, ..., αk satisfying 1−α1 = αk, 1−α2 = αk−1, ..., 1−αk/2 = αk

2+1, then Qs(1−2α1, 1−2α2, ..., 1−2α

k

2) = Q(α1, α2, ..., αk). This allows us to

compare estimators of Q(α1, ..., αk) that are constructed by empirical quantiles

and symmetric quantiles. For αi < αj, we need to develop the asymptotic

variances and covariances for Fsn(−)(1 − 2αi) and F

(+)

sn (1 − 2αj) for i 6= j. Let’s

denote σ_α2₁, σ_1−α2 ₁, σα11−α2 as the asymptotic variance of n

1/2_(F(−)

sn (1 − 2α1) −

F−1(α1)) and n1/2(F

(+)

sn (1 − 2α1) − F−1(1 − α1)) and asymptotic covariance of

n1/2(Fsn(−)(1 − 2α1) − F−1(α1)) and n1/2(F (+)

sn (1 − 2α2) − F−1(1 − α2)). Since

y1, ..., yn is a random sample from distribution F , we may derive the followings

from (3.10) and (3.11), σ_α2₁ = σ_1−α2 ₁, σα11−α2 = σ1−α2α1, σα1α2 = σα2α1, σ1−α11−α2 = σ1−α21−α1, σ_α2₁ = ( 1 2f (µ)) 2_{+ 2α} 1(1 − 2α1)( 1 2f (F−1(1 − α1)) )2, σα1α2 = ( 1 2f (µ)) 2_{+ 2α} 1(1 − 2α2)( 1 4f (F−1_{(1 − α} 1))f (F−1(1 − α2)) ), (3.12) σα11−α2 =      (_{2f (µ)}1 )2− 2α1(1 − 2α2)(_{4f (F}−1_(1−α₁_{))f (F}1 −1_(1−α₂₎₎) if α1 < α2 (_{2f (µ)}1 )2− 2α2(1 − 2α1)(_{4f (F}−1_(1−α 1 1))f (F−1(1−α2))) if α1 > α2 , σ1−α11−α2 = ( 1 2f (µ)) 2 + 2α1(1 − 2α2)( 1 4f (F−1(1 − α1))f (F−1(1 − α2)) ).

With careful arrangement and derivations from (3.8) - (3.12) and denoting αij =

αi(1 − αj) and fij = f (F−1(1 − α₂i))f (F−1(1 −

αj

(39)

matrix for the symmetric quantile vector Qsn(1 − 2α1, 1 − 2α2, ..., 1 − 2α`) is stated in the following theorem.

Theorem 3.3.1 When the distribution F is symmetric about its median µ, then the asymptotic covariance matrix for the symmetric quantile vector is

Σs = 1 4f (µ)211 0₊ 1 4              A1 A0₁       A2 A3       A0₃ A0₂       A4 A0₄              where A1 =               α11 f11 α12 f12 . . . α1(`−1) f1(`−1) α1` f1` α22 f22 . . . α2(`−1) f2(`−1) α2` f2` ... ... ... α(`−1)(`−1) f(`−1)(`−1) α(`−1)` f(`−1)` α`` f``               A2 =               −α1` f1` − α1(`−1) f1(`−1) . . . − α12 f12 − α11 f11 −α2(`−1) f2(`−1) . . . − α22 f22 − α21 f21 ... ... ... −α(`−1)2 f(`−1)2 − α(`−1)1 f(`−1)1 −α`1 f`1              

(40)

A3 =               −α`2 f`2 −α`3 f`3 − α(`−1)3 f(`−1)3 ... ... −α`(`−1) f`(`−1) − α(`−1)(`−1) f(`−1)(`−1) . . . − α3(`−1) f3(`−1) −α`` f`` − α(`−1)` f(`−1)` . . . − α3` f3` − α2` f2`               A4 =               α`` f`` α(`−1)` f(`−1)` . . . α2` f2` α1` f1` α(`−1)(`−1) f(`−1)(`−1) . . . α2(`−1) f2(`−1) α1(`−1) f1(`−1) ... ... ... α22 f22 α12 f12 α11 f11               , and ` = k 2.

3.4 A Comparison between Empirical and Symmetric

Quantile Control Charts

When the interest is of constructing a symmetric quantiles based control chart, we wish to establish evidences for supporting the use of this new control chart. Since the symmetric quantile and the empirical quantile are both asymp-totically normal and both are consistent for a same population quantile vector

Q(α1, ..., αk), one interesting criterion in comparing these two control charts

is to study the efficiencies in terms of the ratio of sizes of their asymptotic covariance matrices.

(41)

Grimshaw and Alt (1997) pointed out that for effective use of a quantile chart in detection of distributional shift we should select quantile percentages

αi’s so that their corresponding quantile differences F_O−1(αi) − F_I−1(αi), with

FO and FI respectively representing the distribution function of in-control and

likely out-of-control processes, are large. For verification of this concern, Chiang

et al. (2006) showed that α0_is should be lie away of 0.5. We consider comparing

the estimator of quantile pair    F−1(α) F−1(1 − α)   . (3.13)

The asymptotic covariance matrices of the estimators of empirical quantile and symmetric quantile provide us to define the efficiency as

Ef f = Det(Cov    F_n−1(α) F_n−1(1 − α)   ) Det(Cov    Fsn(−)(1 − 2α) Fsn(+)(1 − 2α)   ) = α(1 − α)/(f (F −1_{(α)f (F}−1_{(1 − α))} σ4 α− σα1−α2 (3.14) where the covariance matrix of empirical quantiles is decribed in subsection 3.3

and that of the symmetric quantiles and notations of σ_α4 and σ2_α1−α are listed in

(3.12). When Ef f > 1, the estimation of quantile vector (3.13) is more efficient by the symmetric quantiles. On the other hand, if Ef f < 1, it prefers to estimate quantile vector by empirical quantiles. For computation of efficiencies of (3.14), we consider the distributions including normal distribution, Cauchy

(42)

distribution, Laplace distribution (Lap(b)) with pdf

f (y) = 1

2be

−|y|_b _{, y ∈ R}

and the Cauchy distribution with pdf

f (x, δ, γ) = 1

π[

γ

(x − δ)2_{+ γ}2], x ∈ R.

From (3.14), the computation of efficiency requires only the density function and its corresponding distribution function. We display the resulted efficiencies in Table 3.1.

Table 3.1. The efficiencies, Ef f , for symmetric quantiles based estimator

Dist/1 − 2α 0.6 0.7 0.8 0.9 0.95 0.98 Normal 0.85 0.79 0.81 1.03 1.50 2.72 Laplace 1.25 1.67 2.5 5 10 25 Cauchy 0.98 1.4 3.35 21.93 166.98 2573 t(r) r = 1 0.98 1.4 3.35 21.93 166.98 2573 r = 2 0.9 1.01 1.53 4.29 14.58 83.01 r = 5 0.87 0.86 1.02 1.75 3.54 10.24

The estimation of quantile vector (3.13) by empirical quantiles is more effi-cient when the quantile percentage is 0.6. However, it is impressed that it gains more precision to use symmetric quantile to construct the quantile vector esti-mator when percentage is equal or more than 0.8. In fact, the case that when

(43)

the underlying distribution is the Laplace one the estimator of quantile vec-tor constructed by symmetric quantiles totally dominate the one by empirical quantiles.

3.5 Comparison of Average Run Length for Two

Quan-tile Charts

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Σ−1₀ (Qn(α) − Q0(α)) → χ2(k) in distribution.

Hence we actually consider the following hypothesis

(44)

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Σ₀−1(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_Σ−1

0 (Qn(α) − Q0(α)) ≥ χ2α(k))

= P_Q_aX+b_(α),Σ_aX+b(n(Qn(α) − QaX+b(α) + (−(1 − a)Q0(α) + b))0Σ−1_aX+b

(Qn(α) − QaX+b(α) + (−(1 − a)Q0(α) + b)) ≥ 1 a2χ 2 α(k)) = P (χ2(k, n a2(−(1 − a)Q0(α) + b) 0_Σ−1 0 (−(1 − a)Q0(α) + b)) ≥ 1 a2χ 2 α(k)) (3.16) (3.17)

since√nΣ−1/2_aX+b0(Qn(α)−QaX+b(α)+(−(1−a)Q0(α)+b))

app ∼ Nk( √ nΣ−1/2_aX+b0(−(1− a)Q0(α) + b), Ik) = Nk( √ n a Σ −1/20

(45)

noncentral chi-square distribution with non-centrality parameter λ. The ARL

for quantile vector Qn(α) is

ARL = 1 P (χ2_k(λ = _an2(−(1 − a)Q0(α) + b)0Σ −1 0 (−(1 − a)Q0(α) + b)) ≥ 1 a2χ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.

(46)

(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 ARL0_ss

(47)

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.

(48)

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.

(49)

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, σ2_N), ∼ 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+ (µ 2 M + σ 2 M) eµN[(xi−1)2+(xj−1)2]+ σ2_{N [(xi−1)}2+(x_{j −1)}2]2 2 −µ2_MeµN(xi−1)2+ σ2_{N (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

(50)

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

(51)

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 y_i and the usual sample covariance matrix

of {y_i, 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 ≡ v

0

ryi 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

(52)

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

T_i2 = (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)2T 2 i ∼ 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

(53)

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 v0_rµ₀ 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

v0_rµ₀± Zα/2 √

(54)

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

max

1≤r≤K|

Sr− v_√ 0rµ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 χ2_K) when the process is in control. Thus, the upper control limit is the

100(1 − α) percentile of χ2_K.

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

(55)

mean of the profile has been shifted from µ₀ to µ₀ + δ. The probability of detecting the shift by the r-th PC-score chart is

p = 1 − P (|Sr − v 0 rµ0 √ λr | ≤ Zα/2) = 1 − P (− v0_rδ √ λr − Zα/2 ≤ Z ≤ − v0_rδ √ λr + Zα/2) = 1 − Φ(−v 0 rδ √ λr + Zα/2) + Φ(− v0_rδ √ λr − Zα/2),

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 1≤r≤K| Sr− v_√ 0rµ0 λr | ≤ Zα0_/2) = 1 − K Y r=1 P (|Sr− v 0 rµ0 √ λr | ≤ Zα0_/2) = 1 − K Y r=1 [Φ(−v 0 rδ √ λr + Zα0_/2) − Φ(− v0_rδ √ λr − Zα0_/2)],

where α0 = 1 − (1 − α)1/K is the individual false-alarm rate.

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 > χ2_K,α) = P (χ2_K(ξ) > χ2_K,α),

where χ2_K,α denotes the 100(1 − α) percentile of the central chi-square

(56)

4.2 Simulation and Comparative Studies

4.2.1 Settings for Simulation

In our simulation study, we generate profiles from M V N (µµµ₀, ΣΣΣ), where µµµ₀ is

given in (4.4) and ΣΣΣ is given in (4.2) with µI = 1, σI = 0.2, µM = 15, σM = 1,

µN = −1.5, σN = 0.3, and x = 0.64, 0.8, . . . , 3.52. Both of x and y values are

scaled variables, not the actual temperature levels and the amount of aspartame

dissolved in the dissolving process. Denote the in-control ARL by ARL0. All

charts are designed to have the same ARL0 = 370.3704, which corresponds to

the false-alarm rate of α = 0.0027.

4.2.2 A Study on the Number of Principal Components

To study how the choice of the number of effective principal components affects the detecting power of the monitoring scheme, we conduct a simulation study. In this study, the detecting power is measured by the ability of the monitoring scheme in detecting the real out-of-control profiles. For example, in a data set of fifty simulated profiles with three out-of-control profiles, if the scheme catches two of the three, then the detecting power measured is 2/3. The false-alarm rate can be measured in a similar way. Let the number of the principal components used be k. Then for each data set, compute the detect-ing power and the percentage of the total variation explained by k principal components for various values of k.

包含區間, 多變量管制圖 及構面性資料統計品質管制

國 立 交 通 大 學

統計學研究所

統計學研究所

統計學研究所

統計學研究所

博

博

博

博士論文

士論文

士論文

士論文

包含區間

包含區間

包含區間

包含區間,

,

,

, 多變量管制圖

多變量管制圖

多變量管制圖

多變量管制圖

及

及

及

及

構面性資料統計品質管制

構面性資料統計品質管制

構面性資料統計品質管制

構面性資料統計品質管制

Coverage Interval, Multivariate Control Chart

Coverage Interval, Multivariate Control Chart

Coverage Interval, Multivariate Control Chart

Coverage Interval, Multivariate Control Chart

and Profile Monitoring

and Profile Monitoring

and Profile Monitoring

and Profile Monitoring

研

研

研

研

究

究

究

究

生

生 :

生

生

:

:

: 林碩慧

林碩慧

林碩慧

林碩慧

指導教授

指導教授

指導教授

指導教授 :

:

:

: 陳鄰安

陳鄰安

陳鄰安

陳鄰安

教授

教授

教授

教授

指導教授

指導教授

指導教授

指導教授 :

:

:

: 洪志真

洪志真

洪志真

包含區間, 多變量管制圖及構面性資料統計品質管制

國立交通大學