<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Noise model diagnostics
Step1:  Visualisation of the data
After obtaining these parameters, it is useful to visualise the data and the fit.
Step2:  Plotting autocorrelation of the residuals
Next we look at the autocorrelation plot of the residuals to evaluate the noise model (using pints.residuals_diagnostics.plot_residuals_autocorrelation).
Step3:  The figure shows no significant autocorrelation in the residuals. Therefore, the assumption of independent noise may be valid.
Case 2
Step4:  Visualisation of the data
As before we plot the data and the inferred trajectory.
Step5:  Now the autocorrelation plot of the residuals shows high autocorrelation at small lags, which is typical of AR(1) noise. Therefore, this visualisation suggests that the assumption of independent Gaussian noise which we made during inference is invalid.
Case 3
Step6:  Visualisation of the data
As before we plot the data and the inferred trajectories.
<ASSISTANT_TASK:>
Python Code:
import pints
import pints.toy as toy
import pints.plot
import numpy as np
import matplotlib.pyplot as plt
# Use the toy logistic model
model = toy.LogisticModel()
real_parameters = [0.015, 500]
times = np.linspace(0, 1000, 100)
org_values = model.simulate(real_parameters, times)
# Add independent Gaussian noise
noise = 50
values = org_values + np.random.normal(0, noise, org_values.shape)
# Set up the problem and run the optimisation 
problem = pints.SingleOutputProblem(model, times, values)
score = pints.SumOfSquaresError(problem)
boundaries = pints.RectangularBoundaries([0, 200], [1, 1000])
x0 = np.array([0.5, 500])
found_parameters, found_value = pints.optimise(
    score,
    x0,
    boundaries=boundaries,
    method=pints.XNES,
    )
print('Score at true solution: ')
print(score(real_parameters))
print('Found solution:          True parameters:' )
for k, x in enumerate(found_parameters):
    print(pints.strfloat(x) + '    ' + pints.strfloat(real_parameters[k]))
Explanation: Noise model diagnostics: autocorrelation of the residuals
This example shows how to use the autocorrelation plots of the residuals to check assumptions of the noise model.
Three cases are shown. In the first two, optimisation is used to obtain a best-fit parameter vector in a single output problem. In the first case the noise is correctly specified and in the second case the noise is misspecified. The third case demonstrates the same method in a multiple output problem with Bayesian inference.
Case 1: Correctly specified noise
For the first example, we will use optimisation to obtain the best-fit parameter vector. See Optimisation First Example for more details. We begin with a problem in which the noise is correctly specified: both the data generation and the model use independent Gaussian noise.
End of explanation
fig, ax = pints.plot.series(np.array([found_parameters]), problem, ref_parameters=real_parameters)
fig.set_size_inches(15, 7.5)
plt.show()
Explanation: Visualisation of the data
After obtaining these parameters, it is useful to visualise the data and the fit.
End of explanation
from pints.residuals_diagnostics import plot_residuals_autocorrelation
# Plot the autocorrelation
fig = plot_residuals_autocorrelation(np.array([found_parameters]),
                                     problem)
plt.show()
Explanation: Plotting autocorrelation of the residuals
Next we look at the autocorrelation plot of the residuals to evaluate the noise model (using pints.residuals_diagnostics.plot_residuals_autocorrelation).
End of explanation
import pints.noise
# Use the toy logistic model
model = toy.LogisticModel()
real_parameters = [0.015, 500]
times = np.linspace(0, 1000, 100)
org_values = model.simulate(real_parameters, times)
# Add AR(1) noise
rho = 0.75
sigma = 50
values = org_values + pints.noise.ar1(rho, sigma, len(org_values))
# Set up the problem and run the optimisation 
problem = pints.SingleOutputProblem(model, times, values)
score = pints.SumOfSquaresError(problem)
boundaries = pints.RectangularBoundaries([0, 200], [1, 1000])
x0 = np.array([0.5, 500])
found_parameters, found_value = pints.optimise(
    score,
    x0,
    boundaries=boundaries,
    method=pints.XNES,
    )
print('Score at true solution: ')
print(score(real_parameters))
print('Found solution:          True parameters:' )
for k, x in enumerate(found_parameters):
    print(pints.strfloat(x) + '    ' + pints.strfloat(real_parameters[k]))
Explanation: The figure shows no significant autocorrelation in the residuals. Therefore, the assumption of independent noise may be valid.
Case 2: Incorrectly specified noise
For the next case, we generate data with an AR(1) (first order autoregressive) noise model. However, we deliberately misspecify the model and assume independent Gaussian noise (as before) when fitting the parameters.
End of explanation
fig, ax = pints.plot.series(np.array([found_parameters]), problem, ref_parameters=real_parameters)
fig.set_size_inches(15, 7.5)
plt.show()
# Plot the autocorrelation
fig = plot_residuals_autocorrelation(np.array([found_parameters]),
                                     problem)
plt.show()
Explanation: Visualisation of the data
As before we plot the data and the inferred trajectory.
End of explanation
import numpy as np
import matplotlib.pyplot as plt
import pints
import pints.toy
model = pints.toy.LotkaVolterraModel()
times = np.linspace(0, 3, 50)
parameters = model.suggested_parameters()
model.set_initial_conditions([2, 2])
org_values = model.simulate(parameters, times)
# Add noise
sigma = 0.1
values = org_values + np.random.normal(0, sigma, org_values.shape)
# 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([0, 0, 0, 0, 0, 0], [6, 6, 6, 6, 1, 1])
log_likelihood = pints.GaussianLogLikelihood(problem)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Run MCMC on the noisy data
x0 = [[1.0, 1.0, 1.0, 1.0, 0.1, 0.1]]*3
mcmc = pints.MCMCController(log_posterior, 3, x0)
mcmc.set_max_iterations(4000)
print('Running')
chains = mcmc.run()
print('Done!')
Explanation: Now the autocorrelation plot of the residuals shows high autocorrelation at small lags, which is typical of AR(1) noise. Therefore, this visualisation suggests that the assumption of independent Gaussian noise which we made during inference is invalid.
Case 3: Multiple output Bayesian inference problem
The plot_residuals_autocorrelation function also works with Bayesian inference and multiple output problems. For the final example, we demonstrate the same strategy in this setting.
For this example, the Lotka-Volterra model is used. See the Lotka-Volterra example for more details. As in Case 1, the true data is generated with independent Gaussian noise.
End of explanation
# Get the first MCMC chain
chain1 = chains[0]
# Cut off the burn-in samples
chain1 = chain1[2500:]
fig, ax = pints.plot.series(chain1, problem, ref_parameters=parameters)
fig.set_size_inches(15, 7.5)
plt.show()
# Plot the autocorrelation
fig = plot_residuals_autocorrelation(chain1, problem)
plt.show()
Explanation: Visualisation of the data
As before we plot the data and the inferred trajectories.
End of explanation
<END_TASK>