One-Sample Bayesian Predictive Analyses for a Nonhomogeneous Poisson Process with Delayed S-Shaped Intensity Function Using Non-Informative Priors

Otieno Collins, Orawo Luke Akong’o, Matiri George Munene

Department of Mathematics, Egerton University, Egerton, Kenya.

DOI: 10.4236/ojs.2023.135034                

Abstract

The delayed S-shaped software reliability growth model (SRGM) is one of the non-homogeneous Poisson process (NHPP) models which have been proposed for software reliability assessment. The model is distinctive because it has a mean value function that reflects the delay in failure reporting: there is a delay between failure detection and reporting time. The model captures error detection, isolation, and removal processes, thus is appropriate for software reliability analysis. Predictive analysis in software testing is useful in modifying, debugging, and determining when to terminate software development testing processes. However, Bayesian predictive analyses on the delayed S-shaped model have not been extensively explored. This paper uses the delayed S-shaped SRGM to address four issues in one-sample prediction associated with the software development testing process. Bayesian approach based on non-informative priors was used to derive explicit solutions for the four issues, and the developed methodologies were illustrated using real data.

Keywords

Failure Intensity, Non-Informative Priors, Software Reliability Model, Bayesian Approach, Non-Homogeneous Poisson Process

1. Introduction

Software reliability assessment has become indispensable for all software developers as it involves reliability testing and debugging processes. The main objective of reliability testing is to identify potential defects that could lead to system failures, crashes, or malfunctions during real-world usage. Debugging is the process of detecting and correcting errors in software [1] . Software end-users are often concerned with the reliability of software products they acquire from the market. A single software error can cause a failure, which can be avoided by producing reliable software [2] . Defective software may not only damage the reputation of the producer but also attract legal procedures in case of a lawsuit and may expose the developer to preventable costs. As such, software developers are concerned with the reliability of their software products before releasing them into the market. Many researchers have developed methods to test software to detect, isolate, and remove faults during development. The reliability of software products is ensured by running tests that emulate the end-user environment.

Although the approaches to testing software reliability have been proposed, a number of issues still arise with testing. Such may include determining the optimal release time of software, the optimal cost of producing software, and the reliability of the software. Moreover, it is difficult and time-consuming to emulate the end-user environment [2] . These issues are well addressed through software reliability modeling by performing predictive analyses using historical data of software failures. The analyses often involve determining the optimal release time of software and expected costs at the time of release and constructing a prediction interval for a future observable failure.

Over the past decades, several nonhomogeneous Poisson process (NHPP) software reliability growth models (SRGMs) have been developed and used in reliability assessment. Such models are classified as perfect or imperfect debugging NHPP SRGMs [1] . Perfect debugging SRGMs assume that errors are immediately removed with certainty when detected, without introducing new faults. On the other hand, imperfect debugging NHPP SRGMs assume that when faults are detected, they may not be removed with certainty, and new errors may also be introduced into the software during the correction process [1] . The delayed S-shaped software reliability model developed by [3] is one of the perfect NHPP SRGMs. The model is distinctive because of its mean value function, which is S-shaped, reflecting the delay in failure reporting [1] . The mean value function is given by;

$m\left(t\right)=\alpha \left(1-\left(1+\beta t\right){\text{e}}^{\beta t}\right)$ (1)

and the model has an intensity function given by;

$\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ (2)

The delayed S-shaped NHPP SRGM is based on the following assumptions [4] :

1) Errors are removed immediately when detected, without introducing new errors into the software.

2) The current number of faults in a software and the probability of failure detection are proportional.

3) All the faults in a software are mutually independent from the failure detection point of view.

4) Errors in a software leads to failures at random times.

5) The time between ${\left(j-1\right)}^{\text{th}}$ and $j^{\text{th}}$ failures depend on the time to the ${\left(j-1\right)}^{\text{th}}$ failure.

6) The initial error content of the software is a random variable.

Many researchers have considered the use of the delayed S-shaped software reliability model in software reliability testing. [5] performed Bayesian predictive analysis on the model using gamma-distributed informative prior, determining the optimal release time, expected costs, and the estimated reliability of the software at the time of release. [6] performed Bayesian interval estimation and compared Bayesian credible intervals with Wald confidence intervals and found that Bayesian approach yielded shorter intervals with higher coverage probability. The results imply that Bayesian method yields better parameter estimates, which can then be used to enhance accuracy in prediction. Although the estimation of the parameters of the model is essential in software reliability assessment, the main goal is to obtain accurate predictions to provide adequate information to help software developers in planning tests and inform them about when to terminate the testing process. [7] performed interval estimation on the delayed S-shaped model and obtained 90% prediction intervals for the reliability prediction at any

future time, t. However, none of these studies used the $\frac{1}{\alpha \beta }$ non-informative

prior, and predictive issues addressed in this paper have not been explored using the delayed S-shaped model.

This paper focuses on single-sample predictive analyses on the delayed S-shaped software reliability model using Bayesian approach. We first outline four issues in software reliability testing. The issues have been addressed by [8] and [9] using the Power law Process (PLP), [2] using the Goel-Okumoto (1979) software reliability model, and [10] using Musa-Okumoto SRGM. The study used Bayesian approach with non-informative priors and the Yamada delayed S-shaped SRGM to develop and derive predictive distributions presented in section 2.2, and applied the developed methodologies to secondary software data to address the four issues, as discussed in section 4.

2. Predictive Issues and Bayesian Method

In this analysis, we assume that a reliability growth testing is performed on a software and the cumulative number of failures, denoted by $N\left(t\right)$ , is observed in the time interval $\left(0,t\right]$ . Another assumption is that the cumulative number of failures and failure times ( $0<t_1<t_2<\cdots$ ) follow the NHPP with the intensity function given in Equation (2). When testing stops after a predetermined number of failures, n, the failure data is said to be failure-truncated, and the n failure times are denoted by $Y_{obs}^f={\left[t_i\right]}_{i=1}^n$ . However, if testing stops at a predetermined time, t, the failure data is said to be time-truncated, and the corresponding observed failure data is denoted by $Y_{obs}^t=\left[n,t_1,\cdots ,t_n;t\right]$ .

Let $\underset{\_}{t}$ be the vector of observed failure times and $\pi \left(\theta \right)$ be the prior density of $\theta =\left(\alpha ,\beta \right)$ . Then according to Bayes’ rule, the posterior density $\pi \left(\theta /\underset{\_}{t}\right)$ is obtained using the formula:

$\pi \left(\theta |\underset{\_}{t}\right)=\frac{g\left(\theta ,\underset{\_}{t}\right)}{f\left(\underset{\_}{t}\right)}=\frac{f\left(\underset{\_}{t}|\theta \right)\pi \left(\theta \right)}{f\left(\underset{\_}{t}\right)}$ (3)

where $g\left(\theta ,\underset{\_}{t}\right)$ is the joint density of $\theta$ and $\underset{\_}{T}$ , and $f\left(\underset{\_}{t}\right)$ is the marginal density of $\underset{\_}{T}$ . The posterior predictive distribution of t⁺ is given as;

$\begin{array}{c}f\left(t^{+}|\underset{\_}{t}\right)={\displaystyle \int f}\left(t^{+},\theta |\underset{\_}{t}\right)\text{d}\theta \\ ={\displaystyle \int f}\left(t^{+}|\theta ,\underset{\_}{t}\right)\pi \left(\theta |\underset{\_}{t}\right)\text{d}\theta \\ ={\displaystyle \int f}\left(t^{+}|\theta \right)\pi \left(\theta |\underset{\_}{t}\right)\text{d}\theta \end{array}$ (4)

2.1. Issues in Single-Sample Software Reliability Prediction

The study addressed the following four issues in single-sample software reliability testing:

1) Suppose that the predetermined target value, ${\lambda }_{tv}$ , for the software failure rate is not achieved at time T, what is the probability that the target value will be achieved at time $\tau$ , $\tau >T$ ?

2) Suppose that the target value, ${\lambda }_{tv}$ , for the software failure rate is not achieved at time T, how long will it take so that the software failure rate will be attained at ${\lambda }_{tv}$ ?

3) What is the probability that at most k software failures will occur in the failure time period $\left(T,\tau \right]$ , $\tau >T$ ?

4) What is the upper prediction limit (UPL) of $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ with level $\lambda$ , $\tau$ being a predetermined value greater than T?

2.2. Prior, Posterior, and Predictive Distributions

Let $Y_{obs}$ represent $Y_{obs}^t$ or $Y_{obs}^f$ . The joint density of $Y_{obs}$ is obtained as:

$f\left(Y_{obs}|\alpha ,\beta \right)={\alpha }^n{\beta }^{2n}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}$ (5)

and the log-likelihood function is given by:

$\begin{array}{c}l\left(\alpha ,\beta |\underset{\_}{t}\right)=n\log \alpha +2n\log \beta +\log \left({\displaystyle {\prod }_{i=1}^n{t}_i}\right)\\ -\beta {\displaystyle \sum t_i}-\alpha \left(1-\left(1+\beta T\right)\exp \left(-\beta T\right)\right)\end{array}$ (6)

Case 1: When the shape parameter, β, is known, we adopt the following non-informative prior distribution for α;

$\pi \left(\alpha \right)\propto \frac{1}{\alpha },\alpha >0$ (7)

The posterior distribution of α can be obtained from Equation (3) as:

$\pi \left(\alpha |Y_{obs}\right)=\frac{f\left(Y_{obs}\alpha ,\beta \right)\pi \left(\alpha ,\beta \right)}{{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }f}}\left(Y_{obs}\alpha ,\beta \right)\pi \left(\alpha ,\beta \right)\text{d}\alpha \text{d}\beta }$

$\pi \left(\alpha |Y_{obs}\right)={\left[\Gamma \left(n\right)\right]}^{-1}{\alpha }^{n-1}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n$ (8)

Let t⁺ be the random variable being predicted. The posterior predictive density of t⁺ is obtained using Equation (4) as:

$f\left(t^{+}|Y_{obs}\right)={\displaystyle {\int }_0^{\infty }f}\left(t^{+}|Y_{obs},\alpha \right)\pi \left(\alpha |Y_{obs}\right)\text{d}\alpha$ (9)

The Bayesian UPL for t⁺ with level $\gamma$ satisfies:

$\gamma ={\displaystyle {\int }_{-\infty }^{y_U^{\left(\beta \right)}}p\left(t^{+}|Y_{obs}\right)\text{d}{t}^{+}}$ (10)

Case 2: When the shape parameter, β, is unknown, we consider the following non-informative density of α and β, assuming they are mutually independent.

$\pi \left(\alpha ,\beta \right)\propto \frac{1}{\alpha \beta },\alpha ,\beta >0$ (11)

The corresponding posterior joint density is obtained using Equation (3) as:

$\pi \left(\alpha ,\beta |Y_{obs}\right)={\left[k\Gamma \left(n\right)\right]}^{-1}{\alpha }^{n-1}{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}$ (12)

where;

$k={\displaystyle {\int }_0^{\infty }\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left(1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right)}^n}\text{d}\beta }$ (13)

Similar to Equation (9) and Equation (10), if t⁺ is the random variable being predicted, the posterior predictive distribution becomes:

$f\left(t^{+}|Y_{obs}\right)={\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }f}}\left(t^{+}|Y_{obs},\alpha ,\beta \right)\pi \left(\alpha ,\beta |Y_{obs}\right)\text{d}\alpha \text{d}\beta$ (14)

and the Bayesian UPL is:

$\gamma ={\displaystyle {\int }_{-\infty }^{yU}p}\left(t^{+}|Y_{obs}\right)\text{d}{t}^{+}$ (15)

3. Main Results for Prediction

In this section, we present the main results of the four issues stated in section 2.1 as propositions, and their proofs are given in the Appendix.

Proposition 1:

The probability that the target value ${\lambda }_{tv}$ will be achieved at time $\tau$ ( $\tau >T$ ) is:

$\gamma =\{\begin{array}{ll}1-\underset{h=0}{\overset{n-1}{{\displaystyle \sum }}}\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}},\hfill & \text{if}\beta \text{isknown}\hfill \\ 1-\frac{1}{k}\underset{h=0}{\overset{n-1}{{\displaystyle \sum }}}\underset{0}{\overset{\infty }{{\displaystyle \int }}}\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta ,\hfill & \text{if}\beta \text{isunknown}\hfill \end{array}$

Proposition II:

Let ${\tau }^{\ast }$ denote the time required to attain ${\lambda }_{tv}$ . For a specified level $\gamma$ ;

${\tau }^{\ast }=\{\begin{array}{ll}\left[-\frac{1}{\beta }{W}_n\left(\frac{-2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\beta {\chi }^2\left(2n;\gamma \right)}\right)\right]-T\hfill & \text{if}\beta \text{isknown}\hfill \\ \tau -T\hfill & \text{if}\beta \text{isunknown}\hfill \end{array}$

Remark 1:

$\gamma =1-\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}$ (16)

Proposition III:

The probability that at most k failures will occur in the time interval $\left(T,\tau \right]$ , $\tau >T$ is given by:

${\gamma }_k=\{\begin{array}{ll}\frac{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^n}{\displaystyle {\sum }_{j=n}^{n+k}\left(\begin{array}{c}j-n\\ n-1\end{array}\right)}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^j}{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}\hfill & \text{if}\beta \text{isknown}\hfill \\ {\displaystyle {\sum }_{j=n}^{n+k}\frac{\Gamma \left(j\right)}{c\left(j-n\right)!\Gamma \left(n\right)}}{\displaystyle {\int }_0^{\infty }{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}\text{d}\beta \hfill & \text{if}\beta \text{isunknown}\hfill \end{array}$

Proposition IV:

The Bayesian UPL of $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ with level $\gamma$ is obtained as:

${\lambda }_U^{\left(\beta \right)}=\{\begin{array}{ll}\frac{\left({\beta }^{2}t{\text{e}}^{-\beta t}\right){\chi }^2\left(2n;\gamma \right)}{2\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}\hfill & \text{if}\beta \text{isknown}\hfill \\ {\lambda }_{tv}\hfill & \text{if}\beta \text{isunknown}\hfill \end{array}$

Remark 2:

$\gamma =1-\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}$ (17)

4. Real Data Application

In this section, real data in Table 1 was used for single-sample Bayesian prediction, and the results are presented and discussed. The study used failure times (cumulative time between failures) $0<t_1<t_2<\cdots <t_{22}$ , where n = 22. The data is given by [11] , while [12] argued that it has been widely used in assessing software reliability models. For the case where β is assumed known, the study performed maximum likelihood estimation (MLE), fixing initial guess of the parameters arbitrarily at $\alpha =18$ and $\beta =0.00342$ , where the values were chosen such that, they were closer to the ML estimates [11] obtained using their S-shaped model. However, parameter estimation methods can be used to get the initial set of these parameters. The initial parameter values were used with the log-likelihood function of the delayed S-shaped model given by Equation (6) to obtain ML estimates of α and β. The MLE for β was obtained as $\beta =0.007609807$ . The study used this value for the cases where β is assumed known.

Table 1. Time between failures data.

Proposition I: Suppose the target value is given by ${\lambda }_{tv}=0.02$ . At time $T=100$ , the MLE of the achieved software rate for this software is $\lambda =\alpha \beta {\text{e}}^{-100\beta }=0.0640$ , which is greater than 0.02. In this regard, the target value cannot be achieved at the initial time, $T=100$ . Assuming a time greater than 100, $\tau =500$ , and we want to determine the probability of achieving ${\lambda }_{tv}$ at this time: 1) when β is known ( $\beta =0.007609807$ ), from the first formula in Proposition I, we obtain $\gamma =9.0\times {10}^{-9}$ . Therefore, it is almost unlikely that the target value will be achieved. 2) When β is unknown, from the second formula in Proposition I, we obtain $\gamma =3.0\times {10}^{-7}$ . Thus, it is less likely that the target software failure intensity of 0.02 will be achieved at time $\tau =500\text{h}$ .

Since the target value was not achieved at $\tau =500$ , we want to assess the relationship between the probability and time by varying τ, while holding ${\lambda }_{tv}=0.02$ constant. This was illustrated in the case when β is unknown. Table 2 shows the probabilities and time, indicating that it is almost unlikely to achieve the target value at any time between 130 h and 630 h. However, it is almost certain that at 830 h and above, the target software failure intensity will be achieved. It is crucial to note that these values do not indicate the exact time at which the failure intensity is achieved but the time at which it shall have been achieved with some probability. Figure 1 displays the relationship, indicating that the probability of achieving a software failure intensity of 0.02 increases with time. It can be observed in the figure that the probability increases rapidly between 630 h and 780 h, suggesting that between these time periods and/or above, the target software failure intensity is achievable.

