# EM Algorithm for Estimating the Parameters of Quasi-Lindley Model with Application

Authors: M. Kayid; Nassr S. Al-Maflehi
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The quasi-Lindley distribution is a flexible model useful in reliability analysis, management science, and engineering analysis. In this paper, an expectation-maximization (EM) algorithm was applied to estimate the parameters of this model for uncensored and right-censored data. Simulation studies show that the estimates of EM perform better than maximum likelihood estimates (MLEs) for both uncensored and censored data. In an illustrative example, the waiting times of a bank’s customers are analyzed and the estimator of the EM algorithm is compared with the MLE. The analysis of the data can be useful for the management of the bank.


---

## Body

## 1. Introduction

The quasi-Lindley distribution proposed by Shanker and Mishra [1] is a generalization of the Lindley distribution introduced by Lindley [2] and is quite useful in reliability theory and survival analysis. The probability density function (pdf) of the quasi-Lindley distribution is given by(1)fx=λα+1α+λxe−λx,x>0,α>0,λ>0,and is a mixture of gamma distributions G1,λ and G2,λ with weights α/α+1 and 1/α+1, respectively. The hazard rate function of the quasi-Lindley model is(2)hx=λα+λx1+α+λx,x>0,α>0,λ>0,which is an increasing function.An important feature of the quasi-Lindley model is that, unlike the Lindley model and its many other generalizations, it is scale invariant. Nevertheless, it is not complicated but sufficiently flexible. Shanker and Mishra [1] studied some of its basic properties and dynamic reliability measures. They also discussed the maximum likelihood estimator (MLE) for its parameters. The MLE is theoretically consistent and efficient, but in practice, it strongly depends on the initial values and the computational approach, which can be achieved by directly maximizing the log-likelihood function or by solving the likelihood equations. Moreover, the simulation results for the MLE of the quasi-Lindley distribution (especially for α) show extremely large values for the mean square error (MSE) (see Tables 1 and 2). This motivates us to investigate the EM algorithm for estimating the parameters.Table 1
Simulation results for MLE and EM estimators of theQL distribution for uncensored data (in every cell, the first and second lines are related to α and λ, respectively).
Methodn100200α,λBMSECPCILMBMSECPCILMMLE0.4, 0.020.6081131.030.891829.480.2018113.330.90923.200.000050.0000050.9540.00990.000080.0000020.94840.00650.5, 0.52.99552514.50.880457.310.125026.720.914.8490.00260.00370.94720.23820.001250.001810.95540.16490.7, 0.17.99125877.700.8742215.570.6875251.7280.901819.1970.000410.000190.94540.05400.000450.000090.95740.0372EM0.4, 0.02−0.000070.012000.98681.19150.003040.005990.99940.80470.000170.0000020.9790.009170.000110.0000010.99360.00620.5, 0.50.003940.016730.98721.44740.001870.008630.99920.98880.003570.001640.98020.24200.001550.000800.99520.16640.7, 0.10.007670.028470.97622.09850.000210.013510.99761.43790.000830.000060.9740.05480.000510.000030.99440.0378Table 2
Simulation results for MLE and EM estimators of theQL distribution when the censorship rate is p=0.2 (in every cell, the first and second lines are related to α and λ, respectively).
Methodn100200α,λBMSECPCILMBMSECPCILMMLE0.4, 0.023.09912996.3210.8758171.840.2813258.77140.929.9925−0.000130.0000080.94520.0107−0.000090.0000040.95880.00750.5, 0.56.44872928.100.867439.530.6962126.66130.905461.3642−0.003630.00630.9390.3021−0.002350.002890.96040.20620.7, 0.115.6256471.680.853958.732.6733556.140.8947215.307−0.001390.000310.91820.0699−0.000990.000160.95120.0498EM0.4, 0.020.005380.011080.98121.32330.002310.005570.99760.87620.000080.0000030.9780.011130.000090.00000150.9940.00750.5, 0.50.002930.015060.96981.62630.002020.007520.99661.10440.003270.001930.9670.29930.001010.000960.99440.20470.7, 0.10.007850.024330.94522.39500.003800.011930.98721.65840.000290.000070.94340.069180.000380.000040.98560.0482In statistics, when the data are collected from a mixture or competing risks model, the EM algorithm is an effective tool for estimating parameters of latent variable models. In their foundational work, Dempster et al. [3] introduced the EM algorithm. Many authors after them have used the idea of EM in their work to provide better estimation of the parameters of the models they are considering. For example, Elmahdy and Aboutahoun [4] and Almhana et al. [5] used the EM algorithm to estimate the parameters for a mixture of Weibull models and a mixture of gamma models, respectively. In addition, Bee et al. [6] applied the EM algorithm to estimate the parameters of Pareto mixture models, Ghosh et al. [7] used the EM algorithm for a mixture of Weibull and Pareto (IV) models, and Balakrishnan and Pal [8] used the EM-based likelihood inference in their work. For detailed discussions of the EM algorithm, we refer the readers to McLachlan and Krishnan [9] and Mengersen et al. [10]. In addition, Wu [11] proved some results related to the convergence of the EM algorithm.In this paper, a specific EM algorithm is developed to obtain a more reliable estimate for the parameters of the quasi-Lindley distribution for uncensored and right-censored data. The paper is organized as follows. Section2 discusses the EM algorithm for the quasi-Lindley distribution when the data were uncensored. In Section 3, the EM algorithm is extended to right-censored data. Section 4 examines the behavior of the MLE and EM estimates and compares them through simulations. In Section 5, both MLE and EM estimates are computed for a real dataset. Finally, Section 6 draws the conclusions for the paper.


