A Study on the X ¯ and S Control Charts with Unequal Sample Sizes

Abstract

The control charts based on $\overline{X}$ and S are widely used to monitor the mean and variability of variables and can help quality engineers identify and investigate causes of the process variation. The usual requirement behind these control charts is that the sample sizes from the process are all equal, whereas this requirement may not be satisfied in practice due to missing observations, cost constraints, etc. To deal with this situation, several conventional methods were proposed. However, some methods based on weighted average approaches and an average sample size often result in degraded performance of the control charts because the adopted estimators are biased towards underestimating the true population parameters. These observations motivate us to investigate the existing methods with rigorous proofs and we provide a guideline to practitioners for the best selection to construct the $\overline{X}$ and S control charts when the sample sizes are not equal.

1. Introduction

Control charts, also known as Shewhart control charts [1,2,3], have been used to determine if a manufacturing process is in a state of control. In particular, the $\overline{X}$ and S charts have been widely used to monitor or detect the mean and variability of variables. Here, a variable is a quality characteristic measured on a numerical scale. For example, variables include continuous measurement process data such as length, pressure, width, temperature, and volume, in a time-ordered sequence.

Due to their importance and usefulness in real life applications, these traditional types of univariate and control charts have still received much attention in the literature. We observe that these control charts are usually adopted for continuously monitoring numerous data and for solving the problem of process control in the Industry 4.0 framework; see, for example [4,5,6], among others. Of particular note is that [7] developed a nice qcr package in R to generate Shewhart-type charts and obtained numerical results of interest for a process quality control. More recently, based on the concept of data depth, Ref. [8] proposed a novel alternative way for constructing control charts when the critical to quality (CTQ) variables of the process are functional and also developed the Phase I and II control charts for stabilizing and monitoring the processes, respectively.

It deserves mentioning that these traditional Shewhart control charts mentioned above consist of the upper and lower control limits (for short, UCL and LCL) and the center line (CL). It is noteworthy that the American Standard is based on $\mathrm{CL}\pm 3·\mathrm{SE}$ control limits with an ideal false alarm rate (FAR) of 0.27% while the British Standard is based on $\mathrm{CL}\pm 3.09·\mathrm{SE}$ with an ideal FAR of 0.20%, where $\mathrm{SE}$ denotes the standard error.

The usual requirement behind these control charts is that the sample sizes from the process are all equal. In practice, however, it is often the case that this requirement can not be met due to wrong or missing observations in collecting them. In such setting, the three conventional approaches below are widely used to deal with unequal sample sizes:

(i)

A weighted average approach in calculating $\overline{\overline{X}}$ and ${\overline{S}}^2$.

(ii)

Control limits based on an average sample size.

(iii)

A weighted average in calculating $\overline{S}$.

For more details, see Subsection 6.3.2 of Montgomery [9] and Subsection 3.8.B of ASTM [10]. The first approach uses variable-width control limits which are determined by the sample-specific values such as $n_i$, $A_3$, $B_3$, and $B_4$. To estimate the scale parameter, a weighted average of sample variances is calculated first and then its square-root is taken to estimate the scale parameter. The second approach uses fixed-width control limits which is based on the average of the sample sizes. For more details on these two methods, see Subsection 6.3.2 of Montgomery [9]. The third approach is very similar to calculating ${\overline{S}}^2$ in the first approach. However, it uses a weighted average of sample standard deviations directly to estimate the scale parameter. For more details, see Subsection 3.8.B of ASTM [10]. It is known that these three approaches may be satisfactory when the sample sizes are not very different. Given that the average of the sample sizes is not necessarily an integer in general, a practical alternative to the second approach is the use of a modal sample size.

However, when using these ad hoc approaches above, the parameter estimators are biased and they actually underestimate the true population parameters as will be shown in Remark 2 in Section 2 and Remarks 3 and 4 in Section 3. These underestimating ad hoc approaches could result in degraded performance of the control charts. Nonetheless, these biased methods are widely covered in many popular textbooks. These observations motivate us to clarify these conventional methods and investigate other existing methods, especially when the samples are not equal in size. Through the rigorous proofs, we provide a guideline for the best selection of the methods to construct the $\overline{X}$ and S control charts.

This paper is organized as follows. In Section 2, we provide two location estimators and four scale estimators with unequal sample sizes and show that they are all unbiased. In Section 3, we provide the variances of the estimators considered in this paper and show the inequality relations among them. In Section 4, we provide the relative efficiency of the methods and conduct simulation results to compare the performance of the location and scale estimators. In Section 5, we illustrate how to construct various Shewhart-type control charts (i.e., S, $S^2$, and $\overline{X}$ charts) using the provided estimators. In Section 6, we provide the empirical estimates of the average run length (ARL) and the standard deviation of the run length (SDRL) through using the extensive Monte Carlo simulations. Three real-data examples are presented in Section 7 for illustrative purposes. Some concluding remarks are given in Section 8.

2. Estimation of Process Parameters with Unequal Sample Sizes

In this section, we provide two location estimators and three scale estimators for the process parameters under the assumption that each sample has different sample sizes. In parametric statistical quality control, the underlying distribution is used to construct the control charts. A quality characteristic is assumed to be normally distributed, which is most widely used in most practical cases. Under this assumption, we show that the estimators provided in this section are all unbiased.

We assume that we have m samples and that each sample has different sample sizes. Let $X_{ij}$ be the ith sample (subgroup) of size $n_i$ from a stable manufacturing process, where $i=1,2,\dots ,m$ and $j=1,2,\dots ,n_i$. We also assume that $X_{ij}$ are independent and identically distributed as normal with mean $\mu$ and variance ${\sigma }^2$.

2.1. Location Parameter

Using the ith sample above, the sample mean and the sample variance are given by

$${\overline{X}}_i=\frac{1}{n_i}\sum _{j=1}^{n_i}{X}_{ij},\mathrm{and}{S}_i^2=\frac{1}{n_i-1}\sum _{j=1}^{n_i}{(X_{ij}-{\overline{X}}_i)}^2,$$

where $i=1,2,\dots ,m$. Montgomery [9] provides two location estimators of the population mean parameter $\mu$ in Equations (6.2) and (6.30) in his book, which are given by

$$\begin{array}{cc}\hfill {\overline{\overline{X}}}_A\hfill & =\frac{{\overline{X}}_1+{\overline{X}}_2+\cdots +{\overline{X}}_m}{m}=\frac{1}{m}\sum _{i=1}^m{\overline{X}}_i\hfill \end{array}$$

(1)

and

$$\begin{array}{cc}\hfill {\overline{\overline{X}}}_B\hfill & =\frac{n_1{\overline{X}}_1+n_2{\overline{X}}_2+\cdots +n_m{\overline{X}}_m}{n_1+n_2+\dots +n_m}=\frac{1}{N}\sum _{i=1}^m{n}_i{\overline{X}}_i,\hfill \end{array}$$

(2)

where $n_i\ge 2$ and $N={\sum }_{i=1}^m{n}_i$.

These grand averages can be used as the CL on the $\overline{X}$ chart. Since $E\left({\overline{X}}_i\right)=\mu$ for $i=1,2,\dots ,m$, it is easily seen that $E\left({\overline{\overline{X}}}_A\right)=\mu$ and $E\left({\overline{\overline{X}}}_B\right)=\mu$, showing that these two estimators are unbiased. In addition, the variances of ${\overline{\overline{X}}}_A$ and ${\overline{\overline{X}}}_B$ are given by $\mathrm{Var}\left({\overline{\overline{X}}}_A\right)={\sigma }^2{\sum }_{i=1}^m{n}_i^{-1}/m^2$ and $\mathrm{Var}\left({\overline{\overline{X}}}_B\right)={\sigma }^2/N$, which results in $\mathrm{Var}\left({\overline{\overline{X}}}_A\right)\ge \mathrm{Var}\left({\overline{\overline{X}}}_B\right)$. Thus, ${\overline{\overline{X}}}_B$ is preferred to ${\overline{\overline{X}}}_A$. It should be noted that, for the case of an equal sample size, we have ${\overline{\overline{X}}}_A={\overline{\overline{X}}}_B=\overline{\overline{X}}$, where $\overline{\overline{X}}={\sum }_{i=1}^m{\overline{X}}_i/m$.

2.2. Scale Parameter

It is well known that $S_i^2$ is an unbiased estimator of ${\sigma }^2$. However, $S_i$ is not an unbiased estimator of $\sigma$ as below. Since $(n_i-1)S_i^2/{\sigma }^2$ is distributed as the gamma with shape $(n_i-1)/2$ and scale 2, we have

$$E\left[{\left\{\frac{(n_i-1)S_i^2}{{\sigma }^2}\right\}}^{1/2}\right]=\frac{2^{1/2}·\mathsf{\Gamma }(n_i/2)}{\mathsf{\Gamma }((n_i-1)/2)}.$$

Then, we have

$$E\left[S_i\right]=\sqrt{\frac{2}{n_i-1}}\frac{\mathsf{\Gamma }(n_i/2)}{\mathsf{\Gamma }((n_i-1)/2)}·\sigma =c_4\left(n_i\right)·\sigma ,$$

(3)

where

$$c_4\left(n_i\right)=\sqrt{\frac{2}{n_i-1}}\frac{\mathsf{\Gamma }(n_i/2)}{\mathsf{\Gamma }((n_i-1)/2)}.$$

(4)

This shows that $S_i/c_4\left(n_i\right)$ is actually an unbiased estimator of $\sigma$. Note that $c_4(·)$ is the normal-consistent unbiasing factor, which is a function of the sample size and this $c_4$ notation was originally used in ASQC [11] to the best of our knowledge. For more details on $c_4(·)$, the interested reader is referred to Vardeman [12]. In Appendix A, we also provide an approximate calculation of $c_4(·)$ which can be used for practically easier calculation.

Thus, we can estimate $\sigma$ using $S_i/c_4\left(n_i\right)$ which are are all unbiased. Analogous to Equation (1), one can use

$${\overline{S}}_A=\frac{S_1/c_4\left(n_1\right)+S_2/c_4\left(n_2\right)+\cdots +S_m/c_4\left(n_m\right)}{m}=\frac{1}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)},$$

(5)

which is clearly unbiased for $\sigma$. Since ${\sum }_{i=1}^{m}E\left[S_i\right]={\sum }_{i=1}^m{c}_4\left(n_i\right)·\sigma$ from Equation (3), the estimator below

$${\overline{S}}_B=\frac{S_1+S_2+\cdots +S_m}{c_4\left(n_1\right)+c_4\left(n_2\right)+\cdots +c_4\left(n_m\right)}=\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}$$

(6)

is also unbiased for $\sigma$. In addition, we consider the following unbiased estimator proposed by Burr [13]

$${\overline{S}}_C={\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}.$$

(7)

This estimator is actually the best linear unbiased estimator (BLUE) and we provide the proof below.

Theorem 1.

The estimator ${\overline{S}}_C$ in Equation (7) is the BLUE.

Proof. 

First, we consider a linear unbiased estimator in the form of ${\sum }_{i=1}^m{w}_i{S}_i$. Then, its variance and expectation are given by

$$\mathrm{Var}\left(\sum _{i=1}^m{w}_i{S}_i\right)=\sum _{i=1}^m{w}_i^2\{1-c_4{\left(n_i\right)}^2\}{\sigma }^2\mathrm{and}E\left(\sum _{i=1}^m{w}_i{S}_i\right)=\sum _{i=1}^m{w}_i{c}_4\left(n_i\right)\sigma .$$

To obtain the BLUE, we need to minimize $\mathrm{Var}\left({\sum }_{i=1}^m{w}_i{S}_i\right)$ with the unbiasedness condition $E\left({\sum }_{i=1}^m{w}_i{S}_i\right)=\sigma$. Thus, our objective is to minimize

$$\sum _{i=1}^m{w}_i^2\{1-c_4{\left(n_i\right)}^2\}\mathrm{subject}\mathrm{to}\sum _{i=1}^m{w}_i{c}_4\left(n_i\right)=1,$$

which can be easily solved by using the method of Lagrange multipliers. The auxiliary function with the Lagrange multiplier $\lambda$ is given by

$$\mathsf{\Psi }=\sum _{i=1}^m{w}_i^2\{1-c_4{\left(n_i\right)}^2\}-\lambda \left\{\sum _{i=1}^m{w}_i{c}_4\left(n_i\right)-1\right\}.$$

It is immediate from $\partial \mathsf{\Psi }/\partial w_k$ that $2w_k\{1-c_4{\left(n_k\right)}^2\}-\lambda c_4\left(n_k\right)=0$, which results in

$$w_k=\frac{\lambda c_4\left(n_k\right)}{2\{1-c_4{\left(n_k\right)}^2\}}.$$

(8)

Multiplying $c_4\left(n_k\right)$ to Equation (8) and then making the sum of the two sides, we have

$$\sum _{i=1}^m{w}_i{c}_4\left(n_i\right)=\frac{\lambda }{2}\sum _{i=1}^m\frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}.$$

Since ${\sum }_{i=1}^m{w}_i{c}_4\left(n_i\right)=1$, we first solve the above for $\lambda$ and then substitute it into Equation (8), which provides

$$w_k={\displaystyle \frac{{\displaystyle \frac{c_4\left(n_k\right)}{1-c_4{\left(n_k\right)}^2}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}},$$

which results in ${\overline{S}}_C={\sum }_{i=1}^m{w}_i{S}_i$. This completes the proof.  □

Remark 1.

It should be noted that, for the case of an equal sample size ($n_1=n_2=\cdots =n_m=n$), we can easily show that ${\overline{S}}_A={\overline{S}}_B={\overline{S}}_C=\overline{S}/c_4\left(n\right)$, where $\overline{S}={\sum }_{i=1}^m{S}_i/m$.

One can also incorporate the pooled sample variance in estimating $\sigma$ given by

$$S_p^2=\frac{{\displaystyle \sum _{i=1}^m(n_i-1)S_i^2}}{N-m},$$

(9)

where $N={\sum }_{i=1}^m{n}_i$ again. However, $S_p$ is not unbiased for $\sigma$ although $S_p^2$ is. This is because

$$E\left[S_p\right]=c_4(N-m+1)\sigma .$$

(10)

Based on this, Burr [13] suggested the following unbiased estimator of $\sigma$,

$${\overline{S}}_D=\frac{S_p}{c_4(N-m+1)}.$$

(11)

Remark 2.

The weighted average approach introduced in Subsection 6.3.2 of Montgomery [9] uses the pooled sample standard deviation $S_p$ from Equation (9) to estimate σ. However, since $c_4\left(x\right)<1$, which will be shown in Lemma 1, it is immediate from Equation (10) that $S_p$ clearly underestimates the true parameter σ.

We have introduced the four unbiased scale estimators of $\sigma$ which are denoted by ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, and ${\overline{S}}_D$. A natural question appears: which of the four estimators should be recommended for estimating $\sigma$ in practical applications? In the following section, we clarify this question by providing a guideline in terms of inequalities of their variances of the estimators under consideration.

3. Inequalities of the Variances of the Scale Estimators

We first obtain the variance of ${\overline{S}}_A$ in Equation (5), which is obtained as

$$\mathrm{Var}\left({\overline{S}}_A\right)=\frac{1}{m^2}\sum _{i=1}^m\frac{1}{c_4{\left(n_i\right)}^2}\mathrm{Var}\left(S_i\right).$$

(12)

Using Equation (3) and the unbiasedness property of $S_i^2$, we have

$$\mathrm{Var}\left(S_i\right)=E\left(S_i^2\right)-E{\left(S_i\right)}^2={\sigma }^2-c_4{\left(n_i\right)}^2{\sigma }^2={\sigma }^2\left\{1-c_4{\left(n_i\right)}^2\right\}.$$

(13)

Substituting Equation (13) into (12), we have

$$\mathrm{Var}\left({\overline{S}}_A\right)=\frac{{\sigma }^2}{m^2}\sum _{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}.$$

(14)

Similarly, the variance of ${\overline{S}}_B$ in Equation (6) is easily obtained as

$$\mathrm{Var}\left({\overline{S}}_B\right)=\frac{{\sum }_{i=1}^m\mathrm{Var}\left(S_i\right)}{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}={\sigma }^2·\frac{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}$$

(15)

and that of ${\overline{S}}_C$ in Equation (7) is also obtained as

$$\mathrm{Var}\left({\overline{S}}_C\right)=\frac{{\sigma }^2}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}.$$

(16)

Finally, for the case of ${\overline{S}}_D$, we have

$$\mathrm{Var}\left({\overline{S}}_D\right)=E\left({\overline{S}}_D^2\right)-E{\left({\overline{S}}_D\right)}^2.$$

(17)

Using Equation (11) and the unbiasedness property of ${\overline{S}}_D$, we can rewrite Equation (17) as

$$\mathrm{Var}\left({\overline{S}}_D\right)=\frac{1}{c_4{(N-m+1)}^2}E\left(S_p^2\right)-{\sigma }^2.$$

Since $S_p^2$ is also unbiased for ${\sigma }^2$, we have

$$\mathrm{Var}\left({\overline{S}}_D\right)={\sigma }^2\left\{\frac{1}{c_4{(N-m+1)}^2}-1\right\}.$$

(18)

Next, as aforementioned, we here answer the question of how to choose the best one among the four estimators by comparing their variances of these four scale estimators. To be more specific, we prove the inequality relations among the four scale estimators as follows.

Lemma 1.

The function $c_4\left(x\right)$ defined in Equation (4) is monotonically increasing and

$$\sqrt{\frac{2x-3}{2x-2}}<c_4\left(x\right)<1.$$

Proof. 

Using $(x-1)/2·\mathsf{\Gamma }\left(\right(x-1)/2)=\mathsf{\Gamma }\left(\right(x+1)/2)$, we can rewrite $c_4(·)$ in Equation (4) as

$$c_4\left(x\right)=\sqrt{\frac{2}{x-1}}\frac{\mathsf{\Gamma }(x/2)}{\mathsf{\Gamma }\left(\right(x-1)/2)}={\left[\frac{\mathsf{\Gamma }{(x/2)}^2}{\mathsf{\Gamma }\left(\right(x-1)/2)\mathsf{\Gamma }\left(\right(x+1)/2)}\right]}^{1/2},$$

(19)

where $x\ge 2$. Watson [14] showed that the function

$$\theta \left(x\right)=-x+x\frac{\mathsf{\Gamma }\left(x\right)\mathsf{\Gamma }(x+1)}{\mathsf{\Gamma }{(x+1/2)}^2}$$

(20)

is monotonically decreasing. It is noteworthy that the function $\theta \left(x\right)$ was motivated by Wallis’ famous infinite fraction for $\pi$ [15]. We can rewrite $c_4\left(x\right)$ using $\theta \left(x\right)$ in Equation (20) as

$$c_4\left(x\right)=\frac{1}{\sqrt{{\displaystyle \frac{\theta \left(\right(x-1)/2)}{(x-1)/2}+1}}}.$$

(21)

Since both $\theta \left(\right(x-1)/2)$ and $1/\left\{\right(x-1)/2\}$ are positive and monotonically decreasing,

$$\frac{1}{c_4\left(x\right)}=\sqrt{\theta ((x-1)/2)·\frac{1}{(x-1)/2}+1}$$

is also decreasing. Thus, $c_4\left(x\right)$ is monotonically increasing. Watson [14] and Mortici [16] also showed that

$$\sqrt{x+\frac{1}{4}}<\frac{\mathsf{\Gamma }(x+1)}{\mathsf{\Gamma }(x+\frac{1}{2})}\le \sqrt{x+\frac{1}{\pi }}.$$

Multiplying $\sqrt{2/(2x+1)}$ on each of the above terms, some algebra shows that

$$\sqrt{\frac{4x+1}{4x+2}}<c_4(2x+2)\le \sqrt{\frac{2x+2/\pi }{2x+1}}<1.$$

For convenience, we let $x^{*}=2x+2$. Then, we have

$$\sqrt{\frac{2x^{*}-3}{2x^{*}-2}}<c_4\left(x^{*}\right)<1,$$

which completes the proof.  □

Lemma 2

(Chebyshev’s sum inequality). If $a_1\ge a_2\ge \cdots \ge a_m$ and $b_1\ge b_2\ge \cdots \ge b_m$, then we have

$$m\sum _{i=1}^m{a}_i{b}_i\ge \sum _{i=1}^m{a}_i·\sum _{i=1}^m{b}_i.$$

Similarly, if $a_1\le a_2\le \cdots \le a_m$ and $b_1\ge b_2\ge \cdots \ge b_m$, then we have

$$m\sum _{i=1}^m{a}_i{b}_i\le \sum _{i=1}^m{a}_i·\sum _{i=1}^m{b}_i.$$

Proof. 

First, we will consider the case where both $a_i$ and $b_i$ are increasing. In this case, both $(a_i-a_j)$ and $(b_i-b_j)$ have the same sign, or at least one of them can have the zero value. Then, the value of $(a_i-a_j)(b_i-b_j)$ is positive or zero for any i and j, which results in

$$\sum _{i=1}^m\sum _{j=1}^m(a_i-a_j)(b_i-b_j)\ge 0.$$

After the tedious algebra of the above, we have

$$m\sum _{i=1}^m{a}_i{b}_i\ge \sum _{i=1}^m{a}_i·\sum _{i=1}^m{b}_i.$$

Next, we consider the case when $a_i$ is increasing and $b_i$ is decreasing. In this case, $(a_i-a_j)$ and $(b_i-b_j)$ have different signs, or at least one of them can have the zero value. Then, the value of $(a_i-a_j)(b_i-b_j)$ is negative or zero for any i and j. Similar to the above approach, we have ${\sum }_{i=1}^m{\sum }_{j=1}^m(a_i-a_j)(b_i-b_j)\le 0,$ which results in $m{\sum }_{i=1}^m{a}_i{b}_i\le {\sum }_{i=1}^m{a}_i·{\sum }_{i=1}^m{b}_i$. This completes the proof.  □