Since the software failure intensity of 0.02 ( ${\lambda }_{tv}=0.02$ ) was not achieved at $\tau =500\text{h}$ , we want to determine the failure intensity that is most likely at this time period. An illustration was performed for the case when β is unknown to assess the relationship between failure intensity and probability while holding the time, τ, at 500 h. Suppose the interest is to determine the failure intensity that is most achievable at $\tau =500\text{h}$ . The results displayed in Figure 2 were obtained. It can be observed that it is unlikely that the software will have a failure intensity of less than 0.05 at a time period of 500 h. It can also be observed that the most probable failure intensity is between 0.06 and 0.13. Since the delayed S-Shaped model assumes that errors are immediately removed after they are detected, and no further errors are introduced into the software [4] , a failure intensity that corresponds to a probability that tends to one (1) was selected as the highest, while the failure intensities at which the probability is exactly one were ignored.

Table 2. The probability that the software failure intensity of 0.02 is achieved at different time periods, τ.

Figure 1. The probability of achieving ${\lambda }_{tv}=0.02$ at different time periods.

Figure 2. The probability of achieving ${\lambda }_{tv}=\lambda$ at $\tau =500\text{h}$ .

Proposition II: The next step involved determining how long it will require to achieve the target value ${\lambda }_{tv}$ since it was not attained at time $T=100$ . 1) When β is known ( $\beta =0.0076098073$ ), and γ is 0.90, we obtain ${\tau }^{\ast }=583.365\text{h}$ using the first formula in Proposition II. The result implies that it will take another 583.365h to achieve the desired failure rate of 0.02 (i.e., the target failure intensity will be achieved at $\tau =683.365\text{h}$ ). 2) When β is unknown, we obtain ${\tau }^{\ast }=657\text{h}$ using the second formula in Proposition II and Remark 1, where ${\tau }^{\ast }$ is the value that satisfies Equation (16). Thus, it will require 657 additional hours to achieve the target failure rate of 0.02. The value implies that the target software failure intensity will be achieved at $\tau =757\text{h}$ .

Proposition III: Since the study has established that the probable software failure intensity at a future time interval $\left(100,\tau =500\right]$ , is high, we want to obtain the probability that at most k failures will occur at a future time interval $\left(T,{\tau }^{\ast }\right]$ , where $T<{\tau }^{\ast }<\tau$ . Suppose the interest is to determine the probability that at most k failures will occur in the future time interval $\left(T,{\tau }^{\ast }\right]=\left(100,130\right]$ , the following results were obtained for $k=25$ :

1). When β is known ( $\beta =0.007609807$ ), from the first formula in Proposition III, we obtain ${\gamma }_0=0.000213$ , ${\gamma }_1=0.00171$ , ${\gamma }_2=0.00720$ , ${\gamma }_3=0.02122$ , ${\gamma }_4=0.04916$ , ${\gamma }_5=0.09552$ , ${\gamma }_6=0.1621$ , ${\gamma }_7=0.2470$ , ${\gamma }_8=0.3452$ , ${\gamma }_9=0.4497$ , ${\gamma }_{10}=0.5529$ , ${\gamma }_{11}=0.6488$ , ${\gamma }_{12}=0.7329$ , ${\gamma }_{13}=0.8031$ , ${\gamma }_{14}=0.8590$ , ${\gamma }_{15}=0.9019$ , ${\gamma }_{16}=0.9335$ , ${\gamma }_{17}=0.9560$ , ${\gamma }_{18}=0.9716$ , ${\gamma }_{19}=0.9821$ , ${\gamma }_{20}=0.9889$ , ${\gamma }_{21}=0.9933$ , ${\gamma }_{22}=0.9960$ , ${\gamma }_{23}=0.9977$ , ${\gamma }_{24}=0.9987$ , and ${\gamma }_{25}=0.9992$ .

2). When β is unknown, using the second formula in Proposition III, we get ${\gamma }_0=0.000024$ , ${\gamma }_1=0.000230$ , ${\gamma }_2=0.00113$ , ${\gamma }_3=0.00391$ , ${\gamma }_4=0.01054$ , ${\gamma }_5=0.02375$ , ${\gamma }_6=0.04651$ , ${\gamma }_7=0.08137$ , ${\gamma }_8=0.1298$ , ${\gamma }_9=0.1916$ , ${\gamma }_{10}=0.2649$ , ${\gamma }_{11}=0.3466$ , ${\gamma }_{12}=0.4327$ , ${\gamma }_{13}=0.5189$ , ${\gamma }_{14}=0.6014$ , ${\gamma }_{15}=0.6772$ , ${\gamma }_{16}=0.7444$ , ${\gamma }_{17}=0.8018$ , ${\gamma }_{18}=0.8495$ , ${\gamma }_{19}=0.8880$ , ${\gamma }_{20}=0.9182$ , ${\gamma }_{21}=0.9413$ , ${\gamma }_{22}=0.9586$ , ${\gamma }_{23}=0.9713$ , ${\gamma }_{24}=0.9804$ , and ${\gamma }_{25}=0.9868$ .

Figure 3. The graph of the probabilities ${\gamma }_k$ that at most k failures (k = 25) will occur in the time interval $\left(100,130\right]$ for the cases of known and unknown β.

