Count Random Variables

Definition

The observation of randomness patterns serves as guidance for the craft of probabilistic modelling. The most used count models—Binomial, Poisson, Negative Binomial—are the discrete Morris’ natural exponential families whose variance is at most quadratic on the mean, and the solutions of Katz–Panjer recurrence relation, aside from being members of the generalised power series and hypergeometric distribution families, and this accounts for their many advantageous characteristics. Some other basic count models are also described, as well as models with less obvious but useful randomness patterns in connection with maximum entropy characterisations, such as Zipf and Good models. Simple tools, such as truncation, thinning, or parameter randomisation, are straightforward ways of constructing other count models.

1. Introduction

For any $\mathcal{S}={\left\{x_k\right\}}_{k\in \mathbb{K}}$, with $\mathbb{K}\subseteq {\mathbb{N}}_0=\left\{0,1,\dots \right\}$, and for any sequence ${\left\{p_k\right\}}_{k\in \mathbb{K}}$ such that $p_k\ge 0$ for any $k\in \mathbb{K}$ and ${\sum }_{x_k\in \mathcal{S}}{p}_k=1$, $X=\left\{\begin{array}{c}{x}_k,x_k\in \mathcal{S}\hfill \\ p_k=\mathbb{P}[X=x_k]\hfill \end{array}\right.$ is a discrete lattice random variable with support $\mathcal{S}$ and probability mass function ${\left\{p_k\right\}}_{k\in \mathbb{K}}$. If $x_k=k\in {\mathbb{N}}_0$, X is a count random variable.

In most cases, the probability mass function ${\left\{p_k\right\}}_{k\in \mathbb{K}}$ is not interesting, since it is difficult to deal with and there is no clear interpretation of the pattern of randomness it describes. The craft of probabilistic modelling (Gani (1986) [1]) uses a diversity of criteria to describe and select models, namely, those arising from randomness patterns (such as counts in Bernoulli trials, sampling with or without replacement, and random draws from urns). Another source of the rationale description of count models are characterisation theorems based on structural properties (e.g., a power series distribution with mean = variance, or maximum Shannon entropy with prescribed arithmetic and/or geometric mean). Recurrence relationships (for instance, $p_{k+1}=(a+{\displaystyle \frac{b}{k+1}})p_k,k=\nu ,\nu +1,\dots$) or mathematical properties (for instance, the variance being at most a quadratic function of the expectation) also define interesting families of discrete random variables. On the other hand, asymptotic properties such as arithmetic properties, namely, infinite divisibility, discrete self-decomposability, and stability, serve as guidance in model choice.

Section 2 describes the discrete uniform random variables, modelling equiprobability patterns resulting from the principle of insufficient reason, of which the Bernoulli random variable with parameter ${\displaystyle \frac{1}{2}}$ is the simplest example. Section 3 is a detailed overview of count models—Binomial, Negative Binomial—arising from the observation of random patterns in Bernoulli trials (including the Poisson random variable, as a limit of Binomial random variables under a mean stability restriction, and the Hypergeometric sampling without replacement model, herein in the context of conditioning on the sum of two independent random variables). In Section 4, the recurrence relation holding for the probability mass function of Binomial, Poisson, and Negative Binomial random variables investigated by Katz [2] and by Panjer [3] is extended to describe Hess et al. [4]’s family of basic count distributions. Section 5 briefly discusses alternative organisations of count models, namely, via Power Series distributions or Kemp’s [5] generalised hypergeometric probability distributions. Section 6 constrasts the equilibrium pattern of Zipf’s [6] law with the equiprobability modelled by the discrete uniform random variables. Section 7 and Section 8 discuss ways of transforming random models, respectively, by randomising parameters and via the discretisation of countinuous random variables. Section 9 briefly discusses the role of characterisations in the craft of modelling count data. Further issues are briefly described in Section 10.

2. Bernoulli Random Variables, Principle of Insufficient Reason and Discrete Uniform Random Variables

Let $p\in [0,1]$ be the probability of the event $\mathcal{E}$ occurring as the outcome of an experiment. Either $\mathcal{E}$ (sometimes referred to as success) occurs once, or its complementary event $\overline{\mathcal{E}}$ (referred to as failure) occurs, so the number of occurrences of $\mathcal{E}$ in a single trial is either 1, with probability p, or 0, with probability $1-p$. We shall use the notation

$$B=\left\{\begin{array}{cc}0& 1\\ 1-p& p\end{array}\right.,$$

where the first line indicates the support of the count random variable B, and the second line the probabilities of the support points.

The above random variable is called “Bernoulli”, with parameter p, in honour of the genial probabilist Jacques Bernoulli, author of the fascinating Ars Conjectandi [7], published posthumously in 1713 by his nephew Niklaus Bernoulli (also a probabilist).

We shall use the notation $X\sim Bernoulli\left(p\right)$; if $p={\displaystyle \frac{1}{2}}$, meaning that $\mathcal{E}$ and $\overline{\mathcal{E}}$ are equiprobable, it should be assumed that there is insufficient reason to assign different probabilities to $\mathcal{E}$ and $\overline{\mathcal{E}}$.

The principle of insufficient reason (renamed principle of indifference by Maynard Keynes [8]) was used in the foundation texts of Jacques Bernoulli [7] and Laplace [9] for assigning epistemic probabilities to equiprobable events. The natural extension of $X\sim Bernoulli\left({\displaystyle \frac{1}{2}}\right)$ is the equiprobable count model

$$Y=\left\{\begin{array}{c}k,k=0,\dots ,n\hfill \\ p_k={\displaystyle \frac{1}{n+1}}\hfill \end{array}\right.,$$

named Discrete Uniform with parameter n, that we denote as $Y\sim DiscreteUniform\left(n\right)$.

3. Count Random Variables in Bernoulli Trials

Let $\left\{T_1,T_2,\dots \right\}$ be a sequence of independent random experiments (trials) whose outcomes are either $\mathcal{E}$, with probability $p\in (0,1)$, or $\overline{\mathcal{E}}$, with probability $1-p$. The assumption of independence means that the outcome of trial $T_k,k=1,2,\dots$ has no influence whatsoever on the outcome of any other trial.

In this setting, the counts that matter are the following:

• The number k of outcomes $\mathcal{E}$ in n trials (n fixed, k random). From the definition of Bernoulli trials, it is easy to conclude that, with X denoting such a count random variable, $\mathbb{P}[X=k]=p_k=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)p^k{(1-p)}^{n-k},k=0,\dots ,n$. This random variable, with probability mass function ${\left\{p_k=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)p^k{(1-p)}^{n-k}\right\}}_{k=0}^n$, is called a binomial random variable with parameters $n\in {\mathbb{N}}^{+}$ and $p\in (0,1)$, and we shall denote it as $X\sim Binomial(n,p)$ or

  $$X=\left\{\begin{array}{c}k,k=0,\dots ,n\hfill \\ p_k=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)p^k{(1-p)}^{n-k}\hfill \end{array}\right..$$

  The expectation is $\mathbb{E}\left(X\right)=np$ and the variance $\mathrm{var}\left(X\right)=np(1-p)$. Hence, the dispersion index is $\mathcal{D}\left(X\right)={\displaystyle \frac{\mathrm{var}\left(X\right)}{\mathbb{E}\left(X\right)}}=1-p<1$. For that reason, we say that the random variable $X\sim Binomial(n,p)$ is underdispersed.
• The number n of trials needed to observe k times the outcome $\mathcal{E}$ (k fixed, n random). The simple case is $k=1$. As in this case, with $W_1^{\ast }$ denoting such a count random variable, $\mathbb{P}[W_1^{\ast }=n]=\mathbb{P}[{\overline{\mathcal{E}}}_1\dots {\overline{\mathcal{E}}}_{n-1}{\mathcal{E}}_n]={(1-p)}^{n-1}p$. $W_1^{\ast }$ is called geometric (or sometimes Pascal) random variable, and we shall use the notation $W_1^{\ast }\sim Geometric\left(p\right)$, or

  $$W_1^{\ast }=\left\{\begin{array}{c}n,n=1,2,\dots \hfill \\ p_n=p{(1-p)}^{n-1}\hfill \end{array}\right..$$

  Its expectation is $\mathbb{E}\left(W_1^{\ast }\right)={\displaystyle \frac{1}{p}}$ and its variance is $\mathrm{var}\left(W_1^{\ast }\right)={\displaystyle \frac{1-p}{p^2}}$.
  More generally, let $W_k^{\ast }$ be the number of trials needed for ${\mathcal{E}}_n$ being the k-th $\mathcal{E}$ of outcomes. Obviously, due to the independence of the Bernoulli trials,

  $$\mathbb{P}(W_k^{\ast }=n)=\mathbb{P}(k-1{\mathcal{E}}^{\prime }\mathrm{s}\mathrm{in}n-1\mathrm{trials})\times \mathbb{P}\left({\mathcal{E}}_n\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{k-1}}\right)p^{k-1}{(1-p)}^{n-k}\times p.$$

  $W_k^{\ast }$ is a Negative Binomial random variable with parameters k and p, that we denote $W_k^{\ast }\sim NegativeBinomial(k,p)$, or

  $$W_k^{\ast }=\left\{\begin{array}{c}n,n=k,k+1,\dots \hfill \\ p_n=\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{k-1}}\right)p^k{(1-p)}^{n-k}\hfill \end{array}\right..$$

  Obviously, ${\displaystyle W_k^{\ast }\stackrel{\mathrm{d}}{=}{W}_{1,1}^{\ast }+\dots +W_{1,k}^{\ast }}$, with the $W_{1,i}^{\ast }\sim Geometric\left(p\right),i=1,\dots ,k$, independent. So $\mathbb{E}\left(W_k^{\ast }\right)={\displaystyle \frac{k}{p}}$ and $\mathrm{var}\left(W_k^{\ast }\right)={\displaystyle \frac{k(1-p)}{p^2}}$.
  It is sometimes convenient to shift the Negative Binomial random variables to start at 0, i.e., to count the number of $\overline{\mathcal{E}}$s that precede the k-th $\mathcal{E}$. In other words, to define

  $$W_k=W_k^{\ast }-k\sim NegativeBinomial(k,p,0),W_k=\left\{\begin{array}{c}n,n=0,1,\dots \hfill \\ p_n=\left({\displaystyle \genfrac{}{}{0pt}{}{n+k-1}{k-1}}\right)p^k{(1-p)}^n\hfill \end{array}\right..$$

  Its dispersion index is $\mathcal{D}\left(W_n\right)={\displaystyle \frac{\mathrm{var}\left(W_n\right)}{\mathbb{E}\left(W_n\right)}}={\displaystyle \frac{\frac{k(1-p)}{p^2}}{k(\frac{1}{p}-1)}}={\displaystyle \frac{1}{p}}>1$, and for that reason the Negative Binomial random variables are considered overdispersed.
  Note also that using the gamma function extension of factorials, $\Gamma \left(\alpha \right)={\int }_0^{\infty }{x}^{\alpha -1}{\mathrm{e}}^{-x}\mathrm{d}x$, $\alpha >0$, satisfying the recurrence relation $\Gamma (\alpha +1)=\alpha \Gamma \left(\alpha \right)$, so that $\Gamma (n+1)=n!$, we may consider Negative Binomial random variables with parameter $\alpha >0,p\in (0,1)$: $W_{\alpha }\sim NegativeBinomial(\alpha ,p,0)=\left\{\begin{array}{c}n,n=0,1,\dots \hfill \\ p_n={\displaystyle \frac{\Gamma (n+\alpha )}{n!\Gamma \left(\alpha \right)}}{p}^{\alpha }{(1-p)}^n\hfill \end{array}\right..$

In many situations, asymptotic results are paramount in modelling decisions or simplifications. The first central limit theorem (a name coined by Pólya [10] in 1920), establishing that if $X\sim Binomial(n,p)$ and n is large, then ${\displaystyle \frac{X-np}{\sqrt{np(1-p})}}$ can be approximated by $Z\sim Gaussian(0,1)$, appeared in the second edition of Abraham de Moivre’s [11] The Doctrine of Chances.

The other important asymptotic result about Binomial sequences is Poisson’s [12] law of rare events:

Let ${\left\{X_n\sim Binomial(n,p_n)\right\}}_{n\in \mathbb{N}}$, mean-stable in the sense that ${\displaystyle \mathbb{E}\left(X_n\right)=n{p}_n\underset{n\to \infty }{⟶}\mu >0}$ (observe that this implies that ${\displaystyle p_n=O\left(\frac{1}{n}\right)\underset{n\to \infty }{⟶}0}$, and this is the rationale for the name “law of rare events”). Then,

$${\displaystyle X_n\underset{n\to \infty }{\overset{\mathrm{d}}{⟶}}Y=\left\{\begin{array}{c}k,k=0,1,\dots \hfill \\ p_k={\mathrm{e}}^{-\mu }\frac{{\mu }^k}{k!}\hfill \end{array}\right.;}$$

Y is a Poisson random variable with parameter $\mu >0$, that we denote $Y\sim Poisson\left(\mu \right)$.

$\mathbb{E}\left(Y\right)=\mathrm{var}\left(Y\right)=\mu$, and hence, in what concerns dispersion, the Poisson random variable is a yardstick, in the sense that its dispersion index is $\mathcal{D}\left(Y\right)=1.$ McCabe and Skeels (2020) [13] and Di Noia et al. (2024) [14] extensively investigated testing Poissoness vs. overdispersion or underdispersion; cf. also Mijburgh and Visagie (2020) [15]’s overview of goodness-of-fit tests for the Poisson distribution.

The Binomial, the Negative Binomial, and the Poisson random variables are the discrete members of Morris’ [16] Natural Exponential Family (NEF) with quadratic variance function in the mean value (QVF). (Recall that X is a member of NEF if its probability density function can be written as $f_X\left(x\right|\theta )=h\left(x\right)exp[\theta x-A\left(\theta \right)]$, and therefore its cumulant generating function has the simple form ${\mathcal{K}}_X\left(t\right)=A(\theta +t)-A\left(\theta \right)$.) In the sequel, Morris [17] treated “topics that can be handled within this unified NEF-QVF formulation, including unbiased estimation, Bhattacharyya and Cramér-Rao lower bounds, conditional distributions and moments, quadratic regression, conjugate prior distributions, moments of conjugate priors and posterior distributions, empirical Bayes and $G_2$ minimax, marginal distributions and their moments, parametric empirical Bayes, and characterisations”, and this shows the relevance of Binomial, Negative Binomial, and Poisson count models in Statistical Inference.

It is also worth mentioning that Binomial, Negative Binomial, and Poisson random variables have relevant additive properties, in the sense that

• If $X\sim Binomial(n,p)$ and $Y\sim Binomial(m,p)$ are independent, then $X+Y\sim Binomial(n+m,p)$.
• If $X\sim NegativeBinomial(k,p)$ and $Y\sim NegativeBinomial(j,p)$ are independent, then $X+Y\sim NegativeBinomial(k+j,p)$.
• If $X\sim Poisson\left(\mu \right)$ and $Y\sim Poisson\left(\theta \right)$ are independent, then $X+Y\sim Poisson(\mu +\theta )$.

Furthermore, if $X\sim Poisson\left(\mu \right)$ is subject to Binomial filtering or thinning, the resulting Y has probability mass function

$$\begin{array}{cc}\hfill \mathbb{P}(Y=k)& =\sum _{j=0}^{\infty }\mathbb{P}(Y=k|X=k+j)\mathbb{P}(X=k+j)\hfill \\ \hfill & =\sum _{j=0}^{\infty }\left({\displaystyle \genfrac{}{}{0pt}{}{k+j}{k}}\right)p^k{(1-p)}^j{\mathrm{e}}^{-\mu }{\displaystyle \frac{{\mu }^{k+j}}{(k+j)!}}={\mathrm{e}}^{-p\mu }{\displaystyle \frac{{\left(p\mu \right)}^k}{k!}},k=0,1,\dots \hfill \end{array}$$

i.e., $Y\sim Poisson\left(p\mu \right)$.

Many natural phenomena are subject to filtering; for instance, using a Poisson model for the number of eggs laid by turtles, the number of hatched eggs, the number of surviving newborns until they reach the ocean, the number of turtles surviving the first year, and so forth, can be modelled by a chain of filtered Poisson random variables, a very useful tool in population dynamics modelling.

On the other hand, the Negative Binomial is a “random Poisson”, obtained when its parameter is a random variable with a Gamma distribution: Let $X|\Lambda \sim \mathit{Poisson}(\Lambda )$ with $\Lambda \sim \mathit{Gamma}(\alpha ,\delta )$. Then,

$$\begin{array}{cc}\hfill \mathbb{P}[X=k]& ={\int }_0^{\infty }{\displaystyle \frac{{\lambda }^k{\mathrm{e}}^{-\lambda }}{k!}}{\displaystyle \frac{{\lambda }^{\alpha -1}{\mathrm{e}}^{-\lambda /\delta }}{{\delta }^{\alpha }\Gamma \left(\alpha \right)}}\mathrm{d}\lambda ={\displaystyle \frac{1}{k!{\delta }^{\alpha }\Gamma \left(\alpha \right)}}{\int }_0^{\infty }{\left({\displaystyle \frac{\delta y}{1+\delta }}\right)}^{\alpha +k-1}{\displaystyle \frac{{\delta \mathrm{e}}^{-y}}{1+\delta }}\mathrm{d}y\hfill \\ \hfill & ={\displaystyle \frac{\Gamma (\alpha +k)}{k!\Gamma \left(\alpha \right)}}{\left({\displaystyle \frac{1}{1+\delta }}\right)}^{\alpha }{\left({\displaystyle \frac{\delta }{1+\delta }}\right)}^k=\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha +k-1}{k}}\right){\left({\displaystyle \frac{1}{1+\delta }}\right)}^{\alpha }{\left({\displaystyle \frac{\delta }{1+\delta }}\right)}^k,\hfill \end{array}$$

