This example shows you how to perform Bayesian inference on a Normal distribution and a time-series problem, using MALA for MCMC.
First, we create a simple normal distribution
import pints
import pints.toy
import numpy as np
import matplotlib.pyplot as plt
# Create log pdf
log_pdf = pints.toy.GaussianLogPDF([2, 4], [[1, 0], [0, 3]])
# Contour plot of pdf
levels = np.linspace(-3,12,20)
num_points = 100
x = np.linspace(-1, 5, num_points)
y = np.linspace(-0, 8, num_points)
X, Y = np.meshgrid(x, y)
Z = np.zeros(X.shape)
Z = np.exp([[log_pdf([i, j]) for i in x] for j in y])
plt.contour(X, Y, Z)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Now we set up and run a sampling routine using MALA MCMC
# Choose starting points for 3 mcmc chains
xs = [
[2, 1],
[3, 3],
[5, 4],
]
# Create mcmc routine
mcmc = pints.MCMCController(log_pdf, 3, xs, method=pints.MALAMCMC)
# Add stopping criterion
mcmc.set_max_iterations(2000)
# Set up modest logging
mcmc.set_log_to_screen(True)
mcmc.set_log_interval(100)
# # Update step sizes used by individual samplers (which is then scaled by sigma0)
for sampler in mcmc.samplers():
sampler.set_epsilon([1.5, 1.5])
# Run!
print('Running...')
full_chains = mcmc.run()
print('Done!')
Running... Using Metropolis-Adjusted Langevin Algorithm (MALA) Generating 3 chains. Running in sequential mode. Iter. Eval. Accept. Accept. Accept. Time m:s 0 3 0 0 0 0:00.0 1 6 0.5 0.5 0.5 0:00.0 2 9 0.667 0.667 0.667 0:00.0 3 12 0.5 0.5 0.5 0:00.0 100 303 0.683 0.703 0.683 0:00.2 200 603 0.681592 0.711 0.692 0:00.3 300 903 0.698 0.714 0.724 0:00.4 400 1203 0.698 0.731 0.738 0:00.6 500 1503 0.681 0.737 0.738523 0:00.7 600 1803 0.702 0.74 0.73 0:00.8 700 2103 0.709 0.743224 0.729 0:00.9 800 2403 0.710362 0.745 0.734 0:01.1 900 2703 0.716 0.74 0.736 0:01.2 1000 3003 0.711 0.734 0.737 0:01.3 1100 3303 0.715713 0.736 0.738 0:01.4 1200 3603 0.715 0.734388 0.741 0:01.6 1300 3903 0.718 0.739 0.741 0:01.7 1400 4203 0.719 0.739 0.738 0:01.8 1500 4503 0.725 0.736 0.732 0:02.0 1600 4803 0.723 0.737664 0.732 0:02.1 1700 5103 0.725 0.739565 0.731 0:02.2 1800 5403 0.729 0.741 0.733 0:02.3 1900 5703 0.729616 0.741 0.732 0:02.5 2000 6000 0.7285 0.739 0.7325 0:02.6 Halting: Maximum number of iterations (2000) reached. Done!
# Show traces and histograms
import pints.plot
pints.plot.trace(full_chains)
plt.show()
# Discard warm up
chains = full_chains[:, 1000:]
# Check convergence and other properties of chains
results = pints.MCMCSummary(chains=chains, time=mcmc.time(), parameter_names=['mean_x', 'mean_y'])
print(results)
# Look at distribution in chain 0
pints.plot.pairwise(chains[0], kde=True)
plt.show()
# Check Kullback-Leibler divergence of chains
print(log_pdf.kl_divergence(chains[0]))
print(log_pdf.kl_divergence(chains[1]))
print(log_pdf.kl_divergence(chains[2]))
param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ess per sec. ------- ------ ------ ------ ----- ----- ----- ------- ------ ------ -------------- mean_x 2.00 1.01 0.04 1.30 1.98 2.67 3.99 1.00 924.67 355.47 mean_y 4.08 1.65 1.00 2.95 4.06 5.22 7.25 1.01 441.72 169.81
0.012208212988188416 0.003970589498094368 0.006245525296897592
We now try the same method on a time-series problem
First, we try it in 1d, using a wrapper around the LogisticModel to make it one-dimensional.
import pints.toy as toy
# Create a wrapper around the logistic model, turning it into a 1d model
class Model(pints.ForwardModel):
def __init__(self):
self.model = toy.LogisticModel()
def simulate(self, x, times):
return self.model.simulate([x[0], 500], times)
def simulateS1(self, x, times):
values, gradient = self.model.simulateS1([x[0], 500], times)
gradient = gradient[:, 0]
return values, gradient
def n_parameters(self):
return 1
# Load a forward model
model = Model()
# Create some toy data
real_parameters = np.array([0.015])
times = np.linspace(0, 1000, 50)
org_values = model.simulate(real_parameters, times)
# Add noise
np.random.seed(1)
noise = 10
values = org_values + np.random.normal(0, noise, org_values.shape)
plt.figure()
plt.plot(times, values)
plt.plot(times, org_values)
plt.show()
# Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, values)
# Create a log-likelihood function
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, noise)
# Create a uniform prior over the parameters
log_prior = pints.UniformLogPrior(
[0.01],
[0.02]
)
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Choose starting points for mcmc chains
xs = [
real_parameters * 1.01,
real_parameters * 0.9,
real_parameters * 1.15,
]
# Create mcmc routine
mcmc = pints.MCMCController(log_likelihood, len(xs), xs, method=pints.MALAMCMC)
# Add stopping criterion
mcmc.set_max_iterations(2000)
# Set up modest logging
mcmc.set_log_to_screen(True)
mcmc.set_log_interval(100)
# Run!
print('Running...')
chains = mcmc.run()
print('Done!')
Running... Using Metropolis-Adjusted Langevin Algorithm (MALA) Generating 3 chains. Running in sequential mode. Iter. Eval. Accept. Accept. Accept. Time m:s 0 3 0 0 0 0:00.0 1 6 0.5 0.5 0.5 0:00.0 2 9 0.667 0.667 0.667 0:00.0 3 12 0.75 0.75 0.75 0:00.0 100 303 0.990099 0.990099 0.990099 0:00.2 200 603 0.995 0.995 0.99 0:00.3 300 903 0.993 0.993 0.99 0:00.4 400 1203 0.995 0.993 0.993 0:00.6 500 1503 0.992016 0.988024 0.99002 0:00.7 600 1803 0.992 0.99 0.99 0:00.9 700 2103 0.993 0.991 0.991 0:01.0 800 2403 0.99 0.993 0.993 0:01.1 900 2703 0.988 0.993 0.992 0:01.3 1000 3003 0.988012 0.994006 0.993007 0:01.4 1100 3303 0.988 0.995 0.993 0:01.6 1200 3603 0.989 0.994 0.993 0:01.7 1300 3903 0.99 0.995 0.992 0:01.8 1400 4203 0.99 0.995 0.992 0:02.0 1500 4503 0.99 0.995 0.991 0:02.1 1600 4803 0.991 0.995 0.992 0:02.2 1700 5103 0.991 0.995 0.991 0:02.4 1800 5403 0.992 0.995 0.991116 0:02.5 1900 5703 0.992 0.995 0.992 0:02.7 2000 6000 0.992 0.9955 0.992 0:02.8 Halting: Maximum number of iterations (2000) reached. Done!
# Show trace and histogram
pints.plot.trace(chains)
plt.show()
# Show predicted time series for the first chain
pints.plot.series(chains[0, 200:], problem, real_parameters)
plt.show()
Finally, we try MALA MCMC on a 2d logistic model problem.
import pints
import pints.toy as toy
import pints.plot
import numpy as np
import matplotlib.pyplot as plt
# Load a forward model
model = toy.LogisticModel()
# Create some toy data
real_parameters = np.array([0.015, 500])
org_values = model.simulate(real_parameters, times)
# Add noise
np.random.seed(1)
noise = 10
values = org_values + np.random.normal(0, noise, org_values.shape)
# Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, values)
# Create a log-likelihood function
log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, noise)
# Create a uniform prior over the parameters
log_prior = pints.UniformLogPrior(
[0.01, 400],
[0.02, 600]
)
# Create a posterior log-likelihood (log(likelihood * prior))
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Choose starting points for 3 mcmc chains
xs = [
real_parameters * 1.01,
real_parameters * 0.9,
real_parameters * 1.1,
]
# Create mcmc routine
mcmc = pints.MCMCController(log_posterior, len(xs), xs, method=pints.MALAMCMC)
# Add stopping criterion
mcmc.set_max_iterations(4000)
# Set up modest logging
mcmc.set_log_to_screen(True)
mcmc.set_log_interval(100)
# Run!
print('Running...')
chains = mcmc.run()
print('Done!')
Running... Using Metropolis-Adjusted Langevin Algorithm (MALA) Generating 3 chains. Running in sequential mode. Iter. Eval. Accept. Accept. Accept. Time m:s 0 3 0 0 0 0:00.0 1 6 0.5 0.5 0.5 0:00.0 2 9 0.333 0.667 0.667 0:00.0 3 12 0.5 0.75 0.75 0:00.0 100 303 0.970297 0.970297 0.980198 0:00.1 200 603 0.975 0.98 0.99 0:00.3 300 903 0.98 0.983 0.986711 0:00.5 400 1203 0.975 0.985 0.985 0:00.6 500 1503 0.976 0.984 0.986 0:00.8 600 1803 0.975 0.987 0.985025 0:00.9 700 2103 0.977 0.989 0.984 0:01.0 800 2403 0.974 0.99 0.985 0:01.2 900 2703 0.977 0.99 0.984 0:01.3 1000 3003 0.979021 0.991009 0.983017 0:01.5 1100 3303 0.98 0.992 0.981 0:01.6 1200 3603 0.981 0.992 0.981 0:01.7 1300 3903 0.982 0.991545 0.982 0:01.9 1400 4203 0.983 0.991 0.982 0:02.0 1500 4503 0.983 0.99 0.983 0:02.2 1600 4803 0.984 0.989 0.983 0:02.3 1700 5103 0.985 0.99 0.984 0:02.4 1800 5403 0.985 0.991 0.984 0:02.6 1900 5703 0.986 0.991 0.983 0:02.7 2000 6003 0.986 0.99 0.983 0:02.9 2100 6303 0.985 0.989 0.982 0:03.0 2200 6603 0.985 0.989 0.982 0:03.1 2300 6903 0.985 0.989 0.982 0:03.3 2400 7203 0.985 0.99 0.982 0:03.4 2500 7503 0.984 0.990004 0.981 0:03.5 2600 7803 0.984 0.99 0.981 0:03.7 2700 8103 0.984 0.99 0.981 0:03.8 2800 8403 0.984 0.99 0.981 0:03.9 2900 8703 0.984 0.991 0.981 0:04.1 3000 9003 0.984 0.991 0.981 0:04.2 3100 9303 0.984 0.99 0.981 0:04.3 3200 9603 0.984 0.991 0.982 0:04.5 3300 9903 0.985 0.991 0.982 0:04.6 3400 10203 0.985 0.990885 0.983 0:04.8 3500 10503 0.985 0.991 0.982862 0:04.9 3600 10803 0.985 0.991 0.983 0:05.0 3700 11103 0.985 0.991 0.982167 0:05.2 3800 11403 0.986 0.991 0.982 0:05.3 3900 11703 0.986 0.991 0.982 0:05.4 4000 12000 0.986 0.99125 0.98225 0:05.6 Halting: Maximum number of iterations (4000) reached. Done!
# Check convergence and other properties of chains
results = pints.MCMCSummary(chains=chains[:,1000:], time=mcmc.time(), parameter_names=['growth rate', 'capacity'])
print(results)
# Show traces and histograms
pints.plot.trace(chains)
plt.show()
param mean std. 2.5% 25% 50% 75% 97.5% rhat ess ess per sec. ----------- ------ ------ ------ ------ ------ ------ ------- ------ ------ -------------- growth rate 0.01 0.00 0.01 0.01 0.01 0.02 0.02 1.00 362.16 64.88 capacity 500.53 2.06 496.33 499.21 500.55 501.88 504.54 1.00 510.86 91.52
Chains have converged!