Prediction in trend-renewal processes for repairable systems

Some problems of point and interval prediction in a trend-renewal process (TRP) are considered. TRP’s, whose realizations depend on a renewal distribution as well as on a trend function, comprise the non-homogeneous Poisson and renewal processes and serve as useful reliability models for repairable systems. For these processes, some possible ideas and methods for constructing the predicted next failure time and the prediction interval for the next failure time are presented. A method of constructing the predictors is also presented in the case when the renewal distribution of a TRP is unknown (and consequently, the likelihood function of this process is unknown). Using the prediction methods proposed, simulations are conducted to compare the predicted times and prediction intervals for a TRP with completely unknown renewal distribution with the corresponding results for the TRP with a Weibull renewal distribution and power law type trend function. The prediction methods are also applied to some real data.


Introduction
We consider a class of point and interval prediction problems for stochastic models determined by the trend-renewal process (TRP) which is defined to be a time-transformed renewal process (RP), where the time transformation is given by a trend function λ(·). Most commonly used reliability models for repairable systems, such as non-homogeneous Poisson process (NHPP) and RP are special cases of the TRP.
The TRP was introduced and investigated first by Lindqvist (1993) and by Lindqvist et al. (1994) (see also Doksum 2003 andLindqvist 2006). Parametric inference on the parameters of the TRP was considered in the paper of Lindqvist et al. (2003), where the authors also proposed corresponding models, called heterogeneous trend-renewal processes, that extend the TRP to cases involving unobserved heterogeneity.
General approaches to predictive inference are considered by Barndorff-Nielsen and Cox (1996), Smith (1999) and Lawless and Fredette (2005). Prediction problems were considered for some special cases of a TRP. Exact prediction intervals, based on maximum likelihood estimation, for the kth future observation of the non-homogeneous Poisson process with power law intensity (the Weibull process) were derived by Engelhardt and Bain (1978). In the paper of Calabria and Pulcini (1996), prediction limits on the current system lifetime modeled by power law process have been derived both in the maximum likelihood and Bayesian context. The prediction problem for the Weibull process with incomplete observations is considered in the paper of Yu et al. (2008). Recently, Helmers and Mangku (2012) construct and investigate a (1 − α)-upper prediction bound for a future observation of a cyclic Poisson process.
In the present paper we take up the problem of point and interval prediction for the next failure time of a TRP in the case when the renewal distribution of the process is unknown. If the renewal distribution is unknown, then the likelihood function of the TRP is unknown and consequently the ML method of estimating unknown parameters cannot be used. We present some ideas of predicting the next failure time and constructing the prediction intervals for the next failure time in such models. The main idea of prediction consists in finding at the first step the estimates of unknown trend parameters by minimizing the sample variance of the transformed working times (one of the methods of estimating proposed by Jokiel-Rokita and Magiera (2012) for the TRP model) and then in exploiting the inverse transformation of the cumulative trend function of the TRP. The method proposed allows to predict the next failure time and construct the prediction interval for the next failure time in the case when the renewal distribution of the process is completely unknown. Using the prediction methods proposed, the predicted times and prediction intervals for a TRP with completely unknown renewal distribution are compared in a simulation study with the corresponding results for the TRP with a Weibull renewal distribution and power law type trend function. The behavior of the predictions proposed is also examined in the TRP model with a Gompertz renewal distribution and the same power law type trend function. The prediction methods are applied to some real data too.
The article is organized as follows. Some methods determining the predicted next failure time are presented in Sect. 3. Prediction intervals for the next failure time are constructed in Sect. 4. The method of constructing the predicted times and prediction intervals in the case when the form of the renewal distribution function F is unknown is proposed in Sects. 3.3 and 4.2, respectively. Using various methods proposed for determining the predicted next failure time and the prediction intervals, the performance of the predicted times and intervals is illustrated in Sect. 5 for simulated TRP models. The predicted next failure times and the prediction intervals constructed for a TRP with completely unknown renewal distribution function F and power law trend function (power law TRP) are compared with the predictors in the reference to the TRP models with F corresponding to the Weibull (Sect. 5.1) and Gompertz (Sect. 5.3) renewal distributions and the same trend function. In particular, for the power law process, the results are illustrated in Sect. 5.2 for five prediction intervals constructed by various methods. In Sect. 6 the results concerning the determination of the predicted next failure times and the prediction intervals in a TRP model are applied to some real data. Section 7 contains concluding remarks.

