A Bivariate Beta from Gamma Ratios for Determining a Potential Variance Change Point: Inspired from a Process Control Scenario

Abstract

Within statistical process control (SPC), normality is often assumed as the underlying probabilistic generator where the process variance is assumed equal for all rational subgroups. The parameters of the underlying process are usually assumed to be known—if this is not the case, some challenges arise in the estimation of unknown parameters in the SPC environment especially in the case of few observations. This paper proposes a bivariate beta type distribution to guide the user in the detection of a permanent upward or downward step shift in the process’ variance that does not directly rely on parameter estimates, and as such presents itself as an attractive and intuitive approach for not only potentially identifying the magnitude of the shift, but also the position in time where this shift is most likely to occur. Certain statistical properties of this distribution are derived and simulation illustrates the theoretical results. In particular, some insights are gained by comparing the newly proposed model’s performance with an existing approach. A multivariate extension is described, and useful relationships between the derived model and other bivariate beta distributions are also included.

1. Introduction

1.1. Problem Contextualisation

The monitoring of the variance of independent and identically distributed (i.i.d) normal random variables over time by taking successive, independent samples of measurements over time remains an interesting and valuable research consideration within quality control environment. In this case, when the variance ${\sigma }^2$ changes to ${\sigma }_1^2=\lambda {\sigma }^2$ for some $\lambda \ne 1$, the practitioner needs to investigate the scope of such a change (a value of $\lambda$), and ideally, the position within the successive measurements where such a change could’ve taken place. Suppose that $X_{ij}$ are i.i.d. $N(\mu ,{\sigma }^2)$, $i=0,1,2,\cdots ,\kappa -1$ and $X_{ij}\sim$i.i.d. $N\left(\mu ,{\sigma }_1^2=\lambda {\sigma }^2\right),i=\kappa ,\kappa +1,\cdots ,m$ where $j=1,2,\cdots ,n_i\ge 2$ and $\lambda >0$, as outlined in Figure 1. The values of $\kappa$ and $\lambda$ are assumed to be unknown, but deterministic in nature. The order of these samples is important and cannot be re-ordered; in other words, the samples have a set sequence corresponding to the order in which they were obtained.

Thus, inspired by a practical objective, this paper aims to present a theoretically motivated framework to

• present and contextualise this problem within the quality control environment;
• follow a systematic approach to build up the distributional foundations from this practical perspective;
• exploratively focus on the development of the (new) resulting bivariate beta distribution;
• compare this model with an existing approach; and
• determine whether $\lambda \ne 1,$ and if this is indeed the case, to determine $\kappa$, the location of where the shift in the variance occurred.

Therefore, from sample $\kappa$ onwards, the process would be considered out-of-control. Note that it is assumed that the shift occurs between two samples.

Assume that both the process mean $\left(\mu \right)$ and variance $\left({\sigma }^2\right)$ are unknown, and that they are estimated by their respective minimum variance unbiased estimators (MVUE), given by:

$$\begin{array}{ccc}\hfill {\overline{X}}_i& =& \frac{{\sum }_{j=1}^{n_i}{X}_{ij}}{n_i},i=0,1,2,\cdots ,m,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill S_i^2& =& \frac{1}{n_i-1}\sum _{j=1}^{n_i}{\left(X_{ij}-{\overline{X}}_i\right)}^2,i=0,1,2,\cdots ,m,\hfill \end{array}$$

(2)

where ${\overline{X}}_i$ and $S_i^2$ denote the mean and variance of sample i respectively. Some particular notes on Equation (1) and (2) include:

• The index variable i ranges from 0 to m: a total of $m+1$ independent rational subgroups or samples.
• At least two samples are needed for a potential shift between them to be possible, therefore we assume $m\ge 1$.
• The sample size $n_i$ can vary between different samples.
• $n_i\ge 2$ is necessary since the process mean and variance are both assumed to be unknown and have to be estimated.
• The pooling approach here is to use $m-r+1$ and r sample means and variances in the construction of the test statistic in Section 1.2. Alternatively one can consider a single mean/variance in this construction, which would result in additional information ${\sum }_{i=r}^m{n}_i-1$ and ${\sum }_{i=0}^{r-1}{n}_i-1$ such that $n_i\ge 1$ and probability density functions are valid. In this case, the approach would reduce to a two sample comparison testing for a change in the variance.

The problem of determining if a shift in the process variance has occurred can be divided into two stages, namely before the potential shift and after, as indicated below ($Gamma(·,·)$ denotes the usual gamma distribution with suitable shape and scale parameters [1]).

1.2. A Solution: Sequential Statistic Framework

The proposed distribution compares all the sample variances before a certain point (where the potential shift occurs), with all sample variances after the time of the shift. In essence, the multi-sample hypothesis testing problem is approached by using m sequential two sample tests; described as:

$$\begin{array}{ccc}{S}_0^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_1^2,S_2^2,\cdots ,S_m^2\\ S_0^2,S_1^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_2^2,S_3^2,\cdots ,S_m^2\\ & \mathrm{and}\mathrm{so}\mathrm{forth}\mathrm{until}\\ S_0^2,S_1^2,\cdots ,S_{m-1}^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_m^2.\end{array}$$

(5)

Note that, due to the assumption that we have $m+1$ samples, the procedure to identify the possible change in the variance requires m comparisons. Of these m comparisons the one that leads to the largest disparity between the sample variables on the left and the right of (5) will indicate the most likely position where the process experienced a change in variance. In essence, what is then needed is to quantify the difference between the sample variances of the left and right of (5), and then determine which of the m different comparisons has the largest difference between the sample variances, and finally to use some measure to determine if this maximum difference is within some set tolerance range. Our proposed mathematical construct on how to achieve this is discussed for the remained of this introduction. Assuming that no shift in the process variance has occurred, it is possible to construct a series of two sample statistics that correspond to the general procedure described in (5). Each statistic corresponds to whether at sample $r=\kappa$ the two independent samples (the sample variances before time r and the sample variances after and including time r) are from normal distributions with the same unknown variance ${\sigma }^2$. This can alternatively be viewed as testing whether ${\sigma }^2={\sigma }_1^2$, which is similar to testing $\lambda =1$. As such, it follows that detecting a shift in the process variance can be reduced to the following hypothesis test:

$$H_0:{\sigma }^2={\sigma }_1^2\mathrm{vs}.H_A:{\sigma }^2\ne {\sigma }_1^2\mathrm{or}\mathrm{alternatively}{H}_0:\lambda =1\mathrm{vs}.H_A:\lambda \ne 1.$$

Suppose that a shift of size $\lambda$ occurs in the process variance, then the corresponding random variables after the shift would be $W_i$∼$Gamma(\frac{n_i-1}{2},2\lambda ),i=\kappa ,\kappa +1,\cdots ,m$ distributed (see (4)). If there is no shift in the variance, $\lambda =1$, and it follows that $W_i$∼$Gamma(\frac{n_i-1}{2},2),i=0,1,\cdots ,m$. Thus the hypothesis being investigated can be changed depending on the choice of the scale parameter of the specific gamma random variables. From (5) it follows that the series of statistics that forms the building blocks of the process are given by

$$U_r^{*}=\frac{\left(\frac{{\sum }_{i=r}^m\left(n_i-1\right)S_i^2}{\lambda {\sigma }^2{\sum }_{i=r}^m\left(n_i-1\right)}\right)}{\left(\frac{{\sum }_{i=0}^{r-1}\left(n_i-1\right)S_i^2}{{\sigma }^2{\sum }_{i=0}^{r-1}\left(n_i-1\right)}\right)}\equiv \frac{\frac{{\sum }_{i=r}^m{W}_i}{{\sum }_{i=r}^m\left(n_i-1\right)}}{\frac{{\sum }_{i=0}^{r-1}{W}_i}{{\sum }_{i=0}^{r-1}\left(n_i-1\right)}},r=1,2,\cdots ,m-1,m,$$

(6)

where $W_i$∼$Gamma(\frac{n_i-1}{2},2)$ for $i=1,2,\cdots ,r-1$ and $W_i$∼$Gamma(\frac{n_i-1}{2},2\lambda )$ for $i=r,r+1,\cdots ,m$.

In essence, (6) implies that the sample variances before the potential shift are pooled together, and the sample variances after the potential shift are pooled together. The numerator of the statistic at time r is the average, weighted by each statistic’s degrees of freedom, of all the sample variances between and including samples r and m, while the denominator is the corresponding weighted average of all the sample variances between and including samples 0 and $r-1$; this is graphically presented in Figure 2. Thus $U_r^{*}$ is the test statistic that is typically used to test the equality of two population variances from a normal distribution.

From (3), (3) and (6), it follows that if no shift has occurred in the process variance, i.e., $\lambda =1$, each statistic $U_r^{*}$ is univariate F distributed with ${\sum }_{i=r}^m\left(n_i-1\right)$ and ${\sum }_{i=0}^{r-1}\left(n_i-1\right)$ degrees of freedom, respectively.

The reasoning behind the critical values that would indicate whether $\lambda \ne 1$ is justified by inspecting the sequence of statistics in (6). Suppose that an increase $\left(\lambda >1\right)$ in the process variance has occurred from sample $r=\kappa$ onward, then:

• The statistic $U_r^{*}$’s numerator will contain only sample variances that come from a $N\left(\mu ,\lambda {\sigma }^2\right),\lambda >1$ distribution, whereas the denominator will contain only sample variances that come from a $N\left(\mu ,{\sigma }^2\right)$ distribution.
• If $k_1$ is some integer value such that $1\le r-k_1<r$, then the statistic $U_{r-k_1}^{*}$ will contain $k_1$ sample variances in its numerator that are from a $N\left(\mu ,{\sigma }^2\right)$ distribution. This will reduce the weighted average of the sample variances of $U_{r-k_1}^{*}$’s numerator in comparison to the numerator of $U_r^{*}$.
• Similarly, if $k_2$ is some integer value such that $r<r+k_2\le m$, then the statistic $U_{r+k_2}^{*}$ will contain $k_2$ sample variances in its denominator that are from a $N\left(\mu ,\lambda {\sigma }^2\right)$ distribution. This will increase the weighted average of sample variances of $U_{r+k_2}^{*}$’s denominator in comparison to the denominator of $U_r^{*}$.
• Thus, any statistic other than the one immediately following the shift in the process variance, will contain either smaller sample variances in its numerator, or larger sample variances in its denominator (on average). Either of these scenarios result in a high probability that all other statistics are smaller relative to $U_r^{*}$.
• This leads to the conclusion that the most probable place where an upwards shift in the process variance will be detected is at the statistic immediately following the shift. The value that this statistic assumes also has a high likelihood of being the maximum value of all the $U_r^{*},r=1,2,\cdots ,m-1,m$ statistics.
• As such, the most reasonable method of calculating the critical value to detect an upwards shift in the process variance is to calculate the maximum order statistic of the charting statistics $U_r^{*},r=1,2,\cdots ,m-1,m$, (under the null hypothesis) and to set the critical value equal to some percentile of the distribution of the maximum order statistic.

Using a similar but inverted argument, it can be justified that the critical value of the control chart should be set equal to some percentile of the minimum order statistic of the charting statistics, under the null hypothesis of no shift having occurred, if the detection of a downward shift in the process variance is of concern. Due to space limitations, the minimum order statistics are not presented in this article.

To simplify matters going forward and for notational purposes we omit the factors ${\sum }_{i=r}^m\left(n_i-1\right)$ and ${\sum }_{i=0}^{r-1}\left(n_i-1\right)$ in (6), and drop the superscript *, and therefore the statistics of interest become

$$U_r=\frac{{\sum }_{i=r}^m{W}_i}{{\sum }_{i=0}^{r-1}{W}_i},r=1,2,\cdots ,m-1,m,$$

(7)

where $W_i\sim Gamma\left({\alpha }_i>0,{\beta }_i>0\right),i=0,1,\cdots ,m$ and independent. Since ${\alpha }_i=\frac{n_i-1}{2}$ and ${\beta }_i=2\lambda$, the shape parameter is related to the sample size of the $i_{th}$ sample, and the scale parameter is related to the underlying distribution’s variance. The theoretical focus here is based on the statistics in (7), whereas the statistics in (6) are those that are practically applicable and is the basis in the simulation study in Section 3.

Constructing statistics using ratios of random variables as in (7) is of practical interest in many areas of science. Ref. [2] studied and derived the joint density functions of ratios of Rayleigh, Rician, Nakagami-m, and Weibull random variables; [3] approached the ratios of generalised gamma variables via exact- and near exact solutions, and [4] derived closed-form expressions for the ratio of independent non-identically distributed variables from an $\alpha$-$\mu$ distribution which have applications in the performance analysis of wireless communication systems.

The proposed model in this paper will be compared to the model of [5] in Section 3, and is described here for the convenience of the reader. If $r=2$, that the bivariate joint probability density function of the statistics in (10) is given by

