Show code cell source
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks");
The point-predictive Distribution - Separating Epistemic and Aleatory Uncertainty#
We will demonstrate how we can separate epistemic and aleatory uncertainty.
Example (Quadratic)#
Let’s repeat what we did above with a quadratic example. Here are some synthetic data:
np.random.seed(12345)
num_obs = 10
x = -1.0 + 2 * np.random.rand(num_obs)
w0_true = -0.5
w1_true = 2.0
w2_true = 2.0
sigma_true = 0.1
y = (
w0_true
+ w1_true * x
+ w2_true * x ** 2
+ sigma_true * np.random.randn(num_obs)
)
fig, ax = plt.subplots()
ax.plot(x, y, 'x', label='Observed data')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);
Let’s also copy-paste the code from previous hands-on activities for generating the design matrix and fitting the models.
import scipy
def get_polynomial_design_matrix(x, degree):
"""Return the polynomial design matrix of ``degree`` evaluated at ``x``.
Arguments:
x -- A 2D array with only one column.
degree -- An integer greater than zero.
"""
assert isinstance(x, np.ndarray), 'x is not a numpy array.'
assert x.ndim == 2, 'You must make x a 2D array.'
assert x.shape[1] == 1, 'x must be a column.'
cols = []
for i in range(degree+1):
cols.append(x ** i)
return np.hstack(cols)
def get_fourier_design_matrix(x, L, num_terms):
"""Fourier expansion with ``num_terms`` cosines and sines.
Arguments:
x -- A 2D array with only one column.
L -- The "length" of the domain.
num_terms -- How many Fourier terms do you want.
This is not the number of basis
functions you get. The number of basis functions
is 1 + num_terms / 2. The first one is a constant.
"""
assert isinstance(x, np.ndarray), 'x is not a numpy array.'
assert x.ndim == 2, 'You must make x a 2D array.'
assert x.shape[1] == 1, 'x must be a column.'
N = x.shape[0]
cols = [np.ones((N, 1))]
for i in range(int(num_terms / 2)):
cols.append(np.cos(2 * (i+1) * np.pi / L * x))
cols.append(np.sin(2 * (i+1) * np.pi / L * x))
return np.hstack(cols)
def get_rbf_design_matrix(x, x_centers, ell):
"""Radial basis functions design matrix.
Arguments:
x -- The input points on which you want to evaluate the
design matrix.
x_center -- The centers of the radial basis functions.
ell -- The lengthscale of the radial basis function.
"""
assert isinstance(x, np.ndarray), 'x is not a numpy array.'
assert x.ndim == 2, 'You must make x a 2D array.'
assert x.shape[1] == 1, 'x must be a column.'
N = x.shape[0]
cols = [np.ones((N, 1))]
for i in range(x_centers.shape[0]):
cols.append(np.exp(-(x - x_centers[i]) ** 2 / ell))
return np.hstack(cols)
def find_m_and_S(Phi, y, sigma2, alpha):
"""Return the posterior mean and covariance of the weights
of a Bayesian linear regression problem.
Arguments:
Phi -- The design matrix.
y -- The observed targets.
sigma2 -- The noise variance.
alpha -- The prior weight precision.
"""
A = (
Phi.T @ Phi / sigma2
+ alpha * np.eye(Phi.shape[1])
)
L = scipy.linalg.cho_factor(A)
m = scipy.linalg.cho_solve(
L,
Phi.T @ y / sigma2
)
S = scipy.linalg.cho_solve(
L,
np.eye(Phi.shape[1])
)
return m, S
Fit a
Show code cell source
import scipy.stats as st
# Parameters
degree = 7
sigma2 = 0.1 ** 2
alpha = 5.0
# Weight prior
w_prior = st.multivariate_normal(
mean=np.zeros(degree+1),
cov=np.eye(degree+1) / alpha
)
# Design matrix
Phi = get_polynomial_design_matrix(x[:, None], degree)
# Fit
m, S = find_m_and_S(Phi, y, sigma2, alpha)
# Weight posterior
w_post = st.multivariate_normal(mean=m, cov=S)
As we discussed in the video, it is possible to get the posterior point predictive distribution for
corresponds to the measurement noise. is the epistemic uncertainty induced by limited data.
Here is how to visualize both of these:
Show code cell source
xx = np.linspace(-1, 1, 100)
Phi_xx = get_polynomial_design_matrix(xx[:, None], degree)
# Posterior predictive mean
yy_mean = Phi_xx @ m
# Posterior predictive epistemic variance
yy_var = np.einsum(
'ij,jk,ik->i',
Phi_xx,
S,
Phi_xx
)
# Posterior predictive epistemic + aleatory variance
yy_measured_var = yy_var + sigma2
# 95% posterior predictive credible interval
yy_std = np.sqrt(yy_var)
yy_measured_std = np.sqrt(yy_measured_var)
# Epistemic only
yy_le = yy_mean - 2.0 * yy_std
yy_ue = yy_mean + 2.0 * yy_std
# Epistemic + aleatory
yy_lae = yy_mean - 2.0 * yy_measured_std
yy_uae = yy_mean + 2.0 * yy_measured_std
# The true response for plotting
yy_true = w0_true + w1_true * xx + w2_true * xx ** 2
# Plot
fig, ax = plt.subplots()
ax.plot(xx, yy_mean, 'r', label="Posterior mean")
ax.fill_between(
xx,
yy_le,
yy_ue,
color='red',
alpha=0.25,
label="95% epistemic credible interval"
)
ax.fill_between(
xx,
yy_lae,
yy_le,
color='green',
alpha=0.25
)
ax.fill_between(
xx,
yy_ue,
yy_uae,
color='green',
alpha=0.25,
label="95% epistemic + aleatory credible interval"
)
ax.plot(x, y, 'kx', label='Observed data')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc="best", frameon=False)
sns.despine(trim=True);
Questions#
Rerun the code cells above with a very small
. What happens?Rerun he code cells above with a very big
. What happens?Fix
to and rerun the code cells above with a very small and very big value for . What happens in each case