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

  • 1
    $\begingroup$ Would bootstrapping work? Or do you want an analytical confidence interval? $\endgroup$
    – jld
    Commented Jun 26, 2020 at 16:25
  • $\begingroup$ 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. $\endgroup$
    – Sam
    Commented Jun 26, 2020 at 16:54
  • $\begingroup$ 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. $\endgroup$
    – Sam
    Commented Jun 26, 2020 at 16:56
  • 1
    $\begingroup$ This is essentially the same problem as finding confidence intervals for the mean of a lognormal distribution. $\endgroup$
    – whuber
    Commented Jun 26, 2020 at 17:30

1 Answer 1


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


$$ \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])

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

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:

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.