It should be noted that the above inequality name is coined after Pafnuty Lvovich Chebyshev (1821–1894) who mentioned it in a brief note [17]. He provided it in an integral form and his original proof can be found in Chebyshev [18]. For more details, the readers are also referred to Besenyei [19] and Section 2.17 of Hardy et al. [20].

Theorem 2.

We have

$$\mathrm{Var}\left({\overline{S}}_A\right)\ge \mathrm{Var}\left({\overline{S}}_B\right).$$

Proof. 

For convenience, we rearrange the sample sizes so that $n_1\le n_2\le \cdots \le n_m$. Then, it is immediate from Equations (14) and (15) that it suffices to show

$$\frac{1}{m^2}\sum _{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}\ge \frac{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}.$$

Since $c_4\left(x\right)$ is increasing from Lemma 1, it is easily seen that $1/c_4{\left(x\right)}^2-1$ is decreasing. For convenience, let $a_i=c_4\left(n_i\right)$ and $b_i=1/c_4{\left(n_i\right)}^2-1$. Then, we have $a_1\le a_2\le \cdots \le a_m$ and $b_1\ge b_2\ge \cdots \ge b_m$. Thus, we observe from Lemma 2 that

$$\sum _{i=1}^m{c}_4\left(n_i\right)·\sum _{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}\ge m·\sum _{i=1}^m\left\{\frac{1}{c_4\left(n_i\right)}-c_4\left(n_i\right)\right\}.$$

(22)

Applying Lemma 2 again with $a_i=c_4\left(n_i\right)$ (increasing) and $b_i=1/c_4\left(n_i\right)-c_4\left(n_i\right)$ (decreasing), we have

$$\sum _{i=1}^m{c}_4\left(n_i\right)·\sum _{i=1}^m\left\{\frac{1}{c_4\left(n_i\right)}-c_4\left(n_i\right)\right\}\ge m·\sum _{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\},$$

which results in

$$\sum _{i=1}^m\left\{\frac{1}{c_4\left(n_i\right)}-c_4\left(n_i\right)\right\}\ge \frac{m{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}.$$

(23)

Comparing Equations (22) and (23), we have

$$\sum _{i=1}^m{c}_4\left(n_i\right)·\sum _{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}\ge \frac{m^2{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\sum }_{i=1}^m{c}_4\left(n_i\right)},$$

which results in

$$\frac{1}{m^2}\sum _{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}\ge \frac{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}.$$

This completes the proof.  □

Theorem 3.

We have

$$\mathrm{Var}\left({\overline{S}}_B\right)\ge \mathrm{Var}\left({\overline{S}}_C\right).$$

Proof. 

We have the variances of ${\overline{S}}_B$ and ${\overline{S}}_C$ from Equations (15) and (16), which are given by

$$\mathrm{Var}\left({\overline{S}}_B\right)=\frac{{\sigma }^2·{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}}{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}\mathrm{and}\mathrm{Var}\left({\overline{S}}_C\right)=\frac{{\sigma }^2}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}.$$

Thus, it suffices to show

$${\left\{\sum _{i=1}^m{c}_4\left(n_i\right)\right\}}^2\le \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}·\sum _{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}.$$

(24)

For convenience, we let $a_i=c_4\left(n_i\right)/\sqrt{1-c_4{\left(n_i\right)}^2}$ and $b_i=\sqrt{1-c_4{\left(n_i\right)}^2}$. Then, it is immediate from the Cauchy–Schwarz inequality, ${\left\{\sum a_i{b}_i\right\}}^2\le \left\{\sum a_i^2\right\}\left\{\sum b_i^2\right\}$ that the inequality in Equation (24) holds. This completes the proof.  □

Lemma 3.

The function

$$\frac{c_4{\left(x\right)}^2}{1-c_4{\left(x\right)}^2}$$

is concave.

Proof. 

It is immediate from Equation (21) that we have

$$\frac{c_4{\left(x\right)}^2}{1-c_4{\left(x\right)}^2}=\frac{(x-1)/2}{\theta \left(\right(x-1)/2)}.$$

(25)

For convenience, we let $y=(x-1)/2$. Then, it suffices to show $y/\theta \left(y\right)$ is concave. Bustoz and Ismail [21] showed that $\theta \left(y\right)$ is completely monotonic on $[-1/2,\infty )$ using the representation by Watson [14]. For more details on complete monotonicity, refer to Section XIII.4 of Feller [22]. It is well known that completely monotonic functions are log-convex. For example, see Lemma 4.3 of Merkle [23], Theorem 1 of Fink [24], Exercise 6 in Section 2.1 of Niculescu and Persson [25], and Equation (3.4) of van Haeringen [26] with $n=0$ and $m=1$, among others.

Thus, $log\theta \left(y\right)$ is convex so $-log\theta \left(y\right)$ is concave. Since $logy$ is also concave, $log(y/\theta \left(y\right))=logy+\left\{-log\theta \left(y\right)\right\}$ is concave. This implies that $y/\theta \left(y\right)$ is log-concave. The log-concavity of $y/\theta \left(y\right)$ guarantees that it is concave. This completes the proof.  □

Theorem 4.

We have

$$\mathrm{Var}\left({\overline{S}}_C\right)\ge \mathrm{Var}\left({\overline{S}}_D\right).$$

Proof. 

From Equations (16) and (18), we have

$$\mathrm{Var}\left({\overline{S}}_C\right)=\frac{{\sigma }^2}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}\mathrm{and}\mathrm{Var}\left({\overline{S}}_D\right)={\sigma }^2\left\{\frac{1}{c_4{(N-m+1)}^2}-1\right\}.$$

Thus, it suffices to show

$$\frac{c_4{(N-m+1)}^2}{1-c_4{(N-m+1)}^2}\ge \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}.$$

It is immediate from Lemma 3 that $c_4{\left(x\right)}^2/\{1-c_4{\left(x\right)}^2\}$ is concave. Then, using Jensen’s inequality, we have

$$\frac{c_4{\left(\overline{n}\right)}^2}{1-c_4{\left(\overline{n}\right)}^2}\ge \frac{1}{m}\sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}},$$

where $\overline{n}={\sum }_{i=1}^m{n}_i/m$. Thus, it suffices to show

$$\frac{c_4{(N-m+1)}^2}{1-c_4{(N-m+1)}^2}\ge \frac{m{c}_4{\left(\overline{n}\right)}^2}{1-c_4{\left(\overline{n}\right)}^2}.$$

(26)

Using Equation (25) in Lemma 3, we have

$$\begin{array}{c}\hfill \frac{c_4{(N-m+1)}^2}{1-c_4{(N-m+1)}^2}=\frac{(N-m)/2}{\theta \left(\right(N-m)/2)}\hfill \end{array}$$

(27)

and

$$\begin{array}{c}\hfill \frac{m{c}_4{\left(\overline{n}\right)}^2}{1-c_4{\left(\overline{n}\right)}^2}=\frac{m(\overline{n}-1)/2}{\theta ((\overline{n}-1)/2)}=\frac{(N-m)/2}{\theta ((\overline{n}-1)/2)},\hfill \end{array}$$

(28)

where $\theta \left(x\right)$ is defined in Equation (20). Comparing Equations (27) and (28), we need to show

$$\theta ((N-m)/2)\le \theta ((\overline{n}-1)/2).$$

Since $(\overline{n}-1)/2\le (N-m)/2$ and $\theta \left(x\right)$ is decreasing from Watson [14], the above inequality in Equation (26) holds. This completes the proof.  □

In combination with the inequalities in Theorems 2–4, we have the following result:

$$\mathrm{Var}\left({\overline{S}}_A\right)\ge \mathrm{Var}\left({\overline{S}}_B\right)\ge \mathrm{Var}\left({\overline{S}}_C\right)\ge \mathrm{Var}\left({\overline{S}}_D\right).$$

(29)

Lemma 4.

The $c_4\left(x\right)$ defined in Equation (4) satisfies

$$\frac{1}{m}\sum _{i=1}^m{c}_4\left(n_i\right)\le c_4\left(\overline{n}\right).$$

Proof. 

First, we will show that $c_4\left(x\right)$ is concave. Taking the logarithm of $c_4\left(x\right)$ in Equation (19), we have

$$log{c}_4\left(x\right)=log\mathsf{\Gamma }\left(\frac{x}{2}\right)-\frac{1}{2}log\mathsf{\Gamma }\left(\frac{x-1}{2}\right)-\frac{1}{2}log\mathsf{\Gamma }\left(\frac{x+1}{2}\right),$$

where $x\ge 2$. It is well known that the second derivative of $log\mathsf{\Gamma }\left(x\right)$ can be expressed as the sum of the series

$$\frac{d^2}{d{x}^2}log\mathsf{\Gamma }\left(x\right)=\sum _{k=0}^{\infty }\frac{1}{{(x+k)}^2}.$$

(30)

For more details, see Merkle [27] and Section 11.14 (iv) of Schilling [28].

Using Equation (30), we have

$$\begin{array}{cc}\hfill \frac{d^2}{d{x}^2}log{c}_4\left(x\right)\hfill & =\sum _{k=0}^{\infty }\frac{1}{{(x/2+k)}^2}-\frac{1}{2}\sum _{k=0}^{\infty }\frac{1}{{(x/2-1/2+k)}^2}-\frac{1}{2}\sum _{k=0}^{\infty }\frac{1}{{(x/2+1/2+k)}^2}\hfill \\ \hfill \hfill & =\sum _{k=0}^{\infty }\left[\frac{1}{{(x/2+k)}^2}-\frac{1}{2{(x/2-1/2+k)}^2}-\frac{1}{2{(x/2+1/2+k)}^2}\right]\hfill \\ \hfill \hfill & =-\sum _{k=0}^{\infty }\frac{\frac{3}{4}{(x/2-1/2+k)}^2+\frac{3}{4}(x/2-1/2+k)+\frac{1}{8}}{{(x/2+k)}^2{(x/2-1/2+k)}^2{(x/2+1/2+k)}^2}<0.\hfill \end{array}$$

Thus, $c_4\left(x\right)$ is log-concave which guarantees that $c_4\left(x\right)$ is concave. Then, using Jensen’s inequality, we have

$$\frac{1}{m}\sum _{i=1}^m{c}_4\left(n_i\right)\le c_4\left(\overline{n}\right),$$

where $\overline{n}={\sum }_{i=1}^m{n}_i/m$. This completes the proof.  □

Remark 3.

It is worth mentioning that one can argue an alternative way based on an average sample size $\overline{n}={\sum }_{i=1}^m{n}_i/m$. For example, see Section 6.3.2 of Montgomery [9]. In such setting, one may consider the following quotient for the estimator of σ

$${\overline{S}}^{*}=\frac{\overline{S}}{c_4\left(\overline{n}\right)},$$

(31)

where

$$\overline{S}=\frac{1}{m}\sum _{i=1}^m{S}_i.$$

(32)

However, this estimator is not unbiased because

$$E\left[{\overline{S}}^{*}\right]=\frac{\frac{1}{m}{\sum }_{i=1}^{m}E\left[S_i\right]}{c_4\left(\overline{n}\right)\sigma }=\frac{\frac{1}{m}{\sum }_{i=1}^m{c}_4\left(n_i\right)}{c_4\left(\overline{n}\right)}·\sigma .$$

Because $\frac{1}{m}{\sum }_{i=1}^m{c}_4\left(n_i\right)\le c_4\left(\overline{n}\right)$ from Lemma 4, ${\overline{S}}^{*}$ can underestimate the true parameter σ. However, for the case of an equal sample size $(n_1=n_2=\cdots =n_m=n)$, we have ${\overline{S}}^{*}=\overline{S}/c_4\left(n\right)$ so that ${\overline{S}}_A={\overline{S}}_B={\overline{S}}_C={\overline{S}}^{*}$. Thus, in this special case of an equal sample size, ${\overline{S}}^{*}$ is unbiased.

Remark 4.

There is another alternative in Subsection 3.8.B of ASTM [10], which is based on the weighted average of sample standard deviations given by

$${\overline{S}}_w=\frac{n_1{S}_1+n_2{S}_2+\cdots +n_m{S}_m}{n_1+n_2+\cdots +n_m}=\frac{{\sum }_{i=1}^m{n}_i{S}_i}{N}.$$

(33)

Then, we have

$$E\left[{\overline{S}}_w\right]=\frac{{\sum }_{i=1}^m{n}_i{c}_4\left(n_i\right)}{N}·\sigma .$$

It deserves mentioning that the estimator above still underestimates the true parameter σ because $c_4\left(x\right)<1$ from Lemma 1. In addition, for the case of an equal sample size, we have ${\overline{S}}_w=\overline{S}$, where $\overline{S}=\frac{1}{m}{\sum }_{i=1}^m{S}_i$.

Theorem 5.

Let $S_N^2={\sum }_{i=1}^m{\sum }_{j=1}^{n_i}{(X_{ij}-\overline{\overline{X}})}^2/(N-1)$ and $\overline{\overline{X}}=\frac{1}{N}{\sum }_{i=1}^m{\sum }_{j=1}^{n_i}{X}_{ij}$. Then, $\overline{\overline{X}}$ and ${\overline{S}}_E=S_N/c_4\left(N\right)$ are the uniform minimum variance unbiased estimators of μ and σ, respectively. Thus, we have

$$\mathrm{Var}\left({\overline{S}}_D\right)\ge \mathrm{Var}\left({\overline{S}}_E\right).$$

Proof. 

One can obtain the uniform minimum variance unbiased estimator (UMVUE) using complete sufficient statistics as described in Theorem 7.3.23 of Casella and Berger [29]. For more details on the UMVUE, one can refer to Definition 1.6 of Lehmann and Casella [30].

The control charts we have developed are under the assumption that $X_{ij}$ are independent and identically distributed as normal with mean $\mu$ and variance ${\sigma }^2$ which leads to the joint complete sufficient statistics for $\mu$ and ${\sigma }^2$ given by $\overline{\overline{X}}=\frac{1}{N}{\sum }_{i=1}^m{\sum }_{j=1}^{n_i}{X}_{ij}$ and ${\sum }_{i=1}^m{\sum }_{j=1}^{n_i}{(X_{ij}-\overline{\overline{X}})}^2$, respectively. Since $E\left[\overline{\overline{X}}\right]=\mu$ and $E\left[{\overline{S}}_E\right]=\sigma$, $\overline{\overline{X}}$ and ${\overline{S}}_E$ are the UMVUEs of $\mu$ and $\sigma$, respectively. This completes the proof.  □

Remark 5.

It is noteworthy that, analogous to Equation (18), we have

$$\mathrm{Var}\left({\overline{S}}_E\right)={\sigma }^2\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\},$$

(34)

which also results in $\mathrm{Var}\left({\overline{S}}_D\right)\ge \mathrm{Var}\left({\overline{S}}_E\right)$ since $c_4\left(x\right)$ is increasing from Lemma 1.

Remark 6.

Although ${\overline{S}}_E$ can attain the minimum variance, we do not adopt this estimator to construct the control charts. The main reason can be discussed as follows. Consider the case that $X_{ij}$ are independent and identically distributed as normal with different means ${\mu }_i$ for $j=1,2,\dots ,n_i$ and variance ${\sigma }^2$. Then, the pooled sample variance, $S_p^2$, is a complete sufficient statistic. See Example 2.3 in Section 2.2 of Lehmann and Casella [30]. This implies that ${\overline{S}}_E$ is better under the null hypothesis $H_0:{\mu }_1={\mu }_2=\cdots ={\mu }_m$ (in control), whereas ${\overline{S}}_D$ is better under the alternative (out of control). Thus, if the process is out of control, the control charts using ${\overline{S}}_E$ could have wider control limits, which may result in an increase in the rate incorrectly signaling that the process is in-control when the process is actually out of control.

In addition, one can think of the case of non-constancy of σ. In this heteroscedasticity case, Burr [13] mentioned that ${\overline{S}}_C$ seems preferable to ${\overline{S}}_D$. We think that this case should be investigated more thoroughly in a sequel paper.

4. Comparison of the Performance

In this section, we provide the relative efficiency of the methods to compare their performance. We also carried out Monte Carlo simulations to compare the empirical biases and variances.

4.1. Relative Efficiency

When we compare the performance of unbiased estimators (say, ${\widehat{\theta }}_1$ and ${\widehat{\theta }}_0$), the relative efficiency (RE) is widely used in the statistics literature. See Section 2.2 of Lehmann [31] for more details. The RE of ${\widehat{\theta }}_1$ with respect to ${\widehat{\theta }}_0$ is given by

$$\mathrm{RE}({\widehat{\theta }}_1\mid {\widehat{\theta }}_0)=\frac{\mathrm{Var}\left({\widehat{\theta }}_0\right)}{\mathrm{Var}\left({\widehat{\theta }}_1\right)},$$

where ${\widehat{\theta }}_0$ is a reference estimator. In general, the estimator with the smaller variance of the two estimators is used as a reference estimator so that $\mathrm{RE}({\widehat{\theta }}_1\mid {\widehat{\theta }}_0)\le 1$.

To estimate the location parameter, we considered ${\overline{\overline{X}}}_A$ in Equation (1) and ${\overline{\overline{X}}}_B$ in Equation (2). Then, the RE of ${\overline{\overline{X}}}_A$ with respect to ${\overline{\overline{X}}}_B$ is easily obtained as

$$\mathrm{RE}({\overline{\overline{X}}}_A\mid {\overline{\overline{X}}}_B)=\frac{m^2}{\left({\sum }_{i=1}^m{n}_i\right)·\left({\sum }_{i=1}^m{n}_i^{-1}\right)}.$$

Note that $\mathrm{RE}({\overline{\overline{X}}}_A\mid {\overline{\overline{X}}}_B)\le 1$ where the equality holds if and only if $n_1=n_2=\cdots =n_m$ due to the inequality of the arithmetic mean and the harmonic mean [32].

For the case of the scale parameter, we considered the five estimators. Among them, ${\overline{S}}_E$ has the smallest variance. Thus, it can be used as a reference to compare the performance of the scale estimators and the RE is then given by

$$\mathrm{RE}({\overline{S}}_j\mid {\overline{S}}_E)=\frac{\mathrm{Var}\left({\overline{S}}_E\right)}{\mathrm{Var}\left({\overline{S}}_j\right)},$$

where $j=A,B,C,D$. For notational brevity, we denote

$$\mathrm{RE}\left({\overline{S}}_j\right)=\mathrm{RE}({\overline{S}}_j\mid {\overline{S}}_E).$$

(35)

It is immediate from Equations (14)–(16), (18), and (34) that we have

