Phase-II monitoring of the lower truncated proportional hazard rate process based on progressive Type-II censoring

TL;DR

This paper introduces new control charts for monitoring processes with censored data, useful in reliability engineering and lifetime experiments.

Contribution

Four new control charts are proposed for monitoring a lower truncated proportional hazard rate model under progressive Type-II censoring.

Findings

Three EWMA charts using likelihood ratio and maximum likelihood estimators were developed.

A novel weighted log-likelihood ratio chart was introduced and evaluated.

A simulation study and practical example demonstrated the effectiveness of the proposed methods.

Abstract

The problem of monitoring statistical processes using complete data has been extensively studied by researchers. However, in fields such as reliability engineering and lifetime experiments, complete samples are often not available. To address this gap, we introduce four control charts designed for monitoring both parameters of a family of distributions known as the lower truncated proportional hazard rate model, specifically under progressively Type-II censoring. Three of these control charts are exponentially weighted moving average (EWMA) charts that utilize the likelihood ratio statistic and maximum likelihood estimators. The fourth chart is based on a novel weighted log-likelihood ratio statistic. We conduct a Monte Carlo simulation study to assess the performance of the proposed control charts. Finally, we present a practical example to illustrate the application of our methods.

Figures

PDF and CDF of two populations based on the estimated parameters.

Monitoring of the lifetimes of mice after receiving a radiation dose of 300 roentgens at an age of 5-6 weeks.

Equations

Full text

Introduction

Evaluating the lifetime of a product is a critical aspect of product reliability. Research on product lifetime has been conducted across various fields, including industrial manufacturing, medicine, and reliability engineering.

To make inferences about a product’s lifetime under standard manufacturing conditions and to swiftly identify any defective items, it is essential to conduct life tests and reliability experiments. Statistical Process Monitoring (SPM) serves as a widely used tool for making inferences regarding process performance and has been implemented extensively in industrial settings. Within SPM, control charts play a pivotal role as graphical tools for monitoring process status and determining whether the process is in control (IC). The development of efficient control charts is crucial for the continuous monitoring of product lifetime, enabling real-time alerts when signs of deterioration are detected.

The EWMA control chart was first introduced by [1]. It was designed to detect small shifts that allows for quicker detection of changes in the process, making it a valuable tool in quality control and monitoring. Due to ability of this control chart, many researchers have extended to several models in the recent years. For instance, see [2–4].

In reliability engineering and life testing, we usually do not observe complete data and we come across the censored data to save time and cost involving the experiment. Some researchers focused on developing control charts based on censored data. For instance, [5] proposed an EWMA control chart to monitor the Weibull mean based on random censored data. [6] proposed a control chart to monitor the shape parameter of a Weibull process under Type-II censoring. [7] introduced an EMWA control chart for the Poisson exponential lifetime model under Type-I censored data. [8] proposed a cumulative sum (CUSUM) control chart to monitor the mean of the Weibull lifetime model with the fixed shape parameter. [9] discussed a weighted likelihood control chart based on the Weibull distributed observations under random censoring. [10] proposed a weighted EWMA chart for Weibull distribution in Type-I censored samples. [11] proposed some control charts for monitoring the lower Weibull percentiles under Type-II censoring. [12, 13] proposed some control charts under Type-II censored reliability tests by assuming that the lifetimes follow the Weibull distribution with fixed and stable shape parameter.

Unlike conventional censoring, progressive censoring offers more flexibility by permitting the removal of certain living subjects during the experiment. Progressive Type-II censoring has been introduced by [14, 15]. In recent years, a lot of researches has been made based on progressively Type-II censored samples; for comprehensive discussions in the theory of this topic, one may refer to [16, 17]. In progressive Type-II censoring, it is assumed that the removals of still operating units are carried out at observed failure times and that the censoring scheme $[eqn]$ is known in advance. Moreover, the number of units (n) and the number of observed failure times (m) are prefixed. Starting all n units at the same time, the first progressive censoring step takes place at the observation of the first failure time X1:m:n, at this time, r1 units are randomly chosen from the still operating units and withdrawn from the experiment. Then, the experiment continues with the reduced sample size $[eqn]$ . After observing the next failure at time X2:m:n, r2 units are randomly removed from $[eqn]$ active units. This process continues until the mth failure is observed. Then, the experiment ends, and $[eqn]$ are said to be progressively Type-II censored order statistics. This scheme includes as special cases the complete sample situation (when m = n and $[eqn]$ ) and the conventional Type-II right censoring situation (when $[eqn]$ and $[eqn]$ ). The conventional Type-II right censoring scheme arises in a life-testing experiment whenever the experimenter can not observe the objects with the high lifetimes and one is forced to only record the first m smallest lifetimes of n objects for inference. But in the setting of the progressive Type-II censoring, the experimenter may choose different censoring schemes $[eqn]$ with fixed m and n to get the different inferences. This property of progressive Type-II censoring allows the experimenter to choose a censoring scheme to get the optimal inference. This subject provides valuable design ideas to reliability practitioners and evidence for the usefulness and efficiency of progressive Type-II right censoring as compared to conventional Type-II right censoring. For more details, we refer to [16].

Recently, developing control charts based on progressively Type-II censored samples has been considered by some researchers. For instance, [18] proposed a control chart to monitor the shape parameter of the Weibull distributed process under progressive Type-II censoring. [19] introduced some control charts for two-parameter exponential distribution for monitoring the location and scale parameters, simultaneously. [20] proposed some control charts for monitoring the scale parameter and also joint shape and scale parameters in the Weibull distributed process.

Modern products typically exhibit high reliability, and models that include small values in their domain often fail to accurately describe the lifetimes of such products. To address this issue, [21] introduced the lower truncated proportional hazard rate (LTPHR) model. The proportional hazard rate (PHR) model is a significant framework in reliability theory and various other fields. A random variable is said to follow the PHR model, if its survival function can be expressed as $[eqn]$ where $[eqn]$ is a baseline absolutely continuous survival function and $[eqn]$ . For a more comprehensive discussion about the PHR model, we refer to [22–24]. A prevalent way to describe a random variable having a certain lower endpoint is a lower truncation of the well-known distributions. Let $[eqn]$ , $[eqn]$ be the survival function of a random variable X. The lower truncated version of X at point $[eqn]$ is defined as $[eqn]$ with the corresponding survival function $[eqn]$ , $[eqn]$ . By this idea, a random variable X is said to follow the lower truncated proportional hazard rate model if its survival function can be expressed as $[eqn]$ , $[eqn]$ , $[eqn]$ , where $[eqn]$ is a baseline absolutely continuous survival function with the corresponding density function f(x). In this setting, we denote $[eqn]$ . The LTPHR model has been recently considered in some research. For instance, see [25, 26].

In this paper, we focus on designing control charts for Phase-II joint monitoring of the parameters $[eqn]$ and $[eqn]$ in a LTPHR process, utilizing progressively Type-II censored data. Sect 1 introduces four distinct control charts tailored for this purpose. We discuss about the implementation design in Sect 2. A comprehensive simulation study is presented in Sect 3, evaluating the performance of our proposed control charts under various scenarios. Sect 4 demonstrates the application of our methods through an analysis of a real dataset. In Sect 5 we make a conclusion and the R codes used in simulation study are provided in Sect 6.

1 Control charts

In this section, we present control charts designed to identify potential shifts in both parameters of the lower truncated proportional hazard rate model under a progressive Type-II censoring scheme.

In Phase II, the batch of progressive Type-II censored samples $[eqn]$ from a LTPHR model are collected sequentially over time with the same censoring scheme $[eqn]$ . Let $[eqn]$ be the IC value of the parameter vector $[eqn]$ .

1.1 EWMA control chart based on the likelihood ratio

In this subsection, we provide an EWMA control chart based on the likelihood ratio (EWMA-LR) to monitor the LTPHR process of progressively Type-II censored data.

Let $[eqn]$ be the ith batch of progressively Type-II censored sample from $[eqn]$ with the observation $[eqn]$ . According to [16] the likelihood function corresponding to $[eqn]$ can be written as

$$[eqn]$$

with the corresponding log-likelihood

$$[eqn]$$

where c is a free of $[eqn]$ , $[eqn]$ , and $[eqn]$ denotes the indicator function. The MLEs of $[eqn]$ and $[eqn]$ based on the ith batch sample, denoted by $[eqn]$ and $[eqn]$ respectively, are obtained by maximizing the function Li. Based on Theorem 1 of [21], MLEs of the parameters are obtained as the following:

$$[eqn]$$

Let $[eqn]$ be the IC value of $[eqn]$ . The monitoring statistic, denoted by LRi, based on the likelihood ratio, is given by

$$[eqn]$$

Therefore, the charting statistic of the EWMA-LR control chart is given by

$$[eqn]$$

where $[eqn]$ is a smoothing parameter and $[eqn]$ . The smoothing parameter $[eqn]$ plays a crucial role in the EWMA control chart, as it determines the weight assigned to recent observations. This parameter directly influences the smoothness and sensitivity of the control chart in response to process changes. A higher value of $[eqn]$ gives more weight to recent observations, making the control chart more responsive to immediate changes in the process. While this heightened sensitivity can be beneficial for quick detection of shifts, it also increases the rate of false alarms. Conversely, a smaller value of $[eqn]$ places greater emphasis on past observations, resulting in a smoother control chart. This approach reduces sensitivity to immediate changes, which may decrease the rate of false alarms but can also delay the detection of small shifts in the process. In summary, selecting an appropriate value for $[eqn]$ is essential for balancing sensitivity and false alarm rates in EWMA control charts. In most practical applications, one may consider $[eqn]$ .

Note that EL0 does not have a closed-form expression and can be computed using Monte Carlo simulations. The EWMA-LR control chart generates a signal when ELi>h , where the threshold value h is selected to achieve a specified expected IC average run length ( $[eqn]$ ).

Based on Theorem 2.1 of [25], $[eqn]$ follows the chi-square distribution with 2m–2 degrees of freedom and $[eqn]$ follows the beta distribution with the parameter vector $[eqn]$ . Thus, for an IC process, the distribution of LRi only depends on n and m and is free of the censoring scheme $[eqn]$ . Hence, the EWMA-LR control chart for the LTPHR model under progressive Type-II censoring does not depend on the censoring scheme.

1.2 EWMA-Max-MLE control chart

Here we utilize the idea of [27] to construct EWMA-Max-MLE charting statistic. Based on Theorem 2.1 of [25], we have

$$[eqn]$$

where $[eqn]$ denotes the beta distribution with the parameter vector $[eqn]$ and $[eqn]$ denotes the chi-squared distribution with 2m–2 degrees of freedom.

By the probability integral transform, $[eqn]$ and $[eqn]$ independently follow the standard normal distribution when the process is IC, where $[eqn]$ denotes the cumulative distribution function of the beta distribution with parameter $[eqn]$ , $[eqn]$ denotes the cumulative distribution function of the chi-squared distribution with 2m–2 degrees of freedom, and $[eqn]$ denotes the quantile function of the standard normal distribution.

Define $[eqn]$ . An EWMA-Max-MLE control charting statistic can be obtained as

$$[eqn]$$

where $[eqn]$ such that according to [19], we have EM0 = 1.128379. The EWMA-Max-MLE control chart generates signals when ELi>h, where the value of h is selected to attain a specified expected $[eqn]$ .

It is clear that the distribution of Si only depends on n and m and is free of the censoring scheme $[eqn]$ . Hence, to construct the EWMA-Max-MLE control chart for the LTPHR model under progressive Type-II censoring, the censoring scheme is not required.

1.3 EWMA control chart based on the Manhattan distance

In this subsection, we provide an EWMA control chart based on the Manhattan distance (EWMA-MD) to monitor the LTPHR model process of progressively Type-II censored data. For this purpose, we use the Manhattan distance between the points $[eqn]$ and $[eqn]$ , of the form $[eqn]$ , in ith batch sample. The charting statistic of the EWMA-MD control chart is given by

$$[eqn]$$

where $[eqn]$ is a smoothing parameter and $[eqn]$ . The EWMA-MD control chart generates signals when EMDi>h, where the value of h is selected to get a specified expected $[eqn]$ . Also, by a similar argument to the previous subsections, we conclude that the EWMA-MD control chart for the LTPHR model under progressive Type-II censoring does not depend on the censoring scheme and only depends on m and n.

1.4 WL control chart

Hypothesis testing based on the weighted log-likelihood was proposed by [28]. Then, [29] used this idea to design the control charts for monitoring the Poisson count process. Recently, [20] proposed a WL control chart for monitoring a Weibull process under Type-II censoring by the weighted log-likelihood function.

