A Testing Coverage Model Based on NHPP Software Reliability Considering the Software Operating Environment and the Sensitivity Analysis

Abstract

We have been attempting to evaluate software quality and improve its reliability. Therefore, research on a software reliability model was part of the effort. Currently, software is used in various fields and environments; hence, one must provide quantitative confidence standards when using software. Therefore, we consider the testing coverage and uncertainty or randomness of an operating environment. In this paper, we propose a new testing coverage model based on NHPP software reliability with the uncertainty of operating environments, and we provide a sensitivity analysis to study the impact of each parameter of the proposed model. We examine the goodness-of-fit of a new testing coverage model based on NHPP software reliability and other existing models based on two datasets. The comparative results for the goodness-of-fit show that the proposed model does significantly better than the existing models. In addition, the results for the sensitivity analysis show that the parameters of the proposed model affect the mean value function.

1. Introduction

Software combines technologies such as artificial intelligence (AI), Internet of Things (IOT), and big data, which are important elements in the fourth industrial revolution, to create new forms of value. Software is a very important factor in the fourth industrial revolution [1]. Software reliability is defined as the probability that a software will run error free for a period of time. It is very difficult and complex to develop the new technologies and theories necessary to improve it. Therefore, the focus of software development is to improve the reliability and stability of software systems. In addition, unpredictable results occur during the software development process. Generally, the software development process consists of four stages: specification, design, coding, and testing [1]. Software faults are detected and corrected in the testing phase, which is the final stage of software development. Because the number of software faults and the interval between faults has a significant impact on software stability, software failure prediction is an important area of study for software developers, enterprises, and research institutions. The software reliability model makes it easier to evaluate software reliability using the fault data collected in a test or live environment. Furthermore, the software reliability model can measure the number of software faults, fault interval, and reliability, and the fault detection rate can be estimated and variously predicted. Conversely, testing coverage is one of the most important issues not solved in the software development process, and it is important for both software development and the customers of software products. Testing coverage is a measure that enables software developers to evaluate the quality of the software tested and determine how much additional effort is needed to improve the reliability of the software [2,3]. Testing coverage can provide customers with a quantitative confidence criterion when they plan to buy or use software products [4,5]. In addition, the software is used in a variety of operating environments. Therefore, the uncertainty or randomness of the software operating environment must be considered.

In the past, various statistical models have been proposed for the purpose of evaluating software reliability. The model based on the non-homogeneous Poisson process (NHPP) proved to be a very successful approach for practical software reliability [6]. On the basis of the NHPP, the mean value function can be used to obtain the expected number of faults to a certain point in time. Various NHPP software reliability models have been proposed to date. In the early days, most NHPP software reliability models were developed based on the assumption that faults detected in the testing phase were removed immediately with no debugging time delay, no new faults were introduced, and software systems used in the field environments were the same as or close to those used in the development-testing environment [7,8,9,10]. In the mid-1990s, studies began on a software reliability model consistent with a variety of software operating environments due to rapid changes in the industrial structure and environment. In the early 2000s, researchers began to explore new approaches such as the application of calibration factors to the software reliability model considering the uncertainty of the operating environment [11,12,13]. Recently, an NHPP software reliability model considering the uncertainty of the software operating environment has been proposed [14,15,16,17,18,19]. In addition, many testing coverage functions have been proposed in terms of different distributions, and software reliability models based on different testing coverage function have also been developed [2,20,21,22].

In this paper, we discuss a new testing coverage model based on NHPP software reliability with the uncertainty of operating environments and sensitivity analysis in order to study the impact of each parameter of the proposed model. We examine the goodness-of-fit of a new testing coverage model based on NHPP software reliability and other existing NHPP models based on two datasets. The explicit solution of the mean value function for a new testing coverage model is derived in Section 2. The various criteria for comparative model analysis are discussed in Section 3. Model analysis and results based on two actual datasets are discussed in Section 4. The impact of each parameter of the proposed model based on sensitivity analysis is discussed in Section 5. Finally, Section 6 presents the conclusions and remarks.

2. Testing Coverage Model based on NHPP Software Reliability

In this study, the basic assumption is to utilize the NHPP to describe the failure phenomenon during the testing phase, and the counting process $\mathrm{N}\left(\mathrm{t}\right)$ of the NHPP represents the cumulative number of failures up to the point of execution time t.

$$\mathrm{Pr}\left\{\mathrm{N}\left(\mathrm{t}\right)=\mathrm{n}\right\}=\frac{{\left\{\mathrm{m}\left(\mathrm{t}\right)\right\}}^{\mathrm{n}}}{\mathrm{n}!}\exp \left\{-\mathrm{m}\left(\mathrm{t}\right)\right\},\mathrm{n}=0,1,2,\dots ,\mathrm{t}\ge 0.$$

The mean value function $\mathrm{m}\left(\mathrm{t}\right)$ is

$$\mathrm{m}\left(\mathrm{t}\right)={{\displaystyle \int }}_0^{\mathrm{t}}\mathsf{\lambda }(\mathrm{s})\mathrm{d},$$

where, $\mathsf{\lambda }\left(\mathrm{s}\right)$ is the intensity function.

2.1. A General Testing Coverage Model based on NHPP Software Reliability

A general mean value function $\mathrm{m}\left(\mathrm{t}\right)$ of the testing coverage models based on NHPP software reliability using the differential equation is as follows [2]

$$\frac{\mathrm{d}\mathrm{m}\left(\mathrm{t}\right)}{\mathrm{dt}}=\frac{{\mathrm{c}}^{\prime }\left(\mathrm{t}\right)}{1-\mathrm{c}\left(\mathrm{t}\right)}\left[\mathrm{a}\left(\mathrm{t}\right)-\mathrm{m}\left(\mathrm{t}\right)\right],$$

(1)

where, $\mathrm{a}\left(\mathrm{t}\right)$ is the total number of faults detected in the software by time $\mathrm{t}$, $\mathrm{c}\left(\mathrm{t}\right)$ is the testing coverage function, i.e., the percentage of the code coverage by time $\mathrm{t}$, ${\mathrm{c}}^{\prime }\left(\mathrm{t}\right)$ is the derivative of the testing coverage function.

Solving Equation (1) using different functions $\mathrm{a}\left(\mathrm{t}\right)$ and $\frac{{\mathrm{c}}^{\prime }\left(\mathrm{t}\right)}{1-\mathrm{c}\left(\mathrm{t}\right)}$ yields the following mean value function $\mathrm{m}\left(\mathrm{t}\right)$ [2]:

$$\mathrm{m}\left(\mathrm{t}\right)={\mathrm{e}}^{-\mathrm{B}\left(\mathrm{t}\right)}\left[{\mathrm{m}}_0+{{\displaystyle \int }}_{{\mathrm{t}}_0}^{\mathrm{t}}\mathrm{a}\left(\mathsf{\tau }\right)\frac{{\mathrm{c}}^{\prime }\left(\mathsf{\tau }\right)}{1-\mathrm{c}\left(\mathsf{\tau }\right)}{\mathrm{e}}^{\mathrm{B}\left(\mathsf{\tau }\right)}\mathrm{b}\mathsf{\tau }\right],$$

(2)

where, $\mathrm{B}\left(\mathrm{t}\right)={{\displaystyle \int }}_{{\mathrm{t}}_0}^{\mathrm{t}}\frac{{\mathrm{c}}^{\prime }\left(\mathrm{s}\right)}{1-\mathrm{c}\left(\mathrm{s}\right)}\mathrm{ds}$ and $\mathrm{m}\left({\mathrm{t}}_0\right)={\mathrm{m}}_0$ is the marginal condition of Equation (2), with ${\mathrm{t}}_0$ representing the start time of the testing process.

