Annulus distribution

In this notebook we create an annulus distribution, which has the following density,

$$p(x|r_0, \sigma) \propto \text{exp}(-\frac{(|x|-r_0)^2}{2\sigma^2}),$$

where $\beta > 0$ and $|x|$ is the Euclidean norm of the d-dimensional $x$. This distribution was created due to its similarity in geometry to the distribution of a high-dimensional normal (albeit in lower dimensions), for its use in testing MCMC algorithms.

Plotting a two-dimensional version of this function with $r0=10$ and $\sigma=1$.


In [2]:
import pints
import pints.toy
import numpy as np
import matplotlib.pyplot as plt

# Create log pdf (default is 2-dimensional with r0=10 and sigma=1)
log_pdf = pints.toy.AnnulusLogPDF()

# Contour plot of pdf
num_points = 100
x = np.linspace(-15, 15, num_points)
y = np.linspace(-15, 15, 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()


Generate independent samples from this distribution and plot them


In [3]:
samples = log_pdf.sample(100)

num_points = 100
x = np.linspace(-15, 15, num_points)
y = np.linspace(-15, 15, 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.scatter(samples[:, 0], samples[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.show()


Use adaptive covariance MCMC to sample from this (un-normalised) pdf.


In [4]:
# Create an adaptive covariance MCMC routine
x0 = np.random.uniform([2, 2], [8, 8], size=(4, 2))
mcmc = pints.MCMCController(log_pdf, 4, x0, method=pints.HaarioBardenetACMC)

# Set maximum number of iterations
mcmc.set_max_iterations(4000)

# Disable logging
mcmc.set_log_to_screen(False)

# Number of chains
num_chains = 4

# Run!
print('Running...')
chains = mcmc.run()
print('Done!')

# Discard warm-up
chains = [chain[1000:] for chain in chains]


Running...
Done!

Scatter plot of the samples. Adaptive covariance MCMC seems to do ok at sampling from this distribution.


In [5]:
stacked = np.vstack(chains)
plt.contour(X, Y, Z, colors='k', alpha=0.5)
plt.scatter(stacked[:,0], stacked[:,1], marker='.', alpha=0.2)
plt.xlim(-15, 15)
plt.ylim(-15, 15)
plt.xlabel('x')
plt.ylabel('y')
plt.show()


Try Hamiltonian Monte Carlo on same problem.


In [8]:
# Create an adaptive covariance MCMC routine
x0 = np.random.uniform([2, 2], [8, 8], size=(4, 2))
sigma0 = [2, 2]
mcmc = pints.MCMCController(log_pdf, 4, x0, method=pints.HamiltonianMCMC, sigma0=sigma0)

# Set maximum number of iterations
mcmc.set_max_iterations(500)

# Disable logging
# mcmc.set_log_to_screen(False)

# Number of chains
num_chains = 4

# Run!
print('Running...')
chains = mcmc.run()
print('Done!')

# Discard warm-up
chains = [chain[200:] for chain in chains]


Running...
Using Hamiltonian Monte Carlo
Generating 4 chains.
Running in sequential mode.
Iter. Eval. Accept.   Accept.   Accept.   Accept.   Time m:s
0     4      0         0         0         0          0:00.0
1     84     0.333     0.333     0.333     0.333      0:00.0
2     164    0.5       0.5       0.5       0.5        0:00.1
3     244    0.6       0.6       0.6       0.6        0:00.1
20    1604   0.909     0.909     0.909     0.909      0:00.4
40    3204   0.952381  0.952381  0.952381  0.952381   0:00.7
60    4804   0.968     0.968     0.968     0.968      0:01.1
80    6404   0.976     0.976     0.976     0.976      0:01.5
100   8004   0.98      0.98      0.98      0.98       0:01.8
120   9604   0.984     0.984     0.984     0.984      0:02.2
140   11204  0.986     0.986     0.986     0.986      0:02.5
160   12804  0.988     0.988     0.988     0.988      0:02.8
180   14404  0.989011  0.989011  0.989011  0.989011   0:03.2
200   16004  0.990099  0.990099  0.990099  0.990099   0:03.5
220   17604  0.990991  0.990991  0.990991  0.990991   0:03.9
240   19204  0.992     0.992     0.992     0.988      0:04.2
260   20804  0.992     0.992     0.992     0.989      0:04.7
280   22404  0.993     0.993     0.993     0.989      0:05.0
300   24004  0.993     0.993     0.993     0.99       0:05.4
320   25604  0.994     0.994     0.994     0.991      0:05.7
340   27204  0.994152  0.994152  0.994152  0.991      0:06.1
360   28804  0.994     0.994     0.994     0.992      0:06.4
380   30404  0.995     0.995     0.995     0.992      0:06.7
400   32004  0.995     0.995     0.995     0.993      0:07.1
420   33604  0.995     0.995     0.995     0.992891   0:07.4
440   35204  0.995     0.995     0.995     0.993      0:07.7
460   36804  0.995671  0.995671  0.995671  0.994      0:08.1
480   38404  0.996     0.996     0.996     0.994      0:08.4
500   39924  0.996008  0.996008  0.994012  0.994012   0:08.7
Halting: Maximum number of iterations (500) reached.
Done!

A single chain of HMC moves much more naturally around the annulus.


In [9]:
plt.contour(X, Y, Z, colors='k', alpha=0.5)
plt.plot(chains[0][:, 0], chains[0][:, 1])
plt.xlim(-15, 15)
plt.ylim(-15, 15)
plt.xlabel('x')
plt.ylabel('y')
plt.show()


3-dimensional annulus

Now creating a 3 dimensional annulus with $r0=20$ and $\sigma=0.5$, then using Adaptive covariance MCMC to sample from it.


In [16]:
log_pdf = pints.toy.AnnulusLogPDF(dimensions=3, r0=20, sigma=0.5)

# Create an adaptive covariance MCMC routine
x0 = np.zeros(log_pdf.n_parameters()) + np.random.normal(0, 1, size=(4, log_pdf.n_parameters()))
mcmc = pints.MCMCController(log_pdf, 4, x0, method=pints.HaarioBardenetACMC)

# Set maximum number of iterations
mcmc.set_max_iterations(4000)

# Disable logging
mcmc.set_log_to_screen(False)

# Number of chains
num_chains = 4

# Run!
print('Running...')
chains = mcmc.run()
print('Done!')

# Discard warm-up
chains = [chain[1000:] for chain in chains]
stacked = np.vstack(chains)


Running...
Done!

The samples are near to the surface of a sphere of radius 20.


In [17]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(stacked[:, 0], stacked[:, 1], stacked[:, 2], '.', alpha=0.1)
ax.legend()

plt.show()


No handles with labels found to put in legend.

We can see that the mean of the samples is a long way from the true value (0, 0, 0)


In [18]:
a_mean = np.mean(stacked, axis=0)

print("True mean = " + str(log_pdf.mean()))
print("Sample mean = " + str(a_mean))


True mean = [0. 0. 0.]
Sample mean = [-1.71392093 -3.99819504  1.70288497]

10-dimensional annulus

Now creating a 10 dimensional annulus with $r0=15$ and $\sigma=2$, then using Adaptive covariance MCMC to sample from it.


In [19]:
log_pdf = pints.toy.AnnulusLogPDF(dimensions=10, r0=15, sigma=2)

# Create an adaptive covariance MCMC routine
x0 = np.zeros(log_pdf.n_parameters()) + np.random.normal(0, 1, size=(4, log_pdf.n_parameters()))
mcmc = pints.MCMCController(log_pdf, 4, x0, method=pints.HaarioBardenetACMC)

# Set maximum number of iterations
mcmc.set_max_iterations(8000)

# Disable logging
mcmc.set_log_to_screen(False)

# Number of chains
num_chains = 4

# Run!
print('Running...')
chains = mcmc.run()
print('Done!')

# Discard warm-up
chains = [chain[1000:] for chain in chains]


Running...
Done!

Compare the theoretical mean and variance of the normed distance from the origin with the sample-based estimates. Does ok!


In [20]:
chain = np.vstack(chains)
d = list(map(lambda x: np.linalg.norm(x), chain))
a_mean = np.mean(d)
a_var = np.var(d)

print("True normed mean = " + str(log_pdf.mean_normed()))
print("Sample normed mean = " + str(a_mean))

print("True normed var = " + str(log_pdf.var_normed()))
print("Sample normed var = " + str(a_var))


True normed mean = 17.12813200595798
Sample normed mean = 16.12186750819023
True normed var = 3.5490740758476136
Sample normed var = 3.7994461819482375

Less good at recapitulating the actual mean.


In [21]:
a_mean = np.mean(chain, axis=0)

print("True mean = " + str(log_pdf.mean()))
print("Sample mean = " + str(a_mean))


True mean = [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Sample mean = [ 0.17639465 -0.10537406 -0.59015238  0.38752625  0.14268302  0.4043969
 -0.09440046  0.38721041  0.64282244  0.22301249]