This notebook contains a longer version of the example shown on the PINTS repository landing page. In it, we quickly showcase some PINTS functionality, but without going into detail. For a more thorough introduction, please see the full set of examples, here.
import matplotlib.pyplot as plt
import numpy as np
import pints
# Define a custom model, in this case a variant of the Hodgkin-Huxley
# model of ionic currents, after a single voltage-step
class MyModel(pints.ForwardModel):
def n_parameters(self):
return 3
def simulate(self, parameters, times):
tau_m, tau_h, g = parameters
V = 20
m0 = 0.001
h0 = 0.999
m = 0.999 - (0.999 - m0) * np.exp(-times / tau_m)
h = 0.001 - (0.001 - h0) * np.exp(-times / tau_h)
return g * m**3 * h * (V - 50)
# Instantiate the model
model = MyModel()
# Load experimental data
# For this example, we'll generate synthetic data with added Gaussian noise
def load_experimental_data():
times = np.linspace(0, 20, 1000)
values = model.simulate([0.5, 3, 4], times)
values += np.random.normal(size=values.shape) * 3
return times, values
times, values = load_experimental_data()
# Show the noisy data
plt.figure()
plt.plot(times, values)
plt.show()
# Combine the forward model and experimental data into an inference problem
problem = pints.SingleOutputProblem(model, times, values)
# Select one of PINTS' error measures, and run an optimisation to infer the parameters
error_measure = pints.SumOfSquaresError(problem)
initial_point = [1, 2, 3]
optimisation = pints.OptimisationController(
error_measure, initial_point, method=pints.XNES)
optimisation.set_log_to_screen(False)
parameters, error = optimisation.run()
/usr/lib/python3.6/site-packages/ipykernel_launcher.py:15: RuntimeWarning: overflow encountered in exp from ipykernel import kernelapp as app /usr/lib/python3.6/site-packages/ipykernel_launcher.py:17: RuntimeWarning: overflow encountered in power /usr/lib/python3.6/site-packages/ipykernel_launcher.py:17: RuntimeWarning: overflow encountered in multiply /home/michael/dev/pints/pints/_error_measures.py:297: RuntimeWarning: overflow encountered in square return np.sum((np.sum(((self._problem.evaluate(x) - self._values)**2),
# Show the quality of fit
plt.figure()
plt.plot(times, values, alpha=0.5)
plt.plot(times, model.simulate(parameters, times), lw=2)
plt.show()
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, sigma=3)
optimisation = pints.OptimisationController(
error_measure, initial_point, method=pints.XNES)
optimisation.set_log_to_screen(False)
parameters, error = optimisation.run()
plt.figure()
plt.plot(times, values, alpha=0.5)
plt.plot(times, model.simulate(parameters, times), lw=2)
plt.show()
# Set up a sampling routine with 3 MCMC chains, all starting from the same point
mcmc = pints.MCMCController(
log_likelihood, 3, [initial_point]*3, method=pints.AdaptiveCovarianceMCMC)
mcmc.set_max_iterations(2000)
mcmc.set_log_to_screen(False)
samples = mcmc.run()
# Use a diagnostic plot to check if the chains have converged
import pints.plot
pints.plot.trace(samples)
plt.show()
# Plot several predictions that are all likely sources of the experimental data
pints.plot.series(np.vstack(samples[:,1000:]), problem)
plt.show()
pints.plot.series(np.vstack(samples[:,1000:]), problem)
plt.xlim(1, 2)
plt.ylim(-64, -55)
plt.show()