From 0bc6c0c371ce5eb909150f0cafdbb4ec48f1e5a7 Mon Sep 17 00:00:00 2001
From: Kees van Kempen <k.vankempen@student.science.ru.nl>
Date: Thu, 8 Sep 2022 10:36:26 +0200
Subject: [PATCH] Initial commit: fetch 01

---
 Exercise sheet 1/exercise_sheet_01.ipynb | 627 +++++++++++++++++++++++
 1 file changed, 627 insertions(+)
 create mode 100644 Exercise sheet 1/exercise_sheet_01.ipynb

diff --git a/Exercise sheet 1/exercise_sheet_01.ipynb b/Exercise sheet 1/exercise_sheet_01.ipynb
new file mode 100644
index 0000000..51540b2
--- /dev/null
+++ b/Exercise sheet 1/exercise_sheet_01.ipynb	
@@ -0,0 +1,627 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "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",
+    "* 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,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "NAME = \"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "42a51f89f2fde456dc33fa324c98aed9",
+     "grade": false,
+     "grade_id": "cell-f051f1e1b3b3e7be",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__Exercise sheet 1__\n",
+    "\n",
+    "Code from the lecture:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "b99ee0ea4e023a0a82dc338192363b57",
+     "grade": false,
+     "grade_id": "cell-9b66485da15cd77c",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt\n",
+    "\n",
+    "rng = np.random.default_rng()\n",
+    "%matplotlib inline\n",
+    "\n",
+    "def random_in_square():\n",
+    "    \"\"\"Returns a random position in the square [-1,1)x[-1,1).\"\"\"\n",
+    "    return rng.uniform(-1,1,2)  \n",
+    "\n",
+    "def is_in_circle(x):\n",
+    "    return np.dot(x,x) < 1\n",
+    "\n",
+    "def simulate_number_of_hits(N):\n",
+    "    \"\"\"Simulates number of hits in case of N trials in the pebble game.\"\"\"\n",
+    "    number_hits = 0\n",
+    "    for i in range(N):\n",
+    "        position = random_in_square()\n",
+    "        if is_in_circle(position):\n",
+    "            number_hits += 1\n",
+    "    return number_hits\n",
+    "\n",
+    "def random_in_disk():\n",
+    "    \"\"\"Returns a uniform point in the unit disk via rejection.\"\"\"\n",
+    "    position = random_in_square()\n",
+    "    while not is_in_circle(position):\n",
+    "        position = random_in_square()\n",
+    "    return position\n",
+    "\n",
+    "def is_in_square(x):\n",
+    "    \"\"\"Returns True if x is in the square (-1,1)^2.\"\"\"\n",
+    "    return np.abs(x[0]) < 1 and np.abs(x[1]) < 1\n",
+    "    \n",
+    "def sample_next_position_naively(position,delta):\n",
+    "    \"\"\"Keep trying a throw until it ends up in the square.\"\"\"\n",
+    "    while True:\n",
+    "        next_position = position + delta*random_in_disk()\n",
+    "        if is_in_square(next_position):\n",
+    "            return next_position\n",
+    "    \n",
+    "def naive_markov_pebble(start,delta,N):\n",
+    "    \"\"\"Simulates the number of hits in the naive Markov-chain version \n",
+    "    of the pebble game.\"\"\"\n",
+    "    number_hits = 0\n",
+    "    position = start\n",
+    "    for i in range(N):\n",
+    "        position = sample_next_position_naively(position,delta)\n",
+    "        if is_in_circle(position):\n",
+    "            number_hits += 1\n",
+    "    return number_hits\n",
+    "\n",
+    "def naive_markov_pebble_generator(start,delta,N):\n",
+    "    \"\"\"Same as naive_markov_pebble but only yields the positions.\"\"\"\n",
+    "    position = start\n",
+    "    for i in range(N):\n",
+    "        position = sample_next_position_naively(position,delta)\n",
+    "        yield position\n",
+    "\n",
+    "def sample_next_position(position,delta):\n",
+    "    \"\"\"Attempt a throw and reject when outside the square.\"\"\"\n",
+    "    next_position = position + delta*random_in_disk()\n",
+    "    if is_in_square(next_position):\n",
+    "        return next_position  # accept!\n",
+    "    else:\n",
+    "        return position  # reject!\n",
+    "\n",
+    "def markov_pebble(start,delta,N):\n",
+    "    \"\"\"Simulates the number of hits in the proper Markov-chain version of the pebble game.\"\"\"\n",
+    "    number_hits = 0\n",
+    "    position = start\n",
+    "    for i in range(N):\n",
+    "        position = sample_next_position(position,delta)\n",
+    "        if is_in_circle(position):\n",
+    "            number_hits += 1\n",
+    "    return number_hits\n",
+    "\n",
+    "def markov_pebble_generator(start,delta,N):\n",
+    "    \"\"\"Same as markov_pebble but only yields the positions.\"\"\"\n",
+    "    position = start\n",
+    "    for i in range(N):\n",
+    "        position = sample_next_position(position,delta)\n",
+    "        yield position"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "48539a8a121b164023c052bb11a71ca3",
+     "grade": false,
+     "grade_id": "cell-99420e260f8510fe",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "## Empirical convergence rate in the pebble game\n",
+    "\n",
+    "__(a)__ Write a function `pi_stddev` that estimates the standard deviation $\\sigma$ of the $\\pi$-estimate using $n$ trials by running the direct-sampling pebble game $m$ times. Store this data for $n=2^4,2^5,\\ldots,2^{14}$ and $m=200$ in an array. *Hint*: you may use the NumPy function [`np.std`](https://numpy.org/doc/stable/reference/generated/numpy.std.html). __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "6896cc4d2b06bc6cd9db88bbd93fb488",
+     "grade": false,
+     "grade_id": "cell-0e3910cb420ef878",
+     "locked": false,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def pi_stddev(n,m):\n",
+    "    \"\"\"Estimate the standard deviation in the pi estimate in the case of n trials,\n",
+    "    based on m runs of the direct-sampling pebble game.\"\"\"\n",
+    "    # YOUR CODE HERE\n",
+    "    raise NotImplementedError()\n",
+    "    \n",
+    "stddev_data = np.array([[2**k,pi_stddev(2**k,200)] for k in range(4,12)])\n",
+    "stddev_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "9eec8b91bf40df9b97932746f8682097",
+     "grade": true,
+     "grade_id": "cell-90b5022604409d60",
+     "locked": true,
+     "points": 12,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# If done correctly, your code should pass the following tests\n",
+    "from nose.tools import assert_almost_equal\n",
+    "assert_almost_equal(pi_stddev(2**3,1000),0.582,delta=0.03)\n",
+    "assert_almost_equal(pi_stddev(2**4,1000),0.411,delta=0.03)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "bea3c12c013efeb34a936fe5be28149c",
+     "grade": false,
+     "grade_id": "cell-7103a31a2719c5a0",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__(b)__ Write a function `fit_power_law` that takes an array of $(n,\\sigma)$ pairs and determines best-fit parameters $a,p$ for the curve $\\sigma = a n^p$. This is best done by fitting a straight line to the data on a log-log scale ($\\log\\sigma=\\log a+p\\log n$). *Hint*: use [`curve_fit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) from SciPy. __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "7ef0a128432a026f4caa85746fd729b6",
+     "grade": false,
+     "grade_id": "cell-8656545dbb5ec948",
+     "locked": false,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def fit_power_law(stddev_data):\n",
+    "    \"\"\"Compute the best fit parameters a and p.\"\"\"\n",
+    "    # YOUR CODE HERE\n",
+    "    raise NotImplementedError()\n",
+    "    return a_fit, p_fit"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "1983b7014d8bac9a28f6fd36fddb3589",
+     "grade": true,
+     "grade_id": "cell-311260afe5628598",
+     "locked": true,
+     "points": 12,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from nose.tools import assert_almost_equal\n",
+    "assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[0],3.36,delta=0.05)\n",
+    "assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[1],-0.31,delta=0.03)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "63cac93ff9a1399db5b3eea45c7c6ff4",
+     "grade": false,
+     "grade_id": "cell-5c9e6c8197ea6ab5",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__(c)__ Plot the data against the fitted curve $\\sigma = a n^p$ in a log-log plot. Don't forget to properly label your axes. *Hint*: use [`loglog`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.loglog.html) from matplotlib. __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "df285cc68ad2e9230aee61ae26cda962",
+     "grade": true,
+     "grade_id": "cell-1755ac9d1fb8402f",
+     "locked": false,
+     "points": 12,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# YOUR CODE HERE\n",
+    "raise NotImplementedError()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "efacdbc9ad60113c9fd39e149f09e955",
+     "grade": false,
+     "grade_id": "cell-8be1d8340bd0c065",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "## Volume of a unit ball in other dimensions\n",
+    "\n",
+    "In this exercise you will adapt the direct-sampling Monte Carlo method of the pebble game to higher dimensions to\n",
+    "estimate the $d$-dimensional volume of the $d$-dimensional ball of radius $1$.\n",
+    "\n",
+    "__(a)__ Adapt the direct-sampling Monte Carlo method to the pebble game in a $d$-dimensional hypercube $(-1,1)^d$ with an inscribed $d$-dimensional unit ball. __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "71be71d20de0f6ee438696d00c9306e8",
+     "grade": false,
+     "grade_id": "cell-120de066cebda898",
+     "locked": false,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def number_of_hits_in_d_dimensions(N,d):\n",
+    "    \"\"\"Simulates number of hits in case of N trials in the d-dimensional direct-sampling pebble game.\"\"\"\n",
+    "    # YOUR CODE HERE\n",
+    "    raise NotImplementedError()\n",
+    "    return number_hits"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "521b242a339c4a797d346f192ac1fe2d",
+     "grade": true,
+     "grade_id": "cell-e8e97677ad30041e",
+     "locked": true,
+     "points": 12,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from nose.tools import assert_almost_equal\n",
+    "assert_almost_equal(number_of_hits_in_d_dimensions(100,1),100,delta=1)\n",
+    "assert_almost_equal(number_of_hits_in_d_dimensions(2000,3),1045,delta=50)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "fb64e8b2afc2e451908b8743d83a7ea6",
+     "grade": false,
+     "grade_id": "cell-2b9ff58ee25d374f",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__(b)__ Compare your estimates for the volume $V_d$ of the $d$-dimensional unit ball for $N=10000$ trials and $d=1,2,\\ldots,7$ to the [exact formula](https://en.wikipedia.org/wiki/Volume_of_an_n-ball) $V_d = \\frac{\\pi^{d/2}}{\\Gamma(\\tfrac{d}{2}+1)}$ in a plot. *Hint*: the Gamma function is available in `scipy.special.gamma`. __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "68784cc26ebb01389c1ad1675f0e5ef8",
+     "grade": true,
+     "grade_id": "cell-82ea663199a69ee2",
+     "locked": false,
+     "points": 12,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# YOUR CODE HERE\n",
+    "raise NotImplementedError()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "0eff565d5e6284f83691c38b751d3ce8",
+     "grade": false,
+     "grade_id": "cell-beac1a536ba9895a",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "## Efficiency of the Metropolis algorithm and the 1/2-rule\n",
+    "\n",
+    "In the lecture we mentioned the \"1/2 rule\" for Metropolis algorithms: if moves are rarely rejected then typically you do not explore the domain efficiently, while if most moves are rejected you are wasting efforts proposing moves. A good rule of thumb then is to aim for a rejection rate of $1/2$. In this exercise you are asked to test whether this rule of thumb makes any sense for the Markov-chain pebble game by varying the throwing range $\\delta$. \n",
+    "\n",
+    "__(a)__ Estimate the mean square deviation $\\mathbb{E}[(\\tfrac{4\\text{hits}}{\\text{trials}} - \\pi)^2 ]$ from $\\pi$ for different values of $\\delta$ ranging between $0$ and $3$, but fixed number of trials (say, 2000). For a decent estimate of the mean you will need at least $100$ repetitions. __(14 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "836a5749c3515998c0bf7a7b524886b4",
+     "grade": true,
+     "grade_id": "cell-e3762269c892c3df",
+     "locked": false,
+     "points": 14,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# YOUR CODE HERE\n",
+    "raise NotImplementedError()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "1819e4933a0845a4463f559ed8613bb2",
+     "grade": false,
+     "grade_id": "cell-2b9b2164827f2cb9",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__(b)__ Measure the rejection rate for the same simulations as in (a). __(14 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "ec4f5855aa1a96d41a81fe64542ff847",
+     "grade": true,
+     "grade_id": "cell-9726c572b9859251",
+     "locked": false,
+     "points": 14,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# YOUR CODE HERE\n",
+    "raise NotImplementedError()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "deletable": false,
+    "editable": false,
+    "nbgrader": {
+     "cell_type": "markdown",
+     "checksum": "bfc0b7ec0fcf893a10443ba881f99ff5",
+     "grade": false,
+     "grade_id": "cell-6b4d8983ff5d3d37",
+     "locked": true,
+     "schema_version": 3,
+     "solution": false,
+     "task": false
+    }
+   },
+   "source": [
+    "__(c)__ Plot both the mean square deviation and the rejection rate as function of $\\delta$. How well does the 1/2 rule apply in this situation? __(12 pts)__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "deletable": false,
+    "nbgrader": {
+     "cell_type": "code",
+     "checksum": "9ef86bd5fb28c1a681ea44586c49477a",
+     "grade": true,
+     "grade_id": "cell-0677a8cb2d72d30f",
+     "locked": false,
+     "points": 12,
+     "schema_version": 3,
+     "solution": true,
+     "task": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# 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"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}