On monitoring of linear profiles using Bayesian methods

Computers & Industrial Engineering 94 (2016) 245–268

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

On monitoring of linear profiles using Bayesian methods Tahir Abbas a,b,⇑, Zhengming Qian a, Shabbir Ahmad c, Muhammad Riaz d a

Department of Statistics, School of Economics, Xiamen University, Xiamen 361005, China Department of Statistics, Government College University Lahore, Lahore 54000, Pakistan c Department of Mathematics, COMSATS Institute of Information Technology, Wah Cantt 47040, Pakistan d Department of Mathematics and Statistics, King Fahad University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia b

a r t i c l e

i n f o

Article history: Received 21 June 2015 Received in revised form 30 January 2016 Accepted 5 February 2016 Available online 12 February 2016 Keywords: Linear profile Prior distribution Posterior distribution Run length properties Performance Comparison Index (PCI)

a b s t r a c t In many industrial applications, the quality of a process or product is distinctly illustrated by a linear profile. Mostly, a linear profile is shaped or modeled through the usage of linear regression model. Classical and Bayesian set-ups are two possible ways of defining the design structure of a linear profile. This study proposes a fresh Bayesian approach by using different priors for monitoring the linear profiles of processes. In this research we constructed three novel univariate Bayesian EWMA control charts for the Y-intercepts, the slopes and the errors variance under phase II methods by using non-conjugate and conjugate priors. These control charts are used to monitor the Y-intercepts, the slope coefficients and increase in process standard deviations, respectively. This study confirmed that the Bayesian methods distinguish sustainable shifts in the process parameters superior than the competing methods. Moreover, the Bayesian control charting structures with conjugate priors provide better performance for monitoring the Y-intercepts and slopes than the one of non-conjugate priors, while both priors perform almost equivalently in case of errors variance. The individual and overall performance of control charts are evaluated by using (i.e., ARL, SDRL, and MDRL) and (i.e., EQL, RARL, and PCI), respectively. The practical example is considered as an illustration to justify the supremacy of proposed approach and recommendations are given for future perspective. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Since the inception of 1924, the significance of Statistical Process Control (SPC) has been repeatedly asserted by its extensive use in manufacturing industry. Considering the applicability of SPC, the basic premise states that the quality of product with a single feature is well represented by the distribution of univariate quality characteristic and, in the presence of multiple quality characteristics then multivariate distribution is considered as the most efficient approach. Two different approaches can be considered while constructing control charts: first, distribution function can be used to determine process position for single or numerous quality features. Second, a profile approach can be used to monitor a product or process. In the latter, quality of product or process is well represented by relationship among dependent variable and one or several regressors. Profiles occur in many areas of manufacturing industries, i.e., Kang and Albin (2000), provided the case when the characterization of ⇑ Corresponding author at: Department of Statistics, School of Economics, Xiamen University, Xiamen 361005, China. E-mail address: [email protected] (T. Abbas). http://dx.doi.org/10.1016/j.cie.2016.02.007 0360-8352/Ó 2016 Elsevier Ltd. All rights reserved.

amount of dissolved aspartame (an artificial sweetener) per liter of water at different levels of temperature is under study. In literature many practical situations can be seen when the process is positively managed by using profile i.e., in simple as well as complex ways. Likewise, authors highlighted the profile monitoring in both situations (e.g., Amiri, Eyvazian, Zou, & Noorossana, 2012; Keramatpoura, Niakib, & Amiric, 2014; Zhang, He, Zhang, & Woodall, 2013). Usually, control charting in SPC is carried out using either phase I or phase II applications, however both exhibit distinct goals and manifestations. These two methods of control charting structures and their properties are discussed in detail by Montgomery (2009, pp. 476–478). Kim, Mahmoud, and Woodall (2003), Zhu and Lin (2010), Zhang et al. (2013), Keramatpoura et al. (2014). Kang and Albin (2000), presented an efficient and extensive use of phase I and phase II methods for control charting structures. On the other hand Brill (2001) went for more general perspective than linear one. Moreover, Kim et al. (2003), used coded values of regressor and constructed three univariate EWMA (Exponentially Weighted Moving Average) control charts using phase II methods for profile monitoring. They used both logarithms and natural scale of process standard deviations for profile monitoring.

246

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Using simulation, they tested that deviation from mean method confer better performance, compared with that of Kang and Albin (2000) method. Moreover, a range of control charting structures are available in scientific literature which are designed under classical and Bayesian setups. The technological and scientific advancement in industrial field made the machines very costly and upgrading also requires high capital. This makes essential to use improved methods to monitor the statistical variations of process parameters. An efficient SPC techniques can timely identifies these variations. The scope of efficient SPC tools is not limited to production process but also cover the areas such as manufacturing and refining, food and drink, pharmaceuticals etc. These efficient SPC techniques interact with more information about the current and past behaviors of the processes, and produce better diagnostics. Bayesian methods are considered an efficient choice to address such diverse and complex problems. Bayesian control charting structures for the distribution functions of particular quality characteristics are widely available in the literature and can be used to enhance the performance of process parameter (e.g., Ali & Riaz, 2014; Matthias, Tana, & Jianjun, 2012; Marcellus, 2007; Roberts, 1966). In previous studies the authors focused the distribution function approach to determine process position for simple and multivariate quality characteristics and tried to interpret the performance of mean shifts for the posterior distributions of processes or products. From the Bayesian perspective, it could be argued that since probability is artistically elegant to communicate uncertainty and hence, considering vagueness of the parameters of interest in a probabilistic manner seems more logical. From a parametrical vagueness standpoint, it is more appealing to aggregate most credible values of the parameter rather than looking for one value of parameter (based on sampled data). Also, Bayesian philosophy encourages examination of probability distribution for parameter under study, in addition to observing distribution of sampled data. A plethora of scholarly literature testifies the wide application of Bayesian approach (e.g., Durham & Geweke, 2014; Florens & Simoni, 2012; Zellner, Ando, Basturk, Hoogerheide, & Van Dijk, 2014). Thus, in the Bayesian approach, a combine model of response and parameters is specified by using the Bayes theorem, named after Thomas Bayes, a famous English statistician:

PðcjDÞ ¼

LðDjcÞPðcÞ ; PðDÞ

ð1Þ

where, ‘‘D” is the sample data and ‘‘ c ” refers to the parameter under study. Let LðDjcÞ be the likelihood function or joint probability distribution of sample data, PðcÞ prior distribution of the parameter, PðDÞ normalizing quantity (most strict rule of Bayesian Approach) and simply scale the posterior distribution. An observable benefit of the Bayesian approach is the ability to pack into a single corridor for finding a good estimator. The posterior distribution is the combination of joint probability distribution (the likelihood function) of sample data and prior distribution or adding prior knowledge into likelihood function. Referring to previous formula:

pðcjDÞ / LðDjcÞ  pðcÞ:

ð2Þ

The assortment of priors has proved to be a daunting experience for researchers, based on unyielding priors information. The expert’s opinions may vary from one to another. Priors can be informative and non-informative, depending upon the amount of information contains. During the selection of informative priors, conjugate priors are attractive due to mathematical simplification of posterior distributions; however non-conjugate priors are also valuable choice in many applications. Choosing conjugate priors for regression models endows convenient close form of the posterior distributions (e.g., Gelman, Carlin, Stern, & Rubin, 2003; Paciorek, 2006).

This study further enhances the SPC methodology by introducing new Bayesian methods to construct the Bayesian EWMA control charts for the process parameters defined by linear profiles. We consider different non-informative (uniform priors) and informative priors (non-conjugate and conjugate priors) and then constructing the posterior distributions of the process parameters. The idea is based on the premise that a deep and watchful analysis of a production process parameters can lead to improve the SPC system. The production process parameters can be analyzed by using the information such as the time series of measurement, the maintenance and repairing history and judgments and observations of engineers about the process. These informations about the process parameters are well represented in the form of probability distributions. These probability distributions mostly called prior distributions weights the production process to search for assignable causes when the process is in control or out of control. We tried to refine the efficiency of the EWMA control charts for profile monitoring by accumulating these prior informations. The flow chart in Fig. 1 provides the better understanding of computational algorithm adopted for these situations. Our simulation study reflects that there are production process (defined by a linear profile) where Bayesian approach has better expected overall performance than that of classical approach in terms of different individual and overall performance measures. The outlines of this research article are as: Section 2 presents control charting structures in phase II methods. Section 3 describes comparative performance measures analysis. Section 4 shows an example of EWMA control charts for profile monitoring under classical and Bayesian setups and finally, Section 5 presents conclusion and recommendations. 2. Linear profile and control charting structures Let us consider set of observations ðxi ; yij Þ, where i ¼ 1; 2; . . . ; n and j ¼ 1; 2; . . . ; k. For in-control situation of statistical process, simple linear regression model is assumed as:

Y ij ¼ A0 þ A1 X i þ eij :

ð3Þ

Let’s take a straightforward case by assuming fixed X-values and known in-control values of parameters. Kim et al. (2003) takes deviation form of X-values, results in independent least square estimates Ryan (1997, pp. 38–39), then customized version of Eq. (3) is:

Y ij ¼ B0 þ B1 X i þ eij ;

ð4Þ

 B1 ¼ A1 and where, B0 ¼ A0 þ A1 X;

X i

 ¼ X i  X.

2.1. Classical estimates of linear profile The least squares estimates of B0 and B1 in above model of Eq. (4) for the jth sample are:

j ; boj ¼ y

And b1j ¼

SxyðjÞ ; Sxx

ð5Þ

P P P j ¼ 1n ni¼1 yij ; SxyðjÞ ¼ ni¼1 ðxi  xÞyij , Sxx ¼ ni¼1 ðxi  xÞ2 . where y The ordinary least squares estimators are normally distributed with mean vector lc (subscript ‘‘c” stands for classical) and the variance–covariance matrix Rc :

!

lc ¼ ðB0 ; B1 ÞT and Rc ¼ 2

2

r20 r01 ; r10 r21

where r20 ¼ rn , r21 ¼ Srxx , and r01 ¼

Pn



i¼1

ðxi xÞr2

Sxx

ð6Þ ¼ 0: For every succes-

sion of random variables, a separate EWMA control chart can be constructed for the Y-intercept, the slope of variable X and variance of error terms to detect sustainable shifts under phase II methods.

247

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Generate sample of size “n” from standard normal probability distribution

Construct the Posterior distribution of linear profiles parameters

Noninformative

Informative Priors

Conjugate

Nonconjugate

Calculate the posterior estimates of the linear profiles parameters

Elicitating the hyper-parameters of prior distributions Determine the values of control limits coefficients of LI, LS, LE for Intercept, Slope and errors variance respectively, for pre-specified ARL0

For non-symmetrical distribution use non-symmetric control limits, (LIU, LIL), (LSU, LSL), (LEU, LEL)

For symmetrical distribution use symmetric control limits, LI, LS, LE

Generate the single observation from standard normal probability distribution and construct three univariate Bayesian EWMA statistics for profiles parameters Replicate the above process “m” times to for three univariate EWMA statistics simultaneously

Count the values that fall outside the control limits for intercept, slope and errors variance Repeat the above process 10,000 times and calculate the average run length then use these values to determine other performance Fig. 1. Procedural flow chart of Bayesian EWMA charts.

2.2. Posteriors estimates of linear profile

2.3. Posterior estimates using uniform priors

In light of the Bayesian approach, we consider the different prior (i.e., uniform prior, conjugate prior, non-conjugate prior) distributions to construct the posterior distributions for the three unknown(s) in linear regression model of Eq. (4). We assume that error terms are identically and independently normally distributed for the jth sample with mean zero and fixed variance.

Consider the independent flat priors (uniform priors) for the Y-intercept, the slope and the errors variance:

pðB0 ; B1 Þ / 1;

pðr2 Þ / r2 ;

ð7Þ

The uniform priors above are improper as their integration results are undefined. Even though the literature suggests that the priors are improper, as long as the posterior distributions are valid the

248

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

justifiable statistical inferences can be conducted. The likelihood functions for the jth sample of the sampling distribution of dependent variable Yj are specified as:

 Lðyj ; B0 ; B1 ; r2 Þ ¼

1

n=2

) n   1 X  2 :  exp  2 y  B0  B1 xi 2r i¼1 ij

ð8Þ

Now, following Eq. (2) the joint posterior distribution for Y-intercept and slope is given as:

(

n  2 1 X y  B0  B1 xi 2 2r i¼1 ij

)

ð9Þ

r ¼

spread out the square and compare the coefficients of B20 and B21 (which results in posterior variances), while matching the coefficients of B0 and B1 results in posterior means of posterior distributions for the Y-intercept and slope. The algebraic simplifications provide the similar posterior estimates and the joint posterior distributions for the Y-intercept and slope as given in Eqs. (5) and (6). This shows that profile monitoring with Bayesian approach using uniform priors and classical approach provide identical results. Constructing the posterior distributions for errors variance, again following Eq. (2) the algebraic simplification results in the posterior distribution for variance of errors as inverse gamma dis  ðn2Þr2 ; 2 0 . This again shows the tributions i.e., pðr2 jdataÞ  IG n2 2 similarity with classical approach. This also indicates that the Bayesian analysis with uniform priors will provide the identical results with that of classical approach. 2.4. Posterior estimates using non-conjugate priors Now, we consider the non-conjugate priors for the Y-intercept and slope in Eq. (8) above to construct the posterior distributions. We consider the Bramwell, Holdsworth, Pinton (BHP) distribution as prior knowledge of process parameters, introduced by Bramwell, Holdsworth, and Pinton (1998). This probability distribution is continuous, univariate, and unimodel with infinite range. This two parameters (e.g., location, scale) probability distribution is well known for monitoring rare fluctuations in systems or processes. The functional form is straightforward and given as:

np x  t x  to 1 exp  exp : jkj 2 k k 2

BHPðx; t; kÞ ¼ p

) n   1 X  2 y  B  B x 0 1 i 2r2 i¼1 ij      p B0  t0 B0  t0  exp  exp 2 k0 k0      p B1  t1 B1  t1  exp ;  exp 2 k1 k1

(

pðB0 ; B1 jdataÞ / exp

00

r2 k40 exp

  t0

  t0

nk40 exp

k0

k0

and

þ r2

r

2 1nj

¼

r2 k41 exp

  t1 k

 1 : Sxx k41 exp tk11 þ r2

ð14Þ

The distributions of Bayes estimates bonj and b1nj are normally distributed with mean lbn (subscript ‘‘bn” stands for Bayes estimates using non-conjugate priors) and variance–covariance matrix with zero covariance are Rbn . The informative (non-conjugate) prior for r2 selected as Levy distribution with one specified hyper-parameter ‘‘g0 ”. The Levy distribution is stable continuous probability distribution for nonnegative random variables. The Levy distribution is frequently used to represent process variation (e.g., Jonathan & Roger, 2014). The posterior distribution of errors variance can be constructed by using Eq. (2) above as:

pðr

2

(

n  2 1 X jdataÞ / r exp y  B0  B1 xi 2r2 i¼1 ij ng o   3 0 :  r2 2 exp 2r2



 n 2 2

)

ð15Þ

After some algebra of Eq. (15) above, the posterior distribution for errors variance obtain as the inverse gamma distribution i.e., r2  IGðanj ; bnj Þ. For the jth sample, the posterior parameters of the posterior inverse gamma distribution are (‘‘see Appendix A for derivation”).

anj ¼

nþ1 2

and bnj ¼ Pn

When

r20 ¼

i¼1

g0 2

ðyij B0 B1 xi Þ n2

þ

ðn  2Þr20 ; 2

ð16Þ

2

.

For the case of conjugate Bayesian analysis for the normal linear regression model the priors for the Y-intercept and slope are con    sidered as normal priors e.g., B0  N b0 ; j20 and B1  N b1 ; j21 respectively, where (b0 , b1 , j20 , j21 ) are the specified hyperparameters. The above priors are considered conjugate priors as the resultant posterior distributions are also normal. Now, following Eq. (2) the joint posterior distribution of the Y-intercept and slope is given as:

(

ð11Þ

where ðt0 ; t1 ; k0 ; k1 Þ are the specified hyper-parameters. The special case of Taylor series, the Maclaurin series has been used to simplify the Eq. (11) above. 000

f ðxÞx f ðxÞx2 f ðxÞx3 þ þ þ  1! 2! 3!

2 onj

ð13Þ

2.5. Posterior estimates using conjugate priors

ð10Þ

Then, following Eq. (2) the joint posterior distribution of Y-intercept and slope is given as:

0

     j k0 þ p2 r2 k30 exp tk0 þ t0  k20 r2 ny 0   bonj ¼ and nk40 exp tk00 þ r2       k1 SxyðjÞ þ p2 r2 k31 exp tk11 þ t1  k21 r2   : b1nj ¼ Sxx k41 exp tk11 þ r2 While their variances are:

 1;

For Eq. (9) above with uniform priors on the Y-intercept and slope,

f ðxÞ ¼ f ðxÞ þ

B1  Nðb1nj ; r21nj Þ. Designed for jth sample, the Bayes estimator of B0 and B1 are (‘‘see Appendix A for derivation”).



2pr2 (

pðB0 ; B1 jdataÞ / exp

and slope are normal for the jth sample e.g., B0  Nðb0nj ; r20nj Þ and

ð12Þ

Here, ‘‘x” is the real or complex number that is centered on zero. After some algebra, the posterior distributions of the Y-intercept

n  2 1 X pðB0 ; B1 jdataÞ / exp y  B0  B1 xi 2r2 i¼1 ij  1 2 ð B  b Þ  exp 0 0 2j20  1 2  exp : ð B  b Þ 1 1 2j21

)

ð17Þ

For the Eq. (10) above, we expand the square terms and compare the coefficients of the Y-intercept and slope. After simplification we acquire joint posterior distributions for the Y-intercept and

249

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

slope as normal e.g., B0  Nðb0nj ; r20nj Þ and B1  Nðb1nj ; r21nj Þ. Designed for jth sample, the Bayes estimator of B0 and B1 are:

bonj ¼

j j20 þ b0 r2 ny nj20 þ r2

and b1nj ¼

SxyðjÞ j21 þ b1 r2 : Sxx j21 þ r2

ð18Þ

While their variances are:

r2onj ¼

pðr2 jdataÞ /

j20 r2 j21 r2 2 and r ¼ : 1nj 2 nj0 þ r2 Sxx j21 þ r2

ð19Þ

The distributions of Bayes estimates bonj and b1nj are also normally distributed with mean lbc (subscript ‘‘bc” stands for Bayes estimates using conjugate priors) and variance–covariance matrix Rbc with zero covariance. The conjugate prior of inverse gamma distributions i.e., r2  IGða0 ; b0 Þ, (when a0 and b0 are the specified hyperparameters) for errors variance ensure that the posterior distribution for variance of errors is again inverse gamma distribution that can easily be derived. The posterior distribution for errors variance is given as:

(

n  2 1 X y  B0  B1 xi 2r2 i¼1 ij   a 1 1 b :  r2 0 exp 2 0



r2

n2

)

exp

r

ð20Þ

After some algebra of Eq. (20) above the posterior distribution for errors variance that is the inverse gamma distribution i.e., r2  IGðanj ; bnj Þ. For the jth sample, the posterior parameters of the posterior inverse gamma distribution are:

anj ¼ a0 þ

n 2

and bnj ¼ b0 þ

ðn  2Þr20 : 2

ð21Þ

2.6. EWMA control structures We propose three univariate Bayesian EWMA control charts with the objective of monitoring the Y-intercepts first, then slopes and at the end errors variance. It seems logical to construct three univariate EWMA control charts (when the course of action is characterized by profile), while the persistent

Table 1 ARL, SDRL and MDRL values for sensitivity analysis of NIG priors for shifts in Y-intercept at ARL0 = 200. dI

Sensitivity Analysis of Hyper-parameters (Location) EWMA-BCP-1 b0 ¼ 15, b1 ¼ 1:5

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

EWMA-BCP-2 b0 ¼ 20, b1 ¼ 2:5

EWMA-BCP-3 b0 ¼ 25, b1 ¼ 3:5

EWMA -BCP-4 b0 ¼ 30, b1 ¼ 4:5

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

199.95 54.66 14.49 6.90 4.28 2.96 2.22 1.80 1.51 1.32 1.17

200.79 51.54 11.32 4.45 2.44 1.50 1.06 0.81 0.64 0.50 0.38

135.00 39.00 12.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.22 44.78 12.94 6.35 3.99 2.87 2.17 1.77 1.49 1.31 1.17

197.97 40.72 9.38 3.91 2.18 1.44 1.02 0.79 0.62 0.50 0.39

137.00 32.50 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.27 40.45 12.39 6.26 3.98 2.88 2.20 1.79 1.53 1.30 1.18

195.19 36.16 8.94 3.65 2.07 1.40 1.03 0.78 0.63 0.49 0.40

137.00 29.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.98 39.08 12.32 6.35 4.08 2.96 2.27 1.85 1.57 1.36 1.20

194.61 34.54 8.66 3.54 2.06 1.41 1.02 0.80 0.63 0.52 0.41

138.00 29.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

Sensitivity Analysis of Hyper-parameters (Scale)

j20 ¼ 45, j21 ¼ 16 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

200.15 46.90 13.24 6.49 4.02 2.84 2.21 1.79 1.50 1.30 1.17

199.35 42.56 9.90 4.07 2.25 1.43 1.03 0.80 0.62 0.50 0.39

j20 ¼ 35, j21 ¼ 12 137.00 34.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.22 44.78 12.94 6.35 3.99 2.87 2.17 1.77 1.49 1.31 1.17

197.97 40.72 9.38 3.91 2.18 1.44 1.02 0.79 0.62 0.50 0.39

j20 ¼ 25, j21 ¼ 8 137.00 32.50 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.61 42.27 12.58 6.27 3.98 2.82 2.19 1.80 1.49 1.31 1.16

196.64 38.02 9.02 3.73 2.13 1.41 1.04 0.78 0.62 0.50 0.37

j20 ¼ 20, j21 ¼ 6 137.00 31.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

200.03 40.48 12.41 6.24 4.00 2.86 2.20 1.81 1.52 1.31 1.19

196.02 36.14 8.91 3.64 2.09 1.40 1.03 0.79 0.62 0.50 0.40

138.00 29.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

Sensitivity Analysis of Hyper-parameters (IG)

a0 ¼ 0:1, b0 ¼ 0:1 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

199.92 44.78 12.88 6.39 3.97 2.88 2.16 1.78 1.49 1.30 1.17

198.11 40.62 9.34 3.96 2.19 1.45 1.01 0.79 0.62 0.50 0.38

a0 ¼ 0:68, b0 ¼ 0:3 137.00 32.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.22 44.78 12.94 6.35 3.99 2.87 2.17 1.77 1.49 1.31 1.17

Note: NIG = Normal Inverse Gamma, IG = Inverse Gamma.

197.97 40.72 9.38 3.91 2.18 1.44 1.02 0.79 0.62 0.50 0.39

a0 ¼ 1:5, b0 ¼ 0:9 137.00 32.50 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.88 44.77 12.87 6.40 3.98 2.89 2.16 1.77 1.49 1.30 1.17

198.11 40.62 9.31 3.96 2.19 1.45 1.02 0.79 0.62 0.50 0.38

a0 ¼ 2:5, b0 ¼ 2:0 137.00 32.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

200.24 44.80 12.90 6.38 3.96 2.88 2.16 1.77 1.50 1.30 1.17

198.38 40.67 9.37 3.98 2.17 1.45 1.02 0.78 0.62 0.49 0.39

137.00 32.00 10.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

250

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kj20 r2  ; and LCLB ¼ B0  LI ð2  kÞ nj20 þ r2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kj20 r2  : UCLB ¼ B0 þ LI ð2  kÞ nj20 þ r2

revolutionization in process expected that at least one will signal. When each parameter corresponds to a control chart, it will leads to easier detection of possible process change. We integrate the noted control charts by the preceding rationale. First, Crowder and Hamilton (1992) stressed the use of logarithmic transformation of sample variances for monitoring increase in variance of process and then, the comparative use of natural scale for monitoring of linear profile. It is provided here with the Bayesian control charting structure under conjugate priors. The Bayesian control charting structure under non-conjugate priors and classical control charting structures can also be constructed on the similar lines. First, monitoring the Y-intercept estimates, consider the posterior estimates b0nj to compute Bayesian EWMA statistic:

The initial value chosen here as B1 . Then upper and lower control limits are:

EWMAI ðjÞ ¼ kb0nj þ ð1  kÞEWMAI ðj  1Þ

LCLB ¼ B1  LS

_ j ¼ 1; 2; . . .

ð22Þ

Here, kð0 < k 6 1Þ is the smoothing constant and B0 is chosen as preliminary value. The control chart will endow with out-ofcontrol signal whenever EWMAI ðjÞ < LCLB or EWMAI ðjÞ > UCLB , where the upper and lower control limits are:

ð23Þ

The subscript ‘‘I” and ‘‘B” in Eqs. (22) and (23) above stand for intercept and Bayesian respectively. Second, for monitoring the slope estimates, consider the posterior estimates b1nj to compute Bayesian EWMA statistic:

EWMAS ðjÞ ¼ kb1nj þ ð1  kÞEWMAS ðj  1Þ

_ j ¼ 1; 2; . . .

ð24Þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kj21 r2  ; and ð2  kÞ Sxx j21 þ r2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kj21 r2  : UCLB ¼ B1 þ LS ð2  kÞ Sxx j21 þ r2

ð25Þ

The subscript ‘‘S” in Eqs. (24) and (25) above stands for slope.

Table 2 ARL, SDRL and MDRL values for sensitivity analysis of NIG priors for shifts in slope at ARL0 = 200. dS

Sensitivity Analysis of Hyper-parameters (Location) EWMA-BCP-1 b0 ¼ 15, b1 ¼ 1:5

0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

EWMA-BCP-2 b0 ¼ 20, b1 ¼ 2:5

EWMA-BCP-3 b0 ¼ 25, b1 ¼ 3:5

EWMA-BCP-4 b0 ¼ 30, b1 ¼ 4:5

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

199.95 95.42 34.74 15.92 9.10 6.18 4.44 3.52 2.85 2.38 2.06

200.79 95.33 30.39 12.21 6.28 3.80 2.55 1.86 1.46 1.15 0.96

135.00 67.00 26.00 13.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.22 77.58 29.43 13.97 8.51 5.79 4.31 3.40 2.81 2.32 2.01

197.97 76.30 25.21 10.34 5.54 3.42 2.33 1.75 1.38 1.09 0.91

137.00 54.00 22.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.27 70.22 27.29 13.36 8.22 5.76 4.33 3.42 2.80 2.38 2.05

195.19 67.78 23.09 9.65 5.19 3.30 2.37 1.74 1.34 1.11 0.91

137.00 48.00 21.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.98 67.96 26.69 13.41 8.35 5.90 4.49 3.52 2.90 2.47 2.12

194.61 64.51 22.23 9.48 5.22 3.26 2.37 1.74 1.37 1.10 0.92

138.00 48.00 20.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

Sensitivity Analysis of Hyper-parameters (Scales)

j20 ¼ 45, j21 ¼ 16 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

200.15 81.24 30.65 14.50 8.52 5.89 4.36 3.36 2.78 2.34 2.04

199.35 80.44 26.42 10.66 5.64 3.47 2.39 1.76 1.36 1.11 0.95

j20 ¼ 35, j21 ¼ 12 137.00 57.00 23.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.22 77.58 29.43 13.97 8.51 5.79 4.31 3.40 2.81 2.32 2.01

197.97 76.30 25.21 10.34 5.54 3.42 2.33 1.75 1.38 1.09 0.91

j20 ¼ 25, j21 ¼ 8 137.00 54.00 22.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

200.03 71.07 27.49 13.38 8.25 5.83 4.29 3.42 2.82 2.38 2.06

196.02 68.80 23.10 9.76 5.21 3.39 2.32 1.76 1.35 1.09 0.92

j20 ¼ 20, j21 ¼ 6 138.00 49.00 21.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.55 67.81 26.42 13.45 8.29 5.84 4.46 3.50 2.90 2.44 2.13

195.09 64.45 21.95 9.48 5.18 3.25 2.36 1.74 1.37 1.10 0.92

138.00 48.00 20.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

Sensitivity Analysis of Hyper-parameters (IG)

a0 ¼ 0:1, b0 ¼ 0:1 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

199.92 77.73 29.50 14.08 8.46 5.79 4.33 3.37 2.81 2.32 2.03

198.11 76.25 25.25 10.36 5.56 3.40 2.33 1.74 1.38 1.09 0.91

a0 ¼ 0:68, b0 ¼ 0:3 137.00 54.00 22.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.22 77.58 29.43 13.97 8.51 5.79 4.31 3.40 2.81 2.32 2.01

197.97 76.30 25.21 10.34 5.54 3.42 2.33 1.75 1.38 1.09 0.91

a0 ¼ 1:5, b0 ¼ 0:9 137.00 54.00 22.00 11.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.88 77.67 29.48 14.07 8.47 5.79 4.32 3.36 2.81 2.32 2.02

198.11 76.24 25.24 10.37 5.56 3.40 2.33 1.73 1.38 1.10 0.91

a0 ¼ 2:5, b0 ¼ 2:0 137.00 54.00 22.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

200.24 77.67 29.50 14.11 8.47 5.78 4.33 3.38 2.80 2.31 2.03

198.38 76.17 25.25 10.36 5.55 3.39 2.33 1.75 1.39 1.09 0.92

137.00 54.00 22.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

251

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Third, the Bayesian EWMA control chart for errors variance is constructed. The posterior inverse gamma distribution of errors variance is used to derive the formula for variance of the logarithm values of mean square error for the jth sample (‘‘see Appendix B for derivation”).

Var½lnðMSEj Þ ¼

2 2 4 16 þ þ  : Q 1 Q 2 3Q 3 15Q 4

ð26Þ

With distinction to previous cases of the Y-intercept and slope, Bayesian EWMA control chart for errors variance is one-sided. The Bayesian EWMA test statistic for monitoring the increase in errors variance is:

   EWMAE ðjÞ ¼ max k lnðMSEj Þ þ ð1  kÞEWMAE ðj  1Þ; ln r20 _ j ¼ 1; 2; . . .

ð27Þ

  Here the initial value is EWMAE ð0Þ ¼ ln r20 . It is in following of Kim et al. (2003) that r20 ¼ 1, the in-control value of r2 that results EWMAE ð0Þ ¼ 0. When EWMAE > UCLB ,

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  k Var lnðMSEj Þ ; UCLB ¼ 0 þ LE ð2  kÞ

ð28Þ

The subscript ‘‘E” in Eqs. (27) and (28) above stands for errors under log scale.

Finally, in contrast to above Bayesian EWMA control chart for monitoring the variance of errors under log scale, it is demonstrated the Bayesian EWMA control chart for variance of errors under natural scale in order to evaluate the effectiveness and possible differences of natural scale as compare to log scale. Further, the following Bayesian EWMA statistic for errors variance was adopted based on natural scale:

 EWMAE ðjÞ ¼ max kðMSEj  1Þ þ ð1  kÞEWMAE ðj  1Þ; 0 _ j ¼ 1; 2; . . .

ð29Þ

The initial value under natural scale again equal zero i.e., EWMAE ð0Þ ¼ 0, as the in-control value of process variance is one, and then the upper control limit is defined as:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k VarðMSEj Þ; UCLB ¼ 0 þ LE ð2  kÞ

ð30Þ



The subscript ‘‘E⁄” in Eqs. (29) and (30) above stands for errors under natural scale. 2.7. Elicitation of hyper-parameters For the simulation study of previous construction, now incorporate elicitation of hyper-parameters values for the prior

Table 3 ARL, SDRL and MDRL values for sensitivity analysis of BHP levy priors for shifts in Y-intercept at ARL0 = 200. dI

Sensitivity Analysis of Hyper-parameters (Location) EWMA-BNCP-1 t0 ¼ 15, t1 ¼ 1:5

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

EWMA-BNCP-2

EWMA-BNCP-3

EWMA-BNCP-4

t0 ¼ 20, t1 ¼ 2:5

t0 ¼ 25, t1 ¼ 3:5

t0 ¼ 30, t1 ¼ 4:5

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

199.95 52.45 13.97 6.56 4.13 2.89 2.20 1.73 1.47 1.28 1.16

199.38 49.48 10.69 4.21 2.30 1.49 1.03 0.79 0.62 0.50 0.38

137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.95 52.39 13.98 6.57 4.12 2.89 2.20 1.74 1.47 1.28 1.16

199.41 49.40 10.69 4.22 2.30 1.49 1.03 0.79 0.62 0.50 0.38

137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

200.03 52.39 13.98 6.56 4.13 2.89 2.20 1.73 1.47 1.28 1.16

199.46 49.43 10.69 4.21 2.30 1.49 1.03 0.79 0.62 0.49 0.38

137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

200.03 52.39 13.98 6.56 4.13 2.89 2.20 1.73 1.47 1.28 1.16

199.46 49.43 10.69 4.21 2.30 1.49 1.03 0.79 0.62 0.49 0.38

137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

Sensitivity Analysis of Hyper-parameters (Scale) k0 ¼ 45, k1 ¼ 16 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

199.27 53.41 14.11 6.60 4.14 2.90 2.21 1.74 1.48 1.29 1.16

198.75 50.75 10.84 4.26 2.32 1.48 1.05 0.79 0.63 0.50 0.37

k0 ¼ 35, k1 ¼ 12 136.00 38.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.95 52.39 13.98 6.57 4.12 2.89 2.20 1.74 1.47 1.28 1.16

199.41 49.40 10.69 4.22 2.30 1.49 1.03 0.79 0.62 0.50 0.38

k0 ¼ 25, k1 ¼ 8 137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.39 51.41 13.81 6.56 4.07 2.87 2.23 1.77 1.47 1.29 1.16

200.06 48.56 10.40 4.22 2.23 1.48 1.04 0.79 0.61 0.49 0.37

k0 ¼ 20, k1 ¼ 6 136.00 36.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

198.46 50.12 13.54 6.47 4.04 2.86 2.19 1.78 1.48 1.29 1.16

198.23 46.76 10.22 4.17 2.26 1.46 1.06 0.80 0.61 0.49 0.37

136.00 36.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

198.75 48.64 10.47 4.22 2.27 1.49 1.02 0.79 0.61 0.48 0.37

138.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

Sensitivity Analysis of Hyper-parameters (Levy)

g0 ¼ 0:01 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

200.66 51.75 13.85 6.54 4.10 2.88 2.19 1.74 1.46 1.28 1.15

g0 ¼ 0:3 199.36 48.70 10.57 4.20 2.28 1.49 1.02 0.79 0.61 0.49 0.37

138.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.89 51.82 13.79 6.54 4.11 2.85 2.19 1.75 1.46 1.27 1.15

g0 ¼ 1:5 198.71 48.65 10.47 4.22 2.27 1.48 1.02 0.79 0.61 0.48 0.37

138.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.27 51.64 13.83 6.55 4.05 2.86 2.21 1.76 1.46 1.29 1.15

g0 ¼ 2:5 198.13 48.63 10.44 4.23 2.24 1.48 1.03 0.79 0.61 0.49 0.37

137.00 37.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.95 51.79 13.80 6.54 4.11 2.86 2.19 1.75 1.46 1.27 1.15

252

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

from expert to expert as the knowledge and expertise levels are not same. Many methods are available in the literature for elicitation of hyper-parameters (i.e., Garthwaite, Kadane, & O’Hagan, 2005; O’Hagan et al., 2006).

distributions. Elicitation is a technique used to add up expert opinion about parameters under study. In a more formal Bayesian way, it is the mode of selecting prior distributions for posterior models. The possible values of the hyper-parameters for priors may vary

Table 4 ARL, SDRL and MDRL values for sensitivity analysis of BHP levy priors for shifts in slope at ARL0 = 200. dS

Sensitivity Analysis of Hyper-parameters (Location) EWMA-BNCP-1 t0 ¼ 15, t1 ¼ 1:5

0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

EWMA-BNCP-2

EWMA-BNCP-3

EWMA-BNCP-4

t0 ¼ 20, t1 ¼ 2:5

t0 ¼ 25, t1 ¼ 3:5

t0 ¼ 30, t1 ¼ 4:5

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

199.95 87.27 32.23 14.69 8.72 5.84 4.31 3.40 2.78 2.32 1.98

199.38 87.02 28.18 11.31 5.71 3.64 2.44 1.83 1.39 1.13 0.92

137.00 61.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.95 87.20 32.20 14.65 8.75 5.84 4.31 3.40 2.77 2.31 1.99

199.41 87.02 28.16 11.27 5.78 3.64 2.44 1.82 1.40 1.13 0.91

137.00 61.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

200.03 87.14 32.20 14.65 8.75 5.84 4.31 3.41 2.77 2.32 1.98

199.46 86.96 28.17 11.27 5.77 3.64 2.44 1.82 1.39 1.13 0.91

137.00 61.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

200.03 87.14 32.18 14.65 8.75 5.85 4.31 3.40 2.78 2.31 1.99

199.46 86.96 28.16 11.26 5.77 3.64 2.44 1.82 1.40 1.13 0.91

137.00 61.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

Sensitivity Analysis of Hyper-parameters (Scale) k0 ¼ 45, k1 ¼ 16 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

199.27 89.40 32.72 14.91 8.86 5.83 4.39 3.44 2.80 2.30 1.99

k0 ¼ 35, k1 ¼ 12

198.75 89.15 28.63 11.53 5.84 3.56 2.50 1.83 1.42 1.13 0.92

136.00 62.00 25.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.95 87.20 32.20 14.65 8.75 5.84 4.31 3.40 2.77 2.31 1.99

k0 ¼ 25, k1 ¼ 8

199.41 87.02 28.16 11.27 5.78 3.64 2.44 1.82 1.40 1.13 0.91

137.00 61.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.39 84.00 31.52 14.48 8.60 5.88 4.29 3.36 2.77 2.33 2.01

200.06 82.32 27.43 11.00 5.65 3.50 2.39 1.78 1.39 1.12 0.92

k0 ¼ 20, k1 ¼ 6 136.00 59.00 24.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

198.46 79.97 30.37 14.29 8.56 5.74 4.30 3.39 2.71 2.30 1.99

198.23 77.49 26.01 10.64 5.76 3.46 2.37 1.77 1.37 1.09 0.90

136.00 56.00 23.00 12.00 7.00 5.00 4.00 3.00 2.00 2.00 2.00

198.75 84.72 28.10 11.21 5.67 3.60 2.39 1.78 1.39 1.12 0.93

138.00 60.00 24.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

Sensitivity Analysis of Hyper-parameters (Levy)

g0 ¼ 0:01 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

200.66 86.39 31.88 14.50 8.65 5.82 4.26 3.37 2.75 2.31 1.98

g0 ¼ 0:3 199.36 85.04 27.86 11.11 5.72 3.59 2.43 1.80 1.38 1.12 0.93

138.00 60.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.89 85.81 31.96 14.59 8.60 5.86 4.25 3.35 2.75 2.32 1.98

g0 ¼ 1:5 198.71 84.70 28.04 11.22 5.68 3.60 2.40 1.78 1.38 1.13 0.93

138.00 60.00 24.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.27 85.55 31.89 14.62 8.59 5.82 4.25 3.37 2.75 2.33 1.98

g0 ¼ 2:5 198.13 84.45 27.81 11.21 5.66 3.54 2.42 1.79 1.39 1.11 0.93

137.00 59.00 24.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

199.95 85.86 32.02 14.57 8.59 5.85 4.25 3.36 2.75 2.31 1.98

Table 5 Overall and individual in-control ARL for EWMA control charts of classical and Bayesian setups. Overall In-Control ARL of 200

Overall In-Control ARL of 370

Overall In-Control ARL of 500

Classical

Bayesian-BNCP

Bayesian-BCP

Classical

Bayesian-BNCP

Bayesian-BCP

Classical

Bayesian-BNCP

Bayesian-BCP

Log scale LI In-Control ARL LS In-Control ARL LE In-Control ARL

3.0156 587.5 3.0109 577.3 1.3723 701.8

3.0344 585.7 3.0307 548.1 1.8500 700.9

3.5872 494.2 3.2209 1023.9 1.9623 498.3

3.2046 1073.8 3.2051 1074.5 1.4623 1186.7

3.2454 1125.9 3.2807 1198.2 1.942 1064.1

3.8572 1065.3 3.3209 1438.5 2.1223 958.2

3.2586 1281.5 3.2636 1302.1 1.5623 2267.4

3.3214 1429.5 3.3107 1322.3 2.0573 1850.2

3.9382 1377.2 3.4039 1911.5 2.2063 1375.8

Natural scale LI In-Control ARL LS In-Control ARL LE In-Control ARL

3.0156 587.5 3.0109 577.3 4.200 695.1

3.0344 585.7 3.0307 548.1 5.6620 673.5

3.5872 493.2 3.2209 1023.9 4.3079 480.9

3.1950 1041.1 3.1800 990.9 4.6540 1346.5

3.2464 1128.7 3.2807 1198.3 5.8920 1028

3.8572 1065.3 3.3209 1438.5 4.8250 955.8

3.2550 1266.3 3.2800 1379.4 4.9490 2077

3.3214 1429.5 3.3107 1322.3 6.2185 1854.3

3.9382 1377.2 3.4039 1911.5 5.0700 1366.3

253

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

The method applied here is the one introduce by Garthwaitea, Al-Awadhib, Elfadalya, and Jenkinsonc (2013). They focus on models which entail variables that follow Normal, Binomial, Poisson or Exponential distributions. They exploit Piece-wise linear function to sculpt association between every continuous independent and dependent variable as it acquires elastic fallout that can correspond to broad diversity of opinions. Their method seems appealing as they use interactive computer program to elicit beliefs. The elicitation method is based on conditional mean priors by Bedrick, Christensen, and Johnson (1996). Various locations in predictor space are specified and expert assesses prior distributions for mean of potential observations at each location. The assessment of median and two quartiles is divided into three stages: first, asking question about median of expected response at reference knot. Knot is the point of predictor for which we are asked to assess medians of expected response. Then, asking question to assess the upper and lower quartiles at reference knot. These assessments are formulated via method of bisection, used and described by Huber (1974) and Garthwaite and Dickey (1985). The following steps can be employed for assessment.  Characterize the response variables, its minimum and maximum values and finally identify type of regression model i.e. ordinary linear regression.  Nature and amount of explanatory variables, their minimum and maximum value, followed by reference value and reference knots for all predictors.

 Assess median and two quartiles of expected response at reference value of the explanatory variable.  Evaluate median values for the expected response at each reference knot of the explanatory variable.  We assumed median assessment for the expected response at reference values of the predictor, and then asses lower and upper quartile of the expected response at reference values of the explanatory variable and repeated the process up to the last reference knot. Lastly, for the elicitation of hyper-parameters of the inverse gamma prior many methods are available in literature (i.e., Garthwaite & Dickey, 1988; Garthwaite & Dickey, 1992). The methods applied in the literature mentioned above have some drawbacks, that is why for this study the method of Elfadaly and Garthwaite (2014) is used, which is an extension of Garthwaite and Dickey (1988) techniques. The method calls for estimating an experiment and acquiring two mean responses of the expected response for fixed value of the regressor. However, randomness may trigger mean response to vary for any given value of the regressor. Repeat the above experiment several times and assess median differences of expected responses of dependent variable. Repeating and guessing median of differences of mean responses of dependent variable triumph over the possibility of over or under estimation of median assessment. After ending up with assessment, resultant hyper-parameters for the prior distributions of the Y-intercept, the slope and the

Table 6 ARL, SDRL and MDRL comparisons for classical and Bayesian EWMA control charts under shifts in. Y-intercept. dI

ARL0 = 200 EWMA-C

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

EWMA-BNCP

EWMA-BCP

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

207.43 59.90 16.22 7.92 5.16 3.82 3.08 2.60 2.28 2.06 1.88

201.38 54.09 11.52 4.17 2.20 1.34 0.97 0.73 0.56 0.46 0.43

145.00 44.00 13.00 7.00 5.00 4.00 3.00 2.00 2.00 2.00 2.00

199.67 50.56 13.66 6.53 4.07 2.85 2.23 1.76 1.49 1.28 1.16

198.82 47.63 10.25 4.18 2.23 1.48 1.06 0.78 0.62 0.48 0.38

137.00 36.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

199.62 38.69 12.39 6.47 4.18 3.04 2.36 1.89 1.63 1.39 1.23

195.03 33.97 8.49 3.54 2.06 1.41 1.03 0.80 0.66 0.53 0.43

138.00 28.00 10.00 6.00 4.00 3.00 2.00 2.00 2.00 1.00 1.00

361.81 80.75 13.93 4.77 2.39 1.47 1.01 0.76 0.60 0.46 0.41

257.00 62.00 15.50 8.00 5.00 4.00 3.00 3.00 2.00 2.00 2.00

370.08 77.08 17.53 7.71 4.68 3.16 2.43 1.92 1.60 1.38 1.21

362.47 72.65 13.39 4.72 2.52 1.58 1.14 0.84 0.66 0.55 0.43

257.00 54.00 14.00 7.00 4.00 3.00 2.00 2.00 2.00 1.00 1.00

369.61 62.94 16.57 7.93 5.00 3.52 2.72 2.17 1.81 1.56 1.36

358.86 56.69 11.63 4.29 2.44 1.60 1.13 0.88 0.71 0.61 0.51

260.00 46.00 14.00 7.00 5.00 3.00 3.00 2.00 2.00 2.00 1.00

491.04 96.19 15.18 5.02 2.44 1.48 1.03 0.77 0.62 0.49 0.41

349.00 70.00 16.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00 2.00

500.64 88.43 18.93 8.22 4.86 3.27 2.47 1.98 1.65 1.39 1.25

505.90 84.77 14.59 5.18 2.59 1.61 1.13 0.88 0.68 0.55 0.45

347.00 62.00 15.00 7.00 4.00 3.00 2.00 2.00 2.00 1.00 1.00

500.14 71.23 17.74 8.37 5.15 3.62 2.79 2.24 1.86 1.60 1.38

498.03 65.24 12.60 4.74 2.40 1.61 1.16 0.92 0.72 0.62 0.52

350.00 51.00 14.00 7.00 5.00 3.00 3.00 2.00 2.00 2.00 1.00

ARL0 = 370 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

369.95 86.73 19.39 8.96 5.61 4.11 3.29 2.76 2.41 2.16 1.99 ARL0 = 500

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

499.92 100.63 20.79 9.23 5.77 4.19 3.35 2.79 2.43 2.19 2.01

254

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

errors variance are as: B0  Nð20; 35Þ, B1  Nð2:5; 12Þ, and r2  IGð0:68; 0:3Þ:

3.2. Extra Quadratic Loss (EQL) As far as the overall performance and usefulness of control charts is concerned, EQL is the substitute that straightforwardly evaluates expected loss because of poor quality. Mathematically, it is computed as:

3. Performance measures and comparison This section highlights some measures used for efficiency comparison of performance in EWMA control charts for profile monitoring, under classical and Bayesian setups. The formal description of these measures is as follows:

3.1. Average Run Length (ARL) The set of sample points that required detecting out-ofcontrol signal is called the ‘‘run”. The amount of samples in use for the run is called ‘‘run length”. Run length is computed with respect to mean, median and standard deviation as ARL (average of Run Length distribution), MDRL (Median of Run Length distribution) and SDRL (Standard Deviation of Run Length distribution), respectively. These are regularly used measures of control charts on performance comparisons, when shifts take place in process parameters at some discrete point. The in-control ARL are calculated by 1=u where u is the probability of type-I error or false alarm rate. The out-of-control ARL are calculated by 1=ð1  /Þ where / is the probability of type-II error.

EQL ¼

1 dmax  dmin

Z

dmax

d2 ARLðdÞdd;

ð31Þ

dmin

where, the mathematical formulation of EQL defined it as the weighted average of ARLs for the entire domain of shifts (i.e., dmin < d < dmax ), while the square of the shifts (d2 ) as a weight. 3.3. Relative Average Run Length (RARL) It is another measure that checks the overall efficiency of control charts comparative to yardstick for each shift in process by viewing the intimacy of ARL to that yardstick. The decision for the yardstick chart is based on the minimum value of EQL. Mathematically it is defined as:

RARL ¼

1 dmax  dmin

Z

dmax

dmin

ARLðdÞ dd; ARLbmk ðdÞ

ð32Þ

where, ARLbmk stand for average run length of the bench mark (bmk) control chart. Bench mark control chart is the chart with smaller ARL values. The RARL values of the bench mark chart equals 1, while

Table 7 ARL, SDRL and MDRL comparisons for classical and Bayesian EWMA control charts under shifts in slope. dS

ARL0 = 200 EWMA-C

0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

EWMA-BNCP

EWMA-BCP

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

207.43 104.12 36.72 17.20 10.17 7.10 5.49 4.44 3.75 3.26 2.91

201.38 98.46 31.29 12.25 5.94 3.43 2.33 1.71 1.31 1.03 0.85

145.00 74.00 27.00 14.00 9.00 6.00 5.00 4.00 4.00 3.00 3.00

199.67 81.01 30.66 14.45 8.50 5.83 4.31 3.35 2.72 2.31 1.99

198.82 79.15 26.32 10.85 5.63 3.47 2.39 1.79 1.34 1.10 0.92

137.00 57.00 23.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00

199.62 67.58 26.64 13.75 8.53 6.05 4.62 3.65 3.00 2.55 2.19

195.03 63.63 21.65 9.44 5.23 3.28 2.42 1.74 1.38 1.12 0.94

138.00 48.00 20.00 12.00 8.00 5.00 4.00 3.00 3.00 2.00 2.00

361.81 160.77 44.00 16.00 7.23 4.02 2.56 1.83 1.40 1.10 0.91

257.00 118.00 37.00 16.00 10.00 7.00 6.00 5.00 4.00 3.00 3.00

370.08 140.29 43.10 18.29 10.15 6.93 4.91 3.86 3.08 2.57 2.18

362.47 134.97 38.57 14.32 6.69 4.11 2.73 2.00 1.51 1.21 0.99

257.00 98.00 32.00 15.00 9.00 6.00 4.00 4.00 3.00 2.00 2.00

369.61 117.55 37.04 17.58 10.38 7.02 5.31 4.20 3.44 2.85 2.46

358.86 111.50 31.33 12.56 6.32 3.75 2.59 1.98 1.53 1.20 1.01

260.00 84.00 28.00 14.00 9.00 6.00 5.00 4.00 3.00 3.00 2.00

491.04 194.98 49.19 17.00 7.65 4.12 2.72 1.92 1.42 1.11 0.94

349.00 141.00 42.00 18.00 10.00 7.00 6.00 5.00 4.00 3.00 3.00

500.64 169.82 48.82 19.53 10.96 7.12 5.19 3.98 3.23 2.68 2.27

505.90 171.00 44.13 15.36 7.43 4.19 2.78 2.02 1.59 1.26 1.03

347.00 118.00 36.00 15.00 9.00 6.00 5.00 4.00 3.00 3.00 2.00

500.14 141.77 43.09 19.12 11.06 7.50 5.48 4.31 3.55 2.96 2.54

492.04 137.59 36.89 13.99 6.77 3.99 2.68 1.99 1.56 1.25 1.04

354.00 100.00 32.00 15.00 10.00 7.00 5.00 4.00 3.00 3.00 2.00

ARL0 = 370 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

369.95 166.97 50.30 21.22 11.88 7.99 6.03 4.84 4.07 3.50 3.11 ARL0 = 500

0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250

499.92 202.10 56.35 22.81 12.49 8.22 6.19 4.95 4.16 3.58 3.16

255

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

for the other control charts it is greater than 1. The RARL value greater than 1 for a particular chart shows the inferiority of control chart in performance compared to the bench mark control chart. 3.4. Performance Comparison Index (PCI) Incorporates the performance and ranking based on average of EQL (AEQL). The PCI values are computed by taking the ratio of the AEQL of a control chart and the AEQL of the best chart under identical conditions. The PCI is equal to one when it has lowest AEQL while all other charts have greater than one. The mathematical form of PCI is as:

PCI ¼

AEQL AEQLbestchart

ð33Þ

These performance measures with some modifications may be seen in literature (e.g., Ou, Wu, & Goh, 2011; Ou, Wu, & Tsung, 2012; Riaz, Abbasi, Ahmad, & Zaman, 2014; Ryu, Wan, & Kims, 2010; Wu, Jiao, Yang, Liu, & Wang, 2009; Zaman, Riaz, Ahmad, & Abbasi, 2015). The next subsections will present a comprehensive comparison of this proposed Bayesian approach with the approach existing in the literature. We have designed three EWMA control charts for the Y-intercept, the slope and the errors variance to compare the performance of proposed Bayesian approach with the mentioned classical approach. The value of smoothing constant ‘k’ is set equivalent to 0.2 as it is well known that smaller the value of ‘k’ better

the performance of control chart (Lucas & Saccucci, 1990). We consider fixed X-values of 2, 4, 6, and 8, and then take the deviation  ¼ 5). The known values of Y-intercept and slope from mean (X are set equal to 13 and 2, respectively under the model in Eq. (4). The error terms are assumed to follow standard normal distribution. 3.5. Sensitivity analysis of hyper-parameters In this section, it is presented the sensitivity analysis of the hyper-parameters, while performing the Bayesian analysis using informative priors. Our study comprises different informative (non-conjugate and conjugate) priors for process parameters when the process is defined by a linear profile. For the values of hyperparameters (location, scale) of normal and BHP prior distributions proceed with the values provided above by using the method of Garthwaitea et al. (2013), and for the values of hyper-parameters of inverse gamma and levy prior distribution initiate with values obtained by using the method of Elfadaly and Garthwaite (2014) above. These different values of hyper-parameters calculated above give some ideas about the possible location and scale values of the prior distributions. Now, we consider the different values of hyper-parameters for sensitivity analysis and check the possible impact on the values of different performance measures. We monitor the impact of increase and decrease in the values of hyperparameters on the performance of control charting structure under Bayesian scheme of linear profile.

Table 8 ARL, SDRL and MDRL comparisons for classical and Bayesian EWMA control charts under shifts in errors variance at log scale. dE

ARL0 = 200 EWMA-C

EWMA-BNCP

EWMA-BCP

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

207.43 35.09 13.09 7.43 5.17 4.02 3.34 2.87 2.54 2.30 2.09

201.38 30.96 10.40 5.13 3.18 2.24 1.79 1.43 1.25 1.08 0.96

145.00 26.00 10.00 6.00 4.00 3.00 3.00 3.00 2.00 2.00 2.00

199.67 32.19 10.48 5.39 3.39 2.54 2.01 1.73 1.53 1.40 1.31

198.82 32.83 10.14 5.08 2.93 2.07 1.46 1.16 0.93 0.80 0.68

137.00 22.00 7.00 4.00 2.00 2.00 1.00 1.00 1.00 1.00 1.00

199.62 32.13 10.74 5.73 3.68 2.73 2.18 1.81 1.65 1.47 1.37

195.03 31.64 9.90 5.14 3.13 2.22 1.63 1.23 1.07 0.86 0.72

138.00 22.00 8.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

369.95 49.10 16.41 8.60 5.88 4.49 3.69 3.13 2.76 2.48 2.26

361.81 44.91 13.24 6.11 3.70 2.59 1.99 1.56 1.31 1.15 1.03

257.00 35.00 12.00 7.00 5.00 4.00 3.00 3.00 2.00 2.00 2.00

370.08 45.72 13.77 6.45 4.08 2.89 2.26 1.91 1.69 1.51 1.38

362.47 46.87 13.32 5.79 3.59 2.37 1.70 1.36 1.12 0.93 0.75

257.00 31.00 10.00 5.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

369.61 44.09 13.75 6.91 4.36 3.18 2.42 2.09 1.81 1.62 1.45

358.86 43.48 12.73 6.03 3.68 2.55 1.83 1.49 1.24 1.03 0.83

260.00 30.00 10.00 5.00 3.00 2.00 2.00 2.00 1.00 1.00 1.00

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

499.92 63.56 19.12 9.75 6.46 4.87 3.98 3.35 2.94 2.65 2.41

491.04 59.50 16.07 7.03 4.05 2.85 2.11 1.66 1.37 1.24 1.08

349.00 45.00 14.00 8.00 5.00 4.00 3.00 3.00 3.00 2.00 2.00

500.64 57.54 16.47 7.46 4.45 3.10 2.41 2.00 1.77 1.56 1.46

505.90 56.03 15.71 6.87 3.83 2.59 1.80 1.45 1.17 0.94 0.81

347.00 41.00 12.00 5.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

500.14 56.33 16.73 7.71 4.79 3.33 2.61 2.13 1.86 1.66 1.52

498.03 54.91 15.59 6.84 3.95 2.73 2.00 1.54 1.24 1.06 0.88

350.00 40.00 12.00 6.00 4.00 2.00 2.00 2.00 1.00 1.00 1.00

ARL0 = 370

ARL0 = 500

256

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

All the control charts are designed to have in-control ARL of 200 with sample size n ¼ 4. The right combinations of LI ; LS ; LE ; LE (the control limit coefficients for the Y-intercept, the slope, and the errors variance under log and natural scale that determine the width of the control limits) are selected to have in-control ARL of 200. We have evaluated the performance measures using the simulations of 10,000 iterations for shifts in the Y-intercept, slope and errors variance under log and natural scales. Tables 1–4 summarize the results of ARL, SDRL and MDRL for different values of hyper-parameters of non-conjugate and conjugate priors for the impact on the Y-intercept and slope. The impact it

showed on the performance structure of Bayesian EWMA control charts of process parameters as follows:  The values of ARL, SDRL and MDRL for the Y-intercept and slope increase for normal priors, while remain same for BHP priors with the decrease in the location hyper-parameters.  There is continuous decrease in the values of ARL, SDRL and MDRL for the Y-intercept and slope with increase in the location hyper-parameters for normal priors. After certain increase in the values of the location hyper-parameters the values of performance measures remain almost the same.

Table 9 ARL, SDRL and MDRL comparisons for classical and Bayesian EWMA control charts under shifts in errors variance at natural scale. dE

ARL0 = 200 EWMA-C

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

EWMA-BNCP

EWMA-BCP

ARL

SDRL

MDRL

ARL

SDRL

MDRL

ARL

SDRL

MDRL

207.02 33.32 12.00 6.69 4.61 3.44 2.83 2.38 2.09 1.87 1.71

200.49 29.72 9.59 4.80 3.08 2.16 1.72 1.40 1.21 1.05 0.93

145.00 24.00 9.00 5.00 4.00 3.00 2.00 2.00 2.00 2.00 1.00

199.62 31.79 11.33 6.01 3.93 2.89 2.27 1.91 1.67 1.53 1.40

196.79 28.74 9.37 4.96 3.14 2.22 1.63 1.31 1.08 0.92 0.77

139.00 24.00 9.00 5.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

199.52 29.92 10.05 5.35 3.50 2.57 2.10 1.79 1.56 1.43 1.33

196.13 29.20 9.26 4.80 3.00 2.00 1.61 1.21 0.97 0.82 0.69

139.00 21.00 7.00 4.00 3.00 2.00 1.00 1.00 1.00 1.00 1.00

363.51 43.69 12.23 5.75 3.49 2.49 1.95 1.59 1.32 1.13 0.99

257.00 34.00 12.00 6.00 4.00 3.00 3.00 2.00 2.00 2.00 2.00

370.08 43.38 14.05 7.25 4.58 3.27 2.53 2.11 1.84 1.62 1.48

366.84 39.27 11.35 5.68 3.61 2.50 1.88 1.51 1.25 1.03 0.87

255.00 32.00 11.00 6.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

369.68 41.32 12.72 6.42 4.05 2.90 2.31 1.96 1.68 1.52 1.42

364.76 40.67 11.70 5.58 3.38 2.32 1.78 1.39 1.12 0.90 0.80

253.50 29.00 9.00 5.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

495.59 55.29 14.60 6.38 3.78 2.60 2.06 1.64 1.40 1.20 1.06

348.00 42.00 13.00 7.00 5.00 4.00 3.00 2.00 2.00 2.00 2.00

499.56 53.58 16.23 8.09 4.95 3.49 2.64 2.13 1.87 1.64 1.51

502.15 48.65 13.37 6.51 3.99 2.72 2.00 1.51 1.26 1.04 0.89

347.00 39.00 13.00 7.00 4.00 3.00 2.00 2.00 1.00 1.00 1.00

499.37 50.95 14.32 7.05 4.42 3.17 2.41 2.03 1.76 1.55 1.44

507.87 48.79 13.32 6.10 3.63 2.64 1.88 1.46 1.22 0.95 0.83

345.00 36.00 10.00 5.00 3.00 2.00 2.00 1.00 1.00 1.00 1.00

ARL0 = 370 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

369.68 47.49 15.15 7.88 5.15 3.90 3.12 2.62 2.27 1.98 1.81

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00

499.94 59.02 17.34 8.61 5.62 4.10 3.31 2.71 2.37 2.09 1.89

ARL0 = 500

Table 10 EQL, RARL and PCI values for classical & Bayesian EWMA control charts for linear Profile at ARL0 = 200, 370 and 500. ARL0

Charts

Intercept

Slope

Error (log Scale)

Error (Natural Scale)

EQL

RARL

PCI

EQL

RARL

PCI

EQL

RARL

PCI

EQL

RARL

PCI

200

EWMA-C EWMA-BNCP EWMA-BCP

4.0767 2.8956 2.9814

1.3521 1.0000 1.0099

1.4079 1.0000 1.0296

7.2431 5.5818 5.7554

1.2727 1.0000 1.0114

1.2976 1.0000 1.0311

4.2071 2.6972 2.8565

1.4766 1.0000 1.0528

1.5598 1.0000 1.0591

3.5621 2.9698 2.7360

1.2700 1.0876 1.0000

1.3019 1.0855 1.0000

370

EWMA-C EWMA-BNCP EWMA-BCP

4.4793 3.3068 3.5322

1.3078 1.0000 1.0545

1.3545 1.0000 1.0681

8.4108 6.7726 6.9459

1.2219 1.0000 1.0153

1.2419 1.0000 1.0256

4.7008 3.0832 3.2785

1.4459 1.0000 1.0537

1.5246 1.0000 1.0633

3.9831 3.3492 3.0725

1.2663 1.0897 1.0000

1.2964 1.0901 1.0000

500

EWMA-C EWMA-BNCP EWMA-BCP

4.6186 3.4571 3.6730

1.2945 1.0000 1.0518

1.3360 1.0000 1.0625

8.8710 7.2536 7.4124

1.2062 1.0000 1.0160

1.2230 1.0000 1.0219

5.1346 3.3466 3.5214

1.4606 1.0000 1.0462

1.5342 1.0000 1.0522

4.2628 3.5408 3.2723

1.2716 1.0843 1.0000

1.3027 1.0820 1.0000

257

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

 The increase in the location hyper-parameters for BHP priors does not show considerable impact on the ARL, SDRL and MDRL values for the Y-intercept and slope shifts.  We find an increase in the values of individual performance measures for the Y-intercept and slope shifts with increase in the scale hyper-parameters for normal priors with low magnitude of difference.  We observe decrease in ARL, SDRL and MDRL values with decrease in the scale hyper-parameters for normal priors when there are shifts in Y-intercept and slope, while this decrease is negligible after certain decrease in the scale hyper-parameters values.  There is a slight increase or decrease in the values of ARL, SDRL and MDRL for shifts in the Y-intercept and slope with increase or decrease in the scale hyper-parameters for BHP priors.

(a)

 The performance of Bayesian EWMA control charts for the Y-intercept and slope shifts in terms of individual performance measures shows a potent behavior for increase or decrease in the values of hyper-parameters for inverse gamma prior.  There is no increase and decrease in the values of ARL, SDRL and MDRL with the increase or decrease in hyper-parameters under levy prior for shifts in the Y-intercept and slope. The sensitivity analysis for shifts in the errors variance (log and natural scales) is also performed for hyper-parameters of nonconjugate and conjugate priors. The values are not reported here as the behavior of the different performance measures for increase and decrease in the values of hyper-parameters is almost identical. The magnitude of difference in the values of individual performance measures for shifts in errors variance is negligible for difference values of hyper-parameters. The hyper-parameters of

70

EWMA-C

60

EWMA-BNCP EWMA- BCP

ARL

50 40 30 20 10 0

0.20

0.40

0.60

0.80

1.00

δI

(b)

100

EWMA-C

90

EWMA-BNCP

80

EWMA- BCP

ARL

70 60 50 40 30 20 10 0

0.20

0.40

0.60

0.80

1.00

δI

(c)

120

EWMA-C EWMA-BNCP

100

EWMA- BCP

ARL

80 60 40 20 0

0.20

0.40

0.60

0.80

δI Fig. 2. ARL curves for Y-intercept at (a) ARL = 200; (b) ARL = 370; (c) ARL = 500.

1.00

258

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

non-conjugate and conjugate priors with lowest values of performance measures in Tables 1–4 are considered for further analysis of Bayesian EWMA control charting structures for in-control ARL of 200, 370 and 500 under log and natural scales. 3.6. Comparisons In this subsection, we provide comprehensive comparisons of proposed Bayesian EWMA control charting structure of process parameters, with that of existing classical approach. The performance of the control charting structures is compared in terms of individual performance measures (e.g., ARL, SDRL, and MDRL) and overall performance measures (e.g., EQL, RARL, and PCI). The set of control charting structures for linear profile monitoring considered for the comparisons consists of three univariate EWMA control charts for classical approach as EWMA-C and that of Bayesian approach by EWMA-BNCP under non-conjugate priors and EWMA-BCP under conjugate priors. The EWMA-C, EWMA-BNCP,

(a)

and EWMA-BCP control charts are constructed to have same incontrol ARL when measured independently under log and natural scales. The EWMA control charts of process parameters under classical and Bayesian setups are designed in such a way that the combination of three univariate EWMA control charts has same in-control ARLs of 200, 370, and 500. Table 5 presents the right combinations of LI ; LS ; LE ; LE for individual and overall in-control ARL under classical and Bayesian setups. The Y-intercept shifts, the slope shifts and the errors variance shifts in units of standard deviations for log and natural scales are incorporated in our simulation study. The Tables 6–9, show the behavior of ARL, SDRL, and MDRL at range of values of shifts based on in-control ARL values of 200, 370 and 500. The shifts in the Y-intercept from B0 To B0 þ dI r, shifts in the slope from B1 To B1 þ dS r and shifts in the errors variance from r To dE r are incorporated in this study. In Tables 6 and 7, the in-and out-of-control values of ARL, MDRL and SDRL for shifts in the Y-intercept and slope under EWMA-C,

120

EWMA-C EWMA-BNCP

100

EWMA-BCP

ARL

80 60 40 20 0

0.025

0.050

0.075

0.100

0.125

δS

(b)

180

EWMA-C

160

EWMA-BNCP

140

EWMA-BCP

ARL

120 100 80 60 40 20 0

0.025

0.050

0.075

0.100

0.125

δS

(c)

250

EWMA-C EWMA-BNCP

200

EWMA-BCP

ARL

150

100

50

0

0.025

0.050

0.075

0.100

δS Fig. 3. ARL curves for slope at (a) ARL = 200; (b) ARL = 370; (c) ARL = 500.

0.125

259

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

EWMA-BNCP, and EWMA-BCP control charts reflect that proposed Bayesian approach shows quick detection of out-of-control signals compared to classical approach. It is also observed that EWMA-BCP control charts perform better than EWMA-BNCP control charts for small shifts, while EWMA-BNCP control charts have slight advantage over EWMA-BCP control charts for larger shifts. This indicates that Bayesian analysis under conjugate priors for profile monitoring is preferred over non-conjugate priors used in this study for small to moderate values of shifts. The tendency of poor performance for larger shifts appears for monitoring the linear profile as EWMA control charts are not sensitive to larger shifts. In this situation, Shewhart type control charts can be used with EWMA control charts as an alternative to enhance the performance ability of the control charts. The parallel performance measures of SDRL and MDRL shows the same pattern with that of ARL. The out-of-control SDRL and MDRL values for shifts in the Y-intercept and slope show quicker detection in sustainable shifts than that of ARL values.

(a)

The above results indicate that reducing significance levels increases the efficiency of EWMA-BNCP and EWMA-BCP control charts. There is an increase in the efficiency of EWMA-BCP control charts with the decrease in significance level. This shows that Bayesian profile monitoring is more effectual and more statistically considerable. Pragmatically, the overall significance in the performance of EWMA-BCP control charts improves as we enlarge incontrol value of ARL or reduce the probability of the type-I error. Tables 8 and 9 provide the in-and out-of-control ARL, SDRL, and MDRL values when there are shifts in errors variance under log and natural scales in units of standard deviation. The out-of-control ARL, SDRL, and MDRL values reveal that EWMA-BNCP and EWMA-BCP control charts perform better than EWMA-C control charts but the magnitude of difference is low. It is also observed that EWMA-BNCP control charts show slightly better performance over EWMA-BCP control charts under log scales, while EWMA-BCP control charts are better in case of natural scales. This indicates

40 EWMA-C

35

EWMA-BNCP

30

EWMA-BCP

ARL

25 20 15 10 5 0

1.20

1.40

1.60

1.80

2.00

δE

(b)

60

EWMA-C EWMA-BNCP

50

EWMA- BCP

ARL

40 30 20 10 0

1.20

1.40

1.60

1.80

2.00

δE

(c)

70

EWMA-C

60

EWMA-BNCP EWMA- BCP

ARL

50 40 30 20 10 0

1.20

1.40

1.60

1.80

2.00

δE Fig. 4. ARL curves for Var(MSE) at log scale at (a) ARL = 200; (b) ARL = 370; (c) ARL = 500.

260

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

(a)

35 EWMA-C

30

EWMA-BNCP EWMA-BCP

25

ARL

20 15 10 5 0

1.20

1.40

1.60

1.80

2.00

δE

(b)

50 EWMA-C

45

EWMA-BNCP

40

EWMA-BCP

35

ARL

30 25 20 15 10 5 0

1.20

1.40

1.60

1.80

2.00

δE

(c)

70 EWMA-C

60

EWMA-BNCP EWMA-BCP

50

ARL

40 30 20 10 0

1.20

1.40

1.60

1.80

2.00

δE Fig. 5. ARL curves for Var(MSE) at natural scale at (a) ARL = 200; (b) ARL = 370; (c) ARL = 500.

Table 11 ARL values for Shewhart control charts schemes (Roberts Vs Bayesian) at ARL0 = 740. Roberts schemes

Bayesian schemes

Shifts

Rm/n

RSh

Mk

Gr

Ch

GRDd

BS-BHP

BS-N

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

740 79 19 7.9 4.4 2.7 1.9 1.5

740 50 11.9 5.7 3.8 2.9 2.4 2.2

740 40 10.2 6 4.6 3.8 3.3 2.9

740 40 10.1 5.1 3.5 2.7 2.2 1.9

740 34 10 5.8 4.3 3.5 2.9 2.5

740 32 9.2 5.2 3.7 2.9 2.4 2.1

740.58 25.28 2.78 1.18 1.01 1.00 1.00 1.00

740.33 19.98 2.49 1.13 1.01 1.00 1.00 1.00

that change of scale has improved the performance of EWMA-BCP control charts for monitoring the shifts in errors variance. The proposed structure exhibits limited viability in performance as the standard deviation shifts are usually monitored by the R-charts. It is due to the fact that increase in terms of standard deviation

increases the variation of the estimators of the Y-intercept and slope, then smoothing these estimators with EWMA control charts are not appropriate. The individual performance measures comparisons in Tables 8 and 9 shows that the natural scale performs slightly better than log scale for EWMA-BNCP, EWMA-BCP and EWMA-C control charts. The change of scale from log to natural provides similar values of individual performance measures for shifts in the Y-intercept and slope. It means change of scale has no effect on the monitoring of the Y-intercept and slope. The values of ARL, SDRL, and MDRL for the shifts in Y-intercept and slope under natural scale are not provided here because they are identical with that of log scale (see Tables 6 and 7). The in-and out-of-control values of ARL, SDRL, and MDRL in Tables 6–9 reflect the superiority of proposed Bayesian methods with faster detection of sustainable shifts in the parameters of the underlying model in Eq. (4). It can be concluded that EWMA-BCP control charts dominate all other control charts to detect positive smaller-moderate value shifts

261

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

ARL

in the Y-intercept and the slope parameters of the regression model. The above performance measures of ARL, SDRL, and MDRL in Tables 6–9 are the individual performance measures of EWMA-C, EWMA-BNCP, and EWMA-BCP control charts under classical and Bayesian setups. Table 10 provides the values of RARL, EQL, and PCI for in-control ARL of 200, 370 and 500 to measure the overall performance of EWMA-C, EWMA-BNCP, and EWMA-BCP control charts under classical and Bayesian setups. The performance evaluations based on the ARL are at specific points for given control charting structures, whereas the alternative approach of EQL illustrates the overall efficiency of a control chart for given sequences of shifts. The control charts with smaller EQL shows better performance. It can be seen in Table 6 that the values of EQL for the Y-intercepts, the slopes and the error variances (log and natural scale) for the Bayesian approach are smaller than that of classical approach. This shows the superiority of Bayesian methods over classical methods. It is also observed that under Bayesian setup, the EQL values of EWMA-BNCP and EWMA-BCP control charts are roughly the same. This indicates that overall performance of EWMA-BNCP, and EWMA-BCP control charts are almost similar for the Y-intercept, slope and errors variance under log and natural scales. Another overall performance measure the Relative ARL (RARL) describes the overall effectiveness of a control chart relative to bench mark control chart. The control charts with lowest ARL values (see Tables 6–9) are considered the benchmark chart with the RARL value equal to 1. The RARL values greater than 1 for EWMA-C control charts indicate the inferiority of the classical approach compare to Bayesian approach, while the RARL values of EWMABNCP and EWMA-BCP control charts are almost same. The control charts with the lowest EQL have PCI value equal to one. The results in Table 6 indicate that the control charts under Bayesian approach for the Y-intercepts, slopes and error variances (log and natural scale) have PCI values almost equal to one. This also confirms that control charts under Bayesian approach are better than control charts under classical approach. Overall performance measures values are shown with the results of the individual performance measures of EWMA-C, EWMA-BNCP, and EWMA-BCP control charts under classical and Bayesian setups. This indicates that the overall performance of EWMA-BNCP and EWMA-BCP control charts is better than that of EWMA-C control charts for monitoring simple linear profiles. In brief, the inclusion of more information about the process parameters enhances the detection ability of the EWMA control structures. The selection of prior distributions and the values of hyper-parameters require careful consideration for an effective application of Bayesian control structures and the performances

of their corresponding performance measures i.e., ARL, SDRL, MDRL, RARL, EQL and PCI. We have presented the overall comparative view of EWMA-C, EWMA-BNCP, and EWMA-BCP control charts by some graphical displays of ARL curves using ARL values given in Tables 6–9. This graphical display shows the superiority of proposed EWMA-BNCP and EWMA-BCP control charts over the EWMA-C control charts, while EWMA-BCP control charts show better performance than the EWMA-BNCP control charts. Fig. 2(a–c) presents ARL values for in-control ARL 200, 370, and 500 respectively, when there are shifts in the Y-intercept in units of standard deviations. Fig. 2(a–c) indicates that EWMABCP control charts are at the lower side reflecting quicker detection of sustainable shifts in comparison to EWMA-C and EWMABNCP control charts. The performance of EWMA-BNCP and EWMA-BCP control charts is significant for smaller-moderate value shifts while almost negligible for larger value shifts. The insignificant behavior of EWMA-type control charts for large value shifts is a general phenomenon and can be handled with the induction of Shewhart-type control charts with each of the EWMA control charts. Fig. 3(a–c) provides the ARL curves for same in-control ARL of 200, 370, and 500 under the slope shifts in units of standard deviations. This tendency and behavior of ARL curves are very much similar with that of Fig. 2(a–c) for EWMA-C, EWMA-BNCP, and EWMA-BCP control charts. This shows that overall, EWMA-BCP control charts perform better than EWMA-C and EWMA-BNCP control charts for monitoring the slopes of X. The graphs in Figs. 4(a–c) and 5(a–c) provide the ARL curves for in-control ARL values of 200, 370, and 500 while considering increase in error variances under log and natural scales. These graphical representations indicate that the control charting structures under natural scale have slight advantage over the control charting structure under log scale. The display of the curves show that overall EWMA-BNCP and EWMA-BCP control charts perform equivalently, while EWMA-BCP control charts have slight advantage over EWMA-C control charts. It is also noted that change of scale has more affect on EWMA-BCP control charts compare to EWMA-BNCP control charts. Overall, control charting structures under Bayesian setup are efficient than that of competing setup. 3.7. Comparison with Roberts schemes of tests In this subsection we will see how well our proposed Bayesian scheme performs against the other Bayesian schemes available in literature. Consider the scheme of Roberts (1966) in which he tried to extend the different test procedures to monitor the sensitivity of process mean. These efforts to supplement the Shewhart test,

90

Rn/m

80

RSh

70

Mk

60

Gr

50

Ch

GRd

40

BS-BHP

30

BS-N

20 10 0 0.5

1

1.5

2

2.5

δ Fig. 6. ARL Comparisons for Shewhart Control Charts Schemes at ARL0 = 740.

3

262

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

various types of test procedures have been introduced i.e., Run Test (Rm/n), Run Sum Test (RSr), Moving Average Test (Mk), Geometric Moving Average Test (Gr), Cumulative Sum Test (Ch), and Girshick–Rubin Test (GRDd). These test procedures are based on

(a)

transformations applied at standard average control chart to monitor the sensitivity of process mean. This subsection considers a special case to compare the performance of proposed Bayesian schemes (using non-conjugate and EWMA-C EWMA-BNCP

14.1 13.9 13.7

UCLBNCP UCLC

Values

13.5 13.3 13.1 12.9 12.7

LCLC LCLBNCP

12.5 12.3 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 Subgroups

(b)

EWMA-C EWMA-BCP

14.1

Values

13.9 13.7

UCLBCP

13.5

UCLC

13.3 13.1 12.9 12.7 LCLC LCLBCP

12.5 12.3 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 Subgroups

Fig. 7. Control chart display of EWMA-C with (a) BNCP, (b) BCP charts for shift at Y-intercept.

2.5

(a)

EWMA-C EWMA-BNCP

2.4 UCLBNCP UCLC

2.3

Values

2.2 2.1 2.0 1.9 LCLC LCLBNCP

1.8 1.7 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Sub groups 2.5

(b)

EWMA-C EWMA-BCP

2.4 UCLBCP UCLC

2.3

Values

2.2 2.1 2.0 1.9

LCLC LCLBCP

1.8 1.7 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Sub groups Fig. 8. Control chart display of EWMA-C with (a) BNCP, (b) BCP charts for shift at slope.

263

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

conjugate priors) with the test procedures studied by Roberts (1966). Let us assume that the X-values are zero then the model in Eq. (4) becomes:

yij ¼ B0 þ eij :

ð34Þ

Then the expectation over Eq. (34) above results in Eðyij Þ ¼ B0 , which shows that the process monitoring of the Y-intercept is equivalent to monitoring of process mean. The Bayesian schemes under nonconjugate prior (BS-BHP) and conjugate prior (BS-N) are constructed to monitor the process mean. We consider the posterior estimates of the Y-intercept to construct the simple average control chart. The sample size (n ¼ 8) is taken for this study to compute the ARL values. For the average control chart to monitor the Y-intercept under Bayesian scheme using non-conjugate and conjugate prior, the values of the control limit coefficients are chosen 3.2084 and 3.3266, respectively to give the in-control ARL of 740. Table 11 and Fig. 6 summarize the simulation results of in- and out-of-control ARL values for shifts in the Y-intercept in units of standard deviations. The smaller out-of-control ARL values for the Bayesian schemes of BS-BHP and BS-N in Table 11 show superiority over all the schemes studied by Roberts (1966). The performance of BS-N scheme is better than BS-BHP scheme for small shifts while equivalent for larger values of shifts. Fig. 6 indicates that the BS-BHP and BS-N schemes are at the lower side reflecting quicker detection of sustainable shifts in process parameter compare to the schemes of Roberts (1966). The performance of BS-N scheme is significant and more practical. 4. Illustrative Example Let us consider an example to clarify our findings based on different performance measures for classical and Bayesian setups. Suppose that variables Y and X are defined as, Y = density of wood

board, and X = depths across the thickness of board. Some other real life example may be seen in literature (i.e., Keramatpoura et al., 2014; Zhang et al., 2013; Zhu & Lin, 2010). We have simulated 50 samples by using model in Eq. (4) above for three unknowns after considering non-conjugate and conjugate priors. We have in-control conditions where dI ¼ dS ¼ 0, for Yintercept and slope, and m0 = 30 (sample values). The remaining m1 = 20 sample values are obtained from out-of-control situation for Y-intercept shifts as dI ¼ 0:4, for slope shifts dS ¼ 0:125;. The in-control condition for errors variance where dE ¼ 1, while the out-of-control situation for errors variance shifts as d ¼ 1:8. We compute values of control charting statistics of EWMA-C, EWMABNCP, and EWMA-BCP. The control charts are constructed by plotting sample index at X-axis and the values of EWMA-C, EWMA-BNCP, and EWMA-BCP on Y-axis for the Y-intercepts in Fig. 6(a and b), for the slopes in Fig. 7(a and b) and for errors variance in Fig. 8(a and b). The bold line represents the values of EWMA-C and control limits, while the other line represents the control limits and values of EWMA-BNCP and EWMA-BCP. Fig. 7(a) indicates that for Y-intercepts, the EWMA-C control chart detects three out-of-control signals after 30th sample and the EWMA-BNCP control chart detects 12. Fig. 7(b) indicates that the EWMA-C control chart detects three out-of-control signals and that the EWMA-BCP control chart detects 14. The above results indicate that the EWMA-BCP control charts are superior to the EWMA-BNCP and EWMA-C control charts with regards to monitoring the Y-intercepts. As Fig. 8(a) shows, the EWMA-C control chart detects three outof-control signals for slope, whereas the EWMA-BNCP control charts detects nine. Fig. 8(b) shows that the EWMA-C control chart detects three out-of-control signals, whereas the EWMA-BCP control charts detects 12. The above results show that the EWMA-BCP control charts are superior to the EWMA-BNCP and EWMA-C control charts for monitoring the slopes.

(a)

1.4

EWMA-C EWMA-BNCP

1.2

Values

1.0 0.8 UCLC UCLBCNP

0.6 0.4 0.2 0.0 1

3

5

7

9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Sub groups

(b)

1.4

EWMA-C EWMA-BCP

1.2

Values

1.0 0.8 UCLC UCLBCP

0.6 0.4 0.2 0.0 1

3

5

7

9

11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

Sub groups Fig. 9. Control chart display of EWMA-C with (a) BNCP, (b) BCP charts for shift at V(MSE).

264

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

As Fig. 9(a) shows, the EWMA-C control chart detects two outof-control signals for variance of errors, whereas the EWMABNCP control chart detects 10. Fig. 9(b) shows that the EWMA-C control chart detects two out-of-control signals while the EWMA-BCP control chart detects nine. These results show that the EWMA-BCP and EWMA-BNCP control charts are better than the EWMA-C control charts for monitoring the errors variance. It can be observed from the above finding that acquiring more knowledge about the process structures in the form of prior distributions can yield tangible benefits. These benefits should be comparing with the cost on such informations. The detection ability of the EWMA control structures after shift has been introduced improves with the inclusion of prior informations. The weights of

conjugate priors enhance the efficiency of EWMA control structures more as compare to non conjugate priors. This indicates that normal priors are more compatible when the samples are drawn from normal distribution then that of non normal priors. The differences above in detection abilities conclude that the EWMABCP control charts are efficient at linear profile monitoring. 4.1. A case study The significance of propose Bayesian approach is further justified with application to real life data. For our study proposes, we used the data set based on the non isothermal continuous stirred tank chemical reactor model named as CSTR process originally

Ewma -C Ewma-BCP

1.3 1.2

UCLEB UCLEC

1.1

Values

1.0 0.9 0.8 0.7 0.6 LCLEB LCLEC

0.5

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101

0.4

Sub groups Fig. 10. Control chart display of EWMA-C with BCP charts for shift at Y-intercept.

Ewma -C

2.0

Ewma -BCP

1.5 UCLEB UCLEC

Values

1.0 0.5 0.0 -0.5

LCLEC LCLEB

-1.0

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101

-1.5

Sub groups Fig. 11. Control chart display of EWMA-C with BCP charts for shift at slope.

Ewma A-C

0.6

Ewma -BCP

0.5

0.3 UCLEB UCLEC

0.2 0.1 0.0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99 101

Values

0.4

Sub groups Fig. 12. Control chart display of EWMA-C with BCP charts for shift at V(MSE).

265

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

introduced by Marlin (2000) and can be widely seen in literature as benchmark in fault detection and diagnosis (i.e., Ahmad, Abbasi, Riaz, & Abbas, 2014; Shi, Lv, Fei, & Liang, 2013). The CSTR process data consist of nine process variables, from which, for our example of profile monitoring we have considered Inlet Concentration of Solvent Flow (CAS in Kmole/m3) as explanatory variables X while Outlet Concentration of the Product (CA in Kmole/m3) as response variable Y. The two variables originally contain N = 1024 observations collected over half minute sampling interval in which first N0 = 512 observations shows in-control scenario, while remaining N1 = 512 observations shows out-of-control scenario when the shift is introduced to monitor process parameters. For out-of-control scenario the shift of 0.3, 0.175 and 1.2 has been introduced for Yintercept, slope and errors variance, respectively. The bivariate data set is further divided into m = 102 profiles each of size n = 10 when each profile is collected over sampling interval of 5 min. The first m0 = 51 profiles shows in-control state, while the remaining m1 = 51 profiles shows out-of-control state when the shift is introduced to monitor the process parameters. The control limits for designed control charts are calculated with fixed ARL0 = 200 based on the in-control data, while the values of control limits coefficients are provided in Table 5. The value of smoothing constant ‘k’ is set equivalent to 0.2. The prior distributions for the Y-intercept, slope and errors variance with specified hyperparameters are selected as: B0  Nð16; 80:2Þ, B1  Nð1:05; 490Þ, and r2  IGð1:5; 0:3Þ, respectively. Using the aforementioned information the results are demonstrated in the form of three control charts for the Y-intercept, slope and errors variance by plotting the profiles on the horizontal axis and the values of statistics for EWMA-C and EWMA-BCP charts on the vertical axis as shown in Figs. 10–12. In figures, the solid lines represent the control limits and EWMA statistics for EWMA-C charts, while the dotted lines represent the control limits and EWMA statistics for EWMA-BCP charts. Fig.10 indicates that the EWMA-C control chart detects seven out-of-control signals, whereas the EWMA-BCP control chart detects 12 when there are shifts in the Y-intercept. Fig. 11 shows that the EWMA-C control chart detects nine out-of-control signals while the EWMA-BCP control chart detects 16 when there are shifts in slope. Fig. 12 indicates that the EWMA-C control chart detects 22 out-of-control signals while the EWMA-BCP control chart detects 25 when there are shifts in errors variance. This concludes that the better detection ability is offered by EWMA-BCP charts while monitoring the profile parameters of Y-intercepts, slopes and errors variance. The above finding indicates that the out-of-control sample points shown in-control by EWMA-C charts has been identified by EWMA-BCP charts which concludes that acquiring more knowledge about the process structures in the form of prior distributions are beneficial that enhance the detection potentials of the EWMA control charting structures.

5. Conclusions and recommendations This article introduced Bayesian approach by using noninformative (uniform) priors and informative (non-conjugate and conjugate) priors for linear profiles monitoring under phase II method. The results under non-informative (uniform) priors are identical to that of classical approach, while results are improved with the use of informative priors. The results of different individual and overall performance measures demonstrate that the EWMA-BCP control charts are more efficient than that of EWMABNCP and EWMA-C. It is noted that the performance of the EWMA-BCP control charts are better than the EWMA-BNCP control charts for monitoring the Y-intercept and slope, while it is almost

identical for monitoring the errors variance. The results under natural scale are slightly improved for monitoring increase in errors variance for both classical and Bayesian EWMA control charts than that of log scale. The change of scale does not change the performance results, while monitoring the Y-intercept and slope. The efficiency of the EWMA-BNCP and EWMA-BCP control charts is high for small shifts, while its efficacy decreases as the values of shifts increase. Referring to the results of ARL, MDRL, and SDRL, the values of EQL, RARL and PCI also affirms that EWMA-BCP control charts are more practical than the EWMA-BNCP and EWMA-C control charts. Overall, our proposed Bayesian control charting structures with non-conjugate and conjugate priors are superior to existing counterparts for monitoring the Y-intercepts, slopes, and error variances. From a future prospective, our proposed approach entails productive implications and vast potentials for applications. Further, the noted approach bears substantial depth and hence calls for extension. The author initiated an experimental design using simple linear regression model. The next part on this direction might be the use of non linear regression models, multiple regression models and polynomial regression models etc. Appendix A. Let us consider the Bramwell, Holdsworth, Pinton (BHP) distribution given in Eq. (10) above as prior knowledge of process parameters of the profile model given in Eq. (4) above. Now following Eq. (2) above, the joint posterior distribution for Y-intercept and slope are obtained as:

(

) n   1 X  2 y  B  B x 0 1 i 2r2 i¼1 ij      p B0  t0 B0  t0  exp  exp 2 k0 k0      p B1  t1 B1  t1  exp  exp 2 k1 k1

pðB0 ; B1 jdataÞ / exp

Expand the square term and rewrite the above equation as:

( pðB0 ;B1 jdataÞ / exp

1 2r2

nB20 2B0

n n n X X X yij þB21 x2 xi yij i 2B1 i¼1

i¼1

!)

i¼1

n    o n    o exp p2 B0k0t0 exp B0k0t0 exp p2 B1k1t1 exp B1k1t1 )

(

pðB0 ;B1 jdataÞ / exp ( exp

1 2r2

1 2r2

nB20 2B0

! ) n     X B0 t0 B0 t0 p yij þ 2 k0 exp k0 i¼1

! ) n n     X X 2 B1 t1 B1 t1 2  p B1 xi 2B1 xi yij þ 2 k1 exp k1 i¼1

i¼1

pðB0 ; B1 jdataÞ / pðB0 jdataÞ  pðB1 jdataÞ

ðA:1Þ

Where,

(

n X 1 pðB0 jdataÞ / exp nB20  2B0 yij 2 2r i¼1  ) B0  t0  exp k0

(

!



þ

p B0  t0 2

n n X X 1 pðB1 jdataÞ / exp B21 x2 xi yij i  2B1 2 2r i¼1 i¼1  ) B1  t1  exp k1



k0 ðA:2Þ !



þ

p B1  t1 2



k1 ðA:3Þ

266

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Let us first consider Eq. (A.2) above and simplify by using the second order Maclaurin series approximation. We know that the Taylor series can be represented as: 0

k30

exp

  k0

pðxÞ ¼ f ðxÞ þ

ðA:4Þ

i¼1



bonj ¼

00

f ð xÞ f ðxÞ ðx  xÞ þ ðx  xÞ2 1! 2! 000 f ðxÞ þ ð x  xÞ 3 þ . . . 3!

n X k0 yij þ p2 r2

t0

r2onj ¼ 

r2 k40 exp

nk40 exp

nk40 exp  

  t0 k0

! þr

þ r2

2



t0 

k20





! ;

and

t0

  t0 k0

k0

þ r2

: ðA:8Þ

Maclaurin series is a special case of Taylor series when x ¼ 0 i.e. 0

pðxÞ ¼ f ð0Þ þ

00

Solving Eq. (A.3) on the similar lines we obtained:

000

f ð0Þ f ð0Þ 2 f ð0Þ 3 ðxÞ þ ð xÞ þ ðxÞ þ    1! 2! 3!

Let pðxÞ ¼ exp



B0 t0 k0

ðA:5Þ



    t0 1 t0 0 ; f ð0Þ ¼ ; f ð0Þ ¼ exp  exp  k0 k0 k0     1 t0 1 t0 00 000 ; f ð0Þ ¼ 3 exp  ;... f ð0Þ ¼ 2 exp  k0 k0 k0 k0 Now, putting the values together in Eq. (A.5) of Maclaurin series:

     t0 1 t0 B0  t0 1 þ exp  þ 2 pðxÞ ¼ exp  k0 k0 k0 k0 2k0       t0 B0  t0 2 1 t0 B0  t0 3  exp  þ 3 exp  ;... k0 k0 k0 k0 6k0 Using the second order approximation



  P    n t1  p r2 þ r2 t  k2 x y þ k 1 1 ij 1 i¼1 i k1 2   P  b1nj ¼ ; n 2 2 k41 exp tk11 i¼1 xi þ r   2 4 r k1 exp tk11  P : r21nj ¼  4 n 2 2 k1 exp tk11 i¼1 xi þ r k31 exp

and

ðA:9Þ

Now consider the one parameter Levy distribution as prior knowledge for errors variance of the profile model. Again following Eq. (2) above, then the posterior distribution for errors variance are obtained as: ( ) n ng o   n    3 1 X  2 0 pðr2 jdataÞ / r2 2 exp y  B  B x :  r2 2 exp 0 1 i ij 2 2r i¼1 2r2     n3 1 g0 ðn  2Þr20 pðr2 jdataÞ / r2 2 2 exp þ : 2 r 2 2    1 a 1 b : pðr2 jdataÞ / r2 nj exp 2 nj

r

   (    2 ) B0  t0 t0 1 B0  t0 1 B0  t0 ¼ exp  1þ þ 2 exp k0 k0 k0 k0 k0 2k0 ðA:6Þ

ðA:10Þ where,

anj ¼

nþ1 g ðn  2Þr20 ; bnj ¼ 0 þ : 2 2 2

ðA:11Þ

Replace the approximation in Eq. (A.6) in Eq. (A.2).

"

( ! n X 1 p 2 2 pðB0 jdataÞ / exp r nB0  2B0 yij  2 B0 2 2k 0 i¼1      2 !)# t0 1 B0  t0 1 B0  t0 1þ þ 2 þ 2 exp  k0 k0 k0 k0 k0 Simplifying the above equation we get. 8 < 1    pðB0 jdataÞ / exp :2 nr2 þ k4 exp  t0 0 k0 !1 9 0 n   X  > pk1 2 t0  2 2 0 > r yij þ 2 þ k0 exp  k0 t0 k0  1 C> > B = B 2 C i¼1 C      2B B B 0 0 B C> 4 t 0 2 @ A> nr þ k0 exp  k0 > > ; Now, completing the square of above Equation to make it a proper normal distribution as:

pðB0 jdataÞ / exp

  1 2  ronj B0  bonj 2 2

Where,

r

)