## 2. Uncensored Data

Assume thatx1,x2,…,xn be an independent and identically distributed (iid) random sample from quasi-Lindley distribution with parameters α,λ, briefly QLα,λ. The log-likelihood function of the parameters is(3)lα,λ;x=nlnλ−nlnα+1+∑i=1nlnα+λxi−λ∑i=1nxi.The likelihood equations can be obtained by partial differentiation of this log-likelihood function with respect toα and λ as follows:(4)∂∂αlα,λ;x=∑i=1n1α+λxi−nα+1=0,∂∂λlα,λ;x=nλ+∑i=1nxiα+λxi−∑i=1nxi=0.The MLE can be calculated by directly maximizing the log-likelihood function (3) directly or by solving the likelihood equations. Let l=lnfX, and the Fisher information matrix of the quasi-Lindley distribution is(5)K=E−∂2l∂α2E−∂2l∂α∂λE−∂2l∂λ∂αE−∂2l∂λ2.When we have an iid random sampleXi, i=1,2,…,n from QLα0,λ0, the MLE, α^,λ^, weakly converges to bivariate normal Nα0,λ0,n−1K−1 where K−1 is the inverse of the information matrix.


### 2.1. The EM Algorithm for the Complete Data

SinceQLα,λ is a mixture of two gamma distributions G1,λ and G2,λ, the EM algorithm can be used to estimate its parameters. Let Xi, i=1,2,…,n, be an iid random sample of QLα,λ. In the EM approach, for each Xi, we consider one latent random variable Zi which determines that Xi belongs to G1,λ or G2,λ. In other words, Xi|Zi=j∼IGj,θ, j=1,2, PZi=1=α/α+1, and PZi=2=1/α+1. For a brief representation, take θ=α,λ. The likelihood function can be written in the following form.(6)Lθ;x,z=∏i=1n∏j=12gjxi;θPZi=jIzi=j,where the indicator Izi=j equals 1 when zi=j and otherwise equals to 0. Also,(7)gjxi;λ=λjxij−1exp−λxi,is the pdf of the underlying gamma distribution and(8)PZi=j=1α+1αIj=1.Then, the log-likelihood function is(9)lθ;x,z=lnLθ;x,z=∑i=1n∑j=12Izi=jlnλjxij−1exp−λxi1α+1αIj=1.The EM algorithm goes through two steps: the expectation step (E) and the maximization step (M). In each iteration, the E step constructs the expected value of the log-likelihood with respect to the current estimate of the conditional latent variable. In the M step, the constructed in the E step is maximized to provide the estimates. The iterative process can be terminated when the improvement in the expectation function falls below a predetermined small value.


#### 2.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, the conditional distribution of Zi is obtained by Bayes theorem:(10)pij,t=PZi=j|Xi=xi,θt=fXi=xi|Zi=j,θtPZi=j|θtfXi=xi|θt=λtjxij−1e−λtxiαt+1−1αtIj=1∑j=12λtjxij−1e−λtxiαt+1−1αtIj=1,i=1,2,…,n,j=1,2,and after simplification, we have(11)pi1,t=αtαt+λtxi,i=1,2,…,n,and pi2,t=1−pi1,t. These probabilities are known as membership probabilities at iteration t and are used to construct the expectation function Qθ|θt as follows:(12)Qθ|θt=EZ|X,θtlθ;x,Z=∑i=1nEZi|Xi,βt∑j=12Izi=jlnλjxij−1exp−λxi1α+1αIj=1=∑i=1nPZi=1|Xi=xi,θtlnλe−λxiαα+1+PZi=2|Xi=xi,θtlnλ2xie−λxi1α+1=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi+∑i=1npi1,tlnα−lnα+1.The last expressions of (12) show that expectation can be expressed as the sum of two statements, one of which depends only on λ and the other only on α, i.e.,(13)Qθ|θt=Q1λ|θt+Q2α|θt,where(14)Q1λ|θt=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi,Q2α|θt=∑i=1npi1,tlnα−lnα+1.


#### 2.1.2. TheM Step

To estimate the parameters att+1 iteration, we maximize the Qθ|θt in terms of θ. So, we have(15)θt+1=argmaxθQθ|θt,which, by (13), reduces to the following separate maximization problems.(16)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,where Q1λ|θt and Q2α|θt are determined by (14) and (15), respectively. By solving the equation ∂/∂λQ1λ|θt=0, the estimation of λ at t+1 iteration is obtained.(17)λt+1=n+∑i=1npi2,t∑i=1nXi.On the other hand, solving the equation∂/∂αQ2α|θt=0, we have(18)αt+1=∑i=1npi1,t∑i=1npi2,t.The sequenceθt will converge to θ, and the iterative process can be concluded when for some predefined small ϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ. This means that further iterations do not improve the objective function considerably. For detailed information about convergence of the EM algorithm, see Wu [11].


## 2.1. The EM Algorithm for the Complete Data

SinceQLα,λ is a mixture of two gamma distributions G1,λ and G2,λ, the EM algorithm can be used to estimate its parameters. Let Xi, i=1,2,…,n, be an iid random sample of QLα,λ. In the EM approach, for each Xi, we consider one latent random variable Zi which determines that Xi belongs to G1,λ or G2,λ. In other words, Xi|Zi=j∼IGj,θ, j=1,2, PZi=1=α/α+1, and PZi=2=1/α+1. For a brief representation, take θ=α,λ. The likelihood function can be written in the following form.(6)Lθ;x,z=∏i=1n∏j=12gjxi;θPZi=jIzi=j,where the indicator Izi=j equals 1 when zi=j and otherwise equals to 0. Also,(7)gjxi;λ=λjxij−1exp−λxi,is the pdf of the underlying gamma distribution and(8)PZi=j=1α+1αIj=1.Then, the log-likelihood function is(9)lθ;x,z=lnLθ;x,z=∑i=1n∑j=12Izi=jlnλjxij−1exp−λxi1α+1αIj=1.The EM algorithm goes through two steps: the expectation step (E) and the maximization step (M). In each iteration, the E step constructs the expected value of the log-likelihood with respect to the current estimate of the conditional latent variable. In the M step, the constructed in the E step is maximized to provide the estimates. The iterative process can be terminated when the improvement in the expectation function falls below a predetermined small value.


### 2.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, the conditional distribution of Zi is obtained by Bayes theorem:(10)pij,t=PZi=j|Xi=xi,θt=fXi=xi|Zi=j,θtPZi=j|θtfXi=xi|θt=λtjxij−1e−λtxiαt+1−1αtIj=1∑j=12λtjxij−1e−λtxiαt+1−1αtIj=1,i=1,2,…,n,j=1,2,and after simplification, we have(11)pi1,t=αtαt+λtxi,i=1,2,…,n,and pi2,t=1−pi1,t. These probabilities are known as membership probabilities at iteration t and are used to construct the expectation function Qθ|θt as follows:(12)Qθ|θt=EZ|X,θtlθ;x,Z=∑i=1nEZi|Xi,βt∑j=12Izi=jlnλjxij−1exp−λxi1α+1αIj=1=∑i=1nPZi=1|Xi=xi,θtlnλe−λxiαα+1+PZi=2|Xi=xi,θtlnλ2xie−λxi1α+1=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi+∑i=1npi1,tlnα−lnα+1.The last expressions of (12) show that expectation can be expressed as the sum of two statements, one of which depends only on λ and the other only on α, i.e.,(13)Qθ|θt=Q1λ|θt+Q2α|θt,where(14)Q1λ|θt=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi,Q2α|θt=∑i=1npi1,tlnα−lnα+1.


### 2.1.2. TheM Step

To estimate the parameters att+1 iteration, we maximize the Qθ|θt in terms of θ. So, we have(15)θt+1=argmaxθQθ|θt,which, by (13), reduces to the following separate maximization problems.(16)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,where Q1λ|θt and Q2α|θt are determined by (14) and (15), respectively. By solving the equation ∂/∂λQ1λ|θt=0, the estimation of λ at t+1 iteration is obtained.(17)λt+1=n+∑i=1npi2,t∑i=1nXi.On the other hand, solving the equation∂/∂αQ2α|θt=0, we have(18)αt+1=∑i=1npi1,t∑i=1npi2,t.The sequenceθt will converge to θ, and the iterative process can be concluded when for some predefined small ϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ. This means that further iterations do not improve the objective function considerably. For detailed information about convergence of the EM algorithm, see Wu [11].


