Analysis of Software Reliability Data using Exponential Power Model

In this paper, Exponential Power (EP) model is proposed to analyze the software reliability data and the present work is an attempt to represent that the model is as software reliability model. The approximate MLE using Artificial Neural Network (ANN) method and the Markov chain Monte Carlo (MCMC) methods are used to estimate the parameters of the EP model. A procedure is developed to estimate the parameters of the EP model using MCMC simulation method in OpenBUGS by incorporating a module into OpenBUGS. The R functions are developed to study the various statistical properties of the proposed model and the output analysis of MCMC samples generated from OpenBUGS. A real software reliability data set is considered for illustration of the proposed methodology under informative set of priors. KeywordsEP model, Probability density, function, Cumulative density function, Hazard rate function, Reliability function, Parameter estimation, MLE, Bayesian estimation.


INTRODUCTION
Exponential models play a central role in analyses of lifetime or survival data, in part because of their convenient statistical theory, their important 'lack of memory' property and their constant hazard rates.In circumstances where the oneparameter family of exponential distributions is not sufficiently broad, a number of wider families such as the gamma, Weibull and lognormal models are in common use.Adding parameters to a well-established family of models is a time honoured device for obtaining more flexible new families of models.The Exponential Power model is introduced by [14] as a lifetime model.This model has been discussed by many authors [4], [9] and [12].; ( , ) 0, x 0 The probability density function (pdf) associated with eq (1) is given by x e e xp 1 e ; ( , ) 0, x 0 We shall write EP(α, λ) to denote Exponential Power model with parameters α λ. α as "shape parameter" by [4] and [14].The R functions dexp.power( ) and pexp.power( ) given in SoftreliaR package can be used for the computation of pdf and cdf, respectively.Some of the typical EP density functions for different values of  and for λ = 1 are depicted in Figure1.It is clear from the Figure 1 that the density function of the Exponential Power model can take different shapes.

1) Mode
The mode can be obtained by solving the non-linear equation 2) The quantile function For a continuous distribution F(x), the p percentile (also referred to as fractile or quantile), x p , for a given p, 0 < p <1, is a number such that pp P(X x ) F(x ) p    .
(4) http://ijacsa.thesai.org/ The quantile for p=0.25 and p=0.75 are called first and third quartiles and the p=0.50 quantile is called the median(Q 2 ).The five parameters are often referred to as the five-number summary or explanatory data analysis.Together, these parameters give a great deal of information about the model in terms of the centre, spread, and skewness.Graphically, the five numbers are often displayed as a boxplot.The quantile function of Exponential Power model can be obtained by solving The computation of quantiles the R function qexp.power( ), given in SoftreliaR package, can be used.In particular, for p=0.5 we get 3) The random deviate generation Let U be the uniform (0,1) random variable and F(.) a cdf for which F -1 (.) exists.Then F -1 (u) is a draw from distribution F(.) .Therefore, the random deviate can be generated from EP() by where u has the U(0, 1) distribution.The R function rexp.power( ), given in SoftreliaR package, generates the random deviate from EP(α, λ).

4) Reliability function/survival function
The reliability/survival function The R function sexp.power( ) given in SoftreliaR package computes the reliability/ survival function.

5) The Hazard Function
The hazard function of Exponential Power model is given by and the allied R function hexp.power( ) given in SoftreliaR package.Since the shape of h(x) depends on the value of the shape parameter .When  ≥ 1, the failure rate function is increasing.When  < 1, the failure rate function is of bathtub shape.Thus the shape parameter  plays an important role for the model.
Since differentiating equation ( 9) w.r.to x, we have Setting h (x)  = 0 and after simplification, we obtain the change point as It easily follows that the sign of   hx  is determined by which is negative for all x ≤ x 0 and positive for all x ≥ x 0 .

6) The cumulative hazard function
The cumulative hazard function H(x) defined as ) can be obtained with the help of pexp.power( ) function given in SoftreliaR package by choosing arguments lower.tail=FALSEand log.p=TRUE.i.e.

-pexp.power(x, alpha, lambda, lower.tail = FALSE, log.p = TRUE) 7) Failure rate average (fra) and Conditional survival function(crf)
Two other relevant functions useful in reliability analysis are failure rate average (fra) and conditional survival function (crf) The failure rate average of X is given by where H(x) is the cumulative hazard function.http://ijacsa.thesai.org/ The survival function (s.f.) and the conditional survival of X are defined by respectively, where F(•) is the cdf of X. Similarly to h(x) and FRA(x), the distribution of X belongs to the new better than used (NBU), exponential, or new worse than used (NWU) classes, when R ( The R functions hra.exp.power( ) and crf.exp.power( ) given in SoftreliaR package can be used for the failure rate average (fra) and conditional survival function(crf), respectively.

