This example shows how PINTS handles multiplicative Gaussian noise. In a multiplicative Gaussian noise process, the variance of the noise term scales with the magnitude of the output:
$$X(t) = f(t;\theta) + f(t;\theta)^\eta \epsilon(t)$$where
$$\epsilon(t) \overset{\text{iid}}{\sim} N(0, \sigma)$$Therefore, portions of the time-series with greater magnitude will also see a higher standard deviation of the error terms.
In this example, we generate data according to this noise model, and use optimisation to infer the parameters.
In [7]:
from __future__ import print_function
import pints
import pints.toy as toy
import pints.plot
import numpy as np
import matplotlib.pyplot as plt
import pints.noise
# Use the toy logistic model
model = toy.LogisticModel()
# Set some parameters
real_parameters = [0.1, 50]
# Create fake data
times = np.linspace(0, 100, 500)
values = model.simulate(real_parameters, times)
# Add noise according to multiplicative Gaussian
noisy_values = values + pints.noise.multiplicative_gaussian(0.8, 0.5, values)
# Create an inference problem
problem = pints.SingleOutputProblem(model, times, noisy_values)
# Show the generated data
plt.figure(figsize=(15,10))
plt.xlabel('Time')
plt.ylabel('Values')
plt.plot(times, noisy_values, label='data')
plt.plot(times, values, label='true values')
plt.legend()
plt.show()
In this case, the way the noise magnitude increases with the magnitude of the function is visible in the plot of the data above. However, if it were not so clear, we might initially make the incorrect assumption of independent Gaussian noise. In the code below, we estimate the parameters under this incorrect assumption.
In [8]:
# Make a log likelihood assuming independent Gaussian noise
log_likelihood = pints.GaussianLogLikelihood(problem)
# Set up the optimisation problem
y0 = np.array([.1, 10, 1])
boundaries = pints.RectangularBoundaries([0, 5, 1e-3], [1, 500, 10])
opt = pints.OptimisationController(log_likelihood,
y0,
boundaries=boundaries,
method=pints.XNES)
# Run the optimisation
y1, g1 = opt.run()
print('Estimated parameters:')
print(y1)
After performing the fit based on the incorrect assumption of iid Gaussian noise, it is useful to look at a diagnostic plot of residuals vs. function magnitude. This is possible using the plot_residuals_vs_output function in pints.residuals_diagnostics.
In [9]:
import pints.residuals_diagnostics
fig = pints.residuals_diagnostics.plot_residuals_vs_output(
np.array([[y1[0], y1[1]]]),
problem)
plt.show()
This plot shows higher magnitudes in the residuals as the function output increases, suggesting that multiplicative Gaussian noise may be applicable.
In [10]:
# Show the generated data
simulated_values = problem.evaluate(y1[:2])
# Get the standard deviation (which is constant in this noise model)
sigma = y1[2]
plt.figure()
plt.xlabel('Time')
plt.ylabel('Values')
plt.plot(times, noisy_values, label='data', alpha=0.5)
plt.fill_between(times,
simulated_values - sigma,
simulated_values + sigma,
alpha=0.2)
plt.plot(times, simulated_values, label='best fit parameters')
plt.legend()
plt.show()
Similarly, this plot of the data and the fit also highlights the inaccuracy of the independent Gaussian assumption. Although the best fit model parameters are plausible, the measurement noise variance is being overestimated at small function values and underestimated at large function values.
The next step is to use the correct likelihood , which is available as pints.MultiplicativeGaussianLogLikelihood.
In [11]:
# Make a log likelihood assuming the correct multiplicative noise
log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
# Set up the optimisation
y0 = np.array([.1, 10, 0.25, .1])
boundaries = pints.RectangularBoundaries([0, 5, 1e-3, 1e-3], [1, 500, 10, 10])
opt = pints.OptimisationController(log_likelihood,
y0,
boundaries=boundaries,
method=pints.XNES)
# Run the optimisation
y1, g1 = opt.run()
print('Estimated parameters:')
print(y1)
In [12]:
# Show the generated data
simulated_values = problem.evaluate(y1[:2])
# Get the inferred noise model parameters
eta = y1[2]
sigma = y1[3]
plt.figure()
plt.xlabel('Time')
plt.ylabel('Values')
plt.plot(times, noisy_values, alpha=0.5, label='data')
plt.fill_between(times,
simulated_values - sigma*simulated_values**eta,
simulated_values + sigma*simulated_values**eta,
alpha=0.2)
plt.plot(times, simulated_values, label='best fit parameters')
plt.legend()
plt.show()
With the correct noise model, we can see that the standard deviation of the measurement process scales appropriately with the magnitude of the trajectory.