A Computer Program to Calculate Two-Stage Short-Run Control Chart Factors for ( X, MR ) Charts

This paper is the second in a series of two papers that fully develops two-stage short-run ( X, MR ) control charts. This paper describes the development and execution of a computer program that accurately calculates ﬁrst- and second-stage short-run control chart factors for ( X, MR ) charts using the equations derived in the ﬁrst paper. The software used is Mathcad . The program accepts values for number of subgroups, α for the X chart, and α for the MR chart both above the upper control limit and below the lower control limit. Tables are generated for speciﬁc values of these inputs and the implications of the results are discussed. A numerical example illustrates the use of the program.


Introduction
Hillier's two-stage short-run theory of control charting has been applied to (X, R) charts (see Hillier 1969), (X, v) and (X, √ v) charts (see Yang and Hillier 1970), (X, s) charts (see Elam and Case 2005a), and (X, MR) charts (see Elam and Case 2006), where R, v, s, and MR are the range, variance, standard deviation, and moving range (size two) of a subgroup, respectively. Elam and Case (2001) and Elam and Case (2003) develop computer programs to address issues concerning the limited tabulated results provided by Hillier (1969) and Yang and Hillier (1970), respectively. Elam and Case (2005b) develop a computer program using the equations in Elam and Case (2005a) to calculate the factors required to determine twostage short-run control limits for (X, s) charts so that the application of these control charts would not be limited. This paper describes the development and execution of a computer program that uses the equations in Elam and Case (2006) to accurately calculate first-and second-stage short-run control chart factors for (X, MR) charts regardless of the number of subgroups, α for the X chart, and α for the MR chart both above the upper control limit and below the lower control limit (α is the probability of a false alarm). The software used is Mathcad (Mathsoft 2005). The computer program also uses functions provided by the software. This paper first presents the computer program. Tables generated by the program are then presented and compared with legitimate results in the literature. Also, implications of the tabulated results are discussed. Following a numerical example that illustrates the use of the program, final conclusions describing the contributions of the program are given. It should be noted that results from the program are for processes generating parts with independent measurements that follow a Normal distribution.

The computer program
This section of the paper presents the computer program, which is available at http:// www.jstatsoft.org/v15/i11/ or starting at http://program.20m.com/. The program has seven pages, each of which is further divided into sections. When executing the program, it is possible for a section of code in the program to turn red and have the error message "Unknown Error." To correct this, delete one character in the red code and type it back in. Click the mouse arrow outside of the code. The code should turn black, indicating that the error has been eliminated. If not, repeat the procedure (it will eventually correct the problem).

Page 1
The first page of the program begins with the data entry section. The program requires the user to enter the following values: alphaInd (α for the X chart), alphaMRUCL (α for the MR chart above the UCL), alphaMRLCL (α for the MR chart below the LCL), and m (number of subgroups (i.e., the number of MRs plus one)). If no lower control limit on the MR chart is desired, the entry for alphaMRLCL should be left blank (do not enter zero). Before a value can be entered, the cursor must be moved to the right side of the appropriate equal sign. This may be done using the arrow keys on the keyboard or by moving the mouse arrow to the right side of the equal sign and clicking once with the left mouse button. The program is activated by paging down once the last entry is made. Mathcad allows the user to immediately page down to the output section of the program (explained later in subsection 2.7) after the last entry is made. The next part of page 1 is section 1.1 of the program. The value TOL is the tolerance. The calculations that use this value will be accurate to ten places to the right of the decimal. The functions dnorm(x, 0, 1) and pnorm(x, 0, 1) in Mathcad are the pdf and cdf, respectively, of the standard Normal distribution. The equations for the pdf and cdf are also given in case the dnorm or pnorm function fails to calculate a result. In Mathcad, an equation turns red if it does not calculate a result due to an error. If the dnorm function gives an error, type f(x) on the left side of the equal sign in Equation 1: If the pnorm function gives an error, type F(x) on the left side of the equal sign in Equation 2: The equation for P(W), the probability integral (or cumulative distribution function (cdf)) of the range for subgroups of size two sampled from a standard Normal population, is Equation 1 in Elam and Case (2006). As stated in Elam and Case (2006), W represents the (standardized) range w/σ, where w is the range of a subgroup and σ is the population standard deviation. The equation for d2, the mean of the distribution of the range W = w/σ for subgroups of size two sampled from a Normal population with mean µ and variance equal to one, is Equation 2 in Elam and Case (2006).

