Hermite distribution

Let X₁ and X₂ be two independent Poisson variables with parameters a₁ and a₂. The probability distribution of the random variable Y = X₁ + 2X₂ is the Hermite distribution with parameters a₁ and a₂ and probability mass function is given by [8]

${\displaystyle p_{n}=P(Y=n)=e^{-(a_{1}+a_{2})}\sum _{j=0}^{\lfloor n/2\rfloor }{\frac {a_{1}^{n-2j}a_{2}^{j}}{(n-2j)!j!}}}$

where

• n = 0, 1, 2, ...
• a₁, a₂ ≥ 0.
• (n − 2j)! and j! are the factorials of (n − 2j) and j, respectively.
• ${\textstyle \lfloor n/2\rfloor }$ is the integer part of  n/2.

The probability generating function of the probability mass is,[8]

${\displaystyle G_{Y}(s)=\sum _{n=0}^{\infty }p_{n}s^{n}=\exp(a_{1}(s-1)+a_{2}(s^{2}-1))}$

When a random variable Y = X₁ + 2X₂ is distributed by an Hermite distribution, where X₁ and X₂ are two independent Poisson variables with parameters a₁ and a₂, we write

${\displaystyle Y\ \sim \operatorname {Herm} (a_{1},a_{2})\,}$

The moment generating function of a random variable X is defined as the expected value of e^t, as a function of the real parameter t. For an Hermite distribution with parameters X₁ and X₂, the moment generating function exists and is equal to

${\displaystyle M(t)=G(e^{t})=\exp(a_{1}(e^{t}-1)+a_{2}(e^{2t}-1))}$

The cumulant generating function is the logarithm of the moment generating function and is equal to [4]

${\displaystyle K(t)=\log(M(t))=a_{1}(e^{t}-1)+a_{2}(e^{2t}-1)}$

If we consider the coefficient of (it)^rr! in the expansion of K(t) we obtain the r-cumulant

${\displaystyle k_{n}=a_{1}+2^{n}a_{2}}$

Hence the mean and the succeeding three moments about it are

Order	Moment	Cumulant
1	${\displaystyle \mu _{1}=k_{1}=a_{1}+2a_{2}}$	${\displaystyle \mu }$
2	${\displaystyle \mu _{2}=k_{2}=a_{1}+4a_{2}}$	${\displaystyle \sigma ^{2}}$
3	${\displaystyle \mu _{3}=k_{3}=a_{1}+8a_{2}}$	${\displaystyle k_{3}}$
4	${\displaystyle \mu _{4}=k_{4}+3k_{2}^{2}=a_{1}+16a_{2}+3(a_{1}+4a_{2})^{2}}$	${\displaystyle k_{4}}$

The skewness is the third moment centered around the mean divided by the 3/2 power of the standard deviation, and for the hermite distribution is,[4]

${\displaystyle \gamma _{1}={\frac {\mu _{3}}{\mu _{2}^{3/2}}}={\frac {a_{1}+8a_{2}}{(a_{1}+4a_{2})^{3/2}}}}$

• Always ${\displaystyle \gamma _{1}>0}$, so the mass of the distribution is concentrated on the left.

The kurtosis is the fourth moment centered around the mean, divided by the square of the variance, and for the Hermite distribution is,[4]

${\displaystyle \beta _{2}={\frac {\mu _{4}}{\mu _{2}^{2}}}={\frac {a_{1}+16a_{2}+3(a_{1}+4a_{2})^{2}}{(a_{1}+4a_{2})^{2}}}={\frac {a_{1}+16a_{2}}{(a_{1}+4a_{2})^{2}}}+3}$

The excess kurtosis is just a correction to make the kurtosis of the normal distribution equal to zero, and it is the following,

${\displaystyle \gamma _{2}={\frac {\mu _{4}}{\mu _{2}^{2}}}-3={\frac {a_{1}+16a_{2}}{(a_{1}+4a_{2})^{2}}}}$

• Always ${\displaystyle \beta _{2}>3}$, or ${\displaystyle \gamma _{2}>0}$ the distribution has a high acute peak around the mean and fatter tails.

In a discrete distribution the characteristic function of any real-valued random variable is defined as the expected value of ${\displaystyle e^{itX}}$, where i is the imaginary unit and t ∈ R

${\displaystyle \phi (t)=E[e^{itX}]=\sum _{j=0}^{\infty }e^{ijt}P[X=j]}$

This function is related to the moment-generating function via ${\displaystyle \phi _{x}(t)=M_{X}(it)}$. Hence for this distribution the characteristic function is,[1]

${\displaystyle \phi _{x}(t)=\exp(a_{1}(e^{it}-1)+a_{2}(e^{2it}-1))}$

The cumulative distribution function is,[1]

${\displaystyle {\begin{aligned}F(x;a_{1},a_{2})&=P(X\leq x)\\&=\exp(-(a_{1}+a_{2}))\sum _{i=0}^{\lfloor x\rfloor }\sum _{j=0}^{[i/2]}{\frac {a_{1}^{i-2j}a_{2}^{j}}{(i-2j)!j!}}\end{aligned}}}$

• This distribution can have any number of modes. As an example, the fitted distribution for McKendrick’s [2] data has an estimated parameters of ${\displaystyle {\hat {a}}_{1}=0.0135}$, ${\displaystyle {\hat {a}}_{2}=0.0932}$. Therefore, the first five estimated probabilities are 0.899, 0.012, 0.084, 0.001, 0.004.

Example of a multi-modal data, Hermite distribution(0.1,1.5).

• This distribution is closed under addition or closed under convolutions.[9] Like the Poisson distribution, the Hermite distribution has this property. Given two Hermite-distributed random variables ${\displaystyle X_{1}\sim \operatorname {Herm} (a_{1},a_{2})}$ and ${\displaystyle X_{2}\sim \operatorname {Herm} (b_{1},b_{2})}$, then Y = X₁ + X₂ follows an Hermite distribution, ${\displaystyle Y\sim \operatorname {Herm} (a_{1}+b_{1},a_{2}+b_{2})}$.
• This distribution allows a moderate overdispersion, so it can be used when data has this property.[9] A random variable has overdispersion, or it is overdispersed with respect the Poisson distribution, when its variance is greater than its expected value. The Hermite distribution allows a moderate overdispersion because the coefficient of dispersion is always between 1 and 2,

${\displaystyle d={\frac {\operatorname {Var} (Y)}{\operatorname {E} (Y)}}={\frac {a_{1}+4a_{2}}{a_{1}+2a_{2}}}=1+{\frac {2a_{2}}{a_{1}+2a_{2}}}}$

The mean and the variance of the Hermite distribution are ${\displaystyle \mu =a_{1}+2a_{2}}$ and ${\displaystyle \sigma ^{2}=a_{1}+4a_{2}}$, respectively. So we have these two equation,

${\displaystyle {\begin{cases}{\bar {x}}=a_{1}+2a_{2}\\\sigma ^{2}=a_{1}+4a_{2}\end{cases}}}$

Solving these two equation we get the moment estimators ${\displaystyle {\hat {a_{1}}}}$ and ${\displaystyle {\hat {a_{2}}}}$ of a₁ and a₂.[6]

${\displaystyle {\hat {a_{1}}}=2{\bar {x}}-\sigma ^{2}}$

${\displaystyle {\hat {a_{2}}}={\frac {\sigma ^{2}-{\hat {x}}}{2}}}$

Since a₁ and a₂ both are positive, the estimator ${\displaystyle {\hat {a_{1}}}}$ and ${\displaystyle {\hat {a_{2}}}}$ are admissible (≥ 0) only if, ${\displaystyle {\bar {x}}<\sigma ^{2}<2{\bar {x}}}$.

Given a sample X₁, ..., X_m are independent random variables each having an Hermite distribution we wish to estimate the value of the parameters ${\displaystyle {\hat {a_{1}}}}$ and ${\displaystyle {\hat {a_{2}}}}$. We know that the mean and the variance of the distribution are ${\displaystyle \mu =a_{1}+2a_{2}}$ and ${\displaystyle \sigma ^{2}=a_{1}+4a_{2}}$, respectively. Using these two equation,

${\displaystyle {\begin{cases}a_{1}=\mu (2-d)\\[4pt]a_{2}={\dfrac {\mu (d-1)}{2}}\end{cases}}}$

We can parameterize the probability function by μ and d

${\displaystyle P(X=x)=\exp \left(-\left(\mu (2-d)+{\frac {\mu (d-1)}{2}}\right)\right)\sum _{j=0}^{[x/2]}{\frac {(\mu (2-d))^{x-2j}\left({\frac {\mu (d-1)}{2}}\right)^{j}}{(x-2j)!j!}}}$

Hence the log-likelihood function is,[9]

${\displaystyle {\begin{aligned}{\mathcal {L}}(x_{1},\ldots ,x_{m};\mu ,d)&=\log({\mathcal {L}}(x_{1},\ldots ,x_{m};\mu ,d))\\&=m\mu \left(-1+{\frac {d-1}{2}}\right)+\log(\mu (2-d))\sum _{i=1}^{m}x_{i}+\sum _{i=1}^{m}\log(q_{i}(\theta ))\end{aligned}}}$

where

• ${\displaystyle q_{i}(\theta )=\sum _{j=0}^{[x_{i}/2]}{\frac {\theta ^{j}}{(x_{i}-2j)!j!}}}$
• ${\displaystyle \theta ={\frac {d-1}{2\mu (2-d)^{2}}}}$

From the log-likelihood function, the likelihood equations are,[9]

${\displaystyle {\frac {\partial l}{\partial \mu }}=m\left(-1+{\frac {d-1}{2}}\right)+{\frac {1}{\mu }}\sum _{i=1}^{m}x_{i}-{\frac {d-1}{2\mu ^{2}(2-d)^{2}}}\sum _{i=1}^{m}{\frac {q_{i}^{'}(\theta )}{q_{i}(\theta )}}}$

${\displaystyle {\frac {\partial l}{\partial d}}=m{\frac {\mu }{2}}-{\frac {\sum _{i=1}^{m}x_{i}}{2-d}}-{\frac {d}{2\mu (2-d)^{3}}}\sum _{i=1}^{m}\sum _{i=1}^{m}{\frac {q_{i}^{'}(\theta )}{q_{i}(\theta )}}}$

Straightforward calculations show that,[9]

• ${\displaystyle \mu ={\bar {x}}}$
• And d can be found by solving,

${\displaystyle \sum _{i=1}^{m}{\frac {q_{i}^{'}({\tilde {\theta }})}{q_{i}({\tilde {\theta }})}}=m({\bar {x}}(2-d))^{2}}$

where ${\displaystyle {\tilde {\theta }}={\frac {d-1}{2{\bar {x}}(2-d)^{2}}}}$

• It can be shown that the log-likelihood function is strictly concave in the domain of the parameters. Consequently, the MLE is unique.

The likelihood equation does not always have a solution like as it shows the following proposition,

Proposition:[9] Let X₁, ..., X_m come from a generalized Hermite distribution with fixed n. Then the MLEs of the parameters are ${\displaystyle {\hat {\mu }}}$ and ${\displaystyle {\tilde {d}}}$ if only if ${\displaystyle m^{(2)}/{\bar {x}}^{2}>1}$, where ${\displaystyle m^{(2)}=\sum _{i=1}^{n}x_{i}(x_{i}-1)/n}$ indicates the empirical factorial momement of order 2.

• Remark 1: The condition ${\displaystyle m^{(2)}/{\bar {x}}^{2}>1}$ is equivalent to ${\displaystyle {\tilde {d}}>1}$ where ${\displaystyle {\tilde {d}}=\sigma ^{2}/{\bar {x}}}$ is the empirical dispersion index
• Remark 2: If the condition is not satisfied, then the MLEs of the parameters are ${\displaystyle {\hat {\mu }}={\bar {x}}}$ and ${\displaystyle {\tilde {d}}=1}$, that is, the data are fitted using the Poisson distribution.

A usual choice for discrete distributions is the zero relative frequency of the data set which is equated to the probability of zero under the assumed distribution. Observing that ${\displaystyle f_{0}=\exp(-(a_{1}+a_{2}))}$ and ${\displaystyle \mu =a_{1}+2a_{2}}$. Following the example of Y. C. Patel (1976) the resulting system of equations,

${\displaystyle {\begin{cases}{\bar {x}}=a_{1}+2a_{2}\\f_{0}=\exp(-(a_{1}+a_{2}))\end{cases}}}$

We obtain the zero frequency and the mean estimator a₁ of ${\displaystyle {\hat {a_{1}}}}$ and a₂ of ${\displaystyle {\hat {a_{2}}}}$,[6]

${\displaystyle {\hat {a_{1}}}=-({\bar {x}}+2\log(f_{0}))}$

${\displaystyle {\hat {a_{2}}}={\bar {x}}+\log(f_{0})}$

where ${\displaystyle f_{0}={\frac {n_{0}}{n}}}$, is the zero relative frequency, n > 0

It can be seen that for distributions with a high probability at 0, the efficiency is high.

• For admissible values of ${\displaystyle {\hat {a_{1}}}}$ and ${\displaystyle {\hat {a_{2}}}}$, we must have

${\displaystyle -\log \left({\frac {n_{0}}{n}}\right)<{\bar {x}}<-2\log \left({\frac {n_{0}}{n}}\right)}$

The likelihood-ratio test statistic [9] for hermite distribution is,

${\displaystyle W=2({\mathcal {L}}(X;{\hat {\mu }},{\hat {d}})-{\mathcal {L}}(X;{\hat {\mu }},1))}$

Where ${\displaystyle {\mathcal {L}}()}$ is the log-likelihood function. As d = 1 belongs to the boundary of the domain of parameters, under the null hypothesis, W does not have an asymptotic ${\displaystyle \chi _{1}^{2}}$ distribution as expected. It can be established that the asymptotic distribution of W is a 50:50 mixture of the constant 0 and the ${\displaystyle \chi _{1}^{2}}$. The α upper-tail percentage points for this mixture are the same as the 2α upper-tail percentage points for a ${\displaystyle \chi _{1}^{2}}$; for instance, for α = 0.01, 0.05, and 0.10 they are 5.41189, 2.70554 and 1.64237.

The score statistic is,[9]

${\displaystyle S_{2}=2m\left[{\frac {m^{(2)}-{\bar {x}}^{2}}{2{\bar {x}}}}\right]^{2}={\frac {m({\tilde {d}}-1)^{2}}{2}}}$

where m is the number of observations.

The asymptotic distribution of the score test statistic under the null hypothesis is a ${\displaystyle \chi _{1}^{2}}$ distribution. It may be convenient to use a signed version of the score test, that is, ${\displaystyle \operatorname {sgn} (m^{(2)}-{\bar {x}}^{2}){\sqrt {S}}}$, following asymptotically a standard normal.