$$\begin{array}{cc}\hfill \mathrm{RE}\left({\overline{S}}_A\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{m^2}{{\sum }_{i=1}^m\left\{\frac{1}{c_4{\left(n_i\right)}^2}-1\right\}},\hfill \\ \hfill \mathrm{RE}\left({\overline{S}}_B\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)\right\}}^2}{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}},\hfill \\ \hfill \mathrm{RE}\left({\overline{S}}_C\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\sum _{i=1}^m\frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2},\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \mathrm{RE}\left({\overline{S}}_D\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{c_4{(N-m+1)}^2}{1-c_4{(N-m+1)}^2}.\hfill \end{array}$$

It can be easily seen from Equations (29) and (34) that we have $\mathrm{RE}\left({\overline{S}}_A\right)\le \mathrm{RE}\left({\overline{S}}_B\right)\le \mathrm{RE}\left({\overline{S}}_C\right)\le \mathrm{RE}\left({\overline{S}}_D\right)\le 1$. In particular, when $n_1=n_2=\cdots =n_m$, we have $\mathrm{RE}\left({\overline{S}}_A\right)=\mathrm{RE}\left({\overline{S}}_B\right)=\mathrm{RE}\left({\overline{S}}_C\right)\le \mathrm{RE}\left({\overline{S}}_D\right)$.

We have considered the RE to compare the performance of the above unbiased estimators. However, ${\overline{S}}^{*}$ in Equation (31) and ${\overline{S}}_w$ in Equation (33) are not unbiased as mentioned in Remarks 2 and 3, respectively. In addition, $\overline{S}$ in Equation (32) is not unbiased as easily seen from Equation (3). Thus, it is reasonable to consider the mean square error (MSE) to obtain the RE since the MSE can be regarded as a overall measure of bias and dispersion. For other utilization and modification of the RE, one can refer to Park et al. [33,34] and Ouyang et al. [35] which considered the ratio of the determinants of the covariance matrices, that is, the generalized variances [36,37]. Since the MSE is the same as the variance for unbiased estimators (${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, and ${\overline{S}}_D$), the RE based on the MSE is the same as that based on the variance in this unbiased case. In addition, the variance of ${\overline{S}}_E$ is the same as its MSE. Thus, we consider the RE of $\overline{S}$, ${\overline{S}}^{*}$, and ${\overline{S}}_w$ as follows. We denote them by

$$\mathrm{RE}\left(\overline{S}\right)=\frac{\mathrm{Var}\left({\overline{S}}_E\right)}{\mathrm{MSE}\left(\overline{S}\right)},\mathrm{RE}\left({\overline{S}}^{*}\right)=\frac{\mathrm{Var}\left({\overline{S}}_E\right)}{\mathrm{MSE}\left({\overline{S}}^{*}\right)},\mathrm{and}\mathrm{RE}\left({\overline{S}}_w\right)=\frac{\mathrm{Var}\left({\overline{S}}_E\right)}{\mathrm{MSE}\left({\overline{S}}_w\right)}.$$

We next obtain their biases which are easily obtained using Equation (3)

$$\begin{array}{cc}\mathrm{Bias}\left(\overline{S}\right)\hfill & =\left[\frac{1}{m}\sum _{i=1}^m{c}_4\left(n_i\right)-1\right]\sigma ,\hfill \\ \mathrm{Bias}\left({\overline{S}}^{*}\right)\hfill & =\left[\frac{{\sum }_{i=1}^m{c}_4\left(n_i\right)}{m{c}_4\left(\overline{n}\right)}-1\right]\sigma ,\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \hfill & \mathrm{Bias}\left({\overline{S}}_w\right)\hfill & =\left[\frac{{\sum }_{i=1}^m{n}_i{c}_4\left(n_i\right)}{N}-1\right]\sigma .\hfill \end{array}$$

In addition, the variances are also easily obtained by using $\mathrm{Var}\left(S_i\right)={\sigma }^2\left\{1-c_4{\left(n_i\right)}^2\right\}$ so that we have

$$\begin{array}{cc}\mathrm{Var}\left(\overline{S}\right)\hfill & =\frac{{\sigma }^2}{m^2}\sum _{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\},\hfill \\ \mathrm{Var}\left({\overline{S}}^{*}\right)\hfill & =\frac{1}{c_4{\left(\overline{n}\right)}^2}\frac{{\sigma }^2}{m^2}\sum _{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\},\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \hfill & \mathrm{Var}\left({\overline{S}}_w\right)\hfill & =\frac{{\sigma }^2}{N^2}\sum _{i=1}^m{n}_i^2\left\{1-c_4{\left(n_i\right)}^2\right\}.\hfill \end{array}$$

Considering that the MSE is the variance plus the squared bias, we can obtain the RE based on the MSE

$$\begin{array}{cc}\mathrm{RE}\left(\overline{S}\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{m^2}{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}+{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)-m\right\}}^2},\hfill \\ \mathrm{RE}\left({\overline{S}}^{*}\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{m^2{c}_4{\left(\overline{n}\right)}^2}{{\sum }_{i=1}^m\left\{1-c_4{\left(n_i\right)}^2\right\}+{\left\{{\sum }_{i=1}^m{c}_4\left(n_i\right)-m{c}_4\left(\overline{n}\right)\right\}}^2},\hfill \end{array}$$

and

$$\begin{array}{cc}\mathrm{RE}\left({\overline{S}}_w\right)\hfill & =\left\{\frac{1}{c_4{\left(N\right)}^2}-1\right\}·\frac{N^2}{{\sum }_{i=1}^m{n}_i^2\left\{1-c_4{\left(n_i\right)}^2\right\}+{\left\{{\sum }_{i=1}^m{n}_i{c}_4\left(n_i\right)-N\right\}}^2}.\hfill \end{array}$$

4.2. Empirical Biases and Variances

The RE is a useful statistical tool to compare the performance of unbiased estimators. However, the conventional methods such as in Equations (31)–(33) are biased. To compare these with the methods provided here, we obtain the empirical biases and variances by carrying out Monte Carlo simulations.

We consider two cases: an equal sample size and unequal sample sizes. We first generated samples of equal size with $n_1=n_2=n_3=3$, with $n_1=n_2=n_3=10$, and  with $n_1=n_2=n_3=20$. We next generated samples of unequal sizes with $n_1=3$, $n_2=5$, $n_3=7$, with $n_1=5$, $n_2=10$, $n_3=15$, and with $n_1=10$, $n_2=20$, $n_3=30$. Again, let $X_{ij}$ be the ith sample (subgroup) of size $n_i$. Then, $X_{ij}$ were generated from the normal distribution with mean ${\mu }_0=100$ and standard deviation ${\sigma }_0=10$ and we obtained the scale estimates including the unbiased estimators (${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, ${\overline{S}}_E$) and the conventional methods ($\overline{S}$, ${\overline{S}}^{*}$, ${\overline{S}}_w$). In order to obtain the empirical biases and variances of these estimates, we repeated this simulation ten million times ($I={10}^7$) and the results are summarized in

Table 1. The empirical bias, variance, MSE, and RE with an equal sample size.

	${\overline{\mathit{S}}}_{\mathit{A}}$	${\overline{\mathit{S}}}_{\mathit{B}}$	${\overline{\mathit{S}}}_{\mathit{C}}$	${\overline{\mathit{S}}}_{\mathit{D}}$	${\overline{\mathit{S}}}_{\mathit{E}}$	$\overline{\mathit{S}}$	${\overline{\mathit{S}}}^{*}$	${\overline{\mathit{S}}}_{\mathit{w}}$
	($n_1=n_2=n_3=3$)
Bias	0.0005	0.0005	0.0005	0.0004	0.0008	$-1.1372$	0.0005	$-1.1372$
Var	9.1087	9.1087	9.1087	8.6503	6.4354	7.1539	9.1087	7.1539
MSE	9.1087	9.1087	9.1087	8.6503	6.4354	8.4473	9.1087	8.4473
RE	0.7065	0.7065	0.7065	0.7439	1.0000	0.7618	0.7065	0.7618
	($n_1=n_2=n_3=10$)
Bias	0.0001	0.0001	0.0001	0.0002	0.0002	$-0.2733$	0.0001	$-0.2733$
Var	1.9007	1.9007	1.9007	1.8689	1.7388	1.7982	1.9007	1.7982
MSE	1.9007	1.9007	1.9007	1.8689	1.7388	1.8729	1.9007	1.8729
RE	0.9148	0.9148	0.9148	0.9304	1.0000	0.9284	0.9148	0.9284
	($n_1=n_2=n_3=20$)
Bias	0.0001	0.0001	0.0001	0.0001	0.0002	$-0.1305$	0.0001	$-0.1305$
Var	0.8885	0.8885	0.8885	0.8811	0.8511	0.8655	0.8885	0.8655
MSE	0.8885	0.8885	0.8885	0.8811	0.8511	0.8825	0.8885	0.8825
RE	0.9578	0.9578	0.9578	0.9660	1.0000	0.9644	0.9578	0.9644

for the case of an equal sample size and

Table 2. The empirical bias, variance, MSE, and RE with unequal sample sizes.

	${\overline{\mathit{S}}}_{\mathit{A}}$	${\overline{\mathit{S}}}_{\mathit{B}}$	${\overline{\mathit{S}}}_{\mathit{C}}$	${\overline{\mathit{S}}}_{\mathit{D}}$	${\overline{\mathit{S}}}_{\mathit{E}}$	$\overline{\mathit{S}}$	${\overline{\mathit{S}}}^{*}$	${\overline{\mathit{S}}}_{\mathit{w}}$
	($n_1=3$, $n_2=5$, $n_3=7$)
Bias	0.0008	0.0008	0.0009	0.0009	0.0007	$-0.7140$	$-0.1211$	$-0.6164$
Var	5.4627	5.2936	4.3850	4.2500	3.6337	4.5640	5.1653	3.8867
MSE	5.4627	5.2936	4.3850	4.2500	3.6337	5.0738	5.1800	4.2667
RE	0.6652	0.6864	0.8287	0.8550	1.0000	0.7162	0.7015	0.8517
	($n_1=5$, $n_2=10$, $n_3=15$)
Bias	0.0006	0.0006	0.0003	0.0003	0.0002	$-0.3496$	$-0.0783$	$-0.2792$
Var	2.5013	2.4512	1.8990	1.8685	1.7388	2.2825	2.4126	1.7990
MSE	2.5013	2.4512	1.8990	1.8685	1.7388	2.4047	2.4188	1.8770
RE	0.6952	0.7094	0.9157	0.9306	1.0000	0.7231	0.7189	0.9264
	($n_1=10$, $n_2=20$, $n_3=30$)
Bias	$-0.0001$	$-0.0001$	0.0001	0.0001	0.0002	$-0.1634$	$-0.0331$	$-0.1319$
Var	1.1229	1.1138	0.8885	0.8811	0.8511	1.0777	1.1064	0.8657
MSE	1.1229	1.1138	0.8885	0.8811	0.8511	1.1044	1.1075	0.8831
RE	0.7579	0.7641	0.9579	0.9659	1.0000	0.7706	0.7684	0.9637

for the case of unequal sample sizes. In addition, the MSEs along with the squared empirical biases (red) and variances (light blue) are plotted in Figure 1 for the case of an equal sample size and in Figure 2 for the case of unequal sample sizes. In the figures, ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, $\overline{S}$, ${\overline{S}}^{*}$, and  ${\overline{S}}_w$ are denoted by A, B, C, D, S, S*, and Sw, respectively.

Comparing the simulation results in Table 1, we can observe that the empirical results of ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, and ${\overline{S}}^{*}$ are the same for the case of an equal sample size. These results are quite reasonable as pointed out in Remark 3. We can notice that the empirical biases are not exactly zero, although quite negligible, because these biases are due to a random phenomenon of Monte Carlo simulation. In addition, $\overline{S}$ and ${\overline{S}}_w$ have the same results as pointed out in Remark 4. On the other hand, for the case of unequal sample sizes, they all have different results.

For the case of an equal sample size, the empirical variances and biases are noticeably different for a smaller sample size, but all of them are getting closer as the sample size is increasing as also shown in Figure 1. For the case of unequal sample sizes, the empirical biases are getting smaller as the sample sizes are increasing. However, the variances are still noticeably different even with large sample sizes.

Comparing the empirical variances only, ${\overline{S}}_w$ performs very well, whereas it is seriously biased. Thus, it is reasonable to compare the MSEs along with the biases and we can conclude that ${\overline{S}}_D$ is overall the best. Note that we do not recommend the use of ${\overline{S}}_E$, even though ${\overline{S}}_E$ always has the best results in all the measures. This is because it can lead to degraded performance when the process is out of control, as aforementioned in Remark 6.

5. Construction of the Control Charts with Unequal Sample Sizes

We briefly introduce how to construct the control charts and then discuss how to implement the estimators provided here to construct the S chart in Section 5.1 and improve the S chart using the probability limits in Section 5.2. We also discuss the construction of the $\overline{X}$ chart in Section 5.3.

In general, we construct statistical quality control charts based on two phases [9,38], usually denoted by Phase-I and Phase-II. Then, one can establish control limits with a set of stable manufacturing process data in Phase-I. Then, we monitor the process in Phase-II using the control limits obtained in Phase-I. We assume that we have m samples from a stable manufacturing process (Phase-I) and each sample has different sample sizes, denoted by $n_i$ where $i=1,2,\dots ,m$. Then, we monitor the process with a sample of size $n_k$ (Phase-II).

5.1. The S Chart

From the statistical asymptotic theory (for example, Corollarys 6–10 of Arnold [39]), we can have an approximate distribution

$$\frac{S_k-E\left(S_k\right)}{\mathrm{SE}\left(S_k\right)}\stackrel{•}{\sim }N(0,1),$$

where $S_k$ is the sample standard deviation with sample size $n_k$. In order to construct the $\mathrm{CL}\pm 3·\mathrm{SE}$ control limits, we can set up $(S_k-E\left(S_k\right))/\mathrm{SE}\left(S_k\right)=\pm 3$. Solving this for $S_k$, we have

$$E\left(S_k\right)\pm 3·\mathrm{SE}\left(S_k\right)=c_4\left(n_k\right)\sigma \pm 3\sqrt{1-c_4{\left(n_k\right)}^2}\sigma .$$

Since $\sigma$ is generally unknown, we need to estimate $\sigma$. One can use ${\overline{S}}_A$ in Equation (5), ${\overline{S}}_B$ in Equation (6), ${\overline{S}}_C$ in Equation (7), and ${\overline{S}}_D$ in Equation (11). Using ${\overline{S}}_A$, we can construct the S chart with sample size $n_k$ as follows:

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_A\hfill & =\frac{c_4\left(n_k\right)}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)}+3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{1}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)}\hfill \\ \hfill {\mathrm{CL}}_A\hfill & =\frac{c_4\left(n_k\right)}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)}\hfill \\ \hfill {\mathrm{LCL}}_A\hfill & =\frac{c_4\left(n_k\right)}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)}-3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{1}{m}\sum _{i=1}^m\frac{S_i}{c_4\left(n_i\right)},\hfill \end{array}$$