Proposition IV: The Bayesian Upper Prediction Limit (UPL) of $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ given $\tau =700\text{h}$ and the level, $\gamma =0.90$ : 1) when β is known ( $\beta =0.0076098073$ ), the Bayesian UPL was obtained using the first formula in (Proposition IV) as ${\lambda }_U^{\left(\beta \right)}\left(\tau \right)=0.01805$ ; 2) when β is unknown, using the second part of the formula in Proposition IV and Remark 2, the Bayesian UPL of $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ with level $\gamma =0.90$ was obtained as ${\lambda }_U^{\left(\beta \right)}\left(\tau \right)=0.02883$ (Figure 3).

5. Conclusion

Software reliability remains a priority for software developers. Defective software may damage the reputation of the developers, expose them to preventable costs, and cause significant damage to end-users. To address these issues, software reliability modeling comes in handy. Predictive analyses have been performed in software testing to provide information for making decisions regarding the optimal software production costs, release time, and whether and when the desired reliability would be achieved. This study used non-informative priors to derive explicit solutions for four predictive issues in software testing using the Yamada delayed S-shaped model, and applied the derived methodologies with secondary software failure data. Bayesian approach was used because it incorporates prior information about the parameters of the model and is appropriate even when historical data is insufficient. The obtained solutions can be adequately applied in software quality assessment.

Appendices: Proof of Propositions I-IV

We begin by stating the following identity without proof:

${\displaystyle {\int }_{D\left(m;a,b\right)}\text{d}F\left(t_1\right)\cdots \text{d}F\left(t_m\right)}={\left[F\left(b\right)-F\left(a\right)\right]}^m/m!$ (A1)

where m is any positive integer, a and b are two real numbers such that a < b, $F\left(t\right)$ is an increasing and differentiable function, and $D\left(m;a,b\right)=\left[\left(t_1,\cdots ,t_m\right):a<t_1<\cdots <t_m<b\right]$ .

Proof of Proposition I:

The probability is given by:

Let $\pi \left(Y_{obs}|{\lambda }_{\tau }\right)$ denote the posterior density of $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ . Thus, the probability that ${\lambda }_{tv}$ will be achieved at time $\tau$ ( $\tau >T$ ) is given by:

$\gamma =\mathrm{Pr}\left({\lambda }_{\tau }\le {\lambda }_{tv}|Y_{obs}\right)={\displaystyle {\int }_0^{{\lambda }_{tv}}p}\left(Y_{obs}|{\lambda }_{\tau }\right)\text{d}{\lambda }_{\tau }$ (A2)

When β is known, we have $\alpha =\frac{\lambda \left(t\right)}{{\beta }^{2}t{\text{e}}^{-\beta t}}$ , obtained from the intensity function, and $\frac{\text{d}\alpha }{\text{d}{\lambda }_{\tau }}=\frac{1}{{\beta }^{2}t{\text{e}}^{-\beta t}}$ .

The posterior distribution of ${\lambda }_{\tau }$ is $\pi \left({\lambda }_{\tau }|Y_{obs}\right)=\pi \left(\alpha |Y_{obs}\right)\left|\frac{\text{d}\alpha }{\text{d}{\lambda }_{\tau }}\right|$ . Thus, the posterior density becomes:

$\pi \left({\lambda }_{\tau }|Y_{obs}\right)=\frac{1}{\Gamma \left(n\right)}{\left(\frac{\lambda \left(t\right)}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)}^{n-1}{\text{e}}^{-\left(\frac{\lambda \left(t\right)}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}$

Which reduces to:

$\pi \left({\lambda }_{\tau }|Y_{obs}\right)=\frac{{\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}^n}{\Gamma \left(n\right)}{\lambda }_{\tau }^{n-1}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}$ (A3)

From Equation (A3), ${\lambda }_{\tau }$ follows a gamma distribution with parameters n and $\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}$ . However, there is a relationship between gamma and Poisson distribution defined by:

$\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^{\lambda }{x}^{\alpha -1}{\text{e}}^{-\beta x}\text{d}x}=1-{\displaystyle {\sum }_{h=0}^{\alpha -1}\frac{{\left(\beta \lambda \right)}^h}{h!}{\text{e}}^{-\beta \lambda }}$ (A4)

From Equation (A2), ( $\gamma ={\displaystyle {\int }_0^{{\lambda }_{tv}}p}\left(Y_{obs}|{\lambda }_{\tau }\right)\text{d}{\lambda }_{\tau }$ ), and Equation (A4), it follows that:

$\gamma =1-{\displaystyle {\sum }_{h=0}^{n-1}\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}}$ (A5)

Equation (A5) shows the first formula in proposition I.

When β is unknown: from the intensity function, $\alpha =\frac{\lambda \left(t\right)}{{\beta }^{2}t{\text{e}}^{-\beta t}}$ and let $\beta =\beta$ . The Jacobian is $\frac{d\left(\alpha ,\beta \right)}{d\left({\lambda }_{\tau },\beta \right)}=\frac{1}{{\beta }^{2}t{\text{e}}^{-\beta t}}$ . The joint posterior density of $\left({\lambda }_{\tau },\beta \right)$ is given by:

$\pi \left({\lambda }_{\tau },\beta |Y_{obs}\right)=\pi \left(\alpha ,\beta |Y_{obs}\right)\left|\frac{d\left(\alpha ,\beta \right)}{d\left({\lambda }_{\tau },\beta \right)}\right|$

From Equation (12):

$\begin{array}{l}\pi \left({\lambda }_{\tau },\beta |Y_{obs}\right)={\left[k\Gamma \left(n\right)\right]}^{-1}{\alpha }^{n-1}{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\\ ={\left[k\Gamma \left(n\right)\right]}^{-1}{\left(\frac{{\lambda }_{\tau }}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)}^{n-1}{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\left(\frac{{\lambda }_{\tau }}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\\ ={\left(\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)}^{n-1}\frac{{\beta }^{2n-1}}{k\Gamma \left(n\right)}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\lambda }_{\tau }^{n-1}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\\ ={\left(\frac{1}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right)}^n\frac{{\beta }^{2n-1}}{k\Gamma \left(n\right)}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\lambda }_{\tau }^{n-1}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}\end{array}$

