- Email: [email protected]
Poisson Process Models in Software Reliability
Copyrighl © IFAC Induslrial Process Control Syslems. Bruges. Belgium. 1988
POISSON PROCESS MODELS IN SOFTWARE RELIABILITY I. D. AI-Ayoubi* and M.
J.
Newby**
*Postgraduate School of Mathematical Scientes , Bradford, UK **Faculty of Industrial Engineering, Eindhoven University of Techn ology, Eindhoven, The Netherlands.
Abstract . A simple generalisation of the Goel-Okumoto models for software reliability is developed . A classical likelihood approach demonstrates that parameter estimation may be unsatisfactory and unstable . Further, although methods based on goodness of fit are more satisfactory they still produce widely differing predictions . It is suggested that Bayesian ljIethods offer a solution to the problems . The models are equally applicable in reliability growth/degradation situations . Keywords . Computer reliability theory .
software ;
data
reduction
and
analysis ;
random
processes ;
INTRODUCTION Software reliability modelling and reliability growth modelling make extensive use of Poisson and non-homogeneous Poisson processes (Ascher and Feingold, 1984). Poisson models are characterised by their mean value functions, or by their rate of occurrence of failures (Ascher and Feingold, 1984) . Poisson models are not strictly correct for software because the assumption of a fixed initial number of unrevealed failures suggests an order statistic process . However, in practice the difference appears to be unimportant .
This process in the time defined by T - G(u;k) is a homogeneous Poisson process with rate er defined on the interval [0,1] and with an absorbing barrier at 1 .
Goel and Okumoto (1979, 1985) defined by mean value functions
inter-failure times {Xi}t_, with ti -
exhibit
Data
may
failure
be
recorded
times {ti}t_"
as
a
sequence
of
known
as
a
sequence
of
known
L Xj'
or as
i
models
j-,
grouped data {(fi,ti)}!_, where fi events occur in
a,b > 0 ,
Mo(t) - a{l - exp(-bt)}
ESTIMATION
t;> 0
(ti-, ,ti)
and
the
total
number
observed
is
s
and M,(t) - a{l - exp(-bt k )}
a , b,k > 0 ,
L
~(w)
- m(ti_, + w) exp (-[M(ti_, + w) - M(ti_,)]}'
fi ' The density function for inter-failure i-, interval , w, given a failure at ti-, is (Cox and Lewis, 1966)
t;> O.
The main features of the models are that the functions Mi are bounded and that Mi(~) - a is the initial number of unrevealed faults . This suggests models of the form M(t)
n -
The likelihood functions are then easily found . For data known exactly
aC(u ;k)
where u - t/P with er,f3,k > 0 and t ;> 0 and with G(O) - 0, G(~) - 1 . In other words , any statistical G: [O,~) ... [0,1), will do the job . the rate function can be written m(t)
where t* is the last event time, tn, if the observation stopped then, or t* is the time at which observation ceased.
distribution, In this case
For grouped data
~ g(u;k)
where g is the density of G. Thus er represents the number of faults, f3 a change of scale , k a shape parameter and G describes the passage of time .
With our specialisation
When the above a •• umptions are made the failures form a sequence of events in a non-homogeneous Poi.son proce •• with probability function Pn(t) -
and
[M~~)]neXP[_M(t)] _ [aC~I)]neXp[_aC(U)].
L,
123
s { }fi . ern e -aC(u) S K n G(ui) - G(ui_,) i-,
124
1. D. Al-:houbi and \1.
In both cases the max imum likelihood estimator of o is n
where u* - t*/ ~ recorded .
and
t*
is
the
latest
J.
:\e"by
which is discussed later. Convergence was good in that the gradient norms were of the order 10- 6 at the final value, and the Hessian matrix was sign definite .
time
Three simple forms were chosen for G. The second, the Weibull, includes the Goel-Okumoto models. The distributions are : 1
Pare to
G(u) - 1
Weibull
k G(u) - 1 - exp[-u ] ,
Extreme value
( l+u)k'
k G(u) - exp[-u- ],
k > 0, u
>
k > 0, u
> 0;
TABLE 1
Pareto
Weibull
Ext. value
- 143 . 7
- 143.9
-143 . 8
0
377 . 5
160 . 1
502 . 7
~
1.9
10 . 7
72.5
k
0 .2
0.7
0.2
0
346 .6
21. 0
197 . 6
s .d. ~
.2
3.5
100.7
k
0.2
0.1
0.1
corr (o ,~ )
-0.8
0.7
0 . 97
(o,k)
-0.9
-0 . 6
-0 . 9
corr ( ~,k)
0.9
-0.7
-0.96
coefficient of variation of 0
0.9
0.1
0 .4
0 . 1889
0.2004
0.1999
Expected faults t- 20
127
127
127
Expected faults remaining, t-20
250
33
376
0 . 15
0 . 15
0.11
0;
log-like 1 ihood
k > 0, u
> O.
To examine the behaviour of the likelihood function a profile likelihood which is a function of only ~ and k is constructed by replacing 0 by O. The likelihood surface and likelihood contours for the profile likelihood can now be plotted . The profile likelihoods are
s.d .
s . d.
corr
Kn
n n
{g(u .) } ~
i-, ~G(u*)
and
* , k) L,W
n -n n e
Summary of Results
Kolmogorov-Smirnov stat .
G(u*)n
Problems in estimation arise from numerical instabilities in calculating the likelihoods, and from the shape of the likelihood function itself. The form of L~ shows that the product term will be rather flat since G increases and so each [G(ui) - G(ui_,)] / G(u*) is between zero and one . The estimation was carried out by using a modified Newton-Raplison technique implemented as the NAG routine E04KBF ( NAG , 1988 ) to find the maximum of the negative of the log-likelihood . The results of the estimation routine were unstable in that the estimates o , ~ , k of o,~,k were strongly dependent on the initial values and the direction of approach to the extreme. Clearly the loglikelihood had local maxima and possibly saddle points which diverted the maximisation routine from the global maximum. Moreover, for some extremes the Hessian matrix was not positive definite . RESULTS The approach is illustrated by applying it to Musa's data set one (Musa , 1980) and the results are typical of the performance of the models . The instabili ty has already been mentioned on the other hand goodness of fit appeared to be largely unaffected by the parameter values. Lastly, the different models produced large differences in the value of 0, the estimate of the total number of errors . The results are given in Table 1 . The results were obtained by choosing a parameter region after examining the graphical evidence
Estimated reliability, t-20
The choice of well behaved regions of the profile like lihood produces superfic ially satisfac tory results . Certainly all three models are accepted by the Kolmogorov-Smirnov test at the 5% level and the Pare to appears to be a marginally better fit than the other two . Further evidence of good fit is gained from the v alue of the mean value functions at t - 20. Closer examination" shows large differences in the values of 0 and consequently in estimates of the time and cost o~ fault removal. The coefficient of variation of 0 is very large in the Pare to model, moderate in the extreme value and small for the Weibull showing that in two of the models there is high variability in the estimates . Lastly the correlation coefficients are high indicating a need for a better parametric form (Smith and others , 1987) . The graphical evidence in Fig. 1 to Fig. 6 shows that the profile likelihoods have several local maxima and saddle points and the correlation between ~ and k can be seen from the existence of ridges. Neverthe less all three forms have distinct global maxima so that the graphs play a valuable part in guiding the estimation process .
BAYESIAN APPRO.
Poisson Process Models in Software Reliability
Since the likelihoods are known and may be written in the form for i - 0, 1 the posterior distribution of may be found from Bayes result g(Cl',~,k/data)
Cl',~,k
with prior
Li(data;a,~,k)
where the integral is over an appropriate space . Furthermore , in both cases n -aG(u*) Cl' e ~(data/~,k) where u* is defined in a way appropriate to the data. This form for Li suggests the use of gamma conjugate priors for Cl' chosen so that Cl' is independent of ~ and k. To obtain results an effective numerical technique (Smith and others, 1987) is needed in conjunction with the graphical techniques described above . The assumption of an indep~ndent gamma prior for Cl' still leaves considerable problems with the integration .
REFERENCES Ascher, H. , H. Feingold (1984) . Repairable Systems Reliability. Marcel Dekker, New York . Goel , A.L., K. Okumoto (1979). A time-dependent error detection rate model for software reliability and other performance measures . IEEE Trans . on Reliability, R-28, no . 3, 206-211. Goel, A.L . (1985) . Software reliability models : assumptions limitations and applicability . IEEE Trans. on Software Eng., SE-ll, no. 12. Cox, D. R ., P.A. Lewis (1966) . The Statistical Analysis of Series of Events. Methuen, London . Numerical Algorithms Group (1988) . NAG Manual. Musa, D., (1980) . Software Reliability Data . Data and Analysis Center for Software, Rome, New York. Smith, A.J.M., A. M. Skene, J.E.M . Shaw, and J .C. Naylor (1987). Progress with numerical and graphical methods for practical Bayesian statistics . The Statistician, 36, 75-82 .
125
I. D. AI-Ayoubi and M. J. Newby
126
Fig. 1.
Profile likelihood for Pareto model .
Fig. 3.
Profile likelihood for Weibull model.
3.Q
2.5
2. 0
k
1 S
l . fJ
0. '
7. J
i Q. Q'
12
:l
i ~d
I i I J
20 Q
B
Fig. 2.
Contours of profi le li ke lihood for P areto mode 1.
Fig. 4.
Contours of profile likelihood for Weibull model.
Poisson Process Models in Softwar e Reliability
Profile likeliho o d for ex treme value model.
Fig. S.
LO 0. 9 0. 8 O. I
O. b k
~. )'
~.'
0.3'
O. "
~
O. J 10
Fig. 6.
30
50' B
f pro f 1· le like lihood for Contour s 0 1 e mode l. extreme va u
127
POISSON PROCESS MODELS IN SOFTWARE RELIABILITY I. D. AI-Ayoubi* and M.
J.
Newby**
*Postgraduate School of Mathematical Scientes , Bradford, UK **Faculty of Industrial Engineering, Eindhoven University of Techn ology, Eindhoven, The Netherlands.
Abstract . A simple generalisation of the Goel-Okumoto models for software reliability is developed . A classical likelihood approach demonstrates that parameter estimation may be unsatisfactory and unstable . Further, although methods based on goodness of fit are more satisfactory they still produce widely differing predictions . It is suggested that Bayesian ljIethods offer a solution to the problems . The models are equally applicable in reliability growth/degradation situations . Keywords . Computer reliability theory .
software ;
data
reduction
and
analysis ;
random
processes ;
INTRODUCTION Software reliability modelling and reliability growth modelling make extensive use of Poisson and non-homogeneous Poisson processes (Ascher and Feingold, 1984). Poisson models are characterised by their mean value functions, or by their rate of occurrence of failures (Ascher and Feingold, 1984) . Poisson models are not strictly correct for software because the assumption of a fixed initial number of unrevealed failures suggests an order statistic process . However, in practice the difference appears to be unimportant .
This process in the time defined by T - G(u;k) is a homogeneous Poisson process with rate er defined on the interval [0,1] and with an absorbing barrier at 1 .
Goel and Okumoto (1979, 1985) defined by mean value functions
inter-failure times {Xi}t_, with ti -
exhibit
Data
may
failure
be
recorded
times {ti}t_"
as
a
sequence
of
known
as
a
sequence
of
known
L Xj'
or as
i
models
j-,
grouped data {(fi,ti)}!_, where fi events occur in
a,b > 0 ,
Mo(t) - a{l - exp(-bt)}
ESTIMATION
t;> 0
(ti-, ,ti)
and
the
total
number
observed
is
s
and M,(t) - a{l - exp(-bt k )}
a , b,k > 0 ,
L
~(w)
- m(ti_, + w) exp (-[M(ti_, + w) - M(ti_,)]}'
fi ' The density function for inter-failure i-, interval , w, given a failure at ti-, is (Cox and Lewis, 1966)
t;> O.
The main features of the models are that the functions Mi are bounded and that Mi(~) - a is the initial number of unrevealed faults . This suggests models of the form M(t)
n -
The likelihood functions are then easily found . For data known exactly
aC(u ;k)
where u - t/P with er,f3,k > 0 and t ;> 0 and with G(O) - 0, G(~) - 1 . In other words , any statistical G: [O,~) ... [0,1), will do the job . the rate function can be written m(t)
where t* is the last event time, tn, if the observation stopped then, or t* is the time at which observation ceased.
distribution, In this case
For grouped data
~ g(u;k)
where g is the density of G. Thus er represents the number of faults, f3 a change of scale , k a shape parameter and G describes the passage of time .
With our specialisation
When the above a •• umptions are made the failures form a sequence of events in a non-homogeneous Poi.son proce •• with probability function Pn(t) -
and
[M~~)]neXP[_M(t)] _ [aC~I)]neXp[_aC(U)].
L,
123
s { }fi . ern e -aC(u) S K n G(ui) - G(ui_,) i-,
124
1. D. Al-:houbi and \1.
In both cases the max imum likelihood estimator of o is n
where u* - t*/ ~ recorded .
and
t*
is
the
latest
J.
:\e"by
which is discussed later. Convergence was good in that the gradient norms were of the order 10- 6 at the final value, and the Hessian matrix was sign definite .
time
Three simple forms were chosen for G. The second, the Weibull, includes the Goel-Okumoto models. The distributions are : 1
Pare to
G(u) - 1
Weibull
k G(u) - 1 - exp[-u ] ,
Extreme value
( l+u)k'
k G(u) - exp[-u- ],
k > 0, u
>
k > 0, u
> 0;
TABLE 1
Pareto
Weibull
Ext. value
- 143 . 7
- 143.9
-143 . 8
0
377 . 5
160 . 1
502 . 7
~
1.9
10 . 7
72.5
k
0 .2
0.7
0.2
0
346 .6
21. 0
197 . 6
s .d. ~
.2
3.5
100.7
k
0.2
0.1
0.1
corr (o ,~ )
-0.8
0.7
0 . 97
(o,k)
-0.9
-0 . 6
-0 . 9
corr ( ~,k)
0.9
-0.7
-0.96
coefficient of variation of 0
0.9
0.1
0 .4
0 . 1889
0.2004
0.1999
Expected faults t- 20
127
127
127
Expected faults remaining, t-20
250
33
376
0 . 15
0 . 15
0.11
0;
log-like 1 ihood
k > 0, u
> O.
To examine the behaviour of the likelihood function a profile likelihood which is a function of only ~ and k is constructed by replacing 0 by O. The likelihood surface and likelihood contours for the profile likelihood can now be plotted . The profile likelihoods are
s.d .
s . d.
corr
Kn
n n
{g(u .) } ~
i-, ~G(u*)
and
* , k) L,W
n -n n e
Summary of Results
Kolmogorov-Smirnov stat .
G(u*)n
Problems in estimation arise from numerical instabilities in calculating the likelihoods, and from the shape of the likelihood function itself. The form of L~ shows that the product term will be rather flat since G increases and so each [G(ui) - G(ui_,)] / G(u*) is between zero and one . The estimation was carried out by using a modified Newton-Raplison technique implemented as the NAG routine E04KBF ( NAG , 1988 ) to find the maximum of the negative of the log-likelihood . The results of the estimation routine were unstable in that the estimates o , ~ , k of o,~,k were strongly dependent on the initial values and the direction of approach to the extreme. Clearly the loglikelihood had local maxima and possibly saddle points which diverted the maximisation routine from the global maximum. Moreover, for some extremes the Hessian matrix was not positive definite . RESULTS The approach is illustrated by applying it to Musa's data set one (Musa , 1980) and the results are typical of the performance of the models . The instabili ty has already been mentioned on the other hand goodness of fit appeared to be largely unaffected by the parameter values. Lastly, the different models produced large differences in the value of 0, the estimate of the total number of errors . The results are given in Table 1 . The results were obtained by choosing a parameter region after examining the graphical evidence
Estimated reliability, t-20
The choice of well behaved regions of the profile like lihood produces superfic ially satisfac tory results . Certainly all three models are accepted by the Kolmogorov-Smirnov test at the 5% level and the Pare to appears to be a marginally better fit than the other two . Further evidence of good fit is gained from the v alue of the mean value functions at t - 20. Closer examination" shows large differences in the values of 0 and consequently in estimates of the time and cost o~ fault removal. The coefficient of variation of 0 is very large in the Pare to model, moderate in the extreme value and small for the Weibull showing that in two of the models there is high variability in the estimates . Lastly the correlation coefficients are high indicating a need for a better parametric form (Smith and others , 1987) . The graphical evidence in Fig. 1 to Fig. 6 shows that the profile likelihoods have several local maxima and saddle points and the correlation between ~ and k can be seen from the existence of ridges. Neverthe less all three forms have distinct global maxima so that the graphs play a valuable part in guiding the estimation process .
BAYESIAN APPRO.
Poisson Process Models in Software Reliability
Since the likelihoods are known and may be written in the form for i - 0, 1 the posterior distribution of may be found from Bayes result g(Cl',~,k/data)
Cl',~,k
with prior
Li(data;a,~,k)
where the integral is over an appropriate space . Furthermore , in both cases n -aG(u*) Cl' e ~(data/~,k) where u* is defined in a way appropriate to the data. This form for Li suggests the use of gamma conjugate priors for Cl' chosen so that Cl' is independent of ~ and k. To obtain results an effective numerical technique (Smith and others, 1987) is needed in conjunction with the graphical techniques described above . The assumption of an indep~ndent gamma prior for Cl' still leaves considerable problems with the integration .
REFERENCES Ascher, H. , H. Feingold (1984) . Repairable Systems Reliability. Marcel Dekker, New York . Goel , A.L., K. Okumoto (1979). A time-dependent error detection rate model for software reliability and other performance measures . IEEE Trans . on Reliability, R-28, no . 3, 206-211. Goel, A.L . (1985) . Software reliability models : assumptions limitations and applicability . IEEE Trans. on Software Eng., SE-ll, no. 12. Cox, D. R ., P.A. Lewis (1966) . The Statistical Analysis of Series of Events. Methuen, London . Numerical Algorithms Group (1988) . NAG Manual. Musa, D., (1980) . Software Reliability Data . Data and Analysis Center for Software, Rome, New York. Smith, A.J.M., A. M. Skene, J.E.M . Shaw, and J .C. Naylor (1987). Progress with numerical and graphical methods for practical Bayesian statistics . The Statistician, 36, 75-82 .
125
I. D. AI-Ayoubi and M. J. Newby
126
Fig. 1.
Profile likelihood for Pareto model .
Fig. 3.
Profile likelihood for Weibull model.
3.Q
2.5
2. 0
k
1 S
l . fJ
0. '
7. J
i Q. Q'
12
:l
i ~d
I i I J
20 Q
B
Fig. 2.
Contours of profi le li ke lihood for P areto mode 1.
Fig. 4.
Contours of profile likelihood for Weibull model.
Poisson Process Models in Softwar e Reliability
Profile likeliho o d for ex treme value model.
Fig. S.
LO 0. 9 0. 8 O. I
O. b k
~. )'
~.'
0.3'
O. "
~
O. J 10
Fig. 6.
30
50' B
f pro f 1· le like lihood for Contour s 0 1 e mode l. extreme va u
127