<figure><img src="image_1.png" /></figure>

HAL
open science

# On evaluating the efficiency of the delta-lognormal mean estimator and predictor

Aubry Philippe

<figure><img src="image_2.png" /></figure>

## To cite this version:

Aubry Philippe. On evaluating the efficiency of the delta-lognormal mean estimator and predictor. MethodsX, 2025, 9, 10.1016/j.mex.2022.101830. hal-05169181

HAL Id: hal-05169181

https://hal.science/hal-05169181v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Specifications table

<table><tr><td>Subject area:</td><td>Probability distributions</td></tr><tr><td>More specific subject area:</td><td>Relative efficiency of statistical estimators or predictors</td></tr><tr><td>Method name:</td><td>UMVU-based estimation or prediction of the mean for delta-lognormal data</td></tr><tr><td>Name and reference of original method:</td><td><ul><li>J. Aitchison, J. Brown, The lognormal distribution, Cambridge University Press, Cambridge, UK, 1957.</li><li>K. Shimizu, Point estimation, in: E. L. Crow, K. Shimizu (Eds.), Lognormal distributions: theory and applications, Marcel Dekker, New York, USA, 1988, pp. 27-86.</li><li>S. Smith, Evaluating the efficiency of the Δ-distribution mean estimator, Biometrics 44 (1988) 485-493.</li></ul></td></tr><tr><td>Resource availability:</td><td>not applicable</td></tr></table>

Method details

Let $\mathcal{U}$ be a finite population of sampling units, unambiguously identifiable by integer labels $i = 1, 2, \dots, N$. Let $y$ be a variable of interest measured or observed on the sampling units, and the total $t_{\mathcal{U}} = N\bar{y}_{\mathcal{U}}$, with $\bar{y}_{\mathcal{U}}$ the mean of $y$ defined over $\mathcal{U}$. A sample $s \subseteq \mathcal{U}$ of size $n$ is drawn from $\mathcal{U}$ by an ignorable selection mechanism. Classical (frequentist) statistics assume that the $y_i$ ($i \in \mathcal{U}$) are random variables of joint distribution $\xi$. Said another way, $\mathcal{U}$ is itself a random sample drawn from an infinite set of populations sharing the same general statistical properties (i.e., a superpopulation), described by stochastic model $\xi$. In this article, we refer to the situation where, from the sample $s$ at hand, the purpose is either to estimate the expectation of $y$ in the model $\xi$, or to predict the finite population mean $\bar{y}_{\mathcal{U}}$ (or equivalently, the total $t_{\mathcal{U}}$).

We consider here a variable of interest $y$ taking nonnegative values ($y \ge 0$). The sample $s$ can be partitioned into $s = s_0 \cup s_1$, $s_0 \cap s_1 = \emptyset$, with $s_0$ of size $n_0$ having zero values ($y = 0$) and $s_1$ of size $n_1$ having positive values ($y > 0$). If $n_1 = 1$, we note the unique positive value $y_{s_1}$. A characteristic of such data is that they may exhibit a high proportion of zero values. To take this into account in a sufficiently flexible manner, one approach is to use a two-component mixture model. There are two possibilities: (i) increasing the probability of zero values from a distribution defined for $y \ge 0$ (zero-inflated distributions); (ii) introducing a dichotomy between $y = 0$ and $y > 0$ in a mixture model with two separately estimable parts (hurdle-at-zero, conditional, two-part, and delta distributions designate the same thing). If the second possibility is adopted, then a suitable model for nonnegative values is written as:

$$G(y; p_0, \boldsymbol{\theta}) = \begin{cases} p_0 & y = 0 \\ p_0 + (1 - p_0)F(y; \boldsymbol{\theta}) & y > 0 \end{cases} \quad (1)$$

where $F(y; \boldsymbol{\theta})$ is a cumulative distribution with parameters $\boldsymbol{\theta}$, corresponding to a positive distribution, either discrete or continuous; here we consider the lognormal distribution.

Lognormal distribution

Let $z$ be a random variable distributed according to the standard normal distribution (i.e., $z \sim \text{Norm}(0, 1)$) of probability density:

$$\phi(z) = \frac{1}{\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} z^2 \right\} \quad z \in \mathbb{R} \qquad (2)$$

Then, $y = \exp(\mu + \sigma z)$ follows a lognormal distribution that is completely specified by $\mu$ and $\sigma^2$. The probability density function of the lognormal distribution is written as [1, p. 8, Eq. (2.5)]:

$$f(y; \mu, \sigma^2) = \frac{1}{y\sigma} \phi((\ln y - \mu)/\sigma) = \frac{1}{y\sigma\sqrt{2\pi}} \exp \left\{ -\frac{1}{2} \left( \frac{\ln y - \mu}{\sigma} \right)^2 \right\} \quad y \in \mathbb{R}_+^* \qquad (3)$$
## Delta-lognormal distribution

The lognormal distribution is no longer appropriate when zero values must be accounted for. This leads to using the delta-lognormal distribution [2], often also called the Δ-distribution [1] and occasionally the Bernoulli-lognormal two-part model (e.g., [3, p. 703]).

The delta-lognormal distribution results from a mixture of a Dirac mass at 0 with probability $p_0$ and a lognormal distribution with probability $(1 - p_0)$, that is:

$$g(y; p_0, \mu, \sigma^2) = p_0 \delta(y) + (1 - p_0) f(y; \mu, \sigma^2) \quad (4)$$

where $\delta(y)$ is a Dirac distribution that concentrates a unit mass at 0.

The first three cumulants (i.e., expectation, variance, and third central moment) of the distribution (4) are written as [1, p. 95, Eq. (9.43)-(9.45)]:

$$\kappa_1 = (1 - p_0) \alpha \quad (5)$$

$$\kappa_2 = (1 - p_0) \alpha^2 \{\beta - (1 - p_0)\} \quad (6)$$

$$\kappa_3 = (1 - p_0) \alpha^3 \{\beta^3 - 3(1 - p_0)\beta + 2(1 - p_0)^2\} \quad (7)$$

with:

$$\alpha = \exp \left( \mu + \frac{1}{2} \sigma^2 \right) \quad (8)$$

$$\beta = \exp(\sigma^2) \quad (9)$$

Even if it is not necessarily the best possible definition, in this article, skewness is classically defined as:

$$\gamma_1 = E \left\{ \left( \frac{y - \kappa_1}{\kappa_2^{1/2}} \right)^3 \right\} = \frac{\kappa_3}{\kappa_2^{3/2}} \quad (10)$$

The skewness $\gamma_1$ of the delta-lognormal distribution increases dramatically as $\sigma^2$ increases. For small values of $\sigma^2$, as $p_0$ increases, $\gamma_1$ first decreases and then increases. For $\sigma^2 > 1$ approximately, as $p_0$ increases, $\gamma_1$ only increases (Fig. 1).

## Relative efficiency assessement

In the context of a finite population, one can indifferently consider the mean $\bar{y}_{it}$ or the total $t_{it}$ as the quantity of interest. For the finite population that has actually been sampled, these quantities have fixed values that one would be able to know exactly if $s = U$ (ignoring possible measure or observation errors). Under a superpopulation model $\xi$, these statistics are random variables whose values one wants to predict. In an infinite population (superpopulation), one is interested in estimating the expectation $E(y) = \kappa_1$.

For delta-lognormal data, depending on the level of skewness of the distribution, the question arises of the gain in precision that can be achieved by relying on the uniformly minimum-variance unbiased estimator (UMVUE) compared to using the sample mean $\bar{y}_s$, either for estimating $\kappa_1$ or for predicting $\bar{y}_{it}$ (or equivalently, $t_{it}$). In other words, one may compare the situation where the shape of the distribution is known, to the situation where it is unknown (or known but not taken into account), first in the case of an infinite population (estimation context), then in that of a finite population (prediction context).

In the estimation context, Aitchison and Brown [1, p. 98, Fig. 9.1] provided relative efficiency results only for $p_0 = 0.5$ and for the degenerate case of the lognormal distribution ($p_0 = 0$), using a variance approximation (see the validation section). By doing so, the sample size $n$ is disregarded in
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Skewness $\gamma_1$ for the delta-lognormal distribution as a function of $p_0$ and $\sigma^2$. For maximum bounds of the parameter space considered on the figure ($p_0 = 0.9, \sigma^2 = 2$), the skewness is $\gamma_1 \approx 64.5$.</figcaption></figure>

the relative efficiency assessment. Shimizu [2] did not document the relative efficiency in the case of the delta-lognormal distribution. Smith [4] considered the relative efficiency for $p_0 = 0.1$ and $p = 0.5$, for very small sample sizes, using exact or approximate variance. To our knowledge, the relative efficiency assessment in the prediction context has not been documented yet.

In this technical article, after providing a compendium of fundamental formulas for UMVU estimation in the case of the delta-lognormal distribution, we document the relative efficiency more thoroughly than in the past by considering both the estimation and prediction contexts, taking into account the sample size (and the finite population size in the prediction context), and varying the probability of getting a zero value up to $p_0 = 0.9$. In all cases we use the exact expression of the variance of the estimator (or predictor).

## Estimation context

### Unknown shape of the distribution

Let $s$ be a sample of size $n$ drawn by random sampling from an infinite population of unknown shape. The unbiased estimator of $\kappa_1$ is the sample mean $\bar{y}_s$ and its sampling variance is written as:

$$V(\bar{y}_s) = \frac{\kappa_2}{n} \qquad (11)$$

where $\kappa_2$ is estimated without bias by:

$$S_s^2 = \frac{1}{n-1} \sum_{i \in s} (y_i - \bar{y}_s)^2 \qquad (12)$$

The sampling variance is then estimated without bias by:

$$\hat{V}(\bar{y}_s) = \frac{S_s^2}{n} \qquad (13)$$
Known shape of the distribution

When y is distributed according to a delta-lognormal distribution, $κ_1$ can be estimated by the UMVUE [1, p. 97, Eq. (9.54)] (typo corrected); [4]:

$$ \hat{\kappa}_1 = \begin{cases} (1 - \hat{p}_0) \exp(\bar{y}_{\ln}) g_{n_1-1} \left( \frac{1}{2} s_{\ln}^2 \right) & n_1 > 1 \\ \frac{y_{s_1}}{n} & n_1 = 1 \\ 0 & n_1 = 0 \end{cases} \qquad (14) $$

with:

$$ \bar{y}_{\ln} = \frac{1}{n_1} \sum_{i \in s_1} \ln y_i \qquad (15) $$

$$ s_{\ln}^2 = \frac{1}{n_1 - 1} \sum_{i \in s_1} (\ln y_i - \bar{y}_{\ln})^2 \qquad (16) $$

and $g_m(t)$ an infinite series introduced by Finney [5, Eq. (10)], which can be written as [6, Eq. (1.2)], [7, Eq. (3)]:

$$ g_m(t) = \sum_{j=0}^{\infty} \frac{1}{j!} \frac{m+2j}{m} \left( \frac{m}{m+1} t \right)^j \prod_{i=1}^{j} \frac{m}{m+2i} \qquad (17) $$

The exact variance of $\hat{\kappa}_1$ (14) was provided by Smith [4, Eq. (6)]. The function $g_m(t)$ belongs to the class of generalized hypergeometric functions and can be written as a particular instance of the confluent hypergeometric limit function (here denoted as ${}_0F_1$) as [6, Eq. (2.1)]:

$$ g_m(t) = {}_0F_1 \left( \frac{m}{2}; \frac{m^2 t}{2(m+1)} \right) \qquad (18) $$

with [8, Eq. (4)], [9, p. 333, Eq. (16.3.1)]:

$$ {}_0F_1(a; z) = \sum_{j=0}^{\infty} \frac{z^j}{j! (a)_j} \qquad (19) $$

where $(a)_j$ is the notation used in special function theory for the rising factorial:

$$ (a)_j = \begin{cases} 1 & j=0 \\ \prod_{i=0}^{j-1} (a+i) & j \ge 1 \end{cases} \qquad (20) $$

Note that the numerical evaluation of special functions ${}_0F_1(a; z)$ and $g_m(t)$ is addressed later in the article.

Using the confluent hypergeometric limit function, according to Shimizu [2, p. 50], estimators $\hat{\kappa}_1$ and $\hat{\kappa}_2$ can be written, respectively, as:

$$ \hat{\kappa}_1 = \begin{cases} (1 - \hat{p}_0) \exp(\bar{y}_{\ln}) {}_0F_1 \left( \frac{n_1-1}{2}; \frac{(n_1-1)^2}{4n_1} s_{\ln}^2 \right) & n_1 > 1 \\ \frac{y_{s_1}}{n} & n_1 = 1 \\ 0 & n_1 = 0 \end{cases} \qquad (21) $$
and

$$
\hat{\kappa}_2 = \begin{cases} (1 - \hat{p}_0) \exp(2\bar{y}_{\ln}) \left[ 0F_1\left(\frac{n_1-1}{2}; \frac{(n_1-1)^2}{n_1}s_{\ln}^2\right) \right. \\ \left. - \left(\frac{n_1-1}{n-1}\right)0F_1\left(\frac{n_1-1}{2}; \frac{(n_1-1)(n_1-2)}{2n_1}s_{\ln}^2\right) \right] & n_1 > 1 \\ \frac{y_{s_1}^2}{n} & n_1 = 1 \\ 0 & n_1 = 0 \end{cases} \quad (22)
$$

The exact variance of $\hat{\kappa}_1$ is [10, Remark 3.1], [2, pp. 50-51]:

$$
V(\hat{\kappa}_1) = \alpha^2 \left[ \frac{1}{n^2} \sum_{j=2}^{n} \binom{n}{j} (1-p_0)^j p_0^{n-j} j^2 \exp\left(\frac{\sigma^2}{j}\right) \times_0 F_1\left(\frac{j-1}{2}; \frac{(j-1)^2}{4j^2}\sigma^4\right) - (1-p_0)^2 \right] (23)
$$

Equivalent approximations of $V(\hat{\kappa}_1)$ were given by Shimizu [2, p. 51] and Aitchison and Brown [1, p. 99, Eq. (9.58)] (see the validation section); in this article we use the exact expression (23), which translates into Algorithm 1, using Algorithm 2 [11].

**Algorithm 1** Exact variance $V(\hat{\kappa}_1)$ (23) with $n$ the sample size, $p, m, v$ the values for parameters $p_0, \mu, \sigma^2$, respectively.

1: **function ExactVariance(n, p, m, v)**

2: $q \leftarrow 1 - p$

3: $x \leftarrow q$

4: $y \leftarrow p^{(n-1)}$

5: $S \leftarrow 0$

6: **for** $j = 2$ **to** $n$ **do**

7: $x \leftarrow x \times q$

8: $y \leftarrow y/p$

9: $z \leftarrow ((j-1)/(2 \times j) \times v)^2$

10: $S \leftarrow S + \textbf{Binomial}(n, j) \times x \times y \times j^2 \times \exp(v/j) \times \textbf{Hypergeometric0F1}((j-1)/2, z)$

11: **end for** $j$

12: $r \leftarrow \exp(2 \times m + v) \times (S/n^2 - q^2)$

13: [ return $r$ ]

14: **end function**

Pennington [12] has provided the UMVUE for the variance $V(\hat{\kappa}_1)$ expressed with the function $g_m(t)$ [12, Eq. (4)]; for the sake of homogeneity, we present it here expressed using the confluent hypergeometric limit function:

$$
\hat{V}(\hat{\kappa}_1) = \begin{cases} (1 - \hat{p}_0) \exp(2\bar{y}_{\ln}) \left[ (1 - \hat{p}_0) {}_0F_1\left(\frac{n_1-1}{2}; \frac{(n_1-1)^2}{4n_1}s_{\ln}^2\right) \right. \\ \left. - \left(\frac{n_1-1}{n-1}\right) {}_0F_1\left(\frac{n_1-1}{2}; \frac{(n_1-1)(n_1-2)}{2n_1}s_{\ln}^2\right) \right] & n_1 > 1 \\ \left(\frac{y_{s_1}}{n}\right)^2 & n_1 = 1 \\ 0 & n_1 = 0 \end{cases} \quad (24)
$$

Relative efficiency

We compare the precision of the estimator for $\kappa_1$ according to whether the distribution shape is unknown or known. It is expected that taking into account the knowledge of the distribution shape will lead to a gain in precision. The relative efficiency is defined in the same way as that used by
<caption>**Algorithm 2** Iterative computation of the binomial coefficient.</caption>

1: function Binomial(n, k)
2: p ← n − k
3: if k < p then
4:     p ← k
5:     k ← n − p
6: end if
7: if p < 0 then
8:     [ return 0 ]
9: end if
10: if p = 0 then
11:     [ return 1 ]
12: end if
13: r ← k + 1
14: for i = 2 to p do
15:     r ← r × (k + i) / i
16: end for i
17: [ return r ]
18: end function

▷ error in the arguments

▷ trivial case

Aitchison and Brown [1, p. 99, Eq. (9.62)]:

$$ \text{eff}_1 = \frac{V(\hat{\kappa}_1)}{V(\hat{\gamma}_s)} \qquad (25) $$

that is, taking into account the shape of the distribution leads to a gain in precision, which is higher when eff₁ is low ($V(\hat{\kappa}_1) < V(\hat{\gamma}_s)$). Hence, we quantify the gain in precision by expressing it as a function of $p_0$, $\sigma^2$ (or, equivalently, $\sigma$) and $n$. The parameter $\mu$ vanishes through the elimination of the $\alpha^2$ term that appears in the numerator and denominator of the relative efficiency.

We vary $p_0 = 0(0.025)0.9$, $\sigma = 0.05(0.05)2$ and $n = 50, 100, 500, 1000$. The results show that the gain in precision increases substantially (eff₁ decreases) as $\sigma$ increases and increases more weakly – for a fixed ($p_0, \sigma$) point – as $n$ increases (Fig. 2). For a fixed $\sigma$ value, the gain decreases (eff₁ increases) as $p_0$ increases, and this is all the more pronounced when $n$ is small (see Fig. 2a). This phenomenon diminishes as $n$ increases (see Fig. 2b–d). Finally, we note a small difference between $n = 500$ (Fig. 2c) and $n = 1000$ (Fig. 2d), which suggests a convergence of eff₁ toward its asymptotic limit.

To illustrate the speed of convergence of eff₁ toward its asymptotic limit, we set $\sigma^2 = 2$, and we repeat the calculations for $p_0 = 0.0(0.1)0.90$ and $50 \le n \le 1000$. The asymptotic limit and the speed of convergence of eff₁ depend on $p_0$; the lower $p_0$ is, the higher the efficiency gain and the speed of convergence (Fig. 3).

## Prediction context

By denoting $r = U - s$, the total defined on the population can be written as:

$$ t_{U} = \sum_{i \in U} y_{i} = \sum_{i \in S} y_{i} + \sum_{i \in R} y_{i} \qquad (26) $$

that is the sum of the totals defined over the sample ($t_s$) and the remaining part in the population ($t_r$). A predictor of $t_{U}$ can be written as:

$$ \tilde{t}_{U} = \sum_{i \in S} y_i + \sum_{i \in R} \tilde{y}_i \qquad (27) $$
Consider a simple mean model where the random variables $y_i$ ($i = 1, \dots, N$) have the same expectations and variances and are not correlated, that is:

$$E(y_i) = \kappa_1 \qquad (28a)$$

$$V(y_i) = \kappa_2 \qquad (28b)$$

$$\mathrm{Cov}(y_i, y_j) = 0 \quad (i \neq j) \qquad (28c)$$

The empirical predictor is written as:

$$\tilde{t}_{il} = \sum_{i \in S} y_i + \underbrace{(N - n)\hat{E}(y_i)}_{\tilde{t}_r} \qquad (29)$$

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Relative efficiency $\tilde{t}_1$ as a function of $p_0$ and $\sigma$ for different sample sizes ($n$). (a) $n = 50$. (b) $n = 100$. (c) $n = 500$. (d) $n = 1000$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Relative efficiency $\tilde{e}_1$ as a function of sample size ($n$) for $\sigma^2 = 2$ and $p_0 = 0.0(0.1)0.9$.</figcaption></figure>

Predictor (29) is model-unbiased; that is, if the model is correct, then $E(\tilde{t}_{il} - t_{il}) = 0$. The variance of the prediction error is obtained as:

$$V(\tilde{t}_{il} - t_{il}) = V(\tilde{t}_r - t_r) = V(\tilde{t}_r) + V(t_r) - 2 \underbrace{\text{Cov}(\tilde{t}_r, t_r)}_{0} \quad (30)$$

The $\xi$-covariance between $\tilde{t}_r$ and $t_r$ is zero since: (i) $\tilde{t}_r$ is a function of the set of values in $s$ ($\{y_i, i \in s\}$), not of the set of values in $r$ ($\{y_i, i \in r\}$) which we do not know; and (ii) the two sets of values $\{y_i, i \in s\}$ and $\{y_i, i \in r\}$ are uncorrelated under the model (see (28c)).

### Distribution of unknown shape

Predictor (29) can be written as:

$$\tilde{t}_{il}^{\text{exp}} = \sum_{i \in s} y_i + (N - n)\bar{y}_s = n\bar{y}_s + (N - n)\bar{y}_s = N\bar{y}_s \quad (31)$$

which we can designate as an expansion predictor [13]. From relation (30), the prediction error variance of $(\tilde{t}_{il}^{\text{exp}} - t_{il})$ is obtained as:

$$V(\tilde{t}_{il}^{\text{exp}} - t_{il}) = (N - n)^2 [V(\bar{y}_s) + V(\bar{y}_r)] \quad (32)$$

$$= (N - n)^2 \left( \frac{1}{n} + \frac{1}{N - n} \right) \kappa_2 \quad (33)$$

$$= N^2 \left( 1 - \frac{n}{N} \right) \frac{\kappa_2}{n} \quad (34)$$
The predictor of the mean is $\bar{y}_s$, and its prediction error variance is $V(\bar{y}_s - \bar{y}_{t_l}) = (1 - n/N)\kappa_2/n$. It follows the limit:

$$ \lim_{N \to \infty} V(\bar{y}_s - \bar{y}_{t_l}) = V(\bar{y}_s) \qquad (35) $$

For prediction error variance (34), an unbiased estimator is obtained by substituting estimator $\hat{\kappa}_2$ (22) for parameter $\kappa_2$.

### Distribution of known shape

Predictor (29) can be written as:

$$ \tilde{\hat{\kappa}}_{t_l}^{\text{mvu}} = \sum_{i \in S} y_i + (N - n)\hat{\kappa}_1 \qquad (36) $$

From relation (30), the prediction error variance $(\tilde{\hat{\kappa}}_{t_l}^{\text{mvu}} - t_{t_l})$ is obtained as:

$$ V(\tilde{\hat{\kappa}}_{t_l}^{\text{mvu}} - t_{t_l}) = (N - n)^2 [V(\hat{\kappa}_1) + V(\bar{y}_r)] \qquad (37) $$

$$ = (N - n)[(N - n)V(\hat{\kappa}_1) + \kappa_2] \qquad (38) $$

where the variance $V(\hat{\kappa}_1)$ is given by expression (23). The mean predictor is $\tilde{\bar{y}}_{t_l} = N^{-1}\tilde{\hat{\kappa}}_{t_l}^{\text{mvu}}$ (e.g., [14, Eq. (3)]), and its prediction error variance is $V(\tilde{\bar{y}}_{t_l} - \bar{y}_{t_l}) = N^{-2}V(\tilde{\hat{\kappa}}_{t_l}^{\text{mvu}} - t_{t_l})$. It follows the limits:

$$ \lim_{N \to \infty} \tilde{\bar{y}}_{t_l} = \lim_{N \to \infty} \left[ \frac{1}{N} \sum_{i \in S} y_i + \frac{(N - n)}{N} \hat{\kappa}_1 \right] = \hat{\kappa}_1 \qquad (39) $$

$$ \lim_{N \to \infty} V(\tilde{\bar{y}}_{t_l} - \bar{y}_{t_l}) = \lim_{N \to \infty} \left[ \frac{(N - n)^2}{N^2} V(\hat{\kappa}_1) + \frac{(N - n)}{N^2} \kappa_2 \right] = V(\hat{\kappa}_1) \qquad (40) $$

For prediction error variance (38), an unbiased estimator is obtained by substituting estimators $\hat{V}(\hat{\kappa}_1)$ (24) and $\hat{\kappa}_2$ (22) for $V(\hat{\kappa}_1)$ and $\kappa_2$, respectively.

### Relative efficiency

We compare the precision of the $t_{t_l}$ prediction depending on whether the shape of the distribution is unknown or known. In addition to $p_0$, $\sigma^2$ and $n$, we must also vary the finite population size $N$. The relative efficiency is defined as:

$$ \text{eff}_2 = \frac{V(\tilde{\hat{\kappa}}_{t_l}^{\text{mvu}} - t_{t_l})}{V(\tilde{\hat{\kappa}}_{t_l}^{\text{exp}} - t_{t_l})} = \frac{V(\tilde{\bar{y}}_{t_l} - \bar{y}_{t_l})}{V(\bar{y}_s - \bar{y}_{t_l})} \qquad (41) $$

From (35) and (40), we obtain:

$$ \lim_{N \to \infty} \text{eff}_2 = \text{eff}_1 \qquad (42) $$

First, we examine the effect of the population size $N$ for a fixed sample size $n$, which is equivalent to examining the effect of the sampling fraction $f = n/N$. We set $n = 50$ and vary the sampling fraction as $f = 0.4, 0.2, 0.1, 0.05$ ($N = 125, 250, 500, 1000$). As before, we vary $p_0 = 0(0.025)0.9$ and $\sigma = 0.05(0.05)2$. The obtained results (Fig. 4) show a smaller gain in precision than in the case of an infinite population (Fig. 2a). This finite population effect is less pronounced as $f$ tends toward 0 ($N \to \infty$) since in that scenario the situation tends toward the asymptotic result of the case considered here, which corresponds to Fig. 2a.

To illustrate the speed of convergence of $\text{eff}_2$ toward its asymptotic limit $\text{eff}_1$, we set $\sigma^2 = 2$, and we repeat the calculations for $p_0 = 0.0(0.1)0.90, n = 50$ and $500 \le N \le 10000$. As in the infinite population case, the asymptotic limit and the speed of convergence of $\text{eff}_2$ depend on $p_0$. The lower $p_0$ is, the higher the efficiency gain ($\text{eff}_2$ decreases). In contrast, the speed of convergence becomes less important when $p_0$ decreases (Fig. 5). In practice, we can consider that we are almost at convergence
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Relative efficiency $\text{eff}_2$ as a function of $p_0$ and $\sigma$ for a fixed sample size ($n = 50$) and different finite population sizes ($N$) and hence different sampling fractions ($f$). (a) $N = 125$, $f = 0.4$. (b) $N = 250$, $f = 0.2$. (c) $N = 500$, $f = 0.1$. (d) $N = 1000$, $f = 0.05$.</figcaption></figure>

as soon as $N = 5000$ (or $N = 10000$ if one is more conservative). Note that limit values correspond to the values for $n = 50$ in Fig. 3 and are represented by dotted half lines in Fig. 5.

Finally, we examine the gain in precision when $n$ and $N$ increase jointly, keeping the sampling fraction constant. For $f = 0.1$, we increase the population size as $N = 500, 1000, 5000, 10000$ ($n = 50, 100, 500, 1000$). The obtained results (Fig. 6) show a gain in precision that increases ($\text{eff}_2$ decreases) as $N$ and $n$ jointly increase. Note that Fig. 6a is the same as Fig. 4c ($N = 500, n = 50$). There is only a small difference between the case $N = 5000, n = 500$ (Fig. 6c) and the case $N = 10000, n = 1000$ (Fig. 6d), which suggests the convergence of $\text{eff}_2$ toward its asymptotic limit (in the sense that $n$ and $N$ jointly increase and for $f = 0.1$).

## Numerical evaluation of the generalized hypergeometric functions

In this article, we need to numerically evaluate ${}_0F_1(a; z)$ ($a \in \mathbb{R}_+^*$ and $z \in \mathbb{R}_+^*$) — or equivalently $g_m(t)$ — to calculate $V(\hat{k}_1)$ (23), which is directly involved in the relative efficiency $\text{eff}_1$ (25) and through formula (38) in the relative efficiency $\text{eff}_2$ (41).
<figure><img src="image_8.png" /><figcaption>Fig. 5. Relative efficiency $\text{eff}_2$ as a function of population size (N) for a fixed sample size (n = 50), for $\sigma^2 = 2$ and $p_0 = 0.0(0.1)0.9$. Details in the text.</figcaption></figure>

### Evaluating $\text{O}_1 F_1(a; z)$ by recurrence

The confluent hypergeometric limit function (19) can be written as:

$$ \text{O}_1 F_1(a; z) = 1 + \frac{z}{a} + \sum_{j=2}^{\infty} \frac{1}{\prod_{i=0}^{j-1} (a+i)} \times \frac{z^j}{j!} \quad (43) $$

which gives the recurrence relation:

$$ \begin{aligned} A_1 &= \frac{z}{a} & S_1 &= 1 + A_1 \\ A_j &= A_{j-1} \times \frac{1}{a+j-1} \times \frac{z}{j} & S_j &= S_{j-1} + A_j \quad \text{for } j = 2, 3, \dots, k \end{aligned} \quad (44) $$

and translates into Algorithm 3.

The series is infinite, but in practice, it is only required to compute the sum until it reaches a level of convergence considered sufficient; for example, with the convergence criterion $|S_j - S_{j-1}| < \epsilon$, we used $\epsilon = 10^{-10}$. Another way to proceed is to determine the number of terms necessary to reach the precision allowed by the computer at hand (see [15, pp. 88-89]). Computing successive terms of the series by means of the recurrence relation is a computational method often used by default (e.g., [16, p. 99]). For other methods, the interested reader is referred to [17,18].

### Evaluating $g_m(t)$ by recurrence

One can proceed in the same way as previously for the function $g_m(t)$ (17), which can be written as (see, for example, [4, Eq. (1)], with $m = n_1 - 1$):

$$ g_m(t) = 1 + \frac{m}{m+1} \times t + \sum_{j=2}^{\infty} \frac{m^{2j-1}}{(m+1)^j \prod_{i=2}^{j-1} (m+2(i-1))} \times \frac{t^j}{j!} \quad (45) $$
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Relative efficiency $\text{eff}_2$ as a function of $p_0$ and $\sigma$ for a fixed sampling fraction $f = 0.1$ and increasing finite population ($N$) and sample ($n$) sizes. (a) $N = 500$ ($n = 50$). (b) $N = 1000$ ($n = 100$). (c) $N = 5000$ ($n = 500$). (d) $N = 10000$ ($n = 1000$).</figcaption></figure>

**Algorithm 3** Confluent hypergeometric limit function ${}_0F_1$.

1: **function Hypergeometric0F1(a, z, ϵ)**

2: A ← z/a

3: R ← 1 + A

4: j = 2

5: **loop**

6: A ← A × 1/(a + j − 1) × z/j

7: S ← R + A

8: **if** |S − R| < ϵ **then**

9: [ return S ]

10: **end if**

11: j ← j + 1

12: R ← S

13: **end loop**

14: **end function**
which gives the recurrence relation:

$$ 
\begin{aligned}
A_1 &= \frac{m}{m+1} \times t \\
A_j &= A_{j-1} \times \frac{m^2}{(m+1)(m+2(j-1))} \times \frac{t}{j}
\end{aligned}
\qquad S_1 = 1 + A_1 \\
\qquad S_j = S_{j-1} + A_j \quad \text{for } j = 2, 3, \dots, k \qquad (46)
$$

and translates into Algorithm 4.

<caption>**Algorithm 4** Finney's generalized hypergeometric function $g_m(t)$</caption>

1: **function Finney(m, t, ϵ)**
2: &nbsp;&nbsp;&nbsp;&nbsp;A ← m/(m + 1) × t
3: &nbsp;&nbsp;&nbsp;&nbsp;B ← A × m
4: &nbsp;&nbsp;&nbsp;&nbsp;R ← 1 + A
5: &nbsp;&nbsp;&nbsp;&nbsp;j ← 2
6: &nbsp;&nbsp;&nbsp;&nbsp;**loop**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;A ← A × B/(m + 2 × (j − 1))/j
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;S ← R + A
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** |S − R| < ϵ **then**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;[ return S ]
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
12: &nbsp;&nbsp;&nbsp;&nbsp;j ← j + 1
13: &nbsp;&nbsp;&nbsp;&nbsp;R ← S
14: &nbsp;&nbsp;&nbsp;&nbsp;**end loop**
15: **end function**

## Validation

In this section: (i) we check the calculation accuracy of the confluent hypergeometric limit function ${}_0F_1(a; z)$ required for computing $V(\hat{\kappa}_1)$ using Algorithm 1; (ii) we assess the correctness of Algorithm 1 by comparing its results against those of Monte Carlo simulations; (iii) we justify the use of the exact expression for the variance $V(\hat{\kappa}_1)$ (23) for accurately documenting the relative efficiency when considering $\hat{\kappa}_1$ in the estimation or prediction contexts.

*Precision of ${}_0F_1(a; z)$ evaluated by recurrence*

We verify that the calculations performed in this article — using double-precision floating-point arithmetic — by implementing the computation of ${}_0F_1(a; z)$ through Algorithm 3 with $\epsilon = 10^{-10}$ have sufficient numerical accuracy. To do this, we first examine the ranges of variation of parameter $a$ and argument $z$ and then compare the results of our procedure with reference values.

