The COM-negative binomial distribution: modeling overdispersion and ultrahigh zero-inflated count data

TL;DR

This paper introduces the COM-negative binomial distribution, a flexible model for overdispersed and zero-inflated count data, with theoretical properties and practical applications demonstrated through parameter estimation and goodness-of-fit tests.

Contribution

It develops a new three-parameter COM-type negative binomial distribution with comprehensive theoretical properties and applies it to real-world overdispersed, zero-inflated datasets.

Findings

Distribution models overdispersion and zero-inflation effectively.

Maximum likelihood estimation achieves good parameter fit.

Distribution shows strong goodness-of-fit in applied datasets.

Abstract

In this paper, we focus on the COM-type negative binomial distribution with three parameters, which belongs to COM-type $(a,b,0)$ class distributions and family of equilibrium distributions of arbitrary birth-death process. Besides, we show abundant distributional properties such as overdispersion and underdispersion, log-concavity, log-convexity (infinite divisibility), pseudo compound Poisson, stochastic ordering and asymptotic approximation. Some characterizations including sum of equicorrelated geometrically distributed random variables, conditional distribution, limit distribution of COM-negative hypergeometric distribution, and Stein's identity are given for theoretical properties. COM-negative binomial distribution was applied to overdispersion and ultrahigh zero-inflated data sets. With the aid of ratio regression, we employ maximum likelihood method to estimate the parameters and the goodness-of-fit are evaluated by the discrete Kolmogorov-Smirnov test.

Tables

Table 1. Table 1: Fitting results of simulated data from NB ( 1 , 0.5 ) 1 0.5 (1,0.5)

Sample values	Frequency	Fitted Values
		CMNB	NB	CMP
0	5060	5054	5049	5019
1	2480	2494	2492	2502
2	1199	1237	1236	1248
3	638	615	614	622
4	318	306	306	310
5	165	152	152	155
6	74	76	76	77
7	33	38	38	38
8	20	19	19	19
9	8	9	9	10
10	4	5	5	5
11	1	2	2	2
Total	10000	10007	9998	10007
par1		0.97	0.99	0.50
par2		1.00	0.50	0.00
par3		0.51
${\chi }^2$		5.29	5.44	5.83
d.f. of ${\chi }^2$		8.00	9.00	9.00
p value of ${\chi }^2$		0.73	0.79	0.76
K-S		0.002199	0.003976	0.004451
p value of K-S		1.000000	0.997435	0.988831

Table 2. Table 2: Fitting results of simulated data from CMNB(0.005, 0.1 0.5)

Sample values	Frequency	Fitted Values
		CMNB	NB	CMP
0	6442	6449	6421	5937
1	1866	1878	1951	2413
2	874	869	837	981
3	435	415	394	399
4	188	201	194	162
5	101	98	98	66
6	55	48	50	27
7	19	24	26	11
8	7	12	14	4
9	8	6	7	2
10	2	3	4	1
11	2	1	2	0
12	1	1	1	0
Total	10000	10005	9999	10003
par1		0.01	0.55	0.41
par2		0.12	0.45	0.00
par3		0.50		0.44
${\chi }^2$		7.26	16.38	281.94
d.f. of ${\chi }^2$		9.00	10.00	10.00
p value of ${\chi }^2$		0.61	0.09	0.00
K-S		0.001484	0.006484	0.050678
p value of K-S		1.000000	0.794579	0.000000

Table 3. Table 3: Fit of Willmot2 data

No. of claims	Frequency	Fitted Values
		CMNB	NB	CMP
0	3719	3720	3719	3681
1	232	231	230	294
2	38	39	40	23
3	7	8	8	2
4	3	2	2	0
5	1	1	0	0
Total	4000	4001	3999	3991
par1		0.57	0.22	0.08
par2		3.06	0.71	0.00
par3		0.35
${\chi }^2$		1.01	1.56	173.28
d.f. of ${\chi }^2$		2.00	3.00	3.00
p value of ${\chi }^2$		0.60	0.67	0.00
K-S		0.000250	0.000500	0.009500
p value of K-S		1.000000	1.000000	0.863178

Table 4. Table 4: Fit of Car insurance claim data

No. of Claims	Frequency	Fitted Values
		CMNB	NB	CMP
0	27141	27177	27166	26599
1	5789	5666	5664	6430
2	1443	1554	1563	1554
3	457	466	467	376
4	155	146	145	91
5	56	47	46	22
6	27	15	15	5
7	2	5	5	1
8	1	2	2	0
9	1	1	1	0
Total	35072	35079	35074	35078
par1		0.95	0.61	0.24
par2		10.40	0.66	0.00
par3		0.36
${\chi }^2$		24.16	27.88	300.60
d.f. of ${\chi }^2$		6.00	7.00	7.00
p value of ${\chi }^2$		0.00	0.00	0.00
K-S		0.002667	0.002905	0.015584
p value of K-S		0.964199	0.928729	0.000000

Equations

$$\mathrm{P}(X=k)=\frac{{\lambda }^k}{{(k!)}^{\nu }}\cdot \frac{1}{Z(\lambda ,\nu )},(k=0,1,2,\cdots ),$$

$$\mathrm{P}(X=k)=\frac{\Gamma (r+k)}{k!\Gamma (r)}{p}^k(1-p{)}^r,(k=0,1,2,\dots ),$$

$$\mathrm{P}(X=k)=\frac{{(\frac{\Gamma (r+k)}{k!\Gamma (r)})}^{\nu }{p}^k{(1-p)}^r}{\sum _{i=0}^{\infty }{(\frac{\Gamma (r+i)}{i!\Gamma (r)})}^{\nu }{p}^i{(1-p)}^r}={\left(\frac{\Gamma (r+k)}{k!\Gamma (r)}\right)}^{\nu }{p}^k{(1-p)}^r\frac{1}{C(r,\nu ,p)},(k=0,1,2,\dots ),$$

$$\mathrm{P}(X=k)=\frac{{[\frac{\Gamma (r+k)}{k!\Gamma (r)}{\tilde{p}}^k{(1-\tilde{p})}^r]}^{\nu }}{\sum _{i=0}^{\infty }[{\frac{\Gamma (r+i)}{i!\Gamma (r)}}^{\nu }{\tilde{p}}^i{(1-\tilde{p})}^r{]}^{\nu }}={\left(\frac{\Gamma (r+k)}{k!\Gamma (r)}\right)}^{\nu }{\tilde{p}}^k(1-\tilde{p}{)}^r\frac{1}{C(r,\nu ,\tilde{p})},(k=0,1,2,\dots ),$$

$$\mathrm{P}(X=k)=\frac{\Gamma {(r+k)}^{\beta }}{{(k!)}^{\alpha }}{\theta }^k/\sum _{i=1}^{\infty }\frac{\Gamma {(r+i)}^{\beta }}{{(i!)}^{\alpha }}{\theta }^i(k=0,1,2,\dots ),$$

$$\frac{\mathrm{P}(X=k)}{\mathrm{P}(X=k-1)}=p{\left(\frac{\Gamma (k+r)}{k!\Gamma (r)}\right)}^{\nu }\text{raise}3.01385pt/\text{lower}3.01385pt{\left(\frac{\Gamma (k-1+r)}{(k-1)!\Gamma (r)}\right)}^{\nu }=p\cdot {\left(\frac{k-1+r}{k}\right)}^{\nu }.$$

$$p_k=(a+\frac{b}{k}{)}^{\nu }{p}_{k-1},(k=1,2,\cdots )$$

$$\frac{p_k}{p_{k-1}}=p{\left(1+\frac{r-1}{k}\right)}^{\nu }={\left(p^{1/\nu }+\frac{(r-1)p^{1/\nu }}{k}\right)}^{\nu }.$$

$$\mathrm{P}(X=k)=\frac{{\left(\genfrac{}{}{0ex}{}{m}{k}\right)}^v{p}^k{(1-p)}^{m-k}}{\sum _{i=0}^m{\left(\genfrac{}{}{0ex}{}{m}{i}\right)}^v{p}^i{(1-p)}^{m-i}},k=0,1,\cdots ,m,$$

$$\mathrm{P}(X=k)=\frac{{\left(\genfrac{}{}{0ex}{}{m}{k}\right)}^2{p}^k{(1-p)}^{m-k}}{\sum _{i=0}^{\infty }{\left(\genfrac{}{}{0ex}{}{m}{i}\right)}^2{p}^i{(1-p)}^{m-i}},k=0,1,2,\cdots ,$$

