# Confidence of Mean Squared Error

Let $$X \sim \mathcal{N}(\mu, \sigma^2)$$, where $$\mu$$ and $$\sigma$$ are unknown. I would like to estimate

$$\mathbb{E}\left[(X - a)^2\right]$$

where $$a$$ it is a known constant. For this purpose, I have $$n$$ realizations of $$X$$, $$x_1, \dots, x_n$$, and my estimator is

$$e = n^{-1}\sum_{i=1}^n (x_i - a)^2.$$

I would like now to bound my estimator using confidence intervals. I can recognize that

$$e = n^{-1}\sum_{i=1}^n (x_i - a)^2 = n^{-1}\sum_{i=1}^n (\sigma z_i + \mu - a)^2 = n^{-1}\sum_{i=1}^n (\sigma z_i + \mu - a)^2 \\ = n^{-1}\sum_{i=1}^n \sigma^2 z_i^2 + 2\sigma z_i \mu - 2a\sigma z_i - 2a\mu + a^2 + \mu^2 \\ \stackrel{D}{=} \left(\sigma^2\frac{\chi^2_n}{n}\right) + 2\sigma \frac{\sum_iz_i}{n} (\mu-a) +(\mu - a)^2$$ where $$\chi_n^2$$ is a Chi-Squared with $$p=n$$. My confusion comes from the fact that, while I am able to compute the confidence interval of a chi-squared or of a normal distributed variable, here we have a sum over the two variable. Moreover, to be precise, the two variables are not independent, since they have been generated using the same dataset.

How can I compute the confidence interval of my estimator $$e$$ given my observations $$x_1, \dots, x_n$$?

• Would bootstrapping work? Or do you want an analytical confidence interval?
– jld
Jun 26, 2020 at 16:25
• Bootstrap would work (or even better, sampling M different datasets), however, I don't like the assumption of usual C.I., which is that data are normally distributed. Also because It can happened that my C.I. are lower than zero, which is clearly an underestimate (MSE cannot be lower than zero). I would like to use a similar C.I. to the one it is used for the empirical variance.
– Sam
Jun 26, 2020 at 16:54
• Therefore, yes, would be nice to derive an analytical solution (but I fear it is not possible); otherwise some reasonable bounds of the usual C.I. If also this last solution will not be possible, I'll just use classic C.I. with Normal distribution assumption.
– Sam
Jun 26, 2020 at 16:56
• This is essentially the same problem as finding confidence intervals for the mean of a lognormal distribution.
– whuber
Jun 26, 2020 at 17:30

I found this solution. I don't know if it is statistically sound, but on different numerical trials, it seems to work well (the bound holds and it is tight).

Let's reduce the problem to the following:

we want to estimate $$l_c, u_c$$ such that

$$p(\hat{z}_n > z - l_c) \leq \gamma \\ p(\hat{z}_n < z + u_l) \leq 1 - \gamma$$

where

$$\hat{z}_n = n^{-1}\sum_{i=1}^nx_i^2 \\ x_i \sim \mathcal{N}(\mu, \sigma^2) \\ z = \mu^2 + \sigma^2.$$

We notice that

$$\hat{z}_n \stackrel{D}{=} n^{-1}\sum_{i=1}^n \sigma^2\left(z_i^2 + 2\frac{\mu}{\sigma}z_i + \frac{z_i}{\sigma_i}^2\right) \\ = \frac{\sigma^2}{n} \sum_{i=1}^n \left(\mu + z_i\right)^2\\ \stackrel{D}{=} \frac{\sigma^2}{n} \chi_{n, n\mu^2/\sigma^2}^2$$

where $$\chi_{n, \mu/\sigma}$$ it is an noncentral chi-squared with parameters $$k=n, \lambda=n\mu^2/\sigma^2$$.

at this point, we have access to the pdf, the cdf, and the ppf via known numerical heuristic. An implementation of the noncentered chi-squared is on scipy.

The parameters $$\mu$$ and $$\sigma$$ can be estimated in the usual way, since we assume $$x_i$$ to be normal distributed.

I attach a spinnet of the program to estimate the confidence intervals at $$\gamma=0.05$$:

import numpy as np
from scipy.stats import ncx2
import matplotlib.pyplot as plt

mu = 2.
sigma = 2.

ground_truth_z = mu**2 + sigma**2

n = 1000
support = np.arange(1, n+1)

# Values
x = np.random.normal(mu, sigma, size=n)

# Online estimate
x_cum_d = np.cumsum(x**2)/support

def estimate_interval(x):
estimated_sigma = np.std(x)
estimated_mu = np.mean(x)
k = x.shape[0]
mu_chi_2 = k*(estimated_mu/estimated_sigma)**2

l_ci = estimated_sigma**2 * ncx2.ppf(0.05, k, mu_chi_2)/k
u_ci = estimated_sigma**2 * ncx2.ppf(0.95, k, mu_chi_2)/k

return mu_chi_2, l_ci, u_ci

lower_interval = []
upper_interval = []
for i in range(1, n+1):
m, li, ui = estimate_interval(x[:i])
lower_interval.append(li)
upper_interval.append(ui)

plt.plot(support, x_cum_d, label="Online Estimate")
plt.fill_between(support, lower_interval, upper_interval, alpha=0.5, label="Confidence Interval")
plt.hlines(ground_truth_z, 0, n+1, label="Ground truth")
plt.legend(loc='best')
plt.show()


Since $$\mu$$ and $$\sigma$$ are only estimated, I fancy that the bound might be less correct for small $$n$$. However, I think this bound should still be unbiased.

An example of the estimation