{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a routine to return the Gaussian PDF value given and input, but add an option to return the natural log as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(x,mu,sig,log=False) :\n",
    "    if log :\n",
    "        return np.log(1./sig/np.pi/2) - (x-mu)**2/2/sig**2\n",
    "    else :\n",
    "        return 1/sig/np.pi/2 *np.exp(-(x-mu)**2/2/sig**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write a routine to return the likelihood of a full data set, where each point is drawn from a Gaussian with known mean and standard deviation. Include an option to return the log(likelihood) instead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "def likelihood(sample,mu,sig,log=False) :\n",
    "    p=1\n",
    "    sum=0\n",
    "    for i,s in enumerate(sample) :\n",
    "        #print(i)\n",
    "        if log :\n",
    "            sum += np.log(gauss(s,mu,sig))\n",
    "        else :\n",
    "            p *= gauss(s,mu,sig)\n",
    "        \n",
    "    if log : return sum \n",
    "    else : return p    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a random sample from a normal distribution with some mean and standard devation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu=10\n",
    "sig=2\n",
    "n=100\n",
    "sample=np.random.normal(mu, sig, size=n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now numerically try a range of \"guesses\" for the mean, calculate the likelihood for each, and determine the maximum likelihood value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/holtz/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \"\"\"\n",
      "/Users/holtz/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  \"\"\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.5\n",
      "1.0\n",
      "1.5\n",
      "2.0\n",
      "2.5\n",
      "3.0\n",
      "3.5\n",
      "4.0\n",
      "4.5\n",
      "max:  10.559999999999967\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3yb1b348c/RsCRLlvdOnMRk7+E4hBVKKA1QSCDstkDLLaUt7eXy66CXlku53LZ00VJoKS2rQMsIUEIJo+yVOHsnTpzheMeOh2xZW+f3hyTHcTzkkcQk3/fr5Zek5zmPnqOAvzr+Puf5HqW1RgghxMnLcKI7IIQQ4tiSQC+EECc5CfRCCHGSk0AvhBAnOQn0QghxkjOd6A50lZGRoUePHn2iuyGEEJ8p69ata9BaZ3a3b9gF+tGjR7N27doT3Q0hhPhMUUqV97RPUjdCCHGSk0AvhBAnOQn0QghxkpNAL4QQJ7m4Ar1SapFSqlQpVaaUuqOb/Ral1HPR/SVKqdHR7Wal1JNKqS1KqR1KqR8NbfeFEEL0pc9Ar5QyAg8BFwKTgWuVUpO7NLsJaNJajwXuB+6Lbr8SsGitpwFzgG/EvgSEEEIcH/GM6IuBMq31Xq21H3gWWNylzWLgyejzZcBCpZQCNGBXSpkAG+AHXEPScyGEEHGJJ9DnAxWdXldGt3XbRmsdBFqAdCJB3w3UAAeAX2utG7ueQCl1s1JqrVJqbX19fb8/hBDHVeVaWP8UtNae6J4IEZd4bphS3WzrWsS+pzbFQAjIA1KBj5RSb2ut9x7RUOtHgEcAioqKpEC+GN5eux1qNoEywNfegpFzT3SPhOhVPCP6SmBkp9cjgOqe2kTTNMlAI3Ad8IbWOqC1Pgh8AhQNttNCnDABD9RtgymXgw7DgU9PdI+E6FM8gX4NME4pNUYplQBcAyzv0mY5cEP0+RXAuzqydNUB4DwVYQdOB3YOTdeFOAFqt0A4CFOXgnNE5LUQw1yfgT6ac78VeBPYATyvtd6mlLpHKXVptNmjQLpSqgy4HYhNwXwIcABbiXxhPK613jzEn0GI46dqXeQxfw7kTJNALz4T4ipqprVeAazosu2uTs+9RKZSdj2urbvtQnxmVa2HpFxw5kYC/e43I+kcs+1E90yIHsmdsUL0R9W6yGgeIoFeh+Hg9hPbJyH6IIFeiHh5mqBxD+TPjrzOmRZ5lPSNGOYk0AsRr1hAz50ZeUwZBRanBHox7EmgFyJeLZWRx9TRkUeDAbKnQu3WE9YlIeIhgV6IeLmqIo/OvMPb0guhuceFfYQYFiTQCxEvVw3YUo+cYZOUC211EA6duH4J0QcJ9ELEq7UGnF3KPCXlRGbeuKVGkxi+JNALES9XVWQE31nsdWvN8e+PEHGSQC9EvFw1R+bnITKij+0TYpiSQC9EPIJ+cB/sJtBHX8uIXgxjEuiFiEdbtPZ810Bvz4yUK5ba9GIYk0AvRDxc0crcSV0CvdEE9iwZ0YthTQK9EPGIBfquI3qI5OllRC+GMQn0QsSjI9DnHr0vKVdG9GJYk0AvRDxaa8CcCNaUo/c5JdCL4U0CvRDxiM2hV90sj5yUC+2HIOg7/v0SIg5xBXql1CKlVKlSqkwpdUc3+y1Kqeei+0uUUqOj27+klNrY6SeslJo5tB9BiOOguzn0MbG59G11x68/QvRDn4FeKWUksiTghcBk4Fql1OQuzW4CmrTWY4H7gfsAtNbPaK1naq1nAl8B9mutNw7lBxDiuGhvAEdW9/tid8fKTVNimIpnRF8MlGmt92qt/cCzwOIubRYDT0afLwMWKnXU37jXAv8YTGeFOGHcDZCY3v2+2Ihe8vRimIon0OcDFZ1eV0a3ddsmuph4C9D1t+Jqegj0SqmblVJrlVJr6+ulOJQYZkJB8Db3HOjtmZHH9obj1ych+iGeQN/N1Sd0f9oopeYB7Vrrbldo0Fo/orUu0loXZWZmxtElIY4jT1PksadAb0uLPLY3Hp/+CNFP8QT6SmBkp9cjgOqe2iilTEAy0Pn/+muQtI34rGo/FHlMTOt+vykBEpIk0IthK55AvwYYp5Qao5RKIBK0l3dpsxy4Ifr8CuBdrbUGUEoZgCuJ5PaF+OzpCPQ9jOgh8iUQayfEMGPqq4HWOqiUuhV4EzACj2mttyml7gHWaq2XA48CTymlyoiM5K/p9BbnAJVa671D330hjoN4A71HRvRieOoz0ANorVcAK7psu6vTcy+RUXt3x74PnD7wLgpxgsUV6NNlRC+GLbkzVoi+xAK4rYccPUQDvYzoxfAkgV6IvrQ3QoIDzNae29jSJNCLYUsCvRB9aT/U84ybmMR08LdGVqISYpiRQC9EX9oP9Z6fB0hMjTzKBVkxDEmgF6IvcQX69MNthRhmJNAL0Zd4Ar3cHSuGMQn0QvSlvVFG9OIzTQK9EL0J+iIXWeO5GAuSoxfDkgR6IXoTS8X0OaKPpW5kRC+GHwn0QvQmnrtiAUyWyFz79qZj3ych+kkCvRC9iTfQQ/SmKRnRi+FHAr0QvYmn/EGMFDYTw5QEeiF6422JPNpS+m4rpYrFMCWBXojeeJsjj9Z4Ar1UsBTDkwR6IXrjaQaDGcy2vtvaUg8vOyjEMCKBXojeeJsjaRvV3bLIXVhTwOuCcPjY90uIfogr0CulFimlSpVSZUqpO7rZb1FKPRfdX6KUGt1p33Sl1Eql1Dal1BalVC+1XoUYZjzN8aVtIDKiR4Ov5Zh2SYj+6jPQK6WMwEPAhcBk4Fql1OQuzW4CmrTWY4H7gfuix5qAp4FbtNZTgHOBwJD1XohjzdsS34VYONzO03zs+iPEAMQzoi8GyrTWe7XWfiKLfC/u0mYx8GT0+TJgoVJKARcAm7XWmwC01oe01qGh6boQx4G3HyP6WDvJ04thJp5Anw9UdHpdGd3WbRutdRBoAdKB8YBWSr2plFqvlPpBdydQSt2slFqrlFpbX1/f388gxLHjaQZrcnxtYyN6r4zoxfAST6Dv7iqUjrONCTgL+FL08TKl1MKjGmr9iNa6SGtdlJmZGUeXhDhOYhdj42GLLT4igV4ML/EE+kpgZKfXI4DqntpE8/LJQGN0+wda6watdTuwApg92E4LcVxoHcnR9zd1IyN6MczEE+jXAOOUUmOUUgnANcDyLm2WAzdEn18BvKu11sCbwHSlVGL0C2ABsH1oui7EMeZrBR0ewMVYydGL4cXUVwOtdVApdSuRoG0EHtNab1NK3QOs1VovBx4FnlJKlREZyV8TPbZJKfVbIl8WGlihtX7tGH0WIYZWx12xcebozTYwWiR1I4adPgM9gNZ6BZG0S+dtd3V67gWu7OHYp4lMsRTis8XTj/IHMbZUSd2IYUfujBWiJ7GAHW/qJtZWUjdimJFAL0RPYpUr+zOit6ZI6kYMOxLoheiJZ4AjekndiGFGAr0QPenvxViIVrCUWjdieJFAL0RPPM2gDJCQFP8xVsnRi+FHAr0QPfG2REbzhn78mthSwN8KoeCx65cQ/SSBXoie9KegWUzH3bGSvhHDhwR6IXrSn4JmMbF6N3JBVgwjEuiF6El/CprFSBkEMQxJoBeiJ/1ZXSqmoya9jOjF8CGBXoie9Gd1qRipSS+GIQn0QnRH64FdjO2oSS+pGzF8SKAXojsBD4T8/b8YK6kbMQxJoBeiOwMpaAZgSgBzoqRuxLAigV6I7gykoFmMFDYTw0xcgV4ptUgpVaqUKlNK3dHNfotS6rno/hKl1Ojo9tFKKY9SamP05+Gh7b4Qx8hACprFSE16Mcz0ufCIUsoIPAR8nsgasGuUUsu11p2XBLwJaNJaj1VKXQPcB1wd3bdHaz1ziPstxLE1kIJmMVKTXgwz8Yzoi4EyrfVerbUfeBZY3KXNYuDJ6PNlwEKllBq6bgpxnA1kdakYSd2IYSaeQJ8PVHR6XRnd1m0brXUQaAHSo/vGKKU2KKU+UEqdPcj+CnF8xHL0semS/SE16cUwE8+asd2NzHWcbWqAAq31IaXUHOCfSqkpWmvXEQcrdTNwM0BBQUEcXRLiGIsFaouz/8faUmVEL4aVeEb0lcDITq9HANU9tVFKmYBkoFFr7dNaHwLQWq8D9gDju55Aa/2I1rpIa12UmZnZ/08hxFDzNEfq0BvjGQt1YU2BgBuC/qHvlxADEE+gXwOMU0qNUUolANcAy7u0WQ7cEH1+BfCu1lorpTKjF3NRShUC44C9Q9N1IY6hgRQ0i5EyCGKY6XO4orUOKqVuBd4EjMBjWuttSql7gLVa6+XAo8BTSqkyoJHIlwHAOcA9SqkgEAJu0Vo3HosPIsSQGkhBs5jOd8c6soauT0IMUFx/l2qtVwArumy7q9NzL3BlN8e9CLw4yD4KcfwNpKBZjNSkF8OM3BkrRHe8A1h0JEZq0othRgK9EN0ZqtSNEMOABHohujOoi7GSuhHDiwR6IboK+iHQPogRfTTlIyN6MUxIoBeiq467YgcY6I2myBx8ydGLYUICvRBdDaagWYyUQRDDiAR6IboaTEGzGJsUNhPDhwR6Iboa6OpSnVllRC+GDwn0QnQ1mNWlYqQmvRhGJNAL0VUsQA92RC+pGzFMSKAXoqshuRgrywmK4UMCvRBdeZrBZAOTZeDvYUuBoBcCnqHrlxADJIFeiK4GU9AsRsogiGFEAr0QXXkHUecmRmrSi2FEAr0QXXkGUbkyJlbvRkb0YhiQQC9EV4MpaBbTkbqRKZbixIsr0CulFimlSpVSZUqpO7rZb1FKPRfdX6KUGt1lf4FSqk0p9b2h6bYQx5CnRVI34qTSZ6CPrvn6EHAhMBm4Vik1uUuzm4AmrfVY4H7gvi777wdeH3x3hTgO5GKsOMnEM6IvBsq01nu11n7gWWBxlzaLgSejz5cBC5VSCkAptYTIguDbhqbLQhxD4RD4hmBEb00GVJ8jeq314M4jRBziCfT5QEWn15XRbd220VoHgRYgXSllB34I/LS3EyilblZKrVVKra2vr4+370IMvY7yB4O8GGswgtXZY46+rKmM29+/ndlPz+bCFy/koY0PEQqHBndOIXoQT6BX3WzrOgzpqc1Pgfu11m29nUBr/YjWukhrXZSZmRlHl4Q4RoaioFlMD2UQttRv4boV17GyeiWXj72cMcljeHjTw3zz7W/iCcoNVmLomeJoUwmM7PR6BFDdQ5tKpZQJSAYagXnAFUqpXwIpQFgp5dVaPzjongtxLAxFQbOYbmrSV7ZWcuu7t5JmTeOpC58iMzEysHlx14vcvfJufr3m1/xk/k8Gf24hOokn0K8BximlxgBVwDXAdV3aLAduAFYCVwDv6kjy8exYA6XU3UCbBHkxrHmGcERvSz1iRK+15t5V9+IP+Xly0ZMdQR5g6fil7Hft54ltT3BW/ll8ruBzgz+/EFF9pm6iOfdbgTeBHcDzWuttSql7lFKXRps9SiQnXwbcDhw1BVOIz4ShKGgWYz2yVPFHVR/xSfUnfHPGNxmdPPqo5t+Z9R3Gp47nF6t/QSAUGPz5hYiKZ0SP1noFsKLLtrs6PfcCV/bxHncPoH9CHF9DsbpUTKfUTSgc4ldrfsWY5DFcO+nabpsnGBO4bfZtfOudb/HPPf/kyvG9/koJETe5M1aIzga7MHhnsYuxWvN+xfvsd+3n1pm3YjaYezzkrPyzmJ4xnb9s/ouM6sWQkUAvRGfeZjCYwZw4+PeypUI4AIF2ntrxFHn2PM4rOK/XQ5RS3DLjFmrcNbyx/43B90EIJNALcaRYQTPV3Yzhfor+VbC9Zg3r6tZx3aTrMBn6zpaemX8mI5NG8uLuFwffByGQQC/EkYaioFlMNM//fNnL2Ew2Lht3WVyHGZSBpeOWsq5uHXtb9g5NX8QpTQK9EJ15hqAWfYwtBZ+Ct2pWsrBgIc4EZ9yHLhm7BJPBxLJdy4amL+KUJoFeiM6GoqBZjC2VD202WkMeLim8pF+HptvS+dzIz7Fi7wopjSAGTQK9EJ0NxepSMdYUXnXYyTDZmZc7r9+HLxq9iEPeQ6ytWzs0/RGnLAn0QnQ2FKtLRbmMRj5KtHGRoxCjwdjv488ecTY2k4039785JP0Rpy4J9ELEhMOREX1i2pC83QcH1xNUigvMWQM63mayce6Ic3m7/G2C4eCQ9EmcmiTQCxHjawEdBtvQBPp3K94jMxRmWnjgv2ZfGPMFmnxNrK5dPSR9EqcmCfRCxLQ3Rh5jC3sPgjfo5ZPqTzgvaMQwiOUEz8w7E6vRyvsV7w+6T+LUJYFeiJhYAbIhSN18Wv0pnqCH81TSoJYTtJqsnJ57Oh9UfCCrUYkBk0AvREws0A9B6ub9ivdJMicx15I56AXCF4xcQLW7mt3NuwfdL3FqkkAvRMwQpW601nxc9THz8+ZjTkwb9ALhC0YsAOCDig8G9T7i1CWBXogYTzTQDzJ1U9pUSr2nnrPyzzqqJv1AZCZmMiV9Ch9USqAXAyOBXogYTxOgBj2P/uOqj4FIyeFITfoWGGR+/az8s9jSsAWX3zWo9xGnprgCvVJqkVKqVClVppQ6avUopZRFKfVcdH+JUmp0dHuxUmpj9GeTUiq+qk5CnAjtjZEgP4Cbmzr7qPIjJqVNiiwVaEsFHQJf66Dec37efMI6zJqaNYN6H3Fq6jPQK6WMwEPAhcBk4Fql1OQuzW4CmrTWY4H7gfui27cCRVrrmcAi4M/RxcOFGH48jYNO27j8LjbVb4qM5uFwOYVBpm+mZ04n0ZTIp9WfDup9xKkpnhF9MVCmtd6rtfYDzwKLu7RZDDwZfb4MWKiUUlrr9uiaswBWQOaHieHL0zToC7Fra9cS0iHm582PbIgVSBvkzBuzwczcnLmsrFk5qPcRp6Z4An0+UNHpdWV0W7dtooG9BUgHUErNU0ptA7YAt3QK/B2UUjcrpdYqpdbW19f3/1MIMRTaGwc9tXJ17WqsRiszMmdENnSM6AcX6CGSvqloraCytXLQ7yVOLfEE+u6W2uk6Mu+xjda6RGs9BZgL/EgpZT2qodaPaK2LtNZFmZmZcXRJiGPA0zTo1E1JTQmzsmaRYEyIbIj9hTDIET3A/NzIXwkyqhf9FU+grwRGdno9AqjuqU00B58MNHZuoLXeAbiBqQPtrBDH1CBTNw2eBsqayyjOLT680TY0OXqAMcljyErMYmW1BHrRP/EE+jXAOKXUGKVUAnANsLxLm+XADdHnVwDvaq119BgTgFJqFDAB2D8kPRdiKIUC4HMNKnWzpjYyI+b03NMPbxzC1I1SijPyzqCkpkQWIxH90megj+bUbwXeBHYAz2uttyml7lFKXRpt9iiQrpQqA24HYlMwzwI2KaU2Ai8D39JaNwz1hxBi0GKBeBCpm5KaEpLMSUxMm3h4Y4IdDKYhSd1AJH3j8rvY0bhjSN5PnBrimuqotV4BrOiy7a5Oz73Ald0c9xTw1CD7KMSx5xl8+YOSmhLm5MzBZOj0a6VU5D2HYEQPdKxUtbJ6JVMzJAsq4iN3xgoBg65zU9VWRWVb5ZFpm5ghKIMQk25LZ2LaRLkgK/pFAr0QMOgSxatrIguDFOcUH73TljJkqRuIpG82HtyIJ+gZsvcUJze5S1UIGFTqxl9Zie/hJ/j5ZoV++AZ2er2YMjKwTBhP0vnn4zQ6MXgODVlXi3OLeXzb42w4uIEz8s4YsvcVJy8J9EJAp9RN/CP6wMGDHPzlr3C99hpTlebQqFSSLrgAQ2Iiwbo62jdsoO3tdzjoMJMxK0zqf4RQxsHV0QGYnTUbkzKxuma1BHoRFwn0QgC468FoAUtSXM1db71FzX/fifb7MX5lKd9wvsx3L7idBeOv6Gijw2Ha166l4ae3U/fRIVpv/Cp5v7wPc27uoLqaaE5kWuY0WUdWxE1y9EIAuBvAnhmZJdMLrTX1DzxA1Xf/k4TCQgqXv8K6yyfTlKSYlzPviLbKYMBeXEzB/7uQ3OJmvNu2se+qq/Bs3Tbo7s7Lnce2Q9to9Q+uKqY4NUigFwIiI3p7Rq9NdDhM3b3/R8Mf/0Ty0ssZ9fRTJIwezeqa1eTZ8xiRNKLb45Q9k5TCdkY//ieU2Uz59dfTvm7doLpbnFNMWIdZVze49xGnBgn0QkA00PdeZ+ngr35N0zPPkPbVr5J7770YEhII6zCra1dTnFuM6umvgegXiCXHwehnn8WcnU3F12/Gs2nTgLs7I3MGFqOFkpqSAb+HOHVIoBcCDqduenDoscdpfPxxUr/8ZbJ+8P2OoL6zcScuv6vjRqZuObKi5ziIOSuLgicex5ieTsUt38RfUdHzcb1IMCYwK2sWJbUS6EXfJNALoXWvqZvW99/n4K9+RdKiRWT/94+OGLn3On8+JvYF4o5U/zBnZzPykT+jw2EqvnELoba2AXV7Xu48djft5tAQTt0UJycJ9EL4WiHk63ZE7y8vp/p738cyaSJ5P/8ZynDkr0xJbUlHVckexd637WDHJsuYMYx44AH8+/dTe/dP0QNYUzZ28XdNnSwvKHongV4Id3Sxmy6BXgcCVH3v+2A0MvIPf8Bgsx2xPxAKsK5u3VGzbY5iTYkUNnMfuaiOfV4xmd+5Fde//kXLSy/3u9uT0ifhMDskTy/6JIFeiGhKpWugb/jTn/Bu2ULuT+/GnN91UTXYemgrnqCn9/w8gMEQeW/30aunpd98M4nz5lF777349uzpV7dNBhNF2UUd6SMheiKBXoiOEf3hHH37+vU0PPxnkpcswbloUbeHldSUoFDMzZnb9znsGd0GemU0kvfLX2Kw2aj6r9sJ+3z96npxbjEHWg9Q01bTr+PEqUUCvRBdUjehtjaqf/BDzHl5ZP/4zh4PK6kpYWLaRJItyX2fo4cRPYA5O4u8n/8M365dNPzxT/3qeuwisNwlK3ojgV6IjtRNZER/8De/IVBdTd4v78PocHR7iCfoYVP9pr7TNjH2LGjreeF7x4IFJC9ZwqFHH8W7c2fcXR+XOo5US6oEetGruAK9UmqRUqpUKVWmlLqjm/0WpdRz0f0lSqnR0e2fV0qtU0ptiT6eN7TdF2IIuOvBkgwmC+0bNtD87HOkfvlLJM6e3eMhGw5uIBAO9D6tsrNY6qaX2TVZP/wBRqeTmh//BB0MxvW2BmWgOLeYVTWrBjRzR5wa+gz0Sikj8BBwITAZuFYpNblLs5uAJq31WOB+4L7o9gbgEq31NCJryspqU2L4ic6h14EAtf9zN6bsbDK/+5+9HrKqZhUmg4k52XPiO4cjC4Ie8Pc8Z96UmkrOj+/Eu3UrjU89HXf3i3OKOdh+kHJXedzHiFNLPNUri4EyrfVeAKXUs8BiYHunNouBu6PPlwEPKqWU1npDpzbbAKtSyqK17t8VJyGOpWj5g8Ynn8S3axcjHnoQo8Pe6yGrqlcxM3MmiebE+M7RcdNUfa8VMpMuvBDHq/+i/ve/J+nz55Mwovv6OZ3F0kera1czyjmKtiYfzXXttDV58Xsii4gbzQbsyQmkZCeSnJWIwdB78TZxcokn0OcDne/TrgS6JiY72mitg0qpFiCdyIg+Zimwobsgr5S6GbgZoKCgIO7OCzEk3A34DSOo//NDOM5fSNLChb02b/I2saNxB9+Z9Z34z2GP3lDVVg9phT02U0qRc9dP2HPxFzl4332M+MMf+nzrfOsIZrUuYO9LPp5o/oT2Fn+v7U0WI3ljkxkzPYOxRdlY7eb4P4f4TIon0Hf31d81GdhrG6XUFCLpnAu6O4HW+hHgEYCioiJJNIrjy11PXYkFDAZyfvzjPpvHblDqdn3YnsSmbvYw86Yzc24uGd/4BvX330/bJ5/gOPPMbtu11HvY9lEVOz6tYV7b5XjNbvJnppBTmEJ6np2kdCsJNhNKQcAXxt3so6nWTd1+F5U7m/jgH7v4+IUyxszMYOo5+eSPH/jC6GJ4iyfQVwIjO70eAVT30KZSKWUCkoFGAKXUCOBl4Hqtdf/uCBHiWAsFaNvbStu2g2TefjvmnJw+D1lVs4okcxJT0qfEf56OwmZ9B3qAtK/eSPOLL1L3s59j/+fLKPPhUbe72cfqV/ey49MaUIoxMzJoPe0A91beyQuXvsCEtKPTPZZEcKRayB7jZOL8yMIn9Qda2bGyhl2raylbe5ARE1OZt7iQnDFxTBcVnynxBPo1wDil1BigCrgGuK5Lm+VELrauBK4A3tVaa6VUCvAa8COt9SdD120hhoZurqZuvRNzVgppN1zfd3utWVm9kuLcYoyGfiwLmBj/iB7AkJBA9o/uoPKb36LxmWdIv/FGgv4Q694oZ+O/DxAOa6Z/biSzLijAnmKh1p2JXqYpqSlhQtqEuM6RWZBEZkESZ1x2Gls/rGL9m+W8eN86Cmdlcs7V47GnWOL/fGJY6zPQR3PutwJvAkbgMa31NqXUPcBarfVy4FHgKaVUGZGR/DXRw28FxgI/UUr9JLrtAq31QYQYBpr+/jR+l5kRty7FYOk7sFW0VlDtruZrU7/WvxOZEiLr0bbWxn2I49xzsZ9zNg0PPoRn+uf44JUamuvaGTc3m3mXFpKcebj2To49h1HOUayuXc31U/r+wjqya0Zmnl/A5LPy2PROBeveKOfvO0s4c+lYJp2Z23OdffGZEdeasVrrFcCKLtvu6vTcC1zZzXH3AvcOso9CHBPBpibqn1iGPduH4/zPx3XMqppVAJye14/8fIwzH1xds549U0qR8f07eOe/HuXAI2U40qxcettMRk7sfgHz4pxiVuxbQSAcwGzo/wXWBKuJuRePYdzcbN57aifvPb2T3WvrOP+rk7Eny+j+s0zujBWnrPoHHiDs8ZE9uwWV3Pc0RoCV1SvJtedSkDSA2WHOPHBVxd28tdHLipcaOTBiIXnVn7DkssQegzzAGXln4A642XhwY//71klKViJL/msWC66bQO2eFp7/vzVUlTYN6j3FiSWBXpySvKWlND/3PKlnFWJJM0JizwE0JhQOUVJbwvy8+QNLZzjz4h7RV+xo5PmfraGxxs0F149jatPbNP7y5+hwuMdjTs89HZMy8VHVR/3vWxfKoJh6Tj5X3FFEgs3EK7/bwNrX98vdt59REujFKUdrTd3Pfo7R6STzrDRIyoE4Avf2Q9tp9bf2b1plZ858aG+AgLfXvm16p4JXH5F1rHYAACAASURBVNiILSmBK+8oYtwZI8n63v/Du2ULLa8s7/FYR4KD2dmz+bjq44H1rxvp+Q6u/FERY4uyKXllL/9+bDvBQGjI3l8cH3Hl6IU4mbS+9W/aS0rI+Z+7MHqfioy04/BB5QcYlIFZSVOpKSvF3dREW1MjPncb4XCIcCiEUgasDgeWRDv25BRScvJwZmZhMBoPn6e1BtLGHPX+4VCYj18oY8v7lYyZkcH5X51MgjXyK+q85BIa//53Dv72NyR9/vM93rl7dv7Z/Gbdb6h115Jj73uqaDwSrCY+/7XJpOfbWfXPvbQe8nDhLdNJdCYMyfuLY08Ntz/FioqK9Nq1a090N8RJKuz1svfiL2Kw2xnz0ouoh4ogfzZc8Vi37YN+PzW7d3Jg22ZeX/k8jmaNub1/vzMGo4nU3Dxyc5zkVb1E/pd+TWrRJUekf/zeIG89uo3yLYeYef5I5l8+9qgyBZ7Nm9l/1dWkf/3rZP2/27s9157mPSx5ZQl3zb+LK8cfNT9i0PasP8jbj2/H5kzg0u/OJCU7zhIQ4phTSq3TWhd1t09G9OKU0vj44wSqqih44gmU0RjJmU+8+Ig27a4WytasYnfJJ1Ru30ow4EcphXb4sJ82ltNnnE9qbj5JaenYU1KxOpIwmIwoZUCHw/ja3XjdbbibGmmqraa5ppr6A/sp276Nre3j4dePkJL9KoVzijltzjwyx0xgxUNbqdvXwoJrxzN1QfcXhm3Tp5O8ZAmNTzxByhVLSRg16qg2hcmF5Nnz+Kjyo2MS6E+bnYUjzcq/HtzES79exyXfmUlmQc+1e8TwICN6ccoI1NSw58KLcJxzDiMe+D20N8Ivx8AXfkZo7s3sWVvClnffonzLRnQ4THJ2DqfNLqZg2gzWmHbzfxvv45XFr1CY0nOtmt5oj4vGe8ZTUfAl9jY7OLBtM6FAAKPZiTJOZMGXLmXmBTN7/wwHD7J30YUkzp/PyIce7LbNvavuZfme5Xx8zcckGI9NeqWp1s3yBzbibw9y8benkzdOyiecaDKiFwI4+KtfgdZk/eAHkQ2ualoDCWxcU8XWx79Ke0szSemZzL10KeNPP4us0YUd6ZXfvP0kIxwjGJN8dG49XsrmJD05gfRRBmZ+86e01Lt48b5luOo2EvKt5Z1HV7O7ZDqzFl1K4Zy5GLq589aclUX6LbdQ/9vf9lgH5+z8s3mu9DnW1a1jft78Afe3N6k5di7/3hxefWAjyx/YxBe+PpUx0zP6PlCcEBLoxSnBvXo1rhWvk3HrrSSMyKfhwH7W/v1JdpTNRe/ZQuGcecw4fxGjZsw6KsC2B9opqSnhqglXDf4u0ehNU21NPl57aDsBXyFL/3sJabmKre+/zaa3VvDKr+8lOSub4sVXMuXchRhNR978lHbjDTQvW0bdz3+O/Z//RJmO/DWemzOXBEMCH1V9dMwCPUBSmpXLvjebf/1hE68/vIWFN0xiwryhuQAshpakbsRJTweD7Ft6BeHWVlIefYRVryyjdOVHmMxGpjkOMOd7j5I8tufVpN478B7ffe+7/OWCvwx8amXM00txNQV5pe6HeFoDfPHWGeSNS+nYHQ6FKFu7irXLX6KmrJSkjEzmLbmKqZ87/4iA3/rOO1R++1ayf/xj0r78paNOc8u/b6HaXc3yJT1Pxxwqfm+QFX/aQtWuJs69bgJTzs4/5ucUR5PUjTip6LAm1OQlcLCdUIuPUIs/8tjsI+wJogNhtD9EOBCGYBgdCmIs/DraYuLgg1so0IWMnzGTVEcjloPvYWzMxre/BXOOHYP16F+JDyo/wGF2MCcrztWketFsHMcrpUUEEoJcetvMoypFGoxGxs87k3HFZ1C+aT2fLvs7b//1IUpefp7TL7+aqZ/7PAajEcd552E/Yz71f/gDzosvwpR6ZI787BFn84vVv6DcVc4o59EXbYdSgtXEF789nTce2cr7z5QS9IeZsXBk3weK40ZG9GJY02FN8GA7vv0tBKrcBGrdBOrcaH+nO0QNYEyyYExOwGA3oxKMKLMBZTYQ9LVT/sYbuK1WjMYE0nNHkpqeB+1hwg1NhENH1nAxJlsw5yRiyrFjKXBiGuXgghWLmJ01m9+c+5tBfZbGajev/PIjwj4Pl/5wAZmj+74bV2tN+eYNfPrCM9TsLiV9RAELvnITY2bOwbd7N3uXXEbq1VeRc9ddRxxX01bDBS9ewG2zb+OmaTcNqt/xCgXD/PvRbezZUM+8Swspumj0cTmviOhtRC+BXgwrWmsCNW58Zc349rXg2+9CeyILZRsSTZhz7B0/puxETKlWDA4zqsucc601u1Z9zDsP/hZPwM/kufM566ZbSErrdMHwLwvRJiehS/5OoMFDoNZNsK498mVysB1Ckd+NAwk1JBamMX7mdKzjUzEk9r9gWH1FK8t/vxFD2Multv8i/fZXIP20fv27lK1eyYfPPE5zXQ2jps9iwVduIvS3Z2j6+98Z/cLz2KYcWR//2n9dC8A/vviPfvd3oMKhMO/8bQe7SuqYs2gU8xYXSvXL40RSN2JYC/tD+PY0493ZiHdnI6HoUnimDBu2KelYRidjGePEmGaNK2g0Vlfx7uMPU755A852HwumFzHl+92sHNW0DzXpEkwZtsi5OhUM04Ew/qpW3vnkNXz7Wpi7P5/GnaWgIKHAiXViKtYJaZhz7X32qW6fi1f/sBGzxcjiq+2k/LMCDu3pV6BXSjFu3hkUzpnLxjdfY9WLz/LUD77LtAULyU1Po/bunzL62X9E7g2IWjhqIb9f/3tq2mrIdeTGfa7BMBgNnH/DZEwJRta9UU7AH+KsK8dJsD/BJNCLE0IHwnhLG2nfXI93RyM6EEYlGLGMS8F5fhrWCakYnb2XxtVa4/f78Xg8eL1e2lpdbP3wPXaVfIrRbCbXkUmqduO/ZDFbt27FaDRiNBoxm83YDEGs7X5szjEkaH1UIFJmAwmjnDyw7m+MPn00S867BX9la+TLqLQJ15vluN4sx5RuxTY9E9v0TMw5iUe9T3VZM/96cBM2h5nFt83CaW2L7Ggc2GJrRpOZORcvYfKChXz6/DNsemsFuwpzGL9rH85nnyX9S4cvzJ5fcD6/X/973j7wNl+Z/JUBnW8glEFx7nUTMJuNbHq3gqA/zILrJsiC5CeQpG7EcaPDGl9ZM+0bD+LZdgjtC2Gwm7FNy4iM3Mcko0yH6+yFw2FcLhcNDQ0cOnSIlpYWWlpacLlcuFwuWltbCfdSzTFeBoMBh8NBcnIyTqcTp9NJcnIyreZWfrzhx9yx4A4uG3fZEceEWv14dzTSvqUe355mCIMp00bijEwS52RjSrVSsbORFX/cjCPVyuLbZuFItYDW8IsCmH41XPzrQfe9bt8e3v7rH6ktKyXN4+fC/72PnBmzOvYvXb4Uq8nKMxc9M+hz9ZfWmpLle1n3ejnji7NZeMMkDEapo3isDDpHr5RaBPyeyApTf9Va/6LLfgvwN2AOcAi4Wmu9XymVDiwD5gJPaK1v7etcEuhPPsFGL+51dbSvqyPU7ENZTdimppM4IxNLYQrKqPB6vdTW1lJTU0NNTQ11dXUcOnSIYDDY8T5Go/GIYJyUlIRBh9m/bjX1e0pJTk1j3qWXk5mSTs23b8U2dQr5996L1ppQKNTxEwgE8Ox6D8/KR/GefjvtRidtbW1HfIl0Pq/JZCI9PZ2srCxyc3M7fqxWKwChNj+erYfwbK7Ht68FgHC2nY37W/GkWbjkttlHFgD784JIWeSvvDwk/746HGbDC//g4+efJmgyMueSy5l/xbUkWG08uuVRfrf+d6y4bAUjncd3JkxYa9pCYT7+dzmr3q8gZ1IqMy8ZQ1CBLxzGrzWB8OH4oxQYUCjAoMBiMGAzGEg0GrAZDdgMCpvRQLLJSIJBvjC6GlSOXillBB4CPk9kEfA1SqnlWuvtnZrdBDRprccqpa4B7gOuBrzAT4Cp0R9xitCBMJ7tDbjX1EVGvIBlbArJF43BMjGVuoaDlJbvpHJjJTU1NTQ2NnYcm5SURHZ2NmPGjCEjI4P09HQyMjJwOBwdqZFQMMj6Fa+wctk/0FrzuaXXUPTFJRhNZiq+fSvO1lYK//tOEtJ6mNlS8SKwFc5bBAlHVoLUWuNqdXHjSzdSaCrkgswLaGhooLy8nC1btnS0S0tLIzc3l/z8fAoKCsidOwXtClD52l78mxuYbTWgwmF87x3AVJRDQr4jcmD6aVB59GBG6xChUDvBYBuhkJtgyE047AcdRusQoNE6hCaMQmEwWDAYrBgMFiZ8cT4pB3ez7q1/s/ZfL7Lz0w8576vf4KJJF/G79b/jtX2vccuMWwb+HzTKHw5T5Q1Q5fNT7w9y0B/gYPSx3hekIRCkJRjCFf3RAFZgUTIQhs0DS1l1lWg0kGoykmY2kWI2kmo2kWk2kWMxk2Mxkxt9zEkw4zD1Y23fk1Q8OfpioExrvRdAKfUssBjoHOgXA3dHny8DHlRKKa21G/hYKTV26LoshrNgsxf3qhrcq2sJtwcxplhI/FweTXkh9jbWUr5pIxX/qsDvj1xwTU5OJi8vj5kzZ5Kbm0tOTg5JSb0XyarcvpW3H/0jhyoPcFrR6Zx34804M7MAcK1YQds775D1/e+RMKKXG3ca94Ej56ggD5ELn9vatlGmyrjlrFv4wugvdOxra2vr+KujpqaGyspKtm3bBoDZbCbNmU1buYmstBwuuXgm4R1NtGzcRcP2j9Ej3BjHhtCpTfj9jfg33Ig/2Ijff4hg0EUo1N7ff+4jTYG0KZBGLeFAGbsrV2GqdfLD7ASaax5ht70FqzUXqyUHiyUHizWPBHMaSh0eHWutaQgE2e32ccDr44DXT4XXT4XHzwGvnxpfgK45gASlyEwwkZUQCbAT7VaSTUac0Z9kkxGHycjBHY3seK+SnLwkFlx2GnabGbOKjODDgEajNWgifw34whpPOEx7KIwnHMYTijx3BUM0BUI0BYORx0CQaq+H+kAAV/DoVF6yycgoawKjbBZG2xIYbbMwKvqYZzFjOAUuFMcT6POBik6vK4F5PbWJLibeAqQDDfF0Qil1M3AzQEHBAJZoEyeU1hrfnhbaVlbj3X4IjaZ9rJm6TD8HWsrZv3p/RyokKyuLGTNmUFBQwKhRo3A6nXGfp72lmQ+efoztH76LMzObJT/4CafNOfy/YqDuILU/vQfrjOmk3XBD72/WuBfSei5O9krZKzgTnJw78twjtjscDsaNG8e4ceM6trlcLvbt28SOzR/R2voR2TPqsVrb+KTCjTXDjSErcPgNPIA2Yk5xYnFXY0kagcM+EbM5BaPRjtFkx2S0YzTaMZkcGAwWUAYUBpQyglIojGhChMN+wmFf5CfkIxz2463aQ8Nzj2GcPAr3CCeNtWU4LUHsDj8HKh4DHfnvEEbRQCbVajT1pinUGk+jUmdzIJiMK3x4BKyAPIuZkdYEzkx1MNKaQIE1gRHWBLISzGQlmEg2GeObVZOVQqk1kXee3EFF826++J0ZWGxDOx/EHQpR5wtS4/NT6wtQ4wtQ6QtQ7vGxta2d1xuaCXb6pkpQisJECxPs1o6fiXYro20WjCfRF0A8/8rdfdquX+rxtOmR1voR4BGI5OjjPU6cWGFfiPYNB2lbWU17nYsqWzM1I9qoaK+jtaINKiAjI4M5c+YwZswYCgoKSEzsf/3ycDjElnfe5KN/PEnA62PeZVcx77KrMFusHW201tT85MeEfT7yfv6Lo+q/HKVxL4xd2O2uVn8r7xx4hyVjl2AxHjnzJxRqp7VtB62t22lr3Y67vQy3u4xg0EVmPmQCBoMdozEbny+D5iYjTU0GvF4HBkMWucnjyWlwkFOXikUlYJuSgX1+HpbC5KGZgpgHSSttNPzyT4z/4x/RF03htcf+xPoDVdSOKiD5nEXsVSZ2eDSecPR8QUgOtpJLJUV6DflUkUclecY2RjkySU0ah8M+EUfSJBz2cRiNA69BP2FeDiazgbce3cYr92/g0u/OxOro/30JPbEbjRQmGilM7H7GVjCsqfL5Kff42e/xsc/jY3e7j/Wudl452NzRzmJQjE20MNFuY6rDxrQkG9McNpLNn82JivH0uhLofBVnBNB14ctYm0qllAlIBhoRJ6Vgg4e2ldXUrztAeaCOA4mNVCU2EAqHsbZaOe200zjttNMoLCwkJSWl7zfsRd3eMt7+60PU7tlNwdTpnPe1b5Kef/RFxeYXXsD94Udk33knlsI+Kkz6WqGttttVngBe3/c6vpCPS0YvpLHxE1rbttPauo3W1u20t+8lNoYxm1OxJ44j3HoWdTudZOdN4YzF52K15R4RtFtaWti7dy9lZWWU7d3LFk85yhIm35rEyD0jGLmtkrSsdBxn5JE4KwtDwsBzyu5QiF3XfoX3Gz3s2rSb/bYsyooXE5ob2Z/Q7GVMyM3VowqYlupkXKKFsXYraWYTWofx+Wppb9+Hu30vbvcu2lp3UFPzEqGQO3oGhc02iiTHJJzO6TidM0hKmorJ1P2KV905bXYWF5oNvPHnrbz82/Vc+p8zsSf3PpV2qJgMilE2C6NsFs7hyBShOxhiV7uPUreHUreXUreXT5vbeLHu8MLoo6wJTE9KZHpSLPgnkp4w/IN/n7NuooF7F7AQqALWANdprbd1avNtYJrW+pboxdjLtdZXddp/I1Aks24+u7SOTI2sen8XpeW7KTfUU2eIzDBJSUlh0qRJTJw4kZEjR2IYghkRXncbnzz3FBvfWkGiM5lzr/8PJp65oNtRr7+8nH2XXY51+nQKHnsU1df5y1fC44vguudh/Beiny+M272b5pb1vLz1N+QaPaQaDq/tarHkkpQ05fCPYzImQxbvP1NKaUkt088bwVlXjDvqDt2uwuEwVVVV7HryP9lpmky9NxIk0gxJjPSnM9qUzei540man48pzdrre2mtqfIFWNviZk2LmzUuN9vaPLEbeslobmKCu4Wi+cUk6YP8+ZPvcWvzHNyrdmB1JHHu9f/BpLPO7fMvCa3DeL2VtLXtpLWtNPLYug2vN5bRNWC3j8XpnEGycwZO5wzs9vEYDL0HwNj0U3uKhcW3zSKpj897otT7A2xt9bClzcPm1na2tHoo9/o79udbzMxyJjLHaWeOM5FpSYnYTsA00qGYXnkR8Dsi0ysf01r/n1LqHmCt1nq5UsoKPAXMIjKSv6bTxdv9gBNIAJqBC7rM2DmCBPrhJeQLsv+DHWxfu5m93hqaDJEbfrIzs5k0JRLcs7Ozh+zORx0Os/2j9/jwmcfxuFzM/MLFnHn1l7Ekdj9iDHu97L/mWgI1NRS+/BLmvDjWf131MIF/34HrxsdpCZTT0rKBFtdGQqHIZ2sLgcU+iQl5F0ZGrI7JJCQcOXvH6w7w+sNbqN7dzLxLC5lz4aj+/Rs8eQn4Wmm84iVKS0spLS2lvLwcrTU2nUBBOINx+YVMPG8m9nHpKKXwh8NsbfWwxhUJ7Otc7dT4Ivl/m8HAbGcic5MjwWamMxHjsheou+d/yfr+90j72te4fPnlGJWRh6bdxzt//SM1ZaUUTJ3Owpu+RVpe96ta9cbvP4TLtTny07oJl2szgUBk9Gs0JpLsnEVyylxSUopIds7EaLQd9R410RvKEmwmLv72DDJGOPrdjxOhORBka5uHLa0eNrW2s97VzoFo8DcpmOKwdQT+omQ7BdaEY353sNS6Ef0SCoXYv2MPWz7eQFntPtrwooD8tFwmz5nGpMmTSE0d+hWF6vbt4d3HHqZ61w5yxo7n/Ju+RXZh7xO2qu+8k5YXX2LEw38i6dxze2wXCDTR3LyGpqYSmspfpM3QGr2yZMDhmEBy8iySnbN4aOcbvF+3jbeveAerqfsRZkt9O/96cDOuQx4WXj+J8cUDqMH+77tg1Z/gR1Vgisyx93g87N69m51bd1C2p4wWg6IhKQN32ggOZmexx2zGF/19HWlNYG6ynaJoIJlst2Hqpt5P1X/eRus77zDq6ad4w7GP//n0f3j0gkcpyp7D5rff5ON/PEnQ72Pu4iuZt+RKTAkDX5FKa43HcwCXaxMtLetpbllDW1spoFHKRFLSVFKS55ASDf5mc+T/ofoDrbz20Cb83hAX3DSF0Z/RBUzq/QHWu9pZG/0S3tjaTnsoMgso3WxiTnTUPzfZzkxnIolDPOqXQC/6FAgE2LNnD9vXbmHXvjK8IR9GbWCkI4fJM6Yw+YwZOBzHZrTlaXXxyXNPsentN7AlOTnnuhuZsmBhnymY5hdfpObOH5N+yzfIuu22I/b5/Ycigb15Fc1Nq2lzlwJgMFhIbtWkhFJJOefXOJ3TMZkin+uA6wCX/PMSbph8A7cXdb/4ds2eFlb8aTNaay66ZfoRteT7ZetLsOyrcPMHkDeTkNbscns7UjBrmt3sj44QDeEwmW3N5LQ0MjWoWVg4grPnTO1zGipAyOVi3+VL0eEQ+c89w4XvXsO0jGk8uDCyDKG7uYn3//ZXdn7yASk5uZx7/dcpnD13yEafgYCLlpZ1NLeso7l5DS7XZrSOfK7ExLGkpBSRmlKMWc3knUcPUl/RyplLxzJj4cjPfH2cYFhT2u5lXTTwr3e52d3uAyKj/mmORIqTI4G/ONlOlmVwF6Ul0ItuxUaQO7bvoGz3bgKhIAnaxEiVycSx45l8wWzsWcl9v9EAxWbTfPzsU/ja3cz6wheZf+V1WO19f6F4tm6j/EtfwjZ7FgV//Sv+UBPNTSU0Na+mubkEt3s3AAaDLTKKTC0mNWUezsQJGH4xGs74Lpz/P0e8592f3s2re17ljaVvkJmYedQ5d66s4f1nSnGkWvjirTNIyR747BNX/R7W//1m1hR9n7X2cax3uWmNjv4yzKbIaD3ZzlxnIlMdVmrW7Wbbqk2UNZbjUh4A8jJymDxzKhMnTiQjo+dRsGfLVsq//GWsU6bwxn+dzkPbH+HlS19mbOrhv5b2b97Au489TFNNFQXTZnLu9f9BZsHoAX++noRCPlpbt9DcvIbmlrW0tKwjGGwFwGoZgbt+PAd3FpA36izOufIcTCfZzU6NgeDhayotbja2tuON3h08yprAVTlp/L8xA1ulSwK96OByudi5cyc7d+5k//79hMNhErFQEMzgNOcIxp89jaQ5uYOa+dEXrTX7N63nw2cep+HAfkZOnsZ5X/0GGXEGFn9lJXu+fhW+sSES/uNsWjybaW+P3HFpNNpJTp5Naso8UlPnkZQ0FYOhUzqiegM8ci5c+QRMOVy/ps5dx6KXFrF03FJ+fPqRlS5DgTAfvbCbbR9WkT8hhS98fSo2R/wpDq01+z1+1rjcHb/kO91eNGDQYSYlJVIU/ZN+bh/53ECTl4r3drJ9y1b2h+poMESCZGZGBhOjF8Tz8vKOOt71+utU/dftWC/+AtfM+pSzRpx9VH39UDDAprdWsHLZP/C1tzPtvAs446ovYU85dgt/ax2irW1n9At6Nc1NawgEI3n+sD+NzOz5ZGSdTmrKPBITT76Sx/5wmC2tHlZH/7+Y4rBJoBf9p7Wmvr6e0tJSdu7cSVVVFQApJgcF3nRGhzMZOX40SWfmYxmbcsx/ker2lvHhM49zYOsmkrNzOPvaGxh/+ll9ntfrraapeTWNdR/RsPs1gmmRi5BGo6Pjz/+UlHkkJU3BYOjlT+B1T8Kr34XvrD+iTPD/rvxfXtr9Eq9e9iojkg5fmGxt9PLmX7ZSt8/FrAsKOH1xYZ+FuTyhMJtb21nT4maty82alnYOBSI3KiUZDRQl2yOBveQXzGrdTtLX3+rrn+0oOhCifWM9tR/voayhnHJzAzU0odE4nU4mTpzIxIkTGTVqFMZo6eKGhx+m/ne/p+yKufz3uA0su2QZE9ImHN3/tlZWLfsHG996DVNCAsVLrmL2RZdiTjj2UyAjs5/K2LXpHQ7s+RBbRikma2R2l9mcHvnvnFpMSkoxDvv4I+7qPdVJoD/FhEIhKioqOoJ7U1NkhJTtyKDAk0ZBeypp9lQc83Kxz83BlHLsf4Gbaqr49IW/s/OTD7AlOTl96TXM+PyFRy18HePxVNLcHE3FNK3G4z0AgMFnJGE3ZM+5lqxpS3E4Jvc5je8Iy78D216BH+6H6DWAvS17ufyVy7ly/JXcefqdHU13r63j/WdK0Vqz8PpJnDY7q9u3rPb6WRu9CLemxc3WNg+B6O9Voc1CUXJi9MKpnQl26+Fb7ru5INtfWmv8+120fVpN07YaKlQDB5zNVPgPEgwFsVqtjB8/nkmTJnHaaafR8JO7aHnlFf52kY3Wi8/oyNV3p7G6kg+efoy961bjSE1j3mVXM/W8CzCZh+4Gp9401rh545EttLXtY9LnDpE6ah8tzavx+iK38ZhMyaSkzI1+yc/t//8LJxkJ9KcAn8/Hnj17KC0tZdeuXXg8HoxGIwUZ+YwKZZBbbcMetmIZn4qjOAfrpDTUcZjr21hdyaqXnmPnxx9gNJuZc/Fi5l669IjpklprvN6KyIyY5hKam1fj9Ub+8oj9MqfYZ+P749sEV2xhxP2/w7loUf87ozXcPxXyZsI1h8v2fufd77Cmdg2vXfYa6bZ0/J4gHz23i52rasn+/+2de3QdR5ngf19337d0r96ynpZkO46dxK+8nBAISSYJCTAh2cwAB3YgPAYWGOCcHRYyO8POMofhNYcMsMyZyYHZkAkHCAOBQBIygYQkE/KwncSOH7Et25JsWdb7Svd9+1H7R7fka1myJVmSHW//zqnT3dVVt7/73arvVlVX1dce58YPXUSi1p0amLMdXktl2TaeZet4hpdLpjiGNWFDeXRyCGZTPEbNqRbT7PkV/OT98MFHoe1Nc/8+U7CSeTJb+sluPUZhLMfRaJLDlSm60r3kCnkMw2BFRwf1L22h6sknue96izs/9y9c03TNKT/38K4dPPfgA/S+vpvymlo23/EeLrr2BvTTrT5eAMyCzTM/3svrzx+jW1jn1gAAHZpJREFUtrWcP7prLZGK5PFGQPIlcrluwOvdJTZRUXElFZWXEy+/5MRhu/Mc39CfhyilGB4eprOzk/3799PV1YVt24TDYVY0tbPcrqGuO4iRA608QHRTPWVXLMOoPnku82Iw3HuYFx96cNLAb7j57Vz2jtuJVVSWLE7axlhyK6PJFykUjgEQCFRRUXGF10W/krLYBTiZHIc//jFyL79Cw5e/TMXt75qfUIN74btXwDv+ES67C4Cnep7i0099ms9s+gwfueQjdO0Y4ukf7SWTLLDpluXUX9/EK5nc5MyJXSWt9dZwkMsSMTZ50+YuKgvPbfvc/Dh8vR2u+hTc+L/n952mQTmKwv5RMluOkdszgmPbDC8r0hNPcmD0MOOpcUQpagcGGKwe4+Ofvoe62rpTDp9N+K597sEHONa5j0RdPZe94w4ueusNJ2xFsVgcfGWQp374OmbB5sp3drDuhmZ0r6GSLxxzX+4mXyKZ3FLyIj5MIrGRioorqay4nHh8A7p+bi7KWgh8Q3+eUCgU6OrqYv/+/XR2dpJMuntzVFdXs6K5jVZVR+VBQQ0VwNCIXFRNbFMdoZWViL74L7GUUhzetYNtj/yCgy9vwQiF2HDT29l4y01Y0k1y7GXGxrYxPv7q5EyLQKCaysorqay4koqKK4jFTnQ7ZyeT9HzsY+R37qLpG18nfuut8xfw+X+Cx++Gz74GFa2ki2lu++VtxINx7nvzAzzxy4M8150kuTxCbk2cnVaRIW9sPaJp3upHbwVkIkptcAGGMO57B2SH4RPPn/lnTYOdLpJ9eYDMtn6s/ixKU6RahZ6yYXbtfpGUN+5eUVHBihUrWLlyJe3t7ZN77U9FKcXBl7fw4s9/Ql/nXiLlcTa+7Z1suPntRMpnv0HdfMiMFfj9D/fStWOI6qYy3vq+1SzrOHlWmDu1dqvXO9xCOr0Hdy5/kER8vdtDrLySRHzjnLZuONfxDf0bFNu26e3tpauri0OHDtHT04Nt2wQCAdrb2+loWE5jNkGos4jZ5+5FEmyLE9tUT2RdDVp4acYrrWKRvc8/y7ZHfsFg90HKG8JceN1KqjqETG4X6fTruBvRCrHYKioSl5JIbCKR2EQkMvOK0vy+fRz51F9g9fXR9I/3UH7D9JuQzZoH/guMdsNfbGWwaPLfX/hn/jA0yMrYbRy0hPHo8db4ikiITQl3Nsyl8SgXTrMgaUH4w3fgP/4aPrsTKhbXMYh5LEP21QGyrw66DmACGodTLzM8uIcjbQmGY+UUi0VEhJaWlsk9ixoaGiZf6E6glKL39V1sefhnk3/qa998HetvvJW6tpl3BT1TlFIcfHWQZ3+yn8xYgTVXN3DlOzuIneI9k2mOkRzb6s3qeYlUehdK2ZOLuBKJjcTL3X17IpHWN+zMHt/Qv0GwbZu+vr4TDLtpuuO/9fX1dHR00F7ZTE0yQnHXKFa/u395cHncdcd3cc2SvFidYKD7EK898zOOdP4Wo2yERIsQrcujcLcS0PUY8fj6ScMej28gEJhdq2/8N49z9K/+Ci0Wpflb3yK6adO8ZHSU4nC+yJ7kGLsf/yo7mq9nR7Sdo4XjWwdXjdusUjrXXlDD5oYKLimPUL5U87cH98F3L4e3fxMu//CSPFI5imLPONlXB8nuGEBlbZRtYhZ7yb2lnaORLAcPH6Kvrw9w99lvbW2lra2N5cuX09jYiFEyPj90uJutv36Ivc89g2UWaVi5mnU33sLqq65ZtGGdYt7ipV8f4rWnjqDpwoY/amXjTa0EZ9G4say0u3I3+RKjyS2kUq/hON5CJqOChLdZmxsuIRh8Y6zU9Q39OUo+n+fIkSOToaenZ9IhR21tLW1tbbQ1tlJvJdAO5cjvG8VJmyCucY96xl1fgp3/lFIUiv0M92+lZ9/jjAxvQ48OEYjZXgqNWGwV8fglxMsvIZHYNKuNrabiZDIMfPMeRn/4QyLr19P07W8TqJ9+tstURk2LPek8ezK5yePrmTwZ+7gzihWGQ11Roe8eomMgwtpIiBvuWEXLhTN4olpslIJvb4R4E9z1yNI/3lG8vuNV9vz4l1yaXYcerQUUwdY4TluYY9EUR9L9dHd3Mzg4CLiGv7m5mba2NlpaWmhsbCQcDpNLp9jzzJNsf+IxRo4eIRSNserKN7HmmmtpXnsxmrbwf55jg1le+OVBOrcOEIoZrL++hUve2kw4NvthNccxyWT2Mz6+3du3Zwfp9D7cXiiEw83eRnZrKPNCOHzyWoWzjW/ozwEcx2FoaIjDhw9PGvaJigOuQ46WlhbaWpbTqFVhHLMoHByj2DMODmhRg9CqSsIXVhFeVYE+hwU7c8WyUqQz+0in95JJ72NsbBfp1OsocXsQSoGVjhGLrqGp/Xqqai+jvGzttJtWzYX0fz7HsS9+EbOvj8r3v5+6z/0l2pS9V5RS9BVMOrMF9mfzHMgW6Mzm2Z8tTM5+Aag0dNaURVgTC7M2FmbFb/4H6mgde7K3kBu1GS8b5NY/uZx1l7efdsfJRee5b7lTLT/2DDSsPysi/OrAr/juQ3fzv55dSR0rCV14LYi7vYOEdUIdFTjLQ/SHUxwePkp3dzf9/f2T+Wtra2lqapoM5sggu3//W/ZveQEznyNWWcXqq97MBZuvoWHVBQtu9Pu7xtn6aBddO4YIhHQuenMjF1/bRKJ2fquXbTvLeGqXZ/y3k0rt9mb3uPbSMOKe0b/Q/QMoX0MsugpdX7oe9VR8Q7/EFAoFBgYGOHbs2GTo7++f9LIUDodpbm52W0M1y6h14siRPIVDYxSPpMFRIBBoKiN8QSXh1VUEW8oX1CAppTDNYbLZLrLZQ+4e5Jn9pDN7J6c2AjiWQW4oQG4khFh1LGu5mlUb72BZ+9oFa9EUDhxg8NvfIfX44wQ7Olj2d39H9pJL6MkV6c67DiI6swU6M3k6c4XJjaLAXYC0IhpmZTR03LCXRaj3pjX2HRhj72+3s297GktFGEh0caB1C195/1/RklhaZ9kzkkvCN9fCmnfCHf9y1sS4f9f9fPPFr/OXu9u59LFDaPEaKu78GIGWDRQOprCT7vCGVh4ktLwcpzHEUDhDf36Y3r6j9Pb2ks26jQHDMGhoaKCuthYtn2Xs0H76X3sFxywSiSfo2HgZHZsuZ/m6TYTm4YxmJoaOpHn5N10ceHkQx1G0rq1i7ZsbWX5xNUbgzP5cLCtDJrOPVHoP6fQe0qk9pDN7S1xACpFwC9HYCmLRDu+4gmi046TdTxcD39AvEoVCgeHhYYaGhiZDf38/w8PDk2nC4TDLli2jvr6eZTV11OmVlI8HMHvTFHvT2CPenueaEGwuI9SRINSeILg8fsYvU11jPkI+30su1zNp0LPZQ2RzhyZnvgAIBhq1mONljHQXGO91yA2HqaheRfvGy7lg8zXUtS3MEnSlFCnbobvzIHt+9SidBw/RV9/I0MZLObaske6CeYIxF6ApHGCVZ9BXRsOsiIZYFQ1TFzROkEkpxWBPikPbh9i3pZ/xwRyGViRYvoWftG9Faovce+O9tMTPESM/wWOfhy3fh09tmdEhylJw/677+Yet/8C11go+81IV5tN/wFi2jKq77qLsulsxjxYpdo1T6EkdL7uGEGwuJ9BURjbhMCBjHEsNcrTvKAMDAxQKhcnPj0XCBCyT4lA/TmoMwzZpaG6hde3FtKy5mMYL1xIMn/kU4EyywO7njrLrmV4yY0WCEYMVG2tZeVkdTasq0QMLs4ZEKYdcrptUeg+ZjOttLJs9SDZ7cHLcH1wnNdHoCqLRdqKRVsLhZiKRViKRFgKBqgWpV76hnyfuQp48yWRyMoyMjEwa9VSqxFCKUFFRQX19PfU1ddRGKqlyyomMa9iDOcz+DPbo8R9erwxNVo5gcznB1vI57y/jOCbF4iD5/FEv9E6GnHftOLmSHEI43Ego2IIqJCiMhUgeydO3e4DxY3lQQrisnNZLNtC+4VLa1m2krKp6TjJZjmLYtBgyrUmfnb2eU2n3vMDRbJ7slKXrkUnPP65P0uWR0OSxJRw85ZaumWSB3v2j9O5L0rNzmPRoARFovKCSxtqd/PvI3/NYwuCy+sv4xrXfoCZyDr5cGzsC/3QV1K2FDz4C+tlb4fm7nt9x97N3o4vOFyN3cuFPt5F/+RW0aJTE7beTuON2wmvX4qRMij3jFLrHKXaPY/ZlUKb3B60LgbooRkOMfMJhRM8wZI4xlBqmf2CA4eFhSm2PWEW0Qh7NLJIoi1Fbv4ymtjbaLlxLY8dKNH1+rXHHdjiyd5R9L/Vz8JVBzIJNIKTTsraK5RdV07AyQUV9dMHH211nLUfJZDvJZg6SyR6Y/AMoFk90pa3rMSLhZsKRFqqr3kJz8/vm9Uzf0E+D4zjkcjlSqRTpdJp0Ok0qlSKVSp1g2Cdejk4QCoWoqa6mKlFFVSRBQotRYUcpywUhaWKP5LHHS/LoQqA2glEfI1AfJdhURqC5HH2Gl0WOY2KaSUwriVkcpVgcpFAcoFjwjsUhioUBCsVBTPNkb42BQCXhcBPhcBMaVTj5CIUxg8yQYrQnw8ChHtIjx3scifplNK1eS9OFa2lafRFVjU0nbA9ccBzGTJukZZM0LUYtm6GixWDRZMi0GCxa3rXFkGkyYtonyaQBtY5FbXKE6sM91AwPssyxabtoDavetJkVDXXUBIxZVTazYDN8NM3Q4TSDPSmO7k+SnJh9FNZpWl1J+/pa4quEB3fdw48PPYIlwkc3fIKPrPsoxrm8RH7Hg/Dzj7o7a974JTiLL/sOpw7zN8/9Ddv6t7G6cjWfDN/Mqt91kn70MZRpEmhpIX7zTcTe/BYiG9ajhUIoR2EN5TD70phHMxT7MphH0+4Eggk0cT1nVQdJxUzG9SxjKsNILsnASD+j42OYJb05AGyboEAkFCReXkZldQ11DY3UNTUTTyQoKysjGo2e1rOZVbQ58vooXa8N0e01CAAi5QGWdSRoWFlBXWs5VU2xOW1aN1dsO0sud4Rc/jC5XA/5kvPKis2sXv238/rchfAw9TbgW7gepr6nlPrqlPsh4H7gUmAYeLdSqsu7dzfwYcAGPq2UevxUz5qvobcsi7GxMXK5HLlcjmw2O3k+9XrCsDuOc9LnhAJB4rE48XAZ8WCMci1KmR2izAwRywUwUgqVtU7KpyeC6JVh9CoDqVZIjYNUORCzsVUay85gW+7RMsdcY26OHjfqphsmvBxNRSRAMFhDKFhLMFiLrlUgTjl2IYSZCVIY08gOO6SGx0kNDTLSf4yMozADQYqBEHYkSrC+kXBDE8H6RgI1dUhFFUUjQNqeMOYWSdNmzJow7Da5aXQ0QbmuURsMUBM0qA0a1AQMqiyTytFh4sf6SOzfS3zbFuK7d2M4NkZdHWXXX0f8ppuIbt48437zZtEmkywwPpRjfCjP+GCOsaEcI0czJAeyk27nQ1GDhhUJt+W+qoJcxQgvHnuRZ3uf4T+PPIujHG42hU/f+n1amq6YRSk6B/jVZ2DbfbD+vXDL1yC8eNtEnw5HOTze9TjfeeU7HE4dpiZSw40Vm7mhu5z6Fw+Qe+ElsCwkECC8fh2RdesJX7ia0OrVBNva0ELui0k7Y2IN5bAGc94xizmUwx7JH+8BeEhAw4oL46ECY2QYLY6SLCQZN9OkrQIFsXFmaN0HdI1wMEQsFqWsrIzyeILyRJxYWTmRSIRwOHxCKIw7DHRlOHZgjL7OMcYGj/d+o/EgVY0xqhpixGsilFeHKa8OE68OE4ouzV4/c+WMDL2I6Lg+Y2/EdQK+BXhvqTtAEfkEsK7EZ+ztSql3i8ha4EfAFUAj8FvgAqXUyc0+j/ka+u79nTzw03vRRCHiHA+aQ0jXCRk6QU0jpGmERScoGkElBG0h6CgMhIBoiChs3cTWTRzdxNEsnLCDCtrYQRsnYOMYthuvW1iaiYOF5RQo2gUsZWNjYKNNOerY6DjoWBgoiaKkDIcIjopgqxCWE8J2Aph2ENMOUHQMckWdnKmRt4SCrSgoRRHBMgwsPYClG9iGgakHcAIBbCOAaQQw5zCrIaprVBg6CUOnIqCT0DUSIlQIxHFIKIeEY1NeyBPPZqhMp6gYT2KMjmD29VPsH3DDwDB2JoujGSjRkcpqAmsvIrDqIrRVFyL1TZhFhVmwKOZtzIJNIWuRGy+SHS+STRXJjRcxCycWD90QyqpClC8LEavX0WpMrKo0o4FBulNddI11sW/kdfpz7iymBke4ZTzJHZFWlr/3Z1BeP+fydNZQCp7+Ovz+7yGUgA3vhY7roKoD4g0QOr2zkYXGdmyePvI0vz74a57rfY6slcUQgzWhVq4erGR1l0XN3n5CXccQ83gjSK+uJtDURKCxkUBDA3plJXpFAj1RgZ5IIGUxRII4BQNVACcHdtrBSZnYaRMnXcROm6jciQ0rE5uM5MlRJCsFMpIjQ56MypGjQF5MCppNUXNwTtMp0hB0TcPQdHRNR0NDlICtoUywi4KyNQQNHAPQ0DSDQDBEMBwgGA4RjgQIR0OEIiGCoYAXDPcYDhIKBwiGDQxDRw8Y6IZGIKCjGzpGUMcwdDRNQ9M0RGTew0hnauivAv5WKXWzd303gFLqKyVpHvfSPO85Ez8G1AJfKE1bmm6m583X0N93/z18tXk9CikJGg4CCE5JnALvKJP31Tm63alhmhi27QXHOyovQMAWDFswHJ2AA4bl3gtZiqDlEDQVIdMhaB0/Bk01eR2wFMe/+dQCJtOeKk9fSnTUGQyFOGJT1HPkA+PkAykvpMkbKXLBMdKhYUbDw4yExkGmL6dhR9FmWbQXi1yWL7A5l6elejWy+ZOw7t1ndaz7jOjbDs9+E/b9BqzjTsoJloERBj0AWsD9fqIvwjDP9J9XRLFVd3hJd9ivO+zTFMc097fRbUXDCCwfUCwbhZpxRd2Ye6weh+DJHeFpccQNSvPONR0Jl0M4joTjSCCKBCLu0YigGVE0I4KmT8RF0fQQogVRmo6lgak5FDApikUBi6KYFLEoiIWFjYmNJbZ77h0tbCxx3HvY2DJz73ahqCtE+MRXPj+vvKcy9LOpBU3A4ZLrI8CVM6VRSlkiMgZUe/EvTMnbNI2Afw78OUBra+ssRDqZqmiCtf2H3G69ctwpio6DUq6NEKUQB1AgjnLtxuR1ydEBHEEcz7YoN61QevQCJZ+NQncUmqPQlXvUHDdOVwpNMXltKIWmFLoC3RE0JegKDCXojqA7uNeO41U3xYlGd6Y/56nxpfmm3JPSo3dPBOVVMLeSCWjgaKBEcHRQhmAbgmMorIADugOagxLbVaBmozQHxEFpNoiN0k3QiyijCLqJ8s5FHBBXhAgQQU6yVyGpJkwtYdEJi05EdKq10GSo08JowSjEatxW77J1kDipiL3xaFgPf/oDMPOus5Sxw5Dqg1Q/WDmwTXAs9zhzB3l+nKLxFwSu9sIEWWUzpEwGlMlQo0Wu0SaDQ1bZ9GJzSCm3qpkORs7GyDluKDqIqdBNt4GiWwrDVOiWV6dK6qo4eUTl0ZyBkjig6IVJ2U+WWQARA5EAIgF0CaBJkIgWICoBEAOFDuL+aUrpueiAjoiGEreniug4mtuInKgnDoIStwHkiObGaaXx4onm1TEoiWMyziiZqbOQzMbQT/fXPlWdM6WZTV6UUvcC94Lbop+FTCfxx3d+iD+eT0Yfn3OZQBiWXwVcdbYlmZEo0OoFn3OT2YxXHAFKJx03A0dnSuMN3SSAkVnm9fHx8fFZRGZj6LcAq0SkXUSCwHuAh6ekeRj4gHd+J/Ckcgf/HwbeIyIhEWkHVgEvLYzoPj4+Pj6z4bRDN96Y+6eAx3GnV/6rUmqXiHwJ2KqUehj4PvBvItKJ25J/j5d3l4g8COwGLOCTp5px4+Pj4+Oz8Px/u2DKx8fH53ziVLNuzs05hT4+Pj4+C4Zv6H18fHzOc3xD7+Pj43Oe4xt6Hx8fn/Occ+5lrIgMAt1n8BE1wNBpUy09vlxzw5dr7pyrsvlyzY35yrVcKVU73Y1zztCfKSKydaY3z2cTX6654cs1d85V2Xy55sZiyOUP3fj4+Pic5/iG3sfHx+c853w09PeebQFmwJdrbvhyzZ1zVTZfrrmx4HKdd2P0Pj4+Pj4ncj626H18fHx8SvANvY+Pj895zhvG0IvIZ0Rkp4jsEpHPTnNfROTbItIpIjtEZFPJvQ+IyH4vfGBq3kWW632ePDtE5A8isr7kXpeIvCYir4rIgu/kNgvZ3ioiY97zXxWRL5bce5uI7PX0+YUllutzJTLtFBFbRKq8ewumMxH5VxEZEJGdJXFVIvKEV1aeEJHKGfJOW6ZE5FJPvk6vPM7Lx998ZRORDSLyvKfbHSLy7pJ794nIoRLdblgqubx0dsmzHy6JbxeRF738P/G2Q18SuUTkuhKZXhWRvIi8y7u3WPr6E+/3cURkxmmUM9XBeelLKXXOB+BiYCeuMxsD18n4qilpbgUew/VqtRl40YuvAg56x0rvvHIJ5bp64nnALRNyedddQM1Z1NlbgV9Pk1cHDgAduN7jtgNrl0quKenfievfYMF1BrwF2ATsLIn7OvAF7/wLwNemyTdjmcL1t3CVVw4fA25ZYtkumNAn0Aj0ARXe9X3AnWdDZ9699AzxDwLv8c7/GfhvSynXlN91BIgusr7WAKuB3wOXzZBvxjo4H329UVr0a4AXlFJZpZQFPA3cPiXNbcD9yuUFoEJEGoCbgSeUUiNKqVHgCeBtSyWXUuoP3nPB9Z/bvEDPPmPZTsEVQKdS6qBSqgj8GFe/Z0Ou9wI/WqBnn4BS6hncil3KbcAPvPMfAO+aJuu0Zcorb3Gl1PPKrYX3z5B/0WRTSu1TSu33zo8CA8C0qyWXUq6Z8Ho81wP/Pp/8CyzXncBjSqnsXJ8/F7mUUnuUUntPk3XaOjhffb1RDP1O4C0iUi0iUdzWe8uUNNM5MW86RfxSyVXKh3FbeRMo4D9EZJu4DtIXktnKdpWIbBeRx0TkIi/unNCZd/9twM9KohdTZwD1Sqk+AO9YN02aU5W1I9PEL6Vsk4jIFbitwQMl0V/2hnTuEZHQEssVFpGtIvLCxPAIUA0kvT99WFidzUlfuA6TpjYqFkNfs2GmMjYvfc3GOfhZRym1R0S+httySuN2Y6wpyc7IQfkiyuUKJ3IdrqG/piT6TUqpoyJSBzwhIq97LYClku1l3P0x0iJyK/ALXHeP54TOcIdtnlNKlbaIFk1nc2DJy9pc8XoX/wZ8QCnleNF3A8dwjf+9wOeBLy2hWK3eb9cBPCkirwHj06Rbcp15+roE15PeBGdTXwtaxt4oLXqUUt9XSm1SSr0Ftyu0f0qSmRyRL6qD8lnIhYisA74H3KaUGi7Je9Q7DgAP4XbXFozTyaaUGldKpb3zR4GAiNRwDujM46QW1mLrDOj3Kv1E5R+YJs2pylrzNPFLKRsiEgceAf7aG8YE3FatN7RZAP4vC6e7WclV8tsdxB2f3oi7eVeFiEw0OhdSZ7OSy+NPgYeUUmaJvIulr9kwUxmbl77eMIbea8EhIq3AHZzcxXoY+DNx2QyMed21x4GbRKTSe+t+Eyf+ay+qXF78z4H/qpTaVxIfE5HyiXNPrp0sILOQbZk35jfRzdeAYWbnEH7R5PLuJYBrgV+WxC26zjjR0f0HSp9fwrRlyitvKRHZ7On1z2bIv2iyeb/XQ7jvq3465d6E0RPccd2F0t1s5KqcGPrwGhNvAnZ77zKewh0fnzH/YslVwknvghZRX7Nh2jo4b32d7m3tuRKAZ3GdjG8HbvDiPg583DsX4Lu445GvUfI2G/gQ0OmFu5ZYru8Bo8CrXtjqxXd4ebYDu4D/eRZ09inv2dtxXxRfXZL3VmCfp88Fle10cnnXHwR+PCXfguoMt2L3ASZuC+rDuGOgv8PtZfwOqPLSXgZ873Rlyku309Pb/8Fbfb5UsgHv9/K8WhI2ePee9OrGTuABoGwJ5brae/Z27/jhKb/rS54ufwqElvi3bAN6AW3KZy6Wvm73zgtAP24jAdxZUo+erg7OR1/+Fgg+Pj4+5zlvmKEbHx8fH5/54Rt6Hx8fn/Mc39D7+Pj4nOf4ht7Hx8fnPMc39D4+Pj7nOb6h9/Hx8TnP8Q29j4+Pz3nO/wNHS5eixU+KaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.arange(9,11,0.01)\n",
    "for sig in np.arange(0.,5,0.5) :\n",
    "    mu0=10\n",
    "    sample=np.random.normal(mu0, sig, size=n)\n",
    "    lik=[]\n",
    "    for mu in x :\n",
    "      lik.append(likelihood(sample,mu,sig))\n",
    "    print(sig)\n",
    "    lik=np.array(lik)\n",
    "    plt.plot(x,lik/lik.sum())\n",
    "imax=np.argmax(lik)\n",
    "print('max: ',x[imax])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate that the maximum likelihood value is equal to the mean of the sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max:  9.399999999999991\n",
      "mean:  9.40355278188661\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXRU93338fdXO0ISkpAEQixiFWA2A8bY2Aa8gp3YjpcEn8ZxHKeu89hN06Y5T/IkJ03T5zknafs0fVKnySGOYztO7DhesYt3g/GGYyDsIPZFgDYWLYCElt/zx1zRiSyhbWbuLJ/XOTqM7ty58/HV6OOr3/zmXnPOISIisS/J7wAiIhIaKnQRkTihQhcRiRMqdBGROKFCFxGJEyp0EZE44Wuhm9mjZlZtZltDtL3XzOyUmb3SafmvzGyTmW02s2fNLMtbnm5mvzezPWb2sZmVhiKHiIgf/D5CfwxYEsLt/QtwdxfL/9Y5N9M5NwM4BDzkLb8POOmcmwD8BPhxCLOIiESUr4XunFsDnAheZmbjvSPt9Wb2nplN7sP23gYaulhe723bgEFAx6epbgEe924/C1zjrSMiEnP8PkLvynLgr51zc4C/B/4zFBs1s18DlcBk4D+8xSXAYQDnXCtQBwwNxfOJiERait8Bgnlj25cDfwg6UE737rsN+GEXDzvinLuhp2075+41s2QCZf4F4NdAV0fjOheCiMSkqCp0An8xnHLOzep8h3PueeD5gWzcOddmZr8HvkWg0CuAUUCFmaUAQ+g0BCQiEiuiasjFG+veb2Z3QmDM28xmDmSb3jYmdNwGPgvs9O5eAdzj3b4DeMfpbGUiEqPMz/4ys6eARUABUAX8A/AO8HOgGEgFnnbOdTXU0tX23iMwRp4FHCcwi+VN4D0gh8AQyybga865ejPLAH4DXEzgyHyZc25fqP77REQiyddCFxGR0ImqIRcREek/394ULSgocKWlpX49vYhITFq/fn2tc66wq/t8K/TS0lLWrVvn19OLiMQkMzvY3X0achERiRMqdBGRONFjoZvZKDNbZWY7zGybmf1NF+uYmf3UO2vhZjObHZ64IiLSnd6MobcC33TObTCzbGC9mb3pnNsetM5SYKL3dSmBeeSXhjytiIh0q8cjdOfcMefcBu92A7CDwEmtgt0CPOEC1gK5ZlYc8rQiItKtPo2hexeAuBj4uNNd589a6Kng06UvIiJh1OtC986E+BzwjY7ziwff3cVDPvURVDO738zWmdm6mpqaviUVEZEL6lWhm1kqgTL/rXfWw846zlrYYSRwtPNKzrnlzrm5zrm5hYVdzouX7rz67cCXiEg3ejPLxYBfATucc//WzWorgC95s13mA3XOuWMhzCmVWwJfIiLd6M0slwUErtO5xcw2esv+FzAawDn3C2AlcCOwBzgD3Bv6qCIiciE9Frpz7n26HiMPXscBD4YqlIiI9J0+KSoiEidU6CIicUKFLiISJ1ToIiJxQoUuIhInVOhRrqmljTW7aqg4dYaKU2d4fkMFdWdb/I4lIlHItysWyYW1tzueWXeYf3tzF9UNzTyddhaAv3tmE6nJxj2XlfK3101icLp+hCISoDaIQmfPtfHNP2xk5ZZKZo/O5Ue3T+eSD/MxgxevXcBTHx/ikff388b2Kh7/yjzGFgz2O7KIRAENuUSZs+fa+OKvPubVrZV898YpPPe1y7l68jCSzUjCmDUqlx/fMYNn/uoyGptbuePnH7LjWOdzpYlIIlKhR5HWtnYe+t0GNhw6ycN3zeYvrxpH4FQ6nzZvbD5/eOAyUpOTuO+xT6iub4pwWhGJNir0KPLvb+3m7Z3V/PCWadw0o+frg4wvzOKRe+Zy8kwLf/mb9ZxrbY9AShGJVir0KPHR3uP8bPUe7pwzkrvnj+n146aVDOFf75zJpsOnePid3WFMKCLRToUeBc6ea+Pv/7CJ0qGD+cHNF/X58TfNKOa22SX8bPVeNh4+FYaEIhILVOhR4OFVuzly6iw/um16v6ch/uDmiyjMSud7L26hvf1TF4sSkQSgQvfZvppGlq/Zx22zS7h03NB+bycnI5Xv3DiZrUfqeXZ9RQgTikisUKH77F/fKCctOYnvLJ0y4G3dPHMEc8fk8c+v76SxuTUE6UQklqjQfbSloo6VWyq578pxFGanD3h7ZsZ3b5pCbeM5Hv/wwMADikhMUaH76F/eKCcvM5W/vHJsyLZ58eg8rp5cxPI1+2ho0jlfRBKJCt0nWyrqWLOrhvuvGk92RmpIt/23106i7mwLj31wIKTbFZHopkL3yS/e3Ut2RgpfnD865NuePnIIi8sKeezDAzS1tIV8+yISnVToPthfe5qVW49x9/wxIT8673D/VeM5fvocz284Epbti0j0UaH74LEP9pOalMSXF5SG7Tnmj8tneskQHnlvn+aliyQIFXqENTS18Oz6Cj4zo5ii7IywPY+Z8dUrx7Kv9jTv7q4J2/OISPRQoUfYc+srOH2ujXsuLw37cy2dVkxBVjpPfnQw7M8lIv5ToUeQc44n1h5k5qhcZo7KDfvzpaUkseySUbxTXs3hE2fC/nwi4i8VegR9cuAk+2pO88VLQz+zpTt3XToaA57646GIPaeI+EOFHkFPf3KI7PSUXp3rPFRKcgexqKyI5zccoU1vjorENRV6hNSdbWHllmPcPGsEmWmRvZTr7bNHUlnfxId7ayP6vCISWSr0CHl501GaWtpZdknkhls6XDOliJyMFM1JF4lzKvQIeWnjESYWZTGtJCfiz52RmsxnZo7gta2VOgujSBxToUdAxckzfHLgJLdeXNLtRZ/D7fbZJZxtaePVLcd8eX4RCT8VegSs2HQUCJyv3C+zR+dROjST5zbo4hci8UqFHgErNh5lzpg8RuVn+pbBzLht9kjW7jtBxUnNSReJRyr0MNtZWc/OygZumeXf0XmHz11cAsALenNUJC6p0MNsxcajJCcZN06P3Nzz7ozKz2Te2PzzQ0AiEl9U6GHU3u54aeNRrphQQEHWwC8xFwo3TS9md3Uju6sa/I4iIiGmQg+jDYdOcuTUWW692P/hlg5Lpg3HDF7dWul3FBEJMRV6GL2y+RjpKUlcN3W431HOG5aTwdwxeazU9EWRuKNCD5P2dsdrWytZOKmQrPTIftS/J0unFbOzsoF9NY1+RxGREOqx0M3sUTOrNrOt3dy/yMzqzGyj9/X90MeMPZuP1FFZ38QNF0XP0XmHJdMCmTTsIhJfenOE/hiwpId13nPOzfK+fjjwWLHv9W2VpCQZ10wp8jvKp4zIHcTFo3M17CISZ3osdOfcGuBEBLLEDecCwy2XjR9Kbmaa33G6dNP0YrYdrefg8dN+RxGREAnVGPplZrbJzF41s4u6W8nM7jezdWa2rqYmfq9zubu6kf21p7k+CodbOmjYRST+hKLQNwBjnHMzgf8AXuxuRefccufcXOfc3MLCwhA8dXR6fWslZnDD1GF+R+nWyLxMZo4cokIXiSMDLnTnXL1zrtG7vRJINbOCASeLYa9tq+TiUbkU5WT4HeWCrps6jE2HT1Fd3+R3FBEJgQEXupkNN++csGY2z9vm8YFuN1YdPnGGbUfrzw9pRLNrvb8g3t5Z7XMSEQmF3kxbfAr4CCgzswozu8/MHjCzB7xV7gC2mtkm4KfAMudcwl688vVtgSGMaJyu2FnZsGxG5g3ire1VfkcRkRDo8RMvzrm7erj/YeDhkCWKce/srGbSsCzGDB3sd5QemRnXThnGU388xJlzrRG/1qmIhJY+KRpCjc2tfHLgBIvLom/ueXeumzqM5tZ23t+tC0iLxDoVegh9sKeWljbHwrLYmcEzb2w+2RkpvLVDwy4isU6FHkKry2vISk9h7ph8v6P0WmpyEovKinh7RzVt7Qn71odIXFChh4hzjtXl1SyYMJS0lNjarddOKeL46XNsPHzK7ygiMgCx1TxRbFdVI8fqmmJq/LzDoklFJCeZhl1EYpwKPURWlQfmcsfS+HmHIZmpzCvN520VukhMU6GHyOryaiYPz6Z4yCC/o/TL1ZOL2FXVyJFTZ/2OIiL9pEIPgYamFtYdOMmiGBxu6bDI+8tidbk+NSoSq1ToIfDBnlpa2935UoxFE4qyKMkdxKqd8XsWTJF4p0IPgdXlNWSnpzBnTJ7fUfrNzFhUVsiHe2tpbm3zO46I9IMKfYAC0xVruGJiAanJsb07F5cVceZcG5/sP+l3FBHph9huoCiws7KByvqmmB5u6XD5hKGkJSdpHF0kRqnQB2h1eWDMeeGk2H1DtENmWgqXjss/PwVTRGKLCn2AVpdXM6U4h+FDovtiFr21qKyIvTWnOXzijN9RRKSPVOgDUN/UwrqDJ+NiuKWDpi+KxC4V+gB8sLuWtnYXkx/37864gsGMzs9kVbmmL4rEGhX6AKwqryY7I4XZo3P9jhIywdMXm1o0fVEklqjQ+8k5x7u7arhyYgEpMT5dsbPFZUU0tbTz8f4TfkcRkT6IryaKoB3HGqiqb47pj/t3Z/64wCmANY4uEltU6P3UMbVv0aT4eUO0w6C0ZC4bN/T8lEwRiQ0q9H56t7yGqcU5FOXEx3TFzhaVFbK/9jQHak/7HUVEekmF3g91Z1tYf+gkiyfH39F5h46ZOxp2EYkdKvR+eN+brhiP4+cdSgsGUzo0k9W7NOwiEitU6P2wuryanIwULh4VP9MVu7KorIi1+45r+qJIjFCh95FzjtW7arhyUmHcTVfsbOGkQk1fFIkh8d1IYbDtaD01Dc1xObulM01fFIktKvQ+etcbU47Fi0H31aC0ZOaPG8q7mr4oEhNU6H20amc100pyKMqOz+mKnS2aVMi+2tMcOq6zL4pEOxV6H9SdaWHDoZMsioNzn/fW+bMv7tKwi0i0U6H3wXt7amh3xPX8887GFgxmVP4gDbuIxAAVeh+s2lnDkEGpzBoVuxeD7iszY9GkIj7cq+mLItFOhd5L7e3/fXbF5CTzO05ELSor5GxLG58c0PRFkWimQu+l7cfqqW1sjquLWfTWZeM7Lh6tYReRaKZC76VVOwNvCl6VAPPPO+u4eLTmo4tENxV6L63eVcP0kiEUZqf7HcUXCycV6uLRIlFOhd4Lp86c40+HTrI4AT5M1J2O6Yvv6mRdIlFLhd4La3bX0u5gYQKOn3cYX5hFSe4gjaOLRDEVei+sLq8mNzOVWXF+dsULCb54dHOrpi+KRCMVeg/a2x3vltdw1cTChJuu2NmisiLOnGtj3YGTfkcRkS70WOhm9qiZVZvZ1m7uNzP7qZntMbPNZjY79DH9s/VoHcdPnzs/hpzILvemL2ocXSQ69eYI/TFgyQXuXwpM9L7uB34+8FjRY3V5DWaJOV2xs8HpKVwyNk/TF0WiVI+F7pxbA1zoI4K3AE+4gLVArpkVhyqg31aVVzOjZAgFWYk5XbGzhZMK2VXVyNFTZ/2OIiKdhGIMvQQ4HPR9hbfsU8zsfjNbZ2bramqi/8/2k6fPsfHwqYSe3dLZovMXj47+n59IoglFoXf1TqHrakXn3HLn3Fzn3NzCwugfwlizuwbnSOj5551NLMpixJAMDbuIRKFQFHoFMCro+5HA0RBs13ery2vIy0xlxsjEna7YmZmxsCxw9sVzre1+xxGRIKEo9BXAl7zZLvOBOufcsRBs11ft7Y41u2pYOEnTFTtbVFZIY3Mr6w9q+qJINEnpaQUzewpYBBSYWQXwD0AqgHPuF8BK4EZgD3AGuDdcYSNpy5GO6YoaP+/s8vFDSUkyVu+q5rLxQ/2OIyKeHgvdOXdXD/c74MGQJYoS7+ys1nTFbmRnpDK3NI93y2v4ztIpfscREY8+KdqNd3ZWM3t0HvmD0/yOEpUWlRWxs7KBY3WavigSLVToXaiqb2LLkTqunqzhlu6cP/uipi+KRA0Vehc6LmZxzRQVenfKhmUzPCdDpwEQiSIq9C68vbOaktxBlA3L9jtK1Oo4++L7u2tpadP0RZFooELvpKmljfd313L15CLMNF3xQhZOKqShuZUNmr4oEhVU6J2s3Xecsy1tXK3hlh4tmFhASpKxSuPoIlFBhd7JOzurGZSazGXjNL+6Jzne9MV3dlb5HUVEUKH/Geccb++oZsGEAjJSk/2OExOunTKMXVWNuni0SBRQoQfZVdXIkVNnNbulD66dMgyAt3boKF3Ebyr0IG97QweL9XH/XistGMz4wsEqdJEooEIP8s6OaqaV5DB8SIbfUWLKtVOH8fG+E9Q3tfgdRSShqdA9J06fY8Ohk1w9eZjfUWLOtVOG0eqdnVJE/KNC97y7q5p2B9dq/LzPZo/OIy8zlbe2a9hFxE8qdM/bO6opzE5n2oghfkeJOclJxuLJRawqr6FVnxoV8Y0KHTjX2s67u2pYXFZIki5m0S/XThlG3dkWXfRCxEcqdOCjfcdpaGrl+qnD/Y4Ss66cWEBqsmm2i4iPVOjA69sqyUxL5oqJBX5HiVnZGalcPr6AN7ZXEbjmiYhEWsIXelu7441tVSwuK9KnQwdoybThHDx+hp2VDX5HEUlICV/ofzp0ktrGZq6/SNMVB+q6qcMwg9e2VvodRSQhJXyhv76tktTkwCwNGZiCrHQuKc3n9W0qdBE/JHShO+d4fVsVl48vICcj1e84ceGGi4azs7KB/bWn/Y4iknASutB3VjZw6MQZbrhIs1tC5QZv6EpH6SKRl9CF/vq2SswCY78SGiPzMpleMkTj6CI+SOhCf21rJXPH5FGYne53lLiyZNpwNh4+xbG6s35HEUkoCVvo+2tPs7OyQcMtYdCxT9/Ypg8ZiURSwhb6K5uOAnDTjGKfk8SfCUVZTCjK4r+2HPM7ikhCSdhCf3nzUS4pzaN4yCC/o8Slz8wo5pMDJ6isa/I7ikjCSMhC31XVwK6qRj4zY4TfUeLWZ2aMwDl0lC4SQQlZ6K9sOkqSwdLpGj8PlwlFWUwtzuFlb2hLRMIv4QrdOccrm48xf9xQirJ1qblwunnWCDYePsXhE2f8jiKSEBKu0Lcfq2df7WkNt0TATdMDbzi/vFlH6SKRkHCF/srmY6QkGUumabgl3EblZzJ7dC4rNqrQRSIhoQo9MNxylAUTCsgfnOZ3nIRw88wR7KxsYHeVTqkrEm4JVejrD57k8ImzfHamhlsi5cYZxSQZvLxZs11Ewi2hCv25DRVkpiWzVMMtEVOUncH8cUN5aeMRXclIJMwSptCbWtp4ZdMxlkwbzuD0FL/jJJTbZ4/k4PEzfHJAF5AWCaeEKfQ3tlfR0NzKHbNH+h0l4SydPpzBack8u/6w31FE4lrCFPpz6ysoyR3E/HFD/Y6ScDLTUrhpRjH/tfkYZ861+h1HJG4lRKFX1Tfx3u4aPndxCUlJ5nechHTHnFGcPtfGq1t0nnSRcOlVoZvZEjMrN7M9ZvbtLu7/spnVmNlG7+uroY/afy/+6QjtDm6bXeJ3lIR1SWkeY4Zm8uz6Cr+jiMStHgvdzJKBnwFLganAXWY2tYtVf++cm+V9PRLinP3mnOO5DRXMHp3LuMIsv+MkLDPjjtkj+WjfcZ0KQCRMenOEPg/Y45zb55w7BzwN3BLeWKGz8fApdlU1cvscvRnqt9vmjMQMnt9wxO8oInGpN4VeAgRPT6jwlnV2u5ltNrNnzWxUVxsys/vNbJ2ZraupqelH3L57cu0hBqclc8ssDbf4rSR3EJePH8oz6w7T1q456SKh1ptC7+pdxM6/jS8Dpc65GcBbwONdbcg5t9w5N9c5N7ewsLBvSfvh1JlzvLL5KLdeXEKW5p5Hhb+4dAxHTp1ldXm131FE4k5vCr0CCD7iHgn82dmWnHPHnXPN3re/BOaEJt7APLu+gubWdr44f4zfUcRz3dRhDMtJ54mPDvodRSTu9KbQPwEmmtlYM0sDlgErglcws+ALc94M7AhdxP5pb3f89uNDzBmTx5TiHL/jiCc1OYm75o3m3V01HDx+2u84InGlx0J3zrUCDwGvEyjqZ5xz28zsh2Z2s7fa181sm5ltAr4OfDlcgXvrw73H2V97mi/OH+13FOnkrnmjSUkynlyro3SRUOrVwLJzbiWwstOy7wfd/g7wndBGG5gn1x4kLzOVpdOKe15ZImpYTgY3XDScZ9ZV8HfXlTEoLdnvSCJxIS4/KXqs7ixv7qji83NHkZGqsohGd182hrqzLbqakUgIxWWh//qDAwB6MzSKXTo2n0nDsnjiowM6ra5IiMRdodc3tfC7jw9x0/RiRuVn+h1HumFm3LtgLFuP1PPh3uN+xxGJC3FX6L/7+BCNza3cf9U4v6NID26bXUJRdjq/eHev31FE4kJcFXpzaxuPvr+fKyYUMK1kiN9xpAfpKcl85YqxvLe7li0VdX7HEYl5cVXoL208SnVDs47OY8hfXDqa7IwUHaWLhEDcFHp7u2P5mn1MKc7hyokFfseRXsrOSOXu+WNYufUY+2v1QSORgYibQn9jeyV7qhv5q6vGYaaLWMSSexeMJTU5ieVrdJQuMhBxUeht7Y7/+8YuxhUO5jMz9EGiWFOYnc7n547k2fUVOle6yADERaGv2HSE3dWNfPO6MlKS4+I/KeE8uHgCSWb85K1dfkcRiVkx334tbe385M3dTC3OYem04X7HkX4qHjKIL19eygt/OkJ5ZYPfcURiUswX+jPrDnPoxBm+dUOZLgAd4762aDxZ6Sn8y+vlfkcRiUkxXehNLW389O3dzBmTx6Ky8F8wQ8IrNzONBxaO560dVaw/eMLvOCIxJ6YL/Zdr9lFV38y3bijTzJY4ce+CUgqy0vnxq+U6x4tIH8VsoR8+cYaHV+3hpunFzB831O84EiKZaSl849qJ/PHACV7efMzvOCIxJWYL/Z9e2U6SGd+9aYrfUSTE7po3mhkjh/BPr2ynvqnF7zgiMSMmC31VeTVvbK/i69dMZETuIL/jSIglJxn/+9Zp1DY2829vaBqjSG/FXKE3t7bxjyu2Ma5wMPddMdbvOBImM0bmcvf8MTzx0QG2HtGJu0R6I+YK/YUNRzhw/Az/ePNFpKXEXHzpg29eX0b+4HS+++JW2tr1BqlIT2KuET8/dxRP3ncpV07UNMV4N2RQKt//7FQ2HT7Ff67a43cckagXc4WelGRcobMpJoybZ47gllkj+Pe3d7Ph0Em/44hEtZgrdEk8/3TrNIbnZPCNpzfSoFkvIt1SoUvUy8lI5d+XzaLi5Bn+YcU2v+OIRC0VusSES0rzeWjxBJ7fcITffHTA7zgiUUmFLjHjb66dxNWTi/jBy9tZs6vG7zgiUUeFLjEjOcn46V0XM7Eoiwd/u4E91TrNrkgwFbrElKz0FB65Zy7pqUl85bF1VNc3+R1JJGqo0CXmjMzL5JdfmkttYzPLfrlWpS7iUaFLTLp4dB6Pf2UeVXVNLFu+liqVuogKXWLXJaX5gVKvD5S6LjAtiU6FLjFtbmk+T9w3j+ONzdz88Pt8tPe435FEfKNCl5g3Z0w+Lz10BfmD07j7Vx/z5NqDfkcS8YUKXeLC2ILBvPDgAq6cWMD3XtzKg7/dwPHGZr9jiUSUCl3iRk5GKo/ccwnfuqGMN7ZXcv1P1vDqFl3GThKHCl3iSnKS8eDiCbz811dQnJvB1367gS89+ke2H633O5pI2KnQJS5NHp7DC/9jAd+7aQqbDp/ipv94j7/7/Ub21jT6HU0kbFL8DiASLqnJSXz1ynHcOXcUP1+9l19/sJ/n/3SEhZMK+fKCUhZOLCQpyfyOKRIyKnSJe0MGpfLtpZP56pVj+d3Hh/jN2oPc++tPGJ6TwY3Ti/nszGJmjcrFTOUusU2FLgmjICudr18zkQcWjue1bZWs2HiUJ9ce5NEP9lOQlc6CCUO5YkIBl5TmM2ZopgpeYk6vCt3MlgD/D0gGHnHO/ajT/enAE8Ac4DjwBefcgdBGFQmNtJQkbp45gptnjqDubAtvba/ivd01vL/nOC9tPApATkYKM0bmUjY8m7EFgxlXMJixhYMZlp2hYRqJWj0WupklAz8DrgMqgE/MbIVzbnvQavcBJ51zE8xsGfBj4AvhCCwSSkMGpXL7nJHcPmckzjl2VTXyp0Mn2VRRx+aKU/z244M0tbSfXz89JYlhORkUZqdTlJ1OYXY6eZlpZGekkJORSlZGCtkZKWRnpDI4LZm0lCRSk5PO/5vu/Zus/ylIGPTmCH0esMc5tw/AzJ4GbgGCC/0W4Afe7WeBh83MnHMuhFlFwsrMKBueTdnwbJbNCyxrb3dU1jexv/Y0+2tPc/D4aaobmqlpaGZ3dSMf7j1O3dm+X+c0OclITTZSk5NIMiPJIMkMs0COJAPD+9eMpKT//j7JDDrWD/E+6BCO0aZwpI3VUbE7547ivivGhny7vSn0EuBw0PcVwKXdreOcazWzOmAoUBu8kpndD9wPMHr06H5GFomcpCRjRO4gRuQOYsGEgi7XaWt3NDa30tDU4v3bSmNTK6fPtXKutZ2WtnbOtbbT3NpOS5vjXGs759razt92zuGAdudod+AcOOdodw7n8JZ53xP4PnBfeI6XwrHZsGyT2D1ezMtMDct2e1PoXf0/sPOe7M06OOeWA8sB5s6dG7s/DZEgyUnGkEGpDBkUnl9Skd7qzQeLKoBRQd+PBI52t46ZpQBDgBOhCCgiIr3Tm0L/BJhoZmPNLA1YBqzotM4K4B7v9h3AOxo/FxGJrB6HXLwx8YeA1wlMW3zUObfNzH4IrHPOrQB+BfzGzPYQODJfFs7QIiLyab2ah+6cWwms7LTs+0G3m4A7QxtNRET6QifnEhGJEyp0EZE4oUIXEYkTOjlXrBg+3e8EIhLlVOixYumPel5HRBKahlxEROKECl1EJE6o0EVE4oQKXUQkTqjQRUTihApdRCROqNBFROKECl1EJE6YX6ctN7Ma4GA/H15Ap8vbRZFozaZcfROtuSB6sylX3/Q31xjnXGFXd/hW6ANhZuucc3P9ztGVaM2mXH0TrbkgerMpV9+EI5eGXERE4oQKXUQkTsRqoS/3O8AFRGs25eqbaM0F0ZtNufom5LlicgxdREQ+LVaP0EVEpBMVuohInIi6QjezvzGzrWa2zcy+0cX9ZmY/NbM9ZrbZzGYH3XePme32vsT4ifIAAAXISURBVO6JcK6/8PJsNrMPzWxm0H0HzGyLmW00s3URzrXIzOq8595oZt8Pum+JmZV7+/LboczVy2zfCsq11czazCzfuy9k+8zMHjWzajPbGrQs38ze9F4rb5pZXjeP7fI1ZWZzvHx7vNejRSqXmc0ys4+8/brZzL4QdN9jZrY/aL/O6muugWTz1msLev4VQcvHmtnH3uN/b2ZpkcplZouDMm00syYzu9W7b8D7rJtcd3o/o3Yz63Z6Yne/h/3aX865qPkCpgFbgUwCV1N6C5jYaZ0bgVcBA+YDH3vL84F93r953u28COa6vOP5gKUdubzvDwAFPu2vRcArXTw2GdgLjAPSgE3A1Ehm67T+Z4F3wrHPgKuA2cDWoGX/DHzbu/1t4MddPK7b1xTwR+Ay73X4KrA0grkmdexLYARwDMj1vn8MuMOvfebd19jN8meAZd7tXwBfi2SuTj/XE0BmqPZZN7mmAGXAamBuN4/r9vewP/sr2o7QpwBrnXNnnHOtwLvA5zqtcwvwhAtYC+SaWTFwA/Cmc+6Ec+4k8CawJFK5nHMfes8LsBYYGaLnHlCuC5gH7HHO7XPOnQOeJrBv/cp2F/BUCJ//POfcGgK/wMFuAR73bj8O3NrFQ7t8TXmvtxzn3Ecu8Nv2RDePD0su59wu59xu7/ZRoBro8pOD/TWAfdYl7y+Yq4Fn+/P4EOe6A3jVOXemr8/fl1zOuR3OufIeHtrl72F/91e0FfpW4CozG2pmmQSOxkd1WqcEOBz0fYW3rLvlkcoV7D4CR20dHPCGma03s/tDlKkvuS4zs01m9qqZXeQtC+f+6ks2vPuXAM8FLQ7XPuswzDl3DMD7t6iLdS70WqvoYnmkcp1nZvMIHNntDVr8f7yhmJ+YWXqIcvUlW4aZrTOztR3DGsBQ4JT3P3fwcZ8By/j0wUO49llPunuN9Wt/RdVFop1zO8zsxwSOhBoJ/PnR2mm1rsYq3QWWRypXIJzZYgKFfkXQ4gXOuaNmVgS8aWY7vf+jRyLXBgLnfmg0sxuBF4GJhHF/9SFbh88CHzjngo9wwrLP+ijir7W+8P5S+A1wj3Ou3Vv8HaCSQMkvB/4n8MMIRxvt/ezGAe+Y2Ragvov1/Npn04HXgxb7uc9C+hqLtiN0nHO/cs7Nds5dReBPmN2dVqngz4/0RgJHL7A8UrkwsxnAI8AtzrnjQY896v1bDbxA4M+siORyztU75xq92yuBVDMrIMz7qzfZgnzqiCmc+8xT5f1yd/ySV3exzoVeayO7WB6pXJhZDvBfwPe8oUcgcITqDUc2A78mtPutV9mCfnb7CIwfX0zgJFS5ZtZxEBnxfeb5PPCCc64lKG8491lPunuN9Wt/RV2he0dkmNlo4DY+/afRCuBLFjAfqPP+zHoduN7M8rx3ua/nz/8vHNZc3vLngbudc7uClg82s+yO216urYRIL3IN98bjOv48TwKOA58AE7130tMIlOoKQqgXP0vMbAiwEHgpaFlY95lnBdAxa+We4OcP0uVrynu9NZjZfG/ffqmbx4cll/fzeoHAe0l/6HRfR7EZgTHXUO633mTL6xiy8A4cFgDbvfcaVhEYv+728eHKFeRT79WEeZ/1pMvfw37vr57eNY30F/AesJ3An+jXeMseAB7wbhvwMwJjhlsIevcY+Aqwx/u6N8K5HgFOAhu9r3Xe8nHeYzYB24DvRjjXQ97zbiLwZu3lQY+9Edjl7cuQ5upNNu/7LwNPd3pcSPcZgV/gY0ALgSOi+wiMUb5N4K+Gt4F8b925wCM9vaa89bZ6++5hvE9dRyIX8EXvMRuDvmZ5973j/V5sBZ4EsiK5zwjM9tri/ey2APd1+rn+0duXfwDSI/yzLAWOAEmdtjngfdZNrs95t5uBKgIHAxCYmbSyp9/D/uwvffRfRCRORN2Qi4iI9I8KXUQkTqjQRUTihApdRCROqNBFROKECl1EJE6o0EVE4sT/Bz4aMVbSR3a7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,lik)\n",
    "mean=sample.mean()\n",
    "plt.plot([mean,mean],plt.ylim())\n",
    "print('max: ',x[imax])\n",
    "print('mean: ',mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2\n"
     ]
    }
   ],
   "source": [
    "uncertainty=sig/np.sqrt(n)\n",
    "print(uncertainty)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample=np.array([9.75138801,  8.72848018, 10.47220725, 10.52206875, 10.94855483,\n",
    "       10.70213686,  9.42958808,  9.39819465,  8.08670111,  9.79906428,\n",
    "        9.89062253,  8.87365529,  7.81200978, 12.00335652, 11.98713523,\n",
    "        9.77634694,  8.83557441,  9.64112154,  8.89143572,  9.5988839 ,\n",
    "        6.30222694,  9.79593013,  9.89055665,  9.29519617,  6.62671069,\n",
    "        9.69356416,  9.63405782, 15.16734247, 10.79155629, 11.56763008,\n",
    "       10.44287662, 10.64927055, 11.07173284, 12.21214659,  6.38186678,\n",
    "        5.1656226 , 10.39657475, 10.50973562,  7.81338867,  8.20394825,\n",
    "        8.02506472,  9.68805718, 10.29317738, 10.99694181, 10.99182911,\n",
    "        8.24196036,  9.05744792,  7.95655252,  9.06387778,  6.48429789,\n",
    "       11.92729567,  9.06172082, 10.91036909,  8.62835056,  9.96073153,\n",
    "       10.86420525,  6.92702787,  9.03256025,  9.23834502,  9.12067597,\n",
    "       12.68244086, 10.48175634,  5.98863088,  8.79955396, 10.56155435,\n",
    "        7.58495139,  9.74763   ,  6.07843654, 11.02994454,  5.52875172,\n",
    "        7.4839738 ,  7.32670021, 10.6497431 ,  8.12136434, 14.0732171 ,\n",
    "       10.49015059,  8.53914883,  9.8874647 ,  6.33685379, 10.59453499,\n",
    "       10.58877273, 10.38488386,  8.37622205,  9.94371813,  7.02291363,\n",
    "       11.8620325 ,  5.11087953, 10.64512124, 10.78919949,  7.46748207,\n",
    "        9.01515847,  9.03683969,  5.38810793, 12.06649875,  6.84187233,\n",
    "       11.6799807 , 11.70097749,  6.4195111 , 10.73156769, 10.06578757])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 9.75138801,  8.72848018, 10.47220725, 10.52206875, 10.94855483,\n",
       "       10.70213686,  9.42958808,  9.39819465,  8.08670111,  9.79906428,\n",
       "        9.89062253,  8.87365529,  7.81200978, 12.00335652, 11.98713523,\n",
       "        9.77634694,  8.83557441,  9.64112154,  8.89143572,  9.5988839 ,\n",
       "        6.30222694,  9.79593013,  9.89055665,  9.29519617,  6.62671069,\n",
       "        9.69356416,  9.63405782, 15.16734247, 10.79155629, 11.56763008,\n",
       "       10.44287662, 10.64927055, 11.07173284, 12.21214659,  6.38186678,\n",
       "        5.1656226 , 10.39657475, 10.50973562,  7.81338867,  8.20394825,\n",
       "        8.02506472,  9.68805718, 10.29317738, 10.99694181, 10.99182911,\n",
       "        8.24196036,  9.05744792,  7.95655252,  9.06387778,  6.48429789,\n",
       "       11.92729567,  9.06172082, 10.91036909,  8.62835056,  9.96073153,\n",
       "       10.86420525,  6.92702787,  9.03256025,  9.23834502,  9.12067597,\n",
       "       12.68244086, 10.48175634,  5.98863088,  8.79955396, 10.56155435,\n",
       "        7.58495139,  9.74763   ,  6.07843654, 11.02994454,  5.52875172,\n",
       "        7.4839738 ,  7.32670021, 10.6497431 ,  8.12136434, 14.0732171 ,\n",
       "       10.49015059,  8.53914883,  9.8874647 ,  6.33685379, 10.59453499,\n",
       "       10.58877273, 10.38488386,  8.37622205,  9.94371813,  7.02291363,\n",
       "       11.8620325 ,  5.11087953, 10.64512124, 10.78919949,  7.46748207,\n",
       "        9.01515847,  9.03683969,  5.38810793, 12.06649875,  6.84187233,\n",
       "       11.6799807 , 11.70097749,  6.4195111 , 10.73156769, 10.06578757])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}