This example builds on the adaptive covariance MCMC example, and shows you a different way to plot the results.
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 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 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 can now run the MCMC routine and plot the histograms of the inferred parameters.
print('Running...')
chains = mcmc.run()
print('Done!')
Running... Done!
# Select chain 0 and discard warm-up
chain = chains[0]
chain = chain[3000:]
plt.figure()
plt.xlabel('Parameter 1')
plt.ylabel('Frequency')
bins = np.linspace(np.min(chain[:,0]), np.max(chain[:,0]), 50)
plt.hist(chain[:,0], bins=bins)
plt.show()
Plotting a histogram of two variables (showing their correlation) is a bit more involved. We can use Matplotlib's hist2d method, but we need to do some extra work to get a fixed aspect ratio.
fig, ax = plt.subplots(1, 1)
ax.set_xlabel('parameter 1')
ax.set_ylabel('parameter 2')
xbins = np.linspace(np.min(chain[:,0]), np.max(chain[:,0]), 25)
ybins = np.linspace(np.min(chain[:,1]), np.max(chain[:,1]), 25)
histplot = ax.hist2d(chain[:,0], chain[:,1], bins=(xbins, ybins), cmap=plt.cm.Blues)
# Fix the aspect ratio
ax.set_aspect((xbins[-1] - xbins[0]) / (ybins[-1] - ybins[0]))
plt.show()
We now have all the ingredients to create a matrix of parameter distribution plots.
# Create figure with several subplots
n_bins = 25
n_param = log_likelihood.n_parameters()
fig, axes = plt.subplots(n_param, n_param, figsize=(12, 12))
for i in range(n_param):
for j in range(n_param):
ax = axes[i, j]
# Create subplot
if i == j:
# Diagonal: Plot a 1d histogram
bins = np.linspace(np.min(chain[:,i]), np.max(chain[:,i]), n_bins)
ax.hist(chain[:,i], bins=bins, normed=True)
ax.axvline(real_parameters[i], c='g', label='True value')
ax.legend()
elif i < j:
# Upper right: No plot
ax.axis('off')
else:
# Lower left: Plot a 2d histogram
xbins = np.linspace(np.min(chain[:,j]), np.max(chain[:,j]), n_bins)
ybins = np.linspace(np.min(chain[:,i]), np.max(chain[:,i]), n_bins)
histplot = ax.hist2d(chain[:,j], chain[:,i], bins=(xbins, ybins), normed=True, cmap=plt.cm.Blues)
# Fix the aspect ratio
ax.set_aspect((xbins[-1] - xbins[0]) / (ybins[-1] - ybins[0]))
ax.axhline(real_parameters[i], c='g')
ax.axvline(real_parameters[j], c='g')
# Customise tick labels
if j > 0:
# Only show y tick labels for the first column
axes[i,j].set_yticklabels([])
if i < n_param-1:
# Only show x tick labels for the last row
ax.set_xticklabels([])
else:
# Rotate the x tick labels to fit in the plot
for tl in axes[i,j].get_xticklabels():
tl.set_rotation(45)
# Add labels for subplots at the edges
axes[i,0].set_ylabel('parameter %d'%(i+1))
axes[-1,i].set_xlabel('parameter %d'%(i+1))
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")