diff --git a/Exercise sheet 1/exercise_sheet_01.ipynb b/Exercise sheet 1/exercise_sheet_01.ipynb index d98c688..f3f0c9d 100644 --- a/Exercise sheet 1/exercise_sheet_01.ipynb +++ b/Exercise sheet 1/exercise_sheet_01.ipynb @@ -537,13 +537,17 @@ " return np.pi**(d/2)/scipy.special.gamma(d/2 + 1)\n", "\n", "N = 10000\n", - "for d in range(1, 8):\n", - " print(\"N = {}\".format(N),\n", - " \"\\td = {}\".format(d),\n", - " \"\\tV_d = {:.3f}\".format(volume_d_ball(d)),\n", - " \"\\tdirect-sampled V_d = {:5.3f}\".format(number_of_hits_in_d_dimensions(N, d)/N*2**d)\n", - " )\n", - "# TODO: Plot the results as is asked." + "d_range = np.array(range(1, 8))\n", + "V_d_mc = [number_of_hits_in_d_dimensions(N, d)/N*2**d for d in d_range]\n", + "\n", + "## Plotting\n", + "plt.figure()\n", + "plt.plot(d_range, V_d_mc, label=\"direct-sampling\")\n", + "plt.plot(d_range, volume_d_ball(d_range), label=\"exact\")\n", + "plt.xlabel(\"dimension $d$\")\n", + "plt.ylabel(\"volume $V_d$\")\n", + "plt.legend()\n", + "plt.show()" ] }, { @@ -572,7 +576,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 11, "metadata": { "deletable": false, "nbgrader": { @@ -592,37 +596,37 @@ "name": "stdout", "output_type": "stream", "text": [ - "N = 2000 \tdelta = 0.000 \tpi_estimate = 4.000\n", - "N = 2000 \tdelta = 0.333 \tpi_estimate = 3.149\n", - "N = 2000 \tdelta = 0.667 \tpi_estimate = 3.139\n", - "N = 2000 \tdelta = 1.000 \tpi_estimate = 3.135\n", - "N = 2000 \tdelta = 1.333 \tpi_estimate = 3.146\n", - "N = 2000 \tdelta = 1.667 \tpi_estimate = 3.150\n", - "N = 2000 \tdelta = 2.000 \tpi_estimate = 3.150\n", - "N = 2000 \tdelta = 2.333 \tpi_estimate = 3.141\n", - "N = 2000 \tdelta = 2.667 \tpi_estimate = 3.146\n", - "N = 2000 \tdelta = 3.000 \tpi_estimate = 3.146\n" + "N delta 4*hits/trials E[(4*hits/trials - pi)^2]\n", + "2000 0.0000 4.0000 0.737\n", + "2000 0.3333 3.1409 0.014\n", + "2000 0.6667 3.1351 0.005\n", + "2000 1.0000 3.1381 0.006\n", + "2000 1.3333 3.1406 0.004\n", + "2000 1.6667 3.1316 0.007\n", + "2000 2.0000 3.1377 0.007\n", + "2000 2.3333 3.1311 0.011\n", + "2000 2.6667 3.1547 0.011\n", + "2000 3.0000 3.1417 0.015\n" ] } ], "source": [ "N = 2000 # number of trials\n", "m = 200 # number of repetitions\n", + "delta_range = np.linspace(0, 3, 10) # throwing ranges\n", "\n", "p0 = np.array([0., 0.]) # starting position in origin\n", + "hits = np.zeros((len(delta_range), m), dtype=int)\n", "\n", - "for delta in np.linspace(0, 3, 10):\n", - " pi_estimate = 0\n", + "for index_delta, delta in enumerate(delta_range):\n", " for i in range(m):\n", - " hits = markov_pebble(p0, delta, N)\n", - " pi_estimate += hits/N*4\n", - " pi_estimate /= m\n", - " print(\"N = {}\".format(N),\n", - " \"\\tdelta = {:.3f}\".format(delta),\n", - " \"\\tpi_estimate = {:.3f}\".format(pi_estimate)\n", - " )\n", + " hits[index_delta][i] = markov_pebble(p0, delta, N)\n", "\n", - "# TODO: Investigate the rejection rate and the mean square deviation. Maybe only the latter, as the next question concerns the former." + "print(\"{:4} {:6} {:8} {:<10}\".format('N','delta','4*hits/trials','E[(4*hits/trials - pi)^2]'))\n", + "for index_delta, delta in enumerate(delta_range):\n", + " pi_estimates = hits[index_delta]/N*4\n", + " print(\"{:4} {:1.4f} {:1.4f} {:1.3f}\"\n", + " .format(N, delta, np.mean(pi_estimates), np.mean((pi_estimates - np.pi)**2)))" ] }, {