Software reliability growth models with normal failure time distributions

Reliability Engineering and System Safety 116 (2013) 135–141

Contents lists available at SciVerse ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Software reliability growth models with normal failure time distributions Hiroyuki Okamura a,n, Tadashi Dohi a, Shunji Osaki b a b

Department of Information Engineering, Graduate School of Engineering, Hiroshima University, 1–4–1 Kagamiyama, Higashi Hiroshima 739–8527, Japan Faculty of Information Sciences and Engineering, Nanzan University, Seto 489–0863, Japan

a r t i c l e i n f o

abstract

Article history: Received 9 August 2011 Accepted 5 February 2012 Available online 6 May 2012

This paper proposes software reliability growth models (SRGM) where the software failure time follows a normal distribution. The proposed model is mathematically tractable and has sufficient ability of fitting to the software failure data. In particular, we consider the parameter estimation algorithm for the SRGM with normal distribution. The developed algorithm is based on an EM (expectationmaximization) algorithm and is quite simple for implementation as software application. Numerical experiment is devoted to investigating the fitting ability of the SRGMs with normal distribution through 16 types of failure time data collected in real software projects. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Software reliability growth model Normal distribution EM algorithm

1. Introduction Software reliability growth models (SRGMs) are widely used to estimate quantitative software reliability measures. The SRGMs are generally defined as stochastic counting processes regarding the number of faults detected or failures experienced in testing phase. Based on the counting process, we compute some of software reliability measures such as the quantitative software reliability which is defined as a probability that no failure occurs during a certain time interval. Since the seminal works of Jelinski and Moranda [1] and Goel and Okumoto [2], a number of SRGMs have been competitively proposed to improve their fitting abilities and accuracy of the future failure prediction [3–7]. In particular, non-homogeneous Poisson processes (NHPPs) have much popularity for the software reliability modeling due to their mathematical tractability. All the NHPP-based SRGMs are characterized by mean value functions that indicate the expected numbers of failures with respect to testing time. In fact, Goel and Okumoto [2], Goel [8], Musa and Okumoto [9], Ohba [10,11], Yamada, Ohba and Osaki j [12], Zhao and Xie [13] and Pham [14] proposed NHPP-based SRGMs whose mean value functions draw typical non-linear curves. The quantitative software reliability assessment generally consists of the following steps: (i) Collect software failure data, (ii) Determine a set of candidates of models, (iii) Estimate parameters for all the candidates, (iv) Choose the best model under a certain criterion, and (v) Calculate the quantitative measures using the best model. In the step (i), we collect the

n

Corresponding author. E-mail addresses: [email protected] (H. Okamura), [email protected] (T. Dohi), [email protected] (S. Osaki). 0951-8320/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ress.2012.02.002

failure time data in software testing or the number of failures per working day. The steps (iii) and (iv) can be solved by statistical arguments such as the maximum likelihood (ML) estimation and the theory of information criterion. In the step (v), the software reliability measures can be computed according to their formulas. However, in the step (ii), there is no definitive method to determine a suitable set of candidates. If the number of candidates is large, it causes a large effort in the parameter estimation procedure of the step (iii). At the same time, since there have been over 200 SRGMs proposed in many papers, it is infeasible to estimate the model parameters for all the existing SRGMs. On the other hand, if the set did not include the sufficient number of candidates, it causes the degradation of fitting ability and prediction accuracy. That is, it is important to determine a suitable set of candidates. The modeling framework is helpful to determine a set of candidates, which involves some of the existing NHPP-based SRGMs. The most well-known and well-defined modeling framework has been proposed by Langberg and Singpurwalla [15]. They revealed the statistical structure of a class of NHPP-based SRGMs. According to their approach, NHPP-based SRGMs consists of only two factors; the number of total software faults before testing (the number of total failures experienced in the software lifetime) and failure time caused by each fault in testing. In particular, when the number of total software faults is given by a Poisson random variable, the framework covers the NHPP-based SRGMs whose mean value function is bounded by a finite value. Then the concrete formulas of NHPP-based SRGMs are dominated by only types of the failure time distribution. In their framework, the candidate selection of NHPP-based SRGMs can be reduced to selecting a set of suitable failure time distributions. In fact, typical NHPP-based SRGMs consist of well-known statistical distributions such as exponential [2], gamma [12,13], Pareto [16,17], logistic [10,18] and extreme-value (Weibull, Frechet and Gompertz)

136

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

distributions [8,19]. Generally, in the area of reliability engineering, the exponential, gamma, Pareto, logistic, extreme-value distributions are also applied to representing the hardware failure time distribution. However, the NHPP-based SRGM whose failure time follows the normal distribution has not been fully discussed. The normal distribution is one of the most prominent distributions in statistics. As is well known, the normal distribution can be derived as the probability distribution of the sum of random variables based on the central limit theorem (CLT). In addition, the normal distribution belongs to the exponential family, likewise exponential distribution and gamma distribution, and the exponential family are analytically tractable. In the reliability engineering for hardware, the normal distribution can also be applied to representing failure time. However, in the software reliability engineering, it has been never assumed that software failure time is normally distributed. One of the reasons is why the normal distribution is defined on the range ð1,1Þ, though the failure time distribution is defined on the positive support. Only the work using normal distribution in the software reliability engineering is in [20]; Achcar et al. used the log-normal distribution as software failure time distribution. One goal of this paper is to develop SRGMs with normal distribution for the software reliability assessment. In particular, we deal with two SRGMs based on truncated normal and log-normal distributions. Another goal is to derive effective parameter estimation algorithms for the SRGMs with normal distribution. The estimation of model parameters is needed to quantitatively assess the software reliability from observed failure data, i.e., in the step (ii). The commonly used method is the ML estimation. The ML estimation is to find the maximum of likelihood function of observed software failure data. Since ML estimates (MLEs) have rational properties like asymptotic efficiency, MLEs are expected to be suitable even in the parameter estimation for NHPPbased SRGMs. Although ML estimation enables us to compute statistically proper estimates for model parameters, we occasionally encounter several difficulties for the parameter estimation. In general, MLEs are obtained by maximizing log-likelihood functions or by solving non-linear equation called likelihood equation which is derived form the first derivative of log-likelihood function. However, even for Goel and Okumoto model, we cannot obtain closed forms of MLEs. In other words, we employ numerical approaches to find MLEs based on log-likelihood functions. The commonly used approaches to find MLEs are Newton’s method, quasi-Newton’s method and Nelder–Mead method. Although Newton’s method is a powerful tool to calculate MLEs, it has the local convergence property and may fail to get the solution due to unsuitable initial guesses. The Nelder–Mead method is one of the direct search methods, which is more stable for initial guesses than Newton’s method, but a few design parameters, such as expansion rate, must be manually adjusted before executing the algorithm. These are serious problems when we develop and implement the software reliability assessment tool. This is one of the reasons why the software reliability assessment tool has not been distributed to practical software development projects. Recently, we developed a parameter estimation algorithm based on the EM (expectation–maximization) principle [21,22] and applied it to the software reliability assessment based on NHPP-based SRGMs [23–25,19]. Our key idea here is to regard the underlying software failure data as incomplete data in the modeling framework by Langberg and Singpurwalla [15]. It was shown that the EM algorithm can be applied to typical NHPPbased SRGMs and can give much advantages on global convergence and reduction of computation efforts. The goal of this paper is to figure out the EM algorithms for the SRGMs with normal distribution.

This paper is organized as follows. In Section 2, we introduce the modeling framework by Langberg and Singpurwalla [15] and develop the SRGMs under normally distributed failure times. In Section 3, we derive the concrete EM algorithms for the normal distribution based SRGMs. Section 4 is devoted to examining the performance of them. Finally, we conclude the paper with remarks in Section 5.

2. NHPP-based SRGM with normal distribution 2.1. Modeling framework Let fXðtÞ,t Z0g denote the number of software failures experienced before time t. According to Langberg and Singpurwalla [15], this paper considers the following model assumptions:

 Assumption A: Software failures occur at mutually indepen-



dent random times. The probability distributions of all failure times are identical. The probability density and cumulative distribution functions (p.d.f. and c.d.f.) are given by f(t) and F(t), respectively. Assumption B: The number of inherent software faults causing failures is given by a Poisson random variable.

Under the assumption that the number of inherent faults is fixed as N, the probability mass function (p.m.f.) of the cumulative number of failures experienced by time t is given by   N ð1Þ PðXðtÞ ¼ nÞ ¼ FðtÞn F ðtÞNn , n where F ðÞ ¼ 1FðÞ. Since N is a Poisson random variable with mean o, the cumulative number of software failures before time t has the following p.m.f.: PðXðtÞ ¼ nÞ ¼

ðoFðtÞÞn oFðtÞ e : n!

ð2Þ

Eq. (2) is equivalent to the probability mass function of NHPP with mean value function oFðtÞ. 2.2. SRGM with normal distribution In the modeling framework of Eq. (2), NHPP-based SRGMs are defined as respective failure time distributions. By substituting statistical distributions defined on positive support into F(t), we obtain corresponding NHPP-based SRGMs. For example, if F(t) is given by an exponential c.d.f., then the corresponding SRGM becomes Goel and Okumoto model [2]. This paper applies normal distributions to the failure time distribution. Although the normal distribution is one of the most prominent distributions in statistics, only Achcar et al. [20] used the log-normal distribution as software failure time distribution. There has not been discussed on the applicability of normal distribution to the failure time distribution. Let FðtÞ and fðtÞ be c.d.f. and p.d.f. of the standard normal distribution with mean 0 and variance 1. Then we have Z t FðtÞ ¼ fðuÞ du, ð3Þ 1

and 1



1 2



fðtÞ ¼ pffiffiffiffiffiffi exp  t2 : 2p

ð4Þ

Based on the standard normal distribution, we develop the SRGMs with truncated normal distribution (TNorm) and log-normal distribution (LNorm). Concretely, the failure time distribution of TNorm

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

marginal log-likelihood function (LLF):

is given by Fðt; m, sÞ ¼ F

  tm

h^ ML ¼ argmax Lðh; DÞ,

F ðm=sÞ,

s

1 o m o 1, s 4 0,

ð5Þ

s

The above is the model based on log-normal order statistics introduced by Achcar et al. [20].

3. EM algorithm for SRGMs with normal distribution 3.1. Maximum likelihood estimation In the software reliability assessment, we should estimate model parameters of NHPP-based SRGMs from observed data. The most commonly used technique to parameter estimation is maximum likelihood (ML) estimation. Let DT ¼ ðt 1 , . . . ,t K ; TÞ be a set of failure times experienced by time T. Without loss of generality, we assume 0 o t1 o    ot K . For the observed data DT , the log-likelihood function (LLF) for NHPP-based SRGMs with failure time data is given by Lðo, h; DT Þ ¼ log pðDT ; o, hÞ, pðDT ; o, hÞ ¼ oK

K Y

f ðt k ; hÞexpðoFðT; hÞÞ,

ð7Þ

ð8Þ

k¼1

where h is a parameter vector for the failure time distribution, f ð; hÞ is a density function of Fð; hÞ and pðÞ is an appropriate probability mass or density function. On the other hand, DG ¼ ðn1 , . . . ,nK Þ denotes the number of failures experienced on the time sequence 0 o t1 o t2 o    o tK . Then the LLF with the grouped data is given by Lðo, h; DG Þ ¼ log pðDG ; o, hÞ, K Y ðoFðtk ÞoFðtk1 ÞÞnk expðoFðtK ; hÞÞ: nk ! k¼1

ð11Þ

h

where F ðÞ ¼ 1FðÞ. On the other hand, the failure time distribution of LNorm is   log tm , 1 o m o 1, s 4 0: ð6Þ Fðt; m, sÞ ¼ F

pðDG ; o, hÞ ¼

137

ð9Þ

ð10Þ

The ML estimation is to find the parameters maximizing the LLF, so-called maximum likelihood estimates (MLEs). In general, MLEs of NHPP-based SRGMs cannot be expressed as closed forms, even for the simplest model, i.e., Goel and Okumoto model. That is, we need to utilize any iterative methods for numerical optimization such as Newton’s method, quasi-Newton’s method, Fisher’s scoring method and Nelder–Mead method. However, it is well known that it is difficult to choose appropriate initial parameters in these iterative methods. In addition, some methods require appropriate design parameters such as reflection and expansion rates in Nelder–Mead method. This property adversely affects making software reliability assessment tools. In the reliability assessment tool with Newton’s or quasi-Newton’s method, the users should change initial parameters and other design parameters depending on observed failure data. 3.2. EM algorithm The EM algorithm is an iterative method for computing ML estimates with incomplete data [21,22]. Let D and Z be observable and unobservable data vectors, respectively, and we generally estimate a model parameter vector h from only the observable data vector D. In the context of ML estimation, the problem corresponds to finding a parameter vector that maximizes a

Lðh; DÞ ¼ log pðD; hÞ ¼ log

Z

pðD,Z; hÞ dZ,

ð12Þ

where pðÞ is any probability density or mass function. Taking account of the posterior distribution of unobservable data vector with the parameter vector h0 and Jensen’s inequality, we have Z pðD,Z; hÞ dZ Lðh; DÞ ¼ log Z pðD,Z; hÞ pðZ9D; h0 Þ dZ ¼ log pðZ9D; h0 Þ Z pðD,Z; hÞ Z pðZ9D; h0 Þlog dZ pðZ9D; h0 Þ  Zðh; h0 Þ:

ð13Þ

The posterior distribution for unobservable data can be obtained from Bayes theorem: pðZ9D; h0 Þ ¼ R

pðD,Z; h0 Þ : pðD,Z; h0 Þ dZ

ð14Þ

Eq. (13) yields Lðh; DÞZðh; h0 Þ ¼ DKL ðpðZ9D; h0 ÞJpðZ9D; hÞÞ,

ð15Þ

where DKL ðPJQ Þ is the Kullback–Leibler distance from the distribution P to the distribution Q. Hence the difference Lðh; DÞLðh0 ; DÞ is given by Lðh; DÞLðh0 ; DÞ ¼ Zðh; h0 ÞZðh0 ; h0 Þ þ DKL ðpðZ9D; h0 ÞJpðZ9D; hÞÞ: ð16Þ Since DKL ðJÞ Z 0, Eq. (16) implies that the maximization of upper bound results in the maximization of marginal LLF. Let Q ðh9h0 Þ denote the conditional expected LLF with respect to the complete data vector ðD,ZÞ using the posterior distribution for unobservable data vector with provisional parameter vector h0 : Z pðZ9D; h0 Þlog pðD,Z; hÞ dZ: ð17Þ Q ðh9h0 Þ ¼ E½log pðD,Z; hÞ9D; h0  ¼ Then Eq. (13) is rewritten in the form: Z Zðh; h0 Þ ¼ Q ðh9h0 Þ pðZ9D; h0 Þlog pðZ9D; h0 Þ dZ:

ð18Þ

Since the second term of the above equation is constant, the maximization of Q ðh9h0 Þ with respect to h is directly reduced to the maximization of Zðh; h0 Þ with respect to h. Based on the above discussion, the EM algorithm consists of E-step and M-step. E-step computes the conditional expected LLF with respect to the complete data vector ðD,ZÞ using the posterior distribution for unobservable data vector with provisional parameter vector h0 , i.e., Q ðh9h0 Þ. In M-step, we find a new parameter vector h00 that maximizes the expected LLF:

h00 :¼ argmax Q ðh9h0 Þ,

ð19Þ

h

and h00 becomes a provisional parameter vector at the next E- and M-steps. These steps surely increase the marginal LLF. The E- and M-steps are repeatedly executed until the parameters converge. 3.3. EM-step formulas for LNorm with failure time data Consider the EM algorithm for LNorm: PðNðtÞ ¼ xÞ ¼

ðoFðt; m, sÞÞx expðoFðt; m, sÞÞ, x!

ð20Þ

138

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

Fðt; m, sÞ ¼ F

  log tm

s

:

ð21Þ

In the modeling framework by Langberg and Singpurwalla [15], failure times after time T are not observable. The complete data are defined by 0 o T 1 o T 2 o    oT N , where N is the number of total faults and Ti is the i-th ordered failure time. That is, unobserved data are Z T ¼ ðT K þ 1 , . . . ,T N ,NÞ. The complete LLF is given by log pðDT ,Z T ; o, hÞ ¼ N log oo þ

N X

where hðÞ is a measurable function. Applying the above formula, we have  " # Z 1  N K X X  E log T k DT ; o, m, s ¼ log t k þ o log t dFðt; m, sÞ  T k¼1

ð22Þ

i¼1

K X

¼

log t k þ oF

ð1Þ

  log Tm ,

s

k¼1

" log f ðT i ; m, sÞ:

k¼1

E

N X

k¼1

 # Z 1  K X  ðlog T k Þ2 DT ; o, m, s ¼ ðlog t k Þ2 þ o ðlog tÞ2 dFðt; m, sÞ  T k¼1

From the standard argument of MLEs, the MLEs of o and h can be derived as

¼

K X

ðlog t k Þ2 þ oF

ð2Þ

  log Tm ,

k¼1

o^ ¼ N,

ð34Þ

s

ð35Þ

ð23Þ where

N 1X m^ ¼ log T k , Nk¼1

ð24Þ

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X 2 s^ ¼ t ðlog T k Þ2 m^ : Nk¼1

ð2Þ

ð25Þ

o :¼ E½N9DT ; o0 , m0 , s0 , 1

o

E½

XN k¼1

ð26Þ

  log T k DT ; o0 , m0 , s0 ,

ð27Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  1 XN  s :¼ E½ k ¼ 1 ðlog T k Þ2 DT ; o0 , m0 , s0 m2 ,

o

ð28Þ

where o0 , m0 and s0 are provisional parameters. On the other hand, by applying Bayes theorem, we have pðN9DT ; o, m, sÞpoK eo

K Y

F ðzÞ ¼

f ðtk ; m, sÞ,

N¼K

ð29Þ

f ðtk ; m, sÞ

k¼1

Z

Z

1

1



 T

N Y

T N1 k ¼ K

ð30Þ

Then the posterior distribution of N becomes the Poisson p.m.f. with mean oF ðT; m, sÞ for N Z K: ðoF ðT; m, sÞÞNK expðoF ðT; m, sÞÞ: ðNKÞ!

ð31Þ

Therefore, the EM-step formula for o can be obtained as

o :¼ K þ o0 F ðT; m0 , s0 Þ:

ð32Þ

On the other hand, we provide the following formula for the expectation on T 1 , . . . ,T N with a general failure time distribution Fðt; hÞ [24]:  " # Z 1  N K X X  E hðT k ÞDT ; o, h ¼ hðt k Þ þ o hðuÞ dFðu; hÞ, ð33Þ  T k¼1

k¼1

1

s

1

ufðuÞ du ¼ sfðzÞ þ mF ðzÞ,

ð36Þ

sz þ m

Z

1

u2 fðuÞ du

sz þ m

ð37Þ

In the case where the data is given by the grouped form, the exact failure times are also regarded as unobserved data. Then we can obtain the similar formulas of the M-step in the case of grouped data as follows:

o :¼ E½N9DG ; o0 , m0 , s0 , m :¼

1

o

" E

N X

k¼1

ð38Þ

 #   log T k DG ; o0 , m0 , s0 , 

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  " # u  N X u1  s :¼ t E ðlog T k Þ2 DG ; o0 , m0 , s0 m2 :  o

ð39Þ

ð40Þ

k¼1

ð41Þ

and

f ðt k ; m, sÞ dT K þ 1    dT N ,

N ZK þ 1:

pðN9DT ; o, m, sÞ ¼

s

Bayes theorem gives us the following formulas: PK ðoF ðtK ; m, sÞÞN k ¼ 1 nk pðN9DT ; o, m, sÞ ¼ expðoF ðtK ; m, sÞÞ, P ðN Kk ¼ 1 nk Þ!

and K Y

Z

3.4. EM-step formulas for LNorm with grouped data

k¼1

pðN9DT ; o, m, sÞpoN eo

1

¼ ðs2 z þ 2msÞfðzÞ þ ðs2 þ m2 ÞF ðzÞ:

Eq. (22) implies that the estimation problem of NHPP-based SRGMs under the complete data can be decomposed into separate data fitting problems for the Poisson distribution and failure time distribution with independent and identically distributed (IID) samples. From Eq. (19), we have the following M-step formulas:

m :¼

ð1Þ

F ðzÞ ¼

o :¼

K X

nk þ o0 F ðtK ; m0 , s0 Þ:

ð42Þ

k¼1

The expectation on T 1 , . . . ,T N with a general failure time distribution Fðt; hÞ in the case of grouped data is also given by [24]:  " # R tk  N K n X X k tk1 hðuÞ dFðu; hÞ  hðT k ÞDG ; o, h ¼ E  Fðtk ; hÞFðtk1 ; hÞ k¼1 k¼1 Z 1 hðuÞ dFðu; hÞ, ð43Þ þo tK

where hðÞ is a measurable function. Using the above formula, we have  " # ð1Þ ð1Þ  N K X X F ðzk1 ÞF ðzk Þ ð1Þ  E þ oF ðzK Þ, log T k DG ; o, m, s ¼ nk  F ðz Þ F ðz Þ k1 k k¼1 k¼1 ð44Þ

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

" E

N X k¼1

 # ð2Þ ð2Þ  K X F ðzk1 ÞF ðzk Þ ð2Þ ðlog T k Þ DT ; o, m, s ¼ nk þ oF ðzK Þ,  F ðzk1 ÞF ðzk Þ 2

k¼1

data. Then we have the following M-step formulas: ~ G; o ~ 0 , m0 , s0 , o~ :¼ E½N þ N9D

ð45Þ 1

Consider the EM algorithm for TNorm

Fðt; m, sÞ ¼

  ðoFðt; m, sÞÞx exp oFðt; m, sÞ , x!

FððtmÞ=sÞ : F ðm=sÞ

ð46Þ

Unlike the EM algorithm for LNorm, samples generated before time 0 are also unobserved. Let T N~ o    oT 1 o 0 be unobserved failure times before time 0, where N~ is the number of failures already experienced by t ¼0. Then unobserved data can be ~ Under the rewritten by Z T ¼ ðT N~ , . . . ,T 1 ,T K þ 1 , . . . ,T N , N,NÞ. assumption, we obtain the M-step formulas

1

"

m :¼ ~ E o

 #  N X  0 0 0 ~ ,m ,s , T k DT ; o  ~

ð48Þ

ð49Þ

ð51Þ

T

It is worth noting that the above equations include a lefttruncated term as well as a right-truncated term. According to the above formula, we have the expected values in the M-step formulas as follows:  " #     N K m

X X ð1Þ ð1Þ Tm  ~ , m, s ¼ ~ mF , T k DT ; o tk þ o  þF E  s s ~ k¼1

ð53Þ "

N X k ¼ N~

~ G; o ~ , m, s ¼ E½N þ N9D

K X

~ ð1F ðz0 ÞÞ þ o ~ F ðzK Þ, nk þ o

"

N X

ð59Þ

 # ð1Þ ð1Þ  K X F ðzk1 ÞF ðzk Þ  ~ , m, s ¼ T k DG ; o nk  F ðzk1 ÞF ðzk Þ k¼1

~ ðmF þo

ð50Þ

Also, for any measurable function hðÞ and a general failure time distribution defined on the range ð1,1Þ, the following equation holds [19,23]:  " # Z 0  N K X X  ~ ~ hðT k ÞDT ; o , h ¼ hðt k Þ þ o hðuÞ dFðu; hÞ E  1 k¼1 k ¼ N~ Z 1 ~ hðuÞ dFðu; hÞ: ð52Þ þo

E

tK

1

According to the above formula, we have the following E-step formulas in the case of TNorm with grouped data:

k ¼ N~

~ is a Poisson parameter for the total number of failures where o ~ into the mean on the range ð1,1Þ. Then we should change o number of failures experienced on the positive support ½0,1Þ after ~ F ð0Þ. the EM algorithm is finished, i.e., o :¼ o Similar to the posterior distribution of N, we have

k ¼ N

For any measurable function hðÞ and a general failure time distribution defined on the range ð1,1Þ, we derive the following equation:  " # R tk  N K n X X k tk1 hðuÞ dFðu; hÞ  ~ ,h ¼ hðT k ÞDG ; o E  Fðtk ; hÞFðtk1 ; hÞ k¼1 k ¼ N~ Z 0 Z 1 ~ ~ hðuÞ dFðu; hÞ þ o hðuÞ dFðu; hÞ: ð58Þ þo

E

k ¼ N

~ T; o ~ 0 , m0 , s0  ¼ K þ o ~ 0 F ðT; m0 , s0 Þ þ o ~ 0 Fð0; m0 , s0 Þ: E½N þ N9D

ð57Þ

k¼1

k ¼ N

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  " # u  N X u1 0 0 0 2 t 2 ~ s :¼ ~ E T k DT ; o , m , s m ,  o ~

ð56Þ

k ¼ N

k ¼ N

ð47Þ

~ T; o ~ 0 , m0 , s0 , o~ :¼ E½N þ N9D

ð55Þ

 #  N X  0 0 0 ~ T k DG ; o , m , s ,  ~

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  " # u  N X u1  ~ 0 , m0 , s0 m2 : s :¼ t ~ E T 2k DG ; o  o ~

3.5. EM-step formulas for TNorm with failure time data

PðNðtÞ ¼ xÞ ¼

"

m :¼ ~ E o

where zk ¼ ðlog tk mÞ=s.

139

 #     K m

X ð2Þ ð2Þ Tm  ~ , m, s ¼ ~ s2 þ m2 F : þF T 2k DT ; o t 2k þ o   s s k¼1

ð54Þ

Similar to the case of EM-step formulas for LNorm with grouped data, the exact failure times are regarded as unobserved

~F ðz0 ÞÞ þ o

ð1Þ

ðzK Þ,

ð60Þ

 " # ð2Þ ð2Þ  N K X X F ðzk1 ÞF ðzk Þ 2 ~ , m, s ¼ T k DG ; o nk E  F ðzk1 ÞF ðzk Þ ~ k¼1

k ¼ N

~ ðs2 þ m2 F þo

ð2Þ

~F ðz0 ÞÞ þ o

ð2Þ

ðzK Þ,

ð61Þ

where zk ¼ ðt k mÞ=s.

4. Numerical experiment We examine the data fitting abilities of LNorm and TNorm through numerical examples. In the experiments, we consider the standard SRGMs based on typical failure time distributions including the normal distribution. Table 1 presents the standard SRGMs and their corresponding failure time distributions. Except for TNorm, they have been proposed in the past literature. Note that all the standard SRGMs have respective stable parameter estimation procedures based on the EM algorithm [24,19,23], Table 1 The standard NHPP-based SRGMs and their failure time distributions. Model abbrv.

Failure time dist.

References

Exp Gamma Pareto TNorm LNorm TLogis LLogis TEvMax

exponential distribution gamma distribution Pareto (type II) distribution truncated normal distribution log-normal distribution truncated logistic distribution log-logistic distribution Gompertz distribution (truncated extreme-value max distribution) Frechet distribution (log-extreme-value max distribution) Gompertz distribution (truncated extreme-value min distribution) Weibull distribution (log-extreme-value min distribution)

[2] [12,13,20] [16,17] – [20] [10,23] [18,23] [26,19]

LEvMax

3.6. EM-step formulas for TNorm with grouped data

ð1Þ

TEvMin LEvMin

[19] [19] [8,19]

140

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

Table 2 Failure data information.

Table 3 AICs for SYS1, SYS2, SYS3, SS2, SS3 and SS4 with failure time data.

Data

Size (instructions)

# of failures

Model

SYS1

SYS2

SYS3

SS2

SS3

SS4

SYS1 SYS2 SYS3 SYS4 SYS5 SYS6 SYS14C SYS17 SYS27 SYS40 SS1A SS1B SS1C SS2 SS3 SS4

21,700 27,700 23,400 33,500 2,445,000 5700 (Hundreds 61,900 126,100 180,000 (Hundreds (Hundreds (Hundreds (Hundreds (Hundreds (Hundreds

136 54 38 53 831 73 36 38 41 101 112 375 277 192 278 196

Exp Gamma Pareto TNorm LNorm TLogis LLogis TEvMax LEvMax TEvMin LEvMin

1954.7 1940.2 1944.3 1961.3 1942.6 1958.3 1940.5 1957.7 1943.9 1963.8 1940.2

903.5 901.0 900.1 907.4 899.6 906.1 900.4 905.9 899.5 908.4 900.7

612.2 604.1 600.7 617.9 600.6 614.9 601.6 614.6 600.3 620.0 602.8

5232.0 5232.2 5235.3 5219.5 5242.0 5219.2 5232.1 5223.1 5254.3 5216.2 5232.3

7313.0 7313.3 7315.3 7313.0 7319.5 7312.9 7314.9 7313.3 7325.0 7313.1 7315.2

5269.7 5270.0 5271.9 5267.2 5273.0 5267.6 5270.3 5268.1 5275.9 5266.3 5273.3

of thousands)

of of of of of of

thousands) thousands) thousands) thousands) thousands) thousands)

Table 4 BICs for SYS1, SYS2, SYS3, SS2, SS3 and SS4 with failure time data.

where the EM algorithms for LNorm and TNorm are presented in the previous section of this paper. In addition, we use benchmark software failure data obtained from DACS (Data & Analysis Center for Software) that opens 16 types of software failure data to the public through the web site.1 Table 2 gives summaries of the software failure data. Furthermore, each failure data set was two data forms; failure time data and grouped data. The failure time data consist of time intervals of failure occurrences which were recorded with wall-clock time. On the other hand, the grouped data are the number of failures per working day. To measure the goodness-of-fit of the SRGMs, this paper uses two information criteria: AIC (Akaike’s information criterion) [27] and BIC (Bayesian information criterion) [28] defined as follows: AIC ¼ 2ðmaximum of log  likelihoodÞ þ 2p,

ð62Þ

BIC ¼ 2ðmaximum of log  likelihoodÞ þ p log m,

