用無母數迴歸方法監控曲線型品質特性之製程

國

立

交

通

大

學

統計學研究所

碩

士

論

文

用無母數迴歸方法監控曲線型品質特性之製程

Profile Monitoring by Nonparametric Regression

研 究 生：翁崢珮

指導教授：洪志眞 教授

用無母數迴歸方法監控曲線型品質特性之製程

Profile Monitoring by Nonparametric Regression

研 究 生：翁崢珮 Student：Zheng-Pei Weng

指導教授：洪志眞 教授 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 Master

in Statistics June 2004

Hsinchu, Taiwan, Republic of China

用無母數迴歸方法監控曲線型品質特性之製程

研究生：翁崢珮 指導教授：洪志眞 教授

國立交通大學統計學研究所

摘 要

此篇論文的重點在，對品質特性為一個函數或一條曲線的製程，

我們提出一些新的控制圖來做監控；除了平滑的假設外，在此不限定

曲線的形式。

我 們 提 出 了 三 個 新 的 監 控 方 法 ： (1) 用 無 母 數 迴 歸 方 法

“B-spline”對曲線做配適，再利用T

控制圖對製程做監控；(2)用典

型的控制圖EWMA、EWMSD或R控制圖對樣本曲線和參考曲線間的殘差做

監控；(3)用事先制定的測度(metrics)及其臨界值對品質特性曲線做

監控。

經過多次的模擬，利用平均連串長度(ARL)來比較所有方法，結

果顯示在曲線型品質特性為指數函數時，我們所提出的新方法都表現

得相當良好。最後在面對製程的各種偏移情形，該使用何種控制圖做

監控，我們也提出了一些建議。

Profile Monitoring by Nonparametric Regression

Student：Zheng-Pei

Weng

Advisor：Dr. Jyh-Jen Horng Shiau

Institute of Statistics

National Chiao Tung University

ABSTRACT

In this paper, we propose and study some control chart methods for monitoring processes for which the quality is characterized by a profile or a function. No assumptions are made on the form of the function except that it represents a smooth curve.

Three approaches of monitoring schemes are proposed: (1) Use the “B-spline” smoothing method to fit the profiles of processes and design a T2 chart to monitor the deviations of the sample profiles from the in-control reference profile. (2) Use typical control charts – the EWMA, or EWMSD, or R chart to monitor the residuals between the sample profiles and the in-control reference profile. (3) Use some metrics defined to measure the deviations of the sample profiles from the in-control reference profile to monitor the sample profiles.

We construct a simulation study using an exponential profile as an example to investigate the effectiveness of the proposed schemes. The performances of these schemes are evaluated and compared in terms of the average run length. The results of the simulation study show that all approaches appear to perform well for exponential profiles. It is found that some control schemes are preferable to others for particular types of process changes.

誌 謝

在統研所的這兩年，除了學習到書本上的知識和一些電腦的專業

技能外，也交到了許多好朋友和球友，讓熱愛活動的我過的很充實，

很快樂。其中我特別感謝我的指導老師 洪志眞教授不厭其煩的指

導，博士班的碩慧學姊細心的指教，使我的論文能夠順利的完成，也

很謝謝我的研究所同學葉小青、寶文、翠英、Ken、Ming、超毅、志

浩、政輝、巧慧、貞佑、忠庭等等，博士班的牛哥和泰賓學長，他們

不僅帶給我生活上的歡樂，在學業上的互相切磋討論也讓我獲益良

多。

此外還要感謝家人給我的支持和鼓勵，老同學的勉勵和加油，最

後將此論文獻給我的家人、朋友們，和有緣相聚的統計所同學及學弟

妹。

翁 崢 珮 謹誌于

國立交通大學統計學研究所

中華民國九十三年六月

Contents

Chinese Abstract --- i

English Abstract --- ii

Acknowledgment --- iii

Contents --- iv

1. Introduction --- 1

2. Literature Review --- 3

3. Methodologies --- 9

3.1 B-splines

---

9

3.2 Residual EWMA and R Charts

---

12

3.3 Metrics ---

---

13

3.4 EWMSD Chart ----

---

15

4. Simulation Studies --- 16

5. Conclusions --- 22

6. References --- 23

1. Introduction

Statistical process control (SPC) has been successfully applied in many domains, especially in industries. We often use control charts to monitor processes in most SPC applications. In the past, it was always assumed that the quality of a product or process can be measured by one or multiple quality characteristics. However, in many practical situations, the quality of a product or process is better characterized by a profile or a function. One example involved aspartame (an artificial sweetener), which is characterized by the amount that dissolves per liter of water at different levels of temperature (Kang and Albin, 2000). For illustration purpose, Figure 1 gives an example of 4 hypothetical sample profiles created based on the Figure 1 of Kang and Albin (2000). Another example is a semiconductor manufacturing application involving the calibration of a mass flow controller in which the performance of the process is characterized by a linear function (Mestek et al., 1994). Some profiles monitoring methods are proposed for this type of processes. Walker and Wright (2002) proposed additive models to assess the sources of variation of the density profile of particleboards. Jin and Shi (2001) used wavelets to monitor and diagnose process faults. Kang and Albin (2000) and Kim, Mahmoud, and Woodall (2003) proposed some methods for monitoring the process in which the performance is characterized by a linear profile.

Kang and Albin (2000) presented two approaches to monitor linear profiles. The first approach uses a multivariate T2 control chart to monitor the profile parameters, slope and intercept, simultaneously. The second approach treats the residuals of the sample profile (from the reference profile) as a subgroup and uses a combined EWMA/R (exponentially weighted moving average/range) chart for profile monitoring. Kim et al. (2003) presented another approach. First, for each profile, code

the X-values so that the X is zero and the estimators of the intercept and the slope of the regression line are independent. Then, construct two two-sided EWMA charts to monitor the intercept and the slope separately, and use a one-sided EWMA chart to monitor the process variation. This scheme is called EWMA3. We remark that these

two papers only addressed liner profiles.

The main objective of this paper is to propose and study profile monitoring schemes for profiles of more flexible shapes such that the schemes can be applied to more general and practical situations. To model a profile with no restriction on the functional form, it is natural to consider the nonparametric regression model in which the function is only assumed to be smooth. In this paper, for the profile modeling, we adopt the spline regression as the curve fitting/smoothing technique for the simplicity and readiness of a direct extension from the simple linear regression used by Kang and Albin (2000) and Kim et al. (2003) in constructing the T2 chart. The second approach we consider is a monitoring scheme based on the residuals of sample profiles from the reference profile, similar to the combined EWMA/R scheme given in Kang and Albin (2000) and the combined EWMA3 scheme by Kim et al. (2003).

But instead of combining several control charts, we consider the EWMA chart for detecting the mean shift, the R chart for variation change, and the exponentially weighted moving standard deviation (EWMSD) chart for the variation increase separately. Since each of these control charts is designed for a particular type of process changes, their detecting powers usually are better than those more general-purposed charts. The 3rd approach is to use some metrics defined to measure the deviations between observed profiles and the reference profile to detect process changes. This method was firstly presented in Gardner et al. (1997), in which spatial signature metrics were used to diagnose the equipment faults. They reported that these metric-related charts are very powerful in detecting standard deviation shifts. Finally,

a simulation study is conducted to compare the performances of all the proposed methods in terms of the average run length (ARL).

The remaining of this paper is organized as follows. Section 2 reviews the monitoring schemes proposed by Kang and Albin (2000) and Kim et al. (2003). Section 3 describes the three proposed approaches in details. Section 4 presents the results of a simulation study with some discussions and recommendations. Finally, Section 5 concludes the paper with a brief summary and discussion.

2. Literature Review

Let {( ,x yi ij),i=1, 2,...,n} be the observations of the j

th

profile, . Kang and Albin (2000) and Kim et al. (2003) assumed that, when the process is in statistical control, the underlying model is

1, 2,...

j=

0 1 , 1, 2,..., ij i ij

