Proof for MSE = Var + Bias2

Question

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

Answer 1

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");
    }
}

Answer 2

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?

Answer 3

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