# Expectation of the product of iid random variables

If we have iid random variables $$X_1,X_2,...,X_N$$ with $$\mathbb{E}X_i=\mu$$, is it true that $$\mathbb{E}\prod X_i=\mu^N$$?

I had no doubt that this is true, until I tried it out with Python, using random.normalvariate() to generate the set of samples, and surprisingly found that the product of all these data points are generally a lot smaller than $$\mu^N$$.

For example, I used that function to generate 2 million (there's absolutely no need to have a dataset of this size, but I went for it anyways) data points that are, supposedly, distributed as $$N(1, 0.2)$$. I was hoping for their product to scatter around 1 as I repeat the trial, but instead I got numbers $$\sim\pm10^{-18500}$$ constantly.

For what it's worth, I've tried sample sizes ranging from 1 to 1000, and they all fell below their respective $$\mu^N$$, significantly -- the difference was visible on a log-scaled plot.

I suspected random.normalvariate() generated something whose PDF is not $$N(1,0.2)$$. But I plotted that 2 million data points and got a perfect bell-shape curve.

I suspected that the data are correlated among themselves. But I tried to compute $$\mathbb{E}\prod X_i$$ with $$\text{corr}(X_i,X_j)=\rho$$ and found that my calculation could not explain it. I'm not hundred percent confident in my calculation though.

And I tried to understand it intuitively, and had the following thought. Say we have a lot of $$N(1,0.2)$$ data points. It is conceivable that they lie symmetrically around 1. So, we can group the data points into pairs that are roughly like $$\{(1-a_n),(1+a_n)\}$$. This should be feasible when the sample size is large enough. But each pair has their product being less than 1. Therefore, the total product is a lot smaller than 1.

So, this seems to me a paradox. I can't dissuade myself from $$\mathbb{E}\prod X_i=\mu^N$$, but neither can I find the loophole in the thought above (or my empirical tests). I feel that I must have made a blatant mistake somewhere, but I can't locate that error. Please help me out if you know the answer!

• If your simulations are giving you $\sim\pm10^{-18500}$ then you really do not know whether their mean is $1$. Try it with products of $2,3,10,100$ terms with a million simulations each and see if the means are about $1$ - they should be Jun 8 at 11:58
• stats.stackexchange.com/questions/3707 is relevant. There is no paradox, though: when $N$ is very large, it is almost certain the product will be close to zero. The mean remains $\mu^N,$ though, because there is a tiny chance--too tiny to be detected by your simulations--that the product will be incredibly large. The error here, if there is any, lies in relying on an inadequate simulation to test your intuition. Sometimes mathematical reasoning is essential.
– whuber
Jun 8 at 13:29
• The pairing argument is interesting, but when you reflect on it, how would it be relevant? How does it apply to products of $N$ values? Consider $N=2,$ for instance, and (with no loss of generality) $\mu=1.$ For every occurrence of $1-a$ in $X_1$ there will be a $1+a$ in $X_1$ -- and likewise for $X_2.$ That gives us four products to track, equal to $(1-a)^2,$ $(1+a)^2,$ $(1-a)(1+a),$ and $(1+a)(1-a).$ They average out exactly to $1.$ It does not suffice to focus on the last two terms. Notice, too, that only one of these four terms, $(1+a)^2,$ exceeds $1.$ (That's a strong hint!)
– whuber
Jun 8 at 14:59
• @whuber That makes so much sense! The hint is very helpful as well. Thank you! Jun 8 at 16:03
• Here is a related question dealing with the same topic. stats.stackexchange.com/questions/480967 Jun 9 at 10:53

First, let's establish the correct identity.

When $$X_1, \ldots, X_N$$ are independent variables with finite expectations $$\mu_i=X_i,$$ then by laws of conditional expectation,

$$E\left[\prod_{i=1}^N X_i\right] =E\left[X_N E\left[\prod_{i=1}^{N-1} X_i \mid X_{N}\right]\right] = E\left[X_N \prod_{i=1}^{N-1} \mu_i\right] = \mu_N\prod_{i=1}^{N-1} \mu_i= \prod_{i=1}^N \mu_i$$

gives a proof by mathematical induction (beginning with the base case $$N=1$$ where $$E\left[\prod_{i=1}^N X_i\right] = E\left[X_1\right] = \mu_1 = \prod_{i=1}^N \mu_i$$ is trivially true).

Now, let's find an explanation for the simulation results.

The pairing argument in the question is an interesting one, because it shows that when multiplications by $$(1-a)$$ and $$(1+a)$$ occur in equal numbers (approximately $$N/2$$ each), the net product is $$(1-a^2)^{N/2}\approx \exp(-Na^2/2).$$ This suggests that when $$N$$ is sufficiently large, it's nearly certain that the product will be tiny--certainly less than the common mean of $$1.$$ The reason this is not a paradox is that there will be a vanishingly small--but still positive--probability of yielding a whopping big number on the order of $$(1+a)^N \approx \exp(aN).$$ This rare chance of a huge product balances out all the tiny products, keeping the mean at $$1.$$

It is not easy to analyze the product of many Normal variables. Instead, we may gain insight from a simpler case. Let $$Y_1, Y_2, \ldots,$$ be a sequence of independent Rademacher variables: that is, each of these has a $$1/2$$ chance of being either $$1$$ or $$-1.$$ Pick some number $$0 \lt a \lt 1$$ and define $$X_i = 1 + aY_i,$$ so that each $$X_i$$ has equal chances of being $$1\pm a.$$ Clearly $$E[X_i] = 1 = \mu_i$$ for all $$i.$$

Consider the product of the first $$N$$ of these $$X_i.$$ Suppose, in a simulation, that $$k$$ of these values equal $$1-a$$ and (therefore) the remaining $$N-k$$ of them equal $$1+a.$$ The product then is $$(1-a)^k(1+a)^{N-k}.$$ How small must $$k$$ be for this product to exceed $$1$$?

Given $$N$$ and $$a,$$ we must solve the inequality

$$(1-a)^k(1+a)^{N-k} \ge 1$$

for $$k.$$ By taking logarithms, this is equivalent to

$$k \le N \frac{\log(1+a)}{\log(1+a) - \log(1-a)}.$$

Because each $$X_i$$ has equal and independent chances of being $$1\pm a,$$ the distribution of $$k$$ is Binomial$$(N, 1/2),$$ which even for moderate sizes of $$N$$ ($$N \ge 10$$ is fine) is nicely approximated by a Normal$$(N/2, \sqrt{N}/2)$$ distribution. Thus, the chance that the product is $$1$$ or greater will be close to the value of the standard Normal distribution at $$Z$$ (the tail area under the Bell Curve left of $$Z$$) where

$$Z = \frac{N \frac{\log(1+a)}{\log(1+a) - \log(1-a)} - \frac{N}{2}}{\sqrt{N}/2} = \text{constant}\times \sqrt{N}.$$

You can see where this is going! As $$N$$ grows large, $$Z$$ is pushed further out to the left, making it less and less likely to observe any product greater than $$1$$ in a simulation.

In the question, where the standard deviation is $$0.2,$$ the value $$a=0.2$$ will closely reproduce the simulation behavior. In this case the constant is

$$\text{constant} = \frac{2\log(1+a)}{\log(1+a)-\log(1-a)} - 1 = -0.100\ldots$$

Taking $$N=2\times 10^6,$$ for instance, as in the question, compute $$Z \approx -142.$$ The chance that $$k$$ is small enough to produce a value this negative is less than $$10^{-10000}.$$ You can't even represent that in double precision floats. It would take far more than the age of the universe to create a simulation that had the remotest chance of producing such an imbalance between the $$1+a$$ and $$1-a$$ values that the product exceeds $$1.$$

In short, for all practical purposes, when $$N$$ is sufficiently large ($$N \gg 5000$$ will do when $$a=1/5$$), you will never observe a value above $$1$$ in this simulation, even though the mean of the product is $$1.$$

• (+1) Beautiful explanation. Jun 8 at 17:08
• I knew the lion by his claw, before checking the name on the answer. Jun 10 at 1:23

As an addendum to @whuber's explanation, consider that \begin{align} \prod_{i=1}^N \vert X_i\vert &= \exp\{ \sum\nolimits_{i=1}^N \overbrace{\ln \vert X_i\vert}^{Y_i} \}\\ &= \exp\{ \sum\nolimits_{i=1}^N Y_i \}\\ &= \exp\{ \sum\nolimits_{i=1}^N (Y_i - \tilde \mu)+N\tilde\mu\}\qquad\qquad\mathbb \ \ E[Y_i]=\tilde\mu\\ &= \exp\left\{\frac{\tilde\sigma\sqrt{N}}{\sqrt{N}\tilde\sigma} \sum\nolimits_{i=1}^N (Y_i - \tilde \mu)+N\tilde\mu\right\}\qquad\text{var}(Y_i)=\tilde\sigma^2\\ &\approx \exp\{\tilde\sigma\sqrt{N}\zeta+N\tilde\mu\}\qquad\quad\qquad\qquad\zeta\sim\mathcal N(0,1)\\ \end{align} which shows that the expression tends to behave more and more erratically as $$N$$ increases. (For $$|\mu|$$ large enough, $$\tilde \mu\approx\log|\mu|$$.) Note that $$\xi=\exp\{\tilde\sigma\sqrt{N}\zeta+N\tilde\mu\}$$ is a log-normal $$\mathcal L\mathcal N(N\tilde\mu,N\tilde\sigma^2)$$ random variable with $$\mathbb E[\xi]=\exp\{N\tilde\mu+N\tilde\sigma^2/2)\qquad \text{var}(\xi)=[\exp\{N\tilde\sigma^2/2\}-1]\exp\{2N\tilde\mu+N\tilde\sigma^2\}$$ As a last remark, not quite related with the spirit of the question, while $$\prod_{i=1}^NX_i$$ is indeed an unbiased estimator of $$\mu^N$$, a more efficient estimator based on the $$N$$ rv's $$X_i$$ would be $$\mathbb E\left[\prod_{i=1}^NX_i\Big\vert\bar X_N\right]$$ thanks to the Rao-Blackwell theorem. I however do not see an easy way out for computing this conditional expectation (which should be a polynomial of degree $$N$$ in $$\bar X_N$$). But since $$\bar X_N$$ is minimal sufficient and complete there exists a single (UMVU) estimator.

• Thank you! Is it essentially saying that the product can be expressed as the exponential of a $\mathcal{N}(N\bar{\mu},N\bar{\sigma}^2)$ random variable by the Central Limit Theorem? Jun 9 at 1:06
• Thanks. Does it imply that, if $X_i$ are such that we can ignore the possibility that they are smaller than zero, then the random variable $\prod X_i$ has the highest probability around $\prod X_i \approx \text{exp}\{N\tilde{\mu}\}$? Jun 9 at 7:05
• The mode of a log-normal distribution $\mathcal L\mathcal N(\mu,\sigma)$ is $\exp\{\mu-\sigma^2\}$. Jun 9 at 8:21

There are several things to note here.

1. Multiplying lots of numbers in computer leads to rounding errors, especially if these numbers vary in magnitude, i.e. small numbers are multiplied by large numbers.
2. You need to sample the $$\prod X_i$$, not $$X_i$$.
3. In simulations you rely on law of large numbers. However the bigger the variance of the random variable you sample, the slower the convergence.

To see the point 3 let us calculate the variance of $$\prod X_i$$:

$$Var(\prod X_i) = E(\prod X_i)^2 - (E\prod X_i)^2 = (EX_i^2)^N - (EX_i)^N.$$

For normal variable $$N(\mu, \sigma^2)$$ this expression turns into:

$$(\sigma^2+\mu^2)^N-\mu^N \approx N\sigma^2$$

if $$\mu=1$$. So the variance scales with $$N$$.

So here is an example in R which showcases this:

> set.seed(666)
> prod_sample <- sapply(1:1000, function(x)prod(rnorm(100, mean = 1, sd = sqrt(0.2))))
> quantile(prod_sample)
0%           25%           50%           75%          100%
-3.001294e+00 -5.287966e-07  5.889490e-11  2.814698e-06  7.951951e+01
> mean(prod_sample)
[1] 0.1500762


This is for $$N=100$$ and the product is sampled 1000 times.

If we chose $$N=3$$, then the result is more aligned with the expectations:

> prod_sample <- sapply(1:1000, function(x)prod(rnorm(3, mean = 1, sd = sqrt(0.2))))
> quantile(prod_sample)
0%        25%        50%        75%       100%
-0.9823311  0.4047946  0.7918063  1.3744908  5.5923739
>
> mean(prod_sample)
[1] 1.000404

• +1 I like the variance argument: it's simple, straightforward, and completely general. To see why, wlg let $\mu=1$ and the common variances of the $X_i$ be $\sigma^2.$ The variance of the product then is $(1+\sigma^2)^N-1$ which grows exponentially (even faster than your linear lower bound) unless the $X_i$ are constant. Coupled with an argument that most values are very close to zero, this quantifies how large the rare larger values must be.
– whuber
Jun 9 at 12:06
• replicate(1000, prod(rnorm(100, mean = 1, sd = sqrt(0.2)))) strikes me as cleaner. Jun 11 at 15:18