diff --git a/Exercise sheet 4/exercise_sheet_04.ipynb b/Exercise sheet 4/exercise_sheet_04.ipynb new file mode 100644 index 0000000..6f2ed37 --- /dev/null +++ b/Exercise sheet 4/exercise_sheet_04.ipynb @@ -0,0 +1,887 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7711bd97", + "metadata": {}, + "source": [ + "# Exercise sheet\n", + "\n", + "Some general remarks about the exercises:\n", + "* For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.\n", + "* 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()`.\n", + "* 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](https://www.python.org/dev/peps/pep-0008/) for the widely adopted coding conventions or [this guide for explanation](https://realpython.com/python-pep8/).\n", + "* 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.\n", + "* 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.\n", + "* 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.\n", + "\n", + "Please fill in your name here:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e39c230e", + "metadata": {}, + "outputs": [], + "source": [ + "NAME = \"\"\n", + "NAMES_OF_COLLABORATORS = \"\"" + ] + }, + { + "cell_type": "markdown", + "id": "2e7caa84", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "0c4606a8", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "1ba146ece4effd3e8dcf33fe63ec69dc", + "grade": false, + "grade_id": "cell-455926422ebe4f85", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "**Exercise sheet 4**\n", + "\n", + "Code from the lectures:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d3bb610", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "ac562d84531ae19e1d920c10a8861b8e", + "grade": false, + "grade_id": "cell-6b2d0465c2abaaa6", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pylab as plt\n", + "import networkx as nx\n", + "\n", + "rng = np.random.default_rng()\n", + "%matplotlib inline\n", + "\n", + "def draw_transition_graph(P):\n", + " # construct a directed graph directly from the matrix\n", + " graph = nx.DiGraph(P) \n", + " # draw it in such a way that edges in both directions are visible and have appropriate width\n", + " nx.draw_networkx(graph,connectionstyle='arc3, rad = 0.15',width=[6*P[u,v] for u,v in graph.edges()])\n", + "\n", + "def sample_next(P,current):\n", + " return rng.choice(len(P),p=P[current])\n", + "\n", + "def sample_chain(P,start,n):\n", + " chain = [start]\n", + " for _ in range(n):\n", + " chain.append(sample_next(P,chain[-1]))\n", + " return chain\n", + "\n", + "def stationary_distributions(P):\n", + " eigenvalues, eigenvectors = np.linalg.eig(np.transpose(P))\n", + " # make list of normalized eigenvectors for which the eigenvalue is very close to 1\n", + " return [eigenvectors[:,i]/np.sum(eigenvectors[:,i]) for i in range(len(eigenvalues)) \n", + " if np.abs(eigenvalues[i]-1) < 1e-10]\n", + "\n", + "def markov_sample_mean(P,start,function,n):\n", + " total = 0\n", + " state = start\n", + " for _ in range(n):\n", + " state = sample_next(P,state)\n", + " total += function[state]\n", + " return total/n" + ] + }, + { + "cell_type": "markdown", + "id": "2f7814d8", + "metadata": {}, + "source": [ + "## Markov Chain on a graph\n", + "\n", + "**(50 points)**\n", + "\n", + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "071e5617", + "metadata": {}, + "outputs": [], + "source": [ + "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)]\n", + "example_graph = nx.Graph(edges)\n", + "nx.draw(example_graph,with_labels=True)" + ] + }, + { + "cell_type": "markdown", + "id": "952c97cd", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "a9ae9f2c0eb7e99b85fd90ae657a9202", + "grade": false, + "grade_id": "cell-a8413209ddb70ab0", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "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.\n", + "\n", + "__(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`](https://networkx.org/documentation/stable/reference/classes/generated/networkx.Graph.neighbors.html). **(10 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "449da919", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "db0642dd7c9f53225d13fb9674c89a56", + "grade": false, + "grade_id": "cell-df840ef2dee49b9f", + "locked": false, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def sample_proposal(graph,x):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "604ce3f4", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "690a774380e6a4383ccb5da2dee8a737", + "grade": true, + "grade_id": "cell-c70fec7da167b7be", + "locked": true, + "points": 10, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "from nose.tools import assert_almost_equal\n", + "assert sample_proposal(nx.Graph([(0,1)]),0)==1\n", + "assert_almost_equal([sample_proposal(example_graph,3) for _ in range(1000)].count(4),333,delta=50)\n", + "assert_almost_equal([sample_proposal(example_graph,8) for _ in range(1000)].count(5),500,delta=60)" + ] + }, + { + "cell_type": "markdown", + "id": "f7a2e658", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "99c1e08eed2976f4625452d9eed6a5d5", + "grade": false, + "grade_id": "cell-4c74e7dae57e50e3", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(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`](https://networkx.org/documentation/stable/reference/classes/generated/networkx.Graph.degree.html). **(15 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f52fcbd5", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "44377fd9fa8cd0c3bc78ab03753db1ce", + "grade": false, + "grade_id": "cell-ca4f5c3685b72c8a", + "locked": false, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def chain_Q_histogram(graph,start,k):\n", + " '''Produce a histogram (a Numpy array of length equal to the number of \n", + " nodes of graph) of the states visited (excluding initial state) by the \n", + " Q Markov chain in the first k steps when started at start.'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "def transition_matrix_Q(graph):\n", + " '''Construct transition matrix Q from graph as two-dimensional Numpy array.'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "# Compare histogram and stationary distribution in a plot\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "812cd2ef", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "c3e08e0577c953ca68fe8159cadc0634", + "grade": true, + "grade_id": "cell-bdea40714e10d0e8", + "locked": true, + "points": 15, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "assert_almost_equal(transition_matrix_Q(example_graph)[3,4],1/3,delta=1e-9)\n", + "assert_almost_equal(transition_matrix_Q(example_graph)[3,7],0.0,delta=1e-9)\n", + "assert_almost_equal(transition_matrix_Q(example_graph)[2,2],0.0,delta=1e-9)\n", + "assert_almost_equal(np.sum(transition_matrix_Q(example_graph)[7]),1.0,delta=1e-9)\n", + "assert chain_Q_histogram(nx.Graph([(0,1)]),0,100)[1] == 50\n", + "assert len(chain_Q_histogram(example_graph,0,100)) == example_graph.number_of_nodes()" + ] + }, + { + "cell_type": "markdown", + "id": "8ebe4946", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "f4338397c60be0fc0069e30bec121faa", + "grade": false, + "grade_id": "cell-5a9debbf863b9e3b", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(c)__ Determine the appropriate Metropolis-Hastings acceptance probability $A(x \\to y)$ for $x\\neq y$\n", + "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)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37804448", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "639522f9ff68136165c0d9daa174b4ba", + "grade": false, + "grade_id": "cell-92e3b11a3af7eb99", + "locked": false, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def acceptance_probability(graph,x,y):\n", + " '''Compute A(x -> y) for the supplied graph (assuming x!=y).'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "def sample_next_state(graph,x):\n", + " '''Return next random state y according to MH transition matrix P(x -> y).'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "058d3728", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "6bf922f8787f02c80681390ce8a8f84c", + "grade": true, + "grade_id": "cell-fbfc2607999d0c99", + "locked": true, + "points": 10, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "assert_almost_equal(acceptance_probability(example_graph,3,4),0.6,delta=1e-9)\n", + "assert_almost_equal(acceptance_probability(example_graph,8,7),0.5,delta=1e-9)\n", + "assert_almost_equal(acceptance_probability(example_graph,7,8),1.0,delta=1e-9)\n", + "assert_almost_equal(acceptance_probability(nx.Graph([(0,1)]),0,1),1,delta=1e-9)" + ] + }, + { + "cell_type": "markdown", + "id": "cbf205bf", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "da69b710a3b1986ac9f9997235c1a472", + "grade": false, + "grade_id": "cell-97f7d289177acf39", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(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)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5de68351", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "bd25c9f6b39133973fd2d1594e210b30", + "grade": false, + "grade_id": "cell-3b7fde395331916d", + "locked": false, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def chain_P_histogram(graph,start,n):\n", + " '''Produce a histogram of the states visited (excluding initial state) \n", + " by the P Markov chain in the first n steps when started at start.'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "def transition_matrix_P(graph):\n", + " '''Construct transition matrix Q from graph as numpy array.'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "# plotting\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b5cc504", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "1b756f8dff589048e791db63f89d8624", + "grade": true, + "grade_id": "cell-8c4b6c60ee96f037", + "locked": true, + "points": 10, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "assert_almost_equal(transition_matrix_P(example_graph)[3,4],1/5,delta=1e-9)\n", + "assert_almost_equal(transition_matrix_P(example_graph)[3,7],0.0,delta=1e-9)\n", + "assert_almost_equal(transition_matrix_P(example_graph)[2,2],0.41666666,delta=1e-5)\n", + "assert_almost_equal(np.sum(transition_matrix_P(example_graph)[7]),1.0,delta=1e-9)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4e9d8aa", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "22ff259ba7ad6066040f3417cc8d4483", + "grade": true, + "grade_id": "cell-a63274314d1b2a2d", + "locked": true, + "points": 5, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "assert len(chain_P_histogram(example_graph,0,100)) == example_graph.number_of_nodes()\n", + "assert_almost_equal(chain_P_histogram(example_graph,0,20000)[8],2222,delta=180)" + ] + }, + { + "cell_type": "markdown", + "id": "89f65139", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "9b277d17c7b55bb15e5708e6ca90c4f9", + "grade": false, + "grade_id": "cell-da50533d4d46c293", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## MCMC simulation of disk model\n", + "\n", + "**(50 points)**\n", + "\n", + "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:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "683d2de3", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "4c7cb2d90103cd790145cae15d46e99e", + "grade": false, + "grade_id": "cell-d661f575ab9f80ea", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "def plot_disk_configuration(positions,L):\n", + " fig,ax = plt.subplots()\n", + " ax.set_aspect('equal')\n", + " ax.set_ylim(0,L)\n", + " ax.set_xlim(0,L)\n", + " ax.set_yticklabels([])\n", + " ax.set_xticklabels([])\n", + " ax.set_yticks([])\n", + " ax.set_xticks([])\n", + " for x,y in positions:\n", + " # consider all horizontal and vertical copies that may be visible\n", + " for x_shift in [z for z in x + [-L,0,L] if -10$ 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)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ccaa8f9f", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "24f09b7137c0184ddb872ef558f756a9", + "grade": true, + "grade_id": "cell-e180b99ca699c610", + "locked": false, + "points": 15, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def MH_disk_move(x,L,delta):\n", + " '''Perform random MH move on configuration x, thus changing the array x (if accepted). \n", + " Return True if move was accepted, False otherwise.'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + "\n", + "# Test run and plot resulting configuration\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}