The Law of the Unconscious Statistician

When I was a teenager, I spent a good chunk of my free time (and in-class) playing around with numbers and ideas. I grew up with a love of puzzles and was extraordinarily lucky to have maths teachers who would nurture that interest. Similarly, my father was (at the time) a programmer, extremely computer-literate, and taught me how to use computers and introduced me to programming at a young age. These combined into a kid who frequently got told off for trying to figure out how to solve problems the clever way, rather than the way we had been taught. It also meant that kid frequently distracted themselves during lessons by playing around with simple programs and toying with numbers, rather than focusing on the task at hand.

Random number generation was one such idea that occupied my brain for a long period of time. In the middle of a class one day, a younger me was busy not doing what they were told and noticed something unusual.

If you take $n$ randomly generated numbers $x_i$ in the interval $(0,1]$ and take their geometric mean and you graph the results, you get something that looks a bit like a normal distribution. To demonstrate this, I’ve written a quick little function in Javascript:

// Returns a random value ~ U(0, 1]
function pollRandom() {
    //Math.random() returns a value in range [0, 1)
    let x = 1 - Math.random();
    return x; 
}

// Returns the geometric mean of n random values, each ~ U(0, 1]
function geometricMeanOfRandoms(n) {
    var value = 1;
    const pow = (1/n);
    for (let i = 0; i < n; i++) {
        value *= Math.pow(pollRandom(), pow);
    }
    return value;
}

Generating a dataset with a sample size of $10,000$ and graphing the result, here is a histogram (bin width $0.05$ over the unit interval) of what we get when we take $n = 5$…

And here’s what we get for $n = 20$…

And here’s what we get for $n = 50$…

Interesting, isn’t it? As we take larger and larger values of $n$, the distribution of the geometric mean of $n$ randomly chosen values seems to be converging on some value — further testing revealed this converged relatively slowly to somewhere around $0.367$. Had I discovered a bug with how Python (my programming language of choice) generated random numbers? Computers can’t truly generate random numbers as they’re deterministic, but one can design an algorithm that makes almost-random numbers, so could it be that an imperfection in the algorithm was the cause of this?

The alternative — a much more exciting prospect to a younger me — was that I had some kind of discovery on my hands. I checked this idea in other programming languages (such as with the Javascript above) and found similar results, so this must be some kind of mathematical fact, I reasoned.

Constants show up in all manners of unlikely places, so I tested a few constants like $\pi$ and $\phi$ in simple expressions to see if I could get something close to $0.367$. It turned out that taking one over Euler’s constant $e \approx 2.718\dots$ came unreasonably close, and so I came up with this hypothesis:

---

Proposition 1.

For each $i \in \N$, let $X_i \sim U(0, 1]$. Then

$$\lim_{n \to \infty} \left( \mathbb{E}{\left[(X_{1} \times \cdots \times X_{n})^{1/n} \right]} \right) = e^{-1}.$$

---

However, at the time of this, I was 15 - I didn’t have the knowledge on how to prove this, and I wouldn’t for many years. After conferring with my maths teachers and getting many shrugs, I stuffed those notes into a drawer and forgot about them for several years…

…until I was in the middle of my masters and discussing problems with a statistician friend of mine when, all of a sudden, I remembered this problem and posed it to him. He had no idea on how you might go about proving this, and so we set out on attempting to prove it once and for all.

It turns out that one only needs a single key theorem to prove this result — one which is likely my favourite name¹ for a theorem in any field:

---

Theorem.

Law of the unconscious statistician (LOTUS). The expected value of a measurable function of $X$, $g(X)$, given that $X$ has probability distribution function $f(x)$, is given by the inner product of $f$ and $g$:

$$\mathbb{E}{\left[ g(X) \right]} = \int_\R g(x) f(x) \; dx.$$

---

In my several statistics modules I took as an undergraduate student, this is something I don’t recall ever being covered — and apparently something my statistician friend hadn’t heard of either. In hindsight, this theorem is actually just a consequence of the definition of the expected value — it’s an immediately corollary with no work required to get there. We spent many hours trying overly clever tricks to solve this problem, when it turns out that the solution was to just go back to the definition and work from there.

Let’s follow this proof through and prove proposition 1.

---

Proof of Proposition 1.

Before we start, we need to argue that this limit exists. Each $X_i ~ U(0,1]$ and so the product of any number of copies must also be in the interval $(0, 1]$. The geometric mean is a continuous, infinitely differentiable function on $\R^{\geq0}$ and so taking the geometric mean of some product $X_1 \times \cdots \times X_n$ should respect limits.

First, examining the expression on the left, we can do some simple tricks to get this closer to a form where we can use LOTUS. Focusing on the expression inside the limit first, the distributive property of exponents lets us get:

$$\mathbb{E}{\left[(X_{1} \times \cdots \times X_{n})^{1/n} \right]} = \mathbb{E}{\left[ \prod^{n}_{i = 1} \left(X_{i}\right)^{1/n} \right]}.$$

As each $X_i$ is independent of one another, we know that the mean of their products is the product of their means, and so

