Parameter Estimation Problem for Doubly Geometric Process with the Gamma Distribution and Some Applications

1. Introduction

In the statistical literature, a data set with occurrence times of successive events generally can be modeled by using a counting process (CP). To determine a suitable stochastic CP model, it has to be tested whether the data set has a trend or not. Pekalp and Aydogdu [32] compared the monotonic trend tests for some counting processes. If the successive interarrival times are independent and identically distributed (iid) (there is no trend), the data set may be modeled by a renewal process (RP). However, in real-life examples, the successive inter-arrival times contain a monotone trend because of the aging effect and the accumulated wear [8]. In this case, this trend can be modeled by a nonhomogeneous Poisson process (NHPP) [2,3,4,5,6,7,8,9,10]. The data set having a monotone trend can also be analyzed by a geometric process (GP). GP is one of the widely used and known models for the monotone trend. Lam [17,18] first introduced the GP process. Furthermore, the process is used as a model in many areas in the reliability context. For details, see [36]. Lam [22] presented the GP theory and its applications. The real datasets with a monotone trend are modeled by GP in [39].

Although the GP is known as the most commonly used model, it has some limitations that can cause difficulties. The two limitations can be given as follows: Using the GP for non-monotone interarrival times with distributions with varying shape parameters is not suitable. The other limitation is that GP only allows logarithmic growth or explosive growth [7]. The GP model causes difficulties in the applications. Therefore, it can be said that the GP could be unsuitable for the mentioned cases.

To overcome the above difficulties, some stochastic models were developed by Wu and Scarf [36], Wu [37], and Wu and Wang [38]. One of the important models is the DGP for such models.

Wu [37] compares the DGP with the GP and exhibits the advantages and the preferability of the DGP. The definitions of CP, RP, GP, and DGP are presented as follows.

Definition 1.1 $N\left(t\right)=sup\left\{n\ge 0:S_n\le t\right\}$is the number of events that occurred in the interval$(0,t],$then$\left\{N\right(t),t\ge 0\}$is called a CP where,$S_0=0,S_n={\sum }_{k=1}^n{X}_{k}n=\mathrm{1,2},\dots$be the occurrence time of$thenth$event,$X_k=S_k-S_{k-1}$inter-arrival time be the time of$the(k-1)th$and$kth$event;

Definition 1.2.If$\{X_k,k=\mathrm{1,2},\dots \}$are (iid) random variables with cumulative distribution function (cdf )F, a CP$\left\{N\right(t),t\ge 0\}$corresponds to an RP.

Definition 1.3.Let’s assume that$\left\{N\right(t),t\ge 0\}$is a CP and$\{X_k,k=\mathrm{1,2},\dots \}$is the set of interarrival times of a CP$\left\{N\right(t),t\ge 0\}$. If there exists a positive constant value$a$defined as a ratio parameter such that$Y_k{=a}^{k-1}{X}_{k}k=\mathrm{1,2},$…,are iid form an RP with a common$F$is the cdf of the first inter-arrival time$X_1$ then the CP corresponds to a GP.

The cdf of $X_k$ is $F_k\left(x\right)=F\left(a^{k-1}x\right),k=\mathrm{1,2},\dots$.Therefore the distribution of $X1$ uniquely determines the distribution of $X_k$. Now, let $E\left(X1\right)=\mu$ and $Var\left(X1\right)={\sigma }^2$. Then the expected value and the variance of the random variable $X_k$ are given as $E\left(X_k\right)={\displaystyle \frac{\mu }{a^{k-1}}},Var\left(X_k\right)={\displaystyle \frac{{\sigma }^2}{a^{2\left(k-1\right)}}},k=\mathrm{1,2},\dots .$.Therefore it is clear that the cdf of $X_1$ uniquely determines the cdf of $X_k\text{'}s$ that is; $F_k\left(x\right)=F\left(a^{k-1}x\right)$, $k=\mathrm{1,2},3,\dots$..

Monotonicity has an important role in the theory of stochastic processes. If $a<(>)1$, then $\{X_k,k=\mathrm{1,2},\dots \}$ is defined as stochastically increasing (decreasing), and if $a=1$ then the GP corresponds to the RP.

Definition 1.4.Let’s assume that$\left\{N\left(t\right),t\ge 0\right\}$is a CP$\left\{X_k,k=\mathrm{1,2},3,\dots \right\}$is the set of interarrival times of this process. If there exists a positive constant$a$ratio parameter such that$Y_k{=a}^{k-1}{X_k}^{h\left(k\right)}k=\mathrm{1,2},$…are iid form an RP with a common F is the CDF of the$X_1$, the process$\left\{N\left(t\right),t\ge 0\right\}$is called a DGP, where$h\left(k\right)$is a positive function of$k$with$k=\mathrm{1,2},\dots h\left(1\right)=1E\left(X1\right)=\mu$and$Var\left(X1\right)={\sigma }^2$.

Wu [37] handled different $h\left(k\right)$ such as ${h\left(k\right)=(1+logk)}^b{,h\left(k\right)=b}^{k-1},{h\left(k\right)=b}^{log\left(k\right)}$, and $h\left(k\right)=1+blogk$ where $b$is a real number, $logk$ denotes the logarithm with base 10. He has shown that the DGP with $h\left(k\right)={\left(1+logk\right)}^b$ is better than the others for ten real data sets in Lam [23]. Then it is taken as $h\left(k\right)={\left(1+logk\right)}^b$ in this study.

Now we assume that the random variable $X_1$ is a continuous random variable with a pdf $f.$ Hence

$$f_{X_k}\left(x\right)=a^{k-1}h\left(k\right)x^{h\left(k\right)-1}f\left(a^{k-1}{x}^{h\left(k\right)}\right),x\ge 0$$

(1)

Moreover, the expected value and the variance of $X_k$ are

$$\begin{gathered} E\left(X_k\right)=a^{\left(1-k\right)h^{-1}\left(k\right)}{\int }_0^{\infty }{x}^{h^{-1}\left(k\right)}f\left(x\right)dx \\ Var\left(X_k\right)=a^{2\left(1-k\right)h^{-1}\left(k\right)}\left({\int }_0^{\infty }{x}^{{2h}^{-1}\left(k\right)}f\left(x\right)dx-{\left({\int }_0^{\infty }{x}^{h^{-1}\left(k\right)}f\left(x\right)dx\right)}^2\right). \end{gathered}$$

Provided the above integral are finite, where $h^{-1}\left(k\right)={\displaystyle \frac{1}{h\left(k\right)}}$.

Wu [37] gives the monotonicity properties of the DGP given as follows.

i) If $P\left(X_1>1\right)=1$ and $b<0$, $0<a<1$ or if $P\left(0<X_1<1\right)=1$ and $0<b<4.898226,0<a<1$, $\{X_k,k=\mathrm{1,2},\dots \}$ increases stochastically.

ii) If $P\left(0<X_1<1\right)=1$ and $b<0,a>1$ or if$P\left(X_1>1\right)=1$ and $0<b<4.898226,a>1$, $\{X_k,k=\mathrm{1,2},3\dots \}$ decreases stochastically.

The parameter estimation problem naturally arises in DGP. The parameter estimation problem for a DGP contains the parameters $a,b,\mu ,$ and${\sigma }^2.$ These parameters determine the mean and the variance of the first inter-arrival time $X_1$. Therefore, the parameter estimation problem for DGP is a very important issue as well as in GP. The some of the studies in GP are given as follows.

Lam and Chan [20], Chan et al. [8], Aydoğdu et al. [3], and Kara et al. [13] used the lognormal, gamma, Weibull, and inverse Gaussian distributions, respectively, for the $X_1$ interarrival time to estimate the parameters for a GP. Kara et al. [16] consider the parameter estimation problem for the gamma geometric process.

Pekalp and Aydogdu [30] considered the power series expansions for the probability distribution, mean value, and variance function of a geometric process with gamma interarrival times. Pekalp and Aydogdu [33] considered the parameter estimation problem for the mean value and variance functions in GP.

Aydogdu and Altındag [5] computed the mean value and variance functions in a geometric process. Altındag [1] evaluated the multiple process data in a geometric process with exponential failures.

Yılmaz et al. [40] used the Bayesian inference with Lindley distribution by a GP.

Pekalp and Aydogdu [26] studied the integral equation for the second moment function in a GP. Pekalp and Aydogdu [27], Pekalp et al. [28] obtained the asymptotic solution of the integral equation for the second moment function of a GP by discriminating the lifetime distributions of the ten real data sets used in [21].