2.2. A General Testing Coverage Model Based on NHPP Software Reliability with the Uncertainty of the Operating Environments

A general mean value function $\mathrm{m}\left(\mathrm{t}\right)$ of the testing coverage models based on NHPP software reliability using the differential equation considering the uncertainty of the operating environments is as follows [14]:

$$\frac{\mathrm{d}\mathrm{m}\left(\mathrm{t}\right)}{\mathrm{dt}}=\mathsf{\eta }\left[\frac{{\mathrm{c}}^{\prime }\left(\mathrm{t}\right)}{1-\mathrm{c}\left(\mathrm{t}\right)}\right]\left[\mathrm{a}\left(\mathrm{t}\right)-\mathrm{m}\left(\mathrm{t}\right)\right].$$

(3)

We find the mean value function $\mathrm{m}\left(\mathrm{t}\right)$ as in Equation (4) using the differential equation of Equation (3) by applying the random variable $\mathsf{\eta }$ to a probability density function in the uncertain operating environments [14].

$$\mathrm{m}\left(\mathrm{t}\right)={\mathrm{e}}^{-\mathsf{\eta }\mathrm{B}\left(\mathrm{t}\right)}\left[{\mathrm{m}}_0+{{\displaystyle \int }}_{{\mathrm{t}}_0}^{\mathrm{t}}\mathrm{a}\left(\mathsf{\tau }\right)\frac{{\mathrm{c}}^{\prime }\left(\mathsf{\tau }\right)}{1-\mathrm{c}\left(\mathsf{\tau }\right)}{\mathrm{e}}^{\mathsf{\eta }\mathrm{B}\left(\mathsf{\tau }\right)}\mathrm{b}\mathsf{\tau }\right].$$

(4)

2.3. A New Testing Coverage Model based on NHPP Software Reliability Considering the Uncertainty of the Operating Environments

In this paper, a new testing coverage model based on NHPP software reliability considering the uncertainty of the operating environments is presented. We apply Equations (3) and (4), which are assumptions of the existing test coverage model, and add the following assumptions [3,14]:

▪

The initial condition of the mean value function $\mathrm{m}\left(\mathrm{t}\right)$ is $\mathrm{m}\left(0\right)=0$;

▪

$\mathrm{a}\left(\mathrm{t}\right)=\mathrm{N}$ is the expected number of faults that exist in the software before testing;

▪

$\mathsf{\eta }$ has a generalized probability density function g with two parameters $\mathsf{\alpha }$ and $\mathsf{\beta }$;

▪

The fault detection rate can be expressed by $\frac{{\mathrm{c}}^{\prime }\left(\mathsf{\tau }\right)}{1-\mathrm{c}\left(\mathsf{\tau }\right)}$.

We can derive the mean value function $\mathrm{m}\left(\mathrm{t}\right)$ based on the assumptions and differential equations [16].

$$\mathrm{m}\left(\mathrm{t}\right)={{\displaystyle \int }}_{\mathsf{\eta }}^{}\mathrm{N}\left(1-{\mathrm{e}}^{-\mathsf{\eta }{{\displaystyle \int }}_0^{\mathrm{t}}\frac{{\mathrm{c}}^{\prime }\left(\mathrm{s}\right)}{1-\mathrm{c}\left(\mathrm{s}\right)}\mathrm{ds}}\right)\mathrm{dg}\left(\mathsf{\eta }\right),$$

$$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}{\left(1-\frac{\mathsf{\beta }}{\mathsf{\beta }+{{\displaystyle \int }}_0^{\mathrm{t}}\frac{{\mathrm{c}}^{\prime }\left(\mathrm{s}\right)}{1-\mathrm{c}\left(\mathrm{s}\right)}\mathrm{ds}}\right)}^{\mathsf{\alpha }},\mathsf{\alpha }\ge 0,\mathsf{\beta }\ge 0.$$

(5)

In this study, we consider a testing coverage function $\mathrm{c}\left(\mathrm{t}\right)$ as follows:

$$\mathrm{c}\left(\mathrm{t}\right)=1-\left(1+\mathrm{dt}\right){\mathrm{e}}^{-\mathrm{bt}},\mathrm{b},\mathrm{d}>0,$$

(6)

where, $\mathrm{b}$ is the failure detection rate and $\mathrm{d}$ represents the shape factor.

We obtain a new mean value function m(t) for the testing coverage model based on NHPP software reliability with uncertainty in the operating environment that can be used to determine the expected number of software failures detected by time $\mathrm{t}$ by substituting the function c(t) above into Equation (5):

$$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}{\left(1-\frac{\mathsf{\beta }}{\mathsf{\beta }+\mathrm{bt}-\ln \left(1+\mathrm{dt}\right)}\right)}^{\mathsf{\alpha }}.$$

3. Various Criteria for Comparative Model Analysis

We estimate the parameters of the NHPP software reliability models in

Table 1. NHPP software reliability models.