y = A +A x +ε i= n, (1) where theεij's are independent and identically distributed (i.i.d.) normal random

variables with mean zero and variance 2 .

0

σ >

In establishing a monitoring scheme, there are usually two phases: Phase I and Phase II. The interests in Phase I are on understanding the process variation and estimating the in-control process parameters from historical data in order to set up appropriate control limits. In Phase II, we then use the control limits established from Phase I to monitor the on-line process data successively. The multivariate T2 and the residual control charts designed for Phase II by the two sets of authors mentioned above are described below.

Kang and Albin (2000) proposed two control strategies to monitor the process. The first strategy is to use the bivariate T2 chart to examine the regression coefficients. Assume that the in-control values of the parametersA₀, A₁, andσ in Equation (1) 2

are known. The least-squares estimators of A₀ and A₁ for the jth sample profile are, respectively, 0j j 1j a = y −a x and _1j xy j( ) xx S a S = , (2) where 1 1 n j _i ij y =n−

∑

∑

∑

∑

. The estimators and jointly follow the bivariate normal distribution with the mean vector

0 j a a1 j 0 1 (A ,A) = µ (3) and the variance-covariance matrix

2 0 01 2 01 1 , σ σ σ σ ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ Σ (4) where 2 2 2 0 2 2 1 1 ( ) 1 , xx xx x n S S σ σ σ σ = + = , and 2 01 . xx x S σ = −σ

It can be seen that and a are dependent. We express the sample slope and the sample intercept in Equation (2) as the vector '

0j a _1j 0 1 ( , ) j = a j a j z )

. Then, the T2 statistic for the jth sample is given by

(5)

2 1

( ) ' (

j j j

T = z −µ Σ z− −µ ,

which follows the chi-square distribution with 2 degrees of freedom. Thus, the upper control limit is UCL =χ_2,2_α, whereχ_2,2_α is the 100(1−α) percentile of the chi-square distribution with 2 degrees of freedom. When there are shifts from the nominal values,

0

A and A₁, Kang and Albin (2000) pointed out that T in Equation (5) follows the _j2

non-central chi-square distribution with the parameter of noncentrality

(6)

2 2

( x) n Sxx

τ = λ β+ +β

when the intercept A₀ shifts to A₀+λσ and the slope A₁ shifts to A₁+βσ . One disadvantage of the T2 chart is that it cannot distinguish the sources of the process change --- mean or variance. Kang and Albin (2000) then proposed another approach to monitor the mean and variance separately with two charts. The second strategy is to use an EWMA chart to monitor the average of the regression residuals and simultaneously an R chart to monitor the range of the residuals. One reason for adding the R chart to the EWMA chart is that the EWMA chart is not sensitive to shifts in the process variation. Another reason is that the EWMA chart based on the residual average is not sensitive to some shifts in and , for which the magnitudes of the residuals tend to be large, but the residual average tends to be very small. Details of these two charts are reviewed below.

0

A A₁

The regression residuals for the jth sample are

(7)

0 1 , 1, 2,..., . ij ij i

e = y −A −A x i= n

The average of the residuals for sample j is

1 . n ij i j e e n = =

∑

Then, the EWMA statistic for sample j is given by

1

(1 ) ,

j j j z =θe + −θ z ₋

where θ ( 0< ≤ ) is a smoothing constant and θ 1 z0 = . The control limits for the 0

(2 ) LCL L n θ σ θ = − − and UCL L (2 )n θ σ θ = − , (8)

where L is a constant selected for the control chart to achieve a prescribed in-control ARL. If is smaller than LCL or greater than UCL, it is regarded as an out-of-control signal. For the R chart, the R statistic for the j

j z th sample is given by max ( ) min ( ). j i ij i ij R = e − e

The control limits are

2

( _R )

LCL=σ d −L d₃ and UCL=σ(d₂+L d_{R 3}), (9) where LR is also a positive constant chosen to achieve a specified false alarm rate, and are constants relating the range and the standard deviation, which depend on the subgroup size n. Commonly used values of and can be found in Montgomery (2001). The process is claimed out of control when any of the EWMA and the R statistics is out of control.

2 d 3 d 2 d d₃

Kang and Albin (2000) used simulations to demonstrate the effectiveness of these two approaches. They reported that the out-of-control ARL of the combined EWMA/R chart for the residuals is smaller than that of the T2 chart for shifts in the intercept and the slope. On the other hand, the multivariate T2 chart requires only one chart.

Kim et al. (2003) first coded the X-values by centering so that the average of the coded X-values is zero. Then, the least-squares estimators of the slope and intercept are independent. See Myers (1990) and Ryan (1997). Accordingly, transform the model in Equation (1) into

(10)

*

0 1 , 1, 2,..., , ij i ij

y =B +B x +ε i= n

estimator of B is ₀ b_{0 j} = yj . The estimator of B is the same as that of1 .

Because and are independent, they are normally distributed with means

1

A

0 j

b b1 j B , 0

1

B and variances , , respectively. Consequently, the authors

recommended three EWMA charts to monitor the slope, the intercept, and the error variance, respectively. The EWMA statistic for the intercept parameter

2 / n σ 2 /S_xx σ 0 B is (11) 0 ( ) (1 ) ( 1) I j I EWMA j =θb + −θ EWMA j−

with θ ( 0< ≤ ) being a smoothing constant and θ 1 EWMAI(0)=B0. The control

limits are 0 (2 ) I LCL B L n θ σ θ = − − and 0 (2 ) I UCL B L n θ σ θ = + − , (12)

where is a constant chosen for achieving a prescribed in-control ARL. When or , an out-of-control signal is alarmed. The EWMA statistic for the slope parameter

I

L

( )

I

EWMA j <LCL EWMA j_I( )>UCL

1 B is (13) 1 ( ) (1 ) ( 1) S j S EWMA j =θb + −θ EWMA j−

with θ ( 0< ≤ ) being a smoothing constant and θ 1 EWMAS(0)= . The lower and B1

upper control limits are, respectively,

1 (2 ) S xx LCL B L S θ σ θ = − − and 1 (2 ) S xx UCL B L S θ σ θ = + − , (14) where is a constant chosen in order to obtain a specified in-control ARL. They also used a one-sided EWMA chart to detect only increases in process variability. The EWMA statistic for that is given by

S L 2 0 ( ) max{ ln( ) (1 ) ( 1) , ln( )} , E j E

EWMA j = θ MSE + −θ EWMA j− σ (15)

where θ ( 0< ≤ ) again is a smoothing constant,θ 1 2 0

(0) ln( )

E

1 * 0 1 1 ( 2) n ( ) j _i ij j 2 j i

MSE = −n −

∑

Then an upper control limit is given by

[ln( )] , (2 )

E j

UCL L θ Var MSE

θ

=

− (16)

where is again a constant selected to give a specified in-control ARL. The process is claimed out of control if any of the three statistics given in (11), (13), and (15), respectively is beyond the control limits.

E

L

Kim et al. (2003) showed by a simulation study that their methods are generally more effective and seem much more interpretable than the methods of Kang and Albin (2000) in Phase II for detecting sustained shifts in either the intercept or the slope or increases in the error variance.

In practice, the in-control values of parameters are usually unknown. Common practice is to use some historical data to estimate these parameters and construct preliminary control charts during Phase I. Kang and Albin (2000) estimated the regression parameters of Equation (1) based on the historical data containing k sample profiles. For each profile, obtain estimates and a of and , respectively. Then the reference line can be estimated by and

0 j a 1 j A0 A1 0 1 ˆ Y =a +a X σˆ2 =MSE, where 0 1 0 k j j a a k = =

∑

∑

∑

2 1 0 ( ) ' ( 1 j j j k T k − = − − ) , − z z S z z (17)

where z=(a₀ ,a₁)' and S is an unbiased estimate of Σ with components

2 11 1 ( ) xx x S MSE n S = + , ₂₂ , xx MSE S S = and ₁₂ ( ). xx x S MSE S = −