Pekalp et.al [29,31,34] estimate the parameters of a DGP by using the ML method under the exponential, Weibull, and lognormal distribution assumptions, respectively, for the first inter-arrival time $X_1$. Eroglu Inan [11] considered the parameter estimation problem for a DGP under the assumption that the first interarrival time has an inverse Gaussian distribution.

Although there are some studies in GP for the gamma distribution, which is used and is an important model in reliability theory, to the best of our knowledge, no study has been done yet about DGP, which removes the lack of GP model. Therefore, the parametric statistical inference problem for DGP under the gamma distribution assumption can be studied.

The nonparametric problem was studied for GP and $\alpha -$ series processes by the some authors [4,15,19,35,37]. Similarly, Jasim and Al-Qazaz [12] proposed inference methods for DGP.

In this study, the parameter estimation problem for a DGP is considered under the assumption that the first inter-arrival time $X_1$ has a gamma distribution. This paper is organized as follows: In Section 2, the gamma distribution and its probabilistic properties are reminded. In Section 3, the log-likelihood function, first derivatives of this function, and the asymptotic joint distribution of the estimators are obtained. In Section 4, an extensive simulation study is performed, and the simulated mean, bias and mean square error (MSE) values are calculated, and the small sample performances of the estimators are evaluated for various values of parameters. In Section 5, the importance of the modelling with DGP is presented by the two illustrative examples. Finally, in Section 6, the results are discussed.

2. Gamma Distribution

The gamma distribution is one of the most commonly used models in the reliability and life testing areas for asymmetric data, for details, [9,23] referred to in [14].

The probability density function (pdf) of the gamma distribution is given by

$$f_X\left(x,\alpha ,\beta \right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right){\beta }^{\alpha }}}{x}^{\alpha -1}{e}^{-{\displaystyle \frac{x}{\beta }}}x>0;\alpha >0,\beta >0$$

(2)

where

$\alpha :$ shape parameter

$\beta$: scale parameter

If $\alpha =1$, the gamma distribution corresponds to an exponential distribution.

Let’s $X$ have a gamma distribution with $\alpha$ and $\beta$ parameters, that is $X~gamma\left(\alpha ,\beta \right)$. Some characteristics of the gamma distribution such as expected value, variance, skewness, and kurtosis are given as $E\left(X\right)=\alpha \beta ,Var\left(X\right)=\alpha {\beta }^2,{\gamma }_1={\displaystyle \frac{2}{\sqrt{\alpha }}}$ and ${\gamma }_2={\displaystyle \frac{6}{\alpha }}$ respectively.

3. Statistical Inference for Gamma Distribution

Let’s assume that {N(t), t ≥ 0} is a DGP with the ratio parameter $a$, $h\left(k\right)={\left(1+logk\right)}^b,$ and the first inter-arrival time distribution $X_1$ has a $gamma\left(\alpha ,\beta \right).$For a realization $\{X_1,X_2,\dots ,X_n\}$of the DGP, the random variables $Y_k=a^{k-1}{X}_k^{h\left(k\right)}$ are iid with the $gamma\left(\alpha ,\beta \right)$ distribution. From the equation (1) we have that the pdf of $X_k$ is

$$f_{X_k}\left(x\right)=a^{k-1}{x}^{{\left(1+logk\right)}^b-1}\mathrm{br}\mathrm{to}\mathrm{break}\times {\left(1+logk\right)}^b{\displaystyle \frac{1}{\Gamma \left(\alpha \right){\beta }^{\alpha }}}{\left(a^{k-1}{x}^{{\left(1+logk\right)}^b}\right)}^{\alpha -1}{e}^{-{\displaystyle \frac{\left(a^{k-1}{x}^{{\left(1+logk\right)}^b}\right)}{\beta }}}$$

(3)

Further, the r. moment of $X_k$ is easily seen that

$$E\left({X_k}^r\right)={\left(1+logk\right)}^b{\beta }^{r{\left(1+logk\right)}^{-b}}{a}^{\left(1-k\right)\left[1+r{\left(1+logk\right)}^{-b}\right]}\mathrm{br}-\mathrm{to}-\mathrm{break}\times {\displaystyle \frac{\Gamma \left(\alpha +{\left(1+logk\right)}^{-b}\right)}{\Gamma \left(\alpha \right)}}r=\mathrm{1,2},\dots .$$

(4)

$Var\left(X_k\right)={\beta }^{2{\left(1+logk\right)}^{-b}}\left({\displaystyle \frac{\Gamma \left(\alpha +{\left(1+logk\right)}^{-b}\right)}{\Gamma \left(\alpha \right)}}\right){\left(1+logk\right)}^b$$\times \left[a^{2\left(1-k\right)\left[1+{\left(1+logk\right)}^{-b}\right]}\left(a^{1-k}-\left({\displaystyle \frac{\Gamma \left(\alpha +{\left(1+logk\right)}^{-b}\right)}{\Gamma \left(\alpha \right)}}\right){\left(1+logk\right)}^b\right)\right]$The likelihood function $L\left(a,b,\alpha ,\beta \right)$ based on the realization $\left\{X_1,\dots ,X_n\right\}$

$$L\left(a,b,\alpha ,\beta \right)={\left({\displaystyle \frac{1}{\Gamma \left(\alpha \right){\beta }^{\alpha }}}\right)}^n{a}^{{\displaystyle \frac{\alpha n\left(n-1\right)}{2}}}{\prod }_{k=1}^n{\left(1+logk\right)}^b\mathrm{br}\mathrm{to}\mathrm{break}\times \prod _{k=1}^n{x_k}^{\alpha {\left(1+logk\right)}^b-1}{e}^{-{\sum }_{k=1}^n\left({\displaystyle \frac{a^{k-1}{x_k}^{{\left(1+logk\right)}^b}}{\beta }}\right)}$$

(5)

Then, log-likelihood function is

$$lnL\left(a,b,\alpha ,\beta \right)=-n\left(\mathit{log}\Gamma \left(\alpha \right)+\alpha \mathit{log}\beta \right)+{\displaystyle \frac{\alpha n\left(n-1\right)}{2}}lna\mathrm{br}\mathrm{to}\mathrm{break}+{\sum }_{k=1}^{n}ln\left({\left(1+logk\right)}^b\right)+{\sum }_{k=1}^n\left(\alpha {\left(1+logk\right)}^b-1\right)ln\left(x_k\right)\mathrm{br}-\mathrm{to}-\mathrm{break}-{\sum }_{k=1}^n\left({\displaystyle \frac{a^{k-1}{x_k}^{{\left(1+logk\right)}^b}}{\beta }}\right)$$

(6)

Firstly, the first partial derivatives of the log-likelihood function defined as Equation (6) are taken with respect to the parameters $a,b,\alpha ,\beta$ respectively, and then the likelihood equations are obtained as follows by equating the derivatives to zero.

$$\begin{gathered} {\displaystyle \frac{\partial lnL}{\partial a}}=\alpha {\displaystyle \frac{n\left(n-1\right)}{2a}}-{\displaystyle \frac{1}{\beta }}\sum _{k=1}^n\left(k-1\right)a^{k-2}{x_k}^{{\left(1+logk\right)}^b}=0 \\ {\displaystyle \frac{\partial lnL}{\partial b}}={\sum }_{k=1}^{n}ln\left(1+logk\right)+{\sum }_{k=1}^n\alpha ln\left(x_k\right){\left(1+logk\right)}^{b}ln\left(1+logk\right) \\ -{\displaystyle \frac{1}{\beta }}\sum _{k=1}^n{a}^{k-1}{x_k}^{{\left(1+logk\right)}^b}\mathit{ln}\left(x_k\right){\left(1+logk\right)}^{b}ln\left(1+logk\right)=0 \\ {\displaystyle \frac{\partial lnL}{\partial \alpha }}={\displaystyle \frac{n\left(n-1\right)}{2}}lna-n\left({\displaystyle \frac{\partial \Gamma \left(\alpha \right)}{\partial \alpha }}\right)-nln\left(\beta \right)+\sum _{k=1}^n{\left(1+logk\right)}^b\mathit{ln}\left(x_k\right)=0 \\ {\displaystyle \frac{\partial lnL}{\partial \beta }}=-{\displaystyle \frac{n\alpha }{\beta }}+\sum _{k=1}^n{\displaystyle \frac{a^{k-1}{x_k}^{{\left(1+logk\right)}^b}}{{\beta }^2}}=0 \end{gathered}$$