II. MAXIMUM LIKELIHOOD ESTIMATION AND INFORMATION MATRIX
Let x=(x 1 , . . ., x n ) be a sample from a distribution with cumulative distribution function (1).The log likelihood function of the parameter L(, λ) is given by Therefore, to obtain the MLE"s of  and λ we can maximize eq.( 15) directly with respect to  and λ or we can solve the following two non-linear equations using iterative procedure [2] and [4]: where I() is the Fisher"s information matrix given by In practice, it is useless that the MLE has asymptotic variance   1 I( ) because we do not know .Hence, we approximate the asymptotic variance by "plugging in" the estimated value of the parameters.The common procedure is to use observed Fisher information matrix Ô( )  (as an estimate of the information matrix I()) given by where H is the Hessian matrix, =() and ˆ= ( , )    .The observed Fisher information is evaluated at MLE rather than determining the expectation of the Hessian at the observed data.This is simply the negative of the Hessian of the loglikelihood at MLE.If the Newton-Raphson algorithm is used to maximize the likelihood then the observed information matrix can easily be calculated.Therefore, the variance-covariance matrix is given by   1 ˆVar( ) cov( , ) H( ) ˆĉov( , ) Var( ) Hence, from the asymptotic normality of MLEs, approximate 100(1-)% confidence intervals for  and can be constructed as where z /2 is the upper percentile of standard normal variate.

III. BAYESIAN ESTIMATION IN OPENBUGS
The most widely used piece of software for applied Bayesian inference is the OpenBUGS.It is a fully extensible modular framework for constructing and analyzing Bayesian full probability models.This open source software requires incorporation of a module (code) to estimate parameters of Exponential Power model.
A module dexp.power_T(alpha,lambda) is written in component Pascal, enables to perform full Bayesian analysis of Exponential Power model into OpenBUGS using the method described in [15] and [16].

A. Implementation of Module -dexp.power_T(alpha, lambda)
The developed module is implemented to obtain the Bayes estimates of the Exponential Power model using MCMC method.The main function of the module is to generate http://ijacsa.thesai.org/MCMC sample from posterior distribution under informative set of priors, i.e.Gamma priors.

1) Data Analysis
We are using software reliability data set SYS2.DAT -86 time-between-failures [10] is considered for illustration of the proposed methodology.In this real data set, Time-betweenfailures is converted to time to failures and scaled.

B. Computation of MLE and Approximate ML estimates using ANN
The Exponential Power model is used to fit this data set.We have started the iterative procedure by maximizing the loglikelihood function given in eq.( 15) directly with an initial guess for  = 0.5 and λ = 0.06, far away from the solution.We have used optim( ) function in R with option Newton-Raphson method.The iterative process stopped only after 7 iterations.We obtain  0.905868898,  0.001531423 and the corresponding log-likelihood value = -592.7172.The similar results are obtained using maxLik package available in R.An estimate of variance-covariance matrix, using eq.( 22), is given by Var( ) cov( , ) 7.265244e-03 -1.474579e-06 ˆ-1.474579e-06 1.266970e-08 ĉov( , ) Var( )  Thus using eq.( 23), we can construct the approximate 95% confidence intervals for the parameters of EP model based on MLE"s.Table I shows the MLE"s with their standard errors and approximate 95% confidence intervals for  and λ.

TABLE I MAXIMUM LIKELIHOOD ESTIMATE, STANDARD ERROR AND 95% CONFIDENCE INTERVAL
An approximate ML estimates based on Artificial Neural Networks are obtained by using the neuralnet package available in R. We have chosen one hidden-layer feedforward neural networks with sigmoid activation function [1].The results are quite close to exact ML estimates.

C. Model Validation
To study the goodness of fit of the Exponential Power model, we compute the Kolmogorov-Smirnov statistic between the empirical distribution function and the fitted distribution function when the parameters are obtained by method of maximum likelihood.For this we can use R function ks.exp.power( ), given in SoftreliaR package.The result of K-S test is D = 0.0514 with the corresponding p-value = 0.9683, Therefore, the high p-value clearly indicates that Exponential Power model can be used to analyze this data set, and we also plot the empirical distribution function and the fitted distribution function in Figure 3. From above result and  The other graphical method widely used for checking whether a fitted model is in agreement with the data is Quantile-Quantile (Q-Q) plots.We run the model to generate two Markov Chains at the length of 40,000 with different starting points of the parameters.The convergence is monitored using trace and ergodic mean plots, we find that the Markov Chain converge together after approximately 2000 observations.Therefore, burnin of 5000 samples is more than enough to erase the effect of starting point(initial values).Finally, samples of size 7000 are formed from the posterior by picking up equally spaced every fifth outcome, i.e. thin=5, starting from 5001.This is done to minimize the auto correlation among the generated deviates.Therefore, we have the posterior sample { 1i , 1i }, i = 1,…,7000 from chain 1 and  2i , 2i }, i = 1,…,7000 from chain 2.
The chain 1 is considered for convergence diagnostics plots.The visual summary is based on posterior sample obtained from chain 2 whereas the numerical summary is presented for both the chains.

