From eb29fae24031f46afb71542f233875685165b341 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Thu, 8 Sep 2022 11:43:59 +0200 Subject: [PATCH] Ignore checkpoints --- .gitignore | 1 + .../exercise_sheet_01-checkpoint.ipynb | 676 ------------------ 2 files changed, 1 insertion(+), 676 deletions(-) create mode 100644 .gitignore delete mode 100644 Exercise sheet 1/.ipynb_checkpoints/exercise_sheet_01-checkpoint.ipynb diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..763513e --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.ipynb_checkpoints diff --git a/Exercise sheet 1/.ipynb_checkpoints/exercise_sheet_01-checkpoint.ipynb b/Exercise sheet 1/.ipynb_checkpoints/exercise_sheet_01-checkpoint.ipynb deleted file mode 100644 index 60d64e8..0000000 --- a/Exercise sheet 1/.ipynb_checkpoints/exercise_sheet_01-checkpoint.ipynb +++ /dev/null @@ -1,676 +0,0 @@ -{ - "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": 1, - "metadata": {}, - "outputs": [], - "source": [ - "NAME = \"Kees van Kempen\"" - ] - }, - { - "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": 2, - "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": 3, - "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": [ - { - "data": { - "text/plain": [ - "array([[1.60000000e+01, 4.14095400e-01],\n", - " [3.20000000e+01, 2.99924470e-01],\n", - " [6.40000000e+01, 2.19101280e-01],\n", - " [1.28000000e+02, 1.49149198e-01],\n", - " [2.56000000e+02, 9.15196946e-02],\n", - " [5.12000000e+02, 7.09785726e-02],\n", - " [1.02400000e+03, 5.33523935e-02],\n", - " [2.04800000e+03, 3.70486705e-02]])" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "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", - " \n", - " return np.std([simulate_number_of_hits(n)/n*4 for i in range(m)])\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": 4, - "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": 5, - "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": [ - "import scipy.optimize\n", - "def fit_power_law(stddev_data):\n", - " \"\"\"Compute the best fit parameters a and p.\"\"\"\n", - " \n", - " # Define the to be fitted function. We use the log, due to the calculation of the uncertainties.\n", - " log_f = lambda n, a, p: np.log(a) + p*np.log(n)\n", - " \n", - " n, 𝜎 = stddev_data.T\n", - " a_fit, p_fit = scipy.optimize.curve_fit(log_f, n, np.log(𝜎))[0]\n", - " \n", - " return a_fit, p_fit" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "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": 7, - "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": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEbCAYAAAAMKCkgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAr4ElEQVR4nO3dd3RU5d7F8e8vBUINvQhIVQRBinSkCgpKU1EQEEQEqYqKioqda8VGF0RBRUBB6dJFeAFpYqGINJGASpMiPeR5/5hwjbkJJCSTM5Psz1pZlzlzyp6Jd3ZOmeeYcw4RERF/CvE6gIiIpH8qGxER8TuVjYiI+J3KRkRE/E5lIyIifqeyERERv1PZiIiI36lsRIKQmZU1sw1mdtzMHoydtsnMGvpr/cHAzDaa2czY7I97nUf+YfpSpwQKM/sVKAhEA+eBzcBHwBjnXEwSl7/fObfIjzEDgpmNA4455x5O5PlfScF7can1ByIzywVsAioBx4A1zrnKXmaSf2jPRgJNS+dcDqA48CrwBDDO20hJZ2ZhabSp4vg+WIN1/f5QEZjknDvonDsLHPI6kPxDZSMByTl31Dk3E2gHdDGzCgBmNtDMdsQe3tlsZrfFTv8YuBKYZWZ/XziEktj88cXONzXetHfNbGjsv68ws2lmdsDMdsU9tGRmv5rZE2b2I3DCzMJiH++N3e5WM7sxzvzOzMrEeTzezAbH/jvR5eLMvwRoBAyPfa1Xx8nRJLH3It46ypnZUjM7Env4rdWl1h9v+XAz+0/sNs/FviZnZj8k9P4mlZkViD0M9qeZHTOzWWaWM/a5nmY2x8xGmNlBM9tnZk3jLF4RyBw7bwdgSUqySCpzzulHPwHxA/wKNElg+m9Ar9h/3wlcge8PpXbACaBwYstfbP548xUHTgI5Yx+HAr8DtWKXXQ88C2QCSgE7gZvjbPd7oBiQBSgL7AGuiH2+BFA6zrYcUCbO4/HA4EstFy/vUnyHyRJ8/xJ7L2OfCwe2A0/Fvp7GwHGg7MXWH28drwHfxr7mbMAi4AugVLz5ZgNHEvmZncB6ywBN8ZVGHmAV8FjscyOBw8DNsb+T54BFcZYdBQwFFgMTgcxe/zetn39+tGcjwWAfvg8enHOfO+f2OedinHNTgG1AjcQWTOr8zrndwHdAm9hJjYGTzrlvgepAfufci865s865ncBYoH2cVQx1zu1xzp3Cd74pM1DezMKdc78653Yk4XVe7nLJVQvIDrwa+3qW4CuFu5OysJnlAB4E7ol9zSeAaUCe2Pfmv5xzLZxzuRL5aRF/3c657c65hc65M865w8BCIHfs09fFZp7vfOfwNsdbvAzQ3zl3o3Ouo3PuTJLfEfE7lY0EgyL4/qLFzDqb2fexh3+OABWAfIktmMz5P+WfD9wOsY/Bt9dzxYV1xK7nKXwXM1yw58I/nHPbgf7A88B+M5tsZldc6kVe7nKX4Qpgj/v3RRe78b3PSVEf2Omc2xZnWm7gj5QGM7M7zWyFme2PfZ8HAr/EPl0RmBVn9gr8u3AyuSRcSCLeUNlIQDOz6vg+BP/PzIrj26PoC+R1zuUCNgIWO7uLt+yl5o/vc6ChmRUFbuOfstkD7Ir3V3kO59wtcZb917adc586527AV1QO32GnC04CWeM8LpTE5ZLjYpeZ7gOKmVnc//9fCexN4rrzA39deGBmhu/9mh1/RjP7Kva8T0I/X8WbtzG+19sfXyHmA/YD35tZSSAM2BpnkSr4Dl8C4JxrkMT84gGVjQQkM8tpZi2AycAnzrmf8J0bcMCB2Hm64vvr9oI/8Z1PueBS8/+Lc+4AvnMVH+Irly2xT60BjsWevM9iZqFmViG2CBPKXtbMGptZZuA0cOHQ2gXfAx1i19MMaJDE5ZIj/nsR12p8564ejz3R3xBoie+9ToqNQFUzq2xmWYBX8L3PU+LP6Jxr7pzLnshP83izV8JX7D/j21P6ACiAb+/lOuCneHsuVYAUXZAgaUdlI4Fmlpkdx/eh8zTwFtAVwDm3GXgT30njP/EdVlkRZ9lXgEGxh7oGJGH+hHwKNOGfvRqcc+fxfRhXBnYBB4H3gchE1pEZ32XbB/EdWiqA77DbBQ/Fru8I0BGYnsTlkuNf70XcJ5zvsuBWQPPYbY0EOjvnfk7Kip1z64D/AHPxXShRCLjFOXfuMrNeMBHfxQt/4NtL2gZsjs17HXH2Yswsb+x2N6Zwm5JG9KVOERHxO+3ZiIiI36lsRETE71Q2IiLidyobERHxO5WNiIj4XVqNUBt08uXL50qUKOF1DBGRoLJ+/fqDzrn88aerbBJRokQJ1q1b53UMEZGgYma7E5quw2giIuJ3KhsREfE7lY2IiPidztmIiKSic+fOERUVxenTp72O4lcREREULVqU8PDwJM2vshERSUVRUVHkyJGDEiVK4Lv7QvrjnOPQoUNERUVRsmTJJC2jw2giIqno9OnT5M2bN90WDYCZkTdv3mTtvalsUtuuZXB456XnE5F0Kz0XzQXJfY0qm9QUcx5m9YeRdWDlcN9jERGPPf/88wwZMiTR56dPn87mzZsTfT41qGxSU0godJkFpRrAgqdhXFP407+/QBGRlFLZBKPIInD3ZLhjHPz1K7xXH5a+CtFnvU4mIhnIf/7zH8qWLUuTJk3YunUrAGPHjqV69epUqlSJO+64g5MnT7Jy5UpmzpzJY489RuXKldmxY0eC86WU7tSZiGrVqrkUD1dz4iDMGwg/fQ4FykPr4VDk+tQJKCIBacuWLZQrVw6AF2ZtYvO+Y6m6/vJX5OS5ltdedJ7169dz7733snr1aqKjo6latSo9e/aka9eu5M2bF4BBgwZRsGBB+vXrx7333kuLFi1o27YtAIcOHUpwvvjivtYLzGy9c65a/Hm1Z+NP2fLBHe/D3VPg1BF4vwnMfxrOpvyvBBGRxCxfvpzbbruNrFmzkjNnTlq1agXAxo0bqVevHhUrVmTixIls2rQpweWTOl9y6Hs2aaFsMyheGxY+B6uGw8+zodUwKFnf62Qi4keX2gPxp4SuFrv33nuZPn06lSpVYvz48SxdujTBZZM6X3JozyatRERCy3egy2zAYEJLmPUQnD7qdTIRSWfq16/Pl19+yalTpzh+/DizZs0C4Pjx4xQuXJhz584xceLE/86fI0cOjh8//t/Hic2XEiqbtFayHvRaCXX6wXcfwYiasPUrr1OJSDpStWpV2rVrR+XKlbnjjjuoV68eAC+99BI1a9akadOmXHPNNf+dv3379rzxxhtUqVKFHTt2JDpfSugCgUSkygUCl7J3PczoB/s3QYW20Pw133keEQlaCZ00T690gUCwKHI99FgKDZ+CzTNgeHX48XPQHwAiks6obLwWlgkaPgE9l0OeUvDF/fBpOzga5XUyEZFUo7IJFAXKQbcFcPMr8OtyGFEL1n0AMTFeJxMRSTGVTSAJCYXavX0XEBSpCrMf9l21dmiH18lERFJEZROI8pSEzjN838X54ycYVQdWvAvno71OJiJyWVQ2gcoMqnaGPquh9I2w8FkY1wT+2Oh1MhGRZFPZBLqchaH9RGj7IRzZA2MawJL/QPQZr5OJSIAaOnQo5cqVI3fu3Lz66qtA2ozsfDEqm2BgBhVuh75rfd/HWfa6bzTpPWu9TiYiAWjkyJHMnTuXv/76i4EDBwIqG0mOrHng9veg41Q487fvfjnznoSzJ7xOJiIBomfPnuzcuZNWrVrx9ttv07dv3wRvI5DWNBBnMLqqKfReBYtfgG9Hws9zoOW7ULqR18lEJK6vBvou8klNhSpC81cTfXr06NHMmzePr7/+mtmzZwNQp04dWrVq9a/bCKQ17dkEq4iccOubcO9cCAmDj9vAjD6+WxmIiAQY7dkEuxJ1odcK391AVw6DbYt8JVSuhdfJROQieyAZjfZs0oPwLND0Bei+GLLlhykd4bMu8Pd+r5OJSICIfxuBtKaySU+uqAI9vobGg2DrXBhRA36YrIE9ReR/biOQ1nSLgUSkyS0G/OnAVpjRF6LWQJkm0OIdyFXM61Qi6Z5uMaBbDGQs+cvCffOg+euwexWMrAVrxmpgTxHxhMomPQsJhZoP+C6TLlod5g6A8bfAwW1eJxORDEZlkxHkLg73fAmtR8L+zTCqLix/SwN7ikiaUdlkFGZQpSP0WQNX3+T7Quj7jeH3H71OJpLuZIRz4cl9jSqbjCZHIWj3Cdz1ERz7HcY0hMUvwrnTXicTSRciIiI4dOhQui4c5xyHDh0iIiIiycvoS50ZVfnWUKIezH8alr8JW2ZBq+FwZU2vk4kEtaJFixIVFcWBAwe8juJXERERFC1aNMnz69LnRAT9pc/JsX0RzOoPR6OgRg+48VnInN3rVCIShHTpsySuTBPfFWs1usOaMTCyNmxf7HUqEUlHVDbikzkH3PIGdP0KwjLDJ7fD9N5w8rDXyUQkHVDZyL8Vrw09/w/qPeob6mZETdg8w+tUIhLkVDbyv8IjfOdteiz1Xb32WWeYcg8c/9PrZCISpFQ2krjC10H3JXDjc/DLfN/AnhsmamBPEUk2lY1cXGg41HvEd8+cAuVgRm/f+Zy/dnudTESCiMpGkibfVb67gt4yBPas8V2xtvq9/w7sGROjvR0RSZzKRpIuJMR3eXTvVXBlLfjqcfiwORz4hbcX/cJDkzdw7PQ5r1OKSABS2Ujy5boSOk2DNqPh4FYYXZe6+yYw78c9NH9nOWt26XJpEfk3lY1cHjOofLdvYM+yt1Dr1xFsKPQK5dhJ+zGrGDJ/K+fO6945IuKjspGUyV4A7poA7T4h69lDjD3zOO8XmcPYrzfTdtRKdh084XVCEQkAKhtJHeVaQp/VWOW7aXxwIuvzPU/ug+u5dehypqz9LV2PgCsil6aykdSTJTe0HgH3TCd7WAzjeZZ3cnzCi9PW0POT9fx14qzXCUXEIyobSX2lG0GvVVCrN01PzObbyKc5t3UBzd5dxv9tO+h1OhHxgMpG/CNzdmj2CtZtATlyRPJB2Gu85IbTd9wiXpq9mdPnznudUETSkMpG/KtYDei5HOo/TtPzy/m/7AP5feUk2gz/P37587jX6UQkjahsxP/CMkPjp7EeS8mevzgjMw3liWOD6TpsFuNX7NLFAyIZgMpG0k6hinD/Ymj6Ig1Df2Rh+AA2zR3JvR+sYf/x016nExE/UtlI2goNg7oPYT1XkOXKSrwRPoYeux/lvrensmizbmEgkl6pbMQb+cpgXebArW9RK/MupsY8woqJLzHoi+85dVYXD4ikNyob8U5ICFTvRmjf1YSXrsdz4R9z+/fd6P3ORDbuPep1OhFJRSob8V5kUUI7TYXbx1IxyyHeO/kwi0Y/ypivf+a8bl0gki6obCQwmMF1dxHeby1c05L+oZ9T7+s7eWbkR+w7csrrdCKSQiobCSzZ85Op/Xhc+08pkeU0Lx14iPnv9GDuhp1eJxORFFDZSECya24lS/+1nKzQga7MpNyXzRjx4QSO6+ZsIkFJZSOBK0suctw5kuhOM8idNYw+ux9k8ZCObNi22+tkIpJMKhsJeGFlGpLrkbX8cW03WkYvpOAnDZk+5QOidXM2kaChspHgkCkbhe58i9OdvyIkIidttjzMitdvZ0/Ub14nE5EkUNlIUMlWqhaFHlvD1mv6UOfMMrKNrcO3M97DxWgvRySQqWwk+IRlpmz7lzncaRGHMxWm1obH+enNWzn6p87liAQqlY0ErYJlqlLyiVWsLP0wV/29jtBRtdg+bzhoFGmRgKOykaAWGhZGnXue57e7FrEttBRlvn2aX9+6kbP7d3gdTUTiUNlIulD22kqUffxrphd9jDzHNhMzsjYHF7wJMRrUUyQQqGwk3ciaORNt7h/E9y3ns4YK5Fv5IgffrY/7c5PX0UQyPJWNpDv1q1XimkfmMDLvk9iR3ZwfVZ8TCwZD9Fmvo4lkWCobSZcK5MxCzz5PML/hDL6KqUG2lW/w97C6ELXe62giGZLKRtKtkBCjQ6Prubr3Zzyb7RmOHzlIzPtNiP7qKTh70ut4IhmKykbSvbKFcvBU/4f5uMoUJkU3Imz1CM4Oqwm7lnkdTSTDUNlIhhARHsrjbWpQrPN79Ax9gd+PnYEJLXEzH4TTuiuoiL+pbCRDqX91fl5+pDdvlPyA96JvxX33MeeH1YCtX3kdTSRdU9lIhpMnWyaGdalLZKtXaXf+JXacyAST2sPU++DEQa/jiaRLKhvJkMyM9jWu5LUH72Vg3qG8ea4t0Ztm4IZXhx8/15A3IqlMZSMZWqn82ZnSuwGu/uPccuZltpzJB1/cD5+2g6NRXscTSTdUNpLhhYeGMODmsgzuficPhL/MS9H3cG7HN7gRtWDtONDtC0RSTGUjEqtGyTzM7t+QgxW60fjUK/xEaZjzCExoCYc0sKdISqhsROKIzBLOu+2r8OhdN9PxzJM843pybu8PMKoOrHgXzkd7HVEkKKlsRBLQpkoR5j5Un58Lt6bu36/yY8T1sPBZGNcE/tjodTyRoKOyEUlEsTxZmdyjNp1vqslth/swKHwA5w7/BmMawJL/QPQZryOKBA2VjchFhIYYfRtfxbRedfm/TDdQ4+jLbMzTFJa9DqPrwZ41XkcUCQoqG5EkqFwsF3MerMdN15enRdQ9PJ/jBc6dPg7jboKvBsLZE15HFAloKhuRJMqWOYzX2l7H6E5VmX6iPLWOvswvxdvB6lEwshbs+NrriCIBS2UjkkzNKhRm3kP1KVf8Cm7a2orXC7/NeQuDj9vAjD5w6ojXEUUCjspG5DIUiozgo/tqMOjWcrz/W2EaHB/Mb+UfgO8nwYiasGW21xFFAorKRuQyhYQY99crxfQ+dcmaLTv1v2vAmGvGEZMtP0zpCJ91gb/3ex1TJCCobERSqPwVOZnZ9wburVOClzdkouXpFzlQ43HYOheGV/ft7WhgT8ngVDYiqSAiPJTnW13Lh12r8+eJGOqurMq0GlNw+crC9J4wsS0c2eN1TBHPqGxEUlGjsgWY178e9a/Kx6Nfn+Jee5HjjV6G3at8V6ytGauBPSVDUtmIpLJ82TMztnM1BrepwOpf/6LBsqtZ1nQWFKsBcwfA+Fvg4DavY4qkKZWNiB+YGZ1qFWd2v3oUjoyg8xd/8FS2FzjTYjjs3wKj6sLyt+D8Oa+jiqQJlY2IH5UpkJ0ve9flgQalmLR2D82XFmPL7Yvg6pth8QswtjH8/oPXMUX8TmUj4meZwkJ4snk5Jt5fk1PnztNy/DZGFnyO83dOgON/wJhGsPhFOHfa66gifqOyEUkjdUrnY95D9bn52kK8Pm8rd/9fQfbd8w1Uag/L34TRN8Bv33odU8QvVDYiaSgyazjDO1RhyJ2V2LT3KDeP/omZJQdBpy98tyz4oBnMfRzO/O11VJFUpbIRSWNmRtvrizL3oXpcVSA7D07awCPr83K82zKo+QCsGQMja8P2xV5HFUk1KhsRjxTPm43PHqhN/yZXMeOHfTQftYF15Z6A++ZBeAR8cjtM7w0nD3sdVSTFVDYiHgoLDaF/k6v57IHahJhx13ureOvn3Jzr/g3UGwA/TPYN7Ll5htdRRVLkkmVjZs+lRRCRjOz64rmZ+1A9bq9alKFLtnPn+xv4tdIj0GMp5CgEn3WGKZ18V6+JBKGk7Nk8Z2avmdlYM+tlZrn9nkokA8qeOYwhd1ZiRIeq7DzwN7cMXc5ne3Pjui+BJs/DLwtgRA3YMFEDe0rQSUrZOOA0MB8oBqw0s0p+TSWSgd16XWHm9a/PdUUjeXzqj/Se9CNHqvaBXiugQHmY0Rs+vg3+2u11VJEkM3eJv5DMbJNz7to4j68GRjvnGvs7nJeqVavm1q1b53UMycBiYhxjl+9kyIKt5M2WmbfuqkSdUnlg3ThY9Lxv76bJc1C9O4To9KsEBjNb75yrFn96Uv4LPWhm11944Jz7BcifmuFE5H+FhBgPNCjNl73rkjVzKB3eX83L87Zypup90HsVFK8NXz0OHzaDA1u9jityUUkpmweBT8zsEzN7wswmArv8nEtEYlUoEsmcfvXoVOtKxizbyW0jVrL9bG7oOBVuew8O/uIbfWDZEA3sKQHrkofRAMwsM9AEqAAcAiY55074OZundBhNAtGizX/yxLQf+ftMNP2bXM19N5Qg8+lDvj2cTV9CwYrQejhcUdnrqJJBJXYYLUllkxGpbCRQ7T9+mkFfbmTB5j8pmS8bz7YsT6OyBWDLbJjzCJw4CHX6QcOBEJ7F67iSwaTknI2IBJACOSIY07ka47tWx4CuH67l/glr2V2gEfRZDZU7wIp3fIfWdq/0Oq4IoLIRCVoNyxZgXv/6PNn8GlbtOETTt5cxZNl+TjZ/B+6ZDufPwofNYc6jcOa413Elg1PZiASxTGEhPNCgNEsGNOTWioUZ/vV2mrz5DbNPlMX1WgW1esPacTCiFmxb6HVcycBUNiLpQMGcEbzdrjJTe9YmV9ZM9P10Ax0mbGRr5aeg2wLIlA0mtoUvHtDAnuIJlY1IOlKtRB5m9buBwW0qsOWPY9wydDnPb8jG0S5LoP7jsHEqDK8OG7/QkDeSplQ2IulMaIjRqVZxvn60IXfXKMZHq36l8TurmJLjHmK6L4XIojC1q29gz2O/ex1XMgiVjUg6lTtbJga3qcjMvjdQMl82npj2E7dNO8r3zaZB0xdh+yLf7Qu++0h7OeJ3KhuRdK5CkUg+71mbd9pV5vejp2kzajWP7WvI4c5LoVAFmNkPPmoNhzUwiPiPykYkAzAz2lQpwpIBDXmgQSmmf7+XBuN+Y1yZYZy/5U3Y+x2MqgOrRkLMea/jSjqkshHJQLJnDuPJ5uWY178+VYvn5qU5P9N8xVWsu3UulLgB5j8J426C/Vu8jirpjMpGJAMqnT8747tWZ2znapw6d562k/bQ2z3B4WYj4PBOGF0Pvnkdos96HVXSCZWNSAZlZjQtX5CFDzfg0aZXs2TrAerMycv7103m/DUt4ev/wJiGsHe911ElHVDZiGRwEeGh9LvxKhY/2pAbrynI4G8O0ujXznxfdxTu1GF4vwkseAbOnvQ6qgQxlY2IAFAkVxZGdKzKp/fXJCI8hDaLI+kVOYpj5e6GlUNhdF3YtdzrmBKkVDYi8i91yuRjzoP1eLZFeVbsOcv1P7Tg02uGExMTAxNawKz+cPqo1zElyKhsROR/hIeGcN8NJVkyoCFtKhfhqe/z0ODvl9leugvuuwm+gT1/me91TAkiKhsRSVT+HJl5485KfNm7DnlyRdJk0808lftNTodlh0/vgmn3+27WJnIJKhsRuaQqV+bmy951ee2Oiiw4WoxKfwxiScGuuE3TYUQN+GmqhryRi1LZiEiShIQY7apfyZIBDbm7dhm677mJO90rHAwvDNO6waT2cHSv1zElQKlsRCRZIrOE83yra5nz4A2EFa5AjT8HMjZLN87vWAoja8G6DyEmxuuYEmBUNiJyWa4plJNJ3WsxtMP1fBBzK41OvsK20NIwuz981AoO7fA6ogSQDFE2ZpbNzCaY2Vgz6+h1HpH0wsxocd0VLH60Aa0a1uXWo4/zXEwPzuzZgBtVF1YO08CeAgRx2ZjZB2a238w2xpvezMy2mtl2MxsYO/l2YKpzrjvQKs3DiqRzWTOFMeDmsix8pAF7S99F/ROvsspVgAWDfCMQ/LnJ64jisaAtG2A80CzuBDMLBUYAzYHywN1mVh4oCuyJnU1/Zon4SfG82Xi/S3Ve7XozT0c8Tb+zfTn+xw7ce/Xh65ch+ozXEcUjQVs2zrllwOF4k2sA251zO51zZ4HJQGsgCl/hQBC/ZpFg0ahsAeY9XJ/yN93HzdFDmBldC755jZjR9SFqndfxxAPp7YO3CP/swYCvZIoAXwB3mNkoYFZiC5tZDzNbZ2brDhw44N+kIulc5rBQejUszRePtmJJ+cF0PfsYBw4ewL3fBDfvSTh7wuuIkobSW9lYAtOcc+6Ec66rc66Xc25iYgs758Y456o556rlz5/fjzFFMo5CkRG8274Kvbr3pneuUUyMbox9O5Kzw2vDzm+8jidpJL2VTRRQLM7josA+j7KISBw1SuZhSr+muFvf4j6eZ9/RM/BRK85+2QdOHfE6nvhZeiubtcBVZlbSzDIB7YGZHmcSkVhhoSHcU7sEQx7rw4fXfcJ70S0I/eFTTr5TnZgtc7yOJ34UtGVjZpOAVUBZM4sys27OuWigLzAf2AJ85pzTNZciASZPtky8cEd16vYayZN53mb3qcyETOnAXxM6wd86X5oemdPgeQmqVq2aW7dOV82I+Jtzjunrd/HH3Ne57/znRIdlI+bmV8hRvQNYQqdhJZCZ2XrnXLX404N2z0ZE0gcz47Zqpej0+DAmXPcxv5wrQI65vdkzoiXRh3/zOp6kEpWNiASEHBHh9LjjFnL0XszHuXqR98Aazgytwc6vhmpgz3RAZSMiAaVMoUg6PfQK626dw+aQqym1+hm2vdGAP3fp9GswU9mISMAxM+rXqE7FgUtYWOYZCp7cTuT4BqyYMIjTZzTkTTBS2YhIwIrIFEbTTgM43m0FP2evQd1dw/j1tdqsXrnU62iSTCobEQl4Ra4sReUBc9hywzAKuoNUnX87M9/uza4/4w+PKIFKZSMiwcGMck06k/2R79hVuDmtjk7k/Mgb+PjzzzlxJtrrdHIJKhsRCSrhOfJxdc+JHLl9EnkyRdNxY3dmvdaZ2eu2oe8NBi6VjYgEpVzX3UKeAes5WO4e2sfModLM5gweNpLN+455HU0SoLIRkeCVOQcF2g/jfJe5RObIxjOHn2LjqHt4edoqjpw863U6iUNlIyJBL7RkXXL2X83pmg9xR+hyuv3YnhfeeJ1PV//G+RgdWgsEKhsRSR/CI4ho/iKhPZYQma8Ib7s3yDn7froMm8X63X95nS7DU9mISPpyRWUien+Da/wMzcO/Y8RfvZj43qs8MmUD+4+f9jpdhqWyEZH0JzQcqz+A0F4ryV6kPG9lGk2bTQ/RYcg0xi7bybnzGmstralsRCT9yn81od3mQ/M3uCHzNmaHPMpv89+l+dtLWb5N981JS7qfTSJ0PxuRdOav3TC7P+xYwo8h5eh/qhtXla/CoFvLUyxPVq/TpRu6n42IZGy5i0OnL6DNKCpm+p0FEU9x9bb3afbWYt5e+Aunz533OmG6prIRkYzDDCp3wPqsIeyaZjwaMol52V5g0ZKF3PjmN8zb+LtGIfATlY2IZDw5CkK7j+GujygWdpTZEc/Q133KQ598yz3j1rB9/3GvE6Y7GapszKyUmY0zs6leZxGRAFC+NfRZjVVqz91nPmdNnucJifqWZu8sZ/DszRw/fc7rhOmGX8vGzB4ys41mtsnM+qdgPR+Y2X4z25jAc83MbKuZbTezgRdbj3Nup3Ou2+XmEJF0KGseaDMSOn1BZHgME3iOjwpPZdKKLTQa8g1T10cRo1EIUsxvZWNmFYDuQA2gEtDCzK6KN08BM8sRb1qZBFY3HmiWwDZCgRFAc6A8cLeZlTezimY2O95PgVR5YSKSPpW5EXqvwmo+QJ1DX7Ah77O0zLaZAZ//QNvRK9l35JTXCYOaP/dsygHfOudOOueigW+A2+LN0wCYYWYRAGbWHRgaf0XOuWVAQndJqgFsj91jOQtMBlo7535yzrWI97M/KaHNrKWZjTl69GiSX6iIpBOZs0Pz1+C+eWSKyMpzR59haZkp5HR/kydbJq/TBTV/ls1GoL6Z5TWzrMAtQLG4MzjnPgfmAZPNrCNwH3BXMrZRBNgT53FU7LQExWYZDVQxsycTmsc5N8s51yMyMjIZMUQkXbmyFjywHOo9Som9s/nwZF8its32OlVQ81vZOOe2AK8BC/EVyg/A/9xOzzn3OnAaGAW0cs79nYzNWEKbvkimQ865ns650s65V5KxHRHJaMIj4MZnocdSLEch+KwzTOkEx//wOllQ8usFAs65cc65qs65+vgOg22LP4+Z1QMqAF8CzyVzE1H8e2+pKLDvMuOKiPyvwtdB96+hyfPwywIYUQM2TAR9HydZ/H01WoHY/70SuB2YFO/5KsBYoDXQFchjZoOTsYm1wFVmVtLMMgHtgZmpkV1E5L9Cw+CGh6HXCihwLczoDR/f5hsCR5LE39+zmWZmm4FZQB/nXPybSmQF7nTO7XDOxQBdgP/57ZnZJGAVUNbMosysG0DshQd9gfnAFuAz59wm/70cEcnQ8l0F986BW4ZA1FoYWRu+HQ0xGurmUjQQZyI0EKeIXNSRPb6BPbcvgmI1odUwyF/W61Se00CcIiKpKVcx6DgVbnsPDv4Co2+AZW/AeY06kBCVjYjI5TKDSu2hzxq45lZYMhjGNIJ933udLOCobEREUip7AbhzPLSbCCcOwNjGsPA5OKdRBy5Q2YiIpJZyLaDPaqjcAVa8A6Pqwq8rvE4VEFQ2IiKpKUsuaD0cOs+AmGgYfwvMeRROH/M6madUNiIi/lCqIfReBbV6w9pxvsukty30OpVnVDYiIv6SKRs0ewW6LfQN8jmxLXzRA04c8jpZmlPZiIj4W7Hq8MAyaPAEbJzmG/Jm4xcZasgblY2ISFoIywyNnoIe3/i+ozO1K0zuCMd+9zpZmlDZiIikpUIVoNsiaPoS7FgMI2rCdx+l+70clY2ISFoLDYO6D0KvlVCoIszsBx+1hsO7vE7mNyobERGv5C0NXWZBi3dg73cwqg6sGpEuB/ZU2YiIeCkkBKp19X0ZtEQ9mP8UjLsJ9m/xOlmqUtmIiASCyCLQYQrcMQ7+2gWj68HS1yD6rNfJUoXKRkQkUJhBxba+gT3Lt4alL8OYhrB3vdfJUkxlIyISaLLlg7bj4O7JcOoveL8JLBgEZ096neyyqWxERAJV2ebQ51uo2gVWDoPRdWHXcq9TXRaVjYhIIIuIhJbv+K5acw4mtIBZ/eH0Ua+TJYvKRkQkGJSs7/teTp1+8N0EGFELts7zOlWSqWxERIJFpqxw02DfCARZcsGkdjC1G5w46HWyS1LZiIgEm6LX+8ZYa/gUbJ7hG9jzp6kBPeRNhiobMytlZuPMbKrXWUREUiQsEzR8wjeadO4SMK0bTGoPR/d6nSxBfi0bM3vYzDaZ2UYzm2RmEZe5ng/MbL+ZbUzguWZmttXMtpvZwIutxzm30znX7XIyiIgEpILlfffLufll2PkNjKwF6z6EmBivk/2L38rGzIoADwLVnHMVgFCgfbx5CphZjnjTyiSwuvFAswS2EQqMAJoD5YG7zay8mVU0s9nxfgqkygsTEQk0IaFQu4/vzqBXVIbZ/eGjVnBoh9fJ/svfh9HCgCxmFgZkBfbFe74BMOPCHo+ZdQeGxl+Jc24ZcDiB9dcAtsfusZwFJgOtnXM/OedaxPvZn4qvS0Qk8OQpCZ1nQsuh8PsPvoE9VwyF89FeJ/Nf2Tjn9gJDgN+A34GjzrkF8eb5HJgHTDazjsB9wF3J2EwRYE+cx1Gx0xJkZnnNbDRQxcyeTGSelmY25ujR4LqGXUQE8A15c30X38CepRvDwmdgXFP4c5Onsfx5GC030BooCVwBZDOzTvHnc869DpwGRgGtnHN/J2czCUxL9HIM59wh51xP51xp59wricwzyznXIzIyMhkxREQCTM4roP2n0PYDOPIbvFcfvn4Zos94Esefh9GaALuccwecc+eAL4A68Wcys3pABeBL4LlkbiMKKBbncVH+91CdiEjGZAYV7vAN7FnhDvjmNXivAUStS/Mo/iyb34BaZpbVzAy4EfjXDRrMrAowFt8eUFcgj5kNTsY21gJXmVlJM8uE7wKEmamSXkQkvciWF24fAx0+hzPHfAN7znsKzp5Iswj+PGezGpgKfAf8FLutMfFmywrc6Zzb4ZyLAboAu+Ovy8wmAauAsmYWZWbdYrcRDfQF5uMrss+cc94emBQRCVRX3wS9v4Vq98G3I3wXEOz8Jk02bS6Av3HqpWrVqrl169J+V1NEJE38ugJm9oPDO6BqZ2j6km8InBQys/XOuWrxp2eoEQRERCRWibrQawXUfQg2fAIjasLPc/y2OZWNiEhGFZ4Fmr4I9y/23bBtcgf4vCucTOhrjSmjshERyeiKVIUeS6HRIN8tqC31qyEs1dcoIiLBJzQcGjwGdR+EsMypvnrt2YiIyD/8UDSgshERkTSgshEREb9T2YiIiN+pbERExO9UNiIi4ncqGxER8TuVjYiI+J0G4kyEmR0AjgCXc8vOfMDBVA0kFxPJ5f2eAlmgviavcvl7u6m9/tRaX0rWc7nLpvTzq7hzLn/8iSqbizCzMc65Hpex3LqERj0V/7jc31MgC9TX5FUuf283tdefWutLyXoC7fNLh9EubpbXASRJ0uPvKVBfk1e5/L3d1F5/aq0vJesJqP+GtGfjB9qzEZFgpT2b4BL/jqQiIsHCL59f2rMRERG/056NiIj4ncpGRET8TmUjIiJ+p7JJA2ZWyszGmdlUr7OIiCSHmbUxs7FmNsPMbrrc9ahsLpOZfWBm+81sY7zpzcxsq5ltN7OBAM65nc65bt4kFRH5t2R+fk13znUH7gXaXe42VTaXbzzQLO4EMwsFRgDNgfLA3WZWPu2jiYhc1HiS//k1KPb5y6KyuUzOuWXA4XiTawDbY/dkzgKTgdZpHk5E5CKS8/llPq8BXznnvrvcbapsUlcRYE+cx1FAETPLa2ajgSpm9qQ30URELirBzy+gH9AEaGtmPS935WEpyybxWALTnHPuEHDZvyQRkTSQ2OfXUGBoSleuPZvUFQUUi/O4KLDPoywiIsnh188vlU3qWgtcZWYlzSwT0B6Y6XEmEZGk8Ovnl8rmMpnZJGAVUNbMosysm3MuGugLzAe2AJ855zZ5mVNEJD4vPr80EKeIiPid9mxERMTvVDYiIuJ3KhsREfE7lY2IiPidykZERPxOZSMiIn6nshEREb9T2YiIiN+pbESChJl9aWaDzWy5mf1hZk28ziSSVCobkeBRATjinKsH9AY6epxHJMlUNiJBwMyyApHA27GTwoAjngUSSSaVjUhwuBZY75w7H/v4OmDjReYXCSgqG5HgUAH4Ps7j64AfvYkiknwqG5HgUJF/l00FtGcjQUS3GBAREb/Tno2IiPidykZERPxOZSMiIn6nshEREb9T2YiIiN+pbERExO9UNiIi4ncqGxER8bv/B1sxP6O9OtSXAAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "from matplotlib import pyplot as plt\n", - "\n", - "data = np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]])\n", - "a_fit, p_fit = fit_power_law(data)\n", - "\n", - "plt.figure()\n", - "plt.title(\"Data versus fit of $\\sigma = an^p$\")\n", - "plt.loglog(*data.T, label=\"data\")\n", - "n_range = np.logspace(1, 2, 10)\n", - "plt.loglog(n_range, a_fit*n_range**p_fit, label=\"fit\")\n", - "plt.xlabel(\"$n$\")\n", - "plt.ylabel(\"$\\sigma$\")\n", - "plt.legend()\n", - "plt.show()" - ] - }, - { - "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", - "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 -}