A New Three-Parameter Discrete Distribution with Associated INAR(1) Process and Applications

The aim of this article is to propose a new three-parameter discrete Lindley distribution. A wide range of its structural properties are investigated. This includes the shape of the probability mass function, hazard rate function, moments, skewness, kurtosis, index of dispersion, mean residual life, mean past life and stress-strength reliability. These properties are expressed in explicit forms. The maximum likelihood approach is used to estimate the model parameters. A detailed simulation study is carried out to examine the bias and mean square error of the estimators. Using the proposed distribution, a new first-order integer-valued autoregressive process is introduced for the over-dispersed, equi-dispersed and under-dispersed time series of counts. To demonstrate the importance of the proposed distribution, three data sets on coronavirus, length of stay at psychiatric ward and monthly counts of larceny calls are analyzed.


I. INTRODUCTION
Statistical distributions play an important role in data modeling, inference, and forecasting processes. The occurrence times, frequencies and effects of many events in nature are analyzed by statistical modeling techniques. Most of the events in nature or other scientific fields have their own characteristics. Earthquakes, traffic accidents, counts of landslide or number of people dying from the disease are modeled by discrete probability distributions. Researchers have proposed more flexible distributions to reduce estimation errors in the modeling of these data sets. There are two popular methods used to introduce a new discrete distribution. These are mixed-Poisson type discrete distributions and survival discretization method. The recently introduced discrete distribution based on the survival discretization method can be cited as follows: discrete Lindley (DLi) distribution by Gómez-Déniz and Calderín-Ojeda (2011), discrete inverse Weibull (DIW) distribution by Jazi et al. (2010), discrete Burr type XII (DB-XII) distribution by Para and Jan (2014), discrete Pareto (DPa) distribution by Krishna and Pundir (2009) The associate editor coordinating the review of this manuscript and approving it for publication was Giambattista Gruosso . generalized geometric by Gómez-Déniz (2010), discrete generalized exponential type II by Nekoukhou et al. (2013), discrete Rayleigh (DR) by Roy (2004), two-parameter discrete Lindley by Hussain Das et al. (2018) and Poisson-generalized Lindley distribution by Wongrin and Bodhisuwan (2016). Using the similar approach to Sankaran (1970), several authors have been introduced mixed-Poisson distributions such as Poisson-Bilal distribution by Altun (2020a), a new Poisson-weighted exponential distribution by Altun (2020b), Poisson-xgamma distribution , Poisson-weighted Lindley distribution by Atikankul et al. (2020) and Poissontransmuted exponential distribution by Bhati et al. (2017).
Additionally, the modeling of the time series of counts is an important research area for all applied sciences. For instance, a company in insurance sector cares about predicting the number of claims for next month. Another example is that epidemiologist wants to predict the monthly deaths from a disaster such as coronavirus, bird flu, ebola virus infection. In this case, first-order integer-valued autoregressive process, shortly INAR(1) process, can be used. The The goal of this work is to introduce an alternative discrete distribution to model both over-dispersed and under-dispersed count data sets. The over-dispersion and under-dispersion is a widely studied problem of count data modeling. The over-dispersion appears in the case that the empirical variance is greater than empirical mean. The opposite indicates the under-dispersion. In the real-life data modeling, we encounter quite often with these two problems, however, the over-dispersion problem is seen more. We introduce a new three parameter discrete distribution by using the three-parameter Lindley (Li-3P) distribution introduced by Shanker et al. (2017) based on the survival discretization method. The proposed distribution is called three parameter discrete-Lindley, shortly (DLi-3P), distribution. The contributions of the presented study to statistics literature can be summarized as follows: (i) an alternative model for the over(under)-dispersed data sets is introduced, (ii) the statistical properties of the DLi-3P distributions is studied in detail, (iii) INAR(1) process with DLi-3P innovations is introduced, (iv) three applications to coronavirus, psychiatric ward and larceny calls data sets are analyzed by the models introduced based on the DLi-3P model.
The article is organized as follows. In Section 2, we introduce the DLi-3P distribution. Different statistical and reliability properties are discussed in Section 3. In Section 4, the model parameters are estimated by using the maximum likelihood estimation (MLE) approach. In Section 5, we discuss the computational complexity of the DLi-3P distribution its limitations. In Section 6, INAR(1) process with DLi-3P innovations is introduced and its statistical properties are derived. Simulation study is presented in Section 7. Three distinctive data sets are analyzed to show the importance of the DLi-3P distribution in Section 8. The detail interpretation of the empirical results are given in Section 9. Finally, Section 10 provides some conclusions.

