概似比檢定統計量的指數加權移動平均監控一般線性資料

國

立

交

通

大

學

統計學研究所

碩

士

論

文

概似比檢定統計量的指數加權移動平均監控一般線性

資料

An Exponentially Weighted Moving Average Control Chart

Based on Likelihood Ratio Test Statistics for Monitoring

General Linear Profiles

研 究 生：黃偉振

指導教授：陳志榮 博士

概似比檢定統計量的指數加權移動平均監控一般線性

資料

An Exponentially Weighted Moving Average Control Chart

Based on Likelihood Ratio Test Statistics for Monitoring

General Linear Profiles

研 究 生：黃偉振 Student：Wei-Chen Huang

指導教授：陳志榮 博士 Advisor：Dr. Chih-Rung Chen

國 立 交 通 大 學 理 學 院

統 計 學 研 究 所

碩 士 論 文

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

July 2012

概似比檢定統計量的指數加權移動平均監控一般線性

資料

研究生：黃偉振 指導教授：陳志榮 博士

國立交通大學統計學研究所碩士班

摘要

當製程品質可藉由變數間之線性關係來衡量時，本文提出統計製程 管制可用於監控工業實務上。首先介紹一些概似比檢定統計量的性質。 其次，我們提出概似比檢定統計量的指數加權移動平均來監控一般線性 資料。最後，我們用模擬的方法來說明我們提出的這個方法之表現並討 論它的優缺點。 關鍵字：概似比檢定統計量、指數加權移動平均、一般線性資料

An Exponentially Weighted Moving Average Control Chart

Based on Likelihood Ratio Test Statistics for Monitoring

General Linear Profiles

Student: Wei-Chen Huang Advisor: Dr. Chih-Rung Chen

Institute of Statistics National Chiao Tung University

Abstract

When the quality of a process can be characterized by general linear

profiles, a statistical process control scheme that can be used in industrial

practice is proposed in the paper. First, some properties of the likelihood ratio

test statistics are introduced. Next, an exponentially weighted moving average

control chart based on likelihood ratio test statistics for monitoring general

linear profiles is proposed. Finally, the performance of the proposed

methodology is investigated through a simulation study to show its strength

and weakness.

KEY WORDS: Likelihood ratio test statistics, Exponentially weighted moving average, General linear profiles

誌謝

就讀交大兩年的碩士生活一下子就過去了，時光飛逝，卻令我擁有著充實與 美好的記憶，現在的我即將邁入我人生的另一個階段。 首先，我要特別感謝我的指導教授-陳志榮老師，感謝老師總是不厭其煩地 悉心教導我，鉅細靡遺的批改及修改論文，讓我學習到嚴謹的研究態度。老師的 照顧及付出，點滴在心。此外還要感謝口試委員蔡明田老師、洪慧念老師及許文 郁老師，老師們給予的建議，使我的論文更加完善。 再來感謝陪伴我度過這兩年的同學們，由於我腳的傷勢，一年多來行動都不 太方便，但是大家也都很體諒我，並替我著想，能認識你們真的非常開心。去年 暑假由於我要去動手術，沒辦法和你們一同去畢業旅行，我想這會是我這輩子的 遺憾之一。 然後要感謝交大統計所上的老師們，有老師們的教導讓我在兩年學習到不少。 也感謝所辦的郭姐及怡君姐也一直幫助我們處理一些瑣事，更要感謝家人在生活 上的支持以及鼓勵。最後真的很高興也很感謝交大統研對我們所提供的一切。 黃偉振 謹誌于 國立交通大學統計學研究所 中華民國一 O 一年七月

Contents

1 Introduction 1

1.1 Motivation . . . 1

1.2 Literature Review . . . 1

1.3 Outline . . . 6

2 An EWMA Control Chart Based on LR Test Statistics for Monitoring General Linear Profiles 7 2.1 Model . . . 7

2.2 Known and Unconstrained Case . . . 8

2.3 Known and Constrained Case . . . 10

2.4 Unknown and Unconstrained Case . . . 13

2.5 Unknown and Constrained Case . . . 15

2.6 Proposed Monitoring Scheme . . . 16

3 A Simulation Study 18

5 Future Work 24

List of Tables

5.1 ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in β0 with θj = (β0+ δ0σ, β1, σ)T and τj2 = 4 δ20 for j ≥ 1 . . . 27

5.2 ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in β1 with θj = (β0, β1+ δ1σ, σ)T and τj2 = 120 δ12 for j ≥ 1 . . . . 27

5.3 ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in σ with θj = (β0, β1, δσ)T for j ≥ 1 . . . 28

5.4 ARLs for shifts in θ with θj = (βjT, δσ)T and τj2 = τ2/δ2 for j ≥ 1 . . . 29

1

Introduction

1.1

Motivation

Quality improvement is a key factor for keeping competitiveness in the international market. Statistical process control (SPC) is used for applications of industrial processes by statistical methods that can help practitioner to find out the problems in processes. Roberts (1959) proposed an exponentially weighted moving average (EWMA) control chart to detect a small shift in the process mean, and it nowadays is a widely studied and accepted alternative to the traditional Shewhart ¯X control chart (Shewhart, 1931). The likelihood ratio (LR) test is commonly used to test hypotheses for a parametric family in statistical literature because of its high power in detecting a large shift or for a large sample size. In many situations, the quality of a process may be better characterized and summarized by the relationship between the response variable and one or more explanatory variables. Thus, an EWMA control chart based on LR test statistics for monitoring general linear profiles is proposed in the paper and then some properties of the proposed methodology are investigated.

1.2

Literature Review

SPC refers to statistical methods which are extensively used to monitor and improve the quality of industrial processes and others operations. Most of the studies on SPC are focused on the charting skill which is used to monitor processes in order to distinguish

special significant reasons of variation from general allocation reasons of variation in processes. In the Phase I stage, a set of process data is gathered to construct trial control limits that determine whether or not the process has been in control over the period of time, and then to model the in-control process so that reliable control limits of the control chart can be established for the later Phase II stage. In the Phase II stage, process data are compared with a pre-established standard from the previous Phase I stage to determine whether the process is in control or not.

Shewhart (1931) proposed a Shewhart ¯X control chart which is briefly introduced as follows: At time j (≥ 1), let Xj1, Xj2, . . . , Xjnbe i.i.d. N(µj, σ2) observations, where µj

denotes the unknown process mean at time j and σ the known positive process standard deviation. Set ¯Xj ≡ Pni=1Xji/n for j ≥ 1. Then an out-of-control signal occurs at

time j (≥ 1) if√n| ¯Xj− µ|/σ > L, where µ denotes the known in-control process mean

and L is chosen to achieve a specified in-control average run length (ARL).

Page (1954) proposed a cumulative sum (CUSUM) control chart which is briefly introduced as follows: Let X1, X2, . . . be independent observations such that Xj ∼

N(µj, σ2), where µj denotes the unknown process mean at time j and σ the known

positive process standard deviation. For j ≥ 1, two statistics Cj+and Cj− are iteratively

defined as

C_j+ _{≡ max{0, Cj−1}+ + Xj − (µ + K)}

and

of the specified shift in process mean that should be quickly detected by the scheme. Then an out-of-control signal occurs at time j (≥ 1) if max{Cj+, Cj−} ≥ Hσ, where H

(> 0) is chosen to achieve a specified in-control ARL.

Roberts (1959) proposed an EWMA, originally called geometric moving average, con-trol chart which is briefly introduced as follows: At time j (≥ 1), let Xj1, Xj2, . . . , Xjn

be i.i.d. N(µj, σ2) observations, where µj denotes the unknown process mean at time j

and σ the known positive process standard deviation. The EWMA sequence is

Wj ≡ (1 − λ)Wj−1+ λ( ¯Xj − µ) = λ j−1

X

k=0

(1 − λ)k( ¯X_{j−k − µ), j = 1, 2, . . . ,}

where W0 ≡ 0, µ denotes the known in-control process mean, and λ is a smoothing

parameter chosen in (0, 1]. Observe that the standard deviation of Wj is

σj = s λ[1 − (1 − λ)2j_] n(2 − λ) σ → s λ n(2 − λ)σ

as j → ∞. Then out-of-control signal occurs at time j (≥ 1) if |Wj|/σj > L, where L

is chosen to achieve a specified in-control ARL.

An EWMA control chart is typically designed in the Phase II stage for a manu-facturing process. Stoumbos et al. (2003) investigated the performance of Shewhart, EWMA, and CUSUM control charts for detecting sustained shifts and drifts in the pro-cess mean and variance. In practice, the in-control propro-cess distribution is rarely known exactly, and thus control charts are usually constructed using an approximate in-control process distribution estimated from some historical in-control process data. Jones et al. (2001) utilized a numerical procedure to study the run-length (RL) distribution of

an EWMA control chart using estimated in-control process parameters. Jensen et al. (2006) provided a review of the literature that considered the effect of in-control process parameter estimation on control charts and concluded that the influence of in-control process parameter estimation on control charts should not be ignored.

In most SPC applications, it is assumed that the quality of a process can be suit-ably represented by the joint distribution of quality characteristics. However, in many situations, the quality of a process may be better characterized and summarized by the relationship between the response variable and one or more explanatory variables. Sev-eral methods for monitoring linear profiles have been proposed in literature, e.g., Kim et al. (2003) proposed a control chart for monitoring simple linear profiles in the known in-control process parameter case, and Zou et al. (2006) proposed a control chart based on the LR test statistics for monitoring simple linear profiles to detect a sustained shift in the unknown in-control process parameter case.

Through the modern technology that allows simultaneously monitoring all key qual-ity characteristics during a manufacturing process, the monitored qualqual-ity characteris-tics are generally dependent each other. The purpose of multivariate techniques is to study whether quality characteristics are simultaneously in control or not. Lowry et al. (1992) proposed a multivariate exponentially weighted moving average (MEWMA) con-trol chart giving guidelines for designing this easy-to-implement multivariate procedure. The performance of their control chart is similar to that of a multivariate cumulative sum (MCUSUM) control chart (Crosier, 1988) in detecting a small shift in the process mean.

Zou et al. (2007) [13] proposed an MEWMA control chart for monitoring general lin-ear profiles in the known in-control process parameter case, which is briefly introduced as follows: The process is called in control at time j (≥ 1) if

y_j = Xβ + σεj,

where yj is the n × 1 response vector at time j, X is the known n × p model matrix

of full rank p (< n), β (≡ (β0, β1, ..., βp−1)T) is the known p × 1 vector of real-valued

in-control process regression parameters, σ is the known positive in-control process scale parameter, and εjs are i.i.d. Nn(0n×1, In) standardized error vectors. Set Zj ≡

(ZT j (β), Zj(σ))T, where Zj(β) ≡ ( ˆβj − β)/σ and Zj(σ) ≡ Φ−1(Fχ2 n−p((n − p)ˆσ 2 j/σ2)) with ˆβ_{j ≡ (X}TX)−1_XT_y j, ˆσj2 ≡ (yj − X ˆβj)T(yj − X ˆβj)/(n − p), Φ−1(·) denoting the

inverse of the standard normal cumulative distribution function (c.d.f.), and Fχ2 n−p(·)

the c.d.f. of the chi-squared distribution with n − p degrees of freedom. The MEWMA sequence is Wj ≡ λZj+ (1 − λ)Wj−1 = λ j−1 X k=0 (1 − λ)kZ_j−k, j = 1, 2, . . . ,

where W0 ≡ 0(p+1)×1 and λ is a smoothing parameter chosen in (0, 1]. Observe that

the in-control covariance matrix of Wj is

λ[1 − (1 − λ2j_)]

2 − λ Σ → λ 2 − λΣ

as j → ∞, where Σ =      (XTX)−1 ₀ p×1 0T p×1 1      .

Then an out-of-control signal occurs at time j (≥ 1) if (2−λ)WT

j Σ−1Wj/λ > L, where

L is chosen to achieve a specified in-control ARL.

Zou et al. (2007) [16] proposed a self-starting control chart based on recursive resid-uals for monitoring simple linear profiles to detect a sustained shift in the process intercept, slope, or standard deviation. Although current SPC methods mostly focus on detecting a sustained shift in the process mean, time-varying shifts in the process mean frequently occur in industrial applications. Thus, Zou et al. (2009) investigated several control charts for detecting drifts in the process mean. Recently, Zou et al. (2010) proposed an EWMA control chart based on the LR test statistics for monitoring the process mean and variance.

1.3

Outline

The paper is organized as follows. In Section 2, an EWMA control chart based on LR test statistics for monitoring general linear profiles is proposed and then some prop-erties of the proposed monitoring scheme are investigated. In Section 3, a simulation study is presented to illustrate the proposed methodology. In Section 4, performance comparisons and conclusions are given. Possible future work is discussed in Section 5.

2

An EWMA Control Chart Based on LR Test

Statis-tics for Monitoring General Linear Profiles

In this section, an EWMA control chart based on LR test statistics for monitoring general linear profiles is proposed and then some properties of the proposed methodol-ogy are investigated.

2.1

Model

In this subsection, consider the following model for the purpose of monitoring general linear profiles: Suppose that

yj = Xjβj+ σjεj, j = 1, 2, ..., (2.1)

where yj is the nj× 1 response vector at time j, Xj is the known nj× p model matrix of

full rank p (< nj) at time j, βj (≡ (βj0, βj1, ..., βj,p−1)T) is the unknown p × 1 vector of

real-valued process regression parameters at time j, σj is the unknown positive process

scale parameter at time j, and εjs are independent Nnj(0nj×1, Inj) standardized error

vectors. Set θj ≡ (βTj, σj)T (∈ Rp × (0, ∞) ≡ Θ) for j ≥ 1. Then θj is the process

parameter vector at time j for j ≥ 1. For j ≥ 1, the process is called in control at time j if θj = θ, where θ (≡ (βT, σ)T ≡ (β0, β1, . . . , βp−1, σ)T ∈ Θ) is the in-control process

When θ is unknown, it is assumed that the historical in-control process data {X0, y0}

are available such that

y0 = X0β+ σε0, (2.2)

where y0 is the n0 × 1 response vector, X0 is the known n0 × p model matrix of full

rank p (< n0), and ε0 is an Nn0(0n0×1, In0) standardized error vector independent of

{εj : j ≥ 1}.

2.2

Known and Unconstrained Case

In this subsection, consider the case where the process parameter vector at time j (≥ 1) may be different from the known in-control process parameter vector. For j ≥ 1, the joint probability density function (p.d.f.) of yj is

f (yj; θj) = 1 (2π)nj/2σnj j exp  −kyj − Xjβjk 2 2σ2 j  = 1 (2π)nj/2σnj j exp ( −( ˆβj − βj) T_XT jXj( ˆβj − βj) + njˆσj2 2σ2 j ) ,

where ˆβj ≡ (XTjXj)−1XjTyj and ˆσj ≡ k[Inj− Xj(X

T

jXj)−1XTj]yjk/√nj. Then, for j ≥ 1,

the log-likelihood function for θj is

ℓj(θj) ≡ log[f(yj; θj)] = −n₂j log(2π) −n₂j log(σ2_{j) −} ( ˆβj− βj) T_XT jXj( ˆβj − βj) + njσˆ2j 2σ2 j .

and thus the maximum log-likelihood is ℓj( ˆθj) = − nj 2 − nj 2 log(2π) − nj 2 log(ˆσ 2 j), where ˆθ_{j (≡ ( ˆ}βT

j , ˆσj)T) is the maximum likelihood estimator (MLE) of θj. Hence the

LR test statistic at time j (≥ 1) is

