Distribution of errors

Question

I am struggeling with a basic question and would be happy to get some pointers. I am trying to evaluate an algorithm $f$ which maps some sample $x$ onto a scalar $y\in\mathbb{R}$, i.e. $f(x) = y$. For each each sample I have a ground-truth version of this scalar $\hat{y}$ available and I want to estimate the performance of the algorithm with respect to $z=y-\hat{y}$. Preliminary evaluations have shown that the algorithm is largely unbiased with the mean of $z$ being close to zero and errors are distributed normally with $z \sim \mathcal{N}(\varepsilon, \sigma)$ for $\varepsilon$ close to zero.

However, I am struggling how specifically to evaluate the performance, preferably with an error bound as e.g. a confidence interval. I could

1. Work with the signed errors $z$. However, in this case I only see how to do a static evaluation of the accuracy ("95% between thresholds ...").
2. Work with the unsigned errors $\vert z \vert$. In this case I could estimate the mean absolute error. However, as I understood I would need to work with folded or a half-normal distributions. Here I am struggling to find references on how to compute confidence intervals (apart from bootstrapping), making me question if this is indeed the canonical way to deal with my problem.

What is the recommended procedure to tackle this problem? As mentioned, I would prefer if this method would allow me to bound some unsigned version of the errors (squared, absolute values, ...). Additionally, I want to estimate the necessary sample size to bound the error in a confidence interval of prescribed width.

Answer 1

You say you have sufficient proof that $Z$ is normal $\mathcal{N}(\varepsilon,\sigma^2)$.

Evaluate the goodness of fit

You can use the log-likelihood of the data $\mathcal{D}= \{z_1, z_2, \dots, z_n\}$ as a measure of $f$'s quality

$$\mathcal{L}(\mathcal{D}) = \sum_{i=1}^n\log p(z_i\mid \varepsilon, \sigma)$$ Where $p(z\mid \varepsilon, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(z-\varepsilon)^2}{2\sigma^2})$ the normal pdf.

This would be almost the same as computing the sum of squared residuals.

Confidence region

A confidence region at level $\alpha$ is defined by $P(Z\in C_\alpha) = \alpha$.

Signed errors (Normal case)

In the normal case you can use the normal CDF $\Phi$ as follows: if $Z$ is $\mathcal{N}(\varepsilon,\sigma^2)$ then $\tilde{Z}=(Z-\varepsilon)/\sigma$ is standard normal. A $95\%$ confidence interval can be found by inverting the relationship $$ P(-\delta \leq \tilde{Z} \leq \delta) = 2P(\tilde{Z}\leq \delta) = 2\Phi(\delta) = \alpha$$ because the normal pdf is symmetric about $0$. By inverting, you get $$ \delta = \Phi^{-1}(\alpha/2)$$ This can be done in python3 using scipy.stats.norm.ppf (which is $\Phi^{-1}$).

Once you have $\delta$, your confidence region is simply $$-\delta \leq \tilde{Z}\leq \delta \iff \varepsilon -\sigma\delta \leq Z \leq \varepsilon + \sigma\delta$$

Unsigned errors (FoldedNormal case)

Use the same definition for a confidence region, I'll assume you want an interval of the form $[0,\delta]$. Normalize again to get $\lvert\tilde{Z}\rvert$. Then the relationship is $$ P(0 \leq \lvert\tilde{Z}\rvert \leq \delta) = F_{\mathrm{FN}}(\delta) = \alpha$$ where $F_{\mathrm{FN}}$ is the FoldedNormal CDF. Invert this to get $$ \delta = F_{\mathrm{FN}}^{-1}(\alpha)$$ In order to compute $\delta$ from $\alpha=0.95$ you can use any implementation of the FoldedNormal inverse cdf such as scipy.stats.foldnorm.ppf.

Answer 2

This sounds like a nonlinear regression problem. Your notation for the problem essentially reverses many of the standard statistical conventions (e.g., you are using a hat for the true value rather than the estimator), so I will show you how to write it out in the standard way, and then measure the performance of your method.

