ageron
commented
5 years ago Hi @nicholashub , Great question!
A confidence interval is computed for the mean of a sample from a random distribution. We don't need to know what this actual random distribution is, thanks to the central limit theorem which states that the sampling distribution of the sample mean of any random distribution (provided its mean is finite and its standard deviation is not zero), will tend towards a Gaussian distribution as the sample size increases. So if we have enough instances in the test set (which is basically always the case, since even just 100 instances is sufficient), then the confidence interval should be fairly meaningful (regardless of the actual error distribution).
In the code I used in the book, the random distribution I'm considering is the distribution of the squared errors. I'm not making any assumption about it (I don't need to, as explained above). The stats.t.interval() code will give a confidence interval for the mean of the the squared errors, i.e., the MSE. As you probably know given your question (but I'm writing this for other people who will read this thread), the interpretation of this interval is not trivial: it does not mean that there is 95% chance that the true MSE (i.e., the generalization error) is within this interval. It means that if we were to sample many test sets and compute the confidence interval using this same procedure, then 95% of the time the true MSE would lie within the computed interval. It's a subtle but important difference (to learn more about this, research Frequentistism vs Bayesianism).
Now, if the true MSE lies within an interval [A, B], then the RMSE must lie within the interval [sqrt(A), sqrt(B)], since the RMSE is just the square root of the MSE. However, although the confidence interval [A, B] is centered on the true MSE, on average (i.e., if we were to sample more and more test sets, and compute the mean of the confidence interval centers, it would approach the true MSE), the interval [sqrt(A), sqrt(B)] is not centered on the true RMSE, on average.
Below is a script that creates a population of 1 million values from an arbitrary distribution, as well as some predictions, also from an arbitrary (but different) distribution. Then it samples a test set of 100 instances from that population and computes the confidence interval using the code from the book and checks whether the true RMSE is inside this interval. Then it samples a new test set and does the same thing. And it keeps doing this 10,000 times, and counts the number of times that the confidence interval contained the true RMSE. It ends up very close to 95%, as expected. I also tried the same thing with the MSE, the MAE and just the mean error, and I always get ~95%.
import numpy as np
from scipy import stats
N = 10**6
y_true_pop = 5 * np.random.randn(N) ** 2 + 3 * np.log(100 * np.random.rand(N))
y_pred_pop = y_true_pop + 5 * np.cos(10 * np.random.randn(N)) + 5 * np.random.randn(N) - 0.5
true_mean_error = np.mean(y_pred_pop - y_true_pop)
true_mae = np.mean(np.abs(y_pred_pop - y_true_pop))
true_mse = np.mean((y_pred_pop - y_true_pop) ** 2)
true_rmse = np.sqrt(true_mse)
true_mean_error = np.mean(y_pred_pop - y_true_pop)
def sample(n_samples):
sample = np.random.randint(N, size=n_samples)
y_true = y_true_pop[sample]
y_pred = y_pred_pop[sample]
return y_true, y_pred
def compute_interval(samples, confidence=0.95):
return stats.t.interval(confidence,
len(samples) - 1,
loc=samples.mean(),
scale=stats.sem(samples))
def true_mse_in_interval(confidence=0.95, n_samples=100):
y_true, y_pred = sample(n_samples)
squared_errors = (y_pred - y_true) ** 2
interval = compute_interval(squared_errors, confidence)
return interval[0] < true_mse < interval[1]
def true_rmse_in_interval(confidence=0.95, n_samples=100):
y_true, y_pred = sample(n_samples)
squared_errors = (y_pred - y_true) ** 2
interval = np.sqrt(compute_interval(squared_errors, confidence))
return interval[0] < true_rmse < interval[1]
def true_mae_in_interval(confidence=0.95, n_samples=100):
y_true, y_pred = sample(n_samples)
absolute_errors = np.abs(y_pred - y_true)
interval = compute_interval(absolute_errors, confidence)
return interval[0] < true_mae < interval[1]
def true_mean_error_in_interval(confidence=0.95, n_samples=100):
y_true, y_pred = sample(n_samples)
errors = y_pred - y_true
interval = compute_interval(errors, confidence)
return interval[0] < true_mean_error < interval[1]
def estimate_proba_in_interval(func, confidence=0.95, n_samples=100,
n_experiments=10000):
n_in_interval = 0
for _ in range(n_experiments):
if func(confidence=confidence, n_samples=n_samples):
n_in_interval += 1
return n_in_interval / n_experiments
print("Estimated proba that the MSE is in interval:",
estimate_proba_in_interval(true_mse_in_interval))
print("Estimated proba that the RMSE is in interval:",
estimate_proba_in_interval(true_rmse_in_interval))
print("Estimated proba that the MAE is in interval:",
estimate_proba_in_interval(true_mae_in_interval))
print("Estimated proba that the Mean Error is in interval:",
estimate_proba_in_interval(true_mean_error_in_interval))The output is:
Estimated proba that the MSE is in interval: 0.9412
Estimated proba that the RMSE is in interval: 0.9369
Estimated proba that the MAE is in interval: 0.9506
Estimated proba that the Mean Error is in interval: 0.9435This is the ratio of test sets where the procedure produced an interval that contained the true metric (over the whole population).
I hope this is all clear?
nicholashub
Lawrenauh
vasili111
In chapter 2, the confidence interval of generalization error (RMSE) is calculated by taking the square root of a t interval. The code is as follows,
It's unintuitive to me, though (isn't the sum of squared errors follow something like a ${\chi}^2$ distribution?). I wonder under what assumptions (such as what distribution of the variables) this result is deduced, thanks!