## 2.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, the conditional distribution of Zi is obtained by Bayes theorem:(10)pij,t=PZi=j|Xi=xi,θt=fXi=xi|Zi=j,θtPZi=j|θtfXi=xi|θt=λtjxij−1e−λtxiαt+1−1αtIj=1∑j=12λtjxij−1e−λtxiαt+1−1αtIj=1,i=1,2,…,n,j=1,2,and after simplification, we have(11)pi1,t=αtαt+λtxi,i=1,2,…,n,and pi2,t=1−pi1,t. These probabilities are known as membership probabilities at iteration t and are used to construct the expectation function Qθ|θt as follows:(12)Qθ|θt=EZ|X,θtlθ;x,Z=∑i=1nEZi|Xi,βt∑j=12Izi=jlnλjxij−1exp−λxi1α+1αIj=1=∑i=1nPZi=1|Xi=xi,θtlnλe−λxiαα+1+PZi=2|Xi=xi,θtlnλ2xie−λxi1α+1=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi+∑i=1npi1,tlnα−lnα+1.The last expressions of (12) show that expectation can be expressed as the sum of two statements, one of which depends only on λ and the other only on α, i.e.,(13)Qθ|θt=Q1λ|θt+Q2α|θt,where(14)Q1λ|θt=∑i=1n1+pi2,tlnλ−λxi+pi2,tlnxi,Q2α|θt=∑i=1npi1,tlnα−lnα+1.


## 2.1.2. TheM Step

To estimate the parameters att+1 iteration, we maximize the Qθ|θt in terms of θ. So, we have(15)θt+1=argmaxθQθ|θt,which, by (13), reduces to the following separate maximization problems.(16)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,where Q1λ|θt and Q2α|θt are determined by (14) and (15), respectively. By solving the equation ∂/∂λQ1λ|θt=0, the estimation of λ at t+1 iteration is obtained.(17)λt+1=n+∑i=1npi2,t∑i=1nXi.On the other hand, solving the equation∂/∂αQ2α|θt=0, we have(18)αt+1=∑i=1npi1,t∑i=1npi2,t.The sequenceθt will converge to θ, and the iterative process can be concluded when for some predefined small ϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ. This means that further iterations do not improve the objective function considerably. For detailed information about convergence of the EM algorithm, see Wu [11].


## 3. Right-Censored Data

Consider an iid random sampleXi, i=1,2,…,n, from QLα,λ which is exposed to right censorship. We say that Xi is censored from the right by a censoring random variable Ci, if Xi>Ci, and in this case, the only information about event time is that it is greater than censoring time Ci. The observations consist of Ti=minXi,Ci and di, where di=1, when the event has not been censored, Xi≤Ci, and di=0, when the event has been censored, Xi>Ci. Given a right-censored sample ti,di, i=1,2,…,n, the log-likelihood function is(19)lα,λ;t,d=∑i=1ndilnfti+∑i=1n1−dilnRti,where f and R show the density and the reliability functions of the quasi-Lindley distribution, respectively. The log-likelihood function simplifies to(20)lθ;t,d=n1lnλα+1+∑i=1ndilnα+λti+∑i=1n1−diln1+λtiα+1−∑i=1nλti,where n1=∑i=1ndi and θ=α,λ.


### 3.1. The EM Algorithm for Right-Censored Data

To implement the EM algorithm, we should include the latent variableZi, i=1,2,…,n, defined in the previous section. Then, the likelihood function for the censored data is(21)Lθ;t,d,z=∏i=1n∏j=12gjti;θPZi=jIzi=jdi∏i=1n∏j=12G¯jti;θPZi=jIzi=j1−di,where gj shows the gamma pdf considered in the previous section and G¯j is its corresponding reliability function. By taking logarithm from (21), the log-likelihood function has the following form:(22)lθ;t,d,z=∑i=1ndiIzi=1lnλe−λtiαα+1+Izi=2lnλ2tie−λti1α+1+∑i=1n1−diIzi=1lne−λtiαα+1+Izi=2lnλti+1e−λti1α+1.Similar to the uncensored data, we should iterate twoE and M steps to find an improved estimation.