Given $n$ sample pairs you can write your problem using the standard model form:

$$y_i = f(x_i) + \varepsilon_i \quad \quad \quad \varepsilon_1,...,\varepsilon_n \sim \text{IID N}(0, \sigma^2).$$

In this model form the function $f$ is the true regression function (that you are trying to estimate) and the error terms $\varepsilon_1,...,\varepsilon_n$ represent deviations of the observed response values from their true expected value. You have developed a method to estimate the function, so let's denote your estimator as $\hat{f}$ and the resulting response estimates and residuals as:

$$\hat{y}_i = \hat{f}(x_i) \quad \quad \quad r_i = y_i - \hat{y}_i.$$

This is the standard notation for the problem, which is essentially the reverse of the notation you are using in your question. Note that the "hats" on the variables refer to the estimates of the underlying true quantities, which are denoted without hats. Consequently, the direction of the residuals is positive when the true value is higher than the estimated value. I strongly recommend that you adopt the standard notational convention, so that you do not confuse people with your analysis.

---

Now that we have framed your problem in standard notation, let's get to the substance of your question. Since this is a nonlinear regression problem, you can use standard measures of "goodness-of-fit" to evaluate the fit of your model. I will not go into detail here, since there are entire books written on regression goodness-of-fit, but I will note that many of the standard statistics make use of the residual sum-of-squares:

$$SS_\text{Res} = \sum_{i=1}^n (y_i - \hat{y}_i)^2.$$

A lower value for the residual sum-of-squares constitutes a better "fit" for the model and the estimation method. If your function estimator $\hat{f}$ is highly accurate then it will tend to give response predictions that are nearer to the true value, which will tend to reduce the magnitude of the residuals and reduce the residual sum-of-squares.

Since there are existing estimation methods for nonlinear regression problems, one obvious thing to do here would be to compare your algorithm with other estimation methods and see how they perform against one another (e.g., as measured by the residual sum-of-squares or some other goodness-of-fit statistic). Ideally, to avoid overfitting you should use a set of training data to train each method and then use a separate set of test data to compare their performance.

Answer 3

In my final solution for my problem I closely followed the strategies detailed in [1]. As described in my question, my targets were two-fold:

1. Calculate a confidence-interval
2. Estimate the necessary sample size to determine a confidence interval of given width.

Among the strategies outlined in [1], the simplest one is based on the central limit theorem which states that for samples $\{x_1, \dots, x_n\}$ from an arbitrary distribution with mean $\Psi$ and variance $\Theta^2$ it holds that $$ \frac{\overline{x} - \Psi}{\Theta/\sqrt{n}} \overset{n\to\infty}{\rightarrow} \mathcal{N}(0, 1), $$ where $\overline{x} := 1/n\sum_{i=0}^n x_i$. If $n$ is large enough, a confidence-interval for the mean of any distribution can thus be created by $$ \left(\overline{x} - z\frac{\sigma(x)}{\sqrt{n}}, \; \overline{x} + z\frac{\sigma(x)}{\sqrt{n}}\right) $$ since sample mean and standard deviation $\overline{x}$ and $\sigma(x)$ are unbiased estimators of $\Psi$ and $\Theta$ respectively. Since this expression additionally easily allows to estimate the sample size $n$ for a prescribed confidence-interval width, this method was theoretically perfectly suitable for my purposes. However, the question left was: Is $n$ large in my application?

To check this, I wrote a small script to estimate the coverage of the so-determined confidence interval for various parameter combinations $(\Psi,\Theta,n)$ in the vicinity of the expected parameters of my application. Turns out, coverage for the 95% confidence interval is between 93% and 97%, which is completely sufficient for what I wanted to achieve.

[1] Pek, Jolynn & Wong, Augustine & Wong, Octavia. (2017). Confidence Intervals for the Mean of Non-Normal Distribution: Transform or Not to Transform. Open Journal of Statistics. 07. 405-421. 10.4236/ojs.2017.73029.