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

01: Plot dim d volumes, calculate variance of Metropolis

Browse Source
This commit is contained in:
Kees van Kempen
2022-09-13 10:43:48 +02:00
parent bc85efdb6f
commit e507977f46
1 changed files with 32 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
Exercise sheet 1/exercise_sheet_01.ipynb
View File

@ -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)))"
]
},
{