$$\begin{array}{ccc}g\left(t_1,t_2\right)\hfill & =\hfill & \frac{\left({\beta }_0^{{\alpha }_1+{\alpha }_2}{\beta }_1^{{\alpha }_0+{\alpha }_2}{\beta }_2^{{\alpha }_0+{\alpha }_1}\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left(t_1\right)}^{{\alpha }_1-1}{\left(t_2\right)}^{{\alpha }_2-1}{\left(1+t_1\right)}^{{\alpha }_2}\hfill \\ \hfill & \times \hfill & {\left({\beta }_1{\beta }_2+{\beta }_0{\beta }_2{t}_1+{\beta }_0{\beta }_1\left(1+t_1\right)t_2\right)}^{-{\alpha }_0-{\alpha }_1-{\alpha }_2}\hfill \end{array}$$

(8)

where $t_1,t_2>0$, ${\alpha }_i,{\beta }_i>0$ for $i=0,1,2,$, $t_1=\frac{w_1}{w_0}$ and $t_2=\frac{w_2}{(w_0+w_1)}$, and $\mathsf{\Gamma }(·)$ is the gamma function [6].

Refs. [5,7] proposed a beta distribution that is used to detect a shift in the process variance which is based on the Q chart that was developed by [8], and as such the series of comparisons of sample statistics they made were

$$\begin{array}{ccc}{S}_1^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_0^2\\ S_2^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_0^2,S_1^2\\ & \mathrm{and}\mathrm{so}\mathrm{forth}\mathrm{until}\\ S_m^2& \mathrm{is}\mathrm{compared}\mathrm{with}& S_0^2,S_1^2,\cdots ,S_{m-1}^2.\end{array}$$

(9)

Using a similar approach to the method described earlier the statistics can be given as

$$T_r=\frac{W_r}{{\sum }_{i=0}^{r-1}{W}_i},r=1,2,\cdots ,m-1,m,$$

(10)

which is graphically presented in Figure 3—with $W_i$ defined in (3) and (4).

Based on the work of [5,9] provided insight for detecting the change in the parameter structure if the underlying process is multivariate normal.

1.3. Outline of Paper

In Section 2, the bivariate joint probability density function which emanates from (7) is derived, accompanied by an exploratory shape analysis. This is followed by the derivation of the marginal probability density functions, the product moment as well as the maximum order statistic of the distribution. In Section 3 the performance of the model that this article proposes is compared to the Q chart studied by [5], this will be conducted through a simulation study. Tables of simulated values for the 95th percentiles of the maximum order statistics of the sets of random variables in (6) and (7) are provided in Section 3, for varying parameter choices, to enable practical application of the proposed model. The values in these tables are corroborated through numerical integration of the derived expressions for the maximum order statistics. The Appendix A and Appendix B contains proofs of the obtained results as well as the positioning of this distribution among other often considered bivariate beta models.

2. Proposed Model

In this section, the joint probability density function of the random variables $U_1$ and $U_2$ (see (7) when $r=2$) is derived, followed by a shape analysis, the derivation of the marginal probability density functions, the product moment as well as the maximum order statistic. A brief review of this new candidate with respect to its partners is provided in the Appendix A and Appendix B—which provides additional insight for modelling as well as expressions with closed form.

2.1. Bivariate Probability Density Function

Theorem 1.

Let $W_i$ be independent gamma random variables with parameters ${\alpha }_i=\frac{n_i-1}{2}>0,$${\beta }_i=2\lambda >0$ for $i=0,1,2.$ Let $U_1=\frac{W_1+W_2}{W_0}$ and $U_2=\frac{W_2}{W_0+W_1}$ then the joint probability density function of $U_1$ and $U_2$ is given by

$$\begin{array}{ccc}f\left(u_1,u_2\right)\hfill & =\hfill & \frac{\left({\beta }_0^{{\alpha }_1+{\alpha }_2}{\beta }_1^{{\alpha }_0+{\alpha }_2}{\beta }_2^{{\alpha }_0+{\alpha }_1}\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left(u_1-u_2\right)}^{{\alpha }_{\mathbf{1}}-1}{u}_2^{{\alpha }_2-1}{\left(1+u_1\right)}^{{\alpha }_2}{\left(1+u_2\right)}^{{\alpha }_0}\hfill \\ \hfill & \times \hfill & {\left({\beta }_1{\beta }_2\left(1+u_2\right)+{\beta }_0{\beta }_2\left(u_1-u_2\right)+{\beta }_0{\beta }_1{u}_2\left(1+u_1\right)\right)}^{-{\alpha }_0-{\alpha }_1-{\alpha }_2}\hfill \end{array}$$

(11)

where $u_1>u_2>0$.

2.2. Shape Analysis

In this section a shape analysis is conducted for the joint probability density function (11). A standard set of parameters has been chosen as a baseline. The parameters are chosen to be ${\alpha }_0={\alpha }_1={\alpha }_2=5$ and ${\beta }_0={\beta }_1={\beta }_2=2$, in other words, a process where all three samples consist of $5\times 2+1=11$ observations, and where no shift has occurred in the process variance. Some of the parameters will then be varied from this baseline in order to investigate the effect that a change in the specific parameters will have on the general shape of the joint probability density function. Note that the change in some parameters will be large - so large that they lose practical realism; this is conducted to emphasise and investigate the general change in the shape, and is not meant to be an indication of the practical applications of the joint probability density function. The functions will only be plotted for $u_1\in [0,3]$, $u_2\in [0,3]$.

In Figure 4 it can be seen what effect increasing the sample sizes (while keeping them equal) has on the joint probability density function. It is seen that increasing the sample sizes also increases the height of the peak of the probability density function. Larger sample sizes also shrink the length and width of the “tails” of the joint probability density function. In essence, the larger the sample sizes ($n_i=11,n_i=51,n_i=101$), the smaller the domain on which the function has non-trivial values.

Figure 5 below demonstrates that a sustained increase in the process variance, irrespective of size, minimally affects the general shape and location of the joint probability density function, but does affect the height of the probability density function’s peak. In the below example, the shift in the process variance occurs at time 1, and as one would hope and expect, the joint probability density function relies heavily on the value of the statistic at time 1, $U_1$. A similar effect, where the joint probability density function relies heavily on the value of $U_2$ is seen when the shift occurs at time 2.

2.3. Marginal Probability Density Functions

Theorem 2.

Assume that $\left(U_1,U_2\right)$ has the joint probability density function in (11), then the marginal probability density function of $U_1$ is given by

$$\begin{array}{ccc}\hfill & \hfill & f\left(u_1\right)\hfill \\ \hfill & \hfill & =\frac{\left({\beta }_0^{{\alpha }_1+{\alpha }_2}{\beta }_1^{{\alpha }_0+{\alpha }_2}{\beta }_2^{{\alpha }_0+{\alpha }_1}\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left({\beta }_1{\beta }_2+{\beta }_0{\beta }_2{u}_1\right)}^{-{\alpha }_0-{\alpha }_1-{\alpha }_2}\hfill \\ \hfill & \hfill & \times {(1+u_1)}^{{\alpha }_2}\sum _{k=0}^{\infty }\left[\left(\frac{{\alpha }_0!}{k!\left({\alpha }_0-k\right)!}\right){\left(u_1\right)}^{k+{\alpha }_1+{\alpha }_2-1}\frac{\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left(k+{\alpha }_2\right)}{\mathsf{\Gamma }\left(k+{\alpha }_1+{\alpha }_2\right)}\right.\hfill \\ \hfill & \hfill & \times \left.{}_2{F}_1\left({\alpha }_0+{\alpha }_1+{\alpha }_2,k+{\alpha }_2;k+{\alpha }_1+{\alpha }_2;-\frac{u_1\left({\beta }_1{\beta }_2-{\beta }_0{\beta }_2+{\beta }_0{\beta }_1+{\beta }_0{\beta }_1{u}_1\right)}{{\beta }_1{\beta }_2+{\beta }_0{\beta }_2{u}_1}\right)\right]\hfill \end{array}$$

(12)

where $u_1>0$, ${\alpha }_0,{\alpha }_1,{\alpha }_2\in \mathbb{Z}$, ${\beta }_1\left({\beta }_2+{\beta }_0\left(1+u_1\right)\right)>{\beta }_0{\beta }_2$, $\left|-\frac{u_1\left({\beta }_1{\beta }_2-{\beta }_0{\beta }_2+{\beta }_0{\beta }_1+{\beta }_0{\beta }_1{u}_1\right)}{{\beta }_1{\beta }_2+{\beta }_0{\beta }_2{u}_1}\right|<1,$ and the marginal probability density function of $U_2$ is given by

$$\begin{array}{ccc}\hfill & \hfill & f\left(u_2\right)\hfill \\ \hfill & \hfill & =\frac{\left({\beta }_0^{{\alpha }_1+{\alpha }_2}{\beta }_1^{{\alpha }_0+{\alpha }_2}{\beta }_2^{{\alpha }_0+{\alpha }_1}\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left(1+u_2\right)}^{{\alpha }_0}\sum _{k=0}^{\infty }\left[\frac{{\alpha }_2!}{k!\left({\alpha }_2-k\right)!}\right.\hfill \\ \hfill & \hfill & \times \sum _{l=0}^{\infty }\left[{\left(-1\right)}^l\left(\genfrac{}{}{0pt}{}{{\alpha }_0+{\alpha }_1+{\alpha }_2+l-1}{l}\right){\left({\beta }_0{\beta }_2+{\beta }_0{\beta }_1{u}_2\right)}^{-{\alpha }_0-{\alpha }_1-{\alpha }_2-l}\right.\hfill \\ \hfill & \hfill & \times \left.\left.{\left({\beta }_1{\beta }_2+{\beta }_1{\beta }_2{u}_2-{\beta }_0{\beta }_2{u}_2+{\beta }_0{\beta }_1{u}_2\right)}^l{u}_2^{k-{\alpha }_0-l-1}\frac{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_2+l-k\right)}{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2+l-k\right)}\right]\right]\hfill \end{array}$$

(13)

where $u_2>0$, ${\alpha }_0,{\alpha }_1,{\alpha }_2\in \mathbb{Z}$, ${\alpha }_i,{\beta }_i>0$ for $i=0,1,2$ and ${}_2{F}_1\left(,\right)$ is the Gauss hypergeometric function ([6] p. 1005). Note that in (12) if ${\alpha }_0\in \left\{1,2,\cdots \right\}$ the sum changes from ${\sum }_{k=0}^{\infty }$ to ${\sum }_{k=0}^{{\alpha }_0}$ (See [6] p. 25, Equations 1.110 and 1.111). Similarly if ${\alpha }_2\in \left\{1,2,\cdots \right\}$ the sum changes from ${\sum }_{k=0}^{\infty }$ to ${\sum }_{k=0}^{{\alpha }_2}$ in (13).

If ${\beta }_0={\beta }_1={\beta }_2$ it follows that:

$$\begin{array}{ccc}\hfill f\left(u_1\right)& =& \frac{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left(1+u_1\right)}^{-{\alpha }_0-{\alpha }_1}\sum _{k=0}^{\infty }\left[\left(\frac{{\alpha }_0!}{k!\left({\alpha }_0-k\right)!}\right){\left(u_1\right)}^{k+{\alpha }_1+{\alpha }_2-1}\right.\hfill \\ & \times & \frac{\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left(k+{\alpha }_2\right)}{\mathsf{\Gamma }\left(k+{\alpha }_1+{\alpha }_2\right)}{}_2{F}_1\left({\alpha }_0+{\alpha }_1+{\alpha }_2,k+{\alpha }_2;k+{\alpha }_1+{\alpha }_2;-u_1\right)\hfill \end{array}$$

where $0<u_1<1$, and

$$\begin{array}{ccc}\hfill f\left(u_2\right)& =& \frac{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_2\right)}{\left(1+u_2\right)}^{{\alpha }_0}\sum _{k=0}^{\infty }\left[\frac{{\alpha }_2!}{k!\left({\alpha }_2-k\right)!}\sum _{j=0}^{\infty }\left[{\left(-1\right)}^j\left(\begin{array}{c}{\alpha }_0+{\alpha }_1+{\alpha }_2+j-1\\ j\end{array}\right)\right.\right.\hfill \\ & \times & \left.\left.{\left(1+u_2\right)}^{-{\alpha }_0-{\alpha }_1-{\alpha }_2}{u}_2^{k-{\alpha }_0-j-1}\frac{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_2+j-k\right)}{\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1+{\alpha }_2+j-k\right)}\right]\right]\hfill \end{array}$$