ð63Þ

where p is the number of free parameters and m is the number of samples. Both AIC and BIC indicate how well the model is fitted to observed samples, and the model with the small AIC or BIC is better than those with the high AIC or BIC. That is, for each set of failure data, we decide the best model by relatively comparing AICs or BICs of the estimated models. In the statical sense, AIC estimates the bias of estimator in the context of ML estimation, and BIC is an approximate value of the (unnormalized) posterior probabilities of mixture parameters for the mixture model of all the candidates. Then we say that there is a significance difference between two models if the difference of the two AICs or BICs is greater than 1. Tables 3 and 4 present AICs and BICs for the standard SRGMs in SYS1, SYS2, SYS3, SS2, SS3 and SS4, when the model parameters are estimated from failure time data. In SYS2 and SYS3, LEvMax can be selected as the best model in terms of AIC and BIC. However, the differences between LEvMax and LNorm in SYS2 and SYS3 are less than 1, i.e., there is no significant difference between them. Thus LNorm is also be selected as the fitted model that evaluates the software reliability measures in SYS2 and SYS3, instead of LEvMax. Similarly, the AICs of TNorm in SS3 and SS4 are close to the smallest AICs provided by TLogis and TEvMin, respectively, and the differences are less than 1. Hence TNorm is also fitted to SS3 and SS4 with respect to AIC. Table 5 shows the number of times that the model is selected as the fitted model in all the 16 types of DACS data with both 1

http://www.dacs.dtic.mil/databases/sled/  swrel.shtml

Model

SYS1

SYS2

SYS3

SS2

SS3

SS4

Exp Gamma Pareto TNorm LNorm TLogis LLogis TEvMax LEvMax TEvMin LEvMin

1960.6 1949.0 1953.1 1970.1 1951.4 1967.0 1949.3 1966.5 1952.6 1972.5 1949.0

907.5 907.0 906.1 913.4 905.6 912.1 906.4 911.9 905.5 914.5 906.7

615.5 609.1 605.7 622.9 605.6 619.9 606.5 619.6 605.3 625.0 607.8

5238.5 5242.0 5245.0 5229.3 5251.8 5229.0 5241.9 5232.9 5264.1 5226.0 5242.1

7320.3 7324.2 7326.2 7323.9 7330.4 7323.8 7325.8 7324.2 7335.9 7324.0 7326.1

5276.3 5279.8 5281.7 5277.1 5282.8 5277.5 5280.1 5277.9 5285.8 5276.2 5283.2

Table 5 The number of times selected as the fitted model. Model

Exp Gamma Pareto TNorm LNorm TLogis LLogis TEvMax LEvMax TEvMin LEvMin

Failure time data

Grouped data

AIC

BIC

AIC

BIC

24 32 7 11 15 13 21 14 15 15 31

41 18 6 8 14 9 12 10 14 10 18

13 32 0 17 7 24 21 18 6 25 24

28 27 0 13 7 22 18 16 6 18 19

failure time and grouped data. As seen in Tables 3 and 4, if the difference from the smallest AIC or BIC to the AIC or BIC of a model is less than 1, the model is also selected as a fitted model. From this result, it is found that Exp, Gamma, TLogis, LLogis and LEvMin are frequently selected as the best model in terms of both AIC and BIC. Except for Pareto, the other models, TNorm, LNorm TEvMax, LEvMax and TEvMin, are almost evenly selected. From the table, it can also be concluded that LNorm has the higher fitting ability than TNorm in the case of failure time data. On the other hand, TNorm is superior to LNorm in the case of grouped data. Empirically, the failure time data tend to draw not S-shaped curves but concave curves. TNorm is better at drawing S-shaped curves than concave ones. Therefore, LNorm and TNorm are better fitted to failure time data and grouped data, respectively. Also, compared with Pareto, both LNorm and TNorm were selected as

H. Okamura et al. / Reliability Engineering and System Safety 116 (2013) 135–141

the best model in the several cases. This implies that LNorm and TNorm can be members of the suitable set of candidates of SRGMs.

[5] [6] [7] [8] [9]

5. Conclusions This paper has developed the NHPP-based SRGMs with normally distributed failure times in the modeling framework of Langberg and Singpurwalla [15]. In addition, we have proposed effective parameter estimation algorithms based on EM principle. The proposed algorithms require only the calculation of c.d.f. of standard normal distribution in the cases of failure time data and grouped data. This is the great advantage of using LNorm and TNorm, and the proposed algorithm can easily be implemented in the application which does not have high computation library such as a spreadsheet application. Moreover, in the numerical experiment, we have investigated the fitting ability of proposed NHPP-based models. As a result, although the fitting abilities of TNorm and LNorm are inferior to several SRGMs such as Gamma, LLogis, TLogis and LEvMin, they were fitted to several data types observed in the real software projects. Hence, TNorm and LNorm should be needed to make the standard SRGMs which provide a suitable set of candidates of SRGMs. In future, we will develop a software reliability tool that can be executed on a spreadsheet application by using the proposed EM algorithms for SRGMs, and will enhance the algorithm for computing the information matrix to evaluate the estimation errors in the context of EM algorithms.