No.	Model	$\mathbf{Mean}\mathbf{Value}\mathbf{Function}\mathbf{m}\left(\mathbf{t}\right)$
1	Goel-Okumoto (GO) [7]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left(1-{\mathrm{e}}^{-\mathrm{bt}}\right)$
2	Yamada et al. (DS) [8]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left(1-\left(1+\mathrm{bt}\right){\mathrm{e}}^{-\mathrm{bt}}\right)$
3	Ohba (IS) [9]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{a}\left(1-{\mathrm{e}}^{-\mathrm{bt}}\right)}{1+{\mathsf{\beta }\mathrm{e}}^{-\mathrm{bt}}}$
4	Yamada et al. (YE) [10]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left(1-{\mathrm{e}}^{-\mathsf{\gamma }\mathsf{\alpha }\left(1-{\mathrm{e}}^{-\mathsf{\beta }\mathrm{t}}\right)}\right)$
5	Yamata et al. (YR) [10]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left(1-{\mathrm{e}}^{-\mathsf{\gamma }\mathsf{\alpha }\left(1-{\mathrm{e}}^{-{\mathsf{\beta }\mathrm{t}}^2/2}\right)}\right)$
6	Yamada et al. (YID 1) [26]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{ab}}{\mathsf{\alpha }+\mathrm{b}}\left({\mathrm{e}}^{\mathsf{\alpha }\mathrm{t}}-{\mathrm{e}}^{-\mathrm{bt}}\right)$
7	Yamada et al. (YID 2) [26]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left[1-{\mathrm{e}}^{-\mathrm{bt}}\right]\left[1-\frac{\mathsf{\alpha }}{\mathrm{b}}\right]+\mathsf{\alpha }\mathrm{at}$
8	Hossain-Dahiya (HDGO) [27]	$\mathrm{m}\left(\mathrm{t}\right)=\log [({\mathrm{e}}^{\mathrm{a}}-\mathrm{c})/({\mathrm{e}}^{{\mathrm{ae}}^{-\mathrm{bt}}}-\mathrm{c})]$
9	Pham et al. (PNZ) [28]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{a}\left[1-{\mathrm{e}}^{-\mathrm{bt}}\right]\left[1-\frac{\mathsf{\alpha }}{\mathrm{b}}\right]+\mathsf{\alpha }\mathrm{at}}{1+{\mathsf{\beta }\mathrm{e}}^{-\mathrm{bt}}}$
10	Pham-Zhang (PZ) [29]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{1}{1+{\mathsf{\beta }\mathrm{e}}^{-\mathrm{bt}}}\ast \left(\left(\mathrm{c}+\mathrm{a}\right)\left[1-{\mathrm{e}}^{-\mathrm{bt}}\right]-\left[\frac{\mathrm{ab}}{\mathrm{b}-\mathsf{\alpha }}\left({\mathrm{e}}^{-\mathsf{\alpha }\mathrm{t}}-{\mathrm{e}}^{-\mathrm{bt}}\right)\right]\right)$
11	Zhang et al. (ZFR) [30]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{a}}{\mathrm{p}-\mathsf{\beta }}[1-{\left(\frac{\left(1+\mathsf{\alpha }\right){\mathrm{e}}^{-\mathrm{bt}}}{1+{\mathsf{\alpha }\mathrm{e}}^{-\mathrm{bt}}}\right)}^{\frac{\mathrm{c}}{\mathrm{b}}\left(\mathrm{p}-\mathsf{\beta }\right)}]$
12	Teng-Pham (TP) [14]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{a}}{\mathrm{p}-\mathrm{q}}[1-{\left(\frac{\mathsf{\beta }}{\mathsf{\beta }+\left(\mathrm{p}-\mathrm{q}\right)\ln \left(\frac{\mathrm{c}+{\mathrm{e}}^{\mathrm{bt}}}{\mathrm{c}+1}\right)}\right)}^{\mathsf{\alpha }}]$
13	Pham (IFD) [1]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\left(1-{\mathrm{e}}^{-\mathrm{bt}}\right)\left(1+\left(\mathrm{b}+\mathrm{d}\right)\mathrm{t}+{\mathrm{bdt}}^2\right)$
14	Pham (PDP) [31]	$\mathrm{m}\left(\mathrm{t}\right)={\mathrm{m}}_0\left(\frac{\mathsf{\gamma }\mathrm{t}+1}{{\mathsf{\gamma }\mathrm{t}}_0+1}\right){\mathrm{e}}^{-\mathsf{\gamma }\left(\mathrm{t}-{\mathrm{t}}_0\right)}+\mathsf{\alpha }\left(\mathsf{\gamma }\mathrm{t}+1\right)(\mathsf{\gamma }\mathrm{t}-1+\left(1-{\mathsf{\gamma }\mathrm{t}}_0\right){\mathrm{e}}^{-\mathsf{\gamma }\left(\mathrm{t}-{\mathrm{t}}_0\right)}$
15	Kapur et al. (KSRGM) [32]	$\mathrm{m}\left(\mathrm{t}\right)=\frac{\mathrm{A}}{1-\mathsf{\alpha }}[1-{\left(\left(1+\mathrm{bt}+\frac{{\mathrm{b}}^2{\mathrm{t}}^2}{2}\right){\mathrm{e}}^{-\mathrm{bt}}\right)}^{\mathrm{p}\left(1-\mathsf{\alpha }\right)}]$
16	Roy et al. (RMD) [33]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{a}\mathsf{\alpha }\left[1-{\mathrm{e}}^{-\mathrm{bt}}\right]-\left[\frac{\mathrm{ab}}{\mathrm{b}-\mathsf{\beta }}\left({\mathrm{e}}^{-\mathsf{\beta }\mathrm{t}}-{\mathrm{e}}^{-\mathrm{bt}}\right)\right]$
17	Chang et al. (CT) [3]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}\left[1-{\left(\frac{\mathsf{\beta }}{\mathsf{\beta }+{\left(\mathrm{at}\right)}^{\mathrm{b}}}\right)}^{\mathsf{\alpha }}\right]$
18	Pham (Vtub) [25]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}\left[1-{\left(\frac{\mathsf{\beta }}{\mathsf{\beta }+{\mathrm{a}}^{\mathrm{t}}{}^{\mathrm{b}}-1}\right)}^{\mathsf{\alpha }}\right]$
19	Song et al. (3PFD) [15]	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}\left[1-\left(\frac{\mathsf{\beta }}{\mathsf{\beta }-\frac{\mathrm{a}}{\mathrm{b}}\ln \left(\frac{\left(1+\mathrm{c}\right){\mathrm{e}}^{-\mathrm{bt}}}{1+{\mathrm{ce}}^{-\mathrm{bt}}}\right)}\right)\right]$
20	Proposed Model	$\mathrm{m}\left(\mathrm{t}\right)=\mathrm{N}{\left(1-\frac{\mathsf{\beta }}{\mathsf{\beta }+\mathrm{bt}-\ln \left(1+\mathrm{dt}\right)}\right)}^{\mathsf{\alpha }}$

using the least squares estimation (LSE) method. We derived the parameters of the mean value function $\mathrm{m}\left(\mathrm{t}\right)$ using Matlab (2016a, MathWorks, Natick, MA, USA) and R (Version 3.3.1) programs based on the LSE method.

3.1. Criteria for Model Comparison

The statistical significance of model comparisons can be confirmed through well-known criteria. We use twelve criteria to estimate the goodness-of-fit of all model and to compare the proposed model with other models in Table 1.

The mean squared error (MSE) measures the average of the squares of the errors that is the average squared difference between the estimated values and the actual data. The root mean square error (RMSE) is a frequently used measure of the differences between values predicted by a model or an estimator and the values observed. Akaike’s information criterion (AIC) is used to compare the capability of each model to maximize the likelihood function (L), while considering the degrees of freedom [23]. The variance measures the standard deviation of the prediction bias [19]. The root mean square prediction error (RMSPE) is a measure of the closeness with which the model predicts the observation [19]. The predictive ratio risk (PRR) measures the distance of the model estimates from the actual data over the model estimate [1]. The Theil statistic (TS) is the average deviation percentage over all periods with regard to the actual values [24]. The predictive power (PP) measures the distance of the model estimates from the actual data [1]. The sum of absolute errors (SAE) measures the absolute distance of the model [16]. The mean absolute errors (MAE) measures the deviation by using the absolute distance of the model [24]. The R² measures how successful the fit is in explaining the variations in the data [4]. The adjusted R² (Adj R²) is a modification to the R² that accounts for the number of explanatory terms in a model relative to the number of data points [4]. The smaller the value of these ten criteria, i.e., MSE, RMSE, AIC, variance, RMSPE, PRR, TS, PP, SAE, and MAE, the better the model fit is (close to 0). Conversely, the higher the value of these two criteria, i.e., ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$, the better the model fit is (close to 1). These criteria are described as follows in Appendix

Table A1. List of criteria for model comparisons.

No.	Criteria	Formula
1	MSE	$\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left(\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{\mathrm{n}-\mathrm{m}}$
2	RMSE	$\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left(\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{\mathrm{n}-\mathrm{m}}}$
3	AIC	$-2\log \mathrm{L}+2\mathrm{m}$
4	Variance	$\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-\mathrm{Bias}\right)}^2}{\mathrm{n}-1}}$
5	RMSPE	$\sqrt{{\mathrm{Variance}}^2+{\mathrm{Bias}}^2}$
6	PRR	${\displaystyle {\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}}{\left(\frac{\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}}{\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)}\right)}^2$
7	TS	$100\ast \sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{y}}_{\mathrm{i}}{}^2}}\%$
8	PP	${\displaystyle {\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}}{\left(\frac{\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}}{{\mathrm{y}}_{\mathrm{i}}}\right)}^2$
9	SAE	${{\displaystyle \sum }}_{\mathrm{i}=0}^{\mathrm{n}}\left|\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}\right|$
10	MAE	$\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}\right|}{\mathrm{n}-\mathrm{m}}$
11	${\mathrm{R}}^2$	$1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=0}^{\mathrm{n}}{\left(\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)-{\mathrm{y}}_{\mathrm{i}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=0}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}$
12	Adj ${\mathrm{R}}^2$	$1-\frac{\left(1-{\mathrm{R}}^2\right)\left(\mathrm{n}-1\right)}{\mathrm{n}-\mathrm{m}-1}$

.

From Table A1, $\hat{\mathrm{m}}\left({\mathrm{t}}_{\mathrm{i}}\right)$ is the estimated cumulative number of failures at ${\mathrm{t}}_{\mathrm{i}}$ for $\mathrm{i}=1,2,\cdots ,\mathrm{n}$; ${\mathrm{y}}_{\mathrm{i}}$ is the total number of failures observed at time ${\mathrm{t}}_{\mathrm{i}}$; $\mathrm{n}$ is the total number of observations; and $\mathrm{m}$ is the number of unknown parameters in the model.

3.2. Distance of the Normalized Criteria

Pham [25] and Li and Pham [4] mentioned the distance of the normalized criteria (NCD) methods for ranking and selecting the best model from various software reliability models according to the characteristics of a set of criteria. We compare the performance of the software reliability models using the twelve criteria in Section 3.1. It is difficult to check the performance of the software reliability model due to the diversity in criteria, so a criteria method that integrates them is needed. It is shown that 10 out of 12 criteria show that the goodness-of-fit (similarity to actual data) improves as the value of the criteria is closer to 0, while for 2 of the 12 criteria, i.e., ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$, being closer to 1 indicates a better goodness-of-fit. In addition, the weight of each criterion is not considered because each criterion is different. We have improved the NCD method according to the characteristics of the criteria proposed in Section 3.1.

The value of NCD is defined as follows:

$${\mathrm{D}}_{\mathrm{k}}=\sqrt{{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{d}}{\left(\frac{{\mathrm{C}}_{\mathrm{kj}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{s}}{\mathrm{C}}_{\mathrm{ij}}}\right)}^2+{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{f}}{\left(\frac{1-{\mathrm{Z}}_{\mathrm{kj}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{s}}\left(1-{\mathrm{Z}}_{\mathrm{ij}}\right)}\right)}^2},\mathrm{k}=1,2,\dots ,\mathrm{s},$$

where, $\mathrm{s}$ is the total number of models; ${\mathrm{C}}_{\mathrm{ij}}$ denotes the value for the ith model of the jth criterion where $\mathrm{i}=1,2,\dots ,\mathrm{s}$; $\mathrm{d}$ denotes the total number of criteria, e.g., MSE, RMSE, AIC, variance, RMSPE, PRR, TS, PP, SAE, and MAE; $\mathrm{f}$ denotes the total number of criteria, e.g., ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$. Table 1 summarizes the mean value functions for NHPP software reliability models.

3.3. Confidence Interval

We use the following (7) to obtain the confidence interval [1] of the proposed model and of the existing NHPP software reliability models.

$$\hat{\mathrm{m}}\left(\mathrm{t}\right)\pm {\mathrm{Z}}_{\mathsf{\alpha }/2}\sqrt{\hat{\mathrm{m}}\left(\mathrm{t}\right)},$$

(7)

where, ${\mathrm{Z}}_{\mathsf{\alpha }/2}$ is $100\left(1-\mathsf{\alpha }\right)$, the percentile of the standard normal distribution. It is possible to confirm whether the value of the mean value function is included in the confidence interval at each time point or not and how much the confidence interval actually contains the value.

4. Data and the Result of Model Analysis

We estimate the parameters on two datasets, find the results of all criteria and the NCD value, and compare the results against each other.

4.1. Dataset 1

Dataset #1 was given by [6] and provided in

Table 2. Failure data-dataset #1.

Week Index	Dataset #1
	Failures	Cumulative Failures
1	2	2
2	11	13
3	2	15
4	4	19
5	3	22
6	1	23
7	1	24
8	2	26
9	4	30
10	0	30
11	4	34
12	1	35
13	3	38
14	0	38

. In dataset #1, the week index ranges from 1 week to 21 weeks, and there are 38 cumulative failures at 14 weeks. Detailed information is recorded in [6].

First, in

Table 3. Results of model parameter estimation and criteria from dataset #1.