Kang and Albin (2000) proved that the T2/2 statistic in Equation (17) has the F distribution. Thus, the upper control limit of the T2 chart in Phase I is

2,( 2) ,

2 _{n k}

UCL= F _{− α}. For EWMA and R charts in Phase I, the control limits are modified by substituting σ with MSE . There are other methods to examine the historical data in Phase I. Stover and Brill (1998) proposed a Hotelling’s T2 approach that is similar to the T2 method of Kang and Albin (2000). The distinction between the two is the estimate of the variance-covariance matrix. Another approach is to use a univariate chart based on the first principal component corresponding to the vectors containing the estimators of the intercept and slope. But Kim et al. (2003) advised against this principal component method since it is unable to detect the shifts in the direction perpendicular to the first principal component.

3. Methodologies

3.1 B-splines

In order to extend linear profiles to any smooth functions, a smooth curve fitting technique is needed for profile smoothing from noisy data. A polynomial function is not very flexible for approximating curves with different degrees of smoothness at different locations. One way to overcome this drawback is to use locally a polynomial approximation of low degree. Another way is to allow the derivatives of the approximating function to have discontinuities at certain locations. This can be accomplished by fitting piecewise polynomials or splines. Frequently, a cubic spline, i.e., a piecewise polynomial with continuous first two derivatives, is used for such

approximation. Consider the nonparametric regression model ( ) , 1,..., ,

i i i

y =m x +ε i= n (18)

where m(x) is a regression curve and ε_i are i.i.d. normally distributed with zero mean and common variance . In this paper, we adopt the B-spline regression method for its popularity and simplicity. See de Boor (1978) for the definition of the B-spline basis. The points where the derivatives of the approximation function could have discontinuities are called knots. Let

2

0

σ >

1,..., b k

t t₊ be the knots, where b is the number of bases and each polynomial is of order k. Then each basis of order k is a k-2 continuously differentiable function. Define B_{l k,} as a B-spline basis that is nonzero only on the interval ( ,t t_{l l k+} ),l=1, 2,..., .b A B-spline curve of order k can be constructed as , 1 ( )_i b _{l l k}( )_{i i} [ , ], 1,..., , l P x =

∑

where ’s are the unknown B-spline coefficients and n is the number of set points in the interval [u, v]. So we can modify Equation (18) to

l c , 1 ( ) , 1,..., . b i _l l l k i i y =

∑

The B-spline basis B_{l k,} of order k can be defined as

1 ,1 1 , , 1 1 1 1 , 1 , ( ) 0 ; 2 , ( ) ( ) ( ). l l l l l l l k l k l k l k l k l l k l for t t t if k B t for t t and t t t t t t if k B t B t B t t t t t + + + − + + − + + ⎧ ⎧ ≤ < = _{= ⎨} ⎪ _{< ≥} ⎪ ⎩ ⎨ − − ⎪ _{≥ = +} ⎪ _{− −} ⎩ 1, −1

The following S-PLUS function “spline.des” generates a matrix of B-spline bases evaluated at points provided by users.

spline.des(knots, x, ord, derives),

where “knots” is the vector of knots for the spline, “x” represents x-coordinates at which to evaluate the spline basis functions, “ord” is the order of the spline, and

“derives” is the order of the derivative to evaluate at each of the points. If this vector is given, it must have the same length as x. The default is a vector of zeros of the same length as x. We show an example of B-spline bases in Figure 2.

After constructing B-spline bases, we use these normalized B-spline functions

, l k

B as the regression basis function. For the given knots, the spline regression method finds the best spline approximation via the following least squares regression:

2 , 1 1 { ( n b i l l k i i l min y c B x = = −

∑

∑

ˆ=( -1 '

c B'B ) B y

wherey=(y₁,...,y_n)' , cˆ=(c1,...,cb) ', and B is the design matrix with the ( i , l )th element B_{l k,} ( ),x_i l=1,...,b , i=1,...,n . Then c has a multivariate normal distribution with the mean vector c and the variance-covariance matrix

ˆ 2 1 ( ) σ − = Σ B'B .

In phase II, it is assumed that the reference profile is known. Denote it by . We now establish a B-spline representation of the reference profile. First obtain n pairs of data,

( )

f x

{( , ( )),x f x_{i i} i=1,..., }n . Let be the least squares solution of (20) with replaced by

c

i

y f x . We then treat this vector c as the true in-control value. For the ( )_i

jth sample profile, compute the vector of sample estimators cj =(c1j,...,cb j) ', where 1j,..., b j

c c are the estimated sample B-spline coefficients. Then the T2 statistic of the

jth sample is given by . But the upper control limit in phase II needs to be modified as UC 2 1 ( j ) ' ( j j T = c −c Σ c−  −c) L=χ_b2_,α 1

. And in Phase I, we also use historical data containing k sample profiles to construct the control chart. Then the reference curve can be estimated by y B cˆ =  with c=( ,...,c c_b) ' and σˆ2 =MSE where

1 1 1 1 ,..., , k k j b j j j b c c c c k k = = =

∑

∑

∑

with MSEj = yj−B cˆj 2/(n b− ) and ‧ denotes the Euclidean norm. The T

2

statistic for the jth sample profile is then modified by

2 1 0 ( ) ' ( 1 j j j k T k − = − − ) , − c c S c  c (22)

where is an unbiased estimate of . Thus, the upper control limit of the T 2 ˆ ( ) σ − = S B'B 1 _Σ 2

chart in Phase I is UCL=bF_{b n b k,( − ) ,α}.

The choice of the number of bases, as in the role of the smoothing parameter in any nonparametric regression methods, is an important issue. The boundary effect is another issue. We will address these two issues in Section 4.

3.2 Residual EWMA and R Charts

The second approach we propose is to use the EWMA and the R chart to monitor the residual average and the range, respectively, instead of the combined EWMA/R chart proposed by Kang and Albin (2000). Kang and Albin (2000) found that the value

L = 3.1151 yields an in-control ARL of approximately 802 for the EWMA chart and

an in-control ARL of approximately 261 for the R-chart. Although ARL0 of their

combined procedure is close to 200, the detecting power of the EWMA chart is unequal to that of the R chart. In order to increase the detecting powers of both the EWMA and the R charts, we decide to monitor the process by using the EWMA chart to detect mean shifts and the R chart to detect standard deviation shifts of the residuals. The regression residual vector for the jth sample is

ˆ { '} .

j = j− j = −

-1

e y y I B(B'B) B y_j

3.3 Metrics

The third approach is to use some metrics to monitor the process. Gardner et al. (1997) presented a new methodology based on some “metrics” for the equipment fault detection. The equipment fault detection was not only in that it incorporates the use of integrated spatial information in a virtual wafer surface that was fitted by thin-plate splines, but also in that it can be used to detect and classify equipment faults at the same time. Their main focus was to use the differences between the observed and the expected virtual wafer surfaces to construct metrics which can be used to detect and diagnose various types of equipment faults. Thus, the following general and specific metrics were designed:

2 1 ( ) R M =

∫

∫

∫

∫

∫

∫

∫

where g is a newly fitted thin-plate spline surface, T is the target surface, and R denotes the wafer surface region. Metrics 1-3 are general metrics used to detect the presence of an equipment fault, and Metrics 4-5 are specific metrics used to detect specific fault patterns. Gardner et al. (1997) proposed an alternate Bayesian (simulation) approach that can be taken to determine the null distribution of the

metrics. According to a procedure given in Green and Silverman (1994), assuming a Gaussian prior distribution, the posterior distribution of the thin-plate spline surface g has the following multivariate normal distribution:

2 _ˆ

ˆ ˆ

~ MVN [ ,σ ( )]λ

g g A ,

where gˆ is the vector of fitted values, 2

ˆ

σ is calculated as the residual sums of squares about the fitted curve divided by an effective degrees of freedom, and A( )λˆ

is the projection matrix which maps the vector of observed values to their predicted values. Simulate M independent sets of n observations from the posterior multivariate normal distribution described above. A spline surface is fitted to each set of observations, and the metrics are calculated. As a result, M independent values are obtained from the null distribution of each metric. The (100×α)th percentile of the M simulated metric values is the critical value for determining if the new curve is far from the target curve. For each newly observed curve, metrics are calculated and compared to their corresponding critical values in order to determine whether or not a fault is detected. In the experimental data, they used the average predicted surface of two wafers as the target surface. 5,000 observations were simulated using parametric bootstrapping from _{ˆ ˆ}2 _ˆ

MVN[ ,g σ A( )]λ where is the average of the predicted values using the thin-plate spline fitting for these two wafers, and

ˆ

g

2 _ˆ

ˆ ( )

σ A λ is the average of the covariance matrices from the thin-plate spline fitting for these two wafers. Then these simulated observations were used to obtain 5,000 metric values and the corresponding empirical α=0.01 “critical value” for each metric.

With the same idea, we select two different and popular metrics to monitor the process in all types of shifts. They are Squared metric and Absolute-value metric. Two metrics are defined in the following:

2 1 1 { ( ) ( )} n i i i M g x T x = =

∑

∑

where g is a newly fitted B-spline curve and T is the reference profile. To construct the control charts, we need to find the critical value of the null distribution of each metric.

Since the distributions of these metrics are difficult to obtain, we use historical data (k sample profiles) to compute the reference profile in Phase I. Fit a B-spline to each of the k sample profiles. Denote the average of k regression estimators for the lth B-spline coefficient by cˆ_l . So the estimated reference profile , and

1 ˆ ( ) b l l k l c B x =

∑

σ = are available by Equation (21). Simulate M sets of n observations from the following model:

, 1ˆ ( ) , 1,..., , b i _l l l k i i y =

∑

ε σ . For each metric, compute the metric value for each of the simulated profiles. Let the (100×α)th percentile of these M values of the metric be the critical value. The process is claimed out-of-control when the metric of the newly observed profile is greater than the critical value.

3.4 EWMSD Chart

The last approach is using exponentially weighted moving standard deviation (EWMSD) to monitor the process. It is very sensitive in detecting shifts in process variation, particularly when shifts are relatively small. Define the sample standard deviation of the residuals for the jth sample as

2 ˆ , 1, 2,.... j j j s j n b − = = − y y

The EWMSD statistic is given by

1

(1 )

j j

v =θs + −θ v_j− ,

where θ ( 0< ≤ ) is a smoothing constant and θ 1 v0 =σ0. In Phase II, assume the

in-control value of σ is known. The control limits for the EWMSD chart are 2

(2 ) LCL L n θ σ σ θ = − − and UCL L (2 )n θ σ σ θ = + − ,

where L is also a constant chosen to give a specified in-control ARL. When the EWMSD statistic is beyond the control limits, the process is claimed out-of-control. In Phase I, the control limits are modified by substituting σ with MSE .

4. Simulation Studies

We assess the performances of these approaches in terms of ARL through simulation studies. For all approaches except the metric method, we assume the underlying reference profile is known. Denote the in-control ARL value by ARL0. All

charts are designed to have the same ARL0 = 200, which corresponds to α=0.005.

For each metric under study, we simulate 50,000 in-control profiles to approximate the distribution of the metric. The critical value is set to be the 99.5th percentile of these simulated metric values such that the ARL0 = 200. The smoothing constant θ

may affect the ARL performances of the EWMA and EWMSD charts. Details can be seen in Lucas and Saccucci (1990). In our study, the smoothing constant is set to 0.2. In Phase II process monitoring, for each control chart, run lengths are generated for the in-control and various out-of-control situations. The run length is the number of profiles generated when the first out-of-control signal occurs. Generate profiles from the underlying model until the first out-of-control signal is alarmed to obtain a run length. Repeat N times to obtain N run lengths and then estimate the ARL by averaging these N run lengths. Denote the estimator byARLˆ . To obtain the standard

error of (the ARL estimator), one simple way is to compute the sample standard deviation of the N simulated run lengths and then divided it by

ˆ ARL

N . Another

estimator of the standard error is

1/ 2 1 ˆ ARL(1 ) ˆ ARL N − ,

since the run length has a geometric distribution with mean 1/ p and standard deviation (1 -p)1/2 / p with p = 1 / ARL. In our simulation, we generate a total of N = 100,000 run lengths to estimate the ARL value. If we construct the control charts such that ARLˆ ₀ is about 200, then the standard error of this estimator is about

1/ 2

1

200(1 ) / 100, 000 0.631 200

− ≈ . We adopt the first approach in this study.

In this paper, we extend the linear profile to a functional-form-free smooth function to adapt to more general cases. To choose an underlying reference profile for our study, we mimic the aspartame curves in Figure 1. As a result, we take an exponential function as the reference profile. In our simulation study, consider the exponential profile of the form 0( 1)2

0 0 N x

Y = +I M e− − +ε , where ε ~N(0,1). The in-control reference profile is 1 15+ e− −(x1)2, which is displayed in Fiugre 3. We choose

i

x -values of 0, 0.08,…, 3.92 (n = 50) in our simulation. Four different types of shifts

are considered in the simulation study: I shift, M shift, N shift, and error variance increase. The curve varies in different ways with different varying coefficients. Figures 4-6 illustrate the effect of the shifts. In these figures, the solid lines represent the true function , while the dotted lines represent the shifted curves. Figure 4 gives the curve that is 2 units upward by an I shift. The dotted line in Figure 5 is the curve of , that is, the y-value is magnified with the same multiple at each point because of an M shift. Figure 6 shows a case of N shift, and the resulting curve seems narrower than the reference curve. Its function is .

2 ( 1) 1 15 x e− − + 2 ( 1) 1 18 x e− − + 2 1.5( 1) 1 15 x e− − +

For T2 chart, we use the sequence (-1.2, -0.8, -0.4, 0, 0.4, 0.8, 1.2, 1.6, 2, 2.4, 2.8, 3.2, 3.6, 4, 4.4, 4.8, 5.2) as the knots with order 4 so that the number of B-spline bases is 13. The knots sequence is equidistant, and the thirteen cubic B-splines are shown in Figure 7. Approximate the true exponential function by spline regression with these 13 B-spline bases. The fitted curve (dashed line) is shown in Figure 8. The circles are the values on the exponential function and the dashed line goes through these values smoothly. So the B-spline approximation is very close to the true function. The average squared error of the approximation is 0.0000099498 and it is computed by

, 1 1 1 ˆ ( ( ) ( )) n b i l l k i f x l c B i n

∑

∑

It is well known that boundary effect is a potential problem in smoothing methods. We use a simple simulation study for T2 chart to illustrate the boundary effect that we encountered in our study. First, we fix the number of set points at 20, and change the number of B-spline bases to 5, 9, and 13 and find that the ARL0 of T2

chart are 197.976, 193.265, and 187.814, respectively. That is, when n = 20, the ARL0

decreases as the number of B-spline bases increases and the boundary effect is more and more obvious. Because the number of set points 20 is small, the B-spline coefficients for the boundary bases are not accurate enough. It leads to a larger T2 statistic and a smaller ARL0. So we increase the number of set points to 50. Compute

the average of 50,000 replications of each B-spline coefficient. The simulation results are shown in Table 1. The first row is the 13 coefficients of the exponential reference curve fitted by B-splines and these values are treated as the true coefficients. The second row gives the differences between the average of 50,000 replications and the true coefficient for the first simulation. The third and fourth rows give the results of the second simulation and the third simulation, respectively. From Table 1, we can see the variability of each coefficient estimate. The two largest differences occur at the

13th and the first coefficient, respectively, in every simulation. In order to achieve ARL0 = 200, we omit the coefficients 1 and 13 in constructing the T2 statistic.

Naturally, the control limit of the T2 chart is adjusted toχ_11,0.0052 .

Instead of using the EWMA and R charts simultaneously to monitor the process, we use the EWMA chart to detect the mean shifts, say, I, M, or N shifts, and use the R chart to detect the shift in the error variance. For the metric method, use Squared metric and Absolute-value metric to detect all types of shifts. All simulation results are presented below.

Table 2 and Figure 9 give the ARL values and curves for shifts in I0 in unit of

σ , respectively. In Table 2, the smallest ARL values are marked with the deep gray and the second smallest are marked with the light gray. The EWMA chart and two metrics perform better than the T2 chart. The EWMA chart detects all shifts faster than others and the Absolute-value metric is the second best.

Table 3 and Figure 10 show the ARL values and curves for M shifts, respectively. The EWMA chart also performs much better than the T2 chart and two metrics over the entire range of shifts considered. We can see that the performances of two metrics and T2 chart are very similar with the Absolute-value metric better for smaller shifts and T2 chart better for larger shifts.

Table 4 and Figure 11 show the ARL values and curves for N shifts, respectively. Again, the detecting power of the EWMA chart is the largest and the T2 chart is the smaller than the other charts. Secondly, the results of two metrics are similar and the Absolute-value metric is faster for smaller shifts, say, 0.14σ or less, and slower for larger shifts.

Table 5 and Figure 12 give the ARLs for shifts in the process standard deviation. The R chart performs slightly better than the T2 chart and two metrics over the entire

range of shifts. And the Squared metric is also slightly faster than the Absolute-value metric over the whole range of shifts. But the EWMSD chart performs much better than the other charts, because it is designed specially for detecting the standard deviation shifts. This demonstrates that the EWMA and EWMSD charts are more sensitive for all shifts than others, especially for smaller shifts.

Table 6 gives the ARL performances when simultaneous M and N shifts are considered. The T2 chart is uniformly slower than other charts over the entire range of shifts. The EWMA chart performs the best with any combination of M and N shifts, and Absolute-value metric is the second best. Tables 7 and 8 present the simulation results with I and M, I and N combinations, respectively. Two tables also show similar results.

To investigate the effect of the number of set points (n) and the number of B-spline bases (b), another simulation study is conducted. Table 9 and Figure 13 give the ARLs of the T2 chart with different number of bases, 5, 9, and 13, for I shifts. The results show that the detecting power of the T2 chart increases as parameter b decreases. Table 10 and Figure 14 give the ARLs of the T2 chart with different n = 30, 40, and 50 for I shifts. We can see that the ARL decreases as n increases with various sizes in I shifts. In order to improve the power of the T2 chart, it seems that a smaller b and larger n is suggested. However, there is another concern – the squared error between the fitted profile estimated by B-splines and the reference profile. Although the detecting power of T2 chart with n = 50 and b = 5 is larger, the average squared error comes to 0.4829616. From Figure 15 we can see that, with b = 5, the B-spline does not approximate the reference profile well. Thus, b should not be too small.

Table 11 and Figure 16 show that the ARL decreases as n increases for the EWMA chart with various shifts in M shifts. Tables 12-13 and Figures17-18 show that the detecting power increases as n increases for all metrics with I shifts or standard

deviation shifts.

Back to linear profiles, we compare the ARL performances of our approaches to that of the methods proposed by Kang and Albin (2000) and Kim et al. (2003). We simulate again by using the same control limits of each chart. The linear profile model used by Kang and Albin (2000) is y_ij = +3 2x_i+ε_ij , where with fixed . . . ~ (0,1 i i d ij N ε ) i

x -values of 2, 4, 6, and 8. Three different types of shifts, intercept shifts, slope

shifts and standard deviation shifts, are considered in their papers. From our simulation study, we notice that our ARL values of EWMA3 are 1 more than the

values given in Kim et al. (2003). We guess the discrepancy may come from the way of counting the run lengths. They probably counted the run lengths as the number of the profiles before the out-of-control signal occurs, so that all values are one less than ours.

Table 14 and Figure 19 give the ARL values and curves for intercept shifts, respectively. The EWMA chart performs much better than othercharts for smaller shifts. The Absolute-value metric detects the larger shifts, say, 1.4σ or more, faster than the other charts. We can clearly see the performance of EWMA3 is better than

that of the methods of Kang and Albin (2000) for smaller shifts, say 0.4σ or less, and the power of detecting larger shifts is the worst.

Table 15 and Figure 20 show the ARL performances for slope shifts. The EWMA chart is much better that the other charts over the entire range of shifts. Two metrics perform well for larger shifts with Absolute-value metric better. And the performance of EWMA3 is worse than that for intercept shifts. So in this case, we can use the

EWMA chart or EWMA3 chart to detect the smaller shifts and use the EWMA chart or

Absolute-value metric to detect the larger shifts of intercept and slope.

EWMSD chart also performs the best for smaller shifts, say 2.4σ or less and EWMA/R chart and the R chart perform better than others for larger shifts. The detecting power of Squared metric, Absolute-value metric, and T2 chart are similar. Our ARL comparisons show that our methods are slightly effective than the methods of Kang and Albin (2000) and Kim et al. (2003) for shifts in either the intercept or slope or increases in the error variance.

5. Conclusions

The profile monitoring is a very useful and promising area of research. In this thesis, we focus on the non-linear profile monitoring. Three profile monitoring approaches are proposed. And from the simulation studies, all approaches appear to perform well for the exponential profile under study.

Comparisons among these approaches are made. For I, M, and N shifts in the exponential profile, the EWMA chart performs the best over entire range of shifts, while the Absolute-value metric has secondly smaller ARL. And for standard deviation shifts, the EWMSD chart is a good choice to detect, while we also can use the R chart to monitor the process. And for linear profiles, our simulation study shows that new methods proposed by us are slightly effective than the methods of Kang and Albin (2000) and Kim et al. (2003). As a design issue, it is observed that increasing n is helpful in reducing the out-of-control ARL values. How to properly choose the number of bases is worth further studies.

This study extends the framework of statistical process control to more general applications. More statistical methods, models, and ideas are needed to extend the framework to a more complete profile monitoring strategy.

References

Crowder, S. V. and Hamilton, M. D. (1992). “An EWMA for Monitoring a Process Standard Deviation”. Journal of Quality Technology, 24, 12-21.

de Boor, C. (1978). A Practical Guide to Splines. Springer, Berlin.

Fan, J. and Gijbels, I. (1996). Local Polynomial Modelling and Its Applications. Chapman & Hall.

Gardner, M. M., Lu, J. C., Gyurcsik, R. S., Wortman, J. J., Hornung, B. E., Heinisch, H. H., Rying, E. A., Rao, S., Davis, J. C., and Mozumder, P. K. (1997).

“Equipment Fault Detection Using Spatial Signatures”. IEEE Transactions on

Components, Packaging, and Manufacturing Technology – Part C, 20, 295-304.

Gisela, E. M. and Frank, U. (1996). Numerical Algorithms with C. Springer-Verlag, New York, 325-330.

Green, P. J. and Silverman, B. W. (1994). Nonparametric Regression and Generalized

Linear Models: A Roughness Penalty Approach. New York: Chapman.

Jin, J. and Shi, J. (1999). “Feature-Preserving Data Compression of Stamping Tonnage Information Using Wavelets”. Technometrics, 41, 327-339.

Jin, J. and Shi, J. (2001). “Automatic Feature Extraction of Waveform Signals for In-Process Diagnostic Performance Improvement”. Journal of Intelligent

Manufacturing, 12, 257-268.

Jones, M. C. and Rice, J. A. (1992). “Displaying the Important Features of Large Collections of Similar Curves”. The American Statistican, 46, 140-145.

Kang, L. and Albin, S. L. (2000). “On-Line Monitoring When the Process Yields a Linear Profile”. Journal of Quality Technology, 32, 418-426.

Kim, K., Mahmoud, M. A., and Woodall, W. H. (2003). “On The Monitoring of Linear Profiles”. To appear in the Journal of Quality Technology.

Lucas, J. M. and Saccucci, M. S. (1990). “Exponentially Weighted Moving Average Control Schemes: Properties and Enhancements”. Technometrics, 32, 1-29. Mestek, O., Pavlik, J., and Suchánek, M. (1994). “Multivariate Control Charts:

Control Charts for Calibration Curves”. Fresenius’ Journal of Analytical

Chemistry, 350, 344-351.

Montgomery, D. C. (2001). Introduction to Statistical Quality Control. 4th Edition, John Wiley & Sons, New York.

Myers, R. H. (1990). Classical and Modern Regression with Applications. 2nd Edition, PWS-Kent Publishing Company, Boston, MA.

Ryan, T. P. (1997). Modern Regression Methods. John Wiley & Sons, New York. Simonoff, J. S. (1996). Smoothing Methods in Statistics. Springer-Verlag, New York. Stover, F. S. and Brill, R. V. (1998). “Statistical Quality Control Applied to Ion

Chromatography Calibrations”. Journal of Chromatography, 804, 37-43. Walker, E. and Wright, S. P. (2002). “Comparing Curves Using Additive Models”.

0 5 10 15

Figure1. Milligrams of aspartame dissolved per liter of water from four samples.

3 4 5 6 7 8 0. 0 0 .2 0. 4 0 .6 1,4 B B2,4 B3,4 B4,4 B5,4 B6,4 7,4 B 0,4 B

x y 0 1 2 3 4 0 5 10 15

Figure 3. The reference profile m x( ) 1 15= + e− −(x1)2.

2 2 ( 1) ( 1) solid line :1 15 dotted line : 3 15 x x e e − − − − + + x y 0 1 2 3 4 0 5 10 15 20

Figure 4. The reference curve (solid line) and the curve affected by I shift (dotted line).

2 2 ( 1) ( 1) solid line :1 15 dotted line :1 18 x x e e − − − − + + x y 0 1 2 3 4 0 5 10 15 20

Figure 5. The reference curve (solid line) and the curve affected by M shift (dotted line). 2 2 ( 1) 1.5( 1) solid line :1 15 dotted line :1 15 x x e e − − − − + + x y 0 1 2 3 4 0 5 10 15 20

Figure 6. The reference curve (solid line) and the curve affected by N shift (dotted line).

x 0 1 2 3 4 0. 0 0 .2 0. 4 0 .6

Figure 7. Thirteen cubic B-splines bases.

x y 0 1 2 3 4 0 5 10 15

Figure 8. Scatter plot is the reference curve and the dashed line is the B-spline estimate.

Table 1. Simulation results to illustrate boundary effects; Underlines are larger deviations. Basis 1 2 3 4 5 Cc 2.874264 6.157306 11.611126 16.178782 16.183834 Z1-C -0.0437711 0.0042971 0.00457416 0.00045094 -0.0013554 Z2-C 0.06077714 0.00684893 -0.0002403 -0.0109680 0.00645003 Z3-C 0.08942159 -0.0136533 0.01052427 -0.0160640 0.01881415 Basis 6 7 8 9 10 C 11.603308 6.176568 2.767467 1.423282 1.070967 Z1-C -0.0031249 -0.0001978 0.00197017 -0.0013328 -0.0054912 Z2-C -0.0079634 0.01199316 -0.0198097 0.0244964 -0.0223153 Z3-C 0.00019966 -0.0030393 0.00455043 0.01870481 -0.0077123 Basis 11 12 13 U 1.008571 1.000276 1.003893 Z1-U 0.00677231 0.00502447 -0.0936737 Z2-U 0.00303271 0.01045864 -0.1506641 Z3-U -0.0079037 0.02057199 -0.2117598

Table 2. ARL comparisons under I shifts fromI₀ ToI₀+ασ α Chart 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 T2 199.86 (0.622) 180.07 (0.569) 138.44 (0.437) 93.67 (0.294) 58.76 (0.183) 35.40 (0.110) 21.27 (0.066) 13.03 (0.040) 8.20 (0.024) 5.38 (0.015) 3.70 (0.010) EWMA 200.42 (0.613) 47.39 (0.135) 14.71 (0.033) 7.62 (0.014) 5.06 (0.007) 3.79 (0.005) 3.06 (0.003) 2.59 (0.003) 2.27 (0.002) 2.04 (0.002) 1.87 (0.002) Squared metric 200.70 (0.622) 179.76 (0.566) 138.92 (0.435) 94.37 (0.297) 58.64 (0.184) 35.15 (0.109) 20.85 (0.064) 12.63 (0.038) 7.86 (0.023) 5.15 (0.015) 3.52 (0.009) Absolute-value metric 199.49 (0.625) 175.20 (0.561) 134.78 (0.424) 88.79 (0.277) 53.08 (0.167) 30.87 (0.096) 17.76 (0.055) 10.53 (0.032) 6.49 (0.019) 4.23 (0.012) 2.93 (0.008) 0 50 100 150 200 0 0.1 0.2 0.3 0.4 0.5

ARL

Absolute value metric

Figure 9. ARL comparisons under I shifts of sizeα. α

Table 3. ARL comparisons under M shifts fromM₀ ToM₀+βσ β Chart 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 T2 199.15 (0.623) 175.79 (0.560) 124.60 (0.391) 77.34 (0.243) 43.92 (0.137) 24.32 (0.075) 13.76 (0.042) 8.05 (0.024) 5.00 (0.014) 3.32 (0.009) 2.35 (0.006) EWMA 199.26 (0.613) 63.21 (0.184) 20.63 (0.051) 10.36 (0.021) 6.65 (0.011) 4.86 (0.007) 3.84 (0.005) 3.20 (0.004) 2.76 (0.003) 2.45 (0.002) 2.21 (0.002) Squared metric 201.30 (0.623) 176.35 (0.554) 128.85 (0.405) 81.84 (0.257) 47.02 (0.147) 26.47 (0.082) 14.93 (0.046) 8.79 (0.026) 5.39 (0.015) 3.55 (0.010) 2.49 (0.006) Absolute-value metric 199.37 (0.620) 173.90 (0.548) 122.30 (0.385) 74.54 (0.233) 42.39 (0.132) 23.79 (0.074) 13.56 (0.041) 8.14 (0.024) 5.10 (0.014) 3.43 (0.009) 2.45 (0.006) 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1

AR

L

Absolute value metric

β

Table 4. ARL comparisons under N shifts fromN₀ ToN₀+γσ γ Chart 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 T2 199.4 (0.628) 173.86 (0.547) 120.61 (0.380) 73.36 (0.231) 41.99 (0.131) 23.82 (0.073) 13.80 (0.042) 8.38 (0.025) 5.37 (0.015) 3.67 (0.010) 2.66 (0.007) EWMA 201.3 (0.608) 50.39 (0.144) 16.21 (0.038) 8.44 (0.016) 5.61 (0.009) 4.23 (0.006) 3.43 (0.004) 2.91 (0.003) 2.55 (0.002) 2.30 (0.002) 2.11 (0.002) Squared metric 201.05 (0.638) 171.84 (0.542) 117.68 (0.369) 69.62 (0.219) 38.53 (0.121) 21.16 (0.065) 11.94 (0.036) 7.17 (0.021) 4.58 (0.013) 3.12 (0.008) 2.28 (0.005) Absolute-value metric 199.35 (0.621) 169.74 (0.537) 114.29 (0.361) 66.73 (0.210) 36.79 (0.115) 20.34 (0.063) 11.67 (0.035) 7.14 (0.021) 4.63 (0.013) 3.19 (0.008) 2.36 (0.006) 0 50 100 150 200 0 0.04 0.08 0.12 0.16 0.2

AR

L

Figure 11. ARL comparisons under N shifts of sizeγ . γ

Table 5. ARL comparisons under standard deviation shifts fromσ Toλσ λ Chart 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 T2 199.71 (0.624) 85.74 (0.270) 42.31 (0.132) 23.53 (0.073) 14.42 (0.044) 9.58 (0.029) 6.76 (0.020) 5.05 (0.014) 3.95 (0.011) 3.22 (0.008) 2.70 (0.007) R 200.69 (0.634) 79.84 (0.253) 36.12 (0.113) 18.50 (0.057) 10.51 (0.032) 6.58 (0.019) 4.45 (0.012) 3.22 (0.008) 2.49 (0.006) 2.01 (0.005) 1.69 (0.003) Squared metric 200.23 (0.632) 81.03 (0.255) 38.73 (0.121) 21.03 (0.065) 12.58 (0.038) 8.26 (0.025) 5.80 (0.017) 4.33 (0.012) 3.41 (0.009) 2.78 (0.007) 2.34 (0.006) Absolute-value metric 199.62 (0.628) 85.42 (0.269) 42.48 (0.133) 23.68 (0.073) 14.39 (0.044) 9.53 (0.029) 6.70 (0.020) 5.03 (0.014) 3.91 (0.011) 3.17 (0.008) 2.65 (0.007) EWMV 200.08 (0.618) 30.69 (0.084) 9.14 (0.018) 5.04 (0.008) 3.50 (0.005) 2.73 (0.003) 2.27 (0.003) 1.95 (0.002) 1.73 (0.002) 1.54 (0.002) 1.40 (0.002) 0 50 100 150 200 1 1.1 1.2 1.3 1.4 1.5

ARL

Figure 12. ARL comparisons under standard deviation shifts of sizeλ λ

Table 6. ARL comparisons under combinations of M and N shifts β T2 EWMA Absolute-value metric 0.1 0.2 0.3 0.4 0.5 132.527 (0.417) 86.171 (0.272) 50.568 (0.158) 28.459 (0.088) 16.038 (0.049) 17.461 (0.042) 9.290 (0.018) 6.166 (0.010) 4.591 (0.006) 3.682 (0.004) -0.02 127.463 (0.401) 82.180 (0.259) 46.314 (0.146) 26.716 (0.083) 14.851 (0.045) 81.769 (0.257) 51.071 (0.160) 29.896 (0.093) 17.262 (0.053) 10.149 (0.030) 8.366 (0.016) 5.720 (0.009) 4.331 (0.006) 3.497 (0.004) 2.969 (0.003) -0.04 75.389 (0.236) 45.172 (0.141) 26.194 (0.081) 15.078 (0.045) 8.938 (0.027) 44.251 (0.138) 27.757 (0.086) 16.727 (0.051) 10.175 (0.031) 6.299 (0.018) 5.294 (0.008) 4.086 (0.005) 3.343 (0.004) 2.849 (0.003) 2.504 (0.002) -0.06 39.093 (0.124) 23.681 (0.074) 14.163 (0.043) 8.560 (0.026) 5.328 (0.015) 22.900 (0.070) 14.721 (0.045) 9.376 (0.028) 5.992 (0.017) 4.021 (0.011) 3.856 (0.005) 3.187 (0.004) 2.745 (0.003) 2.421 (0.002) 2.185 (0.002) -0.08 18.998 (0.059) 12.282 (0.037) 7.735 (0.023) 5.046 (0.014) 3.379 (0.009) 11.786 (0.036) 7.943 (0.024) 5.402 (0.015) 3.728 (0.010) 2.688 (0.007) 3.050 (0.003) 2.635 (0.003) 2.343 (0.002) 2.126 (0.002) 1.957 (0.002) γ -0.10 10.053 (0.030) 6.615 (0.019) 4.439 (0.012) 3.111 (0.008) 2.300 (0.005)

Table 7. ARL comparisons under combinations of I and M shifts. α T2 EWMA Absolute-value metric 0.05 0.10 0.15 0.20 0.25 138.147 (0.436) 95.299 (0.301) 60.914 (0.190) 36.854 (0.115) 22.197 (0.068) 17.308 (0.041) 8.387 (0.016) 5.375 (0.008) 3.961 (0.005) 3.166 (0.004) 0.1 132.874 (0.421) 90.184 (0.283) 55.410 (0.174) 32.321 (0.100) 19.042 (0.058) 89.390 (0.281) 58.511 (0.184) 36.237 (0.113) 22.241 (0.069) 13.733 (0.042) 9.289 (0.018) 5.740 (0.009) 4.157 (0.005) 3.279 (0.004) 2.738 (0.003) 0.2 85.664 (0.272) 53.523 (0.170) 32.597 (0.101) 18.822 (0.057) 11.481 (0.035) 52.159 (0.163) 33.431 (0.104) 20.896 (0.064) 13.239 (0.040) 8.467 (0.025) 6.163 (0.010) 4.357 (0.006) 3.403 (0.004) 2.817 (0.003) 2.425 (0.002) 0.3 49.722 (0.155) 31.011 (0.096) 19.080 (0.058) 11.463 (0.034) 7.294 (0.021) 29.134 (0.090) 18.945 (0.058) 12.252 (0.037) 8.056 (0.024) 5.420 (0.015) 4.601 (0.006) 3.536 (0.004) 2.903 (0.003) 2.489 (0.002) 2.193 (0.002) 0.4 27.553 (0.085) 17.439 (0.053) 11.081 (0.033) 7.096 (0.020) 4.615 (0.013) 16.393 (0.050) 10.946 (0.033) 7.410 (0.022) 5.129 (0.015) 3.637 (0.010) 3.688 (0.004) 2.994 (0.003) 2.549 (0.002) 2.234 (0.002) 2.020 (0.002) β 0.5 15.675 (0.048) 10.151 (0.030) 6.699 (0.020) 4.560 (0.012) 3.201 (0.008)

Table 8. ARL comparisons under combinations of I and N shifts. α T2 EWMA Absolute-value metric 0.05 0.10 0.15 0.20 0.25 133.834 (0.423) 91.698 (0.287) 57.973 (0.182) 35.111 (0.109) 21.181 (0.065) 15.043 (0.034) 7.713 (0.014) 5.088 (0.007) 3.808 (0.005) 3.073 (0.003) -0.02 128.016 (0.404) 85.045 (0.269) 51.029 (0.161) 30.966 (0.097) 17.402 (0.053) 81.250 (0.225) 51.889 (0.163) 32.281 (0.100) 19.822 (0.061) 12.282 (0.037) 7.702 (0.014) 5.074 (0.007) 3.802 (0.005) 3.062 (0.003) 2.596 (0.003) -0.04 74.088 (0.233) 47.078 (0.145) 27.204 (0.084) 16.285 (0.050) 9.805 (0.029) 43.259 (0.135) 27.714 (0.086) 17.410 (0.053) 11.021 (0.033) 7.214 (0.021) 5.013 (0.007) 3.773 (0.005) 3.045 (0.003) 2.582 (0.003) 2.265 (0.002) -0.06 37.935 (0.119) 23.307 (0.070) 14.090 (0.044) 8.931 (0.027) 5.658 (0.016) 22.217 (0.069) 14.523 (0.044) 9.459 (0.028) 6.309 (0.018) 4.375 (0.012) 3.724 (0.005) 3.018 (0.003) 2.565 (0.002) 2.248 (0.002) 2.027 (0.002) -0.08 18.990 (0.059) 12.084 (0.036) 7.822 (0.023) 5.085 (0.015) 3.490 (0.009) 11.466 (0.035) 7.803 (0.023) 5.389 (0.015) 3.822 (0.010) 2.824 (0.007) 2.967 (0.003) 2.531 (0.002) 2.222 (0.002) 2.010 (0.002) 1.841 (0.002) γ -0.10 9.519 (0.028) 6.361 (0.018) 4.420 (0.012) 3.102 (0.008) 2.331 (0.006)

Table 9. ARL comparisons of T2 chart with different number of bases for I shifts. α Chart 0.0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 b = 13 199.6 (0.626) 180.28 (0.569) 138.93 (0.441) 94.34 (0.297) 58.84 (0.184) 35.13 (0.109) 20.89 (0.064) 12.68 (0.038) 7.88 (0.023) 5.14 (0.015) 3.53 (0.009) b = 9 199.22 (0.626) 176.11 (0.556) 128.00 (0.403) 81.01 (0.225) 47.56 (0.147) 27.31 (0.085) 15.80 (0.048) 9.46 (0.028) 5.93 (0.017) 3.92 (0.011) 2.77 (0.007) b = 5 199.61 (0.629) 168.45 (0.529) 108.79 (0.345) 61.25 (0.193) 33.08 (0.103) 17.98 (0.055) 10.26 (0.031) 6.19 (0.018) 4.01 (0.011) 2.76 (0.007) 2.04 (0.005) 0 50 100 150 200 0 0.1 0.2 0.3 0.4 0.5

ARL

Table 10. ARL comparisons of T2 chart with different n’s for I shifts. α Chart 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 n = 50 199.97 (0.624) 180.01 (0.567) 139.09 (0.438) 94.53 (0.296) 58.39 (0.182) 35.19 (0.110) 20.97 (0.064) 12.62 (0.038) 7.88 (0.023) 5.16 (0.015) 3.51 (0.009) n = 40 199.30 (0.621) 184.07 (0.579) 148.98 (0.470) 107.22 (0.338) 72.16 (0.226) 46.37 (0.145) 29.02 (0.090) 18.24 (0.056) 11.73 (0.035) 7.67 (0.023) 5.26 (0.015) n = 30 199.1 (0.622) 187.94 (0.593) 159.24 (0.503) 124.77 (0.393) 90.39 (0.284) 62.39 (0.196) 42.30 (0.133) 28.28 (0.088) 18.86 (0.058) 12.78 (0.039) 8.88 (0.026) 0 50 100 150 200 0 0.1 0.2 0.3 0.4 0.5

ARL

Figure 14. ARL comparisons of T2 chart with different n’s for I shifts.

n=50 n=40 n=30 α x y 0 1 2 3 4 0 5 10 15

Figure 15. The reference curve (solid line) and the fitted B-spline curve with 5 bases and n=50.

Table 11. ARL comparisons of EWMA chart with different n’s for M shifts. β Chart 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 n = 50 201.13 (0.627) 63.98 (0.186) 20.69 (0.051) 10.43 (0.021) 6.66 (0.011) 4.87 (0.007) 3.86 (0.005) 3.21 (0.004) 2.77 (0.003) 2.45 (0.002) 2.22 (0.002) n = 40 199.86 (0.615) 73.48 (0.236) 24.94 (0.064) 12.45 (0.027) 7.82 (0.014) 5.67 (0.009) 4.43 (0.006) 3.65 (0.004) 3.12 (0.003) 2.74 (0.003) 2.46 (0.002) n = 30 200.23 (0.628) 88.39 (0.264) 32.19 (0.087) 15.88 (0.037) 9.83 (0.019) 6.95 (0.012) 5.35 (0.008) 4.36 (0.006) 3.68 (0.004) 3.21 (0.004) 2.85 (0.003) n =20 201.60 (0.637) 109.19 (0.329) 44.66 (0.125) 22.57 (0.057) 13.69 (0.030) 9.45 (0.018) 7.10 (0.012) 5.69 (0.009) 4.73 (0.007) 4.07 (0.005) 3.58 (0.004) 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1

ARL

Figure 16. ARL comparisons of the EWMA chart with different n’s for M shifts.

n=5 n=4 n=3 n=2 0 0 0 0

β

Table 12. ARL comparisons of Squared metric with different n’s for standard deviation shifts. λ Chart 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 n = 50 200.36 (0.632) 81.03 (0.255) 38.73 (0.121) 21.03 (0.065) 12.58 (0.038) 8.26 (0.025) 5.83 (0.017) 4.35 (0.012) 3.41 (0.009) 2.77 (0.007) 2.32 (0.005) n = 40 200.22 (0.639) 89.38 (0.272) 44.15 (0.142) 26.48 (0.079) 17.32 (0.045) 11.64 (0.033) 8.22 (0.023) 5.58 (0.016) 4.17 (0.011) 3.36 (0.008) 2.71 (0.006) n = 30 200.68 (0.631) 96.10 (0.296) 52.55 (0.169) 32.67 (0.090) 22.17 (0.057) 15.11 (0.039) 10.59 (0.031) 6.63 (0.021) 5.14 (0.015) 3.75 (0.010) 3.24 (0.007) n =20 201.09 (0.643) 106.97 (0.324) 62.29 (0.194) 40.09 (0.112) 30.61 (0.083) 19.98 (0.048) 14.77 (0.040) 9.40 (0.029) 6.49 (0.020) 4.45 (0.013) 4.82 (0.009) 0 50 100 150 200 1 1.1 1.2 1.3 1.4 1.5

ARL

Figure 17. ARL comparisons of Squared metric with different n’s for standard

n=5 n=4 n=3 n=2 0 0 0 0 λ λ deviation shifts.

Table 13. ARL comparisons of Absolute-value metric with different n’s for I shifts α Chart 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 n = 50 200.25 (0.631) 176.64 (0.559) 134.26 (0.423) 88.90 (0.280) 54.20 (0.169) 30.86 (0.097) 17.93 (0.056) 10.59 (0.032) 6.58 (0.019) 4.22 (0.011) 2.96 (0.007) n = 40 200.46 (0.632) 179.69 (0.582) 145.41 (0.457) 103.91 (0.328) 66.33 (0.210) 42.20 (0.133) 25.16 (0.076) 15.90 (0.049) 9.83 (0.029) 6.45 (0.019) 4.29 (0.012) n = 30 201.73 (0.642) 184.97 (0.591) 156.89 (0.494) 123.00 (0.387) 85.88 (0.268) 59.14 (0.186) 38.39 (0.120) 25.31 (0.077) 16.33 (0.050) 11.00 (0.033) 7.45 (0.022) n =20 200.42 (0.632) 190.15 (0.598) 174.07 (0.556) 142.36 (0.444) 110.42 (0.346) 85.34 (0.272) 62.67 (0.197) 42.88 (0.134) 31.71 (0.098) 22.46 (0.069) 15.96 (0.049) 0 50 100 150 200 0 0.1 0.2 0.3 0.4 0.5

ARL

Table 14. ARL comparisons under intercept shifts of a linear profile. λ Chart 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 EWMA/R 198.85 (0.619) 65.58 (0.192) 17.54 (0.040) 8.32 (0.013) 5.33 (0.007) 3.95 (0.004) 3.17 (0.003) 2.67 (0.002) 2.33 (0.002) 2.10 (0.001) 1.94 (0.001) T2 198.44 (0.627) 137.21 (0.433) 64.0 (0.198) 27.96 (0.087) 13.20 (0.040) 6.90 (0.020) 4.00 (0.011) 2.57 (0.006) 1.84 (0.004) 1.45 (0.003) 1.23 (0.002) EWMA3 199.36 (0.614) 60.23 (0.170) 17.28 (0.036) 8.91 (0.013) 6.12 (0.007) 4.82 (0.004) 4.08 (0.003) 3.59 (0.002) 3.27 (0.002) 3.05 (0.001) 2.89 (0.001) EWMA 200.6 (0.616) 39.28 (0.109) 12.02 (0.025) 6.37 (0.010) 4.30 (0.006) 3.28 (0.004) 2.68 (0.003) 2.29 (0.002) 2.03 (0.002) 1.84 (0.002) 1.67 (0.002) Squared metric 201.43 (0.638) 137.40 (0.435) 63.31 (0.199) 27.90 (0.087) 13.17 (0.040) 6.89 (0.020) 3.99 (0.011) 2.58 (0.006) 1.84 (0.004) 1.45 (0.003) 1.23 (0.002) Absolute-value metric 199.92 (0.627) 124.64 (0.395) 52.03 (0.163) 21.91 (0.067) 10.31 (0.031) 5.45 (0.016) 3.25 (0.009) 2.19 (0.005) 1.63 (0.003) 1.33 (0.002) 1.16 (0.001) 0 50 100 150 200 0 0.5 1 1.5 2

ARL

λ