The exponential distribution analog of the Grubbs--Weaver method

TL;DR

This paper establishes a new optimal estimator for the exponential distribution's scale parameter, analogous to the Grubbs-Weaver method for normal distributions, and proves the optimality of a 'rule of fours' for subsample size.

Contribution

It introduces and rigorously proves an exponential distribution analog of the Grubbs-Weaver estimator, replacing the 'rule of eights' with a 'rule of fours' for optimality.

Findings

The 'rule of fours' is optimal for the exponential distribution.

The estimator minimizes variance among unbiased estimators.

The result is rigorously proven for all sample sizes.

Abstract

Grubbs and Weaver (JASA 42 (1947) 224--241) suggest a minimum-variance unbiased estimator for the population standard deviation of a normal random variable, where a random sample is drawn and a weighted sum of the ranges of subsamples is calculated. The optimal choice involves using as many subsamples of size eight as possible. They verified their results numerically for samples of size up to 100, and conjectured that their "rule of eights" is valid for all sample sizes. Here we examine the analogous problem where the underlying distribution is exponential and find that a "rule of fours" yields optimality and prove the result rigorously.

Equations

$$X_1,X_2,\dots ,X_n\stackrel{\mathrm{iid}}{\sim }f(x)$$

$$X_{1:n}\le X_{2:n}\le \cdots \le X_{n:n}.$$

$$R_n=X_{n:n}-X_{1:n}$$

$$W_n=\frac{R_n}{\sigma }.$$

$$S=\sqrt{\frac{{\sum }_{i=1}^n(X_i-\bar{X}{)}^2}{n-1}}$$

$$c=\sqrt{\frac{2}{n-1}}\frac{\Gamma (\frac{n}{2})}{\Gamma (\frac{n-1}{2})},$$

$$\Gamma (x)={\int }_0^{\infty }{t}^{x-1}{e}^{-t}dt$$

$$\frac{R_n}{\mathrm{E}(W_n)},$$

$$\hat{\sigma }=\sum _{i=1}^m{a}_i{R}_{n_i}$$

$$(4)(3,1)(2,2)(2,1,1)(1,1,1,1).$$

$$⟨4⟩⟨13⟩⟨2^2⟩⟨1^{2}2⟩⟨1^4⟩,$$

$$P(n)\sim\frac{\pi}{12\cdot 2^{1/2}\ n^{3/2}}\ e^{\pi(2n/3)^{1/2}}\mbox{ as $n\to\infty$}$$

$$p(n)\sim\frac{1}{4\cdot 3^{1/2}\ n}\ e^{\pi(2n/3)^{1/2}}\mbox{ as $n\to\infty$}.$$

$$X_1,X_2,\dots ,X_n\stackrel{\mathrm{iid}}{\sim }\mathrm{Exp}(\theta );$$

$$f(x;\theta )=\frac{1}{\theta }{e}^{-x/\theta },$$

$$\begin{gathered} \{X_1,X_2,\dots ,X_{n_1}\},\{X_{n_1+1},\dots ,X_{n_1+n_2}\}, \\ \dots ,\{X_{n_1+n_2+\cdots +n_{m-1}+1},\dots ,X_{n_1+n_2+\cdots +n_m}\}. \end{gathered}$$

$$d_n=\frac{\mathrm{E}(R_n)}{\sigma }$$

$$k_n^2=\frac{\mathrm{E}(R_n-d_n\sigma {)}^2}{{\sigma }^2}=\frac{\mathrm{Var}(R_n)}{{\sigma }^2}.$$

$$d_n={\theta }^{-1}{H}_{n-1,1}\text{mbox}and{k}_n^2={\theta }^{-2}{H}_{n-1,2}$$

$$H_{n,j}=\sum _{i=1}^n\frac{1}{i^j}$$

$$a_i=\left(\frac{d_{n_i}}{k_{n_i}^2}\right){\left(\sum _{j=1}^m\frac{d_{n_j}^2}{k_{n_j}^2}\right)}^{-1}.$$

$$\hat{\sigma }=\sum _{i=1}^m{a}_i{R}_{n_i},$$

$$\mathrm{E}(\hat{\sigma })=\mathrm{E}\left(\sum _{i=1}^m{a}_i{R}_{n_i}\right)=\sum _{i=1}^m{a}_i\mathrm{E}(R_{n_i})=\sum _{i=1}^m{a}_i{d}_{n_i}\sigma =\sigma .$$

$$\hat{\sigma}=\frac{5880}{27133}\Big{[}\mathrm{Rng}(X_{1},\dots,X_{5})+\mathrm{Rng}(X_{6},\dots,X_{10})+\Big{]}\\ +\frac{8118}{27133}\Big{[}\mathrm{Rng}(X_{11},...,X_{14})+\mathrm{Rng}(X_{15},\dots,X_{18})+\mathrm{Rng}(X_{19},\dots,X_{22})\Big{]},$$

$$C_{j}=\frac{d_{j}^{2}}{k_{j}^{2}}=\frac{\big{(}E(R_{j})\big{)}^{2}}{{\rm Var}\,(R_{j})}=\frac{H_{j-1,1}^{2}}{H_{j-1,2}},$$

$$\mathrm{Var}(\hat{\sigma })=\mathrm{Var}\left(\sum _{i=1}^m{a}_i{R}_{n_i}\right)=\sum _{i=1}^m{a}_i^2\mathrm{Var}(R_{n_i})=\frac{{\sigma }^2}{{\sum }_{i=1}^m{\left(\frac{d_{n_i}}{k_{n_i}}\right)}^2}=\frac{{\sigma }^2}{{\sum }_{i=1}^m{C}_{n_i}}.$$

$$\text{mbox}maximize\left(C_2{f}_2+C_3{f}_3+\cdots +C_n{f}_n\right)$$

$$2f_2+3f_3+4f_4+\cdots +n{f}_n=n$$

$$f_2,f_3,f_4,\dots ,f_n\ge 0$$

$$f_2,f_3,f_4,\dots ,f_n\in \mathbb{Z}.$$

Full text

KEYWORDS: Grubbs–Weaver statistic, exponential distribution, integer partitions, combinatorial optimization

1 Introduction

Suppose

[TABLE]

is a random sample of size $n$ where each $X_{i}$ is a continuous random variable with density function $f(x)$, and (unknown) standard deviation $\sigma$. Denote the order statistics as

[TABLE]

The range of the sample is

[TABLE]

and the standardized sample range is

[TABLE]

In the case where the sample is drawn from a normal population, the minimum variance unbiased estimator of $\sigma$ is known to be the bias corrected sample standard deviation $S/c$, with

[TABLE]

and

[TABLE]

where

[TABLE]

is Euler’s gamma function. On the other hand, an easily calculuated unbiased estimator for $\sigma$ is

[TABLE]

which of course has a comparatively large variance.

Grubbs and Weaver, (1947) study a compromise between $S/c$ and $R_{n}/E(W_{n})$ in the case where the random sample is drawn from a normal population with mean $\mu$ and (unknown) variance $\sigma^{2}$. They partition the sample of size $n$ into $m$ subsamples of sizes $n_{1},n_{2},\dots,n_{m}$ respectively, where each $n_{i}\geq 2$, and $n_{1}+n_{2}+\dots+n_{m}=n$. Then estimate $\sigma$ by $\hat{\sigma}$, where

[TABLE]

with $R_{n_{i}}$ representing the range of the $i$th subsample, which is of size $n_{i}$ (see Remark 2.1 below for clarification), and $a_{1},a_{2},\dots,a_{m}$ are a set of weights chosen to guarantee that $\hat{\sigma}$ will be unbiased. They use the term “group range” to mean the range of the subsample, and thus entitle their paper “The best unbiased estimate of population standard deviation based on group ranges.”

Recall that a partition $\lambda$ of a positive integer $n$ is a representation of $n$ as an unordered sum of positive integers. Each summand is called a part of the partition, and the number of parts in a given partition is called its length. Since the order of the parts is irrelevant, it is often convenient to write the partition $\lambda$ of length $m$ as $\lambda=(n_{1},n_{2},\dots,n_{m})$ where $n_{1}\geq n_{2}\geq\dots\geq n_{m}$ and $n_{1}+n_{2}+\cdots+n_{m}=n$. Thus, e.g., the five partitions of $n=4$ are as follows:

[TABLE]

Alternatively, we may define $f_{i}=f_{i}(\lambda)$ to denote the frequency (or multiplicity) of part $i$ in the partition $\lambda$, i.e. the number of times the part $i$ appears, and employ the “frequency superscript notation” $\langle 1^{f_{1}}2^{f_{2}}3^{f_{3}}\cdots\rangle$. In this notation, the five partitions of $n=4$ are

[TABLE]

where we have followed the convention that $f_{i}=1$ the superscript is omitted and $i^{f_{i}}$ is omitted if $f_{i}=0$.

We are interested in partitions of $n$ where all parts are at least $2$, as a subsample of size $1$ will, by definition, have a range of [math]. Let us call such a partition admissible. If $p(n)$ denotes the number of unrestricted partitions of $n$, it is well-known and easy to prove that the number $P(n)$ of admissible partitions of $n$ is equal to $p(n)-p(n-1)$. The number of admissible partitions of $n$ increases rather rapidly with $n$; e.g., $P(100)=21{,}339{,}417$. To get a rough idea of the size of $P(n)$ for general $n$, we mention in passing that the asymptotic formula

[TABLE]

may be deduced from a theorem of Meinardus, (1954). Equation (1.2) is analogous to the famous asymptotic formula of Hardy and Ramanujan, (1918, p. 79, Eq. 1.41) for the unrestricted partition function

[TABLE]

Associated with each admissible partition $\lambda$ of $n$ is an estimator $\hat{\sigma}=\hat{\sigma}_{\lambda}$ defined in (1.1); for a given $n$, the admissible partition $\lambda$ that corresponds to the $\hat{\sigma}_{\lambda}$ of minimum variance will be called the optimal partition of $n$.

Grubbs and Weaver, (1947) performed extensive computations for $2\leq n\leq 100$, and showed that the optimal partition of $n$, when the underlying distribution is normal, uses as many $8$’s as possible (with occasional $7$’s or $9$’s to adjust for the fact that not every $n$ is a multiple of $8$), except for sporadic exceptions that occur for small $n$.  Grubbs and Weaver, (1947) did not supply a rigorous proof that their assertions held for $n>100$, perhaps owing to the lack of closed form expressions for the expected value and variance of the range when the underlying distribution is normal.

Here we investigate the analogous problem in the case where the underlying distribution is exponential.

2 Precise Statement of the Problem

Let

[TABLE]

i.e. each $X_{i}$ has pdf

[TABLE]

for $x\geq 0$ and [math] otherwise, where $\theta>0$. Since ${\rm Var}\,(X_{i})=\theta^{2}$, by estimating the population standard deviation $\sigma$, we are equivalently estimating the parameter $\theta$.

Remark 2.1**.**

For a given partition $\lambda=(n_{1},n_{2},\dots,n_{m})$ of $n$ with length $m$, we understand that the $m$ subsamples of the random variables $X_{1},X_{2},\dots,X_{n}$ are to be

[TABLE]

For each admissible partition of $n$, we have $f_{1}=0$ since each part must be at least $2$, and $f_{n-1}=0$ because a partition of $n$ with no $1$’s clearly cannot have $n-1$ as a part.

Following the notation in Grubbs and Weaver, (1947), let

[TABLE]

and

[TABLE]

It is well-known (see, e.g., David and Nagaraja, (2003, p. 52, Ex. 3.2.1)) that

[TABLE]

where

[TABLE]

are generalized harmonic numbers.

A linear combination of random variables which gives the minimum variance unbiased estimate has coefficients which are inversely proportional to the variances of the variables. Consequently, let

[TABLE]

Proposition 2.2**.**

The estimator

[TABLE]

where $a_{i}$ is defined in (2.1), is an unbiased estimator of $\sigma$.

Proof.

[TABLE]

∎

Now we state the main theorem:

Theorem 2.3** (The Exponential “Rule of Fours”).**

Fix an integer $n>1$. Write $n=4q+r$, where $r=r(n)=\mathbf{mod}(n,4)$, the least nonnegative residue of $n$ modulo $4$. The optimal partition of $n$ is as follows:

• •

$\langle 4^{q}\rangle$, if $r=0$;

• •

$\langle 4^{q-r}5^{r}\rangle$, if $r=1$ or $2$ and $q\geq r$;

• •

$\langle 3^{1}4^{q}\rangle$, if $r=3$ and $q\geq 1$;

• •

$\langle n\rangle=(n)$, if $2\leq n\leq 5$; and

• •

$\langle 3^{2}\rangle=(3,3)$, if $n=6$.

Example 2.4**.**

Suppose we have the sample size $n=22$ which implies $q=5$ and $r=2$. Thus, according to Theorem 2.3, the optimal partition of $22$ is $(5,5,4,4,4)$ and this in turn gives

[TABLE]

where for $i<j$, $\mathrm{Rng}(X_{i},\dots,X_{j})$ denotes the range of the subsample $X_{i},X_{i+1},\dots,X_{j}$ of the sample $X_{1},X_{2},\dots,X_{n}$.

3 Proof of Theorem 2.3

Next, Theorem 2.3 will be reformulated into an integer linear program. Upon solving the equivalent integer linear program, we will have proved Theorem 2.3.

3.1 Reformulation as an integer linear program

For any given positive integer $j$, let us define

[TABLE]

observing that $C_{j}$ is independent of $\sigma$. Thus, the $C_{j}$ are known absolute constants. For a given random sample of size $n$, we seek a partition $\lambda$ of $n$ such that ${\rm Var}\,\hat{\sigma}$ is minimized when the sample is partitioned according to $\lambda$.

Notice that

[TABLE]

Thus we seek the partition of $n$ which causes the denominator in the rightmost expression of (3.1) to be as large as possible.

For any given $n$, an optimal partition is given by a solution of the following integer linear program:

[TABLE]

subject to

[TABLE]

[TABLE]

[TABLE]

The objective function is the denominator in the far right member of (3.1). The constraints guarantee that we search only over admissible partitions of $n$.

The program (EKP) is an instance of the equality-constrained “knapsack problem.” While integer linear programming problems are in general NP-complete, we will be able to exploit the special structure to provide an optimal solution to (EKP).

3.2 Statement and proof of the key lemma

But first, our solution is contingent upon the following lemma.

Lemma 3.1**.**

The function $C_{n}/n$ attains its maximum value at $n=4$.

Remark 3.2**.**

The significance of maximizing function $C_{n}/n$ is as follows: given a sample of size $N$ that is partitioned into $N/n$ samples of size $n$ each, the denominator in the rightmost member of (3.1) is $NC_{n}/n$. The sample size $N$ is fixed, but we wish to find the part size $n$ that maximizes $C_{n}/n$ in order to maximize the denominator of the rightmost member of (3.1). Of course, this can only be done exactly in the case where the sample size $N$ is a multiple of the part size $n$, accounting for the fact that while $n=4$ gives the optimal part size, for sample sizes that are not multiples of $4$, the optimal partition includes some parts of size $3$ or $5$.

Proof of Lemma 3.1.

We need to show that for $n\neq 4$,

[TABLE]

We will utilize the elementary inequalities

[TABLE]

and

[TABLE]

Inequality (3.2) may be deduced as follows: note that $\int_{x}^{x+1}1/t\ dt=\log(x+1)-\log x$ for $x>0$. Then observe that $1/x>1/t>1/(x+1)$ in the integrand, thus $1/x>\log(x+1)-\log(x)>1/(x+1)$ and thus $\log(n+1)<H_{n,1}<1+\log n$ follows by summation. Inequality (3.3) follows from the fact that

[TABLE]

Applying (3.2) and (3.3), we obtain

[TABLE]

Let $h(n)={(1+\log(n-1))^{2}}/(n-1)$. Thus, $C_{n}/n$ is bounded above by $h(n)$. Let us now temporarily consider $n$ to be a real variable with domain $n>1$. Notice that

[TABLE]

Since $1+\log(n-1)>0$ and $(n-1)^{2}>0$ for $n>1$, $dh/dn$ is negative whenever $\log(n-1)>1$, i.e. when $n>e+1$. Thus $h(n)$ is decreasing for all integers $n\geq 4$. Now $h(34)<121/196$. Thus $C_{n}/n<h(n)<121/196$ for all $n\geq 34$. That $C_{n}/n<121/196$ for $n=2,3$ and $5\leq n\leq 33$ can be verified by direct computation. ∎

3.3 Solving (EKP) by the method of group relaxation

With Lemma 3.1 in hand, we now proceed to solve (EKP). Fix $n=4q+r$, with $r=\mathbf{mod}(n,4)$. We now closely follow the treatment of Lee, (2004, §7.2, p. 184 ff) to solve (EKP).

In general, the first step is to seek an upper bound for

[TABLE]

and then, if necessary, initiate a branch-and-bound procedure. In the case of interest, this initial step will be sufficient to find an optimal solution to (EKP). By Lemma 3.1 we have that

[TABLE]

Relax the nonnegativity restriction on $f_{4}$, and solve for $f_{4}$ in terms of the other $f_{j}$’s, namely

[TABLE]

to obtain the “group relaxation”:

[TABLE]

subject to

[TABLE]

[TABLE]

[TABLE]

The group relaxation (GR) is in turn equivalent to:

[TABLE]

subject to

[TABLE]

[TABLE]

[TABLE]

