# Validating uncertainty quantification

Regression performance is often evaluated by means of cross-validation. However, classical cross-validation only regards the mean of the identified parameters. How can one quantify the quality of the uncertainty estimate, i.e., the validity of the estimated parameter covariance?

An illustrative example: ordinary least squares (OLS) assumes constant error covariance (homoskedasticity), and performing OLS on heteroskedastic data is known to yield unbiased parameter estimates but biased covariance estimates. One solution is to perform weighted OLS with weights $$w_i = 1/\sigma_i$$. Of course many more combinations of statistical properties of the data and different estimators are conceivable, all leading to different parameter covariance estimates. Thus the question arises in practice: how can I validate and compare the uncertainty estimates provided by different estimation procedures, in order to select the estimator with the most reliable uncertainty estimate?

Example code (Matlab): y = 0.1*x1 - 0.2*x2 + 1

rng(1);
x1 = (0:.001:10)-5;
x2 = detrend(x1-sin(x1));  % generate second, uncorrelated covariate
sigma_eps = linspace(1, 10, length(x1));  % heteroskedastic error variance
y = 0.1*x1 +--0.2*x2 + 1 + sigma_eps.*randn(size(x1));
ols = fitlm([x1; x2; ones(size(x1))]', y', 'Intercept', false)
weights = 1./ sigma_eps;
wls = fitlm([weights.*x1; weights.*x2; weights.*ones(size(x1))]', (weights.*y)', 'Intercept', false)


Output:

ols =

Linear regression model:
y ~ x1 + x2 + x3

Estimated Coefficients:
Estimate        SE        tStat        pValue
_________    ________    ________    __________

x1     0.087713    0.021037      4.1694    3.0802e-05
x2    -0.076851    0.085888    -0.89478       0.37092
x3      0.95685    0.060736      15.754    2.9584e-55

...

wls =

Linear regression model:
y ~ x1 + x2 + x3

Estimated Coefficients:
Estimate       SE        tStat       pValue
________    ________    _______    __________

x1    0.081754    0.017311     4.7227    2.3587e-06
x2    -0.15835    0.059521    -2.6604     0.0078181
x3     0.95645    0.060051     15.927     2.021e-56


The true parameters are within 2*SE of the estimated parameters in all cases, hence nothing appears to be severely off. The aim of the question is to find a method to validate the SE estimates empirically.

To understand how to validate the covariance estimate $$\hat\Sigma(\hat \theta)$$, first make sure to really comprehend the meaning of this quantity. To keep things simple, consider the linear regression model $$y_i = \theta \cdot x_i + \varepsilon_i$$. The measurements $$y_i$$ and $$x_i$$ are realizations of random variables $$X$$ and $$Y$$, and hence the parameter estimate $$\hat \theta$$ will also be a random variable which depends on $$X$$ and $$Y$$. What is often just called "the parameter estimate" corresponds to the mean of this variable, while $$\hat\Sigma(\hat \theta)$$ is exactly that: the estimated variance of the random variable $$\hat \theta$$, which describes the spread of parameter estimates we are likely to observe for different realizations of the measurement random variables.

One admittedly rather impractical way to validate the estimated parameter estimate covariance $$\hat\Sigma(\hat \theta)$$ empirically is then to "simply" draw $$B$$ realizations $$j=1,\ldots,B$$ of $$X$$ and $$Y$$, calculate parameter estimates $$\hat \theta^j$$ using these new realizations, and calculate the sample variance of these estimates $$\mathrm{Var}(\hat\theta) = \frac{1}{B-1}\sum_{j=1}^B (\hat\theta^j-\bar\theta)^2$$ which should be approximately equal to the estimated variance $$\hat\Sigma(\hat\theta)$$ Here's a Matlab implementation of this procedure for the example given in the question:

N_iter = 100;
params_ols = zeros(N_iter, 3);
se_ols = zeros(N_iter, 3);
params_wls = zeros(N_iter, 3);
se_wls = zeros(N_iter, 3);
for ii = 1:N_iter
rng(ii);  % generate different realizations of the same RVs
x1 = (0:.001:10)-5;
x2 = detrend(x1-sin(x1));  % generate second, uncorrelated covariate
sigma_eps = linspace(1, 10, length(x1));
y = 0.1*x1 - 0.2*x2 + 1 + sigma_eps.*randn(size(x1));
ols = fitlm([x1; x2; ones(size(x1))]', y', 'Intercept', false);
params_ols(ii, :) = ols.Coefficients.Estimate';
se_ols(ii, :) = ols.Coefficients.SE';
weights = 1./ sigma_eps;
wls = fitlm([weights.*x1; weights.*x2; weights.*ones(size(x1))]', (weights.*y)', 'Intercept', false);
params_wls(ii, :) = wls.Coefficients.Estimate';
se_wls(ii, :) = wls.Coefficients.SE';
end

actual_se_ols = std(params_ols);
actual_se_wls = std(params_wls);
disp('Actual standard deviation of OLS parameter estimates is:')
disp(actual_se_ols);   % 0.0218    0.0834    0.0677
disp('Actual standard deviation of WLS parameter estimates is:')
disp(actual_se_wls);   % 0.0178    0.0606    0.0671


It turns out that (surprisingly?) in this example, both OLS and WLS provided quite accurate SE estimates.

The drawback of this procedure is of course that it requires a lot of additional data. If one has, say, $$N=1000$$ samples available in total and would like to consider $$B=100$$ independent data sets for empirically estimating the variance of the estimator, then the data sets used for training the estimators $$\hat \theta^j$$ will only consist of $$N/B=10$$ samples each. The natural way to solve this problem is of course to use some kind of cross-validation or bootstrapping, i.e., reusing the same data points multiple times. In this case, however, the datasets used for calculating the estimators $$\hat \theta^j$$ are no independent realizations anymore, and I am unsure how this affects the validity of the procedure / how one would need to adjust the variance estimation. I'd be very happy for someone more knowledgable than me to clarify this issue.

Some relevant pointers: