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

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.

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 several software 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 Section 2, 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 time 𝑖=0,1,2,…. Let 𝑌_𝑖 and 𝑁_𝑖=∑^𝑖_𝑘=0𝑌_𝑘 denote the number of software faults detected at testing date 𝑖 and its cumulative value, where 𝑌₀=𝑁₀=0 is assumed without any loss of generality. The stochastic process {𝑁_𝑖∶𝑖=0,1,2,…} is said to be a discrete non-homogeneous Poisson process (D-NHPP) if the probability mass function at time 𝑖 is given byPr𝑁_𝑖=Λ=𝑚^𝑖𝑚𝑚!exp−Λ_𝑖,𝑚=0,1,2,…, (1) where Λ_𝑖=𝐸[𝑁_𝑖] is called the mean value function of a D-NHPP and means the expected cumulative number of software faults detected by testing date 𝑖. The function 𝜆_𝑖=Λ_𝑖−Λ_𝑖−1(𝑖≥1) is called the discrete intensity function and implies the expected number of faults detected at testing date 𝑖, say 𝜆_𝑖=𝐸[𝑌_𝑖].

MLE is one of the most commonly used parametric estimation method. Let 𝜃 denote the vector of parameters in the mean value function Λ_𝑖=Λ_𝑖(𝜃), and 𝑥_𝑖(𝑦_𝑖) denote the realization of 𝑁_𝑖(𝑌_𝑖). When 𝑛 software faults are detected, the log-likelihood function of a D-NHPP is given byℒℒℱ(𝜃)=^𝑛_𝑖=1𝑥_𝑖−𝑥_𝑖−1lnΛ_𝑖(𝜃)−Λ_𝑖−1(𝜃)−Λ_𝑛(𝜃)−^𝑛_𝑖=1ln𝑥_𝑖−𝑥_𝑖−1!, (2) where 𝑥_𝑖=∑^𝑖_𝑘=0𝑦_𝑘 and 𝑦₀=𝑥₀=0. Then, the maximum likelihood estimate of 𝜃, say 𝜃, is given by the solution of argmax_𝜽ℒℒℱ(𝜃). Therefore, the estimate of the software intensity function 𝜆_𝑖=𝜆_𝑖(𝜃)(𝑖=0,1,2,…) can be obtained by 𝜆_𝑖=Λ_𝑖(𝜃)−Λ_𝑖−1(𝜃)(𝑖≥1) with 𝜆₀=0. In the following sections, we consider the problem of estimating the software intensity function from the noise-involved observation 𝑦_𝑖, 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 the Bartlett transform (BT) [18]. Let 𝜂_𝑖 denote the Poisson white noise, that is,𝑌_𝑖=𝜆_𝑖+𝜂_𝑖,𝑖=1,2…,𝑛. (3) Taking the BT, the random variables:𝐵_𝑖=2𝑌_𝑖,𝑖=1,2,…,𝑛 (4) can be approximately regarded as Gaussian random variables with the normal distribution 𝑁(2𝜆_𝑖−1/4,1), so that the realizations:𝑏_𝑖=2√𝑦_𝑖,𝑖=1,2,…,𝑛 (5) can be considered as samples from 𝑁(2𝜆_𝑖−1/4,1). That is, the transformed realizations 𝑏_𝑖(𝑖=1,2,…,𝑛) by the BT are the ones from the normally distributed random variables:𝐵_𝑖=𝜆^′_𝑖+𝜈_𝑖,𝑖=1,2,…,𝑛, (6) where 𝜆^′_𝑖=2𝜆_𝑖−1/4 is the transformed software intensity function, and 𝜈_𝑖 is the Gaussian white noise with unit variance.

Bartlett [18] also showed that𝑏_𝑖=2𝑦_𝑖+12,𝑖=1,2,…,𝑛 (7) is a better transformation since it provides a constant variance more closely to 1, even when the mean of 𝑌_𝑖 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:𝑠_𝑖=2𝑦_𝑖+38,𝑖=1,2,…,𝑛, (8) where 𝑠_𝑖 can be considered as observations of Gaussian random variable 𝑆_𝑖=2𝑌_𝑖+3/8 with the normal distribution 𝑁(2𝜆_𝑖+1/8,1). Freeman and Tukey [19] proposed the following square-root transform (we call it FT), which is also an extension of BT:𝑓_𝑖=𝑦_𝑖+1+√𝑦_𝑖,𝑖=1,2,…,𝑛. (9) They showed that the variance of Gaussian random variable 𝐹_𝑖=2𝑌_𝑖+1+𝑌_𝑖 is the nearest to 1 among BT, AT, and FT if the mean of 𝑌_𝑖 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.

