Analysis of Gumbel Model for Software Reliability Using Bayesian Paradigm

In this paper, we have illustrated the suitability of Gumbel Model for software reliability data. The model parameters are estimated using likelihood based inferential procedure: classical as well as Bayesian. The quasi NewtonRaphson algorithm is applied to obtain the maximum likelihood estimates and associated probability intervals. The Bayesian estimates of the parameters of Gumbel model are obtained using Markov Chain Monte Carlo(MCMC) simulation method in OpenBUGS(established software for Bayesian analysis using Markov Chain Monte Carlo methods). The R functions are developed to study the statistical properties, model validation and comparison tools of the model and the output analysis of MCMC samples generated from OpenBUGS. Details of applying MCMC to parameter estimation for the Gumbel model are elaborated and a real software reliability data set is considered to illustrate the methods of inference discussed in this paper. KeywordsProbability density function; Bayes Estimation; Hazard Function; MLE; OpenBUGS; Uniform Priors.


INTRODUCTION
A frequently occurring problem in reliability analysis is model selection and related issues.In standard applications like regression analysis, model selection may be related to the number of independent variables to include in a final model.In some applications of statistical extreme value analysis, convergence to some standard extreme-value distributions is crucial.
A choice has occasionally to be made between special cases of distributions versus the more general versions.In this chapter, statistical properties of a recently proposed distribution is examined closer and parameter estimation using maximum likelihood as a classical approach by R functions is performed where comparison is made to Bayesian approach using OpenBUGS.
In reliability theory the Gumbel model is used to model the distribution of the maximum (or the minimum) of a number of samples of various distributions.One of the first scientists to apply the theory was a German mathematician Gumbel [1].Gumbel focused primarily on applications of extreme value theory to engineering problems.The potential applicability of the Gumbel model to represent the distribution of maxima relates to extreme value theory which indicates that it is likely to be useful if the distribution of the underlying sample data is of the normal or exponential type.
The Gumbel model is a particular case of the generalized extreme value distribution (also known as the Fisher-Tippett distribution) [2].It is also known as the log-Weibull model and the double exponential model (which is sometimes used to refer to the Laplace model).
It is often incorrectly labelled as Gompertz model [3,4].The Gumbel model's pdf is skewed to the left, unlike the Weibull model's pdf which is skewed to the right [5,6].The Gumbel model is appropriate for modeling strength, which is sometimes skewed to the left.

II. MODEL ANALYSIS
The two-parameter Gumbel model has one location and one scale parameter.The random variable x follows Gumbel model with the location and scale parameter as - <  < and σ > 0 respectively, if it has the following cummulative distribution function(cdf) The corresponding probability density function (pdf) is Some of the specific characteristics of the Gumbel model are: The shape of the Gumbel model is skewed to the left.The pdf of Gumbel model has no shape parameter.This means that the Gumbel pdf has only one shape, which does not change.
The pdf of Gumbel model has location parameter μ which is equal to the mode but differs from median and mean.This is because the Gumbel model is not symmetrical about its μ.
As μ decreases, the pdf is shifted to the left.As μ increases, the pdf is shifted to the right.www.ijarai.thesai.org where x ( , ), ( , ), 0 It is clear from the Figure 1 that the density function and hazard function of the Gumbel model can take different shapes.

Maximum Likelihood Estimation(MLE) and Information Matrix
To obtain maximum likelihood estimators of the parameters  , σ).Let x1, . . ., xn be a sample from a distribution with cumulative distribution function (2.1).The likelihood function is given by Therefore, to obtain the MLE's of  and σ we can maximize directly with respect to  and σ or we can solve the following two non-linear equations using iterative procedure [8, 9, 10 and 11]:

A. Asymptotic Confidence bounds. based on MLE
Since the MLEs of the unknown parameters   σ) cannot be obtained in closed forms, it is not easy to derive the exact distributions of the MLEs.We can derive the asymptotic confidence intervals of these parameters when  is to assume that the MLE ˆ( , )  are approximately bivariate normal with mean(,σ) and covariance matrix ˆˆĉov( , ) var( ) The above approach is used to derive the 100(1 -, ) as in the following forms /2 ˆẑ Var( ) Here, z is the upper ( /2)th percentile of the standard normal distribution.