Wj ≡ 2[ℓj( ˆθj) − ℓj(θ)] = ( ˆβj − β)TXTjXj( ˆβj− β) σ2 + nj  −1 + σˆ 2 j σ2 − log _σˆ2 j σ2  , (2.3) where ( ˆβj− β)TXTjXj( ˆβj− β) σ2 ∼ σ2 j σ2χ 2 p (βj − β)TXTjXj(βj − β) σ2 j ! , ˆ σ2 j σ2 ∼ σ2 j σ2 χ2 nj−p nj ,

and ( ˆβj − β)TXTjXj( ˆβj − β)/σ2 is independent of nj[ˆσj2/σ2− 1 − log(ˆσ2j/σ2)].

When the process is out of control at time j (≥ 1), the distribution of Wj only

depends on p, nj − p, σj/σ, and (βj − β)TXTjXj(βj − β)/σ2j (≡ τj2), and Wj d

→ ∞ as min{nj, the minimum eigenvalue of XTjXj} → ∞, where

Eθj(Wj) = τ 2 j σ2 j σ2 + nj  lognj 2  − ψ nj₂− p  − 1 + σ 2 j σ2 − log _σ2 j σ2  , and Varθj(Wj) = Eθj(W 2 j) − [Eθj(Wj)] 2

some relevant formulas about both ψ(x) and ψ′_{(x). See Appendix A.2 for E} θj(W

2 j) for

j ≥ 1.

When the process is in control at time j (≥ 1), the distribution of Wj only depends

on p and nj − p, and nj  −1 + σˆ 2 j σ2 − log _σˆ2 j σ2  = nj( ˆσ 2 j − σ2)2 2σ4 + Op  1 √_n j  d → Z2 as nj → ∞, where Z ∼ N(0, 1). Then Wj d → χ2 p+1 as nj → ∞, where Eθ(Wj) ≡ Eθj(Wj)|θj=θ = nj  lognj 2  − ψ nj₂− p  = p + 1 + O 1 nj  , and Varθ(Wj) ≡ Varθj(Wj)|θj=θ = nj  njψ′  nj− p 2  − 2  = 2(p + 1) + O 1 nj  .

See Appendix A.3 for both Eθ(Wj) and Varθ(Wj) as nj → ∞.

For j ≥ 1, set ¯ Wj ≡ Wj− Eθ(Wj) pVarθ(Wj) . (2.4)

2.3

Known and Constrained Case

In this subsection, consider the case where the process parameter vector at time j (≥ 1) may be different from the known in-control process parameter vector with the

constraint σj ≥ σ. Then the maximum log-likelihood is ℓj(θj∗) = − nj 2 log(2π) − nj 2 log(σ ∗2 j ) − njσˆj2 2σ∗2 j , where θ∗

j (≡ (β∗Tj , σj∗)T = ( ˆβjT, max{ˆσj, σ})T) is the MLE of θj. Hence the LR test

statistic at time j (≥ 1) is W_j∗ _{≡ 2[ℓ}j(θj∗) − ℓj(θ)] = ( ˆβj − β) T_XT jXj( ˆβj− β) σ2 + nj  −1 + σˆ 2 j σ2 − log _σˆ2 j σ2  · 1{ˆσj>σj}, (2.5)

where ( ˆβj−β)TXTjXj( ˆβj−β)/σ2is independent of nj[−1+ˆσj2/σ2−log(ˆσj2/σ2)]·1{ˆσj>σj}.

When the process is out of control at time j (≥ 1), the distribution of W∗ j only

depends on p, nj − p, σj/σ and τj2, and Wj∗ d

→ ∞ as min{nj, the minimum eigenvalue

of XT_jXj} → ∞, where Eθj(W ∗ j) = (p + τj2) σ2 j σ2 − nj  1 + log  _σ2 j njσ2  E(1{Hj≥aj}) +σ 2 j

σ2E(Hj· 1{Hj≥aj}) − njE(log(Hj) · 1{Hj≥aj}),

and Varθj(W ∗ j) = Eθj(W ∗2 j ) − [Eθj(W ∗ j)]2 with Hj ≡ njσˆ2j/σj2 ∼ χ2nj−p, aj ≡ njσ 2_/σ2

j, and 1{Hj>nj}denoting the indicator function

for {Hj > nj}. See Appendix A.4 for both Eθj(Wj∗) and Eθj(Wj∗2) in detail.

When the process is in control at time j (≥ 1), the distribution of W∗

on p and nj − p, and nj  −1 + σˆ 2 j σ2 − log _σˆ2 j σ2  · 1{ˆσj>σ} = nj(ˆσ 2 j − σ2)2 2σ4 · 1{√nj(ˆσj2−σ 2_)/(√_2σ2 )>0}+ Op  1 √_n j  d → Z2· 1{Z>0} as nj → ∞, where Z ∼ N(0, 1). Then Wj∗ d → χ2 p/2 + χ2p+1/2 as nj → ∞, where Eθ(Wj∗) ≡ Eθj(Wj∗)|θj=θ

