Neyman Type A distribution

In statistics and probability, the Neyman Type A distribution is a discrete probability distribution from the family of Compound Poisson distribution. First of all, to easily understand this distribution we will demonstrate it with the following example explained in Univariate Discret Distributions; we have a statistical model of the distribution of larvae in a unit area of field (in a unit of habitat) by assuming that the variation in the number of clusters of eggs per unit area (per unit of habitat) could be represented by a Poisson distribution with parameter $\lambda$, while the number of larvae developing per cluster of eggs are assumed to have independent Poisson distribution all with the same parameter $\varphi$. If we want to know how many larvae there are, we define a random variable Y as the sum of the number of larvae hatched in each group (given j groups). Therefore, Y = X₁ + X₂ + ... X _j, where X₁,...,X_j are independent Poisson variables with parameter $\lambda$ and $\varphi$.

History

Jerzy Neyman was born in Russia in April 16 of 1894, he was a Polish statistician who spent the first part of his career in Europe. In 1939 he developed the Neyman Type A distribution to describe the distribution of larvae in experimental field plots. Above all, it is used to describe populations based on contagion, e.g., entomology (Beall[1940], Evans[1953]), accidents (Creswell i Froggatt [1963]), and bacteriology.

The original derivation of this distribution was on the basis of a biological model and, presumably, it was expected that a good fit to the data would justify the hypothesized model. However, it is now known that it is possible to derive this distribution from different models (William Feller[1943]), and in view of this, Neyman's distribution derive as Compound Poisson distribution. This interpretation makes them suitable for modelling heterogeneous populations and renders them examples of apparent contagion.

Despite this, the difficulties in dealing with Neyman's Type A arise from the fact that its expressions for probabilities are highly complex. Even estimations of parameters through efficient methods, such as maximum likelihood, are tedious and not easy to understand equations.

Definition

Probability generating function

The probability generating function (pgf) G₁(z), which creates N independent X_j random variables, is used to a branching process. Each X_j produces a random number of individuals, where X₁, X₂,... have the same distribution as X, which is that of X with pgf G₂(z). The total number of  individuals is then the random variable,

$Y=SN=X_1+X_2+...+X_N$

The p.g.f. of the distribution of SN is :

$E[z^{SN}]=E_N[E[z^{SN}|N]]=E_N[G_2(z)]=G_1(G_2(z))$

One of the notations, which is particularly helpful, allows us to use a symbolic representation to refer to  an F1 distribution that has been generalized by an F2 distribution is,

$Y\sim F_1\bigwedge F_N$

In this instance, it is written as,

$Y\sim \mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{\lambda }\mathrm{)}\bigwedge \mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{\varphi }\mathrm{)}$

Finally, the probability generating function is,

$G_Y(z)=\exp (\lambda (e^{\varphi (z-1)}-1))$

From the generating function of probabilities we can calculate the probability mass function explained below.

Probability mass function

Let X₁,X₂,...X_j be Poisson independent variables. The probability distribution of the random variable Y = X₁ +X₂+...X_j is the Neyman's Type A distribution with parameters $\lambda$ and $\varphi$.

$p_x=P(Y=x)=\frac{e^{-\lambda }{\varphi }^x}{x!}\sum _{j=0}^{\mathrm{\infty }}\frac{(\lambda e^{-\varphi }{)}^j{j}^x}{j!}$$x=1,2,...$

Alternatively,

$p_x=P(Y=x)=\frac{e^{-\lambda +\lambda e^{-\varphi }}{\varphi }^x}{x!}\sum _{j=0}^{x}S(x,j){\lambda }^j{e}^{-{\varphi }^j}$

In order to see how the previous expression develops, we must bear in mind that the probability mass function is calculated from the probability generating function, and use the property of Stirling Numbers. Let's see the development

$G(z)=e^{-\lambda +\lambda e^{-\varphi }}\sum _{j=0}^{\mathrm{\infty }}\frac{{\lambda }^j{e}^{-j\varphi }(e^{\varphi z}-1{)}^j}{j|}$

$=e^{-\lambda +\lambda e^{-\varphi }}\sum _{j=0}^{\mathrm{\infty }}(\lambda e^{-\varphi }{)}^j\sum _{x=j}^{\mathrm{\infty }}\frac{S(x,j){\varphi }^x{z}^x}{x!}$

