Proof for MSE = Var + Bias2

I am trying to prove the equality of $$\rm MSE=Var+Bias^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 < 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("\nI=%f Var=%f Bias2=%f MSE=%f", 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


Unless I've misunderstood your code, you appear to be using Ie as both your estimator and as the mean of your estimates.

Ie here should be the mean:

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


But Ie here should be the estimate:

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


If you change your code to something like this you should see the result you expected (apologies for any syntax errors, this is transcribed from the Python code I originally used):

#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;
float ts = 0;

for i in range(nsamp):
float x = rand1();
x = 5 * x*x*x*x;
sum_ += x;
Ie = sum_ / (i + 1);
data[i] = Ie;

float ts += Ie;
float sum2 += Ie * Ie;
float avg = ts / (i + 1);
float var = sum2 / (i + 1) - avg ** 2;

float bias2 = trueValue - avg;
bias2 *= bias2

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

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


I think you should be careful with computing the variance. As far as I understood the situation you consider the estimator

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

and, hence, the variance is

$$\textrm{Var}(\langle I\rangle ) = \int _0^1 \ldots\int _0^1 \left(\frac1N\sum _1^N 5x_i^4 - 1\right)^2 \,\textrm{d}x_1 \ldots \,\textrm{d}x_N$$

which is different from your result, isn't it?

If you want to proof your equality for an estimator (not just checking it for a particular example), you can do this like sketched here: With the formula of the MSE for an estimator of the mean, where $$\bar x$$ denotes the true mean ($$\bar x=\frac 1 N \sum_{i=1}^N x_i$$) and $$\hat x$$ the estimated value for the mean.

$$MSE = \frac 1 N \sum_{i=1}^N (x_i - \hat x)^2$$

Now you can rewrite $$\hat x=\bar x + b$$ where b is the bias and get:

$$MSE = \frac 1 N \sum_{i=1}^N (x_i - \bar x+b)^2$$

That gives us:

$$MSE = \frac 1 N \sum_{i=1}^N (x_i^2 - 2x_i \bar x - 2x_ib + {\bar x}^2 + 2{\bar x}b + b^2)$$

Which we can just rewrite by just reordering the terms as:

$$MSE = \frac 1 N \sum_{i=1}^N (x_i^2 - 2x_i{\bar x} + {\bar x}^2 - 2*x_ib + 2{\bar x}b + b^2)$$

So the first 3 terms are just equivalent to the variance of the $$x_i$$:

$$MSE = Var + \frac 1 N \sum_{i=1}^N (-2x_ib + 2{\bar x}b + b^2)$$

According to the definition of the mean, we know that

$$\frac 1 N \sum_{i=1}^N ({\bar x}b-x_ib) = b(\bar x - \frac 1 N \sum_{i=1}^N x_i)=0$$

So we get what we wanted to show:

$$MSE = Var + b^2$$

If you want to show that for linear regression, you just need to replace the formulas above and e.g. use $$y$$ as a function of $$x_i$$.