= p + [njlog(nj) − nj]E(1{Hj>nj}) + E[(Hj· 1{Hj>nj}) − njE(log(Hj) · 1{Hj>nj}),

and Varθ(Wj∗) ≡ Varθj(W ∗ j)|θj=θ = Eθ(W ∗2 j ) − (Eθ(Wj∗))2

= p(2 + p) + [njlog(nj) − nj]2E(1{Hj>nj}) + 2[njlog(nj) − nj]E(Hj· 1{Hj>nj})

−2nj[njlog(nj) − nj]E(log(Hj) · 1{Hj>nj}) + E(H

2

j · 1{Hj>nj})

+n2_jE([log(Hj)]2· 1_{Hj>nj})

−2njE(Hjlog(Hj) · 1{Hj>nj}) − [Eθ(Wj∗)]2.

See Appendix A.4 for both Eθ(Wj∗) and Varθ(Wj∗).

For j ≥ 1, set ¯ W_j∗ _≡ W ∗ j − Eθ(Wj∗) q Varθ(Wj∗) . (2.6)

2.4

Unknown and Unconstrained Case

In this subsection, consider the case where the process parameter vector at time j (≥ 1) may be different from the unknown in-control process parameter vector. For j (≥ 1), the joint p.d.f. of (yT 0, yTj)T is f (y0, yj; θ, θj) = f (y0; θ) · f(yj; θj) = 1 (2π)n0/2σn0 exp  −ky0− X0βk 2 2σ2  · f(yj; θj) = 1 (2π)n0/2σn0 0 exp ( −( ˆβ− β) T_XT 0X0( ˆβ− β) + n0σˆ2 2σ2 ) · f(yj; θj),

where ˆβ _{≡ (X}T₀X₀)−1XT₀y0 and ˆσ ≡ |[In0 − X0(X

T

0X0)−1XT0]y0k/√n0.

Then, for j ≥ 1, the log-likelihood function for (θT_{, θ}T j)T is ℓ0,j(θ, θj) ≡ log[f(y0; θ)] + ℓj(θj) = −n0 2 log(2π) − n0 2 log(σ 2_{) −}( ˆβ− β)TX T 0X0( ˆβ− β) + n0ˆσ2 2σ2 + ℓj(θj) ≡ ℓ0(θ) + ℓj(θj)

and thus the maximum log-likelihood is

ℓ0,j( ˆθ, ˆθj) = ℓ0( ˆθ) + ℓj( ˆθj) = − n0 2 − n0 2 log(2π) − n0 2 log(ˆσ 2_{) + ℓ} j( ˆθj) = −n0+ n₂ j −n0+ n₂ j log(2π) − n₂0 log(ˆσ2_{) −} nj 2 log(ˆσ 2 j), where ( ˆθT, ˆθ_jT)T _{(≡ ( ˆβ}T_{, ˆσ, ˆθ}T

at time j (≥ 1) is

W0,j ≡ 2[ℓ0,j( ˆθ, ˆθj) − ℓ0,j( ˜θ, ˜θj)] = (n0+ nj) log(˜σ2) − n0log(ˆσ2) − njlog(ˆσ2j), (2.7)

where ˜θ _{≡ ( ˜}βT_{, ˜σ)}T _{and ˜θ} j ≡ ( ˜βTj, ˜σj)T = ˜θ with ˜ β_{≡ (X}T₀X0+ XTjXj)−1(XT0y0 + XTjyj) and ˜ σ ≡ k[In0+nj− (X T 0, XTj)T(XT0X0+ XTjXj)−1(XT0, XTj)](y0T, yjT)Tk √_n 0 + nj = s n0σˆ2+ njσˆj2+ ( ˆβj − ˆβ)T[(X0TX0)−1+ (XTjXj)−1]−1( ˆβj− ˆβ) n0+ nj .

When the process is in control at time j (≥ 1), the distribution of W0,j only depends

on p, n0 − p, and nj − p, and W0,j a.s.→ Wj as min{n0, the minimum eigenvalue of

XT₀X0} → ∞, see Appendix A.5 in detail. Then W0,j d → χ2 p+1 as min{n0/nj, nj} → ∞. For j ≥ 1, set ¯ W0,j ≡ W0,j− Eθ(W0,j) pVarθ(W0,j) , (2.8)

where Eθ(W0,j) ≡ Eθ,θj(W0,j)|θj=θ = n0  logn0 2  − ψ n0₂− p  + nj  lognj 2  − ψ nj₂− p  −(n0 + nj)  log n0+ nj 2  − ψ n0+ n₂j − p  = p + 1 + O  1 min{n0, nj}  , and Varθ(W0,j) ≡ Varθ,θj(W0,j)|θj=θ = n0  n0ψ′  n0− p 2  − 2  + nj  njψ′  nj − p 2  − 2  −(n0+ nj)  (n0+ nj)ψ′  n0+ nj − p 2  − 2  = 2(p + 1) + O  1 min{n0, nj}  .

2.5

Unknown and Constrained Case

In this subsection, consider the case where the process parameter vector at time j (≥ 1) may be different from the unknown in-control process parameter vector with the constraint σj ≥ σ. Then the maximum log-likelihood is

ℓ0,j(θ∗, θj∗) = − n0 + nj 2 − n0 + nj 2 log(2π) − n0 2 log(σ ∗2_{) −} nj 2 log(σ ∗2 j ),

where (θ∗_{, θ}∗ j)T (≡ (β∗T, σ∗, βj∗T, σj∗)T = ( ˆβT, σ∗, ˆβTj, σj∗)T) is the MLE of (θT, θjT)T with σ∗2= ˆσ2_{· 1{ˆ}σj≥ˆσ}+ n0σˆ2+ njσˆ2j n0+ nj · 1{ˆσj<ˆσ} and σj∗2= ˆσj2· 1{ˆσj≥ˆσ}+ n0σˆ2 + njσˆj2 n0 + nj · 1{ˆσj<ˆσ} .

Hence the LR test statistic at time j (≥ 1) is

W_0,j∗ _{≡ 2[ℓ}0,j(θ∗, θ∗j) − ℓ0,j( ˜θ, ˜θj)] = (n0+ nj) log(˜σ 2

) − n0log(σ∗2) − njlog(σj∗2). (2.9)

When the process is in control at time j (≥ 1), the distribution of W∗

0,j only depends

on p, n0 − p, and nj − p, and W0,j∗ a.s.

→ W∗

j as min{n0, the minimum eigenvalue of

XT₀X_{0} → ∞. Then W0,j}∗ _{→ χ}d _p2/2 + χ2p+1/2 as min{n0/nj, nj} → ∞. For j ≥ 1, set

¯ W_0,j∗ _≡ W ∗ 0,j− Eθ(W0,j∗ ) q Varθ(W0,j∗ ) , (2.10)

where both Eθ(W0,j∗ ) ≡ Eθ,θj(W0,j∗ )|θj=θ and Varθ(W0,j∗ ) ≡ Varθ,θj(W0,j∗ )|θj=θ are given

in Appendix A.6.

2.6

Proposed Monitoring Scheme

Based on the standardized LR test statistics, the EWMA sequence is defined as

where U0 ≡ 0, λ is a smoothing parameter in (0, 1], and ˜Wj = ¯Wj in Section 2.2, or ¯Wj∗

in Section 2.3, or ¯W0,j in Section 2.4, or ¯W0,j∗ in Section 2.2. Then

Uj = λ j−1

X

k=0

(1 − λ)kW˜_j−k, j = 1, 2, . . . . (2.12)

Due to we use the LR test statistics, then the larger value is more likely become out-of-control. An out-of-control signal occurs at time j (≥ 1) if

U_j∗ _≡ Uj − Eθ(Uj) pVarθ(Uj) > C, (2.13) where Eθ(Uj) ≡ Eθ,θ1,...,θj(Uj)|θ1=...=θj=θ = 0, Varθ(Uj) ≡ Varθ,θ1,...,θj(Uj)|θ1=...=θj=θ = λ[1 − (1 − λ)2j_] 2 − λ (2.14) for known θ, Varθ(Uj) ≡ Varθ,θ1,...,θj(Uj)|θ1=...=θj=θ = λ[1 − (1 − λ) 2j_] 2 − λ + 2λ 2 X 0≤k1<k2≤j−1 (1 − λ)k1+k2 Covθ( ˜W_j−k1, ˜Wj−k2) (2.15)

for unknown θ with Covθ( ˜W_j−k1, ˜Wj−k2) ≡ Covθ,θ1,...,θ_j−k1( ˜Wj−k1, ˜Wj−k2)|θ1=...=θ_j−k1=θ,

and C is chosen to achieve a specified in-control ARL. When θ is unknown, Covθ( ˜Wj−k1, ˜Wj−k2)

depends on X0, Xj−k1, and Xj−k2 for 0 ≤ k1 < k2 ≤ j − 1. In the paper, each

3

A Simulation Study

In this section, in order to study the performance of the proposed EWMA monitoring scheme based on the LR test statistics, the out-of-control ARLs are compared with those in Kim et al. (2003) and Zou et al. (2007) [13] through a simulation study.

First, consider the case where the in-control parameter vector θ is known. Then C in equation (2.13) can be evaluated as follows:

Step 1 : Choose the desired control ARL and the smoothing constant λ. Here the in-control ARL and λ is chosen as, respectively, 200 and 0.2. Then we can start to construct the EWMA control chart, and in order to conveniently introduce the steps, we consider the case where the in-control parameter vector θ is constrained, and unconstrained case is also suitable for use.

Step 2 : Due to the in-control ARL = 200, so we generate i.i.d. in-control W∗ 1, . . . ,

W∗

200, and we want to standardized the Wj∗, so we evaluate the Eθ(Wj∗) and Varθ(Wj∗),

which are used to evaluate ˜W∗

j. And the data Xj can be random when the p and nj− p

is known.

Step 3 : Evaluate the EWMA sequences {U1, . . . , U200} by equation (2.12), and we want

to standardized the Uj, so we evaluate the Varθ(U1), . . . , Varθ(U200) by equation (2.14),

which use to evaluate U∗

j by equation (2.13).

Step 4 : We want to know the limit value every time, which use to estimate the C. So we find the maximum of {U∗

Step 5 : Due to the in-control ARL is the mean of run length, so we think that the C is between the first sample quartile and third sample quartile of Q(r)_{, r = 1, . . . , 50, 000.}

Therefore we sort the Q(r)_{, r = 1, . . . , 50, 000, and evaluate its first sample quartile}

(≡ q1), second sample quartile (≡ q2) and third sample quartile (≡ q3).

Step 6 : Choose C as each of q1, q2 and q3 to evaluate the corresponding in-control ARL

for 50,000 times. Use the interpolation method to find the value C with the specified in-control ARL= 200.

To evaluate the out-of-control ARL, consider the same shifts in θ in Zou et al. (2007) [13]. The out-of-control ARL can be evaluated as follows:

Step 1 : Due to we want to evaluate the out-of-control ARL, we simulate 50,000 times to average the RL. Then we generate i.i.d. W∗

1 and W2∗ and standardized them, which

use to evaluate U1 and U2.

Step 2 : Standardized both U1 and U2 to U1∗ and U2∗, then to judge U1∗ whether larger

than C. If it is true, then the RL= j = 1, else continue the process: generate i.i.d. W∗

j, j = 1, . . ., until the value Uj∗ > C, j ≥ 1. Then the stopping time ℓ is the RL of

this time.

Step 3 : Repeat Steps 1 and 2 50,000 times to average the RLs, which can obtain the out-of-control ARL.

Next, consider the case where the in-control parameter vector θ is unknown. Then C in equation (2.13) can be evaluated as follows:

Step 1 : Choose the desired in-control ARL and the smoothing constant λ. Here the in-control ARL and λ is chosen as, respectively, 200 and 0.2. Then we can start to

construct the EWMA control chart, and in order to conveniently introduce the steps, we consider the case where the in-control parameter vector θ is constrained, and un-constrained case is also suitable for use.

Step 2 : Due to the in-control ARL = 200, so we generate i.i.d. in-control W∗ 0,1,

. . . , W0,200∗ , and we want to standardized the W0,j∗ , so we evaluate the Eθ(W0,j∗ ) and

Varθ(W0,j∗ ), which are used to evaluate ˜Wj∗. The W0,j∗ and W0,j∗ ′ depends on all of p,

n0− p, n1− p, n2− p, . . . (j 6= j′), so the X0 and Xj must be a design matrix. e.x.,

X_{: n × p X}j =           X .. . X           mjn×p , X0 =           X .. . X           m0n×p m0, mj ≥ 1 for all j ≥ 1

Step 3 : Evaluate the EWMA sequences {U1, . . . , U200} by equation (2.12), and we want

to standardized the Uj, so we evaluate the Varθ(U1), . . . , Varθ(U200) by equation (2.15),

which use to evaluate U∗

j by equation (2.13).

Step 4 : We want to know the limit value every time, which use to estimate the C. So we find the maximum of {U∗

1, . . . , U200∗ } (≡ Q(1)), repeat above steps for 50,000 times

to obtain Q(1)_{, . . . , Q}(50,000)_.

Step 5 : Due to the in-control ARL is the mean of run length, so we think that the C is between the first sample quartile and third sample quartile of Q(r)_{, r = 1, . . . , 50, 000.}

Therefore we sort the Q(r)_{, r = 1, . . . , 50, 000, and evaluate its first sample quartile}

for 50,000 times. Use the interpolation method to find the value C with the specified in-control ARL= 200.

To evaluate the out-of-control ARL, consider the same shifts in θ in Zou et al. (2007) [13]. The out-of-control ARL can be evaluated as follows:

Step 1 : Due to we want to evaluate the out-of-control ARL, we simulate 50,000 times to average the RL. Then we generate i.i.d. W∗

0,1 and W0,2∗ and standardized them, which

use to evaluate U1 and U2.

Step 2 : Standardized both U1 and U2 to U1∗ and U2∗, then to judge U1∗ whether larger

than C. If it is true, then the RL= j = 1, else continue the process: generate i.i.d. W∗

0,j, j = 1, . . ., until the value Uj∗ > C, j ≥ 1. Then the stopping time ℓ is the RL of

this time.

Step 3 : Repeat Steps 1 and 2 50,000 times to average the RLs, which can obtain the out-of-control ARL.

Then we consider the simplest case of model (2.1) to compared with Kim et al. (2003) and Zou et al. (2007) [13] through a simulation study, in which the parameters in the in-control model are β = (3, 2)T_{, σ = 1, and}

X0 =           X .. . X           mn×p X= (1_4×1, x) with x = (2, 4, 6, 8)T_.

4

Comparisons and Conclusions

In this section, first, we compared our propose EWMA control chart with the ZTW and KMW control charts from Kim et al. (2003) and Zou et al. (2007) [13]. We follow the table by Zou et al. (2007) [13], which in term of out-of-control ARL in the simplest case of model (2.1) for shifts in β0, β1 and σ , respectively. Because of our performance

under unconstrained case is not better, therefore we don’t present in this. Then the performance under known and constrained case are given in Table 5-1, 5-2 and 5-3. On Table 5-1 and 5-2, those table show that for β0 and β1, our charts’ performance is not

better than others with small shifts, but it is better than others with large shifts. We think that, the LR test statistic is more powerful with large shift by nature. But on Table 5-3 shows that change for σ, our performance is better than others with any shifts. Above result show that, our propose EWMA chart do not have better performance for shifts in intercept and slope, but have better performance for shifts in σ. We think that it maybe because of the method of ZTW is using MEWMA, the proportion of β and σ are the same, however the proportion of σ is larger than β by our proposed.

In practice, the in-control process parameters are rarely known, and thus control charts usually are constructed using estimates from the historical in-control process data in place of the unknown in-control parameters. Therefore we present the performance under unknown and constrained case are given in Table 5-1, 5-2 and 5-3 left side, those tables show that the performance is similar to known parameter about 125 historical

From the section 2, we can know that the out-of-control ARL is depends on all of p, nj − p, σj/σ and τj2, where the noncentrality parameter

τj2 = (βj − β)TXjTXj(βj − β)/σj2,

Xj = X for all j ≥ 1,

and both p and nj − p are constant. Therefore we present the performance under

combinations of σj/σ and τ2σj2/σ2, are given in Table 5-4 and 5-5. Those tables show

that the performance reduce fast at σj/σ.

In this paper, our propose EWMA chart based on LR test statistics have better performance for shift in σ. And we provide the mean and variance for likelihood ratio test statistic that not presented in any papers. This can be taken as a reference.

5

Future Work

In this section, we provide some idea to extend or improve study. The future research could substitute: Consider the case where only β may change; or only σ may change to observe the performance whether more better for β0 and β1 with small shift. Another

the research could consider possible transformation for w∗

j obtain better performance for

β0 and β1 with small shifts. e.g., T = Φ−1(Fθ(Wj∗)), then in the in-control T ∼ N(0, 1),

then C can be easily found in some literatures, where F (·) is c.d.f., Φ−1_{(·) is inverse}

of the standard normal cumulative distribution faction. Zou et al. (2010) proposed for jointly monitoring the process mean and variance, consider an MEWMA control chart based on LR test statistics for monitoring general linear profiles.

References

[1] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with For-mulas, Graphs, and Mathematical Tables, volume 55. Dover publications., New York, 1964.

[2] R.B. Crosier. Multivariate generalizations of cumulative sum quality-control schemes. Technometrics, 30(3):291–303, 1988.

[3] W.A. Jensen, L.A. Jones-Farmer, C.W. Champ, and W.H. Woodall. Effects of parameter estimation on control chart properties: a literature review. Journal of Quality Technology, 38(4):349–364, 2006.

[4] L.A. Jones, C.W. Champ, and S.E. Rigdon. The performance of exponen-tially weighted moving average charts with estimated parameters. Technometrics, 43(2):156–167, 2001.

[5] K. Kim, M.A. Mahmoud, and W.H. Woodall. On the monitoring of linear profiles. Journal of Quality Technology, 35(3):317–317, 2003.

[6] C.A. Lowry, W.H. Woodall, C.W. Champ, and S.E. Rigdon. A multivariate ex-ponentially weighted moving average control chart. Technometrics, 34(01):46–53, 1992.

[8] S.W. Roberts. Control chart tests based on geometric moving averages. Techno-metrics, 1(03):239–250, 1959.

[9] W.A. Shewhart. Economic Control of Quality of Manufactured Product, volume 509. American Society for Quality, New York, 1931.

[10] Z.G. Stoumbos, M.R. Reynolds, and W.H. Woodall. Ch. 13. control chart schemes for monitoring the mean and variance of processes subject to sustained shifts and drifts. Handbook of Statistics, 22:553–571, 2003.

[11] C. Zou, Y. Liu, and Z. Wang. Comparisons of control schemes for monitoring the means of processes subject to drifts. Metrika, 70(2):141–163, 2009.

[12] C. Zou and F. Tsung. A multivariate sign ewma control chart. Technometrics, 53(1):84–97, 2011.

[13] C. Zou, F. Tsung, and Z. Wang. Monitoring general linear profiles using multivari-ate exponentially weighted moving average schemes. Technometrics, 49(4):395–408, 2007.

[14] C. Zou, J. Zhang, and Z. Wang. A control chart based on likelihood ratio test for monitoring process mean and variability. Quality and Reliability Engineering International, 26(1):63–73, 2010.

[15] C. Zou, Y. Zhang, and Z. Wang. A control chart based on a change-point model for monitoring linear profiles. IIE Transactions, 38(12):1093–1103, 2006.

Table 5.1: ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in β0 with

θj= (β0+ δ0σ, β1, σ)T and τ_j2= 4 δ02 for j ≥ 1

β0

EWMA EWMA EWMA EWMA δ0 τ 2 j m0= 5 m0= 25 m0= 125 m0= ∞ ZTW KMW .1000 0.0400 178.2 174.5 172.9 171.1 131.5 133.7 .2000 0.1600 134.6 120.3 116.3 113.1 59.9 59.1 .3000 0.3600 88.7 72.7 68.6 67.2 29.6 28.3 .4000 0.6400 55.0 41.2 38.9 38.1 17.2 16.2 .5000 1.0000 33.6 24.6 22.7 22.1 11.5 10.7 .6000 1.4400 21.5 15.2 13.9 13.8 8.5 7.9 .8000 2.5600 9.9 6.9 6.5 6.1 5.8 5.1 1.0000 4.0000 5.5 3.9 3.6 3.4 4.1 3.8 1.5000 9.0000 2.2 1.6 1.6 1.5 2.6 2.4 2.0000 16.0000 1.4 1.1 1.1 1.1 2.0 1.9

Table 5.2: ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in β1 with

θj= (β0, β1+ δ1σ, σ)T and τ 2 j = 120 δ 2 1for j ≥ 1 β1

EWMA EWMA EWMA EWMA

δ1 τ_j2 m0= 5 m0= 25 m0= 125 m0= ∞ ZTW KMW .0250 0.0750 168.5 157.6 155.5 153.0 99.0 101.6 .0375 0.1688 133.6 117.8 115.3 112.3 57.4 61.0 .0500 0.3000 102.9 83.5 79.4 77.4 35.0 36.5 .0625 0.4688 74.2 57.8 53.7 53.2 23.1 24.6 .0750 0.6750 53.5 39.2 36.4 35.4 16.4 17.0 .1000 1.2000 27.7 19.2 17.6 17.4 9.8 10.3 .1250 1.8750 15.3 10.6 9.6 9.5 6.9 7.2 .1500 2.7000 9.4 6.4 5.9 5.8 5.3 5.5 .2000 4.8000 4.4 3.1 2.9 2.8 3.7 3.8 .2500 7.5000 2.6 2.0 1.8 1.8 2.9 2.9

Table 5.3: ARL comparisons among EWMA, ZTW, and KMW control charts for shifts in σ with θj= (β0, β1, δσ)T for j ≥ 1

σ

EWMA EWMA EWMA EWMA

δ m0= 5 m0= 25 m0= 125 m0= ∞ ZTW KMW 1.1000 72.6 56.8 51.1 51.1 76.2 72.8 1.1500 48.9 35.1 30.9 30.2 48.7 48.1 1.2000 34.4 23.2 20.5 20.0 33.2 33.5 1.2500 25.1 16.7 14.7 14.6 24.1 24.9 1.3000 19.3 12.6 11.1 11.0 18.4 19.4 1.4000 12.3 7.9 7.0 7.0 12.1 12.7 1.6000 6.5 4.3 4.0 4.0 7.0 7.2 1.8000 4.2 2.9 2.8 2.8 4.9 5.1 2.2000 2.5 1.9 1.8 1.8 3.1 3.2 2.6000 1.8 1.4 1.4 1.4 2.3 2.5 3.0000 1.5 1.3 1.3 1.3 1.9 2.1

Table 5.4: ARLs for shifts in θ with θj= (βjT, δσ)T and τ 2 j = τ 2 /δ2 for j ≥ 1 m0= 5 m0= 25 m0= 125 δ m0= ∞ 1.10 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18 1.19 1.110 1.111 1.112 0.02 200.0 72.6 32.5 16.6 9.8 6.3 4.5 3.3 2.6 2.1 1.8 1.6 1.4 200.0 56.8 20.9 10.4 6.2 4.1 3.0 2.3 1.9 1.6 1.4 1.3 1.2 200.0 51.2 18.9 9.5 5.6 3.8 2.8 2.2 1.8 1.6 1.4 1.3 1.2 200.0 51.1 18.2 9.0 5.4 3.7 2.7 2.1 1.8 1.5 1.4 1.3 1.2 0.22 178.2 70.1 31.4 16.3 9.6 6.3 4.4 3.3 2.6 2.1 1.8 1.6 1.4 174.5 51.3 20.3 10.3 6.1 4.1 3.0 2.3 1.9 1.6 1.4 1.3 1.2 172.9 47.2 18.5 9.3 5.6 3.8 2.8 2.2 1.8 1.6 1.4 1.3 1.2 171.1 43.2 17.1 8.7 5.3 3.6 2.7 2.1 1.8 1.5 1.4 1.3 1.2 0.42 132.9 57.5 27.6 14.9 9.1 6.1 4.3 3.3 2.6 2.1 1.8 1.6 1.4 120.3 41.3 17.9 9.5 5.9 4.0 3.0 2.3 1.9 1.6 1.4 1.3 1.2 116.3 37.7 16.2 8.6 5.4 3.7 2.7 2.2 1.8 1.6 1.4 1.3 1.2 113.1 34.4 15.0 8.2 5.1 3.5 2.6 2.1 1.8 1.5 1.4 1.3 1.2 0.62 88.7 43.1 22.5 13.3 8.4 5.7 4.2 3.2 2.5 2.1 1.8 1.6 1.4 72.7 30.0 14.6 8.5 5.5 3.8 2.9 2.3 1.9 1.6 1.4 1.3 1.2 68.6 27.5 13.3 7.7 5.0 3.5 2.7 2.1 1.8 1.5 1.4 1.3 1.2 67.2 25.3 12.5 7.2 4.7 3.4 2.6 2.1 1.7 1.5 1.4 1.3 1.2 0.82 55.0 30.4 17.9 11.3 7.5 5.4 4.0 3.1 2.5 2.1 1.8 1.6 1.4 41.2 20.8 11.7 7.3 5.0 3.6 2.7 2.2 1.8 1.6 1.4 1.3 1.2 39.0 18.9 10.7 6.7 4.5 3.3 2.6 2.1 1.7 1.5 1.4 1.3 1.2 38.1 17.4 9.9 6.3 4.3 3.2 2.5 2.0 1.7 1.5 1.4 1.3 1.2 τ2 1.02 33.6 21.6 14.0 9.5 6.7 4.9 3.7 3.0 2.4 2.0 1.8 1.6 1.4 24.6 14.4 9.1 6.2 4.4 3.4 2.6 2.1 1.8 1.6 1.4 1.3 1.2 22.7 13.1 8.4 5.7 4.1 3.1 2.5 2.0 1.7 1.5 1.4 1.3 1.2 22.1 12.3 7.8 5.4 3.9 3.0 2.4 2.0 1.7 1.5 1.3 1.2 1.2 1.22 21.5 15.2 10.8 7.8 5.8 4.5 3.5 2.9 2.3 2.0 1.7 1.5 1.4 15.2 10.1 7.2 5.2 3.9 3.1 2.5 2.1 1.8 1.6 1.4 1.3 1.2 13.9 9.4 6.5 4.8 3.6 2.9 2.3 1.9 1.7 1.5 1.4 1.3 1.2 13.8 8.7 6.2 4.5 3.5 2.8 2.2 1.9 1.6 1.5 1.3 1.2 1.2 1.42 14.4 11.0 8.4 6.5 5.1 4.1 3.3 2.7 2.3 1.9 1.7 1.5 1.4 10.0 7.5 5.7 4.4 3.5 2.8 2.3 2.0 1.7 1.5 1.4 1.3 1.2 9.2 6.9 5.2 4.1 3.2 2.6 2.2 1.9 1.6 1.5 1.3 1.3 1.2 8.6 6.5 5.0 3.9 3.1 2.5 2.1 1.8 1.6 1.4 1.3 1.2 1.2 1.62 9.9 8.2 6.7 5.5 4.5 3.7 3.0 2.5 2.2 1.9 1.7 1.5 1.4 6.9 5.6 4.6 3.7 3.1 2.6 2.2 1.9 1.7 1.5 1.4 1.3 1.2 6.5 5.2 4.2 3.5 2.9 2.4 2.1 1.8 1.6 1.4 1.3 1.2 1.2 6.1 4.9 4.0 3.3 2.8 2.4 2.0 1.8 1.6 1.4 1.3 1.2 1.2 1.82 7.2 6.3 5.4 4.6 3.9 3.3 2.8 2.4 2.1 1.8 1.6 1.5 1.4 5.1 4.4 3.8 3.2 2.8 2.4 2.1 1.8 1.6 1.5 1.4 1.3 1.2 4.7 4.1 3.5 3.0 2.6 2.2 2.0 1.7 1.6 1.4 1.3 1.2 1.2 4.5 3.9 3.3 2.9 2.5 2.2 1.9 1.7 1.5 1.4 1.3 1.2 1.2 2.02 5.5 5.0 4.4 3.9 3.4 3.0 2.6 2.3 2.0 1.8 1.6 1.5 1.3 3.9 3.5 3.1 2.8 2.5 2.2 1.9 1.7 1.6 1.4 1.3 1.2 1.2 3.6 3.2 2.9 2.6 2.3 2.1 1.8 1.7 1.5 1.4 1.3 1.2 1.2 3.4 3.1 2.8 2.5 2.2 2.0 1.8 1.6 1.5 1.4 1.3 1.2 1.2

Table 5.5: ARLs for shifts in θ with θj= (βjT, δσ)T and τ 2 j = τ 2 /δ2 for j ≥ 1 m0= 5 m0= 25 m0= 125 δ m0= ∞ 1.1 0 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18 1.19 1.110 1.111 1.112 2.22 4.3 4.0 3.7 3.4 3.0 2.7 2.4 2.2 1.9 1.7 1.6 1.4 1.3 3.1 2.9 2.7 2.4 2.2 2.0 1.8 1.7 1.5 1.4 1.3 1.2 1.2 2.9 2.7 2.5 2.3 2.1 1.9 1.7 1.6 1.5 1.4 1.3 1.2 1.2 2.8 2.6 2.4 2.2 2.0 1.9 1.7 1.6 1.4 1.3 1.3 1.2 1.1 2.42 3.5 3.4 3.1 2.9 2.7 2.5 2.2 2.0 1.8 1.7 1.5 1.4 1.3 2.6 2.5 2.3 2.1 2.0 1.9 1.7 1.6 1.5 1.4 1.3 1.2 1.2 2.4 2.3 2.2 2.0 1.9 1.8 1.6 1.5 1.4 1.3 1.3 1.2 1.2 2.3 2.2 2.1 2.0 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.2 1.1 2.62 2.9 2.8 2.7 2.6 2.4 2.3 2.1 1.9 1.8 1.6 1.5 1.4 1.3 2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.3 1.2 1.2 2.0 2.0 1.9 1.8 1.7 1.7 1.6 1.5 1.4 1.3 1.2 1.2 1.1 2.0 1.9 1.9 1.8 1.7 1.6 1.5 1.4 1.4 1.3 1.2 1.2 1.1 2.82 2.5 2.4 2.4 2.3 2.2 2.1 2.0 1.8 1.7 1.6 1.5 1.4 1.3 1.9 1.9 1.8 1.7 1.7 1.6 1.5 1.5 1.4 1.3 1.3 1.2 1.2 1.8 1.7 1.7 1.7 1.6 1.5 1.5 1.4 1.3 1.3 1.2 1.2 1.1 1.7 1.7 1.7 1.6 1.6 1.5 1.5 1.4 1.3 1.3 1.2 1.2 1.1 3.02 2.2 2.2 2.1 2.1 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.4 1.3 1.6 1.7 1.6 1.6 1.6 1.5 1.5 1.4 1.3 1.3 1.2 1.2 1.1 1.6 1.6 1.6 1.5 1.5 1.5 1.4 1.4 1.3 1.3 1.2 1.2 1.1 1.5 1.5 1.5 1.5 1.5 1.4 1.4 1.3 1.3 1.2 1.2 1.2 1.1 τ2 3.22 1.9 1.9 1.9 1.9 1.8 1.8 1.7 1.7 1.6 1.5 1.4 1.3 1.3 1.5 1.5 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.1 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.2 1.1 1.4 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.3 1.2 1.2 1.1 1.1 3.42 1.7 1.7 1.7 1.7 1.7 1.7 1.6 1.6 1.5 1.4 1.4 1.3 1.2 1.4 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.3 1.2 1.2 1.2 1.1 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.2 1.2 1.2 1.1 1.1 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.2 1.2 1.2 1.1 1.1 3.62 1.6 1.6 1.6 1.6 1.6 1.6 1.5 1.5 1.4 1.4 1.3 1.3 1.2 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.3 1.2 1.2 1.2 1.1 1.1 1.2 1.2 1.3 1.3 1.3 1.3 1.3 1.2 1.2 1.2 1.2 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.1 1.1 3.82 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.2 1.2 1.2 1.3 1.2 1.2 1.2 1.2 1.2 1.2 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.1 1.1 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.1 1.1 1.1 4.02 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.3 1.3 1.2 1.2 1.1 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1

Appendix

A.1

ψ(1) = −γ, ψ(1 2) = −γ − 2 log(2), ψ(x + 1) = ψ(x) + 1 x, ψ(x) = −γ + ∞ X n=1 x − 1 n(n + x − 1), ψ(x) = ln x − 1 2x− 1 12x2 + 1 12x4 − 1 252x6 + . . . , ψ′(1) = π 2 6 , ψ′(1 2) = π2 2 , ψ′(x + 1) = ψ′_{(x) −} 1 x2, ψ′(x) = ∞ X k=0 1 (x + k)2, ψ′_{(x) ∼} 1 x + 1 2x2 + 1 6x3 − 1 30x5 + 1 12x7 − 1 30x9 . . . ,

A.2

Eθj(W 2 j) = _σ2 j σ2 2    2 " p + 2(βj− β) T_XT jXj(βj− β) σ2 j # + " p + (βj− β) T_XT jXj(βj− β) σ2 j #2   +n2_j  log  _σ2 j njσ2  + 1 2 − 2njσ 2 j σ2  log  _σ2 j njσ2  + 1  (nj − p) +n2_j  log  _σ2 j njσ2  + 1   ψ nj − p 2  + log(2)  −2njσ 2 j σ2 (nj − p)  ψ nj− p + 2 2  + log(2)  + _σ2 j σ2 2 (nj− p)(2 + nj− p)2+ n2j ( ψ′ nj − p 2  +  ψ nj − p 2  + log(2) 2)

A.3

nj  lognj 2  − ψ nj₂− p  = nj  lognj 2  − log nj₂− p  + 1 nj − p + O 1 n2 j  = nj  log  nj nj − p  + 1 nj + O 1 n2 j  = nj  log  1 + p nj− p  + 1 nj + O 1 n2 j  = nj  p nj − p + 1 nj + O 1 n2 j   1 

nj  njψ′  nj − p 2  − 2  = nj  nj 2 nj− p + nj 4 2(nj − p)2 − 2 + O  1 n2 j  = nj  2 1 − p/nj + 2 nj(1 − p/nj)2 − 2 + O  1 n2 j  = nj  2 + 2p nj + 2 nj − 2 + O  1 n2 j  = 2(p + 1) + O 1 nj  .

A.4

Eθj(W ∗2 j ) = _σ2 j σ2 2    2 " p + 2(βj− β) T_XT jXj(βj− β) σ2 j # + " p + (βj− β) T_XT jXj(βj− β) σ2 j #2   +n2_j  log njσ 2 0 σ2 j  − 1 2 E(1{Hj≥aj}) + 2nj σ2 j σ2 0  log njσ 2 0 σ2 j  − 1  E(Hj · 1_{Hj≥aj}) −2n2 j  log njσ 2 0 σ2 j  − 1  E(log(Hj) · 1{Hj≥aj}) + _σ2 j σ2 2 E(Hj2· 1{Hj≥aj}) +n2_jE([log(Hj)]2· 1{Hj≥aj}) − 2nj σ2 j

where E(1{Hj≥aj}) = Z ∞ aj xv−1_e−x2 Γ (v) 2v dx = 1 − _{Γ (v) 2}1 _v ∞ X m=0 (−1)m_av+m j 2m_{· m!(v + m)}, with v _≡ nj− p 2 E(Hj · 1{Hj≥aj}) = Z ∞ aj xv_e−x 2 Γ (v) 2vdx = (nj − p) − 1 Γ (v) 2v ∞ X m=0 (−1)m_av+m+1 j 2m_{· m!(v + m + 1)}, E(log(Hj) · 1{Hj≥aj}) = Z ∞ aj log(x)xv−1e −x 2 Γ (v) 2v dx = [ψ (v + log(2))] − _{Γ (v) 2}1 v ∞ X m=0 (−1)m_am j 2m_{· m!} _log(a j)avj v + m − av j (v + m)2  , E(Hjlog(Hj) · 1{Hj≥aj}) =

Z ∞ aj log(x) x v_e−x 2 Γ (v) 2vdx = (nj − p)  ψ nj − p + 2 2 + log(2)  − _{Γ (v) 2}1 v ∞ X m=0 (−1)m_am j 2m_{· m!} G1(aj, v, m) with G1(aj, v, m) ≡ log(aj)av+1j v + m + 1 − av+1_j (v + m + 1)2 !

E(Hj2· 1{Hj≥aj}) = Z ∞ aj xv+1_e−x 2 Γ (v) 2v dx = (nj − p)(nj− p − 2) − 1 Γ (v) 2v ∞ X m=0 (−1)m_av+m+2 j 2m_{· m!(v + m + 2)} E([log(Hj)]2· 1{Hj≥aj}) = Z ∞ aj [log(x)]2x v−1_e−x 2 Γ (v) 2v dx = ψ′(v) + [ψ (v) + log(2)]2₋ 1 Γ (v) 2v ∞ X m=0 (−1)m_am j 2m_{· m!} G2(aj, v, m) with G2(aj, v, m) ≡ _[log(a j)]2avj(v + m) − log(aj)avj (v + m)2 − log(aj)avj(v + m)2− avj2(v + m) (v + m)4 

A.5

W0,j = (n0+ nj) log n0σˆ2+ njσˆj2+ ( ˆβj− ˆβ)T[(XT0X0)−1+ (XTjXj)−1]( ˆβj − ˆβ) (n0+ nj)σ2 ! −n0log  ˆσ2 σ2  − njlog _σˆ2 j σ2  a.s. → " n0(ˆσ2− σ2) σ2 + nj(ˆσj2− σ2) σ2 + ( ˆβj − ˆβ)T[(XT0X0)−1+ (XTjXj)−1]( ˆβj− ˆβ) σ2 # −n0  ˆσ2_{− σ}2 σ2  − njlog _ˆσ2 j σ2  + Op  1 √_n 0  a.s. → Wj+ Op  1 √_n 0 

as min_{n0, the minimum eigenvalue of XT0X0} → ∞.

A.6

W_0,j∗ _{≡ 2[ℓ}0,j(θ∗, θ∗j) − ℓ0,j( ˜θ, ˜θj)] = W0,j− 2[ℓ0,j( ˆθ, ˆθj) − ℓ0,j(θ∗, θ∗j)].

If the process is in control at time j (≥ 1), then

g(Bj) ≡ 2[ℓ0,j( ˆθ, ˆθj) − ℓ0,j(θ∗, θ∗j)]

= 1_{ˆσj<ˆσ0}



−n0log(ˆσ) − njlog(ˆσj) + (n0+ nj) log

_n

0ˆσ02+ njσˆj2

n0+ nj

 = 1_{Bj<bj}{−n0log(1 − Bj) − njlog(Bj) + [n0log(n0) + njlog(nj)

where Bj ≡ Hj/(H0 + Hj) ∼ beta nj₂−p,n0₂−p
, H0 ≡ n0σˆ2/σ2 ∼ χ2n0−p, g(Bj) is

a function of Bj, bj ≡ nj/(n0 + nj), and 1{Bj<bj} denotes the indicator function for

{Bj < bj}. Eθj(W0,j∗ ) = Eθj(W0,j) − Eθj(g(Bj)), where

Eθj(g(Bj)) = [n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)]E(1{Bj<bj})

−n0E(log(1− Bj) · 1_{Bj<bj}) − njE(log(Bj) · 1{Bj<bj})

Varθj(W

∗

0,j) = Varθj(W0,j−g(Bj)) = Varθj(W0,j)+Varθj(g(Bj))−2Covθj(W0,j, g(Bj)),

where Varθj(g(Bj)) = Eθj(g(Bj) 2 ) − [Eθj(g(Bj))] 2 and Covθj(W0,j, g(Bj))

= n2_{0Cov(log(1 − B}j), log(1 − Bj) · 1{Bj<bj}) + n0njCov(log(1 − Bj), log(Bj) · 1{Bj<bj})

−n0[n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)]Cov(log(1 − Bj), 1{Bj<bj})

+n0njCov(log(Bj), log(1 − Bj) · 1{Bj<bj}) + n

2

jCov(log(Bj), log(Bj) · 1{Bj<bj})

−nj[n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)]Cov(log(Bj), 1{Bj<bj})

with

Eθj(g(Bj)

2

) = E({n20[log(1 − Bj)]2+ n2j[log(Bj)]2 + 2n0njlog(1 − Bj) log(Bj)

+[n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)]2

−2n0[n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)] log(1 − Bj)

Let α ≡ (nj− p)/2 and β ≡ (n0− p)/2 E(1{Bj<bj}) = Z bj 0 Γ(α + β)xα−1_{(1 − x)}β−1 Γ(α)Γ(β) dx = Γ(α + β) Γ(α)Γ(β) ∞ X m=0 (−1)mC_mβ−1 b α+m j α + m E(log(Bj) · 1{Bj<bj}) = Z bj 0 log(x)Γ(α + β)xα−1(1 − x)β−1 Γ(α)Γ(β) dx = ∂ ∂αE(1{Bj<bj}) E(log(1− Bj) · 1{Bj<bj}) = Z bj 0 log(1 − x) Γ(α + β)xα−1_{(1 − x)}β−1 Γ(α)Γ(β) dx = ∂ ∂βE(1{Bj<bj}) E(log(Bj) log(1 − Bj) · 1{Bj<bj}) = Z bj 0 log(x) log(1 − x) Γ(α + β)xα−1_{(1 − x)}β−1 Γ(α)Γ(β) dx = ∂2 ∂α∂βE(1{Bj<bj}) E([log(Bj)]2· 1{Bj<bj}) = Z bj 0 [log(x)]2Γ(α + β)xα−1(1 − x)β−1 Γ(α)Γ(β) dx = ∂2 ∂α2E(1{Bj<bj}) E([log(1− Bj)]2· 1{Bj<bj}) = Z bj 0 [log(1 − x)] 2Γ(α + β)xα−1(1 − x)β−1 Γ(α)Γ(β) dx = ∂2 ∂β2E(1{Bj<bj})