4. Exercise sheet#
Some general remarks about the exercises:
For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.
For each part of the exercise a solution box has been added, but you may insert additional boxes. Do not hesitate to add Markdown boxes for textual or LaTeX answers (via
Cell > Cell Type > Markdown). But make sure to replace any part that saysYOUR CODE HEREorYOUR ANSWER HEREand remove theraise NotImplementedError().Please make your code readable by humans (and not just by the Python interpreter): choose informative function and variable names and use consistent formatting. Feel free to check the PEP 8 Style Guide for Python for the widely adopted coding conventions or this guide for explanation.
Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting
Kernel > Restart & Run Allin the jupyter menu.Each sheet has 100 points worth of exercises. Note that only the grades of sheets number 2, 4, 6, 8 count towards the course examination. Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are just for feedback.
Please fill in your name here and, in case you worked together closely with someone, the name of your collaborator(s) as well:
NAME = ""
COLLABORATORS = ""
Exercise sheet 4
Code from the lectures:
import numpy as np
import matplotlib.pylab as plt
import networkx as nx
rng = np.random.default_rng()
%matplotlib inline
def draw_transition_graph(P):
# construct a directed graph directly from the matrix
graph = nx.DiGraph(P)
# draw it in such a way that edges in both directions are visible and have appropriate width
nx.draw_networkx(graph,connectionstyle='arc3, rad = 0.15',width=[6*P[u,v] for u,v in graph.edges()])
def sample_next(P,current):
return rng.choice(len(P),p=P[current])
def sample_chain(P,start,n):
chain = [start]
for _ in range(n):
chain.append(sample_next(P,chain[-1]))
return chain
def stationary_distributions(P):
eigenvalues, eigenvectors = np.linalg.eig(np.transpose(P))
# make list of normalized eigenvectors for which the eigenvalue is very close to 1
return [eigenvectors[:,i]/np.sum(eigenvectors[:,i]) for i in range(len(eigenvalues))
if np.abs(eigenvalues[i]-1) < 1e-10]
def markov_sample_mean(P,start,function,n):
total = 0
state = start
for _ in range(n):
state = sample_next(P,state)
total += function[state]
return total/n
4.1. Markov Chain on a graph#
(50 points)
The goal of this exercise is to use Metropolis-Hastings to sample a uniform vertex in a (finite, undirected) connected graph \(G\). More precisely, the state space \(\Gamma = \{0,\ldots,n-1\}\) is the set of vertices of a graph and the desired probability mass function is \(\pi(x) = 1/n\) for \(x\in\Gamma\). The set of edges is denoted \(E = \{ \{x_1,y_1\}, \ldots,\{x_k,y_k\}\}\), \(x_i,y_i\in\Gamma\), and we assume that there are no edges connecting a vertex with itself (\(x_i\neq y_i\)) and there is at most one edge between any pair of vertices. The neighbors of a vertex \(x\) are the vertices \(y\neq x\) such that \(\{x,y\}\in E\). The degree \(d_x\) of a vertex \(x\) is its number of neighbours. An example is the following graph:
edges = [(0,1),(1,2),(0,3),(1,4),(2,5),(3,4),(4,5),(3,6),(4,7),(5,8),(6,7),(7,8),(5,7),(0,4)]
example_graph = nx.Graph(edges)
nx.draw(example_graph,with_labels=True)
A natural proposal transition matrix is \(Q(x \to y) = \frac{1}{d_x} \mathbf{1}_{\{\{x,y\}\in E\}}\). In other words, when at \(x\) the proposed next state is chosen uniformly among its neighbors.
(a) Write a function sample_proposal that, given a (networkX) graph and node \(x\), samples \(y\) according to transition matrix \(Q(x \to y)\). Hint: a useful Graph member function is neighbors. (10 pts)
def sample_proposal(graph,x):
# YOUR CODE HERE
raise NotImplementedError()
from nose.tools import assert_almost_equal
assert sample_proposal(nx.Graph([(0,1)]),0)==1
assert_almost_equal([sample_proposal(example_graph,3) for _ in range(1000)].count(4),333,delta=50)
assert_almost_equal([sample_proposal(example_graph,8) for _ in range(1000)].count(5),500,delta=60)
(b) Let us consider the Markov chain corresponding to the transition matrix \(Q(x \to y)\). Produce a histogram of the states visited in the first ~20000 steps. Compare this to the exact stationary distribution found by the function stationary_distributions from the lecture applied to the transition matrix \(Q\). Hint: another useful Graph member function is degree. (15 pts)
def chain_Q_histogram(graph,start,k):
'''Produce a histogram (a Numpy array of length equal to the number of
nodes of graph) of the states visited (excluding initial state) by the
Q Markov chain in the first k steps when started at start.'''
# YOUR CODE HERE
raise NotImplementedError()
def transition_matrix_Q(graph):
'''Construct transition matrix Q from graph as two-dimensional Numpy array.'''
# YOUR CODE HERE
raise NotImplementedError()
# Compare histogram and stationary distribution in a plot
# YOUR CODE HERE
raise NotImplementedError()
assert_almost_equal(transition_matrix_Q(example_graph)[3,4],1/3,delta=1e-9)
assert_almost_equal(transition_matrix_Q(example_graph)[3,7],0.0,delta=1e-9)
assert_almost_equal(transition_matrix_Q(example_graph)[2,2],0.0,delta=1e-9)
assert_almost_equal(np.sum(transition_matrix_Q(example_graph)[7]),1.0,delta=1e-9)
assert chain_Q_histogram(nx.Graph([(0,1)]),0,100)[1] == 50
assert len(chain_Q_histogram(example_graph,0,100)) == example_graph.number_of_nodes()
(c) Determine the appropriate Metropolis-Hastings acceptance probability \(A(x \to y)\) for \(x\neq y\) and write a function that, given a graph and \(x\), samples the next state with \(y\) according to the Metropolis-Hastings transition matrix \(P(x \to y)\). (10 pts)
def acceptance_probability(graph,x,y):
'''Compute A(x -> y) for the supplied graph (assuming x!=y).'''
# YOUR CODE HERE
raise NotImplementedError()
def sample_next_state(graph,x):
'''Return next random state y according to MH transition matrix P(x -> y).'''
# YOUR CODE HERE
raise NotImplementedError()
assert_almost_equal(acceptance_probability(example_graph,3,4),0.6,delta=1e-9)
assert_almost_equal(acceptance_probability(example_graph,8,7),0.5,delta=1e-9)
assert_almost_equal(acceptance_probability(example_graph,7,8),1.0,delta=1e-9)
assert_almost_equal(acceptance_probability(nx.Graph([(0,1)]),0,1),1,delta=1e-9)
(d) Do the same as in part (b) but now for the Markov chain corresponding to \(P\). Verify that the histogram of the Markov chain approaches a flat distribution and corroborate this by calculating the explicit matrix \(P\) and applying stationary_distributions to it. Hint: for determining the explicit matrix \(P(x\to y)\), remember that the formula \(P(x\to y) = Q(x\to y)A(x\to y)\) only holds for \(x\neq y\). What is \(P(x\to x)\)? (15 pts)
def chain_P_histogram(graph,start,n):
'''Produce a histogram of the states visited (excluding initial state)
by the P Markov chain in the first n steps when started at start.'''
# YOUR CODE HERE
raise NotImplementedError()
def transition_matrix_P(graph):
'''Construct transition matrix P from graph as numpy array.'''
# YOUR CODE HERE
raise NotImplementedError()
# plotting
# YOUR CODE HERE
raise NotImplementedError()
assert_almost_equal(transition_matrix_P(example_graph)[3,4],1/5,delta=1e-9)
assert_almost_equal(transition_matrix_P(example_graph)[3,7],0.0,delta=1e-9)
assert_almost_equal(transition_matrix_P(example_graph)[2,2],0.41666666,delta=1e-5)
assert_almost_equal(np.sum(transition_matrix_P(example_graph)[7]),1.0,delta=1e-9)
assert len(chain_P_histogram(example_graph,0,100)) == example_graph.number_of_nodes()
assert_almost_equal(chain_P_histogram(example_graph,0,20000)[8],2222,delta=180)
4.2. MCMC simulation of disk model#
(50 points)
Recall that in the disk model with we would like to sample the positions \(x = (x_1,y_1,\ldots,x_N,y_N)\in [0,L)^{2N}\) of \(N\) disks of radius \(1\) in the torus \([0,L)^2\) with uniform density \(\pi(x) = \mathbf{1}_{\{\text{all pairwise distance }\geq 2\}}(x) / Z\), where \(Z\) is the unknown partition function of the model. We will assume \(L > 2\) and \(N\geq 1\). For the purposes of this simulation we will store the state \(x\) in a np.array of dimension \((N,2)\) with values in \([0,L)\). Such a configuration can be conveniently plotted using the following function:
def plot_disk_configuration(positions,L):
fig,ax = plt.subplots()
ax.set_aspect('equal')
ax.set_ylim(0,L)
ax.set_xlim(0,L)
ax.set_yticklabels([])
ax.set_xticklabels([])
ax.set_yticks([])
ax.set_xticks([])
for x,y in positions:
# consider all horizontal and vertical copies that may be visible
for x_shift in [z for z in x + [-L,0,L] if -1<z<L+1]:
for y_shift in [z for z in y + [-L,0,L] if -1<z<L+1]:
ax.add_patch(plt.Circle((x_shift,y_shift),1))
plt.show()
# Example with N=3 and L=5
positions = np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]])
plot_disk_configuration(positions,5)
(a) Write a function two_disks_overlap that tests whether disks at position \(\mathbf{x}_1 = (x_1,y_1) \in [0,L)^{2}\) and position \(\mathbf{x}_2 = (x_2,y_2) \in [0,L)^{2}\) overlap and a function disk_config_valid that checks whether a full configuration is valid (non-overlapping and non-touching).
Hint: The distance between the two disks is \(\sqrt{\delta_x^2+\delta_y^2}\), where \(\delta_x\) and \(\delta_y\) are the minimal horizontal and vertical separation. Note that \(\delta_x = |x_1-x_2 - r|\) where \(r\) is the closest integer multiple of \(L\) to \(x_1-x_2\). (15 pts)
def two_disks_overlap(x1,x2,L):
'''Return True if the disks centered at x1 and x2 (represented as 2-element arrays) overlap in [0,L)^2.'''
# YOUR CODE HERE
raise NotImplementedError()
def disk_config_valid(x,L):
'''Return True if the configuration x (as two-dimensional array) is non-overlapping in [0,L)^2.'''
# YOUR CODE HERE
raise NotImplementedError()
assert two_disks_overlap(np.array([1,1]),np.array([1,1]),5)
assert two_disks_overlap(np.array([0.6,0.6]),np.array([4.1,0.5]),5)
assert two_disks_overlap(np.array([0.3,0.3]),np.array([4.6,4.6]),5)
assert not two_disks_overlap(np.array([1,1]),np.array([3.1,1]),7)
assert not two_disks_overlap(np.array([1,1]),np.array([1,3.1]),7)
assert not two_disks_overlap(np.array([1,1]),np.array([1.01+np.sqrt(2),1.01+np.sqrt(2)]),6)
assert two_disks_overlap(np.array([1,1]),np.array([0.99+np.sqrt(2),0.99+np.sqrt(2)]),6)
assert disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]]),5)
assert not disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4],[4.1,2.3]]),5)
assert disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1]]),6)
assert not disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1],[2.5,2.5]]),6)
(b) Assuming \(N \leq \lceil \frac12 L -1 \rceil^2\) where \(\lceil r\rceil\) is the smallest integer larger or equal to \(r\), write a function generate_initial_positions that produces an arbitrary non-overlapping (and non-touching) initial condition given \(N\) and \(L\). The layout need not be random, any deterministic layout is ok (e.g. grid). (10 pts)
def generate_initial_positions(N,L):
'''Return array of positions of N disks in non-overlapping positions.'''
# YOUR CODE HERE
raise NotImplementedError()
plot_disk_configuration(generate_initial_positions(33,14.5),14.5)
for l in [4,9.2,14.5]:
for n in range(1,int(np.ceil(l/2-1)**2)+1):
assert disk_config_valid(generate_initial_positions(n,l),l), "Failed for n = {}, l = {}".format(n,l)
(c) Write a function remains_valid_after_move that determines whether in a non-overlapping configuration \(x\) moving the \(i\)th disk to next_position results in a valid non-overlapping configuration. (10 pts)
def remains_valid_after_move(x,i,next_position,L):
'''Returns True if replacing x[i] by next_position would yield a valid configuration,
otherwise False.'''
# YOUR CODE HERE
raise NotImplementedError()
assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],0,[4.5,0.5],5)
assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],1,[4.5,0.5],5)
assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,2.5],5)
assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,3.8],5)
(d) Implement the Metropolis-Hastings transition by selecting a uniformly chosen disk and displacing it by \( (\delta\,\mathcal{N}_1,\delta\,\mathcal{N}_2)\) where \(\delta>0\) is a parameter and \(\mathcal{N}_i\) are independent normal random variables (make sure to keep positions within \([0,L)^2\) by taking the new position modulo \(L\)). Test run your simulation for \(L=11.3\) and \(N=20\) and \(\delta = 0.3\) and about \(10000\) Markov chain steps and plot the final state. (15 pts)
def MH_disk_move(x,L,delta):
'''Perform random MH move on configuration x, thus changing the array x (if accepted).
Return True if move was accepted, False otherwise.'''
# YOUR CODE HERE
raise NotImplementedError()
# Test run and plot resulting configuration
# YOUR CODE HERE
raise NotImplementedError()