Page 2
Page 2 of the program begins with section 2.1. The code in this section determines wD4 and wD3, the 1-alphaMRUCL and alphaMRLCL percentage points, respectively, of the distribution of the range W = w/σ for subgroup size two. The values wD4 and wD3 are used to determine the α-based conventional upper and lower control chart constants, respectively, for the MR chart. The roots of the equations DUCL(W) and DLCL(W) are wD4 and wD3, respectively, and are determined using root (a function in Mathcad that uses Ridder's or Brent's method to find the roots of an equation when root bracketing is used). The subprograms Wseed1 and Wseed2 generate seed values seedD4 and seedD3, respectively, for Ridder's or Brent's method.
The subprogram Wseed1 works as follows. Initially, W 0 and W 1 are set equal to 0.01 and 0.02, respectively. A 0 and A 1 result from evaluating DUCL(W) at W 0 and W 1 , respectively. The while loop begins by checking if the product of A 0 and A 1 is negative. If so, the root for DUCL(W) lies between 0.01 and 0.02. If not, W 0 and W 1 are incremented by 0.01. A 0 and A 1 are recalculated and if their product is negative, the root for DUCL(W) lies between 0.02 and 0.03. Otherwise, the while loop repeats. Once a root for DUCL(W) is bracketed, the bracketing values are passed out of the subprogram into the 2 × 1 vector seedD4 to be used by Ridder's or Brent's method to determine wD4. The subprogram Wseed2 works similarly to construct the 2 × 1 vector seedD3 to be used by Ridder's or Brent's method to determine wD3, except the starting value is 0.001.
The next part of page 2 is section 2.2 of the program. The following is a description of the equations in section 2.2 of the program and their corresponding equation numbers in Elam and Case (2006): • d(x): the equation for determining ν, Equation 5 • h(x): the ratio of the variance to the squared mean, both of the χ distribution with x degrees of freedom, Equation 6 (h(x) has been updated for Mathcad, hence "lnΓ" replacing "gammln") • r: the ratio of the variance to the squared mean, both of the distribution of the mean moving rangeMR/σ, Equation 7a • b and c: Equations 7b and 7c, respectively The values rprevm and dprevm(x) have the same meanings as r and d(x), respectively, except they are for m − 1 subgroups. The value ν is the root of the equation d(x) and is determined using root with seed value seedν. The value νprevm is the root of the equation dprevm(x) and is determined using root with seed value seedνprevm. The subprogram dfseed generates the seed values seedν and seedνprevm for Ridder's or Brent's method.
The subprogram dfseed works as follows. Initially, df 0 and df 1 are set equal to 0.9 and 1.1, respectively. A 0 and A 1 result from evaluating y(x) (which is equal to either d(x) or dprevm(x)) at df 0 and df 1 , respectively. The while loop begins by checking if the product of A 0 and A 1 is negative. If so, the root for y(x) lies between 0.9 and 1.1. If not, df 0 and df 1 are incremented by 0.5. A 0 and A 1 are recalculated and if their product is negative, the root for y(x) lies between 1.1 and 1.6. Otherwise, the while loop repeats. Once a root for y(x) is bracketed, the bracketing values are passed out of the subprogram into the 2 × 1 vector is equal to dprevm(x)) to be used by Ridder's or Brent's method to determine ν or νprevm, respectively. Both ν and νprevm are used to determine the two-stage short-run control chart factors for the (X, MR) charts.