#### 3.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, applying the Bayes theorem, we can compute the conditional distribution of Zi as follows:(23)pij,t=diPZi=j|Xi=ti,θt+1−diPZi=j|Xi>ti,θt=difXi=ti|Zi=j,θtPZi=j|θtfXi=ti|θt+1−diPXi>ti|Zi=j,θtPZi=j|θtfXi>ti|θt=diλtjtij−1e−λttiαt+1−1αtIj=1∑j=12λtjtij−1e−λttiαt+1−1αtIj=1+1−diIj=1e−λttiαt/αt+1+Ij=2λtti+1e−λtti1/αt+1e−λttiαt/αt+1+λtti+1e−λtti1/αt+1=diIj=1αt+Ij=2λttiαt+λtti+1−diIj=1αt+Ij=2λtti+1αt+λtti+1,i=1,2,…,n,j=1,2.Specifically, forj=1,(24)pi1,t=diαtαt+λtti+1−diαtαt+λtti+1,i=1,2,…,n,and pi2,t=1−pi1,t. Then, using (22), the expectation function at iteration t can be written in the following form.(25)Qθ|θt=EZ|t,d,θtlθ;t,d,z=∑i=1ndipi1,tlnλe−λtiαα+1+pi2,tlnλ2tie−λti1α+1+∑i=1n1−dipi1,tlne−λtiαα+1+pi2,tlnλti+1e−λti1α+1.Similar to uncensored case, it is straightforward to check thatQθ|θt can be written as two statements in which one of them just depends on α and the other depends on λ. More precisely,(26)Qθ|θt=Q1λ|θt+Q2α|θt,where(27)Q1λ|θt=∑i=1ndi1+pi2,tlnλ+∑i=1ndipi2,tlnti+∑i=1n1−dipi2,tlnλti+1−λ∑i=1nti,(28)Q2α|θt=∑i=1npi1,tlnα−lnα+1


#### 3.1.2. TheM Step

In this step, we should maximize theQθ|θt function to compute the estimations at the t+1 iteration.(29)θt+1=argmaxθQθ|θt,which, by (26), reduces to the following separate maximization problems.(30)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,in which Q1λ|θt and Q2α|θt are determined by (27) and (28), respectively. The likelihood equation ∂/∂λQ1λ|θt=0 which after some algebra can be simplified to(31)λ=∑i=1ndi1+pi2,t∑i=1nti+∑i=1n1−diλti/1+λtipi2,t∑i=1ntiwhich does not yield to an analytical solution for λ, so the solution can be computed by numerical methods. But, clearly (31) implies that the solution of this equation, namely, λt+1, satisfies the inequality(32)λt+1>∑i=1ndi1+pi2,t∑i=1nti.On the other hand, since1+pi2,t>λti/1+λtipi2,t, one upper bound for the solution is(33)λt+1<∑i=1n1+pi2,t∑i=1nti,and in turn, by (32) and (33), we have(34)∑i=1ndi1+pi2,t∑i=1nti<λt+1<∑i=1n1+pi2,t∑i=1nti.These bounds can be applied in numerical processes to find optimized answer. The solution forα can be obtained by solving the equation ∂/∂αQ2α|θt=0 as follows:(35)αt+1=∑i=1npi1,t∑i=1npi2,t.Similar to the uncensored case, the iterative process can be concluded when for some predefined smallϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ.Letθ^ and θ0 be the EM estimator and the real parameter, respectively. Then, θ^−θ0 converges asymptotically to a bivariate normal distribution N0,V, where V can be approximated by the inverse of the observed information matrix with respect to the observed data (see Meng and Rubin [12]). It is computed by evaluating the Hessian matrix of the log-likelihood function with respect to the observed data at the point θ^, and then calculating the inverse of the obtained Hessian matrix, briefly V=Io−1θ^|x. Fortunately, in the case of this study, the log-likelihood function of the observed data is not complicated and can be used to calculate the Hessian matrix and finally the variance approximation. For this purpose, the function “hessian” of the library “pracma” in R is used. Since the asymptotic distribution of the EM estimator is normal, the standard normal quantiles are used to obtain approximate confidence intervals of the parameters.


## 3.1. The EM Algorithm for Right-Censored Data

To implement the EM algorithm, we should include the latent variableZi, i=1,2,…,n, defined in the previous section. Then, the likelihood function for the censored data is(21)Lθ;t,d,z=∏i=1n∏j=12gjti;θPZi=jIzi=jdi∏i=1n∏j=12G¯jti;θPZi=jIzi=j1−di,where gj shows the gamma pdf considered in the previous section and G¯j is its corresponding reliability function. By taking logarithm from (21), the log-likelihood function has the following form:(22)lθ;t,d,z=∑i=1ndiIzi=1lnλe−λtiαα+1+Izi=2lnλ2tie−λti1α+1+∑i=1n1−diIzi=1lne−λtiαα+1+Izi=2lnλti+1e−λti1α+1.Similar to the uncensored data, we should iterate twoE and M steps to find an improved estimation.