n   X  pk1 t0  yij þ 20 þ k2 t0 k2 0 exp  k0 0 1 i¼1

   t0 nr2 þ k4 0 exp  k0    t0 ¼ nr2 þ k4 exp  0 k0

bonj ¼

r2 onj

2

ðA:7Þ !

Appendix B. Successive values of the proposed EWMA for monitoring a process variance at successive time t are obtained as,

   EWMAt ¼ max ð1  kÞEWMAt1 þ kyt ; ln r20 ; 

In sampling theory case, conditional on r2 , probability distribution for 1=s2 is scaled inverse chi-square distribution and probability for r2 conditional on s2 given scale agnostic prior is also scaled inverse chi-square distribution, the marginal posterior distribution of r2 is also scaled inverse chi-square distribution or equivalently s2t values have inverse gamma posterior distribution derived by using Eq. (12) i.e. IGðan ; bn Þ where the values of an and bn given   ðn2Þr2 v d2 ; 02 0 þ 2 0 considin Eq. (13). So it can be written as, IG v20 þ n2 2 ering the unbiased estimator. Here n is the sample size, and r20 underlying in-control process variance. Now, suppose that x random variable follows inverse gamma distribution i.e. s2t ¼ x with parameters an > 0 and bn > 0, then the pdf of x is given as,

ban f ðxÞ ¼ pnffiffiffiffiffi xan 1 eðbn =xÞ ; ;

an

and

ðB:1Þ

 where EWMA0 ¼ ln r , k the smoothing constant ð0 < k 6 1Þ: let  2 us assume that yt ¼ ln st , here s2t are successive sample variances. 2 0

x>0

ðB:2Þ

Let w ¼ lnðxÞ ) g 1 ðwÞ ¼ ew then random variable w has posterior distribution as log-inverse gamma distribution,



d 1

g ðwÞ

: f w ðwÞ ¼ f x ðg 1 ðwÞÞ

dw

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Now use inverse gamma distribution given in Eq. (B.2) we obtain, an

b w ) f w ðwÞ ¼ pnffiffiffiffiffi ean wbn e ; 1 < w < 1

an

ðB:3Þ

The above posterior probability distribution in Eq. (B.3) is loginverse gamma distribution that can be re-expressed as,

1 f w ðwÞ ¼ pffiffiffiffiffi exp fan ðln bn  wÞ  expðln bn  wÞg; 1 < w

an

<1

ðB:4Þ

The posterior estimate bn related with posterior inverse gamma distribution is scale parameter, while lnðbn Þ is location parameter of parallel posterior log-inverse gamma distribution. The design of EWMA chart principally is to follow location fluctuations in process, then using log transformation looks rational for monitoring and detecting possible changes in original process scale parameter. The shape parameter is known and determined by sample size n and prior information a0 the hyper parameter. We know that distribution of log of inverse gamma random variable is equivalent to the distribution of –log of gamma random variable i.e. EðXÞ ¼ EðXÞ and VarðXÞ ¼ VarðXÞ, then mean and variance of variable w ¼ lnðxÞ is EðwÞ ¼ wðan Þ  lnðbn Þ and VarðwÞ ¼ w0 ðan Þ respectively, Lawless (2002), where wðan Þ and w0 ðan Þ are digamma and trigamma functions. The variance only depends on posterior shape parameter an or equivalently the sample size n and prior parameter a0 . Using infinite series expansions, as given by Abramowitz and Stegun (1965).

