{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:00:22.269979",
     "start_time": "2016-03-02T08:00:22.257626"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/robhicks/anaconda/envs/python36/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy.stats import norm,uniform,lognorm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sbn\n",
    "import pandas as pd\n",
    "import pymc3 as pm3\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "np.random.seed(348348756)\n",
    "\n",
    "sbn.set_style('white')\n",
    "sbn.set_context('talk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "For having some data to play with, let's resurrect our simple example yet again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:00:25.762915",
     "start_time": "2016-03-02T08:00:25.411628"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD8CAYAAAB9y7/cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADx5JREFUeJzt3XuQ3eVdx/F3SrgkYAsUQhIjwtTyZQZF6EyhME0RrFaKXOxFTSNailCmQmmVDgyV0ToNAUY6QhXC2EZ0iDhUMJRGCh0rWCxQqEIt4Be8VIsmTSFNLVkukcQ/fied7clefrvnsL/N1/drJsPus8/J+ezJ7mef/V0e5mzfvh1JUj2v6jqAJOmVYcFLUlEWvCQVZcFLUlFzuw4AEBF7Am8E1gMvdxxHknYVuwGLgIcy88X+D86Kgqcp9y91HUKSdlFLgfv6B2dLwa8HWLNmDQsXLuw6iyTtEjZs2MDy5cuh16H9ZkvBvwywcOFClixZ0nUWSdrVjHlo25OsklSUBS9JRVnwklSUBS9JRbU6yRoRbwauBg4HngGuyswbxpi3DFgBLADuAc7OzG8NLa0kqbVJV/ARsR/wWeBaYD/g3cDKiHhr37wjgVXAMuBAYANw3bADS5LaaXOI5keBdZm5JjO3ZeY/AH8LHN83bzlwe2Y+mJnPAxcDp0fEguFGliS1MWnBZ+YjmXnmjvd7K/qlwKN9Uw8HHh/1uGeBzb1xSdIMm9KNThHxGuAO4Ku9/462NzDSNzYCzJ92uik45JJ1037sN644ZYhJJGl2aH0VTUQcCnwZ2AS8IzO39U0ZAeb1jc0HnhsooSRpWloVfES8AXgQuAs4o3eMvd8TQIx6zAHA/r1xSdIMm/QQTUQcBHweuDozr5xg6s3AvRGxGngYWAnc2TsWL0maYW2OwZ9Nc9njZRFx2ajxa4DXAmTmeZn5SEScA6wGFtJs/3vWkPNKklqatOAz83Lg8jZ/WWbeAtwyaChJ0uDcqkCSirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJakoC16SirLgJamouVOZHBHHAGszc/E4H18HnAS8vGMsM/cZKKEkaVpaFXxEzAHOAj4B/O8EU48Clmbmw0PIJkkaQNtDNJcCFwIrxpsQEQuABcDXh5BLkjSgtgW/mmZ1/tAEc44Gvgd8LiK+HRF/HxHHDRpQkjQ9rQo+M9dn5vZJpu0F3E+z0l8C3ATcGRELB4soSZqOKZ1knUhm3g7cPmro+oj4AHAicPOwnkeS1M7QLpOMiHdFxC/2De8FvDCs55AktTe0FTywD3BFRHwdeAr4EDAPuHuIzyFJammgFXxErIqIVQCZeSNwDfB5YDNwGnByZm4ZNKQkaeqmtILPzHuAA0a9f17fx1cCK4eSTJI0ELcqkKSiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKmruVCZHxDHA2sxcPM7HlwErgAXAPcDZmfmtQUNKkqau1Qo+IuZExPuAu4E9xplzJLAKWAYcCGwArhtSTknSFLU9RHMpcCHN6nw8y4HbM/PBzHweuBg4PSIWDJhRkjQNbQ/RrAYuB06YYM7hwP073snMZyNic29847QTFnbIJeum/dhvXHFKZ8/dlUE/50F09Xp1+Tl3pcvvi0HMxtytCj4z1wNExETT9gZG+sZGgPnTSiZJGsgwr6IZAeb1jc0Hnhvic0iSWhpmwT8BfH+JHxEHAPv3xiVJM2xKl0lO4mbg3ohYDTwMrATuzMxnh/gckqSWBlrBR8SqiFgFkJmPAOfQnJDdCCwGzho4oSRpWqa0gs/Me4ADRr1/Xt/HbwFuGUoySdJA3KpAkoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqy4CWpKAtekoqa22ZSRBwN3AAcATwFnJeZD4wx7zHgUGBbb+g/MvOIIWWVJE3BpAUfEXsBdwArgE8BZwK3RcQhmfnSqHnzgAAWZea3X6G8kqSW2hyiORHYlpnXZ+bWzFwNPAuc2jfvJ4ANlrskzQ5tDtEcDjzeN5Y0h2tuHTV2NLA1Iu4Hfgz4R+DCzHxiGEElSVPTZgW/NzDSNzYCzB9j7kPAMuBg4GHgr3uHbiRJM6zNCn4E6C/p+cBzowcy8waaE7EARMRHgd8AjgLuHyymJGmq2qzgn6A5eTpa0HfYJiLOjYi3jhraDdgdeGGghJKkaWmzgv8isGdEXACsormK5iDgrr55i4ELI+LngGeAK4F/Bh4dXlxJUluTruAz80XgZJpj65uAC4DTMnNLRNwZEZf2pq6gKf2vABuB1wFnZOa2Mf5aSdIrrNWNTpn5NeD4McZPHvX2VuA3e38kSR1zqwJJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6SiLHhJKsqCl6Si5raZFBFHAzcARwBPAedl5gNjzPsQ8BHgh4DPAu/PzC3DiytJamvSFXxE7AXcAfwJsC9wLXBbROzRN+/nacr9ROBHgP2B3xt2YElSO20O0ZwIbMvM6zNza2auBp4FTu2bdybw6cx8MjO/C1wGnB0Ruw03siSpjTaHaA4HHu8bS5rDNbf2zfurvjmvAX4Y+M9JnmM3gA0bNrSIM44tm6b90Keffnr6zzuILjMP8Nxd6ezfCTp7vTr9nLuyK34vQye5R3XmmAvpNgW/NzDSNzYCzJ9k3o63++eNZRHA8uXLW0wd257TfiT89N0fH+DR09dl5kGeuytd/TtBd69Xl59zV3bF72XoPPci4F/7B9sU/Agwr29sPvDcJPN2FHv/vLE8BCwF1gMvt5gvSWpW7otoOnQnbQr+CeD8vrEA/nyMedE357vAf0/2BJn5InBfiyySpB+008p9hzYF/0Vgz4i4AFhFczL1IOCuvnk3Aasi4lbgmzRX0KzJzG3TiixJGsikV9H0VtcnA8uATcAFwGmZuSUi7oyIS3vz7gCuBNbRnFTdTHPZpCSpA3O2b9/edQZJ0ivArQokqSgLXpKKsuAlqSgLXpKKarWb5K4mIpbQXNL5FuB/gKsy89puU7UXEcfTbOp2GM3NXx/LzP77DmadiDgGWJuZi3vv7wesBk6iuSfiY5n56Q4jjmuM7EuAP6S5AW8r8Bngot5VZbNOf/5R468C/gb4amZe1Em4FsZ4/fcArqa5em8OzTYoH8jMl7pLObYxsi+m2X33zcALwI3AR7u4ZLzcCj4i5gBraW68ei3wNuB3e6U56/U2Z1sLXJGZrwZ+HfjTiDik02ATiIg5EfE+4G5g9C6jf0xzJ/NBwLuAqyLiyA4ijmuC7DcBT9PspXQU8EaaDfRmlQny7/BbNAudWWmC/Ctp9rs6DHh97+1Z9QNqguyfBP4FOJDm6+aXgF+Z+YQFCx44FlgMXNLb/fIx4Diazc92BfvSfGHM7f2w2ga8xOzewuFS4EJgxY6BiNgHOAP4ncx8ITO/QnP38zndRBzXWNn3ALYAH+9l3wCsAWbjImGn/Dv0fpiexQ9uAjjbjPX67w6cC5yfmZsycxPwTpp/g9lkvNf+MJqjIzv6dRvw/Azm+r6KBf8G4DGa1eKGiHgSeFNmPttxrlZ6Oa8DbqY5NPAlmi/0b3YabGKraVa5o/fDeD2wNTP/bdTYjl1IZ5OdsmfmS5l5Sq/YdzgVeHSmw7Uw1mtPROwJ/BlNUbbZD6or433tzAWOjYinIuK/gA/TYtuTGTbmaw9cRbOQGaG5q/++zPzMDGcDahb8/jR72D8DHAy8F/hkRCztMlRbvWOmI8C7aTZsOxX4g4j4yU6DTSAz12dm/x1ze7PzqmWsXUg7NU727+v9Gn4tzXbYK2cuWTsT5F8J3JWZs3qPp3Hy709zyONUmkMcbwJ+Frh4huNNaILXfg5wOfBqmgXN0oh4/4yG66lY8C8CmzJzZW8l9mWafetP7zhXW+8Ajs3Mv+zlX0ez/cOvdpxrqtruQjprRcQ8mpOrbwNOyMyNHUdqJSJOojmxPevOGbT0Ik03/XZmbu799voJmkN+s1pELKK5wOPKzBzJzMdptnA5t4s8FQs+gb0jYvQVQrvR/FTdFRzMzltLb+392ZU8BeweEQePGgt2/p/HzEoRsT9wL81q8rjM/PeOI03FLwOvAzZGxGbgPcD5EfG5bmO19hTNcet9R43tKt/Di2h++xh90rWz79+Kl0l+AfgOcEVEXAIcA/wC8DOdpmrvC8DKiDiL5vKqt9DkP6nLUFOVmd+LiNtpPpdzaH5VfQ/w9m6TTa53cvs2YAPwzszcpX64Zua5jFoxRsSNwDOz+TLJ0TJzc0SsBS6PiDNoDvd9mObKptnuMZqrr34/Ij5IU/gXAZ/qIky5FXxmPg/8FPDjwEaaKzc+mJkPdJmrrcz8J5pLCi+kuXb8j4Bfy8yHOw02PecAu9N8wd8KfCQzH+w2UivHASfQLAq+ExHP9f78Xce5/j95L/AkzW98X6NZ+FzdZaA2evdJvB04lOYelnuBvwCu6SKPu0lKUlHlVvCSpIYFL0lFWfCSVJQFL0lFWfCSVJQFL0lFWfCSVJQFL0lF/R9PHrLFg1sx1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f66a4233c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = norm(10,3).rvs(10)\n",
    "sigma = 3          # Note this is the std of x\n",
    "mu_prior = 8\n",
    "sigma_prior = 1.5  # Note this is our prior on the std of mu\n",
    "\n",
    "plt.hist(data,bins=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:00:27.024052",
     "start_time": "2016-03-02T08:00:26.986167"
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# a fairly basic mh mcmc random walk sampler:\n",
    "def sampler(data, samples=4, mu_init=9, sigma= 1, proposal_width=3, \n",
    "            mu_prior_mu=10, mu_prior_sd=1):\n",
    "    mu_current = mu_init\n",
    "    posterior = [mu_current]\n",
    "    for i in range(samples):\n",
    "        # suggest new position\n",
    "        mu_proposal = norm(mu_current, proposal_width).rvs()\n",
    "\n",
    "        # Compute likelihood by multiplying probabilities of each data point\n",
    "        likelihood_current = norm(mu_current, sigma).pdf(data).prod()\n",
    "        likelihood_proposal = norm(mu_proposal, sigma).pdf(data).prod()\n",
    "        \n",
    "        # Compute prior probability of current and proposed mu        \n",
    "        prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)\n",
    "        prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)\n",
    "        \n",
    "        p_current = likelihood_current * prior_current\n",
    "        p_proposal = likelihood_proposal * prior_proposal\n",
    "        \n",
    "        # Accept proposal?\n",
    "        c = np.min((1,p_proposal / p_current))\n",
    "        \n",
    "        accept = np.random.rand() <= c\n",
    "                \n",
    "        if accept:\n",
    "            # Update position\n",
    "            mu_current = mu_proposal\n",
    "        \n",
    "        posterior.append(mu_current)\n",
    "        \n",
    "    return np.array(posterior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:00:52.014564",
     "start_time": "2016-03-02T08:00:30.527650"
    }
   },
   "outputs": [],
   "source": [
    "chain = sampler(data,samples=5000,mu_init=9,sigma=sigma,proposal_width=3,\n",
    "                mu_prior_mu=mu_prior,mu_prior_sd = sigma_prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using our Chain for Inference\n",
    "\n",
    "In what follows, I assume we have already assessed chain convergence and are satisfied.\n",
    "\n",
    "In the code below, I ask a number of questions (you might call this inference) using our sampled posterior values, starting with ones most similar to frequentist statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:01:41.299619",
     "start_time": "2016-03-02T08:01:41.295147"
    }
   },
   "outputs": [],
   "source": [
    "# exclude first 100 iterations for burn-in\n",
    "chain_s = chain[100:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:01:41.986904",
     "start_time": "2016-03-02T08:01:41.978087"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean is:  10.180621628544616\n",
      "The standard deviation is:  0.8093876372087345\n",
      "The median is:  10.196002242838691\n",
      "The 95% Credible Interval is:  [ 8.51523653 11.77003799]\n",
      "The probability that the mean is greater than 10:  0.6015098959396041\n"
     ]
    }
   ],
   "source": [
    "print(\"The mean is: \", chain_s.mean())\n",
    "print(\"The standard deviation is: \", chain_s.std())\n",
    "print(\"The median is: \", np.percentile(chain_s,50))\n",
    "print(\"The 95% Credible Interval is: \",np.percentile(chain_s,[2.5,97.5]))\n",
    "print(\"The probability that the mean is greater than 10: \", \\\n",
    "            chain_s[chain_s>10].shape[0]/chain_s.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While the last question isn't one we typically ask in a frequentist setting, we could by imposing normality from our model estimates and calculating areas under the pdf.  Alternatively, we could perform a bootstrapping excercise as we did earlier in the class.  The posterior pdf pictured above does not assume normality of our parameters and in that way is similar to the bootstrap approach for inference. \n",
    "\n",
    "**The Credible Interval** above can be interpreted as\n",
    "> There is a 95% chance that $\\theta$ is in the interval [8.84,9.95]\n",
    "\n",
    "Note: the numbers in brackets will change depending on the credible interval printed above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Classical Hypothesis testing\n",
    "\n",
    "In frequentist statistics, we often perform test as follows:\n",
    "\n",
    "$$\n",
    "H_0 : \\theta = 0 \\\\\n",
    "H_1 : \\theta \\ne 0\n",
    "$$\n",
    "\n",
    "This isn't very easy to do in a Bayesian Context.  We can see that **no values** of our chain is exactly equal to 0 and this would likely be the case even if our mean for this problem was centered on zero. \n",
    "\n",
    "To make a hypothesis test more relevant for our problem, recast to\n",
    "$$\n",
    "H_0 : \\theta = 10 \\\\\n",
    "H_1 : \\theta \\ne 10\n",
    "$$\n",
    "\n",
    "To see this point, note that there are no values of $\\theta$ in our chain exactly equal to 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:01:56.284616",
     "start_time": "2016-03-02T08:01:56.280086"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cases where theta=10:  0\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of cases where theta=10: \", chain_s[chain_s==10].shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So does this mean we have 100% evidence that theta isn't 10?  A tough question to answer.  A better question might be to recast this question into an interval:\n",
    "\n",
    "$$\n",
    "H_0 : \\theta \\in 10 \\pm .1 \\\\\n",
    "H_1 : \\theta \\notin 10 \\pm .1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:02:04.297546",
     "start_time": "2016-03-02T08:02:04.292578"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of H_0:  0.09610283615588655\n"
     ]
    }
   ],
   "source": [
    "print(\"Probability of H_0: \",chain_s[(chain_s<=10.1) & (chain_s>=9.9)].shape[0] / chain_s.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another type of question we might answer are hypothesis tests over intervals.\n",
    "\n",
    "A more common question answered in a Bayesian setting is:\n",
    "$$\n",
    "H_0 : \\theta < 10 \\\\\n",
    "H_1 : \\theta \\text{ not } < 10\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:02:35.932712",
     "start_time": "2016-03-02T08:02:35.926014"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of H_0:  0.3984901040603958\n"
     ]
    }
   ],
   "source": [
    "print(\"Probability of H_0: \",chain_s[chain_s<10].shape[0] / chain_s.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So Bayesians like to attach probabilities to statements.  Additionally note that due to our very small sample size, the prior is definitely affecting our inference here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Highest Posterior Density (HPD) Intervals\n",
    "\n",
    "For skewed or multimodal posteriors, credible intervals may not be as enlightening as they could be.  Suppose we have a MH markov chain and our posterior pdf is plotted below.  The HPD interval, sets a minimum probability threshold for parameters (e.g. posterior value must be greater than .025) and examines the interval of values of $\\theta$ that are  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:03:07.184556",
     "start_time": "2016-03-02T08:03:06.496552"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEgCAYAAABcnHNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGQlJREFUeJzt3X2YnHV97/F32CSbbCCEhxACCEElX5SjwkHrEQELWNtIQYp6FdQAeipwlIjUZwStnlIexCe0AlbxodJqKxwsINJzpCoeRRFrcw7g1yBGiBB5CAGSTZaQbP+4Z5PdYXd2sjub2cnv/bquudi59zuT77Czn/3N777v3z2lv78fSVJZdmh3A5Kkbc/wl6QCGf6SVCDDX5IKZPhLUoEMf0kq0NR2NyBtjYj4HvCKQZv6gfXAr4DPZOYXW/TvPB84IDO/Nc7n+R6wMjNPakVfUqs48lcn+hYwv3bbC3gBcCPwhYh4bYv+jRuBl7XgeU4EzmjB80gt5chfnWh9Zq6s2/bBiHgd8Cbgmhb8G1Na8Bxk5qpWPI/Uaoa/ticbqaaAiIjFwF8CATwGfB04LzPX1b7/KuCvgf8CrAO+C5yTmb+LiOXAfsD7IuKkzFxQe8ypwHuB5wD3AVcDF2bmU7Xv99ee8w3AbOB44EIGTfvUppMuAA4HZgK3Au/LzKW1738ZmEP1u/kK4POZ+a6W/59S8Zz2UceLiNkR8QHgecA3IuKdwFXA3wMvAt4OnAT8U61+V6qpo+8ABwF/DCysPQbgJcAK4LO1r4mI04HPABcDzwfeCZwCfKWunbOoPn38KXB7XZ/7AT8CptX+zcOBp4BbI+LZg0pfA/wUOBj427H9X5Eac+SvTvTaiFhT+3oHqhH0SuBdVKH+IPB3mfmJWs2y2qj8uog4hGon8YzaY36bmb+JiNcDuwFk5sMRsRFYm5kP157jfOCSzPxq7f69tef8dkR8IDOX17b/Y2b+eKDRiBjc99uowv7PM3Nt7fsnAb8Gzq7dANYC/zMzXXhLE8aRvzrRzVSj4oOppm32yMz5tbDfA5gH/KDuMd+v/fdFmfkL4GtUo+qHI+KfgJcD/z7cPxYRc4F9gPMiYs3ADfhmreR5g8qXNej7hcDPB4IfoDYN9VOqTygDfm3wa6I58lcnWpOZ92zlYwYGOn0Ambk4Ij4CLAKOBj4HnBURh2fm+hEe+z6qo4DqPTjo63Vb2dfA8/eN8zmkreLIX9uVzPw98HvgyLpvDZwbcGdEvDAirgDuy8zPZOafAccBh1KNzqGaGhrwUO323My8Z+BGtVP4Y8COTba3FDg0ImYNbIiImVT7Fe5s+kVKLeDIX9uji4CPRcQy4AaqHbSfBb6TmUsjYi+qnbLTI+ISqqA/heqooF/WnuNJ4ICI2Lt2BNBFwCW1I4G+RXXEzxeAu2t/cJrxOapj/r8eER+q/bt/RfXH43Pje8nS1nHkr+1OZn6K6qib06lG1J+lOizztbXvPwC8GjgA+AlwB7A38MrMfKL2NB+nmg5aGhFTM/OTVEcNvQW4C/gS1R+B121FX7+l+kQyjeoQzx/Uvn75GKaxpHGZ4pW8JKk8jvwlqUCGvyQVyPCXpAJNiqN9IqKb6nC3B6nWZ5Ekja6LanXb2zOzb7TiwSZF+FMF/63tbkKSOtQRwA+35gGTJfwfBLj66qvZc889292LJHWElStX8sY3vhGGnmXelMkS/hsB9txzT/bZZ5929yJJnWarp8vd4StJBTL8JalAhr8kFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kq0GQ5yWvMFrx/uEuqjm75Rce2uBNJ6hyO/CWpQIa/JBXI8JekAhn+klQgw1+SCtTU0T4RcQhwJXAQsAw4MzNvG6buTmB/YFNt028z86AW9SpJapFRwz8iZgDXAxcAXwAWA9dGxILMfGpQ3UwggPmZ+fAE9StJaoFmpn2OAjZl5uWZuSEzrwIeBY6rq3sBsNLgl6TJr5lpnwOBu+q2JdUU0DWDth0CbIiIHwPPBf4dODsz725Fo5Kk1mlm5D8L6K3b1gv0DFN7O3AysC/wM+DbtekgSdIk0szIvxeoD/AeYM3gDZl5JdVOYQAi4oPA24GDgR8308zhF/8bzNq1mVJJ0jg0M/K/m2pH7mBB3VRQRJweEa8ctKkLmAasH1eHkqSWa2bkfwvQHRFLgCuojvaZB9xcV7cXcHZE/AnwCHAx8EvgP1rXriSpFUYd+WdmH7CIai5/FbAEOD4z10bETRFxbq30Aqo/CD8FHgKeA5yQmZuGeVpJUhs1dZJXZi4FDhtm+6JBX28A/rJ2kyRNYi7vIEkFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUoKnNFEXEIcCVwEHAMuDMzLytQf1bgEsyc/eWdClJaqlRR/4RMQO4HvgSMAe4DLg2IqaPUP9s4BOtbFKS1FrNTPscBWzKzMszc0NmXgU8ChxXXxgRXcBXgc+3tk1JUis1E/4HAnfVbUuqKaB67wfuBL49zr4kSROomTn/WUBv3bZeoGfwhog4FFgMvLh2kyRNUs2M/HuBmXXbeoA1A3ciYibwFeAvMnMNkqRJrZnwvxuIum3B0KmgFwPPBm6IiNXADcCuEbE6IvZtSaeSpJZpZtrnFqA7IpYAV1BN7cwDbh4oyMxbGTQNFBF/CHzTQz0laXIadeSfmX3AIuBkYBWwBDg+M9dGxE0Rce4E9yhJarGmTvLKzKXAYcNsXzRC/fcAR/2SNEm5vIMkFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUoKau4bs9WvD+G7f6McsvOnYCOpGkbc+RvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFcjwl6QCNbW8Q0QcAlwJHAQsA87MzNvqarqBTwKvB6YD3wPelpm/a2XDkqTxG3XkHxEzgOuBLwFzgMuAayNiel3p+cDzgQDmAo8Cn2lpt5Kklmhm2ucoYFNmXp6ZGzLzKqpgP66u7sPAosxcBcwDZgOPtLRbSVJLNDPtcyBwV922pJoCumbzhsyNwLqI+CvgQ8ADwJGtaVOS1ErNjPxnAb1123qBnhHqL6o95hrg5oiYNvb2JEkToZnw7wVm1m3rAdYMV5yZ6zNzHfAeYD/gBePqUJLUcs2E/91UO3EHC+qmgiLiqoj4H4M2Ta09/+pxdShJarlm5vxvAbojYglwBbCYaofuzXV1PwXeExE3AQ8BnwZuzcx7W9ivJKkFRh35Z2YfsAg4GVgFLAGOz8y1EXFTRJxbK70S+Arwf4HfUk0NvX5CupYkjUtTJ3ll5lLgsGG2Lxr0dT/w0dpNkjSJubyDJBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFcjwl6QCGf6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBpjZTFBGHAFcCBwHLgDMz87Zh6s4DTgdmA78AzsrM/9+6diVJrTDqyD8iZgDXA18C5gCXAddGxPS6utOAU4A/BHYH/g9wY0T46UKSJplmgvkoYFNmXp6ZGzLzKuBR4Li6ut2BCzLz3sx8Gvg0sC+wT0s7liSNWzPTPgcCd9VtS6opoGs2b8i8tK7meKo/EivG0+BksuD9N271Y5ZfdOwEdCJJ49NM+M8Ceuu29QI9Iz0gIo4ErgDOyMxNY29PkjQRmpn26QVm1m3rAdYMVxwRi4EbgSWZ+Q/ja0+SNBGaCf+7gajbFjxzKoiIOB/4FPCazPzyuLuTJE2IZqZ9bgG6I2IJ1VTOYmAecPPgooh4M3AOcFhm/rLVjUqSWmfUkX9m9gGLgJOBVcAS4PjMXBsRN0XEubXSDwA7AT+LiDWDbs+bqOYlSWPT1ElembkUOGyY7YsGfb2whX1JkiaQJ2BJUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUIMNfkgpk+EtSgQx/SSqQ4S9JBTL8JalAhr8kFWhquxvY3i14/41jetzyi45tcSeStIUjf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klQgw1+SCmT4S1KBmjrJKyIOAa4EDgKWAWdm5m0N6j8NbMjMd7ekS0lSS4068o+IGcD1wJeAOcBlwLURMX2Y2t0i4svAO1rcpySphZqZ9jkK2JSZl2fmhsy8CngUOG6Y2h8CTwPXtLBHSVKLNTPtcyBwV922pJoCqg/5YzLzgdroX+MwljWBXA9IUrOaGfnPAnrrtvUCPfWFmflAK5qSJE2sZsK/F5hZt60HWNP6diRJ20Iz4X83EHXbgmdOBUmSOkQzc/63AN0RsQS4AlgMzANunsjGJEkTZ9SRf2b2AYuAk4FVwBLg+MxcGxE3RcS5E9yjJKnFmjrJKzOXAocNs33RCPWnja8tSdJEcnkHSSqQ4S9JBTL8JalAhr8kFcjwl6QCGf6SVKCmDvVUZxjLYnDggnBSiRz5S1KBDH9JKpDhL0kFMvwlqUCGvyQVyPCXpAIZ/pJUII/zlxeLlwrkyF+SCmT4S1KBDH9JKpDhL0kFcoevxsSdxFJnc+QvSQWaVCP/px9/iP6n+ob9XlfPbHbonrX5/obVK6G/f8Tn2mHmTnTN2HHQc/+e/k2bRq6fMYuumbO31D/xMP0bnx65vruHrp6dt9Q/+Qj9T28YuX76DLpm7bKlfs0q+jcM/1oBpkzrZuqOu26+v3HtY2x6av3I9VOnMXWn3bfU9z7Opr7ekeu7pjJ19twt9eueYNP6tSPX79DF1J332FK/fg2b1j05Yj1TpjBtzp6b727qW8veZ3xh5Hpg2i7zB9X3srH3cQB+8N6jhq3ff//92WGHavyydu1aVq5c2fD599tvP6ZOrd7y69at44EHHmhY/6xnPYvp06cD0NfXx4oVKxrW77333syYMQOADRs2cN999zWsnz9/Pj09PQBs3LiR5cuXN6yfN28eO+5Yvaf7+/u59957G9bPnTuX2bO3vKfvvfde+hv8zuy2227MmTNn8/3ly5ezcePGEet32WUXdt11y3v0vvvuY8OGkX8Hdt55Z3bffct7dMWKFfT1jfw7sNNOO7HHHlvecw888ADr1q0bsX7WrFnsueeW99zKlStZu3bk9/TMmTPZa6+9Nt9/6KGHePLJkd/T3d3d7LPPPpvvP/LIIzz++OMj1k+bNo1999138/1Vq1bx2GOPjVjf1dXFggULNt9fvXo1jz766Ij1U6ZM2fz+HJP+/v623xYuXLhg4cKF/VOnTu0Hhr3t+kdn9u/3vhs236ZM7xmxFuifc+QpQ+q7dtq9Yf3sl75uSP203fdtWL/jwYuG1Hfv/byG9T3Pe8WQ+hnPPrRh/cznvGRIfc+BRzSs797n+UPqd3zRnzSsnzZ3wZD62X9wYsP6rtlzh9TvfMSbGtbv0D1rSP0urzyjYT1MGVK/27HnjFJP/xNPPNE/4Lrrrhu1/v77799cf8stt4xaf+edd26uv+OOO0at/9GPfrS5/p577hm1/qabbtpc//DDD49a/41vfGNzfV9f36j1n//85/sHmzFjRsP6j33sY0Pq582b17D+vPPOG1J/wAEHNKx/xzveMaT+xS9+ccP6U089dUj9Mccc07D+hBNOGFJ/4omN39NHH330kPrTTjutYf2hhx46pP7ss89uWP/c5z53SP2HPvShhvVz584dUn/ppZc2rO/u7u6///77+xcuXNi/cOHCBf1bmbtO+0hSgab0N/gYuK1ExALgN2tfejr9M3cetsZpn86f9tnY+8TI9Tjt47TPUE77NDftc8wxxwDsn5nLRywe7vGTKfz7XnUezNp1tHJpVB5ZpBKsWLFizOHvtI8kFWhSHe0jtcpYL2Y/Fn7KUCcy/KVx8oQ3daKmwj8iDgGuBA4ClgFnZuZtw9S9E3gPsBPwL8AZmTnyHhepUGP9ZOIfDbXKqOEfETOA64ELgC8Ai4FrI2JBZj41qO5PqYL/KOD3wD8CHwXeNQF9S0XaVtNZ/pHZ/jUz8j8K2JSZl9fuXxUR5wDHAdcMqlsMfDEzfwUQEecD342I92bmyMeLSZp0tuU+k23JP2pbNBP+BwJ31W1Lqimga+rq/lddzc7A3kDjA56hC4B1q5toR5LGZsGSv293C621bjXd1VddW/vQZsJ/FlB/tlAv0DNK3cDX9XXDmQ/QfetnmyiVJNWZD/x6ax7QTPj3AjPrtvUAa0apGwj9+rrh3A4cATwIOEUkSc3pogr+27f2gc2E/93AWXXbAviHYeqiruZxoPE59EBm9gE/bKIXSdJQWzXiH9BM+N8CdEfEEuAKqh2784Cb6+q+BlwREdcA91Md6XN1Zo68oI4kqS1GXd6hNipfBJwMrAKWAMdn5tqIuCkizq3VXQ9cDNxItYN3NdWhn5KkSWZSLOwmSdq2XNhNkgpk+EtSgQx/SSqQ4S9JBWr7ks4RcRhwGbCQ6iSvj2Rm/TkEHSci/gC4LjP3qt3fBbgKOJrq/IePZOYX29jiuAzz+vYBPkt1st4G4J+Bd9eOFus49a9v0PYdgO8Cd2Tmu9vSXAsM8/ObDnyc6qi+KVRLtbxt8OKNnWKY17YX1arEhwPrgS8DH+y0w9Aj4nCqn9GBwCPAJZl55Vizpa0j/4joAq4DLsrM2cBfAF+pXdaxI0XElIh4C/CvwPRB3/o7qrOd5wGvAy6JiBe2ocVxafD6vgasoFrL6WDgJcD5277D8Wnw+ga8Czhy23bVOg1e34VU63UtBA6ofd1Rf9wavLbPAPcAc6nel38OvGnbdzh2tYD/F6qB8i7A64ELI+KVjDFb2j3tM4fqBzI1IqYAm4Cn6OwlHs4FzqZaAhuAiNgROAH4cGauz8yfUp0h/db2tDguw72+6cBa4K9rr28lcDVwWHtaHJdnvL4BtV+oNzN0AcNOM9zPbxpwOnBWZq7KzFXAa6l+hp1kpJ/dQqpZjoG82wSMfCX4yWk/4MbMvDozN2Xmz4F/o/odG1O2tDX8M/NR4HNUa/9vAG6legPe386+xukqqpHv4LU2DgA2ZOa9g7YNrIzaaZ7x+jLzqcw8thb6A44D/mNbN9cCw/38iIhu4KtUIdnMelWT1Ujvz6nASyNiWUT8DjiHJpZmmWSG/dkBl1CFYS/V6gM/zMx/3sa9jUtm/iIzFw/cr30SOIJqim5M2dLuaZ8dqH4gr6daCO444FMR8aJ29jUemflgZtafOTeLZ440hlsZddIb4fVtVvvofRnVvOSF266z1mjw+i4Ebs7Mjl6DaoTXtyvVNMlxVNMi/w14FfC+bdzeuDT42U0B/gaYTRWKR0TEGdu0uRaKiJ2pLrB1B9Xof0zZ0u5pnxOBl2bmN2ujxxuploc4pc19tVqzK6N2tIiYSbWj94+BV2TmQ21uqSUi4miqnWkdtw+jSX1UWXBeZq6uffL+BNV0QkeLiPlUa5JdnJm9mXkX1TI0p7e3s7GJiP2BH1EttXMiVYaMKVvaHf77wsC1CDbbULttT5YB0yJi30HbgmdeJKdjRcSuwPepRpEvy8zftLmlVjoJeA7wUESsBt4AnBURN7S3rZZZRjUPPmfQti6qEXOnm0/1qWbwDuCOzJiI+K/AT6gW1TwhM9cxjmxp96Ge/5tqj/WbqQ6/OhL4M6pR1nYjM5+MiG9Rvda3Un30fAPw6vZ21hq1nfXXAiuB12Zmx/1iNZKZpzNopBgRXwYe6eRDPQfLzNURcR3wNxFxAtU05TlUR3B1ujupjkK7NCLeQfXH4N1U1yPvGBExD/gO8PHMvHhg+3iypd07fP8f1aFJZ1Mdn/q3wKmZ+bN29jVB3gpMo3ojXgO8JzN/0t6WWuZlwCuAPwIei4g1tdsP2tyXmnca8CuqEeNSqoHZx9vZUCvUzjN5NbA/1XlE3we+Dny6nX2NwX+nOjLy/EG/X2si4gLGmC2u6ilJBWr3nL8kqQ0Mf0kqkOEvSQUy/CWpQIa/JBXI8JekAhn+klSgdp/hK01qEXEq1QJnzwZ+QXUSYra3K2n8HPlLI4iI9wAXUS0HcDDwNNWyzlLH8wxfaRgRcRDV9QhenZn/Wtv2aqpVZ+dtLyuWqlyO/KXhvRtYOhD8NQ/X/rt7G/qRWsrwl+rUri19ItUiWYMNrJv++LbtSGo9p32kOhFxCPBzYD1DryfdRTXvP7vR1cykTuDRPtIzHUh1cZODGRr+lwIzDH5tDwx/6ZnmAI8PPqQzIqYBL6fDrmsrjcQ5f+mZHgZmRcTgS/+9nepazFe3pyWptRz5S890C9U1Xj8aEZcDrwI+DBxfuzKU1PHc4SsNIyJeA3wS2Au4nerSeLe1tyupdQx/SSqQc/6SVCDDX5IKZPhLUoEMf0kqkOEvSQUy/CWpQIa/JBXI8JekAv0njAPkB6K3hZcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f66654d61d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "skewed_chain = lognorm(1,8).rvs(10000)\n",
    "ax = plt.hist(skewed_chain,bins=100,normed=True)\n",
    "ax = plt.axhline(.025, c='k', linestyle='--')\n",
    "plt.xlabel('$\\\\theta$')\n",
    "plt.xlim((8,20))\n",
    "plt.title('Posterior')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:03:09.086059",
     "start_time": "2016-03-02T08:03:09.079366"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean is:  9.643799028459755\n",
      "The standard deviation is:  2.1150266705678\n",
      "The median is:  9.011681077018235\n",
      "The 95% Credible Interval is:  [ 8.15066093 15.03739996]\n"
     ]
    }
   ],
   "source": [
    "print(\"The mean is: \", skewed_chain.mean())\n",
    "print( \"The standard deviation is: \", skewed_chain.std())\n",
    "print(\"The median is: \", np.percentile(skewed_chain,50))\n",
    "print(\"The 95% Credible Interval is: \",np.percentile(skewed_chain,[2.5,97.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:03:09.942654",
     "start_time": "2016-03-02T08:03:09.934056"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The highest posterior density interval is: \n",
      "[ 8.04731833 13.1198713 ]\n",
      "Note: Min Chain Value is:  8.025551305915142\n"
     ]
    }
   ],
   "source": [
    "print(\"The highest posterior density interval is: \")\n",
    "print(pm3.hpd(skewed_chain))\n",
    "print(\"Note: Min Chain Value is: \",skewed_chain.min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:03:11.537711",
     "start_time": "2016-03-02T08:03:10.881613"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEgCAYAAABcnHNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGuZJREFUeJzt3X20XVV57/EvrwlBEXwhiAwMXvBBqEDqaKmgVoRrR5SEGsVKc6MovqASwauUl2uKeqUEClcBBdKSqK2AdZgIhjTG1kiLIr4Uc6MGHyJcpBGCSIqaHBJikvvHWgc3i5Octc/eyT77rO9njD2yM/fc88yV7PM788w111y7bN26FUlSs+za6w5IknY+w1+SGsjwl6QGMvwlqYEMf0lqIMNfkhpo9153QOqWiLgN+NOWoq3ABuAe4OrMnNelr3MEcFhm3tJhO7cBazLzLd3ol9QOR/4aa24Bnl8+DgReCiwGro+IN3bpaywGXt6FdqYD7+lCO1LbHPlrrNmQmWsqZf8rIt4E/A9gQRe+xi5daIPMXNuNdqSRMPzVFJsppoCIiJnA/wQC+C/gi8BHMvPx8vXXAp8A/gB4HPgG8MHM/EVE3A+8EDgvIt6SmZPK97wN+CvgvwEPADcAl2TmE+XrW8s2/xLYB5gGXELLtE85nXQx8ApgL+B24LzMXFG+/jlgX4rv2z8F/i4zP9T1fyk1gtM+GtMiYp+IuAB4CfBPEXEOMB/4R+Bo4P3AW4AvlfWfTTF19DXgSODPgBeX7wH4I2A18OnyORHxbuBq4FLgCOAc4K3A5yvdOYvit4+Tge9X+vlC4A5gj/JrvgJ4Arg9Il7UUvUU4HvAMcBnRvavIjny19jzxohYVz7flWIEvQb4EEWoPwT8fWb+n7LOqnJUfnNETKY4STy+fM/PM/P/RcSpwHMAMvORiNgMrM/MR8o2ZgOXZeY/lH+/r2zznyPigsy8vyy/KTO/M9jRiGjt9/sowv4vMnN9+fpbgHuBs8sHwHrgf2emm3KpI478NdYspRgVH0MxbbN/Zj6/DPv9gYnAv1fe82/ln0dn5nLgCxSj6kci4kvA8cAPh/piEfE84CDgIxGxbvABfLms8pKW6qu20++jgLsGgx+gnIb6HsVvKIPuNfjVDY78Ndasy8yftfmewUHQRoDMnBkRHwOmAK8BrgHOiohXZOaGbbz3PIpVQFUPtTx/vM1+Dba/scM2pKdx5K/GyMyHgYeBV1VeGrw24CcRcVREXAc8kJlXZ+YbgKnAyyhG51BMDQ36Zfk4NDN/NvigOCn8t8AzanZvBfCyiNh7sCAi9qI4r/CT2gcp1eTIX00zB/jbiFgF3EpxgvbTwNcyc0VEHEhxUnbPiLiMIujfSrEq6KdlG78FDouIF5QrgOYAl5UrgW6hWPFzPXB3+QOnjmso1vx/MSL+uvy6H6X44XFNZ4csPZ0jfzVKZn6KYtXNuylG1J+mWJb5xvL1B4HXAYcB3wX+A3gBcFJm/qZs5gqK6aAVEbF7Zn6SYtXQO4CVwGcpfgi8qY1+/ZziN5I9KJZ4/nv5/PgRTGNJw9rFO3lJUvM48pekBjL8JamBDH9JaqBRv9onIsZRLHd7iGJ/FknS8Haj2N32+5m5sfriqA9/iuC/vdedkKQ+9UrgW9XCfgj/hwBuuOEGDjjggF73RZL6wpo1a5gxYwY89SrzJ/VD+G8GOOCAAzjooIN63RdJ6jdDTpd7wleSGsjwl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8JamBDH9JaqB+uMhLkp5m0vm/v2Xy/XNeP2y5nsqRvyQ1kCN/SaOaI/kdw5G/JDWQ4S9JDWT4S1IDGf6S1ECGvyQ1kKt9JPWN1pU/6owjf0lqIMNfkhrI8JekBqo15x8Rk4G5wJHAKuDMzLyzUmcc8EngVGBP4DbgfZn5i/L1zwDvBDa1vO2IzHygw2OQ1HCeC2jfsCP/iBgPLAI+C+wLXAUsjIg9K1VnA0cAATwPeBS4uuX1Y4AZmfmMlofBL0k9UGfa5wRgS2Zem5mbMnM+RbBPrdS7CJiSmWuBicA+wK8AImJX4Chgedd6LkkasTrTPocDKytlSTEFtODJgszNwOMR8VHgr4EHgVeVLx8GTAAuj4jjgdXA7My8taPeS5JGpM7If29goFI2QBHmQ5lTvmcBsDQi9gD2ozgHcBlwIPAx4EsR8dIR9FmS1KE6I/8BYK9K2QRg3VCVM3MDQEScC7wXeGl5cvjElmo3R8Q3gJOBH7XbaUmqw+2gt63OyP9uipO4rYLKVFBEzI+I97YU7V62/1hEnBgR76m0MR7Y0GZ/JUldUGfkvwwYFxGzgOuAmRQndJdW6n0PODcilgC/BK4Ebs/M+yLiYOCKiFgJ3AG8GTgWOL0rRyFJasuwI//M3AhMAU4D1gKzgGmZuT4ilkTEhWXVucDngW8DP6eYGjq1bOM24GxgPvAb4Fxg6uA1AJKknavWRV6ZuQI4bojyKS3PtwIfLx9DtTEPmDeybkqSusntHSSpgQx/SWogw1+SGsjwl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8JamBDH9JaiDDX5IayPCXpAYy/CWpgQx/SWogw1+SGsjwl6QGMvwlqYFq3ckrIiZT3KbxSGAVcGZm3lmpMw74JMWtG/cEbgPeN3irxog4CfgUcAhwF3BGZt7TncOQJLVj2JF/RIwHFgGfBfYFrgIWRsSelaqzgSOAAJ4HPApcXbYxEVgIXADsB/wrcFN3DkGS1K460z4nAFsy89rM3JSZ8ymCfWql3kXAlMxcC0wE9gF+Vb42HViemYsy8wngE8CLIuJlXTkKSVJb6oT/4cDKSllSTAH9viBzc2Y+HhEfBX4O/AkwZ6g2MnMzcG+1DUnSzlEn/PcGBiplA8CEbdSfU75nAbA0IvYYQRuSpB2ozgnfAWCvStkEYN1QlTNzA0BEnAu8F3hpu21IUrdNOn/xk8/vn/P6HvZkdKgz8r+b4iRuq6AyFRQR8yPivS1Fu5ftP1ZtIyJ2Aw6ttiFJ2jnqjPyXAeMiYhZwHTCT4oTu0kq97wHnRsQS4JfAlcDtmXlfRHwFuDQipgO3Uqz6WQ38sDuHIUlqx7Aj/8zcCEwBTgPWArOAaZm5PiKWRMSFZdW5wOeBb1Oc8J1AseafzFwDnEKxIuhR4CRgemZu7e7hSJLqqHWRV2auAI4bonxKy/OtwMfLx1BtfBM4emTdlCR1U63wl6SdqfXkrHYM9/aRpAYy/CWpgQx/SWogw1+SGsjwl6QGMvwlqYEMf0lqIMNfkhrI8JekBvIKX9XmlrjS2OHIX5IayPCXpAYy/CWpgQx/SWogw1+SGsjwl6QGqrXUMyImU9ym8UhgFXBmZt45RL2PAO8G9gGWA2dl5o/L1z4DvBPY1PKWIzLzgY6OQJLUtmHDPyLGA4uAi4HrKW7gvjAiJmXmEy31TgfeCrwaeAA4H1gcEYdk5hbgGGBGZn652wehHcc7KkljU51pnxOALZl5bWZuysz5FDdhn1qp91zg4sy8LzN/B1wJHAwcFBG7AkdR/DYgSeqxOtM+hwMrK2VJMQW04MmCzMsrdaZR/JBYDRwGTAAuj4jjy7LZmXnrCPstSepAnfDfGxiolA1QhPmQIuJVwHXAezJzS0TsB9wGXAZ8H3g98KWIODYzfzSSjmvHcapHGvvqhP8AsFelbAKwbqjKETETuAaYlZk3ApQnh09sqXZzRHwDOBkw/CXtVO5TVS/87wbOqpQFcGO1YkTMBs4BTsnMZS3lJwKHZubclurjgQ1t91ijgt88Un+rE/7LgHERMYtiKmcmMBFY2lopIt4OfBA4LjN/WmljM3BFRKwE7gDeDBwLnN5R7yVJIzJs+GfmxoiYQhH8fwP8DJiWmesjYglwe2b+DXAB8EzgBxHR2sQfZeZtEXE2MB84kOKE8dTM/EV3D0eSVEeti7wycwVw3BDlU1qev3iYNuYB89rtoCSp+9zeQZIayPCXpAbyNo6SRgWvL9m5HPlLUgMZ/pLUQE77qKu8+EvqD478JamBDH9JaiDDX5IayPCXpAYy/CWpgVzto455cY7Ufxz5S1IDGf6S1ECGvyQ1kHP+Apy3l5rG8G+Ynbn9gls9SKOX0z6S1EC1Rv4RMRmYCxwJrALOzMw7h6j3EeDdwD7AcuCszPxx+dpJwKeAQ4C7gDMy855uHIQkqT3Dhn9EjAcWARcD1wMzgYURMSkzn2ipdzrwVuDVwAPA+cDiiDgEeB6wEJgBLKW42ftNwMu6eCzaBufzNVr52eydOtM+JwBbMvPazNyUmfOBR4GplXrPBS7OzPsy83fAlcDBwEHAdGB5Zi4qf2B8AnhRRBj+ktQDdcL/cGBlpSwppoB+X5B5eWZ+vqVoGsUPidXVNjJzM3BvtQ1J0s5RJ/z3BgYqZQPAhG29ISJeBVwHfCAzt4ykDUnSjlMn/AeAvSplE4B1Q1WOiJnAYmBWZt44kjYkSTtWndU+dwNnVcoCuLFaMSJmA+cAp2Tmskobp7bU2w04lKdPJ2kn2pkn21zzL40udcJ/GTAuImZRTOXMBCZSrNp5UkS8HfggcFxm/rTSxleASyNiOnArxWqf1cAPO+u+JGkkhp32ycyNwBTgNGAtMAuYlpnrI2JJRFxYVr0AeCbwg4hY1/J4SWauAU4BLqI4CXwSMD0zt+6AY5IkDaPWRV6ZuQI4bojyKS3PXzxMG98Ejm63g5Kk7nN7B0lqIDd2k9RoTV2M4MhfkhrI8JekBjL8JamBDH9JaiBP+EraqdzGeXRw5C9JDeTIf4xydCVpewx/7XRNXVctjSZO+0hSAxn+ktRAhr8kNZDhL0kNZPhLUgMZ/pLUQIa/JDVQrXX+ETEZmAscCawCzszMO7dT/0pgU2Z+uKXsM8A7gU0tVY/IzAdG0nFJ0sgNG/4RMR5YBFwMXE9xA/eFETEpM5+o1H0OcAXwtvLPVscAMzLzy93ouCRp5OpM+5wAbMnMazNzU2bOp7gJ+9Qh6n4L+B2woLUwInYFjgKWd9hfSVIX1Jn2ORxYWSlLiimgBZXyEzPzwYj4XKX8MGACcHlEHA+sBmZn5q3td1lDGQt7+bjtg7Tz1Bn57w0MVMoGKML8KTLzwW20sR9wG3AZcCDwMeBLEfHS2j2VJHVNnZH/ALBXpWwCsK7uFylPDp/YUnRzRHwDOBn4Ud12JEndUSf87wbOqpQFcGPdLxIRJwKHZubcluLxwIa6bWhsGgvTVVI/qhP+y4BxETELuI5itc9EYGkbX2czcEVErATuAN4MHAuc3lZvJUldMeycf2ZuBKYApwFrgVnAtMxcHxFLIuLCGm3cBpwNzAd+A5wLTM3MX3TQd0nSCNW6yCszVwDHDVE+ZRv1Tx+ibB4wr83+SdJOU52GHMurztzeQZIayPCXpAYy/CWpgbyBex9zmaSkkXLkL0kNZPhLUgMZ/pLUQIa/JDWQ4S9JDWT4S1IDGf6S1ECGvyQ1kOEvSQ1k+EtSAxn+ktRAhr8kNZDhL0kNVGtXz4iYDMwFjgRWAWdm5p3bqX8lsCkzP9xSdhLwKeAQ4C7gjMy8p4O+S+oT7kA7+gw78o+I8cAi4LPAvsBVwMKI2HOIus+JiM8BH6iUTwQWAhcA+wH/CtzUaeclSSNTZ9rnBGBLZl6bmZsycz7wKDB1iLrfAn4HLKiUTweWZ+aizHwC+ATwooh4WQd9l6QdatL5i598jDV1wv9wYGWlLCmmgKpOzMx3Auu210Zmbgbu3UYb0pj+ppNGgzpz/nsDA5WyAWBCtWJmPridNn5Tpw1tn2EoqRvqjPwHgL0qZRN4+uh+R7chSeqSOuF/NxCVsuDpU0G124iI3YBD22xDktQldaZ9lgHjImIWcB0wE5gILG3j63wFuDQipgO3Uqz6WQ38sL3uSuoXTlGObsOO/DNzIzAFOA1YC8wCpmXm+ohYEhEX1mhjDXAKcBHFSqGTgOmZubWTzkuSRqbWRV6ZuQI4bojyKduof/oQZd8Ejm6zf5KkHcDtHSSpgQx/SWogw1+SGqjWnL96y1UTkrrNkb8kNZDhL0kNZPhLUgMZ/pLUQIa/JDWQq31GKVf4SNqRHPlLUgMZ/pLUQE77aNRrnQK7f87re9gTaexw5C9JDeTIfxTxJK/6nZ/h/uHIX5IayJG/+orz/+qVsfbZc+QvSQ1Ua+QfEZOBucCRwCrgzMy8c4h65wDnAs8Evgq8JzPXl68tBl4DbB6sn5nP6PQAJEntGzb8I2I8sAi4GLgemAksjIhJmflES72TKYL/BOBh4Cbg48CHyirHAK/MzB909QjUWGPt13BpZ6oz7XMCsCUzr83MTZk5H3gUmFqpNxOYl5n3ZOavgdnAGRGxW0TsD+wP/LibnZckjUyd8D8cWFkpS4opoO3VS+BZwAuAycBvgVsj4pGI+HZEvHxkXZYkdapO+O8NDFTKBoAJw9QbfD4BGA98BzgbOAj4ArAkIg5ot8OSpM7VOeE7AOxVKZsArBum3uAPh3WZeQtwS8tr10bE+yimlG6q311paM7/S+2pM/K/G4hKWfD0qaBqvQB+DTwYEW+KiDdX6o8HNrTRV0lSl9QZ+S8DxkXELOA6ihO7E4GllXpfAK6LiAXAf1Ks9LkhM7dExDOAORHxY4qloudQ/Jbw9e4cRv/ycnhJvTDsyD8zNwJTgNOAtcAsYFpmro+IJRFxYVlvEXApsBh4AHiMYuknmfk54Erga2X5NGDK4DUAktRPJp2/+MlHv6p1kVdmrgCOG6J8SuXvVwFXbaONS4BLRtBHSVKXubdPD/TzaEGq8vPcn9zbR5IayPCXpAYy/CWpgZzz30mcF915tvdv7QVgUsGRvyQ1kOEvSQ3ktI8kdaBf95Vy5C9JDeTIv8v6dRQgtcMFDP3Pkb8kNZDhL0kN5LTPDuSvxpJGK8NfjeX5GXVbP32mDH81ir+NSQXDX6rop9GbNFKGfxc4mux/2/o/bPoPAj/bY1et8I+IycBc4EiKe/CemZl3DlHvHIpbNz4T+CrwnsFbNUbEacDFwP7AbcAZmflwF46hJ/ymkLQ9o33gMGz4R8R4YBFFcF9PcQP3hRExKTOfaKl3MkXwnwA8DNxEcRP3D0XEURQ3f38tsAK4GrgGeGNXj2YHM/A1aLR/Y2t0GY2flzoj/xOALZl5bfn3+RHxQWAqsKCl3kxgXmbeAxARs4FvRMRfATOAWzLzu+Vr5wEPR8T+mfnLLh1L2+r8hxj4GtSUz0JTjrNXRssPgjrhfziwslKWFFNACyr1vlKp8yzgBeVr33nyhcxHI+Kxsny48N8NYM2aNTW6OrRXXPrNYetMmvWPI25fzdb62fnWeSf0sCfbV+f7QDvX6tWrhyxv/b8a6WeqJTN3G+r1OuG/NzBQKRsAJgxTb/D5hDbaGMrzAWbMmFGj6tDGjfidUntO/Ponet2FbfL7YPTZ1udlXI06bXg+cG+1sE74DwB7VcomAOuGqTcY7OvaaGMo3wdeCTwEbK5RX5JUjPifT5GhT1Mn/O8GzqqUBXDjEPWiUufXwIPV1yLiucCzy/LtysyNwLdq9FOS9FRPG/EP2mXr1q3bfWdEjAPuA+ZQrNiZWT4/ZHAZZ1lvKr9f0fOfFKt97s/M90fEMcC/Aa8HfkCx2ufAzBwdp70lqWGG3dWzHHlPAU4D1gKzgGmZuT4ilkTEhWW9RcClwGLgAeAxiqWfZOZy4F3AfIoTvAcCb+/60UiSahl25C9JGnvcz1+SGsjwl6QGMvwlqYEMf0lqILd07lDdHU/7VUS8AriCYiuOXwGXZebc3vaq+yJiIvAj4B2ZeWuv+9MtEXEQxRLsVwG/ofj/u6q3veqOiDgOuAp4McVFoB/LzOr1R30pIv4YuDkzDyz/vh/FasnXUFw/9bHMnNfJ13Dk34GWHU8/C+xL8UFcGBF79rRjXVJ+4L5KcVz7AacCl0TEST3t2I4xD3hOrzvRTRGxC3AzxcWUzwH+DPhoGZp9LSJ2ozi2OZm5D/BO4PMRMamnHetQROwSEe8Avg605sjfU+yIMBF4E3BZuVvyiBn+nXlyx9PM3JSZ84FHKXY8HQteCCzOzBsyc0tm3gV8E+j78GgVEWcC6ykuThxLjqW4pub88vP5E+DlFJsu9rt9gecBu5c/5LYAT9D/W8BcCJxNsYU+ABHxDODPgYsyc0Nmfo9ih4V3dfKFDP/ObG/H076Xmcszc+bg38vfBF4J/N/e9aq7IuIw4EPAe3vdlx3gD4GfUIwS10TEPcCfZOajPe5Xx8pjuIZiJ4FNwO3AWZnZ7z/A5wPH8NT9eA4DNmXmfS1lHeeM4d+ZTnYr7SsR8SyKKa7/KP/sexGxO/AF4OzMXNvr/uwAz6b47fRXwMHA6cDVEfHKXnaqGyJiV4rvtVMpvt+mAp+KiKN72rEOZeZDmVm98nZv4PFKWcc5Y/h3ppPdSvtGRBwC3EGxvcf0zNzS4y51y2xgeWb+c687soNsBNZm5iWZ+URm3kFxD45TetyvbpgOHJuZXy6PbTHF1jJv7XG/doQdkjOGf2eqO5lS/r06FdS3IuIPge8CS4E/z8zqCKSf/QXwloh4rLy50MHAFyPi/B73q1sS2Lv8DWfQbsAuPepPNx3M029RsKl8jDWrgD0i4uCWso5zxqWenVkGjIuIWfx+x9OJFEHZ98rlj18DrsjMS3vdn27LzMNb/x4R91PMG4+VpZ7/AvwXMKf8gfbHwBuA/97TXnXHv1CsPHs78DmKpaxvoFgKOaZk5m8j4haK430XxVz/XwKv66RdR/4d2N6Opz3tWPecQbGiYnZErGt5XDzcG9V75W9prwb+gGI33RuBD4yF61Ay80cUSx7Pplj3/hngbZn5g552bMd5F7AHsJpi6u7cwXuij5S7ekpSAznyl6QGMvwlqYEMf0lqIMNfkhrI8JekBjL8JamBDH9JaiCv8JU6EBFvA84DXgQsp7jQaCxsmawxzpG/NEIRcS4wB/gwxTa8vwP+oaedkmryCl9pBCLiSIr7GrwuM79elr2OYmfJiZn5y172TxqOI39pZD4MrBgM/tIj5Z/P7UF/pLYY/lKbyvvHTqfYYKvV4J7rv965PZLa57SP1KaImAzcBWzgqfeM3Y1i3n+fIe7GJI0qrvaR2nc4xQ3Dj+Gp4X85MN7gVz8w/KX27Qv8unVJZ0TsARxPsexTGvWc85fa9wjF7RH3bCl7P8W9Vm/oTZek9jjyl9q3jOJesR+PiGuB1wIXUdzFbWNPeybV5AlfaQQi4hTgk8CBwPcpbqvX97dHVHMY/pLUQM75S1IDGf6S1ECGvyQ1kOEvSQ1k+EtSAxn+ktRAhr8kNZDhL0kN9P8BpAIYv3TLjmsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f66654b91d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# suppose our chain was this\n",
    "bimodal_chain = np.append(norm(2,.8).rvs(5000),norm(8,.5).rvs(5000))\n",
    "plt.hist(bimodal_chain,bins=100,normed=True)\n",
    "plt.xlabel('$\\\\theta$')\n",
    "plt.title('Posterior')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:03:19.871711",
     "start_time": "2016-03-02T08:03:19.865312"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The highest posterior density interval is: \n",
      "[0.75956428 8.85741972]\n"
     ]
    }
   ],
   "source": [
    "print(\"The highest posterior density interval is: \")\n",
    "print(pm3.hpd(bimodal_chain))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pymc `hpd` function is not built for multi-modal chains.\n",
    "\n",
    "The following shows what it should be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:08:41.384702",
     "start_time": "2016-03-02T08:08:40.755087"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAJnCAYAAACH/9UKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XlcVPX+x/EXq4iCSu47V+xklmaiqGF4NTNR8UaL27XbpkalltmVLE29mmvlmltldatbWt5U3DJzwVzSrLypHRNDA3MFVBQBgd8fg/NjAHXcmPH4fj4e87iH71nmM9Pne7+f+fqdMx55eXmIiIiIiMjFebo6ABERERGRG4EKZxERERERJ6hwFhERERFxggpnEREREREnqHAWEREREXGCCmcRERERESd4uzoAERGrMwxjLRBRoCkPOAvsAaaZpvneNXqe24H6pmkuusrrrAUOmabZ/VrEJSJiFZpxFhEpGYuAavmP6sCdwFLgXcMwHrpGz7EUaHkNrhMN9LsG1xERsRTNOIuIlIyzpmkeKtT2qmEYDwN/B768Bs/hcQ2ugWmaKdfiOiIiVqPCWUTEtXKwLdvAMIzewCDAAFKBz4DXTNPMyN9/PzAauAPIAFYDL5qmmWwYRiJQBxhiGEZ30zTr5p/zD+CfQD3gAPAJMNY0zaz8/Xn51+wJBAJRwFgKLNXIXwIyBggHSgPxwBDTNHfk7/8AKI9tTIkA5pim+dI1f6dERFxMSzVERFzAMIxAwzBeARoAnxuG8QLwPvBvoDHwHNAdmJ9/fBC25R4rgIZAB+DW/HMAmgFJwPT8bQzD6AtMA8YDtwMvAI8BHxYK53lss96dga2F4qwDbAR88p8zHMgC4g3D+EuBQ7sC3wN3ATOu7F0REXFvmnEWESkZDxmGkZ6/7Ylt5vYQ8BK2gvhPYK5pmm/lH/Nb/mzwV4ZhNMH2hUK//HP2m6b5u2EYjwC3AJimedQwjBzgtGmaR/OvMQyYYJrmR/l/78u/5jLDMF4xTTMxv/0/pmluOh+oYRgF434WW6HczTTN0/n7uwMJwMD8B8Bp4F+maeZd+VskIuLeNOMsIlIyVmKbjb0L21KLyqZpVssvlCsDVYD1hc5Zl/+/jU3T/An4GNts7lHDMOYD9wA/FvdkhmFUAmoCrxmGkX7+AXyRf0iDAof/dpG4GwHbzxfNAPlLR77HNjN+XoKKZhGxOs04i4iUjHTTNPde5jnnJzcyAUzT7G0YxkigI9AWeAd43jCMcNM0z17g3CHY7rZR2J8FtjMuM67z18+8ymuIiNxQNOMsIuJipmkeBg4D9xbadf7ezzsNw2hkGMYs4IBpmtNM03wQ6AI0xTYrDLblHOcdyX+EmKa59/wD2xcIJwJlnQxvB9DUMIwy5xsMwyiNbR31TqdfpIiIBWjGWUTEPYwDJhqG8RsQh+3LfNOBFaZp7jAMozq2L/D5GoYxAVuR/Bi2u2/8mn+NU0B9wzBq5N9pYxwwIf+OG4uw3VnjXWB3frHujHew3dP5M8Mwhuc/7whshfc7V/eSRURuLJpxFhFxA6ZpTsZ2d4u+2GZyp2O7ddxD+fsPApFAfWAL8ANQA7jPNM2T+Zd5E9sSjh2GYXibpvk2trtzPAnsAuZhK6Afvoy49mObCffBdhu69fnb91zB0hMRkRuaR16evsshIiIiInIpmnEWEREREXGCCmcRERERESeocBYRERERcYIKZxERERERJ7j97egMwyiF7X6hfwI5Lg5HRERERKzLC6gGbDVNM7PwTrcvnLEVzfGuDkJEREREbhqtgQ2FG2+EwvlPgE8++YSqVau6OhbLO3XqFIsWLQKga9euBAQEuCyWzJxMzGMmAEZFg1JepVwWi0hJc6e+eE0dOAAR+T+IuG4d1K7t2nhuMpbNKwHcY9y80XPs0KFD9OrVC/Lrz8JuhMI5B6Bq1arUrFnT1bHcFBo0aODqEOzq1ann6hBEXMad+uI1k5kJ587ZtqtUAf3/eomzZF6JnTuMmxbJsWKXB+vLgSIiIiIiTlDhLCIiIiLiBBXO4mDfvn0EBAQQEBDAvn37XBpLakYqb8S/wRvxb5CakerSWERKmjv1RbEO5ZW1ucO4afUcuxHWOEsJysvLIz093b7t0ljIIysny74tcjNxp74o1qG8sjZ3GDetnmOacRYRERERcYIKZxERERERJ6hwFhERERFxggpnEREREREn6MuB4tby8vI4k3LG1WGIiIiIaMa5JG3atIm//e1vNGnShG7duvHzzz/b9+3YsYMGDRrQpEkT+2PWrFkA7N69m44dOxIaGsq4cePs55w+fZqoqCj7t1eLs3DhQqKjo4u0r1mzhrZt2wKQlJSEYRg0adKE6OhoQkJCqFevHqNHj+bw4cP26xSM76677iIqKooFCxZc9DVnZWUxY8YMHnjgAZo0aUJERARjxozh9OnTTr1nP335E3vW7HHq2CVLlhAWFkazZs1ITk5mwIABNG7cmJiYmCLHvvfeezRt2pR77rmH7Oxsp65fnC1bttjfu/PvS5cuXVizZo39mOHDh/P2228Xe/7ChQsxDIOXX365yL6vv/4awzCYNm3aZce1Z88eDMNw6ti2bds6xFvYkiVL6NatG6GhobRs2ZIBAwaQmJjo9PkXsnjx4vM/ayoiInJD0IxzCUlKSiImJoahQ4cSHR3Nhg0b6Nu3L3FxcVSqVIlff/2Ve++9l9mzZxc5d+bMmfTs2ZOuXbvy4IMP8vDDDxMSEsKcOXPo0aMHZcuWvSYxbtiwgaysLGbPns25c+c4dOgQAwYM4PPPPwdsP6G5cOFCAHJzc9m8eTODBg0iOzubnj17FrneuXPneOqpp/D392fmzJkEBweTnJzMsGHDiImJ4aOPPrpoPKW9SxNEEGXLlaW0d+lLxr9w4UJ69uzJwIEDOXz4MCtXruSbb76hVq1aRY79/PPPeeWVV3j44YedeWsuqnz58mzZsgWwzZAvX76cAQMGsHr1aipXrsyoUaMuef7q1as5e/Ysfn5+9vYlS5ZQpkyZq47varz99tssX76cN954gyZNmpCRkcGMGTPo1asXixcv5pZbbrnia0dFRREVFXUNo7WWoKAgxo4da98WuRaUV9ZW2rs09/3lPvu2K1g9xyw545yQkEBCQgJnzpwptr3wbOfvv/9OQkJCkZnbxMREEhISOHXqVLHtl2P9+vXceuutPProo3h7e9OmTRsaNWrEihUrANi1axe33XZbsed6e9s+3+Tl5ZGXl4eXlxcHDx4kPj6eRx999LLiuJQKFSoQGxvLa6+9xsMPP8yePcXP9np6etKqVSuGDBnCtGnTyM3NLXJMXFwc+/fvZ+rUqQQHBwNQo0YNJk2aRGBgIEePHrXPdhf8bxIdHc3ChQv57OPP+H7N96xZtIYhLw0B4LvvviM6Opq7776brl27sm7dOgCefPJJNm/ezNy5c3nmmWfo0KEDYCvOli1b5hBXhw4dOHDgAKNGjWLUqFGcO3eOyZMnc++99xIWFsaAAQMcZtp79uzJI488QlhYGPv377/o++fh4UFkZCS+vr72HImNjWX8+PEXPKdGjRrUrVvXYdb21KlT/PjjjzRv3tzedvr0aUaOHMk999zDPffcw6uvvmrPzdzcXN566y3CwsIIDw9n6dKlDs+xdOlSoqOjadasGc2bN2f48OGXvL9mcnIyc+bMYdq0aYSGhuLl5UXZsmUZMmQIbdq0cegDmzdvJioqiiZNmtC3b19OnDgBQGpqKi+99BJt27alcePGdOnShR9++MH+3p7/15Bp06YxePBg+vXrR5MmTYiMjGTDhg0Xjc/qzvfF2NhYKlSo4OpwxCKUV9ZW2qc04bXDCa8dTmkf1xTOVs8xSxbOISEhhISEFBl4b7/9dkJCQli9erVD+1133UVISAjLly93aG/evDkhISEsWrTIob1169aEhIRcVky5ubkOs4lgKz7PF2K7d+9m+/bttG3bljZt2jB+/Hiysmw3MY+JieHLL7+kY8eOREdHExwczJtvvsnAgQPx8vK65HP/+uuvhIaGOjwGDRp00XOOHDnCZ599RlhY2EWPa926NSkpKfz+++9F9sXHxxMREUGpUqUc2oOCgpg+fTqVKlW66LWfeOIJunTpQu/evZk6dSq//fYbMTExPPPMM3z//fcMGjSIgQMHYpom77//PqGhocTGxjJr1izi4uIA2yx6ZGSkw3VXrlxJ9erVmTJlCsOHD2fq1KmsXr2aTz/9lLVr1xIYGMjAgQPtheX27dsZNGgQ33zzDXXq1LlozLm5ucTFxeHn58edd9550WMLioqKcih2V6xYQdu2bfH19bW3DR8+nH379rFkyRKWLVvGsWPHGD58OAD/+c9/WLlyJQsXLmTp0qX89NNP9vOSkpJ47bXXGDFiBFu3buXTTz8lLi6OzZs3XzSmDRs2ULt27WKXfIwZM8ahqN+0aRPz5s1jzZo1JCUl8emnnwIwceJEAJYtW8bWrVtp2rQpb775ZrHPt2LFCh5//HG2bNlCREQE//rXvy71tomIiJQoLdUoIeHh4UyaNInly5dz3333sWnTJjZv3kzlypUB2ye0sLAwunXrxvHjxxk4cCBTp05l8ODB1K9fn6+++sp+rZ9//pm0tDTuuOMO+vbty7Fjx+jduzcPPvhgsc9922232ZdYnLdmzZoihUlERARgm9n29/enefPmvPLKKxd9XeXKlQMoMisPkJaWRs2aNS/xzjhv6dKltGrVivvvv98eb9u2bVmyZInT63mLs2jRIoYOHWqPdejQoTRr1sz+U6GVKlWiZcuWFzz/xIkThIaGAnD27Fmys7N59tlnL2sJTWRkJJMnTyY9PZ2yZcuyePFiXnjhBT788EP7dVeuXMlnn31m/6evIUOGEBkZyRtvvMGyZcvo1asXNWrUAGDAgAH2wrhy5cosWbKEmjVrkpqaSlpaGuXKlbPPql9IWlqa07MFTz75pH3ZRqtWrUhKSgLgxRdfpHTp0nh5eZGcnExgYOAFn/euu+6yv89dunRh3rx5Tj23iIhISbFk4bx3714AqlWr5tC+a9cuAKpWrerQ/tNPP5Gbm0uVKlUc2r///ntycnLsxe158fHx5OTkXFZMdevWZfLkybz11luMGDGCiIgIunTpQkBAAID9i4AA/v7+9OvXj7feeovBgwcXudaECRMYPnw4s2fPpnnz5nTv3p1OnToRERFxVeuJ1q1bx9GjRwkPDwfg3//+NxUrVrzoOampqUDR9xqgYsWKHDt2rNjzjh8/fsn1sWln09h5ZCd+GX6knU0jJSWF6tWrOxxTvXp1Dh06dNHrXMrx48cdruvv70/58uXtBd6lZsbLlStnX+MMtjwbMGAAgYGBPPHEE/b2bdu20adPH/vfc+fOtW9XrlyZxo0bs2rVKlq2bMmhQ4do2rSpvXA+efIk2dnZDnHWqFGDvLw8jhw5wrFjxxzyt+AHFm9vbxYsWMAXX3yBv78/t99+O9nZ2cUurymoYsWKHD9+vNh9KSkplC9fHk9PT/t7cJ6Pj4992dORI0cYM2YMCQkJBAcHU758+QsuESmYu97e3pb8qdbLkZiYaO+LGzZsoG7duq4NSCxBeWVtaWfTeP/H9wF4ssmTlPcrX+IxWD3HLLlUo169etSrVw9/f/9i2wt/4So4OJh69eoVmSGsW7cu9erVsxe3hdsvR3p6OtWqVWPx4sVs2bKFCRMmsGfPHm6//XZOnDjB+PHjHdZYZ2ZmFlniALZ/8v7LX/6CYRgkJCTQoEEDypYtS7Vq1eyzfFcjJyeH5ORkkpOTnfpwEB8fT6VKlYp8uADbMo7169dz9uxZh/aUlBQiIiLYvHmzfalJwTtbpKWlAZCbl0t2bjZZOVnk5uVSrVo1kpOTHa6VlJR0yeL+UqpXr+5w3dOnT5OamnrFX3y7/fbbad++PRs3bnRoDw0N5ccff7Q/zs9Sn9elSxfi4uKIi4ujc+fODvsqVqyIr6+vQ5xJSUl4enpSoUIFKleuzMGDB+37Cs7qLl26lGXLlvHVV1+xatUqpkyZUmxuFXbPPfeQnJzMr7/+6tCel5fH008/zfTp0y95jUGDBnHfffexefNmPvvsMx544IFLniM2l9sXRZyhvLK23LxcTmae5GTmSXLzLj45cr1YPccsWTi7o7S0NLp3787OnTvJysrik08+4eDBg7Rt25aAgABWrVrF9OnTyc7OZv/+/cyaNavIbeSysrKYOXMmAwYMAKBWrVr8/PPPnDhxgj/++KPYWd/rJScnh3Xr1vHWW2/xwgsv4OHhUeSYjh07UqNGDV544QUOHDgA2L6g+fzzz9O0aVPCwsK45ZZbCAgIYNGiReTk5PDf//7XoQD09PYkO8NWVEdGRrJlyxa+/vpr+/N/++23RdYwX66//e1vzJgxg+TkZDIyMhg7diwhISHceuutV3S9AwcOsHr1apo0aXJZ53Xo0IHt27ezYMGCIneb8PT0JCoqikmTJpGSksKJEyeYMGECERERBAYGEhUVxYcffsi+fftIT09n6tSp9nPT09Px9vbG19eXrKws5s6dS1JSEufOnbtoPFWrVuWJJ55g4MCBbNu2jdzcXFJSUhgxYgRHjx6le/ful3xN6enplC5dGg8PDxISEpg7d+5V3f5PRETElSy5VMMd1axZkxEjRtC/f3/S0tJo2LAh8+bNs8+Kz5o1i9GjR9OiRQv8/Pzo1q0b//jHPxyu8cEHHxAZGWlfOtCnTx/69+/Phx9+SExMzCWXFFyt3bt324tBHx8fateuzdChQ+nUqVOxx3t5efHee+8xZcoUHn/8cVJTU6lQoQIdOnTg+eefx8PDA19fX15//XWmTZvGlClT6NChA23atLFfo9bdtdg4dyMvxLzARx98xIwZM5g0aRJDhgyhRo0avPnmmzRq1OiqXlefPn3IzMykZ8+epKenExYWxpw5c4r9MFCctLQ0hyK5bNmydO7cmX79+l1WHAEBAYSHh/Pnn3/a70JS0CuvvMLEiROJiooiMzOTdu3aMXToUAAefvhhjh49Sq9evcjLy6NHjx7Ex8cD8OCDD7Jp0yb++te/4ufnR7NmzWjfvr1Td4Z5+eWXqVq1KiNHjiQ5ORk/Pz+aN2/OJ598Uuy/MhQ2atQoxo4dy8SJE6lSpQoPPfQQkydPti/xERERuZF4uPs6QsMw6gK/r169+pp+0UyKl5CQYL9jyN69ey97Scq1lJKRwtQttpnTAWEDCCptvftBilyIO/XFayohAc7flWjvXrDK67pBWDavBHCPcfNGz7GkpCTatWsHEGyaZmLh/VqqISIiIiLiBBXOIiIiIiJOUOEsIiIiIuIEfTlQHJQvX57Y2Fj7tiv5efsRXjvcvi1yM3GnvijWobyyNncYN62eY/pyoIiIlBx9OVBcoG7s0mLbE8cVf1couXnpy4EiIiIiIteACmcRERERESeocBYHBw4coGHDhjRs2ND+a3+ucuLsCWZ8P4MZ38/gxNkTLo1FpKS5U18U61BeWZs7jJtWzzF9OVAcZGdns2vXLvu2K+Xk5XD0zFH7tsjNxJ36oliH8sra3GHctHqOaca5BBmGQePGjUlPT3doz87OJiwsjLZt27ooMhERERG5FBXOJczPz4/Vq1c7tMXHx1vyU5mIiIiIlVhvqUZWFvzxR8k8V61a4Ot7Wad06NCBpUuX0rVrV3vbkiVLuP/++/n+++/tbVu3bmXcuHHs37+f4OBghg0bRqNGjQDYtGkTU6ZMITExkaysLMLDwxk/fjylS5emd+/e3H333axdu5Y//viD22+/nXHjxulWfiIiIiJXyVqFc1YWGAYkJpbM89WtC6Z5WcVzZGQk/fr1IzU1lQoVKpCens7WrVsZNmyYvXA+ePAg/fr1Y8KECbRp04ZVq1bRp08fVq5cia+vL88//zwTJkygXbt2HDp0iJ49exIXF8cjjzwCwNKlS5k3bx7ly5fnmWeeYc6cOYwaNep6vAMiIiIiNw0t1ShhQUFBNGvWjK+//hqAVatW0aZNG3wLFN9xcXGEhYVx33334e3tTceOHbn11ltZuXIlpUqV4r///S/t2rXj1KlTHDlyhPLly3P48GH7+VFRUdSqVYuAgADat29PYkl9kBARERGxMGvNOPv62maA3XipBkDnzp358ssv6datG0uWLOGZZ57h9OnT9v0HDx4kPj6e0NBQe9u5c+do2rQpXl5efPvtt3z44YeA7QuHGRkZFPwFyKCgIPu2t7c37v7rkCIiIiI3AmsVzmArZN38J1zbt2/PyJEj2blzJwcOHKBZs2asXbvWvr9SpUpERkYyYcIEe9sff/xBhQoV2L59OzNmzGDBggXUrVsXgMcee+yaxRYYGEhMTIx925VKeZWiWfVm9m2Rm4k79UWxDuWVc27Un+h2h3HT6jnmVOFsGEYTYDbQEPgNeMY0zc0XOf5JYIJpmhULtPUAxgCVgbXAU6ZpHi7+CtZWpkwZ2rRpwz//+U8iIyPx8PBw2N+pUyceeeQRNm3aRIsWLdi+fTtPP/0077zzDtnZ2Xh6euLn50dOTg5Llixh27ZtNGnS5JrEVqlSJd55551rcq2rVca3DJ1ude//kxK5XtypL4p1KK+szR3GTavn2CULZ8Mw/IAl2Ired4HewELDMOqapplVzPF/Ad4CzhVoawTMAu4HdgDTgHeAh67Ba7ghdenShZiYGKZMmVJkX926dZk8eTITJ04kMTGRoKAgXnnlFVq2bElubi4PPPAAXbp0wdPTkzvuuIMHH3yQhIQEF7wKERERkZuHx6XWvxqG0RGYbZpm7QJt/wNGmKb5ZaFjvYB1wEbgyfMzzoZhjAeqmab5WP7ftwCHgeqmaR65xPPXBX5fvXq1bqkmInKjS0iAkBDb9t69br+0Tqzhcpde3KhLNeTqJSUl0a5dO4Bg0zQTC+935q4atwG7CrWZ2JZtFBYL7ASWXewapmkeB9Ly28WNJCUl0apVK1q1akVSUpJLYzmZeZL3tr/He9vf42TmSZfGIlLS3KkvinUor6zNHcZNq+eYM2ucywBnCrWdAfwLNhiG0RTbMo7Q/MdlX0NcLzMzk02bNtm3Xelc7jn+OPmHfVvkZuJOfVGsQ3llbe4wblo9x5yZcT4DlC7U5g+kn//DMIzSwIfA06ZpplPUJa8hIiIiIuLOnJlx3g08X6jNAD4t8Hco8BcgzjCM89f1NwwjDWiUfw3DfrJhVASC8ttFRERERNyeM4Xzt0ApwzD6Y7szRm+gCrDy/AGmacZTYNmFYRhtgC8KfDnwP8A6wzDeB7YBY4Hl+WudRURERETc3iWXapimmQl0BHoAKUB/IMo0zdOGYSw3DGOoE9f4CegDvA8cAaoDT1xN4CIiIiIiJcmpH0AxTXMH0KqY9o4XOH4tULFQ23xg/uWHKCIiIiLies58OVBERERE5Kbn1IyzXBubNm1i/Pjx7N+/n1tvvZWhQ4fSuHFjAHbs2EG3bt3w8/OzH9+vXz+eeeYZdu/ezaBBgzh69CgPP/wwsbGxAJw+fZoePXrw6aefUrZs2Qs+b1ZWFnPnzmXJkiUcPnyYwMBA7r//fl544QXKlCkDQO/evenQoQPt27enV69eABe9Zknw9fKlUZVG9m2Rm0nZsmXdpi+KdSivrM0dxk2r55gK5xKSlJRETEwMQ4cOJTo6mg0bNtC3b1/i4uKoVKkSv/76K/feey+zZ88ucu7MmTPp2bMnXbt25cEHH+Thhx8mJCSEOXPm0KNHj4sm5rlz53jqqafw9/dn5syZBAcHk5yczLBhw4iJieGjjz5yOL5KlSp8/PHH1/z1X4myvmWJbhDt6jBEXMKd+qJYh/LK2txh3LR6jmmpRglZv349t956K48++ije3t60adOGRo0asWLFCgB27drFbbcV/0OK3t62zzd5eXnk5eXh5eXFwYMHiY+P59FHH73o88bFxbF//36mTp1KcHAwADVq1GDSpEkEBgZy9OjRa/gqRURERKzLkjPOKRkpF93v7+OPn/f/L4lIzUglj7wLHl/auzSlff7/91vSzqaRm5dLUOkgp2PKzc11WIYB4Onpyf79+wHYvXs3vr6+tG3bltzcXDp27MiLL76Ir68vMTExvPzyy/aZ5+DgYF566SUGDhyIl5fXRZ83Pj6eiIgISpUq5dAeFBTE9OnTnY5fRETkZlc3dmmx7YnjOpVwJOIqliycp26ZetH9kfUjaV6juf3vWdtmkZlz4Z+FbBfcjtZ1Wtv/fv/H9zmZeZIRbUY4HVN4eDiTJk1i+fLl3HfffWzatInNmzdTuXJlACpUqEBYWBjdunXj+PHjDBw4kKlTpzJ48GDq16/PV199Zb/Wzz//TFpaGnfccQd9+/bl2LFj9O7dmwcffLDI86alpVGzZk2n4zx48CBPPvmk7XW+/z7Vq1d3+txr7VTmKRaZiwDoanQloFSAy2IRKWnu1BfFOpRX1uYO46bVc8yShbM7qlu3LpMnT+att95ixIgRRERE0KVLFwICbEk9a9Ys+7H+/v7069ePt956i8GDBxe51oQJExg+fDizZ8+mefPmdO/enU6dOhEREUFQkOMseMWKFTl27FixMR0/fpxbbrnFoS0jI4OVK1fat10pOzebvSl77dsiNxN36otiHcora3OHcdPqOWbJwnlA2ICL7vf38Xf4+5nQZy65VKOgJ5s8SW5e7mXFlJ6eTrVq1Vi8eLG97dFHH+Wxxx7jxIkTzJo1i+eee87+Rb/MzMwiyysAli1bxl/+8hcMwyAhIYGIiAjKli1LtWrVSEpKKlI4t27dmvHjx3P27FmHpSIpKSlERETw7rvv0qJFi8t6LSIiIiI3I0t+OTCodNBFHwXXNwNUKF3hoscXXN8MUN6v/GWtbwbbkonu3buzc+dOsrKy+OSTTzh48CBt27YlICCAVatWMX36dLKzs9m/fz+zZs0iOtrxm7FZWVnMnDmTAQNsHwxq1arFzz//zIkTJ/jjjz+oVq1akeft2LEjNWrU4IUXXuDAgQMAJCQk8Pzzz9O0aVPCwsIu63WIiIiI3KwsWTi7o5o1azJixAj69+9PixYtWLFiBfPmzcPf3x9PT09mzZo9yW8MAAAgAElEQVTFr7/+SosWLejZsycPPPAA//jHPxyu8cEHHxAZGUmlSpUA6NOnD9988w33338/ffr0sbcX5OXlxXvvvUft2rV5/PHHadKkCX369KFx48a88847eHh4lMjrFxEREbnRWXKphrvq2rUrXbt2LXZfSEgIH3zwwUXP79u3r8PfNWrUYOHChZd83jJlyjB06FCGDh16wWP+/e9/A7bZaBEREREpSjPOIiIiIiJOUOEsIiIiIuIELdUQB/7+/nTu3Nm+7Uo+nj7cesut9m2Rm4k79UWxDuWVtbnDuGn1HFPhLA6qVavGkiVLXB0GAAGlAuh5Z09XhyHiEu7UF8U6lFfW5g7jptVzTEs1REREREScoMJZRERERMQJKpzFwaFDh+jWrRvdunXj0KFDLo0lPSudBTsXsGDnAtKz0l0ai0hJc6e+KNahvLI2dxg3rZ5jKpzFwenTp5k/fz7z58/n9OnTLo0lKyeLnUd3svPoTrJyslwai0hJc6e+KNahvLI2dxg3rZ5jKpzlujt37pwlP3WKiIjIzUWFs4u8/PLL3HHHHRw+fNjVoTiIjo7G19f3ml5z0KBBfPPNN1d0rrnaZPHQxbQPb0///v05duyYfd+7777LHXfcQZMmTeyPbdu2AbBo0SJatmxJeHg4cXFx9nN27NhBnz59LvqcsbGxjB8/vkj7+PHjiY2NBWDhwoU0aNDA4bnDw8MZM2YM2dnZ9usUjK9p06b07t3bHuOFHD58mGHDhnHvvfdy9913ExkZyccff+xwjGEY7Nmz56LXERERkWtLhbMLnDhxgnXr1tGhQwc+++wzV4fjYOHChWRlXdt/3klNTb2i875Z+Q2/xP1Cq6dasXzNckJCQoiJibHv3717Ny+++CI//vij/REaGgrA2LFjmTdvHu+++y5jxoyxnzNx4kT++c9/Xt0LytegQQOH5/7yyy/ZsGEDU6dOtR/Tu3dv+/7vvvuOBx54gKeffpqdO3cWe83Dhw8THR1NuXLl+Oqrr/jhhx8YO3Ys7733HtOnT78mcYuIiMiVseR9nBMSEi66v2LFipQrV87+9++//05ubu4Fjw8KCqJChQr2vxMTE8nJyaFevXpXFN9XX31FaGgovXr1on///sTExODr68u0adPYt28fx48fZ8eOHdSvX59Ro0bRoEEDtmzZwqhRo2jWrBmLFi0iKCiIQYMG0alTJ8A2A9mzZ0/i4uJ4+umn+fvf/86kSZP4+uuvAWjTpg2xsbGULVuWxx57jMqVK/Pmm2+SlZVFdHQ0999/PwMGDCAyMhJfX19yc3Pp3r07zz77LHPmzCE3N5fBgweTkpLCvHnz8PLyIjY2li5dugDw0UcfsWDBAg4ePEipUqXo0aMH/fv3Z8yYMWzbto0ff/yRpKQkYmNj2bp1K+PGjWP//v0EBwczbNgwGjVqVOR9Wrt6LfVa16NivYp4+3jTv39/PvjgA0zTxDAMdu/ezUMPPVTse+zt/f+p7eXlBcDy5csJDg6mfv36V/Tf7VKqVKlCmzZtLjgT7OfnR69evfjll1+YNWsW06ZNK3LM5MmTadq0KYMHD7a3NW7cmDFjxrBixYrrEreIyM2qbuxSV4cgNxhLzjiHhIRc9FH4n73vuuuuix4/c+ZMh+Nbt25NSEjIFce3YMECHnroIe6++26CgoIcCqLly5fTvXt3tm7dSkREBM8++6x9Bnjv3r34+PiwZcsWRo4cSWxsLL/99pv93MzMTL777jt69erF8OHD2bdvH0uWLGHZsmUcO3aM4cOH4+HhwdixY1mzZg3r1q3j7bffxt/fn2effbZInCdPniQ5OZn169fz0ksv8frrr5OSkkJ8fDzPPfcco0ePBmDbtm32QvCHH35g6tSpzJgxg/379/Pqq68SGhpKbGwssbGxHDx4kH79+hETE8PmzZt58skn6dOnD2lpaUWePzc3F2/f/y+APTw88PDwYP/+/WRkZJCYmMhHH33EPffcQ8eOHfniiy/sx7766qvExMTw3HPPMXLkSLKyspg1axYDBgxw6r/Rxx9/TGhoqMOjcN4UjnXPnj2sWrWKsLCwi167devW/PDDD8Xui4+Pp3379kXaW7VqxahRo5yKXURERK4PSxbO7mz79u2cPHmSNm3aANC9e3c++eQT+/6WLVsSGRmJj48PMTExnDlzhu3btwO2n64cPHgwvr6+hIeH07p1a5YvX24/t1OnTvj6+uLt7c3KlSt5+eWXCQoKoly5cgwZMoTly5eTkZFBzZo1GTJkCK+++ioLFixg4sSJDjO0BT3xxBP4+PjQokULcnJy7H+3bt2atLQ0MjIyaNiwIQsXLqRu3bocO3aM7Oxs/Pz8OHLkSJHrxcXFERYWxn333Ye3tzcdO3bk1ltvZeXKlUWODY8IZ2/8XlL/SCU7O5sZM2aQmZlJZmYmx44d4+6776ZHjx6sWbOGf/3rX4wbN45169YB0LFjR9asWcPq1atp164dH374IZGRkSQnJ/PII4/Qo0cP/ve//13wv9Pf//53tm3b5vD4+9//7nDMr7/+ai+qmzVrxvPPP0/Hjh35xz/+ccHrApQrV45Tp04Vuy8tLY2goKCLni8iIiKuYcmlGnv37r3o/ooVKzr8/dNPP11yqUZB8fHx5OTkXFFs8+fPJzU1lXvvvRew3XEiLS2NX375BYDatWvbj/Xy8qJSpUocO3aMSpUqUbVqVUqVKmXfX7VqVYcvy51/XSdPniQ7O5vq1avb99WoUYO8vDyOHDlCnTp16Ny5M+PHj+euu+6iTp06DjE2bdoUDw8Pjh07Zl/S4ulp+4wVEBAA2GZ/IX9W2Nubd955h5UrV3LLLbdwxx132PcVdvDgQeLj4+1rkc+/B02bNi1ybFTXKNbvWs/GWRt5aM5DPNb7MerVq0dgYCC1atVymAEODQ2la9eurF69moiICIfrpKSksHjxYr744gu6devGqFGjyMvL4/XXX2fhwoVFntdZt9122xWdn5qa6vDfpqDz/70Ly8nJ4dSpU5QvX/6yn09uTH5+frRu3dq+LXItKK+szdvTmzrl6ti3XcHqOWbJwvly1x4HBwdf1vF169a9rOPPO3XqFCtWrOCDDz5wKJDHjBnDxx9/TI0aNRxmac+dO8eRI0eoWrUqOTk5HD9+nJycHPua3YMHDzqsDT5fzFasWBFfX1+Sk5PtRX9SUhKenp72tdoTJ07krrvuYs+ePSxbtozIyEj7dT744AP8/f1p166d/ZoXM2/ePPbs2cM333xDQEAA2dnZLFu2rNhjK1WqRGRkJBMmTLC3/fHHHw5ryM87e+Isw58eTs0RNQHbB4IZM2bQoEEDdu7cyXfffUffvn3tx2dmZhbbSadMmULfvn0pVaoUCQkJNGjQgLy8vEuuhb9e4uPjadiwYbH7WrduzapVq+jatatD+9q1axk8eDDx8fGULVu2JMIUF6tRowbr1693dRhiMcorawssFcgTTZ5waQxWzzEt1ShBixYtonbt2jRt2pRKlSrZHw8//DBLly4lNTWV+Ph4Nm7caF+aUKFCBZo0aQLY7sYxZ84csrOzWbduHZs3b7Z/ObAgT09PoqKimDRpEikpKZw4cYIJEyYQERFBYGAgGzduZPHixYwZM4bXXnuNkSNHFjvL6az09HR8fHzw8fHh9OnTjB8/nuzsbM6dOweAr68v6em2XzDq1KkTa9asYdOmTeTl5fHDDz8QFRVV7LKJjRs30q9fP1JSUkhPT2f06NG0atWKypUr4+/vz/Tp01mxYgW5ubls2rSJpUuX8uCDDzpc47fffsM0TTp37gxArVq1+Pnnn9mxYwe1atW64td8JTIyMvjoo4/45ptv6NevX7HHPPvss2zdupW33nqLtLQ0cnJy2Lx5M6+//jpPPfWUimYREREXsuSMs7uaP3++vYArqFWrVlSoUIH58+fTqFEj5s6dy/PPP0/Dhg2ZPXu2fYY5MDCQQ4cOER4ezi233MKUKVOKLLM475VXXmHixIlERUWRmZlJu3btGDp0KOnp6QwdOpQXX3yRatWqUa1aNRYtWsSwYcOKfAnSWU888QSDBw+mZcuWlClThrZt23L33XeTkJDAPffcQ+fOnRk1ahTJycn861//YvLkyUycOJHExESCgoJ45ZVXaNmyZZHrdu3aFdM0iYyMJDc3l4iICPtMdXBwMJMnT+btt98mNjaWKlWqMHbs2CIzuePHj2fw4MH2mfNXX32Vl156CS8vL8aOHXtFr/dy/Pvf/7bfctDPz48777yTDz74AMMwij2+atWqfP7557z99ttERkaSkZFBjRo1ePbZZ+nZs+d1j1dEREQuzCMvL8/VMVyUYRh1gd9Xr15NzZo1XR3OdTVt2jR+++03h/sAn7dlyxYGDBjAli1brmsMR44c4bXXXgNg9OjRVK5c+bo+38WczjrNt79/C0Db4LaU8S3jslhESpo79cVrKiEBzt+VaO9euMLbesqVsWxeXcL1vu1c4rii//rrCu4wbt7oOZaUlES7du0Agk3TTCy8X0s1xMGpU6eYO3cuc+fOveCdH0pKZk4mP/z5Az/8+QOZOZkujUWkpLlTXxTrUF5ZmzuMm1bPMRXOIiIiIiJO0BpnN9K/f/8L7gsLC7vuyzRERERE5MI04ywiIiIi4gQVziIiIiIiTlDhLCIiIiLiBBXOIiIiIiJO0JcDxYGvr6/9Z7x9fX1dGouXhxdVylSxb4vcTNypL4p1KK+szR3GTavnmApncXD+J6ndQTm/csQ0i3F1GCIu4U59UaxDeWVt7jBuWj3HtFRDRERERMQJKpxFRERERJygpRri4NixY4wbNw6A2NhYKlas6LJYzmSfYcOBDQCE1w7H38ffZbGIlDR36otiHcqr66Nu7NJi2xPHdSrRONxh3LR6jqlwFgcnTpzgzTffBCAmJsalCX/23Fk2/rERgNDqoSqc5abiTn1RrEN5ZW3uMG5aPce0VENERERExAkqnEVEREREnKDCWURERETECSqcRUREREScoMJZRERERMQJKpxFRERERJyg29GJA29vb+rWrWvfdiVPD0/K+5W3b4vcTNypL4p1KK+szR3GTavnmPVekVyVOnXq8Pvvv7s6DADK+5XnhRYvuDoMEZdwp74o1qG8sjZ3GDetnmOaxhMRERERcYJTM86GYTQBZgMNgd+AZ0zT3FzomFLA28AjgC+wFnjWNM3k/P0zgKeB7AKn3W6a5oGrfA0iIiIiItfdJQtnwzD8gCXAGOBdoDew0DCMuqZpZhU4dBhwO2AA6cAsYBoQnb//LqCXaZpfXLvw5VpLSUlhxowZADz33HMEBQW5LJaM7Ay+T/4egOY1mlPap7TLYhEpae7UF8U6lFfW5g7jptVzzJkZ578CuaZpzsz/+33DMF4EugBfFjjudcDXNM0MwzBqAYHAMQDDMDyBRsBP1yxyuS5SU1MZPnw4AD179nRt4XwugzWJawC4s8qdKpzlpuJOfVGsQ3llbe4wblo9x5wpnG8DdhVqM7Et27AXzqZp5gAZhmGMAIYDB4F783fXB/yBSYZh3AMkAcNM04y7quhFREREREqIM18OLAOcKdR2BlshXJxx+ed8Caw0DMMHqIBtzfMEoDowEphvGMadVxCziIiIiEiJc2bG+QxQeK7fH9s65iJM0zwLYBjGy0AMcGf+FwnbFTjsK8MwVgOdgf9dbtAiIiIiIiXNmRnn3di+8FeQQaHlG4ZhvG8YRkyBJu/866cZhtHOMIx+ha7hB5y9zHhFRERERFzCmRnnb4FShmH0x3anjN5AFWBloeO+B142DGM5cASYAsSbprnPMIzawJuGYewCNgKPAmHA49fkVYiIiIiIXGeXnHE2TTMT6Aj0AFKA/kCUaZqnDcNYbhjG0PxDZwMfAt8B+7Et53gk/xprgYHA+8BJ4GWgy/l7PIuIiIiIuDunfgDFNM0dQKti2jsW2M4DRuU/irvGe8B7VxamlBRPT08qVKhg33YlDzwo7V3avi1yM3GnvijWobyyNncYN62eY04VznLzCA4OJiUlxdVhAFChdAWGhA9xdRgiLuFOfVGsQ3llbe4wblo9x6z3UUBERERE5DpQ4SwiIiIi4gQt1RAHaWlpfPTRRwA89thjlC9f3mWxnD13lp8P/QxA46qN8fP2c1ksIiXNnfqiWIfyytrcYdy0eo6pcBYHx48fZ+DAgQB06tTJpQl/JvsMy/cuB6D+LfVVOMtNxZ36oliH8sra3GHctHqOaamGiIiIiIgTVDiLiIiIiDhBhbOIiIiIiBNUOIuIiIiIOEGFs4iIiIiIE1Q4i4iIiIg4QbejkyLc6bflPfBwdQgiLuNOfVGsQ3llbe4wblo5x1Q4i4N69eqRk5Pj6jAACCodxOttXnd1GCIu4U59UaxDeWVt7jBuWj3HrPuRQERERETkGlLhLCIiIiLiBC3VEAenTp1i4cKFAERHRxMQEOCyWDLPZbL72G4AGlRsQCnvUi6LRaSkuVNfFOtQXlmbO4ybVs8xFc7i4MiRIzz++OMAhIeHuzThT2ef5qtfvwKgdlhtFc5yU3GnvijWobyyNncYN62eY1qqISIiIiLiBBXOIiIiIiJOUOEsIiIiIuIEFc4iIiIiIk5Q4SwiIiIi4gQVziIiIiIiTlDhLCIiIiLiBN3HWRwEBwdz6tQpAPz9/V0aSwW/CgxtPRQAH08fl8YiUtLcqS+KdSivrM0dxk2r55gKZ3Hg6elJ2bJlXR0GAB4eHvh6+bo6DBGXcKe+KNahvLI2dxg3rZ5jWqohIiIiIuIEzTiLg/T0dFatWgVA+/btXfqpMSsni4SUBADqBdVz+adokZLkTn1RrMPqeVU3dqmrQ3Apdxg3rZ5jKpzFweHDh4mOjgZg7969Lk349Kx0Pt/5OQADwgYQVDrIZbGIlDR36otiHcora3OHcdPqOabCWUREROQ6uNgMeOK4TiUYiVwrWuMsIiIiIuIEFc4iIiIiIk5Q4SwiIiIi4gQVziIiIiIiTlDhLCIiIiLiBBXOIiIiIiJO0O3oxEGdOnVITk4GoHLlyi6NpbxfeV5q+RIAZXzLuDQWkZLmTn1RrEN5ZW3uMG5aPcdUOIsDb29vqlev7uowAPD08CSgVICrwxBxCXfqi2Idyitrc4dx0+o5pqUaIiIiIiJO0IyzODhz5gybN28GoEWLFvj7+7ssluycbJJOJgFQM7AmPl4+LotFpKS5U18U61BeWZs7jJtWzzEVzuLgzz//pF27doDtN+br1avnslhOZZ3iw58/BGBA2ACCSge5LBaRkuZOfVGsQ3llbe4wblo9x7RUQ0RERETECSqcRUREREScoMJZRERERMQJKpxFRERERJygwllERERExAkqnEVEREREnKDCWURERETECbqPszioVasWu3btsm+7UrlS5Xiu2XP2bZGbiTv1RbEO5ZW1ucO4afUcU+EsDnx9fWnQoIGrwwDAy9OLSmUquToMEZdwp74o1qG8sjZ3GDetnmNaqiEiIiIi4gTNOIuDzMxMfvnlFwDuuOMOSpUq5bJYzuWe48jpIwBULlMZb0+lq9w83KkvinUor6zNHcZNq+eYKhFxkJSURGhoKOD635g/mXmSOT/MAWBA2ACCSge5LBaRkuZOfVGsQ3llbe4wblo9x5wqnA3DaALMBhoCvwHPmKa5udAxpYC3gUcAX2At8Kxpmsn5++8DJgPBwHbgKdM091yblyEiIiIicn1dco2zYRh+wBJgHlAemAosNAzDt9Chw4DbAQOoBBwHpuVfowqwEHgFqAB8A/zn2rwEEREREZHrz5kvB/4VyDVNc6Zpmtmmab6PrSjuUui414GOpmmmAFWAQOBY/r5o4CfTNJeYppkFjAb+YhhG02vyKkRERERErjNnCufbgF2F2kxsyzb+v8E0c0zTzDAMYwSwH2gBjCvuGqZp5gAJha8hIiIiIuKunCmcywBnCrWdAfwvcPy4/HO+BFYahuFzBdcQEREREXErzhTOZ4DShdr8gfTiDjZN86xpmhnAy0Ad4M7LvYaIiIiIiLtxpnDeje0LfwUZFFq+YRjG+4ZhxBRo8s6/flrhaxiG4QWEFL6GiIiIiIi7cuZ2dN8CpQzD6A/MAnpj+/LfykLHfQ+8bBjGcuAIMAWIN01zn2EY/wXGG4YRDcRhu7tGEvDjtXkZcq3UqFGDTZs22bddKbBUIE/f/bR9W+Rm4k59UaxDeWVt7jBuWj3HLlk4m6aZaRhGR2xF8xvAXiDKNM3T+UVyvGmab2C7z3Nl4Dts93H+Gts9nTFN85BhGF2x3cf5Q+AnINo0zbzr8JrkKvj5+dGiRQtXhwGAt6c3NQNrujoMEZdwp74o1qG8sjZ3GDetnmNO/QCKaZo7gFbFtHcssJ0HjMp/FHeNNUDjKwtTRERERMS19JPb4iA7O5v9+/cDUKdOHXx8fFwWS05uDicyTwBQrlQ5vDy9XBaLSElzp74o1qG8sjZ3GDetnmPOfDlQbiIHDhygfv361K9fnwMHDrg0lhOZJ5i6ZSpTt0y1/x+ByM3CnfqiWIfyytrcYdy0eo6pcBYRERERcYIKZxERERERJ6hwFhERERFxggpnEREREREnqHAWEREREXGCCmcRERERESeocBYRERERcYJ+AEUcVKtWjZUrV9q3XSnAN4DejXrbt0VuJu7UF8U6lFfW5g7jptVzTIWzOPD39+f+++93dRgA+Hj5UC+onqvDEHEJd+qLYh3KK2tzh3HT6jmmpRoiIiIiIk7QjLM4OHfuHCkpKQAEBQXh7e26FMnNyyUjOwOA0j6l8fTQ5zy5ebhTXxTrUF5ZmzuMm1bPMVUi4mD//v1UqVKFKlWqsH//fpfGknY2jYkbJzJx40TSzqa5NBaRkuZOfVGsQ3llbe4wblo9x1Q4i4iIiIg4QYWziIiIiIgTVDiLiIiIiDhBhbOIiIiIiBNUOIuIiIiIOEGFs4iIiIiIE1Q4i4iIiIg4wVp3pZarVqVKFebPn2/fdqWyvmV55PZH7NsiNxN36otiHcora3OHcdPqOabCWRyULVuWRx55xNVhAODr5UvDyg1dHYaIS7hTXxTrUF5ZmzuMm1bPMS3VEBERERFxgmacxUFeXh5ZWVkA+Pr64uHh4dJYcvJyAPDy8HJpLCIlzZ36oliH8sra3GHctHqOacZZHOzbtw8/Pz/8/PzYt2+fS2NJPZvK6PWjGb1+NKlnU10ai0hJc6e+KNahvLI2dxg3rZ5jKpxFRERERJygwllERERExAla4ywiIiJSwurGLi22PXFcpxKORC6HZpxFRERERJygwllERERExAkqnEVEREREnKDCWURERETECfpyoDioVKkSc+fOtW+7UhmfMkQZUfZtkZuJO/VFsQ7llbW5w7hp9RxT4SwOAgMDefrpp10dBgClvEtxd7W7XR2GiEu4U18U61BeWZs7jJtWzzEt1RARERERcYIKZxERERERJ6hwFgfu9BvzqRmpjF4/mtHrR5OakerSWERKmjv1RbEO5ZW1ucO4afUc0xpncZCXl0dmZqZ926WxkMe53HP2bZGbiTv1RbEO5ZW1ucO4afUc04yziIiIiIgTVDiLiIiIiDhBhbOIiIiIiBNUOIuIiIiIOEGFs4iIiIiIE1Q4i4iIiIg4QbejEwe33HILkyZNsm+7kr+PP/fXu9++LddO3dilF9yXOK5TCUYiF+JOfVGsQ3llbe4wblo9x1Q4i4Py5cvz0ksvuToMAPy8/WhVq5WrwxBxCXfqi2Idyitrc4dx0+o5pqUaIiIiIiJOUOEsIiIiIuIEFc7iIDExkcqVK1O5cmUSExNdGkva2TQmfjeRid9NJO1smktjESlp7tQXxTqUV9bmDuOm1XNMa5zFQU5ODkePHrVvu1JuXi6ns0/bt0VuJu7UF8U6lFfW5g7jptVzTDPOIiIiIiJOUOEsIiIiIuIEp5ZqGIbRBJgNNAR+A54xTXNzMce9BvQFAoGfgOdN0/wlf98M4Gkgu8Apt5umeeCqXoGIiIiISAm45IyzYRh+wBJgHlAemAosNAzDt9BxjwOPAW2AisA3wFLDMM4/x11AL9M0yxZ4qGgWERERkRuCM0s1/grkmqY50zTNbNM03weOA10KHVcRGGOa5j7TNM8BU4DaQM384rkRtlloEREREZEbjjNLNW4DdhVqM7Et2/jS3mCakwodE4WtwE4C6gP+wCTDMO7JbxtmmmbcFcYtIiIiIlKinCmcywBnCrWdwVYIF8swjHuBWUA/0zRzDcOoAKwFJgBbgU7AfMMwwkzT/N+VBC7XR4UKFRg+fLh925VKe5cmok6EfVvkZuJOfVGsQ3llbe4wblo9x5wpnM8Ahd99fyC9uIMNw+gNvAP0N03zU4D8LxK2K3DYV4ZhrAY6Ayqc3UhQUBAjR450dRgAlPYpzV+D/+rqMERcwp36oliH8sra3GHctHqOOVM47waeL9RmAJ8WPtAwjGHAC0BX0zS/LdDeDggxTXN2gcP9gLOXHbGIOK1u7FJXhyAiImIZzhTO3wKlDMPoj235RW+gCrCy4EGGYTwBvAi0Mk3z10LXyAHeNAxjF7AReBQIAx6/quhFRERERErIJe+qYZpmJtAR6AGkAP2BKNM0TxuGsdwwjKH5h74CBADbDMNIL/BoYJrmWmAg8D5wEngZ6GKaZvK1f0lyNQ4cOED9+vWpX78+Bw649m6BJ86eYOqWqUzdMpUTZ0+4NBaRkuZOfVGsQ3llbe4wblo9x5z6ARTTNHcArYpp71hg+9ZLXOM94L3LDVBKVnZ2Nnv37rVvu1JOXg4pGSn2bZGbiTv1RbEO5ZW1ucO4afUcc6pwFhH3cKE1y4njOpVwJCLyf+3de5xdZXno8d/suWRyvxCSgAkksO1LxSqolYqCYjweI6BCzynVo5YeDUoVqKei1GJbrB7B1hv9eMBUBCwFpElcKF8AACAASURBVIJBjFQUREVFKUgBxVeISSBM7veZyVz3Pn/syc6sSUJWksmsNWt+Xz/5ZM0ziz3PNs+732evefe7JI09aW6AIkmSJI15Ns6SJElSCjbOkiRJUgqucZYkSaOKe9QrK15xliRJklLwirMSpk6dyiWXXFI/zlJrUyt/NPeP6scaGe7ckQ95GosqDuuq2PIwbxa9xmyclTBz5ky+8IUvZJ0GABOaJ/Cm8puyTiMT/hpSeRqLKg7rqtjyMG8WvcZcqiFJkiSlYOMsSZIkpWDjrITVq1fzile8gle84hWsXr0601y2d29nycNLWPLwErZ3b880F2mk5Wksqjisq2LLw7xZ9BpzjbMSuru7efjhh+vHWeqr9NG2o61+LI0leRqLKg7rqtjyMG8Wvca84ixJkiSlYOMsSZIkpWDjLEmSJKVg4yxJkiSl4IcDpQyNphudeEdBSdJY5xVnSZIkKQWvOCth8uTJnH/++fXjLI1rHMdJc06qH2vfRtOVa6WTp7Go4rCuii0P82bRa8zGWQmzZs3i+uuvzzoNACa2TORtJ7wt6zSkTORpLKo4rKtiy8O8WfQac6mGJEmSlIKNsyRJkpSCjbMS2traWLhwIQsXLqStrS3TXHZ07+DGR2/kxkdvZEf3jkxzkUZansaiisO6KrY8zJtFrzHXOCth586d3HffffXjLPVWelmxdUX9WBpL8jQWVRzWVbHlYd4seo3ZOEsjwF0vJEka/VyqIUmSJKVg4yxJkiSlYOMsSZIkpeAaZ0mSlEt+PkR54xVnSZIkKQWvOCth4sSJnHvuufXjLLU0tvD7M3+/fiyNJXkaiyoO66rY8jBvFr3GbJyVMGfOHG6//fas0wBgUsskznvxeVmnIWUiT2NRxWFdFVse5s2i15iNsyRJUk7sa133yivPHOFMtDeucZYkSZJSsHFWwtq1aznnnHM455xzWLt2baa5tPe0c+sTt3LrE7fS3tOeaS7SSMvTWFRxWFfFlod5s+g1ZuOshI6ODpYuXcrSpUvp6OjINJee/h5+s/E3/Gbjb+jp78k0F2mk5Wksqjisq2LLw7xZ9BqzcZYkSZJSsHGWJEmSUrBxliRJklKwcZYkSZJSsHGWJEmSUrBxliRJklLwzoFKGD9+PAsXLqwfZ6m51Mxx04+rH0tjSZ7GoorDuiq2PMybRa8xG2clHH300Xz/+9/POg0AJo+bzLtf+u6s05AykaexqOKwrootD/Nm0WvMpRqSJElSCjbOkiRJUgo2zkpYv349559/Pueffz7r16/PNJeOng6W/mYpS3+zlI6e4t22U3o+eRqLKg7rqtjyMG8WvcZsnJWwY8cObrzxRm688UZ27NiRaS7d/d08uvZRHl37KN393ZnmIo20PI1FFYd1VWx5mDeLXmM2zpIkSVIKNs6SJElSCjbOkiRJUgo2zpIkSVIKqW6AEkI4GfgycCLwFPD+GOODeznvcuACYArwKPDBGOMTA997A/AFYAHwCPCeGONvh+NJSJIkSYfbfq84hxBagbuA64FpwNXAHSGEliHnnQ+8G3gdMBP4PrAshFAKIcwG7gD+Gpg+8L1bhu1ZSJIkSYdZmivOZwCVGOM1A19/NYTwIeBs4PZB580EPhVj/B1ACOGLwD8Ac4EzgUdjjHcNfO+TwF+GEF4eY3x4eJ6KhsO4ceN4xSteUT/OUlOpiaMnH10/lsaSPI1FFYd1VWx5mDeLXmNp/l89Afj1kFiktmyj3jjHGP9pyDlvATYBq4c+RoyxP4SwfOAxbJxzZO7cuTz00ENZpwHAlHFTuODlF2SdhpSJPI1FFYd1VWx5mDeLXmNpPhw4EegcEusEJuzrPwghnA5cC1wcY6wczGNIkiRJeZKmce4Exg+JTQDa93ZyCOFdwDLgohjjzQfzGJIkSVLepGmcnwTCkFhgz+UbhBA+Tm3njLfGGG/Y12OEEBqB8t4eQ9nauHEjl1xyCZdccgkbN27MNJfO3k7ufupu7n7qbjp7h/7CQiq2PI1FFYd1VWx5mDeLXmNp1jjfB4wLIVxEbfnFu4DZwHcHnxRC+HPgQ8CpMcbfDHmMbwJXhRDOBb5NbXeN1cAvDy19Dbdt27Zx9dVXA3DxxRczc+bMzHLp6uvi58/9HIBT5p7ChGZX9mjsyNNYVHFYV8WWh3mz6DW23yvOMcZuYBHwdmAzcBHwlhhjRwjh7hDCxwZO/WtgMvCfIYT2QX9+P8a4Fngr8HfUPjD4BuDcGGP1MDwnSZIkadil2qskxvgYcOpe4osGHf/efh7jB8BLDzRBSZIkKQ+85bYkSZKUgo2zJEmSlIKNsyRJkpSCjbMkSZKUQjY3MlduNTc388IXvrB+nKXGhkaOGH9E/VgaS/I0FlUc1lWx5WHeLHqN2Tgr4ZhjjuG3v/1t1mkAMLV1KhedclHWaUiZyNNYVHFYV8WWh3mz6DVm4ywNo/mXLcs6BUmSdJjYOEuSpEx50UGjhY2zEjZv3sznP/95AD70oQ8xY8aMzHLZ2buTn63+GQCvmvsqxjePzywXaaTlaSyqOKyrYsvDvFn0GrNxVsKWLVv45Cc/CcD555+fbePct5MfrfoRACfNOcnGWWNKnsaiisO6KrY8zJtFrzG3o5MkSZJSsHGWJEmSUrBxliRJklJwjbOkQ7KvT8OvvPLMEc5EkqTDy8ZZeh42hZIkaReXakiSJEkpeMVZCY2NjcyePbt+nKVSQ4lJLZPqx9JYkqexqOKwrootD/Nm0WvMxlkJ8+fPZ+3atVmnAcC01ml8+NQPZ52GlIk8jUUVh3VVbHmYN4teY17GkyRJklLwirN0EPb1oUFJklRcNs5K2Lp1K1/5ylcAeO9738u0adMyy6Wrr4tH1jwCwMuOehmtTa2Z5SKNtDyNRRWHdVVseZg3i15jNs5K2LRpE5deeikA55xzTqYF39nbyT3L7wHghJkn2DhrTMnTWFRxWFfFlod5s+g15hpnSZIkKQUbZ0mSJCkFG2dJkiQpBRtnSZIkKQUbZ0mSJCkFG2dJkiQpBbejU0JDQwOtra3140xzoYGmUlP9WBpL8jQWVRzWVbHlYd4seo3ZOCvhuOOOY+fOnVmnAcD08dO5/PTLs05DykSexqKKw7oqtjzMm0WvMZdqSJIkSSnYOEuSJEkpuFRDCdu3b+frX/86AOeddx5TpkzJLJfuvm6eWP8EAC+e9WLGNY3LLBdppOVpLKo4rKtiy8O8WfQas3FWwoYNG7jgggsAeP3rX59pwXf0dnDXb+8CYMH0BTbOGlPyNBZVHNbV6DX/smV7ja+88sz6cR7mzaLXmEs1JEmSpBRsnCVJkqQUbJwlSZKkFGycJUmSpBRsnCVJkqQUbJwlSZKkFGycJUmSpBTcx1kJxx13HN3d3QA0Nzdnmsv01ul8/PSPA1Bq8D3eaJNmz1HtW57GoorDuiq2PMybRa8xG2clNDQ00NLSknUaQC2XxobGrNOQMpGnsajisK6KLQ/zZtFrzMt4kiRJUgpecVZCe3s73/nOdwB485vfzKRJkzLLpae/h6c2PQXAC494IS2NxX0HKw2Vp7Go4rCuii0P82bRa8zGWQnr1q3jvPPOA+Dpp5/OtODbe9r591//OwAXn3IxM8bPyCwXaaTlaSyqOKyrYsvDvFn0GnOphiRJkpSCjbMkSZKUgo2zJEmSlIKNsyRJkpSCjbMkSZKUgo2zJEmSlIKNsyRJkpRCqn2cQwgnA18GTgSeAt4fY3zwec7/ItAbY/zwoNiXgPcCvYNOfVGM8ZmDSVyHx/z589mwYQMA06dPzzSXaa3T+MirPwJAa1NrprlIIy1PY1HFYV0VWx7mzaLX2H4b5xBCK3AX8CngK8C7gDtCCPNjjD1Dzj0C+CzwZwN/D3YS8L9ijN8YjsR1eDQ2NjJz5sys0wCg1FBiQvOErNOQMpGnsajisK6KLQ/zZtFrLM1SjTOASozxmhhjb4zxq8Am4Oy9nPsA0AfcPjgYQigBLwEePcR8JUmSpEykWapxAvDrIbFIbdnG7UPiC2OMbSGEG4bEXwhMAP4phPBqYDXw8Rjjtw88ZR1OnZ2d/OhHPwLg9NNPZ8KE7N659vb3smrbKgCOnXoszY3NmeUijbQ8jUUVh3VVbHmYN4teY2ka54lA55BYJ7VGOCHG2LaPx5gO3A98BngIOBO4LYRwSozx8dTZ6rBbs2YNixYtAmr3mD/++OMzy2VHzw5ueuwmAC4+5WJmjJ+RWS7SSMvTWFRxWFfFlod5s+g1lqZx7gTGD4lNANrT/pCBDxIuHBRaGkK4FzgLsHGWJElS7qVZ4/wkEIbEAnsu39inEMLCEML7hoRbga60jyFJkiRlKc0V5/uAcSGEi4Brqe2qMRv47gH8nH7gsyGEXwM/Bf4EOAU4/4CylSRJkjKy38Y5xtgdQlhErWn+v8DTwFtijB0hhLuBH8cY/+9+HuP+EMIlwFeBo6l9uPDsGONzh/wMJEnSqDD/smVZpyAdklQ3QIkxPgacupf4on2cf/5eYtcB1x1gfpIkSVIupGqcJWm47OuK08orzxzhTCRJOjBpPhwoSZIkjXlecVbCMcccw9NPP10/ztLUcVO5+JSL68fSWJKnsajisK6KLQ/zZtFrzMZZCc3NzbnZrLyx1OhNTzRm5Wksqjisq2LLw7xZ9BpzqYYkSZKUglecldDV1cUvf/lLAE4++WRaW1szy6Wv0seaHWsAOGryUTSVLFeNHXkaiyoO66rY8jBvFr3G7ESU8Nxzz3HqqbWdB7O+x/z27u1c98vaDoYXn3Jx5r9+kkZSnsaiisO6KrY8zJtFrzGXakiSJEkp2DhLkiRJKdg4S5IkSSnYOEuSJEkp2DhLkiRJKdg4S5IkSSm4HZ0EzL9sWdYpSJKknLNxVsLcuXN5+OGH68dZmjJuCu97+fvqx9JYkqexqOKwrootD/Nm0WvMxlkJ48aN42Uve1nWaQDQVGriqMlHZZ2GlIk8jUUVh3VVbHmYN4teYzbOkiRJo9S+lhquvPLMEc5kbLBxVkJPTw9PP/00AOVymZaWlsxy6a/0s3nnZgBmjJ9BY6kxs1ykkZansajisK6KrUo/FdqB2hyaxbxZ9BqzcVbCs88+y4knnghkf4/5bd3b+NJDXwLg4lMuZsb4GZnlouyM1aspeRqLKg7rqtgqdLKj6dsAbOt+UybzZtFrzMZZUi64s4kkKe/cx1mSJElKwcZZkiRJSsHGWZIkSUrBxlmSJElKwcZZkiRJSsHGWZIkSUrB7eiUcPTRR3PffffVj7M0uWUyf/bSP6sfS2NJnsaiisO6KrYS45nYvxDIbt4seo3ZOCth/PjxnHHGGVmnAUBzYzMLpi/IOg0pE3kaiyoO66rYGmiiuToHqM2hWSh6jblUQ5IkSUrBK85K6OvrY+3atQDMmTOHpqbsSqRSrdDe0w7ApJZJlBp8n6exI09jUcVhXRVblQpVuoDaHJrFvFn0GrMTUcKqVauYN28e8+bNY9WqVZnmsrVrK5/72ef43M8+x9aurZnmIo20PI1FFYd1VWwVOtje9E22N30zs3mz6DVm4yxJkiSlYOMsSZIkpWDjLEmSJKVg4yxJkiSlYOMsSZIkpWDjLEmSJKVg4yxJkiSlUKxdqXXI5syZw9KlS+vHWZrUMok/ffGf1o+lsSRPY1HFYV0VW4lWJvafDmQ3bxa9xmyclTBx4kTe+ta3Zp0GAC2NLZww84Ss05AykaexqOKwroqtgWaaq/OA2hyahaLXmEs1JEmSpBS84qyESqVCR0cHUHvXWCpl996qWq3S098D1N45NzQ0HPJjzr9s2SE/hjQS8jQWVRzWVbFVqQJ9teNqdVjmzQNV9Bor1rPRIVuxYgVTpkxhypQprFixItNctnRt4dMPfJpPP/BptnRtyTQXaaTlaSyqOKyrYqvQzram29jWdFtm82bRa8zGWZIkSUrBxlmSJElKwcZZkiRJSsHGWZIkSUrBxlmSJElKwcZZkiRJSsHGWZIkSUrBG6AoYdasWdx444314yxNbJ7IOSecUz+WxpI8jUUVh3VVbCVamdD/KiC7ebPoNWbjrITJkyfz7ne/O+s0ABjXNI6XznnpAf933h1QRZCnsajisK6KrYFmWqrHAbU5NAtFrzGXakiSJEkppLriHEI4GfgycCLwFPD+GOODz3P+F4HeGOOHB8XeAHwBWAA8ArwnxvjbQ8hdkiRJGjH7veIcQmgF7gKuB6YBVwN3hBBa9nLuESGEG4CLh8RnA3cAfw1MB74P3HKoyWv4LV++nFKpRKlUYvny5ZnmsnnnZq64/wquuP8KNu/cnGku0kjL01hUcVhXxdbPDrY23czWppszmzeLXmNplmqcAVRijNfEGHtjjF8FNgFn7+XcB4A+4PYh8XOBR2OMd8UYe4BPAseFEF5+CLnrMKlWq1Sr1azTAKA68D9pLMrTWFRxWFdFVx34k2EGBa6xNI3zCcCvh8QitWUbQy2MMb4XaH++x4gx9gPL9/EYkiRJUu6kaZwnAp1DYp3AhKEnxhjbDvUxJEmSpDxK0zh3AuOHxCaw51Xlw/0YkiRJUmbSNM5PAmFILLDn8o3UjxFCaATKB/gYkiRJUmbSbEd3HzAuhHARcC3wLmA28N0D+DnfBK4KIZwLfJva7hqrgV8eWLqSJElSNvZ7xTnG2A0sAt4ObAYuAt4SY+wIIdwdQvhYisdYC7wV+DtqO3K8ATg3xljMj1xKkiSpcFLdACXG+Bhw6l7ii/Zx/vl7if0AOPD7J2tEzZw5k6uvvrp+nKUJzRNYVF5UP5bGkjyNRRWHdVVsJcYxvv8VQHbzZtFrLFXjrLFj6tSpXHTRRVmnAUBrUyunzD0l6zSkTORpLKo4rKtia6CFcdXaR8pam1ozyaHoNZbmw4GSJEnSmGfjLEmSJKVg46yEFStWMH36dKZPn86KFSsyzWXLzi1c+cCVXPnAlWzZuSXTXKSRlqexqOKwroqtn3a2Nf0725r+PbN5s+g15hpnJVQqFbZu3Vo/zlKVKl19XfVjaSzJ01hUcVhXRVelSk/9KAtFrzGvOEuSJEkpeMVZ0qg0/7Jle42vvPLMEc5EkjRWeMVZkiRJSsErzpIkadjs67dBUhF4xVmSJElKwSvOkgrFtc+SpMPFxlkJM2bM4B/+4R/qx1ka3zSe1y94ff1YGkvyNBZVHNZVsTXQQmvlpUB282bRa8zGWQnTp0/n8ssvzzoNAMY3j+f0Y0/POg0pE3kaiyoO66rYSoyjtfJioDaHZqHoNeYaZ0mSJCkFG2dJkiQpBRtnJaxatYpjjz2WY489llWrVmWay9aurXz+Z5/n8z/7PFu7tmaaizTS8jQWVRzWVbFVaGd741K2Ny7NbN4seo25xlkJfX19PPPMM/XjLFWqFbZ1b6sfS2NJnsaiisO6KrYqVSoNHUB282bRa8wrzpIkSVIKNs6SJElSCjbOkiRJUgo2zpIkSVIKNs6SJElSCjbOkiRJUgpuR6eEadOmcemll9aPs9Ta1Mqr5726fiyNJXkaiyoO66rYGmhmXOVFQHbzZtFrrKFarWadw/MKIcwHVtx7773MnTs363Q0Csy/bFnWKSiHVl55ZtYpCGD5ciiXa8dPPw3HH59tPjpovtbmm695B2f16tUsXLgQYEGMceXQ77tUQ5IkSUrBxlmSJElKwTXOSnj22Wc588zar3eWLVvGvHnzMstlW9c2bn78Zi5f+gQT+19LiYmJ7/trKB2Iff1aOa91lKexqOKwroqtQgcdjT8EYFvXa5jaOnXEcyh6jdk4K6Gnp4fHH3+8fpyl/mo/6zrW0d+whSqVTHORRlqexqKKw7oqtioV+hu2ALU5NAtFrzEbZ0mSpIJ56RX30MjkPeJ5/S3baOEaZ0mSJCkFG2dJkiQpBRtnSZIkKQUbZ0mSJCkFG2dJkiQpBXfVUMKUKVO44IIL6sdZGtc4jpcf9XJaKhUaaM40F2mk5Wksqjisq2JroJmWSrl+nIWi15iNsxKOPPJIvvzlL2edBgATWyZydjibCRV/MaKxJ09jUcVhXRVbiVYmVE7JNIei15iNs3JjX3d2G67zJUmSDoWX8iRJkqQUbJyV8Nxzz3Haaadx2mmn8dxzz2WaS4VO2hvvob3xHip0ZpqLNNLyNBZVHNZVseVh3ix6jblUQwldXV088MAD9eMsVemnr2FD/VgaS/I0FlUc1lWx5WHeLHqNecVZkiRJSsHGWZIkSUrBxlmSJElKwcZZkiRJSsHGWZIkSUrBxlmSJElKwe3olDBp0iTOO++8+nGWGmiiuXJs/VgaS/I0FlUc1lWx5WHeLHqN2Y0oYfbs2dx6661ZpwFAifFMrLwm6zSkTORpLKo4rKtiy8O8WfQac6mGJEmSlIKNsyRJkpSCjbMS1qxZw1lnncVZZ53FmjVrMs2lQicdjffT0Xg/FTozzUUaaXkaiyoO66rY8jBvFr3GXOOshM7OTpYtW1Y/zlKVfnobngOglf5Mc5FGWp7GoorDuiq2NPPm/MuW7TW+8sozhyWHotdYqsY5hHAy8GXgROAp4P0xxgf3ct5fApcCk4FvAe+LMXYMfG8Z8HrY/S8ZYyzexy0lSZJUSPtdqhFCaAXuAq4HpgFXA3eEEFqGnHcWtab5DGAeMAP4xKBTTgJOizFO2vVneJ6CJEmSdPilWeN8BlCJMV4TY+yNMX4V2AScPeS8dwHXxRh/G2PcBnwceE8IoTGEMAuYBTwxnMlLkiRJIyVN43wC8OshsUht2cbznReBqcALgJOBHcC3QwgbQgg/CSG86uBSliRJkkZemsZ5Iuzx0cxOYMJ+ztt1PAFoBX4GXALMBW4C7g4hzDnQhCVJkqQspPlwYCcwfkhsAtC+n/N2NdbtMcY7gTsHfe+aEMJfUFsGckv6dCVJkqRspGmcnwQ+OCQWgJv3cl4Ycs42oC2E8D+AUozxtkHfbwW6DixdHW4TJkxg0aJF9eMsNdBIU/Xo+rE0luRpLKo4rKtiy8O8WfQaS9M43weMCyFcBFxL7UOAs4HvDjnvJuDaEMLtwLPUdtT4txhjJYQwCbgyhPAEte3s/pLa1el7hudpaLgcddRRfOc738k6DQBKTGBS/xlZp6GC29eepjB8+5oejDyNRRWHdVVseZg3i15j+13jHGPsBhYBbwc2AxcBb4kxdoQQ7g4hfGzgvLuAq4BlwDPAVmrb0xFjvAH4IvAfA/G3AIt27fEsSZIk5V2qG6DEGB8DTt1LfNGQr6+mts/z3h7j08CnDyJHSZIkKXNpdtXQGLJu3Tre8Y538I53vIN169ZlmkuFnXSUfkJH6SdU2JlpLtJIy9NYVHFYV8WWh3mz6DVm46yE9vZ2brnlFm655Rba24dunDKyqvTRW1pJb2klVfoyzUUaaXkaiyoO66rY8jBvFr3GbJwlSZKkFFKtcZaksWhfO25kuduGJCk7XnGWJEmSUrBxliRJklKwcZYkSZJSsHGWJEmSUvDDgUpobW3l1FNPrR8fDs93i+PBGmiksXpk/VgaS0ZiLGrssa6KLQ/zZtFrzMZZCS94wQv4yU9+MiyPlbZB3pcSE5jc/8ZhyUUabYZzLEq7WFfFlod5s+g1ZuMsSZL26VAvgkhFYuMsSQfI/Z0laWzyw4FK2LBhAxdeeCEXXnghGzZsyDSXCl10ln5BZ+kXVOjKNBdppOVpLKo4rKtiy8O8WfQas3FWwvbt27n22mu59tpr2b59e6a5VOmlp/QUPaWnqNKbaS7SSMvTWFRxWFfFlod5s+g1ZuMsSZIkpWDjLEmSJKVg4yxJkiSlYOMsSZIkpWDjLEmSJKXgPs46ZG6OL0nS6Ob+9OnYOCuhpaWFE088sX6cpQZKNFan1Y+lsSRPY1HFYV0VWx7mzaLXmI2zEubNm8cTTzyRdRoAlJjI5H7f6WpsytNYVHFYV8WWh3mz6DXmZTxJkiQpBRtnSZIkKQWXaihh06ZN/OM//iMAl156KUcccURmuVToorv0GwDGVU6gRGtmuUgjLU9jUcVhXRVbHubNoteYjbMStm7dylVXXQXA4sWLMy34Kr10l34FQEvleLBx1hiSp7Go4rCuii0P82bRa8zGWZIkSXvlNnVJrnGWJEmSUrBxliRJklJwqYYkSfIusFIKXnGWJEmSUvCKsyQNEz9EI0nFZuOshMbGRubOnVs/zlIDDZSYUD+WxpI8jUUVh3VVbHmYN4teYzbOSpg/fz7PPvvs3r83wuvfSkxiSt85I/ozpbx4vrEoHSzrqtjyMG8WvcZsnCXpMHMJhyQVgx8OlCRJklLwirMStmzZwjXXXAPAhRdeyPTp0zPLpUI3PaWnAGipvJAS4zLLRRppeRqLKg7rqtjyMG8WvcZsnJWwefNm/uZv/gaA8847L9OCr9JDV+m/AGiuHAs2zhpD8jQWVRzWVbHlYd4seo3ZOGufTv/MD2ie/pus05AKy7XPypKv8dKBGzWN8zPPPMOMGTOYMGFCPbZ8+XIA5syZw8SJE+vxFStWUKlUmD17NpMmTarHV65cSX9/P7NmzWLy5Ml7xI888kimTJmS+Jm9vb3MnDmTqVOn1uPPPvssPT09HHHEEUybNq0eX716Nd3d3cyYMSPxDuu5556jq6uL6dOnM2PGjHq8ra2NnTt3Mm3aNI444oh6fM2aNXR2djJ16lRmzpxZj69du5aOjg6mTJnCkUceWY+vW7eO9vZ2Jk+ezKxZs+rx9evXs2PHDiZNmsTs2bPr8Q0bNrB9+3YmTpzInDlz6vGNGzeyatWqPf6/79+5nUpXBw1NLTRNPmJQfAeVrnYampppmrw7z/6udio7d9DQ2EzTlN3xSlc7/Tt30NDYRNOU3flXujvo79xOQ2MjTVN251/p2Ul/XzsNDQ0waEebSncn/Z3baCiVaJo6O3l+x1ZoaKB52pxB8S76O7bsGe/tor99CwDN048a2Ubz3wAAEfBJREFUFO+mv30zAE3T5tR+PlDt66Fvx6aB+GwaGkoD8V76dmysxafOoqFUS7ba30vf9r3F++jbvqEWn3IkDY21YVit9NO3bf1AfCYNjc17xifPpKFpIF6t0Ld13UD8CBqaWgbiVfq2rgWgcdIMSs27rzj0blkzEJ9Oqbl1d3zrWqhWaZw4nVLL3uLTKLWMr8f7tq2jWqnQOGEqpXETBsXXU6307xnfvp5qfz+NE6ZQGjdxUHwD1f4+GsdPptQ6aVB8I9X+XkrjJ9M4OL5jI9W+Xkqtk2gcP3lQfBPVvh5KrRNpHL97DPe1b6ba201p3AQaJ+wew/3tW6j0du0Z79hCpaeLUst4GidOGxTfSqVnJ6WWVhon7h7b/Z3bqHR3UmpupXHSnvGG5nE0Tdo95g90LO2yZs0ajj/++PrXW7duZdOmTbS0tDBv3rx6fNu2bWzcuJHm5maOOeaYenz79u1s2LCBxsZG5s+fX4/v2LGD9evX7xFvb29n3bp1lEolFixYUI93dHSwdm2ttgbn09nZyZo1a/aI79y5k7a2NgCOO+64+ljaZdWqVcydP7++ZVV3dzerV68G4Nhjj6WpqTY2enp66p/SHxzv7e3lmWeeAWDevHm0tNTGQF9fX/21bO7cuYwbVxsD/f39rFy5EoAXvOAFtLbWar1SqbBixQoAjj76aMaPr9V6tVrld7/7HQBHHXVUqvnnd7/7HdVqdY/5Z9e8lHb+WbVqFX19fXvMP7vmpaHzz655aej8s2teGjr/QG1c7jFmdo2lIWNjv2NpjzGzaywNGTO7xlLaMbMrvseY2TWWhsw/u8bSkPln97w0ZP6pz0tD5p9d81Kpkaapg+P7mH8G4nvOP/uZlxg6/+yelwbPP8l5afD8M2hemjqr/sm1aqU2P1Sq7fuelwbPP4l5KcX8U+mvj4GhY2lvfcRoGku7Xsv2ZdR8OPC1r30tDzzwQCL2ohe9iHK5zL333puIn3TSSZTLZe6+++5E/JWvfCXlcpk777wzET/ttNMol8vcdtttifjChQspl8vcdNNNifiiRYsol8tcd911ifjb3vY2yuVyfW3PLn/yJ39CuVzmi1/8YiL+zne+k3K5zGc+85lE/D3veQ/lcplPfvKTifiFF15IuVzmb//2bxPxSy65hHK5zEc/+tFE/CMf+QjlcpkPfehDifjll19OuVzmAx/4QCL+iU98goULFzLU9ge/QduSxWxa9rlEfMdDS2lbspiNdybz3/HIt2lbspgN3/xUMv5f99C2ZDHrv3FFIt7+xH20LVnMulsvT8Q7f/tT2h++i/bHv5+MP/Uz2pYsZs2//lUivnP5Q7X4DZck4l2rHq3Fr/uLRLz72V/RtmQxbUsuSMR72uJAfDFU+nbH1y2vx6s9XbvjG1fV45Wu9nq8d3NbPb7rhRBqzeKu+K4XKqhNQLvivZvb6vFKV3s93rNx9wtStadrd3zd8t1PoNK3O94WE8+tbckFtC1ZTPezv0rE11z3F7QtWUzXqkeT8RsuoW3JYnYufygZ/9e/om3JYjqf+lkivvbmj9K2ZDEdT/4oEV936+W0LVlM+xP3JeLrv3EFbUsWs+O/7knEN3zzU7X4I99OxDfe+Zla/KGlifimZZ+jbclitj/4jUR8891X07ZkMdt+emsy/r1raFuymK0/+tdk/N6v0LZkMVvuvz4R33L/DbQtWczm7y9JxLf++Kba2LjnS4n4tp9+vRa/Oznmt//89tqY+fZnE/Ed//mtgbF0FXtz4YUXJr6+7rrrKJfLLFq0KBG/6aabKJfLe4zj2267jXK5zGmnnZaI33nnnZTLZV75ylcm4nfffTflcpmTTjopEb/33nspl8u86EUvSsQfeOAByuUy5XI5Ef/FL35Rj/f09OzxvF6/cCHt7bvHzJNPPlk/f/PmzfX48uXL6/HBk9qzzz5bjw+erNevX1+PP/XUU/X41q1b6/EnnniiHu/s7KzHH3nkkXq8v7+/Hn/wwQcTuf/e7/0e5XKZH/7wh4n4H/zBH1Aul/ne976XiL/85S+nXC6zbFnytwyvetWrKJfL3HHHHYn46173OsrlMrfemqzdN77xjZTLZb72ta8l4meddRblcpl/+Zd/ScT/+I//mHK5zJe+lKxRgPW3/g3bH0rOh5uWfX7vY+k//rk2ln5ySzJeH0vJfLbcd11tLP3gq8n4D2+sjaV7h4ylB/6tFv/ukLH0s9v2PpZ+ccfex9LDtbG04c5PJ+Ltj95di9+RnFfbH/9+bV667e8T8Y5f3T8wL30sEe+MD9C2ZDFrb/pIMv70z2vxr/2fRHzn7x6uzT/XX5SIdz3zWO11+ivJsd29+snd8w/Verxn7VO755++3t3x9SsGzUs76/HKzm2s/dpf1ualzu31eN/WdYPmpd1jrL990+55aaBRhtobl93z0urdj9+zsz42Hn/88d3Pq6trr31EpVKpn//Tn/408b0TTjiBcrnM/fffn4i/5CUvoVwu893vfjcR/8M//EPK5TJ33XVXIv7qV7+acrnM7bffnoifccYZlMtlbr755kT8TW96E+VymRtuuCERf+9737tH/oONmsZZkiRJylJDtVrd/1kZCiHMB1Zcf/31nHTSSS7VOMxLNR577LH6u8WjL/gXmqcfldlSjd6e9WwvfZOGhgamNr6dRiYPnO9SDXCpRpGXaqy88sz6VVaAH//4x7zmNa+pnz+ql2osXw4Dz2vVvfcy97WvdanGCC7VGFxXs/70U7TMWuBSjQIt1aiUOtnR9C2qlQoTNr+KUnXiYVuqsTs+eP6pcP1/n1jvI55++mmOP/74UTWWHnroId75zncCLIgxrmSIUdM433vvvfVbOOrwGfyiuqtxzko/O9jR9C0AJve9pd44S0U3tHHeNfkUwqDGmaefhqI8r1EiT6/xGn55mDfvXXzCqH7tWr169a7Gf6+N86j5cKBGRqlUoqFl4J1gQzb3ud+tgQaa68fSWFIqlepXTUql519V5+4cSitfr/EaftnPmwfy2jUa2TgrYcGCBRzzodv2f+IIaGQSU/v+JOs0pEwsWLCAbdu2ZZ2GRrF9vaHKy2u8hl8e5s2iv3bZOI8BXo2SJEk6dMW7hi5JkiQdBl5xVsK2bdvqe+dOPPGMxA4II61KDz0NtU+6t1QX0EBLZrlII23btm31PeTf+c53Jj4NLh2sSncHHb/6AZD9a7yGXx7mzaK/dtk4K2Hjxo1s/t61ALQueHmmL6oVutnZ+J8ANPUdTaONs8aQjRs38sEPfhCobdQ/derUfS67ktLq79yem9d4Db88zJt7e+0qEhtnSZIkDYvTP/ODxHHz9N8AxflclY3zGLa3q1e7bpAhqbj8wLAkHZxUjXMI4WTgy8CJwFPA+2OMD+7lvL8ELgUmA98C3hdj7Bj43tuBTwGzgPuB98QY1w3DcxD7nggljT7zL1uWeBM7+KqNxi7f8Gg0O5g+JY+1vd/GOYTQCtxFren9CvAu4I4QwvwYY8+g886i1jSfAawDbgE+AfxVCOElwLXAG4HHgH8G/h/wx8P6bHLsQF/wfIGUpLHpQBsML5xIIyfNFeczgEqM8ZqBr78aQvgQcDZw+6Dz3gVcF2P8LUAI4ePAvSGEjwD/C7gzxvjzge99FFgXQpgVY1w/TM8lF3zBk6TnefO/+IQRzkSShk+afZxPAH49JBapLdt4vvMiMBV4wdDvxRg3AVsH4pIkSVLupbniPBHoHBLrBCbs57xdxxMO4DH2phFg7dq1KU4dfq+56gd7jT/w0TP2/h90bD6M2cD8i/71sD5+w85tNDU11Y9pGXdYf97z5kI7pabu2nHfFqA3s1ykkTYcY/FAXy8O5vXlQF8LX33FUr4+8LzOu2Ipz02d9fyPMxYdxnkkT6/xGn55mDeHs8ZWr149XGmlNqjfbNzb9xuq1erzPkAI4f8A/y3GuGhQ7BvAozHGTw6KPQZ8Ksb49YGvJwE7gHnU1jP/JMZ41aDzNwJvizE+sJ+f/xrgx8+bpCRJkjR8Tttbj5rmivOTwAeHxAJw817OC0PO2Qa0Df1eCGEmMGMgvj8PAacBa4D+FOdLkiRJB6MROIpa/7mHNFecxwG/A66ktjPGuwaOF+zaam7gvLPZvXPGs9R21VgZY/xACOEk4IfAmcB/UttV4+gYo1tESJIkaVTY74cDY4zdwCLg7cBm4CLgLTHGjhDC3SGEjw2cdxdwFbAMeIbah/8uHfjeo8Bi4KvAeuBo4M+H/dlIkiRJh8l+rzhLkiRJSrcdnSRJkjTm2ThLkiRJKdg4S5IkSSnYOEuSJEkppNnHWQUXQjgZ+DK126g/Bbw/xvhgtlnpQA3cLOiz1G5lvxH4TIzxy9lmpYMRQpgNPA787xjjt7PORwcmhDCX2vaspwPbqY3Fq7PNSgcihHAqcDXwe9TuI3FFjHHo/SuUUyGEVwJLY4xHD3w9ndrObq+ndo+RK2KM1x3MY3vFeYwLIbQCdwHXA9OovVDcEUJoyTQxHZCBF4VvUfv3mw78T+DTIYQ3ZJqYDtZ1wBFZJ6EDF0JoAJZSu8HXEcB/B/5+oBHTKBBCaKT2b3hljHEK8F7gxhDC/EwT036FEBpCCP8buAcY3Mf8C9AOzAb+B/CZEMJLDuZn2DjrDKASY7wmxtgbY/wqsAk4O+O8dGCOBZbFGP8txliJMT4C/ABwsh5lQgjvBzqo3UhKo88p1O5VcNnAa+qvgFcBMdu0dACmAUcCTQNvhCpAD969eDT4GHAJ8KldgRDCJOBtwN/FGLtijL+gdvfrxQfzA2ycdQLw6yGxSG3ZhkaJGOOjMcZ37fp64Ar0acB/ZZeVDlQI4YXAXwEXZp2LDtrLgF9Ru6K1NoTwW+CPYoybMs5LKQ38W/0/andA7gV+DHwwxuib2fz7KnASydtlvxDojTH+blDsoPscG2dNBDqHxDqBCRnkomEQQphKbfnNwwN/axQIITQBNwGXxBg3Z52PDtoMar/J2wgcA5wP/HMI4bQsk1J6IYQStXnwf1KbC88GvhBCeGmmiWm/YoxrYoxD7+w3Edg5JHbQfY6NszqB8UNiE6itBdIoE0JYAPwU2AycG2OsZJyS0vs48GiM8TtZJ6JD0g1sjjF+OsbYE2P8KXA78NaM81J65wKnxBi/MfBvuAxYBrw747x0cIa1z7Fx1pNAGBIL7Ll8QzkXQngZ8HPgu8DbYoxD32Er384D/jSEsDWEsJXa1cpbQwiXZZyXDkwEJg78BmGXRqAho3x04I4Bxg2J9Q780ejzFNAcQjhmUOyg+xy3o9N9wLgQwkXUtk96F7VPnX4306x0QAa2L/sP4LMxxquyzkcHLsZ4wuCvQwgrqa2rdDu60eV7wBbgyoE3Pa8EzgH+W6ZZ6UB8j9quRH8O3EBtW8FzqG1lplEmxrgjhHAntX/TxdTWNr8DePPBPJ5XnMe4GGM3sAh4O7Vf718EvCXG2JFpYjpQ76H2KfCPhxDaB/351P7+Q0nDZ+A3Pa8DXgysp/bp/YvdG3/0iDE+Tm3Lskuo7fn7JeDPYoz/mWliOhSLgWZgNbWlU5fGGH9+MA/UUK0OXUMtSZIkaSivOEuSJEkp2DhLkiRJKdg4S5IkSSnYOEuSJEkp2DhLkiRJKdg4S5IkSSnYOEuSJEkpeOdASSqgEMKfAR8FjgMepXYDh5htVpI0unnFWZIKJoRwKXAl8GHgJKAP+FqmSUlSAXjnQEkqkBDCicB/AW+OMd4zEHszsAyYHWNcn2V+kjSaecVZkorlw8Bju5rmARsG/p6ZQT6SVBg2zpJUECGERuBc4PYh3xo/8Pe2kc1IkorFpRqSVBAhhJOBR4AuoH/QtxqprXOeEmP0RV+SDpK7akhScZwAVKh9IHBw4/xPQKtNsyQdGhtnSSqOacC2wdvOhRCagVdT25pOknQIXOMsScWxAZgYQmgZFPsA0An8WzYpSVJxeMVZkorjPqAX+EQI4RrgjcDfAW+JMXZnmpkkFYAfDpSkAgkhvBX4PHA08BBwaYzxwWyzkqRisHGWJEmSUnCNsyRJkpSCjbMkSZKUgo2zJEmSlIKNsyRJkpSCjbMkSZKUgo2zJEmSlIKNsyRJkpSCjbMkSZKUgo2zJEmSlML/BwhLH/6KOnS9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6665320d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(12,10))\n",
    "# suppose our chain was this\n",
    "bimodal_chain = np.append(norm(2,.8).rvs(5000),norm(8,.5).rvs(5000))\n",
    "plt.hist(bimodal_chain,bins=100,normed=True)\n",
    "plt.xlabel('$\\\\theta$')\n",
    "plt.title('Posterior')\n",
    "plt.axhline(.025/2,c='k',linestyle=':',label=\"95% HPD Cutoff for Bi-Modal Chain\")\n",
    "plt.axvline(bimodal_chain.mean(),c='r',label=\"Mean\")\n",
    "plt.axvline(np.percentile(bimodal_chain,2.5),c='g',linestyle='--',alpha=.5,label=\"95% CI\")\n",
    "plt.axvline(np.percentile(bimodal_chain,97.5),c='g',linestyle='--',alpha=.5)\n",
    "\n",
    "plt.axvline(.05,c='k',linestyle='--',label=\"Approximate 95% HPD CI\")\n",
    "plt.axvline(4.05,c='k',linestyle='--')\n",
    "plt.axvline(6.7,c='k',linestyle='--')\n",
    "plt.axvline(9.39,c='k',linestyle='--')\n",
    "\n",
    "plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Chain and MAP\n",
    "\n",
    "How do the moments of our chain (e.g. mean and median) relate to MAP?  To answer that question, we need to solve for MAP given the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:06:17.874461",
     "start_time": "2016-03-02T08:06:17.868861"
    }
   },
   "outputs": [],
   "source": [
    "def log_posterior(mu,data,sigma,mu_mu_prior,mu_sigma_prior):\n",
    "    ln_prior = norm(mu_mu_prior,mu_sigma_prior).logpdf(mu)\n",
    "    ln_like = norm(mu,sigma).logpdf(data).sum()\n",
    "    return -1*(ln_like + ln_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:06:19.282756",
     "start_time": "2016-03-02T08:06:19.220913"
    }
   },
   "outputs": [],
   "source": [
    "startvals = 9\n",
    "from scipy.optimize import minimize\n",
    "res = minimize(log_posterior,startvals, \n",
    "               args=(data,sigma,mu_prior,sigma_prior),options={'disp': True},\n",
    "               method='L-BFGS-B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-03-02T08:06:21.061002",
     "start_time": "2016-03-02T08:06:21.053796"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The MAP is:  10.16031563320772\n",
      "The mean is:  10.180621628544616\n",
      "The median is:  10.196002242838691\n"
     ]
    }
   ],
   "source": [
    "print(\"The MAP is: \", res.x[0])\n",
    "print(\"The mean is: \", chain[100:].mean())\n",
    "print(\"The median is: \", np.percentile(chain[100:],50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For well behaved symmetric posteriors, the mean of the posterior will usually be very close to MAP, but this doesn't necessarily have to be the case.  Note that in the preceding figure (the bimodal posterior), MAP is approximately 8 whereas the mean is approximately 5."
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.6.4"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}