Definitions and preliminaries
Let N(t) denote the number of jumps (failures) in the time interval (0, t] and let T i be the time of the ith failure. De-fine T 0 = 0 and denote by X i = T i − T i−1 , i = 1, 2, . . . , the time between the failure number i − 1 and the failure number i (the so called working time or waiting time). In the context of failure-repair models it is assumed that all repair times are equal to 0. In practice this corresponds to the situation, when repair actions are conducted immediately or the repair times can be neglected with comparison to the times X i between the failures. A review of reliability models for repairable systems is given by Hollander and Sethuraman (2002).
The observed sequence {T i , i = 1, 2, . . .} of occurrence times T 1 , T 2 , . . . (failure times) forms a point process, and {N(t), t ≥ 0} is the corresponding counting process. The two point processes, non-homogeneous Poisson process (NHPP) and RP are widely investigated in the literature on reliability to model minimal repairs and perfect repairs (renewals), respectively. The process {N(t), t ≥ 0} is an RP if the random variables X 1 , X 2 , . . . are independent and identically distributed with cumulative distribution function (cdf) F with F (0) = 0. The RP with the renewal function F will be denoted by RP(F ). If F is the cdf of the exponential distribution E(λ), then RP(F ) is a homogeneous Poisson process with intensity λ, denoted by HPP(λ).
Let λ(t), t ≥ 0, be a nonnegative function, and let Λ(t) = t 0 λ(u)du. The process {N(t), t ≥ 0} is called a TRP with a renewal distribution function F (t) and a trend function λ(t) if the time-transformed process Λ(T 1 ), Λ(T 2 ), . . . is an RP(F ) with the renewal distribution function F , i.e. if the random variables (the so called transformed working times) are i.i.d. with cdf F . The TRP will be denoted by TRP(F, λ(·)).
The class of TRP's is defined by properties of the sequence {T i , i = 1, 2, . . .} and includes the NHPP's and RP's. It also includes modulated power law processes. Equivalently, the corresponding counting process {N(t), t ≥ 0} can be considered, where N(t) = N(Λ(t)) and { N(t), t ≥ 0} represents a RP.
Note that the representation TRP(F, λ(·)) is not unique. For uniqueness it is assumed that the expected value of the renewal distribution defined by F equals 1.
For a TRP(F, λ(·)) the conditional intensity function has the following form where z(t) is the hazard rate corresponding to F : Andersen et al. (1993), see also Lindqvist et al. (2003, formula 2, ch. 2), the likelihood function of a TRP(F, λ(·)) observed in the interval time [0, τ ] with the realizations t 1 , t 2 , . . . , t N(τ ) of the jump (failure) times T 1 , T 2 , . . . , T N(τ ) can be expressed in the form (1) and the log-likelihood function is defined by For a random stopping time τ with respect to the filtration F t = σ {N(u) : u ≤ t}, this formula follows from the fundamental identity of sequential analysis as a consequence of the optional stopping theorem. We consider the case when τ = inf{t ≥ 0 : N(t) = n}, i.e., N(τ ) = n, τ = T n , where n is a given number of failures. This means we suppose that a TRP(F, λ(·)) is observed up to the nth event (failure) appears for the first time, and the values of the jump times T 1 , . . . , T n are recorded. In other words, we consider the so called failure truncation (or inverse sequential) procedure. In this case the log-likelihood function for a TRP(F, λ(·)) takes the form The problem is to predict the next failure time T n+1 and to construct the prediction intervals for T n+1 on the basis of the failure times observed up to time T n .

General ideas
Let us consider a TRP(F, λ(·)) with trend function (intensity) λ(t) and renewal distribution F . Recall, that for the uniqueness of the TRP it is assumed that the renewal distribution defined by F has the expectation value 1. Thus, The condition (2) prompts to the following predicted future time T n+1 : (Jokiel-Rokita and Magiera 2012, Remark 2). Let us note that the prediction formula given by (3) can be used for the TRP model regardless of the renewal distribution is known or it is not known. If the renewal distribution F is known then the conditional density f T n+1 |T n =t n is known and predicting T n+1 in a TRP(F, λ(·)) can be based on the formula where is the conditional density function of the random variable T n+1 , given T n = t n , where z(·) is the hazard rate function corresponding to F . The predicted next failure time defined by (4) is the conditional expectation E(T n+1 |T n = t n ). The point prediction for a TRP(F, λ(·)) can be also based on the following known fact.
We then have the following proposition.

Proposition 1
The conditional next failure times T (i) n+1 , given T n , i = 1, . . . , m, for the TRP(F, λ(·)) with the cumulative conditional intensity function G(t) can be generated according to the formula where V i , i = 1, . . . , m, are the independent random variables simulated from the exponential distribution E(1), and predict the conditional next failure time T n+1 , given T n , by which is an approximation of (4).
Remark 1 Using the fact that Λ(T n+1 ) − Λ(T n ) is a random variable with distribution F , one can generate T (i) n+1 equivalently according to the formula T where Z i , i = 1, . . . , m, are the independent random variables simulated from the distribution F .
In practice, usually the trend function λ(t) and the renewal distribution F are not exactly known, and consequently to predict the next failure time T n+1 formulas (3) and (4) can not be used. In Sect. 3.2 we consider the problem of predicting the next failure time T n+1 in the TRP(F, λ(·)) when the form of trend function λ(t) = λ(t; ϑ) and the form of the renewal function F (t) = F (t; υ) are known, but ϑ and υ are unknown parameters. In Sect. 3.3 we consider the problem of prediction the next failure time T n+1 in the TRP(F, λ(·)) assuming that the parametric form of the trend function λ(t) = λ(t; ϑ) is known, but we make no assumptions on the form of the renewal distribution F .

Point prediction of T n+1 in a TRP when F is known
If a form of the trend function λ(t) = λ(t; ϑ) and a form of the renewal function F (t) = F (t; υ) are known, then consequently the form of the likelihood function is known and we can obtain maximum likelihood estimators (MLE's) ϑ ML and υ ML of the unknown parameters ϑ and υ.
To predict the next failure time T n+1 one can use formula (3) with Λ(t; ϑ ML ) instead of the unknown function Λ. Let us note that applying this method of prediction does not require the knowledge of an estimate of the parameter υ.
Another possibility for predicting T n+1 consists in the mean prediction defined by formula (4) with f T n+1 |T n =t n (t) = f T n+1 |T n =t n (t; θ), where θ = ( υ ML , ϑ ML ), instead of the unknown function f T n+1 |T n =t n (t). Let us notice that applying this method of prediction requires the knowledge of estimates of both the parameters ϑ and υ.
The prediction based on Proposition 1 also requires the knowledge of estimates of both the parameters ϑ and υ, and can be applied in the case when F is known.
For the WPLP(α, β, γ ) the conditional density function defined by (5) takes the following form and the corresponding distribution function is Note that ϕ = α for the PLP(α, β).
In the failure truncation procedure considered the likelihood function is given by and the log-likelihood function is (n; θ) = n(ln ϕ + ln β + ln γ ) We give in Proposition 2 below the formulas determining the ML estimators of α, β and γ of the WPLP(α, β, γ ) in the failure truncation procedure. In general case, for any stopping time, such formulas are given in the paper of Jokiel-Rokita and Magiera (2012).

Proposition 2
The ML estimators ϕ ML , β ML and γ ML of the parameters ϕ, β and γ in the failure truncation procedure for the WPLP(α, β, γ ) are determined as follows: where β ML and γ ML are the solutions of the following system of likelihood equations An estimator α of α is evaluated according to the formula where ϕ and γ are estimators of ϕ and γ .
The predicted next failure time Basing on formula (3) we obtain the formula for the predicted next failure time.

Proposition 3 The point predictor of T n+1 in the WPLP(α, β, γ ) model based on the ML estimates α ML , β ML and γ ML in the failure truncation procedure is determined by
where Let us notice that in the case of the WPLP(α, β, γ ) formula (6) takes the form . . , m, are the random numbers from the distribution E(1). Thus we have the following proposition.

Proposition 4
The conditional next failure time T n+1 , given T n , for the WPLP(α, β, γ ) can be predicted by where Z i , i = 1, . . . , m, are the random numbers from the distribution E(ϕ).
If the parameters α, β and γ of the WPLP(α, β, γ ) model are not known, then basing on Proposition 4 we can predict the next failure time T n+1 by where Z i , i = 1, . . . , m, are the random numbers from the distribution E( ϕ ML ).

Point prediction in the PLP(α, β)
For the WPLP(α, β, 1), i.e. for the PLP(α, β), the ML estimators of α and β can be explicitly determined (see e.g. Rigdon and Basu 2000, pp. 136-137). The ML estimators of α and β for this process are defined by By Proposition 3, using formula (16) we obtain the following form of the predicted next failure time for the PLP(α, β) In the case of PLP(α, β) we have from Proposition 4 the following form of the conditional predicted next failure time where Z i , i = 1, . . . , m, are the random numbers from the distribution E( α P L ML ).

Point prediction in the HPP(α)
For the HPP(α) the ML estimator of α is and the predicted next failure time is determined by where Z i , i = 1, . . . , m, are the random numbers from the distribution E( α P ML ).

Point prediction of T n+1 in a TRP when F is unknown
If a form of the renewal distribution function F of a TRP(F, λ(·)) is unknown, we cannot determine the ML estimators of the process parameters θ and use them in formulas (3) or (4) to determine the predictors of T n+1 . However, in this case to estimate a parameter ϑ of a trend function one can use one of the three methods proposed by Jokiel-Rokita and Magiera (2012): the least squares method, the method of moments and the constrained least squares (CLS) method. The investigations carried out for the WPLP(α, β, γ ) in their paper have demonstrated the advantage of the CLS method over the other methods for most cases of possible range of the process parameters. In the CLS method we define the sum of squares where t i are the realizations of random variables T i , i = 1, . . . , N(τ ), up to time τ , and Λ(t 0 ) := 0. The problem is to find where the restriction set is Remark that The variables W i = Λ(T i ) − Λ(T i−1 ) are the observations from the distribution with the expected value 1. Thus the CLS method consists of deriving such estimate ϑ CLS of the unknown parameter ϑ (of the trend function) which minimizes the sum of squares of deviations of the random variables W i from the expected value 1 (equivalently, which minimizes the sample variance) under the constraint that the sample mean W equals the theoretical expected value of the distribution defined by F . Let us notice that in the case when F is unknown we can predict the next failure time T n+1 using formula (3) only with Λ = Λ(t; ϑ CLS ).
To derive the CLS estimates of the parameters α and β of the PTRP(α, β) one can use the following proposition (Jokiel-Rokita and Magiera 2012).

Proposition 5
The CLS estimators α CLS and β CLS of α and β in the PTRP(α, β) observed in the interval time [0, τ ] are determined by For the PTRP(α, β) we have Λ(T i ) = Λ(T i ; α, β) = αT β i , and using the method defined by (3) we have the following proposition for the predicted next failure time.

Proposition 6
The predicted next failure time obtained by the CLS method in the PTRP(α, β) model is given by In the CLS method, the failure truncation procedure will be applied in which α CLS = n/t

Interval prediction in a TRP
In a TRP(F, λ(·)) with cumulative intensity function Λ(t) = t 0 λ(u)du and renewal distribution function F the transformed working times Then the prediction interval for T n+1 can be determined from the condition , and 1 − ε is the given coverage probability (ε ∈ (0, 1)). Thus we have the following proposition.
Proposition 7 In a TRP(F, λ(·)) the lower and upper bounds T L and T U , respectively, of the prediction interval for the next failure time T n+1 are given by Equivalently, the prediction interval is defined by In particular, in the PTRP(α, β), In the WPLP(α, β, γ ) we can obtain the exact form of the quantiles q F (ε 1 ) and q F (1 − ε 2 ) (F has the Weibull distribution We(γ , 1/Γ (1 + 1/γ )) ≡ We(γ , α/ϕ 1/γ )) and the lower and upper bounds T L and T U are defined by where ϕ is given by (7). In the case when the renewal distribution F is unknown, we propose to estimate the unknown quantiles q F (ε 1 ) and q F (1 − ε 2 ) by the sample quantiles q(ε 1 ) and q(1 − ε 2 ) from the sample (W 1 , . . . , W n ). Thus the estimated lower and upper bounds, T L and T U , are given by Another way of obtaining the prediction interval I = [t L , t U ] for T n+1 , given T n = t n , is to solve the equations with respect to t L and t U , where ε 1 + ε 2 = ε, i.e. t L and t U are quantiles of order ε 1 and 1 − ε 2 of the distribution defined by the density f T n+1 |T n =t n (t) given by (5).

Interval prediction in a TRP when F is known
If a form of the trend function λ(t) = λ(t; ϑ) and a form of the renewal function F (t) = F (t; υ) are known, then consequently the form of the likelihood function is known and we can obtain ML estimators ϑ ML and υ ML of the unknown parameters ϑ and υ, respectively. Thus, by Proposition 7 we have the following result.

Proposition 8
The estimated lower and upper bounds T L and T U , respectively, of the prediction interval for the next failure time T n+1 in a TRP(F, λ(·)) are given by where Λ(t) = Λ(t; ϑ ML ) and F (t) = F (t; υ ML ).
We can also obtain the estimated conditional lower and upper bounds for the next failure time T n+1 by solving equations (23) with respect to t L and t U for f T n+1 |T n =t n (t) = f T n+1 |T n =t n (t; θ), where θ = ( υ ML , ϑ ML ), instead of the unknown function f T n+1 |T n =t n (t).

Interval prediction in the WPLP(α, β, γ )
Basing on formula (22) for the lower and upper bounds for the next failure time T n+1 in the WPLP(α, β, γ ) model with known parameters α, β and γ , we have the following proposition for the WPLP(α, β, γ ) with the unknown parameters.

Proposition 9
The estimated lower and upper bounds T L and T U , respectively, of the prediction interval for the next failure time T n+1 in the WPLP(α, β, γ ) are given by where ϕ ML , β ML , γ ML are estimators of the parameters ϕ, β, γ , and they are defined in Proposition 2.

Interval prediction in the PLP(α, β)
It follows from Proposition 9 that for the PLP(α, β) which is the WPLP(α, β, 1) the estimated prediction interval for the next failure time T n+1 and based on the ML estimators is determined by where α P L ML and β P L ML are defined by (16). For the PLP(α, β) the transformed working times In the case when the parameter β is known, introducing the pivotal function which has Fisher's F -distribution with (2, 2n) degrees of freedom, one can determine the prediction interval for T n+1 from the condition where F 2,2n (q) is the quantile of order q of Fisher's F -distribution F (2, 2n) with (2, 2n) degrees of freedom, i.e., F 2,2n (q) = n[(1 − q) −1/n − 1]. The lower and upper bounds of the confidence interval for T n+1 are then determined by the inequalities i.e., If the parameter β is unknown, one can consider the estimated prediction interval obtained by substituting the parameter β by its ML estimate into formula (26). This is a very satisfactory way when n is sufficiently large, because it can be shown that both Q and have the same exponential E(1) limit distribution. Thus we have the following proposition.

Proposition 10
The estimated lower and upper bounds, based on the ML estimator of β,T P L L andT P L U , of the prediction intervalȊ P L for the next failure time T n+1 in the PLP(α, β) are given by where β P L ML is defined by (16). Engelhardt and Bain (1978) gave the prediction interval for the kth future failure time T n+k , k = 1, 2 . . . , by considering the pivotal function when β is unknown, namely by considering the pivotal function based on T n and the ML estimator of β. It follows from their result for k = 1 that the random variable where β P L ML is the ML estimator of β defined by (16), has the distribution function defined by P (Y ≤ y) = 1 − (1 + y/(n − 1)) −(n−1) . Taking Y as a pivotal function for T n+1 , we have P (q ε 1 ≤ Y ≤ q 1−ε 2 ) = 1 − ε, i.e., where q ε denotes the quantile of order ε of the distribution of the random variable Y , i.e., q ε satisfies the condition 1 − (1 + q ε /(n − 1)) −(n−1) = ε. Thus, the prediction interval I EB based on the idea of Engelhardt and Bain (1978) is defined by Remark 2 Let us notice that the statistic Q can be written in the form Thus, Q is a pivotal function equivalent to Y . The distribution function of Q is given by .
Taking in Proposition 10 the exact quantiles of the distribution determined by F Q instead of quantiles of the F (2, 2n) distribution we obtain the prediction interval with the lower and upper bounds given by formula (28).
In the simulation study of Sect. 5.2 we compare the exact prediction interval I EB with the prediction intervalȊ P L defined in Proposition 10.

Interval prediction in the HPP(α)
For the HPP(α) which is the WPLP(α, 1, 1) one obtains from Proposition 9 the following form of the estimated lower and upper bounds, T P L and T P U , of the prediction interval for the next failure time T n+1 In view of formula (26) with β = 1 we have the following corollary which gives the form of exact prediction intervals for T n+1 .

Corollary 1
The lower and upper bounds, T L and T U , of the prediction interval for the next failure time T n+1 in the HPP(α) are given by

Interval prediction in a TRP when F is unknown
Denote W = ( W 1 , . . . , W n ), where W i , i = 1, . . . , n, are estimators of the transformed working times W i , i.e., If the form of F in a TRP(F, λ(·)) is unknown, to determine the estimated lower and upper bounds for T n+1 we propose to use sample quantiles of W with ϑ = ϑ CLS instead of unknown quantiles q F . Namely, one can use the following proposition.

Proposition 11
The estimated lower and upper bounds T L and T U , respectively, of the prediction interval for the next failure time T n+1 in a TRP(F, λ(·)) with unknown renewal distribution function F are given by where Λ(t) = Λ(t; ϑ) and q n ( W ; ε) denotes the sample quantile of order ε of W .
To obtain an estimate ϑ of the unknown trend parameter ϑ one can apply the CLS method.
In particular, for the PTRP(α, β) model we have the following corollary.

Corollary 2 The estimated lower and upper bounds, T CLS L
and T CLS U , of the prediction interval I CLS for the next failure time T n+1 in the PTRP(α, β) are given by where q n ( W ; ε) denotes the sample quantile of order ε of

Simulation study
The main purpose of the simulation study is to examine how much the predicted next failure times and the prediction intervals constructed by the CLS method for a TRP(F, λ(·)) model with an unknown distribution function F differ from the predictors constructed by the ML method in the TRP(F, λ(·)) with known F and the same trend function λ(·). As a reference TRP(F, λ(·)) model, the WPLP(α, β, γ ) is taken into account. The CLS estimates α CLS and β CLS of the parameters α and β, and consequently the resulting predictors for the corresponding PTRP(α, β) are evaluated on the basis of the realizations of the generated WPLP(α, β, γ ) supposing that we do know nothing about the renewal distribution function F , i.e., that we observe the PTRP(α, β).
Each sample of the WPLP(α, β, γ ) is generated up to a fixed number n + 1 of jumps is reached. For each chosen combination of the parameters α, β and γ the samples of size k of the realizations are generated. In the case γ = 1 and in the case β = 1, γ = 1 the realizations of the PLP(α, β) and HPP(α) are obtained, respectively.
The 'real' last failure time T n and all the predictors are evaluated as the means from the k estimates, such that each of these estimates was derived on the basis of the individual realization up to nth jump (failure) of the process considered. Analogously, the 'real' next failure time T n+1 denotes the mean resulting from the same k repetitions of the realizations involving the (n + 1)-th jump time. The simulation study was carried out for n = 100 and k = 200.
The accuracy of a predictor, say η of η, is measured by the relative error (RE) which is defined by RE( η ) = ( η − η)/η, and by the variability of a predictor η which is determined by the root mean squared error RMSE( η ) = (sd( η )) 2 + (mean( η ) − η) 2 , where sd stands for the standard deviation. In tables the abbreviations re = re( η ) and se = se( η ) is used for the RE( η ) and RMSE( η ) errors, respectively. The relative errors of predictors are given in percentages.

The WPLP model
The WPLP(α, β, γ ) process is generated according to the following formula for the jump times: T 0 = 0, where U i are random numbers from uniform distribution U(0, 1).
The computer program was developed using Mathematica 8.04 package.
The CLS estimates α CLS and β CLS in the PTRP(α, β) model are evaluated on the basis of Proposition 5 and using a constrained global optimization (CGO) method. The ML estimates α ML , β ML and γ ML in the WPLP(α, β, γ ) model are evaluated by finding at first the estimates β ML and γ ML as the solution to the optimization problem ( β ML , γ ML ) = arg max (β,γ ) n; (ϕ, β, γ ) , by using the CGO procedure, where (n; (ϕ, β, γ )) is the log-likelihood function given by (8) with ϕ = ϕ(β, γ ) defined by (11). The ML estimate α ML of α is then evaluated using formula (12) with ϕ defined by (9). In the computer program, to solve the optimization problems the procedure NMaximize of Mathematica package is implemented. We do not attach the tables with values of the parameter estimates for the models considered. The estimators in WPLP and PTRP models were investigated by Jokiel-Rokita and Magiera (2012).
Tables 1, 2, 3, and 4 provide the values related to the predictors in the WPLP(α, β, γ ) model in comparison to a PTRP(α, β) (the TRP with unknown F ) for the number of failures n = 100.
Tables 1 and 2 contain the values of the predicted next failure times and their errors in the WPLP(α, β, γ ) and PTRP(α, β) models compared. The predictor T CLS n+1 in the PTRP(α, β) model, given by formula (21), is compared with the predicted next failure times T n+1 , T ML n+1 , T ML n+1 , given by formulas (14), (13) and (15), respectively, obtained by the other ideas described in the paper and using the ML estimators for the WPLP(α, β, γ ). The prediction T n+1 means the conditional theoretical predicted next failure time in the WPLP(α, β, γ ) evaluated by formula (14) for theoretical values of α, β and γ . The random numbers Z i , i = 1, . . . , m, are generated from the distribution E((αΓ (1 + 1/γ )) γ ) for each realization of the process. There are m = 1000 simulations used for computing the predictions T n+1 and T ML n+1 . In Tables 3 and 4    respectively, are evaluated as the frequency of appearing the next failure time T n+1 in the prediction intervals, are examined too. The cases 6-9 and 10-12 of Tables 1-4 correspond to the PLP and to the HPP model, respectively.
Remarks on the simulation results It follows from Tables 1-4 that the predicted next failure times and the prediction intervals constructed by the CLS method for a TRP(F, λ(·)) model with an unknown distribution function F differ very slightly from the predictors constructed by the ML method in the TRP(F, λ(·)) with known F and the same trend function λ(·). This indicates that in the TRP model considered the lack of knowledge of the renewal distribution has an acceptable vague influence on the accuracy of determining all the predictors of the next failure time T n+1 as well as on the accuracy of all the prediction intervals for this future breakdown time. The predictors based on the CLS estimators have practically the same errors (RE's percentages and RMSE's) as those based on the ML estimators in the model with known F (see Table 2).
The values of Table 4 show that in determining the prediction intervals the CLS method applied to the TRP with unknown renewal distribution function F provides very satisfactory results.
For any fixed γ all the errors of the predictors in the WPLP and PTRP models decrease for greater values of β and smaller values of α.
For a given number n of failures, the RMSE's of all the parameter estimators in the WPLP(α, β, γ ), and consequently of all the predictors become significantly smaller as the parameter γ increases because the variability of the WPLP realizations evidently decreases as γ increases according to the variance formula of the renewal distribution (see a remark in the paper of Jokiel-Rokita and Magiera 2012). Therefore, for the same pairs of (α, β) the RE percentages and the RMSE's of the predicted times become significantly smaller for greater γ (see Table 2).

The PLP model
Tables 5-9 provide the values connected with the predictors for the PLP(α, β) special reference model for which the ML predictors were evaluated by using explicit formulas (16) for the ML estimators of the parameters α, β.
The point predictors T n+1 , T P L n+1 , T P L n+1 and T CLS n+1 are calculated from formulas (14) The lengths L, L P L ,L P L , L EB of the prediction intervals I, I P L ,Ȋ P L , I EB and the coverage frequencies CF(I), CF( I P L ), CF(Ȋ P L ) and CF( I EB ) for the PLP(α, β) are given in Table 8.
The prediction intervals I CLS = [ T CLS L , T CLS U ] for the PLP(α, β), obtained by the CLS method, and their coverage frequencies CF( I CLS ) are presented in Table 9.
Remarks on the simulation results The values of Tables 8 and 9 demonstrate that all the coverage frequencies CF( I P L ), CF(Ȋ P L ), CF( I EB ) and CF( I CLS ) for the PLP(α, β) are very close to the theoretical coverage frequency CF(I), but it is seen from Table 8 that for all the cases of α and β of the PLP(α, β) the lengths of the prediction intervals can be evidently ordered L CLS < L EB and for all the cases of α and β of the PLP(α, β) (except β = 1, α = 0.5): L P L < L <L P L < L EB . The value of Tables 7, 8, and 9 show that the prediction interval I CLS (obtained by the CLS method) corresponds numerically very close to the prediction intervalȊ P L (based on the approximated F (2, 2n) distribution).
Comparing the values of the predictors in the cases 6-12 of Tables 1-4 with the corresponding values in Tables 5-9 we infer that the methods used for constructing the predicted next failure times and the prediction intervals, and based on the CLS estimators in a TRP model with unknown distribution function F (the PTRP(α, β) model), provide the results which in a very satisfactory way match the predictions for the PLP(α, β) as well as for the HPP(α) model with parameter estimators of α and β obtained on the basis of exact formulas.

The GMPLP model
We have considered also the TRP(F, λ(·)) with the same power law trend function λ(t) and with the renewal distribution function where Γ (a, z) = ∞ z u a−1 exp(−u)du denotes the incomplete gamma function. The renewal distribution function F corresponds to the Gompertz distribution Gom(exp(η)Γ (0, η), η) with the parametrization resulting in the expectation 1. The hazard function corresponding to F has the exponential form z(x) = ηΓ (0, η) exp[exp(η)Γ (0, η)x + η]. We call this process the Gompertz power law process and denote it by GMPLP(α, β, η).
The GMPLP(α, β, η) is generated according to the following formula for the jump times: Table 5 The predicted (n + 1)-th failure times T n+1 , T P L n+1 and T P L n+1 for the PLP(α, β), and the predicted next failure time T CLS n+1 for the PTRP(α, β)   199.3015 206.6012 199.3004 206.4344 199.3004 206.5685 199.3009 206.7859 Table 8 The lengths L, L P L , L P L , L EB of the prediction intervals I, I P L ,Ȋ P L , I EB and the coverage frequencies CF(I), CF( I P L ), CF(Ȋ P L ) and CF( I EB ) for the PLP(α, β) ] are evaluated regarding the GMPLP(α, β, η) realization as that of the WPLP(α, β, η) one. We examine, whether in this case one can expect that the CLS predictions should be more accurate than the ML ones. The predictions T MLW n+1 , T CLS n+1 , I MLW and I CLS are given in Tables 10 and 11.