$\begin{array}{l}=\frac{{\left({\beta }^2\tau {\text{e}}^{-\beta \tau }\right)}^{-n}}{k\Gamma \left(n\right)}{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}\left\{\frac{\Gamma \left(n\right)}{{\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}^n}\right\}\frac{{\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}^n}{\Gamma \left(n\right)}{\lambda }_{\tau }^{n-1}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}\\ =\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{k{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\frac{{\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}^n}{\Gamma \left(n\right)}{\lambda }_{\tau }^{n-1}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]}\end{array}$ (A6)

Using Equations (A2), (A4), and (A6), we have;

$\begin{array}{c}\gamma =\frac{1}{k}\underset{0}{\overset{\infty }{{\displaystyle \int }}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\left\{1-{\displaystyle {\sum }_{h=0}^{n-1}\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}}\right\}\text{d}\beta \\ =\frac{1}{k}{\displaystyle {\int }_0^{\infty }\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }\\ -\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}\end{array}$

But $k={\displaystyle {\int }_0^{\infty }\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left(1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right)}^n}\text{d}\beta }$

$\gamma =1-\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}$ (A7)

Equation (A7) is the formula in the second part of proposition I.

Proof of Proposition II:

For given level γ, the time required to attain the target value ${\lambda }_{tv}$ is given by ${\tau }^{\ast }=\tau -T$ , where τ satisfies Equation (A2). When β is known, from Equation (A2), it can be noted that $2\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]{\lambda }_{tv}$ follows a Chi-square distribution with 2n degrees of freedom. Therefore, we have:

$2\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}\right]{\lambda }_{tv}={\chi }^2\left(2n;\gamma \right)$ (A8)

$\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}=\frac{{\chi }^2\left(2n;\gamma \right)}{2{\lambda }_{tv}}$ (A9)

${\beta }^2\tau {\text{e}}^{-\beta \tau }=\frac{2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{{\chi }^2\left(2n;\gamma \right)}$

$\tau {\text{e}}^{-\beta \tau }=\frac{2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{{\beta }^2{\chi }^2\left(2n;\gamma \right)}$

$\begin{array}{c}\tau =-\frac{1}{\beta }{W}_n\left(\frac{-2\beta {\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{{\beta }^2{\chi }^2\left(2n;\gamma \right)}\right)\\ =-\frac{1}{\beta }{W}_n\left(\frac{-2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\beta {\chi }^2\left(2n;\gamma \right)}\right),\text{for}n\in \text{Ζ}\end{array}$ (A10)

where $n\in \text{Ζ}$ denotes the n^th root of the equation $\tau {\text{e}}^{\tau }=\frac{-2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\beta {\chi }^2\left(2n;\gamma \right)}$ and W is the Lambert W function., satisfying $W\left(\tau {\text{e}}^{\tau }\right)=\tau$ .

From Equation (A10), we can obtain the time required to attain ${\lambda }_{tv}$ , ${\tau }^{\ast }$ , as follows:

${\tau }^{\ast }=\left[-\frac{1}{\beta }{W}_n\left(\frac{-2{\lambda }_{tv}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\beta {\chi }^2\left(2n;\gamma \right)}\right)\right]-T$ (A11)

Equation (A11) is the first formula of Proposition II.

When β is unknown, the time required to attain the target value, ${\lambda }_{tv}$ , is τ, which is the solution to:

$\gamma =1-\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^2\tau {\text{e}}^{-\beta \tau }}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}$ (A12)

Proof of Proposition III:

The probability is given by

${\gamma }_k=\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]$

When β is known:

${\gamma }_k={\displaystyle {\int }_0^{\infty }\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]\pi \left(\alpha |Y_{obs}\right)\text{d}\alpha }$ (A13)

But $\pi \left(\alpha |Y_{obs}\right)$ is given by Equation (8) and

$\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]={\displaystyle {\sum }_{j=n}^{n+k}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)/f\left(Y_{obs}|\alpha \right)}$ (A14)

From Equation (5), $f\left(Y_{obs}|\alpha \right)={\alpha }^n{\beta }^{2n}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}$ and

$\begin{array}{l}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)\\ ={\displaystyle {\int }_{D\left(j-n:T,\tau \right)}f}\left(Y_{obs},t_{n+1},\cdots ,t_j,N\left(\tau \right)=j\right){\displaystyle {\prod }_{l=n+1}^j\text{d}{t}_l}\\ ={\displaystyle {\int }_{D\left(j-n:T,\tau \right)}{\alpha }^j{\beta }^{2j}\left({\displaystyle {\prod }_{i=1}^j{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^j{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right)e^{-\beta T}\right]}{\displaystyle {\prod }_{l=n+1}^j\text{d}{t}_l}}\\ ={\alpha }^j{\beta }^{2j}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta \tau \right)e^{-\beta \tau }\right]}{\displaystyle {\int }_{D\left(j-n:T,\tau \right)}\left({\displaystyle {\prod }_{l=n+k}^j{t}_l}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{l=n+k}^j{t}_l}}{\displaystyle {\prod }_{l=n+1}^j\text{d}{t}_l}}\end{array}$ (A15)

We solve the integral part as follows:

${\displaystyle {\int }_0^{t}t{\text{e}}^{-\beta t}\text{d}t}=\frac{1}{{\beta }^2}\left[1-\left(1+\beta t\right){\text{e}}^{-\beta t}\right]$