The examination of $V(\hat{\kappa}_1)$ (23) shows that the minimal value that $a$ can take in this expression is $a = 0.5$ and that for the maximal value $n = 1000$ used in this article, we have $a \simeq 500$. Besides, we may consider $0 < z \le 1000$ for covering the situations addressed in [19] and similar future examples.

For $a = 0.5$, it is possible to use a particular representation of ${}_0F_1(a; z)$ [20, p. 594, Table 7.13.1 (row 5, column 2)]:

$${}_0F_1(0.5; z) = \cosh(2\sqrt{z}) \qquad (47)$$

with $\cosh = (\exp(z) + \exp(-z))/2$. For $z = 0.01(0.01)1$, we obtain an absolute error less than $10^{-12}$ (at most $k = 8$ iterations are required). For $z = 1(1)1000$, we obtain a relative error less than $10^{-12}$. To gain an idea of the order of magnitude, ${}_0F_1(0.5; 1000) \simeq 1.466103955241 \times 10^{27}$ (at most $k = 69$ iterations are required).

For $a > 0.5$, we use as reference values those returned by the function Hypergeometric0F1[a,z] of Mathematica®. For $a = 5(5)500$ and $z = 0.01(0.01)1$, we obtain an absolute error less than $10^{-12}$ (at most $k = 7$ iterations are required). For $a = 5(5)500$ and $z = 10(10)1000$, we obtain a relative error less than $10^{-11}$ (at most $k = 67$ iterations are required).
<caption>**Table 1**</caption>

Results of the Monte Carlo simulations to assess the correctness of Algorithm 1. Details in the text.

<table><thead><tr><th>p<sub>0</sub></th><th>n</th><th>K</th><th>&kappa;<sub>1</sub></th><th>E<sub>MC</sub>(&kappa;<sub>1</sub>)</th><th>V(&kappa;<sub>1</sub>)</th><th>V<sub>MC</sub>(&kappa;<sub>1</sub>)</th></tr></thead><tbody><tr><td>0.1</td><td>50</td><td>10<sup>8</sup></td><td>29.8039</td><td>29.8039</td><td>268.1132</td><td>268.2039</td></tr><tr><td>0.1</td><td>100</td><td>10<sup>8</sup></td><td>29.8039</td><td>29.8027</td><td>126.5944</td><td>126.5823</td></tr><tr><td>0.1</td><td>500</td><td>10<sup>7</sup></td><td>29.8039</td><td>29.8032</td><td>24.1667</td><td>24.1690</td></tr><tr><td>0.1</td><td>1000</td><td>10<sup>7</sup></td><td>29.8039</td><td>29.8028</td><td>12.0127</td><td>12.0134</td></tr><tr><td>0.5</td><td>50</td><td>10<sup>8</sup></td><td>16.5577</td><td>16.5584</td><td>167.2448</td><td>167.2557</td></tr><tr><td>0.5</td><td>100</td><td>10<sup>8</sup></td><td>16.5577</td><td>16.5568</td><td>75.8297</td><td>75.8418</td></tr><tr><td>0.5</td><td>500</td><td>10<sup>7</sup></td><td>16.5577</td><td>16.5574</td><td>13.9909</td><td>13.9922</td></tr><tr><td>0.5</td><td>1000</td><td>10<sup>7</sup></td><td>16.5577</td><td>16.5565</td><td>6.9244</td><td>6.9199</td></tr><tr><td>0.9</td><td>50</td><td>10<sup>8</sup></td><td>3.3115</td><td>3.3113</td><td>61.8555</td><td>61.3674</td></tr><tr><td>0.9</td><td>100</td><td>10<sup>8</sup></td><td>3.3115</td><td>3.3115</td><td>22.3205</td><td>22.3108</td></tr><tr><td>0.9</td><td>500</td><td>10<sup>7</sup></td><td>3.3115</td><td>3.3108</td><td>3.1210</td><td>3.1220</td></tr><tr><td>0.9</td><td>1000</td><td>10<sup>7</sup></td><td>3.3115</td><td>3.3112</td><td>1.4862</td><td>1.4867</td></tr></tbody></table>

## Monte Carlo simulations

We simulate *K* random samples of size *n* from a delta-lognormal distribution of parameters *p*₀, *μ* = 1.5 and *σ* = 2, and we compute the UMVUE (21) for each sample. This results in *K* values *k*₁ from which we calculate Monte Carlo approximations of expectation (denoted as *E*<sub>MC</sub>(*k*₁)) and variance (denoted as *V*<sub>MC</sub>(*k*₁)). Then, we can compare these Monte Carlo approximations to their theoretical values, that is, *κ*₁ — since *k*₁ is unbiased — and *V*(*k*₁), respectively. We vary *n* = 50, 100, 500, 1000 and *p*₀ = 0.1, 0.5, 0.9. For accurate Monte Carlo approximations, we set *K* = 10⁸ for *n* = 50, 100 and *K* = 10⁷ for *n* = 500, 1000. Table 1 shows a very good agreement between the Monte Carlo approximations and the calculated theoretical values. We conclude that, for 50 ≤ *n* ≤ 1000 and *p*₀ ≤ 0.9, one can trust Algorithm 1.

## Variance approximation of V(k̂₁)

An approximation of $V(k̂₁)$ in $O(n^{-2})$ is [2, p. 51]:

$$V_{\text{app}}(k̂_1) = \frac{1}{n}(1 - p_0)\left[p_0 + \frac{1}{2}\sigma^2(\sigma^2 + 2)\right]\exp(2\mu + \sigma^2) \quad (48)$$

Approximation (48) is algebraically equivalent to the approximation given by Aitchison and Brown [1, p. 99, Eq. (9.58)]. Aitchison & Brown indicated that this approximation is given for large *n* and *p*₀ substantially less than 1, without further details. The quality of approximation (48) depends on the sample size *n* and parameters *p*₀ and *σ²*. To specify the domains of use of the approximation, we compute the ratio $V_{\text{app}}(k̂_1)/V(k̂_1)$ for *p* = 0(0.025)0.9, *σ* = 0.05(0.05)2.0 and *n* = 50, 100, 500, 1000 (Fig. 7). The approximation improves when *σ* and *p*₀ are small, and for the same fixed (*p*₀, *σ*) point, the accuracy increases as the sample size *n* increases. Thus, accurate calculations for great *p*₀ (i.e., here up to *p*₀ = 0.9) and *σ²* values requires using formula (23) (translated into Algorithm 1) instead of approximation (48), unless *n* is sufficiently large.

## Additional information

Pennington [12] introduced the use of the delta-lognormal distribution in marine biology to estimate mean abundance more efficiently than can be achieved with the sample mean in the case of highly skewed distributions. The article by Pennington [12] is widely cited in the literature (603 citations according to Google Scholar, at the time of writing this article). Regarding the robustness of this approach, the reader is referred to Myers and Pepin [21,22], Syrjala [23] and Christman [24]. In the context of the mean (or total) prediction, see the recent contribution [19].
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Ratio $V_{\text{app}}(\hat{k}_1)/V(\hat{k}_1)$ as a function of $p_0$ and $\sigma$ for different sample sizes ($n$). (a) $n = 50$. (b) $n = 100$. (c) $n = 500$. (d) $n = 1000$.</figcaption></figure>

**Declaration of Competing Interest**

The author declare that he has no known competing for financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgments**

I thank American Journal Expert (AJE) and my colleague Matthieu Guillemain for checking the English of this article.

**References**

[1] J. Aitchison, J. Brown, The Lognormal Distribution (Reprinted from 1957 Edition), Cambridge University Press, Cambridge, UK, 1969.

[2] K. Shimizu, Point estimation, in: E.L. Crow, K. Shimizu (Eds.), Lognormal distributions: theory and applications, Marcel Dekker, New York, USA, 1988, pp. 27–86.
[3] B. Neelon, A. O'Malley, Two-part models for zero-modified count and semicontinuous data, in: A. Levy, S. Goring, C. Gatsonis, B. Sobolev, E. van Ginneken, R. Busse (Eds.), Health services evaluation, Springer, New York, USA, 2019, pp. 695–716.

[4] S. Smith, Evaluating the efficiency of the Δ-distribution mean estimator, Biometrics 44 (1988) 485–493.

[5] D. Finney, On the distribution of a variate whose logarithm is normally distributed, Supp. J. Roy. Stat. Soc. 7 (1941) 155–161.

[6] D. Ebbeler, A note on estimation in log-normal linear models, J. Stat. Comput. Simul. 2 (1973) 225–231.

[7] X.-H. Zhou, S. Gao, Estimation of the log-normal mean, Stat. Med. 17 (1998) 2251–2264.

[8] G. Taraldsen, A precise estimator for the log-normal mean, Stat. Methodol. 2 (2005) 111–120.

[9] A. Cuyt, V. Petersen, B. Verdonk, H. Waadeland, W. Jones, Handbook of Continued Fractions for Specials Functions, Springer, Dordrecht, The Netherlands, 2008.

[10] K. Shimizu, K. Iwase, Uniformly minimum variance unbiased estimation in lognormal and related distributions, Commun. Stat. - Theory Methods 10 (1981) 1127–1147.

[11] M. Wolfson, H. Wright, Algorithm 160, Combinatorial of M things taken N at a time, Commun. ACM 6 (1963) 161.

[12] M. Pennington, Efficient estimators of abundance, for fish and plankton surveys, Biometrics 39 (1983) 281–286.

[13] H. Bolfarine, S. Zacks, Prediction Theory for Finite Populations, Springer, New York, USA, 1992.

[14] S. Smith, Use of statistical models for the estimation of abundance from groundfish trawl survey data, Can. J. Fish. Aquat. Sci. 47 (1990) 894–903.

[15] R. Forrey, Computing the hypergeometric function, J. Comput. Phys. 137 (1997) 79–100.

[16] J. Pearson, Computation of hypergeometric functions, Master's thesis, University of Oxford, UK, 2009.

[17] F. Johansson, Computing hypergeometric functions rigorously, ACM Trans. Math. Softw. 45 (2019) 1–26.

[18] G. Navas-Palencia, High-precision computation of uniform asymptotic expansions for special functions, Polytechnic University of Catalonia, Barcelona, Spain, 2019 Ph.D. thesis.

[19] P. Aubry, C. Francesiaz, On comparing design-based estimation versus model-based prediction to assess the abundance of biological populations, Ecol. Indic. (2022). (submitted)

[20] A. Prudnikov, Y. Brychkov, O. Marichev, Integrals and Series. Volume 3. More Special Functions, Gordon and Breach Science Publishers, New York, USA, 1990.

[21] R. Myers, P. Pepin, The robustness of lognormal-based estimators of abundance, Biometrics 46 (1990) 1185–1192.

[22] R. Myers, P. Pepin, Rejoinder to the letter to the editors from M. Pennington, "On testing the robustness of lognormal-based estimators", Biometrics 47 (1991) 1623–1624.

[23] S. Syrjala, Critique on the use of the delta distribution for the analysis of trawl survey data, ICES J. Mar. Sci. 57 (2000) 831–842.

[24] M. Christman, Review of Estimation Methods for Parameters of the Delta-Lognormal Distribution, Technical Report, MCC Statistical Consulting LLC, Gainesville, Gainesville, FL, USA, 2019.