where,$log$ and $ln$ denotes the logarithm with base 10 and $e$ , respectively.

It is seen that these likelihood equations can’t be solved analytically. Therefore, the values of the ML estimators of the $a,b,\alpha$ and $\beta$ must be calculated numerically. For this purpose, the estimates are obtained numerically with the Nmaximize subroutine, which is a constrained nonlinear optimization technique of the Mathematica package program.

The Fisher information matrix $I$ and its inverse $I^{-1}$ are derived to find the asymptotic joint distribution of the ML estimators. It is well known that the diagonal elements give the asymptotic variance for each estimator, respectively.

For the information matrix $I=\left(I_{ij}\right)$

$$\begin{gathered} I_{11}=-E\left({\displaystyle \frac{{\partial }^{2}lnL}{\partial a^2}}\right)={\displaystyle \frac{n\left(n-1\right)\alpha }{2a^2}}+{\displaystyle \frac{\alpha }{a^2}}{\sum }_{k=1}^n\left(k^2-3k+2\right) \\ {I}_{12}=-E\left({\displaystyle \frac{\partial lnL}{\partial a}}{\displaystyle \frac{\partial lnL}{\partial b}}\right) \\ ={\sum }_{k=1}^n{\displaystyle \frac{\left(\mathit{ln}\left(1+logk\right)\right)\left(k-1\right)}{a\beta }}\left[E\left(Y_k\mathit{ln}\left(Y_k\right)\right)-\left(k-1\right)lna\alpha \beta \right] \\ {I}_{13}=-E\left({\displaystyle \frac{\partial lnL}{\partial a}}{\displaystyle \frac{\partial lnL}{\partial \alpha }}\right)=-\left({\displaystyle \frac{n\left(n-1\right)}{2a}}\right) \\ {I}_{14}=-E\left({\displaystyle \frac{\partial lnL}{\partial a}}{\displaystyle \frac{\partial lnL}{\partial \beta }}\right)=-\left({\displaystyle \frac{\left(n^2-n\right)\alpha }{2a\beta }}\right) \\ {I}_{22}=-E\left({\displaystyle \frac{{\partial }^{2}lnL}{\partial b^2}}\right)=-\left[\alpha {\left(\mathit{ln}\left(1+logk\right)\right)}^2\left[E\left(\mathit{ln}\left(Y_k\right)\right)-\left(k-1\right)lna\right]\right] \\ +\left[\left({\displaystyle \frac{1}{\beta }}\right){\sum }_{k=1}^n\left(\left(E\left(\left(Y_k\right)ln\left(Y_k\right)\right)-\left(k-1\right)lna\alpha \beta \right){\left(\mathit{ln}\left(1+logk\right)\right)}^2\right)\right] \\ +\left[{\displaystyle \frac{{\left(\mathit{ln}\left(1+logk\right)\right)}^2}{\beta }}\left({\sum }_{k=1}^{n}E\left(\left(Y_k\right){\left(ln\left(Y_k\right)\right)}^2\right)+\alpha \beta {\left(\left(k-1\right)lna\right)}^2\right.\right. \\ \left.-2E\left(\left(Y_k\right)ln\left(Y_k\right)\right)\left(k-1\right)lna\right] \\ {I}_{23}=-E\left({\displaystyle \frac{\partial lnL}{\partial b}}{\displaystyle \frac{\partial lnL}{\partial \alpha }}\right)={\sum }_{k=1}^n\left(\mathit{ln}\left(1+logk\right)\right)\left(E\left(\mathit{ln}\left(Y_k\right)\right)-\left(k-1\right)lna\right) \\ {I}_{24}=-E\left({\displaystyle \frac{\partial lnL}{\partial b}}{\displaystyle \frac{\partial lnL}{\partial \beta }}\right) \\ =-{\displaystyle \frac{1}{{\beta }^2}}{\sum }_{k=1}^n\left(E\left(\left(Y_k\right)ln\left(Y_k\right)\right)-\left(k-1\right)lna\alpha \beta \right)\left(\mathit{ln}\left(1+logk\right)\right) \\ {I}_{33}=-E\left({\displaystyle \frac{{\partial }^{2}lnL}{\partial {\alpha }^2}}\right)=npolygamma\left[1,\alpha -1\right] \\ {I}_{34}=E\left({\displaystyle \frac{{\partial }^2\ln L}{\partial \alpha \partial \beta }}\right)={\displaystyle \frac{n}{\beta }} \\ {I}_{44}=-E\left({\displaystyle \frac{{\partial }^2\ln L}{\partial {\beta }^2}}\right)={\displaystyle \frac{n\alpha }{{\beta }^2}} \end{gathered}$$

where

$polygamma\left[n,x\right]$: $(n+1).$ derivative of log-gamma function at$x$.

Then, by [6] the joint distribution of these ML estimators is asymptotically normal with the mean vector $\left(a,b,\alpha ,\beta \right)$ and the covariance matrix $I^{-1}$ .

$$\left(\begin{array}{c}\widehat{a}\\ \widehat{b}\\ \widehat{\alpha }\\ \widehat{\beta }\end{array}\right)~AN\left(\left(\begin{array}{c}a\\ b\\ \alpha \\ \beta \end{array}\right),I^{-1}\right)$$

(7)

Although the elements of $I^{-1}$ are obtained analytically, the expressions have very large size, so they can’t be given here. Indeed, by the notation $I^{-1}={\left(v_{ij}\right)}_{n\times n}$ it can be written that

$\widehat{a}~AN\left(a,v_{11}\right)$$\widehat{b}~AN\left(b,v_{22}\right)$$\widehat{\alpha }~AN\left(\alpha ,v_{33}\right)$$\widehat{\beta }~AN\left(\beta ,v_{44}\right)$From the theory of ML estimators, it is well known that these estimators are both asymptotically unbiased and consistent.

4. Simulation Analysis

In this section, a simulation study is done to evaluate the performance of the ML estimators. Mathematica is used for the computations. The simulated means, biases and MSE values are calculated for different sample sizes $n=30,50,100$ which are commonly used values in applications. The trend is usually small for $a$ ratio parameter values which are close to 1 [22]. Therefore, $a$ is taken as 0.95, 0.99, 1.01, and 1.05 throughout the study. The monotonic property of the DGP is altered by the parameter $b$’s positive and negative values. The possible effect of positive and negative values of the $b$ is considered by taking -2, 2. The parameter$\alpha$ and$\beta$ values are taken to be 1 and 2. The simulation is evaluated for 1000 repetitions, step by step:

1. Generate a sample$\{Y_1,Y_2,\dots ,Y_n\}$

from a gamma distribution with the parameters $\alpha$and $\beta$.

2. Calculate a realization$\left\{X_1,X_2,\dots ,X_n\right\}$ of the DGP with the parameters $a,b$ by using the equation $X_k={\left({\displaystyle \frac{Y_k}{a^{k-1}}}\right)}^{1/{(1+\mathrm{l}\mathrm{o}\mathrm{g}k)}^b},k=\mathrm{1,2},3\dots ,n$.

3. Calculate the mean, bias, and MSE values for each estimator based on the $\left\{X_1,X_2,\dots ,X_n\right\}$ data set. Table 1, Table 2, Table 3 and Table 4 include the corresponding simulated mean, biases, and MSE values for each estimator.

From Table 1, Table 2, Table 3 and Table 4 the MSE gets small values, and the absolute bias values get closer to zero for increasing $n$ values. This expected case indicates the fact that the estimators $a,b,\alpha$and

$\beta$are asymptotically unbiased and consistent. The diagonal elements of the $I^{-1}$give the lower bounds of the variances MVB of the parameters$a,b,\alpha$ and $\beta$. The simulated variances for the ML estimates and the MVB values are given in Table 5, Table 6, Table 7, Table 8 for all cases. When the MVB values compare the corresponding simulated variance values, and it is seen that the values get close to each other for increasing $n$ values, it can be said that the estimators are highly efficient estimators.

5. Data Analysis

In this section, two real data sets of are presented to illustrate the data analysis and estimation procedure given in Section 3. Further the performance of the ML estimators are compared with the nonparametric MM estimators given by Jasim and Al-Qazaz [12]. These data sets are called the coal mining disasters data and the propulsion diesel engine failure data.

Now we consider a data set $\left\{X_k,k=\mathrm{1,2},\dots \right\}$ comes from a DGP. Jasim and Al-Qazaz [12] obtained the least square (LS) estimates of the parameters $a,b$ and $\mu$ by minimizing the sum of squared error (SSE);

$SSE={\sum }_{k=1}^n{\left[\left(X_k\right)-{\left({\displaystyle \frac{\mu }{a^{k-1}}}\right)}^{{\left(1+logk\right)}^{-b}}\right]}^2$

By taking ${\widehat{a}}_{LSE}$and ${\widehat{b}}_{LSE}$for the parameters $a$ and $b$ in the equation $Y_k=a^{k-1}{X}_k^{{\left(1+logk\right)}^b}$ the prediction of $Y_k$ is

$${\widehat{Y}}_k={\left({\widehat{a}}_{MM}\right)}^{k-1}{X}_k^{{\left(1+logk\right)}^{{\widehat{b}}_{MM}}}$$

(8)

Then the MM estimators of the parameters $\mu ,{\sigma }^2$ are defined by

$${\widehat{\mu }}_{MM}=\overline{{\widehat{Y}}_k}={\displaystyle \frac{1}{n}}{\sum }_{k=1}^n{\widehat{Y}}_k$$

(9)

$${{\widehat{\sigma }}^2}_{MM}=\left\{\begin{array}{c}{\displaystyle \frac{1}{n-1}}{\sum }_{k=1}^n{\left({\widehat{Y}}_k-\overline{{\widehat{Y}}_k}\right)}^2,a\ne 1,h\left(k\right)\ne 1,b\ne 0\\ {\displaystyle \frac{1}{n-1}}{\sum }_{k=1}^n{\left(X_k-\overline{X}\right)}^2,a=1,h\left(k\right)=1,b=0\end{array}\right.$$

(10)

where

$\overline{X}={\displaystyle \frac{1}{n}}{\sum }_{k=1}^n{X}_k$

Note that ${\widehat{a}}_{LSE},$and ${\widehat{b}}_{LSE}$ are called MM estimators by Jasim and Al-Qazaz [12]. The estimators are denoted by ${\widehat{a}}_{MM}$ and ${\widehat{b}}_{MM}$ in respectively.

In the DGP model with the gamma distribution, the prediction of $Y_k$ is given

$${\widehat{Y}}_k={\left({\widehat{a}}_{ML}\right)}^{k-1}{X}_k^{{\left(1+logk\right)}^{{\widehat{b}}_{ML}}}$$

(11)

If it can be showed that the ${\widehat{Y}}_k^{\text{'}}s$ are iid with gamma distribution, the DGP with gamma distribution can be used to model the data sets. The goodness of fitness of the ${\widehat{Y}}_k^{\text{'}}s$ for gamma distribution with parameters $\alpha ,\beta$ are measured by the Kolmogorov-Simirnov (KS) test. The KS test statistic and corresponding $p-$value are given for the data sets in Tables 14,16.

To compare the ML and MM methods, the prediction of $X_k$ for the data set $\left\{X_1,\dots ,X_n\right\}$ is calculated as

$${\widehat{X}}_k=\left\{\begin{array}{c}{\left({{\widehat{a}}_{ML}}^{\left(1-k\right)}\right)}^{{\left(1+logk\right)}^{-{\widehat{b}}_{ML}}}{\left[\widehat{E}\left(X_1\right)\right]}^{{\left(1+logk\right)}^{-{\widehat{b}}_{ML}}},forMLmethod\\ {\left({\displaystyle \frac{{\widehat{\mu }}_{MM}}{{{\widehat{a}}_{MM}}^{k-1}}}\right)}^{{\left(1+logk\right)}^{-{\widehat{b}}_{MM}}},forMMmethod\end{array}\right.$$

(12)

where$\widehat{E}\left(X_1\right)$ is obtained by the Equation 4 as

$\widehat{E}\left(X_1\right)={\left(1+logk\right)}^{\widehat{b}}{\widehat{\beta }}^{{\left(1+logk\right)}^{\widehat{b}}}{\widehat{a}}^{\left(1-k\right)\left[1+{\left(1+logk\right)}^{\widehat{b}}\right]}$$\times {\displaystyle \frac{\Gamma \left(\widehat{\alpha }+{\left(1+logk\right)}^{-\widehat{b}}\right)}{\Gamma \left(\widehat{\alpha }\right)}}$The mean-squared error (${\mathrm{M}\mathrm{S}\mathrm{E}}^{*}$) and maximum percentage error (MPE) criterias are used to evaluate the performance of the DGPs with the ML and MM estimators. ${\mathrm{M}\mathrm{S}\mathrm{E}}^{*}$ and MPE are defined by Lam et al. [21] as

$${MSE}^{*}={\displaystyle \frac{1}{n}}{\sum }_{k=1}^n{\left(X_k-{\widehat{X}}_k\right)}^2$$

(13)

and

$$MPE=\max \left\{{\displaystyle \frac{\left|S_k-{\widehat{S}}_k\right|}{S_k}},k=\mathrm{1,2},\dots ,n\right\}$$

(14)

where

$$S_k=X_1+\dots +X_k,k=\mathrm{1,2},\dots ,n$$

(15)

and

$${\widehat{S}}_k={\sum }_{k=1}^n{\widehat{X}}_k$$

(16)

* The notation MSE* is different from the MSE defined in Section3.

5.1. First Data Set

The first data set, originally studied by Maguire et al. [25] (Data Set No.1). It has 190 observations show the interarrival time between coal mining disasters.

When the data set 1 is modeled by a DGP with the gamma distribution, the ML estimator values of the $a,b,\alpha ,\beta$ parameters, the MM estimates of the parameters $a,b,$MSE, and MPE values, and the estimates of $\mu ,{\sigma }^2$are given in Table 9.

It can be seen from the Table 9 that ML estimators have the smallest MSE and MPE values. The results show that the ML estimators outperform the MM estimators with smaller MSE and MPE values. For this data, the KS test statistic value and the corresponding p value are given in Table 10.

It is seen in Table 10 that the gamma distribution is appropriate for the data set 1. In Figure 1, $S_k$ and ${\widehat{S}}_k$ disaster times are plotted for both ML and MM methods.

The performance of the ML and MM methods can be compared by Figure 1. It is easily seen that the DGP with the ML estimator outperform the DGP with MM estimator. The results are compatible with Table 9. It is clear that the ML method works well as it is expected.

5.2. Second Data Set

The second data set is called the U.S.S. Grampus No. 4 main propulsion diesel engine failure data. (Data set No.2) contains the cumulative operating hours until significant maintenance events occur for one of the thirty engines on nine different submarines. It has 57 observations. The data set was studied by Lee [24].

When the data set 2 is modeled by a DGP with the gamma distribution, the ML estimator values of the $a,b,\alpha ,\beta$ parameters, the MM estimates of the parameters $a,b,$MSE, and MPE values, and the estimates of $\widehat{\mu },{\widehat{\sigma }}^2$ are given in Table 11.

It can be seen from Table 11 that ML estimators have the smallest MSE and MPE values. The results show that the ML estimators outperform the MM estimators with smaller MSE and MPE values. For this data, the KS test statistic value and the corresponding p value are given in Table 12.

It is seen in Table 12 that the gamma distribution is appropriate for the data set 2. In Figure 2, $S_k$ and ${\widehat{S}}_k$ disaster times are plotted for both ML and MM methods.

The performance of the ML and MM methods can be compared by Figure 2. It is easily seen that the DGP with the ML estimator outperform the DGP with MM estimator. The results are compatible with Table 11. It is clear that the ML method works well as it is expected.

6. Conclusion

In this study, the parameter estimation problem for a DGP is considered under the assumption that the first inter-arrival time has a gamma distribution. The estimators of the parameters $a,b,\alpha$ and $\beta$ are calculated by the ML method. The asymptotic joint distribution of the ML estimators are obtained. The unbiasedness and consistency statistical properties of the estimators are investigated. The Monte Carlo simulation study was performed to evaluate the small sample performance of the estimators. It is seen that the bias and the MSE values of the estimators decrease for increasing $n$ sample size values. This expected case supports the fact that the ML estimators are asymptotically unbiased and consistent. Finally, two data sets are considered in applications. The ML and MM estimators’ values, MSE, MPE and the KS values are calculated. As a result, it can be said that the gamma distribution is an appropriate model for these data sets. The DGP with ML estimators outperforms the DGP with MM estimators.