### 3.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, applying the Bayes theorem, we can compute the conditional distribution of Zi as follows:(23)pij,t=diPZi=j|Xi=ti,θt+1−diPZi=j|Xi>ti,θt=difXi=ti|Zi=j,θtPZi=j|θtfXi=ti|θt+1−diPXi>ti|Zi=j,θtPZi=j|θtfXi>ti|θt=diλtjtij−1e−λttiαt+1−1αtIj=1∑j=12λtjtij−1e−λttiαt+1−1αtIj=1+1−diIj=1e−λttiαt/αt+1+Ij=2λtti+1e−λtti1/αt+1e−λttiαt/αt+1+λtti+1e−λtti1/αt+1=diIj=1αt+Ij=2λttiαt+λtti+1−diIj=1αt+Ij=2λtti+1αt+λtti+1,i=1,2,…,n,j=1,2.Specifically, forj=1,(24)pi1,t=diαtαt+λtti+1−diαtαt+λtti+1,i=1,2,…,n,and pi2,t=1−pi1,t. Then, using (22), the expectation function at iteration t can be written in the following form.(25)Qθ|θt=EZ|t,d,θtlθ;t,d,z=∑i=1ndipi1,tlnλe−λtiαα+1+pi2,tlnλ2tie−λti1α+1+∑i=1n1−dipi1,tlne−λtiαα+1+pi2,tlnλti+1e−λti1α+1.Similar to uncensored case, it is straightforward to check thatQθ|θt can be written as two statements in which one of them just depends on α and the other depends on λ. More precisely,(26)Qθ|θt=Q1λ|θt+Q2α|θt,where(27)Q1λ|θt=∑i=1ndi1+pi2,tlnλ+∑i=1ndipi2,tlnti+∑i=1n1−dipi2,tlnλti+1−λ∑i=1nti,(28)Q2α|θt=∑i=1npi1,tlnα−lnα+1


### 3.1.2. TheM Step

In this step, we should maximize theQθ|θt function to compute the estimations at the t+1 iteration.(29)θt+1=argmaxθQθ|θt,which, by (26), reduces to the following separate maximization problems.(30)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,in which Q1λ|θt and Q2α|θt are determined by (27) and (28), respectively. The likelihood equation ∂/∂λQ1λ|θt=0 which after some algebra can be simplified to(31)λ=∑i=1ndi1+pi2,t∑i=1nti+∑i=1n1−diλti/1+λtipi2,t∑i=1ntiwhich does not yield to an analytical solution for λ, so the solution can be computed by numerical methods. But, clearly (31) implies that the solution of this equation, namely, λt+1, satisfies the inequality(32)λt+1>∑i=1ndi1+pi2,t∑i=1nti.On the other hand, since1+pi2,t>λti/1+λtipi2,t, one upper bound for the solution is(33)λt+1<∑i=1n1+pi2,t∑i=1nti,and in turn, by (32) and (33), we have(34)∑i=1ndi1+pi2,t∑i=1nti<λt+1<∑i=1n1+pi2,t∑i=1nti.These bounds can be applied in numerical processes to find optimized answer. The solution forα can be obtained by solving the equation ∂/∂αQ2α|θt=0 as follows:(35)αt+1=∑i=1npi1,t∑i=1npi2,t.Similar to the uncensored case, the iterative process can be concluded when for some predefined smallϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ.Letθ^ and θ0 be the EM estimator and the real parameter, respectively. Then, θ^−θ0 converges asymptotically to a bivariate normal distribution N0,V, where V can be approximated by the inverse of the observed information matrix with respect to the observed data (see Meng and Rubin [12]). It is computed by evaluating the Hessian matrix of the log-likelihood function with respect to the observed data at the point θ^, and then calculating the inverse of the obtained Hessian matrix, briefly V=Io−1θ^|x. Fortunately, in the case of this study, the log-likelihood function of the observed data is not complicated and can be used to calculate the Hessian matrix and finally the variance approximation. For this purpose, the function “hessian” of the library “pracma” in R is used. Since the asymptotic distribution of the EM estimator is normal, the standard normal quantiles are used to obtain approximate confidence intervals of the parameters.


## 3.1.1. TheE Step

