# Estimating $f(\mathbb{E}[X])$ with a guaranteed error performance

Given "black-box" sample access to a random variable $$X$$**, are there results that give an algorithm that approximates $$f(\mathbb{E}[X])$$ with a user-specified error bound, ideally using as few i.i.d. samples from $$X$$ as possible? (See the section "Update (May 11)".)

In this question:

• $$X$$ is a random variable that does not take on a single value with probability 1.
• $$f(x)$$ is a known function belonging to a given class of functions.

The algorithm should—

• Ensure the expected (absolute) error and/or mean squared error of the estimate is within a user-specified error tolerance ($$\epsilon$$), or if that is not possible,
• return an estimate that is within a user-specified tolerance ($$\epsilon$$) on the absolute error or relative error with probability greater than 1 minus $$\delta$$, where $$\delta$$ is user-specified (relative error means $$|\hat\mu/f(\mathbb{E}[X])-1|$$ where $$\hat\mu$$ is the estimate).

## Before the Update

In this section, $$X$$ lies in the interval [0, 1].

The classes of functions $$f$$ that I care about are:

• C1: Continuous functions that map [0, 1] to [0, 1] and are polynomially bounded (meaning that $$f(x)$$ and $$1-f(x)$$ are both bounded from below by $$\min(x^n, (1-x)^n)$$ for some integer $$n$$) (Keane and O'Brien 1994).
• C2: Continuous functions that map [0, 1] to [0, 1] and are not necessarily polynomially bounded.
• C3: Piecewise continuous functions that map [0, 1] to [0, 1] (notably, each piece's domain must be a non-trivial interval and all the pieces must cover all of [0, 1]).

Note the following:

• The following algorithms are different from the algorithms asked for here:
• The Gamma Bernoulli Approximation scheme (Huber 2017) as well as the algorithm of Kunsch et al. 2019 both estimate the mean of $$X$$, namely $$\mathbb{E}[X]$$, rather than a function of that mean, namely $$f(\mathbb{E}[X])$$.
• Some algorithms produce unbiased estimates of $$f(\mathbb{E}[X])$$, but not with a user-specified error (Jacob and Thiery 2015).
• Algorithms like those being asked for here are especially useful because they can help build so-called "approximate Bernoulli factories", or algorithms that approximately sample the probability $$f(\lambda)$$ given a coin with probability of heads of $$\lambda$$.
• I suspect that the algorithm's performance will depend on the "smoothness" of $$f(X)$$, including its modulus of continuity (see also (Holtz et al. 2011)).

## Update (May 11)

The responses so far helped me greatly. I now have the following more focused questions on this topic:

Let $$f(x)$$ be a continuous function.

1. Given the setting at the top of the question, is there an algorithm, besides the algorithm of Kunsch et al. (2019), that can approximate $$\mathbb{E}[X]$$ (or $$f(\mathbb{E}[X])$$) with either a high probability of a "small" absolute error or one of a "small" relative error, when the distribution of $$X$$ is unbounded, and additional assumptions on the distribution of $$X$$ apply, such as—

• being unimodal (having one peak) and symmetric (mirrored on each side of the peak), and/or
• following a geometric distribution, and/or
• having decreasing or nonincreasing probabilities?

Notice that merely having finite moments is not enough (Theorem 3.4, Kunsch et al.). Here, the accuracy tolerances for small error and high probability are user-specified. My article on estimation algorithms already gives a relative-error algorithm for the geometric distribution in a note.

2. How can the method in the answer by "guy" be adapted to discontinuous functions $$g$$, so that the algorithm finds $$g(\mathbb{E}[X])$$ with either a high probability of a "small" absolute error or one of a "small" relative error at all points in [0, 1] except at a "negligible" area around $$g$$'s discontinuities? Is it enough to replace $$g$$ with a continuous function $$f$$ that equals $$g$$ everywhere except at that "negligible" area? Here, the accuracy tolerances for small error, high probability, and "negligible" area are user-specified. Perhaps the tolerance could be defined as the integral of absolute differences between $$f$$ and $$g$$ instead of "negligible area"; in that case, how should the continuous $$f$$ be built? (In this question, $$X$$ lies in the interval [0, 1].)

## Motivation

The answers to this question will help improve my article on randomized estimation algorithms, which describes practical estimation algorithms in a way that programmers can easily implement them.

## References

• Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.
• Huber, M., 2017. A Bernoulli mean estimate with known relative error distribution. Random Structures & Algorithms, 50(2), pp.173-182. (preprint in arXiv:1309.5413v2 [math.ST], 2015).
• Kunsch, Robert J., Erich Novak, and Daniel Rudolf. "Solvable integration problems and optimal sample size selection." Journal of Complexity 53 (2019): 40-67. Also in https://arxiv.org/pdf/1805.08637.pdf.
• Pierre E. Jacob. Alexandre H. Thiery. "On nonnegative unbiased estimators." Ann. Statist. 43 (2) 769 - 784, April 2015.
• Holtz, O., Nazarov, F., Peres, Y., "New Coins from Old, Smoothly", Constructive Approximation 33 (2011).

** And possibly a second source of randomness (such as unbiased random bits). However, the random bits could also come from the random variable $$X$$ itself, assuming $$X$$ doesn't take on the same value with probability 1, via randomness extraction techniques such as the von Neumann extractor (1951). On the other hand, if $$X$$ can be degenerate, this second source will be necessary.

• If you don't know the bounds of $X$ upfront then there's no solution for the identity function, because no matter how many times you sample, you can't rule out the possibility of black swan events that are rare enough not to show up in your sample but big enough to significantly impact the mean. I think my answer here is relevant: stats.stackexchange.com/a/517046 May 2, 2021 at 14:06
• @fblundun : Thank you for the response. I will edit the question, but I don't want to rule out any cases (classes of random variables $X$ and/or classes of functions $f$) for which an algorithm (or at least a practical algorithm) asked for in the question exists. Although such an algorithm might not exist if $X$ is known to be bounded but the bounds are unknown, it might exist under additional assumptions on $X$, such as if the distribution of $X$ is unimodal and symmetric. May 2, 2021 at 14:29
• I feel like I must be missing something because for continuous $f$ and $X$ compactly supported this seems easy. Just use Hoeffding to bound $|\bar X - E(X)|$ with high probability then appeal to the fact that $f$ being continuous means it must be uniformly so to bound $|f(\bar X) - f(E(X))|$ with high probability. How well this works will depend on the modulus of continuity, and you can get (say) bounds in terms of the Lipschitz constant if you are willing to assume the function is Lipschitz.
– guy
May 8, 2021 at 1:42
• After the update, there are well over six questions here, which is too many. Dec 6, 2021 at 19:23
• @MattF: Question 3 was moved to a separate question page. Dec 6, 2021 at 22:00

C1: Yes, by the argument for C2.

C2: Yes, this is possible. The key result here is that any continuous function $$f$$ on $$[0,1]$$ is also uniformly continuous: for any $$\epsilon$$, there is some $$\gamma$$ such that if $$|x_1-x_2|<\gamma$$ then $$|f(x_1)-f(x_2)|<\epsilon$$. So it is enough to have a probability of $$1-\delta$$ of approximating $$E[X]$$ to within $$\gamma$$, and then applying $$f$$ will give a probability of $$1-\delta$$ of approximating $$f(E[X])$$ to within $$\epsilon$$.

C3: No algorithm will work for all functions of the class C3. If we had an algorithm, we could use it on a discontinuous function $$f$$ to tell whether $$E[X]>\frac12$$, but there are always borderline cases that won't work with the desired probability.

In more detail, suppose the user specifies:

• $$f(x) = 1[x>\frac12]$$, i.e. $$1$$ if $$x>\frac12$$, or $$0$$ otherwise.
• $$\epsilon=\frac15$$ as the accuracy -- so getting $$f$$ to an accuracy of $$\epsilon$$ is the same as getting $$f$$ exactly
• $$\delta=\frac14$$, i.e. a probability of $$\frac34$$ of achieving this accuracy

Suppose the algorithm replies that either:

• $$n$$ observations are required, or
• with probability of at least $$1-\frac\delta2$$, at most $$n$$ observations are required.

Now consider a Bernoulli variable $$X$$ which is equally likely to be 0 or 1.

Let $$S$$ be the sample mean of $$X$$ after $$n$$ observations. Then $$P\left[S=\frac12\right]=\binom{n}{n/2}2^{-n}\simeq\sqrt{\frac2{\pi n}}$$ So there is close to a 50% chance that $$S>\frac12$$, and close to a 50% chance that $$S<\frac12$$. There is no way to have a $$\frac34$$ chance of getting the right answer, or even a $$\frac34-\frac18$$ chance of getting the right answer when we allow a $$\frac18$$ possibility of not enough data.

With a slightly longer argument, we could show that the same pattern holds when $$X$$ is Bernoulli and $$E[X]$$ is anything in the range $$\frac12\pm\sqrt{n}/16$$.

So even if we allow a $$\frac\delta2$$ probability that $$n$$ observations are insufficient, there is still no way to get $$f(E[X])$$ which will be accurate to within $$\epsilon$$ the needed $$1-\delta$$ of the time.

• The class C3, as I conceived it, includes only piecewise continuous functions, where each piece's domain is a non-degenerate interval. Your example $f$ is not such a function since one of its pieces lies in $[1/2, 1/2]$, so is not exactly continuous. Does your statement on C3 functions remain valid if each piece's domain is a non-degenerate interval? (Of course, the statement remains valid if we scale the function's range to [0, 1] rather than $[-1, 1]$.) May 8, 2021 at 1:32
• (continuation): Also do you know of algorithms for the classes C1 and C2? May 8, 2021 at 1:39
• I updated to answer for C1 and C2, and to use a function for C3 whose pieces are non-degenerate intervals. May 8, 2021 at 2:30
• Okay, the class C3 of piecewise continuous functions will have to be modified into fully continuous functions, by making each C3 function continuous but "close" enough to the original discontinuous function (in the sense that the new function exactly equals the original function at all but a "negligible" area around its discontinuities). This suggests another accuracy tolerance for the "negligible area", as well as a tradeoff: the smaller the "negligible" area, the higher the new function's modulus of continuity. May 8, 2021 at 2:58
• I believe the proof for the class C3 works because any algorithm for sampling the mean always leaves a window where $f$'s range (difference between maximum and minimum) can be greater than $\epsilon$ (in your example $f$, this window is centered at 1/2 and can be arbitrarily small, and $f$'s range in that window is equal to 1), and this range will never narrow no matter how big $n$ is. Is that correct? May 8, 2021 at 11:26

Elaborating on my comment to the OP, I'll give an explicit answer when $$f$$ is continuous.

Let $$\psi(\epsilon) = \sup\{\psi : |x - y| < \psi \Longrightarrow |f(x) - f(y)| < \epsilon\}$$. Because $$f$$ is continuous on $$[0,1]$$ it is also uniformly continuous, implying that $$\psi(\epsilon) \downarrow 0$$ as $$\epsilon \downarrow 0$$.

By Hoeffding's inequality we have $$\Pr(|f(\bar X) - f(E(X))| \ge \epsilon) \le \Pr(|\bar X - E(X)| \ge \psi(\epsilon)) \le 2\exp\{-2n\psi(\epsilon)^2\},$$ where $$\bar X = \frac{1}{n} \sum_i X_i$$. We can make this less than $$\delta$$ by taking $$n \ge \frac{\log(2/\delta)}{2 \psi(\epsilon)^2}.$$ In the special case where $$f$$ is Lipschitz continuous, i.e., $$|f(x) - f(y)| \le L |x - y|$$, we have $$\psi(\epsilon) \ge \epsilon / L$$, and the bound simplifies to $$n \ge \frac{L^2 \log(2 / \delta)}{2 \epsilon^2}.$$ This can be weakened to allow for non-Lipschitz functions; for example, you get a similar bound by assuming $$|f(x) - f(y)| \le L |x - y|^\alpha$$ for some $$\alpha \le 1$$ (i.e., that $$f$$ is Holder continuous), which should cover just about every continuous function you might care about, at least on $$[0,1]$$.

• The function $\psi$ appears to be different from the modulus of continuity I'm familiar with ($\omega(\epsilon)$; for example, the modulus of continuity for Lipschitz-continuous functions is $\omega(\epsilon) = L\epsilon$ where $L$ is the Lipschitz constant). Can $\omega(\epsilon)$ be transformed easily to $\psi(\epsilon)$? May 8, 2021 at 22:07
• (continuation) Is $\psi(\epsilon)$ derived by solving $\omega(\epsilon)$ for $\epsilon$? In that case, the Hölder case is probably $\psi(\epsilon) = \epsilon^{1/\alpha}/L^{1/\alpha}$ May 8, 2021 at 22:33
• @PeterO. $\psi(\epsilon)$ is a generalized inverse of the modulus of continuity. I don't think I can assume the modulus of continuity is invertible in general, so this seemed like the easiest thing to work with.
– guy
May 8, 2021 at 22:56