We present a WL chart statistic to achieve sequential monitoring. From time 1 to t, by Eq (1), sum of the weighted log-likelihood functions can be written as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the weight of log-likelihood function corresponding to ith batch sample and $[eqn]$ is a smoothing parameter. Similar to MLEs, the weighted maximum log-likelihood estimators (WMLE) can be obtained by maximizing $[eqn]$ . It is clear that $[eqn]$ is an increasing function with respect to $[eqn]$ for $[eqn]$ . Thus, the WMLE of $[eqn]$ , denoted by $[eqn]$ , is obtained as $[eqn]$ , where $[eqn]$ ( $[eqn]$ ) is the MLE of $[eqn]$ based on the ith batch sample. Also, WMLE of $[eqn]$ , denoted by $[eqn]$ is obtained by the equation

$$[eqn]$$

Solving the above equation immediately implies that

$$[eqn]$$

Then, the WL charting statistic can be obtained as

$$[eqn]$$ $$[eqn]$$

The WL control chart generates signals when WLt>h, where the value of h is selected to attain a specified expected $[eqn]$ .

By a simple algebraic calculation, $[eqn]$ can be rewritten as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ were defined in (2). Therefore, $[eqn]$ is a function of $[eqn]$ . On the other hand, based on the discussions in Subsect 1.1, the distribution of $[eqn]$ does not depend on $[eqn]$ which implies that the distribution of $[eqn]$ is free of the censoring scheme $[eqn]$ . Hence, the WL control chart for the LTPHR model under progressive Type-II censoring is free of the censoring scheme and only depends on m and n.

2 Monte Carlo algorithm for implementation of the proposed control charts

As discussed in our proposed control charts, it is essential to determine the control limit h in order to achieve a desired expected value $[eqn]$ . To facilitate this, we propose an algorithm specifically for designing the EWMA-LR control chart. Similar algorithms can be applied to other control charts as well. It is necessary to mention that, after obtaining the suitable control limit by Algorithm 2.1 and fixing h as a desirable value, we can compute the OOC ARL for the given shifted parameters $[eqn]$ and $[eqn]$ using Steps (i)-(vii) of Algorithm 2.1. For this purpose, it is enough to replace $[eqn]$ and $[eqn]$ by the respective values $[eqn]$ and $[eqn]$ in Steps (ii) and (iii).

Algorithm 2.1. Given the life test sample size n, the number of observations m, the baseline survival function $[eqn]$ , the IC values parameters $[eqn]$ and $[eqn]$ , the smoothing parameter $[eqn]$ , the nominal $[eqn]$ , the number of replications of the Monte Carlo simulations M, and a starting value h,

(i) Set i = 1.(ii) Randomly generate B1 and B2 such that $[eqn]$ and $[eqn]$ .(iii) Compute $[eqn]$ and $[eqn]$ .(iv) Calculate ELi by Eq (4).(v) If ELi<h, set i = i + 1 and repeat Steps (ii)-(iv). And if ELi>h, the process is declared to be out of control (OOC) at the ith sample for the first time and set RL = i.(vi) Repeat steps (i)-(v) for M times to obtain $[eqn]$ .(vii) Compute $[eqn]$ . Note that SDRL is also defined as $[eqn]$ .(viii) Compare the obtained value of ARL in Step (vii) to the nominal $[eqn]$ . If the difference between these two values is negligible, record the starting control limit h.(ix) When the obtained ARL is greater (smaller) than $[eqn]$ , decrease (increase) the control limit h to an appropriate value and repeat Step (i)-(viii) until ARL becomes fairly closer to the target nominal value of $[eqn]$ . When we get the control limit under a set of parameters, it can be used to plot the control chart of the phase-II process.

3 Simulation study

In this section, we employ a Monte Carlo simulation to assess the performance of the proposed control charts presented in this paper. All computations are carried out using the $[eqn]$ software.

We mainly evaluate the performance of the proposed control charts based on the ARL and SDRL. We first determined the control limits for the specified values: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , along with various values of n and m. The simulation results are summarized in Table 1. All results presented in this paper are derived from M = 20000 replications.

Table 1: Means and control limits of four proposed control charts for F¯(μ0)=0.8, λ0=1, ARL0=370 and different values of n, m, and γ.

We use the zero-state ARL. The zero-state ARL refers to the expected number of samples collected before a control chart signals an OOC condition, assuming that any process shift occurs at the start of the monitoring period. At a specific time $[eqn]$ , let $[eqn]$ represent the expected number of samples until an OOC signal is triggered, given that no false alarms have occurred prior to this point. In this scenario, a process shift is assumed to happen at time $[eqn]$ . It is important to note that the zero-state ARL is equivalent to CED(1), while the steady-state ARL is defined as the limit of $[eqn]$ as $[eqn]$ tends to infinity, provided that such a limit exists. For more details, we refer to [30–32].