$$\mathbb{E}{\left[ \prod^{n}_{i = 1} \left(X_{i}\right)^{1/n} \right]} = \prod^{n}_{i = 1} \mathbb{E}{\left[ \left(X_{i}\right)^{1/n} \right]}.$$

However, as each $X_i \sim U(0, 1]$ is independent, it should reason that each $\mathbb{E}[(X_{i})^{1/n}]$ is equal to one another, and so

$$\begin{align*} \prod^{n}_{i = 1} \mathbb{E}{\left[ \left(X_{i}\right)^{1/n} \right]} &= \prod^{n}_{i = 1} \mathbb{E}{\left[ \left(X_{1}\right)^{1/n} \right]} \\ &= {\mathbb{E}{\left[ (X_{1})^{1/n} \right]}}^{n}. \end{align*}$$

Thus, by taking $X \sim U(0, 1]$, we can rewrite the left-hand side of the original equation as

$$\lim_{n \to \infty} \left( \mathbb{E}{\left[\left(\prod^{n}_{i = 1} X_{i}\right)^{1/n} \right]} \right) = \lim_{n \to \infty} \left( {\mathbb{E}{\left[ (X_{1})^{1/n} \right]}}^{n} \right).$$

This is now where LOTUS comes into play so we can resolve the expected value $E[X^{1/n}]$. The probability distribution function $f(x)$ for the uniform distribution $U(0, 1]$ is extremely simple:

$$f(x) = \begin{cases} 1 & x \in (0,1] \\ 0 & \textnormal{otherwise} \end{cases}$$

and so we calculate the expected value

$$\begin{align*} \mathbb{E}[X^{1/n}] &= \int_{\R} f(x) x^{1/n} \; dx \\ &= \int^{1}_{0} x^{1/n} \; dx \\ &= 1 - \frac{1}{n+1} \end{align*}$$

and, putting this together with the limit, we obtain

$$\lim_{n \to \infty} \left( \mathbb{E}{\left[(X_{1} \times \cdots \times X_{n})^{1/n} \right]} \right) = \lim_{n \to \infty} \left( \left( 1 - \frac{1}{n+1} \right)^n \right).$$

This final part should look strikingly familiar to anyone who has studied introductory calculus, because we can rewrite this last part as a known identity for powers of $e$:

$$\lim_{n \to \infty} \left( \left( 1 + \frac{k}{n} \right)^n \right) = e^k.$$

Putting these together, we finally obtain the result we were looking for:

$$\begin{align*} \lim_{n \to \infty} \left( \mathbb{E}{\left[(X_{1} \times \cdots \times X_{n})^{1/n} \right]} \right) &= \lim_{n \to \infty} \left( \left( 1 - \frac{1}{n+1} \right)^n \right) \\ &= e^{-1}. \end{align*}$$

---

Even though it’s been nearly five years (oh dear, five years!) since I completed my master’s degree, this is still something I often reflect fondly on: that a curious young student could make an educated guess, and nearly a full decade later only then have the tools to prove themself right. While the key theorem that made it all work is but just a consequence of definition, it was nonetheless a pleasant surprise to discover it.

One can, of course, also solve for the standard deviation using similar techniques to the above. The standard deviation $\sigma_n$ of $\left( \prod X_i \right)^{1/n}$ is defined as follows:

$$\begin{align*} \sigma^2_n &= {\mathbb{E} \left[ \left( \prod {X_{i}} \right)^{2/n} \right]} - {\mathbb{E} \left[\prod {X_{i}}^{1/n} \right]}^2 \\ &= {\mathbb{E} \left[ \prod \left( {X_{i}} \right)^{2/n} \right]} - \left( 1 - \frac{1}{n+1} \right)^{2n} \\ &= {\mathbb{E} \left[ \left( X \right)^{2/n} \right]}^n - \left( 1 - \frac{1}{n+1} \right)^{2n} \\ &= \left(\int^{1}_{0} x^{2/n} \; dx \right)^n - \left( 1 - \frac{1}{n+1} \right)^{2n} \\ &= \left(1 - \frac{2}{n+2}\right)^n - \left( 1 - \frac{1}{n+1} \right)^{2n} \\ \end{align*}$$

Interestingly, $\sigma^2$ turns out to be the difference between (roughly) the $n^{\textnormal{th}}$ estimate for $e^{-2}$ and the square of the $n^{\textnormal{th}}$ estimate for $e^{-1}$. Graphing the standard deviation $\sigma_n$ as $n$ grows obtains us this lovely graph:

After seeing this graph for the first time, I reckoned there was probably a pleasantly coincdentally estimate for what this standard deviation would be. It turns out there is:

---

Proposition 2.

$$\sigma^2_n \approx e^{-2} n^{-1}$$

---

Experimental analysis suggests that the error term is miniscule. To prove this is likely to take a lot of algebraic manipulation of the two products, though I suspect there’s probably a quite pleasant way to do this via integral representations of them. However, a proof of that is outside of the scope of this work. If you’re reading this and you have a fun proof of this: let me know!

Footnotes

1. I could swear that, at the time of first discovering this theorem for myself, the name of it was rather the law of the drunken statistician. However, try as I might, I since have not been able to find a source that uses this name. ↩