II. THE DLi-3P DISTRIBUTION
In this section, we derive the discrete analogous of the Li-3P distribution, shortly DLi-3P distribution, by using the survival discretization method. Assume that the random variable X follows a Li-3P distribution whose probability density function (pdf) and the corresponding survival function (sf) are given, respectively, by where θ > 0, β > 0 and αθ + β > 0. The derivation of Li-3P distribution is similar to the Lindley distribution. The Li-3P distribution is obtained as mixture distribution of exponential (θ ) and gamma (2, θ) with a mixing proportion αθ αθ +β . Using the survival discretization method and survival function of the Li-3P distribution given in (2), the probability mass function (pmf) of the DLi-3P distribution with positive parameter 0 < λ < 1 can be expressed as where λ = e −θ and N 0 = {0, 1, 2, 3, . . . , q} for 0 < q < ∞. The corresponding cdf and sf to (3) are given, respectively, by The pmf in (3) is log-concave, where P x (x+1;λ,α,β) P x (x;λ,α,β) is a decreasing function in x for all values of the model parameters. The several possible shapes of the DLi-3P distribution are displayed in Figure 1. From this figure, we conclude that the DLi-3P distribution could be used to model left-skewed count data sets.
The shapes of the pmf of DLi-3P can be unimodal or decreasing. The hazard rate function (hrf) of the DLi-3P distribution is  where x ∈ N 0 . Figure 2 displays the hrf plots of the DLi-3P distribution for different values of the model parameters. It is observed that the hrf of the DLi-3P distribution has increasing shape.

III. STATISTICAL PROPERTIES
In this section, the statistical properties of the DLi-3P distribution are derived such as mode, raw moments, skewness and kurtosis measures. Additionally, reliability properties of the DLi-3P distribution are derived such as stress-strength, mean residual life (mrl) and mean past life (mpl).

A. MODE
The mode of any discrete distribution shows the value at which the specific discrete distribution takes its maximum value. If X has a DLi-3P distribution, then the mode can be obtained by solving the following non-linear equation Then, the mode of the DLi-3P distribution can be expressed as

B. MOMENTS, SKEWNESS, KURTOSIS AND DISPERSION INDEX
Assume X be a DLi-3P random variable. Then, the probability generating function (pgf) can be expressed in a closed form as where G X (s) = ∞ x=0 s x P x (x; λ, α, β). Replacing s by e s in (9), we get the moment generating function (mgf). The first four derivatives of the mgf, with respect to s at s = 0, give the first four raw moments. Thus, the first four moments of the DLi-3P model are and The variance and dispersion index (DI) of the DLi-3P distribution are given, respectively, by (14) and (15), as shown at the bottom of the next page. The skewness and kurtosis can be derived also in explicit forms by using the below quantities.
The DI can be calculated by dividing the sample variance to sample mean. When the DI is equal to one, it indicates the equi-dispersion. When the DI is greater than one, it indicates the over-dispersion, opposite case indicates the underdispersion. Table 1 presents some numerical results of the mean, variance, DI, skewness and kurtosis for the DLi-3P distribution for different values of the model parameters. From Table 1, the DLi-3P model is appropriate for modeling under(over)-dispersed data sets. Moreover, this model is capable of modeling positively skewed and leptokurtic data sets.

C. MEAN RESIDUAL LIFE AND MEAN PAST LIFE
The mrl and mpl are two commonly used measures to study the ageing behavior of a component or a system of components. The mrl is used to model the burn-in and maintenance of the component. The mrl is defined as (16) where N 0 = {0, 1, 2, 3, . . . , q} for 0 < q < ∞. Let T be a DLi-3P random variable. Then, the mrl, say i , can be expressed in a closed form as The other reliability measure is mpl which measures the time elapsed since the failure of T given that the system has failed sometime before i. The mpl, * i , is given by * where * 0 = 0. If T be a DLi-3P random variable, then the mpl can be represented in a closed form as * The mean of the distribution function can be expressed as The reversed hrf (rhrf) and the mpl are related as If T be a DLi-3P random variable, then the cdf can be recovered by the MPL as

D. STRESS-STRENGTH ANALYSIS
Stress-strength (Str-Sth) has many applications in different scientific fields. Let X Str be a stress and X Sth be a strength of VOLUME 8, 2020 the system. The expected reliability (R Str−Sth ) can be calculated by Let X Str ∼ DLi-3P(λ 1 , α 1 , β 1 ) and X Sth ∼ DLi-3P(λ 2 , α 2 , β 2 ). Then, R Str−Sth can be represented in a closed form as From (23), the value of R Str−Sth depends only on the model parameters. Some numerical results of R Str−Sth are reported in Table 2 by using the DLi-3P distribution for the parameters α 1 = α 2 = 0.01 and β 1 = β 2 = 0.5.

E. GENERATING RANDOM VARIABLES FROM DLi-3P DISTRIBUTION
We introduce an algorithm to generate random variables from the DLi-3P distribution. The below algorithm could be used for this purpose.

3) Compute
4) X = Z To generate random sample of size n from the DLi-3P distribution, the steps 1-4 should be repeated n times. The function W − (·) represents the negative branch of the Lambert-W function.

IV. MAXIMUM LIKELIHOOD ESTIMATION
Assume that the random sample x 1 , x 2 , . . . , x n come from the DLi-3P distribution with unknown parameters λ, α and β. The log-likelihood function of the DLi-3P is By differentiating (25) with respecto the unknown parameters, we have, ∂ ∂λ , ∂ ∂α , and ∂ ∂β , as shown at the bottom of this page.
The simultaneous solutions of these likelihood equations give the MLEs of the model parameters. However, these equations cannot be solved analytically; therefore, an iterative procedure like Newton-Raphson is required to solve it numerically. Here, we use the constrOptim function of R software to maximize the log-likelihood function of the DLi-3P distribution given in (25). The standard errors of the estimated parameters are obtained by means of the squared root of the inverse of the hessian matrix evaluated at estimated model parameters. The fdHess function of R software is used to obtain hessian matrix.

V. THE COMPUTATIONAL COMPLEXITY AND LIMITATIONS OF THE DLi-3P DISTRIBUTION
Proposing a new distribution with adding one or more additional shape parameters increases the model complexity. The proposed distribution, DLi-3P, contains three parameter which requires a good initial parameter vector in estimation step. More importantly, the domain of the parameter λ is (0, 1). The estimation of the parameter λ may yield a value which is outside of its domain.
To overcome these problems one should take into consideration the points given below.
The initial parameter vector should be correctly determined. The generalized simulated annealing method is implemented to obtain a resonable initial vector. The GenSA package of R software is used for this purpose.
Since the domain of the parameter λ is (0, 1), the constrained optimization algorithm should be used to obtain estimated of the parameter λ. The constrOptim function of R software is used for this purpose.

. ESTIMATION OF INAR(1)DLi-3P PROCESS
The three estimation methods are generally used to estimate the unknown parameters of INAR(1) process. Pr These estimation methods are Yule-Walker (YL), conditional least squares (CLS) and conditional maximum likelihood estimation (CMLE). The relative efficiencies of these estimation methods have been discussed in several researches based on the simulation studies (see Bourguignon et al., 2019, andLívio et al., 2018). According to these simulation studies, CMLE method performs better than other two estimation methods for both small and large sample sizes. Based on these facts, we prefer the CMLE method to obtain the unknown parameters of INAR(1)DLi-3P process. The conditional loglikelihood function of the INAR(1)DLi-3P process is where = (p, λ, α, β) represents the parameter vector to be estimated. It is not possible to obtain the explicit formulations of the CMLE of the parameters of the INAR(1)DLi-3P process. Therefore, (34) has to be maximized by using the statistical software such as R, Matlab, S-Plus or Python. Here, we use the constrOptim function of the R software to minimize the minus of the log-likelihood function given in (34). The standard errors of the estimated parameters are obtained by means of the squared roots of the diagonal elements of the Hessian matrix whose elements are numerically calculated by using the fdHess function of the R software. The initial parameter vector of the INAR(1)DLi-3P process is obtained by GenSA package of the R software.

VII. SIMULATION
The finite-sample performance of the MLEs of the parameters of the DLi-3P distribution is investigated by a simulation study. The below simulation procedure is implemented for this purpose.
1) Set the simulation replication number is 1000.
2) Set the parameters of DLi-3P distribution λ = 0.11, α = −0.12 and β = 0.30. 3) Using the given parameter values, generate random variables with sample size n = 5, 10, 15, . . . , 40 from the DLi-3P by repeating N times. 4) For each generated sample size, obtain the λ j , α j and β j , j = 1, 2, . . . , N . 5) Compute the estimated biases and mean-squared errors (MSEs). The required equations can be found in Altun (2020a). The simulation results are summarized graphically in Figure 3. From this figure, we conclude that the estimated biases approach to the desired value, zero, for large

VIII. APPLICATIONS
In this section, we analyze three real data sets by using developed models in the previous sections of the presented study. In the first application, the suitable probability distribution for the numbers of daily deaths from the coronavirus in Iran is investigated. In the second application, the length of stay in a psychiatric ward is analyzed. In the third application, the monthly counts of the larceny calls in Pittsburgh are predicted by INAR(1)DLi-3P process. The fitted models are compared utilizing some criteria, namely, the negative maximized log-likelihood (− ), Akaike information criterion (AIC), corrected Akaike information criterion (CAIC), Bayesian information criterion (BIC), Hannan-Quinn information criterion (HQIC), Chi-square (χ 2 ) test with its corresponding p-value, and Kolmogorov-Smirnov (K-S) test with its corresponding p-value. The below strategy is used to decide best fitted model.
The AIC, CAIC, BIC, HQIC, χ 2 , − and KS test with its p-value are computed for all competitive models as well as DLi-3P distribution The models with p-values greater than 0.05 are identified as potential models. Among the potential models, the best model is determined as the model with the smallest values of the AIC, CAIC, BIC, HQIC, chi 2 , − ell. The computational results are carried out in R software. The used computer features are: Intel Core i5-826U CPU

A. CORONAVIRUS
The firs data set is reported in https://www.worldometers. info/coronavirus/country/iran/ and represents the daily new deaths in Iran from 15 February to 10 March, 2020. We compare the fits of the DLi-3P model with some competitive models such as DLi, DR, DIW, DB-XII, DPa, DBH and Poisson (Poi) models. The all used competitive models, except the Poi distribution, enables to model over-dispersion. To be fair in comparison, we choose the over-dispersed models. The MLEs with their corresponding standard error (SE), confidence interval (C. I) for the parameter(s) and goodness of fit test for the coronavirus data set are listed in Tables 3 and 4, respectively. As seen from the reported values in Table 4, DLi-3P distribution has the lowest values of the goodness-offit statistics which is evidence to conclude that the DLi-3P distribution is more suitable probability distribution than other competitive models for the data used. The execution time of the DLi-3P model for the coronavirus data set is 0.65 seconds. Figure 4 shows the profile log-likelihood functions of the DLi-3P distribution. It is clear that the estimated parameters are maximizers of the log-likelihood function. Figures 5 and 6 display the estimated cdfs and P-P plots of the fitted distributions for the data used. From these figures, we conclude that the DLi-3P distribution provides the best fits among others. Table 5 lists the mean, variance, DI, skewness and kurtosis values of the fitted DLi-3P distribution. As seen from these results, the fitted DLi-3P distribution right skewed and leptokurtic.

B. PSYCHIATRIC WARD
The data presented herein give the length of stay on a psychiatric ward for 67 Male patients (see Chakraborty and Gupta, 2015). The MLEs with their corresponding SE, C. I for the parameter(s) of the fitted distributions and the goodness of fit test results for the data used are listed in Tables 6 and 7, respectively. Since the DLi-3P distribution has the lowest values of the goodness of fit statistics with highest p-value, it could be selected as a best model among others. The execution time of the DLi-3P model for the length of stay data set is 0.05 seconds.
In the Figure 7, the profile log-likelihood functions of the DLi-3P distribution for the fitted data set are plotted to demonstrate that the estimated parameters of the DLi-3P   distribution are the maximizers of the log-likelihood function. As seen from the Figure 7, the estimated parameters of the DLi-3P distribution are the maximizers of the log-likelihood functions of the DLi-3P distribution.      Table 8 lists the mean, variance, DI, skewness and kurtosis values of the fitted DLi-3P distribution. As in first application, the fitted DLi-3P distribution is right-skewed and leptokurtic.

