{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Chapter 6: Bayesian Model fitting II"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise 1 (20 points)\n",
    "\n",
    "You have a set of $n$ data points, $\\{x_i\\}$, which have been drawn from a Gaussian distribution of known mean $\\mu$, and unknown variance, $\\sigma^2$\n",
    "We want to find the posterior PDF over the variance. Show that when you adopt the following prior for the variance\n",
    "\n",
    "$$p(\\sigma^2 | \\alpha, \\beta) = \\frac{1}{Z} \\sigma^{-2 \\alpha - 2} \\exp\\left(-\\frac{\\beta}{\\sigma^2}\\right)$$\n",
    "\n",
    "the posterior is conjugate with this prior ($\\alpha > 0$ and $\\beta>0$). Write down all of your steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercice 2 (35 points)\n",
    "\n",
    "We revisit the exercise of the lighthouse, but this time we will completely solve the problem in 2 dimensions.\n",
    "\n",
    "**Recall of the problem**: A lighthouse is somewhere off a piece of straight coastline at a position $x$ along the shore and a distance $y$ out at sea. It emits short flashes (perfectly collimated) at random times and hence at random azimuths, $\\theta$, from its position, _i.e._, $p(\\theta) = constant$. You see these flashes on the coast while walking along the shore. You record your position $D_k$ along the coastline at the instant you see a flash, but you do not record the direction the flash came from.\n",
    "\n",
    "* the dataset of $\\{D_k\\}$ is in the file `lighthouse.dat`. (Same dataset as before) \n",
    "* All distances are in $km$.\n",
    "\n",
    "We assume (prior information) that the lighthouse is constrained to be located within the rectangle $-2 < x < 2$ , $0 < y < 2$ , with uniform prior, where $x=0$ coincides with $D=0$.\n",
    "\n",
    "Below is a schematic to help you understand the situation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzAAAAFqCAYAAADFgr+jAAAAiXpUWHRSYXcgcHJvZmlsZSB0eXBl\nIGV4aWYAAHjaVY7rDcQwCIP/Z4qOQICYME5VNVI3uPEPmj7uPkVgWZFx2T/HKEtSiYs263CAAnV1\nXkN0mghRZaq5Y06uLTUUv3YRngLejfT9qJd/0wQdw9QMDRs2jnTehYVjZlCmUtbwn2vnS/HvY30K\nzBZfQzcsSfEGfW8AAAoGaVRYdFhNTDpjb20uYWRvYmUueG1wAAAAAAA8P3hwYWNrZXQgYmVnaW49\nIu+7vyIgaWQ9Ilc1TTBNcENlaGlIenJlU3pOVGN6a2M5ZCI/Pgo8eDp4bXBtZXRhIHhtbG5zOng9\nImFkb2JlOm5zOm1ldGEvIiB4OnhtcHRrPSJYTVAgQ29yZSA0LjQuMC1FeGl2MiI+CiA8cmRmOlJE\nRiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMi\nPgogIDxyZGY6RGVzY3JpcHRpb24gcmRmOmFib3V0PSIiCiAgICB4bWxuczpleGlmPSJodHRwOi8v\nbnMuYWRvYmUuY29tL2V4aWYvMS4wLyIKICAgIHhtbG5zOnRpZmY9Imh0dHA6Ly9ucy5hZG9iZS5j\nb20vdGlmZi8xLjAvIgogICBleGlmOlBpeGVsWERpbWVuc2lvbj0iODE2IgogICBleGlmOlBpeGVs\nWURpbWVuc2lvbj0iMzYyIgogICB0aWZmOkltYWdlV2lkdGg9IjgxNiIKICAgdGlmZjpJbWFnZUhl\naWdodD0iMzYyIgogICB0aWZmOk9yaWVudGF0aW9uPSIxIi8+CiA8L3JkZjpSREY+CjwveDp4bXBt\nZXRhPgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAogICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAog\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgCiAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKICAgICAgICAgICAgICAgICAgICAgICAgICAg\nCjw/eHBhY2tldCBlbmQ9InciPz4uWXaiAAAABHNCSVQICAgIfAhkiAAAIABJREFUeNrs3WdAFNff\nhuFn6SCCoqDYC/beS+y995bExMRoNJr4TzPFGKMmpve8Eo2JiVGj2HvvvWPHXlARRUF62YV9P5ig\nhkUhIku5r0/DnrM7M7+dXebZOTNjWOt/3SwAAAAAyAZsKAEAAAAAAgwAAAAAEGAAAAAAEGAAAAAA\ngAADAAAAAAQYAAAAAAQYAAAAACDAAAAAAAABBgAAAEDuYEcJkJW8PWyAju3bSiEAAMATU71uY301\nbSGFIMAAj+/Yvq3yXX2WQgAAgCfCzcVeA5uXohDZGEPIAAAAABBgAAAAAIAAAwAAACDX4hwY/Cfx\n8XFau2SOdm5apYtnTio6KlJ58rqplE9FNWrRXp16DZSjkzOFAgAAAAEG1nXmxGF98s7LuhF05YHH\nw8Nu68j+nTqyf6fm/T5ZH3wxVVVrN6BgAAAAyDAMIUO6XDhzQu8M65sivPxb6K0bendYXx3et4Oi\nAQAAIMNwBAZplpSUqC/HjlJsdJQqVauj7k8PVrXajZSvQEFFht+R/97t+mvad7py6ZwkyWQyatI7\nL2v60l3K656PAgIAAOCxcQQGabZ781pdPHtS/V58Vd//uUItO/ZSwULesrOzV/4CnmrVqZcmz12v\nBs3aJj8nIjxMi2f/QvEAAABAgEHm2rxmsbr2e0Evjfog1T6Ojk4a8/kUeRcrmfzYri1rKR4AAAAI\nMMhcgRfPatCIdx/Zz8nZRX0HjUz++/rVSxQPAAAAGYJzYJBmvyzYkua+dRo1T562d3CgeAAAAMgQ\nHIHBE1HQyzt5upB3MQoCAACADMERmEf49L3h2rp2aartPZ8ZquGjJz7w2MwpX2vW1G8eeKx5++4a\n8/mU/7QM7Wt5Z/h6rfW//kTrFhUZnjxds35TNiQAAABkCI7APMKoMV9o5LuTUgyDqtu4pf5a558i\nvEjSwGFvafaaQypRprzsHRw0+uMf9daE73NV3YKuXEyebtOlLxsSAAAAMgRHYB7B1c1d3QYMlr2j\no76f+Hby483bdVMBz8IWn2MwGFSwkLc8C3mrbZe+j70D/6SPljwJ+7ZvTA56pctVYkMCkOsZE+Ll\nv3Otdq2dr7NH92ryqjMP7b/096+1ZflMtej6nLq/+HYmL2uc/Hf8vazH9j1yWXPa/AEQYHKEjj2f\n1da1y+S/d5sk6cTh/WrXfUCq/ePjYnX+9AmN/WparqtVYqJJm9cslo2NrV763wdsPABytSvnT2rX\n2vnat2mZYqMj0vy8zcv+VEJcrLYsn5lpAebK+ZPauWae9m9enq5lzSnzB5A9MIQsHQaNGH3vH8ua\nxYqKCE+1785Nq1WnUQu55Mmb6+q0ceVCBV8LVPenB6tM+SpsOABytXk/f6yb1y7JxsaQrue17D5I\nDk7OatHt+Uxb1jk/faiQoMsyGAxWqZW15w8ge+AITDpUql5X5avU1JkThxUfF6vVi2er76ARqezE\nz1evgcNzXY1ioiP1589fqmiJMnph5HtsNAByvbe+niNJOn14l354f1Can9f9hbfU/YW3MnVZ3/l+\noSQp4NAO/fTBi5leK2vPH0D2wBGYdOr5zNDk6WV+vyspKTFFn7DbIbp84YxqNWiS6+rz6/cfK+x2\niN7//Gc5ObuwwQDA30qUq5qNlrXaE3ndG1cvakTHclabP4CcgSMw6dSsXVdN+26iQm/d0M3rV7Vz\n02o1bdPlgT5b1ixRs7ZdZWNjmyHzzC6XUd61eY1WLZyl/439SuUqVWdjAYD7ODhmnx91nJzzPJHX\n9d+xxqrzB5AzcAQmvYnPzl5d+t4bj7zkr5Qn6W9YOV+tO+euSwdfPn9aX304Sl37v6iOvZ5lQwGA\nf//DtbXN1ctqNpu1b9OSHFcrAFbYH6cE6de5z/Oa89sPMiYk6Lj/Pp0LOCafSncPdwdePCuT0aiy\nFTLu5PWsfhnl0Fs3NPa1garXpJVGjP6YDQQALMhOJ6Y/iWXdtXaegq9cyHG1ApD5OALzH+TzKKgW\n7Xsk/734vqMwG1csUKtOvXJNLSLD72jMiKdVonQ5vfPJTzLYsEkBAB505dwJLZz2OYUAQICxph73\nncy/Ze0Shd0Okdls1pa1S9SqY+9cUYPoqAiNGfm0XPO6a9y302VnZ8+GAQB4wJHd6/XD+4MUFxNF\nMQBkCIaQ/Uc+Fauqau0GOn5or0xGo1bMn6Ga9ZvIy7uYPAsXyRXh5f1XBshgMOjjH2fK0dHpof3N\nZrOCr12Wd7FSbDwA8AjGhHgd2bVeO9fO15kju9N1J/rw0BBtXT5TJw9uV0hQoIwJcXL38FR+T29V\nqddC9Vp0lYfXf/s/ZTab5b9jjXavW6jLZ4/JmBCnQsXKqHmXgWrULuWPdz+OGaRT/rtSPP7vK5G9\n+vFvqly3WYbPP7XXOHlwmw5sWaELAf66cytYBhsbueUrqOI+lVW9YRvVad451R/lUruKmu/qs4/V\n935XLwRo05I/dPboXoWHhsjJOY+K+1RR9YatVb1ha127dFq+44Y+9HXiYqK0a90Cndi/VUGXzyg6\nPEz2jk4qWLi4KtZqrOZdn/vP2wFAgMnGejw9RMcP7ZUkrVzwp4KvBapN5z45fr0jwsP0/iv9ZZBB\nn0+dL+c8ro98zozJX+iY/x5989sSNhwASMWV8ye1e90C7du0TDFR4el+/q51C+Q3eYKKlamovsPH\nqmS5aoqOuKM9GxZpxawfdO74AS39/esUz/PwKqpPZmx56GtH3rmt3798M0UgCTx7XDO/e09Bl8+o\n99D3H2gb9ekMizvzj9qBz6j5/9u1i6c067sxSkpKVOveL6nH4HeUJ6+7wm4F6/i+LVo3b4r8d6zV\nqtk/acCrE1SpdsrbIfzfilOKuHNLa+b+rG0rZj90fj8uP6mwkOtaN+8X7Vzjl6b13LTkD63+6//U\nYcAIdRk4Sm4eXoqJvKOzx/Zpzdwp8vOd8MjX2LN+kRZO+1TlqjVQ1+dfV9HSFRUdeUeHtq/Wylk/\nav2CX7V5yQz1fOldtewxiA8esh2GkD2Gp1p1lFfhopLu3vtl+/rlatq2S45e51s3g/XW4B6ytbXT\n51PmyTWvW4o+ZrNZxoQEhd0O0YFdWzT+9Rc057cf1LZrPzYaAHiIOT+NU9ClM0pMNKb7udtWzNas\n796XZ5GSev2LWSpbuY7s7B3kXsBL7fsP18DXP0vuazDY6APf5fJdfVa+q88+MrwEXzmvL1/vLXcP\nL7351V/6btERfTl3r3oNeTf5imEbF01X4NnjT6QuGTH/M0f26Ou3+suzSEmN/n6B6rfsJncPT9nZ\nO8jTu4Radn9eH05drTKVayvkeqAmfzhEu9ctTLnjZGurfAUKqcfg0Y9cbjs7e3l6l1CPF99O03oe\n37dFC6ZO0ovvfqfWvQbLw6uo7Ozs5ZbfU3Wadda7Py565JGqpb9/rT+/fVctug/Syx9OVqkKNWTv\n4Kh8BQqpVY8X9Pa38+To7CKTyaj5Uz/R+vnT+OAh2+EIzOOkPxtbde3/on774RNJUoPm7eSSJ2+O\nXd+rly/o/Vf66+b1q5Kk3s0rpvm5jo5Oata2KxsNADzEO98vkCSdPLhd/zd2cJqfd+PqRc2fevd/\nUYcBw2Xv4JiiT4PWPbRzjZ/OHT8gszlJW5fP1jOj0nblyCkTXtEzoyapUu2n7n2vO7uoTe8hiomK\n1Jq5vskhauAbn2V4XR53/rdvXNPUj0coT958euZ/n6Q6PMzF1V3DP5qiSSO6KPz2Tc36fowKFCqq\n8jUapuibnnvV5HHLl6Z+a+dNlSSVLF891UD0wttf6aOX2ig2OjJF+54Ni7R23lSVr95AnZ99zeJr\neJfwUcenR2rJ9K/uBp4/vlG1hq1VuHgZPoDIPvvglODxdOo1UI5OzpKk1p1y7sn7AccO6s0XuiaH\nl/Rq3KpTjg53AJCRSqbzTvQbFv6qRJNJkuRTtV7q38Xt7t2j7MzRvWl+/eEf/fxAeLhf7aYdkqcv\nnTn2ROrxuPP38x2v2OhINWzT65HBw9Utv54eeXeYltmcpD+/fVcmkzFT3vcr505Iknasnpv68rl7\nqEmnASkeT4iL1aK/r/TWovvzD51P1XotkqeTkhK1ZekMPnQgwOQmrm7uatOlj9zzeajuUy1z7Hq+\n+3Jfhd8J/c/Pb9OlLxsLAKSRo0v67kQfcGhn8rRbvoKp9itVsWbydFhI2u8xVrh42VTb8nveOxE8\nPPTGE6nH48w/8OxxHd+3RZJUrUHa/k9Xa9hahYqVliSF3gzS/k3LMud9d3aRJC3741st/eMbGRPi\nLPZr2Lpnisf2blqqqIiwuyG2St2Hzqdg4eIP/H3q8G4+dMhWGEKWAUqVrShbW7scfRnhZbsv8EYD\nQCaxtU3fv+eIsJDkaZMpQQ62zhb75StY6IntdEtSfGxM5oe9R8x/78bF94Ud7zS9psFgUN3mXbRy\n9k+SpMO71qX5KmePw6dqXfnvWCuzOUlr/aZoz/qFatVzsJp07C/n+0YxeJcsl+JCCAGHtidPvzOg\nQbrme+dWMB86EGBym02rF2v42xMoBAAgQ6T3TvR53PIp/PZNSXcvo+zpXcJiPwcHp3TvzKcnbCVm\n0lCr9Mz/9JF7Q+Xy5M2f5tctU7lO8nTgueOZsi7dX3hbp/x3KzY6Ivm9XPzbF1r91/+pUfs+atl9\nUIqjJ/+4cv5k8vQPS49bPA8KyCkYQvaYrl+9pMjwUFWsVptiAACsoniZysnTF08dTrVfVMS9ocAV\nazXO8LBlNputGvYszf/W9cDkaaMxPs2vW+i++5ZFh4dlyrp4FS2ld75foFIVajzweFxstDYvmaHx\nL7XV71++qdCbQSmee/8yxkSG86EAAQap27hyoVrl4JP3AQBZX9POTydP71m/MNV+F0763/3nb2Or\nVrnk/h9JSYn3Alx42s/ldMl778ph9o+4WXNGKlSstEZ/N1+D3/1O3iXLpViX/ZuX65PhnXR0z8YH\n2kzGhOTpqxcD+FCAAAPLzElJ2rRqkVrngptXAgCyrmoNWqlu87v3ITvlv0sHtq5IuSOfmKh18+9e\nprfroDfkVbR0rqhN3nwFkqevnDuZ5ufZ2zs+ECoyk8FgUN0WXTT255Ua/tEUi0dkfvlk5ANH2+4P\nXCcPbONDAQIMLNu8Zomq12mkwkVLUAwAgFU99+Znati2lyTpz6/f0dp5UxUeGiKTyahLp4/qxw9e\n0OUzx9Xluf+pfb9huaYupe+78trpI2m/2lZcbHTydIWajVMJGk92N8pgMKh6w9Z65/sFGv7RFBUo\nXOyBQLpsxnfJf99/H5fd6xYqmmFkIMDkbrOmfqPezStqQJvqmvrNR7oRdEUHd2/VollTNeSNcRQI\nAGB19g5OatXjReX39FaVei20e91CjXuxld7qXVu/f/mmPL1L6gPfZer0zKu5qi61m3ZMnj64beUD\nQ60e5lZwYHKIaGDhssWSku8DJyld94pJSky0+PiIjuV089pFi23VG7bWB5OXq0zle+fcXjh5KHm6\nYq2nHghfs757P03nJB3ds1HfvfMsHyBkK1yFLA2WzPlVURF3f8lYNOsXLZr1i0qUKa+Pf5wp17xu\nFAgAYHV7NizSnJ/G6fm3vlCdZp2z3PIZDDYym5MkSYkmk2ztMmcXpOZT7VSoWGnduHpR0RF3tGXZ\nn2rTe8gjn/fP+UK1mnRI9S71bh5eivs7cNy5FZzqFcKCLp1+4O/4uJgHLot8v32blqrLc69bbHNy\ncdXTr07QpBFdJUkOjveGuTVu30er5/gm3zvmyO71mvHV23r29U9TvSLZtYunNPO797Lk9gI8DEdg\n0mD0xB9VonQ52dnbq1ipsho8aowmz1nH0DEAQJoZEx68Atb9J5c/Tl9J2rd5mf785l0ZE+Ll7JJX\niSbTYy3rv48Q/BM8LPn3r/ypHeFwcnFNnr4VfCXT5m9jY6vn3vxCNja2kqRVf03WjasXHz7/pETt\nXOMnV7f86vPymFT7FS977+pvh3eus9jn6O4NmjpxxAOPBZ49luprbloyQ6E3r6XaXrDwvX2P8jUa\n3gtT+T3VYcArKbaLiS+317YVs3U7+KpMJqOiIsJ0/uRBzZ/6ib58va9sbe1y3VE5ZH8cgUmDBs3a\nqkGzthQCAPCf/fvSt+G3b6Z6L5Z/971zK1geXkVTfe0l079Knv6/D1965LLkK1hYdZp1UueBo+Tk\nnCdFe9i/bmwYHRkuVzfL91CJiYr417Jfs3iBgCKlyuv8iQOSpGUzvtUzoz6RKSFeGxdP163rV/Ty\nh5Of2PzLVKql59/+Un9+/Y7iYqLk+9FQvfrxb/IsUtJiIFo47TOF3gzSyI9/Vb6ChVOtY9NOA3Rw\n20pJ0vI/v5M5KUn1WnaVi6u7Lp0+oi3LZuraxVMa+fFvGj/k3n7Er5+OUulKtRQTGa63v533wGvG\nxURp8riheu2T6RbnfWzfJkmSg6Ozuj7/xgNtHQa8ouAr57R/8/Lkx27fuKa5k8enGiqHjftZbvkL\n8gFFtmL73PC3x1MGZBWzpn6jzgNHUQgAOUaiyaSwm0Fa8Msk3bx26d6O9o1rKlG+mhydXZKPDiSa\nTAoLuf5333tHCW7duKpSFWrIwele3/sd2LJC4aE307xMcTFRuhjgryO71qtu8y5y+PsywcaEeN24\nckGLf/vigaMU4bdvqGjpinJ0cpGN7d35JyUl6s6tG1r06+e6dvH0fWHrhkpWqJ5iWV3z5tOBrXd3\n9oMDz2n9/F+0cdF0xcXG6KX3v5ejc54nOv+ipSuoRPmqOn/ioG4HX9We9QtlMiYob74CcnJ2UUJc\njM4e3as5P43T7RtXNWLCLyrhU+WhdSxQuJiiI+/o0umjSkpM1Cn/ndq4+Het8ftZezYsVoFCxTTy\nk+nKX7CwVs7+Kfl5dvaO8qlaVx2fHiF3D6/kx//pExUeql3rFirRZLq7fC6uigwP1b6NSzR/6iTl\nzVdAIyb88sARIOnu+To1G7eTjY2tLgQceuiRu8LFy2jExF9VslzVXPeZdLS31aI/vtdzw9/mCyqb\nMqz1v26mDMgq2tfylu/qsxQCQI4xomO5R/b553svPX3vFxMVrvlTJmnvxsXpXr6mnQbo6dc+zpRl\n9d+xRitn/59Cgi7Jw6uI6rXsptY9B8vR2SXTapUQF6v9W5bryK71Crp8RhFht2RjY6v8noVVulIt\n1WvRTZVqP5Xm+pnNZu3duFibl/6pG1cuyMHRSSXLV1fj9n1Uq0mH5H4jO5VXuWr11bh9X9Vq0sHi\neSkjOpbT65/PVL6ChXTiwHad8t+haxfPKCIsRHZ29ipSqrxqN+2opzr0T65ZakJvBmnH6rk6cWCb\nbgdfVXxctFzd8qu4TxXVbtJRdVt2lZ2dfa78TLq52Gtg81Ja63+dLygCDECAAQBrmf3DB4qOvKOX\nx05OtY8xIU4xkRG6ffOaLpw8pI2LfpMkfTZ7FwUEAQbZBifxAwCQzR3ctlI718xT9YatH9rP3sFJ\n7gW8VKZSLbXp/ZKGj5+q6Ig7FBAAAQYAAGSO6Mjw5JO0S1esla7nFilZPtULCQAAAQYAAGS4dfOm\nJB9Fye9ZOF3PvXL+pKrWb0ERARBgAABA5rj/krmhN9M+pt9sNmvdvKlq22coRQRAgAEAAJnDZDIm\nT0///HXdDr6apuf89eNY1W7W6aH3OQEAAgwAAMhQPQe/I4PBIEm6eiFAE15urz+/fVeHd61TyPVA\nJcTFKikpUXExUbp64ZQ2L5mhz1/rKZ+qdVW/ZTcKCCDbsaMEAABkX43a9ZabR0HN852okOuBMhkT\ntGf9Iu1Zv8hi/1IVaujZ/32i0hVrUjwABBgAAJD5qtRtrnG/rNXx/Zt1eOd6XT5zVHdu31BCfKzy\n5M0ndw9P+VStr+oNW6tircYUDAABBgAAWJetnZ1qNGqrGo3aUgwAORrnwAAAAAAgwAAAAAAAAQYA\nAABArsU5MAAAIMsyGKSCeRNVOH+iCuRNlLODWY72ZtnZmpVgNCg2waDoeBvdDLfV9TBbRcfx2yxA\ngAEAAMhkBfImqkrxBPl4G+Vob07z88KibHTyqoNOX3NQvNFAIQECDAAAwJOT3zVJTSvHqqiH6T8/\n/6mKcWpYPk5HLjnq0HlHGRMJMgABBgAAIAPZ2EgNysWpeql42WRA3rC1kWqXiVf5IkZtO+GsyyHs\n8gA55vuCEgAAAGtycjCra91o1SydMeHlfq5OSepYO1q1ysRTaCCH4OcIAABgNa5OSepeP1puLklP\nbB4Gg9SwfJzyOidp2wlnig5kcxyBAQAAVuFgZ1anOjFPNLzcr0rxBI7EAAQYAACA/6Z19VgVyJuY\nqfNsWD5OZQoZKT5AgAEAAEg7H2+jSnlZJ0g0qxKbrkszAyDAAACAXMzBzqynKsZabf7ODmY1LB/H\nGwEQYAAAAB6tcvEEuTha9whIxWLWXwYABBgAAJDFGQx3T6a3+g6QQapSnBP6AQIMAADAQxT1MGXa\nVccepVIxTuYHCDAAAAAPCzAFTFlmWfI4JSlfniTeFIAAAwAAYJl3/sQstTxFPEy8KQABBgAAwLKC\nblkrwBTMm8ibAhBgAAAAUrKzNcveNmtd+cuZK5EBBBgAAACLYcEh64UFJ25oCRBgAAAALLHNgnsd\ndrYEGIAAAwAAYEG80ZDllikuwcAbAxBgAAAALISFLBhgYhPYFQIIMAAAABaYzVJkbNba9YiK4wgM\nQIABAACw4OKpw9q4fmeWWqbrYXa8MQABBgAA4J6oiDDN+n6Mvn6zv1av3JxllivJLAWH2fIGAQQY\nAAAAKSkpUdtWzNZHg9to19r5MpuTtHz5NplMWePmkVdv28mYyBAygAADAAByvQsB/vpiVC/NnTxe\nsdERyY9fv35LixZtyhLLePyyI28UkA0x8BMAAGSYyDu3tXj6l9q7YbHM5pT3WLGzs9eeYzHq18+6\nyxkeY6PAW+wGAQQYAACQKyUlJmrbytla/ucPDxxxuV+FGo3Ub8Q4eZfw0bnrMfLxNlpteXcEOMvM\nPSwBAgwAAMh9zp84ID/fCbp64ZTF9nwFC6v30PdVp1mn5Md2nnJW8YImOdpnfoo4H2yvwBB2gQAC\nDAAAyFXCQ0O0ZPqX2rdpaarDxVr1fFEdnx4pR2eXB9pi4g3adMxZHWrFyJCJ59GHx9ho+0ln3jyA\nAAMAAHKLpMREbVk+Uytn/ZTqcLFKtZ9Sv1c+UqFipVN9nUs37bU9wFnNKsdmynLHJhi08kAexSZw\n5TGAAAMAAHKFs8f2yc93goIunbHY7uFVRL2HjlGtJu3T9HonAh1kazCrccW4J3okJjLWRqsPuSg8\nhguwAgQYAACQ44WHhmjRr5/pwJYVqQ4Xa9NniDr0f0UOTukbonX0sqPuxNiqbY0YOdhl/DkxwXds\nteYQR14AAgwAAMjxkhITtWnJH1o1+yfFxUZb7FO5bjP1f2WcPIuU/M/zCQyx05ztedW4YqzKZdDV\nyeKNBu0/56QTgQ5K4opjAAEGAADkbGeO7JGf7wRdDzxnsd3Dq6j6DBujmo3bZcj8YuIN2nDERSev\nmFS9ZIJKehll8x8OmsQmGHTqmoOOXHTkqAtAgAEAADndnds3tGja5zqwdYXFdnsHR7XpPUQdBgyX\nvYNThs8/KNROQaF2yuOYpHJFjPLOn6hC+Uxydkj9MEpYlI2C79jpyi07Xbxpr6Qk3keAAAMAAHI0\nk8mozUtmaNVfPyk+NsZin6r1W6jvsLGPNVwsraLjbXT4oqMOX/x7+aKva/onQ+Th4a7RX81QnNGg\nOKONYuINijdypAUgwAAAgFzj9OFd8vOdqOAr5y22FyxcXH2GfaDqDVtbbRmj4mx0+PDdq591u2KS\ngyP3cwEIMAAAIFcJvRmkRb9+rkPbV1tst3dwUru+Q9Wu38tPZLhYejg63bsZZqLJJDny/gEEGAAA\nkCuYTEZtXDRda+b6pjpcrFqDVuo7fKwKFi5OwQAQYAAAgHUEHNohP9+JunntosV2zyIl1XfYWFWt\n34JiASDAAAAA6wi9GaQFv0zS4Z3rLLbbOzipff/hattniOwdst74LFu7e7st8XExcs6TlzcVIMAA\nAICcxmQyasOCaVo952cZE+Is9qnRqK36Dv9AHl5Fs+x63H/SflJiIm8sQIABAAA5zYkDWzXv548V\nEnTZYrtnkZLq/8o4Va7bLFuszz/n49jasgsDEGAAAECOcfvGNS2YOklHdq+32O7g6KwOA15Rmz5D\nZGdnn23Wa+Lvm3hzAQIMAADIKYwJcVq/4Fet9Zua6nCxWk3aq/fQMfLwKkLBABBgAACAdRzbu1kL\nfpmU6nCxQsVKq98r41SpdhOKBYAAAwAArCPkeqAWTJ2kY3stD69ydHZRx6dHqlXPF7PVcDEAIMAA\nAJCDGBPitG7eL1o3f1qqw8XqNOukXkPeU35PbwoGgAADAACs4+juDZr/yyTdDr5qsd27hI/6vfKh\nKtRsTLEAEGAAAIB1hARd1vwpn+j4/i0W2x2dXdT52dfUovsghosBIMAAAADrMCbEac3cKVq/YJpM\nxgSLfeq26KLeQ96XewEvCgaAAAMAAKzj8K51WjB1kkJvBlls9y7ho/4jx6t89QY5vhYrZn6vfZuW\nqX6rbury3OtsHAABBgAAZBU3r13SvJ8n6uTB7RbbnZzzqPPAUWrZfZBsbG1zRU2iwsN0K/iKosLD\n2EAAAgwAAMgKEuJitXqurzYu/E0mkzFFu8FgUL0WXdVzyHty9/DMVbVp12+YGrfvK1d3DzYUgAAD\nAACszX/HGi2c9lmqw8WKlKqgASM/kk/VermyPh5eReThVYQNBSDAAAAAawq+ckHzfp6gU/67LLY7\n53FTl4Gj1LzrwFwzXAwACDAAAGQx8bExWjXn/7R58R+pDher36qHegweneuGiwEAAQYAgCzk4LaV\nWjjtc925FWyxvViZSuo/8iOVrVyHYgEAAQYAAOu4HnhPY98aAAAgAElEQVROfr4TdObIHovtznnc\n1PX519Ws8zMMFwMAAgwAANYRHxujlbN/1OalM5RoMqVoNxgMati2l3q8OFp58xWgYABAgAEAIPOZ\nzWYd3LpSC6d9qvDQEIt9ipetrAGvTlDpijUpGAAQYAAAsI7rl8/eHS52dK/FdhdXd3Ub9Iaadn5a\nBoMNBXuEI7vX62LAYZWuVFM1GrWlIAABBgAAZITY6EitnP2Ttiz7U0mJiSnaDQYbNWrXWz0Gj5ar\nW34KlkYBB3do28q/1KzzMwQYgAADAAAel9ls1v7Ny7T4ty9SHS5Wslw19R/5kUpVqEHB0qlSnSZy\ncnFV6UoMtQMIMAAA4LFcu3hKfr4TdO74AYvtrm751XXQm2rSsR/Dxf6jGo3acuQFIMAAAIDHERsd\noRUzf9DWFbNTHS72VId+6vbCmwwXAwACDAAA1mE2m7V342Itmf6VIsJuWexTqkJ19R85XiXLVaNg\nAECAAQDAOq6cPyk/3wm6cPKQxXZXt/zq/uLbaty+D8PFAIAAAwCAdcRGR2jZjO+0fdWcVIeLNe00\nQN1eeFMuru4UDAAIMAAAZD6z2azd6xdqyW9fKioizGKf0hVrasCrE1S8bGUKBgAEGAAArOPKuROa\nO3m8Lp46bLHd1d1DPQaPVqO2vWUwGCgYABBgAADIfNGR4Vo+41ttXzVXZnNSinYbW1s16/yMuj7/\nupzzuFGwTBB6M0hR4aFydfeQh1cRCgIQYAAAgNmcpF1r52vp79+kOlysbOU66jdiHMPFMtm6eVO1\nbeVfatb5GQ14dQIFAQgwAADkbpfPHJWf7wRdOn3UYrtbfk/1GPy2GrTuyXAxACDAAABgHVERYVr2\nx7fauWZeqsPFmncdqC4DRzFcDAAIMAAAWIfZnKQdq+dp2YxvFB1xx2Ifn6r11H/EOBUtXZGCAQAB\nBgAA67h46rD8fCco8Oxxi+3uHp7q+dK7qteyG8PFAIAAAwCAdURFhGnJ9K+0e93CVIeLtej2vLoM\nHCUnF1cKBgAEGAAAMp/ZnKTtK+do6R/fKjY6wmKf8tUbqP+Ij+RdshwFAwACDAAA1nEhwF9+k8fr\nyvmTFtvdPTzVe+gY1W3RhWIBAAEGAADriLxzW4unf6m9GxbLbDanaLe1s1PL7i+o87OvydHZhYJl\ncb2GvqduL7wpO3sHigEQYAAAyDmSEhO1beVsLf/zh1SHi1Wo0Uj9RoyTdwkfCpZNODg6y8HRmUIA\nBBgAAHKO8ycOyM93oq5eCLDYnq9gYfUe+r7qNOtEsQCAAAMAgHWEh4ZoyfQvtW/TUovDxezs7NWq\n54vq+PRIhosBAAEGAADrSEpM1JblM7Vy1k+pDherVPsp9R0+ToWLl6FgAECAAQDAOs4e2yc/34kK\nunTaYruHVxH1HjpGtZq0p1gAQIABAMA6wkNDtOjXz3Rgy4pUh4u16TNEHfq/IgcnTvoGAAIMAABW\nkJSYqE1L/tCq2T8pLjbaYp/KdZup3/AP5VW0FAUDAAIMAADWceboXvlNHq/rgecstnt4FVWfYWNU\ns3E7ipVD7VwzTwGHdqhS7SZ6qkM/CgIQYAAAyHrCb9/Uwmmf6cDWFRbb7R0c1ab3EHUYMFz2Dk4U\nLAe7cu6EDm1fLVe3/BQDIMAAAJC1mExGbV4yQ6v++knxsTEW+1St10J9h4+VZ5GSFCwXqNeyq4qV\nrcTNRwECDAAAWcvpw7vk5ztRwVfOW2wvULiY+g4bq+oNW1OsXKRslboqW6UuhQAIMAAAZA1hIde1\ncNpnOrR9tcV2ewcntes7VO36vcxwMQAgwAAAYB0mk1EbF03Xmrm+qQ4Xq9aglfoOH6uChYtTMAAg\nwAAAYB0Bh3bIz3eibl67aLHds0hJ9R02VlXrt6BYAECAAQDAOkJvBmnBL5N0eOc6i+32Dk5q33+4\n2vZ5ieFiAAACDADAOkwmozYsmKY1c6coIT7WYp8ajdqq7/AP5OFVlIIBAAgwAADrOHFgq+b9/LFC\ngi5bbPcsUlL9XxmnynWbUSwAAAEGAGAdt29c04Kpk3Rk93qL7Q6OzuowYLja9BkqOzt7CgaLrpw7\noZDrgfL0LqHiPlUoCECAAQAgYxkT4rR+wW9a6zdFxoQ4i31qNWmv3kPHyMOrCAXDQ+1cM0/bVv6l\nZp2f0YBXJ1AQgAADAEDGObZ3sxb8MinV4WKFipVWv1fGqVLtJhQLaVKoeFlVrNVYhYqXpRgAAQYA\ngIwRcj1QC6ZO0rG9myy2Ozq7qMOAEWrdazDDxZAuLbs/r5bdn6cQAAEGAIDHZ0yI07p5v2jd/Gmp\nDher06yTeg15T/k9vSkYAIAAAwCwjqN7Nmr+1E90O/iqxXbvEj7q98qHqlCzMcUCABBgAADWERJ0\nWfOnfKLj+7dYbHd0dlGnZ15Tyx6DGC4GACDAAACsw5gQpzVzp2j9gmkyGRMs9qnbvIt6D31f7gW8\nKBgAgAADALCOw7vWacHUSQq9GWSx3buEj/qPHK/y1RtQLAAAAQYAYB03r13SvCkf6+SBbRbbnZzz\nqPPAUWrZfZBsbG0pGACAAAMAyHwJcbFaPddXGxf+JpPJmKLdYDCoXouu6jnkPbl7eFIwAAABBgBg\nHf471mjhtM9SHS5WpFQF9R8xTuWq1adYeOJCbwZp7KDmkiTf1WcpCECAAQDgruArFzR/ykQFHNpp\nsd05j5u6DByl5l0HMlwMAECAAQBYR3xsjFbN+T9tXvxHqsPF6rfqoR6DRzNcDFaVaDLJ1o5dGYAA\nAwDItQ5uW6mF0z7XnVvBFtuLlamk/iPGqWyVuhQLVuHg5Jw8nRAfK2e7vBQFIMAAAHKb64HnNM93\nok4f2W2x3TmPm7o+/7qadX6G4WKwKltbdl0AAgwAINeKj43Rytk/avPSGUo0mVK0GwwGNWzbSz1e\nHK28+QpQMAAAAQYAkPnMZrMObl2phdM+VXhoiMU+xctWVv+R41WmUi0KBgAgwAAArOP65bPy852g\nM0f3Wmx3cXVXt0FvqGnnp2Uw2FAwAAABBgCQ+eJiorRi1o/asuxPJSUmpmg3GGzUqF1v9Rg8Wq5u\n+SkYsiRHJxeN+nRG8jQAAgwAIIcxm83av3mZFv/2RarDxUqUq6oBI8erVIUaFAxZmo2trSrWakwh\nAAIMACAnunbxlPx8J+jc8QMW2/O45VO3QW+pScd+DBcDABBgAADWERsdoRUzf9DWFbNTHS72VId+\n6vbCmwwXAwAQYAAA1mE2m7V342Itmf61IsIsDxcrVaG6+o8cr5LlqlEwAAABBgBgHVfOn5Sf7wRd\nOHnIYrurW351f/FtNW7fh+FiAAACDADAOmKjI7RsxnfavmpOqsPFmnYaoK6D3lSevO4UDABAgAEA\nZD6z2azd6xdqyfSvFBUearFP6Yo1NeDVCSpetjIFAwAQYAAA1nHl3AnNnTxeF08dttju6u6hHoNH\nq1Hb3jIYDBQMOUrAoR26cv6kipetrEq1m1AQgAADAMiqoiPDtXzGt9q+aq7M5qQU7Ta2tmra6Wl1\nG/SGnPO4UTDkSEd2rde2lX+pWednCDAAAQYAkBWZzUnatXaBlv7+taIiwiz2KVu5jvqNGMdwMeR4\nZavUVVJSkspWqUsxAAIMACCruXzmqPx8J+jS6aMW293ye6rH4LfVoHVPhoshV6jXsqvqtexKIQAC\nDAAgK4mKCNOyP77VzjXzUh0u1rzLs+ry3P8YLgYAIMAAAKzDbE7SjtXztGzGN4qOuGOxj0/Veuo/\nYpyKlq5IwQAABBgAgHVcOn1EcyePV+DZ4xbb3T081fOld1WvZTeGiwEACDAAAOuIigjTkulfafe6\nhakOF2vR7Xl1GThKTi6uFAwAQIABAGQ+szlJ21fO0dI/vlVsdITFPuWrN1D/ER/Ju2Q5CgYAIMBQ\nAgCwjgsB/vKbPF5Xzp+02O7u4aneQ8eoTvPODBcDAIAAAwDWEXnntpb8/pX2rF8ks9mcot3Wzk4t\nu7+gzs++JkdnFwoG/EtURJjiYqLk5OIqV7f8FAQgwAAAnoSkxERtWzlby//8IdXhYhVqNFK/EePk\nXcKHggGpWPHn99q28i816/yMBrw6gYIABBgAQEY7f+KA/Hwn6uqFAIvt+QoWVu+h76lOs84UCwAA\nAgwAWEd4aIiWTP9K+zYtsThczM7OXq16vqiOT49kuBiQRgNencCRF4AAAwDISEmJidqyfKZWzvop\n1eFilWo/pb7Dx6lw8TIUDAAAAgwAWMfZY/vk5ztRQZdOW2z38Cqi3kPfV60mHSgWAAAEGACwjvDQ\nEC3+9XPt37I81eFibfoMUYf+r8jByZmCAQBAgAGAzJeUmKhNS/7Qqtk/KS422mKfynWbqd/wD+VV\ntBQFAwCAAAMA1nHm6F75TR6v64HnLLZ7eBVRn2EfqGbjdhQLAAACDABYR/jtm1o47TMd2LrCYru9\ng6Pa9B6iDgOGy97BiYIBAECAAYDMZzIZtXnJDK366yfFx8ZY7FO1Xgv1HT5WnkVKUjAAAAgwAGAd\npw/v0ryfP051uFiBwsXUd9hYVW/YmmIBT9DaeVPlv2ONajXpoPb9hlEQgAADALhfWMh1LZz2mQ5t\nX22x3d7BSe36DlW7fi8zXAzIjM/kzSAFnj2uUuWrUwyAAAMA+IfJZNTGRdO1Zq5vqsPFqjVopT7D\nPpCndwkKBgAAAQYArCPg0A7N+3mibly9aLHds0hJ9R02VlXrt6BYAAAQYADAOkJvBmnBL5N0eOc6\ni+32Dk5q33+42vZ5ieFiAAAQYADAOkwmozYsmKY1c6coIT7WYp8ajdqqz7APVKBQUQoGAAABBgCs\n4+SBbfL7eaJCgi5bbPcsUlL9XxmnynWbUSwAAAgwAGAdt29c04Kpk3Rk93qL7Q6OzuowYLja9Bkq\nOzt7CgYAAAEGADKfMSFO6xf8prV+U2RMiLPYp+ZT7dTn5Q/k4VWEggFZTI3GbeVRqKiKl61MMQAC\nDADkbMf3bdH8qZ+kOlysULHS6vfKOFWq3YRiAVlUpdpN+IwCBBgAyNlCrgdqwdRJOrZ3k8V2R2cX\ndRgwQq17DWa4GAAABBgAsA5jQpzWzftF6+ZPS3W4WO2mHdV76PvK7+lNwQAAIMAAgHUc3bNR86d+\notvBVy22e5fwUb9XPlSFmo0pFgAABBgAsI6QoMuaP+UTHd+/xWK7o7OLOj3zmlr2GMRwMQAACDAA\nYB3GhDitmTtFGxb+KmNCvMU+dZt3Ue+h78u9gBcFAwCAAAMA1nF41zotmDpJoTeDLLZ7l/BR/5Hj\nVb56A4oFAAABBgCs4+a1S5o35WOdPLDNYruTcx51evY1terxgmxsbSkYkAMEX7mgO7eCla9gYRUu\nXoaCAAQYAMj6EuJitcbvZ21Y8KtMJmOKdoPBoHotuqrnkPfk7uFJwYAcZMvSGdq28i816/yMBrw6\ngYIABBgAyNr8d6zRwmmfpTpcrEipCuo/YpzKVatPsYAcyNU9vwoWLi5X9/wUAyDAAEDWFXzlguZP\nmaiAQzsttjvncVPnga+pRdfnGC4G5GBdnntdXZ57nUIABBgAyJriY2O0es5kbVr8e6rDxeq36qEe\ng0czXAwAAAIMYFlCQrzWLJqtTasX69K5AMlgUNHipdS0bVd17v288rrno0h4bAe3rdTCaZ/rzq1g\ni+3FylRS/xHjVLZKXYoFAAABBrDs3Knj+mT0UDk6OWnQiHdVs/5TMsigi+cCtODPn7V0znSN/vgH\n1W7YnGLhP7keeE7zfCfq9JHdFtud87ip6/P/U7POzzJcDAAAAgyQuqMHd2vsyGdUvkpNTfKdI0dH\np+S2yjXqadw39TTl63H6cNRzGvf1b2rQrC1FQ5rFx8Zo5eyftHnpH0o0mVK0GwwGNWzbSz1eHK28\n+QpQMAAACDBA6oKvBeqj/w2SvYOjPvz61wfCy/2GvjFOh/Zs06R3Xpbv3A0qVqosxcNDmc1mHdy6\nUgunfarw0BCLfYqXraz+I8erTKVaFAwAgFzKhhIgPb4Z/4ZioiPV67lhcs+f+q/ftrZ2enboG4qP\nj9NX40YpKSmR4iFV1y+f1Q/vPafpX7xhMbw453HTgJHj9d5PiwkvAADkchyBQZrt27FRRw/sksFg\nUPtuAx7Zv3HLjnLO46pTxw5p06pFatOlL0XEA+JiorRi1o/asuxPJSWmDLkGg40ateutHoNHy9WN\n+z0AAACOwCAdFs2cKkkqUaa8ChbyfmR/ewcH1WvcUpK0zO93CohkZrNZ+zYt1YSh7bRp8e8Ww0uJ\nclU1+rt5Gvj6p4QXAA9YMfN7jXuxlVbM/J5iALkQR2CQJrduXNfh/TskSRWr1U7z8ypWq61t65fr\n9HF/XTwboNLlKlHMXO7axVPy852oc8f3W2zP45ZP3Qa9pSYd+8lg4DcWAClFhYfpVvAVRYWHUQyA\nAANYdnD3FpnNZklS0RKl0/y8+wPLsYO7CTC5WGx0hFbM/EFbV8xOdbjYUx36qdsLb3LEBcBDtes3\nTI3b95WruwfFAAgw+LdP3xuurWuXptre85mhGj564gOPzZzytWZN/eaBx5q3764xn0/5T8vQvpZ3\nhq/XWv/r6ep//PC+5GnvYqXS/DzPwkWTpwOOHVS3AYPZqHIZs9msvRsXa8n0rxURZvnqYqUqVFf/\nER+pZPnqFAzAI3l4FZGHVxEKARBgYMmoMV+oas36+uW7CTImJCQ/XrdxS705/lsV8Cyc4jkDh72l\njj2f1fsjBuj61Ut6/cOv1bRt12xdh6uXzidPF0lHgMnv4Zk8fencaTaoXObK+ZOa5ztR508etNju\n6pZf3V98W43b92G4GAAAIMBkBFc3d3UbMFj2jo76fuLbyY83b9fNYniR7t5or2Ahb3kW8lbbLn0f\n++pb6T1a8iTcCLqSPO2WL+3Dexyc7t0nJiriDhtULhEbHaFlM77T9lVzUh0u1rTTAHUd9Kby5HWn\nYAAAIM34yTONOvZ8VrUaNEv++8Th/Q/tHx8Xq/OnT6hLv0E5Yv1joiOTp51c8qT5efZ29snTkZHh\nbEg5nNls1q51C/TRS221dfksi+GldMWaeu/HRRrw6gTCCwAAIMA8SYNGjE6e3rxmsaIiUt8h37lp\nteo0aiGXPHlzxLrHx8UlTzs7pz3AGE3Ge68RG8tGlINdOX9SX7/ZT7O+e19R4aEp2l3dPTTwjc/0\n9rfzVNynCgUDAAAEmCetUvW6Kl+l5t879LFavXh2qn03rpyv1p375Jh1t3dwSJ62tbNN8/NMxnvn\nDTk6O7MR5UDRkeGa+38f6fPXeuriqcMpv2RsbdW860BN+G29GrfrI4PBQNEAAAABJrP0fGZo8vQy\nv9+VlJRyiEzY7RBdvnBGtRo0yTHrff+RpPQcSYmPu9fXxcWVDSgHMZuTtHPNPE0Y0lbbVv4lszkp\nRZ8ylWvr3R8Wqf+Ij+Scx42iAQCAx8ZJ/OnUrF1XTftuokJv3dDN61e1c9NqNW3T5YE+W9YsUbO2\nXWVjY5sh88wKl1Eu4OmlsNs3JUlxcTFyzpO2MHIn9HbytJd3MTagHOLy2WPymzxel04ftdjult9T\nPQa/rQate3LEBUCGiwgL0abFf0iSegweTUGAXIYjMOlNfHb26tL3+eS/l/w1LUWfDSvnq3Xnvjlq\nvUv53LsBZciNtIef8LB7AaZYyTJsQNlcVESY/vrxQ335vz4Ww4uNra1adn9eH01bo4ZtehFeADwR\ncTFRWjf/F61fMI1iALlxf5wSpF/nPs9rzm8/yJiQoOP++3Qu4Jh8KlWTJAVePCuT0aiyFTLuJOWs\ncBnlcpWra8OK+ZKk4KuXVb5yjTQ9L+jKpeTpyjXqsfFkU2ZzknasnqdlM75RdCqXw/apWlf9R3yk\noqUrUjAAT3bnxd7x7+8mM8UAciGOwPwH+TwKqkX7Hsl/L77vKMzGFQvUqlOvHLfODZu1S56+cOZk\nmp8XeOFM8nStBk3ZeLKhS6eP6Iv/9dacnz60GF7cPTz1wuiv9caXfxFeAAAAASar6nHfyfxb1i5R\n2O0Qmc1mbVm7RK069s5x61u4aInkK7AdObArzc87feLuVakq16gn72Il2XCykaiIMM36foy+eqOf\nAs8eT/nlYWurVj1f1EfT1ql+q+4MFwMAAASYrMynYlVVrd1AkmQyGrVi/gwd998rL+9i8ixcJEeu\nc/8XX5UknTlx+IFzW1ITEx2p08f9JUm9Br7MRpNNmM1J2rZitsa/1Fa71s63eHWx8tUb6IPJy9Xn\n5TFy4upyADKZvYNj8nTsfTdaBpA7cA7MY+jx9BAdP7RXkrRywZ8KvhaoNjno3i//9lSrTvKpVE3n\nAo5p7dK56vfCyIf237pumRITTSpfuYaatO7MBpMNXAjwl9/k8bpy3vIwQXcPT/UeOkZ1mnfmiAsA\n6+282DtQBCAX4wjMY+3Qd5RX4aKS7t77Zfv65WratkuOXV+DwaA3xn0rO3t7LZ49TTEP+dUrISFe\nftP/T46OTnprwvfs7GZxkXdua+Z37+mbt/pbDC+2dnZq0/sljf91g+q26ML7CQAACDDZsng2tura\n/8Xkvxs0b/fADR9zIp+KVfXGuG8Udvumfvj4HYs38jSZjPp63P8UcuOa3v9iikr5cGJ3VpWUmKgt\ny2Zq/JB22r1uocUr+lSo0UhjJi9XryHvydHZhaIBsDpbWztVrNVY1Ru2lq0dg0mA3IZP/WPq1Gug\nZk39RvFxsWrdqXeuWOc2XfrKzt5e3098W+8N76+BL78pn0rVFBsTrSP7d2nOr9/LaDTqq2kLuXRy\nFnb+xAH5+U7U1QsBFtvzFSys3kPfU51mDP8DkLU4ODlr1KczKARAgMl+YqOj0nxH+CfF1c1dbbr0\n0Y4NK1X3qZa5ZsNp0b6HqtZqqEWzpuiHSe8o+FqgHByc5FOxqnoNHKY2XfrK3iFrjFGOi4nK8iea\nJ8THysHROVPmFR4aoiXTv9K+TUssHnGxs7NXy54vqNPTr3LEJYu8Z1l5PXJKHdhmWFcgu4qLjZFT\nLvt/na2HkG1ctTBLLEepshXVvH132dnZ56qNp6BXYb385nj9tniHVu4L1OIdZ/TVr4vUsdezWSa8\nhARd1tlj+7J8LQ9tX62E+NgnOo+kxERtXjJDE1/uoL0bF1sML5VqP6UxvivUc/A7hJeHiAoP1Yn9\nW7P9ekSEhejkwe2P9Rr7Ny9TUmIiG0Uaa2Xpqn45ze3gqzp/4hBvOJBJdm9dq+ioCAJMtgkwKxem\n6XK+T9qm1YvVOgdffSw7u3jqsC6eOpLll/P8iUMW77WSUc4e26dPX+2u+VM/UWx0yi85D68iGvrB\nT3pt0h8qXLwMG04atqtLp49k+/W4EOCvS4/5+Th9ZE+qwxDxoJMHdyg48HzO/3ycPqyLp/x5w4FM\ncvqYf/JtKwgwWZzJaNS5gKMKOHrQqstx/eolRYaHqmK12nyCsuiO5uUzWX9H8/KZI7p85liGv254\naIj++PItff/uQAVdOp2i3c7OXu37D9e4qWtUq0kHNpg0unT6aLYIxo8SePb4Y293l88c0+UnGL5z\n2vfRxVOHc8d6BhzmDQcyycmjBxRwLHcd9cy2Aeb86eNKSIjXKSu/YRtXLlSrXHLyfrbc0Tx1RJdO\nH83SwzbiY2MUdOmsLgRk3K8nSYmJ2rhouiYMaat9m5dZHC5WuW4zjZ2ySt1feEsOToxVT1+AOazA\ns8ez/dCpiwGHdenMEYtXE0yLqPBQhQRd5tf2NIgIC9GdW8EZ+jnPutuVvy6dOcKbDmSChIR4nT99\nXGdO5K4fDbLtSfynjt8NLtYMMOakJG1atUif+s7hE5QFGRPide3iKZlMRgUHnpd3yXJZc2f47x3I\njPpl9szRvfKbPF7XA89ZbPfwKqI+wz5Qzcbt2Ej+y+fenKRLp48qIT5WVy8EqES5qtlyPZKSEnXp\nzJHkAF2sTPovd/7PNnsxgADzKP8El5x+BMZkMurq+QCZTEaFBF2WZ5GSvPnAE3Qu4KhMRqPVRyRl\neoC5EXhO1y5mv/HLuzeukCQFHD2gg1uWWeXGeocP7lWRIkUUdO6Igs7xa1NGaNq0qa74r8uQ17oa\neEEmk1GSdHDdbNWo2zhLrvPhLWslSXduBevklnnK657vP71OVES4Nq5aqJOpfInZ2tmpYdM2atyi\ng+zslWF1zm1u3biu2L9v4uq/0U+GqKBsuR43r19TfGyMJOnQ+r9krt8k3a9xZMsySdKNa5d0Ztdi\nOTvnYQNJxdGtSyRJwVfO6/ze5XJwcMyR6xl09VLy9+6BtTNVtWZ93nxkOQ52NmratKkObV2e7ddl\n59YNkqTwsNtav/A3FSjolTsCTGToDV08eSDbLfj5M3fvFh4XFyv/XRuUP3/+JzavgwcP6vjx47Kx\nsZGPj4+qVq2q8PBw7du3T126dMmW9cuqKlWqpNuXMiYMnj16NHn6YsAhFSuYNXeuLp+5t5ynD21R\n6dKl0/X8pKQkHTt2TIcOHZLRaLTYp0SJEmrcuLHc3NwUfu0kG9pjOH3q1APbVclCbtlyPc7etx6X\nTvmrhFf6b8IbeO7E3QmzWWcOblGxYsXYQFL7nP99npA5KUmnD2xSkSJFcuR6nj1+73yoCycOyDvf\nkwtqO3bs0MmTJ1W5cmU1adKEjQzp3t/ICftv9x95Obhzo8qVK5cr3j87zzzZ7xezuLg4RUTcu5LS\njRs3nmiAOX78uOLj4yVJx44d07Fjx5Q/f3516NBBDlnkcsFI6ebNmxans/pypifAXLt2Tbt27VJY\nWJjF9rx586px48YqWZJhHBklJCQkW2xX6dnubty4ke7nm83mFNsuASb1Wt2/3dy4cSPHBpj7t6X/\nsl0B+O/f5SEhIbknwGT3L8h/3ryKFSs+sfm1bNlSe/bsUUREhNzc3FS+fHlVq1ZNtra2fHKyyYc6\nNDRURqNR9vZZ6149kZGRiomJSfcOcVRUlPbs2VhUcqQAACAASURBVKMLFy5Y/mDb2alGjRqqUaOG\n7Ozs2Bie0HZ1584dxcfHy9Ex+w0Huv979M6dO0pISEjXDzL/PIed1UcLDQ2VyWTKEcE3vd+7iYmJ\nT+x/ZcOGDVWvXj3+FyPXiomJUWRkZK74bkmxn1O4cHHVqdM0Wy10UFDIv3bmYp7oOtSpI/Xs+Syf\nlGwkIiJcUVG/JP9tNpuVL19hVaxYOUst5/79ex74+/bt26pVq7FsbCz/QzaZTNqwYbVWrVqm+Pg4\ni32qV6+lfv0GytPTiw0hg8XHx+nXX3994DE3N09VrVojW61HbGyM7tyZ9q/PRyFVqpT2CxLs3Ln1\nXzvpYapdu4lVzkfM6rZt25SiVtnt/25aREZGKjLy3vduUlKSChQoqrJly7ERAE+Av/+Bf323hKp6\n9YZZ7sfaJxJgfHyqyMenSrZa6Nmz/3zg7+vXr6lVqx7Kk8eVrRmSpNWrF6d4zMkpr7p2HZillvPA\ngQdPuE9ISFDp0lVVvXqdFH23bl2nsWNH6fz50xZfq1QpH02Y8J3atu3CBvCE7NixSUlJD16S29Ex\n621Xj96hXp/i0tr29nnStR5btz4YYKKjo1SlSn2VLVuBDeVfNm3a+K8d/QjVrNlExYuXylHruX59\nyhOinZ3ds93nA8guTpx48CJcJpNJJUpUUN0setGijJTt7gOTlJSkw4f3PfBYYmKi/P33sSUj2cGD\ne1I8dvjw/iy3nJa2238v57VrgRo6tI+efrq9xfDi5OSs0aMnatOmo4SXJ+zf3z2SdOjQnmy4Hvsf\n+/OREa+RWxw4sDtHbDf/5Xv34MHdbADAE2Lpe8TS55AAkwWcPRugiIjwXPHPAI8TDPZm+X+kRmOC\njlm4j9E/y240Juj/27vzqKautQ3gD/OMiAKCEyKIyBQGFbVWrVO9zkNxav0carW1tVrb2qvWq71W\nq71a517rgFPhqli1arVaRByqBZRRcURUBFERZYYwfH/QxIQkEDBAQp7fWqwVkp3knJN9znnfffbe\nZ926b/Hmm244fvyg3M94++0RiIi4jrlzv4Yxb0bZIAnn1at/yb1RqKYFmjU5hubl5eLmzUQeh5WQ\nnf0S9+7d1oogQ97vz8ZForqhqPFeWxqSNC6BkReY8iBJlXfqeDn3Qnny5DEePkxRm+VMTIyVO44l\nOvoSwsNPolcvd6xcuQgFBfkyZRwdnbF372/YseNQo+uGot6Bv2wS/OLFcyQn39L4QDMz8ylSUu4o\nnciVlpZqRVD+umJjI2W6HSqqS5qsrKxM7nn40aMHePLkMSsCkYrdvJmIvLxcpY7vTGDUIoGJrFFi\nQ9yp1XXHVhTA3L17ExMnDpIbTJqYmGL+/GWIiLiGt94axB+7Hj18mKIwENOkwD0l5Q4yM5++1noo\n2o+uX4+Tm3BrM3ndxwDg2rU4CIXFjWY9b99OUnjcZTcyItVTdKWlqnMVE5gGDfrknzjVrXWd1K+O\nVJUAN4S4uJrdQOsf/xiFc+eS8OmnC2FgwPsP1beqkt+4OM25ZF/VPqBs1wNF5UpKSuR2i2S9kVVc\nXITExFit2D94ZY6ofmMdbWg00JgEJi3tIY4dC5Xb71okKGgjkpLipebbJ+3x8mUWIiJO4dChYIVl\nzpw5gd9/P4KMjLQGW87c3BxcuhSBS5cilCrv5NQBISG/Y9u2g2jZsg1/6HpWUlKC69fjcPToAYVl\nzp37A5GRF5Cfn6fW63HtWmw163EaUVEXFV5FycnJxsWL4YiKuqjwM0JD9+DWretyu01pk+zsl7hw\n4UyVQcb+/Ttx584Njd5WQqEQiYkxOHYsVGGZs2d/x5Url1BYWKDS7w4O3oYZMwIRHLyNByrSGnl5\nubh8+RwuXTqrsMyRI/uQlBQvt6tvY6GTlqZ+o0+fP3+G2NgoxMVFIzY2CrGxkXj6VPmbpJmamsHD\nwwcCQWd4e/tDIOgMR0dn3p+gEcnPz0NiYgxiY6MQExOJuLhopfvvi7Ro0RICQWepetKkSVOVLmdx\ncRGuXYtDbGwU4uMr6vPt20lKByyGhobw8ekqtZyOjs6sAHWkvLwcKSl3pOpVYmKM0t2i9PT00KFD\nJwgEXcR1ys3Nq97n5C8vL0dy8i3xMTQuLhoJCVeVDiD19fXh7OyG1q0dYWJigry8XKSk3MW9e7eV\nrrvm5hbw9PSFj08XeHtX1N/GOl6rsLAA167Fird1TEwkkpNvKT25g4WFJby8/ODt7Q9v787w8emC\nVq3aqt16lpWVITn51t/n5SjExUUpHMunqF65unrA29sfPj5dIBB0hqurR61vthscvA0REafQq9cA\nTJjwPg9g1OgIhUJcuxYrPpbHx0fj1q3rSicmjTkebvAEJjc3BwkJVyUSlkg8eHBP5d9jaWn194nU\nX/xD2tu34t6hITtwUlK81EmzJjtwTTg6OkvVEw8PH5iamin13tLSUty5cwOxsZHi+lwX/dybNm0m\nVY8Fgs6ws3NgRamF9PRU8YlBlLBkZ79Q6XcYGhrB3d377yS0om45O3eErq6uStdDtPyiuqfq9VAF\na+vmEvW2IoC1sbHTqDpTUlKCW7euiY9HsbFRuHEjQeVX/ps3txXv315e/g2yrVJT70vUqYp6lZOT\nrdLvMDY2gbu7QFwvfHy6wMmpAxscSeuUlZX9HUNU7G+xsVG4di0OxcVFjIcbOoERCotx7VocYmIi\nxa3Rd+7caLBLXHZ29vDyetUS5O3tj6ZNm3EvUoMdWLK1ISEhRuU7sLJELeqS9UTUon7/fjLi4l4F\nvgkJVxUOYq37uuwgVY/r4mqSpsvKypS4qltxgsjISG+QZRFdnaj4vSpa3JW9OiFaj5iYSPE+0lDr\noQoODq2lglcvL39YWjZRi2UTXZETbeuaXpFTtVat2koEHZ3h5eWnsm2VmflUvG+Izs816fmg2gCr\nyd9XpF5deVbHK1JEr+PhwxRxL6PY2CgkJFxFbm5Og8XDov1NU+LhektgcnNzcOTI/xAd/WeDJy6K\nEhmBoDOGDx8LJ6cO3LMaKFiIirqI8+fDGvwEWtmrLh6d4ecXgLy8PJw/fxqxsVG4e/em2vRhNzEx\nhYeHD7y9/TFixHj4+nZlxULFYMfDh0P+btGKVXlf/NrS1dWFs3NH+Ph0QZ8+b2PQoFFVdjeLirqI\nX3/dj9jYKCQmxijddUcTiLo6+Pp2xahRE+Hh4dOgy1NYWIBff92Pc+dOIyYmEikpd9Tmfj9mZubi\nBFgV2+ry5XM4evQAYmMj66TFVxUNSP37D8HAgcOhp6fHAxppNKFQiBMnfkF4+Mkadyuva5JXRIcN\nC0Tnzj2YwMhLaCS7jsXE/FUvs4hZWlqJW7BEl9DYlUx9PXr0QKrFvC66MFRmZGQs02fUyamDwi4/\nOTnZiI+/Iu4CGRsbhdTU+3W+bQwMDODq6iFuMREIOqNDB/da9yfXFiUlJbh5M1GqVf3mzcR6mfyj\nTZt2UuNjPD19YW5uUeuT4M2bieL9IyYmErdvX9eISUwMDAzg7i6Q6iLVoUMntQ5Os7NfyHQ3TE9P\nrYdtZQh3d2+pbnfOzh3rbFuJuuxKXuGrqy67lTk6Oksdz2rShZdIUzVUPCw5Jk20z3Xs6KkxMYRa\nDeKXHrwf+dot8KLWaMm+fhzMr9kkB5GK6klNBpEqswO7unq89qDrZ8+eSE1CERcXjWfPntT683R1\ndeHk1EHi8m5neHgIYGRkzEqhAgUF+UhMjJHqYnbv3u3XanGXvLIrql/W1s3rZT0kx2e87pUDXV3d\n12od1NPTg7NzR6mxHO7u3jA0NNL4epORkY74+GipZDgrK1Ml20qU6Lq7ezf4tOn5+XniAEu0njWd\nNKWy+phEhUhTqToe1tHRgZNTB6kGfHd3AYyNTTR2G6nlLGSS0tIeSrQw/oWLF8MVlm3f3hU9evQR\n/0AuLp3YGq0FJFuiY2IiERZ2XOFNnCwsLNG//1DxSdPDwwcmJqb1spypqfelEpq//jqvsLW8WTMb\nBAS8CR+frvD29oeXlx8sLCz5Y9cjyRb32NgoXLoUgRcvnssta2BgiICAN6XGcqjLld3s7BdSCc3l\ny+eqXI9u3XpJBZb29q3E471EY28iIy8oTGpsbVsgIKCXSq4yaaKUlLviAbhxcdFVbis7OwepeuPp\n6QszM3ONWM8XL55LJTSXLkUonDjC2NhEaj058QhR7ePhiivAUfjrrwsKJwmytm4u3ucEgi4qHS/H\nBKaWBAJ7hcHp8eN/wcenC2u5llu/fjm++26h3NcCA/8Pa9fuVIvlHDWqFy5fPif3tdWrt2H8+Gn8\nMdXIli1rsHTpPLmv9e8/FLt2/aoR67F58yosWzZf7muDBo3E9u2/VPsZb7/tj/j4K3Jf+/HHEAwf\nPo4V5m/9+glw/Xqc3Ne2bg3F4MGjG8V6rl69FKtXL5H72ujR72LDhj2sDEQqNG7cAJw7d1rua99+\nuxFTpsxq1Ouvq2kL7OMjf1CygYEhPDwErNEEX98Aha95e3fWiOWs6jVqGH5+AbV6TZP2D2XXo6py\nrLs12VaNZ5KNqhoP2bBIpHr+/t0axTlJixIY+QdCDw9Bg/cTJvWpI4oGuKrTTq1oWSwsLOHi4sYf\nUs14evoq7JKqSUG7t7f/a6+HonK2ti0a7c0qa0sgkN9oYm/fqlFNIOPr21Xh+FJFDY9E9Dr7nPzj\nsJGRMdzcvJjAaMoP5ufXjbWZAFRMMers3FHmeRMTU3Tq5K1WJ3xFAaYqb3JIqmFkZCy3/ujq6mpU\nC7OJiSk6dHCXeV5PTw+enr6vFZSr0xVO9T9nNa4WUisra7m3IKiY1ZG9I4jqK4bw9PR97YmImMDU\nAUWt69pwuYxeL2jw8vJTq0kd7Owc5LbAsguO+pJ3nHFxcdOYgddVrYerq4fS6+Hk1AFWVtY8DivB\nxcVN7uBZgaDxdauSF1BVBFPsHUFUF40Gjo7OWnsc1rgExszMHB06dJJ5ni1/JEle31BFrcbqFkgy\nCNSsxFgTE055y1yTeqejoyM3WGXyLeckq6srN1lpjDeZ1Zb1JFLvGEI7eiRpZD+VygfJZs1s4OjY\nnjWZqgnQ1G+nlrdM7C+uvuR1FatqIKUmJfg1TT4ql9fT0+NgbYXnLOnGE319fXh7+2tFMMU6QVS/\nsY46NtYygVFw8mWrH1Xm7NxRpjuMOl7ZqHygad3aEc2b2/IHVFPt2rnIdJ3SxK5A7dq5yHRrqmmg\nWbm8vH2O5B97OnXyrrf7T9UneevF8zNR/R1bbG1boFWrtkxgNCXjZJcbqqxya7Ci8SYNrfK4HJ7s\n1VvlrlOKurSq/YG/UrcmS8smcie+qC6BkZx1inVX+XNWY7z6AlRcWZKcCIKz0hHVLTc3TxgZGUvE\nw9ozoZVGJjCVB83yEjVVFzSoaz/syjOjMRnXrHpV1ZTd6k7ySrZA0KXGM99VnnVKE7vS1ZfK3Zwb\nc7IneT5md1iiumVgYCjVaKBNMYRGJjCS05bq6upyAD8pCDS7asROLblsbMXWrABNk38vyeCytsmH\n5PsYrFZN8jzVmJM9yRZgNi4S1W8M4eXlrzXrrbE3mxAFDoqmqCSS7CKjzoGm6IBjYGAId3dv/nAa\nkMCIuk5pcoAmOf6qtl2aRHXXzMycN19VMshQdL+UxhhMsUGGqP7iYW2bSEVjExjRQZIHSFJE1P9a\nT09Prfuci1pjPTwEUn1ZST1JBqCafNWholuT82sdR0XHYU3uSld/56xufyd9fgrvWN8Y2Nu3gp2d\nPWelI6onosaomtzLqzHQ19wfrIv4xElUVcuEhUUTtZ7xR3RTQG0afKfp/P27obi4CLa2LTQ8qA4Q\nJzO1IZp1ig1J1fPwEMDQ0Egrxgr5+gbg3r3bnJWOqB60bu0IW9sWWjeGVmMTGBsbO7Ru7cibZFE1\nJ9Kuat/FUDSzFZNxzeHt3RmFhYWNYv8oLy+v/QlEXx9eXn6su0oQdRHVhmTPzy9AZrpxIqrLfa6b\n1tz/ReMTGADo2bMfXF09WHNJIR+fLrC0tFL75RQIOjMI1LDAv7i4SOPXQyDogtLSUhXUXTYkKRtk\naEOjm69vAJo2bcYfnKgeYx1tuxKuN2/ekiWauvCtW7fT+C4cVLeaNbNBkyZWat8aaGZmgY4dmYxr\nTr2yhampmcYff5o3t4G5uQVsbOxq/RlNmzZD+/YdWCmUYG3dXGo65ca8npaWTWBt3Zw/OlE9MDY2\ngZubZ6MeX1eZTlraa/QfICIiUmNFRYU4fvwgQkK2488/z+LRo7Jq37N8+T8RFLQRU6Z8jAULVnAj\nEhExgSEiIqpbiYkxCAnZjoMHf0Z29gvx82lp1Z/ynJ3NkZ+fBzMzc9y+ncONSUSkZnS5CYiIqLFZ\ntGg2kpNvQ1e35qe5adNmw9TUDFOnfsINSUSkhngFhoiIGq0LF8IQGNhP/L8yV2CIiEi98QoMERFp\nnLt3b8LBofoBq15eftxYRERMYIiIiBrW8eMHlSpnYmLGjUVExASGiIio4ZSXlyM0dI9SZfX19bnB\niIiYwBARETWc4OBtuHPnhlJltem+CERETGCIiIjUTELCVXzzzefcEERETGCIiIjU28mThxEY2A85\nOdncGERETGCIiEidCIVC7N+/C1OnjoS/fxs4OhrDxcUC/foJsGzZl3j06IFSn1NeXo4zZ05g9uxJ\n6N7dBe3amcDZ2RzdurXH9OljEBq6B0JhcY2WLT7+Cr744gP06tUJLi4WcHQ0hr9/GwQG9sOePVtQ\nWFhQ7WdcuxaLuXOnoGvXdnB0NIaHhw3Gjx+IoKBNSEt7iLCw36RmGRs7tj+mTh2Jly+zpD7HwUFH\n6i88/ORrbfeiokIcPhyCwMB+aNlSt9pyY8f2lykXHn4SU6eOhI9PSzg7m6N/fx/s2vUjymtw14Kc\nnGxs3boW48cPhI9PS7Rta4SOHa0wYIBvjX5/IqLGiPeBISJSMxcvhuPzz99H+/aumD17ATw8fFBU\nVIjTp4/hq68+RGFhAUxNzbBq1RaMGjVR4eckJcXjs8+mobS0FDNnzkOPHm+haVNrpKWlIizsODZu\n/A4ZGelwdGyPFSs2o1evAdUu28qVi7Bhwwp07OiJ5cs3wtPTFxkZ6QgN3Y3165ejpKQE7u4CHDx4\nFpaWTeR+xtata/HDD//Gp58uxODBo2FnZ4+srOe4fDkC69YtR1JSvLisvPu2SCY2ytzXRZnyiYkx\nCAnZgV9++VkqSapcvqpyT59m4IsvPsCpU7/K/Y6JE6fj++9/qnZ59+3biaVL56Fbt16YNWs+3N29\nkZWViaNHD2D16iXIzn4JAwNDfP31Krz//qfcYYiICQwRETWckJDtmD9/JsaMmYTVq7fJDEJfs+Yb\n/Oc//6o4gOvoICjoMAYMGCY3CZo8eRj69RuC9et3wcDAUKbMy5dZeO+9IYiO/hN6enr4z3+2YuzY\nKQqXbceODVi0aDYA4PLlZLRp007q9e++W4j165cDAGbOnIfFi/8j8xl//HEckyYNQXDwSfTuPVDm\ndaGwGJMnDxdfSamvBGbIkAAYG5sgLi4aeXm5CssPG9YDhoaGMuXOn7+B8eMHwte3K959dwa8vHyh\nq6uH48dDMX/+hyguLgIAhIaGo3v33gqXdfnyf2Ljxu8wb94SzJv3L5nXb926jsGDu4q/e9Gilfjo\noy+54xARExgiIqp/586dxsSJg9CyZRtERFyHkZGxTJm4uGgMGtRZ/L+rqzvCwxOlyjx8mIL+/QWw\nsGiC8PBEmJtbKPzO58+foW9fb2RkpEFXVxf79v2BHj36yC3bpYsjUlPvw9DQCCkphTKvJyUloG9f\nLwBAu3YuuHjxlkyZESN6IjLyAq5fz4SVlbXc78nMfIoePVyQnf2y3hIYkbNnf8eECW9XW/7cudMY\nN+7VFav27V3x7bcb8Oab/WXKrlr1NdauXQYAePvtEdix45Dcz9y/fxfmzJmM7t17IzQ0XOEybtq0\nEt9++xUAQE9PD+HhiXB27sgdiIi0BsfAEBGpAaGwGF988QFKS0sxZcrHcpMXAPD09EX79q7i/1NT\n78uUWbjwY2Rnv8TYsZOrTF4AwNq6Ob77bjMAoKysDHPmTFY4Jub582cAFN9bpW1bJ/HjJ0/S5ZZJ\nSLgKANi7V3FXqmbNbPDuuzMa5Hfw9vavVbmgoMNykxdR0iKSmBgjt0x+fp54drVp02ZX+d19+w4W\nPy4tLcX27eu5AxERExgiIqpfhw6F4OHDFABQGAgDgK6uLvbu/Q12dg6wtGyCzz9fKvV6fPwV/PHH\ncQBAv35DlPruAQOGiZOiR48e4JdffpZbbvr0OTAzM8fMmfPkvm5q+uqu94oG8puZmQOo6G62YsUC\nheXeeWdSg/wO1SV8r9bVXOr/qq6AtGvnIn789GmG3DIHD+4VJ4hdurxR5XdX7rp3/nwYdyAiYgJD\nRET16/Tpo+LHrVs7Vlm2bVsnxMQ8wo0bLzBjxmdSrx04sFv82MGhtVLfraOjgxEjxon//+03+V2c\n5s9fhtu3c2SSpszMp/jppx/Qr59A/Fxpaancz+jatSeAiqs9GzasQLdu7bF58/fIzn4pVc7V1V2p\n7mGqpq9voGQ5faU/U5S0ARCPhaksIuKU+LGnp63M7GqSf87O0slTenoqdyAiYgJDRET1Kz7+ivix\nou5jyvjzz1djJ5o2bab0+zp37iF3WaoSHf0nZs4ch549O+LOnRtYs2Z7te9ZsGAFLC2txP9nZKRj\n2bIv4e/fGosXz8H9+8kN+jtUnjThdcsBFVfNRBRNpSzZtezevQKkpZUr/Xf3bh53ICJiAkNERPXr\n2bMn4scvXjyv9eekpNwVPy4qKlT6fU5OHcSPRV2ZFLl0KQLDh7+Bd98dDC8vX1y+nIxVq7bAy8uv\n2u9p184Fx49fhq9vV6nnc3NzsG3bOrzxRgfMmjVR6+5zIrnNX+f3JyJiAkNERPVCsmX++vW4Wn+O\nZNetzMynSr9PckYwExMTuWUKCvIxd+5UjB7dG6amZoiIuI6PPvpS4f1eFGnf3hVHj17C5s3BcHV1\nl1n+Q4eC0aePh8L7qTRGRUWvupZduxbHHYKIiAkMEZF6k+zuFRb2m9LvS0yMwZAhAeL/mze3lXpN\nWZLd1iRnOXsVYBdi3LgB2LcvCAEBb2L37qOws7Ov9fpWjLsZjzNnErBz5xG5V2Tef380rly5rCW/\n/6sEUnQPHCIiYgJDpDKFhQUICtqEwMC+8PCwQdu2hvD0tMU777yFbdvWoaAgnxuJaqRjRw/x49DQ\nPTKD2hXZs2eLVPLj5/cqmbl48YzS35+XlyN+/MYbfWVe37BhBaKiLgIAvvrqW7k3xqxtIjNgwDAc\nO3YZO3cekZphq6SkBCtXLtKK319yFrN9+4LYjYyIsQIxgSFSndjYKPTq1QkLF36MCxfO4PnzZxAK\nhcjMfIqLF8OxePEcdO/ujL/+Os+NRUqTnDr5xYvnWLlyYbXvSU9Pxf79u/DWW/8QPzdkyDvix0eO\n7FM461VlosHzOjo6cqcwDg3dI37s6elb6/V0cNBBcvItua8NGDAMp0/Hwt+/u/i56Og/te73z83N\nwWefTVM44F/SqVO/YvTo3tyBiBgrMIEhIvmuX49DYGBf8f06FMnISEdgYF9cuHCGG42UUvmmk0FB\nm7BjxwaF5UtLS/HFFx9AV1cXo0dPFD//j3+MEncBy8rKxPbtG5T6flGiMGTIGLn3NJGcqjcjI12p\nz1Q0lfLBg3sVvsfCwhIrV/4o/t/YWHY8juSsXkKhsFH8/uPGTZVa15MnD+OTT96rciKGpKR4zJ07\nVWYcERExVmACQ0TiYOyTTyYhNzcHfn4B2LTpZ1y58hD37xcjLu4xNm7cKxX4CYVCzJgRyK4gpBQr\nK2uZ+6ssWjQbU6aMwIULYXj5Mgv5+Xm4fTsJu3f/FwMG+OLMmROYPXuB1LTEenp6WLs2CHp6egCA\ntWv/jbt3b1Zbt3/+eSusrZtjyZI1cstYW9uIH+/du0Xm9ZycbHzyyXtSzym6SeW2beuQmnpf4fK0\naeMkftyjRx+Z183NLcWPHzyoetrlygmAoqSqpuWFwmKp/8vKyhR+ZuUrKfKuitnatsDs2Quknvvl\nl5/x5ptu2LlzMx48uAehsBjPnz9DVNRFLF48B4MHB0BfXx9z5y7mDkTEWIEJDBHJ+v33I0hKises\nWfNx9OgljBw5Afb2rWBgYAAbGzuMGjURp05dRf/+r+5+npWVia1b13LjkVI++GAuxo6dIlPvAgP7\nwc3NGs7O5ujVqxO++upDJCXFo3fvgfjooy9lPsfPrxvWrdsFPT095ORk4733hiAl5Y7C4Hrp0nl4\n9OgBgoIOw96+ldxykvV6y5Y1+OGHf+Pp0ww8efIYQUEb0bOnK8rKymBoaCQud+lSBJKTb8kkNhXL\nNFjhDRj/+OMYAMDExBRffvlvmdclxwutXLkIL148R0ZGGr755nNMmzZKqmzlRCkjI63K36ByeUXL\nKDldNVD1jG9ZWZlS/6elPZRb7tNPF2LkyAlSzz18mIIFC2YhIMAJbdsawcPDBsOHv4Ft29bBwMAA\nO3Ycho2NHXceIsYKWkUnLa28nJuBSJng8h00b26L5cs3VVkuPz8Pfft6iccUuLl5ISyM06KScsrL\ny7F27TKsXbtMppVffODW0cGkSTOxZMmaKm96GRb2G7766kM8evQA5uYWmDHjMwwdGog2bdqhuLgI\nV65cxo8/fo+CgnysWbMDLi5uCj8rIyMdgwZ1xuPHj2Rea9bMBsuXb8LQoe/grbc8ceNGovi1Jk2a\nYteuX9GlyxsAKsbAiFhaWmHGjM8wZMgYtG3rhKys5zhx4hBWrPgnLCyaYMuW/VKTEkgGCFOmjJB5\n3s3NEyEhp2Br2wJCoRDp6an4+uvZOH36mLjMoEEj8fXX36Nly9ZSExEIhUI8fvwIixbNxunTR8XP\nDxw4HEuWrIGDQysYGBiipKQE9+/fxeLFjbjXwAAABtBJREFUc6RmCxs1aiIWLVoJW1t7cRe30tJS\npKenYvnyf+Lw4RCpZfjXv1bD3r6lzGQIot9/3bpvqxy/5OzcEf/97//QqZM3dxoixgpMYIhIvj59\nPHDo0Dmp+2UosmfPFsyfPxMAYGpqhjt3crkBqUbu30/G7t0/4syZE0hNvY/S0lI4OLRGz559MWnS\nTLi5eSn1Ofn5eTh8OAQnThzCjRuJePo0A3p6enBwaA1//24YOXKC1ADyqqSnp+K77xYiLOw35Obm\noFWrthg0aCRmzfpSvF+cOXMC8+a9j6KiQgwYMBTz5y+Tuqrj4KCDAwfOwN6+JcLDT+LcudNISkrA\nkyePYWhoCFdXDwwd+g4mTpwOMzNzhcty7FgofvjhGyQn30arVm0xcuQEzJjxmfg9komSImlp5VLL\npUz5mnxuTZdB0qNHD7BnzxaEh5/Egwf3kJ+fC2vr5vD09MXgwWMwatQElc0ER0SMFZjAEBEePLiH\ngICKfvxWVta4fj2TG4WIiIgYK6gAx8AQ1QF7+5bix61ateUGISIiIsYKKqLPTaA9qurOIK8bw/Tp\nY3D8+EGly7/O99dWbZajPrx8+UL8uGfPvqx8RETEGIMxBmMFFeEVGC1y40YWNm7cK9MvU9Ggsc2b\ng7Fu3S7x/y1btsH69buRkPCEG7MakjM+jRkziRuEiIgYYzDGYKygIrwCo0UsLa0watREODl1wIgR\nPcUz3EREnJY7INjAwBBDhozBnDmT0b17b2zdGqrUoDRF1PVqSV3444/jAIA+fd6Gm5snKx8RETHG\nYIzBWEFF9ObNW7KEm0G7tGjRElZW1ggL++3vg8sp6OsbICDgTZmyP/zwDZo3t8WOHYeqnBGIXikp\nKcGXX85Abm4OfvrpAGxsWnCjEBERYwzGGIwVVIRdyLTU5MkfYdiwseL/V636GidOHJIq8+uv+3D/\nfjJ+/PF/nK6zBg4e3IsHD+5h6tRPeI8GIiJijMEYg7GCinEaZS1WWFiAUaN6ITY2CgBgbGyCAwfC\n4OfXDWfP/o6QkO3YtCkY+vrsaaisnJxs9OnjASMjY5w+HQNTUzNuFCIiYozBGIOxAhMYUpXKd9e2\nsrLG4sXf4+TJI/jpp/0wNDTiRqqB+fNn4n//C8LRo3/Cy8uPG4SIiBhjMMZgrMAEhlQtIeEqRozo\niYKCfACAkZExoqMfoFkzG5V+j7pMcajsctT0s0+ePIxp00Zh1aotmDhxOisWERExxtCyGIOxQv3g\nGBiCp6cv1q3bBR2dip2/qKgQCxd+jHLmtkq7efMaPv30/zB58iwekIiIiBhjMFaoQ7wCQwCAhw9T\nMHZsf6k5yefNW4J58/7V6NZV1VdgMjLSMXRoN/j6BmDz5mDo6rJdgIiISBtjDMYK9YNbj5CRkY7P\nPpuKAwfCMHbsZPHza9YsxbFjoY1ufdPSypX6U8aLF88xYcJAuLi4YcOGPTwgERERaXGMwVihfvAK\njJbLysrEjBmBWLnyv2jXzgXFxUUYPboPrly5BAAwMTHF4cPn4enpy41VSXb2S4wb1x/GxiYIDj4J\nY2MTbhQiIiLGGIwV6hhTQC2Wk5ONDz8cj6VLf0C7di4AAENDI+zYcQj29q0AAAUF+Zg8eTiePHnM\nDVbpgDR+/ADo6Ohg9+5j1R6QysvLkZJylxuOiIgYY2hJjMFYgQkMqZhQWIzp08fgiy++gZubl9Rr\nNjZ22LnziHhHS09PxZQpw1FcXMQNh4oWpXfeeQulpaUICfkd5uYW1b5n5cpFmDt3CjceERExxtCC\nGIOxAhMYqgPr1n2LTp284OcXIPd1T09fbNmyX3x33JiYSKxYsUDrt9vjx48wcuSb0NfXx759p2Fp\naSVTpry8HMXFRXj6NANnz/6OyZOHY/365QgM/D9WPCIiYozRyGMMxgp1j2NgtFBcXDSGDeuOyMgU\n2Nk5VFn24MG9+OST9yoqi44Otm07iEGDRmrldktOvoVx4wYgNfV+jd9rbGyCuLjHsLCwZAUkIiLG\nGI00xmCsUD94BUbLJCRcxfvvj4ZQKERw8Ha8fJklt1xZWRmePHmM3NwcqdaCWbMmYteuH5GZ+VSr\nttuVK5cxbFiPWh2QAGDQoJE8IBEREWOMRhxjMFaoP7wCo2Xk3QNF3pTBCxd+jKCgTVV+Vl3coVZd\ntW9vJr6LcG0EB59E794DWQGJiIgxRiONMRgrMIEhIiIiIiKSwS5kRERERETEBIaIiIiIiIgJDBER\nERERMYEhIiIiIiJiAkNERERERMQEhoiIiIiImMAQERERERExgSEiIiIiImICQ0REREREWkJflykM\nERERERFpiP8HSKuc9RxV1VwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image('lh_schema.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the previous episode of this exercise, we worked out the posterior probability distribution $p(y|x, \\{D_k\\}$ of the lighthouse. Now, we work out $p(x, y | \\{D_k\\}$.\n",
    "\n",
    "From the information given, $\\theta_k$ for any $k$-th flash must lie between $-\\pi/2$ and $+\\pi/2$ radians, so the prior for the $k$-th flash is:\n",
    "$$ p(\\theta_k | x, y) = \\frac{1}{\\pi}$$\n",
    "\n",
    "And from elementary trigonometry, we can write:\n",
    "$$ y \\tan(\\theta_k) = D_k − x$$\n",
    "\n",
    "A the variable transformation leads to write the likelihood of a single $D_k$ as\n",
    "$$p( D_k | x,y) = \\frac{y/\\pi}{(D_k - x)^2 + y^2}$$\n",
    "\n",
    "This tells us that the probability that the $k$-th flash will be recorded at position $D_k$ , given the coordinates of the lighthouse $(x, y)$, is given by a Cauchy distribution in $D_k$.\n",
    "\n",
    "Assuming the flashes independent, thus the positions $\\{D_k\\}$ are also independent, the posterior of the position is\n",
    "\n",
    "$$p(x, y | \\{D_k\\}) = \\frac{1}{Z} p( \\{D_k\\} | x,y)p(x, y) = \\frac{1}{Z} \\prod_k \\frac{y/\\pi}{(D_k - x)^2 + y^2}$$\n",
    "\n",
    "with $Z$ is the posterior normalization constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Inferring the position of the lighthouse from the data involves the estimation of both $x$ and $y$. This exercise will show you that from a one dimensional dataset, and a correct model definition, one can infer two dimensional information about the lighthouse.\n",
    "\n",
    "We will use the tools you learned from the lectures so far, however many procedures exist and new methods will be learned in later chapters. In this exercise in particular, you should evaluate the necessary integrals just by summing on a regular, dense grid.\n",
    "\n",
    "1. Gridding the values: Grid the $(x,y)$ space allowed by your prior with a sensible step size and compute at each point the posterior value $p(x, y |\\{D_k\\})$. Be aware of the influence of your stepsize, make a sensible choice (you may need to adjust it in the following steps). Plot the 2D posterior PDF using sensible contours or a grey scale. Make sure you label you contours or shades of grey, otherwise they are meaningless! Make an initial estimate of the position $(x_{max}, y_{max})$ of the lighthouse based on the peak of your map (mode).\n",
    "    \n",
    "2. Marginal distributions: write down the expressions for both marginal distributions $p(x|\\{D_k\\})$, and $p(y|\\{D_k\\})$. Evaluate them numerically on the grid, and plot them. Write down the theoretical expressions for the means (expectation values), $E[x]$ and $E[y]$, and standard deviations, $\\sigma_x$ and $\\sigma_y$, for the two distributions, as well as their numerical approximations using the grid, and evaluate them (give the four numerical values).On top of your both marginal posterior distributions, plot a Gaussian of mean $E[x]$, $E[y]$ and standard deviation $\\sigma_x$, $\\sigma_y$, respectively. How well do these estimates agree with $(x_{max}, y_{max})$, given the width of the distribution? Explain potential source of differences.\n",
    "    \n",
    "3. Do a quadratic approximations of the 2d posterior distribution $p(x, y |\\{D_k\\})$, around the peak value $(x_{max}, y_{max})$. Work out (analytically) the necessary partial derivatives, showing all your steps. Using the covariance matrix of your approximation, plot the approximative 1-$\\sigma$ and $3-\\sigma$ contours on top of the posterior distribution. Compare your result with the ones from the previous question. Why might they differ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise 3 (45 points)\n",
    "\n",
    "In this exercise, we will study a coin flipping problem, where instead of one coin, we will model the use of two.\n",
    "\n",
    "Throughout this exercise, we will assume that a single coin flipping is  a random event that takes only 2 values: head or tail. In addition, we note $p$ the probability of getting a head on a single toss.\n",
    "\n",
    "1) What statistical distribution could represent the single coin flipping process? Generate 1,000 flips of a fair coin  $(p=0.5)$ and for a biased coin ($p=0.2$). Produce an histogram for each case that illustrates the coin tossing behavior. \n",
    "\n",
    "2) Write down an analytic expression of the likelihood of a dataset of $N$ draws resulting of tossing a coin with an unknown $p$. Deduce the posterior probability distribution of $p$, where $p$ is the probability of getting head on one draw and detail your prior hypothesis.\n",
    "\n",
    "You receive a dataset `coinflip.dat` and are told that this dataset contains some results of coin flipping in which $0$ means $tail$ and $1$ means $head$. However, you are given no information about the number of coins nor their probabilities $p_1, p_2, ..., p_k$. The following questions will guide you to whether you can find out the number of coins and their properties.\n",
    "\n",
    "3) First, let's assume a single coin was used. Grid the values of $p$ and assuming a uniform prior, compute the posterior of $p$ on that grid. (make a sensible stepsize choice). Plot your result. Give an estimate of the value of $p$ and its uncertainty, and justify your choice of estimators.\n",
    "\n",
    "4) Now, let's build and test a 2 coin model. For this model, we need to assume the coins with unknown parameters ($p_1$, $p_2$) and that we do not know the ratio of number of tosses $r = n_1/n_2$ in the dataset. Without more information, we assume uniform prior on $r$. However, we can assume the join prior distribution of $p_1, p_2$ as:\n",
    "\\begin{eqnarray}\n",
    "P(p_1, p_2) \\propto & \\\\\n",
    "&1 & if\\ 0 \\leq p_1 < p_2 \\leq 1, \\\\\n",
    "& 0 & otherwise\n",
    "\\end{eqnarray}\n",
    "    \n",
    "    4.1) Explain why this joint prior distribution is more suitable than a uniform prior on each of $p_1$ and $p_2$. Justify why we do not need a prior on the absolute number of tosses per coin. \n",
    "    \n",
    "    4.2) Write down the likelihood of the data with this model. \n",
    "    \n",
    "    4.3) Through gridding, compute the posterior distribution of the model, and plot the 3 joint distributions (each pair) and 3 marginalized ones. Describe what you observe and in particular give the most probable combination of ($p_1, p_2$) and the behavior of the model in the limit of $r$ goes to $0$ or $1$. \n",
    "\n",
    "5) Based on your previous analyses, can you conclude whether one or two coin model is better? Explain.\n",
    "\n",
    "6) You are now told that 2 coins were used, and that $p_2 = 0.8 \\pm 0.01$. Which prior on $p_2$ should be considered given this piece of information and what does that mean for $P(p_1, p_2)$?  With this new prior, plot the new posterior distribution and discuss what you can infer about the two coins properties: $p_1$, $p_2$, $r$.\n",
    "\n",
    "7) Let's suppose that instead of the information on $p_2$ you have learned that the tossing ratio $r$ is almost $1/3$. Therefore, with $N$ is the total number of tosses, we adopt a new prior on $r$: \n",
    "\\begin{equation}\n",
    "P(N \\times r) \\propto Poisson( \\lambda = N/3)\n",
    "\\end{equation}\n",
    "Justify this prior and with this precise information on $r$ and without the previous prior information on $p_2$, plot the new posterior distribution and tell whether you can tell about the two coins properties. Describe the differences with the previous case (question 6). \n",
    "\n",
    "8) Assuming you have now both $P(N\\,r)$ and $P(p_2)$ from question 6 and 7, what value of $p_1$ and uncertainties do you get? Explain your steps.\n",
    "\n",
    "9) Based on your analyses, do you think one piece of information is more important than another? Do you think you could disentangle if 3 coins were used? (no calculation is required, only write down your thoughts on this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}