No.	Model	Parameter Estimation	MSE	RMSE	AIC	Variation	RMSPE	PRR	TS	PP	SAE	MAE	R2	Adj R2
1	GO	$\hat{\mathrm{a}}$ = 46.1404 $\hat{\mathrm{b}}$ = 0.1182	3.6343	1.9064	62.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4374	1.7031	0.9687	0.9630
2	DS	$\hat{\mathrm{a}}$ = 35.9507 $\hat{\mathrm{b}}$ = 0.3988	9.1885	3.0312	68.6375	3.0081	3.0371	0.9541	10.4521	0.3575	32.6041	2.7170	0.9208	0.9064
3	IS	$\hat{\mathrm{a}}$ = 46.1403, $\hat{\mathrm{b}}$ = 0.1182 $\hat{\mathsf{\beta }}$ = 5.93 × 10−9	3.9647	1.9912	64.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4375	1.8580	0.9687	0.9593
4	YE	$\hat{\mathrm{a}}$ = 64.4549, $\hat{\mathsf{\alpha }}$ = 0.0183 $\hat{\mathsf{\beta }}$ = 0.0574, $\hat{\mathsf{\gamma }}$ = 80.1118	4.1162	2.0288	66.7083	1.7794	1.7794	0.5190	6.3861	2.8085	19.6361	1.9636	0.9704	0.9573
5	YR	$\hat{\mathrm{a}}$ = 38.1089, $\hat{\mathsf{\alpha }}$ = 1.6323 $\hat{\mathsf{\beta }}$ = 0.0364, $\hat{\mathsf{\gamma }}$ = 1.4194	16.0972	4.0121	82.0186	3.6707	3.7164	1.9186	12.6289	0.5463	39.5895	3.9590	0.8844	0.8331
6	YID 1	$\hat{\mathrm{a}}$ = 22.4227, $\hat{\mathrm{b}}$ = 0.3115 $\hat{\mathsf{\alpha }}$ = 0.0505	2.9149	1.7073	64.8745	1.5841	1.5883	0.5051	5.6364	4.3732	16.1998	1.4727	0.9770	0.9701
7	YID 2	$\hat{\mathrm{a}}$ = 18.5582, $\hat{\mathrm{b}}$ = 0.3808 $\hat{\mathsf{\alpha }}$ = 0.0952	2.9795	1.7261	64.7603	1.5991	1.6026	0.5092	5.6984	4.3985	16.2150	1.4741	0.9765	0.9694
8	HDGO	$\hat{\mathrm{a}}$ = 46.1403, $\hat{\mathrm{b}}$ = 0.1182 $\hat{\mathrm{c}}$ = 0.00022	3.9647	1.9912	64.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4375	1.8580	0.9687	0.9593
9	PNZ	$\hat{\mathrm{a}}$ = 17.7782, $\hat{\mathrm{b}}$ = 0.4256 $\hat{\mathsf{\alpha }}$ = 0.1008, $\hat{\mathsf{\beta }}$ = 0.0836	3.2378	1.7994	66.6674	1.5908	1.5947	0.5069	5.6639	4.3641	16.1109	1.6111	0.9768	0.9664
10	PZ	$\hat{\mathrm{a}}$ = 4.33 × 10−8, $\hat{\mathrm{b}}$ = 0.1182 $\hat{\mathsf{\alpha }}$ = 0.03204, $\hat{\mathsf{\beta }}$ = 6.35 × 10−8 $\hat{\mathrm{c}}$ = 46.14	4.8458	2.2013	68.8309	1.8327	1.8330	0.5283	6.5734	2.5742	20.4379	2.2709	0.9687	0.9491
11	ZFR	$\hat{\mathrm{a}}$ = 6.8991, $\hat{\mathrm{b}}$ = 0.0132 $\hat{\mathsf{\alpha }}$ = 0.0095, $\hat{\mathsf{\beta }}$ = 0.1242 $\hat{\mathrm{c}}$ = 0.7981, $\hat{\mathrm{p}}$ = 0.2738	5.4529	2.3351	70.8319	1.8327	1.8330	0.5282	6.5743	2.5753	20.4372	2.5546	0.9687	0.9418
12	TP	$\hat{\mathrm{a}}$ = 78.530, $\hat{\mathrm{b}}$ = 0.3843 $\hat{\mathsf{\alpha }}$ = 0.06219, $\hat{\mathsf{\beta }}$ = 0.2910 $\hat{\mathrm{c}}$ = 3.19 × 10−8, $\hat{\mathrm{p}}$ = 0.8530 $\hat{\mathrm{q}}$ = 0.6711	5.2511	2.2915	72.6973	1.6836	1.6843	0.5125	6.0348	3.5656	17.8269	2.5467	0.9736	0.9428
13	IFD	$\hat{\mathrm{a}}$ = 3.2868, $\hat{\mathrm{b}}$ = 0.8524 $\hat{\mathrm{d}}$ = 1.0 × 10−11	13.0533	3.6129	67.8928	3.8336	3.9783	1.2118	11.9274	0.9925	36.9137	3.3558	0.8969	0.8660
14	PDP	$\hat{\mathsf{\alpha }}$ = 1833.3580, $\hat{\mathsf{\gamma }}$ = 0.0119 ${\hat{\mathrm{t}}}_0$ = 28.7597, ${\hat{\mathrm{m}}}_0$ = 143.1294	26.2542	5.1239	138.4862	4.4940	4.4941	0.9330	16.1283	42.1676	43.5907	4.3591	0.8115	0.7277
15	KSRGM	$\hat{\mathrm{A}}$ = 0.9468, $\hat{\mathrm{b}}$ = 10.5138 $\hat{\mathsf{\alpha }}$ = 0.9774, $\hat{\mathrm{p}}$ = 0.6679	5.3204	2.3066	65.7919	2.0616	2.0734	0.5176	7.2604	0.9784	23.2980	2.3298	0.9618	0.9448
16	RMD	$\hat{\mathrm{a}}$ = 473.90, $\hat{\mathrm{b}}$ = 0.3889 $\hat{\mathsf{\alpha }}$ = 1.0380, $\hat{\mathsf{\beta }}$ = 0.003857	3.3303	1.8249	66.7588	1.6042	1.6053	0.5115	5.7443	4.2796	16.7058	1.6706	0.9761	0.9655
17	CT	$\hat{\mathrm{a}}$ = 0.3318, $\hat{\mathrm{b}}$ = 1.110 $\hat{\mathsf{\alpha }}$ = 0.03498, $\hat{\mathsf{\beta }}$ = 0.9627 $\hat{\mathrm{N}}$ = 584.00	4.0547	2.0136	68.0606	1.6777	1.6784	0.4961	6.0130	2.9906	18.1577	2.0175	0.9738	0.9574
18	Vtub	$\hat{\mathrm{a}}$ = 91.140, $\hat{\mathrm{b}}$ = 0.5612 $\hat{\mathsf{\alpha }}$ = 0.04156, $\hat{\mathsf{\beta }}$ = 21.950 $\hat{\mathrm{N}}$ = 74.030	3.7733	1.9425	66.2648	1.6196	1.6206	0.4394	5.8006	2.1000	18.4484	2.0498	0.9756	0.9604
19	3PFD	$\hat{\mathrm{a}}$ = 0.3353, $\hat{\mathrm{b}}$ = 1.46 × 10−7 $\hat{\mathsf{\beta }}$ = 0.18198, $\hat{\mathrm{N}}$ = 68.60676 $\hat{\mathrm{c}}$ = 20.31988	4.3488	2.0854	68.6419	1.7359	1.7361	0.5136	6.2272	3.0687	18.8258	2.0918	0.9719	0.9543
20	New Model	$\hat{\mathrm{b}}$ = 0.2627, $\hat{\mathsf{\alpha }}$ = 0.3487 $\hat{\mathsf{\beta }}$ = 193.7390, $\hat{\mathrm{N}}$ = 184.0433 $\hat{\mathrm{d}}$ = 0.2997	2.4340	1.5601	60.2184	1.2993	1.2997	0.0709	4.6588	0.0599	13.8037	1.5337	0.9843	0.9744

, we obtained the parameter estimates for all twenty models and the values of all twelve criteria using $\mathrm{t}=1,2,\cdots ,21$ of the week index from dataset #1. As shown in Table 3, we can see that the proposed model achieves the best results when comparing the twelve criteria to the other models. In addition, as shown in

Table 4. Results of criteria and NCD values for comparison with dataset #1.

No.	Model	MSE	RMSE	AIC	Variation	RMSPE	PRR	TS	PP	SAE	MAE	R2	Adj R2	NCD Value	Rank
1	GO	3.6343	1.9064	62.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4374	1.7031	0.9687	0.9630	0.1339749	8
2	DS	9.1885	3.0312	68.6375	3.0081	3.0371	0.9541	10.4521	0.3575	32.6041	2.7170	0.9208	0.9064	0.2307686	17
3	IS	3.9647	1.9912	64.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4375	1.8580	0.9687	0.9593	0.1371307	10
4	YE	4.1162	2.0288	66.7083	1.7794	1.7794	0.5190	6.3861	2.8085	19.6361	1.9636	0.9704	0.9573	0.1372862	11
5	YR	16.0972	4.0121	82.0186	3.6707	3.7164	1.9186	12.6289	0.5463	39.5895	3.9590	0.8844	0.8331	0.3432228	19
6	YID 1	2.9149	1.7073	64.8745	1.5841	1.5883	0.5051	5.6364	4.3732	16.1998	1.4727	0.9770	0.9701	0.1245350	2
7	YID 2	2.9795	1.7261	64.7603	1.5991	1.6026	0.5092	5.6984	4.3985	16.2150	1.4741	0.9765	0.9694	0.1254775	3
8	HDGO	3.9647	1.9912	64.8309	1.8327	1.8330	0.5283	6.5734	2.5743	20.4375	1.8580	0.9687	0.9593	0.1371307	9
9	PNZ	3.2378	1.7994	66.6674	1.5908	1.5947	0.5069	5.6639	4.3641	16.1109	1.6111	0.9768	0.9664	0.1275318	5
10	PZ	4.8458	2.2013	68.8309	1.8327	1.8330	0.5283	6.5734	2.5742	20.4379	2.2709	0.9687	0.9491	0.1458965	13
11	ZFR	5.4529	2.3351	70.8319	1.8327	1.8330	0.5282	6.5743	2.5753	20.4372	2.5546	0.9687	0.9418	0.1522740	15
12	TP	5.2511	2.2915	72.6973	1.6836	1.6843	0.5125	6.0348	3.5656	17.8269	2.5467	0.9736	0.9428	0.1481766	14
13	IFD	13.0533	3.6129	67.8928	3.8336	3.9783	1.2118	11.9274	0.9925	36.9137	3.3558	0.8969	0.8660	0.2925481	18
14	PDP	26.2542	5.1239	138.4862	4.4940	4.4941	0.9330	16.1283	42.1676	43.5907	4.3591	0.8115	0.7277	0.6513642	20
15	KSRGM	5.3204	2.3066	65.7919	2.0616	2.0734	0.5176	7.2604	0.9784	23.2980	2.3298	0.9618	0.9448	0.1543406	16
16	RMD	3.3303	1.8249	66.7588	1.6042	1.6053	0.5115	5.7443	4.2796	16.7058	1.6706	0.9761	0.9655	0.1289773	6
17	CT	4.0547	2.0136	68.0606	1.6777	1.6784	0.4961	6.0130	2.9906	18.1577	2.0175	0.9738	0.9574	0.1337059	7
18	Vtub	3.7733	1.9425	66.2648	1.6196	1.6206	0.4394	5.8006	2.1000	18.4484	2.0498	0.9756	0.9604	0.1271222	4
19	3PFD	4.3488	2.0854	68.6419	1.7359	1.7361	0.5136	6.2272	3.0687	18.8258	2.0918	0.9719	0.9543	0.1387252	12
20	New Model	2.4340	1.5601	60.2184	1.2993	1.2997	0.0709	4.6588	0.0599	13.8037	1.5337	0.9843	0.9744	0.0940763	1

, the proposed model achieves the best results when comparing the NCD value to the other models. Looking at Table 3 in detail, the MSE, RMSE, AIC, variance, RMSPE, PRR, TS, PP, and SAE values for the proposed model are the lowest as compared with all models. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values for the proposed model are the largest as compared with all models. First, the MSE value of the proposed model is 2.4340, which is the lowest value compared to the other models; next, the MSE value of the YID 1 model is 2.9149. The RMSE value of the proposed model is 1.5601, which is the lowest value among the models compared. Furthermore, the RMSE value of the YID 1 model is 1.7073. The AIC value of the proposed model is 60.2184, which is also the lowest value among the models compared, and the AIC value of the GO model is 62.8309. The variance value of the proposed model is 1.2993, which is the lowest value among the models compared, and the variance value of the YID 1 model is 1.5841. The RMSPE value of the proposed model is 1.2997, which is the lowest value among models compared, and the RMSPE value of the YID 1 model is 1.5883. The PRR value of the proposed model is 0.0709, which is the lowest value among the models compared, and the PRR value of the Vtub model is 0.4394. The TS value of the proposed model is 4.6588, which is the lowest value among the models compared, and the TS value of the YID 1 model is 5.6364. The PP value of the proposed model is 0.0599, which is the lowest value among the models compared, and the PP value of the DS model is 0.3575. The SAE value of the proposed model is 13.8037, which is the lowest value among models compared, and the SAE value of the PNZ model is 16.1109. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values of the proposed model are 0.9843 and 0.9744, respectively, which are the largest among the models compared. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values of the YID 1 model are 0.9970 and 0.9701, respectively. Although both the YID 1 and YID 2 models are slightly smaller than our proposed model based on the MAE criterion (i.e., the MAE value of YID 1 and YID 2 models are 1.4727 and 1.4741 as compared with our proposed model which is 1.5337), the performance of our proposed model, based on all the other remaining eleven criteria, is not only much better as compared with the YID 1 and YID 2 models but also as compared with all the models as shown in Table 3. Figure 1 shows a graph of the mean value functions for all models based on dataset #1. Figure 2 shows a graph of the relative error value of all models for dataset #1.

From Table 4, looking at the results for the NCD values, the proposed model achieves the lowest values among all models. The NCD value of the proposed model is 0.0940763, which is the lowest value among all the models, and the NCD values of the YID 1 model and the YID 2 model are 0.1245350 and 0.1254775, respectively. Figure 3a shows a 3D plot of the MSE value, model number and NCD value, and Figure 3b shows a 3D plot of the R² value, model number and NCD value for dataset #1.

Table 5. Results of the 95% and 99% confidence intervals for the proposed model for dataset #1.

Time Index	Data Value	95%		99%
		LCL	UCL	LCL	UCL
1	2	0	5.9683	0	5.0607
2	13	2.3023	19.1903	4.3213	17.1714
3	15	5.0574	25.0431	7.4466	22.6539
4	19	7.3770	29.4975	10.0215	26.8531
5	22	9.4082	33.1811	12.2502	30.3391
6	23	11.2268	36.3542	14.2307	33.3503
7	24	12.8800	39.1580	16.0215	36.0165
8	26	14.4001	41.6794	17.6612	38.4182
9	30	15.8101	43.9767	19.1773	40.6094
10	30	17.1273	46.0910	20.5898	42.6285
11	34	18.3651	48.0527	21.9142	44.5036
12	35	19.5340	49.8848	23.1623	46.2565
13	38	20.6424	51.6055	24.3439	47.9040
14	38	21.6972	53.2291	25.4668	49.4595

lists the 95% and 99% confidence intervals (lower confidence limit (LCL) and upper confidence limit (UCL)) for the proposed model for dataset #1, respectively. Figure 4 shows a graph of the 95% and 99% confidence intervals for the proposed model. The relative value and the confidence interval graphs aim to confirm the ability to improve accuracy and whether the number returned by the mean value function is included in the confidence interval of each time point.

4.2. Dataset 2

Dataset #2 was provided by [34] and is shown in

Table 6. Failure data–Dataset #2.

Time Index	Cumulative System Days	Dataset #2
		Failures	Cumulative Failures
1	1249	4	4
2	4721	6	10
3	8786	4	14
4	13,669	3	17
5	19,094	6	23
6	24,750	1	24
7	32,299	2	26
8	40,594	4	30
9	49,476	1	31
10	55,596	0	31
11	58,061	1	32
12	58,588	1	33
13	58,633	0	33

. In dataset #2, the weekly index uses cumulative system days, and there are 33 cumulative failures in 58,633 system days. The detailed information is recorded in [34].

In

Table 7. Results of model parameter estimation and criteria from dataset #2.

No.	Model	Parameter Estimation	MSE	RMSE	AIC	Variation	RMSPE	PRR	TS	PP	SAE	MAE	R2	Adj R2
1	GO	$\hat{\mathrm{a}}$ = 33.0367 $\hat{\mathrm{b}}$ = 5.8 × 10−5	1.6266	1.2754	49.4070	1.2770	1.2937	0.6251	4.6136	0.2420	12.9085	1.1735	0.9839	0.9806
2	DS	$\hat{\mathrm{a}}$ = 31.0119 $\hat{\mathrm{b}}$ = 1.44 × 10−4	7.3713	2.7150	72.2378	2.8666	2.9440	65.1831	9.8214	1.1664	25.6479	2.3316	0.9269	0.9122
3	IS	$\hat{\mathrm{a}}$ = 33.0368, $\hat{\mathrm{b}}$ = 5.8×10−5 $\hat{\mathsf{\beta }}$ = 8.6 × 10−8	1.7892	1.3376	51.4070	1.2770	1.2937	0.6251	4.6136	0.2420	12.9084	1.2908	0.9839	0.9785
4	YE	$\hat{\mathrm{a}}$ = 4438.167, $\hat{\mathsf{\alpha }}$ = 0.2021 $\hat{\mathsf{\beta }}$ = 5.76 × 10−5, $\hat{\mathsf{\gamma }}$ = 0.03696	1.9812	1.4076	53.3655	1.2984	1.3219	0.6363	4.6057	0.2445	12.9425	1.4381	0.9839	0.9759
5	YR	$\hat{\mathrm{a}}$ = 32.7479, $\hat{\mathsf{\alpha }}$ = 2.3803 $\hat{\mathsf{\beta }}$ = 3.87 × 10−9, $\hat{\mathsf{\gamma }}$ = 1.1396	12.8131	3.5795	91.6713	3.4310	3.5266	199.5366	11.7126	1.4270	31.6731	3.5192	0.8960	0.8440
6	YID 1	$\hat{\mathrm{a}}$ = 23.8781, $\hat{\mathrm{b}}$ = 9.83 × 10−5 $\hat{\mathsf{\alpha }}$ = 6.45 × 10−6	1.0219	1.0109	48.4929	0.9407	0.9461	0.2239	3.4867	0.1203	10.0450	1.0045	0.9908	0.9877
7	YID 2	$\hat{\mathrm{a}}$ = 22.6444, $\hat{\mathrm{b}}$ = 0.000106 $\hat{\mathsf{\alpha }}$ = 9.02 × 10−6	0.9979	0.9990	48.4804	0.9221	0.9252	0.1979	3.4455	0.1107	9.8882	0.9888	0.9910	0.9880
8	HDGO	$\hat{\mathrm{a}}$ = 33.0367, $\hat{\mathrm{b}}$ = 5.78 × 10−5 $\hat{\mathrm{c}}$ = 8.6 × 10−8	1.7883	1.3373	51.3850	1.2904	1.3110	0.6357	4.6124	0.2444	12.9353	1.2935	0.9839	0.9785
9	PNZ	$\hat{\mathrm{a}}$ = 22.6389, $\hat{\mathrm{b}}$ = 1.06 × 10−4 $\hat{\mathsf{\alpha }}$ = 0.9 × 10−5, $\hat{\mathsf{\beta }}$ = 0.00096	1.1090	1.0531	50.4830	0.9279	0.9327	0.1993	3.4458	0.1111	9.9067	1.1007	0.9910	0.9865
10	PZ	$\hat{\mathrm{a}}$ = 29.0937, $\hat{\mathrm{b}}$ = 0.00063 $\hat{\mathsf{\alpha }}$ = 4.98 × 10−5, $\hat{\mathsf{\beta }}$ = 0.0006 $\hat{\mathrm{c}}$ = 4.7088	1.3039	1.1419	53.3626	0.9950	1.0134	0.1180	3.5226	0.0801	10.2318	1.2790	0.9906	0.9839
11	ZFR	$\hat{\mathrm{a}}$ = 33.9535, $\hat{\mathrm{b}}$ = 0.01158 $\hat{\mathsf{\alpha }}$ = 0.01502, $\hat{\mathsf{\beta }}$ = 0.01344 $\hat{\mathrm{c}}$ = 5.6 × 10−5, $\hat{\mathrm{p}}$ = 1.0412	2.5591	1.5997	57.3622	1.3106	1.3368	0.6517	4.6163	0.2481	12.9732	1.8533	0.9838	0.9677
12	TP	$\hat{\mathrm{a}}$ = 228.3963, $\hat{\mathrm{b}}$ = 0.000117 $\hat{\mathsf{\alpha }}$ = 0.00517, $\hat{\mathsf{\beta }}$ = 0.03607 $\hat{\mathrm{c}}$ = 0.1539, $\hat{\mathrm{p}}$ = 0.1941 $\hat{\mathrm{q}}$ = 0.0837	1.3344	1.1551	56.0182	0.8245	0.8269	0.0332	3.0861	0.0292	8.2787	1.3798	0.9928	0.9827
13	IFD	$\hat{\mathrm{a}}$ = 3.283, $\hat{\mathrm{b}}$ = 0.000174 $\hat{\mathrm{d}}$ = 2 × 10−15	38.8383	6.2320	65.7427	7.2581	7.6766	23.6805	21.4949	1.8984	64.6863	6.4686	0.6497	0.5330
14	PDP	$\hat{\mathsf{\alpha }}$ = 278.0266, $\hat{\mathsf{\gamma }}$ = 5.8 × 10−6 ${\hat{\mathrm{t}}}_0$ = 4.9046, ${\hat{\mathrm{m}}}_0$ = 15.3182	34.4540	5.8698	134.4592	5.1188	5.1297	1.0557	19.2064	8.5495	49.5436	5.5048	0.7203	0.5805
15	KSRGM	$\hat{\mathrm{A}}$ = 31.8755, $\hat{\mathrm{b}}$ = 0.012 $\hat{\mathsf{\alpha }}$ = 0.0179, $\hat{\mathrm{p}}$ = 0.005	3.2468	1.8019	55.0954	2.2388	2.4094	2.5919	5.8960	0.4886	16.0865	1.7874	0.9736	0.9605
16	RMD	$\hat{\mathrm{a}}$ = 168.8421, $\hat{\mathrm{b}}$ = 5.8 × 10−5 $\hat{\mathsf{\alpha }}$ = 0.1956, $\hat{\mathsf{\beta }}$ = 0.1229	2.0411	1.4287	53.5239	1.3097	1.3312	0.7194	4.6747	0.2609	13.0707	1.4523	0.9834	0.9751
17	CT	$\hat{\mathrm{a}}$ = 0.0068, $\hat{\mathrm{b}}$ = 0.8298 $\hat{\mathsf{\alpha }}$ = 0.7555, $\hat{\mathsf{\beta }}$ = 58.3462 $\hat{\mathrm{N}}$ = 53.3901	0.9229	0.9607	51.9769	0.7878	0.7889	0.0191	2.9636	0.0184	7.8468	0.9809	0.9933	0.9886
18	Vtub	$\hat{\mathrm{a}}$ = 1.0006, $\hat{\mathrm{b}}$ = 0.5024 $\hat{\mathsf{\alpha }}$ = 9.3307, $\hat{\mathsf{\beta }}$ = 2.9469 $\hat{\mathrm{N}}$ = 85.1389	2.0950	1.4474	52.9672	1.2106	1.2193	0.1150	4.4652	0.1989	11.6921	1.4615	0.9849	0.9741
19	3PFD	$\hat{\mathrm{a}}$ = 0.0059, $\hat{\mathrm{b}}$ = 0.000001 $\hat{\mathsf{\beta }}$ = 2.870, $\hat{\mathrm{N}}$ = 41.540 $\hat{\mathrm{c}}$ = 34.600	1.2405	1.1138	52.9005	0.9290	0.9349	0.21+34	3.4360	0.1165	9.8563	1.2320	0.9910	0.9847
20	New Model	$\hat{\mathrm{b}}$ = 0.000108, $\hat{\mathsf{\alpha }}$ = 0.6413 $\hat{\mathsf{\beta }}$ = 2.5635, $\hat{\mathrm{N}}$ = 42.4984 $\hat{\mathrm{d}}$ = 5.4 × 10−5	0.8790	0.9376	52.0859	0.7662	0.7664	0.0149	2.8923	0.0149	7.6843	0.9605	0.9937	0.9891

, we obtained the parameter estimates of all twenty models and the values of all twelve criteria using $\mathrm{t}=1249,4721,8786,\cdots ,58,633$ from the cumulative system days from dataset #2. As shown in Table 7, we can see that the proposed model achieves the best results when comparing the twelve criteria to the other models. As shown in

Table 8. Results for criteria and NCD value for comparison with dataset #2.

