{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Generation of plots for tutorial on good talks\n", "\n", "*This document was created from a Jupyter notebook. You can download the notebook [here](l04_good_plots_for_slides.ipynb).*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", " \n", "
\n", " \n", " BokehJS successfully loaded.\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import warnings\n", "warnings.simplefilter(action='ignore', category=UserWarning)\n", "\n", "import numpy as np\n", "import scipy.stats as st\n", "\n", "import seaborn as sns\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "%config InlineBackend.figure_formats = {'png', 'retina'}\n", "\n", "# Markov chain Monte Carlo package\n", "import emcee\n", "\n", "import bokeh.plotting\n", "import bokeh.io\n", "import bokeh.charts\n", "import bokeh.models\n", "bokeh.io.output_notebook()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data set and necessary calculations\n", "\n", "The data set consists of reversals of worms from three strains, wild type, Channelrhodopsin in the ASH neuron, and Channelrhodopsin in the AVA neuron. We got, respectively, 0/35, 9/35, and 33/36 reversals." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Store data set\n", "r = np.array([0, 9, 33])\n", "n = np.array([35, 35, 36])\n", "strains = ('wild type', 'ASH', 'AVA')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The posterior distribution is\n", "\n", "\\begin{align}\n", "P(p_\\mathrm{rev}) = \\frac{(n+1)!}{r!(n-r)!}\\,p_\\mathrm{rev}^r\\,(1-p_\\mathrm{rev})^{n-r}.\n", "\\end{align}\n", "\n", "We can write a function to compute this, since we'll need it for plotting." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def posterior(p_rev, r, n):\n", " \"\"\"Log posterior of reversal measurements.\"\"\"\n", " return (n + 1) * st.binom.pmf(r, n, p_rev)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Later, we will plot 95% confidence intervals that contain the highest probability density, or HPD. This is easily computed analytically for the case where $r = 0$, since, in this case,\n", "\n", "\\begin{align}\n", "P(p_\\mathrm{rev}) = (n+1)(1-p_\\mathrm{rev})^n\n", "\\end{align}\n", "\n", "This decreases monotonically with $p_\\mathrm{rev}$. Then, the minimum of the HPD interval is $0$. We can find the maximum by evaluating\n", "\n", "\\begin{align}\n", "\\int_0^{p_\\mathrm{rev}^\\mathrm{max}} \\mathrm{d}p_\\mathrm{rev}\\,P(p_\\mathrm{rev}) = \n", "1-(1-p_\\mathrm{rev}^\\mathrm{max})^{n+1} = f_\\mathrm{HPD},\n", "\\end{align}\n", "\n", "where $f_\\mathrm{HPD}$ is the fraction of the total space of values that the HPD covers (in our case, 0.95). This is solved to give\n", "\n", "\\begin{align}\n", "p_\\mathrm{rev}^\\mathrm{max} = 1 - (1-f_\\mathrm{HPD})^{(n+1)^{-1}}.\n", "\\end{align}\n", "\n", "So, we can compute the HPD for the wild type." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Intialize PDs, filling in WT value\n", "hpds = [(0, 1 - (1 - 0.95)**(1 / (n[0] + 1))), None, None]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will compute the HPDs using Markov chain Monte Carlo. First, we need a function for the log posterior." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def log_posterior(p_rev, r, n):\n", " \"\"\"Unnormalized log posterior, used for MCMC calculations.\"\"\"\n", " # Get scalar probability (emcee expects array argument)\n", " p = p_rev[0]\n", "\n", " if p < 0 or p > 1:\n", " return -np.inf\n", " return r * np.log(p) + (n - r) * np.log(1-p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll write a function to set up and run an MCMC calculation." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def sample_posterior_mcmc(r, n, n_walkers=50, n_burn=5000, n_steps=5000):\n", " \"\"\"Get samples of p_rev from the posterior.\"\"\"\n", "\n", " # Specify starting points of walkers\n", " p0 = np.empty((n_walkers, 1))\n", " p0[:,0] = np.random.uniform(0, 1, n_walkers)\n", "\n", " # Instantiate sampler\n", " sampler = emcee.EnsembleSampler(n_walkers, 1, log_posterior, args=(r, n))\n", " \n", " # Do burn-in and sampling\n", " pos, prob, state = sampler.run_mcmc(p0, n_burn, storechain=False)\n", " _ = sampler.run_mcmc(pos, n_steps)\n", " \n", " return sampler" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we need a function to compute the HPD from the samples." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def hpd(sampler, mass_frac):\n", " \"\"\"\n", " Returns highest probability density region given by\n", " a set of samples.\n", " \"\"\"\n", " trace = sampler.flatchain\n", " \n", " # Get sorted list\n", " d = np.sort(np.copy(trace))\n", "\n", " # Number of total samples taken\n", " n = len(trace)\n", " \n", " # Get number of samples that should be included in HPD\n", " n_samples = np.floor(mass_frac * n).astype(int)\n", " \n", " # Get width (in units of data) of all intervals with n_samples samples\n", " int_width = d[n_samples:] - d[:n-n_samples]\n", " \n", " # Pick out minimal interval\n", " min_int = np.argmin(int_width)\n", " \n", " # Return interval\n", " return np.sort(np.array([d[min_int], d[min_int+n_samples]]).flatten())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will compute the mode and HPD for the ASH and AVA strains." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for i in (1, 2):\n", " sampler = sample_posterior_mcmc(r[i], n[i])\n", " hpds[i] = hpd(sampler, 0.95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that all of our calculations are done, we can make our plots." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ugly plots\n", "\n", "We'll begin by making the \"ugly\" plots that are more or less defaults, but not really useful for presentations. We'll start with a bar chart." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAw4AAAIRCAYAAADjteq3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAWJQAAFiUBSVIk8AAAIABJREFUeJzs3XlY1WX+//HXBxTEBIVwQTS3TFxSzIRxJjPTTC0VW9QW\nddS0NMrCMmvSb2aN9i2XNFtsQ00rNZepEbepoawElNRMwfQbKCqaOwixnt8f/WCGBD8INwcOPh/X\nNdc1574/9/1+c3Wd5rzms1kOh8MhAAAAALgEt8puAAAAAEDVR3AAAAAAYIvgAAAAAMAWwQEAAACA\nLYIDAAAAAFsEBwAAAAC2CA4AAAAAbBEcAAAAANgiOAAAAACwRXAAAAAAYIvgAJRTmzZt1KZNm8pu\nA7ji8N0DKgffvSsXwQEAAACALYIDAAAAAFvVJjgcP35cN954o5YsWVLqNefOndOLL76oW2+9VcHB\nwbrrrru0fv36CuwSAAAAcE01KrsBEzIyMvTYY4/pwoULpV6TmZmpUaNGKTExUX379lVAQIA2bdqk\niIgInTlzRg888EAFdgwAAAC4Fpc/43DkyBE9+OCD2r1792WtW7x4sfbt26fnn39es2fP1lNPPaW1\na9eqdevWeu2113T69OkK6hgAAABwPS4dHCIjIzVw4EDt379f3bp1u6y1H3/8sa6++moNGzascKx2\n7dp65JFHlJmZqS+++MJ0uwAAAIDLcungsGTJEjVp0kTLli3TwIED5XA4SrXu8OHDhfdEWJZVZC40\nNFSSFBcXZ7xfAAAAwFW59D0OM2bM0J///GdZlqX/+7//K/W6Q4cOSZKuueaai+b8/f3l6emppKQk\nU20CAAAALs+lg8Nf/vKXMq07e/asJMnb27vY+Tp16igtLa3MfQEAAADVjUtfqlRWOTk5kiQPD49i\n5z08PJSVleXMlgAAAIAqzaXPOJRVrVq1JP0nQPxRdna2vLy8ylXjcl/F7u3trdq1a6tZs2blqovK\nM3z48MpuAbgi8d0DKgffPdeUnJys48ePX9aaxMRESVdocKhbt64klXg5Unp6uvz9/Z3ZktLS0pSW\nlnbZ/yBRdcTGxlZ2C8AVie8eUDn47l15rsjg0Lx5c0lSSkrKRXO//vqrsrKy1KJFi3LVKEhmpXHz\nzTfr+PHj8vb2Vtu2bctVFwAAACjJvn37lJaWpoYNG+rrr7++rLVXZHAICAhQ48aNFR8ff9FcTEyM\nJKlz585O66dZs2Y6fvy42rZtq6VLlzqtLgAAAK4sw4cPV2xsbJkuj78ib46WpIEDB+rYsWP66KOP\nCsfS09P19ttvy8vLSwMHDqzE7gAAAICq5Yo447BgwQJZlqXw8PDCsYceekhRUVF6+eWXFRsbq6ZN\nm2rTpk1KSUnR1KlT5evrW4kdAwAAAFVLtTrj8Me3QBdYuHCh3nzzzSJjderU0ccff6y7775bO3bs\n0PLly+Xj46M5c+bo/vvvd0a7AAAAgMuoNmccBg8erMGDBxc7l5CQUOy4n5+fXnrppYpsCwAAAKgW\nqtUZBwAAAAAVg+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAA\nALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAA\nALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAA\nALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs1ajsBgAAAEx6ZORI\nnTl6rLLbAMrNt3GA3l68uLLbKERwAAAA1cqZo8c0MiCwstsAym3x0SOV3UIRXKoEAAAAwBbBAQAA\nAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAA\nAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAA\nAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAA\nAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAA\nAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAA\nAIAtlw4OeXl5ioyM1B133KFOnTqpd+/eevPNN5Wbm1uq9YmJiRo/frxCQkLUsWNHDRw4UCtWrKjg\nrgEAAADX49LBYfr06Zo1a5b8/Pw0cuRINWrUSPPnz9ekSZNs1yYkJGjYsGH65ptv1KNHD91///3K\nzMzUtGnTNHv2bCd0DwAAALiOGpXdQFnFx8drxYoV6tevn+bOnVs4PmXKFK1bt07R0dHq0aNHievn\nzZun3377TQsXLtStt94qSZo4caIGDx6sDz74QMOGDVNgYGCF/x0AAACAK3DZMw7Lli2TZVkKDw8v\nMh4RESFJWrly5SXX79mzRz4+PoWhQZK8vLx05513Kj8/X7t37zbfNAAAAOCiXDY47NixQ76+vmrV\nqlWR8QYNGqh58+aKi4u75Pp69erpwoULSktLKzKempoqSfLz8zPbMAAAAODCXDI4ZGdnKzU1Vddc\nc02x84GBgTp//rzOnDlT4h7Dhg1TXl6eJk2apEOHDunChQtatWqV1q5dq/bt2yskJKSi2gcAAABc\njkve43Du3DlJkre3d7HzBePp6eny9fUt9pgHH3xQ7u7uevnll9WnT5/C8b/85S+aO3euLMsy3DUA\nAADgulwyOBQ8btXDw6PY+YLxrKysEvfYuXOnFi1aJA8PDw0YMEDe3t767rvv9P333+v111/X1KlT\ny9VjmzZtyrUeAAAAqCixsbGl/r2amJgoyUWDg6enpyQpJyen2Pns7GxJv9/sXJz09HQ9/PDDkqS1\na9cWXvKUm5urSZMmadmyZbr22mt13333mW4dAAAAcEkuGRy8vb3l5uZ20Y3NBQrGS7qU6csvv9S5\nc+cUHh5e5D6JGjVqaNq0adq4caPWrFlTruBQkMxKY/jw4YqNjS1zLQAAAOByhISEaOnSpZe1xiVv\njq5Zs6YaN26slJSUYudTUlLk5+cnHx+fYudTU1NlWZZatmx50dzVV18tX19fHT161GjPAAAAgCtz\nyeAgSV26dNHJkyeVnJxcZPzEiRNKSkpScHBwiWuvvvpqORwOJSUlXTR3/vx5nT17VvXr1zfdMgAA\nAOCyXDY4hIWFyeFwaM6cOXI4HIXjs2fPlmVZGjJkSIlre/bsKS8vL3300Uc6fPhw4Xh+fr5mzpwp\nSbrzzjsrrnkAAADAxbjkPQ6S1K1bN/Xv319RUVEaOnSoQkNDFR8fr/j4ePXt21c9evQoPHbBggVF\n3jLt5+enqVOnaurUqQoLC9Ptt98uHx8fbdu2TYmJiQoJCdGIESMq608DAAAAqhyXDQ6S9Oqrr6p1\n69Zas2aNlixZooCAAE2cOFFjxowpctzChQvl5uZWGBwk6a677lKTJk20aNEibdmyRb/99puaNm2q\nJ554QqNHj1bNmjWd/ecAAAAAVZZLBwd3d3eNHz9e48ePv+RxCQkJxY6HhITwhmgAAACgFFz2HgcA\nAAAAzkNwAAAAAGCL4AAAAADAFsEBAAAAgC2CAwAAAABbBAcAAAAAtggOAAAAAGwRHAAAAADYIjgA\nAAAAsEVwAAAAAGCL4AAAAADAFsEBAAAAgC2CAwAAAABbBAcAAAAAtggOAAAAAGwRHAAAAADYIjgA\nAAAAsEVwAAAAAGCL4AAAAADAFsEBAAAAgC2CAwAAAABbBAcAAAAAtggOAAAAAGwRHAAAAADYIjgA\nAAAAsEVwAAAAAGCL4AAAAADAFsEBAAAAgC2CAwAAAABbBAcAAAAAtggOAAAAAGwRHAAAAADYqvDg\ncOLECf3444/KyMio6FIAAAAAKoix4HDgwAE9++yz2r59e+HYa6+9pp49e2rIkCHq3r27Pv30U1Pl\nAAAAADiRkeBw8OBBDRkyRGvXrlViYqIk6fvvv9d7770nSWrfvr3y8vL0wgsvaNu2bSZKAgAAAHAi\nI8Fh0aJFysjI0JgxYzRgwABJ0qpVq2RZlp5++mmtWrVKn3zyidzd3RUZGWmiJAAAAAAnqmFik5iY\nGLVu3VpPPfWUJCkvL09ff/213NzcFBYWJkkKCgpSly5dtGvXLhMlAQAAADiRkTMOp06dUqtWrQo/\n79y5U2lpaWrbtq3q1atXOF63bl2lpaWZKAkAAADAiYwEB39/f506darwc3R0tCzL0p///Ocix/3y\nyy9FggQAAAAA12AkOAQFBSk+Pl7btm1TUlKS1qxZI0nq3bt34TFLly7Vzz//rBtuuMFESQAAAABO\nZOQeh3Hjxun777/XqFGjJEkOh0PdunVTx44dJUmDBg3S/v37VatWLY0bN85ESQAAAABOZOSMQ+fO\nnRUZGamQkBC1bNlSQ4cO1YIFCwrna9SooXbt2mnJkiXq0KGDiZIAAAAAnMjIGQdJCg4O1uLFi4ud\nW7x4serUqWOqFAAAAAAnM/bm6EshNAAAAACurUxnHKZOnVrmgpZl6cUXXyzzegAAAADOV6bgsHLl\nyjIXJDgAAAAArqdMwWHmzJmm+wAAAABQhZUpOAwePNh0HwAAAACqMKfcHP3fEhMTnV0SAAAAQDkZ\nexzrTz/9pE8//VRHjx5VTk6OHA5H4ZzD4VBWVpZOnjyp1NRU7d2711RZAAAAAE5gJDjs3r1bDz74\nYJHAYFlWkfBgWZYk6brrrjNREgAAAIATGQkO7777rrKzs3X77bfrrrvuUnR0tD755BO9/fbbys/P\n19atW/XJJ5+oRYsW5XoiEwAAAIDKYeQehx9++EENGjTQa6+9ph49euiOO+5Qfn6+cnJy1LNnT02d\nOlUvvviiDhw4oMjISBMlAQAAADiRkeBw9uxZtWvXTjVr1pT0n8uRfvrpp8Jj7r77bjVp0kTr1683\nURIAAACAExkJDrVr15a7u3vhZ29vb9WrV08HDx4sclzbtm2VnJxsoiQAAAAAJzISHFq0aKG9e/cq\nPz+/yNiePXuKHJeRkWGiHAAAAAAnMxIcbrvtNh07dkyTJk3S4cOHJUkhISE6duyYVq9eLen3Jy/F\nxMSoadOmJkoCAAAAcCIjweHBBx9Uu3btFBUVpZdffrlwrFatWvrb3/6m7t27a9iwYcrLy9O9995r\noiQAAAAAJzISHGrVqqWPP/5YTz/9tLp37y5Jql+/vt566y0FBgbq119/lYeHhx566CE98MADJkoC\nAAAAcCJjb4729PTUmDFjioz96U9/0pYtW3T69GnVq1dPbm5GcgoAAAAAJ6vwX/InTpzQkSNH9Ntv\nv1V0KQAAAAAVxFhwOHDggJ599llt3769cOy1115Tz549NWTIEHXv3l2ffvqpqXIAAAAAnMhIcDh4\n8KCGDBmitWvXKjExUZL0/fff67333pMktW/fXnl5eXrhhRe0bds2EyUBAAAAOJGR4LBo0SJlZGRo\nzJgxGjBggCRp1apVsixLTz/9tFatWqVPPvlE7u7uioyMNFESAAAAgBMZuTk6JiZGrVu31lNPPSVJ\nysvL09dffy03NzeFhYVJkoKCgtSlSxft2rXLREkAAAAATmTkjMOpU6fUqlWrws87d+5UWlqa2rZt\nq3r16hWO161bV2lpaSZKAgAAAHAiI8HB399fp06dKvwcHR0ty7L05z//uchxv/zyS5EgAQAAAMA1\nGAkOQUFBio+P17Zt25SUlKQ1a9ZIknr37l14zNKlS/Xzzz/rhhtuMFESAAAAgBMZucdh3Lhx+v77\n7zVq1ChJksPhULdu3dSxY0dJ0qBBg7R//37VqlVL48aNM1ESAAAAgBMZOePQuXNnRUZGKiQkRC1b\nttTQoUO1YMGCwvkaNWqoXbt2WrJkiTp06GCiJAAAAAAnMnLGYf/+/br++uu1ePHiYucXL16sOnXq\nmCgFAAAAoBIYOeMwYcIE3XnnnSXOExoAAAAA12YkOBw/frzI41gBAAAAVC9GgkPLli114MAB5eXl\nmdiu1PLy8hQZGak77rhDnTp1Uu/evfXmm28qNze3VOuzs7P1xhtv6Pbbb1fHjh112223aebMmbxr\nAgAAAPgDI8HhlVdeUUZGhh544AGtXbtWCQkJSk1N1fHjx4v9jynTp0/XrFmz5Ofnp5EjR6pRo0aa\nP3++Jk2aZLs2NzdXY8aM0cKFC9WwYUONGDFCAQEBWrx4scaOHVvq8AEAAABcCYzcHD1q1ChlZWVp\n165d2rVr1yWPtSxLe/fuLXfN+Ph4rVixQv369dPcuXMLx6dMmaJ169YpOjpaPXr0KHH94sWLFRcX\np7FjxxYJGjNmzNDy5cv1z3/+U4MGDSp3nwAAAEB1YCQ4eHl5ycvLy6lvhV62bJksy1J4eHiR8YiI\nCK1bt04rV668ZHBYtmyZmjRpoieeeKLI+OjRo5WRkaFatWpVSN8AAACAKzISHL788ksT21yWHTt2\nyNfX96Kbshs0aKDmzZsrLi6uxLUHDx7U0aNHNXLkSLm7uxeZCwwM1MyZMyukZwAAAMBVGbnHwdmy\ns7OVmpqqa665ptj5wMBAnT9/XmfOnCl2fv/+/bIsS9dee62io6N13333KTg4WN27d9crr7yizMzM\nimwfAAAAcDlGg8Pp06e1aNEijR07VgMGDNArr7wiSXrrrbeMnpU4d+6cJMnb27vY+YLx9PT0YudP\nnDghh8Ohf/3rX3r44YdVt25d3Xfffapfv74+/PBDjR071ulPiAIAAACqMiOXKklSdHS0nn76aaWl\npcnhcMiyLLVt21aStH79es2fP18jRozQs88+W+5aBU888vDwKHa+YDwrK6vY+YIzCtHR0XrppZd0\nzz33SJIcDoeefPJJbdy4UcuXL9fw4cPL3SsAAABQHRgJDgkJCXrsscdkWZZGjRqlm266SaNHjy6c\nf+CBBzR37lwtWbJEXbt2Ve/evctVz9PTU5KUk5NT7Hx2drak32/aLo6b2+8nWtq2bVsYGqTfn/g0\nefJkbdiwQVFRUeUKDm3atCnzWgAAAKAixcbGlvr3amJioiRDwaHgpWvvv/++unXrdtH8sGHD1LFj\nR91zzz1aunRpuYODt7e33NzcSnxRW8F4SZcy1alTR5LUvn37i+YaN24sHx8fHTp0qFw9AgAAANWJ\nkeAQFxenTp06FRsaCrRr105dunTRwYMHy12vZs2aaty4sVJSUoqdT0lJkZ+fn3x8fIqdb968uaSS\nz1jk5uaqbt265eqxIJmVxvDhwxUbG1uuegAAAEBphYSEaOnSpZe1xsjN0RcuXJC/v7/tcd7e3iWe\nJbhcXbp00cmTJ5WcnFxk/MSJE0pKSlJwcHCJazt27KiaNWsqNjZWDoejyNzBgweVkZGhoKAgI30C\nAAAA1YGR4BAQEKC9e/de9CP8v+Xl5Wnv3r1q1KiRiZIKCwuTw+HQnDlzitSdPXu2LMvSkCFDSlxb\np04d9e/fX8eOHdM777xTOJ6bm6tXX31VlmXp7rvvNtInAAAAUB0YuVSpd+/e+uCDDzRv3jw9+eST\nxR4zf/58HT9+XH/9619NlFS3bt3Uv39/RUVFaejQoQoNDVV8fLzi4+PVt2/fIm+NXrBgwUVvmX7m\nmWe0c+dOvf7664U3h2zbtk0JCQnq37+/brnlFiN9AgAAANWBkeDw8MMPa+PGjVq0aJG+/fZbde3a\nVZJ06NAhvf322/r666/1ww8/qH79+ho7dqyJkpKkV199Va1bt9aaNWu0ZMkSBQQEaOLEiRozZkyR\n4xYuXCg3N7ciwcHPz08rVqzQwoULtXnzZu3YsUOBgYGaPHmysXADAAAAVBdGgoOPj4+WLl2qSZMm\nKT4+Xnv27JEk7dy5Uzt37pQkBQUFac6cOfLz8zNRUpLk7u6u8ePHa/z48Zc8LiEhodjxunXr6rnn\nntNzzz1nrCcAAACgOjL2AriAgAAtX75cu3btUkxMjI4dO6a8vDw1aNBAN954o/70pz+ZKgUAAADA\nyYy9AK7gKUSdOnVSp06dTGwLAAAAoIow8lSlsLAwDRgwQO+++65SU1NNbAkAAACgCjESHNq1a6ef\nf/5Zc+bM0a233qqRI0fqs88+U3p6uontAQAAAFQyI5cqrV69WsnJyfriiy+0fv16xcTEKDY2VjNm\nzFDPnj01cOBA3XzzzXJ3dzdRDgAAAICTGbs5ulmzZnr00Uf16KOPKjExsTBEREVFacOGDapXr576\n9++vAQMGXPKtzgAAAACqHiOXKv1RmzZtNGnSJP3rX//SihUrNGbMGHl4eGj58uW6//77K6IkAAAA\ngApUIcGhQFJSkr777jvFxcXpxIkTcjgcqlevXkWWBAAAAFABjF2qVODw4cOKiorS+vXrlZiYKIfD\noVq1aqlfv34aNGiQbrrpJtMlAQAAAFQwI8Hh2LFjhWHhp59+ksPhkJubm0JCQjRo0CD16dNHderU\nMVEKAAAAQCUwEhx69uwpy7LkcDjUunVrDRw4UAMHDlTDhg1NbA8AAACgkhkJDv7+/howYIAGDRpU\n+AZpAAAAANWHkeDw9ddfy82tQu+zBgAAAFCJjASHgtCQm5urTZs2KTY2Vqmpqbrhhhs0btw4rVy5\nUtdffz1nIwAAAAAXZeypSnv27NETTzyhI0eOyOFwyLIs+fj4SJKWLVum/fv365lnntHIkSNNlQQA\nAADgJEauL0pJSdHo0aN15MgR9enTRzNmzJDD4Sicv+WWW1SjRg3NmjVLMTExJkoCAAAAcCIjwWHh\nwoVKS0vTrFmz9Prrr+vee+8tMv/EE09o4cKFcjgc+vDDD02UBAAAAOBERoLD1q1b1bZtWw0aNKjE\nY7p3767g4GDt27fPREkAAAAATmQkOJw9e1ZNmjSxPc7f319nzpwxURIAAACAExkJDvXr19eBAwds\nj/v555/l7+9voiQAAAAAJzISHLp3765ffvlFy5cvL/GY5cuXKzk5WTfddJOJkgAAAACcyMjjWCdM\nmKCNGzdqxowZ2rZtm0JDQyVJp06d0ueff67o6Gj985//lLe3t8aNG2eiJAAAAAAnMhIcGjZsqA8+\n+ECPP/64Nm3apM2bN0uSvvvuO3333XdyOBzy9/fXvHnzSnUvBAAAAICqxdgL4Nq1a6eoqCht3rxZ\n27ZtU2pqqvLy8tSgQQPdeOONuuOOO1SrVi1T5QAAAAA4kZHg8M4776hVq1bq3bu3+vfvr/79+5vY\nFgAAAEAVYSQ4fPDBB2rQoIF69+5tYjsAAAAAVYyRpyplZ2erWbNmJrYCAAAAUAUZCQ59+/bVt99+\nq4SEBBPbAQAAAKhijFyqdMstt+iHH37QPffcoy5duigoKEh169aVm1vxueSRRx4xURYAAACAkxgJ\nDhMnTpRlWXI4HIqJiVFMTIwkybKsIsc5HA5ZlkVwAAAAAFyMkeDw6KOPXhQSAAAAAFQfRoLDY489\nZmIbAAAAAFWUkZujAQAAAFRvBAcAAAAAtggOAAAAAGwRHAAAAADYIjgAAAAAsFWm4DB37lxt2LDB\ndC8AAAAAqqgyBYfly5dr48aNhZ979eqlV155xVhTAAAAAKqWMgWH7Oxs/frrr4Wfjxw5olOnThlr\nCgAAAEDVUqYXwDVv3lw7duzQ6NGj1bBhQ0nSzp079eyzz9qutSxLf//738tSFgAAAEAlKVNwePTR\nRxUREaHvvvtO0u9h4NChQzp06JDtWoIDAAAA4HrKFBz69OmjqKgo/fTTT8rKytIzzzyjG264QcOG\nDTPdHwAAAIAqoEzBQZKaNm2qpk2bSpKeeeYZNWnSRAMHDjTWGAAAAICqo8zB4b8lJCSY2AYAAABA\nFWUkOBT4+eeftXjxYsXFxenXX3+Vh4eH/P39FRoaqnvvvVdBQUEmywEAAABwEmPB4bPPPtMLL7yg\nnJycwrGMjAydPXtWBw4c0IoVK/Q///M/uueee0yVBAAAAOAkZXqPwx/t3r1b06ZNk5ubm8LDwxUV\nFaUff/xRu3bt0hdffKEJEybIzc1NL7zwgn766ScTJQEAAAA4kZHg8N577yk/P18LFixQeHi4WrRo\noZo1a8rT01PXXnutHn/8cS1YsEC5ubmKjIw0URIAAACAExkJDtu3b1enTp108803l3jMzTffrODg\nYMXGxpooCQAAAMCJjASH8+fPq1GjRrbHNWrUSGfOnDFREgAAAIATGQkODRo0UGJiou1xCQkJ8vf3\nN1ESAAAAgBMZCQ7du3dXUlKS3nvvvRKPWbRokZKSknTTTTeZKAkAAADAiYw8jvWRRx7RP//5T82e\nPVsxMTHq27evAgMDJUkpKSnasGGDvv32W/n4+Ojhhx82URIAAACAExkJDgEBAXr//fcVHh6ub775\nRlu3bi0y73A41KBBA73++uuFgQIAAACA6zD2ArhOnTppy5YtioqKUlxcnE6cOFEYGLp27ap+/fqp\nVq1apsoBAAAAcCJjwUGSPD09FRYWprCwMJPbAgAAAKhkRm6OBgAAAFC9ERwAAAAA2CI4AAAAALBF\ncAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAW8be47B7924tXrxYP//8szIzM5Wfn1/scZZlacuW\nLabKAgAAAHACI8Fh+/btGjVqlHJzc+VwOC55rGVZJkoCAAAAcCIjwWHhwoXKycnRoEGDdP/998vf\n3181ahh9KTUAAACASmTk1/2uXbt07bXX6pVXXjGxHQAAAIAqxsjN0W5ubmrZsqWJrQAAAABUQUaC\nQ4cOHZSQkGB7fwMAAAAA12QkODz22GM6cuSIFixYYGI7AAAAAFWMkXscEhIS1L17d7311ltau3at\nOnToIB8fn2KfoGRZll588UUTZQEAAAA4iZHgMGPGDFmWJYfDoaNHj+ro0aMlHktwAAAAAFyPkeAw\nc+ZME9uDNwTUAAAgAElEQVRctry8PC1dulQrV65USkqK6tevr7vuukvjxo277MfB5ufna9iwYdq9\ne7cSEhIqqGMAAADANRkJDoMHDzaxzWWbPn26VqxYoa5du6pXr16Kj4/X/PnzlZiYqNdff/2y9oqM\njNTu3bt5QR0AAABQDONvacvOztaePXt08uRJeXh46Oqrr1bbtm2NvxAuPj5eK1asUL9+/TR37tzC\n8SlTpmjdunWKjo5Wjx49SrVXcnKy5s+fT2gAAAAASmDs13xubq7mz5+vjz76SJmZmUXmvL29NXTo\nUD3++OOqWbOmkXrLli2TZVkKDw8vMh4REaF169Zp5cqVpQ4Ozz//vBo2bCjLspScnGykPwAAAKA6\nMRIc8vLyNH78eG3dulVubm7q1KmTAgMDlZ+fr8OHD2vv3r167733lJCQoHfffddESe3YsUO+vr5q\n1apVkfEGDRqoefPmiouLK9U+H3/8sbZv367Fixfr73//u5HeAAAAgOrGSHD49NNP9c033+j666/X\nnDlz1LRp0yLzhw4dUkREhLZu3arPPvtMd999d7nqZWdnKzU1VcHBwcXOBwYGKikpSWfOnJGvr2+J\n+xw7dkyvvfaa7r33XoWEhJSrJwAAAKA6M/ICuDVr1uiqq67SO++8c1FokKRrrrlGixYtUu3atbVq\n1apy1zt37pyk3y+BKk7BeHp6+iX3mTZtmq666ipNnjy53D0BAAAA1ZmR4HDgwAGFhITIz8+vxGP8\n/PwUEhKigwcPlrtebm6uJMnDw6PY+YLxrKysEvdYu3attm7dqmnTpqlOnTrl7gkAAACozow/VclO\nTk5Ouffw9PS85F7Z2dmSJC8vr2LnT506pZkzZ+q2225T7969y91Pcdq0aVMh+wIAAADlFRsbW+rf\nq4mJiZIMnXFo0aKF4uLiCi8hKs7Zs2cVFxenli1blruet7e33NzclJaWVux8wXhJlzJNnz5dDodD\n06ZNK3cvAAAAwJXAyBmHu+66Sy+99JImTJigOXPmqGHDhkXmU1NTFRERoQsXLmjQoEHlrlezZk01\nbtxYKSkpxc6npKTIz89PPj4+xc5v2rRJlmXppptuumjOsiwFBQUpMDBQ//rXv8rcY0EyK43hw4cr\nNja2zLUAAACAyxESEqKlS5de1hojweG+++7Txo0bFRcXp169eik4OFiBgYGSfv8Rv2vXLuXm5qpr\n1666//77TZRUly5d9I9//EPJyclq1qxZ4fiJEyeUlJSkXr16lbj2j+9+KPDJJ5/o1KlTeuyxx0o8\nWwEAAABciYwEB3d3d73//vuaPXu2Pv30U23fvl3bt28vnPfy8tIDDzygiIgIY2+QDgsL07p16zRn\nzhzNmzev8K3Ps2fPlmVZGjJkSIlrSwoOW7Zs0alTp/Too48a6REAAACoLozdHO3h4aFnn31WERER\n+vHHH3XixAlJv7+QrUOHDqpVq5apUpKkbt26qX///oqKitLQoUMVGhqq+Ph4xcfHq2/fvkXeGr1g\nwYJi3zINAAAAoHSMP1XJ09NTN954o+lti/Xqq6+qdevWWrNmjZYsWaKAgABNnDhRY8aMKXLcwoUL\n5ebmVqrgUHDmAgAAAMB/lCk4fP7555KkW2+9VVdddVXh59IaMGBAWcpexN3dXePHj9f48eMveVxC\nQkKp9lu7dq2JtgAAAIBqp0zB4emnn5ZlWVq/fr1atGhR+Lm0TAUHAAAAAM5RpuAQFhYmy7IKnzxU\n8BkAAABA9VSm4DBr1qxLfgYAAABQvRh5c/TatWu1Y8cO2+O2bNmi119/3URJAAAAAE5kJDhMmTJF\nK1assD1u3bp1+vDDD02UBAAAAOBEZbpU6f3331dmZmaRsYSEBL3xxhslrklPT9c333xj/H0OAAAA\nACpemYLDb7/9pjfeeEOWZcnhcMiyLO3fv1+JiYm2a4cOHVqWkgAAAAAqUZmCw9ixY1WjRg3l5+fL\n4XBo/vz5atu2rfr06VPs8ZZlydPTU82aNVPPnj3L1TAAAAAA5ytTcPDw8NDDDz9c+HnVqlUKDQ21\nfREbAAAAANdUpuDwR19++aUkKScnR7/88ouuu+66wrnDhw9r7969uvnmm+Xl5WWiHAAAAAAnM/JU\nJUnavHmzbrrpJj3//PNFxnfs2KGJEyeqV69e+u6770yVAwAAAOBERoLDjh079Pjjj+vChQtq06ZN\nkbmgoCANHjxY58+f17hx47Rr1y4TJQEAAAA4kZHg8Pbbb8vNzU3vvvuuZsyYUWQuKChIM2fO1Lvv\nvqu8vDy99dZbJkoCAAAAcCIjwWHfvn3q2rWrunXrVuIx3bp1U5cuXUr1hmkAAAAAVYuR4HDhwgV5\ne3vbHufn56fs7GwTJQEAAAA4kZHg0Lx5c23fvv2it0n/t+zsbP3www9q0qSJiZIAAAAAnMhIcLjz\nzjt15swZPfnkkzp79uxF8+np6Zo8ebJOnjyp/v37mygJAAAAwImMvMdh+PDhioqK0r///W/16NFD\nwcHBCggIkCSlpqZq165dyszMVNu2bTV69GgTJQEAAAA4kZHg4OHhocjISM2bN0+rV69WTEzMRfND\nhgzR5MmTeQkcAAAA4IKMBAdJqlOnjp5//nlNnjxZe/bs0a+//qq8vDz5+/urffv2uuqqq0yVAgAA\nAOBkxoJDAQ8PD91www2mtwUAAABQiYwGhxMnTujYsWPKycmRw+EoHM/Pz1dWVpZOnjypr776SgsW\nLDBZFgAAAEAFMxIcsrOz9dRTT2nz5s0mtgMAAABQxRh5HOv777+vTZs2yd3dXe3bt1fjxo0lSaGh\noWrbtq3c3d3lcDjUokULzjYAAAAALshIcNiwYYPc3Ny0fPlyrVq1Sk8++aQkacqUKVq9erW++uor\nde7cWYcOHVKDBg1MlAQAAADgREaCw6FDh9SpUyd17NhRktSxY0c5HA7Fx8dLkurXr6958+bJsiy9\n9957JkoCAAAAcCIjwSE3N1cNGzYs/NykSRPVqFFD+/fvLxxr2LChbrjhBv3www8mSgIAAABwIiPB\noX79+jp16tR/NnVzU2BgYJHgIEl169bVmTNnTJQEAAAA4ERGgkPnzp0VHx+vPXv2FI61adNGe/bs\n0enTpyX9/kjWffv2yc/Pz0RJAAAAAE5kJDiMHDlSDodDDzzwQOFTk8LCwpSdna1HHnlEK1asUHh4\nuFJSUhQcHGyiJAAAAAAnMvIeh44dO+p///d/NWPGDB0+fFiSdOutt+qWW27Rv//9b/34449yOByq\nW7eunnjiCRMlAQAAADiRsTdH33nnnerTp49OnjxZOPbmm29q3bp12r17twICAjRo0KAiN1EDAAAA\ncA1GgsP06dPVokULjRgxovDlb9LvN0kPHjxYgwcPNlEGAAAAQCUxco/D559/rs8//9zEVgAAAACq\nICPBQZKuvvpqU1sBAAAAqGKMBIchQ4Zo69at+ve//21iOwAAAABVjJF7HBo3bqwmTZpo/PjxCgwM\nVFBQkOrWrSs3t4tziWVZevHFF02UBQAAAOAkRoLDSy+9VPjfU1JSlJKSUuKxBAcAAADA9RgJDjNn\nzjSxDQAAAIAqqkzBYcSIEbrxxhv1+OOPS5JCQ0NVu3Zt1atXz2hzAAAAAKqGMt0cvXv3bv3yyy+F\nn3v16sVZBwAAAKAaK1NwqFGjhvbv36+cnBxJksPhkMPhMNoYAAAAgKqjTJcqdejQQTExMQoNDZWv\nr68kacuWLerVq5ftWsuytGXLlrKUBQAAAFBJyhQcnn/+eT3yyCNKSUlRRkaGLMtSRkaGMjIybNda\nllWWkgAAAAAqUZmCw7XXXqstW7bo9OnTys7O1i233KI+ffrob3/7m+n+AAAAAFQB5Xocq5+fnySp\na9euat++vRo2bGikKQAAAABVi5H3OCxdutTENgAAAACqqDI9VQkAAADAlYXgAAAAAMAWwQEAAACA\nLYIDAAAAAFsEBwAAAAC2CA4AAAAAbBEcAAAAANgiOAAAAACwRXAAAAAAYIvgAAAAAMBWjcpuAACq\nq5HjRunYydTKbgMolwD/Rlq86MPKbgNAFUBwAIAKcuxkqhoPuK6y2wDK5ejn+yu7BQBVBJcqAQAA\nALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABsERwAAAAA2CI4AAAA\nALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABsuXRwyMvLU2RkpO64\n4w516tRJvXv31ptvvqnc3NxSrd+zZ48mTJig0NBQdejQQbfddptmz56tzMzMCu4cAAAAcC0uHRym\nT5+uWbNmyc/PTyNHjlSjRo00f/58TZo0yXbttm3bdN9992nr1q3q3r27RowYIV9fX7377rsaOXKk\nsrOznfAXAAAAAK6hRmU3UFbx8fFasWKF+vXrp7lz5xaOT5kyRevWrVN0dLR69OhR4vrp06fL4XDo\n448/VocOHQrHp02bppUrV2r58uX661//WpF/AgAAAOAyXPaMw7Jly2RZlsLDw4uMR0RESJJWrlxZ\n4tqDBw/ql19+Ue/evYuEBkl69NFH5XA49M0335hvGgAAAHBRLnvGYceOHfL19VWrVq2KjDdo0EDN\nmzdXXFxciWvr1Kmjp59+Wq1bt75ormbNmpKkCxcumG0YAAAAcGEuGRyys7OVmpqq4ODgYucDAwOV\nlJSkM2fOyNfX96L5hg0basyYMcWu3bx5syTpuuuuM9cwAAAA4OJc8lKlc+fOSZK8vb2LnS8YT09P\nv6x9T548qfnz58uyLN17773laxIAAACoRlzyjEPB41Y9PDyKnS8Yz8rKKvWe6enpGjdunE6fPq0R\nI0bo+uuvL1ePbdq0Kdd6AAAAoKLExsaW+vdqYmKiJBc94+Dp6SlJysnJKXa+4FGqXl5epdqvICzs\n27dPPXv21DPPPGOmUQAAAKCacMkzDt7e3nJzc1NaWlqx8wXjJV3K9N8OHTqkMWPGKCUlRb169dLc\nuXPl5lb+PFWQzEpj+PDhio2NLXdNAAAAoDRCQkK0dOnSy1rjkmccatasqcaNGyslJaXY+ZSUFPn5\n+cnHx+eS++zbt0/Dhg1TSkqKBg8erPnz5xc+VQkAAADAf7hkcJCkLl266OTJk0pOTi4yfuLECSUl\nJZX4xKUCycnJGj16tM6cOaNRo0bp73//u5EzDQAAAEB15LK/lMPCwuRwODRnzhw5HI7C8dmzZ8uy\nLA0ZMqTEtQ6HQxERETp79qxGjhypyZMnO6NlAAAAwGW55D0OktStWzf1799fUVFRGjp0qEJDQxUf\nH6/4+Hj17dtXPXr0KDx2wYIFRd4yvXnzZv3000/y9PRUrVq19MYbb1y0v7+/v4YNG+a0vwcAAACo\nylw2OEjSq6++qtatW2vNmjVasmSJAgICNHHixIte7rZw4UK5ubkVBoft27fLsixlZ2frnXfeKXbv\noKAgggMAAADw/7l0cHB3d9f48eM1fvz4Sx6XkJBQ5PNzzz2n5557riJbAwAAAKoVl73HAQAAAIDz\nEBwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nERwAAAAA2CI4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtggMAAAAAWwQHAAAAALYIDgAAAABs\nuXRwyMvLU2RkpO644w516tRJvXv31ptvvqnc3NxSrT937pxefPFF3XrrrQoODtZdd92l9evXV3DX\nAAAAgOtx6eAwffp0zZo1S35+fho5cqQaNWqk+fPna9KkSbZrMzMzNWrUKH366afq3LmzHnzwQaWn\npysiIkLLli1zQvcAAACA66hR2Q2UVXx8vFasWKF+/fpp7ty5heNTpkzRunXrFB0drR49epS4fvHi\nxdq3b5+mTZum++67T5I0YcIEDR06VK+99pr69esnPz+/Cv87AAAAAFfgsmccli1bJsuyFB4eXmQ8\nIiJCkrRy5cpLrv/444919dVXa9iwYYVjtWvX1iOPPKLMzEx98cUX5psGAAAAXJTLBocdO3bI19dX\nrVq1KjLeoEEDNW/eXHFxcSWuPXz4sI4fP64bb7xRlmUVmQsNDZWkS64HAAAArjQuGRyys7OVmpqq\na665ptj5wMBAnT9/XmfOnCl2/tChQ5JU7Hp/f395enoqKSnJWL8AAACAq3PJ4HDu3DlJkre3d7Hz\nBePp6enFzp89e/aS6+vUqaO0tLTytgkAAABUGy4ZHAoet+rh4VHsfMF4VlZWsfM5OTm260taCwAA\nAFyJXPKpSp6enpL+EwD+KDs7W5Lk5eVV7HytWrVs15e0trTatGlz2Wv27dun4cOHl6sugKoj9Zej\nOrOs+EsmAVeRdSrT5f636efjqXqjhMuVAVeSmp1l/Pu3b98+SVJsbGypf68mJiZKctHg4O3tLTc3\ntxIvJyoYL+lSpLp16xY57o/S09Pl7+9voNPLk5aWptjYWKfXBVBxsjJ+q+wWgHJzxf9tOsiVA6gm\nqtL3zyWDQ82aNdW4cWOlpKQUO5+SkiI/Pz/5+PgUO9+8efPC4/7o119/VVZWllq0aFGuHguSWWmE\nhYUpJSVFtWvXVrNmzcpVF85X8IUOCQmp5E6AKwvfPaBy8N1zbcnJycrIyFCTJk20du3ay1rrksFB\nkrp06aJ//OMfSk5OLvJj+8SJE0pKSlKvXr1KXBsQEKDGjRsrPj7+ormYmBhJUufOnc03XYLL/YeG\nqqXgNN/SpUsruRPgysJ3D6gcfPeuXC55c7T0+/9L73A4NGfOHDkcjsLx2bNny7IsDRky5JLrBw4c\nqGPHjumjjz4qHEtPT9fbb78tLy8vDRw4sMJ6BwAAAFyNy55x6Natm/r376+oqCgNHTpUoaGhio+P\nV3x8vPr27asePXoUHrtgwYKL3jL90EMPKSoqSi+//LJiY2PVtGlTbdq0SSkpKZo6dap8fX0r488C\nAAAAqiTL8d//d72LycvL06JFi7RmzRodP35cAQEBCgsL05gxY1SzZs3C44KCguTm5qa9e/cWWX/6\n9GnNmTNHX331lTIyMtSyZUs99NBD6tevn7P/FLiwglO2l3NfC4Dy47sHVA6+e1culw4OQFXAv0CB\nysF3D6gcfPeuXC57jwMAAAAA5yE4AAAAALBFcAAAAABgi+AAAAAAwBbBAQAAAIAtnqoEAAAAwBZn\nHAAAAADYIjgAAAAAsEVwAAAAAGCL4AAAAADAFsEBAAAAgC2CAwAAAABbBAcAAAAAtggOAAAAAGwR\nHAAAAADYIjgAAAAAsEVwAAAAAGCL4AAAAADAFsEBAAAAgC2CA64IR44cUVBQkMLDwwvHhg8frqCg\nIKWnp9uunzBhgoKCgnT06FHbY5OSkrRhw4Zy9QtcSR566CEFBQXp4YcfvuRxR44c0fTp09WnTx91\n7NhRXbt21T333KO3335bGRkZFx1f8B2/1Pd2zZo1CgoK0htvvFHuvwNwRZf6/q1bt05BQUF64okn\nbPd56aWXFBQUpM8///yiuczMTHXu3FlBQUH66KOPjPSNykFwwBXBx8dH4eHh6t+/f5Fxy7JKtd6y\nrFIdm5CQoAEDBuiHH34oU5/AlebkyZP6/vvv5eXlpa1bt+r48ePFHrd9+3bdeeedWr16tdq1a6cR\nI0ZowIABys/P17x58zRgwIBi15bme1vafw8A1Y3d969Pnz6qXbu2oqOjlZmZWeI++fn52rBhg7y9\nvdWnT5+L5jdu3KjMzEx5eXlp1apVxv8OOA/BAVcEb2/vYoODaefPn1dOTk6F1gCqk3/84x/Kz8/X\nQw89pLy8PH322WcXHeNwODRlyhR5enrqiy++0Lx58/TUU09p2rRpWr16tSZOnKgjR47opZdeKlMP\nDoejvH8G4JLsvn9eXl66/fbb9dtvv+nLL78scZ9vv/1WJ0+eVL9+/eTp6VlsHV9fX919991KTEzU\njz/+aPxvgXMQHACD+AECXJ61a9fKx8dHY8eOlbe3t1avXn3RMQcOHFBKSop69uyppk2bXjQ/fvx4\nNWzYUF999ZXy8vKc0TZQLZTm+xcWFiaHw6H169eXuM/n/6+9e4+puvwDOP7+clE5iGCoiRfcvMxD\njksatxRwAiakctOVYplN/6gMa5bmtK05pXTEWEOtdIm5MdmxAQqHUo6zHQZHQSMvw4LSUGKazLxQ\nIAe+vz8YJ+l8uWmQ9Pu8tvPHeb7Peb7PV32+ns/3+TzPOXoURVGIj4+3O3b9+nUsFguhoaE899xz\nqKqKwWD4R69DDBwJHMSgkZSUhJ+fH/fv3+9UnpiYiF6vx2KxdCrfvn07er2ea9euaa5x0NLW1sa+\nfftYsGAB/v7+LF68mOPHj/eqf5mZmaxcuRJFUThw4AA+Pj6Ul5cTHR1NQECAZg52ZmYmer2esrIy\nAPR6PRs2bODUqVMsWbIEf39/IiMjycjIsLtugNraWt555x1mz56Nr68vsbGxfP7551it1l71WYh/\n06VLl/jxxx+ZPXs2Q4YMISoqirq6OkpLSzvV6/j3/NNPP3XZVlpaGnv27OnX/grxX9Lb8RccHMy4\nceMwm82aawKbm5sxmUx4e3szc+ZMu+P5+fmoqkpYWBjPPPMMo0ePxmg00tTU1G/XJvqPBA5i0IiI\niKClpYWzZ8/ayu7cucOlS5dQFIWKiopO9UtKSpg2bRoTJkzo9Tk2btxIWloazs7OvPjii3h5eZGS\nksL333/f42eDg4NJSEhAVVUCAgJYu3Yt48ePJy4ujubmZs0ApKCggLFjxxIaGmorq6qqYvXq1eh0\nOpKTk3F3d+fTTz+1W7h28eJFEhMTOXbsGCEhIaxatQoPDw/S09N5/fXXZfZDPPby8vJQFMWWQhgb\nG6v5NHLatGmMHj2ac+fOsWLFCoxGI3fv3u1UJzAwkLCwMBwdHQes/0IMZr0dfwBxcXG0tLRQXFxs\nd8xkMtHY2EhCQoLmeY4cOYKzszNRUVEoikJMTAyNjY3dzmCIx5fTv90BIXorIiKC3bt3U1ZWRkhI\nCACnT5+mra0NV1dXysvLbXXr6uq4fPkyq1ev7nX7FouFo0ePEh4ezq5du3B2dgYgOzubrVu39riA\nMjAwEFVVyc3Nxd/fnzfeeANov+Hu2rWLwsJC4uLibPXPnz/PlStXWLNmTad2ampqSE5OZsuWLUD7\nLMi6desoLi4mLy/PNhX83nvvYbVaycnJwcfHx/b5HTt2kJWVxaFDh1i2bFmvr1+IgdTW1kZBQQGu\nrq6Eh4cD8Oyzz+Lp6YnJZOL333/Hw8MDACcnJ3bu3MnatWs5c+YMFRUVODg4oNfrCQ4OJioqilmz\nZnV5rqysLEaMGKF5rKqqShZHi/87fRl/0J6utGfPHoxGo1060pEjR3BwcNBMU7p48SI1NTVERkbi\n5uYGwMKFC/nyyy8xGAwkJib241WK/iAzDmLQ8PPzY+TIkZ1SkiwWCyNHjiQ6Oppz587ZUhrMZjOK\nohAREdHr9gsLC1EUhbfeessWNAAsX76cyZMnP3S/J06cyKxZsygtLeXWrVu28vz8fBRFYfHixZ3q\n63Q61q1bZ3vv4ODAhg0bUFXVts1dZWUl1dXVLFmypFPQAJCSkoKTk5NmrqoQj4uSkhJu3rxJdHQ0\nQ4YMAcDR0ZEFCxbQ0tJCXl5ep/qhoaEUFBSQnJzMqFGjUFWVqqoq9u/fT3JyMsnJyVy7ds3uPKqq\ncvDgQXbt2qX56m7BpxD/VX0df5MmTSIgIIDS0lJu375tK799+zYlJSWEhIQwduxYu/N0zGo8//zz\ntjI/Pz8mTZpEZWVlt+mH4vEkgYMYNBRFYc6cOVy8eNGWZ2mxWAgMDMTf35+mpibbTg1msxk3N7du\nn0L+3Q8//ICjoyN6vd7u2NNPP/1IfY+Li8NqtVJUVAT8tXWdXq9n2rRpnepOnz7d9mSmw8SJE3F3\nd+fSpUtA+1McgF9++YXMzMxOr3379uHq6mqrK8TjqCNwfvALBbQ/jVRVVXN3JS8vL7Zs2UJJSQm5\nubls3LiR2bNn4+TkxJkzZ1i1ahXNzc2dPqMoCidOnKCqqkrzlZqaKml94v/Ow4y/hIQEWltbOXbs\nmK2sqKiI1tZWzdmG1tZWjEYjw4YNY968eZrnka1ZBx9JVRKDSkREBAUFBZw+fRp/f39qampYtmwZ\nQUFBqKpKeXk5vr6+nDp1ioiICBwceh8b3759m6FDh2p+xt3d/ZH6HRMTw7Zt2ygsLGT58uW2pz1a\nqVRPPvmkZhujR4+mtrYWwJbfXVJSQklJiWZ9RVH4448/0Ol0j9R3If5pjY2NmEwmgC7TCWtqaqis\nrCQgIEDzuF6vR6/X88orr/Dzzz/z2muvUVtbi9FotMu1lsBAiL887PiLjY1l+/btFBYWsnTpUqB9\nN3AMKOAAAAYYSURBVCWdTqf52w1ms5mGhgYURelyHOfn57N+/XqcnOTr6GAhf1NiUJkzZw6KolBW\nVkZTUxOKohAUFMSUKVPw9PSkoqKCmTNncu/evT6lKUF7cHDt2jVaW1vtFlhq7YjUF8OHDycqKoqi\noiJu3LhBUVERTk5OLFq0yK5uVztN3Llzh5EjRwLt6UyKopCamtrlgjQhHldFRUU0NTXh5+fHU089\nZXf88uXLnDp1CoPBQEBAAFu3buXrr7/mq6++wsvLy67+5MmTSUlJYf369Vy5cmUArkCIwauv46+D\nm5sb8+bN4/jx4zQ0NHD//n3Onj1LUlISw4YNs2snNzcXRVGIjo7miSeesDteWlrK1atXKS4uZsGC\nBf/sRYp+I4GDGFQ8PDzw8/PDYrGgqiru7u62VJ+goCDMZjMnT57E0dGRsLCwPrU9Y8YMzp8/T2Vl\npV2KU29/rKa7RZZxcXEYjUZMJhNms5nQ0FA8PT3t6nWkIT2orq6OGzdu2J7qTJ8+HVVVOX/+vF3g\nYLVaSUtLY8KECaxYsaJX/RZiIHWkSWzatEkzDbC+vp7IyEiKiorYvHkzLi4u3Lp1i+LiYl566aVu\n2x4zZkx/dVuI/4S+jr8HZ60TEhL45ptvMJlMtgdqWg+v7t69y8mTJxkxYgTp6emaMwoGg4H3338f\ng8EggcMgImscxKATHh5OdXU1J06cIDAw0FYeFBTEvXv3yMnJwdfX1/Z0vrc6bn4ff/wxjY2NtvLC\nwkLNL/NaOm6OWr8ePWfOHDw9Pdm7dy8NDQ2aOaEAv/32G3v37rW9t1qtfPTRRyiKQlJSEtC+g9OE\nCRM4fPgwlZWVnT7/2WefkZWV1es+CzGQfv31VyoqKhg/fnyXa4e8vLwICQnhzz//pLCwkBdeeAEn\nJycyMjIwm8129W/evMmePXtwcXEhJiamvy9BiEHrYcbfgzr+HzOZTBQXF9s2//g7o9FIc3MzMTEx\nXaYhxcTE4OLigsViob6+/tEvTgwImXEQg05ERASffPIJ9fX1rFq1ylYeFBQEwL1795g7d26f2/Xz\n8+PVV1/liy++ID4+nrlz51JfX4/JZGLSpEm29QXd6VifYDQacXFxISEhgalTpwLtuyMtWrSI/fv3\n4+rqSlRUlGYbOp2OjIwMLBYLU6dOpaysjOrqauLj423pVw4ODuzYsYM1a9awYsUK5s2bh7e3Nxcu\nXMBiseDt7c369ev7/GcgRH/Ly8tDVVXNNL0HJSYmUlpaisFgYOnSpezYsYNNmzaxZs0afH19CQgI\nQKfTUVtby7fffovVaiU9PV0zJUII0e5hx18HR0dHFi5cSHZ2NlartcsfVe3YTam78wwfPpz58+eT\nn5/P4cOHefPNNx/uosSAkhkHMejMmDGDUaNG2dY3dJgyZYqtXCtwUBTFLpXo7+/fffddtm3bhk6n\nw2AwUF1dzfbt23u9XmLcuHG8/fbbODg4kJ2dbZfi1PE0dP78+QwdOlSzDW9vb3bv3k1DQwM5OTm0\ntbWxefNmPvzww071Zs2aZZviPXv2LAcPHqS+vp6VK1dy6NAhRo0a1as+CzGQOvZ87+mLS3R0NG5u\nbly4cIHq6mpiY2MpKCjg5Zdfpqmpifz8fPbv38+5c+dYuHAheXl5msF4b36jQeveIMR/0cOOvwfF\nx8fT0tKCoiiaM+dXr17lu+++Y/z48T3ubJiYmIiiKOTm5vb9YsS/QlFluwkhBkxOTg4ffPABWVlZ\nBAcH2x3X6/X4+PjITVQIIYQQjx2ZcRBigNy9e5cDBw7g7e2tGTQIIYQQQjzOZI2DEP2svLyc1NRU\nrl+/zq1bt9i5c+e/3SUhhBBCiD6TGQch+tmYMWO4efMmbW1trFu3rtvcUsm1FkIIIcTjStY4CCGE\nEEIIIXokMw5CCCGEEEKIHkngIIQQQgghhOiRBA5CCCGEEEKIHkngIIQQQgghhOiRBA5CCCGEEEKI\nHkngIIQQQgghhOiRBA5CCCGEEEKIHkngIIQQQgghhOiRBA5CCCGEEEKIHkngIIQQQgghhOiRBA5C\nCCGEEEKIHkngIIQQQgghhOiRBA5CCCGEEEKIHv0PhTxUtKsq104AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": { "image/png": { "height": 264, "width": 391 } }, "output_type": "display_data" } ], "source": [ "# Make a bad bar plot using Seaborn with its defaults\n", "sns.set_style('ticks')\n", "sns.barplot(x=strains, y=r/n)\n", "plt.ylabel('fraction of reversals')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll use Bokeh to make an ugly plot of the posterior probabilities." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# p_rev for plotting\n", "p_rev = np.linspace(0, 1, 300)\n", "\n", "# Make plots\n", "p = bokeh.plotting.figure(plot_width=650, plot_height=450)\n", "p.xaxis.axis_label = 'reversal probability'\n", "p.yaxis.axis_label = 'probability'\n", "\n", "colors = ['blue', 'green', 'red']\n", "line_style = ['solid', 'dashed', 'dotted']\n", "\n", "for i, r_val in enumerate(r):\n", " n_val = n[i]\n", " p.line(p_rev, posterior(p_rev, r_val, n_val), line_width=4, \n", " color=colors[i], line_dash=line_style[i], legend=strains[i])\n", "\n", "bokeh.io.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pretty plots\n", "\n", "Now, we'll use Bokeh to make pretty plots for our presentation. We will not label the axes, since we will label these on the slide with well-placed horizontal text. We also do not annotate, beyond tick labels.\n", "\n", "Note that I use the font Concourse T3 for the tick labels. This font is a font I bought and installed on my system. You will not have it (unless you bought it), so you will have to change that font selection to make the following code work.\n", "\n", "We'll start by specifying the colors for our bars and curves." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": true }, "outputs": [], "source": [ "colors_bar = ('gray', '#a6cee3', '#fdbf6f')\n", "colors_plot = ('gray', '#1f78b4', '#ff7f00')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll start with our horizontal bar chart. We use Bokeh's `quad()` function to get the bars." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set up canvas on which to paint plot\n", "p = bokeh.plotting.figure(plot_width=550, plot_height=290, x_range=[0, 1],\n", " y_range=[0.5, 3.25], \n", " border_fill_alpha=0, background_fill_alpha=0, )\n", "p.axis.major_tick_in = 0\n", "p.axis.minor_tick_line_color = None\n", "p.axis.major_label_text_font = 'Concourse T3'\n", "p.grid.grid_line_color = None\n", "p.axis.axis_line_color = '#a6aaa9'\n", "p.axis.major_tick_line_color = '#a6aaa9'\n", "p.axis.major_label_text_color = '#53585f'\n", "p.axis.axis_line_width = 2\n", "p.axis.major_tick_line_width = 2\n", "p.yaxis.visible = None\n", "p.outline_line_alpha = 0\n", "p.xaxis[0].ticker = bokeh.models.FixedTicker(ticks=[0, 1])\n", "\n", "# Grid lines\n", "for i in range(1, 4):\n", " p.line([0, 1], [i]*2, color='#A6AAA9', line_width=2)\n", "\n", "# Populate bars\n", "for i in range(3):\n", " p.quad(0, r[i]/n[i], 3-i+0.25, 3-i-0.25, color=colors_bar[i])\n", "\n", "bokeh.io.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll make a plot of the posterior probabilities." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Values of p for plot\n", "p_rev = np.linspace(0, 1, 300)\n", "\n", "# Set up plot\n", "p = bokeh.plotting.figure(plot_width=650, plot_height=450, x_range=[-0.01,1.01],\n", " border_fill_alpha=0, background_fill_alpha=0, )\n", "p.axis.major_tick_in = 0\n", "p.axis.minor_tick_line_color = None\n", "p.axis.major_label_text_font = 'Concourse T3'\n", "p.grid.grid_line_color = None\n", "p.axis.axis_line_color = '#a6aaa9'\n", "p.axis.major_tick_line_color = '#a6aaa9'\n", "p.axis.major_label_text_color = '#53585f'\n", "p.axis.axis_line_width = 2\n", "p.axis.major_tick_line_width = 2\n", "p.outline_line_alpha = 0\n", "\n", "# Populate glyphs\n", "for i, r_val in enumerate(r):\n", " n_val = n[i]\n", " p.line(p_rev, posterior(p_rev, r_val, n_val), line_width=4, \n", " color=colors_plot[i])\n", "\n", "bokeh.io.show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we will make a plot of the confidence intervals." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set up plot\n", "p = bokeh.plotting.figure(plot_width=550, plot_height=290, x_range=[-0.02,1.02],\n", " y_range=[0.7, 3.3], \n", " border_fill_alpha=0, background_fill_alpha=0, )\n", "p.axis.major_tick_in = 0\n", "p.axis.minor_tick_line_color = None\n", "p.axis.major_label_text_font = 'Concourse T3'\n", "p.grid.grid_line_color = None\n", "p.axis.axis_line_color = '#a6aaa9'\n", "p.axis.major_tick_line_color = '#a6aaa9'\n", "p.axis.major_label_text_color = '#53585f'\n", "p.axis.axis_line_width = 2\n", "p.axis.major_tick_line_width = 2\n", "p.yaxis.visible = None\n", "p.outline_line_alpha = 0\n", "\n", "# Grid lines\n", "for i in range(1, 4):\n", " p.line([0, 1], [i]*2, color='#A6AAA9', line_width=2)\n", "\n", "# Plot data\n", "for i in range(3):\n", " p.circle(r[i]/n[i], 3-i, color=colors_plot[i], size=10)\n", " p.line(hpds[i], [3-i]*2, color=colors_plot[i], line_width=6, line_cap='round')\n", "\n", "bokeh.io.show(p)" ] } ], "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }