Show code cell source
MAKE_BOOK_FIGURES=Trueimport numpy as npimport scipy.stats as stimport matplotlib as mplimport matplotlib.pyplot as plt%matplotlib inlineimport matplotlib_inlinematplotlib_inline.backend_inline.set_matplotlib_formats('svg')import seaborn as snssns.set_context("paper")sns.set_style("ticks")def set_book_style(): plt.style.use('seaborn-v0_8-white') sns.set_style("ticks") sns.set_palette("deep") mpl.rcParams.update({ # Font settings 'font.family': 'serif', # For academic publishing 'font.size': 8, # As requested, 10pt font 'axes.labelsize': 8, 'axes.titlesize': 8, 'xtick.labelsize': 7, # Slightly smaller for better readability 'ytick.labelsize': 7, 'legend.fontsize': 7, # Line and marker settings for consistency 'axes.linewidth': 0.5, 'grid.linewidth': 0.5, 'lines.linewidth': 1.0, 'lines.markersize': 4, # Layout to prevent clipped labels 'figure.constrained_layout.use': True, # Default DPI (will override when saving) 'figure.dpi': 600, 'savefig.dpi': 600, # Despine - remove top and right spines 'axes.spines.top': False, 'axes.spines.right': False, # Remove legend frame 'legend.frameon': False, # Additional trim settings 'figure.autolayout': True, # Alternative to constrained_layout 'savefig.bbox': 'tight', # Trim when saving 'savefig.pad_inches': 0.1 # Small padding to ensure nothing gets cut off })def set_notebook_style(): plt.style.use('seaborn-v0_8-white') sns.set_style("ticks") sns.set_palette("deep") mpl.rcParams.update({ # Font settings - using default sizes 'font.family': 'serif', 'axes.labelsize': 10, 'axes.titlesize': 10, 'xtick.labelsize': 9, 'ytick.labelsize': 9, 'legend.fontsize': 9, # Line and marker settings 'axes.linewidth': 0.5, 'grid.linewidth': 0.5, 'lines.linewidth': 1.0, 'lines.markersize': 4, # Layout settings 'figure.constrained_layout.use': True, # Remove only top and right spines 'axes.spines.top': False, 'axes.spines.right': False, # Remove legend frame 'legend.frameon': False, # Additional settings 'figure.autolayout': True, 'savefig.bbox': 'tight', 'savefig.pad_inches': 0.1 })def save_for_book(fig, filename, is_vector=True, **kwargs): """ Save a figure with book-optimized settings. Parameters: ----------- fig : matplotlib figure The figure to save filename : str Filename without extension is_vector : bool If True, saves as vector at 1000 dpi. If False, saves as raster at 600 dpi. **kwargs : dict Additional kwargs to pass to savefig """ # Set appropriate DPI and format based on figure type if is_vector: dpi = 1000 ext = '.pdf' else: dpi = 600 ext = '.tif' # Save the figure with book settings fig.savefig(f"{filename}{ext}", dpi=dpi, **kwargs)def make_full_width_fig(): return plt.subplots(figsize=(4.7, 2.9), constrained_layout=True)def make_half_width_fig(): return plt.subplots(figsize=(2.35, 1.45), constrained_layout=True)if MAKE_BOOK_FIGURES: set_book_style()else: set_notebook_style()make_full_width_fig = make_full_width_fig if MAKE_BOOK_FIGURES else lambda: plt.subplots()make_half_width_fig = make_half_width_fig if MAKE_BOOK_FIGURES else lambda: plt.subplots()
Diagnostics for Posterior Predictive#
We introduce measures that quantify how good the posterior predictive distribution is.
Example (Quadratic)#
We start with our quadratic synthetic example because we know that our model is adequate. You will see how the standarized errors will turn out perfect.
np.random.seed(12345)
num_obs = 20
x = -1.0 + 2 * np.random.rand(num_obs)
w0_true = -0.5
w1_true = 1.0
w2_true = 2.0
sigma_true = 0.1
true_func = lambda x: (
w0_true
+ w1_true * x
+ w2_true * x ** 2
)
observe_func = lambda x: (
true_func(x)
+ sigma_true * np.random.randn(x.shape[0])
)
y = observe_func(x)
num_valid = 50
x_valid = -1.0 + 2 * np.random.rand(num_valid)
y_valid = observe_func(x_valid)
fig, ax = plt.subplots()
ax.plot(x, y, 'x', label='Observed data')
ax.plot(x_valid, y_valid, 'o', label='Validation data')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc='best', frameon=False)
sns.despine(trim=True);
Let’s copy-paste again our previous code.
Note
It is generally a bad idea to copy-paste your code from notebook to notebook.
Copy-pasting can introduce bugs that are very hard to fix.
The best approach would be to put all our code in a text file, say regression.py
and then import everything from the file.
We will try this out in the homework assignment.
Show code cell source
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 plot_posterior_predictive(
model,
xx,
phi_func,
phi_func_args=(),
y_true=None
):
"""Plot the posterior predictive separating
aleatory and espitemic uncertainty.
Arguments:
model -- A trained model.
xx -- The points on which to evaluate
the posterior predictive.
phi_func -- The function to use to compute
the design matrix.
Keyword Arguments:
phi_func_args -- Any arguments passed to the
function that calculates the
design matrix.
y_true -- The true response for plotting.
"""
Phi_xx = phi_func(
xx[:, None],
*phi_func_args
)
yy_mean, yy_measured_std = model.predict(
Phi_xx,
return_std=True
)
sigma = np.sqrt(1.0 / model.alpha_)
yy_std = np.sqrt(yy_measured_std ** 2 - sigma**2)
yy_le = yy_mean - 2.0 * yy_std
yy_ue = yy_mean + 2.0 * yy_std
yy_lae = yy_mean - 2.0 * yy_measured_std
yy_uae = yy_mean + 2.0 * yy_measured_std
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')
if y_true is not None:
ax.plot(xx, y_true, "--", label="True response")
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc="best", frameon=False)
sns.despine(trim=True)
def plot_posterior_samples(
model,
xx,
phi_func,
phi_func_args=(),
num_samples=10,
y_true=None,
nugget=1e-6
):
"""Plot posterior samples from the model.
Arguments:
model -- A trained model.
xx -- The points on which to evaluate
the posterior predictive.
phi_func -- The function to use to compute
the design matrix.
Keyword Arguments:
phi_func_args -- Any arguments passed to the
function that calculates the
design matrix.
num_samples -- The number of samples to take.
y_true -- The true response for plotting.
nugget -- A small number to add the covariance
if it is not positive definite
(numerically).
"""
Phi_xx = phi_func(
xx[:, None],
*phi_func_args
)
m = model.coef_
S = model.sigma_
w_post = st.multivariate_normal(
mean=m,
cov=S + nugget * np.eye(S.shape[0])
)
fig, ax = plt.subplots()
for _ in range(num_samples):
w_sample = w_post.rvs()
yy_sample = Phi_xx @ w_sample
ax.plot(xx, yy_sample, 'r', lw=0.5)
ax.plot([], [], "r", lw=0.5, label="Posterior samples")
ax.plot(x, y, 'kx', label='Observed data')
ax.plot(xx, yy_true, label='True response surface')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
plt.legend(loc="best", frameon=False)
sns.despine(trim=True);
Let’s fit the model using automatic relevance determination:
from sklearn.linear_model import ARDRegression
# Parameters
degree = 4
# Design matrix
Phi = get_polynomial_design_matrix(x[:, None], degree)
# Fit
model = ARDRegression(
fit_intercept=False
).fit(Phi, y)
Let’s visualize the resulting model:
xx = np.linspace(-1, 1, 100)
plot_posterior_predictive(
model,
xx,
get_polynomial_design_matrix,
phi_func_args=(degree,),
y_true=true_func(xx)
)
And now let’s make predictions on the validation data:
Phi_valid = get_polynomial_design_matrix(
x_valid[:, None],
degree
)
y_predict, y_std = model.predict(
Phi_valid,
return_std=True
)
First, let’s do the observations vs predictions plot:
fig, ax = plt.subplots()
ax.plot(y_predict, y_valid, 'o')
yys = np.linspace(
y_valid.min(),
y_valid.max(),
100)
ax.plot(yys, yys, 'r-')
ax.set_xlabel('Predictions')
ax.set_ylabel('Observations')
sns.despine(trim=True);
It’s okay. However, notice that the predictions do not fall on the red line because there is a little noise in the observations.
Let’s now compute the standardized errors:
eps = (y_valid - y_predict) / y_std
Remember that the standardized errors must follow a standard normal if the model is correct. There are various ways to check this. First, let’s plot them:
idx = np.arange(1, eps.shape[0] + 1)
fig, ax = plt.subplots()
ax.plot(idx, eps, 'o', label='Standarized errors')
ax.plot(idx, 1.96 * np.ones(eps.shape[0]), 'r--')
ax.plot(idx, -1.96 * np.ones(eps.shape[0]), 'r--')
ax.set_xlabel('$i$')
ax.set_ylabel('$\epsilon_i$')
sns.despine(trim=True);
Notice that most standardized errors fall within the 95% central credible interval for \(N(0,1)\), indicating that the model is good.
The other plot we can do is the histogram of the standardized errors compared to the probability density of the standard normal:
What if you wanted to get the samples from the posterior? You would have to do a little bit of manual work to translate the posterior weights and their variance back to the original values…
import scipy.stats as st
fig, ax = plt.subplots()
ax.hist(eps, alpha=0.5, density=True)
ee = np.linspace(eps.min(), eps.max(), 100)
ax.plot(ee, st.norm.pdf(ee), 'r')
sns.despine(trim=True);
It’s not perfect, but it’s a good fit.
The final diagnostic is the so-called quantile-quantile plot or q-q plot. The q-q plot compares the empirical quantiles of the standardized errors (computed by building and inverting the empirical cumulative distribution function; see Lecture 9). Here is how to do it:
fig, ax = plt.subplots()
st.probplot(eps, dist=st.norm, plot=ax)
sns.despine(trim=True);
The q-q plot above also shows a good fit.
Questions#
Rerun the code blocks above with many training observations, say 100. Did the diagnostics improve, or are they the same?
Keep the number of training observations to 100 and change the validation points to 1000. What do the diagnostics look like now?
Let’s now try a model that is inadequate for the data. Set the polynomial degree to 1 and rerun everything. What do the diagnostics look like now?