Robust estimation of the mean with bounded relative standard deviation

Mark Huber

REVIEW 4 major objections 4 minor 9 references

Robust estimation of the mean with bounded relative standard deviation

T0 review · 4 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read A modified M-estimator proves an (ε,δ)-relative-error guarantee for the mean using at most 2(c²ε⁻²+1)/(1−ε²) ln(2/δ) draws.

desk verdict A real constant-factor improvement in sample complexity for bounded-relative-variance mean estimation, but the printed proof has a definitional gap that needs a mechanical fix. read the letter →

arxiv 1908.05386 v1 pith:OO7HZCKB submitted 2019-08-15 stat.CO math.PR

Mark Huber This is my paper · ORCID

classification stat.COmath.PR MSC 62F3562G3565C05

keywords meanestimationrelativeerrorrandomizedapproximationschemeM-estimatorChernoffboundsamplecomplexityrobuststatistics

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper asks how many independent draws of a nonnegative random variable $X$ are needed to estimate its mean $\mu$ to within a relative error $\epsilon$ with failure probability at most $\delta$, when all that is known is an upper bound $c$ on the relative standard deviation. The main result is a new estimator, a modified Catoni $M$-estimator, that uses at most $\lceil 2(c^2\epsilon^{-2}+1)(1-\epsilon^2)^{-1}\ln(2/\delta)\rceil$ draws. The previous best used the same setup but carried an extra $(1+\epsilon)^2$ factor; the new bound is within a factor $1+O(\epsilon^2)$ of the information-theoretic minimum, so it nearly reaches the performance that normal samples would allow. The improvement comes from choosing the estimator's scale parameter $\lambda$ so that the two tail probability bounds are equal, balancing the worst case instead of letting one tail dominate. The same estimator is also fast in practice and can be computed exactly with square and cube roots.

What carries the argument

The central object is the Catoni-type $\Psi_\lambda$ estimator, a robust $M$-estimator whose influence function $d_h$ is close to the identity for small deviations from $\mu$ but saturates at $\pm 5/6$ for large ones, interpolating between sample mean and median. Lemma 1 sandwiches $d_h(u)$ between $-\ln(1-u+u^2/2)$ and $\ln(1+u+u^2/2)$, so exponentiating the score converts tail probabilities into products of factors that depend only on $\mathbb{E}[X/m-1]$ and $\mathbb{E}[(X/m-1)^2]$; that is, only on the mean and variance. Choosing $\lambda$ to balance the two tail bounds produces the gap $1-\frac{1}{2}\epsilon^2(c^2+\epsilon^2)^{-1}(1-\epsilon^2)$, and the sample count follows as $(1/\text{gap})\ln(2/\delta)$.

What would settle it

Take $X=1+Z$ where $Z$ is $\pm 1$ with equal probability (so $\mu=1$ and $\mathrm{SD}(X)=1$), choose $\epsilon=0.1$ and $\delta=0.01$, compute $n$ from Theorem 1, and simulate many batches; if the observed fraction of roots of $\Psi_\lambda$ falling outside $(0.9,1.1)$ is substantially larger than $\delta$, the claimed bound is false.

Watch

Extended reading notes

Core claim

Under the assumptions $X\ge 0$, $\mathbb{E}[X]=\mu$, and $\mathrm{SD}(X)\le c\mu$, there is an $(\epsilon,\delta)$-randomized approximation scheme whose estimate is any root of the $M$-estimator equation $\Psi_\lambda(m)=\sum_i \lambda^{-1}d_h(\lambda(X_i/m-1))=0$, with $d_h(u)=(5/6)\mathbf{1}(u>1)+(u-u^3/6)\mathbf{1}(|u|\le 1)-(5/6)\mathbf{1}(u<-1)$, that uses at most $\lceil 2(c^2\epsilon^{-2}+1)(1-\epsilon^2)^{-1}\ln(2/\delta)\rceil$ draws. The proof bounds $\Pr(\Psi_\lambda((1+\epsilon)\mu)\ge 0)$ and $\Pr(\Psi_\lambda((1-\epsilon)\mu)\le 0)$ separately by exponentiating $d_h$, then chooses $\lambda=\epsilon(1-\epsilon^2)/(c^2+\epsilon^2)$ to equalize the two quadratic upper bounds. This replaces the $(1+\epsilon)^2$ factor of the previous scheme with $(1+\epsilon^2/c^2)/(1-\epsilon^2)$, a second-order correction.