$$\frac{\mathrm{P}(X=k)}{\mathrm{P}(X=k-1)}=\frac{{\lambda }_{k-1}}{{\mu }_k},(k=1,2,\dots ).$$

$$H_{\alpha }^R(X)=\frac{1}{1-\alpha }\ln \sum _{i=0}^{\infty }[P(X=i){]}^{\alpha },(\alpha \ne 1).$$

$$H_{\alpha }^T(X)=\frac{1}{1-\alpha }\left(\sum _{i=0}^{\infty }[P(X=i){]}^{\alpha }-1\right),(\alpha \ne 1).$$

$$\frac{p_{k+1}{p}_{k-1}}{p_k^2}=\frac{p_{k+1}}{p_k}/\frac{p_k}{p_{k-1}}\le (\ge )1,k\ge 1.$$

$$M=\frac{p_{k+1}}{p_k}/\frac{p_k}{p_{k-1}}=p{\left(\frac{r+k}{k+1}\right)}^{\nu }/p{\left(\frac{r+k-1}{k}\right)}^{\nu }={\left(\frac{k^2+kr}{k^2+kr+r-1}\right)}^{\nu }.$$

$$G(z)=\sum _{k=0}^{\infty }{p}_k{z}^k=e^{\sum _{i=1}^{\infty }{\alpha }_i\lambda (z^i-1)},(\left|z\right|\le 1),$$

$$G(z)=\sum _{k=0}^{\infty }{p}_k{z}^k=\exp \left\{\sum _{i=1}^{\infty }{\alpha }_i\lambda (z^i-1)\right\},$$

$$\left|(1-z)G(z)\right|$$

$$\ge p_0-\left|(p_0-p_1)|z|+(p_1-p_2)|z^2|+\cdots \right|$$

$$>p_0-\left|(p_0-p_1)+(p_1-p_2)+\cdots \right|=p_0-|p_0|=0.$$

$$P_{n+1}=\frac{\lambda }{n+1}[{\alpha }_1{P}_n+2{\alpha }_2{P}_{n-1}+\cdots +(n+1){\alpha }_{n+1}{P}_0],(P_0=e^{-\lambda },n=0,1,\cdots )$$

$$(\frac{\Gamma (r+n+1)}{(n+1)!}{)}^{\nu }{p}^{n+1}=\frac{\lambda }{n+1}[{\alpha }_1(\frac{\Gamma (r+n)}{n!}{)}^{\nu }{p}^n+2{\alpha }_2(\frac{\Gamma (r+n-1)}{(n-1)!}{)}^{\nu }{p}^{n-1}+\cdots +(n+1){\alpha }_{n+1}],(n=0,1,\cdots ),$$

$$\mathrm{P}(X=k)=\frac{w(k)}{\mathrm{E}[w(X)]}\cdot \frac{{\theta }^k}{k!}{e}^{-\theta },$$

$$\mathrm{P}(X=k)=\frac{{(1-p)}^r{e}^p}{C(r,\nu ,p)}[\Gamma (1+k){]}^{1-\nu }[\Gamma (r+k){]}^{\nu }\frac{p^k}{k!}{e}^{-p},(k=0,1,2,\dots ).$$

$$f(k,r,\nu )=w(k)=[\Gamma (1+k){]}^{1-\nu }[\Gamma (r+k){]}^{\nu }.$$

$$\frac{d^2\log f(k,r,\nu )}{d{k}^2}=\frac{(1-\nu )d^2\log \Gamma (k+1)}{d{k}^2}+\frac{\nu d^2\log \Gamma (k+r)}{d{k}^2}=\sum _{i=0}^{\infty }(\frac{1-\nu }{{(i+k+1)}^2}+\frac{\nu }{{(i+k+r)}^2}),(\forall k\in \mathrm{N}).$$

$$\frac{\mathrm{P}(Y=n)}{\mathrm{P}(X=n)}={\left(\frac{\Gamma (r+n)}{n!\Gamma (r)}\right)}^{{\nu }_2-{\nu }_1}\frac{C(r,{\nu }_1,p)}{C(r,{\nu }_2,p)},(n=0,1,2,\dots ).$$

$$\underset{r\to \infty }{\lim }\mathrm{P}(X=k\mid r,\nu ,p)=\frac{{\lambda }^k}{{(k!)}^{\nu }}\cdot \frac{1}{Z(\lambda ,\nu )},(k=0,1,2,\cdots ).$$

$$\mathrm{P}(X=k\mid r,\nu ,p)=\frac{{\lambda }^k}{{(k!)}^{\nu }}\cdot {\left(\frac{\Gamma (r+k)}{\Gamma (r)r^k}\right)}^{\nu }\cdot \frac{1}{{(1+\lambda /r^{\nu })}^k}/\frac{C(r,\nu ,p)}{(\frac{r^{\nu }}{r^{\nu }+\lambda }{)}^r}.$$

$$\underset{r\to \infty }{\lim }\mathrm{P}(X=k\mid r,\nu ,p)=\frac{{\lambda }^k}{{(k!)}^{\nu }}/\frac{C(r,\nu ,p)}{(\frac{r^{\nu }}{r^{\nu }+\lambda }{)}^r}=\frac{{\lambda }^k}{{(k!)}^{\nu }}\cdot \frac{1}{Z(\lambda ,\nu )}$$

$$\underset{r\to +\infty }{\lim }\frac{C(r,\nu ,p)}{(\frac{r^v}{r^v+\lambda }{)}^r}=\underset{n\to +\infty }{\lim }\underset{r\to +\infty }{\lim }\sum _{i=0}^n\frac{{\lambda }^i}{{(i!)}^{\nu }}\cdot {\left(\frac{\Gamma (r+i)}{\Gamma (r)r^i}\right)}^{\nu }\cdot \frac{1}{{(1+\lambda /r^{\nu })}^i}=\underset{n\to +\infty }{\lim }\sum _{i=0}^n\frac{{\lambda }^i}{{(i!)}^{\nu }}=Z(\lambda ,\nu ).$$

Full text

††journal: Frontiers of Mathematics in China

1 Introduction

Before 2005, Conway-Maxwell-Poisson distribution (denoted as COM-Poisson distribution) had been rarely used since Conway and Maxwell (1962) briefly introduced it for modeling of queuing systems with state-dependent service time, see also Wimmer and Altmann (1999), Wimmer et al. (1995). About ten years ago, the COM-Poisson distribution with two parameters was revived by Shmueli et al. (2005) as a generalization of Poisson distribution. More recently, there has been a fast growth of researches on COM-Poisson distribution in terms of related statistical theory and applied methodology, see Sellers et al. (2012) and the references therein. The probability mass function (p.m.f.) of the COM-Poisson random variable (r. v) $X$ is given by

[TABLE]

where $\lambda,\nu>0$ and $Z(\lambda,\nu)=\sum\limits_{i=0}^{\infty}{\frac{{{\lambda^{i}}}}{{{{(i!)}^{\nu}}}}}$. We denote (1) as $X\sim{\rm{CMP}}(\lambda,\nu)$.

Kokonendji et al. (2008) proved that COM-Poisson distribution was overdispersed when $\nu\in{\rm{[0,1)}}$ and underdispersed when $\nu\in{\rm{(1,+}}\infty{\rm{)}}$. Another extension of Poisson is negative binomial, which is a noted discrete distribution with overdispersion property and is widely applied in actuarial sciences(see Denuit et al. (2007), Kaas (2008)). The p.m.f. of the negative binomial r.v $X$ is

[TABLE]

where $r\in(0,\infty)$ and $p\in(0,1)$.

In this paper, we propose a COM-negative binomial (denoted by CMNB) distribution , which extends the negative binomial distribution and depends on three parameters by replacing $\frac{{\Gamma(r+k)}}{{k!{\mkern 1.0mu}\Gamma(r)}}$ in (2) with $({\frac{{\Gamma(r+k)}}{{k!{\mkern 1.0mu}\Gamma(r)}})^{\nu}}$ and divide the normalization constant $C(r,\nu,p)=\sum\limits_{i=0}^{\infty}{(\frac{{\Gamma(r+i)}}{{i!{\mkern 1.0mu}\Gamma(r)}})}^{\nu}{p^{i}}{(1-p)^{r}}$.

Definition 1.1**.**

A r.v. $X$ is said to follow CMNB distribution $({\rm{CMNB}}(r,\nu,p))$ with three parameters $(r,\nu,p)$ if the p.m.f. is given by

[TABLE]

where $r,\nu\in(0,\infty)$ and $p\in(0,1)$.

In another point of view, the alternative form of (3) can be written as

[TABLE]

where $\tilde{p}={p^{1/\nu}}$.

When $r\leq 1$, we will show that CMNB (3) is discrete compound Poisson, which has wide application in risk theory (includes non-life insurance) as well, see Zhang et al. (2014) and the references therein. For illustrating the p.m.f of $X$ defined in (3), we plot 12 cases of CMNB in Figture 1.

It is easy to see that our CMNB distribution belongs to the COM-type extension of $(a,b,0)$ class. The $(a,b,0)$ class distribution is a famous family of distributions which sometimes refers to Katz class (see remarks in section 2.3.1 of Johnson et al. (2005)). It has significant applications in non-life insurance mathematics, especially for modelling claim counts ( loss models, collective risk models), see Denuit et al. (2007). A classic result in non-life insurance textbooks states that the $(a,b,0)$ class distribution only contains degenerate, binomial, Poisson and the negative binomial distribution. After adding a new parameter $\nu\in{{\rm{R}}^{+}}$, we define the COM-type extension of $(a,b,0)$ distribution, and it is convenient to see that degenerate, COM-Poisson, COM-binomial and CMNB belong to this class of distributions.

Shmueli et al. (2005) firstly proposed the COM-binomial distribution which is presented as a sum of equicorrelated Bernoulli variables. Borges et al. (2014) studied some properties and an asymptotic approximation (e.g. COM-binomial approximates to COM-Poisson under some conditions) of this family of distributions in detail. We will show that some results of COM-Poisson can be extended in our CMNB distribution. Kadane (2016) gave the exchangeably properties, sufficient statistics and multivariate extension of COM-binomial distribution. Another variant of CMNB distribution has been studied by Imoto (2014), it just replaced the term ${\Gamma(r+k)}$ in (2) by ${\Gamma(r+k)}^{\nu}$ and then divided the normalization constant. We will give adequate reasons to support our extension in succeeding sections. Chakraborty and Imoto (2016) considered the extended COM-Poisson distribution (ECOMP$(r,\theta,\alpha,\beta)$):

[TABLE]

where the parameter space is $(r\geq 0,\theta>0,\alpha>\beta)\cup(r>0,0<\theta<1,\alpha=\beta)$. ECOMP distribution combines Imoto (2014)’s extension ($\alpha=1$), Chakraborty and Ong (2016)’s extension (COMNB$(r,\theta,\alpha,\beta)$ with $\beta=1$) and the CMNB distribution (3) with $\alpha=\beta=\nu$. The COM-Poisson is a special case of ECOMP when $\beta=0$. ECOMP distribution (5) has the queuing systems characterization (birth-death process with arrival rate ${{\lambda_{k}}={{(r+k)}^{\nu}}}$ and service rate ${{\mu_{k}}={k^{\nu}}\mu}$ for $k\geq 1$ ), see also Brown and Xia (2001) for arbitrary birth-death process characterization.

The rest of the article is organized as the follows. In section 2, we propose the COM-type $(a,b,0)$ class distributions and demonstrate some example of COM-type $(a,b,0)$ class which includes the CMNB. Further more, some properties of CMNB, such as Renyi entropy and Tsallis entropy representation, overdispersion and underdispersion, log-concavity, log-convexity (infinite divisibility), pseudo compound Poisson, stochastic ordering and asymptotic approximation are studied. In section 3, some conditional distribution characterizations and Stein identity characterization are presented by using related lemmas, and we also show that COM-negative hypergeometric can approximate to CMNB. In section 4, inverse methods were introduced to generate CMNB distributed random variables. Section 5 then estimates the parameters by maximum likelihood method, in which the initial values are provided by recursive formula. In section 6, two simulated data sets and two applications to actuarial claim data sets are given as examples. In section 7, we provide some potential and further research suggestions based on the properties and characterizations of CMNB distribution.

2 Properties

2.1 Recursive formula and ultrahigh zero-inflated property

The recursive formula (or ratios of consecutive probabilities) is given by

[TABLE]

We say that a zero-inflated count data $X$ following some discrete distribution is ultrahigh zero-inflated if $\mathrm{P}(X=0)/\mathrm{P}(X=1)\gg 1$. For CMNB case with two parameters, we have $\frac{\text{P}_{\text{CMNB}}(X=0)}{\text{P}_{\text{CMNB}}(X=1)}=\frac{1}{p}(\frac{1}{r})^{\nu}$; and for negative binomial case, we get $\frac{\text{P}_{\text{NB}}(X=0)}{\text{P}_{\text{NB}}(X=1)}=\frac{1}{{pr}}.$ If we choose ${r<1}$ and ${\nu>1}$ in CMNB case, then this CMNB distribution is more flexible to deal with the ratio $\frac{{\mathrm{P}(X=0)}}{{\mathrm{P}(X=1)}}$ comparing to the NB distribution, as $\frac{1}{p}{\left({\frac{1}{r}}\right)^{\nu}}\gg\frac{1}{{pr}}$ when $\nu>1$. (For examples, the plots of CMNB distribution with ${\nu=2}$ in Figure 1; $\nu=10.4$ in Table 4 of section 6.2). In the insurance company, the more zero insurance claims, the less risk to bankrupt. The following definition provides a generalization of $(a,b,0)$ class distribution.

Definition 2.1**.**

Let $X$ be a discrete r.v., if ${p_{k}}=\mathrm{P}(X=k)$ satisfies the recursive formula

[TABLE]

for some constants $a,b\in{\rm{R}}$ and $\nu\in{{\rm{R}}^{+}}$, then we call it COM-type $(a,b,0)$ class distribution. We denote this class as $COM(a,b,\nu,0)$ .

The COM-Poisson distribution ${\rm{CMP}}(\lambda,\nu)$ satisfies the case $a=0$, since ${p_{k}}=\frac{\lambda}{{{k^{\nu}}}}{p_{k-1}}$ and $b={\lambda^{1/\nu}}$.

From (6), it is easy to see that CMNB distribution belongs to COM-type $(a,b,0)$ class distribution with ${p_{0}}=1/C(r,\nu,p)$.

[TABLE]

The COM-binomial distribution (CMB), see Shmueli et al. (2005), Borges et al. (2014), with p.m.f.

[TABLE]

where $\nu\in{{\rm{R}}^{+}},m\in{{\rm{Z}}^{+}},p\in(0,1)$. We denote (9) as $X\sim{\rm{CMB}}(m,p,\nu)$. Since the ratio of consecutive probabilities is $\frac{{\mathrm{P}(X=k)}}{{\mathrm{P}(X=k-1)}}=\frac{p}{{1-p}}{\left({\frac{{m+1-k}}{k}}\right)^{\nu}}$, COM-binomial distribution belongs to COM type $(a,b,0)$ class distribution.

Remark 1: As we know, the $(a,b,0)$ class distribution only contains degenerate distribution, binomial distribution, Poisson distribution and the negative binomial distribution. But there are other distributions belongs to COM type $(a,b,0)$ class. For example, $m\in{{\rm{Z}}^{+}}$ can be replaced by $m\in{{\rm{R}}^{+}}$ in (9), and the p.m.f. is given by

[TABLE]

where $p/(1-p)<1$ such that $\sum\limits_{i=0}^{\infty}{{{{m}\choose{i}}^{2}}{p^{i}}{{(1-p)}^{m-i}}}<\infty$.

Remark 2: Brown and Xia (2001) considered the very large class of stationary distribution of birth-death process with arrival rate ${\lambda_{k}}$ and service rate ${\mu_{k}}$ by the recursive formula:

[TABLE]

Thus we can construct a birth-death process with arrival rate ${\lambda_{k}}=c{[a(k+1)+b]^{\nu}}$ and service rate ${\mu_{k}}=c{k^{\nu}}$, and $c$ is a positive constant, which characterizes the COM-type $(a,b,0)$ class distribution.

2.2 Related to Rényi entropy and Tsallis entropy

Notice the Rényi entropy (see Rényi (1961)) in the information theory, which generalizes the Shannon entropy. The Renyi entropy of order $\alpha$ of a discrete r.v. $X$:

[TABLE]

Let $X$ be negative binomial distributed in (2) and $X_{\nu}$ be CMNB distributed in (4). Then the normalization constant $C(r,\nu,\tilde{p})$ in (4) has Rényi entropy representation $C(r,\nu,\tilde{p})={e^{(1-\nu)H_{\nu}^{R}(X)}}$, so $P({X_{\nu}}=x)={P^{\nu}}(X=x)/{e^{(1-\nu)H_{\nu}^{R}(X)}}$.

Another generalization of Shannon entropy in physic is the Tsallis entropy. For r.v. $X$, its Tsallis entropy of order $\alpha$ is defined by

[TABLE]

This entropy was introduced by Tsallis (1988) as a basis for generalizing the Boltzmann-Gibbs statistics. Also, the normalization constant $C(r,\nu,\tilde{p})$ in (4) has Tsallis entropy representation $C(r,\nu,\tilde{p})=1+(1-\nu)H_{\nu}^{T}(X)$, then $P({X_{\nu}}=x)={P^{\nu}}(X=x)/[1+(1-\nu)H_{\nu}^{T}(X)]$.

2.3 Log-concave, Log-convex, Infinite divisibility

This subsection deals with log-concavity and log-convexity of the CMNB distribution. A discrete distribution with ${p_{k}}=\mathrm{P}(X=k)$ is said to have log-concave (log-convex) p.m.f. if

[TABLE]

Lemma 2.1**.**

The CMNB distribution is log-concave if $r\geq 1$ and log-convex if $r\leq 1$.

Proof.

In fact, using the ratio of consecutive probabilities (6), we have

[TABLE]

Then $M\leq 1$ iff $r\geq 1$(log-concave) and $M\geq 1$ iff $r\leq 1$(log-convex). ∎

Remark 3: Ibragimov (1956) called a distribution strongly unimodal if it is unimodal and its convolution with any unimodal distribution is unimodal. He showed that the strongly unimodal distributions is equivalent to the log-concave distributions. So CMNB distribution has strong unimodality(see Figure 1 for example) when $r\geq 1$. Steutel (1970) showed that all log-convex discrete distributions are infinitely divisible, the background and detailed proof can be found in Steutel and van Harn (2003). Then we obtain infinite divisibility of CMNB distribution when $r\leq 1$.

Corollary 2.1**.**

The CMNB distribution (3) is discrete infinitely divisible (discrete compound Poisson distribution) if $r\leq 1$.

Feller’s characterization of the discrete infinite divisibility showed that a non-negative integer valued r.v. $X$ is infinitely divisibleiff its distribution is a discrete compound Poisson distribution with probability generating function (in short, p.g.f.):

[TABLE]

where $\sum_{i=1}^{\infty}{\alpha_{i}}=1,{\alpha_{i}}\geq 0,\lambda>0.$

For a theoretical treatment of discrete infinite divisibility (or discrete compound Poisson distribution), we refer readers to section 2 of Steutel and van Harn (2003), section 9.3 of Johnson et al. (2005), Zhang and Li (2016).

Considering some ${\alpha_{i}}$ being negative in (10), it turns into a generalization of the discrete compound Poisson distribution:

Definition 2.2** (Discrete pseudo compound Poisson distribution).**

If a discrete r.v. $X$ with $\mathrm{P}(X=k)=p_{k}$, $k\in\mathbb{N}$, has a p.g.f. of the form

[TABLE]

where $\sum\limits_{i=1}^{\infty}\alpha_{i}=1$, $\sum\limits_{i=1}^{\infty}\left|\alpha_{i}\right|<\infty$, $\alpha_{i}\in\mathbb{R}$, and $\lambda>0$, then $X$ is said to follow a discrete pseudo compound Poisson distribution, abbreviated as DPCP.

Next, we will give two lemmas on the non-vanishing p.g.f. characterization of DPCP, see Zhang et al. (2014) and Zhang et al. (2017).

Lemma 2.2**.**

Let $p_{k}=\mathrm{P}(X=k)$, for any discrete r.v. $X$, its p.g.f. $G(z)=\sum\limits_{k=0}^{\infty}{{p_{k}}}{z^{k}}$ has no zeros in $-1\leq z\leq 1$iff $X$ is DPCP distributed.

The proof of Lemma 2.2 is based on Wiener-Lévy theorem, which is a sophisticated theorem in Fourier analysis, see Zygmund (2002).

Lemma 2.3**.**

(Lévy-Wiener theorem) Let $F(\theta)=\sum\limits_{k=-\infty}^{\infty}{{c_{k}}{e^{ik\theta}}},\theta\in[0,2\pi]$ be a absolutely convergent Fourier series with $\left\|F\right\|=\sum\limits_{k=-\infty}^{\infty}{\left|{{c_{k}}}\right|}<\infty$. The value of $F(\theta)$ lies on a curve $C$, and $H(t)$ is an analytic (not necessarily single-valued) function of a complex variable which is regular at every point of $C$. Then $H[F(\theta)]$ has an absolutely convergent Fourier series.

Lemma 2.4**.**

For any discrete r.v. $X$ with p.g.f. $G(z)=\sum\limits_{k=0}^{\infty}{{p_{k}}}{z^{k}},(\left|z\right|\leq 1)$. If ${p_{0}}>{p_{1}}>{p_{2}}>\cdots$ , then $X$ is DPCP distributed.

Proof.

First, we show that $G(z)$ has no zeros in $\left|z\right|<1$, since

[TABLE]

And notice that $G(1)=1,G(-1)={p_{0}}-{p_{1}}+{p_{2}}-{p_{3}}+\cdots>0$, so $z=\pm 1$ are not zeros point. ∎

The condition in the next corollary is weaker than that of Corollary 2.1, and the result (DPCP) is also weaker than Corollary 2.1 (DCP).

Corollary 2.2**.**

The CMNB distribution (3) is discrete pseudo compound Poisson distribution if $(p{r^{\nu}}<1,r>1)$ or $(r\leq 1)$.

Proof.

On the one hand, $r\leq 1$ deduces that CMNB belongs to discrete compound Poisson by Corollary 2.1, hence CMNB is discrete pseudo compound Poisson. On the other hand, by using Lemma 2.4, we need to guarantee that $\frac{{\mathrm{P}(X=k)}}{{\mathrm{P}(X=k-1)}}=p{(\frac{{k-1+r}}{k})^{\nu}}<1$ for $k=1,2,\cdots$. The $p(\frac{{k-1+r}}{k})^{\nu}$ is a decreasing function with respect to $k$ when $r>1$, $\frac{{\mathrm{P}(X=k)}}{{\mathrm{P}(X=k-1)}}$ reaches its maximum $p{r^{\nu}}$ as $k=1$. So $p{r^{\nu}}<1$ is the other case. ∎

Applying the recurrence relation (Lévy-Adelson-Panjer recursion) of p.m.f. of DPCP distribution, see Remark 1 in Zhang et al. (2014)

[TABLE]

and ${P_{k}}=\mathrm{P}(X=k)={(\frac{{\Gamma(r+k)}}{{k!\Gamma(r)}})^{\nu}}{p^{k}}{P_{0}}$, then the DPCP parametrization $(\lambda,{\alpha_{1}},{\alpha_{2}},\cdots)$ of CMNB is determined by the following system of equations:

[TABLE]

where $\lambda=\log{P_{0}}$.

2.4 Overdispersion and underdispersion

In statistics, for a given random sample $X$, overdispersion means that $\mathrm{E}[X]<\mathrm{Var}(X)$. Conversely, underdispersion means that $\mathrm{E}[X]>\mathrm{Var}(X)$. Moreover, equal-dispersion means that $\mathrm{E}[X]=\mathrm{Var}(X)$. Gómez-Déniz (2011) summerized the phenomena of insurance count claims data, which were characterized by two features: (i) Overdispersion, i.e., the variance is greater than the mean; (ii) Zero-inflated, i.e. the presence of a high percentage of zero values in the empirical distribution.

The CMNB distribution belongs to the family of weighted Poisson distribution (see Kokonendji et al. (2008)) with p.m.f.

[TABLE]

where $w(k)$ is a non-negative weighted function.

Then weighted Poisson representation of CMNB distribution is

[TABLE]

Therefore, CMNB distribution in (3) can be seen as a weighted Poisson distribution with weighted function

[TABLE]

Theorem 3 and its corollary in Kokonendji et al. (2008) provide an “iff” condition to prove overdispersion and underdispersion of the weighted Poisson distribution.

Lemma 2.5**.**