B. Data Analysis
In this section we present the analysis of one real data set for illustration of the proposed methodology.The data set contains 36 months of defect-discovery times for a release of Controller Software consisting of about 500,000 lines of code installed on over 100,000 controllers.The defects are those that were present in the code of the particular release of the software, and were discovered as a result of failures reported by users of that release, or possibly of the follow-on release of the product.[13] First we compute the maximum likelihood estimates.

C. Computation of MLE and model validation
The Gumbel model is used to fit this data set.We have started the iterative procedure by maximizing the loglikelihood function given in (3.1) directly with an initial guess for  = 202.0 and  = 145.0,far away from the solution.We have used optim( ) function in R with option Newton-Raphson method.The iterative process stopped only after 1211 iterations.We obtain  212.1565,  151.7684 and the corresponding log-likelihood value = -734.5823.The similar results are obtained using maxLik package available in R.An estimate of variance-covariance matrix, using (3.4) Thus using (3.5), we can construct the approximate 95% confidence intervals for the parameters of Gumbel model based on MLE's.Table 1 shows the MLE's with their standard errors and approximate 95% confidence intervals for  and σ.To check the validity of the model, we compute the Kolmogorov-Smirnov (KS) distance 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.gumbel( ), given in [7].The result of K-S test is D =0.0699 with the corresponding p-value = 0. 0.6501, therefore, the high p-value clearly indicates that Gumbel model can be used to analyze this data set.We also plot the empirical distribution function and the fitted distribution function in Fig. 2. Therefore, it is clear that the estimated Gumbel model provides excellent good fit to the given data.

D. Bayesian Estimation in OpenBUGS
A module dgumbel(mu, sigma) is written in component Pascal, given in [13] enables to perform full Bayesian analysis of Gumbel model into OpenBUGS using the method described in [14,15].

1) Bayesian Analysis under Uniform Priors
The developed module is implemented to obtain the Bayes estimates of the Gumbel model using MCMC method.The main function of the module is to generate MCMC sample from posterior distribution for given set of uniform priors.Which is frequently happens that the experimenter knows in advance that b] but has no strong opinion about any subset of values over this range.In such a case a uniform distribution over [a, b] may be a good approximation of the prior distribution, its p.d.f. is given by 1 We run the two parallel chains for sufficiently large number of iterations, say 5000 the burn-in, until convergence results.Final posterior sample of size 7000 is taken by choosing thinning interval five i.e. every fifth outcome is stored.
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.

E. Convergence diagnostics
Before examining the parameter estimates or performing other inference, it is a good idea to look at plots of the sequential(dependent) realizations of the parameter estimates and plots thereof.We have found that if the Markov chain is not mixing well or is not sampling from the stationary www.ijarai.thesai.orgdistribution, this is usually apparent in sequential plots of one or more realizations.The sequential plot of parameters is the plot that most often exhibits difficulties in the Markov chain.Fig. 3 shows the sequential realizations of the parameters of the model.In this case Markov chain seems to be mixing well enough and is likely to be sampling from the stationary distribution.The plot looks like a horizontal band, with no long upward or downward trends, then we have evidence that the chain has converged.

 Running Mean (Ergodic mean) Plot
In order to study the convergence pattern, we have plotted a time series (iteration number) graph of the running mean for each parameter in the chain.The mean of all sampled values up to and including that at a given iteration gives the running mean.In the Fig. 4 given below, a systematic pattern of convergence based on ergodic averages can be seen after an initial transient behavior of the chain.

G. Visual summary by using Box plots
The boxes represent in Fig. 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.

2) Bayesian Analysis under Gamma Priors
The developed module is implemented to obtain the Bayes estimates of the Gumbel model using MCMC method to generate MCMC sample from posterior distribution for given set of gamma priors, which is most widely used prior distribution of  is the inverted gamma distribution with parameters a and b (>0) with p.d.f.given by We also run the two parallel chains for sufficiently large number of iterations, say 5000 the burn-in, until convergence results.Final posterior sample of size 7000 is taken by choosing thinning interval five i.e. every fifth outcome is stored and same procedure is adopted for analysis as used in the case of uniform priors.