Load-bearing premise

The theorem is only usable when the user already knows a finite constant $c$ such that the standard deviation of $X$ is at most $c\mu$; if $c$ is unknown, too small, or infinite, the stated sample bound does not guarantee the $(\epsilon,\delta)$ property.

Editorial extensions

If this is right

Any randomized approximation scheme that simulates a nonnegative random variable with a known relative standard deviation $c$ can now meet the $(\epsilon,\delta)$ guarantee with $\lceil 2(c^2\epsilon^{-2}+1)(1-\epsilon^2)^{-1}\ln(2/\delta)\rceil$ draws, improving the previous $(1+\epsilon)^2$ overhead.
Because the normal-distribution lower bound is $2c^2\epsilon^{-2}\ln(1/\delta)$ plus lower-order terms, the new scheme is within a $1+O(\epsilon^2)$ factor of the best possible sample count for worst-case distributions carrying only mean and variance information.
The estimator is implementable: an approximate root by binary search on sorted data costs $O(n\ln n)$ time, and an exact root is available from the cubic formula over at most $2n$ intervals, so the improved bound is not merely asymptotic.
When $c$ is known only approximately or the relative standard deviation is infinite, the estimator still converges to a fixed centrality value rather than diverging like the sample mean, so the same code behaves as a robust location estimator.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

One natural extension the paper does not prove: if $c$ is unknown, a pilot sample could estimate it and the same balanced choice of $\lambda$ might yield an adaptive $(\epsilon,\delta)$ bound; the numerical experiments on unknown and infinite relative standard deviations are consistent with this working in practice.
The tail-balancing trick is not tied to this particular influence function; other bounded influence functions whose logarithms are controlled by quadratics could receive the same optimization, potentially transferring the constant improvement to related estimators.
Because the exact root reduces to a cubic equation, the method suggests a general recipe for exactly solvable $M$-estimators with piecewise polynomial influence functions, where roots of higher-degree polynomials would be computed numerically but the same variance-only tail analysis would still apply.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

This is a short paper from Huber with a genuine improvement. The main theorem says: if X ≥ 0 has mean μ and standard deviation at most cμ, then for any ε, δ there is an (ε,δ)-ras using at most ⌈2(c²ε⁻²+1)/(1−ε²) ln(2/δ)⌉ draws. That is better than the previous bound 2c²ε⁻²(1+ε)² ln(2/δ), and the gain is a true second-order improvement, not just noise. The mechanism is clean: a Catoni-type M-estimator with a new score function d_h and a tail-balancing choice of λ that equalizes the two Chernoff upper bounds. I re-did the algebra in Lemmas 3–7 and it works, provided the estimator is read as a sum. The exact root-finding idea is also nice, and the paper is honest about the mean/median interpolation behavior.

Now the soft spots. (1) The timing code in Section 6 uses (1+ε²)⁻¹ instead of (1−ε²)⁻¹ in the sample count, so Table 3 does not measure the theorem's algorithm. Reproducibility bug, not a math flaw. (2) Section 3.2's exact method is not correct as printed: the interval m(i,j) definition, the claim about each point entering/leaving once, and the cubic equation's coefficients all need checking. I did not fully verify the cubic, but you cannot trust that description as written. (3) The load-bearing issue: Lemma 2 defines Ψλ as the average (1/n)Σ λ⁻¹ d_h(λ(X_i/m−1)), while Lemma 3's proof and the code use the unnormalized sum Σ d_h(...). Under the average definition, the Chernoff product does not factor into [1+E[u]+E[u²]/2]^n; you get a product involving E[d_h]/n. Also, in Lemma 3 the variable u already contains λ, but Lemma 4 is written with an extra λ, and the comparison drops it. None of this is fatal. The estimator in the code and the bound in Theorem 1 are consistent once you redefine Ψλ as the sum and track λ explicitly. But as printed, the proof-to-theorem link fails. This is exactly what the stress-test note says, and I agree.