Given the estimate of the parameters at iterationt, θt, applying the Bayes theorem, we can compute the conditional distribution of Zi as follows:(23)pij,t=diPZi=j|Xi=ti,θt+1−diPZi=j|Xi>ti,θt=difXi=ti|Zi=j,θtPZi=j|θtfXi=ti|θt+1−diPXi>ti|Zi=j,θtPZi=j|θtfXi>ti|θt=diλtjtij−1e−λttiαt+1−1αtIj=1∑j=12λtjtij−1e−λttiαt+1−1αtIj=1+1−diIj=1e−λttiαt/αt+1+Ij=2λtti+1e−λtti1/αt+1e−λttiαt/αt+1+λtti+1e−λtti1/αt+1=diIj=1αt+Ij=2λttiαt+λtti+1−diIj=1αt+Ij=2λtti+1αt+λtti+1,i=1,2,…,n,j=1,2.Specifically, forj=1,(24)pi1,t=diαtαt+λtti+1−diαtαt+λtti+1,i=1,2,…,n,and pi2,t=1−pi1,t. Then, using (22), the expectation function at iteration t can be written in the following form.(25)Qθ|θt=EZ|t,d,θtlθ;t,d,z=∑i=1ndipi1,tlnλe−λtiαα+1+pi2,tlnλ2tie−λti1α+1+∑i=1n1−dipi1,tlne−λtiαα+1+pi2,tlnλti+1e−λti1α+1.Similar to uncensored case, it is straightforward to check thatQθ|θt can be written as two statements in which one of them just depends on α and the other depends on λ. More precisely,(26)Qθ|θt=Q1λ|θt+Q2α|θt,where(27)Q1λ|θt=∑i=1ndi1+pi2,tlnλ+∑i=1ndipi2,tlnti+∑i=1n1−dipi2,tlnλti+1−λ∑i=1nti,(28)Q2α|θt=∑i=1npi1,tlnα−lnα+1


## 3.1.2. TheM Step

In this step, we should maximize theQθ|θt function to compute the estimations at the t+1 iteration.(29)θt+1=argmaxθQθ|θt,which, by (26), reduces to the following separate maximization problems.(30)λt+1=argmaxλQ1λ|θt,αt+1=argmaxαQ2α|θt,in which Q1λ|θt and Q2α|θt are determined by (27) and (28), respectively. The likelihood equation ∂/∂λQ1λ|θt=0 which after some algebra can be simplified to(31)λ=∑i=1ndi1+pi2,t∑i=1nti+∑i=1n1−diλti/1+λtipi2,t∑i=1ntiwhich does not yield to an analytical solution for λ, so the solution can be computed by numerical methods. But, clearly (31) implies that the solution of this equation, namely, λt+1, satisfies the inequality(32)λt+1>∑i=1ndi1+pi2,t∑i=1nti.On the other hand, since1+pi2,t>λti/1+λtipi2,t, one upper bound for the solution is(33)λt+1<∑i=1n1+pi2,t∑i=1nti,and in turn, by (32) and (33), we have(34)∑i=1ndi1+pi2,t∑i=1nti<λt+1<∑i=1n1+pi2,t∑i=1nti.These bounds can be applied in numerical processes to find optimized answer. The solution forα can be obtained by solving the equation ∂/∂αQ2α|θt=0 as follows:(35)αt+1=∑i=1npi1,t∑i=1npi2,t.Similar to the uncensored case, the iterative process can be concluded when for some predefined smallϵ>0, Qθt+1|θt+1<Qθt|θt+ϵ.Letθ^ and θ0 be the EM estimator and the real parameter, respectively. Then, θ^−θ0 converges asymptotically to a bivariate normal distribution N0,V, where V can be approximated by the inverse of the observed information matrix with respect to the observed data (see Meng and Rubin [12]). It is computed by evaluating the Hessian matrix of the log-likelihood function with respect to the observed data at the point θ^, and then calculating the inverse of the obtained Hessian matrix, briefly V=Io−1θ^|x. Fortunately, in the case of this study, the log-likelihood function of the observed data is not complicated and can be used to calculate the Hessian matrix and finally the variance approximation. For this purpose, the function “hessian” of the library “pracma” in R is used. Since the asymptotic distribution of the EM estimator is normal, the standard normal quantiles are used to obtain approximate confidence intervals of the parameters.


## 4. Simulations