$=\frac{e^{-\lambda +\lambda e^{-\varphi }}{\varphi }^x}{x!}\sum _{j=0}^{x}S(x,j){\lambda }^j{e}^{-{\varphi }^j}$

Another form to estimate the probabilities is with recurring successions,

$p_x=P(Y=x)=\frac{\lambda \varphi e^{-\varphi }}{x}\sum _{r=0}^{x-1}\frac{{\varphi }^r}{r!}{p}_{x-r-1}$,$p_0=\exp (-\lambda +\lambda e^{-\varphi })$

Although its length varies directly with n, this recurrence relation is only employed for numerical computation and is particularly useful for computer applications.

where

• x = 0, 1, 2, ... , except for probabilities of recurring successions, where x = 1, 2, 3, ...
• $j\le x$
• $\lambda$, $\varphi >0$.
• x! and j! are the factorials of x and j, respectively.
• one of the properties of Stirling numbers of the second kind is as follows:

$(e^{\varphi z}-1{)}^j=j!\sum _{x=j}^{\mathrm{\infty }}\frac{S(x,j){\varphi z}^x}{x!}$

Notation

$Y\sim \mathrm{NA}(\lambda ,\varphi )$

Properties

Moment and cumulant generating functions

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 $\mathrm{N}\mathrm{A}\mathrm{(}\mathrm{\lambda },\mathrm{\varphi }\mathrm{)}$, the moment generating function exists and is equal to

$M(t)=G_Y(e^t)=\exp (\lambda (e^{\varphi (e^t-1)}-1))$

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

$K(t)=\log (M(t))=\lambda (e^{\varphi (e^t-1)}-1)$

In the following table we can see the moments of the order from 1 to 4

Order	Moment	Cumulant
1	${\mu }_1=\lambda \varphi$	$\mu$
2	${\mu }_2=\lambda \varphi (1+\varphi )$	${\sigma }^2$
3	${\mu }_3=\lambda \varphi (1+3\varphi +{\varphi }^2)$	$k_3$
4	${\mu }_4=\lambda \varphi (1+7\varphi +6{\varphi }^2+{\varphi }^3)+3{\lambda }^2{\varphi }^2(1+\varphi {)}^2$	$k_4+3{\sigma }^4$

Skewness

The skewness is the third moment centered around the mean divided by the 3/2 power of the standard deviation, and for the $\mathrm{NA}$ distribution is,

${\gamma }_1=\frac{{\mu }_3}{{\mu }_2^{3/2}}=\frac{\lambda \varphi (1+3\varphi +{\varphi }^2)}{(\lambda \varphi (1+\varphi ){)}^{3/2}}$

Kurtosis

The kurtosis is the fourth moment centered around the mean, divided by the square of the variance, and for the $\mathrm{NA}$ distribution is,