where we assign zero to $\mathrm{LCL}$ if it is negative. Next, using ${\overline{S}}_B$, we construct the S chart as follows:

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_B\hfill & =c_4\left(n_k\right)·\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}+3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}\hfill \\ \hfill {\mathrm{CL}}_B\hfill & =c_4\left(n_k\right)·\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}\hfill \\ \hfill {\mathrm{LCL}}_B\hfill & =c_4\left(n_k\right)·\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}-3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{{\sum }_{i=1}^m{S}_i}{{\sum }_{i=1}^m{c}_4\left(n_i\right)}.\hfill \end{array}$$

In addition, using ${\overline{S}}_C$, we can construct the S chart as follows:

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_C\hfill & =c_4\left(n_k\right)·{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}+3\sqrt{1-c_4{\left(n_k\right)}^2}·{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}\hfill \\ \hfill {\mathrm{CL}}_C\hfill & =c_4\left(n_k\right)·{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}\hfill \\ \hfill {\mathrm{LCL}}_C\hfill & =c_4\left(n_k\right)·{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}-3\sqrt{1-c_4{\left(n_k\right)}^2}·{\displaystyle \frac{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4\left(n_i\right)S_i}{1-c_4{\left(n_i\right)}^2}}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \frac{c_4{\left(n_i\right)}^2}{1-c_4{\left(n_i\right)}^2}}}}}.\hfill \end{array}$$

Of particular note is that, for the case of an equal sample size ($n_1=n_2=\cdots =n_m=n$) in Phase-I, it is easily seen that ${\mathrm{UCL}}_A={\mathrm{UCL}}_B={\mathrm{UCL}}_C$, ${\mathrm{CL}}_A={\mathrm{CL}}_B={\mathrm{CL}}_C$, and  ${\mathrm{LCL}}_A={\mathrm{LCL}}_B={\mathrm{LCL}}_C$. Thus, we have the S chart below with sample size $n_k$ to use in Phase-II

$$\begin{array}{cc}\hfill \mathrm{UCL}\hfill & =\frac{c_4\left(n_k\right)}{c_4\left(n\right)}\overline{S}+\frac{3\sqrt{1-c_4{\left(n_k\right)}^2}}{c_4\left(n\right)}\overline{S}\hfill \\ \hfill \mathrm{CL}\hfill & =\frac{c_4\left(n_k\right)}{c_4\left(n\right)}\overline{S}\hfill \\ \hfill \mathrm{LCL}\hfill & =\frac{c_4\left(n_k\right)}{c_4\left(n\right)}\overline{S}-\frac{3\sqrt{1-c_4{\left(n_k\right)}^2}}{c_4\left(n\right)}\overline{S},\hfill \end{array}$$

(36)

where $\overline{S}={\sum }_{i=1}^m{S}_i/m$. We assign zero to $\mathrm{LCL}$ if it is negative. Furthermore, if we assume $n_k=n$, we have

$$\begin{array}{cc}\hfill \mathrm{UCL}\hfill & =\overline{S}+\frac{3\sqrt{1-c_4{\left(n\right)}^2}}{c_4\left(n\right)}\overline{S}=B_4\left(n\right)\overline{S}\hfill \\ \hfill \mathrm{CL}\hfill & =\overline{S}\hfill \\ \hfill \mathrm{LCL}\hfill & =\overline{S}-\frac{3\sqrt{1-c_4{\left(n\right)}^2}}{c_4\left(n\right)}\overline{S}=B_3\left(n\right)\overline{S},\hfill \end{array}$$

where $B_3\left(n\right)=max\left\{1-3\sqrt{1-c_4{\left(n\right)}^2}/c_4\left(n\right),0\right\}$, and $B_4\left(n\right)=1+3\sqrt{1-c_4{\left(n\right)}^2}/c_4\left(n\right)$. This is a well-known S chart introduced in the quality control literature. For example, see the chart in Equation (6.27) of Montgomery [9]. This indicates that the proposed S chart includes the existing S chart as a special case.

Using ${\overline{S}}_D$, we construct the S chart as follows:

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_D\hfill & =c_4\left(n_k\right)·\frac{S_p}{c_4(N-m+1)}+3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{S_p}{c_4(N-m+1)}\hfill \\ \hfill {\mathrm{CL}}_D\hfill & =c_4\left(n_k\right)·\frac{S_p}{c_4(N-m+1)}\hfill \\ \hfill {\mathrm{LCL}}_D\hfill & =c_4\left(n_k\right)·\frac{S_p}{c_4(N-m+1)}-3\sqrt{1-c_4{\left(n_k\right)}^2}·\frac{S_p}{c_4(N-m+1)},\hfill \end{array}$$

where $S_p={\left\{{\sum }_{i=1}^m(n_i-1)S_i^2/(N-m)\right\}}^{1/2}$ and $N={\sum }_{i=1}^m{n}_i$. Unlike the previous cases, this control chart is not the same as the one in Equation (36) even for the case of an equal sample size.

5.2. The S and $S^2$ Charts with Probability Limits

We can improve the above S charts by using probability limits as mentioned in Section 4.7.4 of Ryan [40]. It follows from the following Chi-square distribution result

$$\frac{(n_k-1)S_k^2}{{\sigma }^2}\sim {\chi }_{n_k-1}^2$$

that

$$P\left[{\chi }_{1-\alpha /2,n_k-1}^2\le \frac{(n_k-1)S_k^2}{{\sigma }^2}\le {\chi }_{\alpha /2,n_k-1}^2\right]=1-\alpha ,$$

(37)

where ${\chi }_{\gamma ,\nu }^2$ is the $\gamma$th upper quantile of the Chi-square distribution with $\nu$ degrees of freedom. Rewriting Equation (37) about $S_k$, we then have

$$P\left[\sigma ·\sqrt{\frac{{\chi }_{1-\alpha /2,n_k-1}^2}{n_k-1}}\le S_k\le \sigma ·\sqrt{\frac{{\chi }_{\alpha /2,n_k-1}^2}{n_k-1}}\right]=1-\alpha .$$

Thus, we can construct the S chart with probability limits such that $\mathrm{UCL}=\sigma {\{{\chi }_{\alpha /2,n_k-1}^2/(n_k-1)\}}^{1/2}$, $\mathrm{CL}=\sigma$, and $\mathrm{LCL}=\sigma {\{{\chi }_{1-\alpha /2,n_k-1}^2/(n_k-1)\}}^{1/2}$. In practice, since $\sigma$ is unknown, we need to estimate $\sigma$. Thus, with the estimator $\widehat{\sigma }$, we can obtain

$$\begin{array}{cc}\hfill \mathrm{UCL}\hfill & =\widehat{\sigma }·\sqrt{\frac{{\chi }_{\alpha /2,n_k-1}^2}{n_k-1}}\hfill \\ \hfill \mathrm{CL}\hfill & =\widehat{\sigma }\hfill \\ \hfill \mathrm{LCL}\hfill & =\widehat{\sigma }·\sqrt{\frac{{\chi }_{1-\alpha /2,n_k-1}^2}{n_k-1}}.\hfill \end{array}$$

In the above, one construct the S chart with probability limits by using ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$ or ${\overline{S}}_D$ instead of $\widehat{\sigma }$.

Next, we consider the construction of the $S^2$ chart. We rewrite Equation (37) about $S_k^2$ and we then have

$$P\left[{\sigma }^2·\frac{{\chi }_{1-\alpha /2,n_k-1}^2}{n_k-1}\le S_k^2\le {\sigma }^2·\frac{{\chi }_{\alpha /2,n_k-1}^2}{n_k-1}\right]=1-\alpha .$$

Using the above along with ${\widehat{\sigma }}^2=S_p^2$ where $S_p^2$ is the pooled sample variance denoted in Equation (9), we can also construct the $S^2$ chart as follows:

$$\begin{array}{cc}\hfill \mathrm{UCL}\hfill & =S_p^2·\frac{{\chi }_{\alpha /2,n_k-1}^2}{n_k-1}\hfill \\ \hfill \mathrm{CL}\hfill & =S_p^2\hfill \\ \hfill \mathrm{LCL}\hfill & =S_p^2·\frac{{\chi }_{1-\alpha /2,n_k-1}^2}{n_k-1}.\hfill \end{array}$$

Note that none of ${\overline{S}}_A^2$, ${\overline{S}}_B^2$, ${\overline{S}}_C^2$, or ${\overline{S}}_D^2$ is unbiased for ${\sigma }^2$, whereas $S_p^2$ is unbiased. Thus, it is not recommended to use any of the four estimators to construct the $S^2$ chart.

5.3. The $\overline{X}$ Chart

From the statistical asymptotic theory, we have

$$\frac{{\overline{X}}_k-E\left({\overline{X}}_k\right)}{\mathrm{SE}\left({\overline{X}}_k\right)}\stackrel{•}{\sim }N(0,1),$$

where ${\overline{X}}_k$ is the sample mean with sample size $n_k$. In order to construct the $\mathrm{CL}\pm 3·\mathrm{SE}$ control limits, we can set up $({\overline{X}}_k-E\left({\overline{X}}_k\right))/\mathrm{SE}\left({\overline{X}}_k\right)=\pm 3$. Solving this for ${\overline{X}}_k$, we have

$$E\left({\overline{X}}_k\right)\pm 3·\mathrm{SE}\left({\overline{X}}_k\right)=\mu \pm \frac{3\sigma }{\sqrt{n_k}}.$$

Since $\mu$ and $\sigma$ are unknown in practice, we need to estimate them. With the estimates $\widehat{\mu }$ and $\widehat{\sigma }$, we have

$$\begin{array}{cc}\hfill \mathrm{UCL}\hfill & =\widehat{\mu }+\frac{3\widehat{\sigma }}{\sqrt{n_k}}\hfill \\ \hfill \mathrm{CL}\hfill & =\widehat{\mu }\hfill \\ \hfill \mathrm{LCL}\hfill & =\widehat{\mu }-\frac{3\widehat{\sigma }}{\sqrt{n_k}}.\hfill \end{array}$$

To estimate $\mu$, one can use ${\overline{\overline{X}}}_A$ defined in Equation (1) or ${\overline{\overline{X}}}_B$ defined in Equation (2). To estimate $\sigma$, one can use any of ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, and ${\overline{S}}_D$.

