kees-study/cds-monte-carlo-methods
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

04: Fetch week 4

Browse Source
This commit is contained in:
Kees van Kempen
2022-09-29 15:41:01 +02:00
parent 24b3f3202b
commit b3eb03d239
1 changed files with 887 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

887
Exercise sheet 4/exercise_sheet_04.ipynb Normal file
View File

@ -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 -1<z<L+1]:\n",
" for y_shift in [z for z in y + [-L,0,L] if -1<z<L+1]:\n",
" ax.add_patch(plt.Circle((x_shift,y_shift),1))\n",
" plt.show()\n",
" \n",
"# Example with N=3 and L=5\n",
"positions = np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]])\n",
"plot_disk_configuration(positions,5)"
]
},
{
"cell_type": "markdown",
"id": "efa2f1c9",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "d29e40f9d0e8d304d266f016150aeab7",
"grade": false,
"grade_id": "cell-527056f298dedf7c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"__(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)**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19c010eb",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "bbf8ae650d9d4a66a3e4ac4722a3f60e",
"grade": false,
"grade_id": "cell-1b2a61bf719003e0",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def two_disks_overlap(x1,x2,L):\n",
" '''Return True if the disks centered at x1 and x2 (represented as 2-element arrays) overlap in [0,L)^2.'''\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" \n",
"def disk_config_valid(x,L):\n",
" '''Return True if the configuration x (as two-dimensional array) is non-overlapping in [0,L)^2.'''\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e0f2024a",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "05fc821fed0690aabc13fd66da800694",
"grade": true,
"grade_id": "cell-9f544deda0526691",
"locked": true,
"points": 10,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert two_disks_overlap(np.array([1,1]),np.array([1,1]),5)\n",
"assert two_disks_overlap(np.array([0.6,0.6]),np.array([4.1,0.5]),5)\n",
"assert two_disks_overlap(np.array([0.3,0.3]),np.array([4.6,4.6]),5)\n",
"assert not two_disks_overlap(np.array([1,1]),np.array([3.1,1]),7)\n",
"assert not two_disks_overlap(np.array([1,1]),np.array([1,3.1]),7)\n",
"assert not two_disks_overlap(np.array([1,1]),np.array([1.01+np.sqrt(2),1.01+np.sqrt(2)]),6)\n",
"assert two_disks_overlap(np.array([1,1]),np.array([0.99+np.sqrt(2),0.99+np.sqrt(2)]),6)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "68e57e57",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "8705cb3a08fc636963fe61806348d093",
"grade": true,
"grade_id": "cell-699454de327d56d5",
"locked": true,
"points": 5,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]]),5)\n",
"assert not disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4],[4.1,2.3]]),5)\n",
"assert disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1]]),6)\n",
"assert not disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1],[2.5,2.5]]),6)"
]
},
{
"cell_type": "markdown",
"id": "fe7c5736",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "db023e65fc2613d624f29977a82bf1f1",
"grade": false,
"grade_id": "cell-b0e00f11b0f85525",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"__(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)**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f86a04ed",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "b17579d74c8ff8572d81f18997b49c90",
"grade": false,
"grade_id": "cell-53c1bd894d6fe27d",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def generate_initial_positions(N,L):\n",
" '''Return array of positions of N disks in non-overlapping positions.'''\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()\n",
" \n",
"plot_disk_configuration(generate_initial_positions(33,14.5),14.5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3d3dc3db",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "fa73c21133aa2b913b931b266e7782f0",
"grade": true,
"grade_id": "cell-e1c80d8b59301b93",
"locked": true,
"points": 10,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"for l in [4,9.2,14.5]:\n",
" for n in range(1,int(np.ceil(l/2-1)**2)+1):\n",
" assert disk_config_valid(generate_initial_positions(n,l),l), \"Failed for n = {}, l = {}\".format(n,l)"
]
},
{
"cell_type": "markdown",
"id": "eb5dbe07",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "b4bb37c92263c95c89055362336f894a",
"grade": false,
"grade_id": "cell-1bd76964ec9620ce",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"__(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)**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "56252046",
"metadata": {
"deletable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "04ff82e1b0e3386f2cb15412695c1864",
"grade": false,
"grade_id": "cell-d54b4fa9b2f8eb92",
"locked": false,
"schema_version": 3,
"solution": true,
"task": false
}
},
"outputs": [],
"source": [
"def remains_valid_after_move(x,i,next_position,L):\n",
" '''Returns True if replacing x[i] by next_position would yield a valid configuration,\n",
" otherwise False.'''\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "167bdc3e",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "code",
"checksum": "37fe75b415d4e31369d6b6bb1041609b",
"grade": true,
"grade_id": "cell-e4902309b9869f3d",
"locked": true,
"points": 10,
"schema_version": 3,
"solution": false,
"task": false
}
},
"outputs": [],
"source": [
"assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],0,[4.5,0.5],5)\n",
"assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],1,[4.5,0.5],5)\n",
"assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,2.5],5)\n",
"assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,3.8],5)"
]
},
{
"cell_type": "markdown",
"id": "a8cd0c44",
"metadata": {
"deletable": false,
"editable": false,
"nbgrader": {
"cell_type": "markdown",
"checksum": "8e0cfe7c5bce12cc34e939c5c151ed65",
"grade": false,
"grade_id": "cell-582ca3e158908e6c",
"locked": true,
"schema_version": 3,
"solution": false,
"task": false
}
},
"source": [
"__(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)**"
]
},
{
"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
}