We note that every feasible solution in (GR*′) corresponds to a value of $f_{4}$ in (GR), via Equation (3.6). If we are lucky, the optimal solution we find to (GR′*) yields an $f_{4}\geq 0$, and then we will have also found an optimal solution to (EKP). (If we are not lucky and the corresponding $f_{4}$ in (GR) is negative, then we must embark on a branch-and-bound procedure with potentially many iterations.) Here, however, we will show that for all $n>6$, we are indeed lucky, and the optimal solution to (EKP) will indeed be found directly with no need to initiate branch-and-bound.

In order to solve (GR*′*), we form a weighted directed multigraph $G$ as follows. Let the vertex set be given by $V(G)=\{0,1,2,3\}$. For each $v\in V(G)$ and $j\in\{2,3,5,6,7,8,\dots,n\}$, there is an edge $e(v,j)$ of weight $-C_{j}+jC_{4}/4$ from $v$ to $\mathbf{mod}(v+j,4)$. Note well that in this notation, the “ending vertex” of edge $e(v,j)$ is not $j$ but rather $\mathbf{mod}(v+j,4)$. We seek a minimum weight directed walk from [math] to $r$. Each time we include the edge $e(v,j)$ in the diwalk, we increment $f_{j}$ by $1$. Since for each $j\neq 4$, the edge weight $-C_{j}+jC_{4}/4\geq 0$, we may use the algorithm of Dijkstra, (1959) to find a minimum weight dipath from [math] to $r$.

We first dispense with the trivial case $r=0$. Here we immediately have the optimal solution (via the empty path from [math] to [math])

[TABLE]

We now move on to the remaining cases $1\leq r\leq 3$. Notice that connecting any two vertices in $G$ there are multiple directed edges. For each $j\in\{2,3,4,\dots,n\}\setminus\{4\}$, connecting vertex $v$ to vertex $\mathbf{mod}(v+j,4)$, we have the following edges: $e(v,j),e(v,2j),e(v,3j),\dots,e(v,mj)$, where $m$ is the largest integer such that $v+mj\leq n$. All of these edges have different weights, and since we seek a minimum weight diwalk, for each $j$, we may safely remove all but the one of lowest weight, resulting in a much less “cluttered” digraph $G^{\prime}$. Thus $G^{\prime}$ is a digraph with four vertices and $12$ directed edges (as each of the $4$ vertices now has exactly one directed edge to each of the $3$ other vertices).

[TABLE]

Perform the algorithm of Dijkstra, (1959) on the weighted directed graph $G^{\prime}$ to find that the minimum weight directed paths are as follows:

from $\rightarrow$ to using edges total weight resulting nonzero values

$0\to 1$ $e(0,1)$ $305/8036$ $f_{5}=1,f_{4}=q-1$

$0\to 2$ $e(0,1)$ and $e(1,1)$ $610/8036$ $f_{5}=2,f_{4}=q-2$

$0\to 3$ $e(0,3)$ $51/980$ $f_{3}=1,f_{4}=q$

Thus we arrive at precisely the desired result.

Since the value of $f_{4}$ obtained from the group relaxation (GR) (see Eq. (3.6)) is always nonnegative, then group relaxation solves the original equality knapsack problem (EKP), and thus no branch-and-bound procedure need be undertaken.

With the solution of (EKP), we have proved Theorem 2.3 for $n>6$. For $2\leq n\leq 6$, the theorem can be established by direct calculation. ∎

4 Conclusion

We have proved that an unbiased estimator from an exponential population formed by taking an appropriately weighted sum of the ranges of subsamples obtained by partitioning the original sample of size $n$ has minimal variance when the subsamples are each of size $4$ (or as close as possible when $n$ is not a multiple of $4$). This contrasts with the work of Grubbs and Weaver, (1947), where the optimal subsample sizes are $8$, when the population is normal. A similar analysis could be applied to estimate the standard deviation from populations with other distributions. For example, based on some preliminary calculations, in the case of the Rayleigh distribution, the optimal subsample sizes appear to be $8$, and in the case of the $\chi^{2}$ distribution, the optimal subsample sizes appear to be $3$. Another variant perhaps worth investigating is seeking the optimal partition of the sample where instead of the range of subsamples, we consider various quasi-ranges.

Harter and Balakrishnan, (1996, pp. 9–11) compare the efficiency of estimators of $\sigma$ based on quasi-ranges with the Grubbs–Weaver estimator for $\sigma$ in the case of a normal parent population and remark [p. 11] that the asymptotic efficiency of the Grubbs–Weaver estimator is 75.38 percent. In future work we would like to obtain analogous results for various non-normal parent distributions, including the exponential distribution that was studied in this present paper.

Acknowledgments

The authors thank the Editor and the anonymous referee for helpful suggestions that improved the paper.