1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 02.1\n",
"- Compute the following 1D integral via Monte Carlo \n",
"$$I = \\int_0^1 \\frac{\\pi}{2}\\cos(\\pi x/2) dx = 1$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. sampling a uniform distribution in $[0,1]$\n",
"2. using importance sampling (i.e. sampling a non-uniform probability in $[0,1]$)\n",
"\n",
"<span style=\"color:red\">Show a picture of both your estimations of $I$</span> and their uncertainties with a large number of *throws* $M$ (e.g. $M\\ge 10^4$) as a function of the number of blocks, $N$</font>."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The taylor Expansion of the cosine is:\n",
"\n",
"$$\\cos(x) = 1 - \\frac{x^2}{2} + O(x^3).$$\n",
"\n",
"So I choose to sample the integral with distribution:\n",
"$$d(x) = 1 - x^2.$$\n",
"This must be normalized so that $\\int_0^1 d(x) \\,dx = 1$:\n",
"$$\\int_0^1 \\left(1 - x^2\\right) \\,dx = \\left(x - \\frac{x^3}{3}\\right)\\bigg|_0^1 = \\frac{2}{3} \\quad \\implies \\quad d(x) = \\frac{3}{2}\\left(1 - x^2\\right).$$\n",
"Inverting the cumulative distribution requires solving a cubic equation, so I resort to accept-reject to sample $d(x)$.\n",
"$$y = F(x) = \\int_0^x dx' d(x') \\rightarrow x = F^{-1}(y)$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f16e1d37f40>]"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAAsTAAALEwEAmpwYAAAka0lEQVR4nO3deXRV9b3+8ffnZCQkEIaEOYAQwMhMAEVUnAEHtAUFVERtFRSttbba3lvbXvvrcLXWolJEBGepCiq1Vto64MCYMI8SQCAyzyQhCUm+vz+S5qYYJEDO3uecPK+1shb77J3k2SvhPNnT92vOOURERAACfgcQEZHQoVIQEZFKKgUREamkUhARkUoqBRERqRTtd4Az0bRpU9euXTu/Y4iIhJXs7Oy9zrmU6taFdSm0a9eOrKwsv2OIiIQVM9tyonU6fSQiIpVUCiIiUkmlICIilVQKIiJSSaUgIiKVPCkFM5tmZrvNbNUJ1puZTTSzHDNbYWa9vcglIiL/yasjhReAwd+yfgiQXvFxJ/BnDzKJiMhxPHlOwTn3qZm1+5ZNhgEvufJxvBeYWbKZtXDO7QhGnvU7j/C3FdsxMwJmBAyiooyYQIDoKCMuOoq46ADxMVEkxP77I5rE+GgaxEeTFB9DbLTOvIlI5AmVh9daAduqLOdWvPaNUjCzOyk/miAtLe20vlnO7jye+jiHM5lKol5MFI0SYmiYEEuT+rE0SYylcf1YUpLiSE2KJzUpjuYN42neMJ4G8TGn/41ERDwUKqVg1bxW7Vu2c24KMAUgMzPztN7Wr+regqu6X4VzDueg1DlKyxzHSssoKXUUlZRRVFJK4bEyCopLKCguJb+ohLyiEo4UlnD46DEOHT3GgYJjHCwoZn9BMVu3FrAvr4j84tJvfL/EuGhaJsfTulECrZLr0aZxPdIa16dtkwTaNkkgITZUfgwiUteFyrtRLtCmynJrYHuwv6mZYQYBjJgoiI+JOuOvWVBcwu7DRew+UsTOw4XsPHSU7QcL2X7wKLkHjpL11X4OF5b8x+e0aBhP+6b16ZCSSHqzRDqmJpKemkRKUtwZ5xERORWhUgqzgQlmNgPoDxwK1vWEYEuIjaZd02jaNa1/wm0OFRxjy/58tuwr4Ku9+Wzem8/Gvfm8s/RrjhT9X2E0qR9LlxZJdG7WgIyWDTinZQM6piYSE6XrGSISHJ6Ugpm9DgwCmppZLvALIAbAOTcZeB8YCuQABcBtXuTyS8OEGLonJNO9dfJ/vO6cY/eRIjbsyuPLXUdYv/MI63Yd4bVFWyg8VgZAbFSAs1sk0a11Q7q3SqZHm2Q6piYSFajuDJyIyKkxdyZXW32WmZnp6sIoqaVljs1781m9/RBrth9mRe4hVn19qPKoIjEumu6tG9I7rRF92jWid1ojGtbTxW0RqZ6ZZTvnMqtdp1IIT2Vljs378lm29SDLth1kydYDrNt5hNIyhxl0bpZE//aN6de+Cf3PakzTRF2fEJFyKoU6Ir+ohOXbDpK15QCLNu8ne8sBjh4rvxuqc7MkzuvQhAEdmnBuhya6TVakDlMp1FHHSstY+fUhFmzax/yN+1j81X4Kj5URFTB6tG7IhZ1SuKhTCt1bJ+uahEgdolIQAIpLyliy9QCfb9jLZzl7WZF7EOegUUIMF6SncOnZqVzUKYXkhFi/o4pIEKkUpFoH8ov5LGcvn6zfzdz1e9iXX0zAILNtYy7LSOXyjOa0/5Zba0UkPKkU5KTKyhzLcw/y0brd/GvtbtbuOAxAemoig7s258pzmnNOywaY6TSTSLhTKcgp27a/gH+t3cU/Vu9i4eZ9lDlo07geQ7u14OpuLenaSgUhEq5UCnJG9uUV8a+1u3h/5U6+yNlLSZkjrXEC1/RowbU9WtG5eZLfEUXkFKgUpNYcLCjmH6t38dcV25m3cR+lZY7OzZK4rlcrhvVsScvken5HFJGTUClIUOzNK+L9lTt4d9l2srccwAzObd+E7/RuxdBuLagfFypDa4lIVSoFCbot+/J5Z+l2Zi3NZcu+AhJioxjStQUjMlvTv31jXX8QCSEqBfGMc47sLQd4KzuX91bsIK+ohPZN6zMiszXD+7QmNSne74gidZ5KQXxxtLiU91fu4C9Z21i0eT/RAePyjGaM6pfGwI5NCegpahFfqBTEd5v25DFj8Tbeys5lf34xaY0TuKl/GiMy29C4vp6gFvGSSkFCRlFJKR+s2smrC7ay6Kv9xEYHuKZ7S8YOaEe31g39jidSJ6gUJCSt33mElxd8xawlX1NQXEqfto0YO6AdQ7o2J1qzy4kEjUpBQtrhwmO8mZXLS/O/Ysu+Alo2jOfWAe0Y2TeNhgka4luktqkUJCyUlTk+Wreb5z/fzPxN+0iIjeLGvm24/fz2tGmc4Hc8kYihUpCws2b7YaZ+tonZy7dT5hxDu7Vg3EUd6NpK1x1EzpRKQcLWjkNHeeGLr3h14Vbyikq4IL0p4y/qwHkdmuiBOJHTpFKQsHe48BivLtjKtC82s+dIET3bJDPh4o5cenaqykHkFKkUJGIUHitl5pJcJs/dyLb9R+nSPIn7Lk1n8DnN9TCcSA2pFCTilJSWMXv5dp7+OIdNe/JJT03k3kvTuapbC803LXISKgWJWKVljr+t3MFTH25gw+480lMT+cFl6Qzt2kJHDiIn8G2loCeEJKxFBYxre7Rkzv0X8vToXgBMeG0pQyd+xpzVOwnnP3pE/KBSkIgQCBhXd2/JB/dfyJ9G9qS4pIy7Xs5m2DNf8Mn63SoHkRpSKUhEiQoYw3q24h8/vJDHhndnf34xY6cvZuSUBWRvOeB3PJGQp2sKEtGKS8p4fdFWnvpoA3vzirk8oxkPDe5Mx1TNKy11ly40S52XX1TC9C828+zcTeQXl3Bj3zbcf1knmjXQpD9S96gURCrszy/mqY828MqCLUQHAnz/wrO468KzNJ+01Cm6+0ikQuP6sfzimnP48IFBXHp2KhM/3MCgxz9hxqKtlJaF7x9IIrXFs1Iws8Fmtt7Mcszs4WrWNzSzv5rZcjNbbWa3eZVN6p60Jgk8Pbo3s+4eQJtG9Xh41kqumvgZ8zbu9TuaiK88KQUziwKeAYYAGcAoM8s4brN7gDXOuR7AIOAPZqZ5GiWoeqc1Yub4ATw9uhdHCksY/dxC7nwpiy378v2OJuILr44U+gE5zrlNzrliYAYw7LhtHJBk5aObJQL7gRKP8kkdZlb+jMOHP7qIH1/Zmc9z9nL5E5/yvx+sI79Iv4JSt3hVCq2AbVWWcyteq+pp4GxgO7AS+IFzruz4L2Rmd5pZlpll7dmzJ1h5pQ6Kj4ninos78vGDg7i6ewsmfbKRix//hHeWfq2H36TO8KoUqhuE5vj/ZVcCy4CWQE/gaTNr8I1Pcm6Kcy7TOZeZkpJS2zlFaNYgnidu7MnM8QNo3jCe+/+yjBunLGDdzsN+RxMJOq9KIRdoU2W5NeVHBFXdBsxy5XKAzUAXj/KJfEOfto14++7z+c313fhy1xGumvg5//PXNRwpPOZ3NJGg8aoUFgPpZta+4uLxSGD2cdtsBS4FMLNmQGdgk0f5RKoVFTBG90/j4x8N4sa+bZg+bzOXPTGX91Zs1ykliUielIJzrgSYAMwB1gJvOOdWm9k4MxtXsdmjwAAzWwl8CDzknNP9gRISGtWP5TfXd+Ptu88nJSmOCa8tZcy0RbpLSSKOnmgWOUWlZY5XFmzhsTnrOVZaxn2XpvP9C84iNlrPgkp40BPNIrUoKmDcOqAd/3rgIi49O5XH5qzn6qc+0yisEhFUCiKnqXnDeCbd1IdpYzPJKyxh+OR5PPLuKl2IlrCmUhA5Q5d0acY/HriIW89rx8sLtnDFHz/lo3W7/I4lclpUCiK1IDEuml9eew6zxg8gKT6a21/I4v4ZS9mfX+x3NJFTolIQqUW90hrx3r0XcP9l6fxt5Q4uf2Iuf1uxw+9YIjWmUhCpZbHRAe6/rBPv3XsBrRrV457XlnD3q9nszSvyO5rISakURIKkc/MkZo0fwE8Gd+Zfa3brqEHCgkpBJIiiowLcPagjf7tvIGmNE7jntSVMeG2JrjVIyFIpiHggvVkSM8cP4MdXdmbO6p1c8cdP+XCt7lCS0KNSEPFIdFSAey7uyOwJA2maGMsdL2bx8MwV5GnOBgkhKgURj53dogHvTjifuwd14I2sbQx+8lMWf7Xf71gigEpBxBdx0VH8ZHAX3hx3HgEzbnx2Po9XjKUk4ieVgoiP+rRtzPs/uIDhfVrz9Mc5fPfP89i4J8/vWFKHqRREfJYYF83/Du/B5Jt7s3V/AVdP/JwZi7ZqvgbxhUpBJEQM7tqCD35wIb3Sknl41krGv7KEgwW6dVW8pVIQCSHNG8bzyh39+emQLny4bhdD/vQZCzft8zuW1CEqBZEQEwgYd13UgZnjBxAXHWDUcwv44z+/pEQXocUDKgWRENW9dTLv3XcB1/VsxZ8+3MDoqQvZeajQ71gS4VQKIiEsMS6aJ27syR9G9GDV14cYOvEzPl6/2+9YEsFUCiJh4Lt9WjN7wkBSk+K4bfpifvv3tXqmQYJCpSASJjqmJvLOPeczun8az87dxKgpC9hx6KjfsSTCqBREwkh8TBS/ub4bfxrZkzU7DnPVxM+Z++Uev2NJBFEpiIShYT1b8dd7B5KSGMfY6Yv44z+/pLRMD7vJmVMpiISpDinlp5Ou71V+d9LY6Ys0T4OcMZWCSBirFxvFH0b04Lff6cbCzfu5euJnLNt20O9YEsZUCiJhzswY1S+NWeMHYGbcMHk+ry3U2ElyelQKIhGia6uGvHfvQM7t0ISfvb2Sn7y1gsJjpX7HkjCjUhCJII3qxzJ9bF/uu6Qjb2bncsOz89l+ULetSs2pFEQiTFTAeOCKzky5pQ+b9uRzzVOfM3+jBtWTmlEpiESoK85pzrsTzqdR/Vhufn4hL3yxWdcZ5KQ8KwUzG2xm680sx8wePsE2g8xsmZmtNrO5XmUTiVQdUhJ5++4BXNw5lV/+dQ0/1nUGOQlPSsHMooBngCFABjDKzDKO2yYZmARc65w7BxjhRTaRSJcUH8OUW/pw36XpvJWdy41TFrDrsEZblep5daTQD8hxzm1yzhUDM4Bhx20zGpjlnNsK4JzTUJAitSQQMB64vBOTb+7Nhl1HuPbpz1mu5xmkGl6VQitgW5Xl3IrXquoENDKzT8ws28zGVPeFzOxOM8sys6w9ezTmi8ipGNy1BTPHDyA6EGDEs/N5Z+nXfkeSEONVKVg1rx1/xSsa6ANcBVwJ/NzMOn3jk5yb4pzLdM5lpqSk1H5SkQh3dosGzJ5wPr3aJHP/X5bx2Jx1lGncJKngVSnkAm2qLLcGtlezzQfOuXzn3F7gU6CHR/lE6pQmiXG8fEd/RvVrwzMfb2T8q9nkF5X4HUtCgFelsBhIN7P2ZhYLjARmH7fNu8AFZhZtZglAf2CtR/lE6pzY6AC/ub4bj1ydwT/X7GLEZD3oJh6VgnOuBJgAzKH8jf4N59xqMxtnZuMqtlkLfACsABYBU51zq7zIJ1JXmRm3D2zPtLF92bq/gOue+YIVuQf9jiU+snB+mCUzM9NlZWX5HUMkIqzfeYQ7XlzM3rwi/nhDT4Z0a+F3JAkSM8t2zmVWt05PNIsIAJ2bJ/HOPeeT0aIB419dwrNzN+oJ6DpIpSAilZomxvHa98/l6u4t+O3f1/Gzt1dxrLTM71jioWi/A4hIaImPiWLiyF60bZLAMx9vJPdAAZNu6k1SfIzf0cQDOlIQkW8IBIwfX9mF/x3enfkb9zFi8nx2HNKdSXWBSkFETuiGzDa8cFs/cg8c5bpnvmDN9sN+R5IgUymIyLcamN6Ut8afR8CMEZPn8dkGDS8TyU65FMysfsWopyJSR3Rp3oC37z6fNo0TuG36YmZm5/odSYLkpKVgZgEzG21mfzOz3cA6YEfFnAePmVl68GOKiN+aN4znjXHn0f+sxvzozeU883GOblmNQDU5UvgY6AD8FGjunGvjnEsFLgAWAL8zs5uDmFFEQkSD+Bimj+3HdT1b8tic9fz83VWUajC9iFKTW1Ivc84dO/5F59x+YCYw08x0r5pIHREbHeCJG3rSrGE8z87dxN4jxTw5sifxMTqrHAlOeqTw70IwsyfNrLohsKmuNEQkcgUCxk+HnM3Pr87gg9U7GfP8Ig4V6G0gEpzKheY8YLaZ1QcwsyvM7IvgxBKRcHDHwPZMHNWLpdsOcMOz8zXNZwSocSk45/4beB34xMw+B34EPBysYCISHq7t0ZLpY/uRe6CA70yax6Y9eX5HkjNQ41Iws0uB7wP5QApwn3Pus2AFE5HwMTC9Ka/feS6Fx0oZPnk+K3MP+R1JTtOpnD76L+DnzrlBwHDgL2Z2SVBSiUjY6d46mTfHnUe9mChGTpnPvI17/Y4kp+FUTh9d4pz7vOLfK4EhwK+DFUxEws9ZKYnMHD+Alsn1GDt9MXNW7/Q7kpyimjy8dqI7jnYAl37bNiJS9zRvGM8bd51XPi/DK9m8paefw0qNHl4zs3vNLK3qixVzLZ9nZi8CtwYlnYiEpUb1Y3n1e/0Z0KEpD765nOlfbPY7ktRQTUphMFAKvG5mO8xsjZltBjYAo4A/OudeCGJGEQlD9eOimXprJldkNONXf13DxA83aFiMMFCTJ5ovBt5xzk2qeHK5KXDUOXcwqMlEJOzFx0Qx6abe/GTmCp7455ccKTzGz4aejc44h66alMJ3gUfNrBnlg+EtA5ab2VJgnXOuNIj5RCTMRUcFeHx4D5Lionnus83kF5fy62FdCQRUDKHopKXgnPsegJndD6QDmyk/epgC7AdaBzGfiESAQMD45bXnkBAXzZ8/2cjR4lIeG96d6ChN6RJqTmWO5tuccz3+vWBmk4Af134kEYlEZsZDg7uQGBfNY3PWU1BcwlOjehMbrWIIJafy0zhsZn3+veCcywY61X4kEYlk91zckZ9fncGc1bu46+UsCo/pDHQoOZUjhduBV8xsDZANdAM0LKKInLI7BrYnPibAf729ijteXMxzYzJJiD2VtyMJllN5onkDMAB4H2gGrAWGBimXiES4m/q35fERPZi/cR9jpy0mr6jE70jCqR0pUHGn0ZsVHyIiZ2R4n9bERQe4/y/LGPP8Ql64vR8N4jVnl590hUdEfHVNj5Y8PaoXK3IPccvUhZqsx2cqBRHx3ZBuLfjzzX1Ys+MwNz2/gIMFxX5HqrNUCiISEi7PaMaUWzL5clceo59byIF8FYMfVAoiEjIu7pLKlFv6kLMnj9FTF7JfxeA5z0rBzAab2XozyzGzE07jaWZ9zazUzIZ7lU1EQsegzqlMHZPJpj15jH5ugYrBY56UgplFAc9QPjFPBjDKzDJOsN3vgTle5BKR0HRhpxSev7Uvm/fmqxg85tWRQj8gxzm3yTlXDMwAhlWz3b3ATGC3R7lEJEQNTG/K1Fsz2bw3n5um6hqDV7wqhVbAtirLuRWvVTKzVsD1wORv+0JmdqeZZZlZ1p49e2o9qIiEjgvSU3huTCYbK64x6K6k4POqFKobI/f42TaeBB462VDczrkpzrlM51xmSkpKbeUTkRB1YacUplYUw83P6zmGYPOqFHKBNlWWWwPbj9smE5hhZl8Bw4FJZnadJ+lEJKRd2CmFZ2/uw5c78xgzbSGHC1UMweJVKSwG0s2sfcXcziOB2VU3cM61d861c861A94C7nbOveNRPhEJcRd3SWXSTb1Zs+Mwt05bpLGSgsSTUnDOlQATKL+raC3whnNutZmNM7NxXmQQkfB3WUYznh7dm5W5h7ht+iIKilUMtc3CeSLtzMxMl5WV5XcMEfHYeyu2c9/rSzn3rCZMG9uX+JgovyOFFTPLds5lVrdOTzSLSNi5untL/nBDD+Zv2sddL2dTVKKJemqLSkFEwtL1vVrzu+90Y+6Xe7j3taUcKy3zO1JEUCmISNi6sW8av7r2HP6xZhcPvLGc0rLwPR0eKjT/nYiEtVsHtKOguJTff7COejEBfved7gQC1T0aJTWhUhCRsDd+UAeOHitl4ocbSIiN5hfXZGCmYjgdKgURiQg/vCyd/KISnv98M4lx0Tx4ZWe/I4UllYKIRAQz47+vOpuC4hKe/jiH+nHRjB/Uwe9YYUelICIRw8z49XXdKq8xJMZHc8u5bf2OFVZUCiISUaICxuMjepBfVMIj766iQXw0w3q2OvknCqBbUkUkAsVEBXh6dG/6tWvMj95YzkfrdvkdKWyoFEQkIsXHRDH11kwyWjZg/CtLWLhpn9+RwoJKQUQiVlJ8DC/c1o82jRP43otZrPr6kN+RQp5KQUQiWuP6sbx8Rz8a1Ith7PRFbN6b73ekkKZSEJGI16JhPV6+ox/Owc1TF7LzUKHfkUKWSkFE6oSzUhJ58fZ+HDp6jDHTNN/ziagURKTO6NqqIVPG9OGrvQXc8WIWR4s15PbxVAoiUqcM6NCUJ0f2ZMnWA9zz2hINuX0clYKI1DlDu7Xg0WFd+Wjdbh6auYJwnoGytumJZhGpk24+ty1784p48l8bSE2K5+EhXfyOFBJUCiJSZ/3g0nT2HCli8tyNpCbFcfvA9n5H8p1KQUTqLDPjf4Z1ZV9eMf/z3hqaJsVxbY+Wfsfyla4piEidFhUwnhzZk37tGvPgG8uZt3Gv35F8pVIQkTovPiaK58Zk0rZJAne9lM26nYf9juQblYKICNAwIYYXb+9H/bhobp22iK8PHvU7ki9UCiIiFVom1+OF2/tSUFzK2GmLOFRwzO9InlMpiIhU0aV5A569pQ9f7cvnzpezKCqpW089qxRERI4zoENTHh/Rg4Wb9/PgmysoK6s7D7fpllQRkWoM69mK7QcL+f0H62iZHM9Ph5ztdyRPqBRERE5g3EVn8fXBAp6du4k2jRK4+dy2fkcKOpWCiMgJmBm/vOYcth8s5JF3V9EyOZ5LujTzO1ZQeXZNwcwGm9l6M8sxs4erWX+Tma2o+JhnZj28yiYiciLRUQGeGtWLjJYNmPDaUlbmRvaUnp6UgplFAc8AQ4AMYJSZZRy32WbgIudcd+BRYIoX2URETqZ+XDTTbu1Lo4RYbn9xcUQ/w+DVkUI/IMc5t8k5VwzMAIZV3cA5N885d6BicQHQ2qNsIiInldognum39aWwuJTbpy/mSGFkPsPgVSm0ArZVWc6teO1E7gD+Xt0KM7vTzLLMLGvPnj21GFFE5Nt1apbEn2/uw8Y9edz9amRO0ONVKVg1r1V746+ZXUx5KTxU3Xrn3BTnXKZzLjMlJaUWI4qInNzA9Kb8v+u78tmGvTzy7qqIm6DHq7uPcoE2VZZbA9uP38jMugNTgSHOuX0eZRMROSU39k1jy74CJn2ykbOaJvL9C8/yO1Kt8epIYTGQbmbtzSwWGAnMrrqBmaUBs4BbnHNfepRLROS0PHhFZ4Z2a85v/r6WOat3+h2n1nhSCs65EmACMAdYC7zhnFttZuPMbFzFZo8ATYBJZrbMzLK8yCYicjoCAeOJG3rSvXUy989YxqqvI+NWVQvn82GZmZkuK0vdISL+2X2kkOufmUdJWRnv3jOQ5g3j/Y50UmaW7ZzLrG6dBsQTETkDqUnxPD82k7zCEr730mIKikv8jnRGVAoiImeoS/MGPDW6F2u2H+aHf1kW1qOqqhRERGrBJV2a8V9XZTBn9S4e+8d6v+OcNg2IJyJSS24/vx0b9+Tx50820jElke/2Cb+BGXSkICJSS8yMX117DgM6NOGns1aSvWW/35FOmUpBRKQWxUQFmHRTb1omx3PXy9nkHijwO9IpUSmIiNSy5IRYpt7al6KSMr73Yhb5ReFzR5JKQUQkCDqmJvL06N58uesID7wRPnckqRRERILkok4plXckPfmv8Bi9R3cfiYgE0e3nt2P9zsNM/CiHTs2TuLp7S78jfSsdKYiIBJGZ8eh1Xcls24gH31we8mMkqRRERIIsLjqKybf0oXFCLHe+lMWeI0V+RzohlYKIiAeaJsYxZUwm+wuKGf9KNsUloTlrm0pBRMQjXVs15LHhPcjacoBfzA7NWdt0oVlExEPX9GjJ2h2HmfTJRjJaNuSWc9v6Hek/6EhBRMRjP7qiM5d0SeVXs1ezaHNoDYWhUhAR8VhUwHhyZE/SGidw96vZbD941O9IlVQKIiI+aBAfw5QxmRQeK+POl7MoPFbqdyRApSAi4puOqYk8eWNPVn19mJ+9vTIkLjyrFEREfHRZRjN+eFknZi35mhfnfeV3HJWCiIjf7r2kI5ed3YxH/7aWBZv2+ZpFpSAi4rNAwPjjjT1o2ySBe15d4uuFZ5WCiEgISIqPYcotmRSVlDH+lWzfLjyrFEREQkTH1EQeH9GD5bmH+OXs1b5kUCmIiISQwV2bc8/FHZixeBuvL9rq+fdXKYiIhJgHLu/MhZ1S+MW7q1m69YCn31ulICISYqICxsSRPUltEMfdry5hX553Q22rFEREQlByQiyTb+7D/vxi7n19KSWl3gy1rVIQEQlRXVs15NfXdWXexn08/g9v5nhWKYiIhLARmW24qX8ak+du5INVO4P+/VQKIiIh7pFrMujRJpkH31zOpj15Qf1enpWCmQ02s/VmlmNmD1ez3sxsYsX6FWbW26tsIiKhLC46ikk39SYmyhj/yhIKikuC9r08KQUziwKeAYYAGcAoM8s4brMhQHrFx53An73IJiISDlol1+NPI3vx5e4j/PfbwZvK06sjhX5AjnNuk3OuGJgBDDtum2HAS67cAiDZzFp4lE9EJORd2CmlfETVpV/zWpAebPNqjuZWwLYqy7lA/xps0wrYUXUjM7uT8iMJ0tLSaj2oiEgom3BxR3J259E0MS4oX9+rUrBqXjv+2Kcm2+CcmwJMAcjMzPR/RgoREQ8FAsbEUb2C9/WD9pX/Uy7Qpspya2D7aWwjIiJB5FUpLAbSzay9mcUCI4HZx20zGxhTcRfSucAh59yO47+QiIgEjyenj5xzJWY2AZgDRAHTnHOrzWxcxfrJwPvAUCAHKABu8yKbiIj8H6+uKeCce5/yN/6qr02u8m8H3ONVHhER+SY90SwiIpVUCiIiUkmlICIilVQKIiJSyYI1foYXzGwPsOU0P70psLcW44QD7XPdoH2uG85kn9s651KqWxHWpXAmzCzLOZfpdw4vaZ/rBu1z3RCsfdbpIxERqaRSEBGRSnW5FKb4HcAH2ue6QftcNwRln+vsNQUREfmmunykICIix1EpiIhIpYgvBTMbbGbrzSzHzB6uZr2Z2cSK9SvMrLcfOWtTDfb5pop9XWFm88yshx85a9PJ9rnKdn3NrNTMhnuZLxhqss9mNsjMlpnZajOb63XG2laD3+2GZvZXM1tesc9hPdqymU0zs91mtuoE62v//cs5F7EflA/TvRE4C4gFlgMZx20zFPg75TO/nQss9Du3B/s8AGhU8e8hdWGfq2z3EeWj9Q73O7cHP+dkYA2QVrGc6nduD/b5Z8DvK/6dAuwHYv3Ofgb7fCHQG1h1gvW1/v4V6UcK/YAc59wm51wxMAMYdtw2w4CXXLkFQLKZtfA6aC066T475+Y55w5ULC6gfJa7cFaTnzPAvcBMYLeX4YKkJvs8GpjlnNsK4JwL9/2uyT47IMnMDEikvBRKvI1Ze5xzn1K+DydS6+9fkV4KrYBtVZZzK1471W3Cyanuzx2U/6URzk66z2bWCrgemExkqMnPuRPQyMw+MbNsMxvjWbrgqMk+Pw2cTflUviuBHzjnyryJ54taf//ybJIdn1g1rx1/D25NtgknNd4fM7uY8lIYGNREwVeTfX4SeMg5V1r+R2TYq8k+RwN9gEuBesB8M1vgnPsy2OGCpCb7fCWwDLgE6AD808w+c84dDnI2v9T6+1ekl0Iu0KbKcmvK/4I41W3CSY32x8y6A1OBIc65fR5lC5aa7HMmMKOiEJoCQ82sxDn3jicJa19Nf7f3OufygXwz+xToAYRrKdRkn28DfufKT7jnmNlmoAuwyJuInqv1969IP320GEg3s/ZmFguMBGYft81sYEzFVfxzgUPOuR1eB61FJ91nM0sDZgG3hPFfjVWddJ+dc+2dc+2cc+2At4C7w7gQoGa/2+8CF5hZtJklAP2BtR7nrE012eetlB8ZYWbNgM7AJk9TeqvW378i+kjBOVdiZhOAOZTfuTDNObfazMZVrJ9M+Z0oQ4EcoIDyvzTCVg33+RGgCTCp4i/nEhfGI0zWcJ8jSk322Tm31sw+AFYAZcBU51y1tzaGgxr+nB8FXjCzlZSfWnnIORe2Q2qb2evAIKCpmeUCvwBiIHjvXxrmQkREKkX66SMRETkFKgUREamkUhARkUoqBRERqaRSEBGRSioFERGppFIQEZFKKgWRWlQxX8MKM4s3s/oVY/p39TuXSE3p4TWRWmZmvwbiKR+ELtc591ufI4nUmEpBpJZVjMuzGCgEBjjnSn2OJFJjOn0kUvsaUz7BSxLlRwwiYUNHCiK1zMxmUz4rWHughXNugs+RRGosokdJFfFaxexmJc6518wsCphnZpc45z7yO5tITehIQUREKumagoiIVFIpiIhIJZWCiIhUUimIiEgllYKIiFRSKYiISCWVgoiIVPr/AsrsC+Z951IAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(0, 1,100)\n",
"y = 1 - x**2\n",
"plt.xlabel('x')\n",
"plt.ylabel('$d(x)$')\n",
"plt.plot(x, y)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.0, 'Monte Carlo integration (uniform sampling)')"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiHElEQVR4nO3de5hdVX3/8feHhABBINGMGJKYgcotIAKmSKVqLKgQqVi0mDQKRDDqD6naVuTSClJTaUUFnlrSQBGRyEUqldpYL2jIo5XLJIRICIFwzZiQDBeTaLgFvr8/1ho5OVlnkpnMnpOZ+bye5zxz9lr78t1r79nffTt7KyIwMzOrt0OzAzAzs+2TE4SZmRU5QZiZWZEThJmZFTlBmJlZkROEmZkVOUFYJSQ9KumYHgw3TdKPq4ipL1QZv6SdJN0n6XW9MK5N4pR0lKQHJf1O0vu3dfz9Ue06K+lcSVf20njvlHRQb4yrrzlBVCyvdC9IGlVXvkhSSGrthWnMk3T6No7jPZLmS1ovqUPSbZLet62xdVdEzImId29Nv5JOlfSLqmPqYvqteRkO7SzrTvw9MAOYHxFPbOuICnFeCPxrRLwqIv5rW8ff30XEP0XENv1P1biY1L79jhNE33gEmNrZIemNwC7NC2dTkj4IfBe4BhgL7Al8AfjzHoxr6Jb76h8kDWl2DHU+Dny7onGPB5b0ZMCBtMwrcgvwTkmjmx1It0WEPxV+gEeBvwfuqim7GDgPCKA1l+1B2kB3AI/lYXbIdacCv8jDPUNKOMflupnAS8BzwO9Ie4EABwA/AZ4GlgEnNYhPwOPA57qYhz8CfgY8BTwJzAFG1M3j54HFwPPA0Fx2TK7fCbgEWJk/lwA7NZjWqcAvaroD+ATwYJ73b+SYD8zz/FKe79/WTOviPE+rgVnALjXjOwtYleM4PY//DbnuauByYC7we+AY4L3A3cA6YAVwQc24Hs/D/y5//qQQ/1uBu4C1+e9ba+rmAf8I/BJYD/wYGNWgXV4PPAsMrRv+9O62XX2/wEPAy3n8v8ttuBdpw/Y0sBz4WM14LwBuAq7N7XJ6juVLwP/lcfw38BrSurIuz3trg3nbOY/rKeC3ud89c910YGlun4eBj9cMNwloz8t0TV6u7wcmAw/k2M8txH1DHt9C4E116/ExNf1em7+35rY8JS/zJ4HzaobbBfhWbuOlOZ72unn8CXBKs7dH3d5+NTuAgf7pXOlIG+kDgSGkDc14Nk0Q1wDfB3bLK+QDwGm57lTgReBjefhPkjZwnf/s89h0Q7FrnsZ00sb68LxSH1SI74Acx95dzMMbgHeRNhwtwHzgkrp5XASMI2+M6/7ZLgRuB16bh/8/4B8bTOtUNt/I/QAYQdpIdgDHlvrNZZeQNmyvzm3538CXc92xwBPAQcBw0t54fYJYCxxFOrrembQRemPuPoSUdN6f+2/Nww8txZ9jeAb4SF4OU3P3a2qW20PAfqSNzDzgogbt8l5gSV1Z/XLvcdvVLq/cfRvwb7kNDs3DHp3rLiCtj+/P7dIZ+3LSzsQewH2kdfiYPO/XAN9sMG8fz8tpOGn9fjOwe818/xFpp+AdwAbg8Fw3CdhIOtrdkfT/0QF8Jy/7g0g7EfvUxf3B3P/fkXa2diyssxeweYK4Is/rm0g7Qgfm+otye40kHYEvZvMEcRnwtWZvj7r78SmmvvNt4GTShvZ+4DedFflUxoeAcyJifUQ8CnyVtGHp9FhEXBERL5H2VkaTTgWVHA88GhHfjIiNEbEQ+E/SP0a91+S/qxoFHhHLI+InEfF8RHQAXyP9s9a6LCJWRMSzhVFMAy6MiDV5+C/WzduWXBQRv42Ix4GfkzZYm5Ek0kbisxHxdESsB/4JmJJ7OYm0kVoSERtyHPW+HxG/jIiXI+K5iJgXEb/O3YuB6wrz3sh7gQcj4tt5OVxHWva1p+6+GREP5Ha7sdG8kTby67dyurW2qu1qSRoH/Cnw+dwGi4Ar2XSZ/Soi/iu3S+cy/2ZEPBQRa4EfAg9FxE8jYiPpFOZhDSb5Imk9fENEvBQRCyJiHUBE/E8eZ0TEbaSjrLfVDTszIl4ErgdGAZfm/6MlpNNmh9T0vyAibsr9f42UAI/cUptkX4yIZyPiHuAeUqKAtF79U0Q8ExHtpGRQbz1pGfYrPnfYd75N2vPem7Q3VWsUMIx0aqnTY8CYmu4/XJiMiA1pW8irGkxrPPAWSb+tKRtK+fz1U/nvaNLe1GYkvZa00r+NtGe2A2lPuNaKBrFAOl1RP297ddF/vdqLshtoPN8tpL3QBbl9IO15dl5L2Atoq+m/FPMmZZLeQtpDPJi0jHYibey2Rv18QxfLla7n7RlS23fX1o6/1l5AZ4Lt9Bgwsaa71Hara74/W+huNO1vk44+r5c0gnS66byIeFHSccD5pKOsHUjL99c1wz6Vd5o6p1GKo3a6f4g7Il6W1M7Wr4uN2nIvNm2PUtvsRjp91q/4CKKPRMRjpA3wZOB7ddVPkvaExteUvZ6ao4wtjb6uewVwW0SMqPm8KiI+WRh2We7/A12M/8t5GodExO7Ah0kb3q5iqLWSzedtZRf9b636aT5J2iAcVDPfe0RE5z/yKtIpgE7jtmKc3yGdshoXEXuQrmmoQb/16ucburdcay0G9qm7IPx70gaz0zbf/pqtBF4tqTYh1cfda4+BjogXI+KLETGBdM3meOBkSTuRjnwvJl2TGEG6PlS/7nXHH5a5pB1I68O2rotbs14dSDrq6FecIPrWacCfRcTvawvzHtCNwExJu0kaD/wNaU9qa6wG9qnp/gGwn6SPSNoxf/5Y0oH1A0ZE5Gn9g6TpknaXtIOkP5U0O/e2G/lCsKQxwOe6Mc+QTsv8vaSWfLvvF7oxb11ZDYyVNCzPy8uk88Rfz0c9SBoj6T25/xuB6ZIOlDQ8x7Elu5H2pp+TdATwVzV1HaSLu/sUh0wbs/0k/ZWkoZI+BEwgLZ9uyacuHgSOqCleBJwoabikN5DWr20WEStI14m+LGlnSYfkcc/pjfHXk/ROSW/Mp1rXkXaWXuKVI7YOYGM+mtjWW4jfLOnEnGg/Q7qWcPs2jvNG4BxJI/P/x6dqK3OiezPpQnW/4gTRh/K51LYG1WeS9ggfJt2x9B3gqq0c9aXAByU9I+myfGrg3aRz7ytJh8b/TPpnK8V1E+kayEdz/6tJd6R8P/fyRdKF7rXA/7D5EdCWfIl0amcx6fTAwly2rX5GOsf8hKQnc9nnSRdLb5e0DvgpsD9ARPyQdKrs57mfX+Vhnu9iGv8PuFDSelJCubGzIl/HmAn8UtJvJW1yLjsiniLtDf8t6VTeWcDxEfEkPfPvbHod4OvAC6Tl9S16dwM+lXRxdiVwM3B+RFS1gXsd6e6idaS7gG4jXSBeD/w1qc2fISXnW7ZxWt8nreudNw+cmK9HbIsLSXdTPUJa325i03XqfcC8iOiNo+Y+1XkXjNmgk4+o7iXdcrux2fFsSd4TvZt0N1HDmwqsTNIFpAvhH654Op8EpkTEO3L3HaQ7Eu+tcrpV8BGEDSqS/kLSMEkjSUdV/90fkgNAvotsgpPD9kXS6Pyokh0k7U86Yry5sz4i3tIfkwM4Qdjg83HSOe2HSOe5SxfuzbpjGOn033rSac/vk35D0u/5FJOZmRX5CMLMzIoG1A/lRo0aFa2trc0Ow8ys31iwYMGTEdFSqhtQCaK1tZW2tkZ3kZqZWT1J9b/2/wOfYjIzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMytygjAzs6LKEoSkqyStkVR8SJWSyyQtl7RY0uE1dZ+VtETSvZKuk7RzVXGamVlZlUcQV5NeEt/IccC++TMDuBzSC15Iz4CfGBEHk14XOaXRSMzMBrNJk9KnCpUliIiYDzzdRS8nANfkl5HfDoyQNDrXDQV2yW99Gk7vvJ7SzMy6oZnXIMaw6cu924ExEfEb0jtoHye963VtRPy40UgkzZDUJqmto6Oj0oDNzAaTZiaI0ovHI7/I5QRgb2AvYFdJDd8AFRGzI2JiRExsaSk+b8rMzHqgmQmiHRhX0z2WdCrpGOCRiOjI74r9HvDWJsRnZjaoNTNB3AKcnO9mOpJ0KmkV6dTSkZKGSxJwNOlF5mZm1ocqe9y3pOuAScAoSe3A+cCOABExC5gLTAaWAxuA6bnuDkk3AQuBjaSXtM+uKk4zMyurLEFExNQt1AdwRoO680kJxczMmsS/pDYzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMytygjAzs6LKEoSkqyStkXRvg3pJukzSckmLJR1eUzdC0k2S7pe0VNKfVBWnmZmVVXkEcTVwbBf1xwH75s8M4PKaukuB/42IA4A3AUsritHMzBoYWtWII2K+pNYuejkBuCYiArg9HzWMBn4PvB04NY/nBeCFquI0M7OyZl6DGAOsqOluz2X7AB3ANyXdLelKSbs2I0Azs8GsmQlChbIgHdUcDlweEYeRjijObjgSaYakNkltHR0d1URqZjYINTNBtAPjarrHAitzeXtE3JHLbyIljKKImB0REyNiYktLS2XBmpkNNs1MELcAJ+e7mY4E1kbEqoh4Alghaf/c39HAfU2L0sxskKrsIrWk64BJwChJ7cD5wI4AETELmAtMBpYDG4DpNYOfCcyRNAx4uK7OzMz6QJV3MU3dQn0AZzSoWwRMrCAsMzPbSv4ltZmZFTlBmJlZkROEmZkVOUGYmVmRE4SZmRU5QZiZWZEThJmZFTlBmJlZkROEmZkVOUGYmVmRE4SZmRU5QZiZWZEThJmZFTlBmPWySZPSx6y/c4IwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIicIMzMrqixBSLpK0hpJ9zaol6TLJC2XtFjS4XX1QyTdLekHVcVoZmaNVXkEcTVwbBf1xwH75s8M4PK6+k8DSyuJzMzMtqiyBBER84Gnu+jlBOCaSG4HRkgaDSBpLPBe4Mqq4jMzs6418xrEGGBFTXd7LgO4BDgLeHlLI5E0Q1KbpLaOjo5eD9LMbLBqZoJQoSwkHQ+siYgFWzOSiJgdERMjYmJLS0vvRmhmNog1M0G0A+NquscCK4GjgPdJehS4HvgzSdf2fXhmZoNbMxPELcDJ+W6mI4G1EbEqIs6JiLER0QpMAX4WER9uYpxmZoPS0KpGLOk6YBIwSlI7cD6wI0BEzALmApOB5cAGYHpVsZiZWfdVliAiYuoW6gM4Ywv9zAPm9V5UZma2tfxLamvI7zUwG9ycIMzMrMgJwszMipwgzMysyAnCzMyKnCDMzKzICcLMzIqcICrmW0XNrL9ygjAzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIicIMzMrcoIwM7MiJwgzMyuqLEFIukrSGkn3NqiXpMskLZe0WNLhuXycpJ9LWippiaRPVxWjmZk1VuURxNXAsV3UHwfsmz8zgMtz+UbgbyPiQOBI4AxJEyqM08zMCipLEBExH3i6i15OAK6J5HZghKTREbEqIhbmcawHlgJjqorTzMzKmnkNYgywoqa7nbpEIKkVOAy4o9FIJM2Q1CapraOjo4o4zcwGpWYmCBXK4g+V0quA/wQ+ExHrGo0kImZHxMSImNjS0lJBmGZmg1MzE0Q7MK6meyywEkDSjqTkMCcivteE2MzMBr1mJohbgJPz3UxHAmsjYpUkAf8BLI2IrzUxPjOzQW1oVSOWdB0wCRglqR04H9gRICJmAXOBycByYAMwPQ96FPAR4NeSFuWycyNiblWxmpnZ5ipLEBExdQv1AZxRKP8F5esTZmbWh/xLajMzK3KCMDOzIicIM7N+as4cuP12uO02aG1N3b3JCcLMrB+aMwdmzIDnn0/djz2WunszSThBAJMmpY+ZWX9x3nmwYcOmZRs2pPLe4gRhZtYPPf5498p7wgliK/gIw8y2N69/fffKe6LLBCFpvaR1hc96SQ2fj2RJ1ReQzGzwmjkThg/ftGz48FTeW7r8oVxE7NZ7kxpcGl1AApg2rXlxmdnA0LkdOe20tJ0ZPz4lh97cvvgUU0X64gJSvzFkCBx66Cufiy7qvXEvWgRze+kpLPPmwfHH9864zPrAtGlw5JHwjnfAo4/2/s5nZY/aGOz64gJSv7HLLmlDXoVFi6CtDSZPrmb8ZoOYjyAq0hcXkPq1tWth//1h2bLUPXUqXHFF+v7JT8LEiXDQQXD++a8Mc9dd8Na3wpveBEcckcbxhS/ADTekI5Mbbth0Gm95CyxZ8kr3pEmwYAHceWcaz2GHpb+dMdS64AK4+OJXug8+OO2iAVx7bZr+oYfCxz8OL720TU1htr1ygqhIX1xA6jeefXbTU0w33AB77AH/+q9w6qlw/fXwzDPwsY+l/mfOTEcFixenK/yLF8MLL8CHPgSXXgr33AM//SnsuitceGEqX7Qo/a01ZQrceGP6vmoVrFwJb34zHHAAzJ8Pd9+dhj/33K2fl6VLU/y//GWa5pAhvvvABiyfYqpIX1xA6jcanWJ617vgu9+FM85IG/1ON94Is2fDxo1pw37ffSDB6NHwx3+c+tl99y1P96ST0jS++MU0zr/8y1S+di2ccgo8+GAa74svbv283HprOgrpjOPZZ+G1r9364c36ESeICk2b9spZk3nzmhrK9unll9Me+S67wNNPw9ix8Mgj6dTOXXfByJHpCOO55yAibcy7Y8wYeM1r0hHIDTfAv/97Kv+Hf4B3vhNuvjmdNir9yGXo0BRfp+eeS38jUnL58pd7MMNm/YtPMVnzfP3rcOCBcN118NGPpj35devSqaM99oDVq+GHP0z9HnBAOkV0112pe/36dISx227peyNTpsC//Es6anjjG1PZ2rUpeQBcfXV5uNZWWLgwfV+4MCUugKOPhptugjVrUvfTT6d7mM0GICcIq179NYizz4YHHoArr4SvfhXe9jZ4+9vhS19KF6APOyxdoP7oR+Goo9I4hg1LRwFnnpn6ede70l79O9+ZTkGVLlIDfPCD6RrHSSe9UnbWWXDOOWncjS4wf+ADaeN/6KFw+eWw336pfMKEFOe73w2HHJLiWLWq15rKbHui9GK3gWHixInR1tbW7eE6zzA0Og20pfptGff2as4cXz/pqf66zK1/2tb1TdKCiJhYqvMRhG2mLx4jbGbbPycI24x/BW5mUGGCkHSVpDWS7m1QL0mXSVouabGkw2vqjpW0LNedXVWMVuZfgfecH9BoA0mVRxBXA8d2UX8csG/+zAAuB5A0BPhGrp8ATJU0ocI4rY5/Bd4zPjVnA01lv4OIiPmSWrvo5QTgmkhXyW+XNELSaKAVWB4RDwNIuj73e98WJ7psWbdf3DBn9dE8sOxjPBGvpXXnNczc+wqm7XnrJv1csih/6d6ot3nYZpm5y9HM2OFzbHh55z+UDd/hOWbu8hWYdGsXQw5u591+PRuef90mZRs2wHmnPcG0K6Y0KSob6KrcxjTzGsQYYEVNd3sua1ReJGmGpDZJbS925xexpOQw44HPsSpeR7ADjz3/OmY88DnmrD66W+MZaKbteSuz9/sKo/UE4mXG7/QEs/f7ymaJ0zb1+PPlX1Q3Kjfb3lV6m2s+gvhBRBxcqPsf4MsR8YvcfStwFrAP8J6IOD2XfwQ4IiLO3NL0unuba2tr+TdO48e/8lw2GJy3uUL/jr0ZtnZ9MutNA/U213ZgXE33WGBlF+W9zhdjrTf5AY020DQzQdwCnJzvZjoSWBsRq4C7gH0l7S1pGDAl99vrfDHWetO0aekZgzvtlLrHj0/d/oGh9VeVXaSWdB3psskoSe3A+cCOABExC5gLTAaWAxuA6bluo6RPAT8ChgBXRcSSzSbQC2bOTHeZ1N7z7z0+2xZ+QKMNJFXexTR1C/UBnNGgbi4pgVRqax7J3Xlf+/PPp3PMfuSEmQ0Wg/5x313t8TW6r71zODOzgcyP2uiCHzlhZoOZE0QXfJeTmQ1mThBd8F1OZjaYOUF0wfe1m9lg5gTRBd/XbmaD2aC/i2lLfF+7mQ1WPoIwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIt/FVDHf+WRm/ZWPIMzMrMgJwszMipwgzMysyAnCzMyKnCDMzKzIdzHhO43MzEp8BGFmZkWVJghJx0paJmm5pLML9SMl3SxpsaQ7JR1cU/dZSUsk3SvpOkk7VxmrmZltqrIEIWkI8A3gOGACMFXShLrezgUWRcQhwMnApXnYMcBfAxMj4mBgCDClqljNzGxzVR5BHAEsj4iHI+IF4HrghLp+JgC3AkTE/UCrpD1z3VBgF0lDgeHAygpjNTOzOlUmiDHAipru9lxW6x7gRABJRwDjgbER8RvgYuBxYBWwNiJ+XJqIpBmS2iS1dXR09PIsmJkNXlUmCBXKoq77ImCkpEXAmcDdwEZJI0lHG3sDewG7SvpwaSIRMTsiJkbExJaWll4L3sxssKvyNtd2YFxN91jqThNFxDpgOoAkAY/kz3uARyKiI9d9D3grcG2F8ZqZWY0qjyDuAvaVtLekYaSLzLfU9iBpRK4DOB2Yn5PG48CRkobnxHE0sLTCWM3MrE5lRxARsVHSp4Afke5Cuioilkj6RK6fBRwIXCPpJeA+4LRcd4ekm4CFwEbSqafZVcVqZmabq/SX1BExF5hbVzar5vuvgH0bDHs+cH6V8W0t/9LazAYj/5LazMyKnCDMzKzICcLMzIqcIMzMrMgJwszMipwgzMysyAnCzMyKnCDMzKzICcLMzIqcIMzMrMgJwszMiip9FpP1b34Gldng5iMIMzMrcoIwM7MiJwgzMytygjAzsyInCDMzK3KCMDOzIicIMzMrqjRBSDpW0jJJyyWdXagfKelmSYsl3Snp4Jq6EZJuknS/pKWS/qTKWM3MbFOVJQhJQ4BvAMcBE4CpkibU9XYusCgiDgFOBi6tqbsU+N+IOAB4E7C0qljNzGxzVR5BHAEsj4iHI+IF4HrghLp+JgC3AkTE/UCrpD0l7Q68HfiPXPdCRPy2wljNzKxOlQliDLCiprs9l9W6BzgRQNIRwHhgLLAP0AF8U9Ldkq6UtGuFsZqZWZ0qE4QKZVHXfREwUtIi4EzgbmAj6RlRhwOXR8RhwO+Bza5hAEiaIalNUltHR0dvxW5mNuhVmSDagXE13WOBlbU9RMS6iJgeEYeSrkG0AI/kYdsj4o7c602khLGZiJgdERMjYmJLS0svz4KZ2eBVZYK4C9hX0t6ShgFTgFtqe8h3Kg3LnacD83PSeAJYIWn/XHc0cF+FsZqZWZ3KHvcdERslfQr4ETAEuCoilkj6RK6fBRwIXCPpJVICOK1mFGcCc3ICeRiYXlWsZma2uUrfBxERc4G5dWWzar7/Cti3wbCLgIlVxmdmZo35l9RmZlbkBGFmZkVOEGZmVuQEYWZmRZVepDYbjObNa3YEZr3DRxBmZlbkBGFmZkVOEGZmVuQEYWZmRU4QZmZW5ARhZmZFThBmZlbkBGFmZkVOEGZmVuQEYWZmRU4QZmZW5ARhZmZFThBmZlbkp7mamfVjVT492EcQZmZWVGmCkHSspGWSlks6u1A/UtLNkhZLulPSwXX1QyTdLekHVcZpZmabqyxBSBoCfAM4DpgATJU0oa63c4FFEXEIcDJwaV39p4GlVcVoZmaNVXkEcQSwPCIejogXgOuBE+r6mQDcChAR9wOtkvYEkDQWeC9wZYUxmplZA1UmiDHAipru9lxW6x7gRABJRwDjgbG57hLgLODlriYiaYakNkltHR0dvRC2mZlBtQlChbKo674IGClpEXAmcDewUdLxwJqIWLCliUTE7IiYGBETW1patjVmMzPLqrzNtR0YV9M9FlhZ20NErAOmA0gS8Ej+TAHeJ2kysDOwu6RrI+LDFcZrZmY1qjyCuAvYV9LekoaRNvq31PYgaUSuAzgdmB8R6yLinIgYGxGtebifOTmYmfWtyo4gImKjpE8BPwKGAFdFxBJJn8j1s4ADgWskvQTcB5xWVTxmZtY9iqi/LNB/SeoAHmtQPQp4sg/D6Q7H1jOOrWccW88M1NjGR0TxAu6AShBdkdQWERObHUeJY+sZx9Yzjq1nBmNsftSGmZkVOUGYmVnRYEoQs5sdQBccW884tp5xbD0z6GIbNNcgzMysewbTEYSZmXWDE4SZmRUN+ASxpXdS9HEs4yT9XNJSSUskfTqXXyDpN5IW5c/kJsX3qKRf5xjactmrJf1E0oP578gmxLV/TdsskrRO0mea2W6SrpK0RtK9NWUN20rSOXkdXCbpPU2I7SuS7s/vXrlZ0ohc3irp2Zo2nNWE2Boux+2g3W6oievR/Ny4Pm23LrYb1a9vETFgP6RfcD8E7AMMIz09dkIT4xkNHJ6/7wY8QHrk+QXA320H7fUoMKqu7F+As/P3s4F/3g6W6ROkJ/82rd2AtwOHA/duqa3yMr4H2AnYO6+TQ/o4tncDQ/P3f66JrbW2vya1W3E5bg/tVlf/VeALfd1uXWw3Kl/fBvoRxNa8k6LPRMSqiFiYv68nvQyp/hHo25sTgG/l798C3t+8UAA4GngoIhr9Yr5PRMR84Om64kZtdQJwfUQ8HxGPAMtJ62afxRYRP46Ijbnzdl55rH6fatBujTS93Trlh4meBFxX1fQb6WK7Ufn6NtATxNa8k6IpJLUChwF35KJP5cP/q5pxGicL4MeSFkiakcv2jIhVkFZU4LVNiq3TFDb9J90e2q1To7ba3tbDjwI/rOneW+nVvrdJeluTYiotx+2p3d4GrI6IB2vK+rzd6rYbla9vAz1BbM07KfqcpFcB/wl8JtIjzy8H/gg4FFhFOpRthqMi4nDSa2LPkPT2JsVRpPTk3/cB381F20u7bcl2sx5KOg/YCMzJRauA10fEYcDfAN+RtHsfh9VoOW437QZMZdMdkz5vt8J2o2GvhbIetdtATxBbfCdFX5O0I2khz4mI7wFExOqIeCkiXgauoMLD6K5ExMr8dw1wc45jtaTROfbRwJpmxJYdByyMiNWw/bRbjUZttV2sh5JOAY4HpkU+WZ1PQzyVvy8gna/ery/j6mI5bi/tNpT05ssbOsv6ut1K2w36YH0b6Alii++k6Ev5POZ/AEsj4ms15aNrevsL4N76Yfsgtl0l7db5nXRR815Se52SezsF+H5fx1Zjk7247aHd6jRqq1uAKZJ2krQ3sC9wZ18GJulY4PPA+yJiQ015i6Qh+fs+ObaH+zi2Rsux6e2WHQPcHxHtnQV92W6Nthv0xfrWF1fhm/kBJpOu+j8EnNfkWP6UdKi3GFiUP5OBbwO/zuW3AKObENs+pDsf7gGWdLYV8BrgVuDB/PfVTWq74cBTwB41ZU1rN1KiWgW8SNpjO62rtgLOy+vgMuC4JsS2nHReunO9m5X7/UBe3vcAC4E/b0JsDZdjs9stl18NfKKu3z5rty62G5Wvb37UhpmZFQ30U0xmZtZDThBmZlbkBGFmZkVOEGZmVuQEYWZmRU4QZhWRFJK+WtP9d5IuaGJIZt3iBGFWneeBEyWNanYgZj3hBGFWnY2kdwV/ttmBmPWEE4RZtb4BTJO0R7MDMesuJwizCkV66uY1wF83Oxaz7nKCMKveJaRnDu3a5DjMusUJwqxiEfE0cCMpSZj1G04QZn3jq4DvZrJ+xU9zNTOzIh9BmJlZkROEmZkVOUGYmVmRE4SZmRU5QZiZWZEThJmZFTlBmJlZ0f8HzIt4XEj+BtkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"uniform = np.loadtxt('data/exercise-02.1-uniform.dat')\n",
"plt.errorbar(uniform[:,0], uniform[:,1], yerr=uniform[:,2], fmt='bo')\n",
"plt.axhline(y = 1, color = 'red', linestyle = '-')\n",
"plt.text(60, 1+0.005, 'Exact value', color='red')\n",
"plt.xlabel('N')\n",
"plt.ylabel('I')\n",
"plt.title('Monte Carlo integration (uniform sampling)')"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.0, 'Monte Carlo integration (importance sampling)')"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkgUlEQVR4nO3df7xVVZ3/8dc7EEZMhRQdAuFiYYpmqDd0Ko1SE8lvplMKQ/lbtNEmp2nKH5VmMVmTX398NQmTMScCmcxiGswcG3CmieSS5IBK4k9IRFRCJwy9+Pn+sdbVzeGcC/eyzz3IeT8fj/04Z62199prnbPP/py19z77KCIwMzMrw5sa3QAzM9t+OKiYmVlpHFTMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXqStLjko7qxnITJf28Hm3qCfVsv6S+kh6Q9Oc5PUXSl+qxLiuPpLmSzsrPS9s+JP1I0tgy6iqDg0oPyTvXlyXtXpG/SFJIailhHa9ttFtRxzGS7pH0oqTVkuZJ+sjWtq2rImJ6RHxoS+aVdJqk/6p3mzpZf0t+D3t35HWl/d0wCbgnIp7O6zo3Ir5ap3V1iaSbJX2t0e3Y1pW8fVwBTC6prq3moNKzHgMmdCQkvRPYsXHN2ZikjwH/AtwCDAH2BL4M/J9u1NV783O9MUjq1eg2VDgH+OdGN6LSNvg6NYWIuBfYRVJro9sCQER46oEJeBz4IrCgkPct4BIggJactytpp74aeCIv86ZcdhrwX3m5NaQgdWwumwxsAP4E/C9wXc7fF7gLeB5YCpxUo30CngT+vpM+vA34BfAc8CwwHehf0ccvAPcD64HeOe+oXN4XuBp4Kk9XA31rrOs04L8K6QDOBR7Ofb8+t3m/3OcNud9/KKzrW7lPq4ApwI6F+j4PrMztOCvX//ZcdjNwAzAH+CNwFPBh4D7gBWA5cFmhrifz8v+bp7+o0v73AAuAtfnxPYWyucBXgV8CLwI/B3av8boMBV4Cehfybga+lp+PAVbk/j2T+/hRYBzwu7wdXFxY9jLgh8Cted2/Ad5VKN8vt+8PwBLgIxXrLb5Ok4BXgJfz6/Cveb4LgUdy/Q8AJ1S+z1TZpnP5W4B/yu/TGuDHhbLjgEW5bf8NHNjJtn1Vfj3WkrbPA3JZZ+9rS35fT89la0jb4LtzHX8gf84Kffkl8P/yeh4Cjqx4n8/qyvady3oBV5I+c48B5+f5i9vAjcCljd7PRYSDSo+90HnnStqx75c3lOXAMDYOKrcAPwF2zhv174Azc9lp+UN7dl7+U/nD1rHxvbbR5vROeR2nk3bwB+cNc/8q7ds3t2N4J314O3A0aYc9ELgHuLqij4uAvcg7cDYOKpcD84E98vL/DXy1xrqqfeh+CvQn7VhXA2OrzZvzrgZmk3ZKOwP/Cnw9l40Fngb2B/qRvvVXBpW1wHtJo/k/I+2s35nTB5IC1Ufz/C1VPuSvtSm3YQ3wyfw+TMjp3Qrv2yPAPqSR61zgihqvy4eBJRV5N7NxUGknjTB3IG0rq4Ef5Ndhf1IQ3jvPfxlpm/pYnv9zpB3XDnlaBlwM9AE+SAoM7+jkdXqtLYX2fRx4a57nZFIAGrSF2/S/kQLegNye9+f8g0lB4tC83KmkbW2TLynAMcBC0rbT8UVkUOH12tz7OiX37UP5tfsxaRsenNvw/kJf2oG/zW09Ob8+b6n8fNK17ftcUjAekl+Hf2fT7e2zwI8avZ+LcFDpuRf69aDyReDrpB3bXaSdTOQNuBfpG/7IwnLnAHPz89OAZYWyfnnZP8/p1zbanD4Z+M+KdnyHKt9oSDuGAP6sC336KHBfRR/PqNbv/PwRYFyh7Bjg8Rp1V/vQva+QngVcWGNekXZcbyvk/QXwWH4+jRxgcvrtbBpUbtlM368GrsrPW6p8yF9rEymY3Fux/K+A0wrv2xcLZX8N/KzGeicC8yvybmbjoPIS0Cund85tO7Qw/0Je33FeVqyPtHNdCRyep6fJI+VcPoP8bb7a60SVoFKlD4uA4ze3TQODgFeBAVXquIGKLySkL2zvrzLvB0lfzg4r9qUL7+vgQvlzwMmF9G3ABYW+vBYQc969wCcrP590bfv+BXBOoeyoKtvb2cAvtvSzW89puznu/Qbyz6Rv+MNJo5Ki3UnfCJ8o5D1B+kbU4emOJxGxThLAm2usaxhwqKQ/FPJ6U/14/HP5cRDpm+omJO0BXEva2exM2gGtqZhteY22QPq2Wtm3t3Yyf6WnC8/XUbvfA0k7p4X59YEUaDqO+b8VaCvMX63NG+VJOpR0QvQA0nvUl3T+aUtU9hs6eV/pvG9rSK99Z56LiA35+Uv5cVWh/KWK+l/ra0S8KmkFr78vyyPi1U7a3dn7DYCkU0jfpFty1ptJ23qHWtv0W4DnI6JyG4O0bZ8q6dOFvD5U2Z4i4heSriMdUhoq6XbgcxHxwha+r5WvXWev5e8j7+WzrmzjtbaBt7Lx61ztNd+ZdDiu4XyivodFxBOknfY44EcVxc+SDgUMK+QNBX6/pdVXpJcD8yKif2F6c0R8qsqyS/P8f9lJ/V/P6zgwInYBPkHaWXfWhqKn2LRvT3Uy/5aqXOezpA/7/oV+7xoRHR/SlaRDCR322oI6f0A6nLZXROxKOiSiGvNWquw3dO19Lbof2LvkCyFe67+kN5Fem47zXnvlvA6V7a7s+0ZpScNIx/vPJx3u6w8sZtPtpprlwFsk9a9RNrli2+4XETOqVRQR10bEIaTDf/sAf5+LOntfu2OwCt9kKGcb35LtdT/gt1u5nlI4qDTGmcAHI+KPxcz87XIWMFnSzvkD+Vng+1tY7ypg70L6p8A+kj4paYc8vVvSfpUL5m9XnwW+JOl0SbtIepOk90mammfbmXwyXNJgXv9gbqkZwBclDcyXVn+5C33rzCpgiKQ+uS+vknZkV+XRFZIGSzomzz8LOF3SfpL65XZszs6kb81/kjQa+KtC2WrSYZq9qy6ZTmTvI+mvJPWWdDIwkvT+dElErCCdzB3d1WU7cYikE3OguoB0CHY+8GvSYcTP521nDOlKwJmd1FW5De5ECjSrASSdThoVbFZErATuAL4taUBuwxG5+EbgXEmHKtlJ0oclbTKKy9v8oZJ2yP3puLADOn9fu2MP4G9yWz9O2tnP2co6ZwGfydtwf9LFMJXeT3qtGs5BpQEi4pGIaKtR/GnShv8o6aqYH5DOAWyJa4CPSVoj6dqIeJF0cnE86dvS08A3SEP8au36Iek8zBl5/lXA10gXDgB8hXSCdC3pBGrlSGtzvkY67HQ/8D+kK43K+E3DL0hXJj0t6dmc9wXSSeb5kl4gndx8B0BE3EE6jPcfeZ5f5WXWd7KOvwYul/QiKQjN6iiIiHWkq+9+KekPkg4rLhgRz5GuVPo70mHGzwPHRcSzdM93SOdpyvIT0vvecTHBiRHxSkS8DHwEOJY0+vs2cEpEPNRJXTcBI/Pr8OOIeIB05dKvSNvTO0lXSG2pT5JG7w+RTopfAJA/P2cD1+V2LyOdp6hmF1IQWkM6HPUc6Woz6OR97aZfAyNIr9dk4GP5/d8aN5KuCLyfdKXaHNIFARsgBU3gj5EuLW64jisszJpWHrktJl051N7o9myOpL6kncuR+dv81tR1GekChU+U0bZmJuk00on499V5PccCUyJiWE7fBtwUEVs7IiqFT9RbU5J0Amm0tRNp9Pavb4SAAhAR60mHz6wJSNoR+ABptLIncClwe0d5RHR2HrTH+fCXNatzSMf5HyEdRqh28YLZtkCkQ89rSCPUB9my84AN4cNfZmZWGo9UzMysNE19TmX33XePlpaWRjfDzOwNZeHChc9GxMBqZU0dVFpaWmhrq3Vlr5mZVSOp8g4Rr/HhLzMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWGgcVM7MmM2ZMmurBQcXMzErjoGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWmroFFUnTJD0jaXGNckm6VtIySfdLOrhQNlbS0lx2YSH/VkmL8vS4pEU5v0XSS4WyKfXql5mZ1da7jnXfDFwH3FKj/FhgRJ4OBW4ADpXUC7geOBpYASyQNDsiHoiIkzsWlnQlsLZQ3yMRMarsTpiZ2Zar20glIu4Bnu9kluOBWyKZD/SXNAgYDSyLiEcj4mVgZp73NZIEnATMqE/rzcysOxp5TmUwsLyQXpHzauUXHQ6sioiHC3nDJd0naZ6kw+vRYDMz61w9D39tjqrkRSf5RRPYeJSyEhgaEc9JOgT4saT9I+KFTVYqTQImAQwdOrRbDTczs+oaOVJZAexVSA8BnuokHwBJvYETgVs78iJifUQ8l58vBB4B9qm20oiYGhGtEdE6cODAkrpiZmbQ2KAyGzglXwV2GLA2IlYCC4ARkoZL6gOMz/N2OAp4KCJWdGRIGphP8CNpb9LJ/0d7qiNmZpbU7fCXpBnAGGB3SSuAS4EdACJiCjAHGAcsA9YBp+eydknnA3cCvYBpEbGkUPV4Nj1BfwRwuaR2YANwbkR0dpGAmZnVQd2CSkRM2Ex5AOfVKJtDCjrVyk6rkncbcFvXW2lmZmXyL+rNzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWGgcVMzMrjYOKmZmVxkHFzMxKU7egImmapGckLa5RLknXSlom6X5JBxfKxkpamssuLORfJun3khblaVyh7KI8/1JJx9SrX2ZmVls9Ryo3A2M7KT8WGJGnScANAJJ6Adfn8pHABEkjC8tdFRGj8jQnLzMSGA/sn9f57VyPmZn1oLoFlYi4B3i+k1mOB26JZD7QX9IgYDSwLCIejYiXgZl53s4cD8yMiPUR8RiwLNdjZmY9qJHnVAYDywvpFTmvVn6H8/PhsmmSBmymrk1ImiSpTVLb6tWrt7YPZmZW0Migoip50Uk+pENkbwNGASuBKzdT16aZEVMjojUiWgcOHNilBpuZWed6N3DdK4C9CukhwFNAnxr5RMSqjkxJNwI/3UxdZmbWgxo5UpkNnJKvAjsMWBsRK4EFwAhJwyX1IZ2Anw2Qz7l0OAFYXKhrvKS+koaTTv7f21MdMTOzpG4jFUkzgDHA7pJWAJcCOwBExBRgDjCOdFJ9HXB6LmuXdD5wJ9ALmBYRS3K135Q0inRo63HgnLzMEkmzgAeAduC8iNhQr76ZmVl1iqh66qEptLa2RltbW6ObYWbWo8aMSY9z53ZveUkLI6K1Wpl/UW9mZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmTWT6dJg/H+bNg5aWlC6Tg4qZWZOYPh0mTYL161P6iSdSuszA4qBiZtYkLrkE1q3bOG/dupRfFgcVM7Mm8eSTXcvvDgcVM7MmMXRo1/K7w0HFzKxJTJ4M/fptnNevX8ovi4OKmVmTmDgRpk6Fvn1TetiwlJ44sbx1NPLW92Zm1sMmToQbb0zPu3vvr854pGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK42DipmZlcZBxczMSuOgYmZmpalbUJE0TdIzkhbXKJekayUtk3S/pIMLZWMlLc1lFxby/1HSQ3n+2yX1z/ktkl6StChPU+rVL7MyjBmTJrPtTT1HKjcDYzspPxYYkadJwA0AknoB1+fykcAESSPzMncBB0TEgcDvgIsK9T0SEaPydG6ZHTEzsy1Tt6ASEfcAz3cyy/HALZHMB/pLGgSMBpZFxKMR8TIwM89LRPw8Itrz8vOBIfVqv5mZdV0jz6kMBpYX0ityXq38SmcAdxTSwyXdJ2mepMNrrVTSJEltktpWr17d/dabmdkmGhlUVCUvOsl/fUHpEqAd6Phn5ZXA0Ig4CPgs8ANJu1RbaURMjYjWiGgdOHBgtxtvZmabauSt71cAexXSQ4CngD418gGQdCpwHHBkRARARKwH1ufnCyU9AuwDtNWzA2ZmtrFGjlRmA6fkq8AOA9ZGxEpgATBC0nBJfYDxeV4kjQW+AHwkItZ1VCRpYD7Bj6S9SSf/H+3Z7piZWd1GKpJmAGOA3SWtAC4FdgCIiCnAHGAcsAxYB5yey9olnQ/cCfQCpkXEklztdUBf4C5JAPPzlV5HAJdLagc2AOdGRGcXCZiZWR3ULahExITNlAdwXo2yOaSgU5n/9hrz3wbc1o1mmplZifyLejMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWGgcVMzMrjYOKmZmVxkHFzMxK46Bi1sx69YJRo16frriivLoXLYI5m9zCr3vmzoXjjiunLqurRv6fipk12o47pp1/PSxaBG1tMG5cfeq3bZJHKma2sbVr4R3vgKVLU3rCBLjxxvT8U5+C1lbYf3+49NLXl1mwAN7zHnjXu2D06FTHl78Mt96aRkC33rrxOg49FJYseT09ZgwsXAj33pvqOeig9NjRhqLLLoNvfev19AEHwOOPp+ff/35a/6hRcM45sGHDVr0U1nUOKmbN7KWXNj78deutsOuucN11cNppMHMmrFkDZ5+d5p88OY0+7r8f5s1Ljy+/DCefDNdcA7/9Lfz7v8NOO8Hll6f8RYvSY9H48TBrVnq+ciU89RQccgjsuy/ccw/cd19a/uKLt7wvDz6Y2v/LX6Z19uoF06dvdjErlw9/mTWzWoe/jj4a/uVf4LzzUqDoMGsWTJ0K7e0pGDzwAEgwaBC8+91pnl122fx6TzopreMrX0l1fvzjKX/tWjj1VHj44VTvK69seV/uvjuNdjra8dJLsMceW768lcJBxcw29eqr6Zv/jjvC88/DkCHw2GPpsNOCBTBgQBrJ/OlPEJECQFcMHgy77ZZGOrfeCt/5Tsr/0pfgAx+A229Ph7TGjNl02d69U/s6/OlP6TEiBaSvf70bHbaydHr4S9KLkl6oMr0o6YWeaqSZ9bCrroL99oMZM+CMM9KI4YUX0mGtXXeFVavgjjvSvPvumw5fLViQ0i++mEYyO++cntcyfjx885tpdPLOd6a8tWtTwAG4+ebqy7W0wG9+k57/5jcp2AEceST88IfwzDMp/fzz8MQT3X0FrJs6DSoRsXNE7FJl2jkitmCMa2bbtMpzKhdeCL/7HXz3u3DllXD44XDEEfC1r6WT8AcdlE7Sn3EGvPe9qY4+fdJo49OfTvMcfXQaPXzgA+nwWLUT9QAf+1g6Z3PSSa/nff7zcNFFqe5aJ9n/8i9TwBg1Cm64AfbZJ+WPHJna+aEPwYEHpnasXFnaS2VbRumv4ptTa2trtLW1NboZ1oQ6jurMndvIVliz2trtT9LCiGitVla3q78kTZP0jKTFNcol6VpJyyTdL+ngQtlYSUtz2YWF/LdIukvSw/lxQKHsojz/UknH1KtfZmZWWz0vKb4ZGNtJ+bHAiDxNAm4AkNQLuD6XjwQmSBqZl7kQuDsiRgB35zS5fDywf17nt3M9ZmbWg+oWVCLiHuD5TmY5HrglkvlAf0mDgNHAsoh4NCJeBmbmeTuW+V5+/j3go4X8mRGxPiIeA5bleszMrAc18pLiwcDyQnpFzquWf2h+vmdErASIiJWSOi5CHwzMr1JX55YurX7JolmdXb0oPxnTwEZY06rn9tfIX9RXu7A9OsnvTl2bzihNktQmqe2VrvywyszMNquRI5UVwF6F9BDgKaBPjXyAVZIG5VHKIOCZzdS1iYiYCkyFdPWXL7+xRrhgTHr05meNsNXbXyc/dm3kSGU2cEq+CuwwYG0+tLUAGCFpuKQ+pBPwswvLnJqfnwr8pJA/XlJfScNJJ//v7amOmJlZUreRiqQZpCN2u0taAVwK7AAQEVOAOcA40kn1dcDpuaxd0vnAnUAvYFpEdNzO9ApglqQzgSeBj+dllkiaBTwAtAPnRYRvT2pm1sPqFlQiYsJmygM4r0bZHFLQqcx/DjiyxjKTgcldb6mZmZXFt743M7PSOKiYmVlpHFTMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVp5L2/zMysAep5zzmPVJrUmDG+67+Zlc9BxczMSuOgYtbDpk+H+fNh3jxoaUlps+2Fg4pZD5o+HSZNgvXrU/qJJ1LagcW2Fw4qZj3okktg3bqN89atS/lm2wMHFbMe9OSTXcs3e6NxUDHrQUOHdi3f7I3GQcWsB02eDP36bZzXr1/KN9seOKiY9aCJE2HqVOjbN6WHDUvpiRMb2y6zsvgX9WY9bOJEuPHG9Lyev2w2awSPVMzMrDQOKmZmVpq6BhVJYyUtlbRM0oVVygdIul3S/ZLulXRAoewzkhZLWiLpgkL+rZIW5elxSYtyfouklwplU+rZNzMz21TdzqlI6gVcDxwNrAAWSJodEQ8UZrsYWBQRJ0jaN89/ZA4uZwOjgZeBn0n6t4h4OCJOLqzjSmBtob5HImJUvfpkZmadq+dIZTSwLCIejYiXgZnA8RXzjATuBoiIh4AWSXsC+wHzI2JdRLQD84ATigtKEnASMKOOfTAzsy6oZ1AZDCwvpFfkvKLfAicCSBoNDAOGAIuBIyTtJqkfMA7Yq2LZw4FVEfFwIW+4pPskzZN0eLVGSZokqU1S2+rVq7vbNzMzq6KeQUVV8qIifQUwIJ8X+TRwH9AeEQ8C3wDuAn5GCj7tFctOYONRykpgaEQcBHwW+IGkXTZpQMTUiGiNiNaBAwd2vVdmZlZTPX+nsoKNRxdDgKeKM0TEC8Dp8NrhrMfyRETcBNyUy/4h10dO9yaNcA4p1LUeWJ+fL5T0CLAP0FZyv8zMrIZ6jlQWACMkDZfUBxgPzC7OIKl/LgM4C7gnBxok7ZEfh5ICSHFUchTwUEQUA83AfHEAkvYGRgCP1qVnZmZWVd1GKhHRLul84E6gFzAtIpZIOjeXTyGdkL9F0gbgAeDMQhW3SdoNeAU4LyLWFMrGs+kJ+iOAyyW1AxuAcyPi+Xr0zczMqqvrbVoiYg4wpyJvSuH5r0gjimrLVj3RnstOq5J3G3Bbd9tqZmZbz7+oNzOz0jiomJlZaRxUzMysNA4qZmZWGgcVMzMrjYOKmZmVxkHFzMxK46BiZmalcVAxM7PSOKiYmVlpHFQaZMyYNDXC9Okwfz7MmwctLSltZlYGB5UmM306TJoE69en9BNPpLQDi5mVwUGlyVxyCaxbt3HeunUp38xsazmoNJknn+xavplZVzioNJmhQ7uWb2bWFQ4qTWbyZOjXb+O8fv1SvpnZ1nJQaTITJ8LUqdC3b0oPG5bSEyc2tl1mtn2o6z8/2rZp4kS48cb0fO7chjbFzLYzHqmYmVlpHFTMzKw0DipmZlaaugYVSWMlLZW0TNKFVcoHSLpd0v2S7pV0QKHsM5IWS1oi6YJC/mWSfi9pUZ7GFcouyutaKumYevbNzMw2VbegIqkXcD1wLDASmCBpZMVsFwOLIuJA4BTgmrzsAcDZwGjgXcBxkkYUlrsqIkblaU5eZiQwHtgfGAt8O7fBzMx6SD1HKqOBZRHxaES8DMwEjq+YZyRwN0BEPAS0SNoT2A+YHxHrIqIdmAecsJn1HQ/MjIj1EfEYsCy3wczMekg9g8pgYHkhvSLnFf0WOBFA0mhgGDAEWAwcIWk3Sf2AccBeheXOz4fMpkka0IX1IWmSpDZJbatXr+5+78zMbBP1DCqqkhcV6SuAAZIWAZ8G7gPaI+JB4BvAXcDPSMGnPS9zA/A2YBSwEriyC+sjIqZGRGtEtA4cOLAr/TEzs82o548fV7Dx6GII8FRxhoh4ATgdQJKAx/JERNwE3JTL/iHXR0Ss6lhe0o3AT7d0fWZmVl/1HKksAEZIGi6pD+kk+uziDJL65zKAs4B7cqBB0h75cSjpENmMnB5UqOIE0qEyct3jJfWVNBwYAdxbl56ZmVlVdRupRES7pPOBO4FewLSIWCLp3Fw+hXRC/hZJG4AHgDMLVdwmaTfgFeC8iFiT878paRTp0NbjwDm5viWSZuV62vMyG+rVPzMz21Rd7/2VL/edU5E3pfD8V6QRRbVlD6+R/8lO1jcZ2Obvt9vxd77r16e/85082Td0NLPtg39R38P8d75mtj1zUOlh/jtfM9ueOaj0MP+dr5ltzxxUepj/ztfMtmcOKj3Mf+drZtszB5VuGjMmTV3lv/M1s+2Z/064Afx3vub33bZXHqmYmVlpHFTMzKw0DipmZlYaBxUzMyuNg4qZmZXGQaUbOm4IOW9euiGk79tlZpY4qHSRbwhpZlabg0oX+YaQZma1Oah0kW8IaWZWm4NKF/mGkGZmtTmodNH2ckPIuXN9qxAzK5+DShf5hpBmZrX5hpLd4BtCmplV55GKmZmVpq5BRdJYSUslLZN0YZXyAZJul3S/pHslHVAo+4ykxZKWSLqgkP+Pkh7Ky9wuqX/Ob5H0kqRFeZpSz76Zmdmm6hZUJPUCrgeOBUYCEySNrJjtYmBRRBwInAJck5c9ADgbGA28CzhO0oi8zF3AAXmZ3wEXFep7JCJG5encOnXNzMxqqOdIZTSwLCIejYiXgZnA8RXzjATuBoiIh4AWSXsC+wHzI2JdRLQD84AT8nw/z3kA84EhdeyDmZl1QT2DymBgeSG9IucV/RY4EUDSaGAYKUgsBo6QtJukfsA4YK8q6zgDuKOQHi7pPknzJB1erVGSJklqk9S2evXq7vTLzMxqqOfVX6qSFxXpK4BrJC0C/ge4D2iPiAclfYN0qOt/ScGnvbigpEtyXsddt1YCQyPiOUmHAD+WtH9EvLBRAyKmAlMBWltbK9uzxbb2qi9fNWZm26N6BpUVbDy6GAI8VZwh7/BPB5Ak4LE8ERE3ATflsn/I9ZHTpwLHAUdGROT51wPr8/OFkh4B9gHa6tA3MzOrop6HvxYAIyQNl9QHGA/MLs4gqX8uAzgLuKdjZCFpj/w4lHSIbEZOjwW+AHwkItYV6hqYLw5A0t7ACODROvbPzMwq1G2kEhHtks4H7gR6AdMiYomkc3P5FNIJ+VskbQAeAM4sVHGbpN2AV4DzImJNzr8O6AvclQY3zM9Xeh0BXC6pHdgAnBsRz9erf2Zmtinlo0dNqbW1NdrafHTMzKwrJC2MiNZqZf5FvZmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaZr66i9Jq4EnGt2OBtodeLbRjWgg99/9d/+7Z1hEDKxW0NRBpdlJaqt1WWAzcP/df/e//P778JeZmZXGQcXMzErjoNLcpja6AQ3m/jc3978OfE7FzMxK45GKmZmVxkHFzMxK46DSBCTtJek/JD0oaYmkz+T8t0i6S9LD+XFAo9taT5J65b+b/mlON03/838X/VDSQ3k7+Ism6//f5m1/saQZkv5se+6/pGmSnpG0uJBXs7+SLpK0TNJSScdszbodVJpDO/B3EbEfcBhwnqSRwIXA3RExArg7p7dnnwEeLKSbqf/XAD+LiH2Bd5Feh6bov6TBwN8ArRFxAOn/ncazfff/ZmBsRV7V/uZ9wXhg/7zMtzv+8LA7HFSaQESsjIjf5OcvknYog4Hjge/l2b4HfLQhDewBkoYAHwa+W8huiv5L2oX0J3Y3AUTEyxHxB5qk/1lvYEdJvYF+pL823277HxH3AJV/Ulirv8cDMyNifUQ8BiwDRnd33Q4qTUZSC3AQ8Gtgz4hYCSnwAHs0sGn1djXweeDVQl6z9H9vYDXwT/nw33cl7UST9D8ifg98C3gSWAmsjYif0yT9L6jV38HA8sJ8K3JetzioNBFJbwZuAy6IiBca3Z6eIuk44JmIWNjotjRIb+Bg4IaIOAj4I9vXoZ5O5XMHxwPDgbcCO0n6RGNbtU1Rlbxu/9bEQaVJSNqBFFCmR8SPcvYqSYNy+SDgmUa1r87eC3xE0uPATOCDkr5P8/R/BbAiIn6d0z8kBZlm6f9RwGMRsToiXgF+BLyH5ul/h1r9XQHsVZhvCOnwYLc4qDQBSSIdT38wIv5voWg2cGp+firwk55uW0+IiIsiYkhEtJBOSP4iIj5B8/T/aWC5pHfkrCOBB2iS/pMOex0mqV/+LBxJOq/YLP3vUKu/s4HxkvpKGg6MAO7t7kr8i/omIOl9wH8C/8Pr5xQuJp1XmQUMJX3wPh4RlSf3tiuSxgCfi4jjJO1Gk/Rf0ijSRQp9gEeB00lfKpul/18BTiZdCXkfcBbwZrbT/kuaAYwh3d5+FXAp8GNq9FfSJcAZpNfngoi4o9vrdlAxM7Oy+PCXmZmVxkHFzMxK46BiZmalcVAxM7PSOKiYmVlpHFTMtiGSQtKVhfTnJF3WwCaZdYmDitm2ZT1woqTdG90Qs+5wUDHbtrST/jv8bxvdELPucFAx2/ZcD0yUtGujG2LWVQ4qZtuYfAfpW0h/LGX2huKgYrZtuho4E9ipwe0w6xIHFbNtUL7R3yxSYDF7w3BQMdt2XUm6y6zZG4bvUmxmZqXxSMXMzErjoGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK83/B5cuTsDj/XxjAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"importance = np.loadtxt('data/exercise-02.1-importance.dat')\n",
"plt.errorbar(importance[:5,0], importance[:5,1], yerr=importance[:5,2], fmt='bo')\n",
"plt.axhline(y = 1, color = 'red', linestyle = '-')\n",
"plt.text(60, 1+0.0005, 'Exact value', color='red')\n",
"plt.xlabel('N')\n",
"plt.ylabel('I')\n",
"plt.title('Monte Carlo integration (importance sampling)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 02.2\n",
"- <span style=\"color:red\">3D Random Walks (RW) on a cubic lattice and in the continuum:</span> Repeat many times (e.g. say $10^4$) the simulation of a random walk in 3D always starting at the origin"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"1. on a cubic lattice with lattice constant $a=1$; at each discrete time the walker makes a forward or backward step of length equal to $a$ in one of the 3 principal directions of the lattice: $x$, $y$ or $z$\n",
"2. in the continuum; at each discrete time the walker makes a step of length equal to $a(=1)$ along a **random direction** obtained by sampling **uniformly** the solid angle: $\\theta \\in [0,\\pi]$ and $\\phi \\in [0,2\\pi]$\n",
"\n",
"<span style=\"color:red\">Show a picture of $\\sqrt{\\langle |\\vec{r}_N|^2 \\rangle_{RW}}$ for both RWs</span>, with their statistical uncertainties, as function of the step $i\\in [0,10^2]$.\n",
"\n",
"Suggestion: divide your $M$ throws into a reasonable number of blocks, compute $\\sqrt{\\langle |\\vec{r}_N|^2 \\rangle_{RW}}$ as function of the step $i\\in [0,10^2]$ in each block and use these values (for fixed $i$) to obtain the average value and its statistical uncertainty.\n",
"\n",
"Note that you could try to fit both results with a function like $f(N)=k\\sqrt{N}$. Do your results indicate a diffusive behavior?"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "/gnu/store/3lxr2xg3yscdb3979blgjg0h7xd1n9la-python-3.10.7/bin/python3",
"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.10.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
|
