# 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 says YOUR CODE HERE or YOUR ANSWER HERE and remove the raise 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 All in the jupyter menu.

• For some exercises test cases have been provided in a separate cell in the form of assert statements. When run, a successful test will give no output, whereas a failed test will display an error message.

• 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.

NAME = ""
NAMES_OF_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]


## 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):
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.'''
raise NotImplementedError()

def transition_matrix_Q(graph):
'''Construct transition matrix Q from graph as two-dimensional Numpy array.'''
raise NotImplementedError()

# Compare histogram and stationary distribution in a plot
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).'''
raise NotImplementedError()

def sample_next_state(graph,x):
'''Return next random state y according to MH transition matrix P(x -> y).'''
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.'''
raise NotImplementedError()

def transition_matrix_P(graph):
'''Construct transition matrix Q from graph as numpy array.'''
raise NotImplementedError()

# plotting
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]:
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 \in [0,L)^{2}$$ and position $$\mathbf{x}_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 minimal separation in the $$x$$-direction can be expressed as a function of x1[0]-x2[0] and the minimal separation in the y-direction as a function of x1[1]-x2[1]. Then use pythagoras. (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.'''
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.'''
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.'''
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.'''
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.'''