k = 0, 1, …, so we obtain $X\sim \mathit{NegativeBinomial}\left(\alpha ,{\displaystyle \frac{1}{1+\delta }}\right)$. For this reason, and as the Negative Binomial random variable is overdispersed, the Poisson being the yardstick, some authors in population dynamics modelling consider the Negative Binomial a more dispersed Poisson.

The above results on the sums of Morris’ discrete variable have interesting consequences when conditioning is applied:

If $X\sim Poisson\left(\mu \right)$ and $Y\sim Poisson\left(\theta \right)$ are independent,

$$\begin{array}{cc}\hfill \mathbb{P}(X=k|X+Y=n)& ={\displaystyle \frac{\mathbb{P}(X=k)\mathbb{P}(Y=n-k)}{\mathbb{P}(X+Y=n)}}\hfill \\ \hfill & ={\displaystyle \frac{{\mathrm{e}}^{-\mu }\frac{{\mu }^k}{k!}{\mathrm{e}}^{-\theta }\frac{{\theta }^{n-k}}{(n-k)!}}{{\mathrm{e}}^{-(\mu +\theta )}\frac{{(\mu +\theta )}^n}{n!}}}=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right){\left({\displaystyle \frac{\mu }{\mu +\theta }}\right)}^k{\left({\displaystyle \frac{\theta }{\mu +\theta }}\right)}^{n-k},k=0,1,\dots ,\hfill \end{array}$$