wðxÞ ¼ lnðxÞ  w0 ðxÞ ¼

1 1 1  þ   2x 12x2 120x4

ðB:5Þ

1 1 1 1 þ þ  þ  x 2x2 6x3 30x5

ðB:6Þ

ðB:5Þ ) E½lnðs2 Þ 1 1   ¼ lnðr2 Þ   2 2 a0 þ n2 12 a0 þ n2 2 2

þ

1

 4     120 a0 þ n2 2

ðB:7Þ

ðB:6Þ ) Var½lnðs2 Þ ¼

1

a0 þ n2 2



1 1 þ   þ  3 n2 2 2 a0 þ 2 6 a0 þ n2 2

1  5    30 a0 þ n2 2

ðB:8Þ ) Var½lnðs2 Þ ¼

2 2 4 16 þ þ  ; Q 1 Q 2 3Q 3 15Q 4

ðB:8Þ

ðB:9Þ

where,

Q 1 ¼ 2a0 þ ðn  2Þ

ðiÞ

Q 2 ¼ 2a20 þ ðn  2Þ2 þ 4a0 ðn  2Þ

ðiiÞ

Q 3 ¼ 8a30 þ 12a20 ðn  2Þ þ 6a0 ðn  2Þ2 þ ðn  2Þ3

ðiiiÞ