As an illustration, using ${\overline{\overline{X}}}_B$ and ${\overline{S}}_A$, we have $\overline{X}$ chart as follows:

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_A\hfill & ={\overline{\overline{X}}}_B+\frac{3{\overline{S}}_A}{\sqrt{n_k}}\hfill \\ \hfill {\mathrm{CL}}_A\hfill & ={\overline{\overline{X}}}_B\hfill \\ \hfill {\mathrm{LCL}}_A\hfill & ={\overline{\overline{X}}}_B-\frac{3{\overline{S}}_A}{\sqrt{n_k}}.\hfill \end{array}$$

Similarly, we can construct various $\overline{X}$ charts using a total of eight combinations of $\mu$ and $\sigma$ estimators.

Remark 7.

It should be noted that, when ${\overline{S}}_D$ is used, the $\overline{X}$ chart is somewhat different from the above chart and is given by

$$\begin{array}{cc}\hfill {\mathrm{UCL}}_D\hfill & =\overline{\overline{X}}+\frac{3S_p}{c_4(nm-m+1)\sqrt{n_k}}\hfill \\ \hfill {\mathrm{CL}}_D\hfill & =\overline{\overline{X}}\hfill \\ \hfill {\mathrm{LCL}}_D\hfill & =\overline{\overline{X}}-\frac{3S_p}{c_4(nm-m+1)\sqrt{n_k}}.\hfill \end{array}$$

As shown in Section 4, ${\overline{S}}_D$ performs better than ${\overline{S}}_A$, ${\overline{S}}_B$, and ${\overline{S}}_C$. However, to the best of our knowledge, the $\overline{X}$ chart based on ${\overline{S}}_D$ has not been widely used probably because the calculation of the normal-consistent unbiasing factor $c_4$ is difficult especially with a large argument value. Most textbooks provide the values of $c_4\left(n\right)$ only for $n\le 25$. In Appendix A, for an easy calculation, we provide an approximation for $c_4\left(n\right)$ which is highly accurate within one unit in the ninth decimal place for $n>25$. Thus, using this, one can easily calculate the LCL and UCL of the $\overline{X}$ chart based on ${\overline{S}}_D$.

6. Average Run Length and Standard Deviation of Run Length

To compare the performance of the control charts based on the various estimators, we obtained the empirical estimates of the ARL and the SDRL through using the extensive Monte Carlo simulation method. For this simulation, we considered the $\overline{X}$ chart. In this chart, we only estimated the location with $\widehat{\mu }={\overline{\overline{X}}}_B$ because the RE of ${\overline{\overline{X}}}_B$ is better than that of ${\overline{\overline{X}}}_A$ as shown in Section 4. For the scale, we considered seven different estimates (${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, $\overline{S}$, ${\overline{S}}^{*}$, ${\overline{S}}_w$) and we denoted the charts based on these scale estimates by A, B, C, D, S, S*, and Sw, respectively.

We assume that we have m samples ($n_1,n_2,\dots ,n_m$) in Phase-I. Again, let $X_{ij}$ be the ith sample (subgroup) of size $n_i$. Then, $X_{ij}$’s were generated from the normal distribution with mean ${\mu }_0=10$ and standard deviation ${\sigma }_0=5$ and we obtained the location estimate (${\overline{\overline{X}}}_B$) and the seven scale estimates. Using these estimates, we constructed the seven control charts based on $\mathrm{CL}\pm 3·\mathrm{SE}$ control limits with FAR 0.27%. Then, we monitored the process with a new sample of size $n_k$ from the same normal distribution and obtained the run length in Phase-II. We repeated this simulation one million times ($I={10}^6$) to obtain the run lengths and then estimated the ARL and SDRL based on these run lengths. Note that the simulation results are the same as the ones under different parameter values of ${\mu }_0$ and ${\sigma }_0$. It deserves mentioning that the results are quite reasonable, since the normal distribution is a location-scale family, whereas the results are somewhat dependent on the number of samples and the combination of the sample sizes.

We generated $m=15$ samples with the five different scenarios. The sample sizes of each scenario are given by

$$\mathrm{Scenario}\mathrm{I}{n}_1=n_2=\cdots =n_5=3,n_6=n_7=\cdots =n_{10}=10,n_{11}=n_{12}=\cdots =n_{15}=17,$$

$$\mathrm{Scenario}\mathrm{II}{n}_1=n_2=\cdots =n_5=5,n_6=n_7=\cdots =n_{10}=10,n_{11}=n_{12}=\cdots =n_{15}=15,$$

$$\mathrm{Scenario}\mathrm{III}{n}_1=n_2=\cdots =n_5=7,n_6=n_7=\cdots =n_{10}=10,n_{11}=n_{12}=\cdots =n_{15}=13,$$

$$\mathrm{Scenario}\mathrm{IV}{n}_1=n_2=\cdots =n_5=9,n_6=n_7=\cdots =n_{10}=10,n_{11}=n_{12}=\cdots =n_{15}=11,$$

$$\mathrm{Scenario}\mathrm{V}{n}_1=n_2=\cdots =n_5=10,n_6=n_7=\cdots =n_{10}=10,n_{11}=n_{12}=\cdots =n_{15}=10,$$

We considered Scenarios I–IV for the cases of unequal sample sizes and Scenario V for the case of an equal sample size as a reference. Upon each of these scenarios, we estimated the ARL and SDRL as described above. The simulation results are provided in

Table 3. Estimated ARL and SDRL with $n_k=10$.

	A	B	C	D	$\overline{\mathit{S}}$	${\overline{\mathit{S}}}^{*}$	${\overline{\mathit{S}}}_{\mathit{w}}$
		Scenario I
ARL	475.03	456.02	363.61	361.84	257.78	343.39	270.79
SDRL	1301.18	1184.11	536.61	531.45	499.21	777.03	387.18
		Scenario II
ARL	390.41	388.19	364.59	362.56	269.79	357.12	274.27
SDRL	654.28	643.59	540.99	533.90	421.37	586.10	391.26
		Scenario III
ARL	370.63	370.25	363.84	361.77	273.95	361.84	275.39
SDRL	565.01	563.36	538.59	530.64	400.49	549.63	392.54
		Scenario IV
ARL	364.28	364.23	363.63	361.77	275.28	363.39	275.32
SDRL	539.37	539.16	537.54	531.26	393.40	537.95	392.44
		Scenario V (equal sample size)
ARL	364.36	364.36	364.36	362.58	275.95	364.36	275.95
SDRL	541.45	541.45	541.45	537.31	395.12	541.45	395.12

.

Note that the ideal ARL under the normal distribution is around $1/0.0027\approx 370$ with FAR 0.27%. However, this ideal value is obtained when using the true ${\mu }_0$ and ${\sigma }_0$. In practice, we need to estimate these parameters in Phase-I. Thus, with such an uncertainty due to estimation, the target ARL can be different from the ideal ARL of 370. We observed that the empirical results of A, B, C, and S* have the same value under the case of an equal sample size (Scenario V) and these results have the same tendency when the RE measure is used in Section 4.2. These results are expected as pointed out in Remark 3. The estimated ARL of D (362.58) is very close to that of A, B, C, and S* (364.36). In addition, the estimated SDRL of D (537.31) is very close to that of A, B, C, and S* (541.45). This minor difference may be due to a random phenomenon of the Monte Carlo simulation. Thus, it is quite reasonable to assume that the target ARL is around 360 and the target SDRL is around 540.

In what follows, we analyze and compare our results based on these target values. Note that we did not consider the values from the control charts based on S and Sw because  they have much smaller ARL values than the others, which is mainly due to the fact that they have serious negative bias even with the case of an equal sample size as seen in Table 1.

In Scenario I, there is a serious difference in terms of the sample sizes (3, 10, 17). The ARL values of the charts based on A and B seriously overshoot the target value and their SDRLs are far above the target value. The ARL values of the charts using S and Sw seriously undershoot the target value while the ARL of the chart using S* has a minor underestimate. These results are closely related to the RE. In Table 2, $\overline{S}$ and ${\overline{S}}_w$ have noticeable negative values of the biases while ${\overline{S}}^{*}$ has a decent negative value. This implies that the scale estimates underestimate the true value, which results in a narrower control chart. Thus, using a narrower control chart, one can have a smaller ARL. When it comes to the SDRL, the SDRLs of the charts using A and B seriously overshoot the target and those using S and Sw undershoot the target. The chart using S* overshoots the target. These results are similar to those with the RE.

In Scenarios II–IV, we have similar observations as noticed in Scenario I. Because there is a less serious difference in sample sizes, the observations are not so dramatic as seen in Scenario I, but their tendencies are quite similar. For example, in Scenario IV, the sample sizes are very slightly different, so these results are quite close to those in Scenario V (equal sample size).

In Scenario V, as mentioned earlier, we considered this as a reference. In this case, we can also observe that the charts using S and Sw have the same results and their ARL values seriously undershoot the target value. In Remark 4, we pointed out ${\overline{S}}_w=\overline{S}$ for the case of an equal sample size. We can also observe this in Figure 1a. As mentioned earlier, their underperformance is due to a serious negative bias as seen in Table 1. The results show that the charts based on C and D perform very well. However, when the samples are very small (small number of samples with small sample size), we expect that the chart based on D can be slightly better than that based on C as shown in Figure 1. However, in practice, with a decent size of samples, one can use the charts based on either C or D.

7. Illustrative Examples

Here, we provide three real-data examples to illustrate the applications of the proposed methods into the control charts. All computations were analyzed using the R language [41,42]. The R functions for the $\overline{X}$ and S charts can be obtained in Appendix B.

Example 1.

We consider the data set presented earlier in Table 30 in Section 3.31 of ASTM [10]. The data sets were obtained from ten shipments whose sample sizes were not equal. The sample sizes are 50, 50, 100, 25, 25, 50, 100, 50, 50, 50. The corresponding sample means are given by 55.7, 54.6, 52.6, 55.0, 53.4, 55.2, 53.3, 52.3, 53.7, 54.3 and the corresponding sample standard deviations are 4.35, 4.03, 2.43, 3.56, 3.10, 3.30, 4.18, 4.30, 2.09, 2.67, respectively. Using these, we have ${\overline{\overline{X}}}_A=54.01$, ${\overline{\overline{X}}}_B=53.8$, ${\overline{S}}_A=3.420251$, ${\overline{S}}_B=3.420254$, ${\overline{S}}_C=3.405517$, and  ${\overline{S}}_D=3.491055$.

Using the R functions provided in Appendix B, we can obtain the following results. In the R function, A, B, C, and D denote the control limits based on ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, respectively. The R functionXbarchart()uses ${\overline{\overline{X}}}_B$ as a default for the CL of the $\overline{X}$ chart since ${\overline{\overline{X}}}_B$ performs better as seen in Section 4.1.

> ni = c(50, 50, 100, 25, 25, 50, 100, 50, 50, 50)

> Xbari = c(55.7, 54.6, 52.6, 55.0, 53.4, 55.2, 53.3, 52.3, 53.7, 54.3)

> Si = c(4.35, 4.03, 2.43, 3.56, 3.10, 3.30, 4.18, 4.30, 2.09, 2.67)

> Xbarchart(Xbari, Si, ni=ni, nk=25)

LCL CL  UCL

A 51.74785 53.8 55.85215

B 51.74785 53.8 55.85215

C 51.75669 53.8 55.84331

D 51.70537 53.8~55.89463

> Schart(Si, ni, nk=25)

LCL  CL  UCL

A 1.911697 3.384818 4.857940

B 1.911699 3.384822 4.857945

C 1.903462 3.370238 4.837013

D 1.951272 3.454889 4.958505

The control limits are calculated using Methods A–D and the results are summarized in

Table 4. Control limits for the $\overline{X}$ and S charts for Example 1.

		$\overline{\mathit{X}}$ Chart				S Chart
Method		LCL	CL	UCL		LCL	CL	UCL
				sample size, $n_k=25$
A		51.74785	53.80	55.85215		1.911697	3.384818	4.857940
B		51.74785	53.80	55.85215		1.911699	3.384822	4.857945
C		51.75669	53.80	55.84331		1.903462	3.370238	4.837013
D		51.70537	53.80	55.89463		1.951272	3.454889	4.958505
				sample size, $n_k=50$
A		52.34891	53.80	55.25109		2.369028	3.402846	4.436665
B		52.34891	53.80	55.25109		2.369030	3.402850	4.436669
C		52.35516	53.80	55.24484		2.358823	3.388188	4.417553
D		52.31887	53.80	55.28113		2.418070	3.473290	4.528509
				sample size, $n_k=100$
A		52.77392	53.80	54.82608		2.683351	3.411625	4.139899
B		52.77392	53.80	54.82608		2.683354	3.411629	4.139903
C		52.77834	53.80	54.82166		2.671792	3.396929	4.122065
D		52.75268	53.80	54.84732		2.738900	3.482250	4.225600

. The results using Methods A–C are quite close but those using Method D are slightly different from others.

In addition, the variances of ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, and ${\overline{S}}_E$ are obtained from Equations (14)–(16), (18), and (34), respectively, and we have $\mathrm{Var}\left({\overline{S}}_A\right)=0.0011375146·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_B\right)=0.0011348232·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_C\right)=0.0009301593·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_D\right)=0.0009263542·{\sigma }^2$, and  $\mathrm{Var}\left({\overline{S}}_E\right)=0.0009111612·{\sigma }^2$. Thus, the REs with respect to ${\overline{S}}_E$ are easily obtained from Equation (35) and we have $\mathrm{RE}\left({\overline{S}}_A\right)=80.10\%$, $\mathrm{RE}\left({\overline{S}}_B\right)=80.29\%$, $\mathrm{RE}\left({\overline{S}}_C\right)=97.96\%$, and  $\mathrm{RE}\left({\overline{S}}_D\right)=98.36\%$.