Who is this for? People working with randomized approximation schemes where the relative variance is bounded—DNF counting, convex volume, partition functions—and people who use Catoni-style robust mean estimators. The constant-factor gain is modest but real, and the method is easy to implement once you know the λ formula. The paper deserves a serious referee; desk rejection would be wrong. I would cite this in my own work after the estimator is defined consistently and the code/timing errors are fixed. My recommendation: send to a journal, require a revision that makes Lemma 2 the sum form, redoes Lemma 3 with λ explicit, and corrects Section 3.2 and the timing script.

Referee Report

4 major / 4 minor

Summary. The paper proposes a modified Catoni-Guilini M-estimator for the mean of a nonnegative random variable with known upper bound c on the relative standard deviation. The estimator uses a smoothly trimmed deviation function d_h with a scale parameter λ chosen analytically to balance two Chernoff-type upper bounds, yielding the finite-sample guarantee in Theorem 1: for P(X≥0)=1, E[X]=µ, and SD(X)≤cµ, there is an (ε,δ)-ras using at most ⌈2(c²ε⁻²+1)(1−ε²)⁻¹ ln(2/δ)⌉ draws. This improves the previous constant (1+ε)² to (1+ε²/c²)/(1−ε²). The paper also sketches an approximate binary-search root finder, an 'exact' method based on solving cubic equations, and numerical experiments on exponential, Cauchy-like, and SDE-driven data, with R code included.

Significance. If the main theorem is established, the result is a genuine improvement in the sample-complexity constant for relative-error randomized approximation schemes, moving from a first-order to a second-order gap relative to the Gaussian lower bound. The estimator is parameter-free in the sense that λ is determined explicitly from ε and c, and the sample bound is explicit and non-asymptotic. The inclusion of R code and numerical experiments is a strength, as is the clear statement of the theorem. However, the significance is conditional: the printed proof of Theorem 1 contains definitional and algebraic inconsistencies that must be resolved before the claimed guarantee can be accepted.

major comments (4)

[Section 2, Lemmas 2 and 3] Lemma 2 defines Ψλ(m) = (1/n) Σᵢ λ⁻¹ d_h(λ(Xᵢ/m−1)), but Lemma 3's proof writes exp(λΨλ(m)) as a product of n independent terms exp(λ d_h(Xᵢ/m−1)), dropping both the 1/n and the λ⁻¹ factors. With the printed average definition, λΨλ(m) = (1/n) Σᵢ d_h(λ(Xᵢ/m−1)), so Markov's inequality gives (E[exp(d_h(U)/n)])ⁿ, not (E[exp(d_h(U))])ⁿ; the displayed exponential tail bound does not follow. The proof also sets u=λ(X/m−1) while the lemma statement uses u_ε = X/((1+ε)µ)−1 without λ, so the moment inequalities in Lemma 4 are not the ones being used. Because Theorem 1's sample count is derived from exactly these tail bounds, the proof-to-theorem link is invalid as printed. The intended argument is repairable by defining Ψλ as the unnormalized sum Σᵢ d_h(λ(Xᵢ/m−1)), as in the code, and redoing Lemmas 3 and 5 with λ tracked throughout.
[Section 3.2, exact method] The description of the exact root-finding method contains several errors that make it incorrect as written. The interval defining m(i,j) is given as [m−λ, m+λ], but the condition |λ(X/m−1)|≤1 is equivalent to X∈[m(1−1/λ), m(1+1/λ)], not [m−λ, m+λ]. The contributions of points outside this interval are stated as −(i−1) and n−j, but d_h takes values ±5/6 and the definition of Ψλ includes λ⁻¹, so the outside contributions are ±(5/6)λ⁻¹, not ±1. The cubic equation in r also omits these λ and 5/6 factors. Thus the exact method does not compute a root of Ψλ as defined; this does not affect the sample-complexity theorem but invalidates a stated practical contribution.
[Section 6, time_psi_h function] The R function time_psi_h computes n <- ceiling(2 * (c^2 * epsilon^(-2) + 1) * (1 + epsilon^2)^(-1) * log(2 / delta)), but Theorem 1 and Lemma 7 require (1−ε²)⁻¹, not (1+ε²)⁻¹. Consequently the timing experiments do not implement the proposed estimator's sample count, and the improvements reported in Table 3 are not attributable to the bound claimed in the paper. This is a load-bearing error for the numerical comparison, though not for the theoretical theorem itself.
[Section 2, after Lemma 3] The sentence 'by the union bound the probability that Ψ has a root in [(1−ε)µ,(1+ε)µ] is at least δ/2+δ/2=δ' is incorrect as worded: the union bound gives error probability at most δ/2+δ/2=δ, so the success probability is at least 1−δ. This is a local wording issue, but it obscures the logic of the proof of Theorem 1 and should be fixed.

minor comments (4)

[Section 2, Lemma 5] Lemma 5 says 'Combine Lemmas 3, 4, and 6' but Lemma 6 is stated after Lemma 5; reorder the lemmas or adjust the citation.
[Section 4, first paragraph] The text says that as λ→∞ the values d(λ(Xᵢ/m−1)) evaluate to either 0, 1, or −1; for the function d_h used in the paper, the limiting values are ±5/6, not ±1, and 0 does not occur for Xᵢ≠m.
[Table 3 caption] The caption says 'Behavior of the sample mean, sample median, and M-estimator for 100 draws from the absolute value of a Cauchy distribution', but the table reports timings for the CG and new estimators on an SDE example; the caption should match the content.
[Section 1, introduction] There is a typo, 'appoximation' for 'approximation', and the lower-bound statement 'n ≥ 2c²ε⁻²[ln(1/δ)−o(δ)]' uses ln(1/δ) while the upper bounds use ln(2/δ); this is consistent with the constant-2 convention but should be made explicit.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the sample-complexity bound is derived from Chernoff/Markov bounds and analytic inequalities with lambda chosen by balancing upper bounds, not by fitting or by importing the conclusion.

full rationale

The central claim, Theorem 1, is the sample bound for the robust M-estimator. The derivation chain is self-contained: Lemma 1 bounds d_h by logarithms; Lemma 3 applies Markov's inequality to the exponential of the estimator; Lemma 4 uses only E[X]=mu and Var(X)<=c^2mu^2; Lemmas 5-7 combine these to choose lambda=epsilon(1-epsilon^2)/(c^2+epsilon^2) by equalizing the two quadratic Chernoff upper bounds. No parameter is fitted to data, and no numerical experiment is used to set constants. The self-citation to [5] appears in Lemma 2, which restates the previous algorithm's guarantee, but Theorem 1's improved bound does not depend on Lemma 2; the new lambda and sample count are established inside the present paper. The apparent inconsistency between the average definition of Psi_lambda in Lemma 2 and the product form used in Lemma 3 is a proof-formulation or correctness issue, not a circularity: it does not presuppose Theorem 1's conclusion. Numerical experiments and timings are demonstrations of behavior, not fitted inputs. There is no step where an input is renamed as a prediction, and no load-bearing conclusion is justified solely by an unverified self-citation. Hence the paper has no significant circularity.

Assumptions & free parameters 0 free parameters · 5 assumptions · 0 invented entities

The central claim rests on standard probabilistic inequalities and the explicit moment assumptions of the problem. No free parameters are fitted to data: λ is chosen analytically to equalize two upper bounds. No new entities are introduced.

assumptions (5)

domain assumption The random variables X_i are iid with P(X≥0)=1, E[X]=µ, and Var(X) ≤ c²µ².
This is the problem setting of Theorem 1; nonnegativity is essential for the sign and tail logic.
standard math Markov's inequality and the Chernoff product bound apply to the estimator when expressed as a sum of n independent contributions.
Lemma 3 exponentiates the sum and factorizes the expectation into a product; this requires the sum form rather than the average written in Lemma 2.
standard math The functions d_L(u) = -ln(1-u+u²/2) and d_U(u) = ln(1+u+u²/2) satisfy d_L ≤ d_h ≤ d_U for all real u.
Proved in Lemma 1 using calculus; the central tail bounds depend on this.
standard math Catoni's lower bound n ≥ 2c²ε^{-2}[ln(1/δ)-o(δ)] from [1] is correct and applies to all estimators of this type.
Used only to contextualize near-optimality, not to derive the upper bound.
standard math The function Ψλ(m) is continuous and nonincreasing in m.
Used to conclude a sign change between (1-ε)µ and (1+ε)µ ensures a root in that interval; follows from continuity and monotonicity of d_h.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Robust estimation of the mean with bounded relative standard deviation." pith.science (2026). https://pith.science/paper/OO7HZCKB

@misc{pith2026190805386,
  author       = {Pith},
  title        = {Pith review of: Robust estimation of the mean with bounded relative standard deviation},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/OO7HZCKB}},
  note         = {Machine review of arXiv:1908.05386}
}

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.

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

9 extracted references · 7 canonical work pages

[1]

Catoni, O.: Challenging the empirical mean and empirical variance: A deviation study. Ann. Inst. H. Poincaré Probab. Statist. 48, 1148–1185 (2012)

work page 2012
[2]

Submitted

Catoni, O., Giulini, I.: Dimension free PAC-Bayesian bounds for matrices, vectors, and linear least squares regression with a random design (2017). Submitted. arXiv: 1712.02747

arXiv 2017
[3]

Chernoﬀ, H.: A measure of asymptotic eﬃciency for tests of a hypothesis based on the sum of observations. Ann. of Math. Stat. 23, 493–509 (1952)

work page 1952
[4]

Dyer, M., Frieze, A., Kannan, R.: A random polynomial-time algorithm for approximating the volume of convex bodies. J. Assoc. Comput. Mach.38(1), 1–17 (1991)

work page 1991
[5]

Random Structures Algorithms To appear

Huber, M.: An optimal(ϵ,δ)-approximation scheme for the mean of random variables with bounded relative variance. Random Structures Algorithms To appear

work page
[6]

Huber, M.: Approximation algorithms for the normalizing constant of Gibbs distributions. Ann. Appl. Probab.51(1), 92–105 (2015). arXiv:1206.2689

work page Pith review arXiv 2015
[7]

Huber, P.J.: Robust estimation of a location parameter. Ann. Math. Statist.35(1), 73–101 (1964). DOI 10.1214/aoms/ 1177703732. URL https://doi.org/10.1214/aoms/1177703732

arXiv 1964
[8]

In: Proc

Karp, R.M., Luby, M.: Monte-carlo algorithms for enumerating and reliability problems. In: Proc. FOCS, pp. 56–64 (1983)

work page 1983

Show all 9 references

[9]

Journal de Mathématique Pures et Appliquées2(12), 177–184 (1867)

Tchebichef, P.: Des valeurs moyennes. Journal de Mathématique Pures et Appliquées2(12), 177–184 (1867)

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

Catoni, O.: Challenging the empirical mean and empirical variance: A deviation study. Ann. Inst. H. Poincaré Probab. Statist. 48, 1148–1185 (2012)

work page 2012

[2] [2]

Submitted

Catoni, O., Giulini, I.: Dimension free PAC-Bayesian bounds for matrices, vectors, and linear least squares regression with a random design (2017). Submitted. arXiv: 1712.02747

arXiv 2017

[3] [3]

Chernoﬀ, H.: A measure of asymptotic eﬃciency for tests of a hypothesis based on the sum of observations. Ann. of Math. Stat. 23, 493–509 (1952)

work page 1952

[4] [4]

Dyer, M., Frieze, A., Kannan, R.: A random polynomial-time algorithm for approximating the volume of convex bodies. J. Assoc. Comput. Mach.38(1), 1–17 (1991)

work page 1991

[5] [5]

Random Structures Algorithms To appear

Huber, M.: An optimal(ϵ,δ)-approximation scheme for the mean of random variables with bounded relative variance. Random Structures Algorithms To appear

work page

[6] [6]

Huber, M.: Approximation algorithms for the normalizing constant of Gibbs distributions. Ann. Appl. Probab.51(1), 92–105 (2015). arXiv:1206.2689

work page Pith review arXiv 2015

[7] [7]

Huber, P.J.: Robust estimation of a location parameter. Ann. Math. Statist.35(1), 73–101 (1964). DOI 10.1214/aoms/ 1177703732. URL https://doi.org/10.1214/aoms/1177703732

arXiv 1964

[8] [8]

In: Proc

Karp, R.M., Luby, M.: Monte-carlo algorithms for enumerating and reliability problems. In: Proc. FOCS, pp. 56–64 (1983)

work page 1983