This example shows how the Goodwin's Oscillator toy model can be used.
Our version of this model has five parameters and three oscillating states as described in [1].
[1] Estimating Bayes factors via thermodynamic integration and population MCMC. Ben Calderhead and Mark Girolami, 2009, Computational Statistics and Data Analysis.
In [1]:
import pints
import pints.toy
import pints.plot
import matplotlib.pyplot as plt
import numpy as np
import time
model = pints.toy.GoodwinOscillatorModel()
We can get an example set of parameters using the suggested_parameters() method:
In [2]:
real_parameters = model.suggested_parameters()
print(real_parameters)
In the same way, we can get a suggested set of sampling times:
In [3]:
times = model.suggested_times()
Now we can run a simulation:
In [4]:
values = model.simulate(real_parameters, times)
This gives us all we need to create a plot of current versus time:
In [5]:
plt.figure()
plt.subplot(3, 1, 1)
plt.plot(times, values[:, 0], 'b')
plt.subplot(3, 1, 2)
plt.plot(times, values[:, 1], 'g')
plt.subplot(3, 1, 3)
plt.plot(times, values[:, 2], 'r')
plt.show()
Now we will add some noise to generate some fake "experimental" data and try to recover the original parameters.
In [6]:
noise1 = 0.001
noise2 = 0.01
noise3 = 0.1
noisy_values = np.array(values, copy=True)
noisy_values[:, 0] += np.random.normal(0, noise1, len(times))
noisy_values[:, 1] += np.random.normal(0, noise2, len(times))
noisy_values[:, 2] += np.random.normal(0, noise3, len(times))
plt.figure()
plt.subplot(3, 1, 1)
plt.plot(times, noisy_values[:, 0], 'b')
plt.subplot(3, 1, 2)
plt.plot(times, noisy_values[:, 1], 'g')
plt.subplot(3, 1, 3)
plt.plot(times, noisy_values[:, 2], 'r')
plt.show()
Now we can try and infer the original parameters:
In [7]:
# Create an object with links to the model and time series
problem = pints.MultiOutputProblem(model, times, values)
# Create a log posterior
log_prior = pints.UniformLogPrior([1, 1, 0.01, 0.01, 0.01], [10, 10, 1, 1, 1])
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, [noise1, noise2, noise3])
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Run MCMC on the noisy data
x0 = [[5, 5, 0.5, 0.5, 0.5]]*3
mcmc = pints.MCMCController(log_posterior, 3, x0)
mcmc.set_max_iterations(5000)
mcmc.set_log_to_screen(False)
start = time.time()
print('Running')
chains = mcmc.run()
print('Done!')
end = time.time()
diff = end - start
Print results.
In [8]:
results = pints.MCMCSummary(chains=chains, time=diff,
parameter_names=["k2", "k3", "m1", "m2", "m3"])
print(results)
Now we can inspect the resulting chains:
In [9]:
pints.plot.trace(chains, ref_parameters=real_parameters)
plt.show()
This is a pretty hard problem!
And what about optimisation?
In [10]:
# Fit to the noisy data
parameters = []
opt = pints.OptimisationController(log_posterior, x0[0], method=pints.XNES)
opt.set_log_to_screen(False)
parameters, fbest = opt.run()
print('')
print(' p1 p2 p3 p4 p5')
print('real ' + ' '.join(['{: 8.4g}'.format(float(x)) for x in real_parameters]))
print('found ' + ' '.join(['{: 8.4g}'.format(x) for x in parameters]))
The Goodwin-oscillator model has sensitivities calculated by the forward sensitivities approach, so we can use samplers that use gradients (although they will be slower per iteration; although perhaps not by ESS per second!), like Relativistic HMC.
In [11]:
problem = pints.MultiOutputProblem(model, times, noisy_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, 0, 0, 0, 0, 0, 0, 0],
[10, 10, 1, 1, 1, 1, 1, 1]
)
# Create a posterior log-likelihood (log(likelihood * prior))
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Choose starting points for 3 mcmc chains
real_parameters1 = np.array(real_parameters.tolist() + [noise1, noise2, noise3])
xs = [
real_parameters1 * 1.1,
real_parameters1 * 0.9,
real_parameters1 * 1.15,
real_parameters1 * 1.2,
]
# Create mcmc routine
mcmc = pints.MCMCController(log_posterior, 4, xs, method=pints.RelativisticMCMC)
# Add stopping criterion
mcmc.set_max_iterations(200)
# Run in parallel
mcmc.set_parallel(True)
mcmc.set_log_interval(1)
for sampler in mcmc.samplers():
sampler.set_leapfrog_step_size([0.1, 0.5, 0.002, 0.002, 0.002, 0.0005, 0.001, 0.01])
sampler.set_leapfrog_steps(10)
# time start
start = time.time()
# Run!
print('Running...')
chains = mcmc.run()
print('Done!')
# end time
end = time.time()
diff = end - start
Print results.
In [12]:
results = pints.MCMCSummary(chains=chains, time=diff,
parameter_names=["k2", "k3", "m1", "m2", "m3",
"sigma_x", "sigma_y", "sigma_z"])
print(results)
In [13]:
pints.plot.trace(chains, ref_parameters=real_parameters1)
plt.show()
Plot posterior predictive distribution.
In [14]:
pints.plot.series(np.vstack(chains), problem)
plt.show()