Example 2.

We consider the data set presented earlier in Table 32 in Section 3.31 of ASTM [10]. The data sets were obtained from 21 tension testing machines whose sample sizes were not equal. The sample sizes are 5, 5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5. The corresponding sample means are given by 73.8, 71.0, 74.2, 71.0, 70.0, 67.0, 73.5, 71.2, 71.2, 71.2, 71.6, 71.2, 74.2, 74.6, 72.4, 75.3, 69.0, 71.8, 72.8, 69.8, 69.00, and the corresponding sample standard deviations are 1.10, 0.71, 0.45, 1.41, 0.00, 2.35, 1.91, 1.79, 0.45, 0.45, 0.55, 0.55, 0.84, 0.55, 0.55, 0.50, 0.71, 0.84, 0.45, 1.30, 0.00, respectively. Using these, we have ${\overline{\overline{X}}}_A=71.70476$, ${\overline{\overline{X}}}_B=71.65243$, ${\overline{S}}_A=0.8869858$, ${\overline{S}}_B=0.8861882$, ${\overline{S}}_C=0.8762927$, and  ${\overline{S}}_D=1.014672$.

Using the R functions in Appendix B with the above data sets, the control limits are easily calculated and we summarized the results in

Table 5. Control limits for the $\overline{X}$ and S charts for Example 2.

		$\overline{\mathit{X}}$ Chart				S Chart
Method		LCl	CL	UCL		LCL	CL	UCL
				sample size, $n_k=4$
A		70.32195	71.65243	72.98291		0	0.8171958	1.851804
B		70.32314	71.65243	72.98171		0	0.8164609	1.850139
C		70.33799	71.65243	72.96687		0	0.8073440	1.829480
D		70.13042	71.65243	73.17444		0	0.9348355	2.118381
				sample size, $n_k=5$
A		70.46241	71.65243	72.84244		0	0.8337539	1.741710
B		70.46348	71.65243	72.84137		0	0.8330041	1.740144
C		70.47676	71.65243	72.82810		0	0.8237026	1.720713
D		70.29110	71.65243	73.01375		0	0.9537773	1.992439

.

As was illustrated in Example 1, the variances of ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, and ${\overline{S}}_E$ are obtained from Equations (14)–(16), (18), and (34), respectively, and we have $\mathrm{Var}\left({\overline{S}}_A\right)=0.006484797·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_B\right)=0.006477515·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_C\right)=0.006434091·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_D\right)=0.006116037·{\sigma }^2$, and  $\mathrm{Var}\left({\overline{S}}_E\right)=0.004913916·{\sigma }^2$. Then, the REs with respect to ${\overline{S}}_E$ are given by $\mathrm{RE}\left({\overline{S}}_A\right)=75.78\%$, $\mathrm{RE}\left({\overline{S}}_B\right)=75.86\%$, $\mathrm{RE}\left({\overline{S}}_C\right)=76.37\%$, and  $\mathrm{RE}\left({\overline{S}}_D\right)=80.34\%$.

Example 3.

We consider the data set presented earlier in Table 6.4 of Montgomery [9] which includes 113 measurements (in millimeters) of the diameters of piston rings for an automotive engine produced by a forging process. The data were obtained from 25 samples and the sample sizes are given by 5, 3, 5, 5, 5, 4, 4, 5, 4, 5, 5, 5, 3, 5, 3, 5, 4, 5, 5, 3, 5, 5, 5, 5, 5. The corresponding sample means are calculated as 74.010, 73.996, 74.008, 74.003, 74.003, 73.996, 73.999, 73.997, 74.004, 73.998, 73.994, 74.001, 73.994, 73.990, 74.008, 73.997, 73.999, 74.007, 73.998, 74.008, 74.000, 74.002, 74.002, 74.005, 73.998, and the corresponding sample standard deviations are 0.0148, 0.0046, 0.0147, 0.0091, 0.0122, 0.0099, 0.0055, 0.0123, 0.0064, 0.0063, 0.0029, 0.0042, 0.0100, 0.0153, 0.0087, 0.0078, 0.0115, 0.0070, 0.0085, 0.0068, 0.0122, 0.0074, 0.0119, 0.0087, 0.0162. Using these, we have ${\overline{\overline{X}}}_A=74.00068$, ${\overline{\overline{X}}}_B=74.00066$, ${\overline{S}}_A=0.01010231$, ${\overline{S}}_B=0.01012067$, ${\overline{S}}_C=0.01030545$, and  ${\overline{S}}_D=0.01032266$.

Similar to the two examples above, we can obtain the control limits using these data sets. The control limits are calculated and summarized in

Table 6. Control limits for the $\overline{X}$ and S charts for Example 3.

		$\overline{\mathit{X}}$ Chart				S Chart
Method		LCL	CL	UCL		LCL	CL	UCL
				sample size, $n_k=3$
A		73.98317	74.00066	74.01816		0	0.008952937	0.02299266
B		73.98313	74.00066	74.01819		0	0.008969207	0.02303445
C		73.98281	74.00066	74.01851		0	0.009132968	0.02345501
D		73.98278	74.00066	74.01854		0	0.009148216	0.02349417
				sample size, $n_k=4$
A		73.98551	74.00066	74.01582		0	0.009307435	0.02109109
B		73.98548	74.00066	74.01584		0	0.009324349	0.02112941
C		73.98521	74.00066	74.01612		0	0.009494595	0.02151520
D		73.98518	74.00066	74.01615		0	0.009510446	0.02155112
				sample size, $n_k=5$
A		73.98711	74.00066	74.01422		0	0.009496024	0.01983717
B		73.98709	74.00066	74.01424		0	0.009513281	0.01987322
C		73.98684	74.00066	74.01449		0	0.009686976	0.02023607
D		73.98681	74.00066	74.01451		0	0.009703148	0.02026986

. In addition, the variances of ${\overline{S}}_A$, ${\overline{S}}_B$, ${\overline{S}}_C$, ${\overline{S}}_D$, and ${\overline{S}}_E$ are obtained as $\mathrm{Var}\left({\overline{S}}_A\right)=0.006472658·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_B\right)=0.006390116·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_C\right)=0.006020000·{\sigma }^2$, $\mathrm{Var}\left({\overline{S}}_D\right)=0.005697867·{\sigma }^2$, and  $\mathrm{Var}\left({\overline{S}}_E\right)=0.004474206·{\sigma }^2$. Using these, the REs with respect to ${\overline{S}}_E$ are calculated as $\mathrm{RE}\left({\overline{S}}_A\right)=69.12\%$, $\mathrm{RE}\left({\overline{S}}_B\right)=70.02\%$, $\mathrm{RE}\left({\overline{S}}_C\right)=74.32\%$, and  $\mathrm{RE}\left({\overline{S}}_D\right)=78.52\%$.

8. Conclusions

In this paper, we have considered several unbiased location and scale estimators for the process parameters of the $\overline{X}$ and S control charts when the sample sizes are not necessarily equal. These estimators are essential for constructing the control limits of the Shewhart-type control charts. A natural question is: among these unbiased estimators, which one should be recommended in practical applications? We clarified this question by providing the inequality relations among the variances of these estimators through the rigorous proofs. We also showed that the conventional ad hoc methods could result in degraded performance of the control charts, mainly because the adopted estimators are all biased and they actually tend to underestimate the true scale parameter.

We also provided the relative efficiency of the methods along with the conventional methods and the empirical estimates of the ARL and the SDRL through using the extensive Monte Carlo simulations. We observed that the chart based on ${\overline{S}}_D$ outperforms the others under consideration from both theoretical and numerical points of view. The only difficulty of using ${\overline{S}}_D$ lies in calculating the normal-consistent unbiasing factor $c_4$ for a large argument value (large sample size) without using a professional software. To resolve this problem, we provided an approximation of the $c_4$ as a function of a sample size which can be easily calculated with a general calculator for the case of a large argument value and is accurate within one unit in the ninth decimal place.

All of the theoretical results revealed an interesting and useful connection between the two fields of statistics and mathematics. For example, the normal-consistent unbiasing factor $c_4$ can be expressed as a function through using the Watson representation, which helps one to understand the behavior of the $c_4$ in depth. We thus expect that the new findings about the $c_4$ can help quality engineers develop more useful results in various engineering statistics fields.

It is noteworthy to mention that all the theoretical and empirical results of this paper require the assumption of the homoscedasticity case of $\sigma$. In an ongoing work, we investigate these estimators more thoroughly for the heteroscedasticity case.