{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Numbers, Histograms, and a Simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Random numbers are often useful both for simulation of physical processes and for generating a collection of test cases.\n",
    "Here we will do a mathematical simulation: approximating $\\pi$ on the basis that the unit circle occupies a fraction $\\pi/4$ of the $2 \\times 2$ square enclosing it.\n",
    "\n",
    "### Disclaimer\n",
    "\n",
    "Actually, the best we have available is **pseudo-random** numbers, generated by algorithms that actually produce a very long but eventually repeating sequence of numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Module `random` within package `numpy`\n",
    "\n",
    "The pseudo-random number generator we use are provided by package Numpy in its module random – full name `numpy.random`.\n",
    "This module contains numerous random number generators; here we look at just a few.\n",
    "\n",
    "We introduce the abbreviation \"npr\" for this, along with the standard abbreviations \"np\" for `numpy` and \"plt\" for module `matplotlib.pyplot` witit package `matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.random as npr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniformly distributed real numbers: `numpy.random.rand`\n",
    "\n",
    "First, the function `rand` (full name `numpy.random.rand`) generates uniformly-distributed real numbers in the semi-open interval $[0,1)$.\n",
    "\n",
    "To generate a single value, use it with no argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6831120389814106\n",
      "0.5548018157238266\n",
      "0.8695626004658878\n",
      "0.3752686546706374\n"
     ]
    }
   ],
   "source": [
    "n_samples = 4\n",
    "for sample in range(n_samples):\n",
    "    print(npr.rand())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrays of random values\n",
    "\n",
    "To generate an array of values all at once, one can specify how many as the first and only input argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.62144348 0.03793401 0.13891253 0.84806369]\n"
     ]
    }
   ],
   "source": [
    "pseudorandomnumbers = npr.rand(n_samples)\n",
    "print(pseudorandomnumbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the first method has an advantage in some situations: neither the whole list of integers from 0 to `n_samples - 1` nor the collection of random numbers is stored at any time: instead, just one value at a time is provided, used, and then \"forgotten\".\n",
    "This can be beneficial or even essential when a very large number of random values is used; it is not unusual for a simulation to require more random values than the computer's memory can hold."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-dimensional arrays of random values\n",
    "\n",
    " We can also generate multi-dimensional arrays, by giving the lengths of each dimension as arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A two-dimensional array of random numbers:\n",
      " [[0.18395217 0.71429007 0.68437103]\n",
      " [0.00242336 0.92888447 0.14129646]]\n",
      "\n",
      "A three-dimensional array of random numbers:\n",
      " [[[0.2511712  0.85518662 0.36048736 0.11214513]\n",
      "  [0.00948166 0.91241279 0.07428897 0.18599767]\n",
      "  [0.75078587 0.31421873 0.00762302 0.15589679]]\n",
      "\n",
      " [[0.78636402 0.48269234 0.08894184 0.33770037]\n",
      "  [0.66866924 0.85747635 0.82185659 0.83088524]\n",
      "  [0.98993293 0.10828591 0.78003614 0.53410965]]]\n"
     ]
    }
   ],
   "source": [
    "numbers2d = npr.rand(2,3)\n",
    "print('A two-dimensional array of random numbers:\\n', numbers2d)\n",
    "numbers3d = npr.rand(2,3,4)\n",
    "print('\\nA three-dimensional array of random numbers:\\n', numbers3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normally distributed real numbers: `numpy.random.randn`\n",
    "\n",
    "The function `randn` has the same interface, but generates numbers with the standard normal distribution of mean zero, standard deviation one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ten normally distributed values:\n",
      " [ 1.01125557 -0.51308499  1.47690102 -0.1352893   0.03629515  1.06888055\n",
      "  0.17867389  0.03297391  1.62309016  0.38036881 -1.45193537 -0.16739935\n",
      "  1.52999647 -0.24956606  0.62335861 -0.23251452  0.4546426  -0.61267693\n",
      "  1.24106299  0.52917226]\n"
     ]
    }
   ],
   "source": [
    "print('Ten normally distributed values:\\n', npr.randn(20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of these 10000000 samples is 0.00032428885458364956\n",
      "and their standard deviation is 1.0002386370861691\n"
     ]
    }
   ],
   "source": [
    "n_samples = 10**7\n",
    "normf_samples = npr.randn(n_samples)\n",
    "mean = sum(normf_samples)/n_samples\n",
    "print('The mean of these', n_samples, 'samples is', mean)\n",
    "standard_deviation = np.sqrt(sum(normf_samples**2)/n_samples - mean**2)\n",
    "print('and their standard deviation is', standard_deviation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** The exact mean and standard deviation of the standard normal distribtion are 0 and 1 respectively,\n",
    "so the slight variations above are due to these being only a sample mean and sample standard deviation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Histogram plots with `matplotlib.random.hist`\n",
    "\n",
    "`matplotlib.random` has a function `hist(x, bins, ...)` for plotting histograms, so we can check what this normal distribution actually looks like.\n",
    "\n",
    "Input parameter `x` is the list of values, and when input parameter `bins` is given an integer value, the data is binned into that many equally wide intervals.\n",
    "\n",
    "The function `hist` also returns three values:\n",
    "- `n`, the number of values in each bin (the bar heights on the histogram)\n",
    "- `bins` (which I prefer to call `bin_edges`), the list of values of the edges between the bins\n",
    "- `patches`, which we can ignore!\n",
    "\n",
    "It is best to assigned this output to variables; otherwise the numerous values are sprayed over the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAX5klEQVR4nO3df6zldZ3f8eerINbWZUG5UDIzdlg7uxHtdpQbnMS4tbLCAMbBRltoKlOXZNRColk3FfUPjLrJ2I1rS6oYlAlDqyD1R5i4UJxSsqYJIBdBfojsXHBWrkxhdFBp3GLGffeP87numeHM996599x7ztz7fCQn53ve38/3ez7fzJ3zOt/P9/s931QVkiQdyd8bdQckSePNoJAkdTIoJEmdDApJUieDQpLU6fhRd2DYTjnllFq/fv2ouyFJx5T77rvvJ1U1MWjeiguK9evXMzU1NepuSNIxJclfH2meQ0+SpE4GhSSpk0EhSepkUEiSOhkUkqROBoUkqZNBIUnqZFBIkjoZFJKkTivuymxpHK2/8i8G1vduv3CZeyIdPfcopBFaf+VfHDFEpHFhUEhLaL5BYFhonBkUkqROBoW0RI52L8FhKI0rg0KS1MmgkMaMexUaN54eKw2ZH/Raaebco0iyLsmdSR5N8kiS97f6y5LsTrKnPZ/c6klydZLpJA8meV3fura29nuSbO2rn5XkobbM1UnS9R6SpOUzn6Gng8AHq+pVwCbg8iRnAlcCd1TVBuCO9hrgfGBDe2wDroHehz5wFfB64Gzgqr4P/mta29nlNrf6kd5DWtE8sK1xMmdQVNW+qvpum34OeBRYA2wBdrZmO4GL2vQW4IbquRs4KcnpwHnA7qo6UFXPAruBzW3eiVV1V1UVcMNh6xr0HpKkZXJUB7OTrAdeC9wDnFZV+6AXJsCprdka4Mm+xWZaras+M6BOx3sc3q9tSaaSTO3fv/9oNkkaKvcCtBLNOyiSvBT4GvCBqvpFV9MBtVpAfd6q6tqqmqyqyYmJiaNZVJI0h3kFRZIX0QuJL1XV11v56TZsRHt+ptVngHV9i68FnpqjvnZAves9pLGyVMcU3EPROJjPWU8BrgMerao/75u1C5g9c2krcEtf/dJ29tMm4Odt2Oh24NwkJ7eD2OcCt7d5zyXZ1N7r0sPWNeg9JEnLZD7XUbwBeBfwUJIHWu0jwHbg5iSXAT8C3tnm3QpcAEwDvwTeDVBVB5J8Ari3tft4VR1o0+8DrgdeAtzWHnS8hyRpmcwZFFX1vxl8HAHgnAHtC7j8COvaAewYUJ8CXjOg/tNB7yFJWj7+hIckqZNBIS2SB5y10hkU0pjzKm2NmkEhSepkUEiSOhkU0gIt95CQw08aFYNCktTJoJAkdTIoJEmdDApJUieDQlqAUR1Y9poKjYJBIUnqZFBIkjoZFJKkTgaFJKnTfO5wtyPJM0ke7qt9JckD7bF39oZGSdYn+Zu+eZ/vW+asJA8lmU5ydbubHUlelmR3kj3t+eRWT2s3neTBJK8b/uZLR2dcDiaPQx+0esxnj+J6YHN/oar+dVVtrKqN9O6l/fW+2Y/Pzquq9/bVrwG2ARvaY3adVwJ3VNUG4I72GuD8vrbb2vKSpGU2Z1BU1beBA4Pmtb2CfwXc2LWOJKcDJ1bVXe0OeDcAF7XZW4CdbXrnYfUbqudu4KS2HknSMlrsMYo3Ak9X1Z6+2hlJ7k/yl0ne2GprgJm+NjOtBnBaVe0DaM+n9i3z5BGWOUSSbUmmkkzt379/cVskSTrEYoPiEg7dm9gHvKKqXgv8MfDlJCcy+J7bNce6571MVV1bVZNVNTkxMTGPbkuS5uv4hS6Y5HjgXwJnzdaq6nng+TZ9X5LHgd+ltzewtm/xtcBTbfrpJKdX1b42tPRMq88A646wjLTqzR7Q3rv9whH3RCvdYvYo/hD4QVX9ZkgpyUSS49r079A7EP1EG1J6LsmmdlzjUuCWttguYGub3npY/dJ29tMm4OezQ1TSKHimkVar+ZweeyNwF/B7SWaSXNZmXcwLD2L/AfBgku8BXwXeW1WzB8LfB3wRmAYeB25r9e3AW5LsAd7SXgPcCjzR2n8B+PdHv3mSpMWac+ipqi45Qv3fDah9jd7psoPaTwGvGVD/KXDOgHoBl8/VP0nS0vLKbElSJ4NCktTJoJCOcR5k11Jb8Omx0mrhB7FWO/coJEmdDApJUieDQpLUyaCQJHUyKKQVYFxuqKSVyaCQOvjhKxkUkqQ5GBSSpE4GhSSpk0EhSepkUEgriAfftRT8rSdpAD9wpb8znzvc7UjyTJKH+2ofS/LjJA+0xwV98z6cZDrJY0nO66tvbrXpJFf21c9Ick+SPUm+kuSEVn9xez3d5q8f1kZLkuZvPkNP1wObB9Q/U1Ub2+NWgCRn0rtF6qvbMp9Lcly7j/ZngfOBM4FLWluAT7V1bQCeBWZvtXoZ8GxV/RPgM62dJGmZzRkUVfVt4MBc7ZotwE1V9XxV/ZDe/a7Pbo/pqnqiqn4F3ARsSRLgzfTurw2wE7iob1072/RXgXNae0nSMlrMwewrkjzYhqZObrU1wJN9bWZa7Uj1lwM/q6qDh9UPWVeb//PW/gWSbEsylWRq//79i9gkSdLhFhoU1wCvBDYC+4BPt/qgb/y1gHrXul5YrLq2qiaranJiYqKr35Kko7SgoKiqp6vq11X1t8AX6A0tQW+PYF1f07XAUx31nwAnJTn+sPoh62rzf5v5D4FJq5Y/EKhhW1BQJDm97+XbgdkzonYBF7czls4ANgDfAe4FNrQznE6gd8B7V1UVcCfwjrb8VuCWvnVtbdPvAP5Xay8tKT9kpUPNeR1FkhuBNwGnJJkBrgLelGQjvaGgvcB7AKrqkSQ3A98HDgKXV9Wv23quAG4HjgN2VNUj7S0+BNyU5JPA/cB1rX4d8F+TTNPbk7h40VsrSTpqcwZFVV0yoHzdgNps+z8F/nRA/Vbg1gH1J/i7oav++v8D3jlX/yRJS8uf8JAkdTIoJEmdDApphfKgvIbFHwWUGj9YpcHco5AkdTIoJEmdDApJUieDQpLUyaCQJHUyKKQVzB8I1DAYFBKeGit1MSgkSZ0MCklSJ4NCktTJoJAkdZozKJLsSPJMkof7an+W5AdJHkzyjSQntfr6JH+T5IH2+HzfMmcleSjJdJKrk6TVX5Zkd5I97fnkVk9rN93e53XD33xJ0lzms0dxPbD5sNpu4DVV9fvAXwEf7pv3eFVtbI/39tWvAbbRuz3qhr51XgncUVUbgDvaa4Dz+9pua8tLQ7VaTh9dDduopTNnUFTVt+ndirS/9q2qOthe3g2s7VpHu8f2iVV1V7vv9Q3ARW32FmBnm955WP2G6rkbOOmwe3VLkpbBMI5R/BFwW9/rM5Lcn+Qvk7yx1dYAM31tZloN4LSq2gfQnk/tW+bJIyxziCTbkkwlmdq/f//itkaSdIhFBUWSjwIHgS+10j7gFVX1WuCPgS8nORHIgMVrrtXPd5mquraqJqtqcmJiYn6dlyTNy4JvXJRkK/BW4Jw2nERVPQ8836bvS/I48Lv09gb6h6fWAk+16aeTnF5V+9rQ0jOtPgOsO8IykqRlsqA9iiSbgQ8Bb6uqX/bVJ5Ic16Z/h96B6CfakNJzSTa1s50uBW5pi+0CtrbprYfVL21nP20Cfj47RCVJWj5z7lEkuRF4E3BKkhngKnpnOb0Y2N3Ocr27neH0B8DHkxwEfg28t6pmD4S/j94ZVC+hd0xj9rjGduDmJJcBPwLe2eq3AhcA08AvgXcvZkMlSQuTNmq0YkxOTtbU1NSou6FjwGo9ZXTv9gtH3QWNoST3VdXkoHlemS1J6mRQSJI6GRSSpE4GhSSpk0EhSepkUEirzGo920sLt+Ars6VjlR+U0tFxj0KS1MmgkCR1MigkSZ0MCklSJ4NCWoVWyy1gNRwGhSSpk0GhVcVv0dLRMygkSZ0MCklSp3kFRZIdSZ5J8nBf7WVJdifZ055PbvUkuTrJdJIHk7yub5mtrf2eds/t2fpZSR5qy1zdbpd6xPeQJC2f+e5RXA9sPqx2JXBHVW0A7mivAc6nd6/sDcA24BrofejTu43q64Gzgav6PvivaW1nl9s8x3tIkpbJvIKiqr4NHDisvAXY2aZ3Ahf11W+onruBk5KcDpwH7K6qA1X1LLAb2NzmnVhVd1Xvvqw3HLauQe8haQg8TVbzsZgfBTytqvYBVNW+JKe2+hrgyb52M63WVZ8ZUO96j0Mk2UZvj4RXvOIVi9gkrVR+GEoLtxQHszOgVguoz1tVXVtVk1U1OTExcTSLSpLmsJigeLoNG9Gen2n1GWBdX7u1wFNz1NcOqHe9hyRpmSwmKHYBs2cubQVu6atf2s5+2gT8vA0f3Q6cm+TkdhD7XOD2Nu+5JJva2U6XHrauQe8hSVom8zpGkeRG4E3AKUlm6J29tB24OcllwI+Ad7bmtwIXANPAL4F3A1TVgSSfAO5t7T5eVbMHyN9H78yqlwC3tQcd7yFJWibpnWi0ckxOTtbU1NSou6Ex4UHs+dm7/cJRd0EjluS+qpocNM8rsyUZqOpkUEiSOhkUkqROBoUkqZNBIUnqZFBIAvzdJx2ZQaEVyw89aTgMCklSJ4NCktTJoJAkdTIoJEmdDApJh/AkAB1uMXe4k8aSH3TScLlHIUnqZFBIkjotOCiS/F6SB/oev0jygSQfS/LjvvoFfct8OMl0kseSnNdX39xq00mu7KufkeSeJHuSfCXJCQvfVEnSQiw4KKrqsaraWFUbgbPo3c3uG232Z2bnVdWtAEnOBC4GXg1sBj6X5LgkxwGfBc4HzgQuaW0BPtXWtQF4Frhsof2VJC3MsIaezgEer6q/7mizBbipqp6vqh/Su1Xq2e0xXVVPVNWvgJuALe3+2W8GvtqW3wlcNKT+SpLmaVhBcTFwY9/rK5I8mGRHkpNbbQ3wZF+bmVY7Uv3lwM+q6uBh9RdIsi3JVJKp/fv3L35rdMzyjKfh8AcC1W/RQdGOG7wN+O+tdA3wSmAjsA/49GzTAYvXAuovLFZdW1WTVTU5MTFxFL2XJM1lGNdRnA98t6qeBph9BkjyBeCb7eUMsK5vubXAU216UP0nwElJjm97Ff3tJUnLZBhDT5fQN+yU5PS+eW8HHm7Tu4CLk7w4yRnABuA7wL3AhnaG0wn0hrF2VVUBdwLvaMtvBW4ZQn8lSUdhUXsUSf4B8BbgPX3l/5hkI71hor2z86rqkSQ3A98HDgKXV9Wv23quAG4HjgN2VNUjbV0fAm5K8kngfuC6xfRXknT00vvivnJMTk7W1NTUqLuhZeaB16Wzd/uFo+6ClkGS+6pqctA8r8yWJHUyKCRJnQwKSVIng0KS1Mmg0DHPA9lLy6u0ZVBIkjoZFJKkTgaFJKmTQSFJ6jSMHwWURsIDrNLycI9C0rwYzKuXQSFJ6mRQSJI6GRSS5s2L71Yng0LHJD+spOVjUEiSOi06KJLsTfJQkgeSTLXay5LsTrKnPZ/c6klydZLpJA8meV3fera29nuSbO2rn9XWP92WzWL7LEmav2HtUfyLqtrYd3ekK4E7qmoDcEd7DXA+vXtlbwC2AddAL1iAq4DXA2cDV82GS2uzrW+5zUPqsyRpHpZq6GkLsLNN7wQu6qvfUD13AyclOR04D9hdVQeq6llgN7C5zTuxqu6q3j1bb+hblyRpGQwjKAr4VpL7kmxrtdOqah9Aez611dcAT/YtO9NqXfWZAfVDJNmWZCrJ1P79+4ewSRpXnnUzHvw3WF2G8RMeb6iqp5KcCuxO8oOOtoOOL9QC6ocWqq4FrgWYnJx8wXxJ0sIteo+iqp5qz88A36B3jOHpNmxEe36mNZ8B1vUtvhZ4ao762gF1SdIyWVRQJPmHSX5rdho4F3gY2AXMnrm0FbilTe8CLm1nP20Cft6Gpm4Hzk1ycjuIfS5we5v3XJJN7WynS/vWJUlaBosdejoN+EY7Y/V44MtV9T+S3AvcnOQy4EfAO1v7W4ELgGngl8C7AarqQJJPAPe2dh+vqgNt+n3A9cBLgNvaQ9KIzR6n2Lv9whH3REttUUFRVU8A/2xA/afAOQPqBVx+hHXtAHYMqE8Br1lMP7UyeABVGg2vzJYkdTIoJEmdDAqNPa+dGG/+26x8BoUkqZNBIUnqZFBIkjoZFBprjn8fGzyOtLIZFJKkTgaFJKmTQSFpaBx+WpmG8TPj0tD5gSOND/coJEmdDApJQ+UZUCuPQaGx44eMNF4MCklSpwUHRZJ1Se5M8miSR5K8v9U/luTHSR5ojwv6lvlwkukkjyU5r6++udWmk1zZVz8jyT1J9iT5SpITFtpfSdLCLGaP4iDwwap6FbAJuDzJmW3eZ6pqY3vcCtDmXQy8GtgMfC7JcUmOAz4LnA+cCVzSt55PtXVtAJ4FLltEfzXmHNteWfy3XDkWHBRVta+qvtumnwMeBdZ0LLIFuKmqnq+qH9K7HerZ7TFdVU9U1a+Am4At7R7Zbwa+2pbfCVy00P5KkhZmKMcokqwHXgvc00pXJHkwyY4kJ7faGuDJvsVmWu1I9ZcDP6uqg4fVB73/tiRTSab2798/hC2SNAzuJa4Miw6KJC8FvgZ8oKp+AVwDvBLYCOwDPj3bdMDitYD6C4tV11bVZFVNTkxMHOUWaBz4YSKNr0VdmZ3kRfRC4ktV9XWAqnq6b/4XgG+2lzPAur7F1wJPtelB9Z8AJyU5vu1V9LeXJC2TxZz1FOA64NGq+vO++ul9zd4OPNymdwEXJ3lxkjOADcB3gHuBDe0MpxPoHfDeVVUF3Am8oy2/Fbhlof3VeHJoYnXw3/jYtpg9ijcA7wIeSvJAq32E3llLG+kNE+0F3gNQVY8kuRn4Pr0zpi6vql8DJLkCuB04DthRVY+09X0IuCnJJ4H76QWTJGkZpffFfeWYnJysqampUXdD8+C3zNVp7/YLR90FDZDkvqqaHDTPK7MlSZ0MCknLyj3JY4/3o9Cy84NCOra4R6FlZUgIPNvtWGNQSJI6GRRaFn6D1CD+TRwbDApJUieDQkvOb43q4t7m+POsJy0Z//PraMz+vXhB3vhxj0KS1Mmg0NA5lKDF8G9n/BgUGir/k2sY/LIxXgwKDYX/sbUU/JsaDx7M1oL5n1jLwYPco2dQaEEMCS23/r85Q2N5GRSaN8NB48LQWF5jHxRJNgP/md7d775YVdtH3KVVxXDQuDM0lt5YB0WS44DPAm8BZoB7k+yqqu+Ptmcrj4GgleBIf8cGyOKMdVAAZwPTVfUEQJKbgC307rutOfjhL/Uc7f8Fg+VQ4x4Ua4An+17PAK8/vFGSbcC29vL/JnlsGfo2CqcAPxl1J5bYathGWB3becxuYz51VM2P2e08zD8+0oxxD4oMqNULClXXAtcufXdGK8nUkW5+vlKshm2E1bGdq2EbYXVs57hfcDcDrOt7vRZ4akR9kaRVadyD4l5gQ5IzkpwAXAzsGnGfJGlVGeuhp6o6mOQK4HZ6p8fuqKpHRtytUVrxw2usjm2E1bGdq2EbYRVsZ6peMOQvSdJvjPvQkyRpxAwKSVIng+IYleRPklSSU0bdl2FL8mdJfpDkwSTfSHLSqPs0LEk2J3ksyXSSK0fdn6WQZF2SO5M8muSRJO8fdZ+WSpLjktyf5Juj7stSMiiOQUnW0ftZkx+Nui9LZDfwmqr6feCvgA+PuD9D0feTNOcDZwKXJDlztL1aEgeBD1bVq4BNwOUrdDsB3g88OupOLDWD4tj0GeA/MODiw5Wgqr5VVQfby7vpXT+zEvzmJ2mq6lfA7E/SrChVta+qvtumn6P3QbpmtL0aviRrgQuBL466L0vNoDjGJHkb8OOq+t6o+7JM/gi4bdSdGJJBP0mz4j5A+yVZD7wWuGe0PVkS/4neF7a/HXVHltpYX0exWiX5n8A/GjDro8BHgHOXt0fD17WNVXVLa/NResMYX1rOvi2hef0kzUqR5KXA14APVNUvRt2fYUryVuCZqrovyZtG3Z+lZlCMoar6w0H1JP8UOAP4XhLoDcl8N8nZVfV/lrGLi3akbZyVZCvwVuCcWjkX+6yan6RJ8iJ6IfGlqvr6qPuzBN4AvC3JBcDfB05M8t+q6t+OuF9LwgvujmFJ9gKTVbUSfrnyN9rNqv4c+OdVtX/U/RmWJMfTOzh/DvBjej9R829W2q8NpPctZidwoKo+MOr+LLW2R/EnVfXWUfdlqXiMQuPovwC/BexO8kCSz4+6Q8PQDtDP/iTNo8DNKy0kmjcA7wLe3P79HmjfvHWMco9CktTJPQpJUieDQpLUyaCQJHUyKCRJnQwKSVIng0KS1MmgkCR1+v8+cKrmAkLtRQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Note: the three output values must be assigned to variables, even though we do not need them here.\n",
    "(n, bin_edges, patches) = plt.hist(normf_samples, 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random integers: `numpy.random.randint`\n",
    "\n",
    "One can generate pseudo-random integers, uniformly distributed between specified lower and upper values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "60 random dice rolls:\n",
      " [4 3 2 3 3 3 3 4 1 6 5 4 4 6 6 3 2 6 5 2 5 4 3 1 4 4 1 1 2 2 3 4 3 4 1 6 1\n",
      " 3 6 6 6 5 6 6 2 5 6 4 1 3 2 5 4 5 2 3 6 3 3 5]\n",
      "1 occured 7 times\n",
      "2 occured 8 times\n",
      "3 occured 14 times\n",
      "4 occured 11 times\n",
      "5 occured 8 times\n",
      "6 occured 12 times\n"
     ]
    }
   ],
   "source": [
    "n_dice = 60\n",
    "dice_rolls = npr.randint(1, 6+1, n_dice)\n",
    "print(n_dice, 'random dice rolls:\\n', dice_rolls)\n",
    "# Count each outcome: this needs a list instead of an array:\n",
    "dice_rolls_list = list(dice_rolls)\n",
    "for value in (1, 2, 3, 4, 5, 6):\n",
    "    count = dice_rolls_list.count(value)\n",
    "    print(value, 'occured', count, 'times')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying bin edges for the histogram\n",
    "\n",
    "This time, it is best to explicitly specify a list of the edges of the bins, by making the second argument `bins` a list.\n",
    "\n",
    "With six values, seven edges are needed, and it looks nicest if they are centered on the integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_edges = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAMgklEQVR4nO3db4hlhXnH8e8vriHRKKbsJLWu20lKEIKUKkPaVLBFY9lW0bzoCwWDTYV506amf7BrA5W+s7SkKbSkLGq0xGqLf2iIIVWSiA0Y29lV6581TbBWN5ruiJTEtGBtnr7YW9iM694795zd6zN+P7DM3HPP3PMc2P1y9txz7qSqkCT187ZFDyBJmo8Bl6SmDLgkNWXAJakpAy5JTW07nhvbvn17LS8vH89NSlJ7e/fufamqljYuP64BX15eZm1t7XhuUpLaS/LvR1ruKRRJasqAS1JTBlySmjLgktSUAZekpgy4JDU1NeBJbk5yMMkTR3ju95JUku3HZjxJ0huZ5Qj8FmDXxoVJzgQuAp4beSZJ0gymBryqHgRePsJTfwZcC/iB4pK0AHPdiZnkUuA7VfVYkmnrrgKrADt37pxnc3oTW95976JHGMWzN1y86BGkTdv0m5hJTgI+BfzhLOtX1Z6qWqmqlaWl193KL0ma0zxXofwU8D7gsSTPAjuAfUl+fMzBJElHt+lTKFX1OPCe/388ifhKVb004lySpClmuYzwduAh4KwkB5JcfezHkiRNM/UIvKqumPL88mjTSJJm5p2YktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqalZfiv9zUkOJnnisGV/kuTpJP+S5J4kpx3bMSVJG81yBH4LsGvDsvuBs6vqp4F/Ba4beS5J0hRTA15VDwIvb1h2X1W9Nnn4DWDHMZhNknQU20Z4jV8H/vaNnkyyCqwC7Ny5c4TNSXqrWN5976JHGM2zN1w8+msOehMzyaeA14Db3midqtpTVStVtbK0tDRkc5Kkw8x9BJ7kKuAS4MKqqvFGkiTNYq6AJ9kF/D7wC1X1X+OOJEmaxSyXEd4OPAScleRAkquBvwBOAe5P8miSvzrGc0qSNph6BF5VVxxh8U3HYBZJ0iZ4J6YkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1NcbHyUrt+bGl6sgjcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpqaGvAkNyc5mOSJw5b9WJL7k3xr8vXdx3ZMSdJGsxyB3wLs2rBsN/CVqvoA8JXJY0nScTQ14FX1IPDyhsWXAbdOvr8V+OjIc0mSppj3HPh7q+pFgMnX97zRiklWk6wlWVtfX59zc5KkjY75m5hVtaeqVqpqZWlp6VhvTpLeMuYN+H8kOR1g8vXgeCNJkmYxb8C/AFw1+f4q4O/HGUeSNKtZLiO8HXgIOCvJgSRXAzcAFyX5FnDR5LEk6Tia+jsxq+qKN3jqwpFnkSRtgndiSlJTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYGBTzJbyd5MskTSW5P8o6xBpMkHd3cAU9yBvBbwEpVnQ2cAFw+1mCSpKMbegplG/DOJNuAk4AXho8kSZrFtnl/sKq+k+RPgeeA/wbuq6r7Nq6XZBVYBdi5c+e8m9tSlnffu+gRtIX59+utY8gplHcDlwHvA34CODnJlRvXq6o9VbVSVStLS0vzTypJ+hFDTqF8BPi3qlqvqv8B7gZ+fpyxJEnTDAn4c8DPJTkpSYALgf3jjCVJmmbugFfVw8CdwD7g8clr7RlpLknSFHO/iQlQVdcD1480iyRpE7wTU5KaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoadCfm8eRHZErSj/IIXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYGBTzJaUnuTPJ0kv1JPjzWYJKkoxv6YVZ/Dny5qn41yduBk0aYSZI0g7kDnuRU4Hzg1wCq6lXg1XHGkiRNM+QUyvuBdeBzSR5JcmOSkzeulGQ1yVqStfX19QGbkyQdbkjAtwHnAp+tqnOAHwC7N65UVXuqaqWqVpaWlgZsTpJ0uCEBPwAcqKqHJ4/v5FDQJUnHwdwBr6rvAs8nOWuy6ELgqVGmkiRNNfQqlE8At02uQHkG+PjwkSRJsxgU8Kp6FFgZaRZJ0iZ4J6YkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqanDAk5yQ5JEkXxxjIEnSbMY4Ar8G2D/C60iSNmFQwJPsAC4GbhxnHEnSrIYegX8GuBb44QizSJI2Ye6AJ7kEOFhVe6est5pkLcna+vr6vJuTJG0w5Aj8PODSJM8CdwAXJPn8xpWqak9VrVTVytLS0oDNSZION3fAq+q6qtpRVcvA5cBXq+rK0SaTJB2V14FLUlPbxniRqnoAeGCM15IkzcYjcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqSkDLklNGXBJasqAS1JTBlySmjLgktSUAZekpuYOeJIzk3wtyf4kTya5ZszBJElHt23Az74G/G5V7UtyCrA3yf1V9dRIs0mSjmLuI/CqerGq9k2+/z6wHzhjrMEkSUc3yjnwJMvAOcDDR3huNclakrX19fUxNidJYoSAJ3kXcBfwyar63sbnq2pPVa1U1crS0tLQzUmSJgYFPMmJHIr3bVV19zgjSZJmMeQqlAA3Afur6tPjjSRJmsWQI/DzgI8BFyR5dPLnV0aaS5I0xdyXEVbV14GMOIskaRO8E1OSmjLgktSUAZekpgy4JDVlwCWpKQMuSU0ZcElqyoBLUlMGXJKaMuCS1JQBl6SmDLgkNWXAJakpAy5JTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ1ZcAlqalBAU+yK8k3k3w7ye6xhpIkTTd3wJOcAPwl8MvAB4ErknxwrMEkSUc35Aj8Q8C3q+qZqnoVuAO4bJyxJEnTbBvws2cAzx/2+ADwsxtXSrIKrE4evpLkmwO2eTxsB15a9BAj2Cr7Ae7Lm9FW2Q84TvuSPx704z95pIVDAp4jLKvXLajaA+wZsJ3jKslaVa0seo6htsp+gPvyZrRV9gN678uQUygHgDMPe7wDeGHYOJKkWQ0J+D8DH0jyviRvBy4HvjDOWJKkaeY+hVJVryX5TeAfgBOAm6vqydEmW5w2p3um2Cr7Ae7Lm9FW2Q9ovC+pet1pa0lSA96JKUlNGXBJasqAA0luTnIwyROLnmWoJGcm+VqS/UmeTHLNomeaV5J3JPmnJI9N9uWPFj3TEElOSPJIki8uepYhkjyb5PEkjyZZW/Q8QyQ5LcmdSZ6e/Jv58KJn2gzPgQNJzgdeAf66qs5e9DxDJDkdOL2q9iU5BdgLfLSqnlrwaJuWJMDJVfVKkhOBrwPXVNU3FjzaXJL8DrACnFpVlyx6nnkleRZYqar2N/IkuRX4x6q6cXI13UlV9Z+LnmtWHoEDVfUg8PKi5xhDVb1YVfsm338f2M+hu2bbqUNemTw8cfKn5RFHkh3AxcCNi55FhyQ5FTgfuAmgql7tFG8w4FtakmXgHODhxU4yv8lph0eBg8D9VdV1Xz4DXAv8cNGDjKCA+5LsnXxURlfvB9aBz01Obd2Y5ORFD7UZBnyLSvIu4C7gk1X1vUXPM6+q+t+q+hkO3en7oSTtTnEluQQ4WFV7Fz3LSM6rqnM59EmkvzE5BdnRNuBc4LNVdQ7wA6DVx2Ib8C1ocr74LuC2qrp70fOMYfJf2weAXQseZR7nAZdOzh3fAVyQ5POLHWl+VfXC5OtB4B4OfTJpRweAA4f9r+5ODgW9DQO+xUze+LsJ2F9Vn170PEMkWUpy2uT7dwIfAZ5e7FSbV1XXVdWOqlrm0EdOfLWqrlzwWHNJcvLkzXEmpxt+CWh59VZVfRd4PslZk0UXAq3e7B/yaYRbRpLbgV8Etic5AFxfVTctdqq5nQd8DHh8cu4Y4A+q6ksLnGlepwO3Tn55yNuAv6uq1pfgbQHvBe45dJzANuBvqurLix1pkE8At02uQHkG+PiC59kULyOUpKY8hSJJTRlwSWrKgEtSUwZckpoy4JLUlAGXpKYMuCQ19X8+h8ksDQTcsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(n, bin_edges, patches) = plt.hist(dice_rolls, bin_edges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the above several times, redrawing the histogram each time; you should see a lot of variation.\n",
    "\n",
    "Things average out with more rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 occured a fraction 0.166567 of the time\n",
      "2 occured a fraction 0.166454 of the time\n",
      "3 occured a fraction 0.166568 of the time\n",
      "4 occured a fraction 0.16633 of the time\n",
      "5 occured a fraction 0.167526 of the time\n",
      "6 occured a fraction 0.166555 of the time\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWkUlEQVR4nO3df4xdZ53f8fdnbcKvbXBCBpraVu0tXloTbUtwg7eoiOLdxAGE8weRHO1uLGrJKg2U7S9IdqVGBSKFdrXZRoJIbuyNQ9OYyMDG2jXrtZJQuhL5MSFAcELqqaHxbAIe1k42LILU8O0f93H3Mr7H47nXmWub90u6uud8n+ec+xzJns+cc547J1WFJEmD/MK4ByBJOnMZEpKkToaEJKmTISFJ6mRISJI6LR73AE63iy66qFasWDHuYUjSWeXRRx/9flVNzK6fcyGxYsUKJicnxz0MSTqrJPk/g+pebpIkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1Oue+cS3ppbfi+j8Z9xBOm+/c/O5xD+GMZkj0OZf+4Us6NefS//uXIvC83CRJ6mRISJI6GRKSpE6GhCSp05whkWR7ksNJvjmr/qEkTyXZn+Q/9dVvSDLV2q7oq69vtakk1/fVVyZ5KMmBJJ9Ncl6rv7ytT7X2FafjgCVJp+5UziTuANb3F5L8M2AD8CtV9Sbg91p9NbAReFPb5tNJFiVZBHwKuBJYDVzT+gJ8ErilqlYBR4HNrb4ZOFpVbwBuaf0kSQtozpCoqi8DR2aVPwDcXFU/bn0Ot/oGYGdV/biqvg1MAZe111RVHayqF4GdwIYkAd4J7Grb7wCu6tvXjra8C1jX+kuSFsiw9yR+Gfin7TLQ/0jyj1t9KXCor990q3XVXws8V1XHZtV/Zl+t/fnW/wRJtiSZTDI5MzMz5CFJkmYbNiQWAxcAa4F/D9zTfssf9Jt+DVFnjrafLVZtrao1VbVmYuKE53hLkoY0bEhMA5+vnoeBnwIXtfryvn7LgGdOUv8+sCTJ4ll1+rdp7a/hxMtekqSX0LAh8Uf07iWQ5JeB8+j9wN8NbGwzk1YCq4CHgUeAVW0m03n0bm7vrqoCHgDe1/a7Cbi3Le9u67T2+1t/SdICmfNvNyW5G3gHcFGSaeBGYDuwvU2LfRHY1H6A709yD/AEcAy4rqp+0vbzQWAvsAjYXlX720d8FNiZ5BPAY8C2Vt8GfCbJFL0ziI2n4XglSfMwZ0hU1TUdTb/Z0f8m4KYB9T3AngH1g/RmP82u/wi4eq7xSZJeOn7jWpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVKnOUMiyfYkh9tT6Ga3/bskleSitp4ktyaZSvKNJJf29d2U5EB7beqrvyXJ422bW5Ok1S9Msq/135fkgtNzyJKkU3UqZxJ3AOtnF5MsB34deLqvfCW951qvArYAt7W+F9J77Olb6T2F7sa+H/q3tb7Htzv+WdcD91XVKuC+ti5JWkBzhkRVfZneM6ZnuwX4CFB9tQ3AndXzILAkycXAFcC+qjpSVUeBfcD61nZ+VX2lPSP7TuCqvn3taMs7+uqSpAUy1D2JJO8F/qKqvj6raSlwqG99utVOVp8eUAd4fVU9C9DeX3eS8WxJMplkcmZmZogjkiQNMu+QSPIq4HeB/zCoeUCthqjPS1Vtrao1VbVmYmJivptLkjoMcybx94CVwNeTfAdYBnw1yd+mdyawvK/vMuCZOerLBtQBvtcuR9HeDw8xVknSCOYdElX1eFW9rqpWVNUKej/oL62q7wK7gWvbLKe1wPPtUtFe4PIkF7Qb1pcDe1vbC0nWtllN1wL3to/aDRyfBbWpry5JWiCnMgX2buArwBuTTCfZfJLue4CDwBTwX4F/CVBVR4CPA4+018daDeADwO1tm/8NfLHVbwZ+PckBerOobp7foUmSRrV4rg5Vdc0c7Sv6lgu4rqPfdmD7gPokcMmA+l8C6+YanyTppeM3riVJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1OlUHjq0PcnhJN/sq/3nJN9K8o0kX0iypK/thiRTSZ5KckVffX2rTSW5vq++MslDSQ4k+WyS81r95W19qrWvOF0HLUk6NadyJnEHsH5WbR9wSVX9CvC/gBsAkqwGNgJvatt8OsmiJIuATwFXAquBa1pfgE8Ct1TVKuAocPzJd5uBo1X1BuCW1k+StIDmDImq+jJwZFbtz6rqWFt9EFjWljcAO6vqx1X1bXqPJL2svaaq6mBVvQjsBDa051q/E9jVtt8BXNW3rx1teRewrvWXJC2Q03FP4p/zN8+lXgoc6mubbrWu+muB5/oC53j9Z/bV2p9v/U+QZEuSySSTMzMzIx+QJKlnpJBI8rvAMeCu46UB3WqI+sn2dWKxamtVramqNRMTEycftCTplC0edsMkm4D3AOuq6vgP72lgeV+3ZcAzbXlQ/fvAkiSL29lCf//j+5pOshh4DbMue0mSXlpDnUkkWQ98FHhvVf2wr2k3sLHNTFoJrAIeBh4BVrWZTOfRu7m9u4XLA8D72vabgHv79rWpLb8PuL8vjCRJC2DOM4kkdwPvAC5KMg3cSG8208uBfe1e8oNV9S+qan+Se4An6F2Guq6qftL280FgL7AI2F5V+9tHfBTYmeQTwGPAtlbfBnwmyRS9M4iNp+F4JUnzMGdIVNU1A8rbBtSO978JuGlAfQ+wZ0D9IL3ZT7PrPwKunmt8kqSXjt+4liR1MiQkSZ0MCUlSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktRpzpBIsj3J4STf7KtdmGRfkgPt/YJWT5Jbk0wl+UaSS/u22dT6H2jPxz5ef0uSx9s2t6Y96q7rMyRJC+dUziTuANbPql0P3FdVq4D72jrAlfSea70K2ALcBr0f+PQee/pWek+hu7Hvh/5tre/x7dbP8RmSpAUyZ0hU1ZfpPWO63wZgR1veAVzVV7+zeh4EliS5GLgC2FdVR6rqKLAPWN/azq+qr1RVAXfO2tegz5AkLZBh70m8vqqeBWjvr2v1pcChvn7TrXay+vSA+sk+4wRJtiSZTDI5MzMz5CFJkmY73TeuM6BWQ9Tnpaq2VtWaqlozMTEx380lSR2GDYnvtUtFtPfDrT4NLO/rtwx4Zo76sgH1k32GJGmBDBsSu4HjM5Q2Aff21a9ts5zWAs+3S0V7gcuTXNBuWF8O7G1tLyRZ22Y1XTtrX4M+Q5K0QBbP1SHJ3cA7gIuSTNObpXQzcE+SzcDTwNWt+x7gXcAU8EPg/QBVdSTJx4FHWr+PVdXxm+EfoDeD6pXAF9uLk3yGJGmBzBkSVXVNR9O6AX0LuK5jP9uB7QPqk8AlA+p/OegzJEkLx29cS5I6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSeo0Ukgk+ddJ9if5ZpK7k7wiycokDyU5kOSzSc5rfV/e1qda+4q+/dzQ6k8luaKvvr7VppJcP8pYJUnzN3RIJFkK/CtgTVVdAiwCNgKfBG6pqlXAUWBz22QzcLSq3gDc0vqRZHXb7k3AeuDTSRYlWQR8CrgSWA1c0/pKkhbIqJebFgOvTLIYeBXwLPBOYFdr3wFc1ZY3tHVa+7okafWdVfXjqvo2vedjX9ZeU1V1sKpeBHa2vpKkBTJ0SFTVXwC/BzxNLxyeBx4FnquqY63bNLC0LS8FDrVtj7X+r+2vz9qmq36CJFuSTCaZnJmZGfaQJEmzjHK56QJ6v9mvBP4O8Gp6l4Zmq+ObdLTNt35isWprVa2pqjUTExNzDV2SdIpGudz0a8C3q2qmqv4v8HngnwBL2uUngGXAM215GlgO0NpfAxzpr8/apqsuSVogo4TE08DaJK9q9xbWAU8ADwDva302Afe25d1tndZ+f1VVq29ss59WAquAh4FHgFVtttR59G5u7x5hvJKkeVo8d5fBquqhJLuArwLHgMeArcCfADuTfKLVtrVNtgGfSTJF7wxiY9vP/iT30AuYY8B1VfUTgCQfBPbSmzm1var2DzteSdL8DR0SAFV1I3DjrPJBejOTZvf9EXB1x35uAm4aUN8D7BlljJKk4fmNa0lSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktTJkJAkdRopJJIsSbIrybeSPJnkV5NcmGRfkgPt/YLWN0luTTKV5BtJLu3bz6bW/0CSTX31tyR5vG1za3sCniRpgYx6JvFfgD+tqr8P/EPgSeB64L6qWgXc19YBrqT3aNJVwBbgNoAkF9J7cNFb6T2s6MbjwdL6bOnbbv2I45UkzcPQIZHkfODttMeTVtWLVfUcsAHY0brtAK5qyxuAO6vnQWBJkouBK4B9VXWkqo4C+4D1re38qvpKexb2nX37kiQtgFHOJH4JmAH+MMljSW5P8mrg9VX1LEB7f13rvxQ41Lf9dKudrD49oC5JWiCjhMRi4FLgtqp6M/DX/M2lpUEG3U+oIeon7jjZkmQyyeTMzMzJRy1JOmWjhMQ0MF1VD7X1XfRC43vtUhHt/XBf/+V92y8DnpmjvmxA/QRVtbWq1lTVmomJiREOSZLUb+iQqKrvAoeSvLGV1gFPALuB4zOUNgH3tuXdwLVtltNa4Pl2OWovcHmSC9oN68uBva3thSRr26yma/v2JUlaAItH3P5DwF1JzgMOAu+nFzz3JNkMPA1c3fruAd4FTAE/bH2pqiNJPg480vp9rKqOtOUPAHcArwS+2F6SpAUyUkhU1deANQOa1g3oW8B1HfvZDmwfUJ8ELhlljJKk4fmNa0lSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktTJkJAkdTIkJEmdRg6JJIuSPJbkj9v6yiQPJTmQ5LPtqXUkeXlbn2rtK/r2cUOrP5Xkir76+labSnL9qGOVJM3P6TiT+DDwZN/6J4FbqmoVcBTY3OqbgaNV9QbgltaPJKuBjcCbgPXAp1vwLAI+BVwJrAauaX0lSQtkpJBIsgx4N3B7Ww/wTmBX67IDuKotb2jrtPZ1rf8GYGdV/biqvk3vGdiXtddUVR2sqheBna2vJGmBjHom8QfAR4CftvXXAs9V1bG2Pg0sbctLgUMArf351v//12dt01U/QZItSSaTTM7MzIx4SJKk44YOiSTvAQ5X1aP95QFda462+dZPLFZtrao1VbVmYmLiJKOWJM3H4hG2fRvw3iTvAl4BnE/vzGJJksXtbGEZ8EzrPw0sB6aTLAZeAxzpqx/Xv01XXZK0AIY+k6iqG6pqWVWtoHfj+f6q+g3gAeB9rdsm4N62vLut09rvr6pq9Y1t9tNKYBXwMPAIsKrNljqvfcbuYccrSZq/Uc4kunwU2JnkE8BjwLZW3wZ8JskUvTOIjQBVtT/JPcATwDHguqr6CUCSDwJ7gUXA9qra/xKMV5LU4bSERFV9CfhSWz5Ib2bS7D4/Aq7u2P4m4KYB9T3AntMxRknS/PmNa0lSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktTJkJAkdTIkJEmdhg6JJMuTPJDkyST7k3y41S9Msi/JgfZ+Qasnya1JppJ8I8mlffva1PofSLKpr/6WJI+3bW5NklEOVpI0P6OcSRwD/m1V/QNgLXBdktXA9cB9VbUKuK+tA1xJ7/nVq4AtwG3QCxXgRuCt9J5od+PxYGl9tvRtt36E8UqS5mnokKiqZ6vqq235BeBJYCmwAdjRuu0ArmrLG4A7q+dBYEmSi4ErgH1VdaSqjgL7gPWt7fyq+kpVFXBn374kSQvgtNyTSLICeDPwEPD6qnoWekECvK51Wwoc6ttsutVOVp8eUB/0+VuSTCaZnJmZGfVwJEnNyCGR5BeBzwG/XVV/dbKuA2o1RP3EYtXWqlpTVWsmJibmGrIk6RSNFBJJXkYvIO6qqs+38vfapSLa++FWnwaW922+DHhmjvqyAXVJ0gIZZXZTgG3Ak1X1+31Nu4HjM5Q2Aff21a9ts5zWAs+3y1F7gcuTXNBuWF8O7G1tLyRZ2z7r2r59SZIWwOIRtn0b8FvA40m+1mq/A9wM3JNkM/A0cHVr2wO8C5gCfgi8H6CqjiT5OPBI6/exqjrSlj8A3AG8Evhie0mSFsjQIVFVf87g+wYA6wb0L+C6jn1tB7YPqE8Clww7RknSaPzGtSSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROZ3xIJFmf5KkkU0muH/d4JOnnyRkdEkkWAZ8CrgRWA9ckWT3eUUnSz48zOiSAy4CpqjpYVS8CO4ENYx6TJP3cGPoZ1wtkKXCob30aeOvsTkm2AFva6g+SPLUAYxvFRcD3xz2I0+BcOQ7wWM5E58pxwAIdSz450uZ/d1DxTA+JDKjVCYWqrcDWl344p0eSyapaM+5xjOpcOQ7wWM5E58pxwNl9LGf65aZpYHnf+jLgmTGNRZJ+7pzpIfEIsCrJyiTnARuB3WMekyT93DijLzdV1bEkHwT2AouA7VW1f8zDOh3OmktjczhXjgM8ljPRuXIccBYfS6pOuMQvSRJw5l9ukiSNkSEhSepkSCyQJNuTHE7yzXGPZVRJlid5IMmTSfYn+fC4xzSsJK9I8nCSr7dj+Y/jHtMokixK8liSPx73WEaR5DtJHk/ytSST4x7PKJIsSbIrybfa/5lfHfeY5sN7EgskyduBHwB3VtUl4x7PKJJcDFxcVV9N8reAR4GrquqJMQ9t3pIEeHVV/SDJy4A/Bz5cVQ+OeWhDSfJvgDXA+VX1nnGPZ1hJvgOsqaqz/st0SXYA/7Oqbm+zNF9VVc+Ne1ynyjOJBVJVXwaOjHscp0NVPVtVX23LLwBP0vt2/Fmnen7QVl/WXmflb05JlgHvBm4f91jUk+R84O3ANoCqevFsCggwJDSiJCuANwMPjXckw2uXaL4GHAb2VdXZeix/AHwE+Om4B3IaFPBnSR5tf3bnbPVLwAzwh+0y4O1JXj3uQc2HIaGhJflF4HPAb1fVX417PMOqqp9U1T+i943+y5KcdZcDk7wHOFxVj457LKfJ26rqUnp/Afq6drn2bLQYuBS4rareDPw1cFY98sCQ0FDa9fvPAXdV1efHPZ7ToV0G+BKwfsxDGcbbgPe2a/k7gXcm+W/jHdLwquqZ9n4Y+AK9vwh9NpoGpvvOTnfRC42zhiGheWs3e7cBT1bV7497PKNIMpFkSVt+JfBrwLfGO6r5q6obqmpZVa2g9+dr7q+q3xzzsIaS5NVtQgTt0szlwFk5K7CqvgscSvLGVloHnFUTPM7oP8txLklyN/AO4KIk08CNVbVtvKMa2tuA3wIeb9fyAX6nqvaMcUzDuhjY0R5w9QvAPVV1Vk8fPQe8HvhC73cRFgP/var+dLxDGsmHgLvazKaDwPvHPJ55cQqsJKmTl5skSZ0MCUlSJ0NCktTJkJAkdTIkJEmdDAlJUidDQpLU6f8BgS0Tye9D4jIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_dice = 10**6\n",
    "dice_rolls = npr.randint(1, 6+1, n_dice)\n",
    "# Count each outcome: this needs a list instead of an array:\n",
    "dice_rolls_list = list(dice_rolls)\n",
    "for value in (1, 2, 3, 4, 5, 6):\n",
    "    count = dice_rolls_list.count(value)\n",
    "    print(value, 'occured a fraction', count/n_dice, 'of the time')\n",
    "(n, bin_edges, patches) = plt.hist(dice_rolls, bins = bin_edges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram is now visually boring, but mathematically more satisfying."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"Exercise-A\"></a>\n",
    "#### Exercise A: approximating $\\pi$\n",
    "\n",
    "We can compute approximations of $\\pi$ by using the fact that the unit circle occupies a fraction $\\pi/4$ of the circumscribed square:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc7ea44be10>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAHSCAYAAAAXPUnmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5jTVdrG8e8zIEWUJlUEAUXsdcS1K4oFewXLK1ZWxd67K4iKde3KqqhrAbuoWFDBCupgBUREREUUEcGCdM77x5MscZxhSjI5Seb+XFeuTH5JJk/IMPec8zvFQgiIiIhIfiqKXYCIiIhUn4JcREQkjynIRURE8piCXEREJI8pyEVERPKYglxERCSP1Y1dQHW0aNEidOzYMXYZIiIiWTFu3LifQwgty7ovL4O8Y8eOlJSUxC5DREQkK8zsm/LuU9e6iIhIHlOQi4iI5DEFuYiISB5TkIuIiOQxBbmIiEgeU5CLiIjkMQW5iIhIHlOQi4iI5DEFuYiISB5TkIuIiOQxBbmIiEgeU5CLiIjkMQW5iIhIHlOQi4iI5LGMBLmZ3WdmP5nZ+HLuNzO7xcymmNmnZrZ5yn19zOzLxKVPJuoRERGpLTLVIr8f2GMF9+8JdElc+gJ3AphZc+ByYCugG3C5mTXLUE0iIiIFLyNBHkJ4E/hlBQ/ZD3gwuLFAUzNrC+wOjAwh/BJCmAOMZMV/EIiIiEiKull6nXbAdym3pyeOlXc8azbbDGbNgrXXzuarimRBWAYLFsKypbBsWQ1fAqy0EtSrV/GlSENzpLBNmQItW8JHH2Xn9bIV5FbGsbCC43//BmZ98W55OnTokLHCZs2CP/7I2LcTyb6lS+DPP5df5iWu58+nnP9O5TPzoC3vUqeOB3bp42aweDEsWgQL5sNvv/rtstSpW7nAX2kl/74ieSbbmZKtIJ8OtE+5vQYwI3F8p1LHR5f1DUIIg4HBAMXFxVX87VS+ZEt8dJmvKpIjQoDvv4dJk5ZfPv/cr2fMWP64unWhSxdYd12/rLMONGkCDRpUfKlfP7Ot5cWL/S/lH38s+/LDD8u/Lus3X1ERtGoFbdosv6yxBhQXw9Zb+30iOWinnbL7etkK8uHAKWY2FB/Y9msI4Qczexm4KmWA227AhVmqSST3LFrk/XKlw3rSpL+GXePGsN560KOHXyeDu3Nnb8nmgpVWgtVX90tF/vgDZs4sP/R//BHGj/fwX7rUn9O5swd68rLxxv6HjEgtk5GfejN7FG9ZtzCz6fhI9JUAQgh3ASOAnsAU4E/gmMR9v5jZAOCDxLfqH0JY0aA5kcLw++8wYcLfA/urr5YHFUD79h7Qxx67PKzXXddbp4XU7bzKKn5Za60VP27+fBg3DsaM8ctrr8HDD/t9K68MW24J//jH8nBXq11qgYwEeQjhsAruD0C/cu67D7gvE3WI5KwQPKxfeAFGjIC334YlS/y+evW8O3yjjeDQQ5eHddeuHm6yXMOGsN12fgH/d/322+XBPmYM3HDD8n9btdqlFtBPtEhNmTcPRo3y4B4xAr75xo9vvDGcc44Hy3rrQadOCpfqMoM11/RL795+rDKt9mSw/+MfarVL3tNvD5FMmjJleXCPHg0LF0KjRn4u++KLYc89fcCW1JzKtNqvv355q32ttf7aHa9Wu+QZ/bSKpGPBAnjzzeXh/eWXfnzddaFfP+jZ0wOlfv24ddZm1Wm177EHHHmkf3767CTHKchFquqbb+DFFz24X3vN52w3aADdu8Ppp3uru3Pn2FXKiqyo1f7WW/Dkk/DUU9C0qY9bOOIIf6wWs5EcpCAXqcjixfDOO8tb3RMm+PGOHeGYY7zVttNO3pKT/FS61X7zzctb6Q8/DIMHQ4cOHuhHHgnrrx+7YpH/UZCLlOWHH5a3ukeOhN9+83nRO+zgU8F69vRR5YU0BUyWq1sXdt/dL3feCc8+Cw89BNdeC1dfDZtu6oF+2GGVmycvUoPUTySStGgRPPoobLON/3I+7jgYOxZ69YKnn4bZs+HVV+Gss/wcuEK8dmjUCA4/3P+o+/57b63Xq+czD9ZYwwcy3n+//7EnEoGCXOSnn+DKK72r/PDD4eefYeBA+OQT+O4771bdf39YddXYlUpsrVvDaafBe+/BF1/ApZfC1Kl+iqV1a++Wf/758teZF6kBCnKpvT76CI4+2ldPu/RSX5DlhRd8hbWLLvJpSGp1S3nWWQeuuMKnHL77rvfgvPoq7LMPtG3rsxbGjPGBdCI1SEEutcuSJfDEE7D99rD55v718cfDxInw8st+7lsjk6UqzHz++W23+diK556DXXeF++7z0zRrrw2XXeYteJEaoN9YUjvMng2DBvm0sEMOgenTfSnP6dPh9tt9hTWRdK20Euy9Nwwd6pvA3H+//8xdeaWPq+jWDW65xU/niGSIglwK22efwQkn+KCkCy7wNc2feca7Q886y+cJi9SExo2hTx+f9fDdd8tXkzv9dJ/mdtZZHvYiaVKQS+FZutSnC3Xv7ue5H3oI/u//4NNPfW7wfvtBnTqxq5TapF07OPts+PBD/+PysMO8Zd6pE5x3nu/bLlJNCnIpHHPnwo03eqt7//291X3NNd59PniwD2YTiW3DDf38+eefw0EHeUu9UycfYDl7duzqJA8pyCX/TZoEJ5+8vNWzxhrw+OM+Lej882G11WJXKPJ3XbrAf//rKwXus4//0dmpkw+Mmzs3dnWSRxTkkp+WLfMFOnbf3Qeq3Xuvr4n94Ye+icnBB2sHK8kP663nCxF9+qn/PA8Y4GsaDBigRWakUhTkkl+WLfORwOuuC3vt5ecbBwzwwURDhsBmm8WuUKR6NtzQe5I+/tjX7r/sMg/0q6+GP/6IXZ3kMAW55I+SEp+Xe8wxPiL4kUdg2jS45BJo1Sp2dSKZsckmPrMi+fN+0UXe5X7ddTBvXuzqJAcpyCX3zZrlU8i6dfPgfuABeP99H/lbr17s6kRqxhZb+HKvY8f61+ed53PSb7rJ91MXSVCQS+5assRXy1pnHe9OP/NMXx3rqKO0+prUHlttBS+9BG+/7TMvzjoL1loLbr0VFiyIXZ3kAP02lNz0xhu+hOqpp0JxsQ8EuuEGaNIkdmUicWy7ra/lPnq0j3g/7TS/vusu37lPai0FueSW6dO9y3ynneDXX30t9Fde0RKqIkk77uhh/uqr0KEDnHSSB/o992jXtVpKQS65YeFCn0e77rq+9/dlly1fMEM7kIn8lRnssot3t7/0ErRp4+NIunb1lQy141qtoiCX+EaM8HN/F14IPXp4gF9xBay8cuzKRHKbmc89HzvWB8Y1a+bLEe+9N8yYEbs6yRIFucTz1Ve+otVee/kvpJde8tZ4p06xKxPJL2b+/+iDD3wN91GjfF76I4+odV4LKMgl++bNg4svhvXX93N9117rC7vsvnvsykTyW1GRDxD9+GPvZj/iCN+2V5uyFDQFuWRPCPDYY34e/KqrfEnVL76Ac8/VfHCRTFpnHT9/fs018NxzsMEGvsiMFCQFuWTH+PE+OKdXL2jRwn/J/Pe/sPrqsSsTKUx16vimQSUlvpHQAQf4GgzakKXgKMilZs2dC2ecAZtuCp98Anfe6b9Ytt02dmUitcNGG/lguMsu83PmG24IL78cuyrJIAW51JyHH/Yuvltu8akxkyfDiSd6S0FEsqdePZ8JMnasL6q0xx7+f/H332NXJhmgIJfMW7DAg/vII2HttWHcOG+Ja19wkbiKi/3/4znnwODBvkHLG2/ErkrSpCCXzJo2DbbbzleZuugieOstbS0qkksaNPCd1N5800e577yzr9+ujVjyloJcMuell3yXpilT4NlnYeBAdaOL5KrttvNxKyef7DuqbbaZ7yooeUdBLulbtszPv/Xs6aNjS0pg331jVyUiFWnUyHcYHDkS/vwTtt4aLrlEm7DkGQW5pOeXX3xFqX/9y5eGHDPGz4uLSP7YdVdflKlPH+9J23JLb61LXlCQS/WNG+dbjb7+ug9mu/9+rY8ukq+aNIH77oPhw2HmTA/zgQNhyZLYlUkFFORSPffc43PBly3zAW0nnqhdykQKwT77wIQJcOCB3s2+zTYwaVLsqmQFFORSNfPnw3HH+fSyHXaADz+Ebt1iVyUimbTaajB0qF+++soHwg0dGrsqKYeCXCpv6lRvhd93n/+l/uKLvtyqiBSmXr28db7llnDYYXD99dpNLQfVjV2A5IkRI3wnJfBNGPbeO249IpIdbdrAK6/4Ou3nngvffuvT1TS1NGdkpEVuZnuY2RdmNsXMLijj/pvM7OPEZbKZzU25b2nKfcMzUY9k0NKlvkbzXnvBmmv6ADeFuEjt0qCBd62feSbceqvvXKgFZHJG2i1yM6sD3A70AKYDH5jZ8BDCxORjQghnpjz+VCB1qa/5IYRN061DasDs2XD44f7X+NFHwx13QMOGsasSkRiKiuDGG6F9ezj7bOjRwxd+0tLL0WWiRd4NmBJCmBpCWAQMBfZbweMPAx7NwOtKTfrgA59aNnq0r8l8330KcRHxVvmwYct3Mfz669gV1XqZCPJ2wHcpt6cnjv2Nma0JdAJeTzncwMxKzGysme2fgXokHSF4cG+3nd9+5x0foa6pZSKSdMghvhrczJm+GtyHH8auqFbLRJCX9Ru+vGGNvYEnQghLU451CCEUA4cD/zaztcp8EbO+icAvmTVrVnoVS9nmz4djj4V//tM3UvjwQ98tSUSktO239z/069f3qagvvRS7olorE0E+HWifcnsNYEY5j+1NqW71EMKMxPVUYDR/PX+e+rjBIYTiEEJxy5Yt061ZSps50xd+uP9+H9z2wgs69yUiK7b++suXZd57bxgyJHZFtVImgvwDoIuZdTKzenhY/230uZl1BZoBY1KONTOz+omvWwDbAhNLP1dq2MyZ0L07TJ4Mzz/vG6BoaomIVMbqq/uWqN27e49e//6aa55laQd5CGEJcArwMvA58FgIYYKZ9Tez1C2wDgOGhvCXT3g9oMTMPgFGAdekjnaXLPjpJ/8POG2azxXfa6/YFYlIvmnc2BsBRx0Fl1/up+e0RnvWZGRBmBDCCGBEqWOXlbr9rzKe9y6wUSZqkGqYNWt5iL/wAuy4Y+yKRCRf1avnp+bat/fNVr7/3ke3r7JK7MoKnpZora2SIT51qv8lvdNOsSsSkXxnBldeCXfd5YPfdt7ZT91JjVKQ10Y//wy77OKbITz/vP9nExHJlH/+E555xtdp32YbH38jNUZBXtskQ/zLL33N9O7dY1ckIoVon318QanffvMwHzs2dkUFS0Fem8yeDbvu6n8dP/ecB7qISE3p1s2npzVt6j1/zz4bu6KCpCCvLZIhPmkSDB/uX4uI1LS114Z334WNN4YDD4Q774xdUcFRkNcGv/ziGxx8/rmHeI8esSsSkdqkVSt4/XXo2RNOPhmuuSZ2RQVFQV7okiE+caJ3a+22W+yKRKQ2atQInn7ad1S88EJ46KHYFRWMjMwjlxw1Z44H9/jxPoJ0991jVyQitVndur6M6w8/+Cpw7dpp1kwGqEVeqObO9RD/7DP/K3jPPWNXJCLiC8c89RR06QIHHOBT1CQtCvJClAzxTz7x/zA9e8auSERkuaZNfUnohg3999MPP8SuKK8pyAvNr796F/rHH3uIa+10EclFa67pS0PPnu2/p/74I3ZFeUtBXkiSIf7RR/Dkk76toIhIrtp8c3j8cfj0Uzj0UG20Uk0K8kLx22+wxx7w4Yf+H2OffWJXJCJSsT33hDvugBdfhH79tAVqNWjUeiH47TdviZeUeIjvt1/sikREKq9vX9+F8eqroVMnuOCC2BXlFQV5vps3z1viJSXw2GOw//6xKxIRqborr4RvvvE55h06+HxzqRQFeT4LwXcZGjvWW+IHHBC7IhGR6ikqgvvu833MjznG55jvuGPsqvKCzpHns7vvhocfhv794aCDYlcjIpKe+vV93Yu11vLexc8/j11RXlCQ56uSEjj9dB8octFFsasREcmMZs18jnn9+v777ccfY1eU8xTk+eiXX+Dgg6F1a/jvf71LSkSkUHTs6HPMZ83SHPNKUALkm2XLoE8fmDHDz4uvtlrsikREMm+LLXwA78cfQ+/emmO+AgryfHPttfD883DjjbDVVrGrERGpOXvtBbff7q3z007THPNyaNR6Phk9Gi6+GHr18oUTREQK3Ykn+hzzQYO8y/2882JXlHMU5Pnihx+8e6lLF/jPf8AsdkUiItlx1VU+x/z8832Oee/esSvKKQryfLBkif/g/v47vPYarLpq7IpERLKnqAjuv9/HBvXp43PMt98+dlU5Q+fI88Ell8Cbb8LgwbDBBrGrERHJvuQc806dfBnqL7+MXVHOUJDnuuHD/dzQiSfCEUfErkZEJJ7mzX1zFTNvmS9dGruinKAgz2VTp8JRR/k0jJtuil2NiEh8nTrBrbfCmDH6vZigIM9VCxb4oi9mPl+8QYPYFYmI5IbDDvMlXC+5BCZNil1NdAryXHX66fDRR/Dgg/4XqIiIODO4805o1AiOPrrWd7EryHPRgw/6wLYLLoB99oldjYhI7mnTBm67Dd57D264IXY1USnIc81nn/nAth13hAEDYlcjIpK7evf27Zsvu6xW75SmIM8lv//u58WbNIGhQ6GupvmLiJQr2cW+yirexV5L12NXkOeKEOD44+Grr2DYMO82EhGRFWvd2tdjf//9WtvFriDPFbfd5jv9XHUV7LBD7GpERPLHoYfCQQd5F/vEibGryToFeS4YOxbOPhv23RfOOSd2NSIi+cUM7rgDGjeulV3sCvLYFi6EI4+ENdbwtYSL9JGIiFRZq1bexf7BB3D99bGrySqlRmz//refF7/7bmjWLHY1IiL569BDfcDw5ZfDhAmxq8kaBXlMP/4IV17pXeo9esSuRkQk/91+e63rYleQx3TRRd61XktHWoqIZFyrVn6+vKQErrsudjVZoSCP5YMPYMgQOPNMWHvt2NWIiBSOQw7xy7/+BePHx66mxinIYwjB11Jv3Rouvjh2NSIihef2231xrWOOKfgu9owEuZntYWZfmNkUM7ugjPuPNrNZZvZx4nJ8yn19zOzLxKVPJurJeY8+6lvwXX21n8sREZHMatlyeRf7tdfGrqZGpR3kZlYHuB3YE1gfOMzM1i/jocNCCJsmLvckntscuBzYCugGXG5mhT10e948OO8832O8T+34u0VEJIqDD/aR7AXexZ6JFnk3YEoIYWoIYREwFNivks/dHRgZQvglhDAHGAnskYGactegQfD993DLLZozLiJS0267DZo29VHsixfHrqZGZCJJ2gHfpdyenjhW2kFm9qmZPWFm7av43MIwbZqPojz8cNhmm9jViIgUvpYtfWOVceMKtos9E0FuZRwLpW4/B3QMIWwMvAo8UIXn+gPN+ppZiZmVzJo1q9rFRnXeed4KHzQodiUiIrXHQQdBr15wxRW+VXSByUSQTwfap9xeA5iR+oAQwuwQwsLEzf8AW1T2uSnfY3AIoTiEUNyyZcsMlJ1lb7wBjz8OF1zgy7GKiEj23Habj2I/66zYlWRcJoL8A6CLmXUys3pAb2B46gPMrG3KzX2B5A7wLwO7mVmzxCC33RLHCsvSpT7drEMHbYoiIhJDixZw/vnw6qvw7ruxq8motIM8hLAEOAUP4M+Bx0IIE8ysv5ntm3jYaWY2wcw+AU4Djk489xdgAP7HwAdA/8SxwnLvvfDJJ76Qf8OGsasREamdTjrJA71//9iVZJSFUOYp6ZxWXFwcSkpKMvK9dtrJr0ePzsi3+7u5c6FLF1h/fX8RK2tYgIiIZMWgQX6Kc+xY2GqrGnmJmsgVMxsXQigu6z7Nf6pp/fvD7Nlw880KcRGR2Pr1g9VWK6hWuYK8Jk2aBLfeCiecAJtuGrsaERFZZRUf8DZihK/6VgAU5DXprLOgUSPfqlRERHLDKadAs2YF0ypXkNeUESPgxRd9g/t8nC4nIlKoGjf2nSefew4++ih2NWlTkNeERYv8h6RrVz8fIyIiueW003zp1gJolSvIa8Jtt8HkyXDTTVCvXuxqRESktCZN4Iwz4JlnfHpwHlOQZ9pPP/kygD17wp57xq5GRETKc9pp3s0+YEDsStKiIM+0q6+GP/+EG2+MXYmIiKxIs2a+6uaTT+b1GuwK8kyaNw+GDPH9b7t2jV2NiIhU5IwzYNVV83p2kYI8kx55BH79VQPcRETyRfPmcOqpvqnVxImxq6kWBXmmhAC33w6bbAJbbx27GhERqawzz4SVV87bVrmCPFPGjPGRj/36aSlWEZF80qKFLxIzdKivyJlnFOSZcscdPvrx8MNjVyIiIlV19tm+O+XAgbErqTIFeSb89JOfXzn6aF+SVURE8kvLlnDyyT7WafLk2NVUiYI8E+6911dzO/nk2JWIiEh1nXMO1K+fd61yBXm6li6Fu+6CXXbRlDMRkXzWujWceCI8/DBMmRK7mkpTkKfrhRfg22/VGhcRKQTnngsrrQRXXRW7kkpTkKfrjjugXTvYd9/YlYiISLratoW+feHBB+Hrr2NXUykK8nR8+SW8/DL8859Qt27sakREJBPOP99/p+dJq1xBno677vIP+4QTYlciIiKZsvrq0KePnyv/44/Y1VRIQV5df/7p66ofdBC0aRO7GhERyaQjjoD58+G552JXUiEFeXUNHQpz5miQm4hIIdpuO2+ZDxsWu5IKKcirI7mu+oYbwvbbx65GREQyrajId7J88UWYOzd2NSukIK+O99+HDz/01rjWVRcRKUy9evliX88+G7uSFVKQV8cdd/j+tUceGbsSERGpKVttBWuumfPd6wryqvr5Z/9QjzrKw1xERAqTmXevjxwJs2fHrqZcCvKquu8+WLgQTjopdiUiIlLTeveGJUvgqadiV1IuBXlVJNdV32kn2GCD2NWIiEhN22wzWHvtnO5eV5BXxUsv+ZJ9mnImIlI7mPmgt1GjYObM2NWUSUFeFXfc4evw7r9/7EpERCRbevWCZcvgySdjV1ImBXllff21zyfs29d3xhERkdphww1hvfV8IbAcpCCvrGee8YVgjjkmdiUiIpJNZj7o7e234fvvY1fzNwryyho5Erp29TmFIiJSu/Tq5Y25xx+PXcnfKMgrY+FCeOMN6NEjdiUiIhJD166wySY5OXpdQV4ZY8b4bmcKchGR2qtXLxg7FqZNi13JXyjIK2PkSKhTx+ePi4hI7dSrl18/9ljcOkpRkFfGyJHwj39A48axKxERkVg6d4Ytt8y57nUFeUV++QVKSmDXXWNXIiIisfXq5btfTpkSu5L/UZBX5PXXfaSizo+LiMihh/p1DrXKFeQVefVV3+WsW7fYlYiISGzt28O22+bU4jAZCXIz28PMvjCzKWZ2QRn3n2VmE83sUzN7zczWTLlvqZl9nLgMz0Q9GTVyJOy8s1ZzExER16sXjB8PEyfGrgTIQJCbWR3gdmBPYH3gMDNbv9TDPgKKQwgbA08A16bcNz+EsGnism+69WTU1Kl+Ube6iIgkHXywr/aWI93rmWiRdwOmhBCmhhAWAUOB/VIfEEIYFUL4M3FzLLBGBl635o0c6dcKchERSWrbFnbc0YM8hNjVZCTI2wHfpdyenjhWnuOAF1NuNzCzEjMba2a5ta3YyJF+PmSddWJXIiIiueTAA+GLL+Dbb2NXQt0MfA8r41iZf6KY2ZFAMbBjyuEOIYQZZtYZeN3MPgshfFXGc/sCfQE6dOiQftUVWbrUR6wfcIB3oYiIiCQlB0CPGxd9D45MtMinA+1Tbq8BzCj9IDPbFbgY2DeEsDB5PIQwI3E9FRgNbFbWi4QQBocQikMIxS1btsxA2RUYNw7mzFG3uoiI/N3GG/uKn+PGxa4kI0H+AdDFzDqZWT2gN/CX0edmthlwNx7iP6Ucb2Zm9RNftwC2BXJjGGDy/Pguu8StQ0REck/DhrDBBoUR5CGEJcApwMvA58BjIYQJZtbfzJKj0K8DVgEeLzXNbD2gxMw+AUYB14QQcifIN9sMstH6FxGR/LPFFh7kkQe8ZeIcOSGEEcCIUscuS/m6zPVNQwjvAhtlooaM+uMPePddOPPM2JWIiEiuKi6GIUN8wFvE8+Ra2a0sb74JixdrfXURESnfFlv4deTudQV5WUaOhPr1YbvtYlciIiK5KkcGvCnIyzJyJGy/vQ9mEBERKUuODHhTkJc2YwZMmKBpZyIiUrEcGPCmIC/ttdf8WkEuIiIVKS6Gn3+G776r+LE1REFe2siRPuVsk01iVyIiIrkuOeCtpCRaCQry0t54w7ctLdI/jYiIVCAHBrwprVItWODzATfcMHYlIiKSD3JgwJuCPNW0aX7duXPUMkREJI9EHvCmIE81dapfK8hFRKSyIg94U5CnUpCLiEhVRV7hTUGeaupUWHllaNUqdiUiIpIvkgPeIo1cV5CnmjoVOnUCs9iViIhIvog84E1Bnurrr9WtLiIiVRdxwJuCPCkEb5EryEVEpKoiDnhTkCf9/LPvQ64gFxGRqoo44E1BnqQR6yIiUl0RV3hTkCcpyEVEpLqSA94ijFxXkCclg7xjx6hliIhInkoOeCO7A94U5ElTp0Lbtj6PXEREpKrWXdfHWy1dltWXVZAnJeeQi4iIVEezZn69ZHFWX1ZBnqQ55CIiko7/BfmSrL6sghwgLPO5fwpyERGpLgV5RAsWwrJlCnIREam+ZJAvVpBn34L5fq0gFxGR6lKLPKL5C/xaQS4iItWlwW4RLZgP9ev79DMREZHqaNzYd89UizyC+Qt86lmR/jlERKSaioqgSROdI49iwXzNIRcRkfQ1a6YWeRTzF+j8uIiIpK9ZM50jz7oli2HpEgW5iIikTy3yCBYkRqxrsxQREUmXgjyCZYldarRZioiIpKtZMw12ExERyVtqkYuIiOSxpk19/45lS7P2kgpyERGRTImwTKuCXEREJFMibJyiIBcREckUtchFRETymIJcREQkj0XYAU1BLiIikin52iI3sz3M7Aszm2JmF5Rxf30zG5a4/z0z65hy34WJ41+Y2e6ZqEdERCSKJk38Op8Gu5lZHeB2YE9gfeAwM1u/1MOOA+aEEL+LHU0AACAASURBVNYGbgIGJZ67PtAb2ADYA7gj8f1ERETyT926UKdO3rXIuwFTQghTQwiLgKHAfqUesx/wQOLrJ4BdzMwSx4eGEBaGEL4GpiS+n4iISH5auhRmzcray2UiyNsB36Xcnp44VuZjQghLgF+B1Sr5XADMrK+ZlZhZyaws/gOJiIhU2aKFWXupTAS5lXEsVPIxlXmuHwxhcAihOIRQ3LJlyyqWKCIikkUtspdTmQjy6UD7lNtrADPKe4yZ1QWaAL9U8rkiIiL5o6gIGjTI3stl4Ht8AHQxs05mVg8fvDa81GOGA30SXx8MvB5CCInjvROj2jsBXYD3M1CTiIhI9i1cCMuWwUp1s/aSab9SCGGJmZ0CvAzUAe4LIUwws/5ASQhhOHAv8F8zm4K3xHsnnjvBzB4DJgJLgH4hhOxtGSMiIpJJc+b4dd08CnKAEMIIYESpY5elfL0AOKSc5w4EBmaiDhERkaj+F+QrZe0ltbKbiIhIpkRokSvIRUREMkVBLiIikseSQZ7FwW4KchERkUyZO9ev1SIXERHJQ+paFxERyWNz5kBRHbDsxauCXEREJFPmzMnq+XFQkEPdxK6pye4QERGR6pozJ6vd6qAgX74e7tSpcesQEZH8pyCPoKgO1KunIBcRkfTNmZPVVd1AQe4aNFSQi4hI+tQij6RhAwW5iIikb+5cDXaLokFD+O47WLQodiUiIpKvFi+GP/5QizyKBg0gBPjmm9iViIhIvvrfqm46R559DTVyXURE0hRhVTdQkLsGDf3666/j1iEiIvlLQR5R/XpQv75a5CIiUn0Rdj4DBXmCQadOCnIREak+tcgj69xZQS4iItWnwW6Rde4MX33lo9dFRESqatYsv1aLPJLOneG337R5ioiIVM+nn8Jaa0FRdqNVQZ7UqZNfq3tdRESqY9w42GKLrL+sgjypc2e/VpCLiEhVzZ4N06YpyKNSi1xERKpr3Di/VpBHtOqq0LKlFoUREZGqSwb55ptn/aUV5Kk0BU1ERKpj3DjPkGbNsv7SCvJUCnIREamOceOguDjKSyvIU3Xu7DugLVkSuxIREckXEQe6gYL8rzp3hqVLfW9yERGRyog40A0U5H+lkesiIlJVEQe6gYL8rzSXXEREqiriQDdQkP/VGmv4GrkKchERqaxIK7olKchT1akD6623vJtERERkRZID3SKNWAcF+d/tsgu89RYsWBC7EhERyXWRB7qBgvzvevTwEH/77diViIhIros80A0U5H+3446w0kowcmTsSkREJNdFHugGCvK/a9QIttlGQS4iIhWLPNANFORl69EDPvoIZs2KXYmIiOSqyCu6JSnIy9Kjh1+/9lrcOkREJHclz49HHLEOCvKybbEFNG2q7nURESlfDgx0gzSD3Myam9lIM/sycf23s/1mtqmZjTGzCWb2qZn1SrnvfjP72sw+Tlw2TaeejKlTB7p39yAPIXY1IiKSi3JgoBuk3yK/AHgthNAFeC1xu7Q/gaNCCBsAewD/NrOmKfefG0LYNHH5OM16MqdHD988ZfLk2JWIiEguyoGBbpB+kO8HPJD4+gFg/9IPCCFMDiF8mfh6BvAT0DLN1615yfPk6l4XEZHScmSgG6Qf5K1DCD8AJK5brejBZtYNqAd8lXJ4YKLL/SYzq59mPZmz1lq+G5qCXERESsuBFd2SKgxyM3vVzMaXcdmvKi9kZm2B/wLHhBCWJQ5fCKwLbAk0B85fwfP7mlmJmZXMyta0sB49YNQoWLIkO68nIiL54Z13wCz6QDeoRJCHEHYNIWxYxuVZYGYioJNB/VNZ38PMGgMvAJeEEMamfO8fglsIDAG6raCOwSGE4hBCccuWWeqZ79EDfv8d3n8/O68nIiK5LwR47DHYYQdo3jx2NWl3rQ8H+iS+7gM8W/oBZlYPeBp4MITweKn7kn8EGH5+fXya9WRW9+7+F5e610VEJOmzz2DSJOjVq+LHZkG6QX4N0MPMvgR6JG5jZsVmdk/iMYcCOwBHlzHN7GEz+wz4DGgBXJlmPZnVvLlP9FeQi4hI0tChPk35oINiVwJA3XSeHEKYDexSxvES4PjE1w8BD5Xz/O7pvH5W9OgBgwbBb79B48axqxERkZhCgGHDvMe21QrHd2eNVnarSI8esHQpjB4duxIREYlt3DiYOjVnutVBQV6xrbeGlVdW97qIiHhrvG5dOOCA2JX8j4K8IvXr+x7lCnIRkdpt2TIP8t13z4nR6kkK8sro0QO++MKXbBURkdpp7FjPgRzqVgcFeeVouVYRERk2zHtp96vSemg1TkFeGRtsAG3aKMhFRGqrpUvh8cdhzz1zbgaTgrwyzGDXXeHVV/0ciYiI1C5vvQU//AC9e8eu5G8U5JW1227w889arlVEpDYaNsxnMO29d+xK/kZBXln77guNGsHgwbErERGRbFqyBJ580kO8UaPY1fyNgryymjSB//s/ePRR+OWX2NWIiEi2jBoFs2bl3Gj1JAV5VZx0EixYAEOGxK5ERESyZehQWGUVH+iWgxTkVbHxxrDddnDnnRr0JiJSGyxaBE89BfvvDw0bxq6mTAryqurXD776Cl55JXYlIiJS00aOhLlzc7ZbHRTkVXfggdC6NdxxR+xKRESkpg0bBk2b+sylHKUgr6p69eCEE+D552HatNjViIhITVmwAJ55xjdIqVcvdjXlUpBXR9++vkjM3XfHrkRERGrKiy/C77/n5CIwqRTk1dG+vc8rv+ceWLgwdjUiIlIThg2DFi2ge/fYlayQgry6+vXzld4efzx2JSIikmnz5sFzz8FBB/n+4zlMQV5d3bvDOuto0JuISCF64QX488+cHq2epCCvrqIiOPlkGDMGPvoodjUiIpJJ993nu17usEPsSiqkIE9Hnz6+QIBa5SIiheP99+Hll+GMM6BOndjVVEhBno6mTeGII+Dhh33BABERyX9XXAHNm3uvax5QkKerXz+YPx8eeCB2JSIikq6SEhgxAs4+G1ZdNXY1laIgT9emm8LWW3v3utZfFxHJb/37Q7NmcMopsSupNAV5JvTrB5Mnw+uvx65ERESq66OPfMrZmWdC48axq6k0BXkmHHywLxpw++2xKxERkerq3x+aNIFTT41dSZUoyDOhfn04/ngYPhy++y52NSIiUlWffOLrqp9xhg9kziMK8kw58UQIAQYPjl2JiIhU1YAB3p1++umxK6kyBXmmrLkm7L03/Oc/vhG9iIjkh/Hj4ckn4bTTfKBbnlGQZ1K/fjBzJjz1VOxKRESksgYMgFVW8UFueUhBnkk9esBaa8Ftt8WuREREKmPiRN/86tRTfRGYPKQgz6SiIj+/8s47vqCAiIjktiuvhJVXhrPOil1JtSnIM+2f/4SuXb2LRufKRURy16RJMHSonxZt0SJ2NdWmIM+0evXgppt8gRh1sYuI5K6BA33jq7PPjl1JWhTkNWHPPaFnT194/6efYlcjIiKlTZ4MjzwCJ50ErVrFriYtCvKacuONvin9JZfErkREREq76irvQT333NiVpE1BXlO6dvU5iffc4+v3iohIbvjqK3joIV/Iq3Xr2NWkTUFeky691AdQnH66r/omIiLxXXUV1K0L550Xu5KMUJDXpKZNfTDFW2/5PEUREYnr66/hwQehb19o2zZ2NRmhIK9pxx7re5afe66fMxcRkXiuvtrX/Dj//NiVZExaQW5mzc1spJl9mbguc5FaM1tqZh8nLsNTjncys/cSzx9mZvXSqScn1akDN98M334L118fuxoRkdrrm29gyBDfrbJdu9jVZEy6LfILgNdCCF2A1xK3yzI/hLBp4rJvyvFBwE2J588Bjkuznty0ww5w6KFwzTXa5lREJJZrrgEzuKC8qMpP6Qb5fsADia8fAPav7BPNzIDuwBPVeX7eufZaH/BWQN05IiJ548MPfXfK44+H9u1jV5NR6QZ56xDCDwCJ6/Jm1TcwsxIzG2tmybBeDZgbQliSuD0dKJy+jtLWXNNHSD76KLz9duxqRERqj0WL4OijfarZwIGxq8m4uhU9wMxeBdqUcdfFVXidDiGEGWbWGXjdzD4DfivjceXO0TKzvkBfgA4dOlThpXPIeefBfff5dLQPPvABFyIiUrMGDIDPPoPnn8/L/cYrUmGShBB2DSFsWMblWWCmmbUFSFyXuR5pCGFG4noqMBrYDPgZaGpmyT8m1gBmrKCOwSGE4hBCccuWLavwFnNIo0bexf7hh3D//bGrEREpfOPG+Uj1Pn1gr71iV1Mj0m0SDgf6JL7uAzxb+gFm1szM6ie+bgFsC0wMIQRgFHDwip5fcHr3hm23hQsvhN/K6pQQEZGMWLhweZf6v/8du5oak26QXwP0MLMvgR6J25hZsZndk3jMekCJmX2CB/c1IYSJifvOB84ysyn4OfN706wn95n5dLRZs3wfXBERqRkDBsD48T7IrWnT2NXUmArPka9ICGE2sEsZx0uA4xNfvwtsVM7zpwLd0qkhL22xBRxzjP+FeMIJ0KVL7IpERApLSYlPNzv6aN+NsoBptFUsAwdCgwZ5vw+uiEjOSXapt2kDN90Uu5oapyCPpU0b31Tluefg5ZdjVyMiUjj694cJEwq+Sz1JQR7TaafB2mvDmWfC4sWxqxERyX8lJTBokJ++3HPP2NVkhYI8pvr14cYb4fPP4fbbY1cjIpLfUrvUb7wxdjVZoyCPbe+9/a/GCy+ETz6JXY2ISP664opa1aWepCCPzcwXh2neHA4+GH79NXZFIiL55/33vUv92GNrTZd6koI8F7RqBY895hveH3usb64iIiKVs2CBd6mvvnqt6lJPUpDnim239eVbn3qqoFcgEhHJuH/9y8ca3XMPNGkSu5qsU5DnkjPPhAMO8M1V3nkndjUiIrnvvffguuvguONg991jVxOFgjyXmMGQIb7laa9e8FOZe9CIiAj8tUv9hhtiVxONgjzXNGkCTzwBs2fDEUfA0qWxKxIRyU2XXw6TJtXaLvUkBXku2nRTn1f+6qu+QpGIiPzV2LFw/fVw/PG1tks9SUGeq4491lcmGjAAXnopdjUiIrljwQL//diuXa3uUk9SkOey226DjTbyLvZvv41djYhIbrjssuVd6o0bx64mOgV5Llt5ZT9fvmQJHHIILFoUuyIRkbief95b4SecALvtFruanKAgz3VduvhI9vffh3POiV2NiEg8JSU+o2ezzWrF9qSVpSDPBwceCGedBbfeCsOGxa5GRCT7pk3zvSlatfJWeaNGsSvKGQryfHHNNbDNNj5Cc9Kk2NWIiGTPnDm+fvrChTBihO9uJv+jIM8XK63k67E3bOibq8ybF7siEZGat3Chr3g5dSo8+yyst17sinKOgjyftGsHjzwCEyfCiSdqcxURKWzLlvk0szfe8F0id9ghdkU5SUGeb3bd1ffcfegh33NXRKRQXXIJPPooXH01HHZY7GpyloI8H118sa9kdOqpMG5c7GpERDLv7rs9wP/5Tzj//NjV5DQFeT4qKvIWeevWfr58zpzYFYmIZM6IEXDyydCzpy+MZRa7opymIM9XLVrA44/D99/DUUf5uSQRkXz34Ydw6KG+58SwYVC3buyKcp6CPJ9ttRXceKPPqTztNA1+E5H89s03sNdesNpq/nttlVViV5QX9KdOvuvXz3/4r7/eu59uuUXdUCKSf5JzxefPh9deg7ZtY1eUNxTk+c4Mrr3Wu9ZvvNHPn//73wpzEckfCxf6CpZTpsArr8D668euKK8oyAuBmbfIQ/D1h82WX4uI5LIQ4LjjYPRoH8S7006xK8o7CvJCYeY7Ai1bBjff7LdvvFFhLiK57dJL4eGHYeBA37JZqkxBXkiSLfEQvHu9qGj5uXMRkVzzn/94gJ9wAlx4Yexq8paCvNCYeYiHsLxFft11CnMRyS0vvggnnQR77AF33KHfUWlQkBciM+9eX7bMu9uLimDQIP1HEZHc8NFHPld8o418MyjNFU+L/vUKlZnvXx7C8hb5NdcozEUkrm+/9bnizZrBCy/AqqvGrijvKcgLmZkvbxiCT1ErKoKrrlKYi0gc06d7V/q8efDOO7D66rErKggK8kKXDPNly5a3yAcOVJiLSHaNH+8Lvvz6Kzz3HGy4YeyKCoaCvDYoKvLBJCH4bkJFRTBggMJcRLJj1Cg44ABo1Ajeegs22SR2RQVFQV5bFBXBnXd6yzzZIu/fX2EuIjVr6FDo0wfWXttHqnfoELuigqMgr02KinyP3xDgyiv99hVXxK5KRApRCD5r5txzYYcd4JlnfICbZJyCvLYpKoLBg/0/Wf/+fvvyy2NXJSKFZOlSOPNMnzlz6KHwwAPQoEHsqgqWgrw2KiryFZWWLYN//cu71y+7LHZVIlII5s+HI4+Ep56Cs87y6a9F2jG7JinIa6uiIrjnHm+ZX365h/mll8auSkTy2ezZsN9+8O67vlz0GWfErqhWSCvIzaw5MAzoCEwDDg0hzCn1mJ2Bm1IOrQv0DiE8Y2b3AzsCvybuOzqE8HE6NUkV1KkD997rYX7ZZR7uF18cuyoRyUdff+3Ty6ZNg2HD4JBDYldUa6TbIr8AeC2EcI2ZXZC4fX7qA0IIo4BN4X/BPwV4JeUh54YQnkizDqmuOnXgvvs8zC+5BGbO9I1W6tWLXZmI5Itx43y1tkWLYORI2H772BXVKumeuNgPeCDx9QPA/hU8/mDgxRDCn2m+rmRSnTowZIh3g916K+y8M3z/feyqRCQfvPQS7Lgj1K/vq7UpxLMu3SBvHUL4ASBx3aqCx/cGHi11bKCZfWpmN5lZ/TTrkeqqU8fPaQ0dCp98Aptv7os4iIiUZ8gQ2HtvnyM+Zgyst17simqlCoPczF41s/FlXParyguZWVtgI+DllMMX4ufMtwSaU6pbvtTz+5pZiZmVzJo1qyovLVXRqxe8/z40bw677uprtIcQuyoRySXJ6avHHgvdu8Obb2rd9IgqDPIQwq4hhA3LuDwLzEwEdDKof1rBtzoUeDqEsDjle/8Q3EJgCNBtBXUMDiEUhxCKW7ZsWdn3J9Wx/voe5gceCOefDwcd5Osji4gsWQJ9+/psl6OOguefh8aNY1dVq6XbtT4c6JP4ug/w7AoeexilutVT/ggw/Pz6+DTrkUxZdVXfJ/iGG2D4cNhyS9/0QERqrz/+8Oll99zjM1zuv18DY3NAukF+DdDDzL4EeiRuY2bFZnZP8kFm1hFoD7xR6vkPm9lnwGdAC+DKNOuRTDLzBR1GjYLff4ettoJHHoldlYjEMHOmD4R96SW46y5f5ll7NeSEtKafhRBmA7uUcbwEOD7l9jSgXRmP657O60uWbL89fPihnz8/4ggf1HLDDfpLXKS2mDzZ9xH/8UdfM32ffWJXJCm0bp5UTtu28Npr3kK/7TbYaSeYPj12VSJS0557DrbZxnvlRo9WiOcgBblU3koreUv8scfgs898itrrr8euSkRqwq+/wjHHwL77Qrt23hPXrdzxyBKRglyq7pBDfFR7ixbQowcMGqQpaiKF5NVXYaON4MEHfVDbBx/4XHHJSQpyqZ711vMwP/hguOACOOAATVETyXd//AH9+vkf6Cuv7JufXHmlxsPkOAW5VN8qq/hKcDfdBC+8AMXF8OmnsasSkep4+23YdFO4807fS/yjj3ymiuQ8Bbmkx8zXaB81CubNg3/8Ax56KHZVIlJZCxbAuefCDjvAsmX+f/nGG6Fhw9iVSSUpyCUzttvOp6h16wb/93/ePbdwYeyqRGRFSkp80Or11/tqbZ984hugSF5RkEvmtGnjg2TOOQfuuMN/IXz7beyqRKS0RYvgssu8B+2335Yv8rLqqrErk2pQkEtm1a0L110HTzwBEyf6uu2DBql1LpIrPvvMA3zAADj8cF96effdY1claVCQS8046CD4+GPYZRcf1b7RRvDii7GrEqm9li71P6qLi30xp6ee8ullTZvGrkzSpCCXmtO5Mzz7LIwY4bd79vTFJb76Km5dIrXN5Mk+juWCC3xltgkTfMqoFAQFudS8Pff07rxrrvGV4DbYAC69FP78M3ZlIoVt2TK45RafVjZpEjz8MDz+OGgr6IKiIJfsqF/f9zb/4gvvdr/ySlh3Xf+lolXhRDJv2jTYdVc4/XTfG2HCBD8nrh3LCo6CXLKrXTtvFbz5JjRvDoce6r9sJkyIXZlIYQjB9wvfaCNfWvWee3zBptVXj12Z1BAFucSx/fY+h/X2230FqU028dWktMyrSPV9+y3svTeccIIPavvsMzjuOLXCC5yCXOKpWxdOPtkH4hx3HNx8M6yzDgwZ4uf2RKRyvv8eTjkFunTxldluvtm3He7YMXZlkgUKcomvRQu4+27vBlxrLTj2WN//+IMPYlcmktt+/NGXSF5rLf8/dPTRPqjttNOgSL/eawt90pI7ttjCN2544AEfqLPVVnD88TBrVuzKRHLLrFm+PnrnznDbbXDEEd6zdffd0KFD7OokyxTkkluKiuCoo/yX0plneqivsw7ceissWRK7OpG4Zs+GCy+ETp18Y5ODD4bPP4d77/VjUispyCU3NW4MN9zg26IWF3tX4eabwxtvxK5MJPvmzPG1Fzp18tXZ9t3XZ3o8+KCfF5daTUEuuW299eCVV+DJJ31zh512gt69fYlJkUL366/Qv78H+JVXwh57+Ej0Rx7xdRhEUJBLPjCDAw/0TVguuwyeecYH9/TpA+PGxa5OJPN+/x2uusoD/PLLoXt332L0scd8ZUSRFApyyR8rrwxXXOGjcvv29U0fiot9DenHHoPFi2NXKJKeefPg2ms9wC++2H+2x43zn/WNN45dneQoBbnkn44dffDb9Olw003www/Qq5eP4L36avj559gVilTN/Pk+eK1zZ1/KuFs3eO89GD7cx4aIrICCXPJXkyY+h3byZP+Ft+66cNFF0L69T1v79NPYFYqs2IIF/kdp585w9tm+wuG77/qOgd26xa5O8oSCXPJfnTq+NePIkTB+vJ87f+QR/6W4007w9NO+F7NIrli4EO68E9Ze22dkdO3qMzJeeQW23jp2dZJnFORSWDbYAO66y7vdr70Wvv7aB8qttRZcf71P4xGJZcYMn1a5zjq+PHHHjr617+jRsMMOsauTPKUgl8LUvLmvfPXVVz51rWNHv73GGnDSST4CXiQbfvvNFzbq0cN//s45x0//vPIKvPUW7Lxz7AolzynIpbDVrest8tGj4eOPfQ76kCHect9tN3j+eW3QIpm3eLH/bPXuDW3a+BroU6f6oi5ffOFLEffooV3JJCMU5FJ7bLKJL2U5fToMHOit8n328W7Om2/2lpNIdYUAY8b4LmSrr+4/W6++6psAvfsuTJni0yfXWSd2pVJgFORS+7Ro4aPbv/4ahg6F1q199Hu7dj7waPLk2BVKPpk82Rdt6dLFd+27917YZRd47jmfGnnbbT6ATa1vqSEKcqm9VlrJ55+/845vmXrAAT5QrmtX2HFHHyw3fry3tERSzZwJt9ziU8S6doUBA3wRl/vv9/uGDoW99/afMZEapiAXAV8h7sEH4dtvfW3r337zhTk22gjWXBNOPNHnqv/xR+xKJZZ58+Dhh6FnT++9Of1035Hv+uvhu+98+mOfPr7hj0gW1Y1dgEhOadPGByRdeil8/z289JIvzvHII77Xc7163lrv2dMvXbqoy7SQLVkCr70GDz3k6xHMm+f7fZ93nu8BrnXPJQcoyEXK064dHHecXxYt8i74ESP8cuaZfllrreWhvuOO0LBh7KolXSH4+uYPPeRd5DNnQtOmcPjhcOSRvv55kTozJXcoyEUqo149n++7885w3XUwbRq8+KKH+j33+DKbDRv6LlXJYO/YMXbVUhlLl/pYiDFj/PL22z5VrF49P8995JH+edavH7tSkTIpyEWqo2NHX1jmpJN8vew33vBQf+EFv4DvpZ4M9e2282CQ+H7+GcaOXR7c77/vXeYArVr5CPPzz4dDDoFmzeLWKlIJCnKRdDVoALvv7pebb4Yvv1zeBX/rrb4k5yqr+AIgPXvCnnt6t73UvNKt7TFj/PMBX6N/0019sZatt/ZLp04a8yB5R0EukmlduviI5tNP91Huo0Ytb60//bQ/ZpNNfM7xuut6y33ddX35ToVIembP/ntrOznTINnaPu44vy4u9j3uRfKcglykJq2yiq/wtc8+Pohq4kQP9RdfhEcfhblzlz+2USMP9OQlGfBrr63zs2VJbW0nwzu5mE+dOv7HUp8+am1LwUsryM3sEOBfwHpAtxBCSTmP2wO4GagD3BNCuCZxvBMwFGgOfAj8XwhhUTo1ieQsM5+utMEGvoFLCPDTTzBpEnz+uV9PmuQbaTz88PLn1anj+1WXDvh1160d53CXLfNd63780VfjS4Z2Wa3tY45Z3tpu1Chu3SJZkm6LfDxwIHB3eQ8wszrA7UAPYDrwgZkNDyFMBAYBN4UQhprZXcBxwJ1p1iSSH8x8edjWrX3qWqp587x1mRrwn38OL7/sU+GSWrX6a7Anv27fPvenSM2b5+Fc0WXmTN+EJEmtbZG/SCvIQwifA9iK/wN1A6aEEKYmHjsU2M/MPge6A4cnHvcA3rpXkIs0agSbbeaXVEuX+tS31ICfNAkefxx++WX541Ze2ZcO7dIFmjTxAXmZuNSv70FansWLvZehMgFd1ip5RUX+x0nbtr44z0Yb+XXy0q6dD1BTa1vkf7Jxjrwd8F3K7enAVsBqwNwQwpKU41kfyjtliv8+2WmnbL+ySHXUAdZKXPb2Q3WBjYKH6J9/Lr988ydMmu/hv2zZ8kuFFiUu5e0GZx64qRczf/0lKS1nGi2vtU5dn36XvDSpBy3r/fVYvXq+NnmyYfAnMDVxEckjH3/sw2OypcIgN7NXgTZl3HVxCOHZSrxGWc31sILj5dXRF+gL0KFDh0q8bOW0bJmxbyUSkcFKiYBs0nQFjwt+01qrCwAAB9ZJREFUbj412DNyCR7CpYM5ecn1bn6RDFpllexmS4VBHkLYNc3XmA60T7m9BjAD+BloamZ1E63y5PHy6hgMDAYoLi7O2HZUH32Uqe8kkg8scVGwihSKbPxv/gDoYmadzKwe0BsYHkIIwCjg4MTj+gCVaeGLiIhIQlpBbmYHmNl0YGvgBTN7OXF8dTMbAZBobZ8CvAx8DjwWQpiQ+BbnA2eZ2RT8nPm96dQjIiJS25g3jPNLcXFxKCkpc8q6iIhIwTGzcSGE4rLu04kyERGRPKYgFxERyWMKchERkTymIBcREcljCnIREZE8piAXERHJYwpyERGRPKYgFxERyWMKchERkTymIBcREcljCnIREZE8piAXERHJYwpyERGRPKYgFxERyWMKchERkTyWl/uRm9ks4JsMfssWwM8Z/H4x6b3knkJ5H6D3kosK5X2A3suKrBlCaFnWHXkZ5JlmZiXlbdieb/Reck+hvA/Qe8lFhfI+QO+lutS1LiIikscU5CIiInlMQe4Gxy4gg/Reck+hvA/Qe8lFhfI+QO+lWnSOXEREJI+pRS4iIpLHak2Qm9khZjbBzJaZWbkjCc1sDzP7wsymmNkFKcc7mdl7ZvalmQ0zs3rZqbzMGpub2chELSPNrFkZj9nZzD5OuSwws/0T991vZl+n3Ldp9t9F5d5H4nFLU2odnnI83z6TTc1sTOLn8FMz65VyX/TPpLyf/ZT76yf+nack/t07ptx3YeL4F2a2ezbrLq0S7+MsM5uY+AxeM7M1U+4r82ctlkq8l6PNbFZKzcen3Ncn8fP4pZn1yW7lf1eJ93JTyvuYbGZzU+7Lmc/FzO4zs5/MbHw595uZ3ZJ4n5+a2eYp99XMZxJCqBUXYD2gKzAaKC7nMXWAr4DOQD3gE2D9xH2PAb0TX98FnBTxvVwLXJD4+gJgUAWPbw78AqycuH0/cHAOfCaVeh/AH+Ucz6vPBFgH6JL4enXgB6BpLnwmK/rZT3nMycBdia97A8MSX6+feHx9oFPi+9TJ4fexc8r/hZOS72NFP2s5/F6OBm4r47nNgamJ62aJr5vl8nsp9fhTgfty9HPZAdgcGF/O/T2BFwED/gG8V9OfSa1pkYcQPg8hfFHBw7oBU0IIU0MIi4ChwH5mZkB34InE4x4A9q+5aiu0X6KGytZyMPBiCOHPGq2q6qr6Pv4nHz+TEMLkEMKXia9nAD8BZS7wEEGZP/ulHpP6Hp8Adkl8DvsBQ0MIC0MIXwNTEt8vhgrfRwhhVMr/hbHAGlmusbIq85mUZ3dgZAjhlxDCHGAksEcN1VkZVX0vhwGPZqWyKgohvIk3jMqzH/BgcGOBpmbWlhr8TGpNkFdSO+C7lNvTE8dWA+aGEJaUOh5L6xDCDwCJ61YVPL43f/9PMTDR7XOTmdWviSIrobLvo4GZlZjZ2OTpAfL8MzGzbnjL5KuUwzE/k/J+9st8TOLf/Vf8c6jMc7OlqrUch7eeksr6WYulsu/loMTPzRNm1r6Kz82WSteTONXRCXg95XAufS4VKe+91thnUjcT3yRXmNmrQJsy7ro4hPBsZb5FGcfCCo7XmBW9lyp+n7bARsDLKYcvBH7Eg2QwcD7Qv3qVVvj6mXgfHUIIM8ysM/C6mX0G/FbG4/LpM/kv0CeEsCxxOGufSXlllXGs9L9nzvz/WIFK12JmRwLFwI4ph//2sxZC+Kqs52dBZd7Lc8CjIYSFZnYi3mPSvZLPzaaq1NMbeCKEsDTlWC59LhXJ+v+TggryEMKuaX6L6UD7lNtrADPw9XKbmlndREskebzGrOi9mNlMM2sbQvghEQo/reBbHQo8HUJYnPK9f0h8udDMhgDnZKToMmTifSS6oQkhTDWz0cBmwJPk4WdiZo2BF4BLEt1uye+dtc+kHOX97Jf1mOlmVhdogncxVua52VKpWsxsV/wPsB1DCAuTx8v5WYsVGBW+lxDC7JSb/wEGpTx3p1LPHZ3xCiuvKj8jvYF+qQdy7HOpSHnvtcY+E3Wt/9UHQBfz0dD18B+o4cFHKozCzzUD9AEq08KvKcMTNVSmlr+da0oETfI88/5AmaMvs6DC92FmzZLdzGbWAtgWmJiPn0niZ+pp/PzZ46Xui/2ZlPmzX+oxqe/xYOD1xOcwHOhtPqq9E9AFeD9LdZdW4fsws82Au4F9Qwg/pRwv82cta5X/XWXeS9uUm/sCnye+fhnYLfGemgG78ddeuWyrzM8XZtYVHwg2JuVYrn0uFRkOHJUYvf4P4NfEH+o195nU9Ai/XLkAB+B/ES0EZgIvJ46vDoxIeVxPYDL+197FKcc747+cpgCPA/Ujvpf/b+f+cRKIggCMf1JZWMkFCBWlhZU3oPEENEDDHWjs9AJ2lp6AGvQKEmJhoJRwCBqKHZMXYqIN4U3y/RKy/NlJ3rxZdsLLLm1gAaxjex3v3wIvxX4dYAu0juLfgBVNs3gFrmrNA7iLsS5jO85aE2AA7IGP4nFTS01+O/Zplvfv4/llzPMm5r1bxE4j7gvon6sO/8xjHueAnxrM/jrWKs7lEfiMMb8DvSJ2FLXaAMPac4nXD8DTUVxVdaH5YbSL7/I3zXUWE2ASn18Az5HniuIuqVPVxH92kyQpMZfWJUlKzEYuSVJiNnJJkhKzkUuSlJiNXJKkxGzkkiQlZiOXJCkxG7kkSYkdAPO9bmZCrKrUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[8, 8])\n",
    "angle = np.linspace(0, 2*np.pi)\n",
    "# Red circle \n",
    "plt.plot(np.cos(angle), np.sin(angle), 'r')\n",
    "# Blue square\n",
    "plt.plot([-1,1], [-1,-1], 'b')  # bottom side of square\n",
    "plt.plot([1,1], [-1,1], 'b')  # right side  of square\n",
    "plt.plot([1,-1], [1,1], 'b')  # top side of square\n",
    "plt.plot([-1,-1], [1,-1], 'b')  # left side of square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this fact as follows:\n",
    "- generate a list of $N$ random points in the square $[-1,1] \\times [-1,1]$ that circumscribes the unit circle,\n",
    "by generating successive unifromly distributed random values for both the $x$ and $y$ coordinates.\n",
    "- compute the fraction of these that are inside the unit circle, which should be approximately $\\pi/4$.\n",
    "($N$ needs to be fairly large; you could try $N=100$ initially, but will increase it later.)\n",
    "- Multiply by four and there you are!\n",
    "\n",
    "Do multiple trials with the same number $N$ of samples, to see the variation as an indication of accuracy.\n",
    "\n",
    "Collect the various approximations of $\\pi$ in a list, compute the (sample) mean and (sample) standard deviation of the list, and illustrate with a histogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"Exercise-B\"></a>\n",
    "#### Exercise B\n",
    "\n",
    "It takes a lot of samples to get decent accuracy, so after part (a) is working,\n",
    "experiment with successively more samples;\n",
    "increase the number $N$ of samples per trial by big steps, like factors of 100.\n",
    "\n",
    "For each choice of sample size $N$, compute the mean and standard deviation, and plot the histogram."
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}