${\beta }_2=\frac{{\mu }_4}{{\mu }_2^2}=\frac{\lambda \varphi (1+7\varphi +6{\varphi }^2+{\varphi }^3)+3{\lambda }^2{\varphi }^2(1+\varphi {)}^2}{{\lambda }^2{\varphi }^2(1+\varphi {)}^2}=\frac{1+7\varphi +6{\varphi }^2+{\varphi }^3}{\varphi \lambda (1+\varphi {)}^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,

${\gamma }_2=\frac{{\mu }_4}{{\mu }_2^2}-3=\frac{1+7\varphi +6{\varphi }^2+{\varphi }^3}{\varphi \lambda (1+\varphi {)}^2}$

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

Characteristic function

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

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

This function is related to the moment generating function via ${\varphi }_x(t)=M_X(it)$. Hence for this distribution the characteristic function is,

${\varphi }_x(t)=\exp (\lambda (e^{\varphi (e^{it}-1)}-1))$

• Note that the symbol ${\varphi }_x$ is used to represent the characteristic function.

Cumulative distribution function

The cumulative distribution function is,

$\begin{array}{cc}F(x;\lambda ,\varphi )& =P(Y\le x)\\ & =e^{-\lambda }\sum _{i=0}^x\frac{{\varphi }^i}{i!}\sum _{j=0}^{\mathrm{\infty }}\frac{(\lambda e^{-\varphi }{)}^j{j}^i}{j!}\\ & =e^{-\lambda }\sum _{i=0}^x\sum _{j=0}^{\mathrm{\infty }}\frac{{\varphi }^i(\lambda e^{-\varphi }{)}^j{j}^i}{i!j!}\end{array}$

Other properties

• The index of dispersion is a normalized measure of the dispersion of a probability distribution. It is defined as the ratio of the variance ${\sigma }^2$ to the mean $\mu$,

$d=\frac{{\sigma }^2}{\mu }=1+\varphi$

• From a sample of size N, where each random variable Y_i comes from a $\mathrm{N}\mathrm{A}\mathrm{(}\mathrm{\lambda },\mathrm{\varphi }\mathrm{)}$, where Y₁, Y₂, .., Y_n are independent. This gives the MLE estimator as,

$\sum _{i=1}^N\frac{x_i}{N}=\overline{x}=\overline{\lambda }\overline{\varphi }$ where $\mu$ is the poblational mean of $\overline{x}$

• Between the two earlier expressions We are able to parametrize using $\mu$ and $d$,

$\{\begin{array}{c}\mu =\lambda \varphi \\ d=1+\varphi \end{array}⟶\begin{array}{c}\lambda =\frac{\mu }{d-1}\\ \varphi =d-1\end{array}$

Parameter estimation

Method of moments

The mean and the variance of the NA($\lambda ,\varphi$) are $\mu =\lambda \varphi$ and $\lambda \varphi (1+\varphi )$, respectively. So we have these two equations,

$\{\begin{array}{c}\overline{x}=\lambda \varphi \\ s^2=\lambda \varphi (1+\varphi )\end{array}$

• $s^2$ and $\overline{x}$ are the mostral variance and mean respectively.

Solving these two equations we get the moment estimators $\hat{\lambda }$ and $\hat{\varphi }$ of $\lambda$ and $\varphi$.

$\overline{\lambda }=\frac{{\overline{x}}^2}{s^2-\overline{x}}$

$\overline{\varphi }=\frac{s^2-\overline{x}}{\overline{x}}$

Maximum likelihood

Calculating the maximum likelihood estimator of $\lambda$ and $\varphi$ involves multiplying all the probabilities in the probability mass function to obtain the expression $ℒ(\lambda ,\varphi ;x_1,\dots ,x_n)$.

When we apply the parameterization adjustment defined in "Other Properties," we get $ℒ(\mu ,d;X)$. We may define the Maximum likelihood estimation based on a single parameter if we estimate the $\mu$ as the $\overline{x}$ (sample mean) given a sample X of size N. We can see it below.

$ℒ(d;X)=\prod _{i=1}^{n}P(x_i;d)$

• To estimate the probabilities, we will use the p.m.f. of recurring successions, so that the calculation is less complex.

Testing Poisson assumption

When $\mathrm{N}\mathrm{A}\mathrm{(}\mathrm{\lambda },\mathrm{\varphi }\mathrm{)}$ is used to simulate a data sample it is important to see if the Poisson distribution fits the data well. For this, the following Hypothesis test is used:

$\{\begin{array}{c}{H}_0:d=1\\ H_1:d>1\end{array}$

Likelihood-ratio test

The likelihood-ratio test statistic for $\mathrm{NA}$ is,

$W=2(ℒ(X;\mu ,d)-ℒ(X;\mu ,1))$

Where likelihood $ℒ()$ is the log-likelihood function. W does not have an asymptotic ${\chi }_1^2$ distribution as expected under the null hypothesis since d = 1 is at the parameter domain's edge. In the asymptotic distribution of W, it can be demonstrated that the constant 0 and ${\chi }_1^2$ have a 50:50 mixture. For this mixture, the $\alpha$ upper-tail percentage points are the same as the $2\alpha$ upper-tail percentage points for a ${\chi }_1^2$

Related distributions

The poisson distribution (on { 0, 1, 2, 3, ... }) is a special case of the Neyman Type A distribution, with

$\mathrm{Pois}(\lambda )=\mathrm{NA}(\lambda ,0).$

From moments of order 1 and 2 we can write the population mean and variance based on the parameters $\lambda$ and $\varphi$.

$\mu =\lambda \varphi$

${\sigma }^2=\lambda \varphi (1+\varphi )$

In the dispersion index d we observe that by substituting $\mu$ for the parametrized equation of order 1 and $\sigma$ for the one of order 2, we obtain $d=1+\varphi$. Our variable Y is therefore distributed as a Poisson of parameter $\lambda$ when d approaches to 1.

Then we have that,

$\underset{d\to 1}{\lim }\mathrm{NA}(\lambda ,\varphi )\to \mathrm{Pois}(\lambda )$

Applications

Usage History

When the Neyman's type A species' reproduction results in clusters, a distribution has been used to characterize the dispersion of plants. This typically occurs when a species develops from offspring of parent plants or from seeds that fall close to the parent plant. However, Archibald (1948] observed that there is insufficient data to infer the kind of reproduction from the type of fitted distribution. While Neyman type A produced positive results for plant distributions, Evans (1953] showed that Negative binomial distribution produced positive results for insect distributions. Neyman type A distributions have also been studied in the context of ecology, and the results that unless the plant clusters are so compact as to not lie across the edge of the square used to pick sample locations, the distribution is unlikely to be applicable to plant populations. The compactness of the clusters is a hidden assumption in Neyman's original derivation of the distribution, according to Skellam (1958). The results were shown to be significantly impacted by the square size selection.

In the context of bus driver accidents, Cresswell and Froggatt (1963) derived the Neyman type A based on the following hypotheses:

1. Each driver is susceptible to "spells," the number of which is Poissonian for any given amount of time, with the same parameter $\lambda$ for all drivers.
2. A driver's performance during a spell is poor, and he is likely to experience a Poissonian number of collisions, with the same parameter $\varphi$ for all drivers.
3. Each driver behaves independently.
4. No accidents can occur outside of a spell.

These assumptions lead to a Neyman type A distribution via the $\mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{\lambda }\mathrm{)}$ and $\mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{\varphi }\mathrm{)}$ model. In contrast to their "short distribution," Cresswell and Froggatt called this one the "long distribution" because of its lengthy tail. According to Irwin(1964), a type A distribution can also be obtained by assuming that various drivers have different levels of proneness, or K. with probability:

$\mathrm{P}\mathrm{(}\mathrm{K}\mathrm{=}\mathrm{k}\mathrm{\varphi }\mathrm{)}=\frac{e^{-\lambda }{\lambda }^k}{k!}$

• taking values $0,\varphi ,2\varphi ,3\varphi ,...$

and that a driver with proneness k$\varphi$ has X accidents where:

$\mathrm{P}\mathrm{(}\mathrm{X}\mathrm{=}\mathrm{x}\mathrm{|}\mathrm{k}\mathrm{\varphi }\mathrm{)}=\frac{e^{-k\varphi }(k\varphi {)}^x}{x!}$

This is the $\mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{k}\mathrm{\varphi }\mathrm{)}\wedge \mathrm{P}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{(}\mathrm{\lambda }\mathrm{)}$ model with mixing over the values taken by K.

Distribution was also suggested in the application for the grouping of minority tents for food from 1965 to 1969. In this regard, it was predicted that only the clustering rates or the average number of entities per grouping needed to be approximated, rather than adjusting distributions on d very big data bases.

• The reference of the usage history.

Calculating Neyman Type A probabilities in R

• The code below simulates 5000 instances of $\mathrm{N}\mathrm{A}\mathrm{(}\mathrm{\lambda },\mathrm{\varphi }\mathrm{)}$,

 rNeymanA <- function(n,lambda, phi){
   r <- numeric()
   for (j in 1:n) {
      k = rpois(1,lambda)
      r[j] <- sum(rpois(k,phi)) 
   }
   return(r)
 }

• The mass function of recurring probabilities is implemented in R in order to estimate theoretical probabilities; we can see it below,

dNeyman.rec <- function(x, lambda, phi){
 p <- numeric()
 p[1]<- exp(-lambda + lambda*exp(-phi))
 c <- lambda*phi*exp(-phi)
 if(x == 0){
   return(p[1])
 }
 else{
   for (i in 1:x) {
     suma = 0
     for (r in 0:(i-1)) {
       suma = suma + (phi^r)/(factorial(r))*p[i-r]
     }
     p[i+1] = (c/(i))*suma # +1 per l'R
   }
   res <- p[i+1]
   return(res)
 }
}

We compare results between the relative frequencies obtained by the simulation and the probabilities computed by the p.m.f. . Given two values for the parameters $\lambda =2$ and $\varphi =1$. It is displayed in the following table,

X	0	1	2	3	4	5	6	7	8
Estimated $f_i$	.271	.207	.187	.141	.082	.052	.028	.016	.008
Teoric $P(X=x_i)$	.282	.207	.180	.129	.084	.051	.029	.016	.008