C. LARCENY CALLS
Here, the importance of INAR(1)DLi-3P process is compared with INAR(1)NPWE, INAR(1)P, INAR(1)G and INAR(1)PL processes. The required formulation for the translation probabilities of these competitive models can be found in Altun (2020b). The best fitted model is selected based on the information criteria, AIC and BIC statistics. We analyze the crime data set on the monthly counts of 911 larceny calls which contains 144 observations between Jan. 1990 and Dec. 2001. The data is was available in http:// www.forecastingprinciples.com/index.php/crimedata. Firstly, we investigate whether the data used displays over-dispersion problem. To do this, we calculate the mean, variance and DI of the data set. The following results are obtained, respectively, 19.951, 39.613 and 1.985. Then, the hypothesis test for over-dispersion, proposed by Schweer and Weiß (2014), is applied to decide the whether the observed over-dispersion is statistically significant. The obtained test statistic is 15.936 and its p-value is less than 0.001 which reveal that the data display significant over-dispersion. The fundamental plots of the data used such as autocorrelation function (ACF), partial ACF (PACF), histogram and time series plots are displayed in Figure 10. From the Figure 10 we conclude that the INAR(1) process could be a possible model for this data set since the only first lag is significant in PACF plot. VOLUME 8, 2020 The estimated parameters of the competitive models as well as INAR(1)DLi-3P model and corresponding standard errors, also AIC and BIC statistics, are listed in Table 9. As seen from the results given in Table 9, the INAR(1)DLi-3P model has the lowest values of the AIC and BIC statistics which are evidence to conclude that the INAR(1)DLi-3P model provides better fits than other competitive models for the data used. The execution time of the INAR(1)DLi-3P model for the larceny data set is 1.97 seconds. To check the accuracy of the fitted INAR(1)DLi-3P process, the residual analysis is conducted. We calculate the Pearson residuals of the INAR(1)DLi-3P process by using the following equation where E ( X t | X t−1 ) and Var ( X t | X t−1 ) are given in (32) and (33) (1)DLi-3P process for the data used can be given as follows where the innovation process is The predicted values of the larceny calls obtained by the INAR(1)DLi-3P process and the ACF plot of the Pearson residuals are displayed in 11. The ACF plot of the Pearson residuals confirms that the residuals are uncorrelated.

IX. ANALYSIS OF RESULTS
In this section, we interpret the empirical results more efficiently. In previous section, we analyze three data set to convince the readers in favour of the DLi-3P distribution. The all used data sets are over-dispersed. The competitive models, DLi, DR, DIW, DB-XII, DPa and DBH distributions, except the Poisson distribution, have good properties to model the over-dispersion. However, a few of them achieve to demonstrate acceptable fit to used data sets such as DLi, DIW and DB-XII distributions. The parameter λ of the DLi-3P distribution controls the shape of the distribution and thanks to the parameter λ, the DLi-3P distributions gains much more flexibility than the other competitive models. More importantly, the skewness, kurtosis and DI measures of the DLi-3P distribution has wider range than those of competitive models. Additionally, the profile log-likelihood plots of the DLi-3P distribution reveal the righteousness of the used strategy in estimating the unknown parameter vector of the DLi-3P distribution.

X. CONCLUSIONS
This paper introduces a new three-parameter discrete distribution, shortly DLi-3P distribution. The statistical properties of the DLi-3P distribution are derived in great detail. The maximum likelihood estimation method is used to obtain unknown parameters of the proposed distribution and a brief simulation study is given to discuss the performance of the maximum likelihood estimators of the DLi-3P distribution for both small and large sample sizes. More importantly, a new INAR(1) process with DLi-3P innovations are introduced and studied. The three real data sets are analyzed to convince the readers in favour of the DLi-3P distribution against the other competitive models. We believe that the DLi-3P will increase its popularity and find a wider range of application area in different scientific fields.