Bayesian Approach for Constant-Stress Accelerated Life Testing for Kumaraswamy Weibull Distribution with Censoring

The accelerated life tests provide quick information on the life time distributions by testing materials or products at higher than basic conditional levels of stress such as pressure, high temperature, vibration, voltage or load to induce failures. In this paper, the acceleration model assumed is log linear model. Constant stress tests are discussed based on Type I and Type II censoring. The Kumaraswmay Weibull distribution is used. The estimators of the parameters, reliability, hazard rate functions and p-th percentile at normal condition, low stress, and high stress are obtained. In addition, credible intervals for parameters of the models are constructed. Optimum test plan are designed. Some numerical studies are used to solve the complicated integrals such as Laplace and Markov Chain Monte Carlo methods.


Introduction
Technology advancements are continuously increasing the improvements in manufacturing designs. Therefore, it becomes more difficult to obtain information about lifetime of materials or products or with high reliability at the testing time of that device under normal conditions. Accelerated life tests (ALTs) are often used in such problems in order to shorten the life of test items. This test provides information quickly on the life distribution of the products or materials by testing them at higher than normal levels of stress such as voltage, high temperature, pressure, load or vibration to stimulate early failures. Such testing could save much time, manpower, material and capital. The stress can be applied in different ways: The commonly used methods are constant stress, step stress, and progressive stress (See Nelson (1990) and Bagdonavicius and Nikulin (2002)).
In Bayesian analysis the unknown parameters are considered as random variables. The use of prior information for the unknown parameters concerning engineering facts and material properties is essential. So, this prior information about the characteristics of Weibull distribution is one of the most popular models; it has been extensively used for modeling data in reliability, engineering and biological studies (Murthy et al. (2004)). A distribution was constructed by Kumaraswamy (1980) with two shape parameters (0, 1). Kumaraswamy distribution (Kum distribution) is applicable to many natural phenomena with the outcomes that have upper and lower bounds such as scores obtained in a test, heights of individuals, hydrological data and atmospheric temperatures. A composition between Kum distribution and any distribution can be constructed (See Cordeiro et al.
where are the shape parameters and is a scale parameter. These shape parameters allow a high degree of flexibility of the KumW distribution. It attracts wider applications in engineering, reliability and in other research areas.
The reliability function (rf) of KumW and the hazard rate function (hrf) corresponding to (2) can be written, respectively, as follows: This paper is organized as follows: The Laplace approximation is applied for the estimation of the parameters based on Type I and Type II censored samples in Section 2.
In Section 3, the MCMC method is carried out under censored samples of Type I and Type II.

Laplace Approximation
Laplace approximation was presented to solve complicated integrations. The conditions for Laplace regularity require that the integrals must exist and be finite. Also, the determinant of the Hessians must be bounded away from zero at the optimizers and that the log likelihood must be differentiable (from first to sixth order) on the parameters and all the partial derivatives must be bounded in the neighborhood of the optimizers. The expectation value must be a non-negative function (See Achcar and Rosales (1992)).
Bayesian estimation under Type I censored samples in Subsection (2.1) is derived. In Subsection (2.2) Bayesian estimation under Type II censored samples is obtained. The numerical results are illustrated in Subsection (2.3).

Bayesian Estimation based on Type I censored samples
Suppose that there are k levels of high stress , j=1, 2, ..., k and assume that is the usual condition where < < < ...< , and there are units on test at .When Type I censoring is applied at each stress level, the experiment terminates once all the items fail or when a fixed censoring time is reached. The lifetime at stress level ; , i=1, 2, ..., , j=1, 2, ..., k; is considered to have KumW distribution with the pdf in (2) The stress is assumed to affect only on the shape parameter of the KumW distribution through the log linear model as follows: The corresponding likelihood function is given by: where is the vector of parameters; ( ) , is an indicator variable such that: and Suppose that is unknown. The joint non informative prior distribution of these parameters is as follows: The joint posterior density of can be obtained as follows: where , ( ) is the joint prior distribution of , ∫ ∫ ∫ ∫ ∫ ∫ .
By taking the natural logarithm for (6), then substituting (6) in (9) as follows: one obtains Numerical integration is required to evaluate (12). The Laplace method applied by Achcar (1995), can be used to obtain in as follows: where is the number of the parameters, ̂ ̂ ̂ ̂ ̂ ̂ are the argmax( ); maximize ( ), ( ) is the Hessian matrix of evaluated at ̂ and =5 . and where , , , and .
Then, the joint posterior distribution in (11) for can be written as: where ̃ is given in (13) The marginal posterior densities and the conditional expectations for are obtained using the Laplace method to approximate integrals. Then, the estimators of the parameters based on the squared error loss function are the mean of the posterior distribution and can be obtained as follows: The estimator of is given by and where ̂ ̂ ̂ ̂ are the maximum likelihood estimators (MLEs). Applying Laplace method and using steps analogous to those used for estimating the parameter . One can obtain the Bayes estimators for The Bayes estimator of the rf, ( ) is given by where ( ) and ( | ) are given in (3) and (16) respectively.
The Bayes estimator of the hrf, ( ), is obtained as follows where ( ) is given by (4).

Bayesian Estimation based on Type II censoring
When Type II censoring is applied at each level of stress, once the number of failures out of units are reached the experiment terminates. The lifetime at level of stress , , i=1, 2, ..., , j=1, 2, ..., k, is considered to have KumW distribution with the pdf in (2). It is assumed that the stress effects on the shape parameter of the KumW distribution only through the log linear model as in (5).
The likelihood function of the j-th sample under Type II censoring is presented below: Thus, the likelihood function of the KumW distribution , j=1, 2..., k and i= 1, 2, ..., ; which are identically distributed and independent random variables, can be rewritten as  follows: where ( ) , ( ( ) )-( ) By taking the natural logarithm for (22), then assuming the same joint non informative prior distribution given by (8), the joint posterior density of and can be obtained similarly as given in (9), just by replacing with in (11) .
Applying the Laplace approximation and using steps analogous to those used in (12)- (21), also , ̃, ( ̂ ), the Bayes estimators based on Type II censoring considering the non-informative prior and squared error loss function can be obtained.
The Bayes estimator of a is given by In parallel with the steps used to obtain the Bayes estimator for a the Bayes estimators for ( ) ( ) can be obtained.

Numerical results
This subsection aims to illustrate the theoretical results of estimation problems on the basis of simulated data.

Simulation algorithm
 Several data sets are generated from KumW distribution for a combination of the population parameter and for sample sizes 20, 60 and 100. The transformation between uniform distribution and KumW distribution in step j is  It is assumed that there are only two different levels of stress (k=2), and , which are higher than the stress at usual condition, .


The pre-specified censoring times are and in Type I censored samples, while in Type II censored samples, numbers of test units are allocated to each level of stress , j=1, 2 where ,=0.9( ), j=1, 2 using 1000 replications for each sample size.


Computer program is used depending on MathCad 14 for obtaining the MLEs based on Type I and Type II censored samples. Then, the MLEs are used to solve Laplace approximation in (13) to obtain the Bayes estimates for the model parameters, rf and hrf.


Evaluating the performance of the estimators of , and is considered through some measurements of accuracy. In order to study the precision and variation of estimators, it is convenient to use the relative absolute bias (RAB 1 ) = | | , the mean square error (ER 1 ) and the relative error (RE 1 ) = √ ( ) under Type I censoring and the RAB 2 , ER 2 and RE 2 that are obtained based on Type II censored samples.


The credible intervals are obtained for the parameters using two sided approximate 100(1 credible intervals. The results based on Type I censoring are displayed in Tables 1-3 and Tables 4-6 indicate Type II censoring results.


The reliability function and hazard rate function are estimated for different values of mission times.

Concluding remarks
 It is clear from Tables 1 and 4 that the estimates (E 1 ) and (E 2 ) are close to the population parameters as the sample size increases. Also, as shown in the )%   numerical results the RAB 1 , RAB 2 , ER 1 , ER 2 , RE 1 and RE 2 decrease when the sample size increase. In general, ER 1 is better than ER 2 for all sample sizes. The estimate gives the best results over all the parameters in Type I and Type II censoring.


The two-sided 95% credible intervals for the parameters of KumW are displayed in Tables 2 and 5. These tables contain the standard error (SE), lower bound (L), upper bound (U) and the length (l) of intervals. The interval of the estimates becomes narrower as the sample size increases. The interval length of parameter is the shortest among the other parameters in Type I and Type II censoring. The length of is the shortest among the constant parameters of the acceleration model.
 Tables 3 and 6 indicate that the estimates of rf, hrf and their credible intervals. It is noted that the reliability decreases when the mission time increases while SE and l decrease as increases. The results get better in the sense that the aim of an accelerated life testing experiments is to get large number of failures (reduce the reliability) of the device with high reliability. In other words, when sample size increases, the rf increases. Also, the RAB 1R for the rf, SE and l decrease when the sample size increases. The rf for Type II censoring is better than the rf for Type I censoring. The hrf increases when the mission time increases. The hrf for Type II censoring is better than Type I censoring.

Application
The main aim of this application is to describe how the proposed method in practice might be used. Cordeiro et al. (2010) used Kolmogrov Smirnov goodness of fit test and data points representing failure time. The data was extracted from Murthy et al. (2004). The data that was used had 30 items (n=30) tested with test stopped after 20th failure (r=20). The failure times were arranged in ascending order and divided into two groups. Each group is put under one level of stress. The first group is exposed to level of stress . The second group is exposed to level of stress . The relationship between the shape parameter and the stress is tested through testing "b", the coefficient. Hypothesis test is obtained when with one degree of freedom. As a result, the null hypothesis (b=0) is rejected and the relationship exists.
The initial values of used in this application are =0.5, b=1.5, =2, =2, and =2. The estimated values of , R( ) and h( ) are obtained. The rf and the hrf are estimated for different values of mission times.
Moreover, the precision and variation of MLEs, rf and hrf are studied through some convenient measures such as the RAB 2 . The credible intervals of the parameters, rf and hrf at confidence level 95% are indicated. The results are displayed in Tables 7 and 8. It is concluded that ̂ is the best estimate among all other estimates while the length interval for the parameter has the shortest length over all the parameters. Also when a mission time increases the rf decreases while on the other hand, when the mission time increases the hrf increases. The lengths interval for rf and hrf increase as mission time increases.

MCMC Method
Multiple levels of integration are necessary to obtain the normalizing constant of f(t) and the marginal posterior densities. For complicated models these integrations are often analytically intractable, and sometimes even a numerical integration cannot be directly obtained. In this case, MCMC simulation is the easiest way to get reliable results without evaluating integrals, (See Gelman et al. (2003)). When the number of iterations is large enough, the samples drawn on one parameter can be regarded as simulated observations from its marginal posterior distribution. Functions of the model parameters, such as rf, hrf and p-th percentile of the lifetime distribution at normal conditions ( ), at low stress ( ) and at high stress ( ), can also be conveniently sampled. Posterior inference can be computed using sample statistics. In this study, WinBUGS, (See

Inference based on Type I censored samples
In this subsection, the inference for unknown parameters , , rf , hrf, ( ), ( ) and ( ) are obtained. This subsection aims to discuss the estimation and optimum problems on the basis of simulated data.

Simulation algorithms
The simulation steps are used as follows:  Accelerated life data from KumW distribution are generated using MathCad 14 at different samples of size (30, 60, 100).
 Two accelerated stress levels and are considered and usual levels is taken as .


To avoid posterior dependence on the starting point of a simulation, several chains with over-dispersed starting points in one MCMC simulation should be run. The simulation converges to the target distribution when traces of all chains appear to be mixing together. In this case, three chains with different initials run simultaneously in one simulation. Each chain continues for 50000 iterations. There are two methods to check convergence. One is examining trace plots of the sample values versus iteration. We can be reasonably confident that convergence has achieved since the three chains appear to be well mixed. The initial values of , b, , and for 3 chains are displayed as follows: The first chain is with initial values =0.3, b=0.5, , and the second chain is with initial values =0.4, b=0.6, , and and the third chain is with initial values =0.5, b=0.7, , and . This is an informal approach to convergence diagnosis. A quantitative way of checking convergence is based on an analysis of variance. The Gelman-Rubin convergence statistic, R, is introduced to diagnose convergence. R is defined as the ratio of the width of the central 80% interval of the pooled chains to the average width of the 80% intervals within the individual chains. When a WinBUGS simulation converges, R should be, or close to one, see Gelman and Rubin (1992). The accuracy of a posterior estimate is calculated in terms of Monte Carlo standard error (MC error) of the mean, which is an estimate of the difference between the mean of the sample and the true posterior mean. According to Spiegelhalter et al. (2003), the simulation should be run until the MC error for each estimate is less than 5% of the sample standard deviation. In this case, this rule was achieved.


Assume that the experiment is terminated once the failure of all the items occur or when a fixed censoring time is reached (Type I censoring  (22).


It appears that the distribution of ( ) is quite asymmetric. Therefore, the median value instead of the mean value is chosen to be the point estimated value for ( ).
 Prior distributions are displayed in Table 9, where prior I (P I) is non informative prior, uniform distribution (U) and prior II (P II) is informative prior, Weibull (W), gamma (G) and exponential (exp) distributions. These priors have the same family of the posterior distribution based on Easy fit program 5.5. The posterior mean and variance of , given t can be calculated, respectively, as follows and p is the p-th percentile.
 It is assumed that the values of the parameters are unknown and Bayesian method is applied to determine the optimal stress changing point. The objective function is to minimize the asymptotic variance of the p-th percentile at normal stress, low stress and high stress. Tables 16 and 17 summarize the optimal stress changing times , optimal failure number and the corresponding posterior variance of at p =0.2.  It is noted that when the sample size, n=100, prior II and Type II censoring are applied, the results of , rf and hrf give approximate value of their initial values.

Concluding remarks
 It is clear that increasing the sample size improves the accuracy of the estimates. Although the stress saves the experiment time, the results are better in usual conditions. Also the number of failures increase as stress increases.


In general, it appears that when using the informative prior, the variance of the estimated parameters decreases. This is expected because the prior knowledge was incorporated with data and increased the accuracy of the estimates. It is also reasonable to conclude that the interval length is narrower with informative priors because the posterior depends on information from two sources: the prior knowledge and the data (via the likelihood function).

Application
Considering the data of this application is given in Subsection 2.3.3. The initial parameter values of a, b, in the 3 chains and the priors distribution used in this application are the same as in the simulation study.


The summary of the sampling results with respect to the unknown parameters a, b, ( ) ( ) ( ), ( ) and ( ), where are displayed in Table 18. A simple summary can be generated showing posterior mean, median, MC error, standard error, a 95% posterior credible interval and length.


The iteration runs until the MC error for each estimate is less than 5% of the sample standard deviation. From Table 18, it is noted that the MC error for each estimate is less than 5% of the sample standard deviation and then the rule of MC error was achieved. Also, to check convergence, Gelman-Rubin convergence statistic, R, is introduced. When a WinBUGS simulation converges, R should be one, or close to one.


The two-sided 95% credible intervals for the estimates of parameters, , rf, hrf, ( ), ( ) and ( ) of KumW are displayed in Table 18. The interval length gets narrower as the sample size increases.


In general, it appears that when using the informative prior, the variance of the estimated parameters decrease. This is expected because the prior knowledge was incorporated with data and increased the accuracy of the estimate. It is also reasonable to conclude that the interval length is narrower with informative priors.  The exponentiated exponential distribution if = .  The Weibull distribution if (Khamis (1997) and Liu (2010)).  The Rayleigh distribution if , = .  The exponential distribution if (Ramadan and Ramadan (2012)).