Remarks on the simulation results
We observe that if we apply the ML method, under the WPLP(α, β, η) assumption, to simulated data from the GMPLP(α, β, η) the vast majority of the cases gives more accurate results by the CLS method: the coverage frequencies are almost the same but the interval lengths are smaller in all the cases of parameters.
The simulation study shows that if the ML method does not match the TRP model considered, then the CLS method can give more accurate predictions. For example, if we apply the ML method under the WPLP assumption to data from a PTRP which is not a WPLP, then we can expect more accurate predictions by the CLS method.

Application to some real data set
In this section we consider the problem of point and interval prediction for some real data of failure times, namely for the data set contained in the paper of Lindqvist et al. (2003), and which is given here in Table 12. These data contain 41 failure times of a gas compressor with time censoring at time 7571 (days). Although the data have been observed under the time truncation procedure, we however assume that we only know the first n = 40 failure times of the data set and want to find the point and interval predictions for the next failure time.
Supposing that the set of failure times of Table 12 forms a TRP belonging to the class of WPLP(α, β, γ ), the ML estimates α ML , β ML and γ ML of α, β and γ have been evaluated (on the basis of the previous n = 40 failure times) and presented in Table 13. On the other hand, if no assumptions are made on the renewal distribution function F , the estimates α CLS and β CLS of α and β are given as the parameters of the PTRP(α, β). Consequently, all the predictions for Table 12 The real data  1  4  305  330  651  856  996  1016  1155  1520  1597  1729   1758  1852  2070  2073  2093  2213  3197  3555  3558  3724  3768  4103   4124  4170  4270  4336  4416  4492  4534  4578  4762  5474  5573  5577   5715  6424  6692 6830 6999 T n+1 are evaluated on the basis of only the first n = 40 failure times and then the predictors of T 41 are compared with the observed real last value. In Table 13 the relative errors re( α CLS ) = | α CLS − α ML |/ α ML and re( β CLS ) = | β CLS − β ML |/ β ML are given too. For comparison, in Table 13 there are also given the sum of squares SS CLS := S 2 LS ( α CLS , β CLS ) and SS ML := S 2 LS ( α ML , β ML ), where S 2 LS (ϑ) is defined by (20). The predicted next failure times T ML n+1 , T ML n+1 , T CLS n+1 and the prediction intervals I ML , I CLS for various confidence level (CL) percentages are given in Tables 14 and 15, respectively. The predictions T ML n+1 , T ML n+1 , T CLS n+1 , I ML and I CLS are evaluated according to formulas (13), (15), (21), (24) and (29), respectively. In Table 14 the percentages of the relative errors re( T n+1 ) = ( T n+1 − T n+1 )/T n+1 are given too.
Remarks on the results of the analysis In analyzing the real data of Table 12 we observe in Table 13 a significant relative difference between the estimates α ML and α CLS for a WPLP model and for a PTRP model, respectively, whereas the relative difference between the estimates β ML and β CLS for the respective models is rather small. However, let us note that the sum of squares SS CLS is somewhat smaller then SS ML .
As the results of Table 13 show, the real data of the first n = 40 failure times could be treated as a realization of the WPLP( α ML , β ML , γ ML ) with α ML = 0.049, β ML = 0.76011 and γ ML = 0.83068 or as a realization of the PTRP( α CLS , β CLS ) with α CLS = 0.02586 and β CLS = 0.8295. In both cases, the estimates of β are almost the same.
The percentages of the relative errors of the predicted next failure times are very small (Table 14). However, the CLS method applied to the data set considered forecasts the next failure time T n+1 more accurate: the value T CLS n+1 is closer to the real value T n+1 than the other predictions T ML n+1 and T ML n+1 . The lengths of the prediction intervals I CLS are shorter than the corresponding I ML ones.

Concluding remarks
According to our knowledge the problem of prediction in a TRP has not been engaged a good deal of the literature. In this paper we presented new methods of the prediction for the next failure time in a TRP(F, λ(·)) model in the case when its renewal distribution defined by F is completely unknown. On the basis of the remarks of Sects. 5.1-5.3 we recommend the CLS method for constructing the point and interval predictions in the case of the PTRP(α, β) model.
Besides the methods of interval prediction presented in Sect. 4 one can also consider the method associated with Proposition 4 simulating a large number of T (i) n+1 and taking the sample quantiles of order ε 1 and 1 − ε 2 of the simulated values.
It would be worth to examine the use of the CLS method in constructing the predictions also in the TRP models which are not the PTRP, i.e. which have other than power law type trend function. Considering the time truncation procedure would be also desirable.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.