<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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$.
Step1:  Generate independent samples from this distribution and plot them
Step2:  Use adaptive covariance MCMC to sample from this (un-normalised) pdf.
Step3:  Scatter plot of the samples. Adaptive covariance MCMC seems to do ok at sampling from this distribution.
Step4:  Try Hamiltonian Monte Carlo on same problem.
Step5:  A single chain of HMC moves much more naturally around the annulus.
Step6:  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.
Step7:  The samples are near to the surface of a sphere of radius 20.
Step8:  We can see that the mean of the samples is a long way from the true value (0, 0, 0)
Step9:  10-dimensional annulus
Now creating a 10 dimensional annulus with $r0=15$ and $\sigma=2$, then using Adaptive covariance MCMC to sample from it.
Step10:  Compare the theoretical mean and variance of the normed distance from the origin with the sample-based estimates. Does ok!
Step11:  Less good at recapitulating the actual mean.
<ASSISTANT_TASK:>
Python Code:
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()
Explanation: 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$.
End of explanation
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()
Explanation: Generate independent samples from this distribution and plot them
End of explanation
# 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]
Explanation: Use adaptive covariance MCMC to sample from this (un-normalised) pdf.
End of explanation
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()
Explanation: Scatter plot of the samples. Adaptive covariance MCMC seems to do ok at sampling from this distribution.
End of explanation
# 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]
Explanation: Try Hamiltonian Monte Carlo on same problem.
End of explanation
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()
Explanation: A single chain of HMC moves much more naturally around the annulus.
End of explanation
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)
Explanation: 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.
End of explanation
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()
Explanation: The samples are near to the surface of a sphere of radius 20.
End of explanation
a_mean = np.mean(stacked, axis=0)
print("True mean = " + str(log_pdf.mean()))
print("Sample mean = " + str(a_mean))
Explanation: We can see that the mean of the samples is a long way from the true value (0, 0, 0)
End of explanation
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]
Explanation: 10-dimensional annulus
Now creating a 10 dimensional annulus with $r0=15$ and $\sigma=2$, then using Adaptive covariance MCMC to sample from it.
End of explanation
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))
Explanation: Compare the theoretical mean and variance of the normed distance from the origin with the sample-based estimates. Does ok!
End of explanation
a_mean = np.mean(chain, axis=0)
print("True mean = " + str(log_pdf.mean()))
print("Sample mean = " + str(a_mean))
Explanation: Less good at recapitulating the actual mean.
End of explanation
<END_TASK>