A. Convergence diagnostics
Sequential realization of the parameters  and  can be observed in figure 5.The Markov chain is most likely to be sampling from the stationary distribution and is mixing well.There is ample evidence of convergence of chain as the plots show no long upward or downward trends, but look like a horizontal band, then we has evidence that the chain has converged.

2) Running Mean (Ergodic mean) Plot
The convergence pattern based on Ergodic average as shown in figure 6 is obtained after generating a time series (Iteration number) plot of the running mean for each parameter in the chain.The running mean is computed as the mean of all sampled values up to and including that at a given iteration.
Figure 6 The Ergodic mean plots for 

3) Autocorrelation
The graph shows that the correlation is almost negligible.We may conclude that the samples are independent.
Figure7 The autocorrelation plots for 

4) Brooks-Gelman-Rubin Plot
Uses parallel chains with dispersed initial values to test whether they all converge to the same target distribution.Failure could indicate the presence of a multi-mode posterior distribution (different chains converge to different local modes) or the need to run a longer chain (burn-in is yet to be completed).
Figure 8 The BGR plots for  From the Figure 8, it is clear that convergence is achieved.Thus we can obtain the posterior summary statistics.

B. Numerical Summary
In Table II

C. Visual summary 1) Box plots
The boxes represent inter-quartile ranges and the solid black line at the (approximate) centre of each box is the mean; the arms of each box extend to cover the central 95 per cent of the distribution -their ends correspond, therefore, to the 2.5% and 97.5% quantiles.(Note that this representation differs somewhat from the traditional.)Figure 9 The boxplots for alpha and lambda.

2) Kernel density estimates
Histograms can provide insights on symmetric, behaviour in the tails, presence of multi-modal behaviour, and data outliers; histograms can be compared to the fundamental shapes associated with standard analytic distributions.

D. Comparison with MLE
For the comparison with MLE we have plotted two graphs.In Figure 12, the density functions f(x; , ) using MLEs and Bayesian estimates, computed via MCMC samples under gamma priors, are plotted.http://ijacsa.thesai.org/ Figure 12 The density functions f(x; , )  using MLEs and Bayesian estimates, computed via MCMC samples under gamma priors.
The Figure 13, exhibits the estimated reliability function(dashed line) using Bayes estimate under gamma priors and the empirical reliability function(solid line).It is clear from the Figures, the MLEs and the Bayes estimates with respect to the gamma priors are quite close and fit the data very well.

V. CONCLUSION
In this research paper, we have presented the Exponential Power model as software reliability model which was motivated by the fact that the existing models were inadequate to describe the failure process underlying some of the data sets.
We have developed the tools for empirical modelling, e.g., model analysis, model validation and estimation.The exact as well as approximate ML estimates using ANN of the parameters alpha (α) and lambda () have been obtained.
An attempt has been made to estimate the parameters in Bayesian setup using MCMC simulation method under gamma priors.The proposed methodology is illustrated on a real data set.We have presented the numerical summary and visual summary under different priors which includes Box plots, Kernel density estimates based on MCMC samples.The Bayes estimates are compared with MLE.We have shown that the Exponential Power model is suitable for modeling the software reliability data and the tools developed for analysis can also be used for any other type of data sets.

A
model is said to be an Exponential Power model with shape parameter >0 and scale parameter >0, if the survival function of the model is given by Analysis For α > 0 and, λ > 0 the two-parameter Exponential Power model has the distribution function   x F(x; , ) 1 exp 1 e

Figure 1
Figure 1 Plots of the probability density function of the Exponential Power model for =1 and different values of 

Figure 2
Figure 2 Plots of the hazard function of the Exponential Power model for =1 and different values of  Some of the typical Exponential Power Model hazard functions for different values of  and for = 1 are depicted in Figure 2. It is clear from the Figure 2 that the hazard function of the Exponential Power model can take different shapes.

Figure 3 ,
Figure 3, it is clear that the estimated Exponential Power model provides excellent fit to the given data.

Figure 3
Figure 3 The graph of empirical distribution and fitted distribution function.

Figure 4
Figure 4 Quantile-Quantile(Q-Q) plot using MLEs as estimate.The Q-Q plots show the estimated versus the observed quantiles.If the model fits the data well, the pattern of points on the Q-Q plot will exhibit a 45-degree straight line.Note that all the points of a Q-Q plot are inside the square

Figure 5
Figure 5 Sequential realization of the parameters 

Figure 10
Figure 10 Kernel density estimate and histogram of α based on MCMC samples, vertical lines indicates the corresponding ML and Bayes estimates.

Figure 10
Figure 10 and 11 provide the kernel density estimate of  and λ respectively.It can be seen that  and λ both are symmetric.

Figure 11
Figure 11 Histogram and kernel density estimate of λ based on MCMC samples

Figure 13
Figure 13 The estimated reliability function(dashed line) and the empirical reliability function (solid line).

TABLE II NUMERICAL
SUMMARIES UNDER GAMMA PRIORS