Expectation of the product of iid random variables

Question

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!

Answer 1

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.$

Answer 2

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.

Answer 3

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