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 1/ 2
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ˆ 0 is about 200, then the standard error of this estimator is about
1 1/ 2
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 Y = +I0 M e0 −N0(x−1)2+ε , where ε ~N(0,1). The in-control reference profile is 1 15+ e− −(x1)2, which is displayed in Fiugre 3. We choose
x -values of 0, 0.08,…, 3.92 (n = 50) in our simulation. Four different types of shifts i
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 .
( 1)2
1 15+ e− −x
( 1)2
1 18+ e− −x
1.5( 1)2
1 15+ e− x−
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 in ( ( )f xi bl c Bˆl l k( ))i n
∑
= −∑
= x .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 yij = +3 2xi+εij , where with fixed
. . .
~ (0,1
i i d
ij N
ε )
x -values of 2, 4, 6, and 8. Three different types of shifts, intercept shifts, slope i
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.
Table 16 and Figure 21 give the ARL values for standard deviation shifts. The
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.