A Probability Generating Function Based Goodness-of-Fit Test for the Poisson–Three-Parameter Lindley Distribution

1. Introduction

Goodness-of-fit (GOF) testing plays a fundamental role in statistical modeling, particularly in the analysis of count data, where the selection of an appropriate discrete distribution is essential for reliable statistical inference. Classical goodness-of-fit procedures for discrete models, most notably Pearson’s chi-square test, are known to exhibit several well-documented shortcomings, including strong sensitivity to cell grouping, unreliable behavior in small and moderate samples, and additional difficulties introduced by parameter estimation; see Fisher [1], Cochran [2], Moore [3], Read and Cressie [4], and Cressie [5]. These limitations become more pronounced when dealing with flexible count models designed to accommodate overdispersion or unobserved heterogeneity.

Over the past decades, several authors have proposed goodness-of-fit tests based on probability generating functions (PGFs) for discrete distributions. Foundational contributions are due to Rueda, O’Reilly and Pérez-Abreu [6], who developed Cramér-von Mises type tests for the Poisson distribution based on the empirical PGF, and to Henze [7,8], who studied empirical-distribution-function tests and generating function based tests for a broad class of discrete models. Subsequent developments include tests for the bivariate and multivariate Poisson distributions [9,10] and extensions using minimum-distance PGF estimators [see, e.g., [11,12,13,14]. PGF-based tests are known to possess favorable theoretical and empirical properties, particularly in terms of robustness to parameter estimation and improved power against relevant alternatives, often outperforming classical chi-square based procedures in finite samples.

Among flexible count models, Poisson mixture distributions constitute a well-established class that naturally accommodates overdispersion and unobserved heterogeneity; see Johnson, Kemp and Kotz [15]. In this context, the Poisson–Lindley distribution, originally introduced by Sankaran [16], and several of its extensions have been proposed as flexible alternatives with applications in insurance, reliability, and biological sciences; see Ghitany et al. [17]. Extended versions involving additional parameters have been introduced to enhance modeling flexibility while preserving analytical tractability; see Mahmoudi and Zakerzadeh [18] and Das et al. [19].

Despite its growing use, no formal goodness-of-fit test specifically tailored to the Poisson–Three-Parameter Lindley (PTPL) distribution is currently available. Existing applications typically rely on information criteria or informal diagnostics tools, which do not provide a rigorous inferential framework. This constitutes a clear methodological gap that the present paper aims to address.

The main contribution of this paper is a new goodness-of-fit test for the PTPL distribution based on the Cramér–von Mises distance between the empirical PGF and the theoretical PGF of the PTPL model. The closed-form expression of the PTPL PGF makes this approach computationally and analytically tractable. Parameter estimation is handled via maximum likelihood, and the null distribution of the test statistic is approximated through a parametric bootstrap procedure.

Theoretical properties, including consistency against fixed alternatives and bootstrap validity, are rigorously established. Extensive Monte Carlo simulations assess the empirical size and power of the test against competing alternatives, including the Poisson, Negative Binomial, COM-Poisson [20], and Zero-Inflated Poisson distributions.

The remainder of the paper is organized as follows. Section 2 introduces notation and preliminaries. Section 3 formally defines the PTPL distribution and establishes its basic probabilistic properties. Section 4 describes parameter estimation and its asymptotic properties. The proposed test statistic is presented in Section 5. Asymptotic results are given in Section 6. The bootstrap implementation is described in Section 7. The Monte Carlo study and real data application appear in Section 8 and Section 9, respectively. Section 10 is dedicated to discussion and Section 11 concludes.

2. Notation and Preliminaries

Throughout the paper, the following notation is used. Convergence almost surely, in probability, and in distribution are denoted by $\stackrel{\mathrm{a}.\mathrm{s}.}{\to }$, $\stackrel{p}{\to }$, and $\stackrel{d}{\to }$, respectively. The notation ${\mathbb{N}}_0=\{0,1,2,\dots \}$ is used to denote the set of non-negative integers.

Throughout the paper, a sequence of random variables $X_1,\dots ,X_n$ is said to be independent and identically distributed (i.i.d.) if the $X_i$ are mutually independent and each $X_i$ has the same distribution. The abbreviation i.i.d. is used systematically in what follows.

For a sequence of random functions $\left\{Z_n\left(t\right)\right\}$, uniform convergence over a set $\mathcal{T}$ is understood in the sense that

$$\underset{t\in \mathcal{T}}{sup}|Z_n\left(t\right)-Z\left(t\right)|\stackrel{\mathrm{a}.\mathrm{s}.}{\to }0\mathrm{or}\stackrel{p}{\to }0,$$

almost surely or in probability, depending on the context.

Expectations, variances, and probabilities under the true data-generating process are denoted by $\mathbb{E}\left(·\right)$, $\mathbb{V}\left(·\right)$, and $\mathbb{P}\left(·\right)$, respectively. Bootstrap probabilities are denoted by ${\mathbb{P}}^{*}\left(·\right)$.

The notation ${\mathbb{E}}_{\vartheta }\left(·\right)$, ${\mathbb{V}}_{\vartheta }\left(·\right)$, and ${\mathbb{P}}_{\vartheta }\left(·\right)$ denotes expectation, variance, and probability computed under the distribution $\mathrm{PTPL}\left(\vartheta \right)$. When $\vartheta ={\vartheta }_0$ is the true parameter vector, these are written ${\mathbb{E}}_{{\vartheta }_0}\left(·\right)$, ${\mathbb{V}}_{{\vartheta }_0}\left(·\right)$, and ${\mathbb{P}}_{{\vartheta }_0}\left(·\right)$, respectively; the subscript is omitted when the context is unambiguous.

For a sequence of random variables $\left\{Y_n\right\}$, $Y_n=o_p\left(1\right)$ denotes convergence in probability to zero as $n\to \infty$, and $Y_n=O_p\left(1\right)$ denotes stochastic boundedness, i.e., for every $\epsilon >0$ there exists $M_{\epsilon }<\infty$ such that $\mathbb{P}\left(\right|Y_n|>M_{\epsilon })<\epsilon$ for all n. More generally, $Y_n=o_p\left(r_n\right)$ and $Y_n=O_p\left(r_n\right)$ refer to these properties applied to $Y_n/r_n$ for a deterministic sequence $r_n>0$; see, e.g., van der Vaart [21], Chapter 2.

The symbol ${\mathbb{P}}_n=n^{-1}{\sum }_{j=1}^n{\delta }_{X_j}$ denotes the empirical measure of the sample $X_1,\dots ,X_n$, where ${\delta }_x$ is the Dirac measure at x. Similarly, $P_{\vartheta }$ denotes the probability measure on ${\mathbb{N}}_0$ induced by $\mathrm{PTPL}\left(\vartheta \right)$.

3. The Poisson–Three-Parameter Lindley Distribution

The Lindley distribution, originally introduced as a one-parameter model (see, e.g., Ghitany et al. [17]), has been widely generalized to increase modeling flexibility. Among these generalizations, Shanker et al. [22] proposed a three-parameter Lindley distribution that incorporates additional shape parameters while preserving a relatively simple functional form.

In parallel, Poisson–Lindley type distributions have attracted considerable attention as flexible alternatives to the classical Poisson model capable of capturing overdispersion and unobserved heterogeneity; see Mahmoudi and Zakerzadeh [18] and Das et al. [19].

Motivated by these developments, we study the following Poisson–Three-Parameter Lindley distribution which arises as a reparameterization combining ideas from the three-parameter Lindley framework of Shanker et al. [22] and the Poisson compounding structure of Das et al. [19].

Proposition 1.  (Validity of the PTPL PMF).

Let $\Theta =\{(\theta ,\alpha ,\beta )\in {\mathbb{R}}^3:\theta >0,\beta >0,\alpha \ge -\beta /(1+\theta )\}$. For any $(\theta ,\alpha ,\beta )\in \Theta$, the function

$$p(x;\theta ,\alpha ,\beta )={\displaystyle \frac{{\theta }^2}{(\alpha \theta +\beta ){(1+\theta )}^{x+1}}}\left[\alpha +{\displaystyle \frac{\beta (x+1)}{1+\theta }}\right],x=0,1,2,\dots ,$$

(1)

defines a proper probability mass function on ${\mathbb{N}}_0$.

That is, $p(x;\theta ,\alpha ,\beta )\ge 0$ for all $x\in {\mathbb{N}}_0$, and ${\sum }_{x=0}^{\infty }p(x;\theta ,\alpha ,\beta )=1.$

Proof. 

Non-negativity. For $(\theta ,\alpha ,\beta )\in \Theta$, we have

${\theta }^2>0$,

$(\alpha \theta +\beta )>0$ (since $\alpha \ge -\beta /(1+\theta )$ implies $\alpha \theta \ge -\beta \theta /(1+\theta )>-\beta$),

${(1+\theta )}^{x+1}>0$, and

$$\alpha +\beta (x+1)/(1+\theta )\ge -\beta /(1+\theta )+\beta (x+1)/(1+\theta )=\beta x/(1+\theta )\ge 0\mathrm{for}\mathrm{all}x\ge 0.$$

Hence $p(x;\theta ,\alpha ,\beta )\ge 0$.

Normalization. Set $r=1/(1+\theta )\in (0,1)$.

Using the standard identities ${\sum }_{x=0}^{\infty }{r}^{x+1}=r/(1-r)$ and ${\sum }_{x=0}^{\infty }(x+1)r^{x+1}=r/{(1-r)}^2$, (see, e.g., [23]), we obtain

$${\displaystyle \sum _{x=0}^{\infty }}p(x;\theta ,\alpha ,\beta )={\displaystyle \frac{{\theta }^2}{\alpha \theta +\beta }}\left[\alpha ·{\displaystyle \frac{r}{1-r}}+{\displaystyle \frac{\beta }{1+\theta }}·{\displaystyle \frac{r}{{(1-r)}^2}}\right].$$

Substituting $r=1/(1+\theta )$ so $1-r=\theta /(1+\theta )$, gives

$$={\displaystyle \frac{{\theta }^2}{\alpha \theta +\beta }}\left[{\displaystyle \frac{\alpha }{\theta }}+{\displaystyle \frac{\beta }{{\theta }^2}}\right]={\displaystyle \frac{{\theta }^2}{\alpha \theta +\beta }}·{\displaystyle \frac{\alpha \theta +\beta }{{\theta }^2}}=1.$$

   □

A discrete random variable X following distribution (1) is denoted $X\sim \mathrm{PTPL}(\theta ,\alpha ,\beta )$.

Remark 1.  (Special cases).

When $\alpha =1$ and $\beta =0^{+}$ (in the limit), the PTPL reduces to the Poisson–Lindley distribution of Sankaran [16]. When $\alpha =\beta$, additional simplifications arise. These connections motivate the use of PTPL as a flexible encompassing model.

Probability generating function.

The PGF of $X\sim \mathrm{PTPL}(\theta ,\alpha ,\beta )$, $G\left(t\right)=\mathbb{E}\left(t^X\right)$, is given by

$$G(t;\theta ,\alpha ,\beta )={\displaystyle \frac{{\theta }^2\left[\beta +\alpha (1+\theta -t)\right]}{(\alpha \theta +\beta ){(1+\theta -t)}^2}},t\in [0,1].$$

(2)

This expression follows by summing ${\sum }_{x=0}^{\infty }{t}^{x}p\left(x\right)$ using the same geometric series identities employed in Proposition 1. The PGF is well defined and finite on $[0,1]$ since $1+\theta -t\ge \theta >0$ for all $t\le 1$ and $\theta >0$. The closed-form expression (2) is the key analytical feature that makes the PTPL distribution particularly amenable to PGF-based goodness-of-fit testing.

Remark 2.  (Overdispersion).

A direct calculation using the PGF (2) yields

$$\begin{array}{cc}\hfill \mathbb{E}\left(X\right)& =G^{\prime }(1;\theta ,\alpha ,\beta )={\displaystyle \frac{\alpha \theta +2\beta }{\theta (\alpha \theta +\beta )}},\hfill \\ \hfill \mathbb{V}\left(X\right)& =G^{\prime \prime }(1;\theta ,\alpha ,\beta )+G^{\prime }(1;\theta ,\alpha ,\beta )-{\left[G^{\prime }(1;\theta ,\alpha ,\beta )\right]}^2\hfill \\ \hfill & ={\displaystyle \frac{\alpha \theta +2\beta }{\theta (\alpha \theta +\beta )}}+{\displaystyle \frac{2\beta (\alpha \theta +3\beta )}{{\theta }^2(\alpha \theta +\beta )}}-{\left[{\displaystyle \frac{\alpha \theta +2\beta }{\theta (\alpha \theta +\beta )}}\right]}^2.\hfill \end{array}$$

For $\beta >0$ and $(\theta ,\alpha ,\beta )\in \Theta$, one can verify that

$$\mathbb{V}\left(X\right)-\mathbb{E}\left(X\right)={\displaystyle \frac{2\beta (\alpha \theta +3\beta )}{{\theta }^2(\alpha \theta +\beta )}}-{\left[{\displaystyle \frac{\alpha \theta +2\beta }{\theta (\alpha \theta +\beta )}}\right]}^2>0,$$

confirming that the PTPL is always overdispersed relative to the Poisson model throughout its parameter space Θ.

Figure 1 displays the probability mass function of the PTPL distribution for nine representative parameter configurations, organized by increasing level of overdispersion. The dispersion index $\mathrm{DI}=\mathbb{V}\left(X\right)/\mathbb{E}\left(X\right)$ ranges from approximately $1.5$ (mild overdispersion, top row) to values exceeding 40 (pronounced right skewness, bottom row). The figure illustrates the considerable shape flexibility of the PTPL family: depending on the parameter values, the distribution can resemble a near-Poisson shape, exhibit moderate right skewness with a visible mode away from zero, or display a steeply decreasing profile concentrated near the origin. This breadth of shapes, combined with the closed-form expressions for the PMF and PGF derived above, underscores both the practical relevance of the PTPL model for overdispersed count data and the importance of having a formal goodness-of-fit test that is specifically calibrated to this family.

4. Parameter Estimation for the PTPL Distribution

Let $X_1,\dots ,X_n$ be a random sample from a Poisson–Three-Parameter Lindley distribution with parameter vector $\vartheta ={(\theta ,\alpha ,\beta )}^{\top }$. In this section, we describe two standard estimation procedures for $\vartheta$, namely the method of moments and maximum likelihood estimation, and discuss their asymptotic properties. These results provide the theoretical justification for the plug-in approach adopted in the proposed goodness-of-fit test.

4.1. Method of Moments

Let ${\mu }_r=\mathbb{E}\left(X^r\right)$, $r=1,2,3$, denote the theoretical moments of the PTPL distribution, which are smooth functions of the parameter vector $\vartheta$. Let ${\widehat{\mu }}_r=n^{-1}{\sum }_{i=1}^n{X}_i^r$ be the corresponding sample moments. The method of moments estimator ${\widehat{\vartheta }}_n^{\mathrm{MM}}$ is defined as a solution of the system

$${\widehat{\mu }}_r={\mu }_r\left(\vartheta \right),r=1,2,3,$$

provided that the solution exists and is unique.

Such estimators are widely used for count models due to their simplicity and computational stability; see, for example, Casella and Berger [24].

Proposition 2 

(Consistency of the Method of Moments Estimator). Assume that the mapping $\vartheta \mapsto ({\mu }_1,{\mu }_2,{\mu }_3)$ is injective (i.e., the PTPL distribution is identifiable from its first three moments; see Teicher [25] for general identifiability conditions for Poisson mixtures) and continuous on Θ. Then,

$${\widehat{\vartheta }}_n^{\mathrm{MM}}\stackrel{\mathrm{a}.\mathrm{s}.}{\to }{\vartheta }_0,n\to \infty .$$

Proof. 

By the strong law of large numbers [23], ${\widehat{\mu }}_r\to {\mu }_r\left({\vartheta }_0\right)$ almost surely for $r=1,2,3$. Consistency then follows from the assumed continuity and injectivity of the moment mapping by a standard argument; see van der Vaart [21], Theorem 5.9.    □

Due to the numerical instability and absence of closed-form solutions for the moment equations, the method of moments estimator is used exclusively to provide suitable initial values for the maximum likelihood optimization. Its asymptotic properties are therefore stated for completeness but are not the primary focus of the estimation section.

4.2. Maximum Likelihood Estimation

Given the numerical limitations of the method of moments discussed in the previous section, parameter estimation for the PTPL distribution is primarily conducted via maximum likelihood. Let $\left(p\right(x;\vartheta )$ denote the PTPL probability mass function. The log-likelihood function based on a random sample $\{X_1,\dots ,X_n\}$ is given by

$${\ell }_n\left(\vartheta \right)={\displaystyle \sum _{i=1}^n}logp(X_i;\vartheta ),$$

and the maximum likelihood estimator (MLE), denoted by ${\widehat{\vartheta }}_n^{\mathrm{ML}}$, is defined as the maximizer of ${\ell }_n\left(\vartheta \right)$ over the parameter space $\Theta$, provided that the solution exists and is unique.

Proposition 3 

(Consistency of the MLE). Assume that:

(i) 

The PTPL model is correctly specified, i.e., the true distribution of X belongs to $\left\{p\right(·;\vartheta ):\vartheta \in \Theta \}$.

(ii) 

Identifiability: $\vartheta \ne {\vartheta }^{\prime }$ implies $p(·;\vartheta )\ne p(·;{\vartheta }^{\prime })$ on ${\mathbb{N}}_0$.

(iii) 

The parameter space Θ is compact.

(iv) 

The log-likelihood $logp(x;\vartheta )$ is upper semi-continuous in ϑ for each $x\in {\mathbb{N}}_0$.

(v) 

${\mathbb{E}}_{{\vartheta }_0}[{sup}_{\vartheta \in \Theta }logp(X;\vartheta )]<\infty$.

Then,

$${\widehat{\vartheta }}_n^{\mathrm{ML}}\stackrel{\mathrm{a}.\mathrm{s}.}{\to }{\vartheta }_0,n\to \infty .$$

Proof. 

Under conditions (i)–(v), the result follows from the classical consistency theorem for maximum likelihood estimators of Wald [26]; see also van der Vaart [21], Theorem 5.14, for a modern treatment.    □

Proposition 4 

(Asymptotic Normality of the MLE). In addition to the assumptions of Proposition 3, assume that:

(vi) 

The true parameter ${\vartheta }_0$ lies in the interior of Θ.

(vii) 

The log-likelihood $logp(x;\vartheta )$ is twice continuously differentiable in ϑ on a neighbourhood of ${\vartheta }_0$, for all $x\in {\mathbb{N}}_0$. (This holds for the PTPL PMF (1) by inspection of its explicit form.)

(viii) 

The Fisher information matrix $\mathcal{I}\left({\vartheta }_0\right)={\mathbb{E}}_{{\vartheta }_0}[{\nabla }_{\vartheta }logp(X;{\vartheta }_0){\nabla }_{\vartheta }logp{(X;{\vartheta }_0)}^{\top }]$ is finite and non-singular.

Then,

$$\sqrt{n}\left({\widehat{\vartheta }}_n^{\mathrm{ML}}-{\vartheta }_0\right)\stackrel{d}{\to }\mathcal{N}\left(\mathbf{0},{\mathcal{I}}^{-1}\left({\vartheta }_0\right)\right).$$

Proof. 

Under conditions (i)–(viii), the result is a direct consequence of the asymptotic normality theorem for maximum likelihood estimators; see van der Vaart [21], Theorem 5.41, or Lehmann and Casella [27], Chapter 6. Moreover, under these conditions the MLE satisfies the linear influence-function expansion

$$\sqrt{n}({\widehat{\vartheta }}_n^{\mathrm{ML}}-{\vartheta }_0)={\displaystyle \frac{1}{\sqrt{n}}}{\displaystyle \sum _{j=1}^n}\mathcal{I}{\left({\vartheta }_0\right)}^{-1}{\nabla }_{\vartheta }logp(X_j;{\vartheta }_0)+o_p\left(1\right),$$

which is the expansion required in Assumption A1 below.    □

4.3. Implications for the Goodness-of-Fit Test

Both estimation procedures yield consistent estimators of the PTPL parameters. In particular, the MLE is root-n consistent and asymptotically normal, which justifies its use as the plug-in estimator in the proposed test statistic. These properties also support the validity of the parametric bootstrap procedure, as established in Proposition 8 of Section 6.

In the Monte Carlo study, convergence of the MLE was monitored via the gradient norm and the Hessian condition number. Samples for which the optimization did not converge (convergence code $\ne 0$ in the optim routine) were discarded and replaced, and the corresponding discarded fraction is reported in Table 2 for each sample size and parameter configuration. This approach is standard practice in bootstrap-based simulation studies; see, e.g., Henze [8] and Baringhaus and Henze [13].

5. The Proposed Goodness-of-Fit Test

Let $X_1,\dots ,X_n$ be a random sample of independent and identically distributed observations taking values in ${\mathbb{N}}_0$, drawn from an unknown discrete distribution. Under the null hypothesis, the underlying distribution is assumed to belong to the Poisson–Three-Parameter Lindley family.

The null hypothesis of interest is given by

$$H_0:X\sim \mathrm{PTPL}(\theta ,\alpha ,\beta ),$$

for some unknown parameter vector $\vartheta ={(\theta ,\alpha ,\beta )}^{\top }$ belonging to the admissible parameter space, versus the general alternative

$$H_1:X\mathrm{does}\mathrm{not}\mathrm{follow}\mathrm{a}\mathrm{PTPL}\mathrm{distribution}.$$

Let ${\widehat{\vartheta }}_n$ be a consistent estimator of $\vartheta$, obtained, for instance, by maximum likelihood or by the method of moments. The corresponding theoretical probability generating function under the null hypothesis is then given by

$$G(t;{\widehat{\vartheta }}_n)={\displaystyle \frac{{\widehat{\theta }}_n^2\left({\widehat{\beta }}_n-{\widehat{\alpha }}_{n}t+{\widehat{\alpha }}_n+{\widehat{\alpha }}_n{\widehat{\theta }}_n\right)}{({\widehat{\alpha }}_n{\widehat{\theta }}_n+{\widehat{\beta }}_n){\left(1+{\widehat{\theta }}_n-t\right)}^2}},t\in [0,1].$$

The empirical probability generating function is defined as

$${\widehat{G}}_n\left(t\right)={\displaystyle \frac{1}{n}}{\displaystyle \sum _{j=1}^n}{t}^{X_j},t\in [0,1],$$

which provides a natural nonparametric estimator of the theoretical PGF.

In the goodness-of-fit procedure, the PGF is evaluated on the unit interval $t\in [0,1]$, which is contained in the domain $t<log(1+\theta )$ for all admissible $\theta >0$.

Under the null hypothesis, the empirical PGF ${\widehat{G}}_n\left(t\right)$ is expected to be close to the theoretical PGF $G(t;{\widehat{\vartheta }}_n)$ for all $t\in [0,1]$. Departures from the null model can therefore be detected by measuring discrepancies between these two functions over a suitable range of t-values.

Motivated by classical Cramér–von Mises type statistics and following the general PGF-based goodness-of-fit framework developed in the literature, we define the test statistic

$$T_n=n{\int }_0^1{\left[{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)\right]}^{2}w\left(t\right)dt,$$

(3)

where $w\left(·\right)$ is a nonnegative weight function on $[0,1]$ such that

$${\int }_0^{1}w\left(t\right)dt<\infty .$$

This condition ensures that the integral in (3) is finite for each fixed sample size n.

Typical choices of w used in the literature include polynomial weights such as $w\left(t\right)=t^a$, $a>-1$; see, e.g., Henze [8] and Novoa-Muñoz and Jiménez-Gamero [9].

Large values of the statistic $T_n$ provide evidence against the null hypothesis, indicating that the PTPL model does not adequately describe the observed data. Since the asymptotic null distribution of $T_n$ depends on the unknown parameters and involves complex limiting processes, critical values are obtained via a parametric bootstrap procedure, which is described in detail in Section 7.

Remark 3.  (Domain of evaluation).

The PGF is evaluated on $[0,1]$. This is valid because for all $t\le 1$ and $\theta >0$, we have $1+\theta -t\ge \theta >0$, ensuring that the denominator in (2) is bounded away from zero and $G(t;\theta ,\alpha ,\beta )$ is well defined and continuously differentiable on $[0,1]$. The choice $w\left(t\right)=t^a$, $a>-1$, is motivated by the analytical tractability documented in Appendix A and is standard in the PGF-based testing literature [6,8,9]. In the simulation study we consider $a\in \{0,1,2\}$; the case $a=1$ serves as the primary recommendation based on the power comparisons reported in Section 8.

Remark 4.  (Recommended weight exponent).

Among the values $a\in \{0,1,2\}$ considered in the simulation study, $a=1$ consistently yields the best balance between empirical size control and power across all parameter configurations and sample sizes examined. This is consistent with findings for related PGF-based tests; see Henze [8] and Novoa-Muñoz and Jiménez-Gamero [9]. The choice $a=1$ is therefore recommended as the default in practice.

6. Asymptotic Properties

In this section, we investigate the asymptotic behavior of the goodness-of-fit statistic $T_n$ defined in (3). In particular, we establish the uniform consistency of the empirical probability generating function, the behavior of the test statistic under the null hypothesis, and the consistency of the proposed test against fixed alternatives. Finally, we justify the use of a parametric bootstrap procedure for critical value approximation.

Proposition 5 

(Uniform Consistency of the Empirical PGF). Let $X_1,\dots ,X_n$ be i.i.d. random variables supported on ${\mathbb{N}}_0$ with probability generating function $G\left(t\right)$. Then,

$$\underset{t\in [0,1]}{sup}\left|{\widehat{G}}_n\left(t\right)-G\left(t\right)\right|\stackrel{\mathrm{a}.\mathrm{s}.}{\to }0,n\to \infty .$$

Proof. 

For each fixed $t\in [0,1]$, the random variables $t^{X_1},\dots ,t^{X_n}$ are i.i.d., bounded by 1, and satisfy $\mathbb{E}\left(t^{X_j}\right)=G\left(t\right)$. Hence, by the strong law of large numbers [23,28], ${\widehat{G}}_n\left(t\right)\to G\left(t\right)$ almost surely for every $t\in [0,1]$. Moreover, the class of functions $\{x\mapsto t^x:t\in [0,1]\}$ is Glivenko–Cantelli and uniformly bounded; see van der Vaart and Wellner [29]. Therefore,

$$\underset{t\in [0,1]}{sup}\left|{\widehat{G}}_n\left(t\right)-G\left(t\right)\right|\stackrel{\mathrm{a}.\mathrm{s}.}{\to }0,$$

which completes the proof.    □

Proposition 6 

(Behavior Under the Null Hypothesis). Assume that $H_0$ holds, that is, the data are generated from a PTPL distribution with true parameter vector ${\vartheta }_0$. Let ${\widehat{\vartheta }}_n$ be a consistent estimator of ${\vartheta }_0$. Then,

$$T_n\stackrel{p}{\to }0,n\to \infty .$$

Proof. 

By Proposition 5,

$$\underset{t\in [0,1]}{sup}\left|{\widehat{G}}_n\left(t\right)-G(t;{\vartheta }_0)\right|\stackrel{\mathrm{a}.\mathrm{s}.}{\to }0.$$

Moreover, consistency of ${\widehat{\vartheta }}_n$ and continuity of $\vartheta \mapsto G(t;\vartheta )$ (which is $C^{\infty }$ in $\vartheta$ for each fixed $t\in [0,1]$, as seen from (2)), imply

$$\underset{t\in [0,1]}{sup}\left|G(t;{\widehat{\vartheta }}_n)-G(t;{\vartheta }_0)\right|\stackrel{p}{\to }0.$$

Therefore, by the triangle inequality,

$$\underset{t\in [0,1]}{sup}|{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)|\le \underset{t}{sup}|{\widehat{G}}_n\left(t\right)-G(t;{\vartheta }_0)|+\underset{t}{sup}|G(t;{\widehat{\vartheta }}_n)-G(t;{\vartheta }_0)|\stackrel{p}{\to }0,$$

since both ${\widehat{G}}_n\left(t\right)$ and $G(t;{\widehat{\vartheta }}_n)$ are bounded by one on $[0,1]$.

Since the integrand in (3) satisfies ${[{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)]}^{2}w\left(t\right)\le 4w\left(t\right)$ for all $t\in [0,1]$ (both PGFs are bounded by 1), and ${\int }_0^{1}4w\left(t\right)dt<\infty$ by assumption, the constant function $h\left(t\right)=4w\left(t\right)$ serves as an integrable dominating function independent of n. The dominated convergence theorem [30] then yields $T_n\to 0$ in probability. Similar arguments are used in the analysis of PGF-based goodness-of-fit tests; see, for example, Rueda, O’Reilly and Pérez-Abreu [6].    □

The following theorem establishes the weak limit of the centered and scaled process underlying $T_n$, providing the theoretical foundation for the bootstrap approximation discussed in Proposition 8.

Assumption 1. 

The true parameter vector ${\vartheta }_0={({\theta }_0,{\alpha }_0,{\beta }_0)}^{\top }$ belongs to the interior of Θ. The mapping $\vartheta \mapsto G(·;\vartheta )$ is twice continuously Fréchet-differentiable on $L^2([0,1],w)$, where $L^2([0,1],w)$ denotes the space of square-integrable functions with weight $w\left(t\right)=t^a$, $a>-1$. The Fisher information matrix $\mathcal{I}\left({\vartheta }_0\right)$ is finite and non-singular (conditions (vi)–(viii) of Proposition 4). The estimator ${\widehat{\vartheta }}_n$ satisfies the linear expansion

$$\sqrt{n}({\widehat{\vartheta }}_n-{\vartheta }_0)={\displaystyle \frac{1}{\sqrt{n}}}{\displaystyle \sum _{j=1}^n}{\ell }_{{\vartheta }_0}\left(X_j\right)+o_p\left(1\right),$$

where ${\ell }_{{\vartheta }_0}\left(x\right)=\mathcal{I}{\left({\vartheta }_0\right)}^{-1}{\nabla }_{\vartheta }logp(x;{\vartheta }_0)$ is the influence function of the MLE. Under the conditions of Proposition 4, this expansion holds automatically for the MLE; see van der Vaart [21], Theorem 5.41.

Theorem 1.  (Weak limit of the standardized process under $H_0$).

Under $H_0$ and Assumption 1, the process

$${\mathcal{W}}_n\left(t\right)=\sqrt{n}\left[{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)\right],t\in [0,1],$$

converges weakly in $L^2([0,1],w)$ to a centered Gaussian process $\mathcal{W}$ with covariance kernel

$$K(s,t)=G(min(s,t);{\vartheta }_0)-G(s;{\vartheta }_0)G(t;{\vartheta }_0)-\dot{G}{(s;{\vartheta }_0)}^{\top }\mathcal{I}{\left({\vartheta }_0\right)}^{-1}\dot{G}(t;{\vartheta }_0),$$

where $\dot{G}(t;{\vartheta }_0)={\nabla }_{\vartheta }G(t;{\vartheta }_0)$ denotes the gradient of the theoretical PGF with respect to the parameter vector.

Consequently,

$$T_n=n{\int }_0^1{\left[{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)\right]}^{2}w\left(t\right)\mathrm{d}t\stackrel{d}{⟶}{\int }_0^1\mathcal{W}{\left(t\right)}^{2}w\left(t\right)\mathrm{d}t={\displaystyle \sum _{j=1}^{\infty }}{\lambda }_j{Z}_j^2,$$

where $\left\{Z_j\right\}$ are independent standard normal random variables and $\left\{{\lambda }_j\right\}$ (with ${\sum }_{j=1}^{\infty }{\lambda }_j<\infty$) are the eigenvalues of the integral operator associated with $K(·,·)$ in $L^2([0,1],w)$.

Proof. 

Sketch.

Step 1: Donsker property. The function class $\mathcal{F}=\{x\mapsto t^x:t\in [0,1]\}$ is P-Donsker for any distribution P on ${\mathbb{N}}_0$: since $|t^x|\le 1$ for all $x\in {\mathbb{N}}_0$ and $t\in [0,1]$, the class is uniformly bounded, and its metric entropy satisfies ${logN(\epsilon ,\mathcal{F},\parallel ·\parallel }_{L^2\left(P\right)})=O\left({\epsilon }^{-1}\right)$; see van der Vaart and Wellner [29], Chapter 2.5. Therefore, by the Donsker theorem [29], the empirical PGF process $\sqrt{n}[{\widehat{G}}_n\left(t\right)-G(t;{\vartheta }_0)]$ converges weakly in ${\ell }^{\infty }\left([0,1]\right)$ to a centred Gaussian process $\mathbb{G}$ with covariance $\mathrm{Cov}(\mathbb{G}\left(s\right),\mathbb{G}\left(t\right))=G(min(s,t);{\vartheta }_0)-G(s;{\vartheta }_0)G(t;{\vartheta }_0)$.

Step 2: Linearisation of the plug-in correction. By a first-order Taylor expansion of $G(t;{\widehat{\vartheta }}_n)$ around ${\vartheta }_0$, and Assumption A1,

$$\sqrt{n}[G(t;{\widehat{\vartheta }}_n)-G(t;{\vartheta }_0)]=\dot{G}{(t;{\vartheta }_0)}^{\top }\sqrt{n}({\widehat{\vartheta }}_n-{\vartheta }_0)+o_p\left(1\right),$$

uniformly in $t\in [0,1]$, where $\dot{G}(t;{\vartheta }_0)={\nabla }_{\vartheta }G(t;{\vartheta }_0)$.

Step 3: Combination. The combination of Steps 1 and 2 via the continuous mapping theorem yields ${\mathcal{W}}_n\Rightarrow \mathcal{W}$ in $L^2([0,1],w)$, where $\mathcal{W}$ is the centred Gaussian process with covariance kernel K stated above.

Step 4: Spectral representation. The integral operator $\mathcal{K}$ with kernel $K(·,·)$ on $L^2([0,1],w)$ is compact, self-adjoint, and positive semidefinite: compactness follows from square-integrability of K on ${[0,1]}^2$, self-adjointness from symmetry $K(s,t)=K(t,s)$, and positive semidefiniteness from K being a covariance kernel. Its trace satisfies ${\sum }_{j=1}^{\infty }{\lambda }_j={\int }_0^{1}K(t,t)w\left(t\right)dt\le {\int }_0^{1}w\left(t\right)dt<\infty$, so the eigenvalue series converges. By the spectral theorem for compact self-adjoint operators and Mercer’s theorem [29,30], the Karhunen–Loève decomposition gives ${\int }_0^1\mathcal{W}{\left(t\right)}^{2}w\left(t\right)dt\stackrel{d}{=}{\sum }_{j=1}^{\infty }{\lambda }_j{Z}_j^2$, where $\left\{Z_j\right\}$ are i.i.d. standard normal and $\left\{{\lambda }_j\right\}$ are the eigenvalues of $\mathcal{K}$. See Rueda, O’Reilly and Pérez-Abreu [6] and Novoa-Muñoz and Jiménez-Gamero [10] for analogous derivations in related settings.    □

Remark 5. 

The limiting distribution in Theorem 1 depends on the unknown parameter ${\vartheta }_0$ through the eigenvalues $\left\{{\lambda }_j\right\}$, so critical values cannot be tabulated analytically. This is the primary motivation for using the parametric bootstrap, whose consistency is established in Proposition 8.

Proposition 7 

(Consistency Against Fixed Alternatives). Suppose the true distribution of X does not belong to the PTPL family and has probability generating function $G_{\mathrm{true}}\left(t\right)$. Assume that:

(i) 

${\int }_0^1{\left|G_{\mathrm{true}}\left(t\right)-G(t;\vartheta )\right|}^{2}w\left(t\right)dt>0$ for all $\vartheta \in \Theta$.

(ii) 

The parameter space Θ is compact and $\vartheta \mapsto G(t;\vartheta )$ is continuous for each $t\in [0,1]$.

Then,

$$T_n\stackrel{\mathrm{a}.\mathrm{s}.}{\to }\infty ,n\to \infty ,$$

and the proposed test is consistent against fixed alternatives.

Proof. 

By Proposition 5, ${\widehat{G}}_n\left(t\right)\to G_{\mathrm{true}}\left(t\right)$ uniformly almost surely. Since $\Theta$ is compact and $\vartheta \mapsto {\int }_0^1{|G_{\mathrm{true}}\left(t\right)-G(t;\vartheta )|}^{2}w\left(t\right)dt$ is continuous, its infimum ${\delta }^{*}>0$ is attained at some ${\vartheta }^{*}\in \Theta$. Let ${\widehat{\vartheta }}_n$ converge to ${\vartheta }^{*}$ (by compactness, convergence holds along subsequences, and the argument extends to almost sure convergence of the objective). Therefore,

$$n^{-1}{T}_n={\int }_0^1{\left[{\widehat{G}}_n\left(t\right)-G(t;{\widehat{\vartheta }}_n)\right]}^{2}w\left(t\right)dt\stackrel{\mathrm{a}.\mathrm{s}.}{\to }{\delta }^{*}>0,$$

and hence $T_n=n·n^{-1}{T}_n\stackrel{\mathrm{a}.\mathrm{s}.}{\to }\infty$. Similar consistency arguments can be found in the framework of empirical characteristic function based tests; see Baringhaus and Henze [13].    □

Remark 6. 

The preceding propositions characterize the asymptotic behavior of the proposed test statistic under the null hypothesis and under fixed alternatives. The arguments rely on standard tools from empirical process theory, including uniform laws of large numbers, the continuous mapping theorem, and dominated convergence arguments applied to functionals of the empirical probability generating function. Similar quadratic-type goodness-of-fit statistics based on generating functions have been previously studied in related contexts; see, for example, Rueda, O’Reilly and Pérez-Abreu [6] and Henze [8] for early developments, as well as Baringhaus and Henze  [13,31] for closely related approaches based on empirical characteristic and moment generating functions.

Proposition 8 

(Bootstrap Validity). Under $H_0$ and Assumption 1, the parametric bootstrap consistently approximates the null distribution of $T_n$. Specifically, if $T_n^{*}$ denotes the bootstrap counterpart of $T_n$ computed from a resample $X_1^{*},\dots ,X_n^{*}\stackrel{\mathrm{iid}}{\sim }\mathrm{PTPL}\left({\widehat{\vartheta }}_n\right)$, then

$$\underset{x}{sup}\left|P^{*}(T_n^{*}\le x)-P(T_n\le x)\right|\stackrel{p}{⟶}0,n\to \infty .$$

Proof. 

The statistic $T_n$ can be written as $T_n=\phi (\sqrt{n}({\mathbb{P}}_n-P_{{\vartheta }_0}),\sqrt{n}({\widehat{\vartheta }}_n-{\vartheta }_0))$, where $\phi$ is a Hadamard-differentiable functional of the empirical measure and the normalised estimator: it is a composition of the bounded linear evaluation map $\mu \mapsto \int \mu dw$ and a continuous quadratic functional on ${\ell }^{\infty }\left([0,1]\right)$, hence Hadamard-differentiable tangentially to ${\ell }^{\infty }\left([0,1]\right)\times {\mathbb{R}}^3$; see van der Vaart and Wellner [29], Lemma 3.9.27. Theorem 1 establishes the weak limit of this functional.

Under $H_0$, the bootstrap data-generating distribution $P_{{\widehat{\vartheta }}_n}$ satisfies $\parallel {\widehat{\vartheta }}_n-{\vartheta }_0\parallel =O_p\left(n^{-1/2}\right)$ by Proposition 4. Hence the bootstrap empirical process $\sqrt{n}({\mathbb{P}}_n^{*}-P_{{\widehat{\vartheta }}_n})$ has the same weak limit as $\sqrt{n}({\mathbb{P}}_n-P_{{\vartheta }_0})$ by the parametric bootstrap theorem; see van der Vaart [21], Theorem 23.9, and Bickel and Freedman [32].

The Hadamard-differentiability of $\phi$ and the continuous mapping theorem then imply that the bootstrap distribution of $T_n^{*}$ converges in probability to the same limit as $T_n$, yielding the stated result. Analogous arguments appear in Henze [8] and Novoa-Muñoz and Jiménez-Gamero [10].    □

Remark 7. 

For polynomial weight functions of the form $w\left(t\right)=t^a$ with $a>-1$, the test statistic $T_n$ defined in (3) admits an explicit closed-form representation that does not involve numerical integration. In particular, $T_n$ can be expressed as a finite sum depending on the sample observations and on the model probabilities, together with absolutely convergent series involving the PTPL probability mass function evaluated at the estimated parameter vector.

Although this closed-form expression is theoretically appealing and useful for analytical purposes, its direct implementation requires the evaluation of double sums of order $O\left(n^2\right)$ as well as infinite series that must be truncated in practice. As a consequence, the resulting computational cost becomes prohibitive in Monte Carlo and bootstrap studies. For this reason, and in line with previous contributions on PGF-based goodness-of-fit tests (see, e.g., Henze [8] and Rueda, O’Reilly and Pérez-Abreu [6]), the statistic is computed numerically via quadrature of (3), which proves to be substantially more efficient while yielding stable and accurate results.

7. Bootstrap Implementation

As established in Section 6, the asymptotic null distribution of the test statistic $T_n$ depends on unknown model parameters and involves a nonstandard limiting process. As a consequence, analytical critical values are not readily available. To overcome this difficulty, we adopt a parametric bootstrap procedure, which has proven effective for goodness-of-fit tests based on generating functions; see, for example, Henze [8] and Bickel and Freedman  [32].

The bootstrap procedure is implemented under the null hypothesis $H_0$, using the fitted Poisson–Three-Parameter Lindley model. Let ${\widehat{\vartheta }}_n={({\widehat{\theta }}_n,{\widehat{\alpha }}_n,{\widehat{\beta }}_n)}^{\top }$ denote a consistent estimator of the parameter vector, obtained for instance by maximum likelihood. The bootstrap algorithm proceeds as follows.

7.1. Maximum Likelihood Estimation for the PTPL Model

Let $X_1,\dots ,X_n$ be a random sample from a Poisson–Three-Parameter Lindley distribution with parameter vector $\vartheta ={(\theta ,\alpha ,\beta )}^{\top }$. The likelihood function is given by

$$L_n\left(\vartheta \right)={\displaystyle \prod _{i=1}^n}p(X_i;\vartheta ),$$

where $p(·;\vartheta )$ is defined in (1). The corresponding log-likelihood function is

$${\ell }_n\left(\vartheta \right)={\displaystyle \sum _{i=1}^n}logp(X_i;\vartheta ).$$

The maximum likelihood estimator ${\widehat{\vartheta }}_n={({\widehat{\theta }}_n,{\widehat{\alpha }}_n,{\widehat{\beta }}_n)}^{\top }$ is defined as any solution of

$${\widehat{\vartheta }}_n=arg\underset{\vartheta \in \Theta }{max}{\ell }_n\left(\vartheta \right),$$

provided that the maximizer exists and is unique.

Due to the nonlinear structure of the likelihood equations, closed-form expressions are not available, and numerical optimization methods are required. In practice, the estimator is obtained by maximizing ${\ell }_n\left(\vartheta \right)$ using constrained numerical optimization routines.

Remark 8. 

Two complete algorithmic implementations of the parametric bootstrap procedure described in this section are provided in Appendix A. Algorithm 1 details the Monte Carlo evaluation of the empiricalsize of the proposed GOF-PGF test: data are generated under the null PTPL model, and the outputs are the empirical rejection rates ${\mathrm{err}}_{90}$ and ${\mathrm{err}}_{95}$ at the $10\%$ and $5\%$ nominal levels, respectively. Algorithm 2 details the corresponding evaluation of the empirical power : data are generated from an alternative distribution $\mathcal{D}\left(\psi \right)$ under $H_1$, and the output is the empirical rejection rate ${\widehat{\beta }}_{{\alpha }_0}$ at a prescribed nominal level ${\alpha }_0$. Both algorithms share the same inner bootstrap loop; see Remark A1 in Appendix A.

8. Monte Carlo Study

All simulations were carried out using the R programming language [33]. Given the large number of bootstrap replications required ($M\times B$ evaluations of $T_n$ per experiment), computation time was reduced substantially by exploiting parallel computing. Specifically, the outer Monte Carlo loop was distributed across available processor cores using the parRapply parallelisation routine, which applies a function in parallel over the rows of a matrix via R’s built-in parallel package. This strategy, described in Novoa-Muñoz [34], proved particularly effective in the bootstrap-intensive inner loop of Algorithms 1 and 2, reducing wall-clock time by a factor proportional to the number of cores employed.

8.1. Empirical Size

Table 2 reports the empirical rejection rates ${\widehat{\alpha }}_{0.05}$ and ${\widehat{\alpha }}_{0.10}$ under $H_0$, based on $M=1000$ Monte Carlo replications and $B=500$ bootstrap replications, for $n\in \{50,100,150,200,500\}$, $a\in \{0,1,2\}$, and the six parameter configurations listed in Table 1.

Table 1 summarizes the six PTPL parameter configurations considered in the simulation study, together with the corresponding theoretical mean, variance, and dispersion index.

The results indicate that the bootstrap approximation is accurate even for the smaller sample sizes. In all cases, the empirical size is close to the nominal level, confirming the validity of the parametric bootstrap procedure established in Proposition 8. The configuration ${\mathbf{T}}_6$ ($\theta =3.0$, $\alpha =0.8$, $\beta =0.5$) exhibits slightly elevated rejection rates for $n=50$ (up to $0.060$ at the $5\%$ nominal level), which is consistent with the known sensitivity of bootstrap-based tests to parameter estimation uncertainty in small samples; this inflation diminishes monotonically as n increases, confirming asymptotic size control.

8.2. Empirical Power

To assess the discriminatory ability of the proposed GOF-PGF test, we evaluate its empirical power against four families of count distributions that are commonly used as alternatives in the overdispersed count data literature. All experiments use the same settings as the size study: weight exponent $a=1$, sample size $n=250$, $M=1,000$ Monte Carlo replications, $B=500$ bootstrap replications per sample, and nominal level ${\alpha }_0=0.05$. The empirical power at level ${\alpha }_0$ is estimated as

$${\widehat{\beta }}_{0.05}={\displaystyle \frac{1}{M}}{\displaystyle \sum _{m=1}^M}\mathbf{1}\left\{{\overline{VP}}^{\left(m\right)}\le 0.05\right\},$$

(4)

where ${\overline{VP}}^{\left(m\right)}$ is the bootstrap p-value obtained from the m-th Monte Carlo replicate; see Algorithm 2 in Section 7.

The four alternative distributions considered are:

• D1 — Poisson($\lambda$), a strictly equidispersed model with dispersion index $\mathrm{DI}=1$ for all $\lambda >0$, representing the sharpest possible structural contrast with the PTPL family, which satisfies $\mathrm{DI}>1$ throughout its parameter space.
• D2 — Negative Binomial($r,p$), an overdispersed model with $\mathrm{DI}=1/p$. We fix $r=0.5$ and vary $p\in \{0.10,0.15,0.20,0.25,0.30,0.35\}$, yielding $\mathrm{DI}\in \{10.00,6.67,5.00,4.00,3.33,2.86\}$.
• D3 — COM-Poisson($\lambda ,\nu$), a highly flexible count model that reduces to the Poisson when $\nu =1$ and is strongly overdispersed when $\nu \ll 1$ [20]. We fix $\nu =0.1$ and vary $\lambda \in \{1.47,1.48,1.49,1.50,1.51,1.52\}$. At $\nu =0.1$, the resulting means range from approximately $51.7$ to $70.4$ with $\mathrm{DI}\in [9.07,9.33]$, far beyond the PTPL null range of $\mathrm{DI}\in [1.37,3.18]$.
• D4 — Zero-Inflated Poisson ZIP($\lambda ,\pi$), in which $\mathrm{DI}=1+\lambda \pi \ge 1$. Six configurations are evaluated, covering $\mathrm{DI}\in [1.25,2.40]$ and zero-inflation proportions $\pi \in [0.10,0.70]$.

The results are reported in Table 3.

D1 — Poisson alternatives.

All six Poisson configurations yield high empirical power, ranging from $0.967$ at $\lambda =9.75$ to $0.999$ at $\lambda =10.0$, with a strictly monotone increasing pattern. Since the Poisson distribution is always equidispersed ($\mathrm{DI}=1$) while the PTPL null has $\mathrm{DI}>1$, the two families are structurally distinct regardless of $\lambda$. The high power at $n=250$ confirms that the PGF-based test effectively exploits this difference in dispersion structure, consistent with the theoretical consistency result of Proposition 6.

D2 — Negative Binomial alternatives.

The power decreases as p increases from $0.10$ to $0.35$, corresponding to a decrease in ${\mathrm{DI}}_{\mathrm{NB}}=1/p$ from $10.00$ to $2.86$. This is statistically natural: as ${\mathrm{DI}}_{\mathrm{NB}}$ approaches the PTPL null range $[1.37,3.18]$, discrimination becomes harder. At $p=0.30$ ($\mathrm{DI}=3.33$) and $p=0.35$ ($\mathrm{DI}=2.86$), power drops to $0.867$ and $0.823$, respectively. A minor non-monotonicity at $p=0.25$ ($0.993$) compared to $p=0.20$ ($0.990$) is attributable to Monte Carlo variability and does not affect the overall decreasing trend. Even the lowest value ($0.823$) exceeds the $0.70$ threshold.

D3 — COM-Poisson alternatives.

At $\nu =0.1$, the COM-Poisson is strongly overdispersed ($\mathrm{DI}\in [9.07,9.33]$), far from the PTPL null. The power exhibits a unimodal pattern as a function of $\lambda$, peaking at $\lambda =1.50$ (${\widehat{\beta }}_{0.05}=1.000$) and declining on both sides: $0.876$ at $\lambda =1.47$ and $0.879$ at $\lambda =1.52$. This non-monotone pattern reflects the joint effect of changing mean and distribution shape: as $\lambda$ increases, the empirical PGF shifts in a non-linear manner relative to the PTPL null, producing a unimodal power curve whose peak corresponds to the configuration where the PGF discrepancy $T_n$ is maximised. Although the dispersion indices vary only slightly, the associated means range from $51.7$ to $70.4$, producing differences in the shape of the empirical PGF and hence in $T_n$. For $\lambda \in \{1.49,1.50,1.51\}$, power exceeds $0.99$.

D4 — ZIP alternatives.

Power ranges from $0.901$ to $1.000$. The configuration $\mathrm{ZIP}(2.0,0.70)$, with the highest zero-inflation proportion among the evaluated cases, achieves perfect power. For fixed $\lambda =2.5$, power increases monotonically with $\pi$: $0.901$ at $\pi =0.10$, $0.934$ at $\pi =0.25$, and $0.979$ at $\pi =0.40$, reflecting that larger zero-inflation creates mass concentration at zero that is difficult to reconcile with any PTPL fit.

Overall assessment.

Across all 24 configurations, empirical power ranges from $0.823$ to $1.000$ with mean $0.960$. Twenty out of 24 configurations ($83.3\%$) yield power exceeding $0.90$, and 17 out of 24 ($70.8\%$) exceed $0.95$. Perfect power ($1.000$) is achieved for $\mathrm{NB}(0.5,0.10)$, $\mathrm{CMP}(1.50,0.10)$, and $\mathrm{ZIP}(2.0,0.70)$. The only two configurations below $0.90$ are $\mathrm{NB}(0.5,0.30)$ and $\mathrm{NB}(0.5,0.35)$, both with ${\mathrm{DI}}_{\mathrm{NB}}\in \{3.33,2.86\}$ overlapping with the upper end of the PTPL null dispersion range—precisely the region where any consistent GOF test faces inherently greater difficulty. These results confirm that the proposed test possesses strong discriminatory power against all four alternative families at $n=250$.

9. Real Data Application

To assess the practical utility of the proposed goodness-of-fit test, we apply it to five overdispersed count datasets drawn from industrial safety, biology, reliability, and retail commerce. For each dataset, the PTPL model is fitted by maximum likelihood and the GOF-PGF test is performed with weight exponent $a=1$ and $B=999$ parametric bootstrap replications at the $5\%$ nominal level. For comparison, the Poisson, Negative Binomial, COM-Poisson [20], and Zero-Inflated Poisson (ZIP) models are also fitted; the resulting log-likelihoods, AIC values, and Pearson chi-square goodness-of-fit statistics are reported in Table 4.

The five datasets are described as follows.

• D1 — Industrial accidents [35]. Number of accidents suffered by $n=645$ female munitions workers during a twelve-week period, originally analyzed by Greenwood and Yule [35] and subsequently used as a standard benchmark for overdispersed count models. The sample mean is $\overline{x}=0.4512$, the sample variance is $s^2=0.6300$, and the empirical dispersion index is $\mathrm{DI}=s^2/\overline{x}=1.3964$, confirming mild overdispersion relative to the Poisson model.
• D2 — European corn borer larvae [36]. Number of Pyrausta nubilalis (Hübner) larvae per maize plant recorded by McGuire, Brindley and Bancroft [36] ($n=80$, $\overline{x}=1.8375$, $s^2=4.7707$, $\mathrm{DI}=2.5963$). This dataset exhibits moderate overdispersion and has been widely used to evaluate flexible count distributions; see, e.g., Das et al. [19].
• D3 — European red mites on apple leaves [37]. Number of Panonychus ulmi (Koch) mites per apple leaf observed by Garman [37] ($n=150$, $\overline{x}=1.1467$, $s^2=2.2736$, $\mathrm{DI}=1.9828$).
• D4 — Deaths by horse kicks, Prussian army [38]. The classic dataset due to von Bortkiewicz [38], recording the number of soldiers killed by horse kicks per cavalry corps-year ($n=280$ corps-year units, $\overline{x}=0.7000$, $s^2=0.7627$, $\mathrm{DI}=1.0896$). Although this dataset is only marginally overdispersed, it is included as a near-equidispersed benchmark to examine the calibration of the test under near-Poisson conditions.
• D5 — Products purchased per customer, Black Friday [39]. A contemporary retail dataset freely available on Kaggle [39], comprising $550,069$ transactions recorded at an ABC Pvt Ltd retail store during a Black Friday sales event. The count variable of interest is the number of distinct products purchased per unique customer, derived by grouping the raw transaction file by User_ID and counting the number of distinct Product_ID values per user. This yields $n=5,891$ customer-level observations with $\overline{x}=4.3852$, $s^2=8.3292$, and $\mathrm{DI}=1.8994$, confirming moderate overdispersion. To the best of our knowledge, this is the first application of a PTPL-based model to a Black Friday e-commerce count variable.

The MLE estimates, fit statistics, and test results for all five datasets are summarized in Table 4. The main findings are as follows.

For D1 (industrial accidents), the Poisson model is clearly rejected by the chi-square criterion (${\chi }^2=59.127$, $p<0.001$), reflecting the presence of overdispersion. The Negative Binomial, COM-Poisson, ZIP, and PTPL models all provide substantially better fits, with ZIP achieving the lowest AIC ($1162.34$). The proposed GOF-PGF test does not reject the PTPL model ($T_n^{\mathrm{obs}}=0.000520$, $p_{\mathrm{PGF}}=0.3013$), confirming that the PTPL distribution provides an adequate fit to this dataset.

For D2 (Pyrausta nubilalis larvae), the Poisson model is again strongly rejected (${\chi }^2=65.187$, $p<0.001$). The Negative Binomial achieves the best AIC ($298.10$) and is not rejected by the chi-square test ($p=0.7243$). The PTPL model, with AIC $=300.57$, provides a competitive fit and is not rejected by the proposed test ($T_n^{\mathrm{obs}}=0.015200$, $p_{\mathrm{PGF}}=0.1341$).

For D3 (European red mites), the COM-Poisson achieves the lowest AIC ($448.80$), closely followed by the Negative Binomial ($448.87$) and the PTPL ($450.77$). The Poisson model is strongly rejected. ZIP performs noticeably worse for this dataset ($\mathrm{AIC}=456.87$, $p_{{\chi }^2}=0.0155$). The proposed test does not reject the PTPL fit ($T_n^{\mathrm{obs}}=0.000560$, $p_{\mathrm{PGF}}=0.4785$), further supporting its adequacy as a flexible count model.

For D4 (horse-kick deaths), the dataset is nearly equidispersed ($\mathrm{DI}=1.09$) and all five models provide adequate fits. The Poisson model achieves the best AIC ($630.31$), and none of the competing models offer a meaningful improvement. The GOF-PGF test does not reject PTPL ($T_n^{\mathrm{obs}}=0.000040$, $p_{\mathrm{PGF}}=0.6176$). This result is consistent with the theoretical finding that the test maintains correct size under the null hypothesis even in near-equidispersed scenarios.

For D5 (Black Friday retail counts), the PTPL model achieves the best AIC among all five competing models ($27,391.94$), outperforming the Negative Binomial ($27,841.93$), COM-Poisson ($28,009.34$), and Poisson ($29,237.00$). The ZIP model degenerates to the Poisson in this case ($\widehat{\pi }=0$), as the dataset contains no structural zeros. The estimated parameters are $\widehat{\theta }=0.5218$, $\widehat{\alpha }=-1.1592$, and $\widehat{\beta }=1.8248$, all satisfying the admissibility condition $\alpha \ge -\beta /(1+\theta )$ (here, $-\beta /(1+\theta )=-1.1991$, so $\widehat{\alpha }=-1.1592>-1.1991$). Although the large sample size ($n=5,891$) leads to significant chi-square values for all models, the PTPL chi-square residuals are the smallest across all alternatives, and a visual inspection of the observed versus expected frequencies confirms a markedly better fit at all count values $x\ge 3$. The GOF-PGF test yields $T_n^{\mathrm{obs}}=0.770582$ with a bootstrap p-value of $0.0202$ ($B=999$), suggesting that the PTPL model, while the best available among the candidates, is statistically rejected at the $5\%$ level for this very large dataset, which is expected given the high sensitivity of any consistent test with $n\approx 6,000$ observations. This result illustrates that, in large samples, the GOF-PGF test retains discriminatory power, and that the choice of model should be informed by both the p-value and the AIC criterion.

Taken together, the results across five diverse datasets confirm that the PTPL distribution provides a flexible and competitive fit for overdispersed count data arising in different applied domains, and that the proposed GOF-PGF test correctly identifies adequate fits (D1–D4) while retaining power to detect departures from the null model in large samples (D5).

Remark 9.(Large-sample interpretation of p-values). For very large samples ($n\approx 6,000$), any consistent goodness-of-fit test will eventually reject every parametric model, since no finite-parameter family fits arbitrarily many observations exactly. The rejection of the PTPL model for D5 ($p_{\mathrm{PGF}}=0.0202$) should therefore be interpreted in conjunction with the AIC ranking, which confirms the PTPL as the best-fitting model among all candidates. In practice, the choice of model in large samples should be guided primarily by information criteria and subject-matter considerations, with the p-value serving as a diagnostic tool rather than a definitive criterion [40].

10. Discussion

The proposed PGF-based goodness-of-fit test for the PTPL distribution fills a methodological gap identified in the current literature on flexible count models. The test has three notable strengths: (i) it exploits the closed-form PGF of the PTPL distribution, yielding a computationally efficient statistic evaluated by numerical quadrature; (ii) it is theoretically grounded, with consistency against fixed alternatives (Proposition 6 ) and a rigorous weak limit theorem (Theorem 1); and (iii) its null distribution is consistently approximated by the parametric bootstrap (Proposition 8), avoiding the need for analytically intractable critical values.

The simulation study confirms reliable size control across a wide range of sample sizes and parameter configurations. The choice of weight exponent $a=1$ is recommended as a practical default, as it consistently balances type I error control and power.

Future research directions include:

1.

Local power analysis. Deriving the asymptotic power of $T_n$ against Pitman sequences $G_{1,n}\left(t\right)=G(t;{\vartheta }_0)+n^{-1/2}h\left(t\right)$ would characterise the rate at which the test detects local alternatives and enable comparison with minimax-optimal procedures.

2.

Multivariate and zero-inflated extensions. The PGF approach extends naturally to multivariate PTPL models and to zero-inflated variants, following the framework of Novoa-Muñoz and Jiménez-Gamero [9].

3.

Model selection. Combining the proposed test with information criteria (AIC, BIC) and minimum-distance PGF estimators [41] may yield more robust model selection procedures for overdispersed count data.

4.

Truncated and censored data. Extensions to zero-truncated or right-censored count data, common in reliability and actuarial applications, constitute a natural avenue for future work.

11. Conclusions

This paper has proposed and studied a new goodness-of-fit test for the Poisson–Three-Parameter Lindley distribution based on a Cramér-von Mises type distance between the empirical and theoretical probability generating functions. The test statistic has a closed-form representation, is consistent against fixed alternatives, and its null distribution can be consistently approximated by the parametric bootstrap. Monte Carlo simulations confirm accurate size control and competitive power against several relevant alternatives. The real data application demonstrates the practical utility of the proposed procedure.