Q 4 ¼ 32a50 þ 80a40 ðn  2Þ þ 80a30 ðn  2Þ2 þ 40a20 ðn  2Þ3 þ 10a0 ðn  2Þ4 þ ðn  2Þ5

ðivÞ

The formula in Eq. (B.9) can be used to approximate Var½lnðs2 Þ, which appears as factor in upper control limit when monitoring of error variance is under study.

267

References Abramowitz, M., & Stegun, I. A. (1965). Handbook of mathematical functions. New York, NY: Dover Publications. Ahmad, S., Abbasi, S. A., Riaz, M., & Abbas, N. (2014). On efficient use of auxiliary information for control charting in SPC. Computers and Industrial Engineering, 67 (2014), 173–184. Ali, S., & Riaz, M. (2014). Cumulative quantity control chart for the mixture of inverse rayleigh process. Computers and Industrial Engineering, 73(July 2014), 11–20. Amiri, A., Eyvazian, M., Zou, C., & Noorossana, R. (2012). A parameters reduction method for monitoring multiple linear regression profile. International Journal of Advanced Technology, 58(5–8), 621–629. Bedrick, E. J., Christensen, R., & Johnson, W. (1996). A new perspective on priors for generalized linear models. Journal of the American Statistical Association, 91, 1450–1460. Bramwell, S. T., Holdsworth, P. C. W., & Pinton, J. F. (1998). Universality of rare fluctuations in turbulence and critical phenomena. Nature, 396, 552–554. Brill, R. V. (2001). A case study for control charting a product quality measure that is a continuous function over time. In Presented at the 45th annual fall technical conference, Toronto, Ontario. Crowder, S. V., & Hamilton, M. D. (1992). An EWMA for monitoring a process standard deviation. Journal of Quality Technology, 24, 12–21. Durham, G., & Geweke, J. (2014). Improving asset price prediction when all models are false. Journal of Financial Econometrics, 12(2), 278–306. Elfadaly, F. G., & Garthwaite, P. H. (2014). Eliciting prior distributions for extra parameters in some generalized linear models. Statistical Modeling, XX(X), 1–21. Florens, J., & Simoni, A. (2012). Nonparametric estimation of an instrumental regression: A quasi-Bayesian approach based on regularized posterior. Econometrics, 170, 458–475. Garthwaite, P. H., & Dickey, J. M. (1985). Double- and single-bisection methods for subjective probability assessment in a location-scale family. Econometrics, 29 (1–2), 149–163. Garthwaite, P. H., & Dickey, J. M. (1988). Quantifying expert opinion in linear regression problems. Journal of the Royal Statistical Society, Series B, 50, 462–474. Garthwaite, P. H., & Dickey, J. M. (1992). Elicitation of prior distributions for variable selection problems in regression. Annals of Statistics, 20, 1697–1719. Garthwaite, P. H., Kadane, J. B., & O’Hagan, A. (2005). Statistical methods for eliciting probability distributions. Journal of the American Statistical Association, 100, 680–701. Garthwaitea, P. H., Al-Awadhib, S. A., Elfadalya, F. G., & Jenkinsonc, D. J. (2013). Prior distribution elicitation for generalized linear and piecewise-linear models. Journal of Applied Statistics, 40(1), 59–75. Gelman, A., Carlin, J., Stern, H., & Rubin, D. (2003). Bayesian data analysis.Boca Raton: Chapman & Hall CRC. Huber, G. P. (1974). Methods of quantifying subjective probabilities and multiattribute utilities. Decision Science, 5(3), 430–458. Jonathan, C. B. & Roger, D. C. (2014). Use of a levy distribution for modeling best case execution time variation. In Computer performance engineering. Lecture notes in computer science (vol. 8721 (September 2014), pp. 74–88). Kang, L., & Albin, S. L. (2000). On-line monitoring when the process yields a linear profile. Journal of Quality Technology, 32(4), 418–426. Keramatpoura, M., Niakib, A. T. S., & Amiric, A. (2014). Phase-II monitoring of Ar (1) auto-correlated polynomial profiles. Optimization in Industrial Engineering, 7(2), 53–59. Kim, K., Mahmoud, M. A., & Woodall, W. H. (2003). On the monitoring of linear profiles. Journal of Quality Technology, 35(3), 317–328. Lawless, J. F. (2002). Statistical models and methods for lifetime data (2nd ed.). 978-0471-37215-8New York, NY: John Wiley & Sons. Lucas, J. M., & Saccucci, M. S. (1990). Exponentially weighted moving average control schemes: properties and enhancements. Technometrics, 32(1), 1–12. Marcellus, L. R. (2007). Bayesian monitoring to detect a shift in process mean. Journal of Quality and Reliability Engineering International, 24, 303–313. Marlin, T. E. E. (2000). Process control: Controlling processes and control systems for dynamic performance (2nd ed.)New York: McGraw-Hill. Matthias, H. Y., Tana, Y., & Jianjun, S. (2012). A Bayesian approach for interpreting mean shifts in multivariate quality control. Technometrics, 54(3), 294–307. Montgomery, D. C. (2009). Introduction to statistical quality control (6th ed., pp. 476– 478). New York: John Wiley & sons. O’Hagan, A., Buck, C. E., Daneshkhah, A., Eiser, J. R., Garthwaite, P. H., Jenkinson, D. J., ... Rakow, T. (2006). Uncertain judgements: Eliciting expert probabilities. Chichester: Wiley. Ou, Y. J., Wu, Z., & Goh, T. N. (2011). A new sprt chart for monitoring process mean and variance. International Journal of Production Economics, 132(2), 303–314. Ou, Y. J., Wu, Z., & Tsung, P. (2012). A comparison study of effectiveness and robustness of a control charts for monitoring process mean. International Journal of Production Economics, 135(1), 479–490. Paciorek, J. C. (2006). Misinformation in the conjugate prior for the linear model with implications for free-knot spline modelling. Bayesian Analysis, 1(2), 375–383. Riaz, M., Abbasi, S. A., Ahmad, S., & Zaman, B. (2014). On efficient phase II process monitoring charts. The International Journal of Advanced Manufacturing Technology, 70(9), 2263–2274.

268

T. Abbas et al. / Computers & Industrial Engineering 94 (2016) 245–268

Roberts, S. W. (1966). A comparison of some control chart procedures. Technometrics, 8(3), 411–430. Ryan, T. P. (1997). Modern regression methods (1st ed., pp. 38–39). New York, NY: John Wiley & Sons. Ryu, J. H., Wan, H., & Kims, S. (2010). Otimal control of a cusum chart for a mean shift of unknown size. Journal of Quality Technology, 43(3), 311–326. Shi, X., Lv, Y., Fei, Z., & Liang, J. (2013). A multivariable statistical process monitoring method based on multiscale analysis and principal curves. International Journal of Innovative Computing, Information and Control, 9(4), 1781–1800. Wu, Z., Jiao, J., Yang, M., Liu, Y., & Wang, Z. (2009). An enhanced adaptive cusum control chart. IIE Transactions, 41(7), 642–653.

Zaman, B., Riaz, M., Ahmad, S., & Abbasi, S. A. (2015). On artificial neural networking based process monitoring under bootstrapping using runs rules schemes. The International Journal of Advanced Manufacturing Technology, 76(1), 311–327. Zellner, A., Ando, T., Basturk, N., Hoogerheide, L., & Van Dijk, H. (2014). Bayesian analysis of instrumental variable models: acceptance-rejection within Direct Monte Carlo. Econometric Reviews, 33, 3–35. Zhang, Y., He, Z., Zhang, C., & Woodall, W. H. (2013). Control charts for monitoring linear profiles with within profile correlation using gaussian process models. Quality and Reliability Engineering International, 30(4), 487–501. Zhu, J., & Lin, D. K. J. (2010). Monitoring the slopes of linear profiles. Quality Engineering, 22(1), 1–12.