# A Comparative Study of Data Transformations for Wavelet Shrinkage Estimation with Application to Software Reliability Assessment

Authors: Xiao Xiao; Tadashi Dohi
Journal: Advances in Software Engineering
 (2012)

---

## Abstract

In our previous work, we proposed wavelet shrinkage estimation (WSE) for nonhomogeneous Poisson process (NHPP)-based software reliability models (SRMs), where WSE is a data-transform-based nonparametric estimation method. Among many variance-stabilizing data transformations, the Anscombe transform and the Fisz transform were employed. We have shown that it could provide higher goodness-of-fit performance than the conventional maximum likelihood estimation (MLE) and the least squares estimation (LSE) in many cases, in spite of its non-parametric nature, through numerical experiments with real software-fault count data. With the aim of improving the estimation accuracy of WSE, in this paper we introduce other three data transformations to preprocess the software-fault count data and investigate the influence of different data transformations to the estimation accuracy of WSE through goodness-of-fit test.


---

## Body

## 1. Introduction

In the field of software reliability engineering, the quantitative assessment of software reliability has become one of the main issues of this area. Especially, people are interested in finding severalsoftware intensity functions from the software-fault count data observed in the software testing phases, since the software intensity function in discrete time denotes the number of software faults detected per unit time. This directly makes it possible to estimate the number of remaining software faults and the quantitative software reliability, which is defined as the probability that software system does not fail during a specified time period under a specified operational environment. Moreover, these evaluation measures can be used in the decision making such as allocation of development resources and software release scheduling.Therefore, we are interested in developing a high-accuracy estimation method for the software intensity function.Among over hundreds of software reliability models (SRMs) [1–3], nonhomogeneous Poisson process (NHPP)-based SRMs have gained much popularity in actual software testing phases. In many cases, the NHPP-based SRM is formulated as a parametric model, where the mean value function or its difference in discrete time or derivative in continuous time, called “software intensity function,” can be considered as a unique parameter to govern the probabilistic property. One class of parametric NHPP-based SRMs is concerned with modeling the number of software faults detected in testing phases, initiated by Goel and Okumoto [4]. Afterwards, many parametric NHPP-based SRMs were proposed in the literatures [5–8] from various points of view. However, it is well known in the software reliability engineering community that there does not exist a uniquely best parametric NHPP-based SRM which can fit every type of software-fault count data. This fact implies that nonparametric methods without assuming parametric form should be used to describe the software debugging phenomenon which is different in each testing phase.Apart from the traditional Bayesian framework, some frequentist approaches based on non-parametric statistics were introduced to estimate the quantitative software reliability. Sofer and Miller [9] used an elementary piecewise linear estimator of the NHPP intensity function from the software-fault detection time data and proposed a smoothing technique by means of quadratic programming. Gandy and Jensen [10] applied the kernel-based estimator to estimate the NHPP intensity function in a non-parametric way. Barghout et al. [11] also proposed a kernel-based non-parametric estimation for the order statistics-based SRMs, where the likelihood cross-validation and the prequential likelihood approaches were used to estimate the bandwidth. Wang et al. [12] applied the similar kernel method to the NHPP-based SRMs, where they focused on the local likelihood method with a locally weighted log-likelihood function. By combining the non-parametric estimation with the Bayesian framework, El-Aroui and Soler [13] and Wilson and Samaniego [14] developed non-parametric Bayesian estimation methods. It should be noted that, generally, these non-parametric estimation methods require high computational cost, which, in some cases, may be almost similar to or greater than an effort on model selection in the parametric SRMs.Another class of non-parametric estimation methods for NHPP-based SRMs is the wavelet analysis-based approach, initiated by Xiao and Dohi [15]. They proposed the wavelet shrinkage estimation (WSE), which does not require solving any optimization problem, so that the implementation of estimation algorithms is rather easy than the other non-parametric methods. They compared their method with the conventional maximum likelihood estimation (MLE) and the least squares estimation (LSE) through goodness-of-fit test. It has been shown that WSE could provide higher goodness-of-fit performance than MLE and LSE in many cases, in spite of its non-parametric nature, through numerical experiments with real software-fault count data.The fundamental idea of WSE is to remove the noise included in the observed software-fault count data to get a noise-free estimate of the software intensity function. It is performed through the following three-step procedure. First, the noise variance is stabilized by applying the data transformation to the data. This produces a time-series data in which the noise can be treated as Gaussian white noise. Second, the noise is removed using “Haar-wavelet” based denoising algorithm. Third, an inverse data transformation is applied to the denoised time-series data, obtaining the estimate of the software intensity function. Among many variance-stabilizing data transformations, the Anscombe transform [16] and the Fisz transform [17] were employed in the previous work [15]. The other well-known square root data transformations are Bartlett transform [18] and Freeman transform [19]. Both Anscombe transform and Freeman transform are actually natural extensions of the Bartlett transform.This paper focuses on the first step of WSE and aims at identifying and emphasizing the influence that the data transformation exerts on the accuracy of WSE. The remaining part of this paper is planned as follows. In Section2, we give a preliminary on NHPP-based software reliability modeling. Section 3 is devoted to introduce data transformations. Section 4 describes the WSE for NHPP-based SRMs in details. In Section 5, we carry out the real project data analysis and illustrate numerical examples to examine the effectiveness of the proposed methods. Finally the paper is concluded with future researches in Section 6.


## 2. NHPP-Based Software Reliability Modeling

Suppose that the number of software faults detected through a system test is observed at discrete timei=0,1,2,…. Let Yi and Ni=∑k=0iYk denote the number of software faults detected at testing date i and its cumulative value, where Y0=N0=0 is assumed without any loss of generality. The stochastic process {Ni:i=0,1,2,…} is said to be a discrete non-homogeneous Poisson process (D-NHPP) if the probability mass function at time i is given by(1)Pr{Ni=m}={Λi}mm!exp⁡{-Λi},m=0,1,2,…,
where Λi=E[Ni] is called the mean value function of a D-NHPP and means the expected cumulative number of software faults detected by testing date i. The function λi=Λi-Λi-1(i≥1) is called the discrete intensity function and implies the expected number of faults detected at testing date i, say λi=E[Yi].MLE is one of the most commonly used parametric estimation method. Letθ denote the vector of parameters in the mean value function Λi=Λi(θ), and xi(yi) denote the realization of Ni(Yi). When n software faults are detected, the log-likelihood function of a D-NHPP is given by(2)LLF(θ)=∑i=1n(xi-xi-1)ln⁡[Λi(θ)-Λi-1(θ)]-Λn(θ)-∑i=1nln⁡[(xi-xi-1)!],
where xi=∑k=0iyk and y0=x0=0. Then, the maximum likelihood estimate of θ, say θ̂, is given by the solution of argmaxθℒℒℱ(θ). Therefore, the estimate of the software intensity function λi=λi(θ)(i=0,1,2,…) can be obtained by λ̂i=Λi(θ̂)-Λi-1(θ̂)(i≥1) with λ0=0. In the following sections, we consider the problem of estimating the software intensity function from the noise-involved observation yi, in a non-parametric way.


## 3. Variance Stabilizing Data Transformation

It is very familiar to make use of data transformations (DTs) to stabilize the variance of Poisson data. By using DT, the software-fault count data which follow the D-NHPP are approximately transformed to the Gaussian data. The most fundamental data-transform tool in statistics is theBartlett transform (BT) [18]. Let ηi denote the Poisson white noise, that is,(3)Yi=λi+ηi,i=1,2…,n.
Taking the BT, the random variables:(4)Bi=2Yi,i=1,2,…,n
can be approximately regarded as Gaussian random variables with the normal distribution N(2λi-1/4,1), so that the realizations:(5)bi=2yi,i=1,2,…,n
can be considered as samples from N(2λi-1/4,1). That is, the transformed realizations bi(i=1,2,…,n) by the BT are the ones from the normally distributed random variables:(6)Bi=λi′+νi,i=1,2,…,n,
where λi′=2λi-1/4 is the transformed software intensity function, and νi is the Gaussian white noise with unit variance.Bartlett [18] also showed that(7)bi=2yi+12,i=1,2,…,n
is a better transformation since it provides a constant variance more closely to 1, even when the mean of Yi is not large.The Anscombe transform (AT) [16] is a natural extension of BT and is employed in our previous work [15]. AT is of the following form:(8)si=2yi+38,i=1,2,…,n,
where si can be considered as observations of Gaussian random variable Si=2Yi+3/8 with the normal distribution N(2λi+1/8,1). Freeman and Tukey [19] proposed the following square-root transform (we call it FT), which is also an extension of BT:(9)fi=yi+1+yi,i=1,2,…,n.
They showed that the variance of Gaussian random variable Fi=2Yi+1+Yi is the nearest to 1 among BT, AT, and FT if the mean of Yi is small. Recently, these variance stabilization techniques were used to LSE of the mean value function for the NHPP-based SRMs [20]. Table 1 summaries the DTs mentioned above.Table 1
Representative data transformaions.
Data transformationInverse data transformationDistribution of random variable after DTBT1 [18]Bi=2YiYi=1/4×{(Bi)2}N(2λi-1/4,1)BT2 [18]Bi=2Yi+1/2Yi=1/4×{(Bi)2-2}N(2λi+1/4,1)AT [16]Si=2Yi+3/8Yi=1/4×{(Si)2-3/2}N(2λi+1/8,1)FT [19]Fi=Yi+1+YiYi=1/4×{(Fi)2+(Fi)-2-2}N(2λi,1)As mentioned in Section1, the first step of WSE is to apply the normalizing and variance-stabilizing DTs to the observed software-fault count data. In this paper, we employ BT, AT, and FT in the first and the third steps of WSE. Then, the target of denoising in the second step of WSE is the transformed data bi,si or fi(i=1,2,…,n). Letting bi′,si′, and fi′ denote the denoised bi,si, and fi, respectively, the estimate of the original software intensity function λi can be obtained by taking the inverse DT of bi′,si′ and fi′, as given in Table 1.


## 4. Wavelet Shrinkage Estimation for NHPP-Based SRM

### 4.1. Haar-Wavelet-Based Denoising Procedure