To evaluate the performance of the proposed control charts in terms of OOC ARL and SDRL, we must consider certain shifted parameters. Specifically, when the process goes OOC, the parameters shift from the IC parameters $[eqn]$ and $[eqn]$ to the OOC parameters $[eqn]$ and $[eqn]$ . Therefore, in view of the proposed charting statistics in this paper, we consider shift sizes $[eqn]$ and $[eqn]$ as

$$[eqn]$$

It is clear that $[eqn]$ indicates the IC process. We obtained the ARL and SDRL of the control charts for $[eqn]$ and some values of $[eqn]$ and $[eqn]$ . The results are reported in Tables 2, 3, 4, 5, 6 and 7.

Table 2: ARL and SDRL (in parentheses) values of the proposed control charts for n = 5, m = 3 and different values of γ, δ1 and δ2.

Table 3: ARL and SDRL (in parentheses) values of the proposed control charts for n = 5, m = 4 and different values of γ, δ1 and δ2.

Table 4: ARL and SDRL (in parentheses) values of the proposed control charts for n = 10, m = 4 and different values of γ, δ1 and δ2.

Table 5: ARL and SDRL (in parentheses) values of the proposed control charts for n = 10, m = 6 and different values of γ, δ1 and δ2.

Table 6: ARL and SDRL (in parentheses) values of the proposed control charts for n = 15, m = 6 and different values of γ, δ1 and δ2.

Table 7: ARL and SDRL (in parentheses) values of the proposed control charts for n = 15, m = 12 and different values of γ, δ1 and δ2.

The summary of results derived from Tables 2, 3, 4, 5, 6 and 7 is listed below:

(i) When there is no shift in $[eqn]$ , the simulated $[eqn]$ for all control charts is approximately equal to the nominal value 370. This indicates the stability of all the control charts proposed in this paper. Additionally, the SDRL for the EWMA-LR, EWMA-Max-MLE, and EWMA-MD control charts is close to $[eqn]$ , suggesting that the geometric distribution effectively approximates the IC run length distribution of these charts.

(ii) In certain scenarios, the EWMA-Max-MLE and EWMA-MD control charts exhibit ARL bias; specifically, the OOC ARL can exceed the nominal $[eqn]$ . This outcome is less than ideal for practical applications.

(iii) In most instances, the EWMA-LR control chart outperforms both the EWMA-Max-MLE and EWMA-MD control charts with respect to the ARL and SDRL criteria.

(iv) The WL control chart provides the smallest OOC ARL and SDRL across all cases.

Overall, we recommend using the WL control chart for monitoring the LTPHR process, regardless of $[eqn]$ , n and m.

4 Data analysis

We utilize the data sets presented in [33], which were gathered from a laboratory experiment involving two groups of male RFM strain mice. These mice were exposed to a radiation dose of 300 roentgens when they were 5 to 6 weeks old. The first group was housed in a conventional laboratory environment, while the second group was maintained in a germ-free environment. In his study, [33] identified two primary causes of death: thymic lymphoma and reticulum cell sarcoma, while categorizing all other causes into a single group. For our analysis, we focus on these two main groups: one in a conventional laboratory setting and the other in a germ-free environment.

[25] used these lifetime data by dividing by 1100 and got the data represented in Tables 8 and 9. They showed that the LTPHR models with the baseline survival function

Table 8: Lifetimes of mice lived in a conventional laboratory environment after receiving a radiation dose of 300 roentgens at an age of 5-6 weeks (divided by 1100).

Table 9: Lifetimes of mice lived in a germ-free environment after receiving a radiation dose of 300 roentgens at an age of 5-6 weeks (divided by 1100).

$$[eqn]$$

fit well to the data, such that

$$[eqn]$$

is called an incomplete gamma function. The baseline survival function is called a gamma-uniform distribution introduced by [34]. For the data set in Table 8, the lifetime of mice lived in a conventional laboratory environment divided by 1100, they estimated $[eqn]$ ( $[eqn]$ ) and $[eqn]$ , and for the data set in Table 9, the lifetime of mice lived in a germ-free environment divided by 1100, they estimated $[eqn]$ ( $[eqn]$ ) and $[eqn]$ . Fig 1 presents the probability density function (PDF) and the cumulative distribution function (CDF) of two populations based on the estimated parameters.