where $u_2>0$.

2.4. Product Moment and Order Statistics

Theorem 3.

Assume that $\left(U_1,U_2\right)$ has the joint probability density function (11), then the product moment of $U_1$ and $U_2$ is given by

$$\begin{array}{ccc}\hfill & \hfill & E\left(U_1^r{U}_2^s\right)\hfill \\ \hfill & \hfill & =\sum _{p=0}^r\left(rp\right)\frac{\left({\beta }_0^{{\alpha }_1-p-s}{\beta }_1^{-{\alpha }_1}{\beta }_2^{p+s}\right)\mathsf{\Gamma }\left({\alpha }_2+p+s\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1-p-s\right)\mathsf{\Gamma }\left({\alpha }_1+r-p\right)\mathsf{\Gamma }\left({\alpha }_0-r\right)}{\mathsf{\Gamma }\left({\alpha }_0\right)\mathsf{\Gamma }\left({\alpha }_1\right)\mathsf{\Gamma }\left({\alpha }_2\right)\mathsf{\Gamma }\left({\alpha }_0+{\alpha }_1-p\right)}\hfill \\ \hfill & \hfill & \times {}_2{F}_1\left({\alpha }_0+{\alpha }_1-p-s,{\alpha }_1+r-p;{\alpha }_0+{\alpha }_1-p;1-\frac{{\beta }_0{\beta }_2}{{\beta }_1{\beta }_2}\right)\hfill \end{array}$$

