Maximum Likelihood Estimation for an Inhomogeneous Gamma Process with a Log-linear Rate Function

An inhomogeneous gamma process is a compromise between a renewal process and a nonhomogeneous Poisson process, since its failure probability at a given time depends both on the age of the system and on the distance from the last failure time. The inhomogeneous gamma process with a log-linear rate function is often used in modelling of recurrent event data. In this paper, it is proved that the suitably non-uniform scaled maximum likelihood estimator of the three-dimensional parameter of this model is asymptotically normal, but it enjoys the curious property that the covariance matrix of the asymptotic distribution is singular. A simulation study is presented to illustrate the behaviour of the maximum likelihood estimators in finite samples. Obtained results are also applied to real data analysis.


Introduction
An inhomogeneous gamma process (IGP) was defined by Berman [5] in the following manner. Consider a Poisson process with intensity function (t) . Suppose that an event occurs at the origin, and that thereafter only every th event of the Poisson process is observed. Then, if T 1 , … , T n are the times of the first n events observed after the origin, their joint density is the following where {N(t), t ≥ 0} is the corresponding counting process and z is the hazard function of the gamma distribution G( , 1) with unit scale parameter and shape parameter . The hazard function z(t) of the G( , 1) distribution with < 1 decreases to 1 when t tends to infinity. If > 1 , then the function z(t) increases to 1 when t tends to infinity. Therefore, if the events are thought of as shock, then a value of > 1 indicates that the system is in better condition just after a repair than just before a failure and the larger is, the larger the improvement will be. A value of < 1 indicates that the system is in worse condition just after a repair than just before a failure. When = 1, the IGP reduces to the nonhomogeneous Poisson process. From the above interpretation, one can see that the IGP is an important process especially in the cases when the assumption of minimal repair which characterizes inhomogeneous Poisson process models is violated. Point and interval estimation of the parameter is thus an important task from a practical point of view because allows us to detect whether a system is in better or worse condition after a repair than just before a failure.
A point process {T i , i = 1, 2, … , n} with the joint density (1) for all positive integers n and all real positive is called the IGP with rate function (t) and shape parameter .
An alternative method of deriving the IGP is the following one. Suppose that the random variables i = 1, … , n, are independently and identically distributed according to the gamma G( , 1) distribution. It then follows that (1) is the joint distribution of T 1 , … , T n .
The IGP is a compromise between the renewal process (RP) and the nonhomogeneous Poisson process (NHPP), since its failure probability at a given time t depends both on the age t of the system and on the distance of t from the last failure time. Thus, it seems to be quite realistic model in many practical situations. The IGP for which lim t→∞ (t) = ∞ can be regarded as a special case of a trend renewal process (TRP) introduced and investigated first by Lindqvist [18] and by Lindqvist et al. [20] (see also [11,19,21,22]). The class of the TRP's is a rich family of processes and was considered in the field of reliability [21], finance [34, (1) (2) (t) = (t)z( (t) − (T N(t − ) )), 35], medicine [25], hydrology [9], software engineering [10,29], and to forecasts of volcanoes eruption [4]. The IGP, as a special case of the TRP, can be therefore a relevant model of recurrent events in the above mentioned fields.
In this paper, we consider the maximum likelihood (ML) estimation of the parameters of the IGP with the log-linear intensity function where T 0 = 0, between events tend to get smaller, and the larger is, the larger this trend will be. To give an interpretation of the parameter , let us notice that the intensity function (3) can be written in the following form When < 1 (and > 0 ), the summand (ln )∕ is less than 0, so the intensity increases slowly in the initial phase (for t ∈ (0, −(ln )∕ ) ). When > 1, the summand (ln )∕ is greater than 0, so the intensity function increases very fast from the beginning. Therefore, the summand (ln )∕ can be viewed as the parameter of time translation, and for < 1 we have a translation to the left, for > 1-to the right.
The IGP with the rate function (3) will be denoted by IGPL( , , ). Statistical inference for the IGP was considered by Berman [5] and for modulated Poisson process (a special case of IGP) by Cox [6]. Both papers only seriously addressed questions of hypothesis testing (via the likelihood ratio test), but did not satisfactorily solve the problem of parameter estimation. Inferential and testing procedures for log-linear nonhomogeneous Poisson process (a special case of the IGP considered in this paper) can be found in Ascher and Feingold [1], Lewis [17], MacLean [23], Cox and Lewis [7], Lawless [16] and Kutoyants [14,15]. In the paper of Bandyopadhyaya and Sen [3], the large-sample properties of the ML estimators of the parameters of IGP with power-law form of the intensity function are studied.
The article is organized as follows. In Sect. 2, the log-likelihood equations for the IGPL( , , ) are derived. In Sect. 3, asymptotic properties of ML estimators of the unknown parameters are given. As in the case of IGP with power-law intensity, considered by Bandyopadhyay and Sen [3], in the IGP with log-linear intensity the Hessian matrix of the log-likelihood function converges in probability to a singular matrix. Therefore, to prove the asymptotic normality of ML estimators in the model considered, we used an analogous method as in the paper of Bandyopadhyay and Sen [3]. In Sect. 4, we present the results of simulation study concerning the behaviour of the ML estimators of the model parameters in finite samples. We also illustrate the differences in behaviour of the ML estimator of the parameter compared to ML estimators of and . The asymptotic distribution of the ML estimators, derived in Sect. 3, we apply to obtain the realizations of the pointwise asymptotic confidence intervals for the unknown parameters of IGPL model in the real data analysis contained in Sect. 5. Section 6 contains conclusions and some prospects. Proofs of all theorems formulated in this paper are given in Sect. 7.