No.	Model	MSE	RMSE	AIC	Variation	RMSPE	PRR	TS	PP	SAE	MAE	R2	Adj R2	NCD Value	Rank
1	GO	1.6266	1.2754	49.4070	1.2770	1.2937	0.6251	4.6136	0.2420	12.9085	1.1735	0.9839	0.9806	0.0980582	9
2	DS	7.3713	2.7150	72.2378	2.8666	2.9440	65.1831	9.8214	1.1664	25.6479	2.3316	0.9269	0.9122	0.3193186	17
3	IS	1.7892	1.3376	51.4070	1.2770	1.2937	0.6251	4.6136	0.2420	12.9084	1.2908	0.9839	0.9785	0.1007023	10
4	YE	1.9812	1.4076	53.3655	1.2984	1.3219	0.6363	4.6057	0.2445	12.9425	1.4381	0.9839	0.9759	0.1043514	13
5	YR	12.8131	3.5795	91.6713	3.4310	3.5266	199.5366	11.7126	1.4270	31.6731	3.5192	0.8960	0.8440	0.7414780	18
6	YID 1	1.0219	1.0109	48.4929	0.9407	0.9461	0.2239	3.4867	0.1203	10.0450	1.0045	0.9908	0.9877	0.0783218	4
7	YID 2	0.9979	0.9990	48.4804	0.9221	0.9252	0.1979	3.4455	0.1107	9.8882	0.9888	0.9910	0.9880	0.0773360	3
8	HDGO	1.7883	1.3373	51.3850	1.2904	1.3110	0.6357	4.6124	0.2444	12.9353	1.2935	0.9839	0.9785	0.1010407	11
9	PNZ	1.1090	1.0531	50.4830	0.9279	0.9327	0.1993	3.4458	0.1111	9.9067	1.1007	0.9910	0.9865	0.0800413	5
10	PZ	1.3039	1.1419	53.3626	0.9950	1.0134	0.1180	3.5226	0.0801	10.2318	1.2790	0.9906	0.9839	0.0858796	8
11	ZFR	2.5591	1.5997	57.3622	1.3106	1.3368	0.6517	4.6163	0.2481	12.9732	1.8533	0.9838	0.9677	0.1139585	15
12	TP	1.3344	1.1551	56.0182	0.8245	0.8269	0.0332	3.0861	0.0292	8.2787	1.3798	0.9928	0.9827	0.0824943	6
13	IFD	38.8383	6.2320	65.7427	7.2581	7.6766	23.6805	21.4949	1.8984	64.6863	6.4686	0.6497	0.5330	0.7418298	19
14	PDP	34.4540	5.8698	134.4592	5.1188	5.1297	1.0557	19.2064	8.5495	49.5436	5.5048	0.7203	0.5805	0.8201567	20
15	KSRGM	3.2468	1.8019	55.0954	2.2388	2.4094	2.5919	5.8960	0.4886	16.0865	1.7874	0.9736	0.9605	0.1470449	16
16	RMD	2.0411	1.4287	53.5239	1.3097	1.3312	0.7194	4.6747	0.2609	13.0707	1.4523	0.9834	0.9751	0.1056358	14
17	CT	0.9229	0.9607	51.9769	0.7878	0.7889	0.0191	2.9636	0.0184	7.8468	0.9809	0.9933	0.9886	0.0724187	2
18	Vtub	2.0950	1.4474	52.9672	1.2106	1.2193	0.1150	4.4652	0.1989	11.6921	1.4615	0.9849	0.9741	0.1013477	12
19	3PFD	1.2405	1.1138	52.9005	0.9290	0.9349	0.2134	3.4360	0.1165	9.8563	1.2320	0.9910	0.9847	0.0831768	7
20	New Model	0.8790	0.9376	52.0859	0.7662	0.7664	0.0149	2.8923	0.0149	7.6843	0.9605	0.9937	0.9891	0.0713073	1

, the proposed model also achieves the best results when comparing NCD values to the other models.

Looking at Table 7 in detail, we can observe that the MSE, RMSE, AIC, variance, RMSPE, PRR, TS, PP, and SAE values for the proposed model are the lowest values among all models compared. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values for the proposed model are the largest values among all the models compared. First, the MSE value of the proposed model is 0.8790, which is the lowest value among all the models compared, and the MSE value of the YID 2 model is 0.9229. The RMSE value of the proposed model is 0.9376, which is the lowest among all the models compared, and the RMSE value of the YID 2 model is 0.9607. The variance value of the proposed model is 0.7662, which is the lowest among all the models compared, and the variance value of the CT model is 0.7878. The RMSPE value of the proposed model is 0.7664, which is the lowest among all the models compared, and the RMSPE value of the CT model is 0.7889. The PRR value of the proposed model is 0.0149, which is the lowest among all the models compared, and the PRR value of the CT model is 0.0191. The TS value of the proposed model is 2.8923, which is the lowest among all the models compared, and the TS value of the CT model is 2.9636. The PP value of the proposed model is 0.0149, which is the lowest among all the models compared, and the PP value of the CT model is 0.0184. The SAE value of the proposed model is 7.6843, which is the lowest among all the models compared, and the SAE value of the CT model is 7.8468. The MAE value of the proposed model is 0.9605, which is the lowest among all the models compared, and the MAE value of the CT model is 0.9809. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values of the proposed model are 0.0.9937 and 0.9891, respectively, which are the largest among all the models compared. The ${\mathrm{R}}^2$ and ${\mathrm{Adj}\mathrm{R}}^2$ values of the YID 1 model are 0.9933 and 0.9886, respectively. Although the YID 2 model is smaller than our proposed model based on AIC criterion, the performance of our proposed model based on all the other remaining eleven criteria are not only much better as compared with the YID 2 model but also as compared with all the models, as shown in Table 7. Figure 5 shows a graph of the mean value functions for all the models based on dataset #2. Figure 6 shows a graph of the relative error value of all models for dataset #2.

Moreover, from Table 8, looking at the results for the NCD values, the proposed model achieves the lowest values among all the models compared. The NCD value of the proposed model is 0.0713073, which is the lowest among all the models compared, and the NCD values of the CT model and the YID 2 model are 0.0724187 and 0.0773360, respectively. Figure 7a shows a 3D plot of the MSE value, model number and NCD value, and Figure 7b shows a 3D plot of the R² value, model number and NCD value for dataset #2.

Table 9. Results of the 95% and 99% confidence intervals for the proposed model for dataset #2.

Time Index	Data Value	95%		99%
		LCL	UCL	LCL	UCL
1249	4	0	9.3745	0.1500	8.1220
4721	10	1.6579	17.6748	3.5727	15.7600
8786	14	4.4427	23.8028	6.7572	21.4884
13,669	17	7.1600	29.0935	9.7821	26.4714
19,094	23	9.5560	33.4429	12.4116	30.5873
24,750	24	11.5339	36.8805	14.5640	33.8504
32,299	26	13.5802	40.3256	16.7775	37.1283
40,594	30	15.2819	43.1205	18.6099	39.7925
49,476	31	16.6713	45.3623	20.1012	41.9324
55,596	31	17.4423	46.5924	20.9271	43.1076
58,061	32	17.7190	47.0315	21.2232	43.5273
58,588	33	17.7759	47.1217	21.2841	43.6135
58,633	33	17.7807	47.1293	21.2892	43.6208

lists the 95% and 99% confidence intervals of the proposed model for dataset #2. Figure 8 shows a graph of the 95% and 99% confidence intervals of the proposed model for dataset #2.

5. Sensitivity Analysis

In this study, we investigated the effect of each parameter of the proposed model on the mean value function by performing sensitivity analysis. We conducted the sensitivity analysis by changing one of the parameters of the model and fixing all the other parameters. We examined how the value of the estimated mean function changes when the estimated parameter values obtained from dataset #1 and dataset #2 change from –20% to +20%. Figure 9 shows sensitivity analysis for five parameters in the proposed model based on dataset #1. As can be seen in Figure 9, the cumulative number of detected failures will change with each estimated parameter. Figure 10 shows the sensitivity analysis for five parameters in the proposed model based on dataset #2. From Figure 10, it can be seen that the cumulative number of detected failures will be change with each estimated parameter. Therefore, it appears that the parameters are all influential in the proposed model.

6. Conclusions

Software is used in a variety of areas and environments, driven by the needs of different consumers. That is why we must develop high-quality software and improve its reliability as well. However, in general, software is developed in a controlled testing environment. Therefore, we considered the uncertainty or randomness of the software operating environment, and testing coverage where quantitative confidence criteria for software products can be applied. In this paper, we discussed a new testing coverage model based on NHPP software reliability with the uncertainty of operating environments and provided the sensitivity analysis regarding the impact of each parameter of the proposed model. As a result, we provided twelve criteria and the NCD value to compare the goodness-of-fit for the proposed model with several existing NHPP software reliability models. As can be seen in the results of model analysis, the results show that the proposed model achieves significantly better goodness-of-fit than other models. In addition, we investigated the impact of each parameter of the proposed model on the mean value function by performing sensitivity analysis.