In a simulation study, we investigate the behavior of the MLE and EM estimators and compare them. The fact that the quasi-Lindley model is a mixture of gamma distributions is applied to generate random samples. To generate right-censored sampley1,y2,…,yn, we assume that the censoring random variable Ci follows the degenerate distribution with mean M. Thus, if p is the censoring rate, we can calculate M by solving the equation M=F−11−p where F−1 is the inverse of the distribution function of the quasi-Lindley model. Now, an uncensored sample x1,x2,…,xn is taken from the quasi-Lindley model. Then, the ith instance of the desired right-censored sample is yi=minxi,M.Each cell of Tables1 and 2 shows the results of one run. In every run, r=5000 replicates of samples of size n=100 or 200 were generated by the quasi-Lindley model with selected parameters, and in each run, the MLE and EM estimators were calculated. To calculate the MLE, the log-likelihood function was maximized by using the “optim” function built into R with the standard “Nelder–Mead” optimization method. In both the maximum likelihood method and EM, the initial values are generated from a uniform distribution. Note that checking the termination conditions of the EM process in each EM iteration results in very slow and time-consuming runs. Therefore, the EM algorithm has been tested many times to find a suitable constant for the number of iterations. In this way, we find that 5 iterations are sufficient.Four measures bias (B), mean squared error (MSE), coverage probability (CP), and confidence interval length mean (CILM) forα and λ have been computed. The B and MSE for α are defined to be(36)B=1r∑i=1rα^i−α,MSE=1r∑i=1rα^i−α2,CP=1r∑i=1rIα∈CIαi,ILM=1r∑i=1rlengthofCIαi,where αi^ shows the MLE/EM estimator in the run i and CIαi shows an approximate asymptotic 95 percent confidence interval for α in the ith iteration (see the last paragraph of Section 3). Also, the indicator function I in (30) equals 1 when the real parameter falls inside the confidence interval and otherwise equals zero. These measures are defined for λ similarly. Tables 1 and 2 present the simulation results for uncensored data and censored data with censorship 0.2, respectively. The main observations from these tables are listed in the following:(i)
The MSE decreases as sample size increases, for both MLE and EM estimators and both uncensored data and censored data which indicates that the MLE and EM estimators are consistent.(ii)
The EM estimator outperforms the MLE in terms of the MSE.(iii)
The results show higher CPs and lower CILMs for EM than MLE. Moreover, the CP increases and CILM decreases as sample size increases.


## 5. Application

Table3 shows 100 waiting times of customers of a bank analyzed by Shanker [13]. The quasi-Lindley distribution was fitted to this dataset, and the parameters were estimated using the maximum likelihood method and EM. The “optim” function in the R language was used to calculate the MLE. Table 4 shows the results of the fitting. In terms of the KS, Anderson–Darling (AD), and Cramer–von Mises (CVM) statistics, both methods provide a good fit, but EM outperforms MLE in a close competition. The empirical and the fitted CDFs are shown in Figure 1(a) and also confirm a good fit. The histogram and estimated probability density function are also shown in Figure 1(b). Using the Hessian matrix calculated with the optim function, the variances of the MLE are estimated for the parameters, v^α^=4×10−12 and v^λ^=0.00020. Using these variance estimates and standard normal quantiles, the 95% confidence intervals for α and λ are 0,0.000006 and 0.1744,0.2305, respectively. The left bound of the α confidence interval was a negative value; by the fact that α>0, it was set to 0.Table 3
Waiting times (in minutes) of customers to receive service in a bank.
0.80.81.31.51.81.91.92.12.62.72.93.13.23.33.53.64.04.14.24.24.34.34.44.44.64.74.74.84.94.95.05.35.55.75.76.16.26.26.26.36.76.97.17.17.17.17.47.67.78.08.28.68.68.68.88.88.98.99.59.69.79.810.710.911.011.011.111.211.211.511.912.412.512.913.013.113.313.613.713.914.115.415.417.317.318.118.218.418.919.019.920.621.321.421.923.027.031.633.138.5Table 4
Results of fitting quasi-Lindley model to dataset of Table3.
Methodα^,λ^AICKSADCVMp valuep valuep valueMLE(0.000002, 0.2025)638.60140.04420.19140.03010.98420.97430.9621EM(0.0118, 0.2013)—–0.041980.18580.02790.99460.99380.9828Figure 1
(a) The empirical distribution and fitted quasi-Lindley distribution for dataset of Table3. (b) The histogram and the estimated probability density function for this dataset.
(a)(b)To find the variances of the EM estimator of the parameters, the bootstrap method is used. In this way,r=1000 samples are derived by the function “sample” of R. Then, for each sample, the EM estimates of the parameters are computed. The estimate of the variance of the EM estimators is approximated by the variance of these estimates which are v^α^=0.0015 and v^λ^=0.00023. For each of the parameters, the 2.5% and 97.5% quantiles of the EM estimator can be considered as upper and lower bounds of the 95% confidence intervals. Then, the 95% confidence intervals for α and λ are 0.00024,0.12845 and 0.1750,0.2325, respectively.


## 6. Conclusion

The quasi-Lindley distribution is a scale-invariant version of the Lindley distribution with a shape parameterα and a scale parameter λ and is a simple yet flexible model in reliability theory, survival analysis, management science, and many other fields. The MLE and EM approaches were investigated to estimate the parameters of this model. The simulation results show that the EM algorithm is better than the MLE for estimating the parameters for both uncensored and censored data.



---

*Source: 8467291-2022-01-25.xml*