Bayesian Reliability Analysis of the Power Law Process with Respect to the Higgins-Tsokos Loss Function for Modeling Software Failure Times
BBayesian Reliability Analysis Of The Power LawProcess With Respect To The Higgins-Tsokos LossFunction For Modeling Software Failure Times
Freeh Alenezi and Chris. Tsokos PhD Candidtae, Department of Mathematics and Statistics,
College of Arts andSciences , University of South Florida , Tampa, 33647, USA Distinguished University Professor, Department of Mathematics and Statistics,
College of Arts and Sciences , University of South Florida , Tampa, 33647, USA
The Power Law Process, also known as Non-Homogeneous Poisson Process, has been used in variousaspects, one of which is the software reliability assessment. Specifically, by using its intensity functionto compute the rate of change of a software reliability as time-varying function. Justification ofBayesian analysis applicability to the Power Law Process was shown using real data. The probabilitydistribution that best characterizes the behavior of the key parameter of the intensity function was firstidentified, then the likelihood-based Bayesian reliability estimate of the Power Law Process under theHiggins-Tsokos loss function was obtained. As a result of a simulation study and using real data, theBayesian estimate shows an outstanding performance compared to the maximum likelihood estimateusing different sample sizes. In addition, a sensitivity analysis was performed, resulting in the Bayesianestimate being sensitive to the prior selection; whether parametric or non-parametric.
Index terms—
Reliability growth; intensity function; non-homogeneous Poisson process; kerneldensity; loss function; robustness
Software reliability growth is often tested during the software development process to insure a goodquality product. Repairable software is tested until a failure is detected, then fixed, and tested againuntil a new failure is detected. This reliability improvement of software has been studied for decades.Duane (1964) [1] introduced the "learning curve approach", which is a plot of the failure rate (or1 a r X i v : . [ s t a t . C O ] F e b ayesian Reliability Alenezi & Tsokosthe intensity function) of a system as a function of time. It is used to assess software reliabilityimprovement over time. For example, software reliability has improved when we observe a negativecurve, whereas a positive curve means that reliability is deteriorating. Stability in software reliabilityis achieved when there is no curve, i.e. the graph is a horizontal line. The number of failures in theinterval (0 , t ] , N ( t ) , is considered a Poisson counting process after satisfying the following conditions[9, 6]:1. N ( t = 0) =0.2. Independent increment (counts of disjoint time intervals are independent).3. It has an intensity function V ( t ) = lim ∆ t → P ( N ( t,t +∆ t )=1)∆ t .4. Simultaneous failures do not exist ( lim ∆ t → . P ( N ( t,t +∆ t )=2)∆ t =0).The probability of a random value N ( t ) =n is given by: P ( N ( t ) = n ) = exp (cid:110) − (cid:82) t V ( t ) dt (cid:111) (cid:110)(cid:82) t V ( t ) dt (cid:111) n n ! , t > . (2.1)Crow (1974) proposed a non-homogeneous Poisson process (NHPP) [3], which is a Poisson processwith a time-varying intensity function, given by: V ( t ) = V ( t ; β, θ ) = βθ (cid:18) tθ (cid:19) β − , t > , β > , θ > , (2.2)with β and θ as the shape and scale parameters, respectively. This NHPP is also known as the powerlaw process (PLP).The joint probability density function (PDF) of the ordered failure times T , T , ..., T n from anNHPP with intensity function V ( t ; β, θ ) is given by: f ( t , .., t n ) = n (cid:89) i =1 V ( t i ; β, θ ) exp (cid:26) − (cid:90) w V ( t ; β, θ ) dt (cid:27) , (2.3)where w is the so-called stopping time. Considering the failure truncation case (w = t n ), the conditionalreliability function of the failure time T n given T = t , T = t , T = t ,..., T n − = t n − , T n − = t n − is a function of V ( t ; β, θ ) . 2ayesian Reliability Alenezi & TsokosTo monitor software reliability growth over time, an engineer can use the estimate of the β value,the key parameter in the intensity function, since it plays a significant role during the testing process.For β > , the number of failures would increase because the intensity function is increasing. On theother hand, if the intensity function is decreasing, β < means that the number of failures woulddecrease, indicating improved software reliability. Note that in the case of a homogeneous Poissonprocess pertains when β = 1 , in which case the intensity function will be θ and whatever changes havebeen made have had no effect on the outcome.The NHPP has been used for analyzing software failure times, and for predicting the next failureevent. Several publications show the effectiveness and usefulness of this model in assessing reliabilitygrowth [2, 10, 11, 6]. In addition, NHPP has been used to study drug effectiveness in breast cancertreatment [21] and in the formulation of a software cost model [13].Since the intensity function is driving the NHPP, improving the existing methods to estimate thekey parameter β will certainly improve the accuracy of reliability growth assessment and help thestructuring of maintenance strategies. Molinares and Tsokos [5], obtained a Bayesian estimate of theparameter β and compared it with its approximate maximum likelihood estimate (MLE). The authorsderived the Bayesian estimates with respect to squared-error loss function, using Burr, Jeffreys, andinverted gamma probability distributions as the prior PDFs for β .In performing Bayesian analysis on a real world problem, we need some sort of justification forpursuing this particular type of analysis. Once we have identified the probability distribution thatcharacterized the probabilistic behavior of the failure time, we need to identify the prior PDF of β and a loss function. The squared-error loss function is the most popular loss function used in Bayesiananalysis because of its analytical tractability. It places a small weight on the estimates around thetrue value, but proportionally more weight on estimates far from the true value. Higgins and Tsokos[7] proposed a new loss function that places exponential weight on extreme deviations from the truevalue, while remaining mathematically tractable.In the present study, we investigate the effectiveness of Bayesian analysis in using the Higgins-Tsokos (H-T) loss function (that puts the loss at the end of the process) for modeling software failuretimes. To accomplish this, we use the NHPP as the underlying failure distribution subject to using theBurr PDF as a prior of β . In addition, we utilize the H-T loss function to perform sensitive analysisof prior selections. We employ parametric and non-parametric priors, namely Burr, inverted gamma,Jeffery, and two kernel PDFs. Therefore, the primary objective of the study is to answer the followingquestions within a Bayesian framework: 3ayesian Reliability Alenezi & Tsokos1. What is the performance of the Bayesian estimate of β under the H-T loss function compared toits MLE when modeling software failure times using PLP?2. Is the Bayesian estimate of β , using the H-T loss function in the PLP, sensitive to the selectionof the prior PDF, both parametric and non-parametric?The paper is organized as follows: Section 2 describes the theory and development of the Bayesianreliability model; Section 3 presents the results and discussion; Section are is the conclusions. The probability of achieving n failures of a given system in the time interval (0 , t ] can be written as,[5, 6]: P ( x = n ; t ) = exp (cid:110) − (cid:82) t V ( t ; β, θ ) dt (cid:111) (cid:110)(cid:82) t V ( t ; β, θ ) dt (cid:111) n n ! , t > , (3.1.1)where V ( t ; β, θ ) is the intensity function given by (2.2). The reduced expression P ( x = n ; t ) = 1 n ! exp (cid:26) − tθ β (cid:27) tθ nβ , (3.1.2)is the PLP that is commonly known as the Weibull or NHPP.If the PLP is the underlying failure model of the failure times t , t , t ,... , t n − , and t n , theconditional reliability function of t n given t , t , t ,... , t n − can be written as, [5, 6]: R ( t n | t , t , ..., t n − ) = exp (cid:40)(cid:90) t n t n − V ( t ; β, θ ) dt (cid:41) , tn > t n − > , (3.1.3)since it is independent of t , t , t , ... , t n − .Since the reliability function, equation (3.1.3), is written mathematically as a function of the intensityfunction, estimating the parameter β in the V ( t ; β, θ ) leads to estimation of the reliability function.The (MLE) of β is a function of the largest failure time and the MLE of θ is also a function of the MLEof β . Let T , T , ..., T n denote the first n failure times of the PLP, where T l < T < ... < T n are totaltimes since the initial startup of the system. Thus, the truncated conditional PDF, f i ( t | t , ..., t i − ) , in4ayesian Reliability Alenezi & Tsokosthe Weibull process and is given by, [5, 6]: f i ( t | t , ..., t i − ) = βθ (cid:18) tθ (cid:19) β − exp (cid:26) − tθ β + t i − θ β (cid:27) , t i − < t. (3.1.4)With t = ( t , t , ..., t n ) , the likelihood function for the first n failure times of the PLP T = t , T = t , ..., T n = t n can be written as: L ( t ; β ) = exp (cid:32) − (cid:18) t n θ (cid:19) β (cid:33) (cid:18) βθ (cid:19) n n (cid:89) i =1 (cid:18) t i θ (cid:19) β − . (3.1.5)The MLE for the shape parameter is given by, [2, 3, 5, 6]: ˆ β n = n (cid:80) ni =1 log (cid:16) t n t i (cid:17) , (3.1.6)and for the scale parameter is: ˆ θ n = t n n / ˆ β n . (3.1.7)Note that the MLE of θ depends on the MLE of β using the largest (last) observed failure time. Crow [2, 3] failure data from a system undergoing developmental testing was used, by Molinares &Tsokos [5], to show how β varied depending on the last failure time (largest time), thus they proposeda Bayesian approach to the PLP. The authors also found that the MLE of β follows a four-parameterBurr probability distribution, g ( β ; α, γ, δ, κ ) , known as the four-parameter Burr type XII probabilitydistribution, with a PDF given by: g B ( β ) = g ( β ; α, γ, δ, κ ) = ακ ( β − γδ ) α − δ ( ( β − γδ ) α ) κ +1 , γ ≤ β < ∞ , otherwise , (3.2.1)where the hyperparameters α , γ , δ and κ are being estimated using MLE in the goodness of fit (GOF)test applied to the β estimates. The Crow successive failure data for his system is given in Table 1.5ayesian Reliability Alenezi & Tsokos Table 1: Crow’s failure times of a system under development.
Failure times0.7 3.7 13.2 17.6 54.5 99.2 112.2120.9 151 163 174.5 191.6 282.8 355.2486.3 490.5 513.3 558.4 678.1 688 785.9887 1010.7 1029.1 1034.4 1136.1 1178.9 1259.71297.9 1419.7 1571.7 1629.8 1702.4 1928.9 2072.32525.2 2928.5 3016.4 3181 3256.3 – –According to the reliability growth failure data, the system failed for the first time at . units of time, t = 0 . , and it failed the 40 th time at . units of time, t = 3256 . . The MLE of the parameter β for n = 40 is, [5, 6]: ˆ β = 40 (cid:80) i =1 log (cid:16) . t i (cid:17) (cid:39) . . (3.2.2)If β were treated in a non-Bayesian setting, its MLE would be given by Eq. (3.2.2).In an experimental process, the largest time to failure could occur at any point in the series offailures for a given system. Therefore, consider the case where the largest failure is t = 3181 . Insuch a case, the estimate of β is 0.48.The largest failure time always affects the MLE of β . Thus, it is recommended that β not to bethought of as an unknown constant [5], but rather as an unknown random variable. This recommenda-tion provides the opportunity to study Bayesian analysis in the PLP with respect to various selectionsof loss functions and priors.The Bayesian estimates of β will be derived using H-T loss functions. The H-T loss function (1976) is given by, [5]: L ( ˆ ξ, ξ ) = f exp (cid:110) f ( ˆ ξ − ξ ) (cid:111) + f exp (cid:110) − f ( ˆ ξ − ξ ) (cid:111) f + f − , f , f > . (3.2.1)Higgins and Tsokos [7] showed that it places more weight on the extreme underestimation andoverestimation of the true value when f > f and f < f , respectively. The risk using the H-T lossfunction, where ξ = β represents the estimate of ˆ ξ = ˆ β , is given by: E [ L ( ˆ β, β )] = (cid:90) ∞−∞ [ f exp (cid:110) f ( ˆ β − β ) (cid:111) + f exp (cid:110) − f ( ˆ β − β ) (cid:111) f + f − h ( β | t ) dβ. (3.2.2) E [ L ( ˆ β, β )] with respect to β and setting it equal to zero we solve for ˆ β , theBayesian estimate of β with respect to the H-T loss function, is given by: ˆ β B.T H = 1 f + f ln[ (cid:82) ∞−∞ exp { f β } h ( β | t ) dβ (cid:82) ∞−∞ exp {− f β } h ( β | t ) dβ ] . (3.2.3)The Bayesian estimate of β with respect to the H-T loss function and Burr probability distribution,as the prior, has h ( β | t ) given by: h ( β | t ) = (cid:82) ∞ γ ( βθ ) n exp (cid:110) − (cid:0) t n θ (cid:1) β (cid:111) (cid:81) ni =1 (cid:0) t i θ (cid:1) β − β − γδ ) α − (1+( β − γδ ) α ) κ +1 dβ (cid:82) ∞ γ ( βθ ) n exp (cid:110) − (cid:0) t n θ (cid:1) β (cid:111) (cid:81) ni =1 (cid:0) t i θ (cid:1) β − β − γδ ) α − (1+( β − γδ ) α ) κ +1 dβ . (3.2.4)With the use of Eq. (3.1.3), the conditional reliability of t i , the analytical structure of the con-ditional Bayesian reliability estimate for the PLP that is subject to the above information, is givenby: ˆ R B ( t i | t , t , ..., t i − ) = exp (cid:40) − (cid:90) t i t i − ˆ V (cid:48) B ( t ; β, θ ) dt (cid:41) , t i > t i − > , (3.2.5)where ˆ V (cid:48) B ( t ; β B.T H , θ ) = ˆ β B.T H θ (cid:18) tθ (cid:19) ˆ β B.TH − , θ > , t > , (3.2.6)where ˆ β B.T H is the Bayesian estimate of β using the H-T loss function. We are also interested incomparing the Bayesian estimate, using the H-T loss function, with MLE of the subject parameter fordifferent parametric and non-parametric priors, assuming β has a random behavior and θ is known;and also comparing Eq. (3.1.7) with an adjusted MLE considered as a function of β .7ayesian Reliability Alenezi & Tsokos In this section, we seek the answer to the following question: Is the Bayesian estimate of β , using theH-T loss function in the PLP, sensitive to the selection of the prior, with parametric or non-parametricpriors? Assuming β is a random variable, using simulated data, sensitivity analysis was done for thefollowing parametric and non-parametric priors:1. Jeffreys’ prior [[8]]:Jeffreys’ prior is proportional to the square root of the determinant of the Fisher informationmatrix ( I ( β )). It is a non-informative prior, where the Jeffreys’ prior for the PLP, consideringthat β , I ( β ) is scalar in this case, is given by: g J ( β ) ∝ (cid:112) I ( β ) = (cid:115) − E ( ∂ LogL ( t ; β ) ∂β ) ∝ β , β > . (3.3.1)2. The inverted gamma:The PLP and inverted gamma probability distributions belong to the exponential family ofprobability distributions, which makes the latter a logical choice for an informative parametricprior for β . The inverted gamma probability distribution is given by: g IG ( β ) ∝ (cid:18) µβ (cid:19) v +1 µ Γ( v ) exp (cid:26) − µβ (cid:27) , β > , µ > , v > , (3.3.2)where v and µ are the shape and scale parameters.3. Kernel’ prior:The kernel probability density estimation is a non-parametric method to approximately estimatethe PDF of β using a finite data set. It is given by: g K ( β ) = 1 nh n (cid:88) i =1 K (cid:18) β − β i h (cid:19) , (3.3.3)where K is the kernel function and h is a positive number called the bandwidth.8ayesian Reliability Alenezi & Tsokos Assuming Jeffreys’ PDF, Eq. (3.3.1), as the prior of β and using the likelihood function (3.1.5), theposterior density of β is given by: h J (¯ t | β ) = exp (cid:110)(cid:0) t n θ (cid:1) β (cid:111) β n − θ nβ (cid:81) ni =1 ( t i ) β − (cid:82) ∞ exp (cid:110)(cid:0) t n θ (cid:1) β (cid:111) β n − θ nβ (cid:81) ni =1 ( t i ) β − dβ . (3.3.1.1)Thus, the Jeffreys’ Bayesian estimate of β in V ( t ; β, θ ) under the H-T loss function, using Eq. (3.2.3),is given by: ˆ β J.HT = 1 f + f ln[ (cid:82) ∞ γ exp { f β } h J (¯ t | β ) dβ (cid:82) ∞ γ exp {− f β } h J (¯ t | β ) dβ ] . (3.3.1.2)We cannot obtain a closed analytical form of the Bayesian estimate, ˆ β J.HT , thus we must utilizenumerical method to obtain the subject estimate. Also note that the estimate depends on knowing orbeing able to estimate the scale parameter θ . We proceed with our study with the prior probability density of β given by the inverted gammadistribution Eq. (3.3.2). Using the likelihood Eq. (3.1.5), the posterior density of β is given by: h IG ( t | β ) = β n − v − θ nβ exp (cid:110) − (cid:0) t n θ (cid:1) β − µβ (cid:111) (cid:81) ni =1 ( t i ) β − (cid:82) ∞ β n − v − θ nβ exp (cid:110) − (cid:0) t n θ (cid:1) β − µβ (cid:111) (cid:81) ni =1 ( t i ) β − dβ . (3.3.2.1)Thus, the Bayesian estimate of β under the inverted gamma distribution with respect to the H-T lossfunction, using Eq. (3.2.3) and Eq. (3.3.2.1), is given by: ˆ β IG.HT = 1 f + f ln[ (cid:82) ∞ γ exp { f β } h IG ( t | β ) dβ (cid:82) ∞ γ exp {− f β } h IG ( t | β ) dβ ] . (3.3.2.2)Here as well, we must rely on a numerical estimation of ˆ β IG.HT because we cannot obtain a closedform of the above equation. Also note that the estimate depends on knowing or being able to estimatethe scale parameter θ . 9ayesian Reliability Alenezi & Tsokos Here, we shall assume the non-parametric kernel probability density Eq. (3.3.3) as the prior PDF of β ; using the likelihood Eq. (3.1.5), the posterior density of β is given by: h k (¯ t | β ) = exp (cid:110)(cid:0) t n θ (cid:1) β (cid:111) β n θ nβ (cid:81) ni =1 ( t i ) β − nh (cid:80) ni =1 K (cid:18) β − β i h (cid:19)(cid:82) ∞ exp (cid:110)(cid:0) t n θ (cid:1) β (cid:111) β n θ nβ (cid:81) ni =1 ( t i ) β − nh (cid:80) ni =1 K (cid:18) β − β i h (cid:19) dβ . (3.3.3.1)Thus, the kernel Bayesian estimate of the key parameter β in V ( t ; β, θ ) under the H-T loss function,using Eq. (3.2.3) and Eq. (3.3.1.1), is given by: ˆ β K.HT = 1 f + f ln[ (cid:82) ∞ γ exp { f β } h k (¯ t | β ) dβ (cid:82) ∞ γ exp {− f β } h k (¯ t | β ) dβ ] . (3.3.3.2)We must rely on a numerical estimation because we cannot obtain a closed form solution for ˆ β K.HT .In addition, the kernel function, K ( u ) , and bandwidth, h , will be chosen to minimize the asymptoticmean integrated squared error (AMISE) given by: AM ISE (cid:16) ˆ f ( β ) (cid:17) = (cid:90) E (cid:20)(cid:16) ˆ f ( β ) − f ( β ) (cid:17) (cid:21) dβ, (3.3.3.3)where ˆ f ( β ) and f ( β ) are the estimated probability density of β and the true probability density of β respectively. Below are details of the analysis we conducted using Monte Carlo simulation to generatedata governed by a PLP, followed by using actual data. A Monte Carlo simulation was used to compare the Bayesian (under H-T loss functions) and the MLEapproaches. The parameter β of the intensity function for the PLP was calculated using numericalintegration techniques in conjunction with a Monte Carlo simulation to obtain its Bayesian estimate.Substituting these estimates in the intensity function, we obtained the Bayesian intensity functionestimates, from which the reliability function can be estimated.For a given value of the parameter θ , a stochastic value for the parameter β was generated fromthe Burr PDF. For each pair of values of θ and β , samples of failure times that followed a PLP10ayesian Reliability Alenezi & Tsokoswere generated. This procedure was repeated times for three distinct values of θ . The procedureis summarized in the algorithm (Algorithm 1) given below. StartInitialize the parameter θ and number of iterations p Generate β [ k ] from Burr DistributionGenerate (cid:126)t [ k ] from PLP using β [ k ] Compute MLE of β [ k ] , named ˆ β [ k ] k = , ,..., p Compute Bayesian estimate, ˆ β [ k ] B.HT , of β [ k ] Calculate MSE of ˆ β B.HT
Calculate MSE of ˆ β of β End
Algorithm 1.
Simulation to analyze Bayesian estimates of β for a given θ . For each sample of size , the Bayesian estimates and MLEs of the parameter were calculatedwhen θ ∈ { . , . , } . The comparison is based on the mean squared error (MSE) averaged overthe , repetitions. The results are given in Table 2. Table 2: MSE for Bayesian estimates under the H-T lossfunction and MLE of β , for each assumed θ value. θ MSE of ˆ β MSE of ˆ β B.HT ˆ β B.HT maintains a good accuracy, and is superior to ˆ β in estimating β for thedifferent values of θ . For various sample sizes, the Bayesian estimate under the H-T loss function andthe MLE of the parameter β were calculated and averaged over , repetitions. Table 3 displaysthe simulated result of comparing a true value of β with respect to its MLE and Bayesian estimatesfor n = 20 , , ... , .Again, the Bayesian estimate is uniformly closer to the true value of β than its MLE, even fora very small sample size of n = 20 . A graphical comparison of the true value of β along with the11ayesian Reliability Alenezi & Tsokos Table 3: Bayesian estimates, under H-T loss function, and MLEs forthe parameter β = 0.7054 averaged over 10,000 repetitions n β F ixed ˆ β ˆ β B.HT
20 0.7054 0.784026 0.67526330 0.7054 0.756617 0.69018940 0.7054 0.743982 0.69646750 0.7054 0.73531 0.69915860 0.7054 0.729563 0.70064270 0.7054 0.725977 0.7016980 0.7054 0.723338 0.702382100 0.7054 0.719117 0.703165120 0.7054 0.716315 0.703585140 0.7054 0.714821 0.70398160 0.7054 0.713641 0.704244Bayesian and MLE estimates as functions of sample size is given by Figure 2.
Averagesof β estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆● Fixed ■ MLE ◆ Bayesian.HT
20 40 60 80 100 120 140 1600.680.700.720.740.760.78
Sample size
Figure 2: β estimates versus sample size. Figure 2 shows the MLE of β tends to overestimate whereas the Bayesian estimate tends to under-estimate the true value of β , particularly when considering small sample sizes. The MSEs of the MLEand Bayesian estimates of β is given below by Figure 3.Regardless of sample size, the MSE of the Bayesian estimate of the key parameter β is significantlysmaller than the MSE of the MLE of β ( Figure 3).Since the Bayesian estimate under the H-T loss function for β is better than its MLE, Molinaresand Tsokos proposed to adjust the MLE of the parameter θ using (3.1.6) with a Bayesian estimate of12ayesian Reliability Alenezi & Tsokos MSEsof β estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■● MLE ■ Bayesian.HT
20 40 60 80 100 120 140 1600.010.020.030.04
Sample size
Figure 3: β estimates versus sample size. β instead of its MLE, both of which are needed to estimate the V ( t ; β, θ ) , as given below: ˆ θ B.HT = t n n / ˆ β B.HT . (4.1.1)For various sample sizes and the same β ( β = 0 . ), the Bayesian MLE and MLE of the parameter θ and their corresponding MSEs were computed, averaging over the , repetitions, using the MLEof θ ( θ = 1 . ) of the Crow data. Table 4: MLE and Bayesian estimates under the H-T loss function forthe parameter θ = 1.7441 averaged over 10,000 repetitions. n θ ˆ θ MLE ˆ θ B.HT
20 1.7441 3.17139 1.3642230 1.7441 2.908 1.509740 1.7441 2.73107 1.5811550 1.7441 2.59245 1.6198560 1.7441 2.48865 1.6440670 1.7441 2.41782 1.6608480 1.7441 2.36522 1.67294100 1.7441 2.26774 1.68902120 1.7441 2.20117 1.69923140 1.7441 2.15539 1.70659160 1.7441 2.11872 1.71193Table 4 shows the inferior performance for the MLE of θ and the slow convergence of its averagevalues to θ = 1 . , whereas the adjusted estimate of θ ( ˆ θ B.HT ) using the Bayesian estimate of β under the H-T loss function performed better in estimating the true value of θ . The MLE of andthe Bayesian MLE estimate of θ had a tendency to overestimate and underestimate the parameter θ ,13ayesian Reliability Alenezi & Tsokosrespectively.As expected, based on the Bayesian influence on β , ˆ θ B.HT is a better estimate than the MLE of θ ( ˆ θ ). This can be seen in Figure 4 where the MSEs of θ estimates were ploted against various samplesizes, which demonstrates the excellent performance of ˆ θ B.HT . MSEsof θ estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■● MLE ■ Bayesian.HT
20 40 60 80 100 120 140 16024681012
Sample size
Figure 4: MSE of θ : Bayesian and MLE estimates versus sample size. We also computed the proposed estimate for the parameter θ ( ˆ θ B.HT ) and its MSE over , repetitions for different values of θ (0.5, 1.7441, 4) and sample size n = 40 . The results are given byTable 5. The θ values (including 1.7441) were selected for this simulation are smaller and larger thanthe MLE of θ of the Crow data.Table 5 below shows that the ˆ θ B.HT performed well for the selected θ values. This is particularlytrue for the small and medium value of θ values. Table 5: MSE of θ estimates: Bayesian under the H-T loss function,and MLE of β . θ ˆ θ B.HT
MSE of ˆ θ B.HT
For a fixed value of θ = 1 . and a sample size similar to the size of the collected data, n = 40 ,the estimates of the intensity function ˆ V MLE ( t ) and ˆ V B.HT ( t ) were obtained using ˆ β and ˆ β B.HT ,respectively, in Eq. (2.2). That is, ˆ V (cid:48) MLE ( t ) = ˆ βθ (cid:18) tθ (cid:19) ˆ β − , θ > , t > . (4.1.2)14ayesian Reliability Alenezi & Tsokos V ( t ) V ( t ) V ′ B .HT ( t ) V MLE ( t ) t Figure 5: Graph for θ = 1.7441 and the corresponding β Bayesianestimate and MLE’s used in ˆ V (cid:48) MLE and ˆ V (cid:48) B.HT , estimates of V ( t ; β, θ ) with n = 40. and ˆ V (cid:48) B.HT ( t ) = ˆ β B.HT θ (cid:18) tθ (cid:19) ˆ β B.HT − , θ > , t > . (4.1.3)Their graphs (Figure 5) reveal the superior performance of ˆ V (cid:48) B.HT ( t ) .In order to obtain Bayesian estimates of the intensity function, ˆ V ∗ B.SE and ˆ V ∗ B.HT , we substitutedthe Bayesian estimates of β and its corresponding θ MLE in Eq. (2.2). That is, ˆ V ∗ B.HT ( t ) = ˆ β B.HT ˆ θ (cid:18) t ˆ θ (cid:19) ˆ β B.HT − , t > . (4.1.4)The MLE of the intensity function, ˆ V MLE , is obtained using the MLEs of β and θ . That is, ˆ V MLE ( t ) = ˆ β ˆ θ (cid:18) t ˆ θ (cid:19) ˆ β − , t > . (4.1.5)The Bayesian MLE of the intensity function under the influence of the Bayesian estimate of β ,denoted by ˆ V B.HT , is obtained by substituting ˆ β B.HT and ˆ θ B.HT in Eq. (2.2): ˆ V B.HT ( t ) = ˆ β B.HT ˆ θ B.HT (cid:18) t ˆ θ B.HT (cid:19) ˆ β B.HT − , t > . (4.1.6)To measure the robustness of ˆ V B.HT with respect to ˆ V MLE , we calculated the relative efficiency15ayesian Reliability Alenezi & Tsokos(RE) of the estimate ˆ V B.HT compared to the estimate ˆ V MLE , which is defined as RE ( ˆ V B.HT , ˆ V MLE ) =
IM SE ( ˆ V B.HT ) IM SE ( ˆ V MLE ) = (cid:82) ∞−∞ [ ˆ V B.HT ( t ) − V ( t )] dt (cid:82) ∞−∞ [ ˆ V MLE ( t ) − V ( t )] dt . (4.1.7)If RE = 1 , ˆ V B.HT and ˆ V MLE will be interpreted as equally effective. If
RE < , ˆ V B.HT is moreefficient than ˆ V MLE , contrary to
RE > , in which case ˆ V B.HT is less efficient than ˆ V MLE .Bayesian estimates and MLEs for the parameters β = 0.7054 and θ =1.7441 (Table 6), averaged over10,000 repetitions, were used, for n = 40 , to compare ˆ V B.T H and ˆ V MLE using Eq. (4.1.7).
Table 6: Averages of the Bayesian (under the H-T loss function) and MLEestimates of β and θ β ˆ β ˆ β B.HT θ ˆ θ ˆ θ B.HT β and θ under the H-T loss function, respectively, are closer to the true values than theircorresponding MLE estimates. The results of the comparison among the ˆ V B.T H and ˆ V MLE using Eq.(4.1.7) are given in Tables 7 and 8.The analytical forms of the V ( t ) , ˆ V MLE , ˆ V ∗ B.T H , and ˆ V B.T H (Table 7) were derived by substitutingthe initialized values, MLE estimates of both parameters, Bayesian estimate β and MLE of θ , andBayesian estimates, respectively.Table 8 shows the comparison result of ˆ V B.HT and ˆ V MLE , where the RE ( ˆ V B.HT , ˆ V MLE ) is less than1, which implies that the intensity function using ˆ β B.HT is more efficient than the intensity functionunder ˆ β MLE , establishing the superior relative efficiency of Bayesian estimates under the H-T lossfunction over MLE estimates. The corresponding graph for the intensity functions is given by Figure6. In addition, ˆ V ∗ B.HT computed using a Bayesian estimate for β and MLE estimate for θ , is lessefficient compared to ˆ V MLE , and ˆ V B.HT .The ˆ V B.HT is a better estimate of V ( t ) , compared to the ˆ V MSE and ˆ V ∗ B.HT . Based on the resultsof this section, the Bayesian estimates under the H-T loss function will be used to analyze real datain the following section.
Table 7: Intensity functions with Bayesian and MLE estimates for β and θ V ( t ) ˆ V MLE ˆ V ∗ B.HT ˆ V B.HT . · t − . . · t − . . · t − . . · t − . Table 8: Relative efficiency of ˆ V B.HT compared to ˆ V MLE . RE ( ˆ V B.HT , ˆ V MLE ) RE ( ˆ V B.HT , ˆ V ∗ B.HT ) V ( t ) V ( t ) V * B .HT ( t ) V MLE ( t ) V B .HT ( t ) t Figure 6: Estimates of the intensity function using values in Table 6, n= 40.
Using the software reliability growth data from Table 1, we computed ˆ β B.HT and the adjusted estimateof θ ( ˆ θ B.HT ) in order to obtain a Bayesian estimate of the intensity function under the H-T loss function.We followed the algorithm given below (Algorithm 2) to obtain the Bayesian intensity function for thegiven real data.For the failure data of Crow, provided in Table 1, ˆ β B.HT is . and ˆ θ B.HT is . . There-fore, with the use of ˆ θ B.HT , the Bayesian MLE of the intensity function for the data is given by: ˆ V B.HT ( t ) = 0 . · t − . , t > . (4.2.1)A graphical display of ˆ V B.HT ( t ) is given below.Figure 8 shows the Bayesian MLE estimate of the intensity function ( V ( t ; β, θ ) ) under the H-T lossfunction ( ˆ V ∗ B.HT ), which indicates the improvement of the software reliability over time.To obtain a Bayesian MLE for the reliability function under the H-T loss function, we use thisBayesian estimate for the intensity function. The analytical form for the corresponding Bayesian17ayesian Reliability Alenezi & Tsokos
StartInitialize: (cid:126)t [ k ] = ( t , ...., t n ) as vector of failure times t n is the largest failure time (cid:126) ˆ β [ k ] as vector of MLEs of β Do Goodness of Fit test to fit a PDF g ( β ) for (cid:126) ˆ β [ k ] Compute Bayesian estimate of β under H-T loss function: L ( (cid:126)t | β ) as likelihood function of (cid:126)th ( β | (cid:126)t ) as posterior distribution of β using L ( (cid:126)t | β ) and g ( β )ˆ β B.HT as the Bayesian estimate of β using h ( β | (cid:126)t ) Compute the adjusted MLE of θ , ˆ θ B.HT , using ˆ β B.HT
Obtain the analytical form of the Bayesian MLE ofthe intensity function, ˆ V B.HT , using ˆ β B.HT and ˆ θ B.HT
End
Algorithm 2.
Estimate of the intensity function using Crow data in Table 1 . reliability estimate, based on the real data, is given by: ˆ R B.HT ( t i | t , .., t i − ) = exp (cid:40) − . (cid:90) t i t i − x − . dx (cid:41) , t i > t i − > . (4.2.2)Thus far, we demonstrated not only the applicability of the Bayesian analysis to the PLP, butalso, using real data, the superiority of its performance and influence compared to the MLE of theparameters β and θ , respectively, assuming the Burr PDF is the prior knowledge of the key parameter β . Next section, we study the sensivity of the prior selections, in which an engineer might lack a priorknowledge of parameter β . In the implementation of the simulation procedure we followed Algorithm 1. Random failure times(time to failures) distributed according to the PLP are simulated for a realization of the stochastic scaleparameter β , which follows a Burr type XII probability distribution. Informative parametric priors wereconsidered, such as inverted gamma and Burr probability distributions, whereas the Jefferys prior was18ayesian Reliability Alenezi & Tsokos V B.HT [ t ] t Figure 8: Estimate of the intensity function for the real data in Table1, using ˆ β B.HT and ˆ θ B.HT . chosen as a non-informative prior. In addition, non-parametric priors like kernel density were appliedduring the sensitivity analysis study. Kernel density estimation depends on several variables, includingsample size, bandwidth, and kernel function. In this study, the optimal bandwidth ( h ∗ ) and kernelfunction were chosen such that the asymptotic mean integrated squared error (AMISE) is minimized.The simplified analytical form of AMISE, Eq. (3.3.3.3), is given by: AM ISE (cid:16) ˆ f ( β ) (cid:17) = C ( K ) n · h + ( 14 · h · k · R (cid:16) f (2) ( β ) (cid:17) ) (4.3.1)Where: • C(K)= (cid:82) ( K ( u )) du . • n: sample size. • h: bandwidth. • k = (cid:82) + ∞−∞ u · K ( u ) du . • f (2) ( β ) is the second derivative of Burr PDF. • R ( f (2) ( β )) = (cid:82) ( f (2) ( β )) dβ . 19ayesian Reliability Alenezi & Tsokos • h ∗ = (cid:20) C ( K ) k · R ( f (2) ( β ) ) (cid:21) / · n − / .The minimum AMISE corresponds to the Epanechnikov kernel function ( K ( u ) = (1 − u ) I | u |≤ ), [6].In addition to the Epanechnikov kernel function, the Gaussian kernel function ( K ( u ) = √ π exp (cid:16) − u (cid:17) I I R )was also used in the calculation since it is commonly used for its analytical tractability.Numerical integration techniques were used to compute the Bayesian estimates of β under the H-Tloss function, according to the equations in Section , for each of the five densities and three distinctvalues of θ . Samples of sizes of , , , , , , , , , , and were generated, wherethe parameter θ was assumed to be the MLE of θ ( . ) using Crow’s data. The results, for , repetitions, are presented by Table 9 and Figure 9. Table 9: Averages of the Bayesian estimates (using the subject prior PDFs), under H-T loss functions,and MLEs of the parameter β over 10,000 repetitions n β ˆ β MLE ˆ β B.HT ˆ β IG.HT ˆ β J.HT ˆ β KG.HT ˆ β KE.HT
20 0.7054 0.781303 0.674095 0.799535 0.710047 0.675698 0.67569330 0.7054 0.753465 0.688951 0.762536 0.707341 0.693849 0.69388140 0.7054 0.740665 0.695087 0.745074 0.70651 0.698597 0.69863150 0.7054 0.732738 0.697886 0.734888 0.705885 0.699751 0.6997660 0.7054 0.728669 0.699823 0.728731 0.705847 0.700727 0.70069970 0.7054 0.725499 0.70101 0.724471 0.705776 0.70148 0.70140580 0.7054 0.723539 0.701979 0.721552 0.705882 0.702323 0.7022100 0.7054 0.719621 0.70285 0.717324 0.705664 0.70338 0.703159120 0.7054 0.717465 0.703468 0.71478 0.705632 0.704258 0.703959140 0.7054 0.716224 0.703954 0.713141 0.705692 0.704928 0.70457160 0.7054 0.714739 0.704188 0.711865 0.705629 0.7053 0.704893It can be observed that the Bayesian estimate of β for the Jeffreys prior PDF under the H-T lossfunction converges to the true value faster than other prior PDFs; followed very closely by the Bayesianestimates using the kernel prior PDFs, which tend to converge faster than the inverted gamma priorPDF and are superior to the maximum likelihood approach.Figure 9 presents the graphical behavior of the MLE and Bayesian estimates represented by theiraverage values, where the average values of the Jeffery Bayesian estimates were the closest to the truevalue of β for all sample sizes. Followed by the average values of the kernel Bayesian and the BurrBayesian estimates, which they tend to have similar behavior in estimating for sample sizes of 50 andlarger. Bayesian estimates seems to have insignificant difference for sample sizes of 100 and largerexcept the inverted gamma Bayesian estimate which performed as poor as the MLE in estimating thekey parameter β .Table 10 shows the MSE of the Bayesian estimates of β under the H-T loss function and the subject20ayesian Reliability Alenezi & Tsokos Averagesof β estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □● Fixed ■ MLE ◆ Burr ▲ I. Gamma ▼ Jeffreys ○ Kernel.G □ Kernel.E
20 40 60 80 100 120 140 1600.680.700.720.740.760.780.80
Sample size
Figure 9: Averages of β estimates for different sample sizes. prior PDFs with respect to various sample sizes. Table 10: MSE of the Bayesian estimates, under the H-T loss functions, and MLEs of the parameter β averaged over1000 repetitions for different priors. n MSE of ˆ β MLE
MSE of ˆ β B.HT
MSE of ˆ β IG.HT
MSE of ˆ β J.HT
MSE of ˆ β KG.HT
MSE of ˆ β KE.HT
20 0.0417095 0.00362378 0.0115968
For a sample size of 20, the Jeffery Bayesian estimate of β is the best estimate based on the MSE,followed by the Burr and kernel Bayesian estimates (Table 10). The MSE of the Gaussian KernelBayesian esimate was the lowest for moderate to large sample sizes ( n = 70 , .., ). The invertedgamma Bayesian estimate had the lowest performance compared to other Bayesian estimates. TheMSE of the MLE of the key parameter β indicated its poor performance compared to the Bayesianestimates.Figure 10 shows the MLE and Bayesian estimates of β . It can be observed the Bayesian estimatesunder the H-T loss function are superior to the MLE of the key parameter β for all sample sizesconsidered in this study. Moreover, the Bayesian estimates are presented without MLE in Figure 11to look closely at the performance of Bayesian estimates, under the H-T loss function and the subject21ayesian Reliability Alenezi & Tsokos MSEsof β estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○● MLE ■ Burr ◆ I. Gamma ▲ Jeffreys ▼ Kernel.G ○ Kernel.E
20 40 60 80 100 120 140 1600.010.020.030.04
Sample size
Figure 10: MSEs of β estimates for different sample sizes. prior PDFs, based on their MSEs. MSEsof β estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼● Burr ■ I. Gamma ◆ Jeffreys ▲ Kernel.G ▼ Kernel.E
20 40 60 80 100 120 140 1600.0020.0040.0060.0080.0100.012
Sample size
Figure 11: MSEs of β Bayesian estimates for different sample sizes.
The MSEs of the Jefferys prior, for various sample sizes, were the smallest, followed by MSEs of thekernels and Burr prior PDFs. Bayesian estimates using Gaussian and Epanechnikov kernel probabilitydensities as priors performed similarly, whereas the Bayesian estimate using the inverted gamma PDFhad the lowest performance compared to other Bayesian estimates. It also shows that for a samplesize of n = 30 and larger, the Bayesian estimates using Burr, Jeffery, and kernel PDFs are similar intheir convergence to the true value. For a sample size of n = 70 and larger, all Bayesian estimatestend to converge to the true value of the key parameter β . The Bayesian estimates of β under Jefferys22ayesian Reliability Alenezi & Tsokosand inverted gamma prior PDFs tend to overestimate β , whereas Burr and kernel prior PDFs tend tounderestimate β .For each sample of size 40 based on Monte Carlo simulation, the Bayesian estimates and MLEs of β were calculated when θ ∈ { . , . , } . The comparison is based on the MSE averaged over the , simulated samples. The results are given by Table 11. Table 11: MSE of β Bayesian estimates with Burr, Jeffreys, inverted gamma, and kernel PDFs as priors under the H-Tloss function. MSE of MLE estimate of the parameter β in an NHPP with samples of n = 40 and different values of theparameter θ . θ MSE of ˆ β MLE
MSE of ˆ β B.HT
MSE of ˆ β IG.HT
MSE of ˆ β J.HT
MSE of ˆ β KG.HT
MSE of ˆ β KE.HT
It can be observed that all Bayesian estimates are superior to MLE ( ˆ β ) in estimating β , with samplesize n = 40, while maintaining a consistent behavior for the different values of θ (Table 11). For smallvalue of θ , Jeffrey Bayesian estimate had the lowest MSE value, whereas the MSE of the Gaussiankernel Bayesian estimate was the lowest for moderate and large values of θ .For the case in which we misleadingly assumed the true probability distribution of the key parameter β , we found that the Jeffreys and kernel Bayesian estimates of β had the best performance comparedto the inverted gamma Bayesian estimate of β , indicating that the Bayesian estimate of β is sensitiveto the choice of its prior PDF.As expected, the adjusted MLE of θ produced a better estimate under the mentioned Bayesianinfluence with respect to the H-T loss function. The average values of the MLE and Bayesian estimatesof θ using the MLE and Bayesian estimates of β for various sample sizes, respectively, are presentedby Figure 12. where the average values of the Jeffery Bayesian estimates were the closest to the truevalue of θ for all sample sizes. Followed by the average values of the kernel Bayesian and the BurrBayesian estimates, which tend to have similar behavior in estimating for sample sizes of 40 and larger.When considering sample sizes of 100 and larger, the Bayesian estimates seems to have insignificantdifference except for the inverted gamma Bayesian estimate which which still performs better than theMLE in estimating the parameter θ .The MSE of θ estimates using MLE and Bayesian estimates of β , namely ˆ θ MLE , ˆ θ B.HT , ˆ θ IG.HT , ˆ θ J.HT , ˆ θ KG.HT , and ˆ θ KE.HT , are shown in Table 12 with various sample sizes. For a small sample,moderate, and large sample sizes, n = 20 , , respectively, the Jeffrey Bayesian estimate of θ ( ˆ β J.HT ) performed better than the other estimates, followed by the kernel Bayesian and Burr Bayesianestimates. 23ayesian Reliability Alenezi & Tsokos
Averagesof θ estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○□ □ □ □ □ □ □ □ □ □ □● Fixed ■ MLE ◆ Burr ▲ I. Gamma ▼ Jeffreys ○ Kernel.G □ Kernel.E
20 40 60 80 100 120 140 1601.52.02.53.03.54.0
Sample size
Figure 12: Averages of MLEs of θ using the Bayesian estimates of theparameter β with respect to different priors. Table 12: MSE of the Bayesian estimates, under the H-T loss functions, and MLEs of the parameter θ averaged over1000 repetitions for different priors. n MSE of ˆ θ MLE
MSE of ˆ θ B.HT
MSE of ˆ θ IG.HT
MSE of ˆ θ J.HT
MSE of ˆ θ KG.HT
MSE of ˆ θ KE.HT
20 11.5735 0.175991 1.0789
While the inverted gamma Bayesian estimate of θ outperformed that of the MLE, it did not perform aswell as the other Bayesian estimates based on their MSE values across all sample size. This indicatesthe MLE of θ had the poorest performance compared to the Bayesian estimates.The MSEs of θ estimates using the MLE and Bayesian estimates of the parameter β with respectto different priors is displayed in Figure 13, from which it can be seen that the MLE of θ was extrmelyweak estimator since it has the largest MSEs across the sample sizes. The adjusted θ estimates werewere displayed without the MLE estimates by Figure 14.It can be noted that the θ estimate using the inverted gamma Bayesian estimate of β had the lowestperformance compared to other estimates that used Bayesian estimates of β . In addition, above thesample size n = 40 , the estimate of θ using Burr kernels Bayesian estimates of β performed similarlyin estimating the true value of θ . All θ estimates using Bayesian estimates of β tend to converge tothe true value in similar trajectories, whereas the θ estimate using the Jeffrey Bayesian estimate of β MSEsof θ estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ○ ○ ○ ○● MLE ■ Burr ◆ I. Gamma ▲ Jeffreys ▼ Kernel.G ○ Kernel.E
20 40 60 80 100 120 140 16024681012
Sample size
Figure 13: MSE of the MLEs of θ using MLE and Bayesian estimatesof β with respect to different prior β MSEsof θ estimates ● ● ● ● ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼● Burr ■ I. Gamma ◆ Jeffreys ▲ Kernel.G ▼ Kernel.E
20 40 60 80 100 120 140 1600.20.40.60.81.0
Sample size
Figure 14: MSE of the MLEs of θ using Bayesian estimates of β withrespect to different prior β . converges slightly faster.Therefore, the Bayesian analysis to the PLP under the H-T loss function is sensitive to the priorselection. Based on this finding, an engineer, for example, is recommended to use the Jeffreys or kernelPDFs if he/she lacks any prior knowledge of β . In the present study, we developed the analytical Bayesian form of the key parameter β , under theH-T loss function, in the intensity function, where the underlying failure distribution is the PLP that25ayesian Reliability Alenezi & Tsokosis used for software reliability assessment, among others. The reliability function of the subject modelis written analytically as a function of the intensity function.The behavior of β is characterized by the Burr type XII probability distribution. Real data andnumerical simulation were used to illustrate the efficiency improvement in the estimation of the in-tensity function of PLP under the H-T loss function ( ˆ V B.HT ( t ) ). Based on the , samples ofsoftware failure times, using Monte Carlo simulations and sample size of 40, we found that the Bayesianestimate of β under the H-T loss function ( ˆ β B.HT ) performed better than the MLE of β with respectto three different values of θ ( 0.5, 1.7441, 4 ). Even for different sample sizes ( , , , , , , , , , , and ), similar results were achieved using β =0.7054, θ =1.7441, and averagedover , samples of software failure times.Because the MLE of the second parameter ( θ ) in the intensity function depends on the estimateof β , the adjusted estimate of θ ˆ β B.HT , as expected, performed better than MLE of θ . Moreover,by comparing the relative efficiency metric, mainly using MLEs for both β and θ ( ˆ V MLE ( t ) ), usingBayesian estimates of β under the H-T loss function, and Bayesian MLE of θ ( ˆ V B.HT ( t ) ), we foundthat ˆ V B.HT ( t ) is more efficient in estimating the intensity function V ( t ; β, θ ) .In section , we answered the second research question: Is the Bayesian estimate of the keyparameter, β , using the H-T loss function in the PLP, sensitive to the selection of the prior, whetherparametric or non-parametric? The parametric priors were Burr, Jefferys, and inverted gamma proba-bility distributions. The non-parametric priors were the Gaussian and Epanechnikov kernel densities.The priors’ parameters were estimated using Crow failure times. Additionally, the optimal bandwidthand kernel functions were selected to minimize the asymptotic mean integrated squared error.Using the proposed algorithm, the Bayesian estimate of β under the H-T loss function and BurrPDF as a prior, ˆ β B.HT , performed slightly better than the other Bayesian estimates using differentprior PDFs for small value of θ . The Bayesian estimate of β under H-T loss function and Gaussiankernel PDF as a prior, ˆ β KG.HT , had the smallest MSE comparing to other prior PDFs for moderateand large values of θ . MLE of β continued the poor performance comparing to Bayesian estimatesusing the subject prior PDFs. The Bayesian estimate of β under H-T loss function and inverted gammaPDF as a prior performed the lowest among other Bayesian estimates using Burr, Jeffrey, Gaussian,and Epanechnikov PDFs as priors whereas the latter have slightly different MSEs values. Even fordifferent sample sizes, the MLE of β has the poorest performance comparing to Bayesian estimatesusing different prior PDFs. For small to moderate sample sizes ( n = 20 , , , ), the Bayesianestimate of β under H-T loss function and Jeffrey PDF as a prior, ˆ β J.HT , has the smallest MSE value,26ayesian Reliability Alenezi & Tsokosfollowed closely by the Bayesian estimates using the Gaussian and Epanechnikov kernels, and BurrPDFs as priors. For moderate to large sample sizes ( n = 70 , , .., ), the Bayesian estimate of β under H-T loss function and Gaussian kernel PDF as a prior, ˆ β KG.HT , has the smallest MSE value,followed by the Byesian estimates using the Epanechnikov kernel, Jeffrey, and Burr PDFs as priors.The Bayesian estimates under H-T loss function and the parametric and non-parametric priors wereused the compute the adjusted estimate of θ , where the adjusted θ estimate using the Jeffery Bayesianestimate of β performed the superior estimate comparing to MLE and other Bayesian estimates forvarious sample sizes. Followed by Bayesian estimates using Epanechnikov kernel, Burr, and Gaussiankernel PDFs as priors for moderate to large sample sizes whereas using the Gaussian kernel Bayesianof β to compute the adjusted estimate of θ performed slightly better comparing to the usage of BurrBayesian estimate of β .Thus, based on this aspect of our analysis, we can conclude that the Bayesian analysis approachunder Higgins-Tsokos loss function is superior to the maximum likelihood approach in estimating thereliability function of the Power Law Process. Therefore, the results of this study have the potentialto contribute not only to the reliability analysis field but also to other fields that employ the PowerLaw Process. References [1] J. T. Duane, Learning curve approach to reliability monitoring,
IEEE Trans. Aerospace (1964)563–566.[2] L. H. Crow, Tracking reliability growth, in Proc. Twentieth Conf. Design of Experiments , Report75-2 (US Army Research Office, Research Triangle Park, North Carolina, 1975), pp. 741–754.[3] L. H. Crow, Reliability analysis for complex, repairable systems, in
Reliability and Biometry, eds.F. Proschan and R. J. Serfiing (SIAM, Philadelphia, Pennsylvania, 1974). [4] H. E. Ascher and H. Feingold, Repairable Systems Reliability: Modeling, Inference,
Misconceptionsand Their Causes (Marcel-Dekker, New York, 1984).[5] C. A. Moltnares and C. P. Tsokos, Bayesian Reliability Approach To The Power Law ProcessWith Sensitive Analysis To Prior Selection,
International Journal of Reliability, Quality and SafetyEngineering (World Scientific, 2013). 27ayesian Reliability Alenezi & Tsokos[6] F. N. Alenezi and C. P. Tsokos, The Effectiveness of The Squared Error And Higgins-Tsokos LossFunctions On The Bayesian Reliability Analysis Of Software Failure Times Under The Power LawProcess,
Engineering 11, 272- 299 (Scientific Research Publishing, 2019).[7] J. J. Higgins and C. P. Tsokos, A study of the effect of the loss function on Bayes estimates offailure intensity, MTBF, and Reliability,
Applied Mathematics and Computation (1980) 145–166.[8] Harold Jeffreys, F.R.S., An invariant form for the prior probability in estimation problems, RoyalSociety (1946).[9] D. J. Daley and D. Vere-Jones, An introduction to the theory of point processes,
Vol. I. Probabilityand its Applications (Springer, New York, second edition, 2003).[10] S. Yamada, M. u Ohba and S. Osaki, S-Shaped Reliability Growth Modeling for Software ErrorDetection,
IEEE Trans. Reliab.
R-32 (1983) 475–484.[11] A. L. Goel K. Okumoto, Time-Dependent Error-Detection Rate Model for Software Reliabilityand Other Performance Measures,
IEEE Trans. Reliab.
R-28 (1979) 206–211.[12] M. Engelhardt and J. L. Bain, Statistical analysis of a compound power law model for repairablesystems,
IEEE Trans. Reliab.
R-36 (1987) 392–396.[13] M. Kimura, T. Toyota and S. Yamada, Economic analysis of software release problems withwarranty cost and reliability requirement,
Reliability Engineering and System Safety (1999)49–55.[14] M. Engelhardt and J. L. Bain, Prediction intervals for the Weibull process, Technometrics (1978) 167–169.[15] S. E. Rigdon and A. P. Basu, The effect of assuming a homogeneous Poisson process when thetrue process is a power law process, J. Qual. Tech. (1990) 111–117.[16] S. E. Rigdon and A. P. Basu, The power law process: A model for the reliability of repairablesystems, J. Qual. Tech. (1990) 251–260.[17] S. E. Rigdon and A. P. Basu, Estimating the intensity function of a power law process at thecurrent time: Time truncated case, Commun. Stat.-Simulat. Comput. (1990) 1079–1104.[18] H. E. Ascher, T. T. Y. Lin and D. P. Siewiorek, Modification of error log analysis: Statisticalmodeling and heuristic trend analysis, IEEE Trans. Reliab. (1992) 599–601.28ayesian Reliability Alenezi & Tsokos[19] J. Kyparisis and N. D. Singpurwalla, Bayesian inference for the Weibull process with applicationsto assessing software reliability growth and predicting software failures, in Computer Science andStatistics: Proc. Sixteenth Symp. Interface (1985).[20] W. M. Bassin, Increasing hazard functions and overhaul policy, in
Proc. 1969 Annual Symp.Reliability , Chicago, Vol. (IEEE, 1969), pp. 173–178.[21] Y. Xu and C. P. Tsokos, Non-homogeneous Poisson process for evaluating stage I and II ductalbreast cancer treatment, Journal of Modern Applied Statistical Methods , 10 (2), pp. 646–655, (2011)[22] M. Evans, N. Hastings and B. Peacock, Statistical Distributions (John Wiley and Sons, Inc., NewYork, 2000).[23] C. P. Tsokos, Recent Advances in Life-Testing and Reliability:
A Volume in Honor of AlonsoClifford Cohenm Jr. , ed. N. Balakrishnan (CRC Press, Boca Ratón, 1995).[24] L. J. Bain and M. Engelhardt,
Statistical Analysis of Reliability and Life Testing Models (Marcel-Dekker, New York, 1991).[25] S. K. Bar-Lev, I. Lavi and B. Reiser, Bayesian inference for the power law process,
Ann. Inst.Stat. Math. (1992) 623–639.[26] V. A. Camara and C. P. Tsokos, Bayesian reliability modeling with a new loss function andnumerical simulation, Statistica (2001) 619–630.[27] J. M. Finkelstein, Confidence bounds on the parameters of the Weibull process, Technometrics (1978) 115–117.[28] M. S. Hamada, A. G. Wilson, C. Shane Reese and H. F. Martz, Bayesian Reliability (SpringerScience, New York, 2008).[29] H. Kim, S. Choi and S. Kim, Bayesian inference for the power law process with the power prior,
J. Korean Statist. Soc. (2005) 331–344.[30] K. Kim, S. W. Kim and H. Kim, Intrinsic Bayes factors for model selection in nonhomogeneousPoisson process, Far East J. Theor. Stat. (2003) 15–30.[31] S. W. Kim and D. Sun, Intrinsic priors for model selection using an encompassing model withapplications to censored failure time data,
Lifetime Data Anal. (2000) 251–269.29ayesian Reliability Alenezi & Tsokos[32] J. J. Higgins and C. P. Tsokos, A quasi-Bayes estimate of the failure intensity of a reliability-growthmodel, IEEE Trans. Reliab.
R-30 (1981) 471–475.[33] L. Lee and S. K. Lee, Some results on inference for theWeibull process,
Technometrics (1978)41–45.[34] R. T. Lingham and S. Sivaganesan, Testing hypotheses about the power law process under failuretruncation using intrinsic Bayes factors, Ann. Inst. Stat. Math. (1997) 693–710.[35] H. Raiffa and R. Schlaifer, Applied Statistical Decision Theory (Harvard Graduate School ofBusiness Administration, Boston, 1961).[36] S. K. Bhattacharya, Bayesian approach to life testing and reliability estimation,
J. Amer. Stat.Assoc. (1967) 48–62.[37] W. A. Thompson Jr., On the foundations of reliability, Technometrics (1981) 1–23.[38] C. P. Tsokos and A. N. V. Rao, Estimation of failure intensity for the Weibull process, Reliab.Eng. Syst. Safety45