(14)

where $r,s\in \left\{0,1,2,\cdots \right\},$${\alpha }_0+{\alpha }_1>r+s$, ${\alpha }_0>{\alpha }_1$, ${\alpha }_0+{\alpha }_1>p$, and $\left|1-\frac{{\beta }_0{\beta }_2}{{\beta }_1{\beta }_2}\right|<1$.

The maximum order statistics are of importance in detecting whether a shift in the process variance does indeed occur, as was discussed in Section 1 and will be demonstrated in Section 3. Although a closed form expression for the maximum order statistic is not tractable in a closed form, an expression is provided that can be implemented to calculate numerical values.

Assume that $\left(U_1,U_2\right)$ has the joint probability density function (11), then the maximum order statistic of $U_1$ and $U_2$ can be determined either by

$$\begin{array}{ccc}P\left(max\left(U_1,U_2\right)<z\right)\hfill & =\hfill & P\left(U_1<z\right)\hfill \\ \hfill & =\hfill & {\int }_0^{z}f\left(u_1\right)d{u}_1\hfill \end{array}$$

(15)

when ${\alpha }_0,{\alpha }_1,{\alpha }_2\in \Re$ (since $u_1>u_2$).

3. Comparison Study and Discussion

Deriving the order statistics of the series of statistics in (6) and (7) is a complex task since they are neither independent nor identically distributed, (see [10]). Hence values for the 95th percentiles of the maximum order statistics are simulated for varying numbers of samples (m equal to 1, 4, 9, 14, 19, 24, 29, 49, 99 and 499) and sample sizes (n equal to 2, 5, 10, 15, 20, 25, 30, 50, 100 and 500) in Section 3.1.

In Section 3.1, properties of (7) are also studied when no shift in the variance occurs, this is imperative since the control limits of a control chart are constructed under the null hypothesis that no shift has occurred, or alternatively given that the process is in control. The in control properties that are investigated for both the Q chart and the new proposed model (6) are where the maximum order statistic is most likely to occur when no shift has occurred in the process variance. Practically, this is a very important question since, if the maximum order statistic consistently occurs at roughly the same place in the sequence of samples, it implies that the distribution should be treated with added suspicion if it indicates that a shift in the process occurs at this location.

In Section 3.2, the out of control performance of the newly proposed distribution is compared with the Q chart model investigated by [5]. The probability that the respective control charts will detect a shift in the process variance is investigated for differing sizes of shifts in the process variance. The comparison is made for varying numbers of samples as well as sample sizes.

3.1. Comparison When the Process Is in Control

The location of the maximum order statistic when the process is in control is investigated using graphs. For brevity’s sake only one possible combination ($m=29,n=15$) of the number of samples and sample sizes mentioned above is included in this paper. All of them, however, lead to the same general conclusion.

As can be seen from Figure 6 the [5] distribution’s maximum order statistic occurs most often at the first statistic, with the probability of the maximum occurring at subsequent statistics steadily decreasing. This implies that the Q chart becomes more stable as the process progresses, which makes practical sense since each subsequent statistic includes more of the sample data. The newly proposed distribution’s maximum order statistic occurs most often at the first statistic, and second-most often at the last statistic. This due to the way in which the statistics of the distribution are constructed (see (6)). This implies that while our proposed model may detect shifts at the ends of the samples, signals received at these locations should be treated with a bit of skepticism.

In

Table 1. 95th percentiles of the maximum order statistics of $U_1^{*},U_2^{*},\cdots ,U_m^{*}$.

mn	2	5	10	15	20	25	30	50	100	500
1	214.286	6.652	3.230	2.512	2.189	1.998	1.872	1.615	1.399	1.160
4	242.046	6.229	3.047	2.385	2.093	1.919	1.805	1.569	1.371	1.150
9	254.592	5.992	2.908	2.301	2.023	1.863	1.755	1.534	1.349	1.142
14	251.747	5.874	2.871	2.266	1.997	1.842	1.739	1.522	1.341	1.139
19	255.413	5.836	2.851	2.254	1.987	1.831	1.730	1.517	1.337	1.137
24	257.880	5.813	2.837	2.243	1.980	1.827	1.723	1.514	1.335	1.137
29	259.105	5.815	2.824	2.238	1.977	1.822	1.720	1.512	1.334	1.136
49	259.963	5.799	2.818	2.228	1.968	1.815	1.715	1.508	1.331	1.135
99	258.810	5.797	2.803	2.222	1.959	1.811	1.710	1.504	1.329	1.134
499	263.017	5.772	2.801	2.215	1.957	1.808	1.707	1.501	1.328	1.133

the 95th percentiles of the maximum order statistics of the newly proposed distributions are simulated (to the third decimal) using Monte Carlo simulation. In other words z in the equation $P\left(max\left(U_1^{*},U_2^{*},\cdots ,U_m^{*}\right)<z\right)=0.95.$ Similarly the 95th percentiles of the maximum order statistics of (6) are simulated in

Table 2. 95th percentiles of the maximum order statistics of $U_1,U_2,\cdots ,U_m$.

mn	2	5	10	15	20	25	30	50	100	500
1	399.762	12.091	5.933	4.638	4.072	3.735	3.513	3.062	2.688	2.275
4	1151.103	29.015	14.061	11.045	9.708	8.958	8.434	7.411	6.553	5.618
9	2413.470	57.222	27.649	21.699	19.101	17.623	16.615	14.615	12.984	11.184
14	3691.861	85.496	41.144	32.376	28.469	26.292	24.837	21.864	19.419	16.754
19	4944.966	112.744	54.656	43.071	37.870	34.913	33.020	29.081	25.816	22.318
24	6187.147	141.315	68.325	53.653	47.313	43.549	41.182	36.293	32.248	27.887
29	7528.303	170.174	81.753	64.408	56.621	52.313	49.397	43.531	38.678	33.457
49	12,560.904	281.708	135.824	107.163	94.200	86.900	82.168	72.395	64.349	55.728
99	25,774.144	565.834	270.688	213.065	187.991	173.740	164.160	144.638	128.625	111.365
499	126,434.246	2815.186	1353.458	1064.841	939.259	868.035	819.248	722.038	642.752	556.808

.

Note that:

• Since the critical values of the distribution are derived under the assumption of the null hypothesis, the equations for the maximum order statistic can be simplified to be constructed out of chi-square random variables instead of the more complex gamma case.
• When the sample sizes of all the samples are equal $\left(n_i=n\mathrm{for}i=1,\cdots ,m\right)$ the removal of the constant terms in (7) is superfluous since the series of statistics in (6) and (7) are merely scalar multiples of each other. This can be seen in Table 1 and Table 2 where $P\left(max\left(U_1^{*},U_2^{*},\cdots ,U_m^{*}\right)<z\right)\approx P\left(max\left(U_1,U_2,\cdots ,U_m\right)<mz\right).$
• Instead of using simulation, the maximum order statistics in Table 2 could have been calculated using (15).

  Table 3. Simulated and theoretical 95th percentiles of the maximum order statistics of $U_1,U_2$.

  Methodn	2	5	10	15	20	25	30	50	100	500
  Simulated values	399.762	12.091	5.933	4.638	4.072	3.735	3.513	3.062	2.688	2.275
  Theoretical values	399.000	12.083	5.920	4.640	4.068	3.735	3.515	3.064	2.688	2.276

  demonstrates their equivalence for the bivariate case by numerically integrating (15).

3.2. Comparisons When the Process Is out of Control

In this section, the proposed model’s potential to detect shifts in compared with that of the Q chart form investigated by [5]. The probability of signaling that a shift in the process variance has occurred depends on a few variables. In this paper, these variables are: the number of samples, the sample size of the samples, where in the process the shift in the process variance occurs, and the size of the shift. The figures in this section take all of these parameters into account.

In each of the following figures, the probability of signaling a shift in the variance is displayed as a function of the size of the shift, where the shift size $\lambda$, ranges from $\lambda =1$ (no shift) to $\lambda =5$ (a 500% increase in the process variance). Many different combinations of the number of samples, the sample sizes, as well as the locations of the shift were tested; however, only the graphs that illustrate key findings are included in this paper. The chosen parameters that were simulated are: number of samples $\left(m\right)$ equal to 10, 20 and 30, the sample sizes $\left(n\right)$ equal to 2, 5, 10 and 20, and the location of the shift in the process variance $\left(\kappa \right)$ occurring at (roughly, due to integer sample numbers) 25%, 50% and 75% of the way through the samples.

