uncertainty/equilibration_uncertainty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 2, Morning Session 1, MD/MC workshop\n",
    "### Dr. Michael Shirts, CU Boulder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do we understand the information we get out of simulations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load up some modules.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load up some typical data from a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential = np.loadtxt('potential.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEGCAYAAABcolNbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA28klEQVR4nO3deXhU5fXA8e/JHiABQtgTCZsgiyAERKWKsrhWXKu2bq2te2ur1n3frVVrWzese3H7WS0qIioKLkV2kF12CGsgEEJCyHZ+f9w7k5lkJplAJpNkzud55smd9947c24Cc+Zd7vuKqmKMMcaES0ykAzDGGNO8WaIxxhgTVpZojDHGhJUlGmOMMWFlicYYY0xYxUU6gMYmPT1ds7KyIh2GMcY0KfPmzdupqu0D7bNEU0VWVhZz586NdBjGGNOkiMiGYPus6cwYY0xYWaIxxhgTVpZojDHGhJUlGmOMMWFlicYYY0xYWaIxxhgTVpZojDHGhJUlGmOaAFXl21W5rNtZGOlQjKkzu2HTmEYsv6iU6T/t4IZ3FgKQEBfDTw+dGtmgGgFVRUQiHYYJkdVojGnEBj3wuTfJAJSUVUQumDC6duI8/vTuQu/zfQfKWB+k9lZWXsHhd03h2a9XA5Czu4i9xaUNEaY5SJZojGlE9peU89z01ZSVV7B6x74Gfe8DZeXkFZYc0mus21nIO7M3ArApr4jNe/aHdN6ni7fx4YLNeFb8/dVLPzDqr9O9+7flF5NXWMIN7yzg7TmbKC1Xnvx8JQAjH/+an//ju0OKu7nLLTjgl4y35u9n/sbdDfb+1nRmosKstbtYvnUvlx/XPazvM3XpNu74YDHf33YSSfGxAY/ZW1xKalJ8wH3PTFvFCzPW8MWy7SzYuKfa/vRWifUZrp+r35zH1ytzWf/Y6Qf9Guc9/z92FZZw3tAMfvaXrwG4aHgmY/t15KS+HasdvzZ3H799vXJuwR/W5nFMz3YsyskH4PvVOzmmRzvGPjWDggNlAExauMV7fEWFk5g27Co66Jgbk9U7CogRoUf7VvX6usMe/pK0lgnMv3ssAOOe+oaCA2WH9LeuC6vRmKhwwYQfuO/jZX5lX6/cwbGPTqO4tLze3ueeSUvYVVjC9r3FfuVfrdjOq9+v46ftBRx53+dk3TaZtbn7KC4tp6SsgvIKZfrKHbwwYw1AwCQDsHPfAT5bspU9RU7No7S8ImgTU119vTIX4JB+H7vcGlFhSeVrvD17E/d95P+7/3TxVo577CtmrctjrU/8F730g9/7/+pfs3hw8jJvkvFVobBnf+W39JzdgZNNwSE0q63bWcizX69GVVFVyt3EFi5jnvqGk56cUW+vt2Djbna7fxPf2qrn9+mpQYabJRoTVYpLyyktr+DDBTk88PEytuQXs37XoX1Qf750Gyc88TXLt+5l737nP/Bpz3wLwP/N3UTWbZP5zWtzuf/jZbw/L8d73o85+fS9+zMOv2sKPe/4lMtfnRPS+1397/lkP/QlAI98upxRf51ObsGBoMev21nIuKdnMHd9HoU+H9gzfsoNmFRe+mZtSHEEEhfjdNDv3e//4b4xr4hjHp3Gu3OcZrW7/ruEzXv2c/sHi6u9xqD7P/d7/ur364O+3659ldc98vGvvdvTlm8n+6EvmLRwMwPv+5yFm/Z496kqt7y/iKzbJvPtqtwar+fif83iiakryd9fyqvfr6fnHZ8eUn/Qx4u28OzXq1m1vaDG40rLD70vTlU5+7n/8ct/zfIr/3zpNu/2o1NWNEiysURjosrV/57HXz5bwZ/eXeQdKlxQXPnhu2p7ASMemcaWKn0LB8rKufPDxX41leLScqYu3ca7czaxYVcRpz7zLfvdD+7CEuf4ez9a6vc6E3w+xMvq8O34r+cPYnTfDn7n/uKFmd4P4bW5wftzbv3Pj/y0fR/nvTCTc5//HwDLtuzlsldm8+AnlTWNWDdJPPnFTyHH5VFeoSzatMd7TZ5mM19b84u59T+L+XZVLvn7g39YH6gy4OG8oRlBjx379Dd+z8/4x7dk3TaZK16fy859Jfz7B2fm+iWb873HzN+4m/fmOgl/8o9ba7wuTx/Topx8nnZ/L9vzi8kvKj2oZPD7txfwxNSVjH36m2p9JJ6+LYDed04J6fXmb9ztbT6sylOrXL51r7ds4aY9XPnmPO/zCd+sZc768PfVWKIxUWX6ylw+WrTFr2ynT23gmWmr2La3mBk/+X/Tnbt+NxNnbeTW//wIOM0x2Q99yVVvzmPaih0B32virI0UlQRvhrr5/xaFFPPdZ/TjvKEZvHz5MFKTKrtVZ6/P825fMOEHAL5btZNf/esHyiuUBz9Zxi3vL2L2usrjVmwr4MFPlnHa350a15rcfXy6eCuqSrJPn1JFhbKjoJgNuwp5fvoaXv5uXcDYFm7aw5LN+bz4zRrGP/t9SNdzycuz69QE5VsLrM2SzXv9nns+RBNiY7w1kT1FlUmuQpX8/aVk3TaZs6rE/9vXK2uYl70y29vctKPgAIMe+Jzed07xq1HVZG9xKVe+4b/O1dLN+azc5tRsKiqU26rU7mpr8pu2fDvnPPc/3p6zsdq+igplgtsM62vKkuqJ1TN6L5xsMIBplhZu2sNZz37P5D+MpH+X1n77tu/1/3C4ZuJ8wPlA3+Emnds/WMxFww/zHuO5Y8PT6fy7N+ayL0C/QX0altWWOet3U15R+c05Jqbme0euf3s+e4pKySssCZocfMt/WJvHD2vzGJzZxu967p60hImzNtI6Od5b+7hiZPWBFFU/nIO5/dS+tEiI5e5JS4Me07dTCiu2BW9SeuM3wzkyozWDH/gipPf09dbsjdzynx+57sSe9O2U6i1/b26Ot3azcNMetu8tpqxCadsini+XB/4C8SufpqihD33JPWf0Y03uPv409vCAgzVufG8hH8zfXK3c87t484rhDO3Wttr+7XuLSfEZNPLdqp18syqXYVlp/O6NuZxzVFcANu+uPrJvxqpc/v5V9QTy4ozqzaJVv1SFg9VoTKNXUlZB1m2TuXbivKBDfisqnDvnN+UVMX/jbu8H4CUvz6a0vII+HVNqfZ8HP1nm90378ldn89vX51BQXMpet3lt3c5CFm3aww9r84K9zEE7fWBn7/agjNYM7ZYG+DexTbgkm+FZaQHPf39ejrePpK5J0LcPA5zaGODXxLU1fz8rtvnXGEJ11Qk9Gdk74Cq/vP6b4XRtk8x/rzuOT34/0lvepoX/yLzjD29PmxYJDOjqJIrsAB/OwXiu79mv13BTDTXJox+ZxnGPfcUZfw99uPTDny5n4qyNZD/0Jb94YSb7DpTx5bLtzF6XxwMfLwuYZHxd8vJs9hVX/3u9M3uT3/OLX57FhG/W8vDkZX7X9N8Fm1myOZ97Ji2hokKZv3F3wNerSUlZRVj7aqShRh00FdnZ2WpLOTcueYUlDHmw8lvsk+cP4lyfdvv35mziFrdJC6BdywTv6Ke6io2Ras06V4zszp6iUv4zP/QmnNose+Bk+t0z1a/suhN78uzXa+iUmsTUPx3PS9+s5Z9fr+bGsYfzh9G9/Y69d9ISXp8ZdOVcbhjdm2emrfI+D3RdB+vaUT2ZunQbN47tw3VvzQ963GkDOzF9ZS7Du6fx2q+HU1xaTt+7P/M7ZuE9Y2nTIsGv7NXv13H/x8sY268jXyzb7i33DMU9UFZOflEpG/OKOO+FmcTGCI+fe6RfU+TwrDS/psWD9bPe6Xy7amedzomLkTr1vwG0TIj1G6nnMe+uMbRza0lZt03225eSFOfXv+irfUpijQNEAjl3SAZP/mJQnc7xJSLzVDU70D6r0ZhGr+rd8Df93yLGP/s9RSXOfzLfJBPo+KoeO2cgd51+RMB9gT6MX/5uXa1J5pdHH1bj/qpaJMSx8J6xLLpnHKcN7ARA306pvH/1MXx50wm0To7njEFODefk/p2qnR8XW/lf94TDq9cUfJMMOMk3FJcfm8XdZ/Sr8Zjnpq9hTW5htSRzSpU4Lz+2O0vvP5lXLx8GUO2+osuO6VYtyQAkxsV6Y55yw88C7u+QmuRtVjp3SNdqAwbevWoEL14ytMbr8OXbTOrrJJ8BGKGqa5IBAiYZgPW7CnlxxhoKD5RRdcadYEkG8CaZx84ZGPSYBXeP5SmfxFKfX6SqskRjGp2KCuWHtbsAZ4hmoCG4izbtCTpaJtA9F74uGJbJFSO78+0tJzL7ztGHHjBw/5n9q5WN7JXu3W7TIp4fbh/N5D+MZPrNo9yyBFq3iMfTqCAC2VlptEp0uk77dkpl/WOn06dT9Wa/uNjKT53h3QM3pfnyfCvu2iaZo2s4/r4z+3PqgOqJLRTPXzyEC4dlep+3bRGPiASck2zKDT/j/vEDAr6O59rKKpQOKcFvUO3TKYVXLx/G/Wc6rzMow+mLm33naESEAV1bBz3Xo1/nVGJjxNsc5/HA+P58ffMozs/ODHJm3bx6+TAS4io/brO7tSUpvvaP33Ofn8mjU1bQ/96pHEzj02lHdg66r23LBM4c1KXuL3oQLNGYRueFb9Zw4YQfuPHdhRx+1xS+CXKvw979pVRUKPGxwvlDM3jrt0fX+trDs9K8H36ZaS3okJLEpOuO44qR3WmZEPhOfoDMtGTWP3Y6lx7TzVu2+L5x3u342BjuPM2pJQ3Lasuk647jhjGVzV2dUpPo1DqJ/l1ak5Xe0u+1vYmG0CeJTPCp0fRIb8mie8YFPO62U/sClR35H153bLVv//27+H/IdmmTzEuXZvOz3umE6u4z+iEiPHbukd6yDilJQY8/onNq0H2efqay8gpvggzmxL4dSHb/bq9cPowJlwz1vm/n1Mr3/+cvj6p27rE92zHp+uNY9sDJfrWDr246gUuPyaJ7ektaJcbxwPjKLxHTbjqB7249ka9uOqHGuLq0rnzvxLgYTuzbgXOHOL/3B88awPvXHMuKB09l6f0nc/GIytrUg+Orf2GpzTMXDg66r1VCHG/9zvl/8eeT+5DVrgUAN449HPCvGYeTjTozjcqkhZv5y2fOHFYfLHA6Ue8JMlLp40VbuPejpZSWK4My29A7SIf/g2cN4O7/LuHCYZl+H4QegzLbMCizDXef0Y/i0nISYmNYtnUvE2dt4G23Q/bso5wPifvP7M8bMzcQHyukJMXzxzG9vffj/O74HvTq2IrBGW1o2zLB7/6Fa0b1DHrNw7qn8dnSbXRzPwRCccXI7vzDHVWkQGpyHJ1Sk9hWZUaCq0/oydUnOO/tSTDnDc1gVJ8OFBSXkpwQ671p8sNrj/WeN7ZfR37M2RNy/0Sgb+epydU/Xl799TB+3JRfrdyX58PP0wT16R9+RosavgR4tGuVyDif5ruYGOG+n/ejb+dURvRox/VvLQDgu1tPJKOt/+/6vKEZfLV8B3+/6Cg6tfZPkJeM6Ob9N9jTnRqm6n1W71w5ggvdIeYA5wzJ4IJhmUz/KZdjejg1yBE90nh79ka/gSktE+M4uX8n/v2DM/ji8I4pTPzt0ewuKvHGW5NTB3Ri/OCufhOv+oqJEY7tmc5/rzuOAV1SuXhENyoqlLYhNqXWF0s0ptHILyoN+h8mkM99Oop7dWhVrR8iNkZQVc48sgvnD83wqwUE4+lHGNC1NQ+fNZD9JeX8d+EWUtzmLBFh5u0nefty/jjmcL/zT+xT2aYf777fYWktGD+4a9D3/M1xWYw5ogPd2rUMekxVbVokcNHwTN6evYnEuBhEhG9vPdF7o9+MP48KOiBCRGifkkh7t1nqifMG8c7sjQzKaON33FUn9PQms8y0ZP550ZCg98rE+/xuHxjfny+X7wjYZHZinw5+v6OA15bs9L2kuX/Pfl2C135qE2huu6pJBpw55N67+piArxHoOlol+X90tkzwf56cEEtmWgsuGVFZAx4/uCvZWWl0bZPsd+yQw9pyzlFdGdGjHUf3aAdUT2S++nVOZZn7Jeaxc/y/OA3t1pbFm/O5cezhtPUZtTc4sw0ArZNr/j8wcdYGfjn8sHpfgsESjYmowgNlrNqxj6179ldrUqqLfl1SiYkRXr4smyten8ugjNZMun5k7SfWICZGvO3qyT7fqDu3Tg52ip/O7jfjP47pXeNxIlKnJONxx2lH0K1dS+8Hd3xsDLPvGE1MjJDeKjHk12yfksjvR1ePsVViHKP6tGf6ylzuOaM/g9wPK4/7ft6P8YO7MmvdLk44vDJ5XHpMFpcek1Xn6/H4We90Hjl7IGcd1TD9B6G4+4x+DDmsjfd5alI8s+8YzfBHpgHQIdW/iS8myAd11SQDTq3mqQsG+5clBP5oHtqtLf+55lg+W7KV4tIKWrvJ5NVfDyPOrb1A5SwPdfXZkm38MsjAiENhicY0uJKyCuJjheemr+GD+TmsyXWankIZuXXawE58unibX9kRnVO9syF7ksD+epoo09NvcjBf8FomxoV1dtyUpHhvs5hHh9Tg/SIH4/cn9WLltgLvgIP0Vgns3OfUlDy1hVMGBO9wPhgiUudRfKFonRxf40CImgS6WbVDahLDu6fRr3MqHVOTWHjPWP4+bTWvfL+OQ+36aJnofLG547S+PPLpCgZ0TeXUAZ29nfdVf+e11RJDdcqATmFZUM4SjWlQ2/KLGfHoNO46/QiemLrSb9/2/OIgZ1XyndLl8mOzuOqEHrRrWfltsleHVozq075ak9bB8vyfi9bbzYZ2S2Pm7ZUj87695SQunDAz4Ei4xm7RvYEHTByK966qbG5r0yKBCvcfSrAaTajiYmO8X1LOH5pJckJs0GUnDtX5QzP4P3eaH8/Q8vpmo85MWD39xU8scu9gXrBxNy9960yB8dDk5dWOLa3h/oPj3XtFTvP5Jnffmf3p3DrZb9hoQlwMr/16uLdN+lB5E029vFrTl5wQy6TrR/KX8w7+xr7m7JieTh/LUT7NbIeqbcuEsCUZgMfPPZKHz3aGiB/ROTxfIKxGY8KmrLyCZ6at4plpq1j36Gmc/dz/ajz+myBzLv143zjvhI9xMcKHCzaTs6dhFroa3j2Nt2dvom8T/AZvGt7J/TsFnO2gMYuJEX45/DDGHNGRjvXc9OphicaEjW8/yca84InBt+oOMOuO0aQmxXPEPc50JVVXo3z7yhH1HGlwZx+VwfDu7QJ24hoTSFNKMh4iErYkA9Z0ZuogZ3cRb/usmeFRVFIWcOoW34kdg02GCfjdwf30BYPomJrkHeV1sKNn6pMlGWMOjdVoTMiueG0uK7cXcOqATt5vbapKv3umctHww3i0yrxKZ/6z8p6LK16vPlHp0xcM4siMNmS2bUH+/lJydhd5b4wEZ8LFcIyAMcY0LEs0JmSeGwB37ivxJhrP9Plvz95IZloyvTukMLZfRwqKS2ucPfalS7MZ26+j93nV2YmhaTZBGGOqs0RjQuaZZuT9eTneObTyfVYr9EwdM7JXuvdGMl/nDc3wrpY45oj6GfdvjGn8rI/GhMyzENULM9awbMtenvx8Jcc/UX1t+O9W7/Rbiz0xLob+XVK51me+L2sSMyZ6RKRGIyIPAuOBCmAHcLmqbnH3HQm8CKS6+4eparGIDAVeA5KBT4EbVFVFJBF4AxgK7AIuUNX17mtdBtzlvu1Dqvp6w1xh87Bo0x7unrSEH3Py+dsFg0nzuTHSs+Z8KP7926MZFmRVSGNM8xepGs0Tqnqkqg4GPgHuARCROODfwNWq2h8YBXjaZp4HrgR6u49T3PIrgN2q2gt4Gnjcfa004F7gaGA4cK+IhL72q+H3by/gxxxnpt0/vrvQb1nfUDxz4WBm/HmUJRljolxEEo2q+i483pLKG6/HAT+q6iL3uF2qWi4inYFUVZ2pztrTbwBnueeMBzw1lfeB0eK0y5wMfKGqeaq6G/iCyuRkavHk5yur3fuSk1dEeqvQO+jbB5jY8YwjOzP6IFYtNMY0XRHroxGRh0VkE/Ar3BoNcDigIjJVROaLyC1ueVfAd53RHLfMs28TgKqWAflAO9/yAOdUjeVKEZkrInNzcwPfnR5NVNU7PbyvXYUlnBJg9cWqC2d5pAVISv/85RBedpf2NcZEh7AlGhH5UkSWBHiMB1DVO1U1E5gIXO+eFgeMxEk+I4GzRWQ0BFx60FMLCravpnP8C1UnqGq2qma3b199/fVoU1JeEXTfkRltuPWUvn7J5a3fjQi45nr7WlZHNMZEh7ANBlDVMSEe+hYwGac/JQeYoao7AUTkU2AITr+N7/qzGcAWdzsHyARy3D6e1kCeWz6qyjnTD+JSosrmPfs5vYaO/pTEOH6Rnck1o3qSddtkwJl+/dFzBjKuX0cS42Pokd4KRWtdhtcYEx0i0nQmIr53550JrHC3pwJHikgLN2mcACxT1a1AgYiMcPtfLgUmued8BFzmbp8HfOX240wFxolIW3cQwDi3zNTg1e/WsacoeKd/SpV5x5J9ZpU9sW8Hju2ZTqfWSSEvDmaMaf4idcPmYyLSB2f48gbgagBV3S0iTwFzcJq5PlXVye4511A5vHmK+wB4GXhTRFbj1GQudF8rzx1GPcc97gFVzQv3hTV1ZVXmLPvzyX04c1AXTnpyOqXlSorPErZz7xrjt4SvMcYEEpFEo6rn1rDv3zhNZVXL5wIDApQXA+cHea1XgFcOPtLosnPfAV7733q/suxubclMa0FmWgvW5hb6rf2Sbk1jxpgQ2NdRA0BFhfKzx/3v8k9vlUg/t9P/pUuzOWtwF3q2bxWJ8IwxTZjNdWYAeGv2Rr/1Y646vge3n3aE93nP9q3424VHRSI0Y0wTZzUaA8D6nYV+z32TjDHGHAqr0US5/SXlLN+2l5aJ9k/BGBMe9ukSpXbsLSY5IZYHPl7G/83LYWSvdO++n/VOr+FMY4ypG0s0UWr4I9P8nn+3eicAH18/kr6dUyIRkjGmmbJEE2WunTiPbfnFQfcPzGjdgNEYY6KBJZoosnFXEZ8u3hbpMIwxUcZGnUWRxz5bHukQjDFRyBJNlNhRUBy0NvPiJUMbOBpjTDSxprMoUFRSxvCHp1Urv2h4JqcP7MLI3ukcmdGak/tXX2vGGGMOlSWaKLBhV+VKmbExQnmF0qtDKx4950hv+UfXj4xEaMaYKGBNZ1GgoLjMu33R8EwAstq1iFQ4xpgoY4kmCvzixZne7Y4pSYBTszHGmIZgiaYZW751L69Xmfb/rKO60j29JZeMyIpITMaY6GN9NM3UvZOW8PrMDdXKu7ZJ5uubRzV8QMaYqGU1mmYqUJJ59fJhxFiTmTGmgVmNphkqK6+oVvbdrSeS0dYGABhjGp7VaJqhrQHmMuuYmhSBSIwxxhJNs1NcWs76XYXVyuNj7U9tjIkMazprZvre/Zl3++T+HblpXB/yCksiGJExJtpZommmxhzRgRcvyY50GMYYY01nzYnvIIAh3dpGMBJjjKlkiaYZKSwp9263SU6IYCTGGFOp1kQjIteLiH09bgIKD1TOadamRXwEIzHGmEqh1Gg6AXNE5D0ROUVE7I6/Rmr0kzO823ZfpjGmsag10ajqXUBv4GXgcmCViDwiIj3DHJupg+LScvaXVjad9erQKoLRGGNMpZD6aFRVgW3uowxoC7wvIn8JY2ymDk7663Tv9sJ7xtKrQ0rkgjHGGB+1Dm8WkT8AlwE7gX8Bf1bVUhGJAVYBt4Q3RFObVdsL2OLOBvCvS7Np08IGAhhjGo9Q7qNJB85RVb9ZGlW1QkTOCE9YJlTfrdrJxS/PAuCIzqmM6dcxwhEZY4y/UBLN3wBEJM2nrEBVS1V1eViiMiHzJBmAD689NoKRGGNMYKH00cwHcoGfcJrKcoF1IjJfRIaGMzgTupTEOJLiYyMdhjHGVBNKovkMOE1V01W1HXAq8B5wLfBcOIMzNXPGaDh8R5wZY0xjEkqiyVbVqZ4nqvo5cLyq/gAkhi0yU6t1OytnaS6r0BqONMaYyAmljyZPRG4F3nGfXwDsFpFYoPoKW6bBbMgrinQIxhhTq1BqNL8EMoD/uo9MtywW+EW4AjO12+BTo7nuRLt/1hjTONVYo3FrLX9T1YuDHLK6/kMyofps6TY6t07i+1tPIsbmnDHGNFI11mhUtRxoLyJ2B2AjlFtwgCGHtbUkY4xp1ELpo1kPfC8iHwHethpVfSpcQZnaTVq4mTW5hQzLSqv9YGOMiaBQ+mi2AJ+4x6b4PA6aiDwoIj+KyEIR+VxEurjl8SLyuogsFpHlInK7zzlD3fLVIvJ3zyzSIpIoIu+65bNEJMvnnMtEZJX7uOxQYm5sbnhnIWCjzYwxjV+tNRpVvR9ARFqqamFtx4foCVW9233dPwD3AFcD5wOJqjpQRFoAy0TkbVVdDzwPXAn8AHwKnAJMAa4AdqtqLxG5EHgcuMCdyeBeIBtQYJ6IfKSqu+vpGiJm9ro87/aWPfsjGIkxxtQulIXPjhGRZcBy9/kgETmkGzVVda/P05Y4iQD3Z0sRiQOSgRJgr4h0BlJVdaY7k/QbwFnuOeOB193t94HRbm3nZOALVc1zk8sXOMmpyXt3zibv9qDMNpELxBhjQhBK09nfcD60dwGo6iLg+EN9YxF5WEQ2Ab/CqdGAkygKga3ARuCvqpoHdAVyfE7Pcctwf25yYysD8oF2vuUBzqkay5UiMldE5ubm5h7qpYVdSXnl7Us3jj08gpEYY0ztQl2PZlOVolrnOxGRL0VkSYDHePc171TVTGAicL172nD3tbsA3YGbRKQHEGhYlacWFGxfTef4F6pOUNVsVc1u3759bZcWcaVlTqKJjxXiY0P6ExpjTMSEMupsk4gcC6g7zPkPuM1oNVHVMSHG8BYwGac/5ZfAZ6paCuwQke9x+li+xblp1CMDZ5ACODWVTCDHbXJrDeS55aOqnDM9xJgatS35Tr+MDQQwxjQFoXwdvhq4jsrmq8Hu84MmIr19np4JrHC3NwIniaMlMAJYoapbgQIRGeH2v1wKTHLP+QhnYTaA84Cv3H6cqcA4EWkrIm2BcW5Zk7Y4J58fc/IBOPuogC2BxhjTqIQy6mwnTj9KfXpMRPrgzJW2ASeZATwLvAoswWn6elVVf3T3XQO8hjNIYIr7AHgZeFNEVuPUZC50484TkQeBOe5xD7j9PU3a+Ge/A+Ci4Zk8MH5AhKMxxpjahbKUc3vgd0CW7/Gq+puDfVNVPTdI+T6cIc6B9s0Fqn2yqmpxDee8ArxysHE2Rp7WspYJcdY/Y4xpEkLpo5mE00fyJSEMAjDhs3pHgXfbemeMMU1FKImmhareGvZITK0emlw5BuOaUTZbszGmaQil7eUTETkt7JGYWiXFOUs1f3jtsaS3sjXnjDFNQyiJ5gacZLNfRPaKSIGI7K31LFOv7vhwMZ8t3cbIXukcdVjbSIdjjDEhC2XU2SFNoGnqx1uzNgLQrV2LCEdijDF1E7RGIyIX+2wfV2Xf9dXPMOFS6jPlTGaaJRpjTNNSU9PZjT7b/6iy76CHNpu6Ky6tHOzXtkV8BCMxxpi6qynRSJDtQM9NGO33STStky3RGGOalpoSjQbZDvTchNEni7Z6t5MTQhmRbowxjUdNn1p9ReRHnNpLT3cb93mPsEdmvB74ZJl326qSxpimpqZEc0SDRWGC2rG32LvdISWR43qlRzAaY4ypu6CJRlU3NGQgJrB5GypXnr7rjH7ExlidxhjTtNisjI3cPR8tBeCZCwfz8yM7RzgaY4ypO+tZbsQOlJWTW3CAgV1bM36wrT1jjGmarEbTiOXsdlbS/PVxWZENxBhjDkHQGo2ILCbwMGYBVFWPDFtUBoCte5yBAF3aJEc4EmOMOXg1NZ2d0WBRmIBWuevPdLVEY4xpwmzUWSOVs7uI+z9eRtsW8WS0tURjjGm6au2jEZERIjJHRPaJSImIlNsyAeE30Z2teXdRKSI2pNkY03SFMhjgn8BFwCogGfgt1SfZNPWsU2oSAMOybO0ZY0zTFtLwZlVdLSKxqloOvCoi/wtzXFHPszTA8xcPjXAkxhhzaEJJNEUikgAsFJG/AFuBluENK7p9sWw701fmIgJpLRIiHY4xxhySUBLNJUAscD3wJyATODecQUWz8grld2/MBSAlMY4Ym3LGGNPEhbKUs2f02X7g/vCGYxbl7PFut2lpa88YY5q+mm7YfE9VfxHsxk27YTM85vtMotmltQ1rNsY0fTXVaG5wf9qNmw1oW74zG8Dtp/blxL4dIhyNMcYcuqDDm1XVs6zjtaq6wfcBXNsw4UWf7QUH6NauBVed0JPDO6ZEOhxjjDlkodxHMzZA2an1HYhxhjR/81Mu3drZoD5jTPNRUx/NNTg1lx4+yzgDpADfhzuwaPTZkm3k7y9lzBHWZGaMaT5q6qN5C5gCPArc5lNeoKp5YY0qSv3+7QUAtG+VGOFIjDGm/tQ0qWY+kA9cJCKxQEf3+FYi0kpVNzZQjFEnOSE20iEYY0y9qfU+GhG5HrgP2A5UuMUK2PDmMGmfYjUaY0zzEcrMAH8E+qjqrjDHEvX6dU4lf38p/bu0jnQoxhhTb0IZdbYJpwnNhFFFhbKr8ABDutlszcaY5iWUGs1aYLqITAYOeApV9amwRRWF/vHVarbvPcDe/aWRDsUYY+pVKIlmo/tIcB8mDKat2A7A7qKSCEdijDH1K5RJNe8HEJGWqloY/pCiU3ys04r5twsGRzYQY4ypZ6Es5XyMiCwDlrvPB4nIc2GPLMpsyy/m7KO60qN9q0iHYowx9SqUwQB/A04GdgGo6iLg+DDGFHUqKpRte4vp6C7fbIwxzUkoiQZV3VSlqLw+3lxEbhYRFZF0n7LbRWS1iKwUkZN9yoeKyGJ3399FRNzyRBF51y2fJSJZPudcJiKr3Mdl9RFzOLzwzRrKK5TEuJD+HMYY06SENLxZRI4FVEQSRORm3Ga0QyEimTgTdm70KesHXAj0B04BnnNnJQB4HrgS6O0+TnHLrwB2q2ov4Gngcfe10oB7gaOB4cC9ItIoxw5/ucwZCHBsz3YRjsQYY+pfKInmauA6oCuQAwymfpYJeBq4Bf9F1cYD76jqAVVdB6wGhotIZyBVVWeqqgJvAGf5nPO6u/0+MNqt7ZwMfKGqeaq6G/iCyuTUqCzYtIeOqYkc3cMSjTGm+QlleHMfVf2Vb4GIHMchzOAsImcCm1V1kdsC5tEV+MHneY5bVupuVy33nLMJQFXLRCQfaOdbHuCcqvFciVNb4rDDDju4izpIT3/xE6qQkmTLNhtjmqdQEs0/gCEhlPkRkS+BTgF23QncAYwLdFqAMq2h/GDP8S9UnQBMAMjOzg54TDhsyy/mmWmrALjjtL4N9bbGGNOgalqP5hjgWKC9iNzosysVqHV6YVUdE+R1BwLdAU9tJgOYLyLDcWodmT6HZwBb3PKMAOX4nJMjInFAayDPLR9V5ZzptcXdkF6Ysca7PcDmNzPGNFM19dEkAK1wklGKz2MvcN7BvqGqLlbVDqqapapZOAlhiKpuAz4CLnRHknXH6fSf7S4rXSAiI9z+l0uBSe5LfgR4RpSdB3zl9uNMBcaJSFt3EMA4t6zRSEmqzPMdbGizMaaZqmk9mhnADBF5TVU3iEiKU6z7whWMqi4VkfeAZUAZcJ2qeoZSXwO8BiTjLMg2xS1/GXhTRFbj1GQudF8rT0QeBOa4xz3Q2BZsS4p3KoZXjOwe4UiMMSZ8QumjSRGRBUAagIjsBC5T1SX1EYBbq/F9/jDwcIDj5gIDApQXA+cHee1XgFfqI85wKCguIz5WuOv0IyIdijHGhE0ow5snADeqajdV7Qbc5JaZQ7SnqITWyfFUGXlnjDHNSiiJpqWqfu15oqrTgZZhiyiKrNxeQE+b28wY08yFkmjWisjdIpLlPu4C1oU7sOZu0sLNLNi4hwFdbbSZMaZ5CyXR/AZoD3zgPtKBX4czqGhwwzsLARjQNTWygRhjTJjVdB9NEs70M72AxcBNqmrLP9azgVajMcY0czXVaF4HsnGSzKnAEw0SUZTpnm59NMaY5q2m4c39VHUggIi8DMxumJCav9U7CgD4w0m9iI2xEWfGmOatphqNt5lMVcsaIJaocd3EBQDk7Nkf4UiMMSb8aqrRDBKRve62AMnuc8GZIcB6sQ9SgrvAWY90GyVujGn+apqCptaJM83BSW+VgAhcdULPSIdijDFhZ2sHR8DW/GJG9+1IfKz9+o0xzZ990kXA1vxiurSx2ZqNMdHBEk0DKyopI39/KZ1aW6IxxkQHSzQNbGt+MQBdWidHOBJjjGkYlmga2NY9TqKxGo0xJlpYomlgS7fkA9DdhjYbY6KEJZoGtHPfAR6dsoI2LeLpaEs3G2OihCWaBvTSN2sB2FNkc5MaY6KHJZoGNOOn3EiHYIwxDc4STQPaua8EsKlnjDHRxRJNA9p3oJRhWW1596pjIh2KMcY0GEs0DeRAWTnFpRWM6tOB9imJkQ7HGGMajCWaBlJQ7Ky0kJpU04TZxhjT/FiiaSD5+52RZilJ8RGOxBhjGpYlmgayfa8zI0AHazYzxkQZSzQNZP3OIgC6trU5zowx0cUSTQP55McttEyItTnOjDFRxxJNA1m5rYDTBnYmMc4WLjXGRBdLNA1gb3EpuwpL6NWhVaRDMcaYBmeJpgFscPtnurWzGQGMMdHHEk0DWL+rEICs9BYRjsQYYxqeJZoGsMFNNN3SrEZjjIk+lmgawLqdRXRKTSI5wQYCGGOijyWaBrAoZw99OqVEOgxjjIkISzRhNnd9Hqt37OOkvh0iHYoxxkSEJZowm79xNwA/H9QlwpEYY0xkWKIJs9U79pHeKpG0lgmRDsUYYyLCEk2Ybc0vpmsbm3bGGBO9LNGE2bb8YjqmWqIxxkSviCYaEblZRFRE0t3nY0Vknogsdn+e5HPsULd8tYj8XUTELU8UkXfd8lkikuVzzmUissp9XNbQ11dRoWzes58ubWzGZmNM9IpYohGRTGAssNGneCfwc1UdCFwGvOmz73ngSqC3+zjFLb8C2K2qvYCngcfd108D7gWOBoYD94pI27BdUACb9+ynqKScwzva0GZjTPSKZI3maeAWQD0FqrpAVbe4T5cCSW6NpTOQqqozVVWBN4Cz3OPGA6+72+8Do93azsnAF6qap6q7gS+oTE4NYuW2AgC7h8YYE9UikmhE5Exgs6ouquGwc4EFqnoA6Ark+OzLcctwf24CUNUyIB9o51se4Jyq8VwpInNFZG5ubu5BXFFgK7c7iebwjjZrszEmesWF64VF5EugU4BddwJ3AONqOLc/ThOY5xgJcJjWsq+mc/wLVScAEwCys7MDHnMwftpeQNc2yaQkxdfXSxpjTJMTtkSjqmMClYvIQKA7sMjtz88A5ovIcFXdJiIZwIfApaq6xj0txz3OIwPY4rMvE8gRkTigNZDnlo+qcs70Q7+y0P20fZ/VZowxUa/Bm85UdbGqdlDVLFXNwkkIQ9wk0waYDNyuqt/7nLMVKBCREW7/y6XAJHf3RzgDBwDOA75y+3GmAuNEpK07CGCcW9YgVJWNuwptDRpjTNQLW43mIF0P9ALuFpG73bJxqroDuAZ4DUgGprgPgJeBN0VkNU5N5kIAVc0TkQeBOe5xD6hqXoNcBbC7qJTCknIy02wNGmNMdIt4onFrNZ7th4CHghw3FxgQoLwYOD/IOa8Ar9RLoHW0Kc9ZVTOzrd1DY4yJbjYzQJhs2u0mGqvRGGOinCWaMNmUtx+wRGOMMZZowmRjXhFtW8TTKjHirZPGGBNRlmjCZNnWvRxmI86MMcYSTTgUFJeyaNMeTuzTPtKhGGNMxFmiCYO8whIAMtpa/4wxxlgHQj0pKC7l+elrGH1ER2JjnNlv0lra1DPGGGOJpp6UlSvPTV9DeqtE4mKdRNOmhS3fbIwxlmjqSXJCLAD7S8tZsGYPAIfZ0GZjjLE+mvqSGBdDjEBRSRnrdu5jXL+OpLdKjHRYxhgTcZZo6omI0CIhjoLiMjbsKqJXB5u12RhjwBJNvUpOiGXltgLKKpQe7S3RGGMMWKKpVy0TYlmyOR+AHu3tZk1jjAFLNPUqOSGOwpJyAHqmW43GGGPAEk29aumOPANo3cLuoTHGGLBEU68OlFUA8PDZ1ZbNMcaYqGWJph7l7y8FoIc1mxljjJclmnoU784IYAMBjDGmks0MUI9evCSbqUu30SHFbtQ0xhgPSzT1qFeHVvTq0CvSYRhjTKNiTWfGGGPCyhKNMcaYsLJEY4wxJqws0RhjjAkrSzTGGGPCyhKNMcaYsLJEY4wxJqws0RhjjAkrUdVIx9CoiEgusOEQXiId2FlP4TQV0XbN0Xa9YNccLQ7lmrupavtAOyzR1DMRmauq2ZGOoyFF2zVH2/WCXXO0CNc1W9OZMcaYsLJEY4wxJqws0dS/CZEOIAKi7Zqj7XrBrjlahOWarY/GGGNMWFmNxhhjTFhZojHGGBNWlmjqiYicIiIrRWS1iNwW6Xjqi4hkisjXIrJcRJaKyA1ueZqIfCEiq9yfbX3Oud39PawUkZMjF/3BE5FYEVkgIp+4z5v79bYRkfdFZIX7tz4mCq75T+6/6SUi8raIJDW3axaRV0Rkh4gs8Smr8zWKyFARWezu+7uISJ0CUVV7HOIDiAXWAD2ABGAR0C/ScdXTtXUGhrjbKcBPQD/gL8BtbvltwOPudj/3+hOB7u7vJTbS13EQ130j8Bbwifu8uV/v68Bv3e0EoE1zvmagK7AOSHafvwdc3tyuGTgeGAIs8Smr8zUCs4FjAAGmAKfWJQ6r0dSP4cBqVV2rqiXAO8D4CMdUL1R1q6rOd7cLgOU4/0nH43w44f48y90eD7yjqgdUdR2wGuf302SISAZwOvAvn+LmfL2pOB9ILwOoaomq7qEZX7MrDkgWkTigBbCFZnbNqvoNkFeluE7XKCKdgVRVnalO1nnD55yQWKKpH12BTT7Pc9yyZkVEsoCjgFlAR1XdCk4yAjq4hzWH38XfgFuACp+y5ny9PYBc4FW3ufBfItKSZnzNqroZ+CuwEdgK5Kvq5zTja/ZR12vs6m5XLQ+ZJZr6Eai9slmNGxeRVsB/gD+q6t6aDg1Q1mR+FyJyBrBDVeeFekqAsiZzva44nOaV51X1KKAQp0klmCZ/zW6/xHicJqIuQEsRubimUwKUNalrDkGwazzka7dEUz9ygEyf5xk41fBmQUTicZLMRFX9wC3e7lapcX/ucMub+u/iOOBMEVmP0wR6koj8m+Z7veBcQ46qznKfv4+TeJrzNY8B1qlqrqqWAh8Ax9K8r9mjrteY425XLQ+ZJZr6MQfoLSLdRSQBuBD4KMIx1Qt3dMnLwHJVfcpn10fAZe72ZcAkn/ILRSRRRLoDvXE6EpsEVb1dVTNUNQvn7/iVql5MM71eAFXdBmwSkT5u0WhgGc34mnGazEaISAv33/honP7H5nzNHnW6Rrd5rUBERri/q0t9zglNpEdFNJcHcBrOiKw1wJ2Rjqcer2skTjX5R2Ch+zgNaAdMA1a5P9N8zrnT/T2spI6jUxrTAxhF5aizZn29wGBgrvt3/i/QNgqu+X5gBbAEeBNntFWzumbgbZw+qFKcmskVB3ONQLb7e1oD/BN3VplQHzYFjTHGmLCypjNjjDFhZYnGGGNMWFmiMcYYE1aWaIwxxoSVJRpjjDFhZYnGRA0RKReRhT6PrEjH5EtEzndnTv66SnmMO2PuEncG3TnufQ6IyKci0qae3n9fLfvbiMi1Ps+7iMj79fHepnmz4c0maojIPlVtFWSf4Px/qAi0vyGIyGc4M+lWTTQXAecCv1DVCnfSz0JV3V3P7x/09+Puz8K5r2hAfb6vaf6sRmOilohkuTWI54D5QKaIPC8ic911Su73OXa9iDwiIjPd/UNEZKqIrBGRq32O+7Nb4/jR9/wq73uRWzNZIiKPu2X34Nwc+4KIPFHllM7AVk8SVNUcT5Jx40p3r2WFOyHmEhGZKCJjROR7d92R4e7x94nIzT6xLKlasxORViIyTUTmu3F6ZiJ/DOjp1gafcN9ziXtOkoi86h6/QEROdMsvF5EPROQzN46/1PHPZJqDSN+5ag97NNQDKKdydoMPgSycGZpH+ByT5v6MBaYDR7rP1wPXuNtP49xBnwK0x5mEE2AcMAFnEsIY4BPg+CoxdMGZ/qQ9zmSWXwFnufumA9kB4s5w338h8CRwlM++9UC6ey1lwED3vecBr7ixjAf+6x5/H3Czz/lLgCx3e5/7Mw5nWnjc117tvk4W/uuaeJ8DNwGvutt93WtMwlnjZS3Q2n2+AciM9L8FezTsw2o0JprsV9XB7uNst2yDqv7gc8wvRGQ+sADoj7MYlIdn/rrFwCxVLVDVXKDY7ScZ5z4W4NSQ+uLMF+VrGDBdnckcy4CJOGvBBKWqOUAf4HacxDhNREYHOHSdqi5Wp+azFJimqurGm1XTe1QhwCMi8iPwJc6U8B1rOWckzjQuqOoKnIRyuLtvmqrmq2oxzhxq3eoQi2kG4iIdgDERVujZcDvYbwaGqepuEXkN51u4xwH3Z4XPtud5HM4H9KOq+mIN71e3JXBdqnoAZ2XDKSKyHWfhqWlVDqsak2+8nv/rZfg3mften8evcGpcQ1W11J3JOtBxvmq6Lt+4yrHPnahjNRpjKqXiJJ58EekInFrH86cCvxFn7R5EpKuIdKhyzCzgBLdfJRa4CJhR04u6/UFd3O0Y4EicGsPBWI+zBAAiMgRnPZaqWuM0B5a6fS2eGkgBTnNhIN/gJChE5HDgMJyJGY2xbxbGeKjqIhFZgNPstBb4vo7nfy4iRwAznUFs7AMupnK9D1R1q4jcDnyNUwv4VFVrm3K9A/CSiCS6z2fjzKB7MP4DXCoiC3GWt/gpwDETgY9FZC5Ov9AKN/Zd7uCCJTi1q2d9znkOZyDDYpxa0+WqesD9PZgoZ8ObjTHGhJU1nRljjAkrSzTGGGPCyhKNMcaYsLJEY4wxJqws0RhjjAkrSzTGGGPCyhKNMcaYsPp/Y53HeOlsKXIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(potential)\n",
    "plt.xlabel(\"Frame of Simulation\")\n",
    "plt.ylabel(\"Potential Energy\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question: What should I report as the average potential energy of this simulation?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-35841.17738512487"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(potential)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-35179.34066243115"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(potential[500:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the uncertainty in this average? Does this seem reasonable?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.987736602411432"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(potential[500:])/np.sqrt(500-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reporting distributions of non-normal distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYMElEQVR4nO3deXCc9Z3n8fe3W7dsSZbU8iXJ8g3m8IHAYDtgqCQ4ThiWDANxkkmFZIvJwVQmm8nuMtTu1m7tbKVmtyYJFXa5w6aSOJOd4MAw5GAGE4fDh0x8Ymx8W5YPyYcu69Zv/+iWoziS9cju7ufpfj6vKpWO5+nuj5sfHz/+6Xl+jznnEBGR4Ir4HUBERC5PRS0iEnAqahGRgFNRi4gEnIpaRCTgclLxpJWVla6uri4VTy0ikpW2bt3a4pyLjbQtJUVdV1dHQ0NDKp5aRCQrmdmR0bZp6kNEJOBU1CIiAaeiFhEJOBW1iEjAqahFRAJORS0iEnAqahGRgAtMUQ8MOp5Yv58N+5r9jiIiEiiBKepoxHjqNwd47b1TfkcREQmUwBQ1wIyKYo6cveB3DBGRQAlUUddWFHH0TKffMUREAiVQRT2jvIjGc130Dwz6HUVEJDCCVdQVRfQPOk60dvsdRUQkMAJV1DXlRQAc1Ty1iMhFgSrqGRXFABw5o6IWERkSqKKeUlJAXjTCkbP6haKIyJBAFXU0YlSXF3JUR9QiIhcFqqghfuaHpj5ERH7Pc1GbWdTMfmdmr6QyUG15EUfPXsA5l8qXERHJGOM5ov4asCdVQYbUVhTT0dPPuQt9qX4pEZGM4Kmozawa+DjwbGrjxKc+AI7oCkUREcD7EfV3gH8PjHrJoJk9bGYNZtbQ3HzlK+DNqNC51CIiw41Z1Gb2CeC0c27r5fZzzj3tnKt3ztXHYrErDlRz8YhaRS0iAt6OqJcDf2Jmh4GfAHeZ2Q9TFaggN8rkknwVtYhIwphF7Zx71DlX7ZyrAz4FvO6c+2wqQ80oL+aYpj5ERIAAnkcN8eVOdXWiiEjcuIraOfeGc+4TqQozZEZ5EafaeujuG0j1S4mIBF5gj6hBZ36IiEBQi1pnfoiIXBTIov79cqeapxYRCWRRTyrKZWJ+js78EBEhoEVtZokzP1TUIiKBLGqIX0qudalFRAJc1LXlxRw7d4GBQS13KiLhFuCiLqJvwHGitcvvKCIivgpsUV9cRU/THyIScoEt6qFzqXXRi4iEXWCLelpZIblR05kfIhJ6gS3qaMSonqQzP0REAlvUEL+JgFbRE5GwC3RRzygv4sgZ3ZFcRMIt2EVdUUR7dz+tXbojuYiEV6CLWqvoiYgEvKgvrqKnMz9EJMQCXdQ15YUAHNVypyISYoEu6qK8HGITdUdyEQm3QBc1JM780NSHiIRY4Iu6tqJINxAQkVALfFHPKC/mZFu37kguIqEV+KKurSjEOWg8p6NqEQmn4Bd1+dCNblXUIhJOgS/qoXWpVdQiElaBL+qK4jyK86Jal1pEQivwRR2/I3mxilpEQivwRQ1QW17IYV2dKCIhlRFFPX9KCYdbOuns6fc7iohI2mVEUS+uKWPQwY7GVr+jiIikXUYU9cKaMgC2HTvvaw4RET9kRFGXF+cxo6KIbcfO+R1FRCTtMqKoARbVlOmIWkRCKaOK+lRbDydau/yOIiKSVhlV1ADbjp73NYeISLplTFEvmFZCXjSi6Q8RCZ2MKer8nCjXTivhdypqEQmZMYvazArMbLOZbTez3Wb2X9MRbCSLa8rY2dhK/8CgXxFERNLOyxF1D3CXc24hsAhYZWa3pjTVKBbVlNHVN8C+Ux1+vLyIiC/GLGoXN9SMuYkPl9JUo1ikC19EJIQ8zVGbWdTMtgGngdecc5tG2OdhM2sws4bm5uYkx4ybUVHEpKJcXfgiIqHiqaidcwPOuUVANXCLmV0/wj5PO+fqnXP1sVgsyTHjzIyFuvBFREJmXGd9OOfOA28Aq1IRxotFNWV8cLqD9u4+vyKIiKSVl7M+YmZWlvi6EPgw8H6Kc41qUU0ZzsFOraQnIiHh5Yh6KrDezHYAW4jPUb+S2lijG/qFos6nFpGwyBlrB+fcDmBxGrJ4UlaUx8zKYs1Ti0hoZMyVicMNraTnnC9nCYqIpFXGFnVzew9Nrd1+RxERSbmMLWrQSnoiEg4ZWdTXTi0hLyeiC19EJBQysqjzciJcN61Ev1AUkVDIyKKG+PTHzuOt9GklPRHJchld1N19g+w92e53FBGRlMrYol5cMwnQSnoikv0ytqhrygspL85TUYtI1svYojazixe+iIhks4wtaojPUx9o7qBNK+mJSBbL6KJeXBtfSW/rEZ1PLSLZK6OL+ua6cgpyI/xmb2ruKCMiEgQZXdQFuVFum1XB+r2n/Y4iIpIyGV3UAHdeU8WRMxc41NLpdxQRkZTI+KJeOa8KgPXv66haRLJTxhd1bUURs2LFvLFP89Qikp0yvqgB7pxfxcaDZ+jqHfA7iohI0mVFUa+cH6O3f5B3Drb4HUVEJOmyoqhvmVlOYW6U9e9r+kNEsk9WFHV+TpTlc+Kn6ek+iiKSbbKiqAFWzq+i8VwXB5p1mp6IZJcsKuoYAG/o4hcRyTJZU9TVk4qYWzWBN3Q5uYhkmawpaohfpbj50Fk6e/r9jiIikjRZVdQr58XoHRjk7QNn/I4iIpI0WVXU9XXlFOdFtUiTiGSVrCrqvJwIy+dU8sb7Ok1PRLJHVhU1xOepm1q7+eB0h99RRESSIuuKeug0Pa2mJyLZIuuKemppIddMmajT9EQka2RdUUP8KsUth8/SrpveikgWyNKijtE/6Hhrv07TE5HMl5VFfdOMSUwsyOHX7530O4qIyFXLyqLOjUb42PVT+NWuk7qZgIhkvKwsaoD7FlfT2TvAa3tO+R1FROSqZG1RL51ZztTSAta92+h3FBGRqzJmUZtZjZmtN7M9ZrbbzL6WjmBXKxIx7l00nQ0ftNDS0eN3HBGRK+bliLof+IZz7lrgVuCrZrYgtbGS477F0xkYdLyyvcnvKCIiV2zMonbOnXDOvZv4uh3YA0xPdbBkmD9lIgumlrDud8f9jiIicsXGNUdtZnXAYmDTCNseNrMGM2tobg7OVYH3LZ7O9sZWDjRr7Q8RyUyei9rMJgA/A/7KOdd26Xbn3NPOuXrnXH0sFktmxqvyJ4umETF4SUfVIpKhPBW1meUSL+kfOedeTG2k5JpcUsDyOZWs23ZcS5+KSEbyctaHAc8Be5xzf5/6SMn3bxZN59jZLrYeOed3FBGRcfNyRL0c+HPgLjPblvhYneJcSXX39VMoyI3ol4oikpG8nPXxpnPOnHM3OucWJT5eTUe4ZJmQn8NHF0zhlR0n6O0f9DuOiMi4ZO2ViZe6b8l0Wrv6dD9FEck4oSnqD82ppKI4j59r+kNEMkxoijonGuGehdP41z2nae3SDQVEJHOEpqgBPrlkOr0Dg/xi5wm/o4iIeBaqor5heimzYsW8+K6mP0Qkc4SqqM2M+2+qZvPhs+w92e53HBERT0JV1ABrbq6lIDfC828e8juKiIgnoSvqScV5fHJJNeu2HeeM1qkWkQwQuqIGeGhZHb39g/x401G/o4iIjCmURT138kRunxfjBxuP6EpFEQm8UBY1wBeW19Hc3sMrO3T3FxEJttAW9R3zYsypmsBzbx7S8qciEmihLWoz46HldexuamPLYS1/KiLBFdqiBvjk4mpKC3N1qp6IBFqoi7owL8qnl9by6/dOcuzsBb/jiIiMKNRFDfC522YQMeOFtw/7HUVEZEShL+qppYWsvmEq/7DlGO3dWlVPRIIn9EUN8IUVM+no6ecftzb6HUVE5I+oqIFFNWUsqS3j+28dZmBQp+qJSLCoqBO+sGImR89e4FWtVS0iAaOiTvjY9VOZWzWBb7+2j/4BXVYuIsGhok6IRoxvfHQ+B1s6dWMBEQkUFfUwd183mYXVpXznX/bR0z/gdxwREUBF/QfMjG/efQ1Nrd38aKOWQBWRYFBRX2LF3EqWza7gifX76ezp9zuOiIiKeiR/ffd8znT28v23tAaIiPhPRT2CJbWT+PC1k3lqw0HOX+j1O46IhJyKehR/ffc8Onr6efI3B/2OIiIhp6IexTVTSrh34TReePsQp9u6/Y4jIiGmor6Mr39kHv0Dju+t3+93FBEJMRX1ZcyoKObBm2tYu/mo1qsWEd+oqMfwl3fNJWLGt375vt9RRCSkVNRjmFJawFfvnMM/7zjBG3tP+x1HREJIRe3BX9wxi1mxYv7TS7vo6tWl5SKSXipqD/JzovyP+27g2NkuHn/9A7/jiEjIqKg9unVWBX92UzXPbDjI3pPtfscRkRBRUY/Do6uvZWJBDn+zbieDuhOMiKTJmEVtZs+b2Wkz25WOQEFWXpzHYx9fwNYj5/jJlmN+xxGRkPByRP0CsCrFOTLGny6Zzq2zyvnWL/bQ3N7jdxwRCYExi9o5twE4m4YsGcHM+Nv7bqC7b5D//s/v+R1HREJAc9RXYHZsAl9eOZuXtjWxYV+z33FEJMslrajN7GEzazCzhubm7C+vL6+czczKYh77+U7au/v8jiMiWSxpRe2ce9o5V++cq4/FYsl62sAqyI3yd/ffyPFzXTy2bhfO6SwQEUkNTX1chZvryvn6h+fx8vYmftqgs0BEJDW8nJ63FngHmG9mjWb2xdTHyhxfuXMOy2ZX8F9e3s2+U7oQRkSSz8tZH2ucc1Odc7nOuWrn3HPpCJYpohHjOw8uojgvh0d+/K7WAhGRpNPURxJUlRTw9w8uYt+pDv7bK7v9jiMiWUZFnSR3zIvxpTtms3bzMf5pe5PfcUQki6iok+gbH53HktoyHn1xJ0fP6I4wIpIcKuokyo1GeHzNYiIGj6x9l55+zVeLyNVTUSdZ9aQi/u7+hexobOXRn+3U+dUictVU1Cmw6vop/LuPzOPF3x3n26/t8zuOiGS4HL8DZKu/vGsOjecu8Pjr+5k+qZAHb671O5KIZCgVdYoMrbJ3orWbv1m3iymlhdwxL/svrReR5NPURwrlRiP8788sYW7VBL7yw63sbmr1O5KIZCAVdYpNLMjlhYduoaQwly+8sIWm811+RxKRDKOiToMppQV8/6GbudAzwEPf30KblkUVkXFQUafJNVNK+D+fvYkDzR0qaxEZFxV1Gq2YW8njaxaz/dh5PvvsJs5f6PU7kohkABV1mq2+YSpPfvYm3j/Rzqee3khLh26QKyKXp6L2wYcXTOa5z9dz+EwnDz71Dqfauv2OJCIBpqL2yYfmxvi/D93CydZuHnjqHRrPaREnERmZitpHS2dV8MN/u5Rznb088OQ7HG7p9DuSiASQitpni2snsfbhW+nuH+T+J9+h4fBZvyOJSMCoqAPgumml/PQvbmVCfpQ1z2zkR5uO+B1JRAJERR0Qc6om8tJXV7BsdiWPrdvFoy/upLd/0O9YIhIAKuoAKS3K5fnP38yXV85m7eajrHlmI6d1RohI6KmoAyYaMf7Dqmv43qcX815TG/d87022HTvvdywR8ZGKOqA+ceM0fvblZeRGIzzw5Ds8s+EgA4O6W4xIGKmoA2zBtBL+6ZEV3D4vxt++uocHnnqHA80dfscSkTRTUQfcpOI8nvncTXz7wYXsP93B6u/+VkfXIiGjos4AZsZ9i6t57eu386G5lRePrg/q6FokFFTUGaSqpIBnPld/8ej6Y9/9LU+s309334Df0UQkhVTUGWb40fXK+TH+56/2cuf/eoP/13BM0yEiWUpFnaGqSgp46s/r+YeHb6WqpIBv/uMOPv74b1m/9zTOqbBFsomKOsMtnVXBz7+yjCc+vYSuvvitvj7z7CZ2NJ73O5qIJIml4uirvr7eNTQ0JP155fJ6+wf58aYjPP76fs529nLbrAq+uGImd11TRSRifscTkcsws63OufoRt6mos097dx9rNx/lhbcO09TazczKYh5aXsf9N1VTlJfjdzwRGYGKOqT6Bgb55a6TPPfmIbYdO09JQQ5rltbyQH0Ns2MT/I4nIsOoqIWtR87x3JsH+eWukww6uH56CfcunM49C6cxpbTA73gioaeilotOtXXzyo4TvLztONsbWzGDpTPLuXfRdD66YDIVE/L9jigSSipqGdGhlk5e3tbES9uOczBxG7Drp5ewYk6M2+dWclPdJPJzoj6nFAkHFbVclnOO3U1trH//NL/d38K7R87RP+goyI2wdGYFy+dUsLC6jOunl1Kcr19GiqSCilrGpaOnn40HzvDm/hY2fNDMweb40bYZzIlN4MbqMm6sLuWG6lLmVE2gpCDX58Qime+qi9rMVgHfBaLAs865b11ufxV1dmlu72Hn8fPsaGxNfJynpaP34vbKCfnMqixmZmUxM2Pxz9WTCplcUkB5UZ7O4Rbx4KqK2syiwD7gI0AjsAVY45x7b7THqKizm3OOE63d7DreysGWTg41d3KopZODLZ20dPT8wb45EaNqYj5VJQVMLsmnckI+pYW5lBbmUjL0uSCXksIcivKi5OdEKciNUpAboTA3Sk5UF89KOFyuqL1MON4C7HfOHUw82U+Ae4FRi1qym5kxrayQaWWFf7StrbuPwy2dNJ3v4lRbD6faujnV1sPp9m4OtXSy+dBZ2rr7PS8glRMx8nIiRCNGTsSIRiKJz/EPM7BEJgMY9v2I2a/4Ty0ytklFefz0S7cl/Xm9FPV04Niw7xuBpZfuZGYPAw8D1NbWJiWcZJ6SgtzEHHbZqPs457jQO0Bbdx+tXX20dfXT1tVHV98AXX0D9PQN0N03SFffAN19A/QNDNI/6BgYdPHPA4nPg4M4wDkSnx0O4t+M9LqjbRBJklT9vsZLUY90EPJHI9459zTwNMSnPq4yl2QxM6M4P4fi/Bymlv7xUbmI/CEvE4CNQM2w76uBptTEERGRS3kp6i3AXDObaWZ5wKeAl1MbS0REhow59eGc6zezR4BfET8973nn3O6UJxMREcDbHDXOuVeBV1OcRURERqCTVEVEAk5FLSIScCpqEZGAU1GLiARcSlbPM7Nm4MhldqkEWpL+wsmhbFdG2a5MULMFNRdkb7YZzrnYSBtSUtRjMbOG0RYf8ZuyXRlluzJBzRbUXBDObJr6EBEJOBW1iEjA+VXUT/v0ul4o25VRtisT1GxBzQUhzObLHLWIiHinqQ8RkYBTUYuIBFzSi9rMVpnZXjPbb2b/cYTtZmaPJ7bvMLMlXh+bhmyfSWTaYWZvm9nCYdsOm9lOM9tmZkm/IaSHbCvNrDXx+tvM7D97fWyKc31zWKZdZjZgZuWJbal+z543s9NmtmuU7X6OtbGy+TLWPOTyZZx5zObnWKsxs/VmtsfMdpvZ10bYJ3XjzTmXtA/iy6AeAGYBecB2YMEl+6wGfkH8zjG3Apu8PjYN2ZYBkxJff2woW+L7w0BlMt+vcWZbCbxyJY9NZa5L9r8HeD0d71ni+W8HlgC7Rtnuy1jzmM2vsTZWrrSPM6/ZfB5rU4Elia8nEr/hd9q6LdlH1BdvhOuc6wWGboQ73L3AD1zcRqDMzKZ6fGxKsznn3nbOnUt8u5H43WzS4Wr+7Kl838b73GuAtUl67TE55zYAZy+zi19jbcxsfo01D+/ZaHx/zy6R7rF2wjn3buLrdmAP8fvJDpey8Zbsoh7pRriX/mFG28fLY1OdbbgvEv/bcYgDfm1mWy1+I99k8prtNjPbbma/MLPrxvnYVObCzIqAVcDPhv04le+ZF36NtfFK51jzIt3jbFz8HmtmVgcsBjZdsill483TjQPGwcuNcEfbx9NNdK+C5+c3szuJ/8+zYtiPlzvnmsysCnjNzN5PHAGkK9u7xNcC6DCz1cDPgbkeH5vKXEPuAd5yzg0/Ikrle+aFX2PNMx/G2lj8GGfj5dtYM7MJxP+C+CvnXNulm0d4SFLGW7KPqL3cCHe0fVJ9E11Pz29mNwLPAvc6584M/dw515T4fBpYR/yfM2nL5pxrc851JL5+Fcg1s0ovj01lrmE+xSX/FE3xe+aFX2PNE5/G2mX5NM7Gy5exZma5xEv6R865F0fYJXXjLckT7jnAQWAmv580v+6SfT7OH064b/b62DRkqwX2A8su+XkxMHHY128Dq9KcbQq/v0DpFuBo4j1M2fvm9bmBUuJzi8Xpes+GvU4do/9izJex5jGbL2PNQ660jzOv2fwca4n34AfAdy6zT8rGW1Lf5ESo1cR/I3oAeCzxsy8BXxr2B34isX0nUH+5x6Y527PAOWBb4qMh8fNZiTd3O7Dbp2yPJF57O/FfPi273GPTlSvx/eeBn1zyuHS8Z2uBE0Af8aOWLwZorI2VzZex5iGXL+PMSzafx9oK4tMVO4b9N1udrvGmS8hFRAJOVyaKiAScilpEJOBU1CIiAaeiFhEJOBW1iEjAqahFRAJORS0iEnD/Hzw3dj3KCQXIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "betaf = scipy.stats.beta(0.8,4)\n",
    "x = np.linspace(0,2)\n",
    "plt.plot(x,betaf.pdf(x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD6CAYAAACIyQ0UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAMtUlEQVR4nO3da4zcVR3G8eehpZGbQuhKSAEXjSJKRHAEBUO4iLbUQEx4Ad4iMdmYKIHERKovvL4pbwwYFbIpiEaEKBeDICgREQkWmEK5FgxglYLaqWi4JZLCzxcz0HZ3tnNmO////Hbm+0km7Oz8O316sn16OHPOjCNCAIC8dht2AADAzlHUAJAcRQ0AyVHUAJAcRQ0AyVHUAJBcz6K2fZjt9dvdnrd9fg3ZAACS3M8+atuLJD0j6diI+Ntc1y1dujQmJyd3PR0AjIl169ZtiYiJbo8t7vO5TpH05M5KWpImJyfVbDb7fGoAGF+25+zVfteoz5J01Ry/yZTtpu1mq9Xq82kBAHMpLmrbSySdLumX3R6PiOmIaEREY2Ki6+wdADAP/cyoV0i6LyL+VVUYAMBs/RT12Zpj2QMAUJ2iora9p6RTJV1XbRwAwExFuz4i4mVJ+1ecBQDQBScTASA5ihoAkqOoASC5fk8mVm5y1U1dv79x9cqakwBADsyoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkisqatv72r7G9mO2N9j+cNXBAABtpZ+ZeLGkWyLiTNtLJO1ZYSYAwHZ6FrXtN0s6QdLnJSkiXpH0SrWxAACvK1n6eLuklqQf277f9hrbe828yPaU7abtZqvVGnhQABhXJUW9WNLRki6JiKMkvSRp1cyLImI6IhoR0ZiYmBhwTAAYXyVFvUnSpoi4u3P/GrWLGwBQg55FHRH/lPS07cM63zpF0qOVpgIAvKF018e5kq7s7Ph4StI51UUCAGyvqKgjYr2kRrVRAADdcDIRAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgOYoaAJKjqAEgucUlF9neKOkFSa9K2hoRjSpDAQC2KSrqjpMiYktlSQAAXbH0AQDJlRZ1SPqd7XW2p7pdYHvKdtN2s9VqDS4hAIy50qI+PiKOlrRC0pdsnzDzgoiYjohGRDQmJiYGGhIAxllRUUfEs53/bpZ0vaRjqgwFANim54uJtveStFtEvND5+mOSvlN5shkmV93U9fsbV6+sOQkA1Ktk18cBkq63/fr1P4+IWypNBQB4Q8+ijoinJB1ZQxYAQBdszwOA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5ChqAEiOogaA5Hp+Cnl2k6tu6vr9jatX1pwEAKpRPKO2vcj2/bZvrDIQAGBH/Sx9nCdpQ1VBAADdFRW17YMkrZS0pto4AICZSmfUF0n6qqTX5rrA9pTtpu1mq9UaRDYAgAqK2vYnJG2OiHU7uy4ipiOiERGNiYmJgQUEgHFXMqM+XtLptjdKulrSybZ/VmkqAMAbehZ1RHwtIg6KiElJZ0m6LSI+U3kyAIAkDrwAQHp9HXiJiNsl3V5JEgBAVwv+ZOJcOLEIYFSw9AEAyVHUAJAcRQ0AyVHUAJAcRQ0AyVHUAJAcRQ0AyY3sPuq5sL8awELDjBoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASC5nkVt+02277H9gO1HbH+7jmAAgLaSd8/7n6STI+JF27tLutP2zRGxtuJsAAAVFHVEhKQXO3d379yiylAAgG2K1qhtL7K9XtJmSbdGxN1drpmy3bTdbLVaA44JAOOr6IMDIuJVSe+3va+k620fEREPz7hmWtK0JDUajQU34+YDBQBk1deuj4j4r6TbJS2vIgwAYLaSXR8TnZm0bO8h6aOSHqs4FwCgo2Tp40BJP7G9SO1i/0VE3FhtrDxYEgEwbCW7Ph6UdFQNWQAAXXAyEQCSo6gBIDmKGgCSo6gBIDmKGgCSo6gBIDmKGgCSo6gBIDmKGgCSo6gBILmitznFbLwHCIC6MKMGgOQoagBIjqIGgOQoagBIjqIGgOQoagBIjqIGgOQoagBIjgMvAzbXQRiJwzAA5ocZNQAkR1EDQHI9i9r2wbb/YHuD7Udsn1dHMABAW8ka9VZJX4mI+2zvI2md7Vsj4tGKs40c3sgJwHz0nFFHxD8i4r7O1y9I2iBpWdXBAABtfa1R256UdJSku7s8NmW7abvZarUGFA8AULw9z/bekq6VdH5EPD/z8YiYljQtSY1GIwaWcIyxVAJAKixq27urXdJXRsR11UYaPzvbew0AJbs+LOkySRsi4nvVRwIAbK9kjfp4SZ+VdLLt9Z3baRXnAgB09Fz6iIg7JbmGLACALjiZCADJUdQAkBxFDQDJUdQAkBxFDQDJ8cEBC1C/B2Q4yQgsbMyoASA5ihoAkqOoASA5ihoAkqOoASA5ihoAkqOoASA59lGPMT5BBlgYmFEDQHLMqDELM20gF4p6DPCZjMDCxtIHACRHUQNAchQ1ACRHUQNAchQ1ACTXs6htX257s+2H6wgEANhRyYz6CknLK84BAJhDz6KOiDskPVdDFgBAFwNbo7Y9Zbtpu9lqtQb1tAAw9gZ2MjEipiVNS1Kj0YhBPS/y48g5UC12fQBAchQ1ACTXc+nD9lWSTpS01PYmSd+MiMuqDoZ8hvnmTiyvYJz1LOqIOLuOIBg9Oyt2ChYox9IHACRHUQNAcnxwAIaCNWegHDNqAEiOGTVS4WPDgNmYUQNAcsyoMZJYA8cooaixoLFUgnHA0gcAJEdRA0ByFDUAJMcaNcbKoNa0eVESdaKogQFitwmqwNIHACTHjBqYB7YFok4UNVADlkSwK1j6AIDkKGoASI6lD2CI5rPWPddyCcsro4uiBhaYfsu93+sp9nxY+gCA5JhRA9gBSyj5FBW17eWSLpa0SNKaiFhdaSoA6XD8fnh6FrXtRZJ+KOlUSZsk3Wv7hoh4tOpwAEbPIF9AHRclM+pjJD0REU9Jku2rJZ0hiaIGUIuFchK0qn9QSop6maSnt7u/SdKxMy+yPSVpqnP3RduPzzPTUklb5vlrRxnj0h3jMhtj0l3l4+ILd+mXv22uB0qK2l2+F7O+ETEtabqPUN1/M7sZEY1dfZ5Rw7h0x7jMxph0t5DHpWR73iZJB293/yBJz1YTBwAwU0lR3yvpnbYPtb1E0lmSbqg2FgDgdT2XPiJiq+0vS/qt2tvzLo+IRyrMtMvLJyOKcemOcZmNMeluwY6LI2YtNwMAEuEIOQAkR1EDQHJDKWrby20/bvsJ26u6PG7b3+88/qDto4eRs24F4/Lpzng8aPsu20cOI2fdeo3Ldtd90Parts+sM9+wlIyL7RNtr7f9iO0/1p1xGAr+Hr3F9q9tP9AZl3OGkbMvEVHrTe0XJJ+U9HZJSyQ9IOk9M645TdLNau/h/pCku+vOmXRcjpO0X+frFYzLrOtuk/QbSWcOO3eGcZG0r9oniA/p3H/rsHMnGZevS7qw8/WEpOckLRl29p3dhjGjfuNIekS8Iun1I+nbO0PST6NtraR9bR9Yd9Ca9RyXiLgrIv7TubtW7T3to67k50WSzpV0raTNdYYbopJx+ZSk6yLi75IUEeMwNiXjEpL2sW1Je6td1FvrjdmfYRR1tyPpy+Zxzajp98/8BbX/r2PU9RwX28skfVLSpTXmGraSn5d3SdrP9u2219n+XG3phqdkXH4g6XC1D+49JOm8iHitnnjzM4z3oy45kl50bH3EFP+ZbZ+kdlF/pNJEOZSMy0WSLoiIV9uTpLFQMi6LJX1A0imS9pD0Z9trI+IvVYcbopJx+bik9ZJOlvQOSbfa/lNEPF9xtnkbRlGXHEkfx2PrRX9m2++TtEbSioj4d03ZhqlkXBqSru6U9FJJp9neGhG/qiXhcJT+PdoSES9Jesn2HZKOlDTKRV0yLudIWh3tReonbP9V0rsl3VNPxHkYwmL/YklPSTpU2xb73zvjmpXa8cXEe4a9mJ9kXA6R9ISk44adN9O4zLj+Co3Hi4klPy+HS/p959o9JT0s6YhhZ08wLpdI+lbn6wMkPSNp6bCz7+xW+4w65jiSbvuLnccvVfuV+9PULqWX1f4XcKQVjss3JO0v6Ued2ePWWKDvBlaqcFzGTsm4RMQG27dIelDSa2p/OtPDw0tdvcKfl+9KusL2Q2pPBi+IiNRvC8sRcgBIjpOJAJAcRQ0AyVHUAJAcRQ0AyVHUAJAcRQ0AyVHUAJDc/wETwaG7SiqOFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num = 10000\n",
    "samples = betaf.rvs(num)\n",
    "plt.hist(samples,bins=50,density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pause and do some calculations\n",
    "\n",
    "* What would you report as the mean and standard deviation of this distribution?\n",
    "* What might be a better way to report the behavior of this distribution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean = 0.1680\n"
     ]
    }
   ],
   "source": [
    "mean = np.mean(samples)\n",
    "print(f\"mean = {mean:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "std = 0.1560\n"
     ]
    }
   ],
   "source": [
    "std = np.std(samples)\n",
    "print(f\"std = {std:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.14402158914976904\n",
      "0.47998110650128223\n"
     ]
    }
   ],
   "source": [
    "# what is the 95% confidence interval of this distribution using the formula?\n",
    "print(mean - 2*std)\n",
    "print(mean + 2*std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 2.5% percentile is : 0.0025\n",
      "97.5% percentile is : 0.5743\n"
     ]
    }
   ],
   "source": [
    "xlow = np.percentile(samples,2.5)\n",
    "xhigh = np.percentile(samples,97.5)\n",
    "print(f\" 2.5% percentile is : {xlow:.4f}\")\n",
    "print(f\"97.5% percentile is : {xhigh:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What about the error in the mean of this distribution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distribution_of_means(num,nrepeats):\n",
    "    repeats = np.zeros(nrepeats)\n",
    "    for i in range(nrepeats):\n",
    "        repeats[i] = np.mean(betaf.rvs(num))\n",
    "    return repeats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPMklEQVR4nO3df6zddX3H8edrVMn8sQlrwa5Uy1zNLIuiu3ZOlwXDMhBiKpluZYuSjaRuwUUTXSz+oSZLky6ZuiwZmqpEljhYE0HIQCdDN+L8gRdSsaUyO+mgtqEXMUPdwmx574/77Thebnu+955777m3n+cjubnnfM7ne76vfnPv6377Pd/zPakqJEmnv58ZdwBJ0tKw8CWpERa+JDXCwpekRlj4ktSIVeMOALB69erasGHDuGNI0opy7733PlZVa/rOXxaFv2HDBiYnJ8cdQ5JWlCT/OZf5Qw/pJFmf5EtJ9ifZl+Sd3fgHk3wvyZ7u67KBZa5NciDJg0kumfs/Q5K00Prs4R8D3l1V9yV5PnBvkju7xz5SVX81ODnJJmArcAHwi8A/J3lpVR1fyOCSpLkZuodfVUeq6r7u9g+B/cC6UyyyBbipqp6sqoeAA8DmhQgrSZq/OZ2lk2QD8Erg693QO5Lcn+T6JGd1Y+uARwYWO8QsfyCSbEsymWRyampq7sklSXPSu/CTPA/4DPCuqnoC+CjwEuBC4AjwoRNTZ1n8GRfsqapdVTVRVRNr1vR+kVmSNE+9Cj/Js5gu+09X1c0AVfVoVR2vqqeAj/P0YZtDwPqBxc8DDi9cZEnSfPQ5SyfAJ4H9VfXhgfG1A9OuAPZ2t28DtiY5M8n5wEbgnoWLLEmajz5n6bwOeCvwrSR7urH3AVcmuZDpwzUHgbcDVNW+JLuBB5g+w+caz9CRpPEbWvhV9WVmPy5/xymW2QHsGCGXJGmBLYt32q5UG7bfPpb1Htx5+VjWK2ll8+JpktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNWDXuAJq7DdtvH9u6D+68fGzrljQa9/AlqREWviQ1YmjhJ1mf5EtJ9ifZl+Sd3fjZSe5M8p3u+1kDy1yb5ECSB5Ncspj/AElSP3328I8B766qlwGvAa5JsgnYDtxVVRuBu7r7dI9tBS4ALgWuS3LGYoSXJPU3tPCr6khV3dfd/iGwH1gHbAFu6KbdALypu70FuKmqnqyqh4ADwOYFzi1JmqM5HcNPsgF4JfB14NyqOgLTfxSAc7pp64BHBhY71I3NfK5tSSaTTE5NTc0juiRpLnoXfpLnAZ8B3lVVT5xq6ixj9YyBql1VNVFVE2vWrOkbQ5I0T70KP8mzmC77T1fVzd3wo0nWdo+vBY5244eA9QOLnwccXpi4kqT56nOWToBPAvur6sMDD90GXNXdvgq4dWB8a5Izk5wPbATuWbjIkqT56PNO29cBbwW+lWRPN/Y+YCewO8nVwMPAWwCqal+S3cADTJ/hc01VHV/o4JKkuRla+FX1ZWY/Lg9w8UmW2QHsGCGXJGmB+U5bSWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiNOi8+0HednvErSSuEeviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjTgtzsPX0hnXex4O7rx8LOuVTifu4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjhhZ+kuuTHE2yd2Dsg0m+l2RP93XZwGPXJjmQ5MEklyxWcEnS3PTZw/8UcOks4x+pqgu7rzsAkmwCtgIXdMtcl+SMhQorSZq/oYVfVXcDj/d8vi3ATVX1ZFU9BBwANo+QT5K0QEY5hv+OJPd3h3zO6sbWAY8MzDnUjT1Dkm1JJpNMTk1NjRBDktTHfAv/o8BLgAuBI8CHuvHMMrdme4Kq2lVVE1U1sWbNmnnGkCT1Na/Cr6pHq+p4VT0FfJynD9scAtYPTD0PODxaREnSQphX4SdZO3D3CuDEGTy3AVuTnJnkfGAjcM9oESVJC2HVsAlJbgQuAlYnOQR8ALgoyYVMH645CLwdoKr2JdkNPAAcA66pquOLklySNCdDC7+qrpxl+JOnmL8D2DFKKEnSwvOdtpLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1IihhZ/k+iRHk+wdGDs7yZ1JvtN9P2vgsWuTHEjyYJJLFiu4JGlu+uzhfwq4dMbYduCuqtoI3NXdJ8kmYCtwQbfMdUnOWLC0kqR5G1r4VXU38PiM4S3ADd3tG4A3DYzfVFVPVtVDwAFg88JElSSNYr7H8M+tqiMA3fdzuvF1wCMD8w51Y5KkMVvoF20zy1jNOjHZlmQyyeTU1NQCx5AkzbRqnss9mmRtVR1JshY42o0fAtYPzDsPODzbE1TVLmAXwMTExKx/FKQTNmy/fWzrPrjz8rGtW1pI893Dvw24qrt9FXDrwPjWJGcmOR/YCNwzWkRJ0kIYuoef5EbgImB1kkPAB4CdwO4kVwMPA28BqKp9SXYDDwDHgGuq6vgiZZckzcHQwq+qK0/y0MUnmb8D2DFKKEnSwvOdtpLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUCAtfkhph4UtSIyx8SWqEhS9JjbDwJakRFr4kNcLCl6RGWPiS1AgLX5IasWqUhZMcBH4IHAeOVdVEkrOBfwA2AAeB36uqH4wWU5I0qpEKv/P6qnps4P524K6q2plke3f/vQuwHmksNmy/fSzrPbjz8rGsV6evxTikswW4obt9A/CmRViHJGmORi38Ar6Q5N4k27qxc6vqCED3/ZzZFkyyLclkksmpqakRY0iShhn1kM7rqupwknOAO5N8u++CVbUL2AUwMTFRI+aQJA0x0h5+VR3uvh8FbgE2A48mWQvQfT86akhJ0ujmXfhJnpvk+SduA78D7AVuA67qpl0F3DpqSEnS6EY5pHMucEuSE8/z91X1+STfAHYnuRp4GHjL6DElSaOad+FX1XeBV8wy/n3g4lFCSZIWnu+0laRGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJasRCfIi5pEUwrg9PBz9A/XTlHr4kNcLCl6RGWPiS1AgLX5IaYeFLUiMsfElqhIUvSY2w8CWpERa+JDXCwpekRlj4ktQIC1+SGmHhS1IjLHxJaoSXR5b0DOO6NLOXZV5c7uFLUiMsfElqhIUvSY2w8CWpERa+JDVi0Qo/yaVJHkxyIMn2xVqPJKmfRSn8JGcAfwu8AdgEXJlk02KsS5LUz2Kdh78ZOFBV3wVIchOwBXhgkdYn6TQwrvP/x2kp33uwWIW/Dnhk4P4h4NcHJyTZBmzr7v4oyYPzXNdq4LF5LjtOKzG3mZfGSswMKzP32DPnL+e8yGDmF89lwcUq/MwyVj91p2oXsGvkFSWTVTUx6vMstZWY28xLYyVmhpWZu7XMi/Wi7SFg/cD984DDi7QuSVIPi1X43wA2Jjk/ybOBrcBti7QuSVIPi3JIp6qOJXkH8E/AGcD1VbVvMdbFAhwWGpOVmNvMS2MlZoaVmbupzKmq4bMkSSue77SVpEZY+JLUiBVT+MMu1ZBpf9M9fn+SV40j54xMwzL/SpKvJnkyyXvGkXGmHpn/sNu+9yf5SpJXjCPnTD1yb+ky70kymeQ3x5FzRqZelx9J8uokx5O8eSnznSTLsO18UZL/6rbzniTvH0fOGZmGbucu954k+5L861JnnE2Pbf3nA9t5b/czcvYpn7Sqlv0X0y/8/gfwS8CzgW8Cm2bMuQz4HNPvAXgN8PUVkPkc4NXADuA9K2Q7vxY4q7v9hnFv5znkfh5Pv2b1cuDbyz3zwLwvAncAb17umYGLgH8c98/EHDO/gOmrALyou3/OSsg9Y/4bgS8Oe96Vsof//5dqqKr/BU5cqmHQFuDvatrXgBckWbvUQQcMzVxVR6vqG8BPxhFwFn0yf6WqftDd/RrT77EYtz65f1TdbwbwXGa8EXAM+vxMA/wZ8Bng6FKGO4m+mZeTPpn/ALi5qh6G6d/LJc44m7lu6yuBG4c96Uop/Nku1bBuHnOW0nLL08dcM1/N9P+qxq1X7iRXJPk2cDvwx0uU7WSGZk6yDrgC+NgS5jqVvj8fv5Hkm0k+l+SCpYl2Un0yvxQ4K8m/JLk3yduWLN3J9f5dTPIc4FKmdwxOaaV8iPnQSzX0nLOUlluePnpnTvJ6pgt/7MfC6Zm7qm4BbknyW8BfAL+92MFOoU/mvwbeW1XHk9mmL7k+me8DXlxVP0pyGfBZYONiBzuFPplXAb8GXAz8LPDVJF+rqn9f7HCnMJf+eCPwb1X1+LAnXSmF3+dSDcvtcg7LLU8fvTIneTnwCeANVfX9Jcp2KnPa1lV1d5KXJFldVeO6cFafzBPATV3ZrwYuS3Ksqj67JAmfaWjmqnpi4PYdSa5bAdv5EPBYVf0Y+HGSu4FXAOMs/Ln8TG+lx+EcYMW8aLsK+C5wPk+/gHHBjDmX89Mv2t6z3DMPzP0gy+NF2z7b+UXAAeC14847x9y/zNMv2r4K+N6J+8s184z5n2L8L9r22c4vHNjOm4GHl/t2Bl4G3NXNfQ6wF/jV5b6tu3k/DzwOPLfP866IPfw6yaUakvxJ9/jHmD6L4TKmy+i/gT8aV94u09DMSV4ITAI/BzyV5F1MvxL/xMmed9yZgfcDvwBc1+15HqsxX22wZ+7fBd6W5CfA/wC/X91vzDLOvKz0zPxm4E+THGN6O29d7tu5qvYn+TxwP/AU8Imq2juuzF2uvj8fVwBfqOn/nQzlpRUkqREr5SwdSdKILHxJaoSFL0mNsPAlqREWviQ1wsKXpEZY+JLUiP8DCQSJXCspwZYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(distribution_of_means(2,1000))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAN70lEQVR4nO3dXYwd5X3H8e+vELgIqUrqhbrG7brUVQpS4lQbF5VeEKEWgi8MamiNKopaJKcqVEFKpJhcNJEqS67UJFWlksgpKK6UxrUKFEumL5RGQlEaYEGUYF4aN7iwsYU3LypErWhs/r3YoT4suz5n9+zZYz/+fqSjM/PMMzP/82j129k5M7OpKiRJbfmxcRcgSVp5hrskNchwl6QGGe6S1CDDXZIadO64CwBYs2ZNTU5OjrsMSTqjPPHEE9+tqomFlp0W4T45Ocn09PS4y5CkM0qS/1xsmadlJKlBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQafFHapamskdB8a278O7toxt35IG55G7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1qG+4J1mf5KtJnktyMMlHu/ZPJ/lOkqe613U969yZ5FCSF5JcM8oPIEl6u0EeP3Ac+FhVPZnkXcATSR7qln2uqv60t3OSy4BtwOXATwP/nOQXqurEShYuSVpc3yP3qjpaVU92068BzwHrTrHKVmBvVb1eVS8Ch4DNK1GsJGkwSzrnnmQSeD/waNd0e5Knk9yT5MKubR3wcs9qM5z6l4EkaYUNHO5JLgDuBe6oqleBzwOXApuAo8Bn3uy6wOq1wPa2J5lOMj07O7vUuiVJpzBQuCd5B3PB/uWqug+gql6pqhNV9QbwRU6eepkB1vesfglwZP42q2p3VU1V1dTExMQwn0GSNM8gV8sEuBt4rqo+29O+tqfbDcAz3fR+YFuS85NsADYCj61cyZKkfga5WuZK4Gbgm0me6to+CdyUZBNzp1wOAx8BqKqDSfYBzzJ3pc1tXikjSaurb7hX1ddY+Dz6g6dYZyewc4i6JElD8A5VSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVokKdCSv9vcseBsez38K4tY9mvdKbyyF2SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa1Dfck6xP8tUkzyU5mOSjXfu7kzyU5Fvd+4U969yZ5FCSF5JcM8oPIEl6u0GO3I8DH6uqXwSuAG5LchmwA3i4qjYCD3fzdMu2AZcD1wJ3JTlnFMVLkhbWN9yr6mhVPdlNvwY8B6wDtgJ7um57gOu76a3A3qp6vapeBA4Bm1e4bknSKSzpnHuSSeD9wKPAxVV1FOZ+AQAXdd3WAS/3rDbTtc3f1vYk00mmZ2dnl1G6JGkxA4d7kguAe4E7qurVU3VdoK3e1lC1u6qmqmpqYmJi0DIkSQMYKNyTvIO5YP9yVd3XNb+SZG23fC1wrGufAdb3rH4JcGRlypUkDWKQq2UC3A08V1Wf7Vm0H7ilm74FeKCnfVuS85NsADYCj61cyZKkfs4doM+VwM3AN5M81bV9EtgF7EtyK/AScCNAVR1Msg94lrkrbW6rqhMrXbgkaXF9w72qvsbC59EBrl5knZ3AziHqkiQNwTtUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTp33AVIg5jccWBs+z68a8vY9i0tl0fuktQgj9yHMM6jSUk6FY/cJalBhrskNchwl6QG9Q33JPckOZbkmZ62Tyf5TpKnutd1PcvuTHIoyQtJrhlV4ZKkxQ1y5P4l4NoF2j9XVZu614MASS4DtgGXd+vcleSclSpWkjSYvuFeVY8A3x9we1uBvVX1elW9CBwCNg9RnyRpGYY55357kqe70zYXdm3rgJd7+sx0bW+TZHuS6STTs7OzQ5QhSZpvueH+eeBSYBNwFPhM154F+tZCG6iq3VU1VVVTExMTyyxDkrSQZYV7Vb1SVSeq6g3gi5w89TIDrO/peglwZLgSJUlLtaxwT7K2Z/YG4M0rafYD25Kcn2QDsBF4bLgSJUlL1ffxA0m+AlwFrEkyA3wKuCrJJuZOuRwGPgJQVQeT7AOeBY4Dt1XViZFULklaVN9wr6qbFmi++xT9dwI7hylKkjQc71CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQ33BPck+SY0me6Wl7d5KHknyre7+wZ9mdSQ4leSHJNaMqXJK0uEGO3L8EXDuvbQfwcFVtBB7u5klyGbANuLxb564k56xYtZKkgfQN96p6BPj+vOatwJ5ueg9wfU/73qp6vapeBA4Bm1emVEnSoJZ7zv3iqjoK0L1f1LWvA17u6TfTtUmSVtFKf6GaBdpqwY7J9iTTSaZnZ2dXuAxJOrstN9xfSbIWoHs/1rXPAOt7+l0CHFloA1W1u6qmqmpqYmJimWVIkhay3HDfD9zSTd8CPNDTvi3J+Uk2ABuBx4YrUZK0VOf265DkK8BVwJokM8CngF3AviS3Ai8BNwJU1cEk+4BngePAbVV1YkS1S5IW0Tfcq+qmRRZdvUj/ncDOYYqSJA3HO1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBfe9Qlc52kzsOjGW/h3dtGct+1QaP3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KCh/kF2ksPAa8AJ4HhVTSV5N/A3wCRwGPjNqvrBcGVKkpZiJY7cP1hVm6pqqpvfATxcVRuBh7t5SdIqGsVpma3Anm56D3D9CPYhSTqFYcO9gH9K8kSS7V3bxVV1FKB7v2ihFZNsTzKdZHp2dnbIMiRJvYY65w5cWVVHklwEPJTk+UFXrKrdwG6AqampGrIOSVKPoY7cq+pI934MuB/YDLySZC1A935s2CIlSUuz7HBP8s4k73pzGvh14BlgP3BL1+0W4IFhi5QkLc0wp2UuBu5P8uZ2/rqq/iHJ48C+JLcCLwE3Dl+mJGkplh3uVfVt4H0LtH8PuHqYoiTB5I4DY9v34V1bxrZvrQzvUJWkBhnuktSgYS+FPC2M889XSTodeeQuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBTfybPUkra1z/uvLwri1j2W+LPHKXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWhk4Z7k2iQvJDmUZMeo9iNJeruR3KGa5BzgL4BfA2aAx5Psr6pnR7E/SW0Y152x4zSqu3JHdeS+GThUVd+uqv8F9gJbR7QvSdI8o3q2zDrg5Z75GeCXezsk2Q5s72Z/mOSFedtYA3x3RPWdiRyPkxyLt3I8TjrjxiJ/MtTqP7vYglGFexZoq7fMVO0Gdi+6gWS6qqZWurAzleNxkmPxVo7HSY7FSaM6LTMDrO+ZvwQ4MqJ9SZLmGVW4Pw5sTLIhyXnANmD/iPYlSZpnJKdlqup4ktuBfwTOAe6pqoNL3Myip2zOUo7HSY7FWzkeJzkWnVRV/16SpDOKd6hKUoMMd0lq0NjDvd9jCjLnz7vlTyf5pXHUuVoGGI/3JPnXJK8n+fg4alwtA4zFb3c/E08n+XqS942jztUywHhs7cbiqSTTSX51HHWuhkEfb5LkA0lOJPnwatZ3Wqiqsb2Y+7L1P4CfA84D/g24bF6f64C/Z+7a+SuAR8dZ82kwHhcBHwB2Ah8fd81jHotfAS7spj/kzwYXcPJ7tPcCz4+77nGNRU+/fwEeBD487rpX+zXuI/dBHlOwFfirmvMN4CeSrF3tQldJ3/GoqmNV9Tjwo3EUuIoGGYuvV9UPutlvMHc/RasGGY8fVpdqwDuZd+NgQwZ9vMkfAvcCx1azuNPFuMN9occUrFtGn1acTZ+1n6WOxa3M/YXXqoHGI8kNSZ4HDgC/t0q1rba+Y5FkHXAD8IVVrOu0Mu5w7/uYggH7tOJs+qz9DDwWST7IXLh/YqQVjddA41FV91fVe4DrgT8edVFjMshY/Bnwiao6MfpyTk+jerbMoAZ5TMHZ9CiDs+mz9jPQWCR5L/CXwIeq6nurVNs4LOlno6oeSXJpkjVVdUY9SGsAg4zFFLA3Ccw9TOy6JMer6u9WpcLTwLiP3Ad5TMF+4He6q2auAP6rqo6udqGrxMc2nNR3LJL8DHAfcHNV/fsYalxNg4zHz6dLs+6qsvOAFn/h9R2LqtpQVZNVNQn8LfAHZ1Oww5iP3GuRxxQk+f1u+ReY+6b7OuAQ8N/A746r3lEbZDyS/BQwDfw48EaSO5i7UuDVcdU9CgP+bPwR8JPAXV2mHa9Gnwg44Hj8BnMHQj8C/gf4rZ4vWJsx4Fic9Xz8gCQ1aNynZSRJI2C4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAb9HxiCjZT1CceLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(distribution_of_means(5,1000))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOwklEQVR4nO3dXYxcd3nH8e+vCaS8iqA4qbHdbkDuiyMVB21T2lRVaKomJFUNF1SOWpRKqQxSkKAFtQ69gF5YciVe2ouCZJIUqwWCBaFEhLakKRJCtAmbNA1xTMAQN1nsxstbk/YiEPP0Yo9has96ZufFs/vn+5FGc+Z/zn/O83jt3549e+Y4VYUkqS0/MesCJEmTZ7hLUoMMd0lqkOEuSQ0y3CWpQefOugCACy64oObm5mZdhiStK/fdd983q2pDv3VrItzn5uZYWFiYdRmStK4k+c+V1nlaRpIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGrQmPqGq9WNu950z2e+RvdfOZL/SeuWRuyQ1yHCXpAYZ7pLUIMNdkho0MNyTbEny2SSHkhxM8uZu/J1JvpHkge5xTc+cm5IcTvJIkqum2YAk6XTDXC3zDPDWqro/yQuA+5Lc1a17b1W9q3fjJNuAncAlwEuAf07ys1V1YpKF68fLrK7SAa/U0fo08Mi9qo5V1f3d8lPAIWDTGabsAG6rqqer6lHgMHDZJIqVJA1nVefck8wBlwL3dENvSvJgkluTnN+NbQIe75m2SJ9vBkl2JVlIsrC0tLT6yiVJKxo63JM8H/g48JaqehJ4P/AyYDtwDHj3yU37TK/TBqr2VdV8Vc1v2ND3vwCUJI1oqHBP8iyWg/1DVXU7QFU9UVUnquoHwAf40amXRWBLz/TNwNHJlSxJGmSYq2UC3AIcqqr39Ixv7NnstcBD3fIdwM4k5yW5GNgK3Du5kiVJgwxztczlwOuBLyV5oBt7O3Bdku0sn3I5ArwBoKoOJjkAPMzylTY3eqWMJJ1dA8O9qj5P//Ponz7DnD3AnjHqkiSNwU+oSlKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoIHhnmRLks8mOZTkYJI3d+MvTnJXkq92z+f3zLkpyeEkjyS5apoNSJJON8yR+zPAW6vqF4BXAjcm2QbsBu6uqq3A3d1runU7gUuAq4H3JTlnGsVLkvobGO5Vdayq7u+WnwIOAZuAHcD+brP9wGu65R3AbVX1dFU9ChwGLptw3ZKkM1jVOfckc8ClwD3ARVV1DJa/AQAXdpttAh7vmbbYjZ36XruSLCRZWFpaGqF0SdJKhg73JM8HPg68paqePNOmfcbqtIGqfVU1X1XzGzZsGLYMSdIQhgr3JM9iOdg/VFW3d8NPJNnYrd8IHO/GF4EtPdM3A0cnU64kaRjDXC0T4BbgUFW9p2fVHcD13fL1wCd7xncmOS/JxcBW4N7JlSxJGuTcIba5HHg98KUkD3Rjbwf2AgeS3AA8BrwOoKoOJjkAPMzylTY3VtWJSRcuSVrZwHCvqs/T/zw6wJUrzNkD7BmjLknSGPyEqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNOnfWBWj15nbfOesSJK1xHrlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBnkppDTArC49PbL32pnsV20YeOSe5NYkx5M81DP2ziTfSPJA97imZ91NSQ4neSTJVdMqXJK0smFOy3wQuLrP+Huranv3+DRAkm3ATuCSbs77kpwzqWIlScMZGO5V9Tng20O+3w7gtqp6uqoeBQ4Dl41RnyRpBOP8QvVNSR7sTtuc341tAh7v2WaxGztNkl1JFpIsLC0tjVGGJOlUo4b7+4GXAduBY8C7u/H02bb6vUFV7auq+aqa37Bhw4hlSJL6GSncq+qJqjpRVT8APsCPTr0sAlt6Nt0MHB2vREnSao0U7kk29rx8LXDySpo7gJ1JzktyMbAVuHe8EiVJqzXwOvckHwGuAC5Isgi8A7giyXaWT7kcAd4AUFUHkxwAHgaeAW6sqhNTqVyStKKB4V5V1/UZvuUM2+8B9oxTlCRpPN5+QJIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoYLgnuTXJ8SQP9Yy9OMldSb7aPZ/fs+6mJIeTPJLkqmkVLkla2TBH7h8Erj5lbDdwd1VtBe7uXpNkG7ATuKSb874k50ysWknSUAaGe1V9Dvj2KcM7gP3d8n7gNT3jt1XV01X1KHAYuGwypUqShjXqOfeLquoYQPd8YTe+CXi8Z7vFbuw0SXYlWUiysLS0NGIZkqR+Jv0L1fQZq34bVtW+qpqvqvkNGzZMuAxJ+vE2arg/kWQjQPd8vBtfBLb0bLcZODp6eZKkUYwa7ncA13fL1wOf7BnfmeS8JBcDW4F7xytRkrRa5w7aIMlHgCuAC5IsAu8A9gIHktwAPAa8DqCqDiY5ADwMPAPcWFUnplS7JGkFA8O9qq5bYdWVK2y/B9gzTlGSpPH4CVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaNPB+7pJmY273nTPb95G9185s35oMj9wlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CBvHDaGWd7YSZLOxCN3SWqQ4S5JDRrrtEySI8BTwAngmaqaT/Ji4KPAHHAE+N2q+s54ZUqSVmMSR+6vqqrtVTXfvd4N3F1VW4G7u9eSpLNoGqdldgD7u+X9wGumsA9J0hmMG+4FfCbJfUl2dWMXVdUxgO75wn4Tk+xKspBkYWlpacwyJEm9xr0U8vKqOprkQuCuJF8edmJV7QP2AczPz9eYdUiSeox15F5VR7vn48AngMuAJ5JsBOiej49bpCRpdUYO9yTPS/KCk8vAbwEPAXcA13ebXQ98ctwiJUmrM85pmYuATyQ5+T4frqp/TPJF4ECSG4DHgNeNX6YkaTVGDveq+jrw8j7j3wKuHKcoSdJ4/ISqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaNPJ/kC2pXXO775zJfo/svXYm+22RR+6S1CDDXZIa1MRpmVn9CClJa5VH7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWhq17knuRr4K+Ac4Oaq2jutfUlqg7c9mJypHLknOQf4a+DVwDbguiTbprEvSdLppnXkfhlwuKq+DpDkNmAH8PCU9idJI5vlp9yn9VPDtMJ9E/B4z+tF4Jd7N0iyC9jVvfyfJI+s8F4XAN+ceIVriz2uf633B/Y4FfmLsab/zEorphXu6TNW/+9F1T5g38A3Shaqan5Sha1F9rj+td4f2ON6M62rZRaBLT2vNwNHp7QvSdIpphXuXwS2Jrk4ybOBncAdU9qXJOkUUzktU1XPJHkT8E8sXwp5a1UdHPHtBp66aYA9rn+t9wf2uK6kqgZvJUlaV/yEqiQ1yHCXpAbNNNyTXJ3kkSSHk+zus/7nk/xrkqeTvK1nfEuSzyY5lORgkjef3cqHM2p/PevPSfLvST51dipevXF6TPKiJB9L8uXua/krZ6/y4Y3Z4x91f0cfSvKRJD959iof3hA9/l6SB7vHF5K8fNi5a8Go/a2XrOmrqmbyYPkXrV8DXgo8G/gPYNsp21wI/BKwB3hbz/hG4BXd8guAr5w6d9aPcfrrWf/HwIeBT826n2n0COwH/rBbfjbwoln3NOG/p5uAR4HndK8PAH8w655G7PFXgfO75VcD9ww7d9aPMftb81mz0mOWR+4/vEVBVX0POHmLgh+qquNV9UXg+6eMH6uq+7vlp4BDLP9DWktG7g8gyWbgWuDms1HsiEbuMckLgV8Hbum2+15VffesVL06Y30dWb4i7TlJzgWey9r8vMcwPX6hqr7Tvfw3lj+7MtTcNWDk/tZJ1vQ1y3Dvd4uCVf+hJZkDLgXumUxZEzNuf38J/AnwgwnWNGnj9PhSYAn4m+7U081JnjfpAidg5B6r6hvAu4DHgGPAf1fVZyZe4fhW2+MNwD+MOHcWxunvh9Zw1vQ1y3AfeIuCgW+QPB/4OPCWqnpyIlVNzsj9Jflt4HhV3TfZkiZunK/hucArgPdX1aXA/wJr8XztOF/H81k+QrwYeAnwvCS/P8HaJmXoHpO8iuXw+9PVzp2hcfo7Ob6Ws6avWYb7WLcoSPIslv+wP1RVt0+4tkkYp7/Lgd9JcoTlHyF/I8nfTba8iRinx0VgsapOHgV9jOWwX2vG6fE3gUeraqmqvg/czvK53bVmqB6T/CLLpwl3VNW3VjN3xsbpbz1kTV+zDPeRb1GQJCyfqz1UVe+ZYo3jGLm/qrqpqjZX1Vw371+qai0e8Y3T438Bjyf5uW7oStbmLaHHuZXGY8Arkzy3+zt7JcvnbNeagT0m+WmWvzm9vqq+spq5a8DI/a2TrOlvlr/NBa5h+bfPXwP+rBt7I/DGbvmnWP6u+yTw3W75hcCvsfxj1YPAA93jmln/dnpS/Z3yHlewRq+WGbdHYDuw0H0d/57uaoW19hizxz8Hvgw8BPwtcN6s+xmxx5uB7/T8e1s409y19hi1v/WSNf0e3n5AkhrkJ1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQ/wHsMe03F0DenAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(distribution_of_means(100,1000))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.13662244337810703\n",
      "0.19933707397340616\n"
     ]
    }
   ],
   "source": [
    "# the formula now works.\n",
    "print(mean - 2*std/np.sqrt(100-1))\n",
    "print(mean + 2*std/np.sqrt(100-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determining independent points: the autocorrelation time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important: this is only valid for a _stationary_ timeseries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas implements an autocorrelation function as a function of _lag_ ($\\tau$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.35636797425369826"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start = 500\n",
    "panda_pot = pd.Series(potential[start:]-np.mean(potential[start:]))\n",
    "panda_pot.autocorr(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Programming Exercises\n",
    "\n",
    " * Use the `autocorr` function to calculate the entire autocorrelation function for all $\\tau$.  \n",
    " * Bonus points: Compute the autocorrelation function from scratch and compare.\n",
    " * Can you use your estimate the lag time $\\tau$ at which the system becomes uncorrelated (i.e. the autocorrelation function) to show points spaced more than $\\tau$ are independent? \n",
    " * Note: the formulas assume that the average if the timeseries is zero, so you should subtract it off first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = 700\n",
    "stationary = pd.Series(potential[start:])-np.mean(potential[start:])\n",
    "n_autocorr = np.shape(stationary-1)[0]\n",
    "nlim = int(n_autocorr/2)\n",
    "acf = np.zeros(nlim)\n",
    "for i in range(nlim):\n",
    "    acf[i] = stationary.autocorr(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABB4klEQVR4nO3dd3hcZ5X48e+Zpt675CK519hOnMROw4QUJ4EUYBdCSYDNhixkNywsEGCX/mMXWDphQ2BDDQmBhCRAKinOprrFcS9yU7Gs3vvMvL8/7r2jmdFIGtmSRrbP53n8WHPnau5rWXPPnLecV4wxKKWUUiNxJboBSimlpjcNFEoppUalgUIppdSoNFAopZQalQYKpZRSo/IkugGTIT8/35SXlye6GUopdcrYsmVLkzGmINZzp2WgKC8vZ/PmzYluhlJKnTJE5OhIz2nXk1JKqVFpoFBKKTUqDRRKKaVGpYFCKaXUqDRQKKWUGlVCA4WI3CsiDSKyc4TnRUR+KCKVIrJdRM6e6jYqpdSZLtEZxS+B9aM8fxUw3/5zK/A/U9AmpZRSYRIaKIwxLwIto5xyHfBrY3kNyBaRkslqzw+fPcCG/Y2T9fJKKXVKSnRGMZYyoDrscY19bBgRuVVENovI5sbGE7vZ373hIC8d0EChlFLhpnugkBjHYu60ZIy5xxiz2hizuqAg5ir0MXlcwmBAN3JSSqlw0z1Q1AAzwx7PAI5N1sV8HheDgeBkvbxSSp2SpnugeAy4yZ79tAZoN8bUTdbFPC4NFEopFS2hRQFF5H5gHZAvIjXAlwAvgDHmbuBx4GqgEugBPjyZ7fF6BL92PSmlVISEBgpjzI1jPG+Aj09Rc/C6XAxoRqGUUhGme9fTlPK6XZpRKKVUFA0UYTxu0TEKpZSKooEijNftYjCoGYVSSoXTQBHG6xYG/ZpRKKVUOA0UYbxuF/6gBgqllAqngSKMx+1iQAezlVIqggaKMD7telJKqWE0UITxuLTrSSmlommgCOP1uLQooFJKRdFAEcbr0nUUSikVTQNFGK9biwIqpVQ0DRRhtCigUkoNp4EijEeLAiql1DAaKML4PFoUUCmlommgCOPRwWyllBpGA0UYq4SHwdoGQymlFGigiOB1C4CupVBKqTAaKMJ43daPQ1dnK6XUEA0UYTx2oBj0a0ahlFIODRRhfE7Xk2YUSikVooEiTCij0JlPSikVooEiTGiMQgezlVIqRANFGGfWk67OVkqpIQkNFCKyXkT2iUiliNwZ4/ksEfmziLwpIrtE5MOT2R7NKJRSariEBQoRcQN3AVcBS4AbRWRJ1GkfB3YbY1YA64DviIhvstrkcTnrKDSjUEopRyIzivOASmPMIWPMAPAAcF3UOQbIEBEB0oEWwD9ZDfJ6dDBbKaWiJTJQlAHVYY9r7GPhfgwsBo4BO4A7jDGTdhf3upxAoV1PSinlSGSgkBjHou/QVwLbgFJgJfBjEcmM+WIit4rIZhHZ3NjYeEINcgaz/ZpRKKVUSCIDRQ0wM+zxDKzMIdyHgYeNpRI4DCyK9WLGmHuMMauNMasLCgpOqEHOOgqd9aSUUkMSGSg2AfNFpMIeoH4v8FjUOVXA2wBEpAhYCByarAb5dNaTUkoN40nUhY0xfhG5HXgKcAP3GmN2icht9vN3A18DfikiO7C6qj5rjGmarDZ53DrrSSmloiUsUAAYYx4HHo86dnfY18eAK6aqPc46isGgZhRKKeXQldlhQvtR+DWjUEophwaKMLofhVJKDaeBIownVOtJu56UUsqhgSLM0KwnzSiUUsqhgSKM7kehlFLDaaAIExrM1q4npZQK0UARZqjWk2YUSinl0EARxuUS3C7RldlKKRVGA0UUr1s0o1BKqTAaKKJ4XS4do1BKqTAaKKJ4PS7NKJRSKowGiigel+jKbKWUCqOBIorX7WLAr11PSinl0EARxevWjEIppcJpoIjidesYhVJKhdNAEcXj1llPSikVTgNFFJ+uo1BKqQgaKKJ43C5dma2UUmE0UETxuoUBzSiUUipEA0UUr9ul+1EopVQYDRRRvDqYrZRSETRQRPG4dDBbKaXCaaCIorWelFIqkgaKKF6XaNeTUkqFSWigEJH1IrJPRCpF5M4RzlknIttEZJeIbJjsNulgtlJKRfIk6sIi4gbuAi4HaoBNIvKYMWZ32DnZwE+A9caYKhEpnOx2edwuBjSjUEqpkERmFOcBlcaYQ8aYAeAB4Lqoc94HPGyMqQIwxjRMdqN8WhRQKaUiJDJQlAHVYY9r7GPhFgA5IvKCiGwRkZtGejERuVVENovI5sbGxhNulMftYtCvgUIppRyJDBQS41h0n48HOAe4BrgS+A8RWRDrxYwx9xhjVhtjVhcUFJxwo7xuF4NB7XpSSilHwsYosDKImWGPZwDHYpzTZIzpBrpF5EVgBbB/shrl1aKASikVIZEZxSZgvohUiIgPeC/wWNQ5jwIXi4hHRFKB84E9k9kor9uFMRDQrEIppYAEZhTGGL+I3A48BbiBe40xu0TkNvv5u40xe0TkSWA7EAR+bozZOZnt8ritHrHBQBC3yz2Zl1JKqVNCIrueMMY8DjwedezuqMffBr49VW3yua0kazAQJNmrgUIppXRldhSPy8kotOtJKaVAA8UwXo/1I9HV2UopZdFAEcVrdz3p5kVKKWXRQBHFaw9m63aoSill0UARxRs2mK2UUkoDxTAelxMoNKNQSinQQDGMzzO0jkIppZQGimGcjEIryCqllEUDRZTQrCe/dj0ppRRooBgmNOtJMwqllAI0UAyjs56UUiqSBoooQ0UBtetJKaVAA8UwPs0olFIqggaKKB63U+tJMwqllAINFMM4g9la60kppSwaKKJ4NaNQSqkIGiii6KwnpZSKpIEiSvhWqEoppTRQDDM060m7npRSCjRQDONshao73CmllEUDRRS3SxDRriellHJooIgiInhdLgaD2vWklFKggSImr1sY9GtGoZRSMEagEJG/s/+umIyLi8h6EdknIpUicuco550rIgERefdktCNaktdNnz8wFZdSSqlpb6yM4nP23w9N9IVFxA3cBVwFLAFuFJElI5z3TeCpiW7DSPLSfDR1DkzV5ZRSalrzjPF8s4g8D1SIyGPRTxpjrj2Ja58HVBpjDgGIyAPAdcDuqPP+GStQnXsS1xqXwswkGjr7pupySik1rY0VKK4BzgZ+A3xngq9dBlSHPa4Bzg8/QUTKgBuASxkjUIjIrcCtALNmzTqphhVmJLPpSMtJvYZSSp0uRg0UxpgB4DURucAY0zjB15ZYl4x6/H3gs8aYgEis08O+0Zh7gHsAVq9efVJTlgoykmjs7McYw1jXVUqp0128s55+JyLZzgMRyRGRkx0zqAFmhj2eARyLOmc18ICIHAHeDfxERK4/yeuOqTAjiX5/kI4+/2RfSimlpr2xup4cBcaYNueBMaZVRApP8tqbgPn2jKpa4L3A+8JPMMaEZluJyC+BvxhjHjnJ646pICMJgMbOPrJSvJN9OaWUmtbizSgCIhLq+BeRcoZ3E42LMcYP3I41m2kP8KAxZpeI3CYit53Ma58sJ1A0dPYnshlKKTUtxJtRfAF4SUQ22I8vwR44PhnGmMeBx6OO3T3CuR862evFqzAjGYBGDRRKKRVfoDDGPCkiq7GCwzbgUaB3EtuVUIWZdkbRoYFCKaXiChQicgtwB9aA8zZgDfAq1rTV005Gkockj4vGLg0USikV7xjFHVjrGI4aY94KrAImerrstCEi1qK7Dl10p5RS8QaKPmNMH4CIJBlj9gILJ69ZiVeYkayD2UopRfyBosZeR/EI8IyIPMrwNQ+nlUJ70V24m+/dyLef2pugFimlVGLEO5h9g/3ll+3aT1nAk5PWqmmgICOJVw42hx43d/WzYX8j3f26CE8pdWaJd3psiDFmw9hnnfoKM5Jo7x2kbzBAstfNq4esoFHd2hM65+GtNSR73Vy9vCRRzVRKqUk37kBxpnDWUjR19TMjJ5WXK61AUd/RHwoeP3j2AP2DQdYvLcbl0ppQSqnTk+5wN4Lo1dmvHGzC67aCQW1bL/3+ANUtPRzv6OON6rZENVMppSadBooRhAJFRz81rT0cbe7hiiXFAFS39FDV3IOzrfYTO+oS1UyllJp0GihG4KzObuzq5xW72+k951rFbqtbeznU1A1AaVYyT+w8jjEnVfpKKaWmLQ0UI8hLS8Il8MbRVv66o478dB8XzsvH53FR3dLDoUYrUPzjJXOobetlR217gluslFKTQwPFCNwuYUZOKg+/UcuG/Y1csqDAOpadYgeKLgoykrhhVRkel/D4juOJbrJSSk0KnfU0ij/ctpYaezrsouJMAGbmplLd2kOSx82c/DSyU32cPStHt05VSp22NFCMoigzmaLM5IhjM3NT2Fbdhktg/TJr/URxVjJv1rQloIVKKTX5NFCM08ycVNp7BwGYW5AGQH56Ek1aF0opdZrSMYpxmpmbGvp6jh0oCjKS6B4I0DOg5T2UUqcfDRTjNDMnLFDkpwOQn+4DoKlzICFtUkqpyaSBYpxm5qYA4HULM3Ksr/MzhtZcKJUou46188kHt+EPBBPdFHWa0UAxTlkpXjKSPMzKTcXjtn58Bel2oNBxCpVAf9lex8Nbazlor/FRaqJooBgnEWFxaSYrZmSHjjnlPppGyCgONnbx2qHmmM8pNVEqG7oA2F/fGfP56pYePvyLjaHJGErFS2c9nYB7P3QubhmqFpubZo9RjBAovvaX3ew73smrn3vblLRPnZkO2oHiwAiB4rm9DTy/r5E3qlpZt7BwKpumTnEaKE5AelLkj83rdpGb5ovZ9RQIGrYcaaVrwM9gIIjXrUmcmnj9/gBHW6zFofvru2Kes88OIEebe2I+r9RI9K41QfLTfTEzin3HO+ns92MMHG/vS0DL1JngSFMPgaDB6xb2N8TOKJxM43CTjmGo8UlooBCR9SKyT0QqReTOGM+/X0S2239eEZEViWhnPPLTk2jqGj49dvPRodIedRoo1CRxxicunl/A0eYe+v2BiOeNMew7bgWKI80aKNT4JCxQiIgbuAu4ClgC3CgiS6JOOwy8xRhzFvA14J6pbWX8rEAxPKPYdKQVj737XV1771Q3S50hKhu6EIH1S4sJBM2wrKGhs5+OPj8i2vWkxi+RGcV5QKUx5pAxZgB4ALgu/ARjzCvGmFb74WvAjCluY9wKMpKGjVEYY9h0uIVLFhQAcKwtMqP42+56/vPxPVPWRnX6qmzsoiw7heUzsoDh4xRONrFqZjbVLT261kKNSyIDRRlQHfa4xj42kn8AnhjpSRG5VUQ2i8jmxsbGCWpi/PLTk+iJKuNR29bL8Y4+1i0sICvFG5FR7DrWzu33b+WnLx4asfSHMYbOPp3KqCIZY4ZtlFXZ0MW8wnTmFKThdsmwmU/OlNkrlhbjDxpqWjW7VfFLZKCQGMdibhMnIm/FChSfHenFjDH3GGNWG2NWFxQUTFAT4xerjMfmI1YytHp2LiVZyRxrs96cbT0D3PbbLQz4rU91Iw0u/teTezn/G8+G5r0Hgoadte30DQZinq/ODPdvrGbV157h0W21gPV7caixi3kF6SR53MzOSx22luJAfRd5aT7OmZ0D6DiFGp9EBooaYGbY4xnAseiTROQs4OfAdcaYabtqrSBUxmOoe+mlyiYykjwsLM6gNDsl1PX0w2crOd7exzduWA4Q2i0v3NO7jvPTDYfoGQiw65i1e96f3qjl7T96iRVfeZqP/HITvQMaMM5Em4600NYzyB0PbOOTv9/G3uMd9PuDzC+yao8tKMzgQHTXU30nC4oyKM+zClke0ZlPahwSGSg2AfNFpEJEfMB7gcfCTxCRWcDDwAeNMfsT0Ma45YfKeFgZxY6adh7eWsO1K0txu4SSrORQ19PGI82cW57L9avKELFWboerbunhU394k/mF1ht/V20HAJsOt5CZ7OGa5SU8t7eBbdVtU/SvU8Gg4aO/2cyL+6e+WzNaVUsP55bncMfb5vPItlr+/u5XAZhn/74sKErnSHN3KPM0xnCgvpMFRenkp/tI87k5ogPaahwSFiiMMX7gduApYA/woDFml4jcJiK32ad9EcgDfiIi20Rkc4KaO6bwMh6DgSCffWg7+elJfGb9IgBKs1No7RmkrWeAvXWdrJyZTbLXzYyclGEZxQObqugdCPC/N59LSVYyO+2M4s2aNlbNyuHOq63X3Hu8Ywr/hWe22rZentpVz1O7Er/l7dHmHubkp/Ovly/gwY+uJTvVh9ctzCvIAKA8P42gIdTVWdvWS/dAgAXFGYgI5flp2vWkxiWhK7ONMY8Dj0cduzvs61uAW6a6XSfCKePR2NnPPS8eYnddB3d/4ByyUrwAlGZbO+U9vbsef9CwcmY2YJUqj84o3qxuZ2FxBrPyUllamsWuYx109/vZX9/JlUuLKUhPIi/Nx966oX5oYwwisYZ91EQ4kTUIL+xrIDfNx1lhdcFOVs+An6aufmblWeXuV5fn8uQnLqauvY+sVOt3rdjelbG+o585BemhbqgFRUOBZFdt+4S1SZ3+dGX2BHHKeDy16zjfeXofbz+rhPXLikPPl2RZJcmf2FEHwMpZ2YC1+dHhpu7QLBZjDNtr2kI3l2VlmRxs7GLj4RaCBlbOzEZEWFSSEcootte0seprz/CnN2qm6F975nHKXxxpir/L5vMP7+C/n57YHtPqFitLmBW2gVZGsjcUBAAKQ4HCGhOrtvd9n20Hl/K8VGpaexmMc4psV7+fTz64bcQaUur0p4FiAuWn+9h73Bo0/Oa7zop4rtQOFC9VNlGWnUJhhvVmnluQTs9AgOP2m/pocw8dfX7OsufDLy3Nwhi4f2MVACvsTGRxcSb76jsJBA1/3VFHW88gn3zwTR7cVI2aeE5Gcay9N65ZZ32DAeo6+jg6wV08zuuFB4poRZlWN6gTKI619eF1C/lp1vHyvDT8QUOtPUXWGMMrlU0EgzEnHfLtJ/fy8NZant3bMGH/DnVq0UAxgWbkpJKb5uNnN60mLapwYFFWEiIwGBjqdoKh7VQPNlg3gO12l4ATKJaVZQLw7N4GZuWmhrq4FpVk0jcY5EhzNy8daGLlzGwumpfPZx7azvaatsn8Z56R9td34nYJxsS3srmmtQdjiPnJ/ZcvH+bn/3fohNpR1RKZHcSSkewlzeemvsNaAFrX3ktxVjIuu0KAs52vs5bijeo23vfz17nv9aPDXmvL0RZ+/Zp13BnzUGceDRQT6L/euZy//PNFEftqO5I87tDMqPBAMbfAmqlyqMnqR95e3UaSxxXqSijOTCYvzUcgGBlgFhVbz79ysJldxzp426JC7nr/2fjcLh7bNmyWsToJg4EgBxu7WDsnD4ivqJ4TTAIxFrf9fnMNP/jbgdA6mgc3VXPnQ9tHfN09dR088oa1ZqKqpYeMZE9o7GskRZnJoYyirq0v1PUJUGhPvGjotJ53Mou7NxxiMBBk17F2Vn31aS777gY+dt9WSjKTKc9L1UBxBtNAMYEKM5MpzU4Z8fnSLKu7yRmfAOtNm57kCe0lsL22nSWlmaFy5CLC0jIru1gRFijmFabjdgm/eOkwABfOzycz2ctF8/N5YufxYSt31Yk73NTNYMBw5dIiIL4B7fCsI3rNQn1HH539fl471MxgIMg3n9zLA5uquey7G/ju0/uGvdaXHt3FJx/cRmv3AFUtPczOSx1z4kJ4oDjW3hv63YOhMYwGu+SMc15tWy8PbKzijge24Xa5mJ2bSorXzbfevYK5Bem6mvsMpoFiCpVkpeB2CctKs0LHRIQ5BWkcauomEDTsqm3nrLKsiO9bWmp1P62cOXQ82etmTr71fRnJntD3rF9WTG1bLzsSPKslGDTTsp7QgD/IL18+PK7Fis74xDmzc8lL88W1WK2qpSdUDDI8U+gbDNDSba21eWZ3PS/sa6S5e4Bvvms5Vy4t4kfPV0a8fmVDFxuPWBMZnt3bQFVzz6jjE46izCTqO/sIBg31HX2UhH2ASU/ykOZz02B3TTV09uPzuFhSkskXH9tFZUMX33vPCv73Q+fywqffykXz8+0FoxoozlQaKKbQu86ZwcfXzSXF5444Prcgnb3HO9lytJXugcCw6ZTXrijlupWlLC+LPL64xAogF8zNC+3ffcWSIjwu4YmdiZ3v//HfbeVj922d8utWNfdwy6820d0fu37WX3cc48t/3s1DW0efIdbvD7D++y/y0JYa9h23xifmFqZRnp8WZ9dTN/OLMshI8kQMaDs3Z5/bxTO76/njlmry0ny88+wZfPkdS/G6XPzi5cOh8x/YWIXHJeTZM+pqWnuZlZs25vWtjKKfxq5+BgMmIqMAK6uot7ueGjr6KMpM4vZL52EM3HJRBRfPjyyDU5aTQkefX2uPTQPVLT2hDy9TRQPFFLp8SRGfvGLhsOOXLiqkqauf99xjrbB1BrIdi0sy+cF7V+HzRP53LSqxxikumpcfOpad6mPt3Dye2FGX0O6nffWd/G1PPQ0dU7sHx4b9DfxtT0NoOmu0R96wxm+eH2MGT21rL3uPd/KFR3bw7N4GKvLTSPK4Kc8bWqzW3e8nMMJMoaMtPczOTWV2fiqHw7qhnNlt15xVwvGOPp7aVc91K8vwul0UZiZz7cpSHtxcQ1vPAP3+AA9treGKpUVctbyY5/Y2MBAIxplRJDPgD7K7zppCXZwV2SVamJFEY4fT9dRPUUYyVy0r5k8fu4A7r1o07PWcLtXoCshqalU193D9XS/z7rtfobV7+P43k0UDxTTwjhWlPPWJS3jHWaWcX5HLHHuAeywXzcunMCOJty0uiji+flkxR5p72FmbuJXbTZ39BA38eXvdlF7XGRtojrGJVFNXPy9VNpHsdfHywaZRp7k6N/QBf5A9dR0stCcXVOSnUt/RT1VzD5d863l+9NyBYd8bCBpqWnqZnZdKeV5aREbhvO77zp+F2+6aetc5Q0WT/+GiCnoHA3z3mf1875kDtPYMcuN5s7hiSXEoKI0248lRZI9DbKtqA6AkRkbhDGbXd/ZRmJmEiLBqVk4oOw1XFgoUE9P91NDRN6U3utNBe+8gH/7lRgYCQbr7/fz4+copu7YGimliQVEGP7xxFb//6NrQDWQsZ83IZuMXLhs2gP725aVkpXj5ryf3JCSrGPAH6eizun6cCqdTxZk+2tI9fBOpv26vIxA0fOryhfQNBnn14Mg1Jp0B3s9fvRiAxSVOoLCC+Md+t4Xm7gFeifEaxzv6GAgEmZ2XRnleWsQU2eN2va8FRRlcPD+fs2ZksTRszGpxSSYXz8/n168e5e4NB1kxM5sL5+azZk4eGfaU63jHKIBQPbDo35HCjCTqO/oxxtDY0R9a1zMSJ1DUjDNQHGrsoq1neEC45debufPh7eN6rTPdvz+yk6qWHn5202refc4Mfv3qEapbpqZmlwaK01BWqpd/vWw+L1c288zu+gl5za5+f9yD085g7ZyCNLbXtHMoqkTJZHICRaxtaR/dVsui4gw+uHY2KV43z43S/XS83Qo07z1vFg9/7AI+dGEFAOX51k16Z20HWSledtS0D/u5OBnE7LxUyvPTIqbIHm/vJ8XrJjPZw13vO5vf3nL+sGt//z0r+d0t5/Pip9/KQ7etxeUSfB4X6xYV4nXLsOwgllBGYU+3zkn1Rj2fRO9ggMbOfjr7/aHzR1KYkYTXLePOKN73s9f55pN7I44ZY9hf38mWo206Oy9Ode29PL6jjo9cWMGaOXl88vKFuF3Ct58aPktuMmigOE29f81s5hWm8/W/7qGyoXPEvvR4BIKGy7+7gR8+F1+q62wJ+6ELyhGBR09gXcd9rx/l7K89w13PV444MB3NGBOWUUQGiuqWHrZWtXHdyjKSvW4ump/Pc3sb2HSkhc//aQfP7K6PuGnVd/SRnuQhPcnD2bNySLc/zTtluufkp/H5qxfROxigMioQVtndX7NyUym3u4mcmUz1HX2UZCUjIqQlechMHr4eIi89iQvm5TMrLzWiG+jOqxZxzwdXx+wailZoZxTtvYOUZqcMm07rZBDO7DhnbcVIXC6hOGxPlXj0DPg53tHHm9WRM/Aau/rpGwzS1NUfWhSoRnf/xmqCxvD+82cDUJyVzE1ry/nrjrpQ9juZNFCcprxuF196xxJqWnu47LsvsvKrT7PzBKfM7qxtp669j9cPDe9mue03W/jRs5H99E6gWFqaydo5eTy6rXbcnxy3VbXR1jPAt5/axw0/eTmu72/qGqDHnvbaHLV/ubOnhzPwf+miQmrbevm7u1/l95uq+cdfb+Y997wWCjD19kygaGlJHr5y7VLuev/ZnFeRF2pruKMtPXjdQml2CuX59v4PdpZxvKNvzE/vIynLTuGtiwrjOjfJ4w5lEbEyECcwOIEinjaVZaeEFufFwxn4rmzoCi0uhKF6VeHXVyMbDAR5YGMV6xYUhIpBArzvvFkEgoY/bJ78sj0aKE5jF88v4NlPreOb71pOZ5+f/zvQdEKv81Kl9X27j3VE1APq7Bvk6d3HeSFqjwZnIDkvLYnrV5ZxpLmHN2vGd0No6OxnWVkWd7xtPvvru+iMI6uoahkaNG6Oyiicrh+nf3/90mIumpfP565axBtfvJyvXLuUjYdb+JO9Avp4Rx/FI3Tx3HxBOYtLMinPSyUrxcubUSVTqpp7mJGTitue1pqe5AllFMfbR37diebc/Euyhi8CdRbd7QwFitEzCmDcaymccwfsle2O8H51DRRje2Z3PQ2d/XxgzeyI4+X5aaydk8fvN1ePWKdromigOM1V5KfxnnNnUZadEpoqOV4v2QGms98fqkQKVjn0oGHYGESzPZCcl+5j/fJifB7XuAe16zv6KMxICpVDiWeGjNPtVJGfNmzWU7VT+sL+lJ2T5uO3t5zPR98yl8xkLzetnU1Wijf0b6lvH/uTv4iwYmY228K6Vlq7B9h1rD0UkESEuQVp7DneGVr8dqIZxXg513FK3IdzuqaGup7iyyiOd/TFPVYVHlT2hP3uOf9Ps3JTTzjLPRN09g3ykxcq+Y9HdlKWncK6hcOzyfeeN5Pqll5ePnhiHwLjpYHiDLG4JDPizRqv3oEAW462smZOLkDElNstR609wVt7BiNu5M1dAyR5XKTbffCXLizkz2/WjWuldmNnP4WZyeSmWTf26AwhlqPNPYhY61CGjVG09jIjZ+TZQiLCvMJ0Khu6CAYNDZ39cd3QV87IYn99Jz0Dfp7YUcfl39tATWsv7zx7aMrr+XPy2FbVRm1bL/6giWsweiI4WUKsjCIjyUOy10V9Rz9JHheZKWNvTVOWnULQDE3xHcuxtl5cAj6PK+J3r7qlh6LMJFbPzpnwjOJoc3eo0vKp7tZfb+FbT+5jSWkm99x0TszZkFcuLSY71csDGye3+0kDxRliSWkmhxq74iqRHW7TkRYGAkE+cmEFHpeE+voBtlS1hr4+FLZaubGrn/z0pNAA6vWrSmnq6o85lTSWAX+Q5u4BCjOSyLVLY8ebURTb9baau/sjxjVqWnuYmTNyHS6AuQVpHGzspqm7H3/QhDYAGs2KmdkEgobPP7yDf7pvK6XZKTx2+0Vct3IoUFw0L5+BQJDH3rQG9acqo3DaXxIjoxCRUDuKMpPj2vRqpEV3g4Eg19/1Mj+OWlNS22ZlTwuLMtgTtslWVUsPM3NSWVaWRWNn/4QOxt73ehWfe3gHPQPxTYCYrqqae3j1UDOfunwBv/mH8yOmUIdL9rp556oZHGzsmtSSORoozhBLSjIIGsZc+l/Z0BlRB+nlyia8buGi+fnML8pg5zHrk2EwaHijqpUL5loDuuHdT81dA+Sl+0KP1y0sJCPZwyNxdj85g+GFGcnkplqvE50hxFLV3MPM3FTy0nwMBkxoXMMYQ3VLb8yqvuHmFqTT1NUf2hEunhu6U27lkW3HuGZ5CX+4bS1L7NpcjnPLc/G5XTy0xSobMmVjFPZ1SmNkFDA0oB3P+ARYZTwAatsi5+5vOtzCtuo2/vvp/RFdjMfaeinNTmFxSQZ76jpCgbu6xapXtdyuQPBGVRvfe2Z/zDLn41XXbgWdWAsuTyWPvWn9HN95zowxz/3M+oU8ccfFcc2GO1EaKM4QS0qsN+Vo4xRvVLVyxfde5BevDNUaeqmyibNn5ZDq87CsNJNdte0YY6hs7KKzz8/1K8vwuiUio2ju7g+VVAfrU8/F8/PZerSVeDhVTYsyk8hNH0egsMtmOEHKuVk0dw/QOxhgxpgZhbWY7mV78D6eG3pBRhJXLi3iwxeW88MbV5HkcQ87J8Xn5pzZOaGfUTyZykR4+/JSvnLtUhYUxV7p74xLxDM+AVbXk8clVDZEjkk9vbueJI+Lc8tz+PQft7PDnrhwrN0JFJk0dw/Q2NnPgD9IXUcfM3NTWVKSiQh8+g9v8oNnD/D//rrnpGtJ1dnjIk1dp+60W2MMj247xrnlOaGFjqNJ9ronfRtkDRRniBk5KWQkeUJ9xdFdUP3+AJ/543aCBnbZ4xA9A35213Vwvr0Pw9JS6w1f39EfGp84ryKXWbmpERlFU+cAeWm+iNefk59OdZzbbzpdEYUZyaT53PjcLlpirO4N1zsQoKGzn9l5qaHuKmd1tjPLZuYoYxQAcwutG6rTRRbvDf2nH1zNl96xdNQV9RfNt6blusTaCXEqZKV6ufmC8hFvIs6AdmGcGUWy121llWHjVMYYntldz8Xz8/npB1eT6nPzi5cPEwwa6tr6KM1ODhWv3F3XQW1bL8ZYmyelJXlYWJRBvz/IRy+ZQ89AILTvxolyMopYCy5PFXuPd3KgoYtrw7ovE00DxRnC5bL22d59rMPaY/urz/Db14ZS/bueP8iBhi6KM5PZbxfU21/fhTFDZc6X2aXMdx1rZ8vRVnLTfMzOS2VOQTqHGq1Py8YYmrv7yUuPvPnMzkslELb95mjCMwoRISfNO+YYhTOTxul6gqGbhTM1dqyup5k5KXjdwvaatgm/oTvrNwozkie1i2A8nExiPGMmy0oz2WlnlTB0879iSTG5aT7WVOSx8UgLTd39DASClGWnsLjY+v3ZU9cZMeMJ4MfvW8Xjd1zE565ezPKyLH77WtUJr9YO2LPKYPg6mlPJo9uO4XEJ1ywvSXRTQqbHb6yaEktKMtl7vJPP/HE7vYMBvvP0Ptp7B9lW3cZPnq/khlVl3HB2GYebukPF8IDQG32x3VVw58M7eHRbLWfPygntp3G0uYdA0NDR52cwYIbdZCvshWeH7YVn++s7eWWEKX2NHX24hFCwyU1LGrPraahsRlqo68n5HmdK71hdTx63i/K8NILG6lKayBv6srIsMpM9oXGD6cAZoxhrVXa4ZWVZNHcPhGY+Pb2rHhG4dLE1dfPcilxqWntD3YylWSlkpXopz0vlyV3HqYra83teYQbzCq06Wu8/fxb76jtD2ep4NXdZkxDg1O56emFfA2vm5IW2PZ4ONFCcQRaXZNLV72fv8U7+9bIFtPUO8q0n93L777ZSlJnMl96xhIVFGfiDhiPN3eyt6yDN5w7dYNOSPNx43iyrXtKacj6z3iqZPic/jYFAkNrW3tAnufxhGYUVKI7a/fTfeHwPt/1mS8zSIvUdVkbidOXkpnnHDBSbjrTgc7uYX5geeoM5balu6SU3zTdsH/NYnHGKiR5HcLuET16+gPefN2tCX/dkOKt8x8q0wjlZpTMO8czuelbPzgn9f59Xbk2jdsq5OwPgt186nzer2/j5S4fxeVwxg9O1K0vJSPZwb9h+HONxrH1o9tSp2vUUDBoON3WHilBOF2O/cyaRiKwHfgC4gZ8bY/4r6nmxn78a6AE+ZIyZ+t1wThPObJxrzirhjsvmc7S5m/tetzbGefC2tWSn+phvD3zur+9kz/FOFpVk4grre//GDcuHva5TFv1gUxdpPutXKi8qo8hPt1coN/dgjGF7TTsdfX521rZHbPEK1l7O4TNxctOS2NHaNuq/7fl9jZw/JzcUDDKSPKG1F/FMjXXMK0yHXZMzhdUpLDhdrJ6dwx9vW8s5s3Pi/p4lJZm4BHYe66AiP43ddR38+zWLQ88vLskgPckTKrjoTKl956oyfvvaUbZVtzGnIC3id8qR6vPwDxdV8P2/HeCRN2q5ftX4+uidyrwip25GUdfRR78/GCr9Ml0kLKMQETdwF3AVsAS4UUSWRJ12FTDf/nMr8D9T2sjTzPKyLL75ruV843rrZv+pKxdSlJnEF65ZzNmzrJvF3IJ0XGJNo91T18Gi4rE/2cyxf6kPNXaHPsXnpUV+YhQRZuelcqS5m2PtfaEMwSkPEq4+qux1buroGUV1Sw+VDV0RK1dz032hWU81rb3MiPNT89xC698yVWsdEklEWF2eO64ZMyk+N/MK09lZ2879G6vxuoUbwm7oHreLs2fnMBAI2gsurcDtcglfuXYpIqNPKrj9rfM4tzyHL/xpx7irDjvrO2KtzD9VOKVeKqZZoJBElfkVkbXAl40xV9qPPwdgjPnPsHN+CrxgjLnffrwPWGeMGXU3nNWrV5vNmzefULvWrVt3Qt93OtlW3YbHJXT1+6nIT4vrprn5SCt56T5SfW4ON3VzzuwcvFF9/Acauuju9zMrN5X99Z24RMhI9rC4JBNjIGAMHpew5WgrOak+5hRYb5aa1l5qWns4vyKPWPe0+o4+Djd1s3JmNslea3rqrmPtuERYXJLJxsMtFGclx7WPQ3e/nx217czMTY1rauKZ6GBDF229gxggK8XL/MLI6be1bb1Ut/SQ6nMP29a3obOfZI+LzJThVXMdA/4g22vbSU/yxPVBxVHV0sPx9j6yU730DgZYEXXtU4Hzu3z2rJxhO1rG44UXXjjha4vIFmPM6ljPJXKMogwIX3deYx8b7zkAiMitIrJZRDY3NjbGOkXFKdXnpsterJbqi693MiPZQ0v3AH2D1vRXj2v4r1ayx0W/P0hXvx8RoSAjic4+P0F7Xcb2mjaCxjAYCEa8SbxuKzr4g7Gn1rb1DJLsdYeChHP9wYBhwB8kaAxJcb7pUnxuslN9ZI1yIzvTpSV5GAwE8QeCFMUYa8iwswhfjDUlhRlJowYJ6/us/TN6x7m6ut9v/d543S78gVNzn4u+wSAukRMKEpMpkWMUsfLd6P/deM6xDhpzD3APWBnFiTbqZCLy6eI7T+/jR/beEy98+QoyYuyZEG17TRvX/vhlkrwuFnvdbPjiFcPO+cPmaj79x+3MyE1lUZKHT1+5kA//chNXvWUOP91wiFzgthuWU/WnHXz5+mWhapl/2X6M23/3Bvd+4hIWFmcw4B8KJH2DAVZ+9Wk+eu4svnzt0tC17nxoO8/tbeA//24FN9+7kbs/ch5vWVAwAT8dtfFwC3//01cpz0vl+X9bN6zrqm8wwFlfeZq/Xz2Dr18/fEwrHt96ci/3vHiI575+VczxjFje/T+v4HEL51Xk8aPnDvC3r181baYix+sjv9zEsbZenvzEJYluSoREBooaYGbY4xlA9A438ZyjJth8e3/ombkpcQUJsEpZXLO8hL/uqBuxD9oZoKtq6eE9q2dyXkUuHpfw0w2HqMhPo61ngJ/93yEgcozAmcXU0j3Aw1tr+OSDb1Jur99wsph1CyODQG6aj5buAb715F4KM5LGNWCrRre0NJPMZA8fvrAi5vhGstfNbz5yXsTeCeNVlJmMP2ho6RkYNoNuJHXtfZxXkUtBug9joKVnIO5V59PFkaZuFk2zGU+Q2K6nTcB8EakQER/wXuCxqHMeA24SyxqgfazxCXXynJIPi4ozxzgz0ievWGDtwTDCQjVndziA5TOySEvysGpWNgBfuXYp71hRymF7MC98+mR4oHhxfyNZKV4WFWdyvL0vNJi61q455chLT8IfNOw61sGX3rE0tEOdOnlpSR5e//xl3LR29ojnnD8nL2bV2ng5s97iLRjoLLYryUoOrb9p6pzcAe0D9Se3c2Q0fyBIVUtPxPtkukjYu8cY4xeR24GnsKbH3muM2SUit9nP3w08jjU1thJreuyHE9XeM8mc/HQykz3j/hQ+tyCdr1y7dFj5DoczRbar389yez7+7ZfOZ09dB5csKCAj2cOvX7VWi4eXlQgVBuwZ4M2ads6vyOXuD54zalucNqxbWMDVy4vH9e9QY0vxDR9/mEjOxkoNHf0sLR37/CZ7sV1JVnIoA3H2RZkMx9p6ufL7L/Kdv1/BDavGLtwXj5pWqwz9dJsaCwleR2GMeRwrGIQfuzvsawN8fKrbdabzeVw8/2/rxhx0jCV6F65wzhTZ/fWdofT6LQsKQmMHK2dmU5GfxpHm7ojuhhz7pn+4sZvDTd28O46KmmfPymHtnDy+dt2ySS+Ypiae0/UYb0bh1HgqyUoJZbSTuZZif30nQUOo0nAsL1c2UZSZbK3NiYNTtWCOBgp1qoiu1TRRVs/OITfNF7PKqojw8bfO4/m9DRFTa71uFxnJHl7YZy3iWhm1QC+WWXmp3H/rmglrt5paBelO11Pkzd4fCOIPmogZbjBUNbY4PKOYxLUUThdpVdi2ruEG/EFuvncjQWN4z7kz+ez6RWSnjl6S47BdL00zCnXG+8p1y0Yt+vbuc2bEzBhy03wcaupGhNA+Bur05fO4yEvzUd8ZmVF8+c+7eGLHcX7/0TWhGlEwlFGUZqeQmezB53bRGEdGYYyhtWdw3HWVnIVx1SMEitaeAfxBw/KyLP6wuYaOPj93ve/s0V+zuZuMZM+IXbeJdGrNHVOnhRPpCnLeyHML0smMcyaWOrUVZibTENb11DcY4E9ba2nuHuADP99ITdj+7YeaukjzuclJ9SJiTaiIJ6N4qbKJ1V9/ho2HW8bVtkNjZBTOtT/+1rnceskcHt9Rx4H60TcNO9zUTUV+2rTsKtVAoU4JzoD2qbjaVp2YosykiK6nZ/c00D0Q4D/evoSeAT+3/GpzKDt9/VAL54SVI8lL98U1RrGjtp2ggf98Ys+4ypsfsccTWnsGY2625Ayk56YlccvFc0jxukNrk2Kpa+9l4+GW0CSP6UYDhTolOBnFypnT842kJl5RRnLEYPaj22opzEjiQxeU85n1i9h7vJP99V3W9rUNXayZkxs6Nz89Ka6MwhkXeKOqjad318fVrn5/gNrW3lDpkuqW4XusOLXJctN85Kb5uGltOX/efixid8Cb7t3IZ/+4nWDQ8J2n92OAf1o3N642TDUNFOqUMBQodOHcmaIoM4mmrn4CQUN7zyAv7GvkHStKcbuEy5cUAfC3PfWhbqM1c4bW0uSlJcXMKPyBILuOtYceO3XJ5hak8a0n9+KPYwfG6pYeggYusWfrxep+coKUM97wjxdX4HW5eHCzVZHIGMNrB5v5/eZq/um+LTy0tYYPX1DOjDF2YUwUDRTqlLBqVg5LSjJZOI4icerUVpiZTNBY+4o8uauOgUCQ61ZaiyqKMpNZXpbFs3vqee1QM6k+d0S3TX6GNUYR3Z301b/s5pofvkRVs3VzP9LczfzCdG6/dB4HG7vZXttOLP5AkBf2NeAPBEO7OTrTumMNaLd0D+B2SahmWF56ErPyUkPXbekeYCAQpDwvlad21ZOV4uVjb513Mj+uSaWzntQpYf2yYtYv04VzZ5KhtRT9/GV7HeV5qRHB4LLFRXz/2f3Utfexujw3Ykr1vIJ0BgJBXjvUElq1/8K+htCCzjeqW8lK9dLUNUBFfhrz7RlUTZ3Ds5DW7gH++f43eKmyiX+/ZjFBO/ismJFNVoo3dkbRPUBOqjeiTlVJVjJ19p4Zziytz65fxJHmHhYWp0/rQpSaUSilpiWnjMf++k5ePdjM+mUlETOC3ra4EGOsm274+ATAO1aUkp/u4+4NBwHrE/yn/7idBUXpJHtdbK9pj9j7wVmk1xy170l7zyDX3fUyGw+3UJqVzH2vV3GosZvcNB9ZqV5m5qbEDBQt3f3DptyWZqVQa++ZcTxsOu8/rZvLpYuKTvjnNBU0UCilpiUno7h/YxX+oOHKpZE306WlmZTYe5CHj0+AVZjwQxeUs2F/I1uOtnLLrzbR3jPI99+ziqWlWWyvaQstmqvITxu2fa7j6d3HqWrp4ec3r+bfrlzI4aZu/rrDym7A2vt7pK6nYYEiO4Wmrn76/QHqOpyV5KdG0UINFEqpaSkvzYdLYPPRVoozk4dNjRYR1i8rJifVG3Na6QfXlJPmc/O+n73Gtuo2fnjjSpaUZrK8LIudtR1UNnQhYq3iT/K4yUj2DNtre8P+Rgoykrh4fj5XLy8hJ9VLZ5+finxrxtPM3FRqWnsJRhUHbO4eGLbLY0m23ZXW3s/x9l48Lpm0CggTTQOFUmpa8rhdoXIcVywtirkvxWfXL+KJOy4ZtpsiQFaql/evmU2/P8jXr1/O+mUlAKyYmUXvYIC/7alnRk5KqJxMfnpSRNdTIGj4vwNNvGVBASJCstfN36+2dj2oyB/KKAYCwWEryGNmFHY13WPtvdS19VGUmYw7zr02Ek0DhVJq2nK6n9YvjT2RIdnrpniU7ptPX7mQx//lYt53/qzQseVl2QDsPd4ZygzAymDCu56217TR3jsYmgYLVtHLrBQv58y2xkSc7XW//dQ+Lvrmczy/15oZ1RajLIiTURxr66WuvW/Udk83GiiUUtNWcVYy2alezqvIHfvkGLxuF0tKI/dVmZOfRoa9P0l4pdbosh8b9jciAhfPyw8dm5mbyrYvXh6aSeUEioe31lLb1suG/Y209gyGXi+ck1HUtfdxvEMDhVJKTYjPXLmQn9+0ekK3NHW5hGX2mEZ52C58eemRi/Q27G9kxYzsUJl7R/jMq1m5qXz9+mU89E9rOWtGNgcaOiNWZYdLsWtR1bb1UtfeS0mmBgqllDpp84syWF1+YtnEaM6yKxBXFAx1PeWn+WjpGSAQNLT1DPBmdVtEt1MsIsIH1szmnNm5LChMZ399V1idp+FVYEuyUthT10HfYJCS7BPfAXCqaaBQSp1x3rKwgPQkD0tKhrql8tKTMMYqEb6tuo2ggbVR025HM78oncbO/tC02+hZT2BNkd1V2wGcOlNjQVdmK6XOQBfMzWfHl6+I6EYKLbrrGgjd7OcXxbc7nXWutbr79UNW7alYGUVpdjIDdj0pHaNQSqlpLnrfBycDaO6ysoLxbiLkVJN97VAzADmpw0tylGSlhH2tgUIppU4p+c5e291WRjFnnJsIlWalkOpz09DZT3aqN+YAfKk9RdYlQ9u9ngo0UCilFEP7xDd39XOo0dptbjxcLgllFSNtrVpqD2AXZiRP6EyuyXbqtFQppSZRdooXt0uobe3lWHtvxGK8eDn7eI/UZeV0NzmL704VGiiUUgorI8hN87G1qhVjoKJgfBkFDA1+j5RRFGUmI3JqjU9AggKFiOSKyDMicsD+e9i2ZSIyU0SeF5E9IrJLRO5IRFuVUmeOvDQfO+zNi+aMs+sJYEEoUMQef/C6Xbx1YSFr5+bHfH66SlRGcSfwrDFmPvCs/TiaH/iUMWYxsAb4uIgsmcI2KqXOMPnpSQwGrEqw5ScQKOaP0fUEcO+HzuWDa2afWAMTJFGB4jrgV/bXvwKujz7BGFNnjNlqf90J7AHKpqqBSqkzj7OWojAjifSk8S8zK8tO4YNrZnPF0um9EdF4JWrBXZExpg6sgCAihaOdLCLlwCrg9Slom1LqDOWspRjvjCeHyyV87fplE9mkaWHSAoWI/A2IVRv4C+N8nXTgIeATxpiOUc67FbgVYNasWSOdppRSI3IyijknMJB9Opu0QGGMuWyk50SkXkRK7GyiBGgY4TwvVpC4zxjz8BjXuwe4B2D16tVmtHOVUioWZ9HdiWYUp6tEjVE8Btxsf30z8Gj0CWItifxfYI8x5rtT2Dal1BnK6Xoqz9NAES5RgeK/gMtF5ABwuf0YESkVkcftcy4EPghcKiLb7D9XJ6a5SqkzwZq5efzjxRVcNP/Umr462cSY06+XZvXq1Wbz5s2JboZSSp0yRGSLMWZ1rOd0ZbZSSqlRaaBQSik1Kg0USimlRqWBQiml1Kg0UCillBqVBgqllFKj0kChlFJqVBoolFJKjeq0XHAnIo3A0RP89nygaQKbM9Gme/tA2zgRpnv7YPq3cbq3D6ZXG2cbYwpiPXFaBoqTISKbR1qdOB1M9/aBtnEiTPf2wfRv43RvH5wabQTtelJKKTUGDRRKKaVGpYFiuHsS3YAxTPf2gbZxIkz39sH0b+N0bx+cGm3UMQqllFKj04xCKaXUqDRQKKWUGpUGCpuIrBeRfSJSKSJ3Jro9ACIyU0SeF5E9IrJLRO6wj+eKyDMicsD+OyfB7XSLyBsi8pdp2r5sEfmjiOy1f5Zrp1MbReRf7f/fnSJyv4gkJ7p9InKviDSIyM6wYyO2SUQ+Z7939onIlQls47ft/+ftIvInEclOVBtjtS/suX8TESMi+WHHpvxnGC8NFFg3OuAu4CpgCXCjiCxJbKsA8AOfMsYsBtYAH7fbdSfwrDFmPvCs/TiR7gD2hD2ebu37AfCkMWYRsAKrrdOijSJSBvwLsNoYswxwA++dBu37JbA+6ljMNtm/k+8Fltrf8xP7PZWINj4DLDPGnAXsBz6XwDbGah8iMhNrC+iqsGOJ+hnGRQOF5Tyg0hhzyBgzADwAXJfgNmGMqTPGbLW/7sS6wZVhte1X9mm/Aq5PSAMBEZkBXAP8POzwdGpfJnAJ8L8AxpgBY0wb06iNgAdIEREPkAocI8HtM8a8CLREHR6pTdcBDxhj+o0xh4FKrPfUlLfRGPO0McZvP3wNmJGoNo7wMwT4HvAZIHwmUUJ+hvHSQGEpA6rDHtfYx6YNESkHVgGvA0XGmDqwgglQmMCmfR/rlz4Ydmw6tW8O0Aj8wu4e+7mIpE2XNhpjaoH/xvp0WQe0G2Oeni7tizJSm6br++cjwBP219OijSJyLVBrjHkz6qlp0b6RaKCwSIxj02besIikAw8BnzDGdCS6PQ4ReTvQYIzZkui2jMIDnA38jzFmFdBN4rvCQux+/uuACqAUSBORDyS2VeM27d4/IvIFrK7b+5xDMU6b0jaKSCrwBeCLsZ6OcWza3IM0UFhqgJlhj2dgpf8JJyJerCBxnzHmYftwvYiU2M+XAA0Jat6FwLUicgSru+5SEfntNGofWP+3NcaY1+3Hf8QKHNOljZcBh40xjcaYQeBh4IJp1L5wI7VpWr1/RORm4O3A+83QQrHp0Ma5WB8I3rTfMzOArSJSPE3aNyINFJZNwHwRqRARH9ag0mMJbhMiIlh963uMMd8Ne+ox4Gb765uBR6e6bQDGmM8ZY2YYY8qxfmbPGWM+MF3aB2CMOQ5Ui8hC+9DbgN1MnzZWAWtEJNX+/34b1ljUdGlfuJHa9BjwXhFJEpEKYD6wMQHtQ0TWA58FrjXG9IQ9lfA2GmN2GGMKjTHl9numBjjb/h1NePtGZYzRP9aHjquxZkkcBL6Q6PbYbboIK/3cDmyz/1wN5GHNOjlg/507Ddq6DviL/fW0ah+wEths/xwfAXKmUxuBrwB7gZ3Ab4CkRLcPuB9rzGQQ64b2D6O1CatL5SCwD7gqgW2sxOrrd94vdyeqjbHaF/X8ESA/kT/DeP9oCQ+llFKj0q4npZRSo9JAoZRSalQaKJRSSo1KA4VSSqlRaaBQSik1Kg0USgEi0jWBr/UJEblplOfXicgFE3W9sNf9W6Ir9arTkwYKpSaQXdjvI8DvRjltHdbq64n2G+Bjk/C66gyngUKpMCKSLiLPishWEdkhIteFPfcf9l4Hz9j7RvxbjJe4FNhq7AqmIvIvIrLb3h/hAbu4423Av4rINhG5WEQKROQhEdlk/7nQ/t4vi8hvROQ5ew+If7SPl4jIi/b37xSRi+1rPwbcOIk/HnWG8iS6AUpNM33ADcaYDntTmddE5DHgHOBdWBV8PcBWIFYxxAujjt8JVBhj+kUk2xjTJiJ3A13GmP8GEJHfAd8zxrwkIrOAp4DF9vefhbUXSRrwhoj8FSsYPGWM+X/2ngWpAMaYVrsERJ4xpnkCfybqDKeBQqlIAnxDRC7BKp1eBhRhlVN51BjTCyAifx7h+0uI3MRpO3CfiDyCVT4klsuAJVapJwAyRSTD/tq5Zq+IPI+1R8Em4F67YOQjxphtYa/VgFWFVgOFmjDa9aRUpPcDBcA5xpiVQD2QTOwy0LH02uc7rsHaPfEcYIs9hhHNBaw1xqy0/5QZa6MqGF5q2hhrQ5xLgFrgN1ED58l2G5SaMBoolIqUhbXHxqCIvBWYbR9/CXiHWPtZp2MFgFj2APMARMQFzDTGPI+1uVM2kA50Ahlh3/M0cLvzQERWhj13nX3NPKxB8E0iMttu48+wqgufbX+fAMVYxeaUmjDa9aRUpPuAP4vIZqzqo3sBjDGb7LGKN4GjWNVo22N8/xNYs4/A2v/6tyKShZWRfM8eo/gz8Ed7oPyfsfbMvktEtmO9J1/EGvAGq9T0X4FZwNeMMcfs/RY+LSKDQBfgZBTnAK+Zoa1AlZoQWj1WqTiJSLoxpsveqexF4FZj72kedd6fgM8YYw6c5PW+TNigdxzn/wB4zBjz7MlcV6lo2vWkVPzuEZFtWDOeHooVJGx3Yg1qT7WdGiTUZNCMQiml1Kg0o1BKKTUqDRRKKaVGpYFCKaXUqDRQKKWUGpUGCqWUUqP6/yW2i1u4UqksAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.arange(nlim)\n",
    "plt.plot(np.arange(nlim),acf)\n",
    "plt.hlines(0,0,nlim,'k')\n",
    "plt.xlabel(\"lag (steps)\")\n",
    "plt.ylabel(\"acf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When does the ACF become essentially zero?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's your guess looking at the ACF?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting ACF to an exponential to estimate correlation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exponential(x, r): #function f(x, r) = e^(r*x)\n",
    "    #return np.e ** (r * x)\n",
    "    return np.exp(r*x) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This function calculates tau by using an exponential fit.\n",
    "#To do this it uses scipy.optimize curve_fit\n",
    "#curve_fit will optimize any function to fit given data\n",
    "\n",
    "def tau_calc(ac_data, function = exponential): #takes in data and a python function\n",
    "    x_data = np.arange(len(ac_data))\n",
    "    pars, cov = curve_fit(f=function, xdata=x_data, ydata=ac_data, p0=[0], bounds=(-np.inf, np.inf))\n",
    "    #curve fit returns an np.array of optimally fit paramters(pars) and their coverience(cov)\n",
    "    #pars in this case will return an optimized value of k to fit the dataset\n",
    "    return -1/pars #tau = -1/k\n",
    "# scipy curve_fit documentation: https://towardsdatascience.com/basic-curve-fitting-of-scientific-data-with-python-9592244a2509"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.7911753567330577\n"
     ]
    }
   ],
   "source": [
    "tau = tau_calc(acf)[0]\n",
    "print(tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This $\\tau$ is the time it takes to go from 1 to $1/e$. So we actually can show we want to go out $2\\tau$ (there's some theory that this far enough!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABH+UlEQVR4nO3dd3icxbX48e9sVe/Fsorl3uUmF2xwodrgYJNAAnEK3BAgQEJILuCQX4Dk5hJCuCkklBBCCaGE0AyEXoyNC7Zs5CoXyZZsVatYve/O7493d62VVtWS5ZXO53n0yLvvu7vzytLZ2TNnZpTWGiGEEP7PNNgNEEII0T8koAshxBAhAV0IIYYICehCCDFESEAXQoghwjJYLxwTE6NTU1MH6+WFEMIv7dixo0xrHevr2KAF9NTUVDIyMgbr5YUQwi8ppfI6OyYpFyGEGCIkoAshxBAhAV0IIYaIQcuhCyG8tbS0kJ+fT2Nj42A3RZwFAgICSEpKwmq19vgxEtCFOEvk5+cTGhpKamoqSqnBbo4YRFprysvLyc/PZ/To0T1+XLcpF6XUU0qpE0qpvZ0cV0qph5VS2Uqp3Uqp2b1otxDCpbGxkejoaAnmAqUU0dHRvf601pMc+jPA8i6OrwDGu75uAB7rVQuEEB4SzIVbX34Xug3oWusNQEUXp6wC/qENW4EIpVRCr1vSQ0f3b2frkz+morRooF5CCCH8Un9UuSQCx9vcznfd14FS6galVIZSKqO0tLRPL1ZxLIsF+U9TXdJpbb0Q4gx444032L9//2A3A4Dc3FymTZvW7TkvvPCC53ZGRgY/+tGPBrppZ1R/BHRfnwt87pqhtX5Ca52utU6PjfU5c7V7gZEAOOvK+vZ4IUS/GOiA7nA4urzdW+0Denp6Og8//PBpPefZpj8Cej6Q3OZ2ElDYD8/rW1AUALrh5IC9hBDD1erVq5kzZw5Tp07liSeeACAkJMRz/JVXXuHaa69l8+bNvPnmm9xxxx3MnDmTnJwcMjMzWbBgAWlpaVxxxRWcPGn8jWZnZ3PhhRcyY8YMZs+eTU5ODlpr7rjjDqZNm8b06dP517/+BcD69etZtmwZ3/zmN5k+fXqH2w6HgzvuuIO5c+eSlpbGX//61w7XkJuby3nnncfs2bOZPXs2mzdvBmDt2rVs3LiRmTNn8oc//IH169ezcuVKACoqKli9ejVpaWksWLCA3bt3A3DffffxX//1XyxdupQxY8ac9W8A/VG2+CZwq1LqJWA+UKW1HrAEt3L10KmXgC6Grl++tY/9hdX9+pxTRoZx71emdnnOU089RVRUFA0NDcydO5evfe1rPs9buHAhl19+OStXruTKK68EIC0tjT//+c8sWbKEe+65h1/+8pf88Y9/ZM2aNaxdu5YrrriCxsZGnE4nr732GpmZmezatYuysjLmzp3L4sWLAdi2bRt79+5l9OjRrF+/3uv2E088QXh4ONu3b6epqYlFixZx8cUXew0gxsXF8eGHHxIQEMDhw4e55ppryMjI4IEHHuChhx7i7bffBow3D7d7772XWbNm8cYbb/DJJ5/wne98h8zMTAAOHDjAp59+Sk1NDRMnTuQHP/hBr2rDz6RuA7pS6kVgKRCjlMoH7gWsAFrrx4F3gEuBbKAeuG6gGgtgDjZ66Ep66EL0u4cffpjXX38dgOPHj3P48OEePa6qqorKykqWLFkCwHe/+12uuuoqampqKCgo4IorrgCMyTIAn3/+Oddccw1ms5n4+HiWLFnC9u3bCQsLY968eV61121vf/DBB+zevZtXXnnF87qHDx9mwoQJnvNbWlq49dZbyczMxGw2c+jQoW7b//nnn/Pqq68CcP7551NeXk5VVRUAl112GXa7HbvdTlxcHCUlJSQlJfXo53KmdRvQtdbXdHNcA7f0W4u6YbUHUa/tqMbKM/WSQpxx3fWkB8L69ev56KOP2LJlC0FBQSxdupTGxkav3m9v66I724S+q83pg4ODO72ttebPf/4zl1xyidc5ubm5nn//4Q9/ID4+nl27duF0Oj1vIr1tp/u67Xa75z6z2Uxra2u3zzdY/G4tF7vVxElCMDd1VUkphOitqqoqIiMjCQoK4sCBA2zduhWA+Ph4srKycDqdnt47QGhoKDU1NQCEh4cTGRnJxo0bAXjuuedYsmQJYWFhJCUl8cYbbwDQ1NREfX09ixcv5l//+hcOh4PS0lI2bNjAvHnzum3jJZdcwmOPPUZLSwsAhw4doq6ursN1JCQkYDKZeO655zyDqW3b297ixYt5/vnnAeONLSYmhrCwsJ7+6M4afhfQbWYTVToEc1PVYDdFiCFl+fLltLa2kpaWxi9+8QsWLFgAwAMPPMDKlSs5//zzSUg4NcXk6quv5ne/+x2zZs0iJyeHZ599ljvuuIO0tDQyMzO55557ACO4P/zww6SlpbFw4UKKi4u54oorSEtLY8aMGZx//vk8+OCDjBgxots2Xn/99UyZMoXZs2czbdo0brzxxg495ptvvplnn32WBQsWcOjQIU8PPy0tDYvFwowZM/jDH/7g9Zj77ruPjIwM0tLSWLt2Lc8+++xp/SwHi+rqo89ASk9P133Z4OJEdSPZv1vGuBg7cbet7/+GCTFIsrKymDx58mA3Q5xFfP1OKKV2aK3TfZ3vdz10u8XMSUKwNksPXQgh2vK7gG6zGCkXe3PlYDdFCCHOKn4Z0E8Sgr21CgYpXSSEEGcjvwvoZpOimlDM2gHNtYPdHCGEOGv4XUAHqDWFGv+ol9JFIYRw88uAXmd21YfKbFEhhPDwy4Be7wno0kMXQgg3vwzoDRbpoQvRnyorK3n00UcHuxleKzv2p4FcC/2ZZ56hsPDUArPXX3/9oK0T758B3Rxu/ENy6EL0i7MloA+UgVwLvX1Af/LJJ5kyZUq/PHdv9cfyuWdcs9XdQ68c1HYIMWDeXQvFe/r3OUdMhxUP+Dy0du1acnJymDlzJhdddBH33nsvq1at4uTJk7S0tPDrX/+aVatWkZuby8qVK9m719gz/qGHHqK2tpb77rvP5/NmZ2dz0003UVpaitls5t///jfx8fE+n7u9Bx98kOeeew6TycSKFSt44AHvtpeWlnLTTTdx7NgxAP74xz+yaNEiPvvsM2677TbAWGBrw4YNrF27lqysLGbOnMl3v/tdZs2a5VlK97777uPo0aMUFRVx6NAhfv/737N161beffddEhMTeeutt7BarfzqV7/irbfeoqGhgYULF/LXv/6VV199lYyMDNasWUNgYCBbtmxhxYoVPPTQQ6Snp/Piiy9y//33o7Xmsssu47e//S1gfBK57bbbePvttwkMDGTdunXEx8f36b+1Lb/soZusdhpUoOTQhegnDzzwAGPHjiUzM5Pf/e53BAQE8Prrr7Nz504+/fRTfvrTn3a5QmJn1qxZwy233MKuXbvYvHkzCQkJPXrud999lzfeeIMvvviCXbt2ceedd3Z47ttuu43bb7+d7du38+qrr3L99dcDxpvMI488QmZmJhs3biQwMJAHHniA8847j8zMTG6//fYOz5WTk8N//vMf1q1bx7e+9S2WLVvGnj17CAwM5D//+Q8At956K9u3b2fv3r00NDTw9ttvc+WVV5Kens7zzz9PZmYmgYGBnucsLCzkrrvu4pNPPiEzM5Pt27d7Fimrq6tjwYIF7Nq1i8WLF/O3v/2t1z9bX/yyh26zmKgxhRIoOXQxVHXSkz5TtNbcfffdbNiwAZPJREFBASUlJb16js7WQm9pafH53G0X5/roo4+47rrrCAoKAiAqKqrD83/00Udeuerq6mpqampYtGgRP/nJT1izZg1f/epXe7R2+YoVK7BarZ5dkZYvXw7A9OnTPUvzfvrppzz44IPU19dTUVHB1KlT+cpXvtLpc27fvp2lS5fi3m5zzZo1bNiwgdWrV2Oz2Ty7Jc2ZM4cPP/yw2zb2hF8GdLvFRI0KJU4CuhAD4vnnn6e0tJQdO3ZgtVpJTU2lsbERi8WC0+n0nNfV+uid9eg7e+72j227DrsvTqeTLVu2ePWKwUgfXXbZZbzzzjssWLCAjz76qLvL9ax5bjKZsFqtntc2mUy0trbS2NjIzTffTEZGBsnJydx3333drg3f1Seatq/Rn2us+2XKxW4xUU2oDIoK0U/arxVeVVVFXFwcVquVTz/9lLy8PMBYG/3EiROUl5fT1NTk2c7Nl87WQu/sudu6+OKLeeqpp6ivrweMPT99nfOXv/zFc9u9ZVxOTg7Tp0/nrrvuIj09nQMHDnS5FnpPuIN3TEwMtbW1nh2ToPN11ufPn89nn31GWVkZDoeDF1980bOj00Dxy4Bus5ioIkTKFoXoJ9HR0SxatIhp06Zxxx13sGbNGjIyMjz54UmTJgFGz/Kee+5h/vz5rFy50nN/Z3ythd7Zc7e1fPlyLr/8ctLT05k5cyYPPfRQh3MefvhhzxrmU6ZM4fHHHweMwdFp06YxY8YMAgMDWbFiRZdrofdEREQE3//+95k+fTqrV69m7ty5nmPXXnstN910EzNnzqShocFzf0JCAr/5zW9YtmyZZ4NsX4O//cnv1kMHuP1fmSw+/BuusG2HO4/0c8uEGByyHrpob8ivhw7GrkWVTlcPvU0+TwghhjP/HBS1mijXwaCd0FQNgRGD3SQhhrVbbrmFTZs2ed132223cd111w1Si4YnvwzoNrOJCkew8fmi4aQEdDFk9KS642z0yCOPDHYThpy+pMP9M+ViMVHqMDZ+lYFRMVQEBARQXl7epz9kMbRorSkvL/fU7veUX/bQ7RYz5U53QJfSRTE0JCUlkZ+fT2lp6WA3RZwFAgICejQpqi2/DOhG2aI7oFcOaluE6C9Wq5XRo0cPdjOEH/PblMtJLbsWCSFEW34Z0O1ePXTJoQshBPhpQLdZTDgw47SFSQ5dCCFcehTQlVLLlVIHlVLZSqm1Po6HK6XeUkrtUkrtU0oNaPGp3WI02xEQIT10IYRw6TagK6XMwCPACmAKcI1Sqv12HLcA+7XWM4ClwP8ppWz93FYPd0BvtUdIDl0IIVx60kOfB2RrrY9orZuBl4D2K8xoIFQZMyJCgAqgf9aD9MHmDui2COmhCyGES08CeiJwvM3tfNd9bf0FmAwUAnuA27TWHRZZUUrdoJTKUEplnE6trc1sBqBZAroQQnj0JKD7mofcfirbJUAmMBKYCfxFKRXW4UFaP6G1Ttdap7t38egLu9VodpM1XAZFhRDCpScBPR9IbnM7CaMn3tZ1wGvakA0cBbpeKPk02MzugB5mTCySFReFEKJHAX07MF4pNdo10Hk18Ga7c44BFwAopeKBicCALVTuzqE3mMMBDY2VA/VSQgjhN7qd+q+1blVK3Qq8D5iBp7TW+5RSN7mOPw78D/CMUmoPRormLq112UA12l3l0mBxZXUaTkJQx01khRBiOOnRWi5a63eAd9rd93ibfxcCF/dv0zrn7qHXmdsEdCGEGOb8dqYoQI0l0rijTlanE0IIvwzodotRtlhtdqVZaksGsTVCCHF28NOAbjS7yhRu3FF7YhBbI4QQZwe/DOjussVGpxkCo6SHLoQQ+GlAN5kUVrOiudUJIfES0IUQAj8N6GD00ptanRASJykXIYTAjwO63WqWHroQQrThtwHdZja5Arqrhy47pQshhjn/DegWE02tDqOH3lIPzbWD3SQhhBhUfhvQ7RYTzQ5XygUkjy6EGPb8NqDbLCaaWlwpF5A8uhBi2PPrgO7dQ5eALoQY3vw2oNst7rJFSbkIIQT4cUC3WcxGQA+MBJNFeuhCiGHPfwO6u2zRZILgOAnoQohhz28Dut1qornVYdyQ2aJCCOHHAd099R9ktqgQQuDHAd1mcaVcQHroQgiBHwd0z8QicPXQT4DTObiNEkKIQeS3Ad0zsQiMgK4d0FAxuI0SQohB5NcB/VQPXWaLCiGE3wZ0u8WMw6lpldmiQggB+HFAt7n2FTWm/7t76DIwKoQYvvw3oLv2FW1ulR66EEKAHwd0u9VoelOrE+whYA2WHroQYljz24Du1UMHVy269NCFEMOX3wZ0u9UM4D1btKZ4EFskhBCDy28DuruH3iTruQghBNDDgK6UWq6UOqiUylZKre3knKVKqUyl1D6l1Gf928yO7Jb2KRdZz0UIMbxZujtBKWUGHgEuAvKB7UqpN7XW+9ucEwE8CizXWh9TSsUNUHs9fAb0xkpobQKLfaBfXgghzjo96aHPA7K11ke01s3AS8Cqdud8E3hNa30MQGs94LkPdx16U2v72aKSdhFCDE89CeiJwPE2t/Nd97U1AYhUSq1XSu1QSn3H1xMppW5QSmUopTJKS0v71mIXm68eOkhAF0IMWz0J6MrHfbrdbQswB7gMuAT4hVJqQocHaf2E1jpda50eGxvb68a2Zbe0r3Jx99Cl0kUIMTx1m0PH6JEnt7mdBBT6OKdMa10H1CmlNgAzgEP90kofTk39d1W5hCcZ36vyB+olhRDirNaTHvp2YLxSarRSygZcDbzZ7px1wHlKKYtSKgiYD2T1b1O9dUi5BMeCJQAqjw3kywohxFmr2x661rpVKXUr8D5gBp7SWu9TSt3kOv641jpLKfUesBtwAk9qrfcOZMPt7QdFlYLwZKg63sWjhBBi6OpJygWt9TvAO+3ue7zd7d8Bv+u/pnWtQw8dICJZeuhCiGFrCMwUbRPQw5OhUnroQojhyW8DeoeUC0BECtSXQXP9ILVKCCEGj98GdKUUNrOpXcolxfgulS5CiGHIbwM6uPYVbZ9yAcmjCyGGJb8O6HaL6dRqi2AMigJUSUAXQgw/fh3QO/TQQxPAZJGBUSHEsOT/Ad3RJqCbzBCWKLXoQohhya8Dut1ioqnF6X1nRIrk0IUQw5JfB/QOPXRwBXTpoQshhh+/Duh2i9l7UBSMSpeaImhtHpxGCSHEIPHrgN6hDh1clS4aqqUWXQgxvPh3QG9f5QJtatEl7SKEGF78OqAbdeg+cugglS5CiGHHrwO6zx56WCKgpIcuhBh2/D6gd+ihW2zGBCMpXRRCDDN+HdCNKhdnxwMRKZJyEUIMO34e0E00ty9bBNnoQggxLPl1QPc5sQiMSpfqAnD6CPZCCDFE+XVAd1e5aK29D0Qkg7PVmGAkhBDDhF8HdJvZhNbQ6mwf0F2li1LpIoQYRvw7oPvahg4g3B3QJY8uhBg+/Dqgu/cV7VCLHjkKlAnKswehVUIIMTj8OqDbLGbAR0C32CEyFcoOnvlGCSHEIPHzgO5OufioZomZCGWHz3CLhBBi8Ph1QO805QIQO8FIuThaz3CrhBBicPh1QO90UBQgZgI4mqEy7wy3SgghBscQDugTje9lh85gi4QQYvD0KKArpZYrpQ4qpbKVUmu7OG+uUsqhlLqy/5rYuVC7BYDqxpaOB2MnGN9LZWBUCDE8dBvQlVJm4BFgBTAFuEYpNaWT834LvN/fjexMUmQQAPknGzoeDAiHkBHSQxdCDBs96aHPA7K11ke01s3AS8AqH+f9EHgVONGP7etSXKgdm8VEfkW97xNixksPXQgxbPQkoCcCbefQ57vu81BKJQJXAI/3X9O6ZzIpkiICOX6yk4Ae6ypdbL/WixBCDEE9CejKx33tI+Qfgbu01l0ub6iUukEplaGUyigtLe1hE7uWFBXE8QofKRcwBkabqqC2pF9eSwghzmY9Cej5QHKb20lAYbtz0oGXlFK5wJXAo0qp1e2fSGv9hNY6XWudHhsb27cWt5Mc2UUPPWa88V3SLkKIYaAnAX07MF4pNVopZQOuBt5se4LWerTWOlVrnQq8AtystX6jvxvrS3JUEJX1LdT4rHSR0kUhxPDRbUDXWrcCt2JUr2QBL2ut9ymlblJK3TTQDexOUmQg0EmlS2gC2EIloAshhgVLT07SWr8DvNPuPp8DoFrra0+/WT2X7CpdPF5Rz+SEMO+DSkmlixBi2PDrmaJgpFwAjvvqoYOr0kV66EKIoc/vA3pkkJVgm5njndaiTzC2omusPrMNE0KIM8zvA7pSiuSoIPK7qkUHWUpXCDHk+X1AB2MJgM5r0V1rushmF0KIIW5IBPTkKKMWXfuaERo5GiwBULz3zDdMCCHOoKER0CODqG92UFHX3PGg2UJNxGRa83ec+YYJIcQZNDQCeheVLrlldbxSHIcuzJTdi4QQQ9oQCejG5CJflS5v7Spkt3MMVmej5NGFEEPa0AjonayLrrVm3a5Cdusxxh0FO72O//njw7y/r/iMtFEIIQbakAjowXYLUcG2Dot0ZRXVkH2ilobQ0dToQK88elOrg4c/Ocwv39xHq8PHFnZCCOFnhkRAB9eqi+1SLm/uKsRiUvzg/AnscY6m5ViG59jB4hpaHJrCqkY+2C/L6woh/N+QCehJkUFeKRetNW/tKuTc8TGcMyaK3Xos9vIsaG0CYHd+FQBRwTae+vzooLRZCCH609AJ6FGBFJxswOk0atF3HjtJQWUDl88YyajoYPYxFpNu9dSj78mvIjLIyi3LxpGRd5Ld+ZWD2HohhDh9QyagJ0cG0exwUlLTCMC6zELsFhMXTx2B1WyiKmq6cWKhMTC6u6CK6UkRfD09iRC7hac35Q5Sy4UQon8MnYDuqkX/1Vv7+fnre1iXWciFk+MJsRsrBEcmjOEkYVCwk8YWB4dLapieGEZogJWr0pN4e3chJ6obB/MShBDitAyZgD5tZBhjYoLZnlvB+/uKCbSa+daCUZ7jE0aE8aVjDI6CHWQVVdPq1ExPjADg2oWptDo1/9yaN0itF0KI09ejDS78QXSInU/+e2mnxyfGh7Jbj2FZ2etk5RlboqYlhQMwKjqYmckRbMutOBNNFUKIATFkeujdmTgilF3OsSg0lTkZxITYSAgP8BxPiQqioLKTFRuFEMIPDJuAnhgRyGHLeACsxV8yPTEcpZTX8aLKRhxOHys2CiGEHxg2Ad1kUkTHJ1FoHsnYukymJ0V4HU+MDKTVqSmRgVEhhJ8aNgEdYFJ8KB83TWG+aT8zRgR6HUuMMG5L2kUI4a+GVUCfMCKUz53TCVZNzDJnex1LinQF9M42mxZC+J2Kumbqm4fPstnDKqBPjA9li3MKDkxEFW/yOjZSeuhCDDlf/+sW7l23b7CbccYMq4A+YUQI1QSTa58IOZ96HQuyGSs2tl+CVwjhn0prmsg+UcvmnPLBbsoZM6wCemyInYVjo2lOWWIsAdBw0ut4YkRgj3voh0tqeG5L7gC0UgjRHzKPVwLGp+6ezgL/x5Zcvv74lgFs1cAaVgFdKcUL31/A5HNXgXbC0Y1exxMjAik42XHXI1/++NFhfrFuH7VNwyc/J4Q/yTx+qsO289jJLs485T+7i9iWW0FpTdNANWtADauA7pE0F2whcMQ77ZIUafTQte66Fr2xxcH6gycAOFpaN2DNFEL0XebxSibEh2Azm9h5rLLb81scTna5Vl3dV1g1sI0bIMMzoJutkHpuhzx6YmQgjS1OKuqau3z4lpxy6podABwpqx2wZgoh+sbp1Ow+XsW80VFMSwxjZ173PfSsomoaW4zdy/YVVg90EwdEjwK6Umq5UuqgUipbKbXWx/E1Sqndrq/NSqkZ/d/UfjZmGZw8CidzPXf1tBb9g/3FBNvMmBTkSA9diLNOTmktNU2tzEiKYHZKJLsLqmhu7XqryR2uoB8WYGH/UA3oSikz8AiwApgCXKOUmtLutKPAEq11GvA/wBP93dB+N3aZ8b1NLz3RVYveVaWLw6n5cH8JyybFkRwVxJFS6aELcbb50jUgOislgtmjImludXabRtmRd5LEiEAWjYsZ0imXeUC21vqI1roZeAlY1fYErfVmrbX7M81WIKl/mzkAYiZAWBIcet9zV1KEsaZ6V5OLMo+fpKy2mYunjmBMTDBHuumhO2VtGCHOuMzjlYQGWBgTE8LslEiAbvPoO/NOMntUJFMSwsgtr6emseUMtLR/9SSgJwLH29zOd93Xme8B7/o6oJS6QSmVoZTKKC0t7XkrB4JSMHU1ZH8E9cayuWGBFkLsli5TLh/sK8FqViydGMuY2BCOlNV2GrT/+lkOM3/1AbllkpYR4kzKPFbJjKQITCbFiPAAEiMCu6x0KaxsoLCqkTkpEUxNDAPgQHHNmWpuv+lJQFc+7vMZwZRSyzAC+l2+jmutn9Bap2ut02NjY3veyoEy/SpwtkDWm4BR1pgYEdhpykVrzfv7ijlnbAxhAVbGxAbT2OKkqF2Nq9aa//vgIL959wDVja1sOdLziQ0FlQ20OLrO9QkhoLy2iVYffysNzQ4OltQwMznCc9+slAi+7GJg1B3s54yKYupIY5+EfQX+l3bpSUDPB5Lb3E4CCtufpJRKA54EVmmt/WNqVsIMiB4He17x3JUY2fnkouwTteSW13PxlHgAxsSEAHjl0Z1OzS/f2s+fP8nmG+nJhAdae7wB9brMAs777Sf8+ePDHY4VVTVwuKRGlvcVAjhZ18yi337Cij9t5KP9JV6lxnsLq3A4tVdAn50SSWFVI8VVvicYZeSeJNBqZlJCKHGhdmJCbH5Z6dKTgL4dGK+UGq2UsgFXA2+2PUEplQK8Bnxba32o/5s5QJQyeum5n0O18R7V1eSi178sAOAiV0AfGxsM4JVHf2t3Ic9szuW6Rak88LXppCWFszu/+3f617/M5/Z/ZeLU8OlB73SU1ppv/HUrF/1hA9PufZ+vPrqJdZkFvb9eIYaIfYVGiWF5XTPX/yODb/x1q6dKJdOVK5+ZEuE5f/Yodx7ddy9957GTzEgOx2o2oZRiysjwoRnQtdatwK3A+0AW8LLWep9S6ial1E2u0+4BooFHlVKZSqmMAWtxf5t2JaBh72uAMbmourG1w4DIJwdKePyzHC5LSyA+zNjpKDbUTojd4tVD/2B/CXGhdu5ZOQWlFNMTwzlYXENji6PTJry6I5+fvLyL+aOjuXHxGPYWVlFVf+r1s0/Ucqyinm+kJ3P1vGSKqxr540cde/FCdKakupFfvbW/29I9f3Gg2Ai27/34PH69ehpHyur42mObuf7ZDD7cX0JSZCAxIXbP+VMSwrBbTJ6g31Z9cyv7CquZ4wr67vMPn6jxu59Xj+rQtdbvaK0naK3Haq3/13Xf41rrx13/vl5rHam1nun6Sh/IRvermHEwchbs+TdwqnSxbdplX2EVt77wJVNHhvO7K9M89yulGBMbzBHXoGerw8nGQ6UsmRDr2Q0pLSmCVqcmq8j3u/2JmkbufHU354yJ5qlr53LB5Hi0hq1HT2Wt1rt67LddOJ57vzKVaxelcrSsjhM1shmH6JlXduTz1KajPZ4Cf7bLKqohPsxOXGgA31owig13LuWOSybyxZFytuVWeKVbAGwWE1NGhrHXR158d76Romkb0KeODKPFoTlU4l8Do8Nzpmh706+CokwoO3xqcpFrYLS4qpHvPZNBeKCVJ7+bTpDNe1/ttqWLXx6vpLqxlWWT4jzH3RtRd5Z2+fxwGQ6n5u5LJxNoMzMzOYJAq5nN2WWecz47VMqE+BDPEr9zU6MA2H50aPxxioHn7pn6Cmj+6EBxNZNGhHluB9ks3LJsHBvvWsZdyyfxg6VjOzxmdHQwxyo6plPdP5tZyd4BHfC7CUYS0AGmfhVQsOffXj30gsoGvvnkVmoaW3jq2rmeVEtbY2JDKKhsoKHZWN/FbFIsGhfjOZ4QHkBMiL3TgL7xcBnRwTamJBi/QDaLibmjozxLftY1tbLtaAVLJpyqCpqWGE6g1cy2o95jzzmltXywr7jbtWjE8OJ0ak/Q6o8AtfqRTT4H7s+UVoeTwyW1TEoI7XAsIsjGD5aO9VSqtDUqOpji6sYO6c89+VWkRgcRGWzz3JcaHUywzcz+Tj5Z91ZtU+sZ+buUgA4QlgBjlsLOfxAToLBZTGw4VMbXHt1MaU0TT183j8kJYT4fOsY1MHq0rI71B0uZkxJJeKDVc1wpxYykcJ+VLk6nZuPhMs4dH4PJdKo6dOHYaA6fqOVETSNbj5TT7HCydOKpXr/VbGL2qAi25Xr30H/26h5ueG4HNz+/0ysHL4a3nNJaqhpasJgUe09zBmRlfTOZxyv57NDgzSM5WlZHs8PJ5BG+/yY7Myo6CK0hv13Rw9GyOsbFhXjdZzIpJieE9cuM0R15FaT/+kNe3Ha8+5NPkwR0t4W3Qk0Rpr2vkBgRyEdZJTi05uUbz2He6KhOH+YuXdx6pJx9hdUsndSxvn56UjjZpbUdlto9UFxDWW0T5433fsyisUYPf0tOOesPlhJkM5OeGul1zrzUaA4UV1PVYATu/JP1bMutYM6oSD7cX8KlD29kR15F738OYsjJcPXOl08bQfaJWhqaOx+g7457sk1WUfWgzYLOcrXBVw+9KynRxkzwvPJTAd3h1Bwtr2NMbEiH86eODGN/YcfrrKxv5vHPcmhq7f7neKy8nu//YweNLU52uZYjGEgS0N3GXgDx02Dzn5k8IoTRMcG89oOFnfbM3UbHGD30Z12bXSydENfhnBlJEWjdcaLCxsNGL+e88TFe908ZGUZYgIVN2WWsP3SChWOjsVvMXufMGx2F1niC9rpMo+zyj9+YySs/WIjZpPjO37dJT12QkXuS6GAbK9NG4tSnKkT6wj24X9fsIM9HPtrtkU+z+enLu/r8Ol05UFSN1aw8nameGhXVMaAXVjbQ3Or0/B23NS4uhLpmB2W13mujv7e3mAfePcAjn2R3eEzbtEpVfQvXPrMNh1MzJiaYnDOw7pMEdDelYOEPoTSLh9PLeP/Hi0l2/QJ0JdBmJjEikLzyeuLD7Ez20WuY3snA6MbDZUyMD+2QmzebFOeMjebdPcUcr2jwyp+7zUqJwGpWfHG0Aq01b3xZQPqoSJKjgpiZHMHj35pDXbODf2Uc681PQQxBO/KMT27TXFPa955GHj2rqBpXAVeX6YgP9pfw6s78DlUinx0q5ZYXdvLe3qIelQQ2tjh4dnOu16eKA8U1jI0NwWbpXfiKCrYRYrd4DYy6K9TG+Ajo7r//9gOp7tuPfZbDwTbLA6w/eILZ//Mh5/zmY9Y8uZWr/7aV4xX1PPHtOcwfE012ae2A59EloLc19asQOhLL1r/06pfFnUdfOiHOU67YVkyIncSIQHa36aE3tjjYllvRoXfutnBsDDWuFM0SH73+AKuZtKQIth2tYH9RNYdP1LJq1qkldqaMDGP+6Cie3Zznc3q0GB5Ka5rILa8nPTWSxIhAIoKs7D+NvHBWUQ1zR0VhMakuJ97klRuB8p9b8zz3tTqc3LNuL//ZXcRN/9zJOb/5mN9/eKjLIPf3z49y75v7eGHbqY5JVlF1t5+cfVFKkRIV5GkbnJrl7SvlktJJQM+rqCcu1E5ogJW7Xt2Nw6n54kg5Nz63g7jQAM4ZE01tk4PqhhYeumoG88dEMzY2mMr6lm73WjhdEtDbsthgwQ8gdyMU7Ozxw9zv7ksndr4+TVq7gdFtRytobnVyno/eNxgDo+7nduf+2pubGsWe/Cpe2nYci0mxcnqC1/HrFo2moLKBj7JKenwtYvBV1DWfVp67LXdKbs6oKJRSTB0Zxt6CvvXQWx1ODpXUkJYUzri4kE4rZirrm6msb8FuMfHazgLP2NG6zELyyut5bM1snr52LjOSI3j448O80cmsZ3euGow3BqdTU1nfTFFVI5NG9C5/7jYqOsgr5XK0rI5Qu4WYEFuHcxMjA1GqY0A/XlHPxBGh/GLlZDKPV3LPur1879kMkiIDeeH78/n9N2ay7pZFbFp7PqtmGp0s96Br9omBTbtIQG9vzrVgD4PP/9Djh8wbHU1sqJ1FnfS2wUi75JXXe3LaGw+XYrOYmJfqe8B1XFwI4+JCWJmW4PM4wPzRUbQ6NS9sO8bSibFeZVdgLFGQFBnIU5tye3wtZ6uy2ib++9+7BryHMxA+O1Ta4wXXGlscXPbwRu58dXefXiu3rI5nN+d6BvIyck9is5g86ZZpI42Zy31ZAC63vI6mVieTE8KY2sXU+FxXwPz+eWOobWrl9S8LcDg1f/k0m8kJYSyfNoJlk+J48jvpzEgK5/53Dvhcqvaxz3KobWrlpiVjOVpWx+accs+g7KQ+9NDBGBg9frLesybSkdI6xsQG+/xkbbeYGREWwPEK77WdjlXUMyo6iNUzE1k8IZbnvzhGZLCV569fQHSb2altjXV9AhjoDXEkoLcXEAYLbjZWYMz9vEcPuSwtgW13X0BYgLXTc9xrMl/9t6089flRPjlwgnmpUQTazD7PV0rxwY8Xc/tFEzp9zjmpkShljNSvntVxRWOzSfHdc1LZdrTC7yeUbM4p55Ud+fx1Q85gN6VX9uRX8d2ntvH+vuIenf/KjnyKqhp5d09RnzYqfv6LPO59cx+/WLcXrTUZeSeZkRTuGVSfMjKMZlcdd2/tLzpVXTJ1ZBhltU2cqO44W9md0lg1cyTTE8N5bksub+0q5GhZHbddMM4TPE0mxa9WTaOstok/tVvKoriqkWc25XLFzER+fOF4ooJt/HNrHgdcg7KT+9hDT40OpsWhKaoygvSR0lqfA6JuyVFBHG/TQ69qaKGyvoWUqCCUUjzw1el8Iz2Z57+3gBHhHeepuCVGBBJoNUsPfVAsug3CU+CdO8DR2v354PMdvq35o6O4/4rpWEyKX729n5zSOs7tokcPxi98V88bFmBl8ogwQuwWLpwc7/Ocr89NJshm5mk/76UXupZi+OeWPE76US/dPXDYvpfnS6vDyRMbjpASFUSrU/Pqzvxev15ueT1mk+L5L45x75v72FdYxZxRpz4FTks0Bujd9eiNLQ6e/yKvRymeA0XVWEyKcXEhnpmUvnrpR8vqUMoIht9eMIpDJbXc++Y+JsaHcvGUEV7nzkiO4Oq5yTy9OddrAPVPHx/GqTW3XzSBAKuZr6cn82FWCZ8eLCUq2EZsqO+ecHfclS7HyutpaHZQWNXoM3/ulhIV5JVycQd3d359ZEQgv70yrdO0qJvJZCwTMtCVLhLQfbEFwfLfwIn9sP3JfnlKpRTfnJ/CWz88lw9uX8z/u2wy18xLOe3nvfvSyTx4ZRoBVt89/fBAK1fOSeLNXQV+3UsvrGzAZjZR1+zg6c25g92cHnOnCNw9wq78Z08RxyrqufvSycxNjeRf24/3uioir7yO8yfFce3CVP6xJY8Whya9zRolo90zIAur0Vrzs9f28PPX9/Lal92/eWQVVTMuLgS7xcxk99R4HzMp88rrGRkeSIDVzFdmjCQ80EpVQws/umC81wQ6tzsumUSI3cJ//3sXv3k3i5/8K5OXM46zZv4oT6XJmvkpOLXms0OlTBoR2m0HqjOeWvSKeo66K1xiO++hp0QFec0udQf3nlTAtTc2NkR66INm0mVGbfqn/wu1J/r1qSfEh3L9eWO8ZpT21bnjY7h0eud5doDbLhhPTIidHzy/w2/r0gsrGxgTG8wlU+N5ZtNRqv1kezB33XZhZdcLqWmteWx9DuPiQrh4SjxXz03haFkdXxzt+eQwp1OTV17P6Jhg7v3KFL45P4XwQKvXpDT3DMi9BVU89lkOr39ZgNmk2NCDmZ9ZRTWe6pKwACspUUE+Sxdzy+sY5QqcgTYzNy4Zw8Kx0ayYNqLDuWCUE/5sxSR251fx9Oe5fHHUqP669fxxnnOSo4JY6iogmNjHdAtAQnggVrMir/xUQO8q5eLuibs3vTnWrofeG+PijGVC6pt79qm/LySgd0YpWPFbaGmAD34x2K05LdEhdv7yzdkUVTby039nntEZfq0OZ7+UTRZWNpIYEcity8ZT3djKc1vyun/QINNae3roxdVd99DXHyzlQHENNy0Zi8mkuHR6AqEBFl7a1vN5BCU1jTS1Oj353fuvmM4Xd19ARJD3YPm0xHB25Vfy4HsH+cqMkXw9PYlN2eVd1oWfrGumuNq7umTqyDCfKZfcsjpS2wTJm5eO44XvL/DZO3e7el4K+355CQd/vZxNa8/nmevmeS1/C/Dtc0YBeNY96guzSZEcGcSxijpPyWLXOXRjbSd3qiWvvJ6oYBuhXYyXdcY9MNrdPsSnQwJ6V2LGw3k/gd0vQeYLg92a0zJnVCQ/v2wyH2Wd4C+fZp+xdZ5veWEnNzy347Sfp7CqgZERgUxPCmfpxFj+/vnRfu/pPPT+Qc5/aD1/+eQwJT4G+3qruLqRqoYWbBYTRd300B/7LIfEiEBWzRwJGD3b1TMTeWdvMZX1zVTUNfPs5lye3nS00zGE3DIj6KRGnwpQvlJxU1xLw85IjuB3V6axZEIctU2tXS6tm+WaXdq2/ntKQhh57TZTrqpv4WR9C6nd5JR9CbZbukylLJsYx9++k87lrp9RX6W4ShePlNWREB7QYQXVttypleOu9V+OV9T3Kd0Cp0oXBzKPLgG9O4vvhNGL4e3boXjPYLfmtFy7MJWVaQn8/sNDTPrFuyz53afc/PwO6poG5iPgiZpGPtxfwsbDpV0G3y055dz9+p5O88X1za1U1reQEGFUEfzw/HFU1DXzUj8vdrThcCmFVQ089MEhFj7wCfeu23taz3fAVRWyYEw05XXNnW5yUt/cyvbcCr42Jwmr+dSf5NXzkmludXLN375g/v0fce+b+/jlW/uZ/5uP+cnLmWSf8J6F6a4uGdVNML1wcjzfnJ/C3749hwCrmUXjorGYVJcLbmW5rqVtQHdvpuw+BpBX4W5D573evlJKcdGU+A7LYPTWqChXQC+t7TJ/DhAbYifAauKYqxTzWEW9Z2C1t1JjgjApyBnAPLoE9O6YLfC1pyAwEv71bWioHOwW9ZlSiv/7+gz+dPVMblk2jrGxIbyzp5jP26y93p/e3lWEU0OLQ5OR23nv728bj/DCF8co7GS/R3f+2b1W/ZxRUcxNjeTvnx/t1w21j1XUc+WcJNb/91LOGRPNqztPb5s/d6/WnfvtrNeffaIWrWFKu2Ujpo4MZ/7oKEqqG/nOOam89+PzeO/H5/H19CTe31vM95713hgst7weq1l51s3vTFSwjfuvmE6ca8mJ0AArc0ZF8tnBzgP6gaJqYkLsXtUlns2U2+TR3TXoXaUxBltKdDC1Ta1kFdV0ux6Me3bpsYp6WhxOCiob+pQ/B6OuPSUqiGzpoQ+ykFi46lmoOg6v3QAO/xiQ88VuMbNqZiI/vXgij31rNnaLia1HBmZP73WZBYyPC8FqVp713dtraHawyfWGsqeTzbTdJYttA9WNi8dSUNnAO3uK+qWt1Y2n6otTY4I5b3wMtU0dtyLsjQNFNSRGBHoG8TobGHWvBzIhvuNg3/PXz2fb3Rfwi5VTmDQijEkjwvj16uncftEE8srrvXatyiuvIzkqCHMXuerOLJkYy/6iap915WC8ObVfp8i9mfKeNmsU5boGGvsa9M4Edzqo2eF7Ua723AG9qLIRh1Of1rWNiwsh54Tk0AdfynxY8SAcfh9euc6vg7qb3WJmzqhIvjjS/8vsHi2rY1d+Fd+Ym8zM5Ai25Pj+FPB5dhlNrnz+nk7KKn0F9PMnxTEuLoTHPzvSLwseuQe9kiONP1b3JJHTyaUfcAXBBNdzdVa6ePhELTaLyWeawmI2YTF3/DOd4dpibffxUz+zvPJ6r/x5b7gXgGufdtmdX8kN/8hgb0E1M5IivI4pZWzm8vGBE54xmdxyIy/dWRnt2aBtSqq7lAucmlzkTif1NYcOxsDo0bK6AVtfSQJ6b8z9HlzyG8h6C/59LbT6zwSXzswfHU1WcXW/lzO+8WUBSsHKtJGcMzaGPQVVnrXb2/o4q4QQu4XxcSGd7upUWNmASUF8m4/7JpPihsVjyCqqZuPh008ZHW9XX+xeAbO4qvezNQGaWh3klNYxaUQYCeHGG1FRJymlg8U1jIsN6VXPeurIMMwmxS7XpxqtNXnldX3uPU5JCCM21O4J6Mcr6rn26W1c/pdNbD1Szo8uGM/Nyzpu67Z6ZiJVDS2sP2iU9p7Om8qZkhQZ5FkxsidL8CZHBlHX7PCsZ97dGEVXxsaF0Oxwesog+5sE9N4652ZY/ls48Da8/G1o8q9NZNubP8ZYV31bbv/10rXWrMss4Jwx0YwID2Dh2Gic2liQrC2nU/PxgRMsmRjLnFGR7Cmo8tnbLqxqJD4soENPdfXMROLD7P2yHICnvtj1x+ruVRf3sYeefaIWh1MzKSGUQJuZiCBr5z30kppe11YH2Yw3wV2uN8Gy2mbqmh19qi4Bo7e9ZEIsGw+X8fL246z400Yyck9y1/JJbFp7Pj+5aILPapBzx8cQFWzzrMdvlCyevekWMCp/RoQFYLOYPFtOdsX9Jvl5dhk2s8nnVpQ9NdCLdElA74sFN8GlD8HhD+CJpVCyb7Bb1GczkyOwWUx80Yc8+qGSGp8TfHbnV5FbXs9q10pzs1IiCLCa2Nwu7bK7oIrSmiYunBzH9KRwKutbfPZcCisbfA702SwmvnfuaDZll5/2LNhjFfVEBFk96/G4/2j7mnJxV7i4NzJOCA/0WbpY3dhCYVUj4+N7t1kDGP93u/MrPb1zgFGnMRi5ZEIsVQ0t3PnqbqaMDOPd287jB0vHdllzbTWbWJmWwEdZJRRWNlBe1zwgFS79bXRMMGNignv0qcj9Jr8zr5KkyMA+jVG4uWvRB2pgVAJ6X837Pnz3LaOH/rcL4Mt/gh9uzhxgNTMrOaJXMxIBNhwqZcWfNvL7Dw51OPZGZgE2s4lLXDMD7RYzc1Oj2NJuYPTjrBJMylhHfnqi701AoPOADvD19GSjPYdPb4/LYxXe1QsB1q571d05UFyN3WLy9JgTwgN8plwOu9YvmehjQLQ7aUkRVNa3cKyi3lNdcjrpjiUTY5k3Ooq7L53Ei99f0ONc8aqZiTS1GuvQnG4bzpT7Lp/KQ1fN6NG57nGVZoez2zVbuhMeaGXBmCiCO1mU73RJQD8dqefCjRshKR3W3QLPfsUve+sLxkSzr7Cqx9PpDxRXc/PzO3E4Ndt9pGo+3F/C4gmxXksbnDM22rOHqttHWSdIT40iMtjGxBGhWM2qw8Co1prCqkZGdrKSXUSQjaTIwNPezd7XhJERYQF9zqEfKDbSKO40kRHQO745HHKteuirwqU7M5KNN8Fd+VXklddhNilPaWdfhAVYefnGc7hh8dhe9UJnp0SQEhXk2YTibE+5gPHzdi9U1p1Am9lTrtkf1Tsv3XAO3z4n9bSfxxcJ6KcrNB6+sw4u+z8o2QuPnwv/+W+o6v1KeYNl/pgonBoyXME5p7SWa57Yyo68jrXjJdWNXPf0doLtZq6ak0RWUbXXxKSCygbyTzawaFy01+MWuja+dpdIFlQ2kFVUzYWTjd2Y7BYzk0aEsaeg0utx5XXNNLc6u6ytdm/m21cOpyb/ZH2HP9b4sIA+p1yyimq8psmPjAjkZH1Lh1UNDxbXEOTaxrC3JsSHYreY2HW8ktzyekZGBPR6W7b+oJRi1cyRnkqXUVFnfw+9t9y/G2dzOSZIQO8fJjPMvR5+uBPS/wsynoI/zYBXr4fCLwe7dd2anRKJzWziiyMVlNc2cd3T29lypJwfvrCTyvpTlTxVDS1c9/R2qhtaeOrauVyWloBT47Wb+XZX6mZuu407po0MI9RuYV1mIS9vP87972QBcEGbZX+nJ4WzO997YNRXyWJ7U0eGc7S8rs8zXkuqG2lxdKwvHhEW0KdB0dKaJspqmzz5c/dzQcdB1sMnahgfF9LlOiedsZpNTEs0dsI6Vl43qKkO9848I8ICOl3j358luwZPT6dk8UyQgN6fgqKMnvptmTDvRjj4njFo+ug5sOEhqDgy2C30KcBqZmZyBBsOl3H9PzIoqW7k16unUVrbxJ2v7EZrTXVjC9/5+xccPlHDI2tmM3VkOLNcm3a07clvy60g1G7psOejxWxi4bhoPtxfwp2v7ubdPUUsnhDrGSQCmJ4YTk1ja4dd2eFU1YkvUxLC0O12s291OLn0Txt7tLhVZyvojQgPoKy2qdezUd3VPO41wwHPsgVFld5pl4PFtX1Kt7ilJYWzp6CKI2V1p1VOd7rGxYUwKyXC5ybpQ4H7d2Mwf8Y90fmqNKLvIlJg+f2w9C7Y/TLseQU++R/jK3ocjFlqfCXNM1I2Z4H5Y6L48yfZKAWPfnM2K6Yn0Nji4Nf/yeKxz3L4YF8J+4uqeWzNHJZONNIk4YFWxseFeC3qtO1oBempkT5zsPdfMZ3vLkwlKSKIhIgAr3VLgFMDowVVntX62k/798W9psi+wmrPZg57CqrYX1TNM5tzubqbdeePtZtU5DYiPACt4URNU69SIi9nHGdkeADpbT6ljHTVordd3qCirpmy2qbTWg52ZnIET2/KpbHFOeiDkc9cOw/6XgByVlsyMZbtuSfP6iUNoIcBXSm1HPgTYAae1Fo/0O64ch2/FKgHrtVa93yX5aEqINyohpn3fag8bmxrd2Q9ZL54auOM0JEwchbETYaYCcYKj5GpxtoxfVzEvy8WT4jlz59k87MVk1jhWl/9vxaNZlN2GQ++dxCLSfHomtlcOMX7DWjOqEje3VuM06k5Wd9M9olavjq743Z4YCzju7CTPRfByAnbLCb2FlRx+QxjRb3CygYCXRUnnRkRFkBUsI19bTY//tw12ehAcQ0Hiqu90h/tHa8wdvlx96LbPi8Y26H1NKAXVjaw4XApP1w2zutNzT3ztLjNwKh7h57xp9VDj/D8e7DLBcO7+D/yd3NGRfHiDQsGuxnd6jagK6XMwCPARUA+sF0p9abWen+b01YA411f84HHXN+FW0QynHOL8dXaDIU7oWCnkWMv/BIOvQe6zYCZLQTCkyA0AULijfVkgmKMQB8UZbxZ2EONDa1tIWANBFuwkc/vg7mpUWxee75XrtpkUjx01QzueGU318xL4aIpHT9NzB4VyUvbj3OkrJZs1xoV80f73vi6OzaLickJYexus6aLsWxuQJfLqiqlmJIQ5rV7zsbsMlKigiiobGBdZiGTlnce0I9VGAOK7T8x9KUW/dUd+WgNV7nKKd0CrGaigm1ePfRDp1Gy6JYaHURYgIXqxtY+TyoSQ0dPeujzgGyt9REApdRLwCqgbUBfBfxDG6NZW5VSEUqpBK11/6ycNNRYbJCywPhya22Gk0eh7BCczDOqZKqOQ00xlOdAbQk4elBCZ7aBJcD1ZQez1bjPbAWT9dR3kxlMFuO7MoPJzEiljH8rk+crWimeCjfBQQWHTK6P1Mr16UFxSWMrLZZiePsNohtaud9Ww8zdH8PedsMzXgFZdXr/XbqMnII6nO++gQnF8vxCLjeb4N33fJx/yo9aKth3ogrHu2/jdGpW5OcxdWQ4J83NVG5vQeskz5uCe9DVffuCvEKWmxW896HXc45pcfD/LMdI2f4OFHRf4qa1Jm7HcR6NsZK8bVOH47+w5BN0xALvGfX5Y3PK+FVALfFbPqevuQoF/DakiILWBkbv2AgmGRbzC6MXw8Tl/f60qruFjZRSVwLLtdbXu25/G5ivtb61zTlvAw9orT933f4YuEtrndHuuW4AbgBISUmZk5d39u86c9bQGlrqob4CGiqMCU1tv1rqjd2VmuugtQlaG43vjmbXVws4W1zfW8HpcH1vNT4ZaG3cp52uL9d92umaMKVP3Xb/G43Wmoq6JuxmRatTY1IQZm//KUF7X4fP+41vLU4nDc0O7BYTdouZmqYWLCYTgVaz9/nttDicNLQ4CLZb0FpT3+wgyGZGa2hocRBks2AxKc8xk0m5nhNqm1qxmFSHBaW01tQ1tWI1qx6twd2qNQ3NDgKsJqw+AmtDiwOn1gS7ptDXtzgATZD19IayWpxOHE5NwGmuEy7OoHNuhmV39+mhSqkdWut0X8d68pvkq+vQ/i+rJ+egtX4CeAIgPT3d/6ZVDialjJSKLdhI35wlFHDHM9s5WFxDUVUDt54/np9cNKHPz2fRmrtf/JL39hbz4g0LuOrxLdx+4QRuu3B8l4/LO1HDhb/fwEOXz2B/YTX//CKP3XdeTKtTM+/XH3LV1GTuu3wq33t2O+sPlqIUfHDjYhIjA5l2z/vcuXwiNy8d5/WcCrj0wU+ZlRLBn66e1W3bf/rSl3x64ATb7roQq4/VBu9/Yy9v7ipk190X0+pwcs6vP+LS6Qn85qvTe/Uzas/q+hKiJwE9H2gbQZKAwj6cI4ao2SkRfHLAWG2vr/lzN6UU/7t6OjvzTnKTa+u69oOVvoyOCSHAamJ/YTWbssuYlxrl6XFfODme/+wpItBmZv3BUu64ZCJ/+SSbR9fncNMSYwXBziaMjAjznrK/t6CKv244QmJEIGNigkmICKDF4aS2ycG7e4u5em5yp0vHJkQEUNXQQn1zK7999wBVDS0snRjbq5+PEF3pSUDfDoxXSo0GCoCrgW+2O+dN4FZXfn0+UCX58+Fj9iijHt1iUsxKiTjt5wsPsvL7b8zkmr9tBbouWXQzmxSTRoTx2aET5JTWcUWbSptVMxN5e3cRT2w4wpr5KdyybByV9c38/fOjnlrxzgJ6fHiA1yDt05tyeWdPESZl7MTUlknBNV2USLpLFx96/xDPbsnj+nNHc8nUEd1emxA91W1A11q3KqVuBd7HKFt8Smu9Tyl1k+v448A7GCWL2Rhli9cNXJPF2WZGUgRmk2JaYniXG+72xoIx0fxgyVge+yzHawf5rkwdGcbzXxgTic4dF+O5f8mEWOJC7YyKDuLer0wF4PvnjeHZLXn88aPDQMcadLcRYXY+qGpEa43WxgYQl05P4A9fn0FBZQPFVY3YrWYCrWYig63EhXb+acJduvjUpqMsmxjLzy6d3KPrEqKnevTXp7V+ByNot73v8Tb/1sAt/ds04S+C7RauW5jKlJGdlwb2xR2XTOQbc5N7XAPu3uMyKtjGlDYzVW0WE+/9eDGhARZPaWJcWABXz03mH1vyCLVbOq1zjw8LoKnVSVWDsbRvWW0TSybEYjEbOwz1pvbbfR3j40J4+JpZp7UMqxC+yExR0S/+38op/f6cSqleBUz3G8rCsdEd1kaJCrZ1OP/GJWN5cdsxkqOCOq1zH9Fmowv3bj6LJ8T4PLc7SZGB3PuVKVw8dUSXa4wL0VcS0MWQMWlEKNMTw/na7KQenZ8YEcjPL53c5WJSp/YDbeSzg6VMHRnWZVqlK0oprls0uk+PFaInJKCLISPAauatH57bq8dc202Adc8WzS6pZcexk9y4eEyf2yfEQJNpZUJ0wd0bf3VnPg6n9ixMJsTZSAK6EF2wWUzEhNg4UFxDqN3SL2WZQgwUCehCdMOddjl3fEyHBbyEOJvIb6cQ3XAvo7tkgszqFGc3CehCdMNdurhEpumLs5xUuQjRjW/MTSYlKoiE8N5v5CzEmSQBXYhupCVFeO0MJMTZSlIuQggxREhAF0KIIUICuhBCDBES0IUQYoiQgC6EEEOEBHQhhBgiJKALIcQQIQFdCCGGCGXsHjcIL6xUKZDXx4fHAGX92Bx/MRyvezheMwzP6x6O1wy9v+5RWmuf61AMWkA/HUqpDK11+mC340wbjtc9HK8Zhud1D8drhv69bkm5CCHEECEBXQghhgh/DehPDHYDBslwvO7heM0wPK97OF4z9ON1+2UOXQghREf+2kMXQgjRjgR0IYQYIvwuoCulliulDiqlspVSawe7PQNBKZWslPpUKZWllNqnlLrNdX+UUupDpdRh1/fIwW5rf1NKmZVSXyql3nbdHg7XHKGUekUpdcD1f37OMLnu212/33uVUi8qpQKG2nUrpZ5SSp1QSu1tc1+n16iU+pkrth1USl3S29fzq4CulDIDjwArgCnANUqpKYPbqgHRCvxUaz0ZWADc4rrOtcDHWuvxwMeu20PNbUBWm9vD4Zr/BLyntZ4EzMC4/iF93UqpROBHQLrWehpgBq5m6F33M8Dydvf5vEbX3/jVwFTXYx51xbwe86uADswDsrXWR7TWzcBLwKpBblO/01oXaa13uv5dg/EHnohxrc+6TnsWWD0oDRwgSqkk4DLgyTZ3D/VrDgMWA38H0Fo3a60rGeLX7WIBApVSFiAIKGSIXbfWegNQ0e7uzq5xFfCS1rpJa30UyMaIeT3mbwE9ETje5na+674hSymVCswCvgDitdZFYAR9IG4QmzYQ/gjcCTjb3DfUr3kMUAo87Uo1PamUCmaIX7fWugB4CDgGFAFVWusPGOLX7dLZNZ52fPO3gK583Ddk6y6VUiHAq8CPtdbVg92egaSUWgmc0FrvGOy2nGEWYDbwmNZ6FlCH/6cZuuXKG68CRgMjgWCl1LcGt1WD7rTjm78F9Hwguc3tJIyPaUOOUsqKEcyf11q/5rq7RCmV4DqeAJwYrPYNgEXA5UqpXIxU2vlKqX8ytK8ZjN/pfK31F67br2AE+KF+3RcCR7XWpVrrFuA1YCFD/7qh82s87fjmbwF9OzBeKTVaKWXDGEB4c5Db1O+UUgojp5qltf59m0NvAt91/fu7wLoz3baBorX+mdY6SWudivH/+onW+lsM4WsG0FoXA8eVUhNdd10A7GeIXzdGqmWBUirI9ft+AcZY0VC/buj8Gt8ErlZK2ZVSo4HxwLZePbPW2q++gEuBQ0AO8PPBbs8AXeO5GB+1dgOZrq9LgWiMUfHDru9Rg93WAbr+pcDbrn8P+WsGZgIZrv/vN4DIYXLdvwQOAHuB5wD7ULtu4EWMMYIWjB7497q6RuDnrth2EFjR29eTqf9CCDFE+FvKRQghRCckoAshxBAhAV0IIYYICehCCDFESEAXQoghQgK6EEIMERLQhRBiiPj/OBM/PGMS4yAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_data = np.arange(100)\n",
    "plt.plot(x_data,acf[0:100], label='autocorrelation') \n",
    "plt.plot(exponential(x_data, -1/tau), label = 'tau_calc estimation')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Key point*: If the data fits an exponential well, we can treat it as uncorrelated samples if they samples are $2\\tau$ apart. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-35179.34066243115\n"
     ]
    }
   ],
   "source": [
    "start = 500\n",
    "mean = np.mean(potential[start:])\n",
    "print(mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nsamples = (len(potential)-start)/(2*tau)\n",
    "stderr_mean = np.std(potential[start:])/np.sqrt(nsamples-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What should we use instead of `nsamples` above?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean = -35179.34066243115, sderr_mean = 16.580879179923137\n"
     ]
    }
   ],
   "source": [
    "print(f\"Mean = {mean}, sderr_mean = {stderr_mean}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we are probably interested in is the average potential energy per molecule, not the total potential energy. There were 900 water molecules in the simulation above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can estimate the heat of vaporization $H_{vap}$ by:\n",
    "    \n",
    "\\begin{eqnarray}\n",
    "H_{vap} &=& H_{gas}-H_{liquid} \\\\\n",
    "        &=& U_{gas} + PV_{gas} - (U_{liquid} + PV_{liquid}) \\\\ \n",
    "        &=& \\left(U_{gas} - U_{liquid}\\right) - P \\left(V_{gas}-V_{liquid}\\right) \\\\ \n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the kinetic energy of liquid and vapor is the same, then this is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray}\n",
    "U_{pot,gas} - U_{pot,liquid} + P(V_{gas}- V_{liquid})\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assume ideality of gas ($PV=nRT$), and zero internal energy (which is valid for rigid water, like TIP3P or SPC/E), then we get:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{eqnarray}\n",
    "       &=& - U_{pot,liquid} + P\\left( \\frac{nRT}{P} - V_{liquid}\\right) \\\\\n",
    "        &=& - U_{pot,liquid} + nRT - PV_{liquid} \\\\\n",
    "\\mathrm{Molar} &=& -\\frac{U_{pot,liquid}}{N} + RT - \\frac{PV_{liquid}}{N}\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "What is the $H_{vap}$ of water predicted by this simulation? What is the uncertainty in the estimate? How does it compare to the experimental $H_{vap}$ of water is 40.7 kJ/mol?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hvap = 41.581 +/- 0.018 kJ/mol\n"
     ]
    }
   ],
   "source": [
    "print(f\"Hvap = {-mean/900 + (300)*(0.00831)-(8.31)*(18.02)*(0.101) /(900*1000):.3f} +/- {stderr_mean/900:.3f} kJ/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the $\\frac{PV_{liquid}}{N}$ term is almost zero; we actually hould have thrown that out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0222444444444443e-06"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.101)*(18.02)/(900*1000)  # 18.02 g/mol x 1 L*cm^3 / 1000 g  x 1 atm x 0.101  kJ / L*atm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automated tool for equilibration detection and correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use scipy statistical tests to see if two parts of the distributions are within the uncertainties of each other.  NOTE: have to use uncorrelated samples, or it will erroneously say that they are NOT within uncertainties of each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A higher T-statistic means more difference between the data sets\n",
    "\n",
    "A lower P-value indicates a low probability that the difference in means was by chance\n",
    "\n",
    "P-value < 0.05 suggest they are two different data sets! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise:  Use this code to estimate better than \"eyeballing\" what fraction is equibrated  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 613,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=-14.984452443412309, pvalue=6.184533203981602e-46)"
      ]
     },
     "execution_count": 613,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start = 0\n",
    "cut = 500\n",
    "stop = 1000\n",
    "per_ind_sample = 1  # the number of indices between independent samples\n",
    "# compare \n",
    "scipy.stats.ttest_ind(potential[start:cut:per_ind_sample],potential[cut:stop:per_ind_sample],equal_var=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definitely not independent!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=2.025269013602619, pvalue=0.05124877861250221)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start = 600\n",
    "cut = 800\n",
    "stop = 1000\n",
    "per_ind_sample = 12  # What SHOULD we use here?\n",
    "# compare \n",
    "scipy.stats.ttest_ind(potential[start:cut:per_ind_sample],potential[cut:stop:per_ind_sample],equal_var=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pymbar.timeseries tools "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymbar import timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "[t0, g, Neff_max] = timeseries.detectEquilibration(potential)\n",
    "# t0 is the initial point detected as starting the stationary point.\n",
    "# g is the estimate of the correlation time (approximately 2tau)\n",
    "# Neff_max is an estimate of the _effective_ number of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "745\n",
      "5.662475\n",
      "45.386513\n"
     ]
    }
   ],
   "source": [
    "print(t0)\n",
    "print(g)\n",
    "print(Neff_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can subsample the data, picking every $n$th particle so that all the points are statistically independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "84"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newdata = timeseries.subsampleCorrelatedData(potential[500:],g=6)\n",
    "len(potential[500:])\n",
    "len(newdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}