Page 3
Page 3 of the program begins with section 3.1. The following is a description of the equations in section 3.1 of the program and their corresponding equation numbers in Elam and Case (2006): • P3(z), the probability integral of the studentized range for subgroups of size two sampled from a Normal population, Equation 3a • cν, P1(z), and P2(z), components of P3(z), Equations 3b, 3c, and 3d, respectively (cν has been updated for Mathcad, hence "lnΓ" replacing "gammln") As stated in Elam and Case (2006), the variable z is equal to 5 × Q. Q represents the studentized range w/s, where w is the range of a subgroup and s is an independent estimate (based on ν degrees of freedom) of the population standard deviation. The function lnΓ in the equation for cν calculates the natural logarithm of the gamma (Γ) function.
Section 3.2 contains the calculations required to determine qD4, the 1-alphaMRUCL percentage point of the distribution of the studentized range Q = w/s for subgroup size two with ν degrees of freedom. The value qD4 is used to determine the second-stage short-run upper control chart factor for the MR chart. The subprogram Zseed1 generates the seed value seed1 for Ridder's or Brent's method (the first root function) or for the Secant or Mueller's method (the second root function -a function in Mathcad that finds the roots of an equation when root bracketing is not used). Either root-finding method determines the root of D(x). The result of dividing this root by five is qD4. Both root functions are given because one may not work when the other one does. If the first root function fails, type qD4 on the left side of the equal sign in The subprogram Zseed1 begins by generating values for Z 0 and Z 1 . A 0 and A 1 result from evaluating P3(z) at Z 0 and Z 1 , respectively. The while loop continually increments Z 0 and Z 1 by 5.0 and evaluates P3(z) at these two values until A 1 becomes greater than 1-alphaMRUCL for the first time, at which point A 0 will be less than 1-alphaMRUCL. When this occurs, P3(z) is equal to 1-alphaMRUCL for some value z between Z 0 and Z 1 . An initial guess for this value is determined using linterp (a function in Mathcad that performs linear interpolation) and stored in Zguess. The initial guess is passed out of the subprogram as seed1.

Page 4
Page 4 of the program is section 4.1. The code in this section is used to determine qD3, the alphaMRLCL percentage point of the distribution of the studentized range Q = w/s for subgroup size two with ν degrees of freedom. The value qD3 is used to determine the secondstage short-run lower control chart factor for the MR chart. The subprogram Zseed2 generates the value seed2 that is used to determine an initial value for qD3. An improved value for qD3 is then calculated by determining the root of the equation P3(z)-alphaMRLCL using the Secant or Mueller's method with the seed value seed2 and dividing this root by five.
The ability of the Secant or Mueller's method to calculate a result depends upon the values for alphaMRLCL and m (Ridder's or Brent's method should not be used). It is not a problem if it does not calculate a result because the initial value for qD3 and the improved value match to several places to the right of the decimal. This phenomenon is discussed in more detail when the tabulated results of the program are presented later in section 3. The Monitor Results area in the bottom right hand corner of section 4.1 indicates how closely the two values for qD3 match until the root routine fails. This will dictate the number of decimal places that can be used to display qD3 and the second-stage short-run lower control chart factor for the MR chart.
The code in the subprogram Zseed2 that begins with the first line of code and includes the while loop and the two for loops constructs 21 × 1 vectors Zv for z and Av for P3(z). The first row of each vector is zero. The while loop determines the first value Z where P3(Z) is greater than alphaMRLCL. This Z and the corresponding value P3(Z) are stored in the second rows of Zv and Av, respectively. The two for loops generate values for the remaining rows of Zv and Av. Two different for loops are used because P3(z) may encounter an error for some i (i: 1, 2, . . . , 20). The value for i where the error occurs can be skipped using the dual for loop construction. When the execution of this section of code is complete, P3(z) is equal to alphaMRLCL for some value z between Zv 0 and Zv 1 .
The code in the subprogram Zseed2 that starts in the line where the variable Zguess first appears to the last line of the subprogram is derived from Harter, Clemm, and Guthrie (1959). This code searches for and estimates the value z where P3(z) is equal to alphaMRLCL. Zguess is the initial guess for this value z. It is determined using linterp, the 21×1 vectors for P3(z) and z previously determined, and alphaMRLCL. Aguess, the estimated value for P3(Zguess), is determined using the Mathcad functions lspline and interp. The function lspline returns a vector vs which interp uses to create a cubic, piecewise polynomial that passes through all of the (x, y) data points, and assumes the resultant spline curve is linear at the endpoints. The function interp returns a spline interpolated value of Av at a point Zguess and stores the result in Aguess. The while loop first checks if Aguess is an accurate estimate (within 10 −15 ) of alphaMRLCL. If so, Zguess is passed out of the subprogram as the value seed2. If not, Aguess and Zguess are entered into the second rows of the previously determined vectors Av and Zv, respectively, if Aguess is more than 10 −15 larger than alphaMRLCL. If Aguess is more than 10 −15 smaller than alphaMRLCL, Aguess and Zguess are entered into the first rows of the vectors Av and Zv, respectively. New values for Zguess and Aguess are determined using the same procedure as before and execution is returned to the beginning of the while loop.

