<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Goodwin's oscillator toy model
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.
We start by creating a toy model
Step1:  The model also provides suggested parameters and sampling times, allowing us to run a simulation
Step2:  This gives us all we need to create a plot of current versus time
Step3:  Now we will add some noise to generate some fake "experimental" data and try to recover the original parameters.
Step4:  Now we can try and infer the original parameters
Step5:  We can use an MCMCSummary to display the results
Step6:  Now we can inspect the resulting chains
Step7:  This is a pretty hard problem!
And what about optimisation?
Step8:  Sampling using relativistic HMC
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.
Step9:  Display the results
Step10:  Plot posterior predictive distribution.
<ASSISTANT_TASK:>
Python Code:
import pints
import pints.plot
import pints.toy
import matplotlib.pyplot as plt
import numpy as np
model = pints.toy.GoodwinOscillatorModel()
Explanation: Goodwin's oscillator toy model
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.
We start by creating a toy model:
End of explanation
real_parameters = model.suggested_parameters()
times = model.suggested_times()
values = model.simulate(real_parameters, times)
Explanation: The model also provides suggested parameters and sampling times, allowing us to run a simulation:
End of explanation
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()
Explanation: This gives us all we need to create a plot of current versus time:
End of explanation
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()
Explanation: Now we will add some noise to generate some fake "experimental" data and try to recover the original parameters.
End of explanation
# 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)
print('Running')
chains = mcmc.run()
print('Done!')
Explanation: Now we can try and infer the original parameters:
End of explanation
results = pints.MCMCSummary(
    chains=chains, 
    time=mcmc.time(),
    parameter_names=['k2', 'k3', 'm1', 'm2', 'm3']
)
print(results)
Explanation: We can use an MCMCSummary to display the results:
End of explanation
pints.plot.trace(chains, ref_parameters=real_parameters)
plt.show()
Explanation: Now we can inspect the resulting chains:
End of explanation
# 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]))
Explanation: This is a pretty hard problem!
And what about optimisation?
End of explanation
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)
# Tune the samplers' hyper-parameters
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)
# Run!
print('Running...')
chains = mcmc.run()
print('Done!')
Explanation: Sampling using relativistic HMC
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.
End of explanation
results = pints.MCMCSummary(
    chains=chains, 
    time=mcmc.time(),
    parameter_names=['k2', 'k3', 'm1', 'm2', 'm3', 'sigma_x', 'sigma_y', 'sigma_z'])
print(results)
pints.plot.trace(chains, ref_parameters=real_parameters1)
plt.show()
Explanation: Display the results:
End of explanation
pints.plot.series(np.vstack(chains), problem)
plt.show()
Explanation: Plot posterior predictive distribution.
End of explanation
<END_TASK>