By simulating graphs such as those in Figure 7, Figure 8 and Figure 9, certain conclusions can be reached about the properties and efficacy of the two competing models ((5) and (9)):

• When $n=2$, irrespective of the number of samples, the newly proposed model outperforms the Q chart. There are some caveats however that should be noted:
  • In the above graphs, each plotted point was simulated 500,000 times during the Monte Carlo process. For all the $n=2$ graphs, the points vary erratically between each 0.05 increases in the shift size, and thus even for a large number of simulations the process cannot be described as “stable”.
  • The Q chart seems to be completely incapable of detecting an increase in the process variance when $n=2$, with the probability of detecting a shift remaining at roughly 5%, irrespective of the size of the shift.
  • While the new model’s probability of detecting a shift does increase as the size of the shift increases, it remains relatively low, at roughly 7% to 10%, just marginally higher than the 5% chance when the process is actually IC. This implies that while it might be theoretically possible to implement the new model for samples sizes of 2, it would likely not be a practically useful technique.
  • The new model’s probability of detecting a shift does not increase as the number of samples increases, as would be expected (and as is the case for the other choices of n).
• From these points above, it can be concluded that using a sample size of 2 does not lead to an effective control chart.
• For a small numbers of samples $\left(m=9\right)$, the newly proposed model outperforms the Q chart for all simulated sample sizes as well as locations of shifts (for all shift sizes).
• When there are 20 samples $\left(m=19\right)$, the newly proposed model outperforms the Q chart in nearly all situations. The Q chart does have a higher probability of detecting a shift in the process variance only when the sample sizes are small $\left(n=5\right)$, and the shift occurs relatively late in the process $\left(\kappa =15\right)$, for shifts in the process variance between $\lambda =3$ and $\lambda =4.75$. Since a 300% to 475% increase in the process variance is unlikely to occur in practice, the newly proposed model would likely be more effective for $m=19$.
• For $m=29$, sweeping statements about the performances of the two methods are more difficult to make since the plotted percentage lines cross often. However it can be said that:
  • For small sample sizes $\left(n=5\right)$, the proposed model outperforms the Q chart for small shifts in the process variance, whereas the Q chart performs better for larger shifts.
  • The Q chart performs at its best when the shift in the process variance occurs late in the series of samples.
  • For larger sample sizes $\left(n=20\right)$ the proposed model outperforms the Q chart when the shift in the process variance occurs early, but when the shift occurs roughly half way through the series of samples, or further, the performance of the two methods are very similar.

4. Concluding Remarks

In this paper a new bivariate beta type distribution is proposed that can be utilised to detect a shift in a process’s variance when the underlying process follows a normal distribution. The proposed model compares favourably with the Q chart model studied by [5] in most situations; especially when the number of samples is small, and when the process variance experiences a change early on in the series of samples. Future work can include focus on (i) when a shift occurs within a sample, (ii) expanding underlying distributional assumptions to that of the class of scale mixtures to consider departures from normality (see [11]), and (iii) the multivariate setup, of which we propose the probability density function in the following theorem and a proof outlined in the Appendix A and Appendix B.

Theorem 4.

Let $W_i$ be independent gamma random variables with parameters $\left({\alpha }_i>0,{\beta }_i>0\right)$ for $i=0,1,2,\cdots ,m.$ Let $U_r=\frac{{\sum }_{i=r}^m{W}_i}{{\sum }_{i=0}^{r-1}{W}_i},r=1,2,\cdots ,m-1,m,$ then the joint probability density function of $U_1,U_2,\cdots ,U_m$ is given by

$$\begin{array}{ccc}\hfill & \hfill & f\left(u_1,u_2,\cdots ,u_m\right)\hfill \\ & \hfill & =\frac{{\prod }_{i=1}^{m-1}\left[{\left(u_i-u_{i+1}\right)}^{{\alpha }_i-1}\right]{\left(u_m\right)}^{{\alpha }_m-1}\mathsf{\Gamma }\left({\sum }_{i=0}^m{\alpha }_i\right)}{{\prod }_{i=0}^m\left[{\beta }_i^{{\alpha }_i}\mathsf{\Gamma }\left({\alpha }_i\right)\right]}\hfill \\ \hfill & \hfill & \times {\left(1+u_1\right)}^{{\sum }_{i=2}^m{\alpha }_i}\prod _{i=2}^m\left[{\left(1+u_i\right)}^{-{\alpha }_{i-1}-{\alpha }_i}\right]\hfill \\ \hfill & \hfill & \times {\left(\frac{1}{{\beta }_0}+\frac{\left(u_1-u_2\right)}{{\beta }_1\left(1+u_2\right)}+\sum _{i=2}^{m-1}\left[\frac{\left(1+u_1\right)\left(u_i-u_{i+1}\right)}{{\beta }_j\left(1+u_i\right)\left(1+u_{i+1}\right)}\right]+\frac{\left(1+u_1\right)u_m}{{\beta }_m\left(1+u_m\right)}\right)}^{-{\sum }_{i=0}^m{\alpha }_i}\hfill \end{array}$$

(16)

where $u_1>u_2>\cdots >u_m>0.$