Substituting the limits T and τ, we get; $\frac{1}{{\beta }^2}\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]-\frac{1}{{\beta }^2}\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]$ , which reduces to $\frac{1}{{\beta }^2}\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]$ . Therefore, the integral part of equation (A15) is obtained as;

$\begin{array}{l}{\displaystyle {\int }_{D\left(j-n:T,\tau \right)}\left({\displaystyle {\prod }_{l=n+k}^j{t}_l}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{l=n+k}^j{t}_l}}{\displaystyle {\prod }_{l=n+1}^j\text{d}{t}_l}}\\ =\frac{1}{{\beta }^{2\left(j-n\right)}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}\end{array}$ (A16)

Substituting Equation (A16) into Equation (A15) we obtain;

$\begin{array}{l}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)\\ ={\alpha }^j{\beta }^{2j}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right)e^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}\frac{1}{{\beta }^{2\left(j-n\right)}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}\end{array}$

From Equation (A14);

$\begin{array}{l}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)/f\left(Y_{obs}|\alpha \right)\\ =\frac{{\alpha }^j{\beta }^{2j}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}\frac{1}{{\beta }^{2\left(j-n\right)}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}}{{\alpha }^n{\beta }^{2n}\left({\displaystyle {\prod }_{i=1}^n{t}_i}\right){\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}}\end{array}$

which reduces to;

$\begin{array}{l}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)/f\left(Y_{obs}|\alpha \right)\\ ={\alpha }^{j-n}{\text{e}}^{-\alpha \left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}\end{array}$

Thus, Equation (A14) becomes;

$\begin{array}{l}\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]={\displaystyle {\sum }_{j=n}^{n+k}f\left(Y_{obs},N\left(\tau \right)=j|\alpha \right)/f\left(Y_{obs}|\alpha \right)}\\ ={\displaystyle {\sum }_{j=n}^{n+k}{\alpha }^{j-n}{\text{e}}^{-\alpha \left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}\end{array}$ (A17)

And hence, Equation (A13) becomes;

$\begin{array}{c}{\gamma }_k={\displaystyle {\int }_0^{\infty }{\displaystyle {\sum }_{j=n}^{n+k}{\alpha }^{j-n}{\text{e}}^{-\alpha \left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{\left(j-n\right)!}\\ \times \left\{{\left[\Gamma \left(n\right)\right]}^{-1}{\alpha }^{n-1}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n\right\}\text{d}\alpha \\ ={\displaystyle {\sum }_{j=n}^{n+k}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}{\left(j-n\right)!\Gamma \left(n\right)}}\\ \times {\displaystyle {\int }_0^{\infty }{\alpha }^{j-n}{\text{e}}^{-\alpha \left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}\text{d}\alpha }\\ ={\displaystyle {\sum }_{j=n}^{n+k}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}{\left(j-n\right)!\text{Γ}\left(n\right)}}\\ \times \frac{\text{Γ}\left(j\right)}{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}{\displaystyle {\int }_0^{\infty }\frac{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}{\Gamma \left(j\right)}{\alpha }^{j-n}{\text{e}}^{-\alpha \left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}\text{d}\alpha }\end{array}$ (A18)

The integral part of Equation (A18) is a gamma distribution with parameters j and $\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]$ , and thus integrates to 1. Hence, Equation (A18) reduces to;

${\gamma }_k={\displaystyle {\sum }_{j=n}^{n+k}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n\text{Γ}\left(j\right)}{\left(j-n\right)!\Gamma \left(n\right){\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}}$ (A19)

Equation (A19) is rearranged to obtain;

${\gamma }_k=\frac{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^n}{\displaystyle {\sum }_{j=n}^{n+k}\left(\begin{array}{c}j-1\\ n-1\end{array}\right)}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^j}{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}$ (A20)

Equation (A20) is the first formula of proposition III.

When β is unknown, we obtain:

${\gamma }_k={\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]P\left(\alpha ,\beta |Y_{obs}\right)\text{d}\alpha \text{d}\beta }}$

where $\mathrm{Pr}\left[N\left(\tau \right)\le n+k|Y_{obs}\right]$ and $\pi \left(\alpha ,\beta |Y_{obs}\right)$ are given by Equations (A17) and (12), respectively.

$\begin{array}{l}{\gamma }_k={\displaystyle {\sum }_{j=n}^{n+k}\frac{1}{c\left(j-n\right)!\Gamma \left(n\right)}}{\displaystyle \int {\displaystyle {\int }_0^{\infty }{\alpha }^{j-n}{\text{e}}^{-\alpha \left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}}}\\ \times {\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}{\alpha }^{n-1}{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}{\text{e}}^{-\alpha \left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}\text{d}\alpha \text{d}\beta \\ ={\displaystyle {\sum }_{j=n}^{n+k}\frac{\Gamma \left(j\right)}{c\left(j-n\right)!\Gamma \left(n\right)}}{\displaystyle {\int }_0^{\infty }{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}\frac{{\left[\left(\beta T{\text{e}}^{-\beta T}+{\text{e}}^{-\beta T}\right)-\left(\beta \tau {\text{e}}^{-\beta \tau }+{\text{e}}^{-\beta \tau }\right)\right]}^{j-n}}{{\left[1-\left(1+\beta \tau \right){\text{e}}^{-\beta \tau }\right]}^j}\text{d}\beta }\end{array}$ (A21)

where c = k as used in Equation (12). Letter c has been substituted for k because the summation in Equation (A21) is from n to (n + k), and the k’s are not the same. Equation (A21) implies the second formula in Proposition III.

Proof of Proposition IV:

When β is known, given a predetermined $\tau$ ( $\tau >T$ ), the Bayesian UPL for $\lambda \left(t\right)$ with level γ, denoted by ${\lambda }_U^{\left(\beta \right)}$ satisfies $\gamma =\mathrm{Pr}\left({\lambda }_t\le {\lambda }_U^{\left(\beta \right)}\left(\tau \right)|Y_{obs}\right)$ . From Equations (A2) and (A8);

$\gamma ={\displaystyle {\int }_0^{{\lambda }_U^{\left(\beta \right)}\left(\tau \right)}f}{\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}\right]}^n\frac{{\lambda }_{\tau }^{n-1}}{\Gamma \left(n\right)}{\text{e}}^{-{\lambda }_{\tau }\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}\right]}\text{d}{\lambda }_{\tau }$ (A22)

From Equation (A14);

$2\left[\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}\right]{\lambda }_U^{\left(\beta \right)}\left(\tau \right)={\chi }^2\left(2n;\gamma \right)$ (A23)

Making ${\lambda }_U^{\left(\beta \right)}\left(\tau \right)$ in Equation (A23) the subject, we get;

${\lambda }_U^{\left(\beta \right)}\left(\tau \right)=\frac{\left({\beta }^{2}t{\text{e}}^{-\beta t}\right){\chi }^2\left(2n;\gamma \right)}{2\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}$ (A24)

Equation (A24) implies the first part of proposition IV.

When β is unknown, the Bayesian UPL for $\lambda \left(t\right)=\alpha {\beta }^{2}t{\text{e}}^{-\beta t}$ with level γ is ${\lambda }_{tv}$ , which is the solution to;

$\gamma =1-\frac{1}{k}{\displaystyle {\sum }_{h=0}^{n-1}{\displaystyle {\int }_0^{\infty }\frac{{\left(\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}{\lambda }_{tv}\right)}^h}{h!}{\text{e}}^{-\frac{1-\left(1+\beta T\right){\text{e}}^{-\beta T}}{{\beta }^{2}t{\text{e}}^{-\beta t}}{\lambda }_{tv}}\frac{{\beta }^{2n-1}{\text{e}}^{-\beta {\displaystyle {\sum }_{i=1}^n{t}_i}}}{{\left[1-\left(1+\beta T\right){\text{e}}^{-\beta T}\right]}^n}\text{d}\beta }}$ (A25)

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Kaur, R. and Panwar, P. (2015) Study of Perfect and Imperfect Debugging NHPP SRGMs Used for Prediction of Faults in a Software. IJCSC, 6, 73-78.
[2]	Akuno, A.O., Orawo, L.A. and Islam, A.S. (2014) One-Sample Bayesian Predictive Analyses for an Exponential Nonhomogeneous Poisson Process in Software Reliability. Open Journal of Statistics, 4, 402-411. https://doi.org/10.4236/ojs.2014.45039
[3]	Yamada, S., Ohba, M. and Osaki, S. (1984) S-Shaped Software Reliability Growth Models and Their Applications. IEEE Transactions on Reliability, R-33, 289-292. https://doi.org/10.1109/TR.1984.5221826
[4]	Hanagal, D.D. and Bhalerao, N.N. (2018) Analysis of Delayed S-Shaped Software Reliability Growth Model with Time Dependent Fault Content Rate Function. Journal of Data Science, 16, 857-878. https://doi.org/10.6339/JDS.201810_16(4).00010
[5]	Lee, T.Q., Yeh, C.W. and Fang, C.C. (2014) Bayesian Software Reliability Prediction Based on Yamada Delayed S-Shaped Model. Applied Mechanics and Materials, 490, 1267-1278. https://doi.org/10.4028/www.scientific.net/AMM.490-491.1267
[6]	Collins, O., Akong’o, O.L., Munene, M.G. and Okenye, J.O. (2023) Bayesian Interval Estimation in a Non-Homogeneous Poisson Process with Delayed S-Shaped Intensity Function. American Journal of Theoretical and Applied Statistics, 12, 43-50.
[7]	Yin, L. and Trivedi, K.S. (1999) Confidence Interval Estimation of NHPP-Based Software Reliability Models. Proceedings 10th International Symposium on Software Reliability Engineering, 6-11.
[8]	Yu, J.W., Tian, G.L. and Tang, M.L. (2007) Predictive Analyses for Nonhomogeneous Poisson Processes with Power Law Using Bayesian Approach. Computational Statistics & Data Analysis, 51, 4254-4268. https://doi.org/10.1016/j.csda.2006.05.010
[9]	Yu, J.W., Tian, G.L. and Tang, M.L. (2011) Bayesian Estimation and Prediction for the Power Law Process with Left Truncated Data. Journal of Data Science, 9, 445-470. https://doi.org/10.6339/JDS.201107_09(3).0009
[10]	Cheruiyot, N., Orawo, L.A.O. and Islam, A.S. (2018) Bayesian Predictive Analyses for Logarithmic Non-Homogeneous Poisson Process in Software Reliability. Open Access Library Journal, 5, 1-13. https://doi.org/10.4236/oalib.1104767
[11]	Ehrlich, W., Prasanna, B., Stampfel, J. and Wu, J. (1993) Determining the Cost of a Stop-Test Decision (Software Reliability). IEEE Software, 10, 33-42. https://doi.org/10.1109/52.199726
[12]	Hung-Cuong, N. and Quyet-Thang, H. (2015) New NHPP SRM Based on Generalized S-Shaped Fault-Detection Rate Function. Nature of Computation and Communication: International Conference, ICTCC 2014, Ho Chi Minh City, 24-25 November 2014, 212-221. https://doi.org/10.1007/978-3-319-15392-6_21