Page 5
Page 5 of the program contains sections 5.1 and 5.2. These sections correspond to sections 3.1 and 3.2, respectively, described earlier in subsection 2.3. The only difference is that the calculations in sections 5.1 and 5.2 use νprevm instead of ν. The calculations are for qD4prevm, the 1-alphaMRUCL percentage point of the distribution of the studentized range Q = w/s for subgroup size two with νprevm degrees of freedom. The value qD4prevm is used to determine the first-stage short-run upper control chart factor for the MR chart.

Page 6
Page 6 of the program is section 6.1. This section corresponds to section 4.1 described earlier in subsection 2.4. The only difference is that the calculations in section 6.1 use νprevm instead of ν. The calculations are for qD3prevm, the alphaMRLCL percentage point of the distribution of the studentized range Q = w/s for subgroup size two with νprevm degrees of freedom. The value qD3prevm is used to determine the first-stage short-run lower control chart factor for the MR chart.

Page 7
Page 7 of the program begins with section 7.1. The value d2starMR is Equation 4 in Elam and Case (2006). The value d2starMRprevm has the same meaning as d2starMR, except it is for m − 1 subgroups. Both d2starMR and d2starMRprevm are used to determine the two-stage short-run control chart factors for the (X, MR) charts. The function qt(adj_alpha, ν) in Mathcad determines the critical value crit_t for a cumulative area of adj_alpha under the Student's t curve with ν degrees of freedom. The value crit_t is used to determine first-and second-stage short-run control chart factors for the X chart. The function norm(adj_alpha, 0, 1) in Mathcad determines the critical value crit_z for a cumulative area of adj_alpha under the standard Normal curve. The value crit_z is used to determine the conventional control chart constant for the X chart.
The next part of page 7 is section 7.2 of the program. The following is a description of the equations in section 7.2 of the program and their corresponding equation numbers in Elam and Case (2006): • E21: the first-stage short-run control chart factor for the X chart, Equation 9 • E22: the second-stage short-run control chart factor for the X chart, Equation 8 • E2: the conventional control chart constant for the X chart, Equation 14 (the equation for E2 is a generalization of the equation for E 2 from Wheeler (1995)'s Tables 3 and 4 to allow for different values of alphaInd) • D41: the first-stage short-run upper control chart factor for the MR chart, Equation 12 • D42: the second-stage short-run upper control chart factor for the MR chart, Equation 10 • D4: the α-based conventional upper control chart constant for the MR chart, Equation 15 • D31: the first-stage short-run lower control chart factor for the MR chart, Equation 13 • D32: the second-stage short-run lower control chart factor for the MR chart, Equation 11 • D3: the α-based conventional lower control chart constant for the MR chart, Equation 16 The last part of page 7 is the output section of the program. The four values entered at the beginning of the program are given. The control chart factors are broken down into first-stage, second-stage, and conventional. The values for ν, d2starMR, νprevm, and d2starMRprevm, the mean of the distribution of the range W = w/σ for subgroup size two and the variance of the distribution of the mean moving rangeMR/σ, and Harter, Clemm, and Guthrie (1959)'s Table  II.2 results for n = 2 (i.e., for subgroup size two) complete the output of the program. To copy results into another software package (like Excel), follow the directions from Mathcad's help menu or highlight a value and copy and paste it into the other software package. When highlighting a value with the mouse arrow, place the arrow in the middle of the value, depress the left mouse button, and drag the arrow to the right. This will ensure just the numerical value of the result is copied and pasted.
The values qD4, qD4prevm, and wD4, as well as qD3, qD3prevm, and wD3, are in  Table 2, the Secant or Mueller's method fails to work for m ≥ 3. As mentioned previously, this is not a serious issue. The reason is that the initial value for qD3 matches the improved value for qD3 (before the Secant or Mueller's method fails) to eight places to the right of the decimal.
Values for E21, D41, D31, E22, D42, D32, E2, D4, and D3 are in Table 3. The E2 value compares favorably to the E 2 value for n = 2 in Wheeler (1995)'s Table 4. It should be noted that the values wD4 and wD3 in Table 2 and D4 and D3 in Table 3 may differ in the ninth or tenth decimal place for different root routines used to calculate wD4 and wD3.
These favorable comparisons validate the program. Consequently, Table 3 Table 3 show some interesting properties. As m increases, the E22 and D42 values converge in a decreasing manner to E2 and D4, respectively. The D32 values also converge in a decreasing manner to D3, though it is not evident from the accuracy shown. This convergence makes sense because more information about the process is at hand for larger m.

