This example builds on the adaptive covariance MCMC example, and shows you a different way to plot the results. The plot shown here is similar to the pairwise scatterplots example, but instead of using histogram plots, we use kernel density estimation (KDE) to represent the pairwise probability distributions of each parameter.
Inference plots:
See the adaptive covariance MCMC example for details.
from __future__ import print_function
import pints
import pints.toy as toy
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
# Load a forward model
model = toy.LogisticModel()
# Create some toy data
real_parameters = [0.015, 500]
times = np.linspace(0, 1000, 100)
org_values = model.simulate(real_parameters, times)
# Add noise
noise = 50
values = org_values + np.random.normal(0, noise, org_values.shape)
real_parameters = np.array(real_parameters + [noise])
# Get properties of the noise sample
noise_sample_mean = np.mean(values - org_values)
noise_sample_std = np.std(values - org_values)
# Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, values)
# Create a log-likelihood function (adds an extra parameter!)
log_likelihood = pints.GaussianLogLikelihood(problem)
# Create a uniform log-prior over both the parameters and the new noise variable
log_prior = pints.UniformLogPrior(
[0.01, 400, noise * 0.1],
[0.02, 600, noise * 100]
)
# Create a posterior log-likelihood (log(likelihood * prior))
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Perform sampling using MCMC, with a single chain
x0 = real_parameters * 1.1
mcmc = pints.MCMCController(log_posterior, 1, [x0])
mcmc.set_max_iterations(6000)
mcmc.set_log_to_screen(False)
We now run the MCMC routine and plot the histogram for a single parameter. In addition, we can use SciPy's Kernel-Density Estimation (KDE) method to obtain a smooth estimate of the probability distribution.
print('Running...')
chains = mcmc.run()
print('Done!')
Running... Done!
# Select chain 0 and discard warm-up
chain = chains[0]
chain = chain[3000:]
def plot_kde_1d(x, ax):
""" Creates a 1d histogram and an estimate of the PDF using KDE. """
xmin = np.min(x)
xmax = np.max(x)
x1 = np.linspace(xmin, xmax, 100)
x2 = np.linspace(xmin, xmax, 50)
kernel = stats.gaussian_kde(x)
f = kernel(x1)
hist = ax.hist(x, bins=x2, normed=True)
ax.plot(x1, f)
fig, ax = plt.subplots(1, 1)
ax.set_xlabel('Parameter 1')
ax.set_ylabel('Probability density')
plot_kde_1d(chain[:, 0], ax)
plt.show()
/home/chonlei/anaconda3/envs/matplotlib_build/lib/python3.7/site-packages/matplotlib/axes/_axes.py:6521: MatplotlibDeprecationWarning: The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead. alternative="'density'", removal="3.1")
We can use KDE again to get very nice plots of 2d distributions.
def plot_kde_2d(x, y, ax):
# Get minimum and maximum values
xmin, xmax = np.min(x), np.max(x)
ymin, ymax = np.min(y), np.max(y)
# Plot values
values = np.vstack([x, y])
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.imshow(np.rot90(values), cmap=plt.cm.Blues, extent=[xmin, xmax, ymin, ymax])
# Create grid
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([xx.ravel(), yy.ravel()])
# Get kernel density estimate and plot contours
kernel = stats.gaussian_kde(values)
f = np.reshape(kernel(positions).T, xx.shape)
ax.contourf(xx, yy, f, cmap='Blues')
ax.contour(xx, yy, f, colors='k')
# Fix aspect ratio
ax.set_aspect((xmax - xmin) / (ymax - ymin))
fig, ax = plt.subplots(1, 1)
ax.set_xlabel('parameter 1')
ax.set_ylabel('parameter 2')
plot_kde_2d(chain[:, 0], chain[:, 1], ax)
plt.show()
Using these functions, we can now create a matrix of scatterplots with KDE.
n_param = log_likelihood.n_parameters()
fig_size = (12, 12)
fig, axes = plt.subplots(n_param, n_param, figsize=fig_size)
for i in range(n_param):
for j in range(n_param):
# Create subplot
if i == j:
# Plot the diagonal
plot_kde_1d(chain[:, i], ax=axes[i, j])
axes[i, j].axvline(real_parameters[i], c='g', label='True value')
axes[i, j].legend()
elif i < j:
# Upper triangle: No plot
axes[i, j].axis('off')
else:
# Lower triangle: Pairwise plot
plot_kde_2d(chain[:, j], chain[:, i], ax=axes[i, j])
axes[i, j].axhline(real_parameters[i], c='g')
axes[i, j].axvline(real_parameters[j], c='g')
# Adjust the tick labels
if i < n_param - 1:
# Only show x tick labels for the last row
axes[i, j].set_xticklabels([])
else:
# Rotation the x tick labels to fit in the plot
for tl in axes[i, j].get_xticklabels():
tl.set_rotation(45)
if j > 0:
# Only show y tick labels for the first column
axes[i, j].set_yticklabels([])
# Add labels to the subplots at the edges
axes[i, 0].set_ylabel('parameter %d' % (i + 1))
axes[-1, i].set_xlabel('parameter %d' % (i + 1))
plt.show()