As mentioned in Section 1, 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 𝑏_𝑖,𝑠_𝑖 or 𝑓_𝑖(𝑖=1,2,…,𝑛). Letting 𝑏^′_𝑖,𝑠^′_𝑖, and 𝑓^′_𝑖 denote the denoised 𝑏_𝑖,𝑠_𝑖, and 𝑓_𝑖, respectively, the estimate of the original software intensity function 𝜆_𝑖 can be obtained by taking the inverse DT of 𝑏^′_𝑖,𝑠^′_𝑖 and 𝑓^′_𝑖, 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𝜙(𝑖)=1()0(),𝜓(𝑖)=110≤𝑖<1otherwise0≤𝑖<21−120(),≤𝑖<1otherwise (10) respectively. By introducing the scaling parameter 𝑗 and shifting parameter 𝑘, the Haar father wavelet and the Haar mother wavelet are defined by𝜙_𝑗,𝑘(𝑖)=2^−𝑗/2𝜙2^−𝑗,𝜓𝑖−𝑘_𝑗,𝑘(𝑖)=2^−𝑗/2𝜓2^−𝑗𝑖−𝑘, (11) respectively. Then the target function, transformed software intensity function 𝜆^′_𝑖(𝑖=1,2,…,𝑛), can be expressed in the following equation:𝜆^′_𝑖=^2𝑗0−1_𝑘=0𝛼_𝑗0,𝑘𝜙_𝑗0,𝑘(𝑖)+^∞_𝑗=𝑗0^2𝑗−1_𝑘=0𝛽_𝑗,𝑘𝜓_𝑗,𝑘(𝑖), (12) where𝛼_𝑗0,𝑘=^𝑛_𝑖=1𝜆^′_𝑖𝜙_𝑗0,𝑘(𝑖), (13)𝛽_𝑗,𝑘=^𝑛_𝑖=1𝜆^′_𝑖𝜓_𝑗,𝑘(𝑖) (14) are called the scaling coefficients and the wavelet coefficients, respectively, for any primary resolution level 𝑗₀(≥0). Due to the implementability, it is reasonable to set an upper limit instead of ∞ for the resolution level 𝑗. In other words, the highest resolution level must be finite in practice. We use 𝐽 to denote the highest resolution level. That is, the range of 𝑗 in the second term of (12) is 𝑗∈[𝑗₀,𝐽]. The mapping from function 𝜆^′_𝑖 to coefficients (𝛼_𝑗0,𝑘,𝛽_𝑗,𝑘) is called the Haar wavelet transform (HWT).

Since 𝑏_𝑖,𝑠_𝑖 or 𝑓_𝑖(𝑖=1,2,…,𝑛) can be considered as the observation of 𝜆^′_𝑖, the empirical scaling coefficients 𝑐_𝑗0,𝑘 and the empirical wavelet coefficients 𝑑_𝑗,𝑘 of 𝜆^′_𝑖 can be calculated by (13) and (14) with 𝜆^′_𝑖 replaced by 𝑏_𝑖,𝑠_𝑖 or 𝑓_𝑖. The noises involved in the empirical wavelet coefficients 𝑑_𝑗,𝑘 can be removed by the thresholding method that we will introduce later. Finally, the estimate of 𝜆^′_𝑖 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:𝛿_𝜏(𝑢)=𝑢1_|𝑢|>𝜏, (15) and the soft thresholding:𝛿_𝜏(𝑢)=sgn(𝑢)(|𝑢|−𝜏)₊, (16) for a fixed threshold level 𝜏(>0), where 1_𝐴 is the indicator function of an event 𝐴, sgn(𝑢) is the sign function of 𝑢 and (𝑢)₊=max(0,𝑢). 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]:𝜏=2log𝑛,𝜏=1−log2log𝑛^−1/2𝜏𝑛2, (17) where 𝑛 is the length of the observation 𝑠_𝑖. 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 𝑑_𝑗,𝑘 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 (𝑛,𝑥_𝑛) 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, whereMSE₁=∑^𝑛_𝑖=1Λ_𝑖−𝑥_𝑖²𝑛,MSE₂=∑^𝑛_𝑖=1𝜆_𝑖−𝑦_𝑖²𝑛. (18) Additionally, we calculate LL (Log Likelihood), which is defined asLL=^𝑛_𝑖=1𝑥_𝑖−𝑥_𝑖−1lnΛ_𝑖−Λ_𝑖−1−Λ_𝑛(𝜃)−^𝑛_𝑖=1ln𝑥_𝑖−𝑥_𝑖−1!, (19) where Λ_𝑖=∑^𝑖_𝑘=1𝜆_𝑘(𝑖=1,2,…,𝑛), and 𝜆_𝑖 is the WSE estimate of the software intensity function 𝜆_𝑖.

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 Tables 3 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 𝛿_𝜏(𝑑_𝑗,𝑘) is the wavelet coefficient 𝑑_𝑗,𝑘 itself. In our numerical experiments, “ut” is relatively large in these 3 software-fault count data sets, which result in 𝛿_ut(𝑑_𝑗,𝑘)=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 Table 2, 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.

Table 3 

Goodness-of-fit test results for different data transformations. (Threshold level: universal threshold.)

Table 4 

Goodness-of-fit test results for different data transformations. (Threshold level: “leave-out-half” cross-validation threshold.)

(a) HBT2-based WSE(h, lht) versus MLE

(b) HBT2-based WSE(h, lht) versus MLE

Figure 1 

Behavior of estimates with MLE and WSE using hard thresholding (DS1).

(a) HBT2-based WSE(s, lht) versus MLE

(b) HBT2-based WSE(s, lht) versus MLE

Figure 2 

Behavior of estimates with MLE and WSE using soft thresholding (DS1).

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 transform 2𝑌_𝑖 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.