The Haar-wavelet-based shrinkage technique can be used as a denoising algorithm for the second step of WSE. In general, the noise removal is performed through the following three steps: (i) expanding the transformed time-series data to obtain the empirical wavelet coefficients, (ii) removing the noise included in the empirical wavelet coefficients using thresholding method, and (iii) making use of the denoised coefficients to calculate the estimate of the transformed software intensity function.


#### 4.1.1. Haar Wavelet Transform

The Haar scaling function and the Haar wavelet function are defined as(10)ϕ(i)={1(0≤i<1)0(otherwise),ψ(i)={1(0≤i<12)-1(12≤i<1)0(otherwise),
respectively. By introducing the scaling parameterj and shifting parameterk, the Haar father wavelet and the Haar mother wavelet are defined by(11)ϕj,k(i)=2-j/2ϕ(2-ji-k),ψj,k(i)=2-j/2ψ(2-ji-k),
respectively. Then the target function, transformed software intensity function λi′(i=1,2,…,n), can be expressed in the following equation:(12)λi′=∑k=02j0-1αj0,kϕj0,k(i)+∑j=j0∞∑k=02j-1βj,kψj,k(i),
where(13)αj0,k=∑i=1nλi′ϕj0,k(i),(14)βj,k=∑i=1nλi′ψj,k(i)
are called the scaling coefficients and the wavelet coefficients, respectively, for any primary resolution level j0(≥0). Due to the implementability, it is reasonable to set an upper limit instead of ∞ for the resolution level j. In other words, the highest resolution level must be finite in practice. We use J to denote the highest resolution level. That is, the range of j in the second term of (12) is j∈[j0,J]. The mapping from function λi′ to coefficients (αj0,k,βj,k) is called the Haar wavelet transform (HWT).Sincebi,si or fi(i=1,2,…,n) can be considered as the observation of λi′, the empirical scaling coefficients cj0,k and the empirical wavelet coefficients dj,k of λi′ can be calculated by (13) and (14) with λi′ replaced by bi,si or fi. The noises involved in the empirical wavelet coefficients dj,k can be removed by the thresholding method that we will introduce later. Finally, the estimate of λi′ can be obtained by taking the inverse HWT with denoised empirical coefficients.


#### 4.1.2. Thresholding

In denoising the empirical wavelet coefficients, the common choices of thresholding method are the hard thresholding:(15)δτ(u)=u1|u|>τ,
and the soft thresholding:(16)δτ(u)=sgn⁡(u)(|u|-τ)+,
for a fixed threshold level τ(>0), where 1A is the indicator function of an event A, sgn(u) is the sign function of u and (u)+=max(0,u). There are many methods to determine the threshold level τ. In this paper, we use the universal threshold [21] and the “leave-out-half” cross-validation threshold [22]:(17)τ=2log⁡n,τ=(1-log⁡2log⁡n)-1/2τ(n2),
where n is the length of the observation si. Hard thresholding is a “keep” or “kill” rule, while soft thresholding is a “shrink” or “kill” rule. Both thresholding methods and both threshold levels will be employed to work on the empirical wavelet coefficients dj,k for denoising.


### 4.2. Wavelet Shrinkage Estimation for NHPP-Based SRM

Since the software-fault count data is Poisson data, the preprocessing is necessary before making use of the Haar-wavelet-based denoising procedure. Xiao and Dohi [15] combined data transformation and the standard denoising procedure to propose the wavelet shrinkage estimation (WSE) for the D-NHPP-based SRMs. They call the HWT combined with AT, the Haar-Anscombe transform (HAT). Similarly, we call the HWT combined with BT and FT, the Haar-Bartlett transform (HBT) and the Haar-Freeman transform (HFT), respectively. In the numerical study, we investigate the goodness-of-fit performance of HBT- and HFT-based WSEs and compare them with the HAT-based WSE [15].


## 4.1. Haar-Wavelet-Based Denoising Procedure

The Haar-wavelet-based shrinkage technique can be used as a denoising algorithm for the second step of WSE. In general, the noise removal is performed through the following three steps: (i) expanding the transformed time-series data to obtain the empirical wavelet coefficients, (ii) removing the noise included in the empirical wavelet coefficients using thresholding method, and (iii) making use of the denoised coefficients to calculate the estimate of the transformed software intensity function.


### 4.1.1. Haar Wavelet Transform

The Haar scaling function and the Haar wavelet function are defined as(10)ϕ(i)={1(0≤i<1)0(otherwise),ψ(i)={1(0≤i<12)-1(12≤i<1)0(otherwise),
respectively. By introducing the scaling parameterj and shifting parameterk, the Haar father wavelet and the Haar mother wavelet are defined by(11)ϕj,k(i)=2-j/2ϕ(2-ji-k),ψj,k(i)=2-j/2ψ(2-ji-k),
respectively. Then the target function, transformed software intensity function λi′(i=1,2,…,n), can be expressed in the following equation:(12)λi′=∑k=02j0-1αj0,kϕj0,k(i)+∑j=j0∞∑k=02j-1βj,kψj,k(i),
where(13)αj0,k=∑i=1nλi′ϕj0,k(i),(14)βj,k=∑i=1nλi′ψj,k(i)
are called the scaling coefficients and the wavelet coefficients, respectively, for any primary resolution level j0(≥0). Due to the implementability, it is reasonable to set an upper limit instead of ∞ for the resolution level j. In other words, the highest resolution level must be finite in practice. We use J to denote the highest resolution level. That is, the range of j in the second term of (12) is j∈[j0,J]. The mapping from function λi′ to coefficients (αj0,k,βj,k) is called the Haar wavelet transform (HWT).Sincebi,si or fi(i=1,2,…,n) can be considered as the observation of λi′, the empirical scaling coefficients cj0,k and the empirical wavelet coefficients dj,k of λi′ can be calculated by (13) and (14) with λi′ replaced by bi,si or fi. The noises involved in the empirical wavelet coefficients dj,k can be removed by the thresholding method that we will introduce later. Finally, the estimate of λi′ can be obtained by taking the inverse HWT with denoised empirical coefficients.


### 4.1.2. Thresholding

In denoising the empirical wavelet coefficients, the common choices of thresholding method are the hard thresholding:(15)δτ(u)=u1|u|>τ,
and the soft thresholding:(16)δτ(u)=sgn⁡(u)(|u|-τ)+,
for a fixed threshold level τ(>0), where 1A is the indicator function of an event A, sgn(u) is the sign function of u and (u)+=max(0,u). There are many methods to determine the threshold level τ. In this paper, we use the universal threshold [21] and the “leave-out-half” cross-validation threshold [22]:(17)τ=2log⁡n,τ=(1-log⁡2log⁡n)-1/2τ(n2),
where n is the length of the observation si. Hard thresholding is a “keep” or “kill” rule, while soft thresholding is a “shrink” or “kill” rule. Both thresholding methods and both threshold levels will be employed to work on the empirical wavelet coefficients dj,k for denoising.


## 4.1.1. Haar Wavelet Transform

The Haar scaling function and the Haar wavelet function are defined as(10)ϕ(i)={1(0≤i<1)0(otherwise),ψ(i)={1(0≤i<12)-1(12≤i<1)0(otherwise),
respectively. By introducing the scaling parameterj and shifting parameterk, the Haar father wavelet and the Haar mother wavelet are defined by(11)ϕj,k(i)=2-j/2ϕ(2-ji-k),ψj,k(i)=2-j/2ψ(2-ji-k),
respectively. Then the target function, transformed software intensity function λi′(i=1,2,…,n), can be expressed in the following equation:(12)λi′=∑k=02j0-1αj0,kϕj0,k(i)+∑j=j0∞∑k=02j-1βj,kψj,k(i),
where(13)αj0,k=∑i=1nλi′ϕj0,k(i),(14)βj,k=∑i=1nλi′ψj,k(i)
are called the scaling coefficients and the wavelet coefficients, respectively, for any primary resolution level j0(≥0). Due to the implementability, it is reasonable to set an upper limit instead of ∞ for the resolution level j. In other words, the highest resolution level must be finite in practice. We use J to denote the highest resolution level. That is, the range of j in the second term of (12) is j∈[j0,J]. The mapping from function λi′ to coefficients (αj0,k,βj,k) is called the Haar wavelet transform (HWT).Sincebi,si or fi(i=1,2,…,n) can be considered as the observation of λi′, the empirical scaling coefficients cj0,k and the empirical wavelet coefficients dj,k of λi′ can be calculated by (13) and (14) with λi′ replaced by bi,si or fi. The noises involved in the empirical wavelet coefficients dj,k can be removed by the thresholding method that we will introduce later. Finally, the estimate of λi′ can be obtained by taking the inverse HWT with denoised empirical coefficients.


## 4.1.2. Thresholding

In denoising the empirical wavelet coefficients, the common choices of thresholding method are the hard thresholding:(15)δτ(u)=u1|u|>τ,
and the soft thresholding:(16)δτ(u)=sgn⁡(u)(|u|-τ)+,
for a fixed threshold level τ(>0), where 1A is the indicator function of an event A, sgn(u) is the sign function of u and (u)+=max(0,u). There are many methods to determine the threshold level τ. In this paper, we use the universal threshold [21] and the “leave-out-half” cross-validation threshold [22]:(17)τ=2log⁡n,τ=(1-log⁡2log⁡n)-1/2τ(n2),
where n is the length of the observation si. Hard thresholding is a “keep” or “kill” rule, while soft thresholding is a “shrink” or “kill” rule. Both thresholding methods and both threshold levels will be employed to work on the empirical wavelet coefficients dj,k for denoising.


## 4.2. Wavelet Shrinkage Estimation for NHPP-Based SRM

Since the software-fault count data is Poisson data, the preprocessing is necessary before making use of the Haar-wavelet-based denoising procedure. Xiao and Dohi [15] combined data transformation and the standard denoising procedure to propose the wavelet shrinkage estimation (WSE) for the D-NHPP-based SRMs. They call the HWT combined with AT, the Haar-Anscombe transform (HAT). Similarly, we call the HWT combined with BT and FT, the Haar-Bartlett transform (HBT) and the Haar-Freeman transform (HFT), respectively. In the numerical study, we investigate the goodness-of-fit performance of HBT- and HFT-based WSEs and compare them with the HAT-based WSE [15].


## 5. Numerical Study

### 5.1. Data Sets and Measures

We use six real project data sets cited from reference [1], where they are named as J1, J3, DATA14, J5, SS1, and DATA8. These data sets are software-fault count data (group data). We rename them for convenience as DS1~DS6 in this paper. Let (n,xn) denote the pair of the final testing date and the total number of detected fault. Then these data sets are presented by (62, 133), (41, 351), (46, 266), (73, 367), (81, 461), and (111, 481), respectively. We employ the MSE (mean squares error) as the goodness-of-fit measures, where(18)MSE1=∑i=1n(Λi-xi)2n,MSE2=∑i=1n(λi-yi)2n.
Additionally, we calculate LL (Log Likelihood), which is defined as(19)LL=∑i=1n(xi-xi-1)ln⁡[Λ̂i-Λ̂i-1]-Λ̂n(θ)-∑i=1nln⁡[(xi-xi-1)!],
where Λ̂i=∑k=1iλ̂k(i=1,2,…,n), and λ̂i is the WSE estimate of the software intensity function λi.


### 5.2. Goodness-of-Fit Test

A total of 16 wavelet-based estimation methods are examined in this paper since the WSE is applied with four thresholding techniques: hard thresholding (h) versus soft thresholding (s); universal threshold (ut) versus “leave-out-half,” cross-validation threshold (lht). Let HBT(·,·), HAT(·,·) and HFT(·,·) denote the WSEs based on Haar-Bartlett Transform, Haar-Anscombe Transform and Haar-Freeman Transform, respectively. HBT1(·,·), and HBT2(·,·) correspond to the transforms in (5) and (7), respectively. Additionally, the result of HAT-based WSE was introduced in [15], but they only showed the results of HAT-based WSE with hard thresholding. Here, we present comprehensively all the results of them for a further discussion.We present the goodness-of-fit results based on different threshold levels in Tables3 and 4, respectively. HBT2 provides smaller MSE and larger MLL than the others when “ut” is used, regardless of the thresholding method employed. It is worth mentioning that “ut” provides the same estimates with hard or soft thresholding in three data sets (DS1, DS2, and DS4). This is due to the relatively large value of “ut”, since when threshold is set to be 0, the output of thresholding δτ(dj,k) is the wavelet coefficient dj,k itself. In our numerical experiments, “ut” is relatively large in these 3 software-fault count data sets, which result in δut(dj,k)=0 whichever hard or soft thresholding is applied. HBT2 also looks better than the others when software thresholding is applied with “lht.” However, when “lht” is combined with hard thresholding, HBT1 (HAT; HFT) possesses the best results in DS3 and DS5 (DS1; DS6), respectively. Since “lht” is considered as a much more proper threshold level than “ut” in analyzing of software-fault count data, we suggest that HBT2 should be selected as an appropriate DT for the WSE.


### 5.3. Comparison with MLE

Our concern in this section is the comparison of WSE with MLE. We estimate the software intensity function by using three parametric D-NHPP-based SRMs listed in Table2, where the MLE is applied to estimate the model parameters for comparison. HBT2-based WSE is selected as a representative one among the 16 wavelet-based estimation methods. Figures 1 and 2 depict the estimation behavior of cumulative number of software faults and its increment per testing date (individual number of software faults detected per testing date) with DS1. The observed data is plotted in bar graph in both figures. Looking at (i) and (iii) in these figures, it is clear that the parametric D-NHPP-based SRMs with maximum likelihood estimator can fit the real data. However, since the parametric D-NHPP-based SRMs assume the software intensity function as smooth function, they can estimate only the average tendency of the individual number of software faults detected at each testing date, but cannot follow the microscopic fluctuated behavior in (ii) and (iv) of Figures 1 and 2. In other words, the estimation accuracy based on the cumulative number of faults is embedded in “cumulative effects” in (i) and (iii). The experimental results performed here give the potential applicability of the wavelet-based estimation methods with different thresholding schemes. Our methods employed here do not need the expensive computational cost comparing with the MLE (within less than one second to get an estimate). This is a powerful advantage in applying the D-NHPP-based SRMs, in addition to the fact that practitioners do not request much time and effort to implement the wavelet-based estimation algorithms.Table 2
Representative discrete NHPP-based SRMs.
SRMMean value functionΛiGeometric (GE)ω{1-(1-p)i}Negative Binomial (NB)ω{∑k=1i((r+k-1)!/(r-1)!k!)pr(1-p)k-1}Discrete Weibull (DW)ω{1-pir}Table 3
Goodness-of-fit test results for different data transformations. (Threshold level: universal threshold.)
DS1MSE1MSE2LLHBT1(h, ut)4.1800.327−153.759HBT2(h, ut)2.4010.316−144.850HAT(h, ut)2.5730.317−145.390HFT(h, ut)9.8420.403−383.785HBT1(s, ut)4.1800.327−153.759HBT2(s, ut)2.4010.316−144.850HAT(s, ut)2.5730.317−145.390HFT(s, ut)3.2310.320−148.023DS2MSE1MSE2LLHBT1(h, ut)12.4491.141−197.109HBT2(h, ut)11.6001.138−196.904HAT(h, ut)11.7251.138−196.953HFT(h, ut)11.8471.138−196.806HBT1(s, ut)12.4501.141−197.109HBT2(s, ut)11.6001.138−196.904HAT(s, ut)11.7251.138−196.953HFT(s, ut)11.8471.138−196.806DS3MSE1MSE2LLHBT1(h, ut)10.1590.719−143.181HBT2(h, ut)9.0550.791−192.903HAT(h, ut)9.2860.793−193.178HFT(h, ut)9.6530.796−167.025HBT1(s, ut)12.9541.037−212.321HBT2(s, ut)10.5281.035−213.978HAT(s, ut)10.8131.035−213.496HFT(s, ut)11.4961.037−213.419DS4MSE1MSE2LLHBT1(h, ut)2.0650.304−171.605HBT2(h, ut)1.7810.304−171.491HAT(h, ut)1.8220.304−171.504HFT(h, ut)1.8640.304−171.516HBT1(s, ut)2.0650.304−171.605HBT2(s, ut)1.7810.304−171.491HAT(s, ut)1.8220.304−171.504HFT(s, ut)1.8640.304−171.516DS5MSE1MSE2LLHBT1(h, ut)14.0130.586−304.198HBT2(h, ut)12.7810.581−301.073HAT(h, ut)12.9410.581−301.447HFT(h, ut)13.2000.582−301.841HBT1(s, ut)14.0540.601−312.551HBT2(s, ut)12.8040.596−309.709HAT(s, ut)12.9680.597−310.005HFT(s, ut)13.2270.598−310.480DS6MSE1MSE2LLHBT1(h, ut)12.7070.521−378.932HBT2(h, ut)11.8990.525−380.497HAT(h, ut)12.0530.526−380.987HFT(h, ut)12.1860.526−381.071HBT1(s, ut)13.5460.585−423.335HBT2(s, ut)12.5050.584−421.480HAT(s, ut)12.6700.584−421.729HFT(s, ut)12.8100.584−421.866Table 4
Goodness-of-fit test results for different data transformations. (Threshold level: “leave-out-half” cross-validation threshold.)
DS1MSE1MSE2LLHBT1(h, lht)0.1120.015−59.704HBT2(h, lht)0.1310.017−60.634HAT(h, lht)0.1140.014−60.324HFT(h, lht)8.0790.343−216.320HBT1(s, lht)0.1590.010−59.581HBT2(s, lht)0.1090.010−60.195HAT(s, lht)0.1150.010−60.118HFT(s, lht)0.1510.010−59.609DS2MSE1MSE2LLHBT1(h, lht)0.4770.102−70.438HBT2(h, lht)0.2250.072−70.227HAT(h, lht)0.2630.077−70.663HFT(h, lht)0.5520.125−71.041HBT1(s, lht)0.3420.027−69.407HBT2(s, lht)0.3630.030−69.781HAT(s, lht)0.3640.030−69.731HFT(s, lht)0.3510.028−69.446DS3MSE1MSE2LLHBT1(h, lht)0.2080.032−55.829HBT2(h, lht)0.2750.041−56.767HAT(h, lht)0.9020.074−57.614HFT(h, lht)0.3230.048−56.411HBT1(s, lht)0.5730.040−55.847HBT2(s, lht)0.4860.037−56.649HAT(s, lht)0.4960.037−56.551HFT(s, lht)0.5370.037−55.922DS4MSE1MSE2LLHBT1(h, lht)0.2880.015−121.008HBT2(h, lht)0.1320.007−120.942HAT(h, lht)0.2320.013−121.031HFT(h, lht)0.2810.014−121.004HBT1(s, lht)0.0970.005−120.893HBT2(s, lht)0.0460.003−120.906HAT(s, lht)0.0500.003−120.907HFT(s, lht)0.0640.004−120.888DS5MSE1MSE2LLHBT1(h, lht)0.0550.012−120.869HBT2(h, lht)0.2210.016−121.205HAT(h, lht)0.2280.017−121.195HFT(h, lht)0.2190.017−120.890HBT1(s, lht)0.1740.007−120.805HBT2(s, lht)0.1490.006−120.913HAT(s, lht)0.1510.006−120.898HFT(s, lht)0.1580.006−120.806DS6MSE1MSE2LLHBT1(h, lht)0.7780.025−166.906HBT2(h, lht)0.5580.017−166.816HAT(h, lht)0.5480.017−166.767HFT(h, lht)0.5400.017−166.673HBT1(s, lht)0.2880.009−166.432HBT2(s, lht)0.2820.009−166.532HAT(s, lht)0.2890.010−166.540HFT(s, lht)0.2940.009−166.440Behavior of estimates with MLE and WSE using hard thresholding (DS1).
(a)
HBT2-based WSE(h, lht) versus MLE(b)
HBT2-based WSE(h, lht) versus MLEBehavior of estimates with MLE and WSE using soft thresholding (DS1).
(a)
HBT2-based WSE(s, lht) versus MLE(b)
HBT2-based WSE(s, lht) versus MLE


## 5.1. Data Sets and Measures

We use six real project data sets cited from reference [1], where they are named as J1, J3, DATA14, J5, SS1, and DATA8. These data sets are software-fault count data (group data). We rename them for convenience as DS1~DS6 in this paper. Let (n,xn) denote the pair of the final testing date and the total number of detected fault. Then these data sets are presented by (62, 133), (41, 351), (46, 266), (73, 367), (81, 461), and (111, 481), respectively. We employ the MSE (mean squares error) as the goodness-of-fit measures, where(18)MSE1=∑i=1n(Λi-xi)2n,MSE2=∑i=1n(λi-yi)2n.
Additionally, we calculate LL (Log Likelihood), which is defined as(19)LL=∑i=1n(xi-xi-1)ln⁡[Λ̂i-Λ̂i-1]-Λ̂n(θ)-∑i=1nln⁡[(xi-xi-1)!],
where Λ̂i=∑k=1iλ̂k(i=1,2,…,n), and λ̂i is the WSE estimate of the software intensity function λi.


## 5.2. Goodness-of-Fit Test

A total of 16 wavelet-based estimation methods are examined in this paper since the WSE is applied with four thresholding techniques: hard thresholding (h) versus soft thresholding (s); universal threshold (ut) versus “leave-out-half,” cross-validation threshold (lht). Let HBT(·,·), HAT(·,·) and HFT(·,·) denote the WSEs based on Haar-Bartlett Transform, Haar-Anscombe Transform and Haar-Freeman Transform, respectively. HBT1(·,·), and HBT2(·,·) correspond to the transforms in (5) and (7), respectively. Additionally, the result of HAT-based WSE was introduced in [15], but they only showed the results of HAT-based WSE with hard thresholding. Here, we present comprehensively all the results of them for a further discussion.We present the goodness-of-fit results based on different threshold levels in Tables3 and 4, respectively. HBT2 provides smaller MSE and larger MLL than the others when “ut” is used, regardless of the thresholding method employed. It is worth mentioning that “ut” provides the same estimates with hard or soft thresholding in three data sets (DS1, DS2, and DS4). This is due to the relatively large value of “ut”, since when threshold is set to be 0, the output of thresholding δτ(dj,k) is the wavelet coefficient dj,k itself. In our numerical experiments, “ut” is relatively large in these 3 software-fault count data sets, which result in δut(dj,k)=0 whichever hard or soft thresholding is applied. HBT2 also looks better than the others when software thresholding is applied with “lht.” However, when “lht” is combined with hard thresholding, HBT1 (HAT; HFT) possesses the best results in DS3 and DS5 (DS1; DS6), respectively. Since “lht” is considered as a much more proper threshold level than “ut” in analyzing of software-fault count data, we suggest that HBT2 should be selected as an appropriate DT for the WSE.


## 5.3. Comparison with MLE

Our concern in this section is the comparison of WSE with MLE. We estimate the software intensity function by using three parametric D-NHPP-based SRMs listed in Table2, where the MLE is applied to estimate the model parameters for comparison. HBT2-based WSE is selected as a representative one among the 16 wavelet-based estimation methods. Figures 1 and 2 depict the estimation behavior of cumulative number of software faults and its increment per testing date (individual number of software faults detected per testing date) with DS1. The observed data is plotted in bar graph in both figures. Looking at (i) and (iii) in these figures, it is clear that the parametric D-NHPP-based SRMs with maximum likelihood estimator can fit the real data. However, since the parametric D-NHPP-based SRMs assume the software intensity function as smooth function, they can estimate only the average tendency of the individual number of software faults detected at each testing date, but cannot follow the microscopic fluctuated behavior in (ii) and (iv) of Figures 1 and 2. In other words, the estimation accuracy based on the cumulative number of faults is embedded in “cumulative effects” in (i) and (iii). The experimental results performed here give the potential applicability of the wavelet-based estimation methods with different thresholding schemes. Our methods employed here do not need the expensive computational cost comparing with the MLE (within less than one second to get an estimate). This is a powerful advantage in applying the D-NHPP-based SRMs, in addition to the fact that practitioners do not request much time and effort to implement the wavelet-based estimation algorithms.Table 2
Representative discrete NHPP-based SRMs.
SRMMean value functionΛiGeometric (GE)ω{1-(1-p)i}Negative Binomial (NB)ω{∑k=1i((r+k-1)!/(r-1)!k!)pr(1-p)k-1}Discrete Weibull (DW)ω{1-pir}Table 3
Goodness-of-fit test results for different data transformations. (Threshold level: universal threshold.)
DS1MSE1MSE2LLHBT1(h, ut)4.1800.327−153.759HBT2(h, ut)2.4010.316−144.850HAT(h, ut)2.5730.317−145.390HFT(h, ut)9.8420.403−383.785HBT1(s, ut)4.1800.327−153.759HBT2(s, ut)2.4010.316−144.850HAT(s, ut)2.5730.317−145.390HFT(s, ut)3.2310.320−148.023DS2MSE1MSE2LLHBT1(h, ut)12.4491.141−197.109HBT2(h, ut)11.6001.138−196.904HAT(h, ut)11.7251.138−196.953HFT(h, ut)11.8471.138−196.806HBT1(s, ut)12.4501.141−197.109HBT2(s, ut)11.6001.138−196.904HAT(s, ut)11.7251.138−196.953HFT(s, ut)11.8471.138−196.806DS3MSE1MSE2LLHBT1(h, ut)10.1590.719−143.181HBT2(h, ut)9.0550.791−192.903HAT(h, ut)9.2860.793−193.178HFT(h, ut)9.6530.796−167.025HBT1(s, ut)12.9541.037−212.321HBT2(s, ut)10.5281.035−213.978HAT(s, ut)10.8131.035−213.496HFT(s, ut)11.4961.037−213.419DS4MSE1MSE2LLHBT1(h, ut)2.0650.304−171.605HBT2(h, ut)1.7810.304−171.491HAT(h, ut)1.8220.304−171.504HFT(h, ut)1.8640.304−171.516HBT1(s, ut)2.0650.304−171.605HBT2(s, ut)1.7810.304−171.491HAT(s, ut)1.8220.304−171.504HFT(s, ut)1.8640.304−171.516DS5MSE1MSE2LLHBT1(h, ut)14.0130.586−304.198HBT2(h, ut)12.7810.581−301.073HAT(h, ut)12.9410.581−301.447HFT(h, ut)13.2000.582−301.841HBT1(s, ut)14.0540.601−312.551HBT2(s, ut)12.8040.596−309.709HAT(s, ut)12.9680.597−310.005HFT(s, ut)13.2270.598−310.480DS6MSE1MSE2LLHBT1(h, ut)12.7070.521−378.932HBT2(h, ut)11.8990.525−380.497HAT(h, ut)12.0530.526−380.987HFT(h, ut)12.1860.526−381.071HBT1(s, ut)13.5460.585−423.335HBT2(s, ut)12.5050.584−421.480HAT(s, ut)12.6700.584−421.729HFT(s, ut)12.8100.584−421.866Table 4
Goodness-of-fit test results for different data transformations. (Threshold level: “leave-out-half” cross-validation threshold.)
DS1MSE1MSE2LLHBT1(h, lht)0.1120.015−59.704HBT2(h, lht)0.1310.017−60.634HAT(h, lht)0.1140.014−60.324HFT(h, lht)8.0790.343−216.320HBT1(s, lht)0.1590.010−59.581HBT2(s, lht)0.1090.010−60.195HAT(s, lht)0.1150.010−60.118HFT(s, lht)0.1510.010−59.609DS2MSE1MSE2LLHBT1(h, lht)0.4770.102−70.438HBT2(h, lht)0.2250.072−70.227HAT(h, lht)0.2630.077−70.663HFT(h, lht)0.5520.125−71.041HBT1(s, lht)0.3420.027−69.407HBT2(s, lht)0.3630.030−69.781HAT(s, lht)0.3640.030−69.731HFT(s, lht)0.3510.028−69.446DS3MSE1MSE2LLHBT1(h, lht)0.2080.032−55.829HBT2(h, lht)0.2750.041−56.767HAT(h, lht)0.9020.074−57.614HFT(h, lht)0.3230.048−56.411HBT1(s, lht)0.5730.040−55.847HBT2(s, lht)0.4860.037−56.649HAT(s, lht)0.4960.037−56.551HFT(s, lht)0.5370.037−55.922DS4MSE1MSE2LLHBT1(h, lht)0.2880.015−121.008HBT2(h, lht)0.1320.007−120.942HAT(h, lht)0.2320.013−121.031HFT(h, lht)0.2810.014−121.004HBT1(s, lht)0.0970.005−120.893HBT2(s, lht)0.0460.003−120.906HAT(s, lht)0.0500.003−120.907HFT(s, lht)0.0640.004−120.888DS5MSE1MSE2LLHBT1(h, lht)0.0550.012−120.869HBT2(h, lht)0.2210.016−121.205HAT(h, lht)0.2280.017−121.195HFT(h, lht)0.2190.017−120.890HBT1(s, lht)0.1740.007−120.805HBT2(s, lht)0.1490.006−120.913HAT(s, lht)0.1510.006−120.898HFT(s, lht)0.1580.006−120.806DS6MSE1MSE2LLHBT1(h, lht)0.7780.025−166.906HBT2(h, lht)0.5580.017−166.816HAT(h, lht)0.5480.017−166.767HFT(h, lht)0.5400.017−166.673HBT1(s, lht)0.2880.009−166.432HBT2(s, lht)0.2820.009−166.532HAT(s, lht)0.2890.010−166.540HFT(s, lht)0.2940.009−166.440Behavior of estimates with MLE and WSE using hard thresholding (DS1).
(a)
HBT2-based WSE(h, lht) versus MLE(b)
HBT2-based WSE(h, lht) versus MLEBehavior of estimates with MLE and WSE using soft thresholding (DS1).
(a)
HBT2-based WSE(s, lht) versus MLE(b)
HBT2-based WSE(s, lht) versus MLE


## 6. Concluding Remarks

In this paper, we have applied the wavelet-based techniques to estimate the software intensity function. Four data transformations were employed to preprocess the software-fault count data. Throughout the numerical evaluation, we could conclude that the wavelet-based estimation methods with Bartlett transform2Yi have much more potential applicability than the other data transformations to the software reliability assessment practice because practitioners are not requested to carry out troublesome procedures on model selection and to take care of computational efficiency such as judgment of convergence and selecting initial guess of parameters in the general purpose of optimization algorithms. Note that, the result obtained here does not mean that the other data transformations are not good because the performance evaluation was executed only through goodness-of-fit test. Although the prediction ability of the proposed methods is out of focus of this paper at the present, the predictive performance should be considered and compared in the future.



---

*Source: 524636-2012-05-13.xml*