Robust estimation of the mean with bounded relative standard deviation

TL;DR

This paper improves the efficiency of randomized algorithms estimating a mean with bounded relative standard deviation by optimizing the sampling process and analyzing practical implementation, supported by numerical experiments.

Contribution

It introduces a modified algorithm that balances tails to reduce sample complexity and discusses practical implementation strategies.

Findings

Improved sample complexity bounds for mean estimation with bounded relative standard deviation.

Demonstrated effectiveness of the new algorithm through numerical experiments.

Provided insights into practical implementation when the relative standard deviation is unknown or infinite.

Abstract

Many randomized approximation algorithms operate by giving a procedure for simulating a random variable $X$ which has mean $\mu$ equal to the target answer, and a relative standard deviation bounded above by a known constant $c$. Examples of this type of algorithm includes methods for approximating the number of satisfying assignments to 2-SAT or DNF, the volume of a convex body, and the partition function of a Gibbs distribution. Because the answer is usually exponentially large in the problem input size, it is typical to require an estimate $\hat \mu$ satisfy $\mathbb{P}(|\hat \mu/\mu - 1| > \epsilon) \leq \delta$, where $\epsilon$ and $\delta$ are user specified nonnegative parameters. The current best algorithm uses $\lceil 2c^2\epsilon^{-2}(1+\epsilon)^2 \ln(2/\delta) \rceil$ samples to achieve such an estimate. By modifying the algorithm in order to balance the tails, it is possible to improve this result to $\lceil 2(c^2\epsilon^{-2} + 1)/(1-\epsilon^2)\ln(2/\delta) \rceil$ samples. Aside from the theoretical improvement, we also consider how to best implement this algorithm in practice. Numerical experiments show the behavior of the estimator on distributions where the relative standard deviation is unknown or infinite.

Tables

Table 1. Table 1: Behavior of the sample average, sample mean, and M 𝑀 M -estimator for 100 exponential random variables with mean 2. here. Repeated five times to show variation.

Mean	Median	$\lambda =0.1$	$\lambda =1$	$\lambda =5$
\svhline 2.34	1.86	2.33	1.93	1.87
1.89	1.35	1.88	1.53	1.40
2.29	1.78	2.28	1.83	1.77
2.02	1.37	2.01	1.48	1.35
2.17	1.37	2.16	1.70	1.39

Table 2. Table 2: Behavior of the sample mean, sample median, and M 𝑀 M -estimator for 100 draws from the absolute value of a Cauchy distribution. Repeated five times to show variation.

Mean	Median	$\lambda =0.1$	$\lambda =1$	$\lambda =5$
\svhline 2.30	0.70	1.95	0.88	0.69
2.70	1.10	2.56	1.27	1.11
10.29	1.01	2.96	1.24	1.02
3.59	1.09	2.58	1.32	1.13
8.07	1.25	4.83	1.68	1.34

Table 3. Table 3: Behavior of the sample mean, sample median, and M 𝑀 M -estimator for 100 draws from the absolute value of a Cauchy distribution. Repeated five times to show variation.

epsilon	delta	CG	New	relative change
\svhline 0.10	1e-06	551.94	461.48	-0.1638946
0.05	1e-06	2013.24	1822.40	-0.0947925

Mean	Median	$\lambda$ = 0.1	$\lambda$ = 1	$\lambda$ = 5
2.34	1.86	2.33	1.93	1.87
1.89	1.35	1.88	1.53	1.40
2.29	1.78	2.28	1.83	1.77
2.02	1.37	2.01	1.48	1.35
2.17	1.37	2.16	1.70	1.39

Mean	Median	$\lambda$ = 0.1	$\lambda$ = 1	$\lambda$ = 5
2.30	0.70	1.95	0.88	0.69
2.70	1.10	2.56	1.27	1.11
10.29	1.01	2.96	1.24	1.02
3.59	1.09	2.58	1.32	1.13
8.07	1.25	4.83	1.68	1.34

epsilon	delta	CG	New	improvement
0.10	1e-06	551.94	461.48	-0.1638946
0.05	1e-06	2013.24	1822.40	-0.0947925

Equations

$$\mathbb{P}\left(\left|\frac{\hat{\mu }}{\mu }-1\right|>ϵ\right)\le \delta .$$

$$S_n=\frac{X_1+\cdots +X_n}{n}.$$

$$\mathbb{P}(|S_n-\mu |>ϵ\mu )\le \frac{\mathrm{Var}(X)}{n{ϵ}^2{\mu }^2}\le c^2{ϵ}^{-2}{n}^{-1}.$$

$$\mathbb{P}(|S_n-\mu |>ϵ\mu )=\mathbb{P}(|Z|\ge ϵ\sqrt{n}/c),$$

$$n=2c^2{ϵ}^{-2}[\ln (2/\delta )-o(\delta )].$$

$$\mathbb{P}(|\hat{\theta }-\mu |>ϵ\mu )\ge (1/2)\mathbb{P}(|Z|\ge ϵ\sqrt{n}/c).$$

$$n\ge 2c^2{ϵ}^{-2}[\ln (1/\delta )-o(\delta )].$$

$$n=\lceil 2c^2{ϵ}^{-2}(1+ϵ{)}^2\ln (2/\delta )\rceil$$

$$\lceil 2(c^2{ϵ}^{-2}+1)(1-{ϵ}^2{)}^{-1}\ln (2/\delta )\rceil$$

$$\frac{1+{ϵ}^2/c^2}{(1+ϵ{)}^2(1-{ϵ}^2)}.$$

$$\Psi (m)=\sum _{i=1}^n\psi (x_i,m).$$

$$f(x,m)=\left\{\begin{array}[]{ll}(x/m)-1&m\neq 0\\ 0&x=m=0.\end{array}\right.$$

$$\sum _{i=1}^{m}f(x_i,m)=0.$$

$$\sum _{i=1}^{n}f(x_i,m)=0\iff m=\frac{{\sum }_{i=1}^n{x}_i}{n}.$$

$$g(x,m)=\text{mathbbold}1\left(f(x,m)>0\right)-\text{mathbbold}1\left(f(x,m)<0\right).$$

$$\sum _{i=1}^{n}g(x_i,m)=0.$$

$$d_f(u)$$

$$d_g(u)$$

$$d_h(u)$$

$$h(x_i,m)=d_h(u(x_i,m)).$$

$$d_L(u)$$

$$d_U(u)$$

$$d_L(u)\le d_h(u)\le d_U(u)$$

$$h_{\lambda }(x_i,m)={\lambda }^{-1}{d}_h(\lambda (x_i/m-1)).$$

$$n=\lceil 2c^2{ϵ}^{-2}\ln (2/\delta )(1+ϵ{)}^2\rceil .$$

$${\Psi }_{\lambda }(m)=\frac{1}{n}\sum _{i=1}^n{\lambda }^{-1}{d}_h(\lambda u(X_i,m)).$$

$$\mathbb{P}({\Psi }_{\lambda }((1+ϵ)\mu )\ge 0)$$

$$\mathbb{P}({\Psi }_{\lambda }((1-ϵ)\mu )\le 0)$$

$$u_{ϵ}=\frac{X}{(1+ϵ)\mu }-1.$$

$$\mathbb{P}({\Psi }_{\lambda }((1+ϵ)\mu )\ge 0)$$

$$\le \mathbb{E}[\exp (\lambda {\Psi }_{\lambda }((1+ϵ)\mu ))]$$

$$\exp (\lambda {\Psi }_{\lambda }(m))=\prod _{i=1}^n\exp (\lambda d_h(X_i/m-1)).$$

Full text

1 Introduction

Suppose we are interested in approximating a target value $\mu$. Then many randomized approximation algorithms work by constructing a random variable $X$ such that ${\mathbb{E}}[X]=\mu$ and $\operatorname{Var}(X)/\mu^{2}\leq c^{2}$ for a known constant $c$. The randomized algorithm then simulates $X_{1},X_{2},\ldots,X_{n}$ as independent identically distributed (iid) draws from $X$. Finally, the values are input into a function to give an estimate $\hat{\mu}$ for $\mu$.

Examples of this type of algorithm include when $\mu$ is the number of solutions to a logic formula in Disjunctive Normal Form (DNF) karpl1983 , the volume of a convex body dyerfk1991 , and the partition function of a Gibbs distribution huber2015a . For all of these problems, the random variable $X$ used is nonnegative with probability 1, and so we shall only consider this case for the rest of the paper.

For these types of problems, generating the samples from a high dimensional distribution is the most computationally intensive part of the algorithm. Hence we will measure the running time of the algorithm by the number of samples from $X$ that must be generated.

The answer $\mu$ for these problems typically grows exponentially quickly in the size of the problem input size. Therefore, it is usual to desire an approximation $\hat{\mu}$ that is accurate when measured by relative error. We use $\epsilon>0$ as our bound on the relative error, and $\delta>0$ as the bound on the probability that the relative error restriction is violated.

Definition 1

Say $\hat{\mu}$ is an $(\epsilon,\delta)$-randomized approximation scheme ($(\epsilon,\delta)$-ras) if

[TABLE]

This is equivalent to considering a loss function $L$ that is $L(\hat{\mu},\mu)={\mathbbold{1}}(|(\hat{\mu}/\mu)-1|>\epsilon)$, and requiring that the expected loss be no more than $\delta$.

Note that this question is slightly different than the problem most statistical estimates are designed to handle. For instance, the classic work of huber1964 is trying to minimize the asymptotic variance of the estimator, not determine how often the relative error is at most $\epsilon$.

The first approach to this problem is often the standard sample average

[TABLE]

For this estimator, only knowing that $\operatorname{Var}(X)\leq c^{2}\mu^{2}$, the best we can say about the probability that the relative error is large comes from Chebyshev’s inequality tchebichef1867 .

[TABLE]

In particular, setting $n=c^{2}\epsilon^{-2}\delta^{-1}$ gives an $(\epsilon,\delta)$-ras. There are simple examples where Chebyshev’s inequality is tight.

While this bound works for all distributions, for most distributions in practice the probability of error will go down exponentially (and not just polynomially) in $n$. We desire an estimate that matches this speed of convergence.

For example, suppose that $X$ is a normal random variable with mean $\mu$ and variance $c^{2}\mu^{2}$ (write $X\sim\textsf{N}(\mu,c^{2}\mu^{2})$. Then the sample average is normally distributed as well. To be precise, $S_{n}\sim\textsf{N}(\mu,c^{2}\mu^{2}/n)$, and it is straightforward to show that

[TABLE]

where $Z$ is a standard normal random variable. This gives

[TABLE]

samples being necessary and sufficient to achieve an $(\epsilon,\delta)$-ras.

We did this calculation for normally distributed samples, but in fact this gives a lower bound on the number of samples needed. For any $X_{1},\ldots,X_{n}$ with mean $\mu>0$ and variance $c^{2}>0$, and an estimate $\hat{\theta}=\hat{\theta}(X_{1},\ldots,X_{n})$,

[TABLE]

See, for instance, Proposition 6.1 of catoni2012 . Hence the minimum number of samples required for all instances is

[TABLE]

The method presented in hubertoappearb comes close to this lower bound, requiring

[TABLE]

samples. So it is larger than optimal by a factor of $(1+\epsilon)^{2}$.

Our main result is to modify the estimator slightly. This has two beneficial effects.

The nuisance factor is reduced from first order in $\epsilon$ to second order. To be precise, the new nuisance factor is $(1+\epsilon^{2}/c^{2})/(1-\epsilon^{2})$. 2. 2.

It is possible to solve exactly for the value of the $M$-estimator (using square and cube roots) rather than through numerical appoximation if so desired.

Theorem 1.1

Suppose that ${\mathbb{P}}(X\geq 0)=1$, ${\mathbb{E}}[X]=\mu$, and the standard deviation of $X$ is at most $c\mu$, where $c>0$. Then there exists an $(\epsilon,\delta)$-ras $\hat{\mu}_{1}$ where at most

[TABLE]

draws from $X$ are used.

Ignoring the ceiling function, the new method uses a number of samples bounded by the old method times a factor of

[TABLE]

For instance, when $\epsilon=0.1$ and $c=2$, this factor is $0.8556\ldots$, and so we obtain an improvement of over 14% in the running time. At first, this might seem slight, but remember that the best improvement we can hope to make based on normal random variables is $1/(1+\epsilon)^{2}=0.8264$, or a bit less than 18%. Therefore, this does not quite obtain the maximum improvement of $1/(1+2\epsilon+\epsilon^{2})$, but does obtain a factor of $(1+O(\epsilon^{2}))/(1+2\epsilon+\epsilon^{2})$.

The remainder of this paper is organized as follows. In Section 2, we describe what $M$-estimators are, and show how to down weight samples that are far away from the mean. Section 3 shows how to find the estimator both approximately and exactly. In Section 4 we consider using these estimators on some small examples and see how they behave numerically.

2 $\Psi$-estimators

The median is a robust centrality estimate, but it is easier to build a random variable whose mean is the target value. P.J. Huber first showed in huber1964 how to build an estimate that has the robust nature of the median estimate while still converging to the mean.

Consider a set of real numbers $x_{1},\ldots,x_{n}$. The sample average of the $\{x_{i}\}$ is the point $m$ where the sum of the distances to points to the right of $m$ equals the sum of the distances to points to the left of $m$.

This idea can be generalized as follows. First, begin with a function $\psi:{\mathbb{R}}\times{\mathbb{R}}^{n}\rightarrow{\mathbb{R}}$. Second, using $\psi$, form the function

[TABLE]

For convenience, we suppress the dependence of $\Psi$ on $(x_{1},\ldots,x_{n})$ in the notation.

Consider the set of zeros of $\Psi$. These zeros form the set of $\Psi$-type $M$-estimators for the center of the points $(x_{1},\ldots,x_{n})$.

For example, suppose our $\psi$ is $f:{\mathbb{R}}^{2}\rightarrow{\mathbb{R}}$ defined as

[TABLE]

Then consider

[TABLE]

For $x$ and $m$ with the same units, the right hand side is unitless. If the $x_{i}$ are nonnegative $\{x_{i}\}$ and at least one $x_{i}$ is positive,

[TABLE]

That is, the unique $M$-estimator using $f$ as our $\psi$ function is the sample average.

Now consider

[TABLE]

Then we wish to find $m$ such that

[TABLE]

Summing $g(x_{i},m)$ adds 1 when $x_{i}$ greater than $m$, and subtracts 1 when $x_{i}$ is $m$. Then the sum of the $g(x_{i},m)$ is zero exactly when there are an equal number of $x_{i}$ that are above and below $m$. When the number of distinct $x_{i}$ values is odd, then $m$ has a unique solution equal to the sample median.

The sample median has the advantage of being robust to large changes in the $x_{i}$ values, but this will converge to the median of the distribution (under mild conditions) rather than the mean. The sample average actually converges to the mean when applied to iid draws from $X$, and this is the value we care about finding. However, the sample average can be badly thrown off by a single large outlier. Our goal is to create a $\Psi$ function that combines the good qualities of each while avoiding the bad qualities.

Both $f$ and $g$ can be expressed using weighted differences. Let $d_{f}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ and $d_{g}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ be defined as

[TABLE]

Let $u(x_{i},m)=x_{i}/m-1$ for $m\neq 0$, and [math] for $x_{i}=m=0$. Then $f(x_{i},m)=d_{f}(u)$ and $g(x_{i},m)=d_{g}(u)$.

Catoni catoni2012 improved upon this $\Psi$-estimator by using a function that approximated $d_{f}$ for $|u|\leq 1$, and approximated $d_{u}$ for $|u|>1$, and could be analyzed using the Chernoff bound approach chernoff1952 . Catoni and Guillini then created an easier to use version of the $\Psi$-estimator in catonig2017 .

The following is a modification of the Catoni and Guillini $\Psi$-estimator. Unlike catonig2017 , the weighted difference here does not have square roots in the constants, which makes them slightly easier to work with from a computational standpoint.

[TABLE]

Then define

[TABLE]

The function $d_{h}$ behaves like the mean weight for values near 0, but like the median weights for values far away from 0. See Figure 1. In order to link this function to the first and second moments of the random variable, the following links $d_{h}$ to weighted differences introduced in catoni2012 . See Figure 2.

Lemma 1

Let

[TABLE]

Then for all $u\in{\mathbb{R}}$,

[TABLE]

Proof

Note $d^{\prime}_{U}(u)=(1+u)/(1+u+u^{2}/2)$. This derivative is positive over $[1,\infty)$ so $d_{U}$ is increasing in this region. Since $d_{U}(1)=\ln(2.5)>5/6=d_{h}(1)$, $d_{U}(u)\geq d_{h}(u)$ over $[1,\infty)$.

Similarly, $d^{\prime}_{U}(u)<0$ for $u\in(-\infty,-1]$, and so it is decreasing in this region, and $d_{U}(-1)=\ln(1/2)\geq-5/6=d_{h}(u)$ so $d_{U}(u)\geq d_{h}(u)$ in $(-\infty,1]$. Finally, inside $[-1,1]$, $d_{h}(u)=u-u^{3}/6$. So a minimum of $d_{U}(u)-d_{h}(u)$ occurs either at $-1$, $1$, or a critical point where $d^{\prime}_{U}(u)-(1-u^{2}/2)=0$. The unique critical point in $[-1,1]$ is at $u=0$, and $d_{U}(0)=d_{h}(0)$. This means $d_{U}(u)-d_{h}(u)\geq 0$ for all $u\in[-1,1]$.

The lower bound follows from $d_{L}(u)=-d_{U}(-u)\leq-d_{h}(-u)=d_{h}(u)$.

It helps to introduce a scale factor $\lambda$ that allows us to extend the effective range where $d(u)\approx u$. Let

[TABLE]

Then $\lambda>0$ is a parameter that can be chosen by the user ahead of running the algorithm based upon $\epsilon$ and $c$.

The following was shown as Lemma 13 of hubertoappearb .

Lemma 2

For given $\epsilon,\delta>0$, let $\epsilon^{\prime}=\epsilon/(1+\epsilon)$ and $\lambda=\epsilon^{\prime}/c^{2}$. Set

[TABLE]

For $X_{1},\ldots,X_{n}$ iid $X$, form the function

[TABLE]

Let $\hat{\mu}$ be any value of $m$ such that $|\Psi_{\lambda}(m)|\leq(\epsilon^{\prime})^{2}/2.$ Then $\hat{\mu}$ is an $(\epsilon,\delta)$-ras for $\mu$.

To improve upon this result and obtain Theorem 1.1, we must be more careful about our choice of $\lambda$.

To meet the requirement of an $(\epsilon,\delta)$-ras, we must show that $\Psi_{\lambda}(m)$ has all its zeros in the interval $((1-\epsilon)\mu,(1+\epsilon)\mu)$ with probability at least $1-\delta$. Given that $\Psi_{\lambda}$ is continuous and decreasing, it suffices to show that $\Psi_{\lambda}((1-\epsilon)\mu)>0$ and $\Psi_{\lambda}((1+\epsilon)\mu)<0$.

Lemma 3

Let $X_{1},\ldots,X_{n}$ be iid $X$. Then for $\lambda>0$,

[TABLE]

where

[TABLE]

Proof

Note that

[TABLE]

by Markov’s inequality. Note

[TABLE]

Each term in the product is independent, therefore the mean of the product is the product of the means (which are identical.)

Let $u=\lambda(X/m-1)$. Then

[TABLE]

From the previous lemma,

[TABLE]

Hence

[TABLE]

Putting $m=(1+\epsilon)\mu$ into this expression then gives the first inequality.

For the second inequality, the steps are nearly identical, but we begin by multiplying by $-\lambda$. This completes the proof.

In particular, if we choose $n$ so that ${\mathbb{P}}(\Psi_{\lambda}((1+\epsilon)\mu)>0)\leq\delta/2$, and ${\mathbb{P}}(\Psi_{\lambda}((1-\epsilon)\mu)<0)\leq\delta/2$, then by the union bound the probability that $\Psi$ has a root in $[(1-\epsilon)\mu,(1+\epsilon)\mu]$ is at least $\delta/2+\delta/2=\delta$.

As is well known, for $g>0$,

[TABLE]

and so if $n\geq(1/g)\ln(2/\delta)$,

[TABLE]

We refer to $g$ as the gap. Since the number of samples needed is inversely proportional to the gap, we wish the gap to be as large as possible. We can lower bound the gap given by the previous lemma in terms of $\lambda$ and $m$.

Lemma 4

For $u=\lambda(X/m-1)$, and $a(m)=1-(m/\mu)$,

[TABLE]

Similarly,

[TABLE]

Proof

From linearity of expectations,

[TABLE]

The second moment is the sum of the variance plus the square of the first moment. Using $\operatorname{Var}(c_{1}X+c_{2})=c_{1}^{2}\operatorname{Var}(X)$, we get

[TABLE]

Since $\operatorname{Var}(X)\leq c^{2}\mu^{2}$,

[TABLE]

giving the first result. The proof of the second statement is similar.

Lemma 5

Let $p_{\ell}={\mathbb{P}}(\Psi_{\lambda}((1-\epsilon)\mu)\leq 0)$ and $p_{r}={\mathbb{P}}(\Psi_{\lambda}((1+\epsilon)\mu)\geq 0)$. Then

[TABLE]

Proof

Note

[TABLE]

Combine Lemmas 3, 4, and 6 with $m$ equal to $(1-\epsilon)\mu$ to get the first inequality. Then set $m=(1+\epsilon)\mu$ and use $a((1+\epsilon)\mu)=-\epsilon$ to get the second.

Our goal is to use as few samples as possible, which means simultaneously minimizing the quantities inside the brackets in Lemma 5. Both of the upper bounds are upward facing parabolas, and they have a unique minimum value.

Lemma 6

The minimum value of $f(\lambda)=1+a_{1}\lambda+(1/2)a_{2}\lambda^{2}$ is

[TABLE]

at $\lambda^{*}=-a_{1}/a_{2}.$

Proof

Complete the square.

However, because the coefficients in the quadratic upper bounds are different, it is not possible to simultaneously minimize these bounds with the same choice of $\lambda$.

What can be said is that these bounds are quadratic with positive coefficient on $\lambda^{2}$, and so the best we can do is to choose a $\lambda$ such that the upper bounds are equal to one another. This gives us our choice of $\lambda$.

Lemma 7

Let

[TABLE]

Then

[TABLE]

Proof

Follows directly from Lemma 5.

This makes the inverse gap

[TABLE]

and immediately gives Theorem 1.1.

3 Computation

Set $\lambda=\epsilon(1-\epsilon^{2})(c^{2}+\epsilon^{2})^{-1}$. Consider how to locate any root of $\Psi_{\lambda}$ for a given set of $X_{1},\ldots,X_{n}$. As before, we assume that the $X_{i}$ are nonnegative and not all identically zero. The function $\Psi_{\lambda}(m)$ is continuous and decreasing, although not necessarily strictly decreasing. Therefore, it might have a set of zeros that form a closed interval.

3.1 An approximate method

Suppose that the points $X_{i}$ are sorted into their order statistics,

[TABLE]

Then $\Psi_{\lambda}(X_{(1)})>0$ and $\Psi_{\lambda}(X_{(n)})<0$, so there exists some $i$ such that $\Psi_{\lambda}(X_{(i)})\geq 0$ and $\Psi_{\lambda}(X_{i+1})\leq 0$. Since any particular value of $\Psi_{\lambda}$ requires $O(n)$ time to compute, this index $i$ can be found using binary search in $O(n\ln(n))$ time.

At this point we can switch from a discrete binary search over $\{1,\ldots,n\}$ to a continuous binary search over the interval $[X_{i},X_{i+1}]$. This allows us to quickly find a root to any desired degree of accuracy.

This is the method that would most likely be used in practice.

3.2 An exact method

Although the approximation procedure is what would be used in practice because of its speed, there does exist a polynomial time exact method for this problem. For $i<j$, and a value $m$, suppose that the points of $\{X_{i}\}$ that fall into the interval $[m-\lambda,m+\lambda]$ are exactly $\{X_{(i)},X_{(i+1)},\ldots,X_{(j)}\}.$ Then say that $m\in m(i,j)$.

That is, define the set $m(i,j)$ for $i<j$ as follows.

[TABLE]

See Figure 3.

Note that

[TABLE]

There are at most $n$ choose 2 such $i<j$ where $m(i,j)$ is nonempty. In fact, since as we slide $m$ from $X_{(1)}$ to $X_{(n)}$, each point can enter or leave the interval $[m-\lambda,m+\lambda]$ exactly once. Therefore there are at most $2n$ pairs $(i,j)$ where $m(i,j)$ is nonempty.

That means that to find a root of $\Psi_{\lambda}$ we need merely check if there is a root $m\in m(i,j)$ for all $i<j$ such that $m(i,j)\neq\emptyset$.

For each $m\in m(i,j)$, the contribution of $X_{(1)},\ldots,X_{(i-1)}$ to $\Psi_{\lambda}(m)$ is $-(i-1)$, and the contribution of $X_{(j+1)},\ldots,X_{(n)}$ is $n-j$. Hence $m\in m(i,j)$ is a zero of $\Psi_{\lambda}$ if an only if $m=1/r$ and

[TABLE]

This last equation is a cubic equation in $r$, and so the value of $r$ that satisfies it (assuming such exists) can be determined exactly using the cubic formula.

4 Numerical Experiments

The $M$-estimator presented here can be thought of as a principled interpolation between the sample mean and the sample median. Because for $u$ small $d(u)\approx u$, as $\lambda\rightarrow 0$, the estimator converges to the sample mean.

At the other extreme, as $\lambda\rightarrow\infty$, all of the $d(\lambda(X_{i}/m-1)$ will evaluate to either 0, 1, or -1. Hence the estimator converges to the sample median.

When $\lambda=\epsilon(1-\epsilon^{2})/(c^{2}+\epsilon^{2})$, we can precisely bound the chance of the relative error being greater than $\epsilon$. However, the estimator can be used for any value of $\lambda$.

For instance, Table 4 records the result of using the estimator for 100 exponential random variables with mean 2 and median $2\ln(2)=1.386\ldots$.