Values in
These results have major implications. A common rule of thumb is that 20 to 30 subgroups of size 4 or 5 (i.e., a combined sample of between 80 and 150 data values) are necessary to use conventional control chart constants for constructing control limits. The results in Table  3 indicate that this may be an incorrect rule when applied to constructing (X, MR) control charts. Consider Figure 1, which is a plot of E22 and E2 for alphaInd=0.0027 and m: 5, 10,15,20,25,30,50,75,100,150. It indicates that if one were to construct X charts using the conventional control chart constant E2 when 150 or less subgroups of size one (i.e., a combined sample of 150 or less data values) are available to estimate the process mean and standard deviation, the upper and lower control limits would not be wide enough, resulting in a higher probability of a false alarm.  Table 3 also indicate that the common rule of thumb, when applied to constructing (X, MR) control charts, may be an incorrect rule. Consider Figure 2, which is a plot of D42 and D4 for alphaMRUCL=0. 005 and m: 5, 10, 15, 20, 25, 30, 50, 75, 100, 150. It indicates that if one were to construct the upper control limit of MR charts using the conventional control chart constant D4 when 150 or less subgroups of size one are available to estimate the process standard deviation, the upper control limit would not be wide enough, resulting in a higher probability of a false alarm.
To the accuracy shown in Table 3, there is no difference between D32 for any m and D3. If increased accuracy is used, then D3 is slightly less than D32 for any m. Consequently, if one were to construct the lower control limit of MR charts using the conventional control chart constant D3 when 150 or less subgroups of size one are available to estimate the process standard deviation, the lower control limit would be slightly too wide, possibly creating a situation in which the probability of detecting a special cause signal is slightly diminished. Quesenberry (1993) also investigated the validity of the common rule of thumb when applied to constructing (X, MR) control charts and concluded that 300 individual values are needed for the X chart before conventional control chart constants may be used. However, for all practical purposes, the computer program presented by this paper eliminates the need for

A numerical example
Consider the data in Table 4 obtained from a process requiring short-run control charting techniques (assume alphaInd=0.0027, alphaMRUCL=0.005, and alphaMRLCL=0.001).  For m=5, the following first-stage short-run control chart factors for the MR chart are obtained from Table 3 The first moving range (MR=0.151) is above UCL(MR). Find, investigate, and remove from the process the special cause (or causes) that created this out of control point, delete the first moving range, recalculate the average moving range (shown as the Revised Average in Table 4), and construct second-stage control limits for the (X, MR) charts (this approach is from Case (1998)). For m=4, the following second-stage short-run control chart factors for the MR chart are obtained from Table 3: D42=13.20218 and D32=0.00157. For m=5, the following second-stage short-run control chart factor for the X chart is obtained from Table  3 These control limits may be used to monitor the future performance of the process.

Conclusions
This paper is the second in a series of two papers that fully develops two-stage short-run (X, MR) control charts. This paper described the development and execution of a computer program that accurately calculates first-and second-stage short-run control chart factors for (X, MR) charts using the equations derived in the first paper. The computer program makes two important contributions. Those involved with quality control in industry will, for the first time, be able to use theoretically precise control chart factors to determine control limits for (X, MR) charts regardless of the number of subgroups and α values. Also, as already mentioned, the computer program eliminates the need for the research question of how many subgroups are enough before conventional control chart constants may be used.