H. Convergence diagnostics
Simulation-based Bayesian inference requires using simulated draws to summarize the posterior distribution or calculate any relevant quantities of interest.We need to treat the simulation draws with care.Trace plots of samples versus the simulation index can be very useful in assessing convergence.The trace indicates if the chain has not yet converged to its stationary distribution-that is, if it needs a longer burn-in period.A trace can also tell whether the chain is mixing well.A chain might have reached stationary if the distribution of points is not changing as the chain progresses.The aspects of stationary that are most recognizable from a trace plot are a relatively constant mean and variance.

 Autocorrelation
The graph shows that the correlation is almost negligible.We may conclude that the samples are independent.www.ijarai.thesai.orgFrom the Fig. 7, it is clear that convergence is achieved.Thus we can obtain the posterior summary statistics.

III. NUMERICAL SUMMARY
In Table 3, we have considered various quantities of interest and their numerical values based on MCMC sample of posterior characteristics for Gumbel model under Gamma priors.

A. Visual summary by using Kernel density estimates
Histograms can provide insights on skewness, 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.

B. Comparison with MLE using Uniform Priors
For the comparison with MLE we have plotted two graphs.

C. Comparison with MLE using Gamma Priors
For the comparison with MLE, we have plotted a graph which exhibits the estimated reliability function (dashed line) using Bayes estimate under gamma priors and the empirical reliability function(solid line).It is clear from Fig. 12, the MLEs and the Bayes estimates with respect to the gamma priors are quite close and fit the data very well.

IV. CONCLUSION
The developed methodology for MLE and Bayesian estimation has been demonstrated on a real data set when both the parameters mu (location) and sigma (scale) of the Gumbel model are unknown under non-informative and informative set of independent priors.The bayes estimates of the said priors, i.e., uniform and gamma have been obtained under squared error, absolute error and zero-one loss functions.A five point summary Minimum (x), Q 1 , Q 2 , Q 3 , Maximum (x) has been computed.The symmetric Bayesian credible intervals and Highest Probability Density (HPD) intervals have been constructed.Through the use of graphical representations the intent is that one can gain a perspective of various meanings and associated interpretations.
The MCMC method provides an alternative method for parameter estimation of the Gumbel model.It is more flexible when compared with the traditional methods such as MLE method.Moreover, 'exact' probability intervals are available rather than relying on estimates of the asymptotic variances.Indeed, the MCMC sample may be used to completely summarize posterior distribution about the parameters, through a kernel estimate.This is also true for any function of the parameters such as hazard function, mean time to failure etc.The MCMC procedure can easily be applied to complex Bayesian modeling relating to Gumbel model

Figure 1 .
Figure 1.Plots of the (a) probability density function and (b) hazard function of the Gumbel model for  =1 and different values of   , σ) is given by

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

Figure 3 .
Figure 3. Sequential realization of the parameters  and .

Figure 4 .
Figure 4.The Ergodic mean plots for mu and sigma.F.Numerical Summary

Figure 5 .
Figure 5.The boxplots for mu and sigma

Figure 6 .
Figure 6.The autocorrelation plots for mu and sigma. Brooks-Gelman-Rubin 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 7 .
Figure 7.The BGR plots for mu and sigma

Figure 8 .
Figure 8. Histogram and kernel density estimate of  based on MCMC samples, vertical lines represent the corresponding MLE and Bayes estimate.

Fig. 8 and
Fig. 8 and Fig. 9 provide the kernel density estimate of  and .The kernel density estimates have been drawn using R with the assumption of Gaussian kernel and properly chosen values of the bandwidths.It can be seen that  and  both are symmetric.

Figure 9 .
Figure 9. Histogram and kernel density estimate of  based on MCMC samples, vertical lines represent the corresponding MLE and Bayes estimate.

Figure 10 .
Figure 10.The density functions ˆf(x; , )  using MLEs and Bayesian estimates, computed via MCMC samples.Whereas, Fig.11represents the Quantile-Quantile(Q-Q) plot of empirical quantiles and theoretical quantiles computed from MLE and Bayes estimates.

Figure 11 .
Figure 11.Quantile-Quantile(Q-Q) plot of empirical quantiles and theoretical quantiles computed from MLE and Bayes estimates.It is clear from the Figures, the MLEs and the Bayes estimates with respect to the uniform priors are quite close and fit the data very well.

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