Confidence Intervals¶
This is a simple demonstration of the algorithms and how well they do with confidence intervals. Each algorithm proposed below are Bayesian in nature and have confidence intervals attached to them. Obviously we won't be able to plot multidimensional confidence intervals so easily but this will at least give us an idea.
Note: This is a heads up but I didn't do any special training procedures to train the algorithms. So the performance may not be optimum.
Below are the Algorithms that have confidence intervals.
- Bayesian Ridge Regression
- ARD Regression
- Gaussian Process Regression
import numpy as np
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, BayesianRidge, ARDRegression
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, ConstantKernel, RBF
from sklearn.multioutput import MultiOutputRegressor
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import time as time
import matplotlib.pyplot as plt
%matplotlib inline
# Reproducibility
rng = np.random.RandomState(0)
# Generate sample data
X = 15 * rng.rand(100, 1)
f = lambda x: np.sin(x)
y = f(np.sin(X).ravel())[:, None]
# print(y.shape)
y += 3 * (0.5 - rng.rand(X.shape[0]))[:, None] # add noise
# Generate plot data
X_plot = np.linspace(0, 20, 10000)[:, None]
Bayesian Ridge Regression¶
# Define GP Model
bayes_model = BayesianRidge(
n_iter=1000,
normalize=True, # good standard practice
verbose=1, # Gives updates
compute_score=True
)
# train GP Model
t0 = time.time()
bayes_model.fit(X, y)
t1 = time.time() - t0
# Predictions
ypred, ystd = bayes_model.predict(X_plot, return_std=True)
fig, ax = plt.subplots(figsize=(10, 5))
plt.scatter(X, y, c='k', label='data')
plt.plot(X_plot, f(X_plot), color='navy', lw=2, label='True')
# Plot standard deviations
plt.fill_between(X_plot[:, 0],
ypred.squeeze() - ystd,
ypred.squeeze() + ystd, color='darkorange',
alpha=0.2)
plt.legend()
plt.show()
The mean looks ok but the variance estimates don't look nearly as good as the GPs below.
ARD Regression¶
# Define GP Model
ard_model = ARDRegression(
n_iter=1000,
verbose=1, # Gives updates
)
# train GP Model
t0 = time.time()
ard_model.fit(X, y)
t1 = time.time() - t0
# Predictions
ypred, ystd = ard_model.predict(X_plot, return_std=True)
fig, ax = plt.subplots(figsize=(10, 5))
plt.scatter(X, y, c='k', label='data')
plt.plot(X_plot, f(X_plot), color='navy', lw=2, label='True')
# Plot standard deviations
plt.fill_between(X_plot[:, 0],
ypred.squeeze() - ystd,
ypred.squeeze() + ystd, color='darkorange',
alpha=0.2)
plt.legend()
plt.show()
Gaussian Process Regression¶
# define kernel function
kernel = ConstantKernel() * RBF() + WhiteKernel()
# define GP model
gp_model = GaussianProcessRegressor(
kernel=kernel, # kernel function (very important)
normalize_y=True, # good standard practice
random_state=123, # reproducibility
n_restarts_optimizer=10, # good practice (avoids local minima)
)
# train GP Model
t0 = time.time()
gp_model.fit(X, y)
t1 = time.time() - t0
# Predictions
ypred, ystd = gp_model.predict(X_plot, return_std=True)
fig, ax = plt.subplots(figsize=(10, 5))
plt.scatter(X, y, c='k', label='data')
plt.plot(X_plot, f(X_plot), color='navy', lw=2, label='True')
# Plot standard deviations
plt.fill_between(X_plot[:, 0],
ypred.squeeze() - ystd,
ypred.squeeze() + ystd, color='darkorange',
alpha=0.2)
plt.legend()
plt.show()
