Generalized Moment Estimators Based on Stein Identities

Abstract

For parameter estimation of continuous and discrete distributions, we propose a generalization of the method of moments (MM), where Stein identities are utilized for improved estimation performance. The construction of these Stein-type MM-estimators makes use of a weight function as implied by an appropriate form of the Stein identity. Our general approach as well as potential benefits thereof are first illustrated by the simple example of the exponential distribution. Afterward, we investigate the more sophisticated two-parameter inverse Gaussian distribution and the two-parameter negative-binomial distribution in great detail, together with illustrative real-world data examples. Given an appropriate choice of the respective weight functions, their Stein-MM estimators, which are defined by simple closed-form formulas and allow for closed-form asymptotic computations, exhibit a better performance regarding bias and mean squared error than competing estimators.

Introduction

For many parametric distributions, so-called Stein identities are available, which rely on moments of functional expressions of a corresponding random variable. These identities are named after Charles Stein, who developed the idea of uniquely characterizing a certain distribution family by such a moment identity [see 21, 22]. Many examples of both continuous and discrete distributions together with their Stein characterizations can be found in Stein et al. [23], Sudheesh [24], Sudheesh and Tibiletti [25], Landsman and Valdez [15], Weiß and Aleksandrov [30], Anastasiou et al. [1] and the references therein. Stein identities are not a mere tool of probability theory. During the last years, there was also a lot of research activity on statistical applications of Stein identities, for example to goodness-of-fit (GoF) tests [4, 31], sequential change-point monitoring by control charts [29], and shrinkage estimation [14]. For references on these and further statistical applications related to Stein’s method, see Betsch et al. [4], Anastasiou et al. [1], Kubokawa [14]. In the present article, a further statistical application of Stein identities is investigated and exemplified, namely to the development of closed-form parameter estimators for continuous or discrete distributions. The idea to construct generalized types of method-of-moments (MM) estimators based on an appropriate type of Stein identity plus weighting function, referred to as Stein-MM estimators, was first explored in some applications by Arnold et al. [2] and Wang and Weiß [28]. Recently, Ebner et al. [6] discussed the Stein-MM approach in a much broader way, and also the present article provides a comprehensive treatment of Stein-MM estimators for various distributions. The main motivation for considering Stein-MM estimation is that the weighting function might be chosen in such a way that the resulting estimator shows better properties (e. g., a reduced bias or mean squared error (MSE)) than the default MM estimator or other existing estimators. Despite the additional flexibility offered by the weighting function, the Stein-MM estimators are computed from simple closed-form expressions, and consistency and asymptotic normality are easily established, also see Ebner et al. [6].

In what follows, we apply the proposed Stein-MM estimation to three different distribution families. We start with the illustrative example of the exponential (Exp) distribution in Sect. 2. This simple one-parameter distribution mainly serves to demonstrate the general approach for deriving the Stein-MM estimator and its asymptotics, and it also indicates the potential benefits of using the Stein-MM approach for parameter estimation. Afterward in Sect. 3, we examine a more sophisticated type of continuous distribution, namely the two-parameter inverse Gaussian (IG) distribution. In Sect. 4, we then turn to a discrete distribution family, namely the two-parameter negative-binomial (NB) distribution. Illustrative real-world data examples are also presented in Sects. 3–4. Note that neither the exponential distribution nor any discrete distribution have been considered by Ebner et al. [6], and their approach to the Stein-MM estimation of the IG-distribution differs from the one proposed here, see the details below. Also Arnold et al. [2], Wang and Weiß [28] did not discuss any of the aforementioned distributions. Finally, we conclude in Sect. 5 and outline topics for future research.

Stein Estimation of Exponential Distribution

The exponential distribution is the most well-known lifetime distribution, which is characterized by the property of being memory-less. It has positive support and depends on the parameter $\lambda >0$, where its probability density function (pdf) is given by $\varphi (x)=\lambda e^{-\lambda x}$ for $x>0$ and zero otherwise. A detailed survey about the properties of and estimators for the $Exp(\lambda )$-distribution can be found in Johnson et al. [10, Chapter 19]. Given the independent and identically distributed (i. i. d.) sample $X_1,\dots ,X_n$ with $X_i\sim Exp(\lambda )$ for $i=1,\dots ,n$, the default estimator of $\lambda >0$, which is an MM estimator and the maximum likelihood (ML) estimator at the same time, is given by $\widehat{\lambda }=1/\overline{X}$, where $\overline{X}=\frac{1}{n}{\sum }_{i=1}^n{X}_i$ denotes the sample mean. This estimator is known to neither be unbiased, nor to be optimal in terms of the MSE, see Elfessi & Reineke [7]. To derive a generalized MM estimator with perhaps improved bias or MSE properties, we consider the exponential Stein identity according to Stein et al. [23, Example 1.6], which states that

2.1

$$\begin{array}{c}\hfill X\sim Exp(\lambda )\text{iff}\mathbb{E}[f^{\prime }(X)]=\lambda \mathbb{E}[f(X)]\end{array}$$for any piecewise differentiable function f with $f(0)=0$ such that $\mathbb{E}[|f^{\prime }(X)|]$, $\mathbb{E}[|f(X)|]$ exist. Solving (2.1) in $\lambda$ and using the sample moments $\overline{f^{\prime }(X)}=\frac{1}{n}{\sum }_{i=1}^n{f}^{\prime }(X_i)$ and $\overline{f(X)}=\frac{1}{n}{\sum }_{i=1}^{n}f(X_i)$ instead of the population moments, the class of Stein-MM estimators for $\lambda$ is obtained as

2.2

$$\begin{array}{c}\hfill {\widehat{\lambda }}_{f,Exp}=\frac{\overline{f^{\prime }(X)}}{\overline{f(X)}}.\end{array}$$Note that the choice $f(x)=x$ leads to the default estimator $\widehat{\lambda }=1/\overline{X}$. Generally, $f(x)\ne x$ might be interpreted as a kind of weighting function, which assigns different weights to large or low values of X than the identity function does. For deriving the asymptotic distribution of the general Stein-MM estimator ${\widehat{\lambda }}_{f,Exp}$, we first define the vectors ${\mathit{Z}}_i$ with $i=1,\dots ,n$ as

2.3

$$\begin{array}{c}\hfill {\mathit{Z}}_i:=(f^{\prime }(X_i),f(X_i){)}^{\top }.\end{array}$$Their mean equals

2.4

$$\begin{array}{c}\hfill {\mathit{\mu }}_Z:=\mathbb{E}[{\mathit{Z}}_i]=\left(\begin{array}{c}\mathbb{E}[f^{\prime }(X_i)]\\ \mathbb{E}[f(X_i)]\end{array}\right)\stackrel{(2.1)}{=}{\mu }_f(0,1,0)\left(\begin{array}{c}\lambda \\ 1\end{array}\right),\end{array}$$where we define ${\mu }_f(k,l,m):=\mathbb{E}[X^k·f{(X)}^l·f^{\prime }{(X)}^m]$ for any $k,l,m\in {\mathbb{N}}_0=\{0,1,\dots \}$. Then, the following central limit theorem (CLT) holds.

Theorem 2.1

If $X_1,\dots ,X_n$ are i. i. d. according to $Exp(\lambda )$, then the sample mean $\overline{\mathit{Z}}$ of ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_n$ according to (2.3) is asymptotically normally distributed as$$\begin{array}{c}\hfill \sqrt{n}(\overline{\mathit{Z}}-{\mathit{\mu }}_Z)\stackrel{\text{d}}{\to }N(\mathbf{0},\mathbf{\Sigma })\text{with}\mathbf{\Sigma }={({\sigma }_{\mathit{ij}})}_{i,j=1,2},\end{array}$$where $N(\mathbf{0},\mathbf{\Sigma })$ denotes the multivariate normal distribution, and where the covariances are given as$$\begin{array}{cc}\hfill {\sigma }_{11}& ={\mu }_f(0,0,2)-{\lambda }^2·{\mu }_f^2(0,1,0),\hfill \\ \hfill {\sigma }_{22}& ={\mu }_f(0,2,0)-{\mu }_f^2(0,1,0),\hfill \\ \hfill {\sigma }_{12}& =\frac{\lambda }{2}({\sigma }_{22}-{\mu }_f^2(0,1,0))=\lambda ({\sigma }_{22}-\frac{1}{2}{\mu }_f(0,2,0)).\hfill \end{array}$$

The proof of Theorem 2.1 is provided by Appendix A.1. In the second step of deriving the asymptotics of ${\widehat{\lambda }}_{f,Exp}$, we define the function $g(u,v):=\frac{u}{v}$. Then, ${\widehat{\lambda }}_{f,Exp}=g(\overline{\mathit{Z}})$ and $\lambda =g({\mathit{\mu }}_Z)$. Applying the Delta method [18] to Theorem 2.1, the following result follows.

Theorem 2.2

If $X_1,\dots ,X_n$ are i. i. d. according to $Exp(\lambda )$, then ${\widehat{\lambda }}_{f,Exp}$ is asymptotically normally distributed, where the asymptotic variance and bias are given by$$\begin{array}{c}\hfill {\sigma }_{f,Exp}^2=\frac{1}{n}·\frac{{\mu }_f(0,0,2)}{{\mu }_f^2(0,1,0)},{\mathbb{B}}_{f,Exp}=\frac{\lambda }{2n}·\frac{{\mu }_f(0,2,0)}{{\mu }_f^2(0,1,0)}.\end{array}$$

The proof of Theorem 2.2 is provided by Appendix A.2. Note that the moments ${\mu }_f(k,l,m)$ involved in Theorems 2.1 and 2.2 can sometimes be derived explicitly, see the subsequent examples, while they can be computed by using numerical integration otherwise.

After having derived the asymptotic variance and bias without explicitly specifying the function f, let us now consider some special cases of this weighting function. Here, our general strategy is as follows. We first specify a parametric class of functions, where the actual parameter value(s) are determined only after a second step. In this second step, we compute asymptotic properties of the Stein-MM estimator such as bias, variance, and MSE, and then we select the parameter value(s) such that some of the aforementioned properties are minimized within the considered class of functions. Detailed examples are presented below. The chosen parametric class of functions should be sufficiently flexible in the sense that by modifying its parameter value(s), it should be possible to move the weight to quite different regions of the real numbers. Its choice may also be guided by the aim of covering existing parameter estimators as special cases. For the exponential distribution discussed here in Sect. 2, we already pointed out that the default ML and MM estimator $\widehat{\lambda }=1/\overline{X}$ corresponds to $f(x)=x$, so the choice of a parametric class of functions including $f(x)=x$ appears reasonable. This leads to our first illustrative example, namely the choice $f_a:[0;\infty )\to [0;\infty )$, $f_a(x)=x^a$, where $a=1$ leads to the default estimator $\widehat{\lambda }=1/\overline{X}$. Here, we have to restrict to $a>0$ to ensure that the condition $f(0)=0$ in (2.1) holds. Using that

2.5

$$\begin{array}{c}\hfill \mathbb{E}[X^a]=\frac{\mathrm{\Gamma }(a+1)}{{\lambda }^a}\text{for}X\sim Exp(\lambda )\text{if}a>-1,\end{array}$$the following corollary to Theorem 2.2 is derived.

Corollary 2.3

Let $X_1,\dots ,X_n$ be i. i. d. according to $Exp(\lambda )$, and let $f_a(x)=x^a$ with $a>\frac{1}{2}$. Then, ${\widehat{\lambda }}_{f_a,Exp}$ is asymptotically normally distributed, where the asymptotic variance and bias are given by$$\begin{array}{c}\hfill {\sigma }_{f_a,Exp}^2=\frac{{\lambda }^2}{n}\left(\begin{array}{c}2(a-1)\\ a-1\end{array}\right),{\mathbb{B}}_{f_a,Exp}=\frac{\lambda }{2n}\left(\begin{array}{c}2a\\ a\end{array}\right).\end{array}$$Furthermore, the MSE equals$$\begin{array}{cc}\hfill {\text{MSE}}_{f_a,Exp}& ={\sigma }_{f_a,Exp}^2+{\mathbb{B}}_{f_a,Exp}^2={\lambda }^2[\frac{1}{n}\left(\begin{array}{c}2(a-1)\\ a-1\end{array}\right)+\frac{1}{4n^2}{\left(\begin{array}{c}2a\\ a\end{array}\right)}^2].\hfill \end{array}$$

The proof of (2.5) and Corollary 2.3 is provided by Appendix A.3. Note that in Corollary 2.3, $\left(\begin{array}{c}r\\ s\end{array}\right)$ denotes the generalized binomial coefficient given by $\mathrm{\Gamma }(r+1)/\mathrm{\Gamma }(s+1)/\mathrm{\Gamma }(r-s+1)$.

Fig. 1 [Images not available. See PDF.]

Plot of (a) $n{\sigma }_{f_a,Exp}^2$ and $n{\mathbb{B}}_{f_a,Exp}$, and (b) MSE${}_{f_a,Exp}$ for $a\in (0.5;1.5)$ and $\lambda =1$. Points indicate minimal MSE values. Dotted line at $a=1$ corresponds to default estimator $\widehat{\lambda }=1/\overline{X}$

In Fig. 1a, the asymptotic variance and bias of ${\widehat{\lambda }}_{f_a,Exp}$ according to Corollary 2.3 are presented. While the variance is minimal for $a=1$ (i. e., for the ordinary MM and ML estimator), the bias decreases with decreasing a (i. e., bias reductions are achieved for sublinear choices of $f_a(x)=x^a$). Hence, an MSE-optimal choice of $f_a(x)$ is obtained for some $a\in (0.5;1)$. This is illustrated by Fig. 1b, where the MSE of Corollary 2.3 is presented for different sample sizes $n\in \{10,25,50,100\}$. The corresponding optimal values of a are determined by numerical optimization as 0.952, 0.978, 0.988, and 0.994, respectively. As a result, especially for small n, we have a reduction of the MSE (and of the bias as well) if using a “true” Stein-MM estimator (i. e., with $a\ne 1$). Certainly, if the focus is mainly on bias reduction, then an even smaller choice of a would be beneficial.

As a second illustrative example, let us consider the functions $f_u(x)=1-u^x$ with $u\in (0;1)$, which are again sublinear, but this time also bounded from above by one. Again, we can derive a corollary to Theorem 2.2, this time by using the moment formula

2.6

$$\begin{array}{c}\hfill \mathbb{E}[u^X]=\frac{\lambda }{\lambda -ln(u)}\text{for}X\sim Exp(\lambda )\text{if}u\in (0;1).\end{array}$$

Corollary 2.4

Let $X_1,\dots ,X_n$ be i. i. d. according to $Exp(\lambda )$, and let $f_u(x)=1-u^x$ with $u\in (0;1)$. Then, ${\widehat{\lambda }}_{f_u,Exp}$ is asymptotically normally distributed, where the asymptotic variance and bias are given by$$\begin{array}{c}\hfill {\sigma }_{f_u,Exp}^2=\frac{\lambda }{n}\frac{{(\lambda -ln(u))}^2}{\lambda -2ln(u)},{\mathbb{B}}_{f_u,Exp}=\frac{{\sigma }_{f_u,Exp}^2}{\lambda -ln(u)}.\end{array}$$Furthermore, the MSE equals$$\begin{array}{cc}\hfill {\text{MSE}}_{f_u,Exp}& =(1+\frac{\lambda }{n(\lambda -2ln(u))})·{\sigma }_{f_u,Exp}^2.\hfill \end{array}$$

The proof of (2.6) and Corollary 2.4 is provided by Appendix A.4.

Fig. 2 [Images not available. See PDF.]

Plot of MSE${}_{f_u,Exp}$ for $u\in (0;1)$, for (a) different n and $\lambda =1$, and (b) different $\lambda$ and $n=25$. Points indicate minimal MSE values

This time, the variance decreases for increasing u, whereas the bias decreases for decreasing u. As a consequence, an MSE-optimal choice is expected for some u inside the interval (0; 1). This is illustrated by Fig. 2a, where the minima for $n\in \{10,25,50,100\}$, given that $\lambda =1$, are attained for $u\approx 0.918$, 0.963, 0.981, and 0.990, respectively. The major difference between the two types of weighting functions in Corollaries 2.3 and 2.4 is given by the role of $\lambda$ within the expression for the MSE. For $f_a(x)$ in Corollary 2.3, $\lambda$ occurs as a simple factor such that the optimal choice for a is the same across different $\lambda$. Hence, the optimal a is simply a function of the sample size n, which is very attractive for applications in practice. For $f_u(x)$ in Corollary 2.4, by contrast, the MSE depends in a more sophisticated way on $\lambda$, and the optimal u differs for different $\lambda$ as illustrated by Fig. 2b. Thus, if one wants to use the weighting function $f_u(x)$ in practice, a two-step procedure appears reasonable, where an initial estimate is computed via $\widehat{\lambda }=1/\overline{X}$, which is then refined by ${\widehat{\lambda }}_{f_u,Exp}$ with u being determined by plugging-in $\widehat{\lambda }$ instead of $\lambda$ (also see Section 2.2 in Ebner et al. [6] for an analogous idea).

We conclude this section by pointing out two further application scenarios for the use of Stein-MM estimators ${\widehat{\lambda }}_{f,Exp}$. First, in analogy to recent Stein-based GoF-tests such as in Weiß et al. [31], ${\widehat{\lambda }}_{f,Exp}$ might be used for GoF-applications. More precisely, the idea could be to select a set $\{f_1,\dots ,f_K\}$ of weighting functions, and to compute ${\widehat{\lambda }}_{f_k}$ for all $k=1,\dots ,K$. As any ${\widehat{\lambda }}_{f_k}$ is a consistent estimator of $\lambda$ according to Theorem 2.2, the obtained values $\{{\widehat{\lambda }}_{f_1},\dots ,{\widehat{\lambda }}_{f_K}\}$ should vary closely around $\lambda$. For other continuous distributions with positive support, such as the IG-distribution considered in the next Sect. 3, we cannot expect that $\overline{f^{\prime }(X)}/\overline{f(X)}$ has an (asymptotically) unique mean for different f, see Remark 3.1, so a larger variation among the values in $\{{\widehat{\lambda }}_{f_1},\dots ,{\widehat{\lambda }}_{f_K}\}$ is expected. Such a discrepancy in variation might give rise for a formal exponential GoF-test. But as the focus of this article is on parameter estimation, we postpone a detailed investigation of this GoF-application to future research.

Table 1. Empirical bias and MSE of ${\widehat{\lambda }}_{f,Exp}$ from ${10}^5$ simulated i. i. d. samples $X_1,\dots ,X_n\sim Exp(1)$, where $\lceil 0.1·n\rceil$ observations randomly selected for additive outlier “$+5$”

	Bias				MSE
$n\backslash f$	$x^{0.9}$	$1-0.9^x$	$ln(1+x)$	x	$x^{0.9}$	$1-0.9^x$	$ln(1+x)$	x
10	$-$0.280	$-$0.277	$-$0.206	$-$0.304	0.107	0.105	0.100	0.114
25	$-$0.343	$-$0.341	$-$0.271	$-$0.366	0.126	0.125	0.091	0.140
50	$-$0.305	$-$0.303	$-$0.238	$-$0.328	0.098	0.097	0.067	0.111
100	$-$0.308	$-$0.306	$-$0.241	$-$0.330	0.098	0.096	0.063	0.111

A second type of application is illustrated by Table 1, which refers to a simulation experiment with ${10}^5$ replications per scenario. For simulated i. i. d. $Exp(1)$-samples of sizes $n\in \{10,25,50,100\}$, about 10 % of the observations were randomly selected and contaminated by an additive outlier, namely by adding 5 to the selected observations. Note that the topic of outliers in exponential data received considerable interest in the literature [10, pp. 528–530]. Then, different estimators ${\widehat{\lambda }}_{f,Exp}$ are computed from the contaminated data, where the first three choices of the weighting function f are characterized by a sublinear increase, whereas the fourth function, $f(x)=x$, corresponds to the default estimator $\widehat{\lambda }=1/\overline{X}$. Table 1 shows that all MM estimators are affected by the outliers, e. g., in terms of the strong negative bias. But comparing the four columns of bias and MSE values, respectively, it gets clear that the novel Stein-MM estimators are more robust against the outliers, having both lower bias and MSE than $\widehat{\lambda }$. Especially the choice $f(x)=ln(1+x)$, a logarithmic weighting scheme, leads to a rather robust estimator. The relatively good performance of the Stein-MM estimators can be explained by the fact that the weighting functions increase sublinearly (which is also beneficial for bias reduction in non-contaminated data, recall the above discussion), so the effect of large observations is damped. To sum up, by choosing an appropriate weighting function f within the Stein-MM estimator ${\widehat{\lambda }}_{f,Exp}$, one cannot only achieve a reduced bias and MSE, but also a reduced sensitivity towards outlying observations.

Stein Estimation of Inverse Gaussian Distribution

Like the exponential distribution considered in the previous Sect. 2, the IG-distribution with parameters $\mu ,\lambda >0$, abbreviated as $IG(\mu ,\lambda )$, has positive support, where the pdf is given by$$\begin{array}{c}\hfill \varphi (x)=\sqrt{\frac{\lambda }{2\pi }}{e}^{\lambda /\mu }{x}^{-3/2}exp(-\frac{\lambda }{2\mu }(\frac{x}{\mu }+\frac{\mu }{x}))\text{for}x>0,\text{and 0 otherwise}.\end{array}$$The IG-distribution is commonly used as a lifetime model (as it can be related to the first-passage time in random walks), but it may also simply serve as a distribution with positive skewness and, thus, as an alternative to, e. g., the lognormal distribution [see ]. Detailed surveys about the properties and applications of $IG(\mu ,\lambda )$, and on many further references, can be found in Folks & Chhikara [8], Seshadri [19] as well as in Johnson et al. [10, Chapter 15]. In what follows, the moment properties of $X\sim IG(\mu ,\lambda )$ are particularly relevant. We have $\mathbb{E}[X]=\mu$, $\mathbb{E}[1/X]=1/\mu +1/\lambda$, and $\mathbb{V}[X]={\mu }^3/\lambda$. In particular, positive and negative moments are related to each other by

3.1

$$\begin{array}{c}\hfill \mathbb{E}[X^{-k}]=\frac{\mathbb{E}[X^{k+1}]}{{\mu }^{2k+1}}\iff \mathbb{E}[{(\mu /X)}^k]=\mathbb{E}[{(X/\mu )}^{k+1}]\text{for}k\in {\mathbb{N}}_0,\end{array}$$see Tweedie [26, p. 372] as well as the aforementioned surveys.

Remark 3.1

At this point, let us briefly recall the discussion in Sect. 2 (p. 7), where we highlighted the property that for i. i. d. exponential samples, the quotient $\overline{f^{\prime }(X)}/\overline{f(X)}$ has an (asymptotically) unique mean for different f. From counterexamples, it is easily seen that this property is not true for $IG(\mu ,\lambda )$-data. The Delta method implies that the mean of $\overline{f^{\prime }(X)}/\overline{f(X)}$ is asymptotically equal to $\mathbb{E}[f^{\prime }(X)]/\mathbb{E}[f(X)]$, which equals

• $1/\mathbb{E}[X]=1/\mu$ for $f(x)=x$, but

• $2\mathbb{E}[X]/\mathbb{E}[X^2]=2\mu /({\mu }^2(1+\frac{\mu }{\lambda }))=2\lambda /(\mu (\lambda +\mu ))$ for $f(x)=x^2$.

From now on, let $X_1,\dots ,X_n$ be an i. i. d. sample from $IG(\mu ,\lambda )$, which shall be used for parameter estimation. Here, one obviously estimates $\mu$ by the sample mean $\overline{X}$, but the estimation of $\lambda$ is more demanding. In the literature, the MM and ML estimation of $\lambda$ have been discussed (see the details below), while our aim is to derive a generalized MM estimator with improved bias and MSE properties based on a Stein identity. In fact, as we shall see, our proposed approach can be understood as a unifying framework that covers the ordinary MM and ML estimator as special cases. A Stein identity for $IG(\mu ,\lambda )$ has been derived by Koudou & Ley [13, p. 172], which states that

3.2

$$\begin{array}{c}\hfill X\sim IG(\mu ,\lambda )\text{iff}\mathbb{E}[f(X)(\lambda X^2-{\mu }^{2}X-\lambda {\mu }^2)]=2{\mu }^2\mathbb{E}[X^2{f}^{\prime }(X)]\end{array}$$holds for all differentiable functions $f:(0;\infty )\to \mathbb{R}$ with ${lim}_{x\to 0}f(x)\varphi (x)={lim}_{x\to \infty }f(x)\varphi (x)=0$. Solving (3.2) in $\lambda$ and using the sample moments $\overline{h(X)}=\frac{1}{n}{\sum }_{i=1}^{n}h(X_i)$ instead of $\mathbb{E}[h(X)]$ (where h might be any of the functions involved in (3.2)), the class of Stein-MM estimators for $\lambda$ is obtained as

3.3

$$\begin{array}{c}\hfill {\widehat{\lambda }}_{f,IG}=\frac{{\overline{X}}^2(2\overline{X^2{f}^{\prime }(X)}+\overline{Xf(X)})}{\overline{X^{2}f(X)}-{\overline{X}}^2\overline{f(X)}}.\end{array}$$Here, the ordinary MM estimator of $\lambda >0$, i. e., ${\widehat{\lambda }}_{\mathit{MM}}={\overline{X}}^3/S^2$ with $S^2$ denoting the empirical variance [27], is included as the special case $f\equiv 1$, whereas the ML estimator ${\widehat{\lambda }}_{\mathit{ML}}=\overline{X}/(\overline{X}·\overline{1/X}-1)$ [26] follows for $f(x)=1/x$. Hence, (3.3) provides a unifying estimation approach that covers the established estimators as special cases.

Remark 3.2

At this point, a reference to Example 2.9 in Ebner et al. [6] is necessary. As already mentioned in Sect. 1, also Ebner et al. [6] proposed a Stein-MM estimator for the IG-distribution, which, however, differs from the one developed here. The crucial difference is given by the fact that Ebner et al. [6] tried a joint estimation of $(\mu ,\lambda )$ based on (3.2), namely by jointly solving two equations that are implied by (3.2) if using two different weight functions $f_1\ne f_2$. The resulting class of estimators, however, does not cover the existing MM and ML estimators, so Ebner et al. [6] did not pursue the Stein-MM estimation of the IG-distribution further. By contrast, as we did not see notable potential for improving the estimation of $\mu$ by $\overline{X}$ (recall the diverse optimality properties of the sample mean as an estimator of the population mean [e. g., 20]), we used (3.2) to only derive an estimator for $\lambda$. In this way, we were able to recover both the MM and ML estimator of $\lambda$ within (3.3).

For deriving the asymptotic distribution of our general Stein-MM estimator ${\widehat{\lambda }}_{f,IG}$ from (3.3), we first define the vectors ${\mathit{Z}}_i$ with $i=1,\dots ,n$ as

3.4

$$\begin{array}{c}\hfill {\mathit{Z}}_i:=(X_i,f(X_i),X_{i}f(X_i),X_i^{2}f(X_i),X_i^2{f}^{\prime }(X_i){)}^{\top }.\end{array}$$Their mean equals

3.5

$$\begin{array}{c}\hfill {\mathit{\mu }}_Z:=\mathbb{E}[{\mathit{Z}}_i]=(\mu ,{\mu }_f(0,1,0),{\mu }_f(1,1,0),{\mu }_f(2,1,0),{\mu }_f(2,0,1){)}^{\top }.\end{array}$$Then, the following CLT holds.

Theorem 3.3

If $X_1,\dots ,X_n$ are i. i. d. according to $IG(\mu ,\lambda )$, then the sample mean $\overline{\mathit{Z}}$ of ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_n$ according to (3.4) is asymptotically normally distributed as$$\begin{array}{c}\hfill \sqrt{n}(\overline{\mathit{Z}}-{\mathit{\mu }}_Z)\stackrel{\text{d}}{\to }N(\mathbf{0},\mathbf{\Sigma })\text{with}\mathbf{\Sigma }={({\sigma }_{\mathit{ij}})}_{i,j=1,\dots ,5},\end{array}$$where $N(\mathbf{0},\mathbf{\Sigma })$ denotes the multivariate normal distribution, and where the covariances are given as$$\begin{array}{ccc}& {\sigma }_{11}={\mu }^3/\lambda ,\hfill & \hfill {\sigma }_{23}={\mu }_f(1,2,0)-{\mu }_f(0,1,0)·{\mu }_f(1,1,0),\\ \hfill & {\sigma }_{12}={\mu }_f(1,1,0)-\mu ·{\mu }_f(0,1,0),\hfill & \hfill {\sigma }_{24}={\mu }_f(2,2,0)-{\mu }_f(0,1,0)·{\mu }_f(2,1,0),\\ \hfill & {\sigma }_{13}={\mu }_f(2,1,0)-\mu ·{\mu }_f(1,1,0),\hfill & \hfill {\sigma }_{25}={\mu }_f(2,1,1)-{\mu }_f(0,1,0)·{\mu }_f(2,0,1),\\ \hfill & {\sigma }_{14}={\mu }_f(3,1,0)-\mu ·{\mu }_f(2,1,0),\hfill & \hfill {\sigma }_{34}={\mu }_f(3,2,0)-{\mu }_f(1,1,0)·{\mu }_f(2,1,0),\\ \hfill & {\sigma }_{15}={\mu }_f(3,0,1)-\mu ·{\mu }_f(2,0,1),\hfill & \hfill {\sigma }_{35}={\mu }_f(3,1,1)-{\mu }_f(1,1,0)·{\mu }_f(2,0,1),\\ \hfill & {\sigma }_{22}={\mu }_f(0,2,0)-{\mu }_f^2(0,1,0),\hfill & \hfill {\sigma }_{45}={\mu }_f(4,1,1)-{\mu }_f(2,1,0)·{\mu }_f(2,0,1),\\ \hfill & {\sigma }_{33}={\mu }_f(2,2,0)-{\mu }_f^2(1,1,0),\hfill & \hfill {\sigma }_{55}={\mu }_f(4,0,2)-{\mu }_f^2(2,0,1),\\ \hfill & {\sigma }_{44}={\mu }_f(4,2,0)-{\mu }_f^2(2,1,0).\hfill \end{array}$$

The proof of Theorem 3.3 is provided by Appendix A.5.

In the second step of deriving the asymptotics of ${\widehat{\lambda }}_{f,IG}$, we define the function $g(x_1,x_2,x_3,x_4,x_5):=x_1^2(2x_5+x_3)/(x_4-x_1^2{x}_2)$. Then, ${\widehat{\lambda }}_{f,IG}=g(\overline{\mathit{Z}})$ and $\lambda =g({\mathit{\mu }}_Z)$. Applying the Delta method [18] to Theorem 3.3, the following result follows.

Theorem 3.4

Let $X_1,\dots ,X_n$ be i. i. d. according to $IG(\mu ,\lambda )$, and define ${\vartheta }_{f,IG}:={\mu }_f(2,1,0)-{\mu }^2{\mu }_f(0,1,0)$. Then, ${\widehat{\lambda }}_{f,IG}$ is asymptotically normally distributed, where the asymptotic variance and bias, respectively, are given by$$\begin{array}{cc}\hfill {\sigma }_{f,IG}^2& =\frac{1}{n}[\frac{{\mu }^4}{{\vartheta }_f^2}[{\mu }_f(2,2,0)-{\mu }_f(1,1,0)(\frac{\lambda {\vartheta }_f}{{\mu }^2}+2{\mu }_f(2,0,1))+4{\mu }_f(3,1,1)+4{\mu }_f(4,0,2)\hfill \\ \hfill & -4{\mu }_f^2(2,0,1)]+\frac{2\lambda \mu }{{\vartheta }_f^2}[{\mu }^3({\mu }_f(1,2,0)+2{\mu }_f(2,1,1)-\frac{\lambda {\vartheta }_f}{{\mu }^2}·{\mu }_f(0,1,0))\hfill \\ \hfill & -\mu ({\mu }_f(3,2,0)+2{\mu }_f(4,1,1))+{\mu }_f(2,1,0)(4{\mu }_f(2,1,0)+4{\mu }_f(3,0,1)-\frac{\lambda {\vartheta }_f}{\mu })]\hfill \\ \hfill & +\frac{{\lambda }^2}{\mu {\vartheta }_f^2}[{\mu }^5·{\mu }_f(0,2,0)+{\mu }^3·{\mu }_f(0,1,0)({\vartheta }_f+{\mu }_f(2,1,0))-2{\mu }^3·{\mu }_f(2,2,0)\hfill \\ \hfill & +4{\mu }_f(2,1,0)({\mu }^2·{\mu }_f(1,1,0)-{\mu }^3·{\mu }_f(0,1,0)-{\mu }_f(3,1,0))\hfill \\ \hfill & +\mu (3{\mu }_f^2(2,1,0)+{\mu }_f(4,2,0))]],\hfill \end{array}$$ and$$\begin{array}{cc}\hfill {\mathbb{B}}_{f,IG}& =\frac{1}{n}[\frac{1}{\mu {\vartheta }_f^2}[{\mu }^2·{\mu }_f(2,1,0)({\mu }_f(2,1,0)+3{\mu }^2·{\mu }_f(0,1,0))+\lambda ({\mu }^5·{\mu }_f(0,2,0)\hfill \\ \hfill & +\mu ·{\mu }_f(4,2,0)-2{\mu }_f(3,1,0)·{\mu }_f(2,1,0)+4{\mu }^2·{\mu }_f(1,1,0)·{\mu }_f(2,1,0)\hfill \\ \hfill & -2{\mu }^3·{\mu }_f(2,2,0))]+\frac{\mu }{{\vartheta }_f^3}[{\mu }_f(2,1,0)(2{\mu }_f^2(2,1,0)-\mu ·{\mu }_f(3,2,0)\hfill \\ \hfill & +4{\mu }_f(2,1,0)·{\mu }_f(3,0,1)-2\mu ·{\mu }_f(4,1,1)+{\mu }^3({\mu }_f(1,2,0)+2{\mu }_f(2,1,1)))\hfill \\ \hfill & -{\mu }^2·{\mu }_f(0,1,0)(2{\mu }_f(3,1,0)(2{\mu }_f(2,0,1)+{\mu }_f(1,1,0))\hfill \\ \hfill & +2{\mu }_f^2(2,1,0)+4{\mu }_f(2,1,0)·{\mu }_f(3,0,1)+{\mu }^3({\mu }_f(1,2,0)+2{\mu }_f(2,1,1))\hfill \\ \hfill & -\mu ({\mu }_f(3,2,0)+2{\mu }_f(4,1,1)))]].\hfill \end{array}$$

The proof of Theorem 3.4 is provided by Appendix A.6.

Before we discuss the effect of f on bias and MSE of ${\widehat{\lambda }}_{f,IG}$, let us first consider the special cases of the ordinary MM and ML estimator. Their asymptotics are immediate consequences of Theorem 3.4. For the MM estimator ${\widehat{\lambda }}_{\mathit{MM}}$, we have to choose $f\equiv 1$ such that $f^{\prime }\equiv 0$. As a consequence,

3.6

$$\begin{array}{c}\hfill {\mu }_1(k,l,m)=\mathbb{E}[X^k]\text{if}m=0,\text{and 0 otherwise}.\end{array}$$This leads to a considerable simplification of Theorem 3.4, see Appendix A.7, which is summarized in the following corollary.

Corollary 3.5

Let $X_1,\dots ,X_n$ be i. i. d. according to $IG(\mu ,\lambda )$, then ${\widehat{\lambda }}_{\mathit{MM}}={\overline{X}}^3/S^2$ is asymptotically normally distributed with asymptotic variance ${\sigma }_{\mathit{MM}}^2=\frac{2}{n}\lambda (\lambda +3\mu )$ and bias ${\mathbb{B}}_{\mathit{MM}}=\frac{3}{n}(\lambda +3\mu )$.

While we are not aware of a reference providing these asymptotics, they can be verified by using Tweedie [27, p. 704]. There, normal asymptotics for the reciprocal $1/{\widehat{\lambda }}_{\mathit{MM}}$ are provided: $\sqrt{n}({\widehat{\lambda }}_{\mathit{MM}}^{-1}-{\lambda }^{-1})\sim N(0,2(1+3\mu /\lambda )/{\lambda }^2)$. Applying the Delta method with $g(x)=1/x$ and $g^{\prime }(x)=-1/x^2$ to it, we conclude that $\sqrt{n}({\widehat{\lambda }}_{\mathit{MM}}-\lambda )$ has the asymptotic variance ${\lambda }^4·2(1+3\mu /\lambda )/{\lambda }^2=2\lambda (\lambda +3\mu )$ like in Corollary 3.5.

Next, we consider the special case of the ML estimator ${\widehat{\lambda }}_{\mathit{ML}}$, which follows by choosing $f(x)=1/x$ such that $f^{\prime }(x)=-1/x^2$. Again, the joint moments ${\mu }_f(k,l,m)$ simplify a lot:

3.7

$$\begin{array}{c}\hfill {\mu }_{1/x}(k,l,m)={(-1)}^m\mathbb{E}[X^{k-l-2m}]\text{for all}k,l,m\in {\mathbb{N}}_0.\end{array}$$Together with Theorem 3.4, see Appendix A.8, we get the following corollary.

Corollary 3.6

Let $X_1,\dots ,X_n$ be i. i. d. according to $IG(\mu ,\lambda )$, then ${\widehat{\lambda }}_{\mathit{ML}}=\overline{X}/(\overline{X}·\overline{1/X}-1)$ is asymptotically normally distributed with asymptotic variance ${\sigma }_{\mathit{ML}}^2=\frac{2}{n}{\lambda }^2$ and bias ${\mathbb{B}}_{\mathit{ML}}=\frac{3}{n}\lambda$.

Comparing Corollaries 3.5 and 3.6, it is interesting to note that the MM estimator has larger asymptotic bias and variance than the ML estimator: ${\sigma }_{\mathit{MM}}^2={\sigma }_{\mathit{ML}}^2+\frac{6\lambda \mu }{n}$ and ${\mathbb{B}}_{\mathit{MM}}={\mathbb{B}}_{\mathit{ML}}+\frac{9\mu }{n}$. To verify the asymptotics of Corollary 3.6, note that the ML estimator ${\widehat{\lambda }}_{\mathit{ML}}$ has been shown to follow an inverted-${\chi }^2$ distribution: ${\widehat{\lambda }}_{\mathit{ML}}\sim n\lambda ·\text{Inv-}{\chi }_{n-1}^2$ [see 26, p. 368]. Using the formulae for mean and variance of $\text{Inv-}{\chi }_{n-1}^2$ [see 3, p. 431], we get$$\begin{array}{c}\hfill \mathbb{E}[{\widehat{\lambda }}_{\mathit{ML}}]=\frac{n}{n-3}\lambda =(1+\frac{3}{n-3})\lambda \approx (1+\frac{3}{n})\lambda ,\mathbb{V}[{\widehat{\lambda }}_{\mathit{ML}}]=\frac{2n^2}{{(n-3)}^2(n-5)}{\lambda }^2\approx \frac{2}{n}{\lambda }^2\end{array}$$for large n, which agrees with Corollary 3.6.

Table 2. Simulated vs. asymptotic mean and standard deviation of estimator ${\widehat{\lambda }}_{f,IG}$ from (3.3)

		Mean		Std. dev.				Mean		Std. dev.
f(x)	n	Sim	Asym	Sim	Asym.	f(x)	n	Sim	Asym	Sim	Asym.
$(\mu ,\lambda )=(1,1)$						$(\mu ,\lambda )=(3,1)$
1	100	1.103	1.120	0.258	0.283	1	100	1.213	1.300	0.365	0.447
	250	1.045	1.048	0.169	0.179		250	1.098	1.120	0.245	0.283
	500	1.023	1.024	0.122	0.126		500	1.053	1.060	0.180	0.200
$x^{-1/3}$	100	1.125	1.138	0.306	0.312	$x^{-1/3}$	100	1.243	1.294	0.442	0.451
	250	1.053	1.055	0.193	0.198		250	1.104	1.118	0.274	0.285
	500	1.027	1.028	0.138	0.140		500	1.055	1.059	0.194	0.202
$x^{-2/3}$	100	1.026	1.025	0.158	0.151	$x^{-2/3}$	100	1.020	1.019	0.160	0.156
	250	1.010	1.010	0.098	0.096		250	1.008	1.008	0.100	0.099
	500	1.006	1.005	0.069	0.068		500	1.004	1.004	0.070	0.070
1/x	100	1.031	1.030	0.149	0.141	1/x	100	1.031	1.030	0.149	0.141
	250	1.012	1.012	0.091	0.089		250	1.012	1.012	0.091	0.089
	500	1.006	1.006	0.064	0.063		500	1.006	1.006	0.064	0.063
$x^{-4/3}$	100	1.032	1.031	0.151	0.143	$x^{-4/3}$	100	1.033	1.032	0.154	0.146
	250	1.013	1.013	0.093	0.091		250	1.013	1.013	0.094	0.093
	500	1.007	1.006	0.065	0.064		500	1.007	1.006	0.066	0.065
$(\mu ,\lambda )=(1,3)$						$(\mu ,\lambda )=(3,3)$
1	100	3.172	3.180	0.595	0.600	1	100	3.309	3.360	0.775	0.849
	250	3.071	3.072	0.378	0.379		250	3.134	3.144	0.506	0.537
	500	3.036	3.036	0.267	0.268		500	3.069	3.072	0.366	0.379
$x^{-1/3}$	100	3.207	3.216	0.677	0.675	$x^{-1/3}$	100	3.374	3.415	0.918	0.937
	250	3.085	3.086	0.427	0.427		250	3.159	3.166	0.580	0.593
	500	3.043	3.043	0.301	0.302		500	3.081	3.083	0.413	0.419
$x^{-2/3}$	100	3.087	3.085	0.462	0.440	$x^{-2/3}$	100	3.078	3.076	0.475	0.454
	250	3.035	3.034	0.284	0.278		250	3.031	3.031	0.293	0.287
	500	3.018	3.017	0.199	0.197		500	3.017	3.015	0.206	0.203
1/x	100	3.093	3.090	0.448	0.424	1/x	100	3.093	3.090	0.448	0.424
	250	3.037	3.036	0.274	0.268		250	3.037	3.036	0.274	0.268
	500	3.019	3.018	0.192	0.190		500	3.019	3.018	0.192	0.190
$x^{-4/3}$	100	3.095	3.092	0.449	0.426	$x^{-4/3}$	100	3.097	3.094	0.453	0.430
	250	3.038	3.037	0.275	0.269		250	3.039	3.038	0.278	0.272
	500	3.020	3.018	0.193	0.191		500	3.020	3.019	0.195	0.192

Remark 3.7

To analyze the performance of the asymptotics provided by Theorem 3.4 (and that of the special cases discussed in Corollaries 3.5 and 3.6), if used as approximations to the true distribution of ${\widehat{\lambda }}_{f,IG}$ for finite sample size n, we did a simulation experiment with ${10}^5$ replications. The obtained results for various choices of $(\mu ,\lambda )$ and f(x) are summarized in Table 2. It can be recognized that the asymptotic approximations for mean and standard deviation generally agree quite well with their simulated counterparts. Only for the case $f(x)\equiv 1$ (the default MM estimator) and sample size $n=100$, we sometimes observe stronger deviations. But in the large majority of estimation scenarios, we have a close agreement such that the conclusions derived from the asymptotic expressions are meaningful for finite sample sizes as well.

In analogy to our discussion in Sect. 2, let us now analyze the performance of the Stein-MM estimator ${\widehat{\lambda }}_{f,IG}$ for the weight functions $f_a:(0;\infty )\to (0;\infty )$, $f_a(x)=x^a$ with $a\in \mathbb{R}\backslash \{-\frac{1}{2}\}$. Recall that this class of weight functions cover the default MM estimator ${\widehat{\lambda }}_{\mathit{MM}}$ for $a=0$ and the ML estimator ${\widehat{\lambda }}_{\mathit{ML}}$ for $a=-1$. The choice $a=-\frac{1}{2}$ (right in the middle between these two special cases) has to be excluded as it leads to a degenerate estimator ${\widehat{\lambda }}_{f_a,IG}$ according to (3.3). For this reason, the subsequent analyses in Figs. 3 and 4 are done separately for $a<-\frac{1}{2}$ (plots on left-hand side, covering the ML estimator) and $a>-\frac{1}{2}$ (plots on right-hand side, covering the MM estimator).

Fig. 3 [Images not available. See PDF.]

Plots of $n{\sigma }_{f_a,IG}^2$ and $n{\mathbb{B}}_{f_a,IG}$, where points indicate minimal variance and bias values. Scenarios $(\mu ,\lambda )=(1,3)$ with a$a\in (-2.5;-0.5)$ and b$a\in (-0.5;0.8)$, and $(\mu ,\lambda )=(3,1)$ with c$a\in (-2.5;-0.5)$ and (d) $a\in (-0.5;0.8)$. Dotted lines at $a=-1$ and $a=0$ correspond to default ML and MM estimator, respectively

Let us start with the analysis of asymptotic bias and variance in Fig. 3. The upper and lower panel consider two different example situations, namely $(\mu ,\lambda )=(1,3)$ and $(\mu ,\lambda )=(3,1)$, respectively, while left-hand and right-hand side are separated by the pole at $a=-\frac{1}{2}$. The right-hand side shows that the default MM estimator is neither (locally) optimal in terms of asymptotic bias nor in terms of variance. In fact, the optimal a for $a>-\frac{1}{2}$ is around $-0.1$ for $(\mu ,\lambda )=(1,3)$, and around $-0.2$ for $(\mu ,\lambda )=(3,1)$. However, comparing the actual values at the Y-axis to those of the plots on the left-hand side, we recognize that the asymptotic bias and variance get considerably smaller for some region with $a<-\frac{1}{2}$. In particular, the ML estimator is clearly superior to the MM estimator, and as shown by Figs. 3a and c, the ML estimator is even optimal in terms of the asymptotic variance. It has to be noted, however, that the curve corresponding to the asymptotic variance is rather flat around $a=-1$, so moderate deviations from $a=-1$ do not have a notable effect on the variance. Thus, it is important to also consider the optimum bias, which is reached for some a around $-0.65$ in both (a) and (c). So it appears to be advisable to choose an $a>-1$ for optimal overall estimation performance.

Fig. 4 [Images not available. See PDF.]

Plots of MSE${}_{f_a,IG}$, where points indicate minimal MSE values. Scenarios $(\mu ,\lambda )=(1,3)$ with a$a\in (-2.5;-0.5)$ and b$a\in (-0.5;0.8)$, and $(\mu ,\lambda )=(3,1)$ with c$a\in (-2.5;-0.5)$ and d$a\in (-0.5;0.8)$. Dotted lines at $a=-1$ and $a=0$ correspond to default ML and MM estimator, respectively

This is confirmed by Fig. 4, where the asymptotic MSE is shown for various sample sizes n and the same scenarios as in Fig. 3. While the ML estimator approaches the MSE-optimum for increasing n, we get an improved performance for smaller sample sizes if choosing $a\in (-1;-0.5)$ appropriately (e. g., $a\approx -0.8$ if $n\le 50$). Generally, an analogous recommendation holds for $a>-\frac{1}{2}$ in parts (b) and (d), with MSE-optima at a around $-0.1$ and $-0.2$, respectively, but much smaller MSE values can be reached for $a<-\frac{1}{2}$. To sum up, the default MM estimator (and more generally, Stein-MM estimators ${\widehat{\lambda }}_{f_a,IG}$ with $a>-\frac{1}{2}$) are not recommended for practice due to their rather large bias, variance, and thus MSE, while the ML estimator constitutes at least a good initial choice for estimating $\lambda$, being optimal in terms of asymptotic variance. However, unless the sample size n is very large, an improved MSE performance can be achieved by reducing a to an appropriate value in $(-1;-0.5)$ in the second step and by computing the corresponding Stein-MM estimate ${\widehat{\lambda }}_{f_a,IG}$.

Table 3. Runoff data from Example 3.8: Stein-MM estimates ${\widehat{\lambda }}_{f_a,IG}$ for different choices of function $f_a(x)=x^a$

a	$-$1.5	–1	$-$0.668	$-$0.125	$-$0.109	0	0.5
Notes		(i)	(ii)	(iii)	(iv)	(v)
${\widehat{\lambda }}_{f_a,IG}$	1.423	1.440	1.429	1.511	1.511	1.512	1.529
(i)   ML estimate, ${\sigma }_{f_a,IG}^2$-optimal				(ii)   ${\mathbb{B}}_{f_a,IG}$-optimal		(v)   MM estimate
(iii)   ${\mathbb{B}}_{f_a,IG}$-optimal given $a>-0.5$				(iv)   ${\sigma }_{f_a,IG}^2$-optimal given $a>-0.5$

Example 3.8

As an illustrative data example, let us consider the $n=25$ runoff amounts at Jug Bridge in Maryland [see $\mu$ is estimated by the sample mean as $\approx 0.803$, and using the ML estimator ${\widehat{\lambda }}_{f_{-1},IG}$ as an initial estimator for $\lambda$, we get the value $\approx 1.440$. As outlined before, this initial model fit might now be used for searching estimators with improved performance. Some examples (together with further estimates for comparative purposes) are summarized in Table 3. The ML estimator ($a=-1$) is also optimal in asymptotic variance, whereas the bias-optimal choice is obtained for a somewhat larger value of a, namely $a\approx -0.668$. The corresponding estimate is slightly lower than the ML estimate, similar to the value for $a=-1.5$, and can thus be seen as a fine-tuning of the initial estimate. By contrast, a notable change in the estimate happens if we turn to $a>-0.5$. The “constrained-optimal” choices (optimal given that $a>-0.5$) as well as the MM estimate lead to nearly the same values (around 1.51) and are thus visibly larger than the actually preferable estimates for $a<-0.5$. Also their variance and bias are about 2.5 times larger than those of the estimates for $a<-0.5$.

Stein Estimation of Negative-binomial Distribution

While the previous sections (and also the research by Ebner et al. [6]) solely focussed on continuous distributions, let us now turn to the case of discrete-valued random variables. Here, the most relevant type are count random variables X, having a quantitative range contained in ${\mathbb{N}}_0$. The probably most well-known distributions for counts are Poisson and binomial distributions, both depending on the (normalized) mean as their only model parameter. But as already discussed in Sect. 3, there is hardly any potential for finding a better estimator of the mean than the sample mean, so we do not further discuss these distributions. Instead, we focus on another popular count distribution, namely the NB-distribution with parameters $\nu >0$ and $\pi \in (0;1)$, abbreviated as $NB(\nu ,\pi )$. Such $X\sim NB(\nu ,\pi )$ has the range ${\mathbb{N}}_0$, probability mass function (pmf) $\mathbb{P}(X=x)=\left(\begin{array}{c}\nu +x-1\\ x\end{array}\right){(1-\pi )}^x{\pi }^{\nu }$, and mean $\mu =\mathbb{E}[X]=\frac{\nu (1-\pi )}{\pi }$. By contrast to the equidispersed Poisson distribution, its variance ${\sigma }^2:=\mathbb{V}[X]=\frac{\nu (1-\pi )}{{\pi }^2}$ is always larger than the mean (overdispersion), which is an important property for applications in practice. A detailed survey about the properties of and estimators for the $NB(\nu ,\pi )$-distribution can be found in Johnson et al. [11, Chapter 5]. Instead of the original parametrization by $(\nu ,\pi )$, it is often advantageous to consider either $(\mu ,\nu )$ or $(\mu ,\pi )$, where $\nu$ or $\pi$, respectively, serve as an additional dispersion parameter once the mean $\mu$ has been fixed. In case of the $(\mu ,\nu )$-parametrization, it holds that $\pi =\frac{\nu }{\nu +\mu }$ and $\mathbb{V}[X]=\frac{\nu +\mu }{\nu }\mu$, whereas we get $\nu =\frac{\pi \mu }{1-\pi }$ and $\mathbb{V}[X]=\frac{1}{\pi }\mu$ for the $(\mu ,\pi )$-parametrization. Besides the ease of interpretation, these parametrizations are advantageous in terms of parameter estimation. While MM estimation is rather obvious, namely $\mu$ by $\overline{X}$ and ${\nu }_{\mathit{MM}}={\overline{X}}^2/(S^2-\overline{X})$, ${\pi }_{\mathit{MM}}=\overline{X}/S^2$, ML estimation is generally demanding as there does not exist a closed-form solution, see the discussion by Kemp & Kemp [12], i. e., numerical optimization is necessary. However, there is an important exception: the NB’s ML estimator of the mean $\mu$ is given by $\overline{X}$ [12, p. 867], i. e., $\overline{X}$ is both the MM and ML estimator with its known appealing performance. So it suffices to find an adequate estimator for $\nu$ or $\pi$, respectively, the ML estimators of which do not have a closed-form expression.

These difficulties in estimating $\nu$ or $\pi$, respectively, serve as our motivation for deriving a generalized MM estimator. For this purpose, we consider the NB’s Stein identity according to Brown & Phillips [5, Lemma 1], which can be expressed as either

4.1

$$\begin{array}{cc}& X\sim NB(\nu ,\frac{\nu }{\nu +\mu })\text{iff}\hfill \\ \hfill & \nu \mathbb{E}[Xf(X)-\mu f(X+1)]=\mu \mathbb{E}[X(f(X+1)-f(X))],\hfill \end{array}$$

4.2

$$\begin{array}{c}\hfill \text{or}\\ \hfill & X\sim NB(\frac{\pi \mu }{1-\pi },\pi )\text{iff}\hfill \\ \hfill & \pi \mathbb{E}[(X-\mu )f(X+1)]=\mathbb{E}[X(f(X+1)-f(X))],\hfill \end{array}$$for any function f such that $\mathbb{E}[|f(X)|]$, $\mathbb{E}[|f(X+1)|]$ exist. Note that the discrete difference $\mathrm{\Delta }f(x):=f(x+1)-f(x)$ in (4.1) and (4.2) plays a similar role as the continuous derivative $f^{\prime }(x)$ in the previous Stein identities (2.1) and (3.2).

Stein-MM estimators are now derived by solving (4.1) in $\nu$ or (4.2) in $\pi$, respectively, and by using again sample moments $\overline{h(X)}=\frac{1}{n}{\sum }_{i=1}^{n}h(X_i)$ instead of the involved population moments $\mathbb{E}[h(X)]$ (with $\mu$ being estimated by $\overline{X}$). As a result, the (closed-form) classes of Stein-MM estimators for $\nu$ and $\pi$ are obtained as

4.3

$$\begin{array}{c}\hfill {\widehat{\nu }}_{f,NB}=\frac{\overline{X}\overline{X\mathrm{\Delta }f(X)}}{\overline{Xf(X)}-\overline{X}\overline{f(X+1)}},{\widehat{\pi }}_{f,NB}=\frac{\overline{X\mathrm{\Delta }f(X)}}{\overline{Xf(X+1)}-\overline{X}\overline{f(X+1)}}.\end{array}$$Note that the choice $f(x)=x$ (hence $\mathrm{\Delta }f(x)=1$) leads to the default MM estimators given above. The ML estimators are not covered by (4.3) this time, because they do not have a closed-form expression at all. Note, however, that the so-called “weighted-mean estimator” for $\nu$ in (2.6) of Kemp & Kemp [12], which was motivated as a kind of approximate ML estimator, is covered by (4.3), namely by choosing $f_{\alpha }(x)={\alpha }^x$ with $\alpha \in (0;1)$. It is also worth pointing to Savani & Zhigljavsky [17], who define an estimator of $\nu$ based on the moment $\overline{f(X)}$ for some specified f; their approach, however, usually does not lead to a closed-form estimator.

For deriving the asymptotic distribution of the general Stein-MM estimator ${\widehat{\nu }}_{f,NB}$ or ${\widehat{\pi }}_{f,NB}$, respectively, we first define the vectors ${\mathit{Z}}_i$ with $i=1,\dots ,n$ as

4.4

$$\begin{array}{c}\hfill {\mathit{Z}}_i:=(X_i,f(X_i+1),X_{i}f(X_i),X_{i}f(X_i+1){)}^{\top }.\end{array}$$Their mean equals

4.5

$$\begin{array}{c}\hfill {\mathit{\mu }}_Z:=\mathbb{E}[{\mathit{Z}}_i]=(\mu ,{\tilde{\mu }}_f(0,0,1),{\tilde{\mu }}_f(1,1,0),{\tilde{\mu }}_f(1,0,1){)}^{\top },\end{array}$$where we define ${\tilde{\mu }}_f(k,l,m):=\mathbb{E}[X^k·f{(X)}^l·f{(X+1)}^m]$ for any $k,l,m\in {\mathbb{N}}_0$. Then, the following CLT holds.

Theorem 4.1

If $X_1,\dots ,X_n$ are i. i. d. according to a negative binomial distribution, then the sample mean $\overline{\mathit{Z}}$ of ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_n$ according to (4.4) is asymptotically normally distributed as$$\begin{array}{c}\hfill \sqrt{n}(\overline{\mathit{Z}}-{\mathit{\mu }}_Z)\stackrel{\text{d}}{\to }N(\mathbf{0},\mathbf{\Sigma })\text{with}\mathbf{\Sigma }={({\sigma }_{\mathit{ij}})}_{i,j=1,\dots ,4},\end{array}$$where $N(\mathbf{0},\mathbf{\Sigma })$ denotes the multivariate normal distribution, and where the covariances are given as$$\begin{array}{ccc}& {\sigma }_{11}={\sigma }^2,\hfill & \hfill {\sigma }_{23}={\tilde{\mu }}_f(1,1,1)-{\tilde{\mu }}_f(0,0,1)·{\tilde{\mu }}_f(1,1,0),\\ \hfill & {\sigma }_{12}={\tilde{\mu }}_f(1,0,1)-\mu ·{\tilde{\mu }}_f(0,0,1),\hfill & \hfill {\sigma }_{24}={\tilde{\mu }}_f(1,0,2)-{\tilde{\mu }}_f(0,0,1)·{\tilde{\mu }}_f(1,0,1),\\ \hfill & {\sigma }_{13}={\tilde{\mu }}_f(2,1,0)-\mu ·{\tilde{\mu }}_f(1,1,0),\hfill & \hfill {\sigma }_{33}={\tilde{\mu }}_f(2,2,0)-{\tilde{\mu }}_f^2(1,1,0),\\ \hfill & {\sigma }_{14}={\tilde{\mu }}_f(2,0,1)-\mu ·{\tilde{\mu }}_f(1,0,1),\hfill & \hfill {\sigma }_{34}={\tilde{\mu }}_f(2,1,1)-{\tilde{\mu }}_f(1,1,0)·{\tilde{\mu }}_f(1,0,1),\\ \hfill & {\sigma }_{22}={\tilde{\mu }}_f(0,0,2)-{\tilde{\mu }}_f^2(0,0,1),\hfill & \hfill {\sigma }_{44}={\tilde{\mu }}_f(2,0,2)-{\tilde{\mu }}_f^2(1,0,1).\end{array}$$

The proof of Theorem 4.1 is provided by Appendix A.9.

In the second step of deriving the Stein-MM estimators’ asymptotics, we define the function

• $g_{\nu }(x_1,x_2,x_3,x_4):=x_1(x_4-x_3)/(x_3-x_1{x}_2)$ for ${\widehat{\nu }}_{f,NB}$,

• $g_{\pi }(x_1,x_2,x_3,x_4):=(x_4-x_3)/(x_4-x_1{x}_2)$ for ${\widehat{\pi }}_{f,NB}$.

Then, ${\widehat{\nu }}_{f,NB}=g_{\nu }(\overline{\mathit{Z}})$, $\nu =g_{\nu }({\mathit{\mu }}_Z)$, and ${\widehat{\pi }}_{f,NB}=g_{\pi }(\overline{\mathit{Z}})$, $\pi =g_{\pi }({\mathit{\mu }}_Z)$, respectively, holds. Applying the Delta method [18] to Theorem 4.1, the following theorems follow.

Theorem 4.2

Let $X_1,\dots ,X_n$ be i. i. d. according to $NB(\nu ,\frac{\nu }{\nu +\mu })$, and define ${\eta }_1:={\tilde{\mu }}_f(1,1,0)-\mu ·{\tilde{\mu }}_f(0,0,1)$. Then, ${\widehat{\nu }}_{f,NB}$ is asymptotically normally distributed, where the asymptotic variance and bias, respectively, are given by$$\begin{array}{cc}\hfill \sigma & {}_{f,NB}^2=\frac{1}{{\eta }_1^{4}n}[{\mu }^4({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(2,0,2)-2{\tilde{\mu }}_f(2,1,1)+{\tilde{\mu }}_f(2,2,0))\hfill \\ \hfill & -2({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(1,0,2)-{\tilde{\mu }}_f(1,1,1)))\hfill \\ \hfill & +{\tilde{\mu }}_f(0,0,2)({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0){)}^2)\hfill \\ \hfill & +2{\mu }^3({\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(1,0,1)({\tilde{\mu }}_f(2,1,1)-{\tilde{\mu }}_f(2,2,0))\hfill \\ \hfill & -{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(2,0,2)-{\tilde{\mu }}_f(2,1,1))+{\tilde{\mu }}_f(1,0,2){\tilde{\mu }}_f(1,1,0))\hfill \\ \hfill & -{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(1,1,1)+{\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(1,0,2)+{\tilde{\mu }}_f(1,1,1)))\hfill \\ \hfill & +{\mu }^2({\tilde{\mu }}_f(1,1,0)(2{\tilde{\mu }}_f(1,0,1)(-{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,0,1)+{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,1,0)\hfill \\ \hfill & +{\tilde{\mu }}_f^2(1,1,0)-{\tilde{\mu }}_f(2,1,1))+{\tilde{\mu }}_f(1,1,0)(2{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,0,1)-2{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,1,0)\hfill \\ \hfill & +{\tilde{\mu }}_f(2,0,2))+2{\tilde{\mu }}_f^3(1,0,1)-4{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(1,1,0))+{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(2,2,0))\hfill \\ \hfill & -2\mu {\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(2,1,0)-{\tilde{\mu }}_f(1,1,0){\tilde{\mu }}_f(2,0,1))\hfill \\ \hfill & +{\tilde{\mu }}_f^2(1,1,0){\sigma }^2({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0){)}^2]\hfill \end{array}$$ and$$\begin{array}{cc}\hfill \mathbb{B}& {}_{f,NB}=-\frac{1}{{\eta }_1^{3}n}[{\mu }^3({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(1,0,2)-{\tilde{\mu }}_f(1,1,1))+{\tilde{\mu }}_f(0,0,2)({\tilde{\mu }}_f(1,1,0)\hfill \\ \hfill & -{\tilde{\mu }}_f(1,0,1)))-{\mu }^2({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(2,1,1)-{\tilde{\mu }}_f(2,2,0))-2{\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(1,1,1)\hfill \\ \hfill & +{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(1,0,2)+{\tilde{\mu }}_f(1,1,1)))\hfill \\ \hfill & +\mu ({\tilde{\mu }}_f(1,0,1)({\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,1,0)+2{\tilde{\mu }}_f^2(1,1,0)-{\tilde{\mu }}_f(2,2,0))\hfill \\ \hfill & +{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,0,1)-2{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,1,0)+{\tilde{\mu }}_f(2,1,1))\hfill \\ \hfill & -2{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(1,1,0))+{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(0,0,1){\sigma }^2({\tilde{\mu }}_f(1,1,0)-{\tilde{\mu }}_f(1,0,1))\hfill \\ \hfill & +{\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(2,1,0)-{\tilde{\mu }}_f(1,1,0){\tilde{\mu }}_f(2,0,1))].\hfill \end{array}$$

The proof of Theorem 4.2 is provided by Appendix A.10.

Theorem 4.3

Let $X_1,\dots ,X_n$ be i. i. d. according to $NB(\frac{\pi \mu }{1-\pi },\pi )$, and define ${\eta }_2:={\tilde{\mu }}_f(1,0,1)-\mu ·{\tilde{\mu }}_f(0,0,1)$. Then, ${\widehat{\pi }}_{f,NB}$ is asymptotically normally distributed, where the asymptotic variance and bias, respectively, are given by$$\begin{array}{cc}\hfill \sigma & {}_{f,NB}^2=\frac{1}{{\eta }_2^{4}n}[{\mu }^2({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(2,0,2)-2{\tilde{\mu }}_f(2,1,1)+{\tilde{\mu }}_f(2,2,0))\hfill \\ \hfill & -2({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(1,0,2)-{\tilde{\mu }}_f(1,1,1)))\hfill \\ \hfill & +{\tilde{\mu }}_f(0,0,2)({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f^2(1,1,0)))\hfill \\ \hfill & -2\mu ({\tilde{\mu }}_f^2(0,0,1)({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(2,0,1)-{\tilde{\mu }}_f(2,1,0))\hfill \\ \hfill & -{\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f^3(1,0,1)-2{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(1,1,0)+{\tilde{\mu }}_f(1,0,1)({\tilde{\mu }}_f^2(1,1,0)\hfill \\ \hfill & +{\tilde{\mu }}_f(2,1,1)-{\tilde{\mu }}_f(2,2,0))+{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f(2,1,1)-{\tilde{\mu }}_f(2,0,2)))\hfill \\ \hfill & +({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(1,1,1)-{\tilde{\mu }}_f(1,0,2){\tilde{\mu }}_f(1,1,0)))\hfill \\ \hfill & +{\tilde{\mu }}_f^2(0,0,1){\sigma }^2({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0){)}^2+{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(2,2,0)\hfill \\ \hfill & -2{\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(1,0,1)-{\tilde{\mu }}_f(1,1,0))({\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(2,1,0)-{\tilde{\mu }}_f(1,1,0){\tilde{\mu }}_f(2,0,1))\hfill \\ \hfill & -2{\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(1,1,0){\tilde{\mu }}_f(2,1,1)+{\tilde{\mu }}_f^2(1,1,0){\tilde{\mu }}_f(2,0,2)]\hfill \end{array}$$ and$$\begin{array}{cc}\hfill \mathbb{B}& {}_{f,NB}=-\frac{1}{{\eta }_2^{3}n}[{\mu }^2({\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(1,0,2)-{\tilde{\mu }}_f(1,1,1))\hfill \\ \hfill & +{\tilde{\mu }}_f(0,0,2)({\tilde{\mu }}_f(1,1,0)-{\tilde{\mu }}_f(1,0,1)))\hfill \\ \hfill & +\mu ({\tilde{\mu }}_f^2(0,0,1)({\tilde{\mu }}_f(2,0,1)-{\tilde{\mu }}_f(2,1,0))\hfill \\ \hfill & +{\tilde{\mu }}_f(0,0,1)({\tilde{\mu }}_f(1,0,1)({\tilde{\mu }}_f(1,1,0)-{\tilde{\mu }}_f(1,0,1))-{\tilde{\mu }}_f(2,0,2)+{\tilde{\mu }}_f(2,1,1))\hfill \\ \hfill & +{\tilde{\mu }}_f(1,0,1)({\tilde{\mu }}_f(1,0,2)+{\tilde{\mu }}_f(1,1,1))-2{\tilde{\mu }}_f(1,0,2){\tilde{\mu }}_f(1,1,0))\hfill \\ \hfill & +{\tilde{\mu }}_f(1,1,0)({\tilde{\mu }}_f^2(0,0,1){\sigma }^2-2{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(2,0,1)+{\tilde{\mu }}_f(2,0,2))\hfill \\ \hfill & +{\tilde{\mu }}_f(0,0,1){\tilde{\mu }}_f(1,0,1)(-{\tilde{\mu }}_f(0,0,1){\sigma }^2+{\tilde{\mu }}_f(2,0,1)+{\tilde{\mu }}_f(2,1,0))\hfill \\ \hfill & -{\tilde{\mu }}_f^3(1,0,1)+{\tilde{\mu }}_f^2(1,0,1){\tilde{\mu }}_f(1,1,0)-{\tilde{\mu }}_f(1,0,1){\tilde{\mu }}_f(2,1,1)].\hfill \end{array}$$

The proof of Theorem 4.3 is provided by Appendix A.11.

Our first special case shall be the function $f_{\alpha }(x)={\alpha }^x$ with $\alpha \in (0;1)$, which is inspired by Kemp & Kemp [12]. For evaluating the asymptotics in Theorems 4.1–4.3, we need to compute the moments$$\begin{array}{c}\hfill {\tilde{\mu }}_{f_{\alpha }}(k,l,m)=\mathbb{E}[X^k·f_{\alpha }{(X)}^l·f_{\alpha }{(X+1)}^m]={\alpha }^m\mathbb{E}[X^k·{\alpha }^{(l+m)X}].\end{array}$$As shown in the following, this can be done by explicit closed-form expressions. The idea is to utilize the probability generating function (pgf) of the NB-distribution,$$\begin{array}{c}\hfill \text{pgf}(z):=\mathbb{E}[z^X]=(\frac{\pi }{1-(1-\pi )z}{)}^{\nu },\end{array}$$together with the following property:$$\begin{array}{c}\hfill \mathbb{E}[X_{(k)}{z}^X]=z^k·\frac{d^k}{d{z}^k}\text{pgf}(z),\end{array}$$where $x_{(r)}:=x\cdots (x-r+1)$ for $r\in {\mathbb{N}}_0$ denote the falling factorials. The main result is summarized by the following lemma.

Lemma 4.4

Let $X\sim NB(\nu ,\pi )$. For the mixed factorial moments, we have$$\begin{array}{c}\hfill \mathbb{E}[X_{(k)}{z}^X]=\frac{{(1-\pi )}^k{(\nu +k-1)}_{(k)}{z}^k}{(1-(1-\pi )z{)}^k}\text{pgf}(z).\end{array}$$

The proof of Lemma 4.4 is provided by Appendix A.12. The factorial moments are easily transformed into raw moments by using the relation $x^k={\sum }_{j=0}^k{S}_{k,j}{x}_{(j)}$, where $S_{k,j}$ are the Sterling numbers of the second kind (see [11], p. 12). Then, $\mathbb{E}[X^k·{\alpha }^{(l+m)X}]$ follows by plugging-in $z={\alpha }^{l+m}$ into Lemma 4.4.

Fig. 5 [Images not available. See PDF.]

Stein-MM estimator ${\widehat{\nu }}_{f,NB}$ for $\mu =2.5$. Plots of $n{\sigma }_{f,NB}^2$ and $n{\mathbb{B}}_{f,NB}$ for parametrization (4.1), where points indicate minimal variance and bias values. Weighting function (a)–(b) $f(x)={\alpha }^x$ with $\alpha \in (0,1)$, and (c)–(d) $f(x)={(x+1)}^a$ with $a\in (-1,1.5)$. The gray graphs in (c)–(d) correspond to the comparative choice $f(x)=x^a$, which leads to the default MM estimator for $a=1$ (dotted lines)

While general closed-form formulae are possible in this way for ${\tilde{\mu }}_{f_{\alpha }}(k,l,m)$ as well as for Theorems 4.1–4.3, the obtained results are very complex such that we decided to omit the final expressions. Instead, we compute ${\tilde{\mu }}_{f_{\alpha }}(k,l,m)$ and, thus, the expressions of Theorems 4.1–4.3 numerically. This is easily done in practice, in fact for any reasonable choice of the function f, by computing$$\begin{array}{c}\hfill {\tilde{\mu }}_{f_{\alpha }}(k,l,m)\approx \sum _{x=0}^M{x}^k·f{(x)}^l·f{(x+1)}^m·P(X=x),\end{array}$$where the upper summation limit M is chosen sufficiently large, e. g., such that $M^k·f{(M)}^l·f{(M+1)}^m·P(X=M)$ falls below a specified tolerance limit. In this way, we generated the illustrative graphs in Figs. 5 (estimator ${\widehat{\nu }}_{f,NB}$) and 6 (estimator ${\widehat{\pi }}_{f,NB}$). There, parts (a)–(b) always refer to the above choice $f_{\alpha }(x)={\alpha }^x$, and clear minima for variance and bias for $f_{\alpha }(x)={\alpha }^x$ can be recognized. To be able to compare with the respective default MM estimator, we did analogous computations for $f_a(x)=x^a$ (where $a=1$ for the default MM estimator), which, however, is only defined for $a>0$ as X becomes zero with positive probability. As can be seen from the gray curves in parts (c)–(d), variance and basis usually do not attain a local minimum for $a>0$. Therefore, parts (c)–(d) mainly focus on a slight modification of the weight function, namely $f_{a,1}(x)={(x+1)}^a$, which is also well-defined for $a<0$.

Fig. 6 [Images not available. See PDF.]

Stein-MM estimator ${\widehat{\pi }}_{f,NB}$ for $\mu =2.5$. Plots of $n{\sigma }_{f,NB}^2$ and $n{\mathbb{B}}_{f,NB}$ for parametrization (4.2), where points indicate minimal variance and bias values. Weighting function (a)–(b) $f_{\alpha }(x)={\alpha }^x$ with $\alpha \in (0,1)$, and (c)–(d) $f_{a,1}(x)={(x+1)}^a$ with $a\in (-1,1.5)$. The gray graphs in (c)–(d) correspond to the comparative choice $f_a(x)=x^a$, which leads to the default MM estimator for $a=1$ (dotted lines)

Table 4. Stein-MM estimators ${\widehat{\nu }}_{f,NB}$ and ${\widehat{\pi }}_{f,NB}$, optimal choices for $\alpha$ or a, respectively (columns “$a_{\mathit{opt}}$”), and corresponding minimal value of variance and bias (columns “min.”)

			${\widehat{\nu }}_{f,NB}$, optimal choice for				${\widehat{\pi }}_{f,NB}$, optimal choice for
$\mu =2.5$			Variance		Bias		Variance		Bias
f(x)	$\nu$	$\pi$	$a_{\mathit{opt}}$	Min	$a_{\mathit{opt}}$	Min	$a_{\mathit{opt}}$	Min	$a_{\mathit{opt}}$	Min
${\alpha }^x$	1.0	0.286	0.751	5.095	0.620	5.093	0.751	0.271	0.544	0.920
	1.5	0.375	0.771	14.271	0.668	10.102	0.771	0.407	0.595	1.144
	2.5	0.500	0.805	59.113	0.727	26.022	0.805	0.641	0.650	1.475
${(x+1)}^a$	1.0	0.286	$-$0.489	5.002	$-$0.990	5.311	$-$0.489	0.267	$-$0.990	0.985
	1.5	0.375	$-$0.332	14.130	$-$0.861	10.440	$-$0.332	0.404	$-$0.990	1.205
	2.5	0.500	$-$0.097	58.908	$-$0.470	26.680	$-$0.097	0.639	$-$0.896	1.544

The optimal choices for $\alpha$ and a, respectively, lead to very similar variance and bias values, see Table 4. While ${\alpha }^x$ leads to a slightly larger variance than ${(x+1)}^a$, its optimal bias is visibly lower. For both choices of f, however, the optimal Stein-MM estimators perform clearly better than the default MM estimator, see the dotted line at $a=1$ in parts (c)–(d) in Figs. 5 and 6. Altogether, also in view of the fact that explicit closed-form expressions are possible for ${\alpha }^x$ (although being rather complex), we prefer to use $f_{\alpha }(x)={\alpha }^x$ as the weighting function, in accordance to Kemp & Kemp [12]. For this choice, we also did a simulation experiment with ${10}^5$ replications (in analogy to Remark 3.7), in order to check the finite-sample performance of the asymptotic expressions for variance and bias. We generally observed a very good agreement between asymptotic and simulated values. Especially for the estimator ${\widehat{\pi }}_{f_{\alpha },NB}$, the asymptotic approximations show an excellent performance, whereas the estimator ${\widehat{\nu }}_{f_{\alpha },NB}$ sometimes leads to extreme estimates if $\nu =2.5$ and $n=100$. But except these few outlying estimates, also ${\widehat{\nu }}_{f_{\alpha },NB}$ is well described by the asymptotic formulae. Detailed simulation results are available from the authors upon request.

Table 5. Counts of red mites on apple leaves from Example 4.5: Stein-MM estimates ${\widehat{\nu }}_{f_{\alpha },NB}$ (upper part) and ${\widehat{\pi }}_{f_{\alpha },NB}$ (lower part) for different choices of function $f_{\alpha }(x)={\alpha }^x$

$\alpha$	0.25	0.5	0.530	0.690	0.75	–
Notes			(i)	(ii)		(iii)
${\widehat{\nu }}_{f_{\alpha },NB}$	0.967	0.963	0.967	1.009	1.032	1.167
(i)   ${\mathbb{B}}_{f_{\alpha },NB}$-optimal		(ii)   ${\sigma }_{f_{\alpha },NB}^2$-optimal		(iii)   MM estimate
$\alpha$	0.222	0.25	0.5	0.690	0.75	—
Notes	(i)			(ii)		(iii)
${\widehat{\pi }}_{f_{\alpha },NB}$	0.459	0.457	0.456	0.468	0.474	0.504
(i)   ${\mathbb{B}}_{f_{\alpha },NB}$-optimal		(ii)   ${\sigma }_{f_{\alpha },NB}^2$-optimal		(iii)   MM estimate

Example 4.5

As an illustrative data example, let us consider $n=150$ counts of red mites on apple leaves (see [16], p. 271), who confirmed “a good fit of the negative binomial” for these data. The parameter $\mu$ is estimated by the sample mean as $\approx 1.147$. In case of the $(\mu ,\nu )$-parametrization, we use the ordinary MM estimator as an initial estimator for $\nu$, leading to the value $\approx 1.167$. Based on this initial model fit, we search for Stein-MM estimators with $f_{\alpha }(x)={\alpha }^x$ having an improved performance. The resulting estimates (together with further estimates for comparative purposes) are summarized in the upper part of Table 5. It can be seen that the initial estimate (last column) is corrected downwards to a value close to 1 (i. e., we essentially end up with the special case of a geometric distribution). Here, it is interesting to note that the numerically computed ML estimate as reported in Rueda & O’Reilly [16], also leads to such a result, namely to the value 1.025. In this context, we also recall Kemp & Kemp [12] who proposed the choice $f_{\alpha }(x)={\alpha }^x$ to get a closed-form approximate ML estimator for $\nu$.

We repeated the aforementioned estimation procedure also for the $(\mu ,\pi )$-parameterization, starting with the initial MM-estimate $\approx 0.504$ for $\pi$, see the lower part of Table 5. Again, the initial estimate is corrected downwards to a value around 0.46.

Conclusions

In this article, we demonstrated how Stein characterizations of (continuous or discrete) distributions can be utilized to derive improved moment estimators of model parameters. The main idea is to first choose an appropriate parametric class of weighting functions. Then, the final parameter values are determined such that the resulting Stein-MM estimator has asymptotically optimal mean or variance properties within the considered class, and certainly also lower variance and bias than existing parameter estimators. Here, the optimal choice from the given class of weighting functions is implemented based closed-form expressions for asymptotic distributions, possibly together with numerical optimization routines. The whole procedure was exemplified for three types of distribution: the continuous exponential and inverse Gaussian distribution, as well as the discrete negative-binomial distribution. For all these distribution families, we observed an appealing performance in various aspects, and we also demonstrated the application of our findings to real-world data examples.

Taking together the present article and the contributions by Arnold et al. [2], Wang & Weiß [28], Ebner et al. [6], Stein-MM estimators are now available for a wide class of continuous distributions. For discrete distributions, however, only the negative-binomial distribution (see Sect. 4 before) and the discrete Lindley distributions [see 28] have been considered so far. Thus, future research should be directed towards Stein-MM estimators for further common types of discrete distribution. Our research also gives rise to several further directions for future research. While our main focus was on selecting the weight function with respect to minimal bias or variance, we also briefly pointed out in Sect. 2 that such a choice could also be motivated by robustness to outliers. In fact, there are some analogies to “M-estimation” as introduced by Huber [9]. It appears to be promising to analyze if robustified MM estimators can be achieved by suitable classes of weighting function. As another direction for future research (briefly sketched in Sect. 2), the performance of GoF-tests based on Stein-MM estimators should be investigated. Finally, one should analyze Stein-MM estimators in a regression or time-series context.

Acknowledgements

The author thank the two referees for their useful comments on an earlier draft of this article.

Author Contributions

SN and CHW contributed to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript.

Funding

Open Access funding enabled and organized by Projekt DEAL. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 437270842.

Availability of data and materials

The datasets analyzed during the current study are available in Folks & Chhikara [8, p. 272] and Rueda & O’Reilly [16, p. 271], respectively.

Declarations

Conflict of interest

The authors declare they have no conflict of interest.

Consent for publication

The authors hereby consent to publication of their work upon its acceptance.

References

1. Anastasiou, A; Barp, A; Briol, F-X; Ebner, B; Gaunt, RE; Ghaderinezhad, F; Gorham, J; Gretton, A; Ley, C; Liu, Q; Mackey, L; Oates, CJ; Reinert, G; Swan, Y. Stein’s method meets computational statistics: a review of some recent developments. Stat. Sci.; 2023; 38, 1 pp. 120-139.MathSciNet ID: 4534646[DOI: https://dx.doi.org/10.1214/22-STS863]

2. Arnold, BC; Castillo, E; Sarabia, JM. A multivariate version of Stein’s identity with applications to moment calculations and estimation of conditionally specified distributions. Commun. Stat. Theory Methods; 2001; 30, 12 pp. 2517-2542.MathSciNet ID: 1877354[DOI: https://dx.doi.org/10.1081/STA-100108446]

3. Bernardo, JM; Smith, AFM. Bayesian Theory; 1994; New York, John Wiley & Sons Inc: [DOI: https://dx.doi.org/10.1002/9780470316870]

4. Betsch, S; Ebner, B; Nestmann, F. Characterizations of non-normalized discrete probability distributions and their application in statistics. Electronic Journal of Statistics; 2022; 16, 1 pp. 1303-1329.MathSciNet ID: 4381061[DOI: https://dx.doi.org/10.1214/22-EJS1983]

5. Brown, TC; Phillips, MJ. Negative binomial approximation with Stein’s method. Methodol. Comput. Appl. Probab.; 1999; 1, 4 pp. 407-421.MathSciNet ID: 1770372[DOI: https://dx.doi.org/10.1023/A:1010094221471]

6. Ebner, B., Fischer, A., Gaunt, R.E., Picker, B., Swan, Y.: Point estimation through Stein’s method. arXiv:2305.19031 (2023)

7. Elfessi, A; Reineke, DM. A Bayesian look at classical estimation: the exponential distribution. J. Stat. Educ.; 2001; 9, 1 5. [DOI: https://dx.doi.org/10.1080/10691898.2001.11910648]

8. Folks, JL; Chhikara, RS. The inverse Gaussian distribution and its statistical application – a review. J. Roy. Stat. Soc. B; 1978; 40, 3 pp. 263-289.MathSciNet ID: 522211[DOI: https://dx.doi.org/10.1111/j.2517-6161.1978.tb01039.x]

9. Huber, PJ. Robust estimation of a location parameter. Ann. Math. Stat.; 1964; 35, 1 pp. 73-101.MathSciNet ID: 161415[DOI: https://dx.doi.org/10.1214/aoms/1177703732]

10. Johnson, N.L., Kotz, S., Balakrishnan, N.: Continuous Univariate Distributions, Volume 1. 2^nd edition, John Wiley & Sons, Inc., New York (1995)

11. Johnson, N.L., Kemp, A.W., Kotz, S.: Univariate Discrete Distributions. 3rd edition, John Wiley & Sons, Inc., New York (2005)

12. Kemp, AW; Kemp, CD. A rapid and efficient estimation procedure for the negative binomial distribution. Biomed. J.; 1987; 29, 7 pp. 865-873.MathSciNet ID: 932069

13. Koudou, AE; Ley, C. Characterizations of GIG laws: a survey. Probab. Surv.; 2014; 11, pp. 161-176.MathSciNet ID: 3264557[DOI: https://dx.doi.org/10.1214/13-PS227]

14. Kubokawa, T.: Stein’s identities and the related topics: an instructive explanation on shrinkage, characterization, normal approximation and goodness-of-fit. Japanese Journal of Statistics and Data Science, in press (2024)

15. Landsman, Z; Valdez, EA. The tail Stein’s identity with applications to risk measures. N. Am. Actuar. J.; 2016; 20, 4 pp. 313-326.MathSciNet ID: 3582061[DOI: https://dx.doi.org/10.1080/10920277.2016.1237879]

16. Rueda, R; O’Reilly, F. Tests of fit for discrete distributions based on the probability generating function. Commun. Stat.-Simul. Comput.; 1999; 28, 1 pp. 259-274.MathSciNet ID: 1702548[DOI: https://dx.doi.org/10.1080/03610919908813547]

17. Savani, V; Zhigljavsky, AA. Efficient estimation of parameters of the negative binomial distribution. Commun. Stat.-Theory Methods; 2006; 35, 5 pp. 767-783.MathSciNet ID: 2256250[DOI: https://dx.doi.org/10.1080/03610920500501346]

18. Serfling, RJ. Approximation Theorems of Mathematical Statistics; 1980; New York, John Wiley & Sons Inc: [DOI: https://dx.doi.org/10.1002/9780470316481]

19. Seshadri, V. The Inverse Gaussian Distribution; 1999; New York, Springer: [DOI: https://dx.doi.org/10.1007/978-1-4612-1456-4]

20. Shuster, JJ. Nonparametric optimality of the sample mean and sample variance. Am. Stat.; 1982; 36, 3 pp. 176-178.MathSciNet ID: 664099[DOI: https://dx.doi.org/10.2307/2683173]

21. Stein, C. A bound for the error in the normal approximation to the distribution of a sum of dependent random variables. Proc. Sixth Berkeley Sympos. Math. Stat. Probab.; 1972; 2, pp. 583-602.MathSciNet ID: 402873

22. Stein, C.: Approximate Computation of Expectations. IMS Lecture Notes, Volume 7, Hayward, California (1986)

23. Stein, C., Diaconis, P., Holmes, S., Reinert, G.: Use of exchangeable pairs in the analysis of simulations. In P. Diaconis & S. Holmes (eds): Stein’s Method: Expository Lectures and Applications, IMS Lecture Notes, Vol. 46, 1–25 (2004)

24. Sudheesh, KK. On Stein’s identity and its applications. Stat. Probab. Lett.; 2009; 79, 12 pp. 1444-1449.MathSciNet ID: 2536504[DOI: https://dx.doi.org/10.1016/j.spl.2009.03.021]

25. Sudheesh, KK; Tibiletti, L. Moment identity for discrete random variable and its applications. Statistics; 2012; 46, 6 pp. 767-775.MathSciNet ID: 2989776[DOI: https://dx.doi.org/10.1080/02331888.2011.555548]

26. Tweedie, MCK. Statistical properties of inverse Gaussian distributions. I. Ann. Math. Stat.; 1957; 28, 2 pp. 362-377.MathSciNet ID: 110132[DOI: https://dx.doi.org/10.1214/aoms/1177706964]

27. Tweedie, MCK. Statistical properties of inverse Gaussian distributions. II. Ann. Math. Stat.; 1957; 28, 3 pp. 696-705.MathSciNet ID: 110132[DOI: https://dx.doi.org/10.1214/aoms/1177706881]

28. Wang, S; Weiß, CH. New characterizations of the (discrete) Lindley distribution and their applications. Math. Comput. Simul.; 2023; 212, pp. 310-322.MathSciNet ID: 4591343[DOI: https://dx.doi.org/10.1016/j.matcom.2023.05.003]

29. Weiß, C.H.: Control charts for Poisson counts based on the Stein–Chen identity. Advanced Statistical Methods in Statistical Process Monitoring, Finance, and Environmental Science, Springer, in press (2023)

30. Weiß, CH; Aleksandrov, B. Computing (bivariate) Poisson moments using Stein-Chen identities. Am. Stat.; 2022; 76, 1 pp. 10-15.MathSciNet ID: 4375551[DOI: https://dx.doi.org/10.1080/00031305.2020.1763836]

31. Weiß, CH; Puig, P; Aleksandrov, B. Optimal Stein-type goodness-of-fit tests for count data. Biomed. J.; 2023; 65, 2 2200073.MathSciNet ID: 4555960