The ML Estimation in the IGPL Model
Let us notice that for the IGPL ( , , ) and Therefore, in contrast to the IGP with power-law intensity considered by Bandyopadhyay and Sen [3], the IGPL can be used to model reliability growth with bounded unknown number of failures. For example in the case < 0 and = 1 , the IGPL( , , ) is known as the Goel-Okumoto software reliability model (see [8]). Maximum likelihood estimation for the class of parametric nonhomogeneous Poisson processes (NHPP's) software reliability models with bounded mean value functions, which contains the Goel-Okumoto model as a special case, was considered by Zhao and Xie [33]. They showed that the ML estimators need not be consistent or asymptotically normal. They also derived asymptotic distribution for a specific NHPP model which is called the k-stage Erlangian NHPP software reliability model (see [13]) (for k = 1 this model is IGPL( , , 1) with < 0 ). Nayak et al. [24] extended the inconsistency results of Zhao and Xie [33] for all estimators of the unknown number of failures (not just the MLE), and for all NHPP models with bounded mean value functions.
From the above-mentioned results, one can see that properties of the ML estimators of IGPL parameters can depend on an assumed model (with decreasing or increasing rate function) and should be considered separately. We will consider the IGPL( , , ) for which > 0, > 0, > 0 . We suppose that the IGPL( , , ) is observed up to the nth event (failure) appears for the first time, and the values t 1 , … , t n of the jump times T 1 , … , T n are recorded. In other words, we consider the so-called failure truncation (or inverse sequential) procedure. It should be noted that the failure truncation procedure cannot be applied to IGPL( , , ) with < 0. Denote = (t 1 , … , t n ) . The likelihood function of the IGPL( , , ), observed until the nth failure occurs, is

Remark 1
The system of likelihood equations given above not always has a solution (̂,̂,̂) ∈ (0, ∞) 3 (see [12]). However, for some realizations of the IGPL, it has more than one solution.

Asymptotic Properties of ML Estimators
From now on, we denote vector of process parameters by = ( , , ) � , and 0 = ( 0 , 0 , 0 ) � indicates the true parameters values. We will use standard symbols o P (⋅) and O P (⋅) for convergence and boundedness in probability, respectively. All limits mentioned in this section will be taken as n → ∞ unless mentioned otherwise. Denote Now, define the scaled matrix obtained from the matrix A n ( ).

Theorem 1
The matrix C n ( 0 ) converges in probability to Note that C is a singular matrix with rank 2. Therefore, we cannot use standard methods (see for example [27,31,32]) to prove asymptotic normality of the ML estimators in the model considered. We proceed by reducing the problem to two-dimensions and appealing to the distributional properties of the IGPL. The parameter will be partitioned as follows � = ( , , ) = ( , � ). Substituting by n [exp( t n ) − 1] −1 into (6) and (7), the score functions n ( , , ; ) and n ( , , ; ) reduce to * Denote where and h are fixed numbers, 0 < < 1 2 , 0 < h < ∞.

Theorem 3 Vector n is asymptotically (singular) normal with mean vector zero and covariance matrix
Corollary 1 Theorem 3 can be applied to construct pointwise asymptotic confidence intervals for the parameters of the IGPL.

Remark 2
As in the case of the MPLP considered by Bandyopadhyay and Sen [3], the asymptotic result of Theorem 3 provides some curious insights into the behaviour of the MLE's of the IGPL parameters. Apart from the singularity and nonuniform scalings of the MLE'e, we have also that the estimator ̂ is asymptotically independent of the estimators ̂ and ̂.

Remark 3
The asymptotic result for the ML estimators of the parameters of nonhomogeneous Poisson process with log-linear intensity, namely, the process IGPL( , , 1) , can be obtained by substituting 0 = 1 in the top 2 × 2 top left submatrix of .

Simulation Study
In this section, we report a simulation study of the finite sample performance of the ML estimators of the IGPL model parameters. For each selected combination of ( 0 , 0 , 0 ) and n number of events, where n ∈ {25, 50, 75, 100} , Monte Carlo simulations with 1000 replications were performed. The IGPL( , , ) realizations (t 1 , … , t n ) were generated according to the formula where G( 0 ) is a random value generated from gamma G( 0 , 1) distribution and t 0 = 0 . For each realization (t 1 , … , t n ) , the MLE of 0 was calculated as a solution to the following equation where S n ( ) and W n ( ; ) are given by (4) and (8), respectively. A solution to Eq. (13) was obtained using Newton-Raphson method, implemented in nleqslv function inside R package nleqslv, with the initial value The initial value start was taken on the basis of an analogous reasoning to the construction of simple estimators, presented in the paper of Bandyopadhyay and Sen [3].
Estimates of and were determined by the formulas and respectively. Performance of the estimates is investigated in terms of bias (Bias) and mean square error (MSE), given by following formulas In Tables 1 and 2, we present the empirical biases and MSE's for = 1.5 and = 0.5 , respectively, for simulated data sets. We have taken values = 1.5 and = 0.5 to consider the case of the intensity function which increases very fast from the beginning ( > 1 implies time translation to the right) and the case of the intensity function which increases slower in the initial phase ( < 1 implies time translation to the left). From the results collected in these tables, we conclude that: • the empirical biases and MSE's of the estimators ̂ and ̂ are not so big even for n = 25 , but the estimator ̂ has rather big empirical MSE's, especially for small n; • the empirical MSE's of the estimator decrease a much slower rate than the empirical MSE's of the estimator ̂ and ̂ as n increases; • the empirical MSE's of the estimator ̂ for n = 100 are almost the same for various values of and when the value of is fixed.

Application to Some Real Data Set
In this section, we apply Theorem 3 to obtain realizations of pointwise asymptotic confidence intervals for the parameters of the IGPL model fitted to a real data set.
For each data set, we calculated the Akaike information criterion (AIC) and Bayesian information criterion (BIC) for five special cases of inhomogeneous gamma process model: the power-law process (PLP), the modulated power-law process (MPLP) considered by Bandyopadyay and Sen [3], the gamma renewal process (GRP), the nonhomogeneous Poisson process with log-linear intensity function (NHPPL), and the IGPL considered in this paper.

Diesel Engine
We consider the failure times (in thousands) in operating hours to unscheduled maintenance actions for the USS Halfbeak No.3 main propulsion diesel engine (see [2]). The data were considered by Rigdon [28], where the author assumed the power-law process model and obtained the ML estimates of the parameters. We assumed that the system was observed until the 71st failure at 25518 h. The values of AIC and BIC are given in Table 3.
The smallest values of the AIC and BIC are for the IGPL model, and therefore, the IGPL model is the best within the class of models considered, regardless of criterion. It is better than the PLP model considered by Rigdon [28].  The estimates (point and interval) of the IGPL parameters are given in Table 4. Realizations of 95% pointwise asymptotic confidence intervals are obtained using Theorem 3. In Table 4, the bootstrap confidence limits are also given for comparison.
The estimated value of is less than 1, what indicates that the system is in worse condition just after a repair than just before a failure.

Air Conditioning
As the second example, we consider the successive failures of the air conditioning system of Boeing 720 jet airplanes nr 7912, presented in work of Proschan [26]. The system was observed till 30th failure at 1788 hours. For numerical reasons, we consider event times in hundreds of hours. The values of AIC and BIC are given in Table 5.
According to AIC and BIC, the most appropriate model for air conditioning failure process is NHPPL. Let us notice that NHPPL is a special case of IGPL process with = 1.
The estimates, pointwise asymptotic confidence intervals and bootstrap confidence limits of the IGPL parameters are given in Table 6.
It can be observed that both (asymptotic and bootstrap) realizations of the confidence intervals for include 1, what suggests the correctness of the model choice based on the previously considered criteria.

Concluding Remarks
Asymptotic properties of ML estimators of the unknown parameter of the IGPL model were given. As in the case of IGP with power-law intensity, considered by Bandyopadhyay and Sen [3], in the IGPL the Hessian matrix of the log-likelihood function converges in probability to a singular matrix. Therefore, to prove the asymptotic normality of ML estimators in the model under study, a non-standard method has been applied. Moreover, the ML estimator enjoys the curious property that the covariance matrix of the asymptotic distribution is singular. The consistency of ML estimators in the IGPL model, as well as in the modulated power-law process considered by Bandyopadhyay and Sen [3], remains as the open problem.

Proof of Theorem 1
To prove Theorem 1, we shall first formulate and prove some lemmas. To simplify the notation of the proof, we define the following random variables and (14) The next lemma provides some necessary results concerning random variables Y i , i = 1, … , n.

Lemma 1
The random variables Y i , i = 1, 2, … , are such that

Proof
(i) Let us note that where i = 1, … , n, are independent gamma random variables with shape parameter 0 and scale parameter 1. Next by application of the Lindeberg-Feller central limit theorem, we obtain the result.
where Denote the first and second term of K n by K n1 and K n2 , , respectively. Using (18), we obtain that Setting a correspondence of e in with a Riemann sum, we observe that as n → ∞ , as n → ∞. These facts enable us to use the Lyapunov's central limit theorem. From the fact that K n2 converges to 0 as n → ∞, part (ii) of the lemma is proved. where as n → ∞. Using Markov inequality, we obtain that The sequence n = (U 1n , U 2n , U 3n ) � , n = 1, 2, … , converges in distribution to a multivariate normal random variable with mean vector zero and covariance matrix Proof From properties of the IGPL, X i = [exp( t i ) − exp( t i−1 )] , for i = 1, … , n, are independent gamma G( , 1) distributed random variables, and U 2n and U 3n can be re-expressed asỸ Then, by an application of bivariate central limit theorem (U 2n , U 3n ) � are asymptotically normal with zero mean vector and covariance matrix Random variable U 1 can be expressed in terms of Y i as Using Stirling's formula, the non-random term on RHS of above equation is o P (1) . From properties of the IGPL and Lemma 1 follows that U 1 is independent of (U 2 , U 3 ) � . Part (iv) of Lemma 1 also entails that U 1 converges in distribution to N(0, −1 ), which ends the proof. ◻

Lemma 3
The random variables T i , i = 1, 2, … , are such that Using the inequality x y ≥ x − 1 y − 1 for x < y and x, y > 1, we have Denote the lower bound of the above expression by B L . It is enough to show that B L = o P (1) . From properties of the IGPL, we know that .
are independent with B( (i − 1), ) distribution. Therefore, Denote the lower and upper bound of the above expression by B L and B U , , respectively. The bound B L by part (i) of this lemma is o P (1). Hence, it will be enough to show that B U = o P (1) . We have that ◻ The elements of C n ( ) = (c ij ), i, j = 1, 2, 3, can be written in the following form where and U 1n , U 2n , U 3n are given by (14), (15), (16), respectively. Now, the theorem follows from Lemma 3 and the fact that n = O p (1).

Proof
(i) We have The first factor of the above expression is O P (1). Using the Taylor series expand of exp function, we obtain Now, it is enough to show that 1∕T n is o P (1). For > 0 , we have Therefore, .
and the part (i) of the lemma is proved. (ii) We have From (i) the first factor of the above expression is o P (1), the second factor where the last convergence follows from the properties of the IGPL and the law of large numbers, and ends proof of (ii). ◻ Lemma 5 (i) The random variable converges in distribution to a bivariate normal random variable with mean vector zero and covariance matrix (ii) C * n ( 0 ) converges in probability to * , as n → ∞.

Proof of Theorem 3
Assuming that the equation * n ( ) = 0 has a solution ̂n = (̂n,̂n) � in the set M n ( 0 ) , we can expand * n (̂n) around 0 and obtain where n is a point on the line segment joining ̂n and 0 . Then, where * n = (Z 2n , Z 3n ) � , Z 2n and Z 3n are given by (10) and (11), respectively. Multiplying both sides of (25) by ( * ) −1 , we have For ∈ M n ( 0 ), the last equality follows from Lemma 5. This implies that * n is asymptotically normal with zero mean and covariance matrix ( * ) −1 . Furthermore, using the equality log b (a − c) = log b a + log b (1 − c a ), we have Using the asymptotic normality of U 2n and * n , and the equality it follows that Using the delta method, we get the expected dependence. ◻

Conflicts of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest. *