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.

  • 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()