PDF and CDF of two populations based on the estimated parameters.

Using such information, we reconstruct realistic data sets adopting a Monte Carlo simulation to illustrate how to use the proposed four control charts. We observe 15 values for $[eqn]$ as the estimated parameters for a set of 15 samples with n = 5 and m = 3 (regardless of the censoring scheme $[eqn]$ ) with $[eqn]$ as data from an apparently IC process, and another 15 observations for $[eqn]$ as the estimated parameters for a set of 15 samples with n = 5 and m = 3 with $[eqn]$ ( $[eqn]$ ) as data from an apparently OOC process. The observed values are reported in Table 10. Based on these values and for $[eqn]$ , the control limits for the EWMA-LR, EWMA-Max-MLE, EWMA-MD, and WL control charts were obtained as 2.6606, 1.3444, 0.8248, and 0.1141, respectively. We display the corresponding control charts in Fig 2. We observe that the OOC signals for the EWMA-Max-MLE monitoring approach begin from the 22nd sample, while the EWMA-MD approach shows OOC signals starting from the 20th sample. In contrast, the WL monitoring method indicates OOC signals from the 18th sample onwards. The EWMA-LR monitoring approach, however, suggests that the process remains IC before the 25th sample, generating OOC signals only between the 25th and 27th observations. Notably, since the actual process shift begins from the 11th subgroup sample, this indicates that the WL control chart is capable of detecting shifts more effectively and earlier than the other monitoring charts.

Monitoring of the lifetimes of mice after receiving a radiation dose of 300 roentgens at an age of 5-6 weeks.

Table 10: Estimated parameters for lifetimes of mice after receiving a radiation dose of 300 roentgens at an age of 5-6 weeks for n = 5, m = 3.

5 Conclusion

In this paper, we introduced four control charts designed for monitoring both parameters of a LTPHR process characterized by the baseline survival function $[eqn]$ and the parameter vector $[eqn]$ . These charts were developed based on progressively Type-II censored data with a censoring scheme denoted as $[eqn]$ . The control charts include: an EWMA chart based on the likelihood ratio statistic (EWMA-LR), an EWMA chart utilizing maximum likelihood estimators (EWMA-Max-MLE), an EWMA chart based on the Manhattan distance (EWMA-MD), and a chart based on a weighted log-likelihood ratio statistic (WL).

To facilitate the implementation of the proposed control charts, we present a Monte Carlo algorithm that allows for the specification of the control chart for a given expected value $[eqn]$ .

We conducted a comprehensive Monte Carlo simulation study to evaluate the performance of the proposed control charts. This study involved comparing the OOC ARL and the SDRL across various scenarios. The simulation results demonstrate that the WL control chart consistently outperforms the others for monitoring the LTPHR process, regardless of the baseline survival function $[eqn]$ , number of units n, and the number of the observed failure times m.

Additionally, we applied our methods to analyze the lifetimes of two groups of male RFM strain mice, further illustrating their practical utility. In this example, the WL chart again showed superior performance.

A key advantage of our proposed control charts is their reliance solely on $[eqn]$ , the IC values parameters $[eqn]$ and $[eqn]$ , the smoothing parameter $[eqn]$ , and the nominal $[eqn]$ . Consequently, they do not depend on the specific censoring scheme employed.

The extension of the proposed control charts to self-starting schemes, along with the development of other censoring mechanisms, such as progressive hybrid censoring and adaptive progressive censoring, may be the focus of future researches.

6 Appendix

6.1 The R code of EWMA-LR chart

####### Computimg EL0

n=5; m=3     # insert n and m

d1=rexp(10000000)

d2=rchisq(rep,2m-2)/(2m)

LR=m*(d2-log(d2)-1)+d1

mean(LR)

####### Computimg ARL0 (IC) or ARL1 (OOC)

n=5; m=3           # insert n and m

EL0=2.0271         # insert EL0

h.EWMA.LR=2.6606   # insert the obtained h based on ARL0

delta1=-.2; delta2=.8  # set dalta1=0 and

dalta2=1 (IC) to obtain ARL0

gamma=.05          # insert the smoothing parameter gamma

Fbmu0=.8; lambda0=1 # inser Fbar(mu0) (IC) and lambda0 (IC)

rep=20000          # number of replication

Fbmu1=Fbmu0*exp(-delta1/lambda0)

lambda1=lambda0/delta2

L=c()

for(k in 1:rep){

i=0

EL=EL0

while(EL<h.EWMA.LR){

Fbmu1hat=Fbmu1exp(-rexp(1,nlambda1))

lambda1hat=2mlambda1/rchisq(1,2*m-2)

d1=nlambda0log(Fbmu0/Fbmu1hat)

d2=lambda0/lambda1hat

if(d1>0) LR=m*(d2-log(d2)-1)+d1 else LR=10000

EL=(1-gamma)EL+gammaLR

i=i+1

}

L[k]=i

}

mean(L)

sd(L)

6.2 The R code of EWMA-Max-MLE chart

#====== Computimg ARL0 (IC) or ARL1 (OOC)

n=5; m=3

EM0=1.128379

h.EWMA.Max.MLE=1.3444

delta1=-.2; delta2=.8

gamma=.05

Fbmu0=.8; lambda0=1

rep=20000

Fbmu1=Fbmu0*exp(-delta1/lambda0)

lambda1=lambda0/delta2

L=c()

for(k in 1:rep){

i=0

EM=EM0

while(EM<h.EWMA.Max.MLE){

Fbmu1hat=Fbmu1exp(-rexp(1,nlambda1))

lambda1hat=2mlambda1/rchisq(1,2*m-2)

Z1=Fbmu1hat/Fbmu0

Z2=2mlambda0/lambda1hat

T1=qnorm(pbeta(Z1,n*lambda0,1))

T2=qnorm(pchisq(Z2,2*m-2))

S=max(abs(T1),abs(T2))

EM=(1-gamma)EM+gammaS

i=i+1

}

L[k]=i

}

mean(L)

sd(L)

6.3 The R code of EWMA-MD chart

#====== Computimg EMD0

n=5; m=3 # insert n and m

lambda0=1 # insert lambda0

d1=log(rbeta(10000000,n**lambda0,1))

d2=rchisq(10000000,2m-2)/(2m)

MD=lambda0*abs(d1)+abs(d2-1)

mean(MD)

#====== Computimg ARL0 (IC) or ARL1 (OOC)

n=5; m=3

EM0=.6992

h.EWMA.MD=.8249

delta1=-.2; delta2=.8

gamma=.05

Fbmu0=.8; lambda0=1

rep=20000

Fbmu1=Fbmu0*exp(-delta1/lambda0)

lambda1=lambda0/delta2

L=c()

for(k in 1:rep){

i=0

EM=EM0

while(EM<h.EWMA.MD){

Fbmu1hat=Fbmu1exp(-rexp(1,nlambda1))

lambda1hat=2mlambda1/rchisq(1,2*m-2)

S=n^2lambda0^2(log(Fbmu0)-log(Fbmu1hat))^2+

m^2/(m-1)*(lambda0/lambda1hat-1)^2-2-m/(m-1)

EM=(1-gamma)EM+gammaS

i=i+1

}

L[k]=i

}

mean(L)

sd(L)

6.4 The R code of WL chart

#====== Computimg ARL0 (IC) or ARL1 (OOC)

n=5; m=3

h.WL=.1141

delta1=-.2; delta2=.8

gamma=.05

Fbmu0=.8; lambda0=1

rep=20000

Fbmu1=Fbmu0*exp(-delta1/lambda0)

lambda1=lambda0/delta2

L=c()

for(k in 1:rep){

K=matrix(,nr=100000,nc=2)

i=1

repeat{

Fbmu1hat=Fbmu1exp(-rexp(1,nlambda1))

lambda1hat=2mlambda1/rchisq(1,2*m-2)

K[i,]=c(Fbmu1hat,lambda1hat)

w=gamma*(1-gamma)^((i-1):0)

W=sum(w)

Fbmu1tilda=max(K[1:i,1])

lambda1tilda=mW/sum(w(m/K[1:i,2]-n*log(K[1:i,1]/Fbmu1tilda)))

d1=nlambda0log(Fbmu0/Fbmu1tilda)

d2=lambda0/lambda1tilda

if(d1>0) WL=W*(m*(d2-log(d2)-1)+d1) else WL=10000

if(WL > h.WL) {AR=i; break} else i=i+1

}

L[k]=AR

}

mean(L)

sd(L)