Approximating likelihood ratios with calibrated classifiers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximating likelihood ratios with calibrated classifiers\n",
    "\n",
    "Gilles Louppe, January 2016.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import theano\n",
    "import theano.tensor as T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Toy problem\n",
    "\n",
    "Let us consider two 1D distributions $p_0$ and $p_1$ for which we want to approximate the ratio $r(x) = \\frac{p_0(x)}{p_1(x)}$ of their densities.\n",
    "\n",
    "- $p_1$ is defined as a mixture of two gaussians;\n",
    "- $p_0$ is defined as a mixture of the same two gaussians + a bump."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from carl.distributions import Normal\n",
    "from carl.distributions import Mixture\n",
    "\n",
    "components = [\n",
    "    Normal(mu=-2.0, sigma=0.75, random_state=0),   # c0\n",
    "    Normal(mu=0.0, sigma=2.0, random_state=1),     # c1\n",
    "    Normal(mu=1.0, sigma=0.5, random_state=2)      # c2 (bump)\n",
    "]\n",
    "\n",
    "bump_coefficient = 0.05\n",
    "g = theano.shared(bump_coefficient) \n",
    "p0 = Mixture(components=components, weights=[0.5 - g / 2., 0.5 - g / 2., g], random_state=10)\n",
    "p1 = Mixture(components=components[:2], weights=[0.5, 0.5], random_state=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: for $p_0$, weights are all tied together through the Theano shared variable `g`. This means that changes to the value stored in `g` also automatically change the weight values and the resulting mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_true = p0.rvs(5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xdc1dX/wPHXYTlARcUt4kbFjSJu3JilZlpqZWWlZVq2\nvza1PWxvq19lQ3NvyZGoKIgDFRXce6KIKIYgnN8fH0BAhAsCn3sv7+fjcR/ezzif+75Xfd9zz+cM\npbVGCCGEfXEwOwAhhBCFT5K7EELYIUnuQghhhyS5CyGEHZLkLoQQdkiSuxBC2CGLkrtSKlApFa2U\n2qeUejmH4wOVUjuUUhFKqXClVOdMx45kPlaYwQshhMiZyqufu1LKAdgH9AJOAZuB4Vrr6EznlNVa\nX0173gKYpbVumrZ9CPDVWl8smrcghBAiO0tq7n7Afq31Ua11MjATGJT5hPTEnsYNSM20rSx8HSGE\nEIXEkqRbCzieaftE2r4slFKDlVJRwGJgdKZDGliplNqslHr8doIVQghhmUKrUWutF6Q1xQwG3sl0\nqLPWui1wB/CUUqpLYb2mEEKInDlZcM5JoE6m7dpp+3KktQ5RStVXSlXSWsdqrU+n7Y9RSs3HaOYJ\nyV5OKSWT3AghRD5prVVO+y2puW8GGiqlvJRSLsBwYFHmE5RSDTI9bwu4aK1jlVJllVJuaftdgb7A\nrlyCNPXx5ptvmh6DtTzks5DPQj4L6/8scpNnzV1rnaKUGg+swPgy+FlrHaWUGmsc1tOAe5RSo4Ak\n4D/g3rTi1YD5abVyJ+BPrfWKvF5TCCHE7bGkWQatdRDgnW3fD5mefwR8lEO5w0Dr24xRCCFEPkkX\nxUwCAgLMDsFqyGdxg3wWN9jKZ1G3bl2UUkX6mDJlSpG/RvZH3bp1Lf4M8hzEVFyUUtpaYhFC2Dal\nVJ5t0rYo+/tK2y7wDVUhhBA2RpK7EELYIUnuQghhhyS5CyFEIbty5Qrh4QWfBPfAgQOcPHnLsaIW\nkeQuhBCFbO7cufj5+eV4bNasWRnPw8PDWbZsGXPmzMlyTsOGDVmzZs1txSDJXQghClFqaioJCQm3\nPH79+vWMPxcvXswdd9yRYyIvU6YMMTExBY7DokFMQgghcrZp0yYWLlxI37592blzJ3369MHT0xOA\nEydOsHfvXlavXo3Wmqeffjqj3LJly1BK8d133xEcHMy6detwdXXF19cXMGrvO3fupFevXgWKS2ru\nImd798LXX8NHH8GqVZCamncZIaycUvl7WMLT05Pk5GQCAgIICwsjLi4ONzc3AJydnenVqxfXrl3j\ntddeo0aNGhnlQkJCGDduHKVKlWLUqFF069aN9evXZxx3c3Pj/PnzBX6vUnMXWV2+DBMmQFAQDBwI\nbm4kPDWB1KRkdj4/hcsNmwDg7u6Cv38Lk4MVIn+KYlxTpUqVKFeuHFprnJ2dKV++fEZSLlu2LHv2\n7KFJkybExcXh6uqaUe7KlStUrVqVsLAw6tWrB8B///2XcTw+Ph53d/cCxyU1d3FDbCx06wZOTnDg\nAEybBp9+yvrPfyfmsbfwe30iDfbHUKWKL3FxSWZHK4RV2Lp1K1evXuXHH3/k7bffxtvbm+PHjfWN\n3n77bXbs2EFiYiJ79+7NUq5///7MnDmT+++/n/LlywPg4uKScXzPnj23vClrCam5C0NSEtxzDwQE\nwKefZv1NqhQXA0eSVLMeDZ4byIEvl0MVC3+zCmHnwsPDeeaZZ7I0uZQtWxaAjz66aT7FjOkD7rrr\nrox9SUlJrFy5kt69e2fsS01NpWLFigWOS2ruwvDqq+DmBlOn3rKxMaFlR45N+p76Lw/F+VJcMQco\nhPWJj49n2bJlxMfHZ9l/5513snnz5hzLODs737SvT58+9OnTh1atWgGwb98++vbte1uxycRhAsLD\njfb1yEioUuWmw0FBW6lSxTdju/YnE7l+dh+1Vi8rziiFsJhMHCY1d5GaCmPGwGef5ZjYc3Jq3LtU\n2rUdVq4s4uCEEAUlyb2kmzEDypSB4cMtLpJaxpU9416Cp5+GtAEZQgjrIsm9JEtKgtdfhw8+sLxT\nb5oYvy5QtSpMn15EwQkhbock95Js9myoWxe6d89/WaXg/fdh8mTjS0IIYVUkuZdUWhtdHp9//pan\npKTAtWu5XKNTJ2jSBP76q/DjE0LcFknuJdW6dZCQAP3733QoJAT69QNXVyhXDsaMacbff9+ief2F\nF4zuk3bYM0EIWybJvaT67DN49llwuPFPQGuYMgXuuw9GjIALF+C//2DixKMEBxudai5ezHadPn2M\nEa3Llxdr+EKI3ElyL4nOnIG1a+GBB7Lsfv11WLgQtm2Dhx82au6OjtCsWQLffAO+vjB2LGQZr6GU\n0bTz2WfF+haEsGY2s1iHUipQKRWtlNqnlHo5h+MDlVI7lFIRSqlwpVRnS8sKE/z+OwwZYmTvNH/9\nZTz++QeqVbu5iIMDjBsHHTrAiy8a7fEZhg2D7dvh4MGij10IG2ATi3UopRyAr4F+gA8wQinVJNtp\nq7TWrbTWbYBHgZ/yUVYUJ63hl1/gkUcydh09anRZX7Ag93FMSsHEiUYrzMyZN+bRoHRpGDUKfvqp\nCAMXwjZYy2IdltTc/YD9WuujWutkYCYwKPMJWuurmTbdgFRLy4piFh4OycnQ2fhxpbXR1PLss9Cy\nZd7FHR3hzTdhyZIqbN2a6cDjjxtfGtItUpQwmzZt4pVXXiE4OJgvv/ySvXv3ZlmsY/Xq1bzyyitM\nmjSJ06dPZ5TLvFiHg4MD8+fPz3Ld9MU6CsqS5F4LOJ5p+0TaviyUUoOVUlHAYmB0fsqKYvTnn/Dg\ngxmDlhYvhv37E2nVaitBQcYjLCwy10tUrQqPPHKS8eMzreHRpAl4exsXFMJaFcFqHYWxWIefnx+r\nV6/Ocl2rWaxDa70AWKCU6gK8A/TJ7zUmT56c8TwgIICAgIDCCk+AkYnnzoW0f0TXr8OkSfDQQ8ep\nUePGxGAxMVtvdYUMvXtfICSkbsZ3BQAPPQR//GFMHSyENSqCLru3u1jHpk2b+PzzzzPmgE93q8U6\nMufJ3FiS3E8CdTJt107blyOtdYhSqr5SqlJ+y1oatCigsDCoWNGoZWNU4itXhvbt4/MoeDMHB6OD\nzP33G10nXVwwkvqzzxqLflSqVMjBC2Gdsi/WUbNmTdauXQsYi3W0adMmY7GOWrVuNFykL9YxcuRI\nypQpc9Mslnv27OGOO+646fUy58kpU6bcMi5LkvtmoKFSygs4DQwHRmQ+QSnVQGt9MO15W8BFax2r\nlMqzrChGs2fD0KGAUYH56CP44ouCz/3VuTM0bmxML/PYY0CFCtC3L8yZY3SKF6IEKIzFOk6cOMGu\nXbs4efJkxhdAkS/WobVOAcYDK4DdwEytdZRSaqxSKv1/8D1KqV1KqW3AV8C9uZUtcLSi4FJTjaQ7\nbBhgLJHq7AwFXFg9wxtvwHvvGfdoAaPv/J9/3t5FhbARhbVYR+3atZkxY0ZGYpfFOoTlNm0yRibt\n2QNK0auX0RvygQduXowjJmYrgYE3trMfz35Oz55Gc/tDD2H0lqlZ0xgJVacOQphBFuuQEaolx4IF\nMGQIYZt28f33u9mxI4kKFbYRFLSVqKjDt3Xp5583mne0xmh8HzIEMg3UEEIUP0nuJcWyZXDnncTF\nJbFunQ+DBrlQs2ZbqlTxJSHh9hbc6N8fLl+GjRvTdgwdavTKEUKYRpJ7SXD8OJw8CX5+JCUpli83\nlkwtLA4OMGECfPll2o6AANi3z3hNIYQpJLmXBMuXG3P4OjoSFuZOo0ZQu3bhvsTDD8OKFXD6NEbT\nzJ13QrYRd0KI4iPJvSRYtgzS+suuWFG5UGvt6fbsiaRDh/O89toJgoK2sq1uSy79IkvwCWEWSe72\n7to1WLMG+vUjJgaiotwoioG/cXFJ3HuvB6tX18bDwxf6jqNs9B64jYmPhLAlv/32G9NzWFP40qVL\nN80bUxwKbfoBYaXWr4dmzcDDg3k/QLt2lyhTpmhGj6ZPPBYZCS1bluG8b0dqLFyYNsJJCHOEhUUS\nF1fwCe3c3V3w929R4PJxcXHMmzePu+++u8DXKAhJ7vYuU5PM339Dt24XgaJJ7koZN2oXLTIS/dnO\nPakxb54kd2GquLikm8Zp5Educy0lJyczbNgwkpKSKFu2LIGBgfTu3Zvr169TpUoV/v77b6ZNm8bK\nlSvp2bMns2fPZty4cZw7d45SpUoxZ86cjEnGCps0y9i5q/MXsrFCXf76awfh4ddxc9tRpK83YIAx\nL1liIpxv18lYkDWXua2FsGULFiygQ4cOLFu2jMqVK+Ps7MzSpUsJDg6mSZMmrFmzhjFjxtC3b1/+\n/fdfKleuzG+//caaNWsYNmwYf//9d5HFJjV3e3bqFE7nz1PKfyQ75jrSrRskJxftfOseHuDjY+T0\nVq3coH17I9sXxV1cIUx26NAh2rRpA4Cvry9Xrlxh9OjRnDx5knPnztG4cWMaNmyYcX5qaiovvPAC\nkZGRXL58uUibaqTmbs/WrCG2pS84OrJqFfTuXTwv27evsVwfYHSJXLKkeF5YiGJWr149tm/fDkBE\nRATXrl3D29ub4OBghgwZkjENcPrqS9u3b+fq1ausXbuWcePGFekUCZLc7dnq1cS2bs+lS7B3r7H+\naXHo0cNY8CkhweFGcrfDeT6EGDx4MBs3bqR///5cunSJxMREFi1axMCBAzl69CgA1atXJzY2lnvv\nvZeqVauyf/9+7rjjjltOLFZYpFnGXmkNq1dz4dVP2LgRfH2NpU4tERV1MNv24XzdkCpXzni90FB3\n7rmnnrEjIgLats3POxCiULi7u1i0AE1u5W/FxcWFRYsWZdn3yiuv3HReUFBQxvP169cXOJb8kORu\nrw4dguRkEjzrsm4OdOtmedGEhNQsyXzLlv03nZPXF0DfvrBgQVqvnPTauyR3YYLb6cZoy6RZxl79\n+y/06sX1FEVYGHTpUriXT/8CSH9kn3ysWzfYs8eN8+eRdnchTCDJ3V6tXg09e7J7txuenkYvluJU\ntqwxYGr+fIxvlv374cyZ4g1CiBJMkrs90tqouffsSXi4O127mhNGp05xRnJ3djbaaZYtMycQIUog\nSe72aPdu4yamlxebNlUwLbm3a3eJkBCIj0eaZkSx8vLyQilldw8vLy+LPwO5oWqP1q+Hbt3Yvx8S\nEx3w9jYnDFfXVLp0MSrsw/v1MyZ9T042avJCFKEjR46YHYLppOZuj9avh65dWbkS2raNR+W4wmLx\nuPtuY4U/qlaFhg0hNNS8YIQoQSS52xutYd26jOTeuvVlU8MZOBCCgoyZhwkMNDaEEEVOkru9OXIE\nUlK4XrchwcHQpk28qeFUqwbNmxudd+jXL9O8BEKIoiTJ3d6kNcls3qKoUwcqVry9xa8LQ0bTjL8/\nHDwIZ8+aHZIQds+i5K6UClRKRSul9imlXs7h+Eil1I60R4hSqmWmY0fS9kcopcILM3iRg0zt7X36\nmB2MYfBgWLgQUhycoVcvY7FVIUSRyjO5K6UcgK+BfoAPMEIp1STbaYeAblrrVsA7wLRMx1KBAK11\nG621X+GELW4pU3Ivrlkg89KggdE8ExqKtLsLUUwsqbn7Afu11ke11snATGBQ5hO01mFa60tpm2FA\nrUyHlYWvI27XuXNw5gyX67YgIiJ/88kUtcGDYfFijHb3lSshNdXskISwa5Yk3VrA8UzbJ8iavLN7\nDFieaVsDK5VSm5VSj+c/RGGxkBDo1Im1IY74+RlTAJgpKuogQUFbCQraSuXKUcyenQh16hhzIWzb\nZm5wQti5Qh3EpJTqATwCZJ6mqrPW+rRSqgpGko/SWofkVH7y5MkZzwMCAggICCjM8OyflbW3Z55d\nsnJlmDw5mUOHoH5600y7diZHKIRtCQ4OJjg42KJzLUnuJ4E6mbZrp+3LIu0m6jQgUGt9MX2/1vp0\n2p8xSqn5GM08eSZ3UQDr18Pnn7P6Cfi//zM7mKwcHKB9+0ssXerBhMBAePtteO01s8MSwqZkr/RO\nmTLlluda0iyzGWiolPJSSrkAw4Ess9MrpeoAc4EHtdYHM+0vq5RyS3vuCvQFdln8ToTlLl+G6Ggu\n1G/PsWPWOXW6n98lY3qZrl1h+3aIizM7JCHsVp7JXWudAowHVgC7gZla6yil1Fil1Ji0014HKgHf\nZuvyWA0IUUpFYNxoXay1ln5wRWHTJmjdmo1bS9GhAzhZ4axBbdrEs3EjXEkpY0wDvHq12SEJYbcs\nSgFa6yDAO9u+HzI9fxy46Wap1vow0Po2YxSWCAuDjh0JCSn8hTkKi6trKv7+sGoVDA4MNEar3nOP\n2WEJYZeki6K9CA1lb+WaLFlyBReXfRm9VKKiDpsdWRYDBqTN/Nuvn3FTVRbOFqJISHK3B1pDWBhH\nqjXj0CE3unRpfMvl78x2552wdCmkNvI27rJGRZkdkhB2SZK7Pdi/H9zciLxQl/r1ze/fnpuGDaFC\nBdgWoWS0qhBFSJK7PQgNBX9/du92o7UV3+FIH9Tk43OGL744xbaqjYj7e47ZYQlhlyS524O0m6nW\nntzTBzX17Vudbdtqono9jlvkDkhIMDs0IeyOJHd7EBpKSnt/oqJcadXK7GDy1ro1HD8O5xLLc6lR\nU7BwxJ0QwnKS3G3dlSuwfz97SrWhQoXrVK5sdkB5c3IypnYPCYHzvh2l3V2IIiDJ3dZt3gytWrE+\nvBTNml0xOxqLdelizJZwvl0nSe5CFAFJ7rYuLAz8/QkJAR8f20nunTvDli1woVZjo839wAGzQxLC\nrkhyt3WhoRkjU20pubu7G4t4RO4qL10ihSgCktxtWdrgpZOe/ly7BrVqXTM7onzp2hU2baogyV2I\nIiDJ3ZYdOgQuLqw95EmXLqCU2QHlT9eusGFDWVbhwfV/g1mxaCNhYZFmhyWEXZDkbssytbdb62Rh\nuWnQAFJSFLH0JLFRK+oeTyAuLsnssISwC5LcbVmm9nZbTO5Kgbf3Kdavh0udAqkQKk0zQhQWSe62\nLDSU+Gb+HD6MVY9MzY239yk2bID4joGUl+QuRKGR5G6rrl6FqCg2/NcWPz9wdjY7oIKpV+8s+/bB\n6Zq+OMWeo/S5M2aHJIRdkORuq7ZsgebNWbe5jE02yaRzdk7F1xdCNzkQ37EfHls2mh2SEHZBkrut\nsoGVlyzVtasxWjW+YyAeW0PNDkcIuyDJ3VaFhpLk25GICGOeFlvWubNxbzi2XV8q79gMSdJjRojb\nJcndFqUNXtpZ1p8mTaBcObMDuj1Vq0LNmrDtRFWu1vQ0Mr0Q4rZIcrdFR48CsPqAl803yaRLb5qJ\nkYnEhCgUktxtUfrgpQ3KbpJ7ly4yBbAQhUmSuy0KDSXVvyMbNhjt1fagaVOIj4eo8m2NXyanTpkd\nkhA2zaLkrpQKVEpFK6X2KaVezuH4SKXUjrRHiFKqpaVlRQGEhXGkmj+VKkGNGmYHUzgcHIwvqvBt\nlaFfP1i61OyQhLBpeSZ3pZQD8DXQD/ABRiilmmQ77RDQTWvdCngHmJaPsiI/EhNh1y7+jW9nN00y\n6YxZIt1h4EBYtMjscISwaZbU3P2A/Vrro1rrZGAmMCjzCVrrMK31pbTNMKCWpWVFPm3bBk2aEBxe\n1u6Su58fREe7crlzIKxda4zCFUIUiCXJvRZwPNP2CW4k75w8BiwvYFmRl0yThXXtanYwhcvVFZo2\nvcLKLRWhXTtYtcrskISwWU6FeTGlVA/gEaBAdcrJkydnPA8ICCAgIKBQ4rIrYWFc6DqYhARo3Njs\nYAqfn98lliypwJD0ppmBA80OSQirERwcTHBwsEXnWpLcTwJ1Mm3XTtuXRdpN1GlAoNb6Yn7Kpsuc\n3MUthIYS2uUDm1ycwxJ+fpf43/8gddJdOHzwAaSmGndbhRA3VXqnTJlyy3Mt+V+zGWiolPJSSrkA\nw4Esd7uUUnWAucCDWuuD+Skr8uHECUhKImhffbtrb09Xo0YSVarApvMNwMMDNm82OyQhbFKeyV1r\nnQKMB1YAu4GZWusopdRYpdSYtNNeByoB3yqlIpRS4bmVLYL3UTKEhtrd4KWc3H03zJ8P3HWX9JoR\nooAsanPXWgcB3tn2/ZDp+ePA45aWFQUUFsZ/rfw58C+0aWN2MEVnyBAYNgw+nD4QNXYMvPuu2SEJ\nYXOkMdOWhIay07Uj7duDi4vZwRSd1q0hJQUiy/jBuXNw+LDZIQlhcyS524pr12DHDoIutLe7LpDZ\nKWXU3uctdIQBA6RpRogCkORuK3bsgEaN+Dfcza7b26OiDhIUtJWaNaP57berbPVqSfyv080OSwib\nU6j93EXhCQuLJC7uxqIVTf6Zi2d7f7bOsP3FOXKTkJBKlSq+dOsG778PJ5s9Scupb8Dp0/YzkY4Q\nxUBq7lYqLi6JKlV8Mx6lIrZzuJo/jRtD+fJmR1f0HByge3dYHVKKc35d0rrPCCEsJcndRrhHRRJ8\nraNdN8lk17MnrFkDZ7v0gtmzzQ5HCJsiyd0GOJ0/jePVBJbsbVSikruvrzG1+16vrhARAWfPmh2S\nEDZDkrsNcN21ibgmzVm/wcFuFuewhLOzsULT+i3VoX9/WLDA7JCEsBmS3G2A285QjtZoQ4UKUKuE\nzanZpw+sW1fRGNU0Z47Z4QhhMyS52wDXXWFsduhg9/3bc+LvD8eOleG4TyCEh0NMjNkhCWETJLlb\nu+vJlI3ayj9xXUpUe3s6Z2fo3PkiMxeVNQY0yY1VISwiyd3KlTkQSVINL8L31iyRyR2gW7eLzJwJ\nPPAA/PGH2eEIYRMkuVs5t8gwzjfsyNWrjjQpoavPtmx5mVOnYJ9XHzh40HgIIXIlyd3Kue4MJdLN\nn2bNrtjl4hyWcHSEe++FmXOd4b77pPYuhAVk+gEr57orjH+bvYyHx1GCgrLWWI8fP4an542FrqKi\nDlOlim9xh1gshg+H0aPh9V8fQD1wP7zxhn0uRSVEIZGauxVzuhiDc+w5lh5uRo0a57JMR1Clii9n\nz17Lsp2QcN3skIuMv78xMWaEU3sjqYeHmx2SEFZNkrsVc921ifgmfhw97kDNmhfzLmDHlIJRo+CX\nX5XcWBXCApLcrZjrzlAOePjTtCk4OaWaHY7pHnoIZsyAa8MegJkzITHR7JCEsFqS3K2Y244NrE/t\nbNdL6uVHvXrQsiUs2lUfWrWSmSKFyIXcULVS6vp1ykZtYUGjjtw3CC6W7FaZDKNHw6efXqZRj954\nfvgpmys2BsDd3QV//xYmRyeE9ZCau5UqfzCaxFoN2LK/Ai0kZ2UYMgR27SrDGf/nKH/8GLWvlqNK\nFd8sC5sIISS5W62Ku7dztHZn6tUDV1ezo7EeZctCly4XWbLChQsDHsJjwU9mhySEVZLkbqXc9+xg\nk1MXWrc2OxLr07fvBRYtgnMDH6Pykl9RSdfMDilPYWGRBAVtzXiEhUWaHZKwcxa1uSulAoHPMb4M\nftZaf5jtuDfwC9AWeEVr/WmmY0eAS0AqkKy19iuc0O2Y1lTcs5PF9TrTOdDsYMwXFZV18JbWhylT\npgnrzjSmbsMWVFw1m3Ptm5oUnWXSl01MFxOz1cRoREmQZ81dKeUAfA30A3yAEUqp7LOcXAAmAB/n\ncIlUIEBr3UYSu4UOHiTVwZEV0XWk5s6NRbPTH1evXufee2HWLDg3YiJV//oMtDY7TCGsiiXNMn7A\nfq31Ua11MjATGJT5BK31ea31ViCnIZLKwtcR6TZs4LiXL1WrKdzdzQ7GOgUGQmQk7Kk3AMerl6m4\nK8LskISwKpYk3VrA8UzbJ9L2WUoDK5VSm5VSj+cnuBIrJIRtpf2lf3suSpeGu+6COfMcODdiInXn\n/2V2SEJYleLo595Za31aKVUFI8lHaa1Dcjpx8uTJGc8DAgIICAgohvCs0IYNrHD5VpJ7Hu65Bx5+\nGJ6c/RCtv/4f6/5vAVdregLS713Yp+DgYIKDgy0615LkfhKok2m7dto+i2itT6f9GaOUmo/RzJNn\nci+xLlxAnzjBQtWB6ZLcc1W7dtqI1dWupLbtQZPlKzj+v28BuWEp7FP2Su+UKVNuea4lzTKbgYZK\nKS+llAswHFiUy/kZ87AqpcoqpdzSnrsCfYFdFrxmybVxI1ebd8CptAPVq5sdjPV7+GH4/XfY2GEA\nlVbMxPnsCbNDEsIq5JnctdYpwHhgBbAbmKm1jlJKjVVKjQFQSlVTSh0HngVeVUodS0vq1YAQpVQE\nEAYs1lqvKKo3Yxc2bCDaowvNm18xOxKb0LIlVKsGm4825/ygR6n+6wdmhySEVbCozV1rHQR4Z9v3\nQ6bnZwHPHIpeAaQzX36EhLCy9GR8fK4Alc2OxiY89BC8915TzvzWgOb3NuXMI5My/X4UomSSLorW\n5L//YPt2ft/XQWru+dCpk/Hn2uhqXLjrEam9C4Ekd+uyaROJjZoTm1wOT0+Zq9xSSkH37nuYNg1O\nP/gSlf75izKnpe1dlGyS3K1JcDD7agQQECDLg+aXj89xrl+HlTurcXbERBr/8rXZIQlhKpnP3ZoE\nBxOkXqbn3WYHYnscHGDcOPjiC+jx6/P43PMNbNx4o83GBoWFRd40lbH03xeWkpq7tUhMRG/Zws/R\nnenRw+xgbFPnzlChAixeXZb9D4+DZ5+FVNtdnjB9srHMD5m3XlhKau7WYtMmEhv48N/F8jRoAAcO\nmB2Q7VEKnnkGXnoJSj/TDM8r8zj19Kscv3MoAMePH8PT88Z4PKkFC3smNXdrERzM3mrd6dFD2ttv\nR4sW4O8PS5f7cOrNP2n850/UpAZVqvhy9uw1qQWLEkOSu7VYu5bliQHSJFMIxo+HiIj6RDn4EDP0\nSTw/Gi9TAosSR5K7NUhMRIeH8+OeLpLcC0HlyhAQsJv33oNTD71C6WP7qbR0utlhCVGsJLlbg/Bw\n/qvbFFWhPF5eZgdjHzp02E9qKsxcUJpD786g9ucvUPHCGbPDEqLYyA1Va7B2LdHVAuhZ3+xA7IeD\ng+bNN+GRR6BTp+aUe+wN7p7xOWeGPYEuVdrs8HK0atVutm93JC7OmeRkRXx8AgMGQI0aZkcmbJEk\nd2sQHMw+R2AmAAAc10lEQVTy/56TJplCVqcOjB0Lr74KP/80Hvflf1Lngyc5+sb/gVI3rc0Kefeg\nKey+5ydOwPTpsGABREY2oUEDR6pWBScnOHDAlVmzwNUVeveGe+8t0EuIEkqSu9nS2tunOXUhTJJ7\noRs2DCIiYOonig5DnmDCzM+oOuMLzo2cmLE2a2Z5zQOffaFrS8rk5Nix0owYAf/8A8OHwwcfwNWr\n26lR48a1ly9fQ2DgcPbtg4UL4b77oG/fWnTpAm5u+X5JUcJIm7vZQkK4Urc55WtXkJ/fRUApeO01\nI8GH7mzGwU8WUu33j3FfNduUeK5cgc8+gxdfbEzbtnD4MHz7LfTsadTWs1MKvL2NvvuzZkFsrDMt\nWsDmzcUfu7AtktzNtnIlEZX70K+f2YHYL1dXmDoVVq9uweqDdTnwxTLqfDSe+vu3F1sMqamweDEM\nHWok+OefX0aLFlsJDd1KUJDxiIo6nOs1PDzgxReP8MknMGAA/PxzMQUvbJI0y5ht5UpmJn/B4L5m\nB2Lf6taFkSPXM2VKX8pNbYXDx/MZNCGQ0226cKnrnUX62ps3w3PPeePkBJ98Aj4+sHz5fzc172zZ\nst+i6w0ZAs2bQ79+EBMDL7+cdeBb9vsCMhK3ZJKau5liYtAHDzLjsD9du5odjP3z9IzlnXfghRdg\n1dVO/P3AS3i98xiVlhRNH/izZ+HRR2HgQBgw4Dy//GIk9sLQuDGEhBhLDL7+etZj2eekkZG4JZMk\ndzOtXs2Zxt3w6+xMmTJmB1My+PsbTTRvvgnLLvRi33f/Uuu716jxw+RCm2QsKcloV2/eHNzdIToa\n+vS5gEMh/2+rVQuCg2HuXOMXgRCZSbOMmVauZH2pPvSVJpli1bq1cRNz3LjmpKaW438/bKLplPtw\n3R1O3MQXCnxdrY0ujS+9BI0awbp10LRpIQaegypVYMUK6NLFeD5qVNG+nrAdktzNojWsXMlP+gU+\nlZupxa5hQ3jyyX/YvHko9z5Tg5efW82Q8P/RadxIcP0/uOMOi6+lNSxZAu+9B5cvw/jxh/H2juXo\nUTh61DgnKurwTW3shcXTE4KCoHt340tFCJBmGfPs28f1pFR2pzQptHZYkT+lS1/nvfdg4kR49yNn\n7tr3CXPv+gg9frxRBT59OtfycXEwcyaMG9eM114zphvevh28vWNvmoc9IeF6kb6Xpk3hl1+M3jgx\nMc5F+lrCNkhyN8uqVeyv14e+/ZRM8Wuyrl1hzhyje+GU9ffQMGEnQdurk9ioOUdHT2HvtgT27DF6\nvQQHV+Tbb2H0aBg0CHbvhjFjjhMRYQwycnQ0730MGABPPw1vv92AJLmHWuJJs4xZgoJYcm2k9G+3\nEs7OcOed0KFDNPXr+7JixUesWP0kfee9SNvfGvCb21hmVBpDSnkX2reHJ56Ali2hdGmIiblsNV/Q\nL70Eixcn8fnnrrz0ktnRCDNZlNyVUoHA5xg1/Z+11h9mO+4N/AK0BV7RWn9qadkS6b//0GvX8rX+\nje2S3K1O48bGg/H1CAqaRMzlMjzy+8c8u7YF26r74dpnKokNm2ecn32OmqJsX8+LUjBx4lGeeaYi\n7doZI19FyZRncldKOQBfA72AU8BmpdRCrXV0ptMuABOAwQUoW/IEBxPr2Yr6VStRsaLZwYi8JNZv\nxtE3f+HUmeNcmzqR1k8HklS1NucHP8bF3sNumqPG0sFIBZH9iySnAUpubim8956xhKy3N7i4FFk4\nwopZ0ubuB+zXWh/VWicDM4FBmU/QWp/XWm8Fst81yrNsibR0KWvdBnDXXWYHIvIjubona3sNI3Lx\nUU4/9gYVNiyj5YA6DPtzKpWW/4nDlfgijyH9iySvAUrNm8NDD8Err0ByspW0GYliZUlyrwUcz7R9\nIm2fJW6nrH3SmsR58/lsb1/Kl99FUNBWwsIizY5K5IejI/Fd7uDQx/PYufQY0T4dqLhiJi0HeFL/\nhbvx2bGhWBJ9Xu6/3xhE9fvvNc0ORZjAqm6oTp48OeN5QEAAAQEBpsVSZKKjSU5M5Vj5NrRpY9So\nCjJlrLAOqW4ViGzdldr9h+N4OQ734AX4zPyKBgNqc7ltdy72GsqlbgNNiU0pYyTuiBGVWL0aevUy\nJQxRiIKDgwkODrboXEuS+0mgTqbt2mn7LJGvspmTu91aupRwjz509bvxUzmnRSPMvCknCialnDsX\n7nqY5U6lGdC1P+7rl+C+eg51Pp6AW836OF1LIC5gECnuHsUWU6VK8OyzR3joocZERBijWIXtyl7p\nnTJlyi3PtSS5bwYaKqW8gNPAcGBELudnbuDLb1n7t3Qpf195KstEYTktGlGUN+XErRVWz5dUtwrE\n9r+f2P7343D1Coe+fpUeYf/g+dlzJPi052w1L5z8enC9crXCCv2WfH0vM2KEMYnZwoVYTbdNUbTy\nTO5a6xSl1HhgBTe6M0YppcYah/U0pVQ1YAtQDkhVSj0DNNNaX8mpbJG9G2sXG0vq1m0sSO3Lo63N\nDkbkpCh6vqSWdSOqRUfq9h+OSrxKhY1B1PnjU5oMbcJ/jVpxsddQLva977ZfJzfvvgudOsE338D4\n8UX6UsJKWNTmrrUOAryz7fsh0/OzgKelZUusxYs50qAXPu7JOa66I+yfLl2WuJ5DWH4tiTt6Dqb8\nppVUXDmLmt+9RinPxji5uXGpY2DOyzIVUPqvkSefLMVzz3nj6LifNm1SZI53OyfTDxSnefP4O2kI\nXbrEmR2JsAK6VGkudbuLI2//TuSSYxxs3Irqv7xHywGe1PriJSrHWHprK3fpv0Zat27Oc885M3Vq\nM86cSS6UawvrJfXH4nLlCqn/ruE79RtftTtkdjTCyqS6lSeiXS+qv/4jpY5E47HoFx78v7dJDVvO\nuXsnGKtF5TBxTX5vxg8YAKGh8OOPngwenOMpwk5Ici8uy5dz0qsTHX3cKVVKmx2NsGLX6jbh5NMf\nsqJec4Y7OVD91/fx/OQZYoaOo4xb1iHN+b0ZrxRMmgTDh5dn3jxjyT5hn6RZprjMm8dcPYShQ80O\nRNiKFCdnYvvfz95fwzj0wWzKHNrNuM+exfPDp3A5UfBff25u8NJLh3nySTh+PO/zhW2S5F4cEhNJ\nXbacr48Nys8aEEJkuOrTniNTfuP7Z6aS4laBpg/5Ue/VkVQ9fbRA12vaNIFnnoEHH4SUlEIOVlgF\nSe7FYeVKzlVrQet+1XB1NTsYYcsS3Nw59dR7RC48xFXvNoz4/QMaPt0ft23r8n2tl182mmnef78I\nAhWmkzb34vDXX8x2HCFNMqLQpLqV5+yoF/m9YnVGpFzD663RJFWvw+mxtx6xmJ2jI/zxB/j6GlMT\ndOxYhAGLYifJvahduULq0mV8rr4i0pwpRoQdS3Fy5sJdD3LhzoepvOx36r45ipFlypFUy4uElnln\n61q14IcfYORIiIgwJhoT9kGaZYraggUcqd2FbkM8KFvW7GCE3XJy4sLAR9g9dy9RzTtQ75XhNHy6\nP2V3hedZdNAgYz3wJ54wFvsW9kGSe1H7809+uno/Dz5odiCiJNDOLkS068XuefuI6zaQBi8NYehf\nn1LqyN5cy02daqwH++uvxRKmKAbSLFOUzp0jZUMos93m8HZ3s4MRJYl2KcX5oU9y4c6HufjWo3R7\nrAtxPe/h1Jg3ue5RI8cVnWbMaEGPHtCmDbSWuY9sntTci9KsWUR6DmDIg645DS4Uosjp0mUI7TqQ\n3XOiSSnrhs99zanx/RskxybctKJT8+bw7bcwcCCcPm125OJ2Sc29COlff+XDU2/Sre5ugoISM/bL\nXO2iuKW4V+bkxKnE3DeBmt+9zriZXxKbdJWYIWPAyTnjvGHDIDraaIdfuxbKlDExaHFbpOZeVCIi\nSDwWQ4RHd/z8fLLUkhISsi81K0TxSKrhxZG3pjNj1P+osG4RPvf6UGHd4ix3Ul97DRo1MtZgTU01\nMVhxWyS5F5Uff2R+pdEEDog1OxIhbnK2Rl0OfP0Px5//glpfvUy7V56CSGMtX6Xg55/hzBl4+mnp\nQWOrJLkXhYQEUv+ayftnRtO9+0WzoxHiluI792fPjB2c6xQAvXsb/SHPnaN0aVi8GDZuhDfeuHX5\nsLBIgoK2ZjxksXfrIcm9KMyezYGqHek60pMyZeR3rbByTs7809CXVV/P5Mi5yyQ18ib6sWeI3rGV\noCCYNQs+/TTnonFxSTfdmBXWQZJ7EdDTfuSj2McZO9bsSISwTEJCKhXr9+DCq3+y75dNVN93mBb3\nDaLqhvmsXKH56iv47DOzoxT5Icm9kO2cPpv4XQcJq9yd06e3EhV12OyQhMiXa16NOfjZInZPmARv\nvEGdh3uy8dvtfPcdvPOOtMHbCknuhazy9N+YXm4iI0ZVlJ4xwqZdaOtvTDgzfDg1Hglke7vHWPn7\nGSZNkgRvCyS5F6YzZ6i8cT3fJI+hd2+zgxGiEDg5wdixsHcvZWtXYs355nj+8T6P3p/ItWtmBydy\nI4OYCtO33/JPxSEEDqtUmIvXC2GKm6YoGPIg/mPGMPb5l4hd2IR3W32Ez+RGVKmSSxl3F/z9WxRH\nuCIbSUGF5epVUr77gbeurubTu80ORojbl3191nXrZhHXtAE8+SruHbfwyIdvcfbBcpx/80s8+rfP\nsUxMzNabrhsWFpmlV418ARQNi5pllFKBSqlopdQ+pdTLtzjnS6XUfqXUdqVUm0z7jyildiilIpRS\nec8/aqumTWNX+c7UDfTAzc3sYIQofOmJu0oVX5z7jCX2nx1sa9WTlm8MwvGRUTidPWnRdaT7ZPHI\nM7krpRyAr4F+gA8wQinVJNs5/YEGWutGwFjgu0yHU4EArXUbrbVfoUVuTRITSfnwY56JeY2hQ8+a\nHY0QxcPRkdjBLdj02142Hvek3uBWdFo5D5V41ezIBJbV3P2A/Vrro1rrZGAmMCjbOYOA6QBa601A\nBaVUtbRjysLXsTnpo/P2PP8m27QPlXrX4syZ/WaHJUSx8mxWjpaL32Vi163EhaXQ+C5vKi37Qyam\nMZklSbcWcDzT9om0fbmdczLTORpYqZTarJR6vKCBWqO4uCSqlveh3uwZvHTlHcaOrSZdH0WJVKYM\nPPWRF3PumciwlL9J/eIrGj/cEfc9O8wOrcQqjhuqnbXWp5VSVTCSfJTWOiSnEydPnpzxPCAggICA\ngGII7/ZUmf0NuxxbU2uwHx4eZkcjhLmaNTuJ/+jhPPF+KM13/sWUt16E8BXw4Yfg5WV2eDYvODiY\n4OBgi861JLmfBOpk2q6dti/7OZ45naO1Pp32Z4xSaj5GM0+eyd0WOF++RNWfP2CYWsd7o82ORgjr\nULEivP+hA6tXP0CLjwfwefSXDGrTFodHR8OkSWaHZ9OyV3qnTJlyy3MtaZbZDDRUSnkppVyA4cCi\nbOcsAkYBKKX8gTit9VmlVFmllFvaflegL7DL8rdi3erP/IUg1yF0erQpFSqYHY0Q1kMpY5LJz6Yd\nYUPvN2mlItkddhnt7U29Wb/KTddikGdy11qnAOOBFcBuYKbWOkopNVYpNSbtnGXAYaXUAeAHYFxa\n8WpAiFIqAggDFmutVxTB+yh++/dTbfkSpjCZYcPMDkYI61S2bCpTp8KMtTV5yvF7BlYMISX8AM2H\nNMZj3jS4nmx2iHbLojZ3rXUQ4J1t3w/ZtsfnUO4wYH9L7WpNyhPjmOr8Mg++VAMXF7MDEsK6NW8O\na9bA4sXe9HtqPr08Ivlg4f/w+eMToh4YDf3aGtV9UWhkhGo+pI+sqxEchPu2kyyocx8/djM7KiFs\ng1LG4tsODrsJCWlP8x9X8WTDlUz84xlYORfeegv69ZMkX0gkuedDXFwS1V3q0XjaNwy6Ppfe/XcC\n9cwOSwib4uRkLMTdv7/izz/70iBqBx82mM8jE57DxaMCTJ4MfftmSfLZpywAmbYgL5Lc80NrPN97\nglnqPnzHd8TVdabZEQlhs9zcjAkn+/TZw44d91Hr26F8UHc2o8ZPJLF0KaKHP86Fth1AKaKiDtOt\n29As5XOat0bcIMn9FnKqKeg/fydx625+azqdj+6BoCCTghPCjpw4sY8uXVJo0cKRuXMDeD7sbp6p\n8R0vfP81qtJfnBo7hS1X5MZrfklyv4X0yY3SuZw4iNfc6Qwq9S+TJpeWZkEhLJB9CmBj3+Es/7du\nTEgGL78Mjz8O77zTh6+2T+DFirMY9+4EHlHJJJZyIS5gMDg65vha0nSTlV3O+VLYHP5LwGvi3byp\nJ3PPO22oXNnsiISwDZlnkkx/5DVFR6VK0KdPJPMXOXK25wgaJ+3mvdRJuH3/MT7DmuIx/0fUtcSb\nymWfbbKkzzgpyT0vWlN7ymOsjvNld9fudOxodkBClAxubvDgg7BgkQMn/Vrgdz2UJxx+JHn2fJoP\nqk+9Wb/CpUtmh2m1JLnnoca3r3Mh/BA/tfmWrt2izQ5HiBLHxQXatTvE7DmKFuO7M7LCMvpcDyIm\n+BSpdevBhAmwd6/ZYVodSe658Pj7a9TsWTzltYRX3ykj7exCmMjREQIC4LvvYNz3LXmvyfe0SN3J\nwuAKJHfqRrtXx1Nh/RKZajiNJPdbqBW0gHLffMD9Hv/w+pdVKFXK7IiEEOkaNoSnnz7G+sO12Tfq\nHZq7HeXjk49Q+sPJNBvciKp/fIrjpVizwzSVJPcc6G++pfq0/2N41TVMmlaPcuXMjkgIkZNKleDF\nF2HPodKUn9CDh5ttZsDFPzkwextN76xPyw9egVWrSmRtXrpCZpaaSsorrxP7zQyGV1/Bqz81ktke\nhbBS2btZli17mA8/akRsrD+Llvnz7NxY7tzxGxNGvUBlx0s4jRkNDz8Mnp45X9DOSM09XUIC1wYO\nY/d3a3mh6ybGf5IkiV0IK5a9m2V6F8tKleCBB+CneZVwf7U/rw+IoG/8HBZ8d5rEpq1J6t4bfv4Z\nLl40+R0ULbusuWcfzJDnQIaICK4OHsmSC/5EjvuLXz4oxYoVx4ohUiFEUVEKmjdP4IUXFNe+9iUo\nyJfHpn+CWr6MJ/f9RbsJz6F69sD5wRFw111Qtuwtr5XvnGIF7DK5Zx9dess5KJKTuf7xZ1x752Oe\nd5xKzZc60dlvFytW3DyKTghhu0qVgkGDYNCgMsTH38OCBffw6e+X8Fi9gEdCptFi1OPEtO1IXPdu\nJPUKoH2fzlnKW5JTrO0LwC6Tu0U2bCD+gSeJPFuN3zqH02N0PI0bN8o4vGXLfhODE0IUlfLlYdQo\nGDWqAvHxD/H++135JNSNKmGLuXf3Avw+/YAr7TviOnIwatBAqF3boutaXKksJiUvue/cSdzEN7ke\nupm3KnxK4OxhTBugCAqSGeaEsDfZb7oeP34MT886WfZVrXqYoZ8MJSnpUTZvfpQvlhzGY1sE/V5a\nQL/nX+d6DU/qtWyN6vcAV1p3QZcqXZxvocBKRnLXGr0+hNg3v8Bh43o+dXqeiyM/ou+gyzg4biMo\nSJphhLBH6Tdd023Zsp+2bbP+P0//le7iAp07Q+PGsfTrN4To6CH8vOw6R2eHU3lZEANWvUHT5F3E\nNOqMS4dmUMsFfHzAwTr7pdh1cne8FEvZZXOo9tsXHLuczHS3p6j81i+0895H7dqNspwrzTBCCMha\n2/f2AW+fUuzc2ZoNrm8xbUMcpTf+S9Pp/+A+YwiVucDllp0pf0cX3MtVQXVqgXa2jnU3bT65Z7+J\n4Xz5EslzFlF+x2RqHFjHKvqwvPYr1Hm6Hb6+V3Bw2EtU1GFq125rYtRCCGuVvbYPkJS0n169wN/f\nHZ4fwvHjdTlQ4QdmLT3NlX82UOXd9XRlBg1Tn+ZUzXb859Mej+aVoZE71K9vytKBNp/c48/GU+dM\nKsnLV1Ih9B88zkexVnflr4YjcH7hT/z7lsdz40wC+99Y31tq6UKI21G6tKZHD+jRowZMHcr160OZ\nNm0PCw7XRm8MpdzmLdRfEUbtL76hnEMCF+r6gm87KvZqS7mOzY35E5ydizRG20ruycnovfu48O8O\n4oLCKBURSteze4hSTdlcrhcxbd/HfUBnzsYtZvDgoXlfTwghCiCnRUiuXTvMHcObwfB+QD9OnNjG\n4aptiQ4+Q8K6rZRev5Xqs/6kmd5NLX2CC5Uacb5qPS7VqYdu6oVu6sWB69epVa9BlusWtEulRcld\nKRUIfI4xovVnrfWHOZzzJdAfSAAe1lpvt7RsFlrDxYskRh3mwpZDXNm6F71zF25Hd+Fx6SDHqEO0\nkw9HqrchvtUrnK3lxt3396FdptGky5fnvhiAEELcjpyabrK3CJQqpenYETp2rA6TBhAWVoeLFwcQ\nEuPMiX2QtPMkqbuiqbf1DDVXraZxyh789Ukulq7JhUoN+a9GA5K9GpJQw5lrjppSjb3Iz7D5PJO7\nUsoB+BroBZwCNiulFmqtozOd0x9ooLVupJTqAHwP+FtSNrMjFdtQ+fJhdKrmEPU5XcqLsxXqEVe7\nK0n9R+HSshZXrh+hf/9B1Eors3z5zEKbJmDLlmDatQsonIvZOPksbpDP4gb5LG7I67PIXrtPX+S7\nalWjkw13t2T58ku07P8sAEfi4eO/ltG8XGOuRx/A6ehBKoQcwCNuP4e++YE66hjawZHYsp5cqeRJ\nco06ObzqDZbU3P2A/VrrowBKqZnAICBzgh4ETAfQWm9SSlVQSlUD6llQNsOKYW+T7OlBhbplqFT5\nOnv33rzi+fLlURaEXDBbt8o/3HTyWdwgn8UN8lnckNdnkVM3zNyULw81veLp0L8h0DBjf0zMVrz7\n+nL2jObk7jgu7jxOQvRxkg8dz/V6liT3WkDmq5zASPh5nVPLwrIZEpom0q2bf8b2tm1y41MIUbLd\nNBCr/DE876lDaarDv7cuV1Q3VAvU7yevhXOFEKKksWQgVk6U1jr3E5TyByZrrQPTtv8H6Mw3RpVS\n3wNrtNZ/p21HA90xmmVyLZvpGrkHIoQQ4iZa6xwr05bU3DcDDZVSXsBpYDgwIts5i4CngL/Tvgzi\ntNZnlVLnLSiba4BCCCHyL8/krrVOUUqNB1ZwoztjlFJqrHFYT9NaL1NK3aGUOoDRFfKR3MoW2bsR\nQggBWNAsI4QQwvZY53RmVkAp9bxSKlUpVcnsWMyilPpIKRWllNqulJqrlCpvdkzFSSkVqJSKVkrt\nU0q9bHY8ZlFK1VZK/auU2q2UilRKPW12TGZTSjkopbYppRaZHcutSHLPgVKqNtAHOGp2LCZbAfho\nrVsD+4FJJsdTbDINwOsH+AAjlFJNzI3KNNeB57TWPkBH4KkS/FmkewbYY3YQuZHknrPPgBfNDsJs\nWutVWuvUtM0wwLIlaexDxuA9rXUykD4Ar8TRWp9Jn05Ea30FiIKMQeIlTlrl7w7gJ7NjyY0k92yU\nUgOB41rrSLNjsTKjgeVmB1GMbjUwr0RTStUFWgObzI3EVOmVP6u+YWlbs0IWEqXUSqBa5l0Yf1Gv\nAa9gNMlkPma3cvksXtVaL04751UgWWv9lwkhCiuhlHID5gDPpNXgSxyl1ADgrNZ6u1IqACvODyUy\nuWut++S0XynVHKgL7FBKKYxmiK1KKT+t9bliDLHY3OqzSKeUehjjJ2jPYgnIepwEMs/MVDttX4mk\nlHLCSOy/a60Xmh2PiToDA5VSdwBlgHJKqela61Emx3UT6QqZC6XUYaCt1vqi2bGYIW265k+Ablrr\nC2bHU5yUUo7AXowZTU8D4cCIkjpOQyk1HTivtX7O7FishVKqO/C81nqg2bHkRNrcc6ex4p9dxeAr\nwA1Ymdbt61uzAyouWusUIH0A3m5gZglO7J2B+4GeSqmItH8LgWbHJXInNXchhLBDUnMXQgg7JMld\nCCHskCR3IYSwQ5LchRDCDklyF0IIOyTJXQgh7JAkdyGEsEOS3IUQwg79P0cY1ferXf6jAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbbb4f5a7b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reals = np.linspace(-5, 5, num=1000)\n",
    "plt.plot(reals, p0.pdf(reals.reshape(-1, 1)), label=r\"$p(x|\\theta_0)$\", color=\"b\")\n",
    "plt.plot(reals, p1.pdf(reals.reshape(-1, 1)), label=r\"$p(x|\\theta_1)$\", color=\"r\")\n",
    "plt.hist(X_true[:, 0], bins=100, normed=True, label=\"data\", alpha=0.2, color=\"b\")\n",
    "plt.xlim(-5, 5)\n",
    "plt.legend(loc=\"best\", prop={'size': 8})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density ratio estimation\n",
    "\n",
    "The density ratio $r(x)$ can be approximated using calibrated classifiers, either directly by learning to classify $x \\sim p_0$ from $x \\sim p_1$, or by decomposing the ratio of the two mixtures as pairs of simpler density ratios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "\n",
    "from carl.ratios import ClassifierRatio\n",
    "from carl.ratios import DecomposedRatio\n",
    "from carl.learning import CalibratedClassifierCV\n",
    "\n",
    "# Classifier\n",
    "# from sklearn.linear_model import ElasticNetCV\n",
    "# clf = ElasticNetCV()\n",
    "\n",
    "from sklearn.neural_network import MLPRegressor\n",
    "clf = MLPRegressor(activation=\"tanh\", hidden_layer_sizes=(10, 10), random_state=0)\n",
    "\n",
    "# from sklearn.ensemble import ExtraTreesRegressor\n",
    "# clf = ExtraTreesRegressor(n_estimators=250, max_leaf_nodes=15)\n",
    "\n",
    "# No calibration\n",
    "cc_none = ClassifierRatio(base_estimator=clf)\n",
    "cc_none.fit(numerator=p0, denominator=p1, n_samples=100000)\n",
    "\n",
    "# Calibration + Direct approximation \n",
    "cv = StratifiedShuffleSplit(n_iter=1, test_size=0.75, random_state=0)  # 25% for training, 75% for calibration\n",
    "cc_direct = ClassifierRatio(base_estimator=CalibratedClassifierCV(clf, cv=cv))\n",
    "cc_direct.fit(numerator=p0, denominator=p1, n_samples=100000)\n",
    "\n",
    "# Calibration + Decomposition of the mixture\n",
    "cc_decomposed = DecomposedRatio(ClassifierRatio(base_estimator=CalibratedClassifierCV(clf, cv=cv)))\n",
    "cc_decomposed.fit(numerator=p0, denominator=p1, n_samples=100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: `CalibratedClassifierRatio` takes three arguments for controlling its execution:\n",
    "- `base_estimator` specifying the classifier to be used,\n",
    "- `calibration` specifying the calibration algorithm (`\"kde\"`, `\"histogram\"`, or a user-defined distribution-like object),\n",
    "- `cv` specifying how to allocate data for training and calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD7CAYAAACPDORaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8TNf7wPHPyR6y2RLEElvt+75F0JaiaHVBKVqlP9Xq\n+i0tRReqrdIVbbWoqlZbWqWWIpba9y1KQkJCgojs28yc3x830iAiZLJNnvfrdV/J3Dn33DOTm2fO\nnHvuc5XWGiGEELbJrrAbIIQQIv9IkBdCCBsmQV4IIWyYBHkhhLBhEuSFEMKGSZAXQggb5lDYDbiR\nUkrmdAohxB3SWqvs1hfJnrzWulCXyZMnF3obisoi74W8F/JeFP33IidFMsgLIYSwDgnyQghhw9Tt\nuvoFTSmls7bJz8+PsLCwQmyRKO6qV69OaGhonusJDAwkICAgz/XYAnkv/lMU3gulFPoWY/JFPshn\nNL4QWySKOzmGhK3LKcjLcI0QQtgwCfJAXFwcXbt2pWvXrnh5edGtWzeefvrpfN3nmTNn2Lx5MwDn\nz5/ngw8+yNf9CSFKJhmuuYG/vz9btmyxWn1aa5S6+VvUhg0b2LZtG5MnT7bavkT2ZLhG2DoZrrkL\nkyZNYuTIkfTo0YNffvmFqVOnAjB//nyWLFkCwOTJk+natSv33Xcf4eHh120/dOhQXnjhBfr06UNE\nRARdu3alc+fOvPjiiwB89dVXLFiwgJ49exISEsKIESMAWLBgAe3bt8ff359jx44V4CsWQtgiCfI5\nqF+/PmvXrqVMmTI39cYPHjxIdHQ0mzZtYtasWbz//vs3bd+lSxdWrVqFt7c3GzduZOvWrURFRREW\nFsaoUaMYMWIEa9asAYxPYpPJxLx589i+fTsLFizgzTffLJDXKYSwXUUurUFuZDP6kaO7/abesmXL\njP39t8NrX/uDgoLYsGED3bp1A6BKlSq33P7SpUuMGTOG2NhYQkNDOX/+fLb7i4qKokaNGiilqFmz\nJjExMXfXcCGEyFAsg3x+Dq9mHbu1szO+6Hh6enLhwgUAjhw5Qtu2balbty69evVi5syZAJjN5pvq\nurb94sWLeeyxxxg8eDADBw5Ea42joyMmk+m68t7e3pw5cwaz2UxoaChlypTJl9cohCg5ZLjmBtd6\n7Vl7782bN+fMmTP07t2b6OhoAFq0aJE5E+fee+9l0aJF2dYD0L17d6ZPn86AAQNISUkBoEmTJmze\nvJkhQ4ZklnN0dGTUqFF06tSJ4cOH88477+Tb6xRClAwyu0bYPDmGhK2T2TVCCFFCSZAXQggbJkFe\nCCFsmFWCvFKqp1LqhFLqpFLq9RzKtVZKpSulHrbGfoUQQuQsz0FeKWUHfA70ABoCg5RS9W5R7n1g\nbV73KYQQInes0ZNvA5zSWodprdOBpUC/bMo9D/wCXLTCPq1q8+bNeHh4EBcXB8CIESM4ffq0Vepe\nuHAh3377LVFRUUyfPh2Azp0733E9hw4d4uDBg5m/f/fdd1ZpnxDCtlkjyPsC57I8Ds9Yl0kpVRno\nr7WeA9zh9aoFo1q1anzzzTf5Vr+Pjw8TJkwAyDZh2TW3mup38OBB9u3bB0DTpk0zc90IIUROCuqK\n19lA1rH6HAP9lClT8rUx2enbty8rV67MTCAGRgriJ554gvj4eJo1a8bs2bOv22bFihXMmDGDUqVK\nMWXKFDw8PBg3bhypqan069eP8ePHZ5YNCwtj4sSJfP/996SnpzNixAiOHTvGlClT6NWrF127dqVN\nmzZcuHCBV1555aZ6vvrqK65cuUJgYCAjR47k77//5p133mHGjBmsXLkSFxcXFixYQJUqVWjSpAmN\nGzfm+PHjLFy4kCZNmhTY+yiEyH+BgYEEBgbmrrAV7hLeDliT5fF44PUbypzOWM4A8UAk0PcW9ems\nbnycHwIDA/XEiRP1zJkz9bJly/SIESN0SEiI/vDDD/XixYu11lo/88wzevfu3ZnbWCwW3bp1a52a\nmpq5LiUlJfP3rl276pSUFL1gwQI9f/58HRoaqocOHaq11rpWrVo6PDxcJyYm6o4dO2qttQ4ICNA7\nd+68bT3X2jtp0iQdGRmpe/bsqbXWetu2bXrMmDFaa60rVqyo09LS9D///KNfeuklq79fxU1BHENC\nFKaMYzzbGG2NnvweoLZSqjpwARgIDLrhg6Tmtd+VUt8BK7XWf9ztDtXUOxvx0ZNvf7WjUoqRI0fy\n6KOP4uvri9aakJAQevfuDRjJxoKDg2ndujVgJB3z8/PDyckps47Tp0/zyiuvkJSUxMmTJ7l4MfvT\nD+XKlcPX1xjRcnD4709wLaFZbusJDQ3N7KW3atWKt99+G4DatWvj6OiIr68vsbGxt33tQgjblecg\nr7U2K6XGAuswxvjna62DlFKjjaf1Vzdukud95iJo3w0PDw/q1q3Lxo0bUUpRu3Zt9u7dS/369dm7\ndy/PPPNMZtkKFSpw9uxZUlNTcXZ2RmvNnDlzGD9+PP7+/nTu3PmW4+tXrlzh/PnzeHl5XZfY7FpC\ns+zqcXR0JDU19bp6/Pz8OHToEAB79uyhVq1aN+3rVm0QQpQMVhmT11qvAeresG7eLco+ZY195pcX\nXniBL7/8EoCRI0cyePBgvvnmG5o0aUKbNm0yyymlmDBhAl26dMHNzY3JkyfTu3dvnnvuORo0aICz\ns/Mt91GhQgWmTJnCwYMHM+8MlfVkbHb1tGvXjuHDh3Ps2DEGDBgAGCdzAwIC6NixI87OzixcuPCm\nuoQQJZskKBM2T44hYeskQZkQQpRQEuSFEMKGSZAXQggbJkFeCCFsmAR5IYSwYRLks9iwYQNdu3Yl\nICCAAQMGEBMTc1OZa8nL1q5dy19//UVYWBhDhw694339/vvvXL16FTCSmB04cCDP7RdCiBsVVO6a\nIu/y5cu88847rF69mlKlShEcHHzTxUdZ9ejRAzBy0twu4Vh2z69YsYJGjRrh5eXFsGHD8v4CRImQ\nng4bNkBQEHh5wX33QZUqxnF2Lu4c+y/s52T0ScLjwgmPCyc+LZ7k9GTSLem4OrhS2qk0ns6eVPGo\nQjXPavh5+dHEpwlVParK9RU2SoJ8htWrV/Pkk09SqlQpwEgNsGDBAhYuXEhiYiLTpk3j3nvvzSy/\ncOFCzGYz3bt35/Tp0/Tr149Lly7x448/Ur16dZo2bUrTpk1p1KgR3t7e19VTt25d1qxZw4kTJ3jk\nkUdISEigc+fOdOnShaFDhxIREUGVKlX4/vvv2bp1KzNmzMDBwYGYmBjWrl2b2UZRsqxYAc8/D9Wq\nQcuWcPZKJM/PX4NPp7+ILbsRBzt7WlZuSf3y9alVphb+1f3xdPbExcEFR3tHktOTSUxPJDYllnNx\n5zh+6Th/nvyTw1GHSTGl0KxiMzpU7UBXv660r9qeUo5ynNkCCfIZLly4cFO2xoEDBzJ8+HDi4uJ4\n9NFHrwvyWcXExLB161b27NnD+++/z5w5czh//jy7du3CxcWFlJSU6+pZu3YtDzzwAJMmTaJGjRpM\nnToVgOXLl9OwYUOWLFnCtGnT+PXXX/H29sbZ2Znly5czffp0NmzYwIMPPpjv74coOrSGGTNg7lz4\ndlEKMd5/8N3B79gZvpP7et3H5Z0P4LBhJpt+r0Llyne3j6iEKA5EHmBr2FYmbZrE4ajDtPFtQ/96\n/elfrz/VPKtZ90WJAlM8g/ydfq3MxdWOlSpVIiIi4rp1f/31F59++ilaay5dunTLbRs3boxSimbN\nmhESEgJA3bp1cXFxuWU9+r+sm5lCQkJo0aIFYCQr279/P97e3jRq1AiAypUrZ47ji5Jj1ixY8GM8\nQ+fMZdj+WTSo0IARzUbw62O/UsqxFHoQvPce9OwJW7YYwzh3ysfNh561e9Kzdk8AEtIS2HB6A8tP\nLOftzW/j5+XH4w0fZ0iTIVRyr2TlVyjyU/EM8vlwiXqvXr149NFHGThwIG5ubgQHB/P6669z5MgR\nUlJS6NSp0y23PXLkCBaLhYMHD2YmCcs6vvn++++zZcuW6+pxdHS8LjkZkJkQ7YEHHmDv3r3UqVPn\nprrk8vySZeXqNKas/RzHQdMJTrqXv574i6YVm15XRil48024eBEGD4ZVq+68H3QjNyc3+tXrR796\n/TBZTGwJ28IPh3+gwZcNaF+lPU82fZL+9frj4uCStx2JfFc8g3w+KF++PJMmTaJPnz4AlC1blqFD\nh9K5c2fatGlDmTJlgOyTf/n4+NC/f38uX77MDz/8cFO5Bx98EH9/f1q3bo1XRjerR48ejBkzhkcf\nfTSzbL9+/Vi2bBldunShcuXKjB8/nm3btuXr6xZF19J9qxmy/iVa9anFdwO3UL9C/VuWVQpmzoRO\nneCzz+CFF6zXDgc7B7rV6Ea3Gt34rNdnLA9azvwD8xm3Zhwjm4/k2VbPUtWzqvV2KKxKEpQJm1fc\njqHYlFjGrRnHr3u20j3tM1Z82CvX2wYHQ9u2sG8f+PnlXxsB/r38L1/u+ZLvD39PtxrdeL7N8/hX\n95dZOoUgpwRlEuSFzSvqx9CuuDgWRUbyYpUqRFzczbAVw2jh3ps9733AyaNu5DSZKiY9nZ8uXiQ0\nJYXpNWuilOLdd40gv3x5wbQ/PjWeRYcW8enuT/Fy8WJCpwn0rdsXOyWX4RQUCfKiRCvqx1Drffuo\n6OjIxitRmK8e5v/8mrBxUgfGPebKU0OvD5TxJhNHExPZGx/PX1eusD02lh5ly7I/IYHPatemZ7ly\npKRAo0Ywbx50715wr8NsMbPixAqmb5tOsimZ8R3HM7DRQBztHQuuESWUBHlRohXlY+hYYiL3HzpI\nh4hPOR0bwZBu8/jtpGJPQhy6QioeDg64ZNwxLCY9HQvQoFQpmrm5cX/ZsvQoWxZPBwd+vXSJ98+e\nZXeLFiilWLIEvvgCtm3L+0nYO6W1Zv3p9UzfNp3Qq6FM6DSBEc1GSLDPRxLkRYlWlI+hF/89yk/H\nf6WHOsPcPnNxsnOhcWP4+GO4737NpfR00iwWNFDWwYHS9vbZjnlbtKbu7t0sqleP9p6emM3QpIlx\nMrZnz4J/XdfsOLeDyYGTCYkJYUqXKQxuPBh7O/vCa5CNkiAvSrSiegwFRwdTf+9+hjue4auu/0Mp\nxe+/w9tvw969d94D//jcOfbHx7O4QQMAli2DDz+EXbsKvjd/o82hm5m4aSLRSdFMDZjKgAYDZMze\niuTOULmQm+RkUDwTlI0YMYL27dtz//3388gjj2ResFUUvfzyy2it2bx5M6GhoQCZ77UtCboURMcf\nH6GUS9nMAK81TJsGb7xxd0F5eMWKrIyOJsFkAmDAAIiNha1brdz4u9DFrwtbhm9hVo9ZfLD9A1rM\na8HqU6uL5Ievzbl25WVRWYwm/efGx/nh0qVLukuXLjoxMVFrrfWpU6d0ZGRktmWHDx+uQ0JCMh+H\nhobqoUOH3rJui8Vyy3qCg4Pz0OrrTZkyRYeFhd1yX9fafOzYMe3v72+1/eaXKVOm6L///tsqdRXE\nMXQnjl88rivPrKyf2f2rHnjsWOb6zZu1vucerc3mu6/7/oMH9c9RUZmP58zRum/fvLTW+iwWi14e\ntFzX/ayu7r6wu95/fn9hN6nYyzjGs42p0pMn++RkPj4+LFiwgK5du9KmTRv+/vvv67ZZuHAh3377\nLUBmgrIOHToQFhYGQNOmTXnyySf58MMPb6rn3LlzrFmzhiFDhjBz5kymTp3Kxo0bMZvNDB48mC5d\nuvDEE09gsVjYvHkzvXr1om/fvnTu3JmkpKRsX0Nu5yY3aNCAmjVrcu7cOS5fvky/fv3o3r07Y8eO\nBSAqKopevXrRrVs33nzzTQBmzJhBp06duPfeewkPDwegSZMmDB06lKZNm7J06VL69OlD69atOX/+\nPGCkenj88cdp3bo1+/btA2Dx4sW0b9+ezp07c+TIEdLT0+nTpw/dunVj4MCBAHTt2pX09HQWLFjA\nq6++ymuvvXbdez1u3Di6dOlC3759iY+PJywsDH9/fx555JHr9l9Unbh8gnu/v5f3u79Pqlt9umbJ\nQfD550YCMrs8/FcOqFCBX7Ok4HjySdi+HU6ezEurrUspRf96/Tnyf0cYUH8AvZb04snlT3I29mxh\nN8023Sr6F9ZCIfTk33//fb169eqb1icnJ2uttY6NjdX333+/1vq/XvGCBQv0/PnzdWhoqG7QoIG2\nWCx6165d+tlnn9Vaa12+fPnM7bOrZ8SIEfr06dNaa6PXumHDBr1s2TL97rvvaq21fu+99/TPP/+s\nAwMDdf/+/bXWWk+bNk3/8ccf2b6GKVOm6NDQ0Gyfu/Hbx4QJE/SuXbv0K6+8onfu3Km11vr111/X\nO3bs0OPGjdPr16/PLBsZGal79uyptdZ627ZtesyYMVprrStWrKjT0tL09u3bdcuWLbXWWi9ZskTP\nmjVLa611uXLldFJSko6IiNB9+vTRZrNZt2nTRpvNZh0aGqr79OmjT58+rUeMGHFdW7t27arNZnPm\ne6K1znyv9+zZo0eOHKm11vqHH37QM2bM0KGhobpFixaZ+//0009vev0FcQzlxpmYM7rKx1X0ggML\ntNZaV9++XZ/I+PYYEaF1mTJax8bmbR+Rqanac8sWnZbl68Cbb2o9dmze6s1PcSlxeuKGibrsjLJ6\n/Prx+mry1cJuUrFDDj35YpnWQAUG3lF5HRCQ4/PZJSeD4pGgbMqUKWzevJmwsDDWrl2Li4sL06dP\np23btrdsb0REBJUrVyYoKIjx48cDkJiYSNu2bTl16hQdOnTILBsaGpqZnbNVq1a8/fbbgPFtx9HR\nkcqVK1O/fv3M9p04cQKAWrVq4erqiqurK3FxcVy+fJnq1atjZ2dH9erViYuLo0aNGjRu3JihQ4fS\nsmVLXnzxxVu2Obv3Z/PmzYDx7eTa/oODg3Oso7BEJ0XTc3FPXu/4OsOaDSM6PZ0Yk4l7XF0B+Ppr\nGDgQPDzyth8fJydqurqyMy6OzhnfEkaNgubNjUyWRTFLtbuzO+90e4dnWz3LpE2TuOfze5jkP4nR\nLUfLtEsrKJZB/nZB+07dmJwsJCSE0qVLZ5tYLDuFmaBsypQpALz99tsMHz6catWyTwl7bbugoCDO\nnj1LlSpVqFevHkOGDKF58+YAmM1mtm7dyo4dO+jevTtaa/z8/Dh06BAAe/bsyXx9WWXXvuDgYJKT\nk4mJicHDw4Py5ctz9uxZTCYT4eHheHp6kp6ezosvvohSih49ejB48ODM7R0dHTFlnEC8platWqxb\ntw6AvXv3ZtuWG9+foiA5PZkHf3yQ/vX6M7aNMSx2KCGBxqVLo5QiPR2++grW/ZkG56IgMhKiMn5G\nRsLVq+DgYCyOjkaayfLloVw542elSuDjYzwP9ChblnUxMZlBvlo1aN8efv4Zhg8vrHfh9nw9fPm2\n37ccjjrMy2tfZs7eOXzS8xPurZl9im+RO8UyyFtbdsnJ5s+fT58+fTITixX3BGVDhgzB3d0dd3d3\n5s+fD8CECRMYNWoUsbGx2Nvb88033/D6668zbNgw3nvvPTp06MC7775LQEAAHTt2xNnZmYULF97y\nfciqatWqPPXUUwQHB/Pll19iZ2fHc889R+fOnbG3t+eLL74gNDSUp556CrPZTO3atfH29s6sNyAg\ngDfeeIPdu3dnfnC1atUKV1dX/P398fDwYMmSJcTExBTpXClmi5nBvw2mjmcNpjn0gIULITKSw6VK\n0dTODl59leRTkQRdjMSjbQJUqAAVKxqLj89/P81mMJkgLQ1CQox5kdHRcOkSXLgAly8b21auTI82\nbfhf9+68s3gxVK4Mvr680qMKb85ryPDhRfe9uqaJTxPWD13P7//+zug/R9PYuzEz759JrbI3f6iL\n25N58iJfdO7cma1FYe4ehXsMvbnhTXac2cLfK8tgFxwCLVpApUoMaN+eB4HhHh48/64P7R+qyOCx\nZe/+rKvJZPT+IyJIO3+eCu7uhGzdSvmwMDh/Hn3iBKsvtsYvcAEN27pZ9TXmpxRTCrN3zuaj7R/x\nTItneKPzG7g7uxd2s4ocuRhKFDh/f3+2bNlS2M0ACu8YWnZsGf9b+ypBO1vhkpQKv/0GTk5cSU+n\n5s6dhLZrR1qMI/fcA+fOgbsVY9fwoCBquLoy+VoqypQU9rcbg0/EPnyPrQdvb+vtrACcjz/PhA0T\n+Pv030zvPp0hTYbIxVRZSJAXJVphHEOHIg9x7/f3cvDqYHw37IJNmyDjJOuXERFsuXqVpQ0bMns2\nHDhgjOJY0+nkZNrs28eJNm0o7+QEQEiwZnnTKbxS/RfUpo3GMFAxszN8J+PWjAPg056f0rbKrScY\nlCRyxasQBSg6KZqHfnqIHyqNxXf+T0Z+gYwAD7AoMpInK1YEYMECGDbM+m2o6erKQG9vPjx3LnNd\nrdqK5c2mcqrZo9CjByQkWH/H+axdlXbseHoHz7V+jod/fphhK4YRmRBZ2M0q0iTIC2FFWmuGrRjG\nI/f05/4Pf4NPPoGq/901KTwlhVPJydxfpgwHD0JMDFh5slimV6tWZf6FC5lpDgCGDoVJpsnQsqXx\noBh+S7ZTdjzZ9ElOPHeCiqUr0nhOY2bvnI3JYrr9xiVQkQ/y1atXRykliyx3vVSvXr3AjtdZO2dx\nOeky009UMca9H3vsuudXRkfTq1w5HOzsWLTIuCI1L1e45sTP1ZUALy8WRkVlrnvsMVi7ThH7/hyI\niIAvv8yfnRcAd2d3Ztw3g60jtrLq1Cqaz2vO5tDNhd2soudWV0kV1kIRuTpRiDu1O3y3rvBBBR0a\nfkxrb2+tDx++qcwDhw7pn6KitNmsdeXKWh8/nr9t2hoTo2vv3KlNWXIoPfSQ1t98o7X+91+ty5fX\n+sSJ/G1EAbBYLHrZsWW66sdV9aBfBumIuIjCblKBQnLXCJG/YlNiGfjrQOb0nkP1ZeuMO2o3bnxd\nmej09Mw7OW3bZkxrr3/re3NbRUdPT8o4OPDH5cuZ64YOhe+/B+65x0h5OW5csRy2yUopxSMNHiHo\nuSD8vPxoMqcJH23/iHRzemE3rdBJkBcij7TWjP5zND1q9WBAvYeMTGOvvXZTuR+iouhdrhyeDg4s\nXWqkMchvSikmVa/Om2fOkG6xANCrFxw9CqGhwNixcPYsrFyZ/40pAKWdSjOt+zS2P72dDWc20HRu\nUzae2VjYzSpUVgnySqmeSqkTSqmTSqnXs3l+sFLqUMayTSnVOLt6hCiOfjjyA8cuHePjHh8bUyVL\nl4ZscgctjopiWMWKmEzwyy83Ddfnmz7lylHN2ZnPM/IzOTsbueZ/+gkjTcKHH8LEicW+N5/VPeXu\nYfXg1UzrPo2n/3iax395nPC48MJuVqHIc5BXStkBnwM9gIbAIKVUvRuKnQb8tdZNgXeBr/O6XyGK\ngvC4cF5e+zKL+i/CxcEFvvsORo686a4fIcnJhKWk0M3Li40boUYNqFmzYNqolOLDWrWYcfYsSRn5\nkh5/PCPIg9G1t7eHVasKpkEFRCkjpfGxMceoV64ezeY24/1t75NmTivsphUoa/Tk2wCntNZhWut0\nYCnQL2sBrfVOrXVsxsOdgK8V9itEodJa8/QfT/N8m+dpXqm5kVdm1Sp49NGbyv588SKPVKiAg50d\nP/1UMEM1WTV2c6OjpydfZeTb79LFSHlz8iTGB9KECcZtqWxQKcdSTO06lV0jd/HPuX9oPKcx60LW\nFXazCow1grwvcC7L43ByDuIjAdu6l5sokebunUtMcgwTOk8wVmzaBA0aGEnFbrA+Jobe5cqRmgor\nVhTcUE1Wr1Stypzz59FaY28PjzySpTc/YIAxpTIfb0NZ2GqVrcXKQSuZef9M/m/V/zHg5wEl4kYl\nBXriVSnVFRgB3DRuL0RxEnwlmEmbJrHooUU42GUkc/39d+jf/6ayJouFvfHxtPfwYN06aNgQfAvh\nu2x7Dw8clGJLrPGl+rohG3t7eOop+Oabgm9YAetzTx+OjTlGU5+mtJjXgve2vEeqKbWwm5VvrJFq\nOALImsS8Ssa66yilmgBfAT211tnfJTvDtRzpYKScDcivSwKFuAtmi5nhK4Yz0X8i9cpnOf20YYOR\ntP0GRxITqersTBlHx0IZqrlGKcXISpX49sIFunh50aGDcaPvo0ehUSOMIN+0KXz8sXF21oa5OLjw\nVpe3GNpkKC+ve5lGcxrxSc9P6FWnV2E3LVcCAwMJzO3Nk241gT63C2APBAPVASfgIFD/hjLVgFNA\nu1zUl0+XCwhhHV/u/lJ3mN9Bmy1Z7rh9/rxx/75s7sL9eXi4HnnihE5M1NrTU+ss99kucOdTUrTX\n1q06yWTSWmv98staT5yYpUDnzlr/+WfhNK4QrT65Wtf5tI5+cMmDOjg6uLCbc8fIz4uhtNZmYCyw\nDjgGLNVaBymlRiulRmUUmwSUBb5USh1QSu3O636FKAyRCZG8FfgW8/rMuz7V7dat0LlztjkKtsfG\n0t7Dg9WroU2bws3yW8nZmRZubqy+cgUwvlX89FOW2ZMDBhjzO0uYB+o8wJH/O0KHqh1o+01b3tr0\nFknpSYXdLKuwypi81nqN1rqu1rqO1vr9jHXztNZfZfz+jNa6nNa6hda6uda6jTX2K0RBe3ntyzzd\n/GkaeTe6/oktW8DfP9tttsfF0cHDg59+MsbBC9sgb29+zMhn06qVcdOpgwcznnz4YfjjD2NlCePs\n4Mz4TuM5MPoAJ6NP0uCLBvwW9FuxT3UuV7wKkUvrQtaxI3wHk/wn3fzk3r1GN/0GF1JTiTOZqGwp\nxbp18NBDBdDQ2xhQoQLrY2KIM5lQypjps3RpxpNVqxq3DNy7t1DbWJiqelZl6SNL+a7fd7y16S16\nLO7BicsnCrtZd02CvBC5kJyezJhVY/j8gc8p7VT6+ifNZuPsZZMmN223Iy7OGKpZpejYEcqWLaAG\n56CMoyNdvLxYnpHPZuBA43xxZof1vvtg/frCa2AR0bVGVw6MPkCvOr3o9G0n/rf+f8Snxhd2s+6Y\nBHkhcmH6tuk0q9iM3vf0vvnJ4GBjoN3T86antsfG0t7Ts8gM1VyTdcimSRNjMs3ua2fKJMhncrR3\n5MV2L3LCMgGIAAAgAElEQVR0zFEuJl6k/hf1WXJkSbEawpEgL8RtnLh8gi/3fMnsnrOzL3DokDH1\nMBvb4+Jo6uDBxo3Qr1+2RQpF3/Ll2REXR2RqKkoZH0CZQzb+/sZwTXJyobaxKKnoVpEF/Rfw86M/\n89H2j+iyoAuHow4XdrNyRYK8EDnQWvN/q/6PSf6TqOJRJftCBw9Cs2Y3rT6UkMDJpCQiAz3o0gW8\nvPK5sXegtL09L1Spwr2HDhGanMzjjxt3KbRYMBKs1atn01e/3q0OVTuw55k9PNH4Ce5ddC/Pr36e\nmOQcL/vJX1rDsWM5FpEgL0QOFh9eTGxKLM+1ee7WhY4dy7ia6D8pZjNDgoL4uHZtfl9qXyhpDG7n\nbT8/nqlcmQ4HDpBUJY6yZWHbtown27WDXbsKtX1Flb2dPaNbjSbouSDSLenU/6I+8/fPx6ItBd+Y\nXbtuOvZuJEFeiFu4knyF19a/xrw+8/5LXZCd06ehVq3rVi25eBFfJycedPZhyxbo2zefG3sXlFKM\nq1KFyX5+TDxzhkGD4McfM55s2xZ27izU9hV15UqVY26fuawavIpvDnxD+/nt2ROxp2AbkYtvWxLk\nhbiF8X+P55EGj9Dat/WtC2ltBPkb8gb/HRPD497e/P67ols38PDI58bmwWMVKvBPXByPDdT88ouR\nTFOCfO61rNySf576hzGtxtB3aV+e+eMZLiVeKpid790LZcrkWESCvBDZ+OfsP6w6tYr3ur2Xc8GL\nF6FUqZuieFhKCrVdXfn558LJOHknyjg64uPoSJp3EnXrZkysueceI7HNpQIKVsWcnbJjWLNhnHju\nBG5ObjT8siFf7P4Ck8WUvzveswfGjMm5bfnbAiGKn3RzOs+uepaP7/8YT5ebp0VeJ5tePMDZ1FQ8\nkl345x948MF8aqgVtXB3Z39CAoMHw5IlGDnmW7SAffsKu2nFiqeLJ7N6zmLjsI0sO76MVl+1YtvZ\nbbff8G4kJRnTd0eOzLGYBHkhbjBr5yx83X15rGEuuuDZBPl0i4WotDR2rXbivvvAzS2fGmpFLdzc\n2B8fz6OPGvc9SUgAWrYs0Ve+5kUj70ZsGraJCZ0mMOjXQTzx2xPWv/3g7t3GzeL9/HIsJkFeiCxC\nr4bywT8f8EWvL1A33MIvW9kE+YjUVCo6OfHLT3ZF6gKonLRwd+dAQgIVKkCHDkb6Glq1kp58Hiil\neLzR4wQ9F0RNr5o0nduUdza/Q3K6la4/2LzZuMXXbUiQFyKD1prn/3qel9q9RK2ytW6/AUBIyE1B\n/mxqKpXsnNm927h9anHQPKMnr7X+b8imZUsJ8lbg5uTGO93eYe8zezl88TD1v6jPsmPL8n7VrAR5\nIe7MihMrCLkSwmsdX8v9RkePGrd6yuJsSgqm8y707WtcV1QcVHBywt3BgTMpKfTrZ2ROvuxZC+Li\n5OSrldQoU4Nljy5jQf8FvLv1XQIWBnAw8uDtN8xOaqoxXNOp022LSpAXAohPjeeFNS8wp/ccnOyd\ncrdRcjIEBRnjolmcTU3l/AFnhgzJh4bmo2vj8u7u8MAD8MuvyujN//NPYTfNpgT4BbB/1H4GNxpM\nj8U9GLVyFBcTL95ZJXv2QN262eZLupEEeSGAyYGTubfmvXTxu/3X30xr1xrzyW/orh+OTCEp1IXu\n3a3cyHx2bYYN8N+QzdChMGdO4TbMBl27avbEcyco7ViaBl804OMdH5NmTstdBZs2QS5viypBXpR4\nBy4c4IcjP/DhfR/e2Ya//mrcSekGe8+l0q2hM/b2VmpgAbnWkwfo2dP4khLSZhAcOWIswurKuJZh\nVs9ZbB2xlXUh62g8pzGrT62+/YYrV0LvbDKiZkOCvCjRzBYzo/8czfTu0ylfqvwdbGiGv/66KV+B\nxWKMyQ+618XKLc1/Ldzd2ZeQgNaapVci8X3/FN/+6ARjx8KsWYXdPJtWv0J9/nriLz6+/2NeWvsS\nvX7odesblVy4AKdOGbebzAUJ8qJEm7dvHs4OzgxvNvzONty5E3x9jTspZXjrzBk6bD2MuXwq9zVx\ntm5DC0BlJyfsgd5HjjD97Fni68TwZegFzCNHw4oVEBlZ2E20aUopet/TmyP/d4R7a95Lp2878dKa\nl7iacvX6gqtWGV+1HB1zVa8EeVFiRSZEMjlwMnN7z73+pty5sXr1dV+X0y0WZoeHcz4Cqqe5USaX\n/4BFiVKKJytWpJW7O/tbtuSX5vVJeOwMK/Z7GrkZvvmmsJtYIjjZO/Fy+5c5/txxEtMTqft5Xebt\nnfdfioQ//rijy6hVUbvDiVJKF7U2Cds06NdB1PCqwbTu0+5842bN4IsvoGNHALZcvcq4k8GE9WzF\nwSMWqvnaRv+p0/J/uRTqwL8BccYNakNCKHYnG4q5AxcO8OLaF7mSfIVP2k6hW9enICzsuhsUKKXQ\nWmd79Z5tHIlC3KF1IevYFb6Lif4T73zj8HBjadcuc9WaK1fwjShL587YTIAH+KaTH6dqX+B4pYZQ\nvrzcFrAQNK/UnMBhgbwd8DYbPvg/ttUrRVD6hVxvbztHoxC5lHlT7l6fU8qx1J1XsHatcR/ULD3a\nNVeucOansjzzjBUbWgTUq+BMzfPleWNDFIwaBV9/XdhNKpGUUjxU/yHeiajL5QEP4L/An7Grx3I5\n6fJtt5UgL0qcaVun0bxSc3rVucucA4GB0K1b5sMEk4mghCSu7vCgZ0/rtLEoebpeedbHR2N5fBBs\n3GjM7hAF7/Rp7E6coP+Lcwl6LgiFov4X9Zm5fWaOm0mQFyXK8UvHmbtvLp/0/OTuKtDaCPJZLkQJ\nSkrC7Wopnh5mh0MON5Aqrl7oWIYUv3iW73SBhx+GxYsLu0kl0xdfwIgR4ORE+VLl+azXZ2wdsZVN\noZty3EyCvCgxLNrCqJWjmBowlcrule+ukjNnjMnwtWtnrtofk0jc4VI89ZSVGlrElHawp5HZk3fX\nXYEnnshyj0BRYBISYMGCm24QUq98Pf4c/GeOm0qQFyXG/P3zMVlMPNvq2buvJDDQyPyXJQ3xL/sS\nqaXcbpfWu1gb2bA8QV7RnKrcxZgvHxRU2E0qWRYtAn//2+aOz44EeVEiRCZE8ubGN/nqwa/ufE58\nVjcM1VgssD0ykaGdikm6ybv0cKVy2LW9widzFAwcKL35gpSWBh98AK/dQXbULCTIixLhpbUv8VTz\np2ji0+TuK9HayOGdJcivWwdpvokMbmfbQd7X2Zna7i4s3BdHXO9BRpCX61kKxsKFRsbJDh3uanMJ\n8sLmrQlew67wXbzV5a28VRQRYaQXrlMnc9UHc9NxdDdTzaX4pTG4Uw9XLIf7mDC6XHKj5rRp/LRj\nR2E3qWSYMwdef/2uN5cgL2xaUnoSY1aNYU7vOXc3Jz6r3buhTZvM8fiB20+y+and1HZzyd2tAou5\n0ZUr076WM6dWufNBWATPx8dzNiWlsJtl+2JjoXr1u95cgrywaZM3TaZdlXb0qN0j75Xt2WMEeSA6\nPZ3f4qN4Iro2s+vUvs2GtqGSszO/tq1HX3MVkkI78PLq1Yz899+838ZO5CwtDZxyeSObbEiQFzZr\nZ/hOFh9ZfPdz4m90rScPfHLgEvb7yzL3cR+6lSljnfqLiTfegNeXteKlX38jIj6e3y/f/qpLkQcS\n5IW4WYophRG/j+CTnp9QoXSFvFdoscDevdC6NQDzj0bzkE95SuVxBKg4atQIut9nx+GKvZl99Civ\nnT5NusVS2M2yXRLkhbjZlMApNKzQkEcbPGqdCv/910jQVa4cO/ZauFDxKtMfLlk9+Kzeew8++Lc/\n/vO/p5qzMwsl13z+kSAvxPX2ROzhu4Pf8UWvL6x3QjRjPF5rGP3FVWqo0lT3uvt/vOKuenW4Z3RX\nzIeP8Z6XF2+HhZFsNhd2s2xTUQjySqmeSqkTSqmTSqls5/oopT5VSp1SSh1USjWzxn6FuFGqKZUR\nv49gdo/Z+Lj5WK/ijPH4FSvgfLUrDKtX1np1F1P/m+TMBnUfFb7ZRBt3dz46d66wm2R7LBYwmchL\nUqQ8B3mllB3wOdADaAgMUkrVu6HMA0AtrXUdYDQwN6/7FSI7UzdPpXbZ2gxsNNC6Fe/fT1L9Frz8\nMrgGRNOnfDnr1l8MeXqCz7Ce/PvZOt6rWovZ4eEypdLa0tONXnwevpFaoyffBjiltQ7TWqcDS4F+\nN5TpBywC0FrvAjyVUlbsZgkBW8K28N3B75jXZ551561bLHD0KO+saEKz3kmYnM00c3OzXv3FWOs3\n7qNjyt8s/diZ53x9+V9ISGE3ybbkcagGrBPkfYGs39PCM9blVCYimzJC3LXYlFieXP4kXz/4tXWH\naQDCwkh19WTB72Vo+/wVHihbFrsScPFTbqjq1ShVrTxbPz1A94vV2B4Xx5arV2+/ocgdKwT5Ipn9\nesqUKZm/BwQEEJAlV4gQ2Xlu9XP0qtOLPvf0sXrdCdsPszu+CXOWwLzUKzxdsaLV91GcOfa6n9nR\n63hwYEve2lCTl4OD2dOyZYm4Cjjf3SLIBwYGEhgYmKsq8nwjb6VUO2CK1rpnxuPxgNZaz8hSZi6w\nSWv9U8bjE0AXrXVUNvXJjbzFHfnxyI9M3TyV/aP35z11wQ20hh8bvot3qQQ67HoPn+3bCW/fHk9b\nvDvI3Vq1Cj76iJeabeLEv5qoqft4o3o1HvH2LuyWFX9hYUaK4bCwHIvl94289wC1lVLVlVJOwEDg\njxvK/AE8mdGYdsDV7AK8EHcq9Goo49aM44eHf7B6gAeYOBF8og7T5fkmbLp6lZZubhLgb9SlC+zZ\nwweTE7GYFT5/1mTimTOY5AKpvCsKY/JaazMwFlgHHAOWaq2DlFKjlVKjMsqsBs4opYKBecCYW1Yo\nRC6lmlJ5bNljvN7xdVpWbmn1+ufOhZ9/hi5lD+PYojGroqPpVU5m1dzEzQ2aN8dx9z/88gtErixD\nUrgzH8qUyrwrKmPyWus1QN0b1s274fFYa+xLiGteXfcqvh6+vNz+ZavX/eWXMH06bPozEYf2Z0mo\nXZuf9u1jb0vrf5jYhIAACAzE/f77Wb9O0fXxekx75QC+Ts48WUnOYdy1hATjQzQP5IpXUSz9dPQn\nVgev5rt+31n1BJ/FAm+9BR99BCs2peF8diM0a8bCK1cI8PKihqur1fZlU3r2hF9+gXPnKF8etv3m\nQr3vmvDM/hC+PhNFYsbVsFprDickMOnMGT4JD8ci599yFhsLHh55qkIGF0WxczL6JGP/GsvaIWvx\ncvGyWr2RkfD000bnKXB9PPcGbSbcyYkuU6eyLzSUPxo1stq+bE6HDtCvHzRrBu3b4/nMM2xf1Jsn\nPmjE2Mhgxtb8l4oujmjAQSn6lS9P4NWrrLtyhcX161PG0THbatMsFjTgbFdC+6OxscZVZ3mQ59k1\n1iaza0ROYlNiaTe/HS+1e4lRLUdZpc5DsYm8+fdF1gcn0sXDi9k1gvh680ZC69Zl4QMP8JdS3OPq\nSnN3d6vsz6YlJcGyZfDVV8aMkLFj+bvWaJ4e70mz+1L5vxcs9KhfCqUU6RYL/zt9mt8uXWJ8tWq4\n2tmRrjXpWhOSnMzOuDgOJiRg0pqKTk40LF2axqVL83CFCrTNY++22PjmG9ixA+bPz7FYTrNrJMiL\nYsNsMfPgjw9Ss0xNPu/1eZ7ri4+Hqb/HMMvzOL5HKjI6wIXg4A1scnWlfLlyrOnUifJ5POlVoh0+\nDDNnwsqVpA0exhznl3h3UTXuvx9GjTJmBioFG2JiWBAZiR3gqBSOdnZUdXamvYcHrdzdcbWzIyw1\nlaOJiRxKSGDu+fN09fLieV9f2np42PZ8/JkzITwcZs3KsZgE+fxkMsGFCxAdbUyq9vMDL6885Zoo\nVkwmCA01enB16kA+jlm/uu5VDkQeYM0Ta3C0z/7r/e1ERsLff8NfazQr0iKwDDvDrCuOPOt2AWbM\nAG9vWLAAykoCMqsJD4dPP4X580nt8zDLyj7LglUVuJjsjn9vd+57wIF27cAnlxcqx5pMfH3+PF9f\nuECKxcJD5cvzUIUKdPL0xL6o/d9FRhpLw4ZwiyGpHP3vf+DuDpMm5VhMgnxepafDqVNw9KixnD5t\nfBUNCzP+gBUqGIvWRsADI9hXrmzkIK9Q4b+fXl5QurRxxvzGxdXVyDZXmAeq2WwE7ORk42dCAly+\n/N9y6RJERUFICAQHw9mzUKkSlCplvC+VK0ODBtC+PXTubNxJycUlz81aeHAh7259l10jd1HW9foA\nrLUx0ywh4b8lNta473Z4hOZ0pJmjoSaOnU2no+UPOlfYwupuDYnycue32bOp4+AAvr7QoweMGVNy\nPqAL2pUrMHs2/PorOj4ey9V4VGI8JhyI0+4k2rljLuWOpbQ7qT7VCR8wDtW2DV5eZC6lSxszCo1L\nFTRBSUn8dukSv12+THhqKl29vGjq5kY1Z2ec7Oxwyvhm4KQUTnZ2OCqV63V3lboiNRXWr4fVq2Ht\nWoiJMT69zp2D5s2hXTvjf6NFC+N/5VbfFFNSYN8+eOwx+P13aNUqx90WuyDvuHRnxiONk07FiRQc\ndTqJdu6kKufbVHDDY53NuozVKJ1tWTs0TjoFZ52Ck07FkTRMyhETjqQrR0w4YFYOmJU9Zuxv2oHC\ngj0m7LGgtBk7LNgpC3bagsKCQmOHRmHJKK8z1huN0MpYYzxD5k+t7Ljxr6VvcSDq69pkPLrxELj2\nUKtrbdCZb4PGLqMdZGmtHRZ17Xd7Ml5ZZhuM12U8Y69N2CszdpiN9w0n0nHI9m9xXWOyY2dC2afg\nmO6ACUcsKpuTcBnbqyw/lQKLvQVHbY8HDvjEXqT8pUjK1KnD/eXLM7JOHRxltkzh0hpSUtBx8USc\niCfkYDyXz8Rjf+QgnXfM4IhHR2ZXmMaxlFrExBj9jrQ0oy/i5PTf4ugIZu9k0uvFYqqWiKVMKjho\ncLCg7TU4aLSDJfOnvV06rvYJ2NunEefogcnBDu2gjbL2FrSjBjMosx3KrFAmO4x/1+wPVDttppQl\nEWedjFk5kKacSVPOmHAEBUpbcCQdB52Oo07DQZuww5zxP2UPWpP1v16hMWFPsl1pUuxuf5GfuX/H\n4hXkd02YgcvRvTidOIS5YhVS6zXD7FMZ9z+WkNKoBbHDX8JU1S/bsTjFzR0xRfZ/GqX+W28XGYHL\n3n9w3bYOlwM7SWvYnNRmbUlt2gZT/SZoF9fr6lVZ6rhxX9mtz3htN5XLWlYpjDl86SbsTOmQload\n2fipTOmQnoZKT888IMj42ymNsU4ZPzPXA8rBHuzt0MoO7OxRdnYoB3u0nR3YGetwdACXUlDKFRyd\n/guYqGxfh10u1ykFKvoyjmtW4vTDd2A2k/rsOMx9H77pq2vWv0XWdYejDjBtXl9++9kOuwf6Y791\nE9rNHfOQ4ZgeHwKlXHF0BCdHsj0enJTCMTYWBg82hpZ+/hlK2D1Zi62kJKPn//HHxt9v0iTj2zBG\nkE9PN5a0NGO5duhnXeC/3x3+PUapNb9SatMqHE//S3K7bqRVrYnH8oVcfOtzEno9nrlrrTVmjJPA\nJjTp2oLppi4WOESG4/Hrd3is/pG4XoOI7/8k6b7Vc/f6TGYcYi5hfzkKHBywOLuiXVyxOLtgcfME\n+9zNKNKAfwOXWwZ5tNZFagG0fu01rVeu1PrKFX2dhAStp0/Xunx5rUeN0vrsWX1XTCatd+3SeuZM\nrR95RGtfX629vbV+/HGtf/hB69jYu6tX3JrFovVff2kdEGD8/Z59VusdO3Lc5OTlk3rA6DI6payH\n1suXGyvNZq3//lvr/v2Nv9m772odHX3rSo4e1bp2ba1feEHrtDQrviBRYC5e1Pr557UuV07rd97R\nOi4ud9uZTMbff8oUrRs00LpqVa1feknrDRu0Tk39r9zevVrXqqX1M8/krm6TSevff9f6gQe0LlvW\nOLYiIu7utVmJEcpvEVNv9URhLRmNzVl0tNb/+5/WZcpo/cQTWh84cPttYmK0XrRI60GDjIOlYUOt\nx4zR+vvvtQ4ONoKQKBinT2s9bZrWNWtq3b698YF+w/sfFhOq3+lfVieW89B648bs6zl2TOsnn9Ta\n3V3rJk20fuwxrSdNMj6o9+3Tetky4wPlu+/y/zWJ/HfqlPH/6+2t9QcfaB0fr3VUlNFhW7rU6ACO\nGqX1ffcZQdvJSWs/P61ffFHr7duNDsKtxMZqPWKE8UHwxx/Zl0lK0nrWLKNMu3ZaL1igdWJi/rzW\nO5RTkC+SwzW5btPVq8Z83E8/hXvugWeegYce+u9En8UCf/0FCxcaJ0ECAqBPH+PqvKpV8+01iFwy\nm2H5cpg82TgJ9dlnUK8eESf3cfRRf5qllsVnzVbjJHZOEhONG22fOPHf8u+/xnf4hQuNk7/Cdhw7\nBlOmwG+/GUNv1atDjRrGUrPmf79Xr37nJ/03boTRo6FJE2O2Ve3axljP4sUwYYJxAnTSJChi6S2K\n3YnXO25TWprxB//2W9i/HwYOhFq1YM4c45LgUaPgkUdkWlxRlZ4OX3wB775LQt+epCxbSsjDXWk7\nb6VVZuYIG2U2g7299etNTjbOBcycCY8/bsysu3TJyFjXtq3192cFth/kswoNNXpvwcFGz75zZ5kS\nV0yEHtnKhhcexHnYUwwZ/nFhN0eUdBcvwocfGt8yx469u3nuBaRkBXlRLO2/sJ8Hf3yQSf6TeLbV\ns4XdHCGKlZyCvCQoE4Vuw+kNDPp1EHP7zOXh+g8XdnOEsClFMrVbujm9sJsgCsjcvXMZ/Ntglj26\nTAK8EPmgSPbky35QlnZV2uFfzR//6v608W2Dq6NcmWhLUk2pvPDXC2w7t41tI7ZRp1ydwm6SEDap\nSI7JRydF88/Zf9gStoUtZ7dw9OJRWlRqkRn0O1TtgLuzpH0trsKuhjHo10H4uPmwqP8i+VsKkUfF\n/sRrQloCO8N3GkE/bAt7z++lfoX6mUG/U7VOlCsl994sDpYeXcoLf73Aqx1e5dUOr2KXXR4aIcQd\nKfZB/kapplT2nN+TGfR3hO+gmmc1OlXtRIeqHehQtQM1y9S07TzTxUx0UjQvrX2JXRG7WPLwkny5\n8bYQJZXNBfkbmSwmDlw4wPZz29kevp1/zv5DmjktM+B3qNqBlpVayrh+IdBas/jwYl5b/xqPNXyM\n6d2nU9qpdGE3SwibYvNBPjvnYs+xI3yHEfjPbefYpWM08m5EhypG0G9ftT1VPKpYocXiVg5GHuSV\nda8QkxzDvD7zaO3burCbJIRNKpFB/kZJ6UnsO78vs7e//dx2XB1caVelHa0qt6J15da0qNQCT5e8\n3TRXQOjVUCZtmsT6kPVM9J/Is62excGuSE7kEsImSJDPhtaa4CvB7I7YzZ7ze9hzfg+HIg9RxaMK\nrX1b06pSK1r7tqZZxWaUcrx90n4BQZeC+Gj7Ryw/sZzn2zzPqx1elZkzQhQACfK5ZLKYOH7pOHsi\njKC/9/xejl86Tp1ydTJ7+k19mtLEp4kErwwWbWHjmY18tvszdobvZGzrsYxpPUZmOwlRgCTI50GK\nKYUjUUfYc34P+y/s51DUIY5fOk5Ft4qZAb+pT1OaVmyKn5dfiZkSGHo1lIUHF/Ldwe8o41qGUS1G\nMazZMPnWI0QhkCBvZSaLieArwRyKPMShKGM5HHWY2JRYGvs0plGFRtQrXy9zqeZZDXu7fEiJWoC0\n1hy9eJQVJ1aw/MRyzsWd47EGj/F0i6dpUalFYTdPiBJNgnwBiU6K5nDUYYIuB3Hi8onM5XLSZeqU\nq2ME/XL1uKfcPdQoUwM/Lz8qulUskr1/k8XEicsn2Bq2lcCwQAJDAyntWJq+dfvSv15/OlXrJCdT\nhSgiJMgXsoS0BE5GnyTokhH8T145SdjVMM5cPUNsSizVvarj5+WHn6cffl5+VHavTEW3ivi4+VDR\nrSIVSlXIt28CSelJnIk5w5mrZwi5EsLRi0c5GHWQ45eO4+vuS8dqHQmoHkAXvy74efnlSxuEEHkj\nQb4IS0pPygz4oVdDCb0aSmRCZOYSlRjFleQrlHMtR4XSFfB09sTD2eO6pbRjaRzsHLC3s8de2Wf+\nTLekk2JKIcWUQnJ6MonpiUQnR3M56TKXEi9xKekS8anx+Hn5UaNMDWp41aCRdyOaV2xOY5/GuDm5\nFfbbI4TIBQnyxZzJYsoMynGpcTctCWkJmCwmzBYzZm3O/Olk74SLg0vmUsqxFOVLlb9u8S7tXSSH\ni4QQuSdBXgghbFhOQV66cEIIYcMkyAshhA2TIC+EEDYsT0FeKVVGKbVOKfWvUmqtUuqm7F5KqSpK\nqY1KqWNKqSNKqRfysk8hhBC5l9ee/Hjgb611XWAjMCGbMibgZa11Q6A98JxSql4e9yuEECIX8hrk\n+wELM35fCPS/sYDWOlJrfTDj9wQgCPDN436FEELkQl6DvLfWOgqMYA5451RYKeUHNAN25XG/Qggh\ncuG2yUeUUusBn6yrAA1MzKb4LSe4K6XcgF+AcRk9eiGEEPnstkFea33frZ5TSkUppXy01lFKqYrA\nxVuUc8AI8N9rrX+/3T6nTJmS+XtAQAABAQG320QIIUqMwMBAAgMDc1U2T1e8KqVmAFe01jOUUq8D\nZbTW47Mptwi4rLV+ORd1yhWvQghxB/ItrYFSqizwM1AVCAMe01pfVUpVAr7WWvdRSnUEtgBHMIZz\nNPCG1nrNLeqUIC+EEHdActcIIYQNk9w1QghRQkmQF0IIGyZBXgghbJgEeSGEsGES5IUQwoZJkBdC\nCBsmQV4IIWyYBHkhhLBhEuSFEMKGSZAXQggbJkFeCCFsmAR5IYSwYRLkhRDChkmQF0IIGyZBXggh\nbJgEeSGEsGES5IUQwoZJkBdCCBsmQV4IIWyYBHkhhLBhEuSFEMKGSZAXQggbJkFeCCFsmAR5IYSw\nYRLkhRDChkmQF0IIGyZBXgghbJgEeSGEsGES5IUQwoZJkBdCCBsmQV4IIWyYBHkhhLBhEuSFEMKG\nSbf9QLcAAARuSURBVJAXQggbJkFeCCFsmAR5IYSwYXkK8kqpMkqpdUqpf5VSa5VSnjmUtVNK7VdK\n/ZGXfQohhMi9vPbkxwN/a63rAhuBCTmUHQccz+P+hBBC3IG8Bvl+wMKM3xcC/bMrpJSqAvQCvsnj\n/oQQQtyBvAZ5b611FIDWOhLwvkW5WcBrgM7j/oQQQtwBh9sVUEqtB3yyrsII1hOzKX5TEFdK9Qai\ntNYHlVIBGdvnaMqUKZm/BwQEEBAQcLtNhBCixAgMDCQwMDBXZZXWd9+5VkoFAQFa6yilVEVgk9a6\n/g1lpgFDABPgCrgDv2mtn7xFnTovbRJCiJJGKYXWOtsOdF6Ha/4Ahmf8Pgz4/cYCWus3tNbVtNY1\ngYHAxlsFeCGEENaV1yA/A7hPKfUv0B14H0ApVUkp9WdeGyeEECJv8jRckx9kuEYIIe5Mfg7XCCGE\nKMIkyAshhA2TIC+EEDZMgrwQQtgwCfJCCGHDJMgLIYQNkyAvhBA2TIK8EELYMAnyQghhwyTICyGE\nDZMgL4QQ/9/e/bxYVcdhHH8/5kawFm1cOIhLYVoMLQZCqFDEYYJpPQQibosERaLsLwjEheHKlYK4\nsEW5CJpoXUE6MkxTuhCx6AdFm3b+eFrcO9yBZq6zkPM5fu/zWt1z7l08fM7l4Xu/nMtpWEo+IqJh\nKfmIiIal5CMiGpaS38R2H6s1CTKLkcxiJLMY6fssUvKb6PtF61JmMZJZjGQWI32fRUo+IqJhKfmI\niIb18vF/1RkiIp43Wz3+r3clHxERz062ayIiGpaSj4hoWEr+KSSdlvRE0svVWapI+kTSmqRlSZ9J\neqk6U5ckzUn6SdIdSR9U56kiaUrSN5JWJa1Ier86UzVJOyTdlPRFdZatpOTHkDQFHAHuV2cp9hUw\nbXsGuAt8WJynM5J2AJ8CR4FpYFHSgdpUZR4Bp2xPA68B707wLNadBH6sDjFOSn6888CZ6hDVbH9t\n+8nw8FtgqjJPx2aBu7bv234IXAPeLs5UwvbvtpeHr/8F1oC9tanqDBeB88Cl6izjpOS3IGkBeGB7\npTpLz5wAvqwO0aG9wIMNx78wwcW2TtJ+YAb4rjZJqfVFYK9vUdxZHaCSpCVgz8ZTDC7Yx8BHDLZq\nNr7XrDGzOGv7xvAzZ4GHtq8WRIyekLQbuA6cHK7oJ46kt4A/bC9LepMe98NEl7ztI5udl/QKsB+4\nLUkMtid+kDRr+88OI3Zmq1msk3ScwU/TQ50E6o9fgX0bjqeG5yaSpJ0MCv6K7c+r8xQ6CCxImgd2\nAS9Kumz7WHGu/8mfobZB0j3gVdv/VGepIGkOOAe8bvvv6jxdkvQC8DNwGPgN+B5YtL1WGqyIpMvA\nX7ZPVWfpC0lvAKdtL1Rn2Uz25LfH9PjnWAcuALuBpeHtYherA3XF9mPgPQZ3GK0C1ya44A8C7wCH\nJN0afhfmqnPFeFnJR0Q0LCv5iIiGpeQjIhqWko+IaFhKPiKiYSn5iIiGpeQjIhqWko+IaFhKPiKi\nYf8B6548zP98pXAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbbaf40d358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(reals, -p0.nnlf(reals.reshape(-1, 1))  \n",
    "                +p1.nnlf(reals.reshape(-1, 1)), label=\"True ratio\")\n",
    "\n",
    "plt.plot(reals, cc_none.predict(reals.reshape(-1, 1), log=True), label=\"No calibration\")\n",
    "plt.plot(reals, cc_direct.predict(reals.reshape(-1, 1), log=True), label=\"Calibration\")\n",
    "plt.plot(reals, cc_decomposed.predict(reals.reshape(-1, 1), log=True), label=\"Calibration + Decomposition\")\n",
    "\n",
    "plt.xlim(-5, 5)\n",
    "plt.ylim(-0.5, 0.5)\n",
    "plt.legend(loc=\"best\", prop={'size': 8})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using density ratios for maximum likelihood fit\n",
    "\n",
    "In the likelihood-free setting, density ratios can be used to find the maximum likelihood estimator $\\theta^* = \\arg \\max_{\\theta} p(D | \\theta)$ by noticing that $\\theta^*$ also maximizes $\\prod_{x \\in D} \\frac{p(x|\\theta)}{p(x|\\theta_1)}$ for some fixed value of $\\theta_1$.\n",
    "\n",
    "As an example, this can be used to find the bump coefficient in $p_1$, as illustrated below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.utils import check_random_state\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "n_trials = 1000  \n",
    "mles = []\n",
    "\n",
    "for i in range(n_trials):  \n",
    "    # Reset\n",
    "    rng = check_random_state(i)\n",
    "    p0.set_params(random_state=rng)\n",
    "    p1.set_params(random_state=rng)\n",
    "    \n",
    "    g.set_value(bump_coefficient)\n",
    "    X_true = p0.rvs(5000)\n",
    "    \n",
    "    # Fit ratio\n",
    "    clf = MLPRegressor(activation=\"tanh\", hidden_layer_sizes=(10, 10), random_state=rng)\n",
    "    cc = DecomposedRatio(ClassifierRatio(base_estimator=CalibratedClassifierCV(clf, cv=cv)))\n",
    "    cc.fit(numerator=p0, denominator=p1, n_samples=100000)\n",
    "\n",
    "    # Inference\n",
    "    def objective(theta):       \n",
    "        g.set_value(theta[0])  # this indirectly updates weights in p1, \n",
    "                               # along with the density ratios computed by cc \n",
    "                               # (without having to retrain the classifiers since \n",
    "                               # g only affects the weights!)\n",
    "\n",
    "        return -np.mean(cc.predict(X_true, log=True))\n",
    "\n",
    "    results = minimize(objective, x0=[0.1], \n",
    "                       constraints=[{'type':'ineq', 'fun': lambda x: x[0]},\n",
    "                                    {'type':'ineq', 'fun': lambda x: 1. - x[0]},])\n",
    "\n",
    "    mles.append(results.x[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.049954620813502319"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(mles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFoVJREFUeJzt3X+MXeV95/H319gQQvGPtGBIbAwJS2gqtJSW4CRb5VIo\ndaiCk0plSasFQ8QikS6oVbIYKsUzClEAZZM0apAmE2BNN23Kj21ws+VHkPeStJJLTCBxw5h1ArYB\nl8ERnUTBSgz4u3/MsXM9nvFczz137j1z3i9p5HvOfe49Xz8e3889zznnOZGZSJLqa16vC5Ak9ZZB\nIEk1ZxBIUs0ZBJJUcwaBJNWcQSBJNTdtEETEHRExGhHfb1l3W0SMRMRTEXF/RCxsee7GiNhWPH9R\ntwqXJJWjnT2Cu4Dfn7DuEeA3MvNsYBtwI0BEvAu4FPh14APA7RER5ZUrSSrbtEGQmf8E/PuEdY9m\n5r5icROwrHh8CfC1zHw9M7czHhLvLq9cSVLZyjhGcBXwj8XjtwHPtzz3YrFOktSnOgqCiPgL4LXM\n/NuS6pEkzbL5M31hRKwBLgZ+t2X1i8DyluVlxbrJXu8kR5I0A5lZ6rHXdvcIovgZX4hYBXwCuCQz\nf9HSbgNwWUQcHRGnAacDj0/1pplZ2Z9169b1vAbr730dday/yrXPhfq7Ydo9goj4G6AB/GpE7ATW\nATcBRwPfLE4K2pSZ12bm0xFxD/A08BpwbXarcklSKaYNgsz840lW33WY9p8BPtNJUZKk2eOVxTPU\naDR6XUJHrL+3qlx/lWuH6tffDdGrkZuIcNRIko5QRJAlHyye8VlDkjSdU089lR07dvS6jEpasWIF\n27dvn5VtuUcgqWuKb6+9LqOSpuq7buwReIxAkmrOIJCkmjMIJKnmDAJJqjmDQFItnXbaaWzcuLHX\nZfQFg0CVNDAw0OsSNMe98cYbvS5h1hgEqqTBwcFel6AKu/zyy9m5cycf/OAHWbhwIZ/97GfZsWMH\n8+bN484772TFihVccMEFPPbYYyxfvvyg17buSWQmt9xyC6effjonnHACl112GWNjY1Nu97bbbuOt\nb30ry5Yt44477mDevHk8++yzXf27tsMgkFQ7d999N6eccgrf+MY3+OlPf8rHP/7xA89961vfYuvW\nrTz88MPA+Hn7U/niF7/Ihg0b+Pa3v82uXbtYsmQJ11577aRtH3roIb7whS+wceNGfvjDH9JsNg/7\n3rPJIJDUMwMDA0TEIT9TDf1N1r6TYcKJF2xFBIODgxx77LEcc8wx075+aGiIT3/605x88sksWLCA\nT37yk9x3333s27fvkLb33nsvV155JWeeeSZvetOb+mp40ykmJPXMwMDAEX0gHmn7mVi2bNn0jQo7\nduzgwx/+MPPmjX+nzkwWLFjA6OgoJ5988kFtd+3axbnnnntgefny5X1z1bVBIKmWphqWaV1/3HHH\nsWfPngPLb7zxBrt37z6wfMopp3DnnXfynve8Z9rtnXzyybzwwgsHlnfu3OnQkNSJdevW9boEVdxJ\nJ510yIHaid/QzzjjDH7+85/z4IMP8vrrr3PzzTezd+/eA89fc8013HTTTezcuROA3bt3s2HDhkm3\nd+mll3LXXXexdetW9uzZw80331zy32jmDAJVUj+Nr6qa1q5dy6c+9Sne8pa38LnPfQ44dC9h4cKF\n3H777Xz0ox9l2bJlHH/88QcNHV1//fWsXr2aiy66iEWLFvHe976Xxx+f/O68q1at4rrrruP888/n\njDPOOLAX0c6xiG5z9lFJXePso1PbunUrZ511Fr/4xS8OHGNoNZuzjxoEmhM2bdrC2NjeadstXnw0\nK1eeNQsVCQyCib7+9a9z8cUX8+qrr7JmzRrmz5/P/fffP2nb2QwCDxZrThgb28sJJ/zWtO12735i\nFqqRJjc0NHQgABqNBl/60pd6XRJgEEjSrHnwwQd7XcKkPFisSvJgsVQeg0CV5FxDUnkMAkmqOY8R\nSOqaFStW9M3Vs1WzYsWKWduWQSCpa7Zv397rEtQGh4YkqeYMAlWScw1J5Zk2CCLijogYjYjvt6xb\nEhGPRMQzEfFwRCxqee7GiNgWESMRcVG3Cle9efqoVJ529gjuAn5/wrq1wKOZ+U5gI3AjQES8C7gU\n+HXgA8Dt4ZEiSepr0wZBZv4T8O8TVq8G1heP1wMfKh5fAnwtM1/PzO3ANuDd5ZQqSeqGmR4jODEz\nRwEy8yXgxGL924DnW9q9WKyTJPWpsg4WO72gJFXUTK8jGI2IpZk5GhEnAS8X618Elre0W1asm1Tr\nAb9Go0Gj0ZhhOaqb2bh3rdQPms0mzWazq9to634EEXEq8A+ZeVaxfCvwSmbeGhE3AEsyc21xsPir\nwHmMDwl9E/gPk914wPsRqBMT52p/6KEn2p6GetWq6dtJ/aon9yOIiL8BGsCvRsROYB1wC3BvRFwF\n7GD8TCEy8+mIuAd4GngNuNZPe0nqb9MGQWb+8RRPXThF+88An+mkKEnS7PHKYkmqOSedU1873L2I\nH3rol7edHBl5rq1jBJIOZRCor011L+Krr1530PrNm7fNZlnSnOLQkCrpmmsGel2CNGcYBJJUcwaB\nJNWcQSBJNWcQSFLNGQSqpKGhgV6XIM0ZBoEqaXh4sNclSHOGQSBJNWcQSFLNGQSSVHMGgSTVnEGg\nSrr66nW9LkGaMwwCVZJzDUnlMQgkqeYMAkmqOYNAkmrOIJCkmjMIVEnONSSVxyBQJTnXkFQeg0CS\nas4gkKSaMwgkqeYMAkmqOYNAleRcQ1J5DAJVknMNSeUxCCSp5joKgoj4s4j414j4fkR8NSKOjogl\nEfFIRDwTEQ9HxKKyipUklW/+TF8YEW8F/htwZmbujYi/Az4CvAt4NDNvi4gbgBuBtaVUq0rYtGkL\nY2N7p223ePHRrFx51ixUJOlwZhwEhaOA4yJiH3As8CLjH/zvL55fDzQxCGplbGwvJ5zwW9O22737\niVmoRtJ0ZhwEmbkrIv4HsBPYAzySmY9GxNLMHC3avBQRJ5ZUq3TA0NDAjA4Yj4z8qK127q2oTjoZ\nGloMrAZWAD8B7o2IPwFyQtOJywcMDAwceNxoNGg0GjMtRzUzPDw4oyB49dV97q2oUprNJs1ms6vb\n6GRo6ELg2cx8BSAi/h54LzC6f68gIk4CXp7qDVqDQJJ0qIlfkgcHy59wsZOzhnYCKyPiTRERwAXA\n08AGYE3R5grggY4qlCR1VSfHCB6PiPuAJ4HXij+/DBwP3BMRVwE7gEvLKFSS1B0dnTWUmYPAxP2U\nVxgfNpIkVYBXFquSnGtIKo9BoEpyriGpPAaBJNWcQSBJNWcQSFLNGQSSVHMGgSppaGig1yVIc4ZB\noEoaHi7/MnuprgwCSao5g0CSas4gkKSa6/QOZdKMtXOTmJGR59q6f4CkmTMI1DPt3CRm8+Ztk67v\n9lxD7YSUdzHTXGEQqJK6PddQOyHlXcw0V3iMQJJqziCQpJozCCSp5gwCSao5DxarbZs2bWFsbO+0\n7WbjlM+hoQFvTiOVxCBQ28bG9rb1AT/VKZ9lGh4eNAikkjg0JEk1ZxBIUs0ZBJJUcwaBJNWcQaBK\n6vZcQ1KdGASqJM8YkspjEEhSzRkEklRzHQVBRCyKiHsjYiQifhAR50XEkoh4JCKeiYiHI2JRWcVK\nksrX6ZXFfwn8Y2b+UUTMB44DbgIezczbIuIG4EZgbYfbkfpOOzevAW9go/434yCIiIXA72TmGoDM\nfB34SUSsBt5fNFsPNDEIVLJ+mGuonZvXgDewUf/rZGjoNODHEXFXRHw3Ir4cEW8GlmbmKEBmvgSc\nWEahUqvh4cFelyDNGZ0MDc0HzgE+lpmbI+LzjH/zzwntJi4fMDAwcOBxo9Gg0Wh0UI4kzT3NZpNm\ns9nVbXQSBC8Az2fm5mL5fsaDYDQilmbmaEScBLw81Ru0BoEk6VATvyQPDpa/NzzjoaFi+Of5iDij\nWHUB8ANgA7CmWHcF8EAnBUqSuqvTs4auA74aEQuAZ4ErgaOAeyLiKmAHcGmH25AkdVFHQZCZ3wPO\nneSpCzt5X2k6zjUklccri1VJvT51VJpLDAJJqjmDQJJqziCQpJozCCSp5gwCVdLQ0ECvS5DmjE6v\nI9AcsGnTFsbG9k7bbmTkubYmWZsNw8ODnjkklcQgEGNje9v6gN+8edssVCNpthkEc1w73/b76Zu+\npNlnEMxx7Xzb95u+VG8eLJakmjMIVEnONSSVxyBQJXnGkFQeg0CSas4gkKSaMwgkqeYMAkmqOYNA\nleRcQ1J5DAJV0vDwYK9LkOYMg0CSas4gkKSaMwgkqeYMAkmqOYNAleRcQ1J5DAJVknMNSeUxCCSp\n5gwCSao571AmddnIyI/aard48dGsXHlWl6uRDtVxEETEPGAz8EJmXhIRS4C/A1YA24FLM/MnnW5H\nqqpXX93X1j2hd+9+YhaqkQ5VxtDQ9cDTLctrgUcz853ARuDGErYhHcS5hqTydBQEEbEMuBj4Ssvq\n1cD64vF64EOdbEOajHMNSeXpdI/g88AngGxZtzQzRwEy8yXgxA63IUnqohkfI4iIPwBGM/OpiGgc\npmlO9cTAwMCBx41Gg0bjcG8jSfXTbDZpNptd3UYnB4vfB1wSERcDxwLHR8RfAy9FxNLMHI2Ik4CX\np3qD1iCQJB1q4pfkwcHyh0VnPDSUmTdl5imZ+XbgMmBjZv4X4B+ANUWzK4AHOq5SktQ13big7Bbg\n9yLiGeCCYlkqlXMNSeUp5YKyzHwMeKx4/ApwYRnvK03FuYak8jjFhCTVnEEgSTXnXENSn3BOIvWK\nQSD1CeckUq84NKRKcq4hqTwGgSrJuYak8jg01Gc2bdrC2Njeads5TiypLAZBnxkb2+s4saRZ5dCQ\nJNWcQSBJNWcQqJKca0gqj0GgSnKuIak8BoEk1ZxBIEk1ZxBIUs0ZBJJUcwaBKsm5hqTyGASqJOca\nkspjEEhSzRkEklRzBoEk1Zyzj1ZUu7c1HBl5rq3ZTCXVl0FQUe3e1nDz5m2zUM3sc64hqTwODamS\nnGtIKo9BIEk1ZxBIUs0ZBJJUcwaBJNXcjIMgIpZFxMaI+EFEbImI64r1SyLikYh4JiIejohF5ZUr\njXOuIak8newRvA78eWb+BvAe4GMRcSawFng0M98JbARu7LxM6WDONSSVZ8ZBkJkvZeZTxeOfASPA\nMmA1sL5oth74UKdFSpK6p5RjBBFxKnA2sAlYmpmjMB4WwIllbEOS1B0dB0FE/ApwH3B9sWeQE5pM\nXJYk9ZGOppiIiPmMh8BfZ+YDxerRiFiamaMRcRLw8lSvHxgYOPC40WjQaDQ6KUeS5pxms0mz2ezq\nNjqda+hO4OnM/MuWdRuANcCtwBXAA5O8Djg4CKQj4VxDqouJX5IHB8s/UWLGQRAR7wP+BNgSEU8y\nPgR0E+MBcE9EXAXsAC4to1CplXMNSeWZcRBk5j8DR03x9IUzfV9J0uzyymJJqjmDQJJqzhvTSHPQ\npk1bGBvbO227xYuPZuXKs2ahIvUzg0CVNDQ04AHjwxgb29vWHex2735iFqpRvzMIZkm739C8x3B7\nhocHaxsE7dyv2t8jHQmDYJa0+w1trt5jWOVp537V/h7pSHiwWJJqziCQpJozCCSp5gwCVZJzDUnl\nMQhUSXU9Y0jqBs8aKkE7p4Z6Op+qzAvU5jaDoATtnBrq6XyqMi9Qm9scGpKkmjMIJKnmDAJV0tDQ\nQK9LkOYMg0CVNDxc/u36pLoyCCSp5jxrSKqxdmYyHW/n6c9zmUEg1Vg7M5mCpz/PdQ4NSVLNGQSq\nJOcaksrj0NBheFex/uVcQ1J5DILD8K5ikurAoSFJqjmDQJJqziCQpJrzGIEqaWhowAPGfaidC9S8\nZ0H/6VoQRMQq4AuM73XckZm3dmtbqp/h4UGDoA+1c4Ga9yzoP10JgoiYB/wVcAGwC/hORDyQmVu7\nsb39ZvMuSps3N/nt32509B69VPX6q67K/V/l2gGazSaNRqPXZfSVbu0RvBvYlpk7ACLia8BqYEZB\nMDY2xs9+9uq07Xbu/DHveMf507Yr4xvJE09U+z9D1euvuir3f5VrB4NgMt0KgrcBz7csv8B4OMzI\nyMguXnnleObPXzBlm3373uC553bzjne0835OtCX1Srv//55/fifLl58ybbsyjznU9d7MlThYvGBB\nELGHzKMO0yqZPz/aej8n2pJ650j+/51zzuzeJ7mu92aOzCz/TSNWAgOZuapYXgtk6wHjiCh/w5JU\nA5nZ3rfeNnUrCI4CnmH8YPG/AY8DH8nMkdI3JknqSFeGhjLzjYj4U+ARfnn6qCEgSX2oK3sEkqTq\nKG2KiYhYFRFbI+L/RcQNU7T5YkRsi4inIuLsYt0xEfEvEfFkRGyJiHUt7ddFxAsR8d3iZ1VZ9ZZR\ne8tz84r6NrSsWxIRj0TEMxHxcEQs6kbtXax/Vvp+hvX/Zsv67RHxveL35/GW9f3c/+3U38/9f3bL\n+kURcW9EjETEDyLivGJ9P/d/O/X39WdPRJxR/M58t/jzJxFxXfHckfd9Znb8w3ig/BBYASwAngLO\nnNDmA8D/KR6fB2xqee7NxZ9HAZuAdxfL64A/L6PGbtVerPsz4H8BG1rW3Qr89+LxDcAtFau/631f\n0u/Os8CSSd63Ev1/mPqr0v//E7iyeDwfWFix/p+q/kp89rS8zy5g2Uz7vqw9ggMXkGXma8D+C8ha\nrQbuBsjMfwEWRcTSYnlP0eaY4h+jdbyq1KPjk+io9ohYBlwMfGWS16wvHq8HPtSd8rtWP3S/76HD\n+osaJ/s9rkT/M3X9+5/rthnXHxELgd/JzLuK517PzJ+2vKav+3+a+qHPP3taXAj8KDNfaHnNEfV9\nWUEw2QVkb5umzYv72xRDE08CLwHfzMzvtLT702KX6Ctd2r3sqHbg88AnODi8AE7MzFGAzHwJOLGs\ngqepraz6oft9P1ltR1p/At+MiO9ExNUtbarS/1PVD/3f/6cBP46Iu4ohii9HxLFFmyr0/+Hqh/7/\n7NnvPwN/27J8xH3fF9NQZ+a+zPxNYBlwXkS8q3jqduDtmXk24yHxuV7VOJmI+ANgNDOfYvzbw+G+\nQfTdUflp6u/rvm/xvsw8h/G9mo9FxH+aol3f9X9hqvqr0P/zgXOALxV/hz3A2uK5if8X+rH/D1d/\nFfqfiFgAXALce5hm0/Z9WUHwItB6LfiyYt3ENssP16bYLfu/wKpieXcWA13AMHBuSfVOrGumtb8P\nuCQinmU8kc+PiLuLNqMtwy8nAS93ofb9tZVe/yz1/f7aZvy7k5n/Vvy5G/h7fjmVSRX6f8r6K9L/\nLwDPZ+bmYv19jH+wArxUgf6fsv4KfPbs9wHgieL3Z78j/t0vKwi+A5weESsi4mjgMmDDhDYbgMuL\n4lYCY5k5GhG/tn+3q9gt+z2KyemKv8R+fwj8a0n1llJ7Zt6Umadk5tuL123MzMtbXrOmeHwF8EAX\nau9a/bPU9x3VHxFvjohfKdYfB1zUUmff9//h6q9C/xfDD89HxBlFuwuAp1tes6Z43Jf9f7j6+/2z\np+X5j3DwsND+16wpHrfX9yUeAV/F+NXE24C1xbprgP/a0uavGD9K/j3gnGLdWcB3GT9i/n3gL1ra\n312sewr4OrC0rHrLqH3Ce7yfg8+6eQvwaPG+jwCLu1F7F+uflb7v8HfntKK+J4Et+19blf6fpv6+\n7/9i/X9k/APtKeB/A4uq0v/T1N/3nz3Am4HdwPET3vOI+94LyiSp5vriYLEkqXcMAkmqOYNAkmrO\nIJCkmjMIJKnmDAJJqjmDQJJqziCQpJr7/xS2H6hJ3qlmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbbaf08cc88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "h = plt.hist(mles, bins=30, normed=1, alpha=0.2)\n",
    "plt.vlines(bump_coefficient, 0, h[0].max()+5, linestyles=\"dashed\", label=\"true g\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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
}