This example shows you how to perform Bayesian inference on a time series, using Adaptive Covariance MCMC.
It follows on from Optimisation: First example
Like in the optimisation example, we start by importing pints:
In [1]:
import pints
Next, we create a model class.
Instead of using a real model, in this example we use the "Logistic" toy model included in pints:
In [2]:
import pints.toy as toy
model = toy.LogisticModel()
In order to generate some test data, we choose an arbitrary set of "true" parameters:
In [3]:
true_parameters = [0.015, 500]
And a number of time points at which to sample:
In [4]:
import numpy as np
times = np.linspace(0, 1000, 400)
Using these parameters and time points, we can now generate some toy data:
In [5]:
org_values = model.simulate(true_parameters, times)
And make it more realistic by adding gaussian noise:
In [6]:
noise = 25
values = org_values + np.random.normal(0, noise, org_values.shape)
We can use matplotlib (or any other plotting package) to look at the data we've created:
In [7]:
import matplotlib.pyplot as plt
plt.figure(figsize=(12,4.5))
plt.xlabel('Time')
plt.ylabel('Values')
plt.plot(times, values, label='Noisy data')
plt.plot(times, org_values, lw=2, label='Noise-free data')
plt.legend()
plt.show()
Now we have enough data (a model, a list of times, and a list of data) to formulate a problem:
In [8]:
problem = pints.SingleOutputProblem(model, times, values)
We now have some toy data, and a model that can be used for forward simulations. To make it into a probabilistic problem, we need to add a noise model. One way to do this is using the GaussianLogLikelihood function, which assumes independently distributed Gaussian noise over the data, and can calculate log-likelihoods:
In [9]:
log_likelihood = pints.GaussianLogLikelihood(problem)
This noise has mean zero, and an unknown standard deviation. How can we find out the standard deviation? By inferring it along with the other parameters. This means we have added one parameter to our problem!
In [10]:
print('Original problem dimension: ' + str(problem.n_parameters()))
In [11]:
print('New dimension: ' + str(log_likelihood.n_parameters()))
(This means we also have to update our vector of true parameters)
In [12]:
true_parameters += [noise]
print(true_parameters)
This log_likelihood represents the conditional probability $p(y|\theta)$, given a set of parameters $\theta$ and a series of $y=$ values, it can calculate the probability of finding those values if the real parameters are $\theta$.
We can use this in a Bayesian inference scheme to find the quantity we're interested in:
$p(\theta|y) = \frac{p(\theta)p(y|\theta)}{p(y)} \propto p(\theta)p(y|\theta)$
To solve this, we now define a prior, indicating our initial ideas about what the parameters should be. Just as we're using a log-likelihood (the natural logarithm of a likelihood), we'll define this using a log-prior. This simplifies the above equation to:
$\log p(\theta|y) \propto \log p(\theta) + \log p(y|\theta)$
In this example we'll assume we don't know too much about the prior except lower and upper bounds for each variable: We assume the first model parameter is somewhere on the interval $[0.01, 0.02]$, the second model parameter on $[400, 600]$, and the standard deviation of the noise is somewhere on $[1, 100]$.
In [13]:
log_prior = pints.UniformLogPrior(
[0.01, 400, 1],
[0.02, 600, 100]
)
With this prior, we can now define the numerator of Bayes' rule -- the unnormalised log posterior, $\log \left[ p(y|\theta) p(\theta) \right]$:
In [14]:
# Create a posterior log-likelihood (log(likelihood * prior))
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
Finally we create a list of guesses to use as initial positions. We'll run three MCMC chains so we create three initial positions:
In [15]:
xs = [
np.array(true_parameters) * 0.9,
np.array(true_parameters) * 1.05,
np.array(true_parameters) * 1.15,
]
And this gives us everything we need to run an MCMC routine:
In [16]:
chains = pints.mcmc_sample(log_posterior, 3, xs)
We can take a further look at the obtained results using Pints's diagnostic plots.
First, we use the trace method to see if the three chains converged to the same solution.
In [17]:
import pints.plot
pints.plot.trace(chains)
plt.show()
Based on this plot, it looks like the three chains become very similar after about 1000 iterations. To be safe, we throw away the first 2000 samples and continue our analysis with the first chain.
In [18]:
chain = chains[0]
chain = chain[2000:]
We can also look for autocorrelation in the chains, using the autocorrelation() method. If everything went well, the samples in the chain should be relatively independent, so the autocorrelation should get quite low when the lag on the x-axis increases.
In [19]:
pints.plot.autocorrelation(chain)
plt.show()
Now we can inspect the inferred distribution by plotting histograms:
In [20]:
fig, axes = pints.plot.histogram([chain], ref_parameters=true_parameters)
# Show where the sample standard deviation of the generated noise is:
noise_sample_std = np.std(values - org_values)
axes[-1].axvline(noise_sample_std, color='orange', label='Sample standard deviation of noise')
axes[-1].legend()
fig.set_size_inches(14, 9)
plt.show()
Here we've analysed each parameter in isolation, but we can also look at correlations between parameters we found using the pairwise() plot.
To speed things up, we'll first apply some thinning to the chain:
In [21]:
chain = chain[::10]
In [22]:
pints.plot.pairwise(chain, kde=True, ref_parameters=true_parameters)
plt.show()
As these plots show, we came pretty close to the original "true" values (represented by the blue line). But not exactly... Worse, the method seems to suggest a normal distribution but around the wrong point. To find out what's going on, we can plot the log-posterior function near the true parameters:
In [23]:
# Plot log-posterior function
fig, axes = pints.plot.function(log_posterior, true_parameters)
# Add a line showing the sample standard deviation of the generated noise
axes[-1].axvline(noise_sample_std, color='orange', label='Sample standard deviation of noise')
axes[-1].legend()
# Customise the figure size
fig.set_size_inches(14, 9)
plt.show()
As this plot (created entirely without MCMC!) shows, the MCMC method did well, but our estimate of the true parameters has become biased by the stochastic noise! You can test this by increasing the number of sample points, which increases the size of the noise sample, and reduces the bias.
Finally, we can look at the bit that really matters: The model predictions made from models with the parameters we found (a posterior predictive check). Thes can be plotted using the series() method.
In [24]:
fig, axes = pints.plot.series(chain, problem)
# Customise the plot, and add the original, noise-free data
fig.set_size_inches(12,4.5)
plt.plot(times, org_values, c='orange', label='Noise-free data')
plt.legend()
plt.show()