MSE & Variance of Monte Carlo Estimator

Question

I am trying to prove the equality of $$\rm MSE(\langle I\rangle)=Var(\langle I \rangle)+Bias(\langle I \rangle)^2$$ but obviously I got something wrong as they don't equal in my calculation:

So here is the example. I use monte carlo to estimate this integral:

$$I = \int_0^1 5x^4~\mathrm dx$$

The value of this integral is 1. Assuming samples are computed from a uniform probability distribution my estimator is:

$$\langle I\rangle = \frac1N\sum _1^N 5x_i^4$$

And the variance of the estimator can be analytically computed as:$$ \textrm{Var}(\langle I\rangle )= \frac1N\int _0^1(5x^4-1) ^2~\mathrm dx=\frac{16}{9N}$$

So Var here is the variance of the estimator, is that right? where, up to each iteration I calculate it as:

$$ \textrm{Var}(\langle I\rangle )= (\mathrm E[\langle I\rangle ^2] - \mathrm E[\langle I\rangle]^2) $$

and in the code:

float Ie = sum / (i + 1); //estimator
float avg2 = sum2 / (i + 1);
float var = avg2 - (Ie * Ie);
var /= i + 1;

and the bias is simply:

float bias2 = (trueValue - Ie);
bias2 *= bias2;

Here is the full code:

#include <iostream>
const int nsamp = 100;
int main()
{
    float trueValue = 1;
    float data[nsamp];
    float sum = 0;
    float sum2=0;
    float mse = 0;
    float sqErr=0;
for (int i = 0; i &lt; nsamp; i++)
{
    float x = rand1();
    data[i] = 5 * x*x*x*x;
    sum += data[i];
    sum2 += data[i] * data[i]; 

    float Ie = sum / (i + 1); //estimator
    float avg2 = sum2 / (i + 1);
    float var = avg2 - (Ie * Ie);
    var /= i + 1;

    float bias2 = (trueValue - Ie);
    bias2 *= bias2;

    sqErr += (trueValue - Ie) * (trueValue - Ie);
    mse = sqErr / (i+1);

    printf(&quot;\nI=%f Var=%f Bias2=%f MSE=%f&quot;, Ie, var, bias2, mse);
}

}

And the output where the mse doesn't equal var+bias2:

I=0.000000 Var=0.000000 Bias2=1.000000 MSE=1.000000
I=0.252220 Var=0.031807 Bias2=0.559176 MSE=0.779588
I=0.170473 Var=0.018592 Bias2=0.688114 MSE=0.749097
I=0.662600 Var=0.192099 Bias2=0.113839 MSE=0.590282
I=0.647206 Var=0.123133 Bias2=0.124464 MSE=0.497118
I=0.583528 Var=0.088888 Bias2=0.173449 MSE=0.443174
I=0.510921 Var=0.069824 Bias2=0.239198 MSE=0.414034
.
.
.
I=0.984586 Var=0.001662 Bias2=0.000238 MSE=0.005417
I=0.983600 Var=0.001659 Bias2=0.000269 MSE=0.005411
I=0.982616 Var=0.001657 Bias2=0.000302 MSE=0.005406
I=0.985189 Var=0.001660 Bias2=0.000219 MSE=0.005401
I=0.984248 Var=0.001658 Bias2=0.000248 MSE=0.005396
I=0.983362 Var=0.001655 Bias2=0.000277 MSE=0.005391

Answer 1

Here's a proof of the equality, it doesn't require choosing a specific function or numerical computations:

\begin{align} E[(I_N - I)^2] &= E[(I_N-E[I_N]+E[I_N]-I)^2] \\ &= E[(I_N-E[I_N])^2 + (E[I_N]-I)^2 + 2(I_N-E[I_N])(E[I_N]-I)] \\ &= E[(I_N-E[I_N])^2] + E[(E[I_N]-I)^2] + E[2(I_N-E[I_N])(E[I_N]-I)] \\ &= Var[I_N] + (E[I_N]-I)^2 + 2(E[I_N]-E[I_N])(E[I_N]-I) \\ &= Var[I_N] + (Bias(I_N))^2 + 0 \end{align}

The 3rd line holds by linearity of $E$, the 4th holds because $E[I_N]-I$ is not a random variable. The last holds since by definition $Bias(I_N) = E[I_N] - I$.

Edit:

In your code you are not computing the variance or the bias, you are computing estimates of those. By definition the bias, variance, and MSE are:

\begin{equation} Bias(\langle I\rangle, I) = E[\langle I \rangle] - I \\ Var(\langle I\rangle) = E[(\langle I \rangle - E[\langle I\rangle])^2] = E[(\langle I \rangle)^2] - (E[\langle I \rangle])^2 \\ MSE(\langle I \rangle, I) = E[(\langle I \rangle - I)^2] = Var[\langle I \rangle] + (Bias(\langle I \rangle, I))^2 \end{equation}

Now let's look at what you are computing in your code. You obviously know the true solution $I=1$. However in code you only compute $\langle I \rangle = \frac{1}{N}\sum_{k=1}^N g(x_k)$ and not $E[\langle I \rangle]$. This means that you don't compute the squared bias $(E[\langle I \rangle] - I)^2$ but instead you computed $(\langle I \rangle - I)^2$. That alone already introduced approximation error. Then for an estimate of the variance you have used $\frac{1}{N}\left(\frac{1}{N}\sum_{k=1}^N (g(x))^2 - (\frac{1}{N}\sum_{k=1}^N g(x))^2\right)$, but note that this is not equal to the variance, it's an estimate of it. And it's even a biased estimate (you need Bessel's correction for example for an unbiased one).

Finally let's address the elephant in the room: I don't think what you are doing for the MSE estimate makes any sense:

\begin{equation} mse = \frac{1}{N}\sum_{k=1}^N(I_k-I)^2 \end{equation}

In fact I am fairly certain that the above "estimator" isn't even consistent (i.e. it doesn't even estimate the MSE), since if I make a large error at the beginning then it will stick around forever.

Either way it should be clear now that you never computed the expectations in the first place, and just approximations of those, and the estimator for the mse that you wrote is wrong.