Let $X$ be a weighted Poisson random variable with mean $\theta>0$, and let $w(k),k\in{\rm{N}}$ be a weighted function not depending on $\theta$. Then, weighted function $k\mapsto w(k)$ is logconvex (logconcave) iff the weighted version $X_{w}$ of $X$ is overdispersed (underdispersed).

Kokonendji et al. (2008) applied it to show that COM-Poisson distribution is overdispersion if $\nu<1$ and underdispersion if $\nu>1$. We employ their methods to get a criterion for overdispersion or underdispersion of CMNB distribution.

Theorem 2.1**.**

Set ${\Delta_{k}}=\sum\limits_{i=0}^{\infty}{(\frac{{1-\nu}}{{{{(i+k+1)}^{2}}}}+\frac{\nu}{{{{(i+k+r)}^{2}}}})}$. The COM-nagative binomial distribution (3) is overdispersion if ${\Delta_{k}}>0$ $(\forall k\in{\rm{N)}}$ and underdispersion if ${\Delta_{k}}<0$ $(\forall k\in{\rm{N)}}$.

Proof.

Function $f(x)$ is logconvex (logconcave) if $\frac{{{d^{2}}\log f(x)}}{{d{x^{2}}}}>0(<0)$. Followed by the formula of logarithmic second derivative of Gamma function (see p54 of Temme (2011)), $\frac{{{d^{2}}\log\Gamma(x)}}{{d{x^{2}}}}=\sum\limits_{i=0}^{\infty}{\frac{1}{{{{(x+i)}^{2}}}}}$, we have

[TABLE]

Applying Lemma 2.5, the proof is complete. ∎

Then, the results of overdispersion can be easily obtained by Theorem 2.1.

Corollary 2.3**.**

In these two cases: 1. $\nu>0,r<1$; 2. $\nu<1,r\in{{\rm{R}}^{+}}$. CMNB distribution is overdispersion.

Remark 4: The result of case 2 ($\nu<1,r\in{{\rm{R}}^{+}}$) can be also obtained from Corollary 2.1 and overdispersion of discrete compound Poisson distribution (equivalently, the discrete infinitely divisible).

2.5 Stochastic ordering

Stochastic ordering is the concept of one r.v. $X$ neither stochastically greater than, less than nor equal to another r.v. $Y$. There are plenty types of stochastic orders, which have various applications in risk theory. Firstly, we present 4 different definitions for discrete r.v.: usual stochastic order, likelihood ratio order, hazard rate order and mean residual life order.

1. $X$ is stochastically less than $Y$ in usual stochastic order (denoted by $X{\leq_{st}}Y$) if ${S_{X}}(n)\geq{S_{Y}}(n)$ for all $n$, where ${S_{X}}(n)=\mathrm{P}(X\geq n)=\sum\limits_{k=n}^{\infty}{{p_{k}}}$ is the survival function $X$ of with p.m.f. $p_{k}$ .

2. $X$ is stochastically less than $Y$ in likelihood ratio order (denoted by $X{\leq_{lr}}Y$) if $\frac{{g(n)}}{{f(n)}}$ increases in $n$ over the union of the supports of $X$ and $Y$, where $f(n)$ and $g(n)$ denotes the p.m.f. of $X$ and $Y$, respectively.

3. $X$ is stochastically less than $Y$ in hazard rate order (denoted by $X{\leq_{hr}}Y$) if ${r_{X}}(n)\geq{r_{Y}}(n)$ for all $n$, where the hazard function of a discrete r.v. $X$ with p.m.f. $p_{k}$ is defined as $r_{X}(n)={p_{n}}/\sum\limits_{k=n}^{\infty}{{p_{k}}}$.

4. $X$ is stochastically less than $Y$ in mean residual life order (denoted by $X{\leq_{MLR}}Y$) if ${\mu_{X}}(n)\geq{\mu_{Y}}(n)$ for all $n$, where the mean residual life function of a discrete r.v. $X$ with p.m.f. $p_{k}$ is defined as ${\mu_{X}}(n)={\rm{E}}(X-n\left|{X\geq n}\right.)=\sum\limits_{k=n}^{\infty}{k{p_{k}}}/\sum\limits_{k=n}^{\infty}{{p_{k}}}-n$.

The relationship among the above four stochastic ordering are $X{\leq_{lr}}Y\Rightarrow X{\leq_{hr}}Y\Rightarrow X{\leq_{MLR}}Y$

(see Theorem 1.C.1 of Shaked and Shanthikumar (2007)) and $X{\leq_{hr}}Y\Rightarrow X{\leq_{st}}Y$ (see Theorem 1.B.1 of Shaked and Shanthikumar (2007)).

Gupta et al. (2014) gave the stochastic ordering between COM-Poisson r.v. $X$ and Poisson distributed r.v. $Y$ with same parameter $\lambda$ in (1), that is $X{\leq_{lr}}Y$, therefore $X{\leq_{st}}Y$, $X{\leq_{hr}}Y$ and $X{\leq_{MLR}}Y$. In the following result we will show that CMNB distribution also has some stochastic ordering properties.

Theorem 2.2**.**

Let $X$ and $Y$ be two r.vs following CMNB distribution with parameters $(r,{\nu}_{1},p)$ and $(r,{\nu}_{2},p)$, respectively. If ${\nu_{1}}\leq{\nu_{2}},r\geq 1$, then $X{\leq_{lr}}Y$, hence $X{\leq_{st}}Y$, $X{\leq_{hr}}Y$ and $X{\leq_{MLR}}Y$.

Proof.

Note that $r\geq 1$, we have $\frac{{\Gamma(r+n+1)}}{{(n+1)!}}/\frac{{\Gamma(r+n)}}{{n!}}=\frac{{r+n}}{{n+1}}\leq 1$. Then

[TABLE]

which is increasing in $n$ as ${\nu_{1}}\geq{\nu_{2}}$. ∎

Especially, assume that $r$ is a positive integer, the CMNB should be called the COM-Pascal distribution. Let $X$ be CMNB distributed and $Y$ be negative binomial distributed with the same parameters $r,p$, it yields to $X{\leq_{lr}}Y$ when ${\nu_{2}}>{\nu_{1}}=1$.

The next theorem is proved in the view of weighted Poisson distribution (12) from weighted function of CMNB distribution. Example 1.C.59 of Shaked and Shanthikumar (2007) states the obvious lemma below:

Lemma 2.6**.**

Define $X_{w}$ as the r.v. with weighted density function ${f_{w}}(x)=\frac{{w(x)}}{{{\rm{E[}}w(X)]}}f(x),(x\geq 0)$, Similarly, for another nonnegative r.v. $Y$ with density function $g$, define $Y_{w}$ as the r.v. with the weighted density function ${g_{w}}(y)=\frac{{w(y)}}{{{\rm{E[}}w(Y)]}}g(y),(y\geq 0)$. If $w(x)$ is an increasing function, then $X{\leq_{hr}}Y\Rightarrow{X_{w}}{\leq_{hr}}{Y_{w}}$.

Theorem 2.3**.**

Let $X$ and $Y$ be two CMNB distributed with parameters $(r,{\nu},p_{1})$ and $(r,{\nu},p_{2})$, respectively. If $({p_{1}}\leq{p_{2}},\nu\leq 1)$ or $({p_{1}}\leq{p_{2}},r\geq 1)$, then $X{\leq_{lr}}Y$, and therefore $X{\leq_{st}}Y$, $X{\leq_{hr}}Y$ and $X{\leq_{MLR}}Y.$

Proof.

For Poisson distributed $X$,$Y$ with mean $p_{1}$,$p_{2}$, if ${p_{1}}\leq{p_{2}}$, then ${\rm{P}}(Y=n)/{\rm{P}}(X=n)={(\frac{{{p_{2}}}}{{{p_{1}}}})^{n}}{e^{-({p_{2}}-{p_{1}})}}$ is increasing for all $n$. So $X{\leq_{hr}}Y$. From section 2.4, we know that CMNB is weight Poisson with weight (13).

On the one hand, when $\nu\leq 1$, we notice that weighted density function $w(x)={[\Gamma(1+x)]^{1-\nu}}{[\Gamma(r+x)]^{\nu}}$ for CMNB distribution is increasing with respect to $x$. On the other hand, $w(x)=\Gamma(1+x){\left({\frac{{\Gamma(r+x)}}{{\Gamma(1+x)}}}\right)^{\nu}}$ is an increasing function with respect to $x$ as $r\geq 1$, that is, ${X_{w}}{\leq_{hr}}{Y_{w}}$. ∎

2.6 Approximate to COM-Poisson distribution

The next theorem enables CMNB distribution to be a suitable generalization since its limit distribution is the COM-Poisson under some conditions. We prove that CMNB distribution converges to the COM-Poisson distribution when $r$ goes to infinity.

Theorem 2.4**.**

Suppose that r.v. $X$ has CMNB distribution with parameters $(r,\nu,p)$, denote the p.m.f. as $\mathrm{P}(\left.{X=k}\right|r,\nu,p)$, and let $\lambda={r^{\nu}}\frac{p}{{1-p}}$. Then

[TABLE]

Proof.

Notice that $p=\frac{\lambda}{{{r^{\nu}}+\lambda}}$, substitute to p.m.f (3), then we obtain

[TABLE]

Hence,

[TABLE]

holds as $\mathop{\lim}\limits_{r\to+\infty}{\left({\frac{{\Gamma(r+k)}}{{\Gamma(r)\;{r^{k}}}}}\right)^{\nu}}\cdot\frac{1}{{{{(1+\lambda/{r^{v}})}^{k}}}}=1,(k=0,1,2,\cdots)$ and

[TABLE]

∎

3 Characterizations

3.1 Sum of equicorrelated geometrically distributed r.v.

It is well known that the binomial r.v. can be seen as the sum of $m$ independent Bernoulli r.v. $Z_{i}$.

[TABLE]

where

[TABLE]

[TABLE]

Shmueli et al. (2005) and Borges et al. (2014) mentioned that the COM-binomial distribution (9) can be presented as a sum of equicorrelated Bernoulli r.vs $\{{Z_{i}}\}_{i=1}^{m}$ with joint distribution

[TABLE]

where $k=\sum_{i=1}^{m}z_{i}$.

As we know, negative binomial distribution can be treated as the sum of $m$ independent geometric r.vs $Z_{i}$ $(i=1,\cdots,m)$:

[TABLE]

[TABLE]

where $\mathrm{P}({Z_{i}}={z_{i}})=p^{z_{i}}(1-p),(z_{i}=1,2,,\cdots).$

It is similar to see that the CMNB distribution can be interpreted as a sum of equicorrelated geometric r.vs $Z_{i}$ $(i=1,\cdots,m)$ with joint distribution

[TABLE]

where $x=\sum_{i=1}^{m}z_{i}$.

The reason is that we assume $\{{Z_{i}}\}_{i=1}^{m}$ is equicorrelated, and $x=\sum_{i=1}^{m}z_{i}$ has ${{m+x-1}\choose{x}}$ feasible positive integer solutions, and each solution has probability $\mathrm{P}({Z_{1}}={z_{1}},\cdots,{Z_{m}}={z_{m}})$. Then, the ${{m+x-1}\choose{x}}$ possible values of random vector $({Z_{1}},\cdots,{Z_{m}})$ such that $S=\sum_{i=1}^{m}Z_{i}$ is CMNB distributed, namely

[TABLE]

Thus we have (15).

3.2 Conditional distribution

In this subsection, two conditional distribution characterizations are obtained for CMNB distribution. For two independent r.vs $X,Y$, what is the form of the conditional distribution of $X$ given $S=X+Y$? Consider the sum of CMNB r.vs with parameters $(r_{x},\nu,p)$ and $(r_{y},\nu,p)$, then

[TABLE]

The conditional distribution $\mathrm{P}(X=k\left|{S=s}\right.)$ is

[TABLE]

Using (16), we naturally define the p.m.f. of COM-negative hypergeometric distribution with parameter $(z,\nu,m,n)$ as follow:

[TABLE]

where $k=0,1,2,\cdots,z$ and $N(m,n,z,\nu)$ is the normalization constant.

When $\nu=1$, COM-negative hypergeometric distribution turns out to be negative hypergeometric distribution, see Wimmer and Altmann (1999), Johnson et al. (2005).

By Patil and Seshadri (1964)’s general characterization theorem for negative binomial, Poisson and geometric distribution, we know that, given $X+Y=x+y$, if the conditional distribution $X\left|{X+Y}\right.$ is negative hypergeometric with parameters $m$ and $n$ for all values of the sum $x+y$,

[TABLE]

then $X$ and $Y$ both are negative binomial distribution, with parameters $(m,p)$ and $(n,p)$,

[TABLE]

respectively (see also Kagan et al. (1973)).

Lemma 3.1**.**

(Patil and Seshadri (1964)) Let $X$ and $Y$ be independent and both discrete (or both continuous) r.vs and suppose $\mathrm{P}(X\left|{X+Y}\right.)$ is the function $c(x,x+y)$. If $\frac{{c(x+y,x+y)c(0,y)}}{{c(x,x+y)c(y,y)}}$ is of the form $h(x+y)/h(x)h(y)$ where $h(\cdot)$ is an arbitrary non-negative function, then

[TABLE]

where $\mathrm{P}(X=x)=f(x)>0,\mathrm{P}(Y=y)=g(y)>0$ and $f(0),g(0)$ are the corresponding normalizer for $f(x)$ and $g(y)$ which make them p.m.f..

Now we apply Lemma 3.1 for characterizing CMNB distribution.

Theorem 3.1**.**

Let X,Y be the independent discrete r.v. with $\mathrm{P}(X=x)=f(x)>0$ and $\mathrm{P}(Y=y)=g(y)>0$. If the $\mathrm{P}(X=x\left|{X+Y=x+y}\right.)$ is the COM-negative hypergeometric distribution (17) with parameters $(z=x+y,\nu,m,n)$ for all $x+y$, then both $X$ and $Y$ have the CMNB distributions with the parameters $(m,\nu,p)$ and $(n,\nu,p)$ respectively.

Proof.

Note that $c(x,x+y)={\left({{{x+y}\choose{x}}\frac{{{\rm{B}}({m}+x,{n}+x+y-x)}}{{{\rm{B}}({m},{n})}}}\right)^{\nu}}/{{N(m,n,z,\nu)}}$. Then

[TABLE]

so

[TABLE]

We have $h(x)={\left({\frac{{\Gamma(m+x)}}{{(x)!\Gamma(m)}}}\right)^{v}},h(y)={\left({\frac{{\Gamma(m+y)}}{{(y)!\Gamma(m)}}}\right)^{v}}$. Apply (18), $\frac{{h(y)c(0,y)}}{{c(y,y)}}={\left({\frac{{\Gamma(n+y)}}{{\Gamma(n)y!}}}\right)^{\nu}}$, let $p={e^{a}}$ and compare with expression (3), hence the proof is complete. ∎

Rao and Rubin (1964) study the following characterization of the Poisson distribution: If $X$ is a discrete r.v. taking only nonnegative integer values and the conditional distribution of $Y$ given $X=x$ is binomial distribution with parameters $\texttt{B}i(x,p)$ ($p$ does not depend on $x$), then $X$ follows the Poisson distribution iff

[TABLE]

Based on the COM-negative hypergeometric distribution above and an extension of Rao and Rubin (1964)’s characterization which established by Shanbhag (1977), the Rao-Rubin characterization for CMNB is obtained.

Lemma 3.2**.**

*(Shanbhag (1977)) Let $X,Y$ be the non-negative r.vs such that $\mathrm{P}(X=z)={p_{z}}$ with

${p_{0}}<1,\,{p_{z}}>0$, and*

[TABLE]

where ${a_{z}}>0$ for all $z\geq 0$, ${b_{0}},{b_{1}}>0$ and ${b_{z}}\geq 0$ for $z\geq 2$, then

[TABLE]

Next, we give the following CMNB extension of Rao-Rubin characterization. The Shanbhag’s result is vital to prove this extension immediately.

Theorem 3.2**.**

Let $X,Y$ be the non-negative integer-valued r.vs such that $\mathrm{P}(X=z)={P_{z}}$ with ${P_{0}}<1,{P_{z}}>0$, and

[TABLE]

then $\mathrm{P}(Y=k)=\mathrm{P}(\left.{Y=k}\right|X=Y),(r=0,1,\cdots)$ iff $X$ follows the CMNB distribution with the parameters $(m+n,\nu,\theta)$ for some $\theta>0$.

Proof.

Since the alternative form of negative hypergeometric distribution can be written as

[TABLE]

From the normalizer in (17), let $L(m,n,z,\nu)=\frac{{N(m,n,z,\nu)}}{{{{[B(m,n)]}^{-\nu}}}}$, we get

[TABLE]

where ${a_{z}},{b_{z}},{c_{z}}$ satisfy the conditions in Lemma 3.2.

Then, ${c_{z}}=\frac{{{{[\Gamma(m+n+z)]}^{\nu}}}}{{{{(z!)}^{\nu}}}}$, compareing with form (3), we conclude that $X$ is the CMNB distributed with the parameters $(m+n,\nu,\theta)$. ∎

3.3 COM-negative hypergeometric approximate to CMNB

We have shown that the conditional distribution of CMNB r.v. $X$ given the sum of two CMNB random variables $X$ and $Y$ follows the COM-negative hypergeometric distribution (17). Next, we show that the COM-negative hypergeometric distribution (17) converges to the CMNB distribution (3).

Theorem 3.3**.**

Let $X$ be the COM-negative hypergeometric r.v. with parameters $(z,\nu,m,n)$ and p.m.f. $\mathrm{P}(X=k|z,\nu,m,n)$ given by (17), Assume $m<\infty$, $\frac{z}{{n+z}}={p^{\frac{1}{\nu}}}$, $p\in(0,1)$, then

[TABLE]

Proof.

First, we multiply the p.m.f of COM-negative hypergeometric distribution by the normalization constant (see (17)):

[TABLE]

Using the Stirling’s formula for the Gamma function, $\Gamma(z)=\sqrt{\frac{{2\pi}}{z}}{\mkern 1.0mu}{\left({\frac{z}{e}}\right)^{z}}\left({1+O\left({\frac{1}{z}}\right)}\right)$, it yields

[TABLE]

From the conditions, we have $\frac{m+n}{m+n+z-k}\rightarrow 1-p^{\frac{1}{\nu}}$ as $n\rightarrow\infty$. Then

[TABLE]

and

[TABLE]

Finally, we obtain

[TABLE]

Note that the normalizing constant $\mathop{\lim}\limits_{n\to\infty}{[N(m,n,z,\nu)]^{-1}}{\left({\frac{{{{(1-{p^{1/\nu}})}^{\nu}}}}{{1-p}}}\right)^{m}}$ exsists and does not depend on $k$. Therefore, $\mathop{\lim}\limits_{n\to\infty}\mathrm{P}(X=k|z,\nu,m,n)$ is the p.m.f. of ${\rm{CMNB}}(m,\nu,p)$ via comparing with (3). ∎

Borges et al. (2014) showed that the COM-Poisson distribution is the limiting distribution of the COM-binomial distribution. Let $X_{m}\sim{\rm{CMB}}(m,p,\nu)$(see (9)) and $\lambda=m^{\nu}p$ for $\nu\geq 0$, then

[TABLE]

namely $\lim_{m\rightarrow\infty}X_{m}\sim{\rm{CMP}}(\lambda,\nu)$.

To sum up, judging from the above mentioned theorems of limiting distribution, we may naturally draw the relationships among some COM-type distributions in Figure 2.

3.4 Functional operator characterization by Stein’s identity

Using recursive formula (6) for CMNB distribution, an extension of functional operator characterization which arises from the negative binomial approximation literatures by Stein-Chen method is presented. The following Stein’s lemma also can be derived from the work of Brown and Xia (2001) who give a large class of approximation distribution which is the equilibrium distribution of a birth-death process with arbitrary arrival rate and service rate.

The functional operator characterization is well-known in the Stein-Chen method literatures. Lemma 3.3 extends the Lemma 1 in Brown and Phillips (1999) for negative binomial approximation.

Lemma 3.3**.**

The r.v. $W$ has distribution ${\rm{CMNB}}(r,\nu,p)$iff the equation

[TABLE]

holds for any bounded function $g:\mathbb{Z}^{+}\to\mathbb{R}$.

Proof.

Sufficiency: When $W$ is CMNB distributed in the form of (3), $f$ is a bounded function. Then we have

[TABLE]

Hence $\mathbb{E}g(W)=0$.

Necessity: $\mathbb{E}g(W)=0$ for all bounded functions. Letting $f(W)={1_{(W=k)}}(W)$, a simple computation shows the recursive formula (6). This verifies that $W\sim{\rm{CMNB}}(r,\nu,p)$. ∎

4 Generating CMNB-distributed random variables

The inverse transform sampling is a basic method in pseudo-random number sampling, namely for generating sample numbers at random from any probability distribution given its cumulative distribution function. The following lemmas are the basis of inverse method.

Lemma 4.1**.**

Suppose $F(x)$ is the cumulative distribution function of random variable $\xi$, then $\theta=F(\xi)$ is uniformly distributed in interval $[0,1]$.

Lemma 4.2**.**

Suppose $\theta\sim U[0,1]$, for any cumulative distribution function $F(x)$, let $\xi=F^{-1}(\theta)$, then $\xi\sim F(x)$, that is, $F(x)$ is the distribution function of random variable $\xi$.

The proof of Lemma 4.1 and 4.2 is trivial and can be found in text books of computational statistics.

Consider generating random sample from CMNB$(r,\nu,p)$, since the p.m.f. is

[TABLE]

It follows easily that

[TABLE]

Multiplying these equations yields

[TABLE]

Thus, for any non-negative $k$, the distribution function $F(k)=\sum_{i=0}^{k}\text{P}(X=i)$ can be easily calculated. By Lemma 4.2, the inverse method alternates between the following steps:

Step 1:

Let $j=0$, $\text{P}(X=j)=(1-p)^{r}/C(r,\nu,p)$;

Step 2:

Generate a random sample $\theta\sim U[0,1]$, and fix its value;

Step 3:

Let $F(j+1)=F(j)+\text{P}(X=j+1)$;

Step 4:

If $F(j)<\theta\leq F(j+1)$, then $x=j+1$; Otherwise, let $j=j+1$, and back to Step 3;

We conclude that, $x$ is the random number drawn from the CMNB distribution with parameters $(r,\nu,p)$.

5 Estimations

5.1 Estimation by three recursive formulas

The method of recursive formula estimation (or refer as ratio regression, see Bohning (2016)) which give crude estimations of the estimated parameters, is originated with Shmueli et al. (2005), and extended by Imoto (2014). However, this crude estimations can be put into the maximum likelihood estimation (MLE) with Newton-Raphson algorithm as initial values.

First, we note that the p.m.f. of CMNB has the following recursive formula:

[TABLE]

Second, applied with the expression of log-concave (log-convex), we have

[TABLE]

Replace “$k$” by “$k+1$” in (23). Then we obtain

[TABLE]

Thus, $(\nu,r,p)$ could be sequentially solved by the system of equations (22), (23) and (24), if we have the sample p.m.f.

5.2 Maximum likelihood estimation

Let r.v. $X$ be distributed as the CMNB distribution with parameters $\theta=(r,\,\nu,\,p)$. We consider the MLE in the case parameters $r$, $\nu$, $p$ are unknown. The likelihood function is

[TABLE]

where $n$ is the sample size, $x_{1},x_{2},\ldots,x_{n}$ are the observed values, and the log-likelihood function is given by

[TABLE]

To find the maximum point, for $\log L(\theta)$, we take the partial derivatives with respect to $r$, $\nu$ and $p$ and set them equal to zero, hence the likelihood equation is given by

[TABLE]

where $\psi(\alpha)=\Gamma^{{}^{\prime}}(\alpha)/\Gamma(\alpha)$ is called the digamma function. The analytical solutions of the above likelihood equations are not tractable, therefore, numerical optimization method is used to obtain the maximum likelihood estimates. Let F(\theta)=\big{(}F_{1}(\theta),\,F_{2}(\theta),\,F_{3}(\theta)\big{)}^{T}, then the Fisher information matrix $I(r,\nu,p)$ would be the Jacobin matrix of $F(\theta)$. Given trial values $\theta_{k}=(r_{k},\,\nu_{k},\,p_{k})$, applied with the scoring method for solving (27), we can update to $\theta_{k+1}=(r_{k+1},\,\nu_{k+1},\,p_{k+1})$ as

[TABLE]

For the initial point $\theta_{0}=(r_{0},\,\nu_{0},\,p_{0})$, we can choose the crude estimated parameters introduced in Subsection 5.1.

5.3 Kolmogorov-Smirnov and Chi-squared goodness-of-fit test

In goodness-of-fit test of discrete distributions, Pearson’s chi-squared test is a popular choice to check distribution model, which is usually better than the others. The larger p-value, the better goodness-of-fit would arise from the assumed models. However, when the data come from an assumed model, with $r$ unknown parameters $\theta=({\theta_{1}},\cdots,{\theta_{r}})$ and the null hypothesis ${H_{0}}$, what we are interested in is

[TABLE]

The $\chi^{2}$ statistic has the following limiting distribution:

[TABLE]

where $n_{i}$ is the number of class $i$, and $m$ classes ${A_{i}}~{}(i=1,2,\cdots,m)$ which are the sub-division of $n$ samples, it satisfies that:

[TABLE]

We also have the another handy expression for calculating $\chi^{2}$ statistic: $\eta=\sum\limits_{i=1}^{m}{\frac{{n_{i}^{2}}}{{n{{\hat{p}}_{0i}}}}-n}$.

One limitation of $\chi^{2}$ statistic is that, if degree of freedom are too small, the approximation to the $\chi^{2}$ distribution would fail. For example, $k=m-1-r=0,1,2$. Since $\chi^{2}$ test of small degree of freedom (denotes it by $k$) did not have enough power (namely ${\mathrm{P}_{{\theta}}}(\eta>{\eta_{0}})$ is not large enough). That is to say, for a $\chi^{2}$ r.v. $\eta$, the p-value is defined by

[TABLE]

which tends to be small when $k$ varies from 3 to 1. An attempt at this has been made in the figure below.

Another drawback is that, with small number of class $m$, Haberman (1988) discussed that: “chi-squared test statistics may be asymptotically inconsistent even in cases in which a satisfactory chi-squared approximation exists for the distribution under the null hypothesis.”

Consequently, a better way of comparing distributions is to use a non-parametric test not depending on the parametric assumption of the models. One of notable non-parametric test is Kolmogorov-Smirnov test by using following statistic:

[TABLE]

which is a nonparametric test of continuous distributions. And it could be modified for discrete distribution:

[TABLE]

where $x_{i}$’s are the discontinuity points ($I$ is the countable index set).

Fortunately, Arnold and Emerson (2011) developed the R package “dgof”, which can calculate Kolmogorov-Smirnov test of the discrete distribution. The non-parametric test avoids the assumption of parametric model, hence it is an effective way to evaluate different distributions’ performance.

6 Applications with simulated and real data

In this section, we will describe four examples of fitting data by CMNB distribution, and we will compare them with those by the negative binomial and COM-Poisson distribution. The notations, CMNB, NB and CMP distribution, in the following tables, are representing the CMNB, negative binomial and COM-Poisson, respectively. By using scoring method in section 5.2, we could estimated the CMNB distribution with three parameters. The original and expected frequencies, parameter estimators (obtained by maximum likelihood method), the $\chi^{2}$ and K-S statistics, and corresponding $p$-values. are all being considered in the tables below.

6.1 Simulated Data

In this subsection, the inverse method mentioned in section 4 were used to generate random variables from CMNB, NB and COM-Poisson distribution with certain parameters. And we compare the performance of the fit by CMNB, NB and COM-Poisson distribution.

Example 1. As our first simulation example, we consider the original data that are following the NB distribution. We use inverse method to draw 10000 random samples from NB distribution with parameter $(r=1,p=0.5)$, and fit the data with above-mentioned distributions, the results are shown in the following table 1.

In this case, although the data are generated from NB distribution, the CMNB distribution can also be recognized as the true distribution, as the estimator of $\nu$ is 1. That is to say, the data can be seen negative binomial distributed, and note that he result of $\chi^{2}$ test and K-S test is nearly the same, all of these indicate that the CMNB distribution is an extension of NB distribution.

Example 2. In this example, we will evaluate the performance of NB distribution, when the original data are generated from CMNB distribution. Assume $X_{1},X_{2},\ldots,X_{n}$ are independent and identically CMNB distributed random variables with parameters $(r,\nu,p)$. Let $n=10000$ and $(r,\nu,p)=(0.005,0.1,0.5)$, we can generate $n$ sample points $x_{1},x_{2},\ldots,x_{n}$ through inverse method introduced in section 4. The fitting results are summarized below:

From table 2, we note first that in $\chi^{2}$ test, the CMNB distribution’s $\chi^{2}$ statistic 7.26 is dramatically smaller than NB distribution’s 16.38, while the corresponding p-values are 0.61 and 0.09 respectively, which implies that NB distribution may not be a reasonable fit. As for the result of the K-S test, it suggests that CMNB distribution performs best, because the p-value of CMNB distribution is approximately 1, while NB’s and COM-Poisson’s are 0.79 and 0.

6.2 Real data analysis

In this subsection, the distributions mentioned above are considered here to analyze real actuarial claim data that have ultrahigh zero-inflated and overdispersion properties, then compare its Kolmogorov-Smirnov test and Chi-squared test.

Example 3. Let us consider the claim counts of the third party liability vehicle insurance (see Willmot (1987) for data set in an Zaire insurance company) which correspond to claims from 4000 vehicle policies. Gómez-Déniz (2014) analyzed the data using negative binomial distribution and found that it is a reasonable fit. We analyze the data by CMNB, NB and CMP distribution, and the results are summarized below:

It’s a well known fact that NB distribution is a popular choice to fit the claim data in actuarial science. However findings in Table 3 suggest that, CMNB might be a better choice for fitting this data. It appears that, for chi-square statistic, fitting the chi-square statistic with the NB distribution, it turns out to be $1.56$, which is dramatically larger than CMNB’s $1.01$, and the p-values correspond to the $\chi^{2}$ statistics is nearly the same. As for Kolmogorov-Smirnov test, although the p-value is both 1, it still can be seen that the K-S statistic of CMNB is slightly smaller than NB’s. All of these show that the CMNB is superior to the fit the data.

Example 4. For this example, we use the car insurance claim data of a Chinese insurance company (see Wang and Lei, (2000)), which were modeled with negative binomial distribution. Total insurance policies are $n=35072$. We analyze the data with the above-mentioned distributions, and the results are shown in table 4.

As can be observed from table 4, CMNB is the clear winner. CMNB distribution outperforms other distributions by either $\chi^{2}$ test or K-S test, as $\chi^{2}$ statistic and K-S statistic of CMNB is slightly smaller than that of NB’s, and dramatically smaller than CMP’s. It shows that their performance is on the whole dominated by CMNB distribution.

7 Further researches

Many discrete distributions in aspects of related statistical models are widely proposed in plenty of research articles. The CMNB distribution can be used in driving count data in some generalized linear models (GLMs). For example, Poisson regression, COM-Poisson regression and negative binomial regression (NBR), whose driving distribution of count data are actually the special case or limiting case of CMNB distribution, see Figure 2 for a visual relations. Even in the popular logistic regression, the assuming distribution lying in this model is Bernoulli distribution which is a limiting case of COM-Poisson distribution. Besides GLMs, high-dimensional GLMs (i.e., the sample size $n$ may be greater than the number of covariates $p$) should also be considered for CMNB, see Zhang and Jia (2017) and references therein for high-dimensional NBR. The discrete frailty item as a random effect in survival analysis model (or Long-term survival models, see Rodrigues et al. (2012)), which could be considered to set some flexible distributions. Besides the researches of regression model, in our paper, both conditional distribution and Stein identity characterizations are obtained. Test of statistic based on the novel characterization is more reasonable, since that the distribution-free and omnibus goodness-of-fit test like Kolmogorov-Smirnov and Chi-squared does not depend on certain types(or families) of distributions, may have their own drawbacks. For example, some new tests of a class of count distributions which includes the Poisson is constructed from Stein identity characterization (see Meintanis and Nikitin (2008)), and some tests of logarithmic series distribution comes from conditional distribution characterization (see Ramalingam and Jagbir (1984)). The goodness-of-fit tests between two continuous probability distributions based on combining Stein s identity is studied by Liu, Lee and Jordan (2016), and we consider study Stein discrepancy goodness-of-fit tests for discrete distribution in the future.

8 Acknowledgments

The proposed COM-negative binomial distribution of this work was as early as conceptualized in Dec. 2014 when the authors saw the online version of Imoto (2014). The authors want to thank Prof. Köhler R. for mailing the valuable encyclopedia of discrete univariate distributions (Wimmer and Altmann (1999)) to us. This work was partly supported by the National Natural Science Foundation of China (No.11201165).

9 References