Based on the standardized LR test statistics, the EWMA sequence is defined as
Uj ≡ λ ˜Wj+ (1 − λ)Uj−1, j = 1, 2, ..., (2.11)
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
Uj∗ ≡ 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+k2Covθ( ˜Wj−k1, ˜Wj−k2) (2.15)
for unknown θ with Covθ( ˜Wj−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 Covθ( ˜Wj−k1, ˜Wj−k2) is approximated by a Monte Carlo estimate.
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 W1∗, . . . , W200∗ , and we want to standardized the Wj∗, so we evaluate the Eθ(Wj∗) and Varθ(Wj∗), which are used to evaluate ˜Wj∗. 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 Uj∗ 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 {U1∗, . . . , U200∗ } (≡ Q(1)), repeat above steps for 50,000 times
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. W1∗ 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.
Wj∗, 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 W0,1∗ , which use to evaluate Uj∗ 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 {U1∗, . . . , 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 (≡ q1), second sample quartile (≡ q2) and third sample quartile (≡ q3).
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. W0,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.
W0,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
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 − β)TXTjXj(β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 wj∗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.
[7] E.S. Page. Continuous inspection schemes. Biometrika, 41:100–115, 1954.
[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.
[16] C. Zou, C. Zhou, Z. Wang, and F. Tsung. A self-starting control chart for linear
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 τj2 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 τj2= 120 δ21for 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 τj2= τ2/δ2for j ≥ 1
Table 5.5: ARLs for shifts in θ with θj= (βjT, δσ)T and τj2= τ2/δ2for j ≥ 1
Appendix
A.2
nj
where
E(Hj2· 1{Hj≥aj}) = Z ∞
aj
xv+1e−x2 Γ (v) 2v dx
= (nj − p)(nj− p − 2) − 1 Γ (v) 2v
X∞ m=0
(−1)mav+m+2j 2m· m!(v + m + 2) E([log(Hj)]2· 1{Hj≥aj}) =
Z ∞
aj
[log(x)]2xv−1e−x2 Γ (v) 2v dx
= ψ′(v) + [ψ (v) + log(2)]2− 1 Γ (v) 2v
X∞ m=0
(−1)mamj
2m· m! G2(aj, v, m)
with
G2(aj, v, m) ≡
[log(aj)]2avj(v + m) − log(aj)avj
(v + m)2 − log(aj)avj(v + m)2− avj2(v + m) (v + m)4
A.5
where Bj ≡ Hj/(H0 + Hj) ∼ beta nj2−p,n02−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(W0,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))
= n20Cov(log(1 − Bj), 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}) + n2jCov(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)
−2nj[n0log(n0) + njlog(nj) − (n0+ nj) log(n0+ nj)] log(Bj)} · 1{Bj<bj}).
Let α ≡ (nj− p)/2 and β ≡ (n0− p)/2