CC BY
88
A New Mixed Poisson Distribution for Over-dispersed Count Data: Theory and Applications
Ramajeyam Tharshan •
Postgraduate Institute of Science, University of Peradeniya, Peradeniya, Sri Lanka Department of Mathematics and Statistics, University of Jaffna, Jaffna, Sri Lanka
tharshan10684@gmail.com
Pushpakanthie Wijekoon •
Department of Statistics and Computer Science, University of Peradeniya, Peradeniya,Sri Lanka
pushpaw@pdn.ac.lk
Abstract
In this paper, an alternative mixed Poisson distribution is proposed by amalgamating Poisson distribution and a modification of the Quasi Lindley distribution. Some fundamental structural properties of the new distribution, namely the shape of the distribution and moments and related measures, are explored. It was noted that the new distribution to be either unimodal or bimodal, and over-dispersed. Further, it has a tendency to accommodate various right tail behaviors and variance-to-mean ratios. Its unknown parameter estimation by using the maximum likelihood estimation method is examined by a simulation study based on the asymptotic theory. Finally, two real-world data sets are used to illustrate the flexibility and potentiality of the new distribution.
Keywords: over-dispersion, mixed Poisson distribution, Lindley distribution, Quasi Lindley distribution, goodness of fit.
1. Introduction
Most of the real-world applications, especially, reliability, actuarial, biomedical, engineering, ecological sciences, and among others, the variable of interest is in the form of count data. The Poisson distribution is a standard tool to model the count data if the empirical and theoretical properties satisfy the related underline assumptions. A random variable X is said to have a Poisson distribution with parameter A if both the E(X) and Var(X) of the distribution equal to the parameter A. This property is commonly known as equidispersion. Even though its probability mass function (pmf) is very flexible to compute its probabilities, in some real-world applications the Poisson distribution fails to match empirical observations. Here the variance of the observed data exceeds the theoretical variance. This phenomenon is explained as over-dispersion or variation inflation (Greenwood and Yule, 1920). The over-dispersion occurs by the failure of the basic assumptions of the Poisson distribution. The reasons might be by phenomena of the clustered structure of the population or population is heterogeneous, and heavy right tail that cannot accommodate by the Poisson distribution (McCullagh et al., 1989; Ridout et al., 1998). The heterogeneity of a population is determined by the Poisson parameter A which differs individual to individual, and then the Var(X) = tE(X); (t > 1), where t is called the index of dispersion.
The mixed Poisson distributions are well-known flexible modeling methods to explain the
heterogeneity of the Poisson parameter as well as heavy right tail behaviors (Feller,1943; Shaked, 1980). The mixed Poisson distribution is a resultant distribution or unconditional distribution by assuming that the Poisson parameter is a random variable that has as a parameterized distribution P. The distribution P and its parameter vector © are called prior distribution and hyperparameter, respectively. Then, the resultant distribution of the random variable X can be expressed mathematically in the following form
!• TO
fx(x) = 0 fx\A (x| A)fP (A)dA, (1)
where X| A has a Poisson distribution with parameter A as
e-AAx
fx\a(x|A) = ——, x = 0,1,2..... A > 0, (2)
x.
A is the random variable of the Poisson parameter A, and fp (A) is the density function of the assumed continuous distribution P to the Poisson parameter A. Hence, the random variable X has the same support of X|A with parameter(s) of the prior distribution. Further, Lynch (1988) showed that the form of the mixing distribution has ascendancy over to the form of the resultant mixed distribution.
In literature, researchers assumed the standard lifetime distributions to model the Poisson parameter A as a classical approach. Greenwood and Yule (1920) used the gamma distribution, and the resultant distribution is negative binomial (NBD). Johnson et al., (1992) assumed the exponential distribution to model the Poisson parameter, and the resultant distribution is geometric distribution (GD). Even though NBD and GD are computationally flexible pmfs, they are not befitting distributions for a higher value of t and long right tail. In this context, researchers have assumed several modifications of the standard lifetime distributions for the Poisson parameter. They were modified to have more flexibility in their shapes and failure rate criteria than the standard lifetime distributions. Confluent Hypergeometric series, Gamma product ratio, Generalized gamma, Shifted gamma, Inverse gamma, and Modified Bessel of the 3rd kind are used by Bhattacharya (1967), Irwin (1975), Albrecht (1984), Ruohonen (1988), Willmot (1993), and Ong and Muthaloo (1995), respectively to model the Poisson parameter. The pmfs of such distributions are derived through the recursive formulas or Laplace transform technique, or by using the special mathematical functions. Hence, computing the probabilities of such distributions is complicated
and they are limited in practice.
The Lindley distribution (LD) is one of the life time distributions introduced by Lindley (1958) having the density function
/a(A) = Tte (1 + A)e-eA, A > 0, e > 0, (3)
where e is the shape parameter, and A is the respective random variable. Equation (3) presents a two-component mixture of two different continuous distributions namely exponential (e) and
q
gamma (2, e) distributions with the mixing proportion, p = ---. Sankaran (1970) introduced
1 + e
the one-parameter discrete Poisson-Lindley distribution (PLD) by combining the Poisson and LD. Its pmf is given as
fx(x) = ^, x = 0,1,2..... e > 0. (4)
Note that the pmf of the PLD is an explicit form. Then, obtaining its probabilities is computationally flexible. However, the PLD flexibility is limited to fit various types of the over-dispersed count data sets since it has only one parameter. Then, as an alternative to PLD, Bhati et al. (2015) have obtained the Generalized Poisson-Lindley distribution (GPLD), where the Poisson parameter is distributed to Two-parameter Lindley distribution (Shanker et al., 2013b); Wongrin
Ramajeyam Tharshan, Pushpakanthie Wijekoon RT&A, No 1 (67) A NEW MIXED POISSON DISTRIBUTION_Volume 17, March 2022
and Bodhisuwan (2016) introduced Poisson-generalized Lindley distribution (PGLD), which was obtained by mixing the Poisson distribution with the generalized Lindley distribution (Elbatal et al., 2013); Grine et al. (2017) have obtained the Quasi-Poisson distribution (PQLD) by modeling Poisson parameter to Quasi Lindley distribution (Shanker et al., 2013a). Table 1 summarizes the mixing proportions, mixing components, and parameters of the above continuous distributions that have been used to model the Poisson parameter. We can see that the mixing proportions of the Two-parameter Lindley and generalized Lindley distributions are incorporated with the scale parameter, 9 of the mixing components. Further, the shape parameter of the mixing component gamma (2,9) is fixed with value 2 for the Two-parameter Lindley and the Quasi Lindley distributions. These settings of such mixing distributions may limit the flexibility of the above-mentioned Poisson mixtures to fit well for the various types of the right tail heaviness and t for an over-dispersed count data (Tharshan and Wijekoon, 2020 a b).
Table 1: Mixing proportions, mixing components, and parameters of some modified-Lindley distributions.
Distribution Mixing proportion Mixing components Parameters
shape scale
Two-parameter Lindley 9 9 + a exponential (9) , gamma (2,9) 9, a
Generalized Lindley 9 1 + 9 gamma (a, 9) , gamma (ft, 9) 9, a, ft
Quasi Lindley a a +1 exponential (9) , gamma (2,9) 9 a
The main contribution of this paper is to propose an alternative mixed Poisson distribution for over-dispersed count data to address the above issues. It is obtained by mixing the Poisson distribution and the Modification of the Quasi Lindley distribution (MQLD) (Tharshan and Wijekoon, 2021). The density function of the MQLD(9, a, 5) is given as
ap-8A / \
fA(A; 9, a, 5) = + 1)r(5^r(5)a3 + (9A)3-1j; A > 0,9 > 0, a3 > -1,5 > 0, (5)
where a, and 5 are shape parameters, 9 is a scale parameter, and A is the respective random variable. Equation (5) presents the mixture of two non-identical distributions, exponential (9) ,
a3
and gamma (5,9) with the mixing proportion, p = 3 1 • We can clearly observe that its mixing
proportion p does not incorporate with scale parameter, 9 of the mixing components. Further, the shape parameter of the mixing component gamma distribution, 5 is not fixed with a value. The authors have shown that these settings of the MQLD provide the capability to capture the various ranges of right tail heaviness measured by excess kurtosis (kurtosis -3), horizontal symmetry measured by skewness, and heterogeneity measured by Fano factor (variance-to-mean ratio) by setting its parameter values. Further, its density function can be either unimodal or bimodal.
The remaining part of this paper is organized as follows: In section 2, we introduce the PMQLD with its explicit forms of the probability mass and distribution functions. Its fundamental structural properties are discussed in section 3. The simulation of its random variables and parameter estimations are discussed in section 4. Finally, a simulation study is done to examine the performance of parameter estimation by using the maximum likelihood estimation method, and some real-world examples are taken to show the applicability of the proposed model by comparing it with some other existing Poisson mixtures, NBD, GD, PLD, GPLD, PQLD, PGLD.
2. Formulation of the new mixed Poisson distribution
In this section, we introduce the new mixed Poisson distribution with its pmf and cumulative distribution function (cdf).
Let the random variable X represent the total counts of a specific experiment with mean A. Then, the traditional distribution to calculate probabilities of such outcomes is the Poisson distribution. The PMQLD is obtained by mixing the Poisson and MQLD (Tharshan and Wijekoon, 2021) for over-dispersed count data. The following theorem gives the pmf of the PMQLD.
Theorem 1. Let X|A is a random variable that follow the Poisson distribution with parameter A, abbreviated as X|A ~ Poisson (A) and the Poisson parameter A ~ MQLD (0, a, 5). Then, the pmf of the PMQLD is defined as
\5-1 , aS-1T
0 ^ r(5)r(x + 1)a3 (1 + 0)5-1 + 05-1 r(x + 5)^ x!(a3 + 1)(1 + 0)x+5 r(5)
fx(x) = —-———-„, --,x = 0,1,2,...,0 > 0,5 > 0,a3 > -1. (6)
x!(a3 1 1 W1 1 ^
Proof. Since X|A ~ Poisson (A) and A ~ MQLD (0, a, 5), the unconditional distribution of X can be obtained by substituting equations (2) and (5) in equation (1) as below
rto P-A\x 0P-0A / \
f(x) = I -->TWTrm(m" + (0A) ) dA
/0 x! (a3 + 1)r(5)
3 ^^ (T(5)a3 (~ e-A(1+0)AxdA + 05-1 H Ax+5-1 e-A(1+0)dA^ x!(a3 + 1)r(5) \ J0 J0 )
= 0 f r(5)a3r(x + 1) 0S-1r(x + 5) = x!(a3 + 1)r(5) V (1 + 0)x+1 + (1 + 0)x+5
0
x!(a3 + 1)r(5)(1 + 0 )x+5
((1 + 0)5-1r(5)a3 r(x + 1) + 05-1 r(x + 5)).
■
Remarks:
1. Equation (6) presents a two-component mixture of GD( 1+) and NBD(5, 1+) with the mixing proportion p = af+T.
2. For a ^ 0, the PMQLD reduces to the NBD(5, ^).
3. For a ^ to, the PMQLD reduces to the GD( ^).
The right tail behaviors of the PMQLD for different values of 0, a, and 5 are illustrated in Figure 1. For fixed a and 5, it is clear that the distribution's right tail approaches to zero at a faster rate when 0 increases. For fixed 0 and 5, and when a is increasing, the distribution's right tail approaches to zero at a slower rate when compared with the changes of 0. Further, for fixed 0 and a, and when 5 is increasing, the distribution captures more right tail. From Figure 2, we may note that the PMQLD may be a bimodal distribution when parameter value 5 is very different (higher value) from the parameter values of 0 and a.
The corresponding cumulative distribution function of PMQLD is given as
x
Fx (x) = E f (x) t=0
Ramajeyam Tharshan, Pushpakanthie Wijekoon RT&A, No 1 (67)
A NEW MIXED POISSON DISTRIBUTION Volume 17, March 2022
_ 5(1 + 9)5-1r(5)a3r(x + 1)((1 + 9)x+1 - 1) + 95r(x + 5 + 1)2f(1,x + 5 + 1;5 + 1; ^) = (a3 + 1)r(5)x!5(1 + 9 )x+5+1
,x = 0,1,2,...,9 > 0,5 > 0,a3 > -1,
where 2F1 (c, d; r; w) is the Gaussian hypergeometric function defined as
rt A \ ^ (c)i (d)iwi 2Fl(c,d;r;w) = g ' ,
i=0 ( )l
which is a special case of the generalized hypergeometric function given by the expression
r I \ ^ (p1 )i...(pa)iWl
aFb(p1,p2,...pa;q1,q2,...qb;w) = g (q1)i...(qfc)ii! ,
and (p)i = rTp+) = p(p + 1)...(p + i + 1) is the Pochhammer symbol.
(7)
8 = 0.75 8 = 1.25
a = 0.75 ^ 0 . a = 0.75
8 = 1.25 £ « 0 " 8 = 1.25
o
1 1 5 10 1 15 1 1 5 10 1 15
x
(a)
10
15
■
8 = 0.1 a = 0.05 8 = 2.75
x
(g)
10 20 30 40 50 60
x (h)
10 20 30 40 50 60
x (i)
Figure 1: The probability mass function of the PMQLD at different parameter values of 9, a,and ,5
10
20
30
40
50
3. Statistical properties of PMQLD
In this section, we present some important statistical properties of the PMQLD such as the shape of the distribution, moments and related measures, probability and moment generating functions, and quantile function.
e = 0.2
a = 0.75 S = 8.5
e = 0.5 a = 1.25 S = 24.5
L
0 50 100 150
T-1-1-r
0 50 100 150 0 50 100 150
Figure 2: Some b/moda/ d/str/but/ons of the PMQLD
3.1. Shape of the distribution
From (6) we can easily derive f (0) = 0((1(+30+1)(1++)5—), and limx^TO f (x) = 0.
The recurrence relation for probabilities is given by
f (x + 1) _ A
f(x) (x + 1)(1 +
x = 0,1,2,...,
where A = (1 + 0)5-1 r(5)a3r(x + 2) + 05-1 (x + 5)r(x + 5), and B = (1 + 0)5-1r(5)a3 r(x + 1) + 05-1 r(x + 5)).
The PMQLD(0, a, 5) has a log-concave probability mass function when
A^(x) = f (x+1) - f (x + 2) > 0, Vx (Gupta et al., 1997)
f (x) f (x + 1)
(x + 2) A2
(x + 1)^(1 + 0)S-1r(S)a3 r(x + 3) + 9S-1r(x + S + 2)^ B
> 1.
(8)
Under this condition, the distribution represents a unimodal distribution. Further, by using (8), it can be shown that
(i) For (1 + 0)S-1 a3 + dSS < (1 + 0)((1 + d)S-1 a3 + 9s-1 ), equation (6) has unique mode at X = 0.
(ii) Equation (6) has a unique mode at X = x0, for
(1 + 0)S-1 r(S)a3 r(x0 + 2) + 9S-1(x0 + S)r(x0 + S)
(x0 + 1)(1 + d)[ (1 + 0)S-1r(S)a3r(x0 + 1) + 9S-1r(x0 + S)
and
(1 + 9)5-1r(5)a3r(x0) + 05-1r(x0 + 5 - 1) < 1
x0 (1 + 9) + 9)5-1 r(5)a3r(x0 + 1) + 95-1 r(x0 + 5)
(iii) Equation (6) has two modes at X = x0 and X = x0 + 1, for
(1 + 9)5-1 r(5)a3r(x0 + 2) + 95-1 (x0 + 5)r(x0 + 5) = (x0 + 1)(1 + 9 )((1 + 9)5-1 r(5)a3r(x0 + 1) + 95-1 r(x0 + 5)). The above facts are also shown in Figures 1 and 2 at different parameter settings. Further, the PMQLD(9, a, 5) has a log-convex probability mass function when A^(x) < 0:
^_(x + 2)A2_ < 1.
(x + 1)^(1 + 9)5-1r(5)a3 r(x + 3) + 95-1r(x + 5 + 2)^ B
3.2. Survival and hazard rate functions
The survival/reliability function is associated with the probability of a system that will survive beyond a specified time. The survival function of the PMQLD is defined as
S(x) = 1 - F(x) = 1 - + 95r(x + 5 + 1)2' + 5 + 1;5 + 1'1++9) (9)
where ft1 = 5(1 + 9)5-1r(5)a3r(x + 1)((1 + 9)x+1 - 1) and ft2 = (a3 + 1)r(5)x!5(1 + 9)x+s+1. The hazard rate function (hrf) is the instantaneous failure rate. The hrf of the PMQLD is defined as
w N .. P(x < X < x + Ax|X > x) h(x) = limAx^0 —-
Ax
59(1 + 9) i (1 + 9)5-1 r(5)a
fM =_
S(x) ft - (ft2 + 9sr(x + 5 + 1)2F1 (1,x + 5 + 1;5 + 1; 1+9)).
59(1 + 9) ^(1 + 9 )5-1 r(5)a3r(x + 1) + 95-1 r(x + 5)^
(10)
Figure 3 provides an illustration of the possible shapes of the PMQLD's hazard rate function at different shape parameter values. According to these illustrations, it is clear that the proposed model has the capability to model the bathtub, monotonic increasing, and decreasing failure rate shapes.
3.3. Moments and related measures
The central tendency, horizontal symmetry, tail heaviness, and dispersion are important characteristics of a distribution. These characteristics can be studied by using the moments. The following theorem provides the rth factorial moment of the PMQLD.
Theorem 2. Let X ^PMQLD(9, a, 5), then the rth factorial moment of X is given as
, = r(5)r(r + 1)a3 + r(5 + r) F(r)= (a3 + 1)r(5)9r . (11)
Proof.
p'(r) = E(nr-1(X - i))
-r~
100
8 = 0.1 a = 0.5 8 = 2.5
T
100
-T~
100
8 = 0.1 a = 0.5 8 = 0.25
Figure 3: The hazard rate function at different values of a, and 5
T
100
to (r) rto e-AAx 0e-0A
xE0 x( -/0 x! (a3 + 1)r(5)
r(5)a3 + (0A)5-1 dA
/■to 0e-0A Ar-0e-
/0 A (a3 + 1)r(5)
r(5)a3+(0A) dA
■
= 0 f r(5)a3r(r + 1) r(5 + r) \ = (a3 + 1)r(5)V 0r+1 + 0r+x J
= r(5)r(r + 1)a3 + r(5 + r) = (a3 + 1)r(5)0r .
Then, the first four raw moments of X can be derived by the following relationship
r
Fr = E(xr) = ES(r,'V« ;r = 1,2,...,
i=0
where S(r, i) is the Stirling numbers of the second kind, and it is defined as
s(r,/) = e: (-1)!-,(')/r, 0 < i < r.
Let
Ki = a3 + 5, K2 = 2a3 + 5(5 + 1), K3 = 6a3 + 5(5 + 1)(5 + 2), k4 = 24a3 + 5(5 + 1)(5 + 2)(5 + 3). Then,
150
150
150
150
150
150
/ _ k1 _ i _ 9k1 + k2 / _ 92k1 + 39k2 + k3 i
p1 = (a3 + 1)9 = P, P2 = (a3 + 1)92' P3 = (a3 + 1)93 ' P4
93k1 + 792k2 + 69k3 + k4 (a3 + 1)94 .
Further, the rth-order moments about the mean can be obtained by using the relationship between moments about the mean and moments about the origin, i.e.
Pr
E
(Y - p)
i=0
g i (-1)
Pr/-l
; r = 1,2,
Therefore, the variance of X, a2 and index of dispersion, 71 are derived as
a
P2
a6(1 + 9) + a3(2 + 9 + 5(9 + 5 - 1)) + 5(9 + 1)
p + p
(a3 +1)2 92 2f a3(a3 + 2 + 5(5 - 1)) + 5 (a3 + 5)2
and
71
p2 _ a3(a3 + 2 + 5(5 - 1)) + 5
p1
(a3 + 1)(a3 + ,
respectively. It is clear that the 71 > 1. Then, the PMQLD is an over-dispersed distribution. Since the mathematical expressions of p3 and p4 are very long, we present the graphical presentations
of the skewness (72) = , P..3-„ and kurtosis (73) = ^ of the PMQLD in Figure 5. The surface
( p 2 ) 3/2 p 22 plots in Figures 4, and 5 show some possible values of the index of dispersion, skewness, and kurtosis that can be accommodated by the PMQLD at different settings of the parameters. Hence, these plots indicate that the PMQLD (9, a, 5) has the capability to accommodate various ranges of the index of dispersion, skewness, and kurtosis at different sets of parameters for over-dispersed count data.
0=0.10
0=0.50
0=1.50
10
10
10
Figure 4: Surface plots for the index of dispersion at different values of 9, a, and 5
Figure 5: Surface plots for the skewness and kurtosis functions at different values of 0, and a
3.4. Probability and moment generating functions
The characteristics of a probability distribution are directly associated with its probability generating function (pgf) and the moment generating function (mgf). The following theorem provides the pgf of the PMQLD.
Theorem 3. The pgf, G(t) = E(tX), X ~PMQLD(0,a,5) is given as
Proof.
0(a3(1 - t + 0)5-1 + 05-1) R G(t) = («3 + 1)(1 - t + 0)5 , t g R.
G(t) = E(tX)
to /-to e-Aii 0e-0X / \
x=0 'X1 it ra^K3+(0A)5-1)
r(5)a3 + (0X)5-1) dX
dX
z>to c\n—0X
r e-A(1-0. 0e
(a3 + 1)r(5)
(a3 + 1)r(5)
^f(5)a3 jf" e-X(1-t+0)dX + 05-1 jf" e-X(1-t+0)X5-1dX
(a3(1 - t + 0)5-1 + 05-1) (a3 + 1)(1 - t + 0)5
(12)
The mgf can be obtained effortlessly from pgf by using the relationship G(ef) = E(etX) = MX(t), and given as
(a3(1 - ef + 0)5-1 + 05-1)
MX(t)
t ^ A V
(a3 + 1)(1 - ef + 3.5. Quantile function
-, t g R.
(13)
The quantile function is a useful function to estimate the quantiles. Let us define the quantiles for random variable X ^PMQLD(0, a, 5). The uth quantile can be derived by solving F(xu) = u for
xu,0 < u < 1. Then, the uth quantile function of the PMQLD is given as
0
01 (xu)+ 05r(xu + 5 + 1)2F (1, Xu + 5 + 1; 5 + 1; ——) - u^u) = 0, 0 < u < 1, (14)
1+0
30
30
0
where ft(x„) = ¿(1 + 9)¿-1 T(S)a3T(xu + 1)((1 + 9)x»+1 - 1) and ft(xu) = (a3 + 1)r(¿)xu!¿(1 +
9 )Xu +¿+1.
Since equation (14) is not a closed-form in xu, the estimates of the quantiles can be evaluated by using any numerical method. Further, the first three quartiles can be calculated by substituting u = 0.25,0.50, and 0.75 in equation (14) and solving the respective equations.
4. Simulation and parameter estimation
4.1. Simulation of the random variables
Here, we provide two different algorithms to simulate the random variables X1, x2,..., xn from the PMQLD(9, a, ¿) with size n based on the inverse transform method.
The first algorithm is obtained by considering the mixing of the PMQLD. Since X|A ^Poisson (A) and A ^MQLD(9, a, ¿), the first algorithm is obtained as follows
Algorithm I:
i Simulate the random variables, u. ^uniform(0,1); i = 1,2,..., n.
ii Solve the non-linear equation for Ai: r(¿)(1 + a3(1 - e-9Ai)) - r(¿, 9Ai) - u (a3 + 1)r(¿) = 0 to simulate the random variables, Ai ^MQLD(9, a, ¿); i = 1,2,..., n.
iii Simulate xi from Poisson (Ai); i = 1,2,..., n.
The second algorithm is obtained from the quantile function of PMQLD discussed in subsection 3.5, and the steps are as follows
Algorithm II:
i Simulate the random variables, ui ^uniform(0,1); i = 1,2,..., n.
ii Solve the non-linear equation for [xu. ];
ft(xu;) + 9¿r(xu,. + ¿ + 1)2F(1,xUi + ¿ + 1;¿ + 1; T+9) - uip2(xu¡) = 0, where ft (x^) and ft(xui) are defined as in section 3.6. [.] denotes the integer part.
4.2. Parameter estimation of PMQLD
In this subsection, we discuss the parameter estimation of the PMQLD by using the method of moment estimation and the maximum likelihood estimation method.
4.2.1 Method of moment estimation (MME)
Given a random samples x\, x2...xn with size n from the PMQLD(9, a, ¿), the method of moment estimators of 9, a, and ¿, abbreviated as 9mme, &mme, and ¿mme, can be evaluated by equating
n
E <
' i=1 the raw-moments, say цr, to the sample moments, say-, r = 1,2,3, i.e. we need to find the
solutions of the following system of non-linear equations:
nx\ - (a3 + 1)9E,n=i x. = 0; n(9xt + K2) - (a3 + 1)92 E,n=i x2 = 0;
n(92K1 + 3K2 + K3) - (a3 + 1)93 E,n=1 x3 = 0,
where K1, k2, and k3 are defined in subsection 3.3. It is clear that these equations are not a closed form. However, the solutions can be derived by using a numerical method.
4.2.2 Maximum likelihood estimation (MLE)
Given a random samples Xi, X2...Xn with size n from the PMQLD(0, a, S), the likelihood function of the i'th sample value Xi is given as
L(0a,S|X) = Xi!(a3 + i)(ig+ e)*+*r(S) + 1)a3(1 + 0^ + ^ r(Xi + S)) ■
Then, the log-likelihood function is given as
log(L(0, a, SlXi )) = l (0, a, S|x)
= n (log(0) - log(a3 + 1) - log(T(S))^ + E,n=1 log^T(3)r(xi + 1)a3(1 + 0)S-1 + 0s-1 r(xi + S)j
- log(Xi!) - (Xi + S)log(1 + 0). (15)
The score functions are
dl(0,a,S|X) = n A r(S)T(Xi + 1)a3(S - 1)(1 + 0)s-2 + r(Xi + S)(S - 1)0S-2 " Xi + S d0 = 0 + i=1 r(S)T(Xi + 1)a3(1 + 0)S-1 + 0S-1 r(Xi + S) i=1 1 + 0 '
dl(0,a,S|x) _£ 3a2r(S)r(Xi + 1)(1 + 0)S-1 3na2
and
da r(¿)r(x¿ + 1)a3 (1 + 9 )ô-1 + 6s-1 r(x + ô) a3 + 1'
dl(9,a,ô|x) _
dô
^ r(xi + 1)a3(r(ô)(1 + 9)s-11og(1 + 9) + (1 + 9)s-1 r(S)^(ô)) + r(x + ô)9s-1 (log(9) + + ô)) ¿Í r(ô)r(xi + 1)a3 (1 + 9)s-1 + 9s-1 r(xi + ô)
-n(log(1 + 9 ) + ^(ô)),
d r' (a)
where = — logr(a) = r( ^. By setting the score functions equal to zero, the maximum da i(a)
likelihood estimators of 9, a, and ô, abbreviated as 9mle, amle, and ¿mle can be derived. These systems of non-linear equations can be solved by a numerical method. Here, the solutions of the parameter estimates will be obtained by using the optzm function in the R package stats.
The asymptotic confidence intervals for the parameters 9, a, and ô are derived by the asymptotic theory. The estimators are asymptotic three-variate normal with mean (9, a, ô) and the observed information matrix
( d2l(9,a,ô|x) d2l(9,a,ô|x) d2l(9,a,ô|x) \
1(9, a, ô) =
„ d92 „ d9da „ d9dô
d l(9,a,ô|x) d2l(9,a,ô|x) d2l(9,a,ô|x)
„ dadO „ da2 „ 3a35
d21 (0, a, 5|x) d21 (O, a, 5|x) d21 (O, a, 5|x)
V ^ 352 )
at O = Omle, a = a mle , and 5 = ¿mle , i.e. (Omle, a mle , 5mle ) ~ N3 ((O, a, 5), 1-1(0, a, 5)). The second order partial derivatives of the log-likelihood function are given in Appendix.
Therefore, a (1 — a)100% confidence interval for the parameters O, a, and 5 are given by
O MLE ± VA^Omle ), a MLE ± Z„/2 JV«r(a MLE ), ¿MLE ± Z«/^ V«r(5MLE ),
wherein, the V«r(OMLE), Var(aMLE), and V«r(5MLE) are the variance of Omle, ^mle, and 5mle, respectively, and can be derived by diagonal elements of 1 (O, a, 5) and za/2 is the critical value at « level of significance.
5. Monte Carlo simulation study and real-world application
This section is devoted to discuss the simulation study and the applicability of PMQLD.
5.1. Monte Carlo simulation study
Here, we examine the accuracy of the MLE method in the unknown parameter estimation of the PMQLD with respect to sample size n. The second algorithm given in subsection 4.1 is used to simulate the random variables from the PMQLD. The sample sizes are taken as 60,100,200, and 300, and the simulation study is repeated 1000 times. The study is designed as follows
(i) Simulate 1000 samples of size n.
(ii) Compute the maximum likelihood estimates for the 1000 samples, say (9, ai, 5>i), i = 1,2, ...1000.
(iii) Compute the average MLEs, biases, and mean square errors (MSEs) by using the following equations
s(n) = 14 E1=T ^, biass(n) = 14 E1=010(s - s), and MSEs(n) = 14 £=00(> - s)2, for s = 9, a, 5, and n = 60,100,200,300.
Tables 2 to 5 represent the average MLEs, biases, and MSEs (in parentheses) of 9, a, and 5 for different values of 9, a, and 5 which are 9 = 0.1, 0.3; a = 0.25, 0.50, 0.75; and 5 = 2.50, 3.50, 4.50. Note that the biases and MSEs decrease as n increases for all parameters. Then, MLE method verifies the asymptotic property for all parameter estimates, and the parameters 9, a, and 5 are over estimated. Further, while the estimation of 9 is good for small value of 9, the estimation of a doses not show a good estimation for small value of a based on average biases and MSEs. However, there is no particular pattern for estimation of 5.
Table 2: Performance ofMLE method for the PMQLD (9 = 0.10, a = 0.50,5)
n = 60 n = 100 n = 200 n = 300
MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE)
5 = 2.50 9 a 5 0.1151 0.0152(0.0014) 0.1086 0.0086(0.0007) 0.1067 0.0067(0.0002) 0.1009 0.0009(0.0001) 0.7519 0.2519(0.0691) 0.7482 0.2482(0.0648) 0.7290 0.2290(0.0557) 0.7195 0.2195(0.0494) 3.0685 0.5685(1.0107) 3.0545 0.5545 (0.9155) 2.8466 0.3466(0.3078) 2.8358 0.3358(0.2186)
5 = 3.50 9 a 5 0.1101 0.0101(0.0010) 0.1051 0.0051(0.0006) 0.1044 0.0044(0.0002) 0.1029 0.0029(0.0001) 0.6733 0.1733(0.0348) 0.6662 0.1662(0.0317) 0.6684 0.1684(0.0300) 0.6509 0.1509(0.0259) 3.8991 0.3991(0.9110) 3.8494 0.3494(0.8928) 3.6791 0.1791(0.2898) 3.6010 0.1010(0.0887)
5 = 4.50 9 a 5 0.1087 0.0087(0.0009) 0.1054 0.0054(0.0006) 0.1037 0.0037(0.0002) 0.1029 0.0021(0.0001) 0.6299 0.1299(0.0213) 0.6268 0.1268(0.0191) 0.6217 0.1217(0.0161) 0.6192 0.1192(0.0156) 4.8664 0.3664(1.3188) 4.7506 0.2506(1.1957) 4.6354 0.1354(0.3871) 4.5805 0.0805(0.0122)
Table 3: Performance ofMLE method for PMOLD(9 = 0.30, a = 0.50,5)
n = 60 n = 100 n = 200 n = 300
MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE)
5 = 2.50 9 a 5 0.3791 0.8392 3.6229 0.0790 (0.0242) 0.3392 (0.1199) 1.1229(3.5560) 0.3635 0.8243 3.4672 0.0635(0.0215) 0.3243(0.1109) 0.9672(2.9115) 0.3378 0.7956 3.1876 0.0378 (0.0075) 0.2956 (0.0933) 0.6876(1.2228) 0.3367 0.7503 3.1509 0.0367(0.0039) 0.2503(0.0853) 0.6509 (0.7488)
5 = 3.50 9 a 5 0.3839 0.7493 4.3284 0.0839 (0.0245) 0.2493 (0.0698) 0.8284(2.5256) 0.3710 0.7215 4.1190 0.0710(0.0216) 0.2215(0.0543) 0.6190(1.9316) 0.3275 0.7170 3.7680 0.0275(0.0042) 0.2170(0.0487) 0.2680(0.4910) 0.3207 0.6983 3.7393 0.0207(0.0018) 0.1983(0.0457) 0.2393(0.2547)
5 = 4.50 9 a 5 0.3638 0.6784 5.1958 0.0638 (0.0210) 0.1784 (0.0373) 0.6958(3.0589) 0.3484 0.6703 4.9240 0.0484(0.0118) 0.1703(0.0330) 0.4240(2.0855) 0.3143 0.6658 4.6154 0.0143(0.0033) 0.1658(0.0299) 0.1154(0.6430) 0.3117 0.6578 4.6713 0.0117(0.0016) 0.1578(0.0277) 0.1713(0.2667)
Table 4: Performance ofMLE method for PMOLD(9 = 0.10, a, 5 = 2.50)
n = 60 n = 100 n = 200 n = 300
MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE)
a = 0.25 9 a 5 0.1128 0.7109 3.0685 0.0128(0.0014) 0.4609(0.2184) 0.5685(1.0447) 0.1079 0.7034 3.0085 0.0079(0.0007) 0.4534(0.2088) 0.5085(0.7195) 0.1080 0.6957 2.9490 0.0080(0.0002) 0.4457(0.2000) 0.4490(0.3138) 0.1027 0.6804 2.8473 0.0027 (0.0001) 0.4304 (0.1888) 0.3473 (0.1865)
a = 0.75 9 a 5 0.1351 0.8487 3.1905 0.0351(0.0065) 0.0987(0.0167) 0.6905(2.0102) 0.1217 0.8414 3.1765 0.0217(0.0017) 0.0914(0.0135) 0.6765(1.4733) 0.1203 0.8318 3.0104 0.0203(0.0023) 0.0818(0.0095) 0.5104(0.7818) 0.1156 0.8190 2.8781 0.0156(0.0004) 0.0690(0.0081) 0.3781(0.3592)
Table 5: Performance ofMLE method for PMOLD(9 = 0.30, a, 5 = 2.50)
n = 60 n = 100 n = 200 n = 300
MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE) MLE Bias(MSE)
a = 0.25 9 a 5 0.3998 0.7823 3.5545 0.0998(0.0329) 0.5323(0.2928) 1.0545(2.7164) 0.3805 0.7795 3.4774 0.0805(0.0240) 0.5295(0.2873) 0.9774(2.2678) 0.3482 0.7611 3.1584 0.0482(0.0060) 0.5111(0.2636) 0.6584(0.7994) 0.3336 0.7437 3.0583 0.0336(0.0030) 0.4937(0.2508) 0.5583(0.5396)
a = 0.75 9 a 5 0.3773 0.8962 2.9598 0.0773(0.0621) 0.1462(0.0280) 0.4598(4.7991) 0.3566 0.8889 2.8278 0.0566(0.0506) 0.1389(0.0221) 0.3278(4.4320) 0.3275 0.8693 2.6089 0.0275(0.0267) 0.1193(0.0186) 0.1089(2.7618) 0.3097 0.8475 2.5723 0.0097(0.0225) 0.0975(0.0175) 0.0723(2.4811)
5.2. Real-world applications
In this subsection, we discuss the real-world applications of the proposed mixed Poisson distribution. Two data sets are considered to illustrate whether the proposed distribution is well fitted compared to some other existing competing Poisson mixtures. The best-fitted distribution was selected based on the negative log-likelihood (-2logL), Akaike Information Criterion (AIC), and chi-square goodness of fit statistic. The unknown parameters of the models are estimated by using the MLE method. Tables 6 and 7 summarize all these statistical measures for each data set, and the standard errors of the parameter estimates are reported in parentheses.
The first data set contains the epileptic seizure counts (Chakraborty, 2010). The sample index dispersion (t) of this data set is 1.867. Since t value is greater than one, the distribution of
this data set is clearly over-dispersed. Also, the skewness and excess kurtosis for this example are 1.239 and 1.680, respectively, which show that the distribution of the data set is positively skewed and leptokurtic. This data set was used to fit the PMQLD, GD, NBD, PLD, GPLD, PQLD, and PGLD. Table 6 presents the estimates of the parameters of distributions and the goodness of fit test. Of all eight distributions, the PMQLD performs well based on the smallest AIC value of 1191.83 and the smallest chi-square value (^2) of 2.93 (p-value=0.71).
The second data set represents the number of roots produced by 270 micro-propagated shoots of the columnar apple cultivar Trajan (Ridout et al., 1998). This is a bimodal data set for which the sample index dispersion, skewness, and excess kurtosis are 3.077, 0.182, and -1.056, respectively. These values indicate that the distribution of the data set is extremely over-dispersed, mild positively skewed, and platykurtic. This data set was also used to fit the same distributions that we used for the first example. Table 7 summarizes the results of parameter estimations and the goodness of fit test. The results show that the PMQLD having AIC=1350.20, = 11.75, p-value=0.47 outperforms clearly than other distributions.
Figure 6 illustrates how the expected values of the proposed distribution adhere with the observed value for the data sets. We can see that the observed values of the first and second data sets are very close to the expected values of the PMQLD, and the observed values of the third data set are very close to the expected values of the ZMPQLD.
Figure 6: Performance of PMQLD for the real-data sets
Table 6. Epileptic seizure counts
Counts Observed Expected
GD NBD PLD GPLD PQLD PGLD PMQLD
0 126 137.95 120.22 128.72 121.93 121.82 122.85 125.65
1 80 83.73 93.00 87.14 90.92 90.95 90.08 80.98
2 59 50.82 59.18 55.26 58.72 58.76 57.85 59.30
3 42 30.85 34.94 33.63 35.20 35.23 35.20 39.25
4 24 18.72 19.84 19.89 20.16 20.17 20.53 22.99
5 8 11.37 10.99 11.52 11.20 11.20 11.54 12.17
6 5 6.90 5.98 6.57 6.08 6.08 6.27 5.94
7 4 4.19 3.22 3.70 3.25 3.25 3.31 2.72
8 3 6.47 3.63 4.57 3.54 3.54 3.37 2.00
Total 351 351 351 351 351 351 351 351
0 = 0.65 ß = 1.00 0 = 0.97 0 = 1.11 0 = 1.12 0a = 1.57 0a = 2.70
(0.04) (0.19) (0.05) (0.13) (0.13) (0.66) (1.26)
a = 1.55 a = 2.76 a = 0.38 a = 1.49 a = 0.82
MLE (0.28) (2.76) (0.33) (0.57) (0.07)
ß = 3.89 5 = 5.89
(2.38) (2.83)
X2 11.42 5.67 5.84 4.85 4.86 4.66 2.93
p-value 0.12 0.46 0.56 0.56 0.56 0.46 0.71
-2logL 1196.79 1189.88 1190.36 1188.96 1188.96 1188.54 1185.83
AIC 1198.79 1193.88 1192.36 1192.96 1192.96 1194.54 1191.83
Table 7. Number of roots
Counts Observed Expected
GD NBD PLD GPLD PQLD PGLD PMQLD
0 64 44.62 36.87 31.09 35.46 35.45 82.81 61.93
1 10 37.25 36.05 32.94 34.00 33.99 17.81 13.92
2 13 31.09 32.16 31.79 31.19 31.19 15.31 8.47
3 15 25.95 27.77 29.06 27.78 27.77 16.46 12.85
4 21 21.66 23.58 25.63 24.20 24.20 17.53 19.30
5 18 18.08 19.83 22.04 20.74 20.74 17.70 24.53
6 24 15.09 16.56 18.60 17.55 17.55 16.95 26.85
7 21 12.60 13.76 15.48 14.69 14.69 15.53 25.94
8 23 10.52 11.39 12.73 12.20 12.20 13.70 22.55
9 21 8.78 9.40 10.37 10.05 10.05 11.72 17.90
10 17 7.33 7.74 8.39 8.23 8.24 9.76 13.13
11 12 6.12 6.37 6.74 6.71 6.71 7.96 8.99
12 5 5.11 5.23 5.38 5.44 5.45 6.36 5.78
13 2 4.26 4.28 4.27 4.40 4.40 5.01 3.52
14 3 3.56 3.51 3.38 3.54 3.54 3.88 2.03
> 15 1 17.98 15.50 12.11 13.82 13.83 11.51 2.31
Total 270 270 270 270 270 270 270 270
0 = 0.12 ß = 0.24 0 = 0.35 0 = 0.37 0a = 0.32 0a = 0.59 0a = 4.23
(0.01) (0.03) (0.02) (0.03) (0.03) (0.08) (1.36)
a = 1.21 a = 0.47 a = 0.68 a = 0.22 a = 0.73
MLE (0.16) (0.26) (0.27) (0.12) (0.04)
ß = 4.25 5a = 29.34
(0.66) (9.28)
X2 121.92 120.76 117.44 110.58 110.58 46.72 11.74
p-value 0.00 0.00 0.00 0.00 0.00 0.00 0.466
-2logL 1464.90 1462.63 1454.10 1451.45 1451.45 1384.21 1344.20
AIC 1466.90 1466.63 1456.10 1455.45 1455.45 1390.21 1350.20
6. Conclusion
This paper proposes an alternative mixed Poisson distribution to model the over-dispersed count data. Explicit expressions of the pmf, hazard rate function, moments, mean, variance, skewness, and kurtosis were derived for the proposed distribution. Its pmf possesses to be either unimodal or bimodal, and hazard rate function presents monotonic increasing, decreasing, and bathtub shapes. The kurtosis and the variance-to-mean ratio functions of the new distribution indicate that
the distribution can capture various ranges of right tail weights as well as the index of dispersions. Further, its structural properties show that the new distribution is much more flexible than its predecessors, negative binomial, geometric, and Poisson-Lindley distributions. The maximum likelihood method was employed to estimate the parameters of the distribution, and the observed information matrix has also been derived. The proposed distribution and some other competing Poisson mixtures have been fitted to two real-world data sets. The results show that the proposed distribution could provide a better fit than a set of common Poisson mixtures considered in these applications.
Acknowledgement
We thank the Postgraduate Institute of Science, University of Peradeniya, Sri Lanka for providing all facilities to do this research.
References
[1] Albrecht, P. (1984). Laplace Transforms, Mellin Transforms and Mixed Poisson Processes. Scandinavian Actuarial Journal, 11: 58-64. https://doi.org/10.1080/03461238.1984.10413753
[2] Bhati, D., Sastry, D. V. S., and Qadri, P. Z. (2015). A New Generalized Poisson Lind-ley Distribution, applications and Properties. Austrian Journal of Statistics, 44: 35-51. https://doi.org/10.17713/ajs.v44i4.54
[3] Bhattacharya, S. K. (1967). A Result in Accident Proneness. Biometrika, 54: 324-325. https://doi.org/10.1093/biomet/54.1-2.324
[4] Cameron, A. C. and Trivedi, P. K. Regression Analysis of Count Data, (2nd Ed.) Cambridge University Press, New York, USA, 2013. https://doi.org/10.1017/CBO9781139013567
[5] Chakraborty, S. (2010). On Some Distributional Properties of the Family of Weighted Generalized Poisson Distribution. Communications in Statistics - Theory and Methods, 39: 2767-2788. https://doi.org/10.1080/03610920903129141
[6] Elbatal, I., Merovci, F., and Elgarhy, M. (2013). A new generalized Lindley distribution, Mathematical Theory and Modeling, 3: 30-47.
[7] Feller, W. (1943). On a Generalized Class of Contagious Distributions. Annals of Mathematical Statistics, 14: 389-400. https://www.jstor.org/stable/2235926
[8] Greenwood, M. and Yule, G. U. (1920). An inquiry into the nature of frequency distributions
representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents, Journal of the Royal Statistical Society, 83: 255-279. https://www.jstor.org/stable/2341080
[9] Grine, R. and Zeghdoudi, H. (2017). On Poisson quasi-Lindley distribution and its applications. J. of Modern Applied Statistical Methods, 16: 403-417.
[10] Irwin, J. (1975). The Generalized Waring Distribution Parts I, II, III. Journal of the Royal Statistical Society A, 138: 18-31 (Part I), 204-227 (Part II), 374-384 (Part III).
[11] Johnson, N. L., Kotz, S., and Kemp, A. W. Univariate Discrete Distributions, (2nd Ed.) John Wiley and Sons, New York, USA, 1992.
[12] Lindley, D. V. (1958). Fiducial Distributions and Bayes' Theorem. Journal of the Royal Statistical Society, 20: 102-107. https://www.jstor.org/stable/2983909
[13] Lynch, J. (1988). Mixtures, Generalized Convexity and Balayages, Scandinavian Journal of Statistics, 15: 203-210.
[14] McCullagh, P. and Nelder, J. Generalized Linear Models, (2nd Ed.) Chapman and Hall, New York, USA, 1989.
[15] Ong, S.H. and Muthaloo, S. (1995). A Class of Discrete Distributions Suited to Fitting Very Long Tailed Data. Communication in Statistics-Simulation and Computation, 24: 929-945. https://doi.org/10.1080/03610919508813285
[16] Ridout, M., Demktrio, C. G. B, and Hinde, J. (1998). Models for count data with many zeros, Invited paper presented at the Nineteenth International Biometric Conference, Capetown, South Africa.
[17] Ruohonen, M. (1988). A Model for the Claim Number Process. ASTIN Bulletin, 18: 57-68. https://doi.org/10.2143/AST.18.L2014960
[18] Sankaran, M. (1970). The discrete Poisson-Lindley distribution, Biometrics, 26: 145-149. https://www.jstor.org/stable/2529053
[19] Shaked, M. (1980). On Mixtures from Exponential Families, Journal of the Royal Statistical Society B, 42: 192-198. https://doi.org/10.1111/j.2517-6161.1980.tb01118.x
[20] Shanker, R. and Mishra, A. (2013a). A Quasi Lindley Distribution. African Journal of Mathematics and Computer Science Research, 6: 64-71.
[21] Shanker, R., Sharma, S., and Shanker, R. (2013b). A Two-Parameter Lindley Distribution for Modeling Waiting and Survival Times Data. Applied Mathematics, 4: 363-368.
[22] Tharshan, R. and Wijekoon, P. (2020a). Location Based Generalized Akash Distribution: Properties and Applications, Open Journal of Statistics, 10: 163-187.
[23] Tharshan, R. and Wijekoon, P. (2020b). A comparison study on a new five-parameter generalized Lindley distribution with its sub-models. Statistics in Transition new Series, 21: 89-117. https://doi.org/10.21307/stattrans-2020-015
[24] Tharshan, R. and Wijekoon, P. (2021). A modification of the Quasi Lindley distribution. Open Journal of Statistics, 11: 369-392.
[25] Willmot, G. E. (1993). On Recursive Evaluation of Mixed Poisson Probabilities and Related Quantities, Scandinavian Actuarial Journal, 18: 114-133.
[26] Wongrin, W. and Bodhisuwan, W. (2016). The Poisson-generalised Lindley distribution and its applications, Songklankarin J.Sci. Technol, 38: 645-66.
Appendix: Elements of the observed information matrix, i(9, a, 5) defined in subsection 4.2.1:
Let us define Ti, T2, T3, T4, T5, T6, T7, Tg, T9, T10, Tu, T12, T13 and T14 as follows T1 = r(5)r(x + 1)a3(5 - 1)(1 + 9)5-2 + r(x + 5)(5 - 1)95-2, T2 = r(5)r(x + 1)a3(1 + 9)5-1 + 95-1 r(xf + 5), T3 = 3a2r(5)r(xi + 1)(1 + 9 )5-1, T4 =
r(xi + 1)a3(r(5)(1 + 9)5-11og(1 + 9) + (1 + 9 )5-1 r (5)^(5)) + r(x + 5)95-1 (log(9) + ^ + 5)), T5 = r(5)r(xi + 1)a3(1 + 9)5-3 + r(xi + 5)95-3, T6 = r(5)r(xi + 1)a3(1 + 9)5-2 + r(xi + 5)95-2, T7 = 1og(1 + 9)(r(5)(1 + 9 )5-11og(1 + 9) + (1 + 9)5-1r (5)^(5)) Tg = (1 + 9)5-1 (r(5)^(5) + (^(5))2 r(5)) + r(5)^(5)(1 + 9 )5-11og(1 + 9), T9 = r(xi + 5)95-11og(9 ) + 95-1r(xi + 5)ty(%i + 5), T10 = r(xi + 5)93-1^1(xi + 5), T11 = r(xi + 1)a3(r(5)(1 + 9)5-11og(1 + 9) + (1 + 9 )5-1 r(5)tf(5)), T12 = 95-1 r(xi + 5)^(xi + 5) + r(xi + 5)95-11og(9), T13 = r(5)((1 + 9)5-2 + 1og(1 + 9)(5 - 1)(1 + 9)5-2) + (5 - 1)(1 + 9)5-2r(5)^(5),
and
T14 = (1og(9) + 0(xi + 5))r(xi + 5)(5 - 1)95-2.
Then, the second order partial derivatives of the log-likelihood function are as follows
d21(0,a, S|x) _ -n E xi + 8 E T2(8 - 1)(S - 2)Ts - Ti(8 - 1)T6 d02 _ 02 + ¿1 (1 + 0)2 + ¿1 T22 ,
da2
T2(6aT(S)T(xi + 1)(1 + 0)8-1) - T3(3a2r(S)r(xi + 1)(1 + 0)8-1) 3na(2(a2 + 1) - 3a3)
d2l(0, a, S|x)
n T./' £a,W8\i 1W1 I 0\S-1\ T.
¿1 T22 (a3 +1)2 '
d2l(0,a,S|x) _ dS2 _
E T2(r(Xi + 1)a3(T7 + T8) + (log(0) + + S))T9 + T10) - T4(Tn + T12) - ^(S) i_1 T22 1
d2l(0,a,S|x) _ E T2(3a2r(S)r(xi + 1)(S - 1)(1 + 0)8-1) - T1T3
^ _E T2 ,
i_1 1
d2l(0,a,S|x) _¿A T2(3a2r(xi + 1)(r(S)(1 + 0)S-1log(1 + 0) + (1 + 0)8-1 T(S)$(S))) - T4T3
dSda ^ T2
i_1 2
and
d21(0,a,S|x) _ " T2(r(xi + 1)a3T13 + r(xi + S)08-2 + T14) - T4T1 n
E
dSd0 i_1 T22 1 + 0'
d2 n 1
where ^ (s) is the trigamma function and defined as ^1(s) _ —-2log(r(s)) _ E 7-T\2.
ds i_1 (s + k)