Acknowledgments This research was supported by Nanzan University Pache Research Subsidy I-A-2 for the 2011 academic year, the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (C), Grant No. 21510167 (2009–2011), Grant No. 23500047 (2011–2013) and Grant No. 23510171 (2011–2013). This paper is also an extended version of the paper [29] presented at International Conference on Quality, Reliability, Risk, Maintenance, and Safety Engineering (ICQR2MSE2011), Xi’an, China, June 17– 19, 2011. References [1] Jelinski Z, Moranda PB. Software reliability research. In: Freiberger W, editor. Statistical computer performance evaluation. New York: Academic Press; 1972. p. 465–484. [2] Goel AL, Okumoto K. Time-dependent error-detection rate model for software reliability and other performance measures. IEEE Transactions on Reliability 1979;R-28:206–211. [3] Lyu MR, editor. Handbook of software reliability engineering. New York: McGraw-Hill; 1996. [4] Musa JD, Iannino A, Okumoto K. Software reliability, measurement, prediction, application. New York: McGraw-Hill; 1987.

[10]

[11] [12] [13]

[14]

[15] [16] [17]

[18]

[19]

[20]

[21]

[22] [23]

[24]

[25]

[26]

[27]

[28] [29]

141

Musa JD. Software reliability engineering. New York: McGraw-Hill; 1999. Pham H. Software reliability. Singapore: Springer; 2000. Xie M. Software reliability modelling. Singapore: World Scientific; 1991. Goel AL. Software reliability models: assumptions, limitations and applicability. IEEE Transactions on Software Engineering 1985;SE-11:1411–1423. Musa JD, Okumoto K. A logarithmic Poisson execution time model for software reliability measurement. In: Proceedings of the 7th international conference on software engineering (ICSE-1084). IEEE CS Press/ACM, 1984. p. 230–8. Ohba M. Inflection S-shaped software reliability growth model. In: Osaki S, Hatoyama Y, editors. Stochastic models in reliability theory. Berlin: SpringerVerlag; 1984. p. 144–165. Ohba M. Software reliability analysis. IBM Journal on Research & Development 1984;28:428–443. Yamada S, Ohba M, Osaki S. S-shaped reliability growth modeling for software error detection. IEEE Transactions on Reliability 1983;R-32:475–478. Zhao M, Xie M. On maximum likelihood estimation for a general nonhomogeneous Poisson process. Scandinavian Journal of Statistics 1996;23: 597–607. Pham H, Zhang X. An NHPP software reliability models and its comparison. International Journal of Reliability, Quality and Safety Engineering 1997;4: 269–282. Langberg N, Singpurwalla ND. Unification of some software reliability models. SIAM Journal on Scientific Computing 1985;6:781–790. Littlewood B. Rationale for a modified Duane model. IEEE Transactions on Reliability 1984;R-33:157–159. Abdel-Ghaly AA, Chan PY, Littlewood B. Evaluation of competing software reliability predictions. IEEE Transactions on Software Engineering 1986;SE12:950–967. Gokhale SS, Trivedi KS. Log-logistic software reliability growth model. In: Proceedings of the 3rd IEEE international high-assurance systems engineering symposium (HASE-1998). IEEE CS Press, 1998. p. 34–41. Ohishi K, Okamura H, Dohi T. Gompertz software reliability model: estimation algorithm and empirical validation. Journal of Systems and Software 2009;82:535–543. Achcar JA, Dey DK, Niverthi M. A Bayesian approach using nonhomogeneous Poisson processes for software reliability models. In: Basu AP, Basu KS, Mukhopadhyay S, editors. Frontiers in reliability. Singapore: World Scientific; 1998. p. 1–18. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. Journal of Royal Statistical Society B 1977;B39:1–38. Wu CFJ. On the convergence properties of the EM algorithm. Annals of Statistics 1983;11:95–103. Okamura H, Dohi T, Osaki S. EM algorithms for logistic software reliability models. In: Proceedings of the 22nd IASTED international conference on software engineering. ACTA Press; 2004. p. 263–8. Okamura H, Watanabe Y, Dohi T. An iterative scheme for maximum likelihood estimation in software reliability modeling, In: Proceedings of the 14th international symposium software reliability engineering, 2003. p. 246–56. Okamura H, Murayama A, Dohi T. EM algorithm for discrete software reliability models: a unified parameter estimation method. In: Proceedings of the 8th IEEE international symposium on high assurance systems engineering, 2004. p. 219–28. Yamada S. A stochastic software reliability growth model with Gompertz curve. Transactions of Information Processing Society of Japan (in Japanese) 1992;33:964–969. Akaike H. Information theory and an extension of the maximum likelihood principle. In: Petrov BN, Csaki F, editors. Proceedings of the 2nd international symposium on information theory. 1973. p. 267–81. Schwarz G. Estimating the dimension of a model. Annals of Statistics 1978;6:461–464. Okamura H, Dohi T, Osaki S. Software reliability growth model with normal distribution and its parameter estimation. In: Proceedings of the international conference on quality, reliability, maintenance and safety engineering (ICQR2MSE 2011). 2011. p. 424–9.