i.e., $X|X+Y=n\sim Binomial\left(n,{\displaystyle \frac{\mu }{\mu +\theta }}\right)$.

Reasoning in the same way, if $X\sim Binomial(n,p)$ and $Y\sim Binomial(m,p)$ are independent,

$$\begin{array}{cc}\hfill \mathbb{P}(X=k|X+Y=\ell )& ={\displaystyle \frac{\mathbb{P}(X=k)\mathbb{P}(Y=\ell -k)}{\mathbb{P}(X+Y=\ell )}}\hfill \\ \hfill & ={\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n}{k}\right)p^k{(1-p)}^{n-k}\left(\genfrac{}{}{0pt}{}{m}{\ell -k}\right)p^{\ell -k}{(1-p)}^{m-\ell +k}}{\left(\genfrac{}{}{0pt}{}{n+m}{\ell }\right)p^{\ell }{(1-p)}^{n+m-\ell }}}\hfill \\ \hfill & ={\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n}{k}\right)\left(\genfrac{}{}{0pt}{}{m}{\ell -k}\right)}{\left(\genfrac{}{}{0pt}{}{n+m}{\ell }\right)}},k=max\{0,\ell -k\},\dots ,min\{\ell ,n\}.\hfill \end{array}$$

This is the probability mass function of $X|X+Y=\ell \sim Hypergeometric(n+m,\ell ,p)$ with $p={\displaystyle \frac{m}{m+n}}$.

If $X\sim Hypergeometric(N,n,p)$, $\mathbb{E}\left(X\right)=np$ and $\mathrm{var}\left(X\right)=np(1-p){\displaystyle \frac{N-n}{N-1}}$.

It is interesting to observe that if $X_P\sim Poisson\left(\mu \right)$, $Y_B\sim Binomial(n,p)$, and $T_H\sim Hypergeometric(N,n,p)$ have the same expectation $np$, then,

$$\mathrm{var}\left(X_P\right)=\mu =np>\mathrm{var}\left(Y_B\right)=np(1-p)>\mathrm{var}\left(T_H\right)=np(1-p){\displaystyle \frac{N-n}{N-1}}.$$

In other words, for these three count models, the number of parameters implies a tradeoff between decreasing simplicity and increasing information. Augmenting the number of parameters supplies more information. However, one should always bear in mind that all models are wrong, but some are useful (as Box [18] judiciously stated), and that the parsimony principle hints that the simplest model that is useful enough should be used. Observe that in model choice, criteria such as the Akaike Information Criterion (AIC) [19] penalise the number of parameters of the model, and this often leads us to choose a model with fewer parameters, although it provides a better fit.

The Hypergeometric random variable is used in simple random sampling without replacement, while the Binomial is appropriate for simple random sampling with replacement. In the Hypergeometric context, we are dealing with exchangeability; in the Binomial setting, with a stronger independence concept.

The Hypergeometric random model is also the basis for estimating the size of a population with the technique of capture–recapture, cf. Seber and Schofield [20]. Suppose that in a controlled investigation, k animals are captured, marked, and released, and that after a while, there is a second instance of capturing, in which j of the n captured are marked (i.e., recaptured). If the unknown size of the population is N, k of which have been marked and $N-k$ are unmarked, j is the observed value of $Y\sim Hypergeometric(N,n,p)$, with $p={\displaystyle \frac{k}{N}}$. The Peterson estimator $\tilde{N}={\displaystyle \frac{nk}{j}}$ of the size of the population is a method of moments estimator (observe that we described a technique with just one recapture of j marked animals, and in a sample $\left(y\right)$ of size one, this is the mean).

4. Katz Count Models and Extensions

The probability mass function of the discrete Morris NEF-QVF Binomial, Poisson, and Negative Binomial random variables share the interesting property

$$p_{k+1}=\left(a+{\displaystyle \frac{b}{k+1}}\right)p_k,k=0,1,\dots$$

(1)

Let $N_{a,b}$ denote the non-degenerate count random variables whose probability mass function satisfies the recursive relation (1).

Multiplying both sides of (1) by $s^{k+1}$ and adding for $k\ge 0$, we find that the corresponding probability generating functions ${\displaystyle {\mathcal{G}}_{a,b}\left(s\right)=\sum _{k=0}^{\infty }{p}_k{s}^k}$ are the solutions of the differential equation ${\displaystyle \frac{{{\mathcal{G}}^{\prime }}_{a,b}\left(s\right)}{{\mathcal{G}}_{a,b}\left(s\right)}}={\displaystyle \frac{a+b}{1-as}}$; for details, cf. Pestana and Velosa (2004) [21]; namely:

• ${\mathcal{G}}_{a,b}\left(s\right)={\left(1-{\displaystyle \frac{a}{a-1}}+{\displaystyle \frac{a}{a-1}}s\right)}^{-{\displaystyle \frac{a+b}{a}}}$, $a<0,{\displaystyle \frac{a}{a-1}}=p\in (0,1),-{\displaystyle \frac{a+b}{a}}=n\in \mathbb{N}$, i.e., $N_{a,b}\sim Binomial\left(-{\displaystyle \frac{a+b}{a}},{\displaystyle \frac{a}{a-1}}\right)$.
• ${\mathcal{G}}_{a,b}\left(s\right)={\mathrm{e}}^{b(s-1)}$, $a=0,b>0$, i.e., $N_{0,b}\sim Poisson\left(b\right)$.
• ${\mathcal{G}}_{a,b}\left(s\right)={\left({\displaystyle \frac{1-a}{1-as}}\right)}^{{\displaystyle \frac{a+b}{a}}}$, $a\in (0,1),a+b>0$, i.e., $N_{a,b}\sim NegativeBinomial\left({\displaystyle \frac{a+b}{a}},a\right)$.

The simple expression (1) for the successive probability atoms has been several times rediscovered in different contexts (McCabe and Skeels, 2020 [13]). Katz (1945, 1965) [2,22] used it to organise a family of discrete models in the same spirit of the continuous Pearson family, and for that reason, the family of random variables $\left\{N_{a,b}\right\}$ whose probability mass function satisfies (1) is referred to as the Katz family. In Risk Theory and Insurance, $\left\{N_{a,b}\right\}$ is known as Panjer’s class, due to Panjer’s (1981) [3] important breakthrough that the recurrence relation (1) implies that the distribution of aggregate claims can be iteratively computed or approximated:

Let $S_1,S_2,\dots$ denote the identically distributed discrete claim sizes, with ${\pi }_k=\mathbb{P}[S_i=k],$$k\in {\mathbb{N}}_0,$ and assume that the $S_i$ are mutually independent and independent of the claim number $N_{a,b}$. Then,

$$\mathbb{E}\left(S_1\left|\sum _{i=1}^n{S}_i=j\right.\right)={\displaystyle \frac{j}{n}},\mathrm{and}\mathbb{P}\left(S_1=k\left|\sum _{i=1}^n{S}_i=j\right.\right)={\displaystyle \frac{{\pi }_k{\pi }_{j-k}^{\ast (n-1)}}{{\pi }_j^{\ast n}}},$$

where $\left\{{\pi }_k^{\ast n}\right\}$ denotes the n-th convolution of $\left\{{\pi }_k\right\}$. Consequently, the probabilities $p_{X,k}=\mathbb{P}[X=j]$ of the aggregate claim amount ${\displaystyle X=\sum _{i=1}^{N_{a,b}}{S}_i}$ (understood as $X=0$ if $N_{a,b}=0$) can be obtained recursively using Panjer’s algorithm:

$$\mathbb{P}[X=j]=\left\{\begin{array}{cc}{\mathcal{G}}_{a,b}\left({\pi }_0\right),\hfill & j=0\hfill \\ \\ {\displaystyle \frac{{\sum }_{k=1}^j\left(a+\frac{bk}{j}\right){\pi }_k{p}_{X,j-k}}{1-a{\pi }_0}},\hfill & j\in \mathbb{N}\hfill \end{array}\right..$$

This requires only $\mathrm{O}\left(j\right)$ computations, while the traditional convolution method would require $\mathrm{O}\left(j^2\right)$. More generally, the cumulative distribution function $F_X={\sum }_{k=0}^{\infty }{p}_k{F}_S^{\ast k}$ of the aggregate claim $X={\sum }_{k=1}^{N_{a,b}}{S}_k$ for an arbitrary distribution function $F_S$ of the individual claim sizes, when the number of claims is $N_{a,b}$ with probability mass $\left\{p_k\right\}$, satisfies the integral equation

$$F_X\left(x\right)=p_0+a{F}_S\ast F_X\left(x\right)+b{\int }_0^{x}v{\int }_0^{x-v}{\displaystyle \frac{\mathrm{d}{F}_X\left(y\right)}{v+y}}\mathrm{d}{F}_S\left(v\right)$$

for $x>0$, if $F_S\left(0\right)=0$. For nonnegative claim sizes with $F_S\left(0\right)>0$,

$$F_X\left(0\right)=\left\{\begin{array}{cc}{\mathrm{e}}^{b[F_S\left(0\right)-1]},& a=0\\ {\left[{\displaystyle \frac{1-a}{1-a{F}_S\left(0\right)}}\right]}^{{\displaystyle \frac{a+b}{a}}},& a\ne 0\end{array}\right..$$

For detailed proofs, cf. Rólski et al. [23] (pp. 118–124). Klugman et al. [24] (pp. 221–224) use an example to discuss in depth the use of the empiric dispersion index as guidance for model choice in the Panjer family.

If we relax the condition on the support, $p_{k+1}=(a+{\displaystyle \frac{b}{k+1}})p_k,k=1,2,\dots$, as investigated by Jewell [25], Willmot [26], or Sundt [27], we also obtain the Logarithmic random variable

$$X_L=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ p_k={\displaystyle \frac{{\theta }^k/k}{-ln(1-\theta )}}\hfill \end{array}\right.,\theta \in (0,1),$$

used by Fisher [28] to model species abundance, and the Engen Extended Negative Binomial (ENB) [29] random variable,

$$X_{ENB}=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ {\displaystyle p_k=\frac{1}{1-{\left(1-\theta \right)}^{\beta }}\left(\genfrac{}{}{0pt}{}{\beta +k-1}{\beta -1}\right){\theta }^k{\left(1-\theta \right)}^{\beta }}\hfill \end{array}\right.,\theta \in (0,1),\beta \in (-1,0).$$

The dispersion indices of the Logarithmic and of the Engen random variables depend on $\theta$. For $X_L\sim Logarithmic\left(\theta \right)$,

$$\mathcal{D}\left(X_L\right)={\displaystyle \frac{-\frac{{\theta }^2+\theta ln(1-\theta )}{{(1-p)}^2{[ln(1-\theta )]}^2}}{-\frac{\theta }{(1-\theta )ln(1-\theta )}}}={\displaystyle \frac{\theta +ln(1-\theta )}{(1-\theta )ln(1-\theta )}}⟹\mathcal{D}\left(X_L\right)\left\{\begin{array}{cc}<1\hfill & \mathrm{if}\theta \in \left(0,1-{\displaystyle \frac{1}{\mathrm{e}}}\right)\hfill \\ =1\hfill & \mathrm{if}\theta =1-{\displaystyle \frac{1}{\mathrm{e}}}\hfill \\ >1\hfill & \mathrm{if}\theta \in \left(1-{\displaystyle \frac{1}{\mathrm{e}}},1\right)\hfill \end{array}\right.,$$

and for $X_{ENB}\sim Engen(\beta ,\theta )$

$$\mathcal{D}\left(X_{ENB}\right)={\displaystyle \frac{1+\beta \theta -\frac{\beta \theta }{1-{(1-\theta )}^{\beta }}}{(1-\theta )}}⟹\mathcal{D}\left(X_E\right)\left\{\begin{array}{cc}<1\hfill & \mathrm{if}\theta \in \left(0,1-{\left({\displaystyle \frac{1}{\beta +1}}\right)}^{1/\beta }\right)\hfill \\ =1\hfill & \mathrm{if}\theta =1-{\left({\displaystyle \frac{1}{\beta +1}}\right)}^{1/\beta }\hfill \\ >1\hfill & \mathrm{if}\theta \in \left(1-{\left({\displaystyle \frac{1}{\beta +1}}\right)}^{1/\beta },1\right)\hfill \end{array}\right..$$

Hess et al. (2002) [4] investigated $\nu$-Panjer classes, $\nu \in {\mathbb{N}}_0$, whose probability mass function satisfies $p_{k+1}=\left(a+{\displaystyle \frac{b}{k+1}}\right)p_k,k=\nu ,\nu +1,\dots$. Excluding the 0-Panjer family, which does not contain the logarithmic and ENB distributions, Hess et al. (2002) [4] established that any $\nu$-Panjer distribution is the left endpoint truncation of an $\left(\nu -1\right)$-Panjer distribution, for $\nu \ge 2$. For this reason, they call the Binomial, Poisson, Negative Binomial, Logarithmic, and ENB distributions basic count models.

Klugman et al.’s [24] notation $N(a,b,n)$ has been used by Fackler (2024) [30] to supply a unified, practical, and intuitive representation of the Panjer distributions and their parameter space, and give an inventory of parameterisations used for Panjer distributions.

Panjer’s recursion has been extended by Tzaninis and Bozikas (2024) [31] to mixed compound distributions, following the use of finite mixtures for modelling dispersion in count data by Ong et al. (2023) [32].

5. Other Families of Count Models

The Katz [2,22] organisation of discrete models using the difference Equation (1), similar to Pearson’s continuous model types description using a differential equation, inspired later essays on the rationale description of discrete models.

5.1. Power Series Distributions

The discrete random variable $X_{\theta },\theta >0$ has a power series distribution (PSD) if

$$X_{\theta }=\left\{\begin{array}{c}k,k\in {\mathbb{N}}_0\hfill \\ p_k={\displaystyle \frac{a_k{\theta }^k}{A\left(\theta \right)}}\hfill \end{array}\right.,a_k\ge 0,A\left(\theta \right)=\sum _{k\in {\mathbb{N}}_0}{a}_k{\theta }^k.$$

A power series distribution is a discrete linear exponential distribution.

If we replace ${\mathbb{N}}_0$ by $\mathcal{S}$, a nonempty enumerable subset of ${\mathbb{N}}_0$, we obtain a generalised power series distribution (GPSD):

$$X_{\theta }=\left\{\begin{array}{c}k,k\in S\hfill \\ p_k={\displaystyle \frac{a_k{\theta }^k}{A\left(\theta \right)}}\hfill \end{array}\right.,a_k\ge 0,A\left(\theta \right)=\sum _{k\in S}{a}_k{\theta }^k.$$

Observe that the probability generating function of $X_{\theta }$ is ${\mathcal{G}}_{X_{\theta }}\left(s\right)={\displaystyle \frac{A\left(\theta s\right)}{A\left(\theta \right)}}$.

The random variables $X_B\sim Binomial(n,\theta )$, with $A\left(\theta \right)={(1+\theta )}^n$, $Y_P\sim Poisson\left(\theta \right)$, with $A\left(\theta \right)={\mathrm{e}}^{\theta }$, $W_{NB}\sim NegativeBinomial(k,\theta )$, with $A\left(\theta \right)={\displaystyle \frac{1}{{(1-\theta )}^k}},k>0$, and $X_L\sim Logarithmic\left(\theta \right),\theta \in (0,1)$, with $A\left(\theta \right)=-ln(1-\theta )$, have GPSD.

Obviously, if a generalised power series distribution is truncated, another generalised power series distribution is obtained. An important issue, seldom discussed in the literature: restricting the support may enlarge the parameter space $\Theta$; for example truncating to the right of $N\in \mathbb{N}$ the support of $X_L\sim Logarithmic\left(\theta \right),\theta \in (0,1)$, the resulting truncated logarithmic random variable $X_{TL}$

$$X_{TL}=\left\{\begin{array}{c}k,k=1,\dots ,N\hfill \\ p_k={\displaystyle \frac{{\theta }^k/k}{{\sum }_{k=1}^N{\theta }^k/k}}\hfill \end{array}\right.,\theta >0,$$

has a parameter space $\Theta =(0,\infty )$, since ${\sum }_{k=1}^N{\displaystyle \frac{{\theta }^k}{k}}<\infty$ for any $\theta >0$.

5.2. Kemp’s Generalised Hypergeometric Probability Distributions

The work by Johnson et al. (2005) [33] is a coordinated thorough description of useful discrete random variables, structured on Kemp’s (1968) [5] use of hypergeometric functions, and the corresponding characterisation sections enhance their randomness patterns. Kemp’s [5] generalised hypergeometric probability distributions, with probability generating function ${\displaystyle \frac{{}_{p}F_q\left(\lambda z\right)}{{}_{p}F_q\left(\lambda \right)}}$, where

$${}_{p}F_q[a_1,\dots ,a_p;b_1,\dots ,b_q;\lambda z]=\sum _{j=0}^{\infty }{\displaystyle \frac{{\left(a_1\right)}_j\dots {\left(a_p\right)}_j}{{\left(b_1\right)}_j\dots {\left(b_q\right)}_j}}{\displaystyle \frac{{\left(\lambda z\right)}^j}{j!}}$$

(${\left(\alpha \right)}_j=\alpha (\alpha +1)\cdots (\alpha +j-1)$ is the Pochhammer’s symbol) is the generalised hypergeometric function with numeratorial parameters $a_1,\dots a_p$ and denominatorial parameters $b_1,\dots ,b_q$. This huge family—Johnson et al. (2005) [33] (pp. 93–94) identify the name of 37 of the more than 70 studied in Patil et al. (1984) [34]—has many useful properties.

6. Uniformity and Power Law Randomness, a Striking Example of Opposite Patterns

The simplest equiprobability pattern of count randomness is modelled by the discrete uniform distribution $X\sim DiscreteUniform\left(n\right)$, $X=\left\{\begin{array}{c}k,k=0,1,\dots n\hfill \\ p_k={\displaystyle \frac{1}{n+1}}\hfill \end{array}\right..$

Equiprobability is, however, the exception; natural phenomena are more prone to present other patterns of equilibrium, such as with Zipf’s model, which is the special case $q=0$ of the Zipf–Mandelbrot family of discrete models given by

$$X_{N,q,\alpha }=\left\{\begin{array}{c}k,k\in \left\{1,2,\dots ,N\right\}\hfill \\ p_k={\displaystyle \frac{1}{H_{N,q,\alpha }}}{\displaystyle \frac{1}{{(k+q)}^{\alpha }}}\hfill \end{array}\right.,$$

where $\alpha >0,q\ge 0$, k is the rank of the data and $H_{N,q,\alpha }={\displaystyle \sum _{j=1}^N}{\displaystyle \frac{1}{{(j+q)}^{\alpha }}}$ is a generalisation of the harmonic number $H_N={\displaystyle \sum _{j=1}^N}{\displaystyle \frac{1}{j}}$. The cumulative distribution function is, for $x\in \{1,\dots ,N\}$, $F\left(x\right)={\displaystyle \frac{H_{\lfloor x\rfloor ,q,\alpha }}{H_{N,q,\alpha }}}$, where $\lfloor x\rfloor$ is the largest integer not greater than x.

Zipf’s law has been introduced to model the equilibrium—in the sense that rank × frequency is approximately constant—in verbal communication, which employs mainly words of ordinary usage (social trend) and scarcely words tied to the user’s vocabulary and thematic preferences (individual trend). In what concerns extra-pair paternity (EPP), there is some evidence that the observed variability of the number of extra-pair offspring (EPO) in the brood, in various species of passerines, results from a delicate equilibrium of needing a social partner to care for raising the brood and the tendency to have offspring sired by a stronger male. This is the reason why the Zipf–Mandelbrot family seems to be a plausible model, cf. Marques et al. (2005) [35].

Notice that Zipf’s law is also a particular case of the extended parameter space of the right-truncated logarithmic distribution where $p_k={\displaystyle \frac{{\theta }^k/k}{ln(1-\theta )}},$ for $k=1,2,\dots$ and $\theta \in (0,1).$ In fact, when considering the truncated logarithmic model where $p_k={\displaystyle \frac{{\theta }^k/k}{{\sum }_{j=1}^{\nu }{\theta }^j/j}}$ for $k=1,2,\dots ,\nu$, the parameter space can be extended to $\theta \in (0,\infty )$, as already observed. In particular, if $\theta =1$, we obtain Zipf’s law. For more information, cf., e.g., Johnson et al. (2005) [33] (Chapter 11) and references therein. It is a power law in the sense that $\mathbb{P}[X=k]\propto k^{-\alpha }$, for some positive $\alpha$, and therefore its log-log plot exhibits a linear signature. When an almost linear signature is observed, it is considered a strong support for choosing a power law model, namely, in the discrete context, a Zipf–Mandelbrot law.

7. Hierarchical Models and Parameter Randomisation

A powerful tool to obtain models is to randomise parameters. For instance, let $X|P\sim \mathit{Binomial}(n,P)$ and $P\sim \mathit{Uniform}(0,1)$. Since

$$\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right){\int }_0^1{p}^k{(1-p)}^{n-k}\mathrm{d}p=\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)B(k+1,n+1-k)={\displaystyle \frac{1}{n+1}},$$

we have $\mathbb{P}[X=k]={\displaystyle \frac{1}{n+1}}$, $k=0,1,\dots ,n$, i.e., $X\sim \mathit{DiscreteUniform}\left(n\right)$. Observe that $P\sim \mathit{Uniform}(0,1)$ can be considered as “null” information about the value of p, thus resulting in discrete equiprobability.

More generally, if $X|P\sim \mathit{Binomial}(n,P)$ with $P\sim \mathit{Beta}(\mu ,\nu )$, we obtain Beta-Binomial distributions, with probability mass function

$$\begin{array}{cc}\hfill p_k=\mathbb{P}[X=k]& =\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right){\int }_0^1{p}^k{(1-p)}^{n-k}{\displaystyle \frac{p^{\mu -1}{(1-p)}^{\nu -1}}{B(\mu ,\nu )}}\mathrm{d}p\hfill \\ \hfill & =\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right){\displaystyle \frac{\Gamma (k+\mu )\Gamma (n-k+\nu )\Gamma (\mu +\nu )}{\Gamma (n+\mu +\nu )\Gamma \left(\mu \right)\Gamma \left(\nu \right)}},k=0,1,\dots ,n.\hfill \end{array}$$

Let $X|P\sim \mathit{Geometric}(P)$ with $P\sim \mathit{Uniform}(0,1)$. As ${\int }_0^{1}p{(1-p)}^{k-1}\mathrm{d}p=B(2,k)$, for $k=1,2,\dots$ we obtain $\mathbb{P}[X=k]={\displaystyle \frac{1}{k(k+1)}}$.

More generally, if $X|P\sim \mathit{Geometric}(P)$ with $P\sim \mathit{Beta}(\mu ,\nu )$,

$$p_k={\int }_0^{1}p{(1-p)}^{k-1}{\displaystyle \frac{p^{\mu -1}{(1-p)}^{\nu -1}}{B(\mu ,\nu )}}\mathrm{d}p={\displaystyle \frac{\mu \Gamma (\mu +\nu )\Gamma (\nu +k-1)}{\Gamma \left(\nu \right)\Gamma (\mu +\nu +k)}}.$$

For example, if $\mu =2$ and $\nu =1$, for $k=1,2,\dots$ we obtain

$$p_k={\displaystyle \frac{2\Gamma \left(k\right)\Gamma \left(3\right)}{\Gamma (k+3)}}={\displaystyle \frac{4}{k(k+1)(k+2)}}.$$

Let $X|P\sim \mathit{Geometric}(P)$ with $P\sim \mathit{BetaBoop}(1,1,1,2)$; cf. Brilhante and Pestana [36]. Using formula 1, ${\int }_0^1{x}^{p-1}{(1-x^r)}^{q-1}(-lnx)\mathrm{d}x={\displaystyle \frac{1}{r^2}}B\left({\displaystyle \frac{p}{r}},q\right)\left[\psi \left({\displaystyle \frac{p}{r}}+q\right)-\psi \left({\displaystyle \frac{p}{r}}\right)\right]$, from subsubsection 4.253 in Gradshteyn and Ryzhik (2007) [37] (p. 538), where $\psi \left(z\right)={\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}ln\Gamma \left(z\right)={\displaystyle \frac{{\Gamma }^{\prime }\left(z\right)}{\Gamma \left(z\right)}}$ is the logarithmic derivative of the gamma function, we conclude that [37],

$$p_k={\int }_0^{1}p{(1-p)}^{k-1}(-lnp)\mathrm{d}p=B(2,k)\left[\psi (2+k)-\psi \left(2\right)\right].$$

Observing that taking the derivative of $\Gamma (z+1)=z\Gamma \left(z\right)$, ${\Gamma }^{\prime }(z+1)=z{\Gamma }^{\prime }\left(z\right)+\Gamma \left(z\right)$, and dividing by $\Gamma (z+1)$, we obtain ${\displaystyle \frac{{\Gamma }^{\prime }(z+1)}{\Gamma (z+1)}}={\displaystyle \frac{z{\Gamma }^{\prime }\left(z\right)}{z\Gamma \left(z\right)}}+{\displaystyle \frac{\Gamma \left(z\right)}{z\Gamma \left(z\right)}},$ and consequently $\psi (z+1)=\psi \left(z\right)+{\displaystyle \frac{1}{z}}$, we deduce that

$$p_k={\displaystyle \frac{1}{k(k+1)}}\left[\psi \left(k\right)+{\displaystyle \frac{1}{k}}+{\displaystyle \frac{1}{k+1}}-\psi \left(2\right)\right].$$

Continuing to apply the recursive expression to $\psi$, we obtain

$$X=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ p_k={\displaystyle \frac{H_{k+1}-1}{k(k+1)}}\hfill \end{array}\right.,$$

where $H_n={\sum }_{j=1}^n{\displaystyle \frac{1}{j}}$ is the n-th harmonic number.

8. Discrete Analogues of Continuous Random Variables

Discrete analogues of continuous random variables may inherit interesting patterns. For instance, if we define

$$W=\left\{\begin{array}{c}k,k=0,1,\dots \hfill \\ \mathbb{P}[W=k]=F_X(k+1)-F_X\left(k\right)\hfill \end{array}\right.,$$

with $X\sim Exponential\left(\lambda \right)$, then $W\sim Geometric(1-{\mathrm{e}}^{-\lambda })$ inherits the memoryless property of the Exponential in the sense that $\mathbb{P}[W\ge k+\ell |W\ge \ell ]=\mathbb{P}[W\ge k]$.

In this context, it is worth citing the discrete lognormal random variable:

$$X=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ {\displaystyle p_k=\frac{1}{kC(\mu ,\sigma )}exp\left[-\frac{1}{2}{\left(\frac{lnk-\mu }{\sigma }\right)}^2\right]}\hfill \end{array}\right.$$

is a model for discrete skew populational data, that is not a power law but imitates a linear signature. For an overview of discrete analogues of continuous probability distributions, cf. Chakraborty (2015) [38].

9. Characterisations of Count Models

Families of discrete distributions have striking structural properties, with relevant consequences for model choice. Characterisation theorems are also helpful in the craft of probabilistic modelling.

As already observed, the dispersion index $\mathcal{D}\left(X\right)={\displaystyle \frac{\mathrm{Var}\left(X\right)}{\mathbb{E}\left(X\right)}}$ of the Poisson distribution is 1. In fact, a PSD random variable X has $\mathcal{D}\left(X\right)=1$ if and only if $X\sim Poisson\left(\mu \right)$, for some $\mu >0$, a characterisation due to Kosambi (1949) [39].

Characterisations are a useful tool in model choice. For instance, Kapur (1989) [40] has shown that a discrete random variable X with fixed arithmetic mean has maximum Shannon entropy $-{\sum }_{k\in S}{p}_{k}ln{p}_k$ if and only if X is a Geometric random variable. Observe that choosing the largest entropy model within a class of distributions amounts to selecting as default the least informative model, and so minimising the prior information.

A discrete random variable X with fixed geometric mean has maximum Shannon entropy if and only if it is a Zipf or Discrete Pareto random variable

$$X_Z=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ p_k={\displaystyle \frac{1}{\zeta \left(\eta \right)k^{\eta }}}\hfill \end{array}\right.,$$

where $\zeta \left(\eta \right)$ is the Riemann zeta function. A discrete random variable X with fixed arithmetic and geometric means has maximum Shannon entropy if and only if X is a Good random variable

$$X_{Go}=\left\{\begin{array}{c}k,k=1,2,\dots \hfill \\ p_k={\displaystyle \frac{q^k}{k^{\eta }\varphi (q,\eta ,1)}}\hfill \end{array}\right.,$$

where $\varphi (q,\eta ,\nu )={\sum }_{k=0}^{\infty }{\displaystyle \frac{q^k}{{(\nu +k)}^{\eta }}},\nu \ne 0,-1,-2,\dots$, is the Lerch function (subsubsection 9.55 in Gradshteyn and Ryzhik (2007) [37]). Observe that the Geometric random variable with support $1,2,\dots$ is the special case when $\eta =0$, and the Logarithmic random variable is the case $\eta =1$.

10. Concluding Remarks

As in the continuous case, asymptotic results also play an important role in modelling; for instance, the Poisson limit in discrete settings is, to a certain extent, comparable to the central limit Gaussian law in general; cf. Steutel and van Harn (1979) [41]. In the same spirit, extremal discrete laws do arise naturally in order statistics contexts; cf. Hall (1996) [42]. In addition, useful models can originate other patterns of randomness via truncation (with the eventual broadening of the parameter space), randomly stopped sums, or randomly stopped extremes, and the randomisation of parameters.

Moreover among all continuous probability distributions with support $[0,\infty )$ and mean $\mu$, the exponential distribution with $\lambda =1/\mu$, modelling interarrival times in the Poisson process, has the largest differential entropy

$${\displaystyle h\left(f\right)=-{\int }_{-\infty }^{\infty }{f}_X\left(x\right)ln{f}_X\left(x\right)\mathrm{d}x,}$$

meaning that the Poisson process has the greatest $\lambda t$-dimensional entropy in the interval $(0,t)$ among all the homogeneous point processes with the same given intensity $\lambda$. For that reason, Poisson streams of events are generally interpreted as representing the pattern of unconstrained randomness. This is a rationale for considering Poisson modelling, since entropy is non-decreasing, and therefore the Poisson random variable is a weak limit in very general frameworks, and in particular of Binomial and of Negative Binomial sequences with a mean stability condition. In a sense, the Poisson distribution models the unconstrained scattering of the type of throwing rice grains when seeding.

In Population Dynamics, dispersion is a clue to model choice. Zhang et al. (2018) [43] investigated dispersion patterns of finite mixtures, and Cahoy et al. (2021) [44] introduced flexible fractional Poisson distributions to model underdispersed and overdisersed count data. Recent advances in what concerns underdispersed count models were made by Huang (2022) [45] and Seck et al. (2022) [46]. Rana et al. (2023) [47] investigated the influence of outliers on overdispersion, and Sengupta and Roy (2023) [48] the role of under-reported count data and zero-inflated models, an issue also investigated by Aswi et al. (2022) [49].

From the additive properties of Binomial random variables shown in Section 3, $X\sim Binomial(n,p)$ is n-divisible, in the sense that it can be decomposed as the sum of n independent and identically distributed $B_j\sim Bernoulli\left(p\right),j=1,\dots ,n,$ random variables; and $X\sim NegativeBinomial(k,p)$ is k-divisible, in the sense that it can be decomposed as the sum of k independent and identically distributed $W_j\sim Geometric\left(p\right),j=1,\dots ,k,$ random variables. On the other hand $X\sim Poisson\left(\mu \right)$ is infinitely divisible, since for any $n=2,3,\dots$, it can be decomposed as the sum of n independent and identically distributed $W_k\sim Poisson(\mu /n)$ random variables.

Moreover, from Raikov’s theorem [50], if $X\sim Poisson\left(\mu \right)$ is decomposed as the sum of two independent random variables, $X=Y+Z$, then Y and Z are also Poisson random variables. This means that the Poisson random variables are extreme points of the convex set of infinitely divisible random variables. A similar result with the Gaussian random variables, conjectured by Lévy [51] and proved by Cramér [52], shows that they are also extreme points of the set of infinitely divisible random variables. From this, Johansen (1966) [53]) supplied a structural proof of the integral representations of infinitely divisible characteristics functions. In fact, the Poisson random variables are the building blocks of infinitely divisible random variables (de Finetti (1929, 1932) [54,55], which is a strong modelling asset, since many observed phenomena are the result of several contributing effects.