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$?
 A: 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(\mu^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 $z_i\sim \mathcal{N}(0,1)$ and $\chi^2_{n, \mu/\sigma}$ is a noncentral chi-squared distribution 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 heuristics. An implementation of the noncentral chi-squared is on scipy.
The parameters $\mu$ and $\sigma$ can be estimated in the usual way, since we assume $x_i$ to be normally distributed.
I attach a snippet 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:

