A randomized selective bootstrap

randomized_selective_bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%load_ext rpy2.ipython\n",
    "%pylab inline\n",
    "import numpy as np\n",
    "from scipy.stats import t as tdist, norm as ndist\n",
    "import scipy.stats\n",
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 50\n",
    "kappa = 0.8\n",
    "sd = 1.5\n",
    "truth = 1. / np.sqrt(n) * sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def logistic(n=1):\n",
    "    U = np.random.sample(n)\n",
    "    return np.log(U / (1 - U))\n",
    "\n",
    "def Gbar(z):\n",
    "    return np.exp(-z) / (1 + np.exp(-z))\n",
    "\n",
    "def selection_event(S, threshold):\n",
    "    return scipy.stats.ttest_1samp(S, 0)[0] + kappa * logistic() > threshold\n",
    "\n",
    "def bootstrap(sample, nsample=500):\n",
    "    boot = []\n",
    "    n = sample.shape[0]\n",
    "    for _ in range(nsample):\n",
    "        bootS = np.random.choice(sample, n, replace=True)\n",
    "        boot.append(bootS.mean())\n",
    "    return np.array(boot)\n",
    "\n",
    "\n",
    "def interval(dbn_func, observed, sd, n, Z, coverage=0.9):\n",
    "\n",
    "    f1 = lambda theta: dbn_func((observed - theta) / sd, np.sqrt(n) * theta, Z) - (1 - coverage) / 2\n",
    "    v1 = scipy.optimize.bisect(f1, observed - 6 * sd, observed + 6 * sd )\n",
    "    f2 = lambda theta: dbn_func((observed - theta) / sd, np.sqrt(n) * theta, Z) - (1 - (1 - coverage) / 2)\n",
    "    v2 = scipy.optimize.bisect(f2, observed - 6 * sd, observed + 6 * sd )\n",
    "    return v2, v1\n",
    "\n",
    "\n",
    "def noise(n):\n",
    "    # this noise has exponential moments...\n",
    "    Z = (np.random.standard_normal(n) + 0.5 * np.random.standard_normal(n) * np.random.binomial(1, 0.5, size=(n,))) \n",
    "    return Z\n",
    "\n",
    "noise_sd = np.std(noise(50000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.713832656894 True 1.2131848998 1\n"
     ]
    }
   ],
   "source": [
    "def simulation(n, threshold=1.65, nsample=2000, coverage=0.90):\n",
    "    count = 0\n",
    "    while True:\n",
    "        count += 1\n",
    "        S = noise(n) / noise_sd * sd + truth\n",
    "        if selection_event(S, threshold):\n",
    "            break\n",
    "    \n",
    "    boot_sample = bootstrap(S, nsample)\n",
    "    bootZ = boot_sample - S.mean()\n",
    "    bootZ /= boot_sample.std()\n",
    "    \n",
    "    gaussZ = np.random.standard_normal(2000)\n",
    "    \n",
    "    def pivot(t, delta, Z):\n",
    "        num = (Gbar((threshold - Z - delta) / kappa) * (Z <= t)).mean()\n",
    "        den = (Gbar((threshold - Z - delta) / kappa)).mean()\n",
    "        return num / den\n",
    "    \n",
    "    bootL, bootU = interval(pivot, S.mean(), boot_sample.std(), n, bootZ, coverage=coverage)\n",
    "    \n",
    "    gaussL, gaussU = interval(pivot, S.mean(), boot_sample.std(), n, gaussZ, coverage=coverage)\n",
    "    boot_cover = (bootL < truth) * (bootU > truth)\n",
    "    gauss_cover = (gaussL < truth) * (gaussU > truth)\n",
    "    return pivot((S.mean() - truth) / boot_sample.std(), np.sqrt(n) * truth, bootZ), count, pivot, boot_cover, (bootU - bootL) / (2 * np.fabs(ndist.ppf((1 - coverage) / 2.)) * np.std(S) / np.sqrt(n))\n",
    "    \n",
    "p, count, pivot, covered, length = simulation(n)\n",
    "print p, covered, length, count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probability of selection: 0.355113636364\n",
      "coverage: 0.893\n",
      "number of root finding fails: 0\n",
      "average ratio of length to T interval: 1.21957241388\n"
     ]
    }
   ],
   "source": [
    "P = []\n",
    "C = []\n",
    "coverage = []\n",
    "num_except = 0\n",
    "length = []\n",
    "for i in range(1000):\n",
    "    try:\n",
    "        result = simulation(n)\n",
    "    except ValueError: # probably root finding failed\n",
    "        num_except += 1\n",
    "    P.append(result[0])\n",
    "    C.append(result[1])\n",
    "    coverage.append(result[-2])\n",
    "    length.append(result[-1])\n",
    "print 'probability of selection:', 1. / np.mean(C)\n",
    "print 'coverage:', np.mean(coverage)\n",
    "print 'number of root finding fails:', num_except\n",
    "print 'average ratio of length to T interval:', np.mean(length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAYAAAB91L6VAAAD8GlDQ1BJQ0MgUHJvZmlsZQAAOI2N\nVd1v21QUP4lvXKQWP6Cxjg4Vi69VU1u5GxqtxgZJk6XpQhq5zdgqpMl1bhpT1za2021Vn/YCbwz4\nA4CyBx6QeEIaDMT2su0BtElTQRXVJKQ9dNpAaJP2gqpwrq9Tu13GuJGvfznndz7v0TVAx1ea45hJ\nGWDe8l01n5GPn5iWO1YhCc9BJ/RAp6Z7TrpcLgIuxoVH1sNfIcHeNwfa6/9zdVappwMknkJsVz19\nHvFpgJSpO64PIN5G+fAp30Hc8TziHS4miFhheJbjLMMzHB8POFPqKGKWi6TXtSriJcT9MzH5bAzz\nHIK1I08t6hq6zHpRdu2aYdJYuk9Q/881bzZa8Xrx6fLmJo/iu4/VXnfH1BB/rmu5ScQvI77m+Bkm\nfxXxvcZcJY14L0DymZp7pML5yTcW61PvIN6JuGr4halQvmjNlCa4bXJ5zj6qhpxrujeKPYMXEd+q\n00KR5yNAlWZzrF+Ie+uNsdC/MO4tTOZafhbroyXuR3Df08bLiHsQf+ja6gTPWVimZl7l/oUrjl8O\ncxDWLbNU5D6JRL2gxkDu16fGuC054OMhclsyXTOOFEL+kmMGs4i5kfNuQ62EnBuam8tzP+Q+tSqh\nz9SuqpZlvR1EfBiOJTSgYMMM7jpYsAEyqJCHDL4dcFFTAwNMlFDUUpQYiadhDmXteeWAw3HEmA2s\n15k1RmnP4RHuhBybdBOF7MfnICmSQ2SYjIBM3iRvkcMki9IRcnDTthyLz2Ld2fTzPjTQK+Mdg8y5\nnkZfFO+se9LQr3/09xZr+5GcaSufeAfAww60mAPx+q8u/bAr8rFCLrx7s+vqEkw8qb+p26n11Aru\nq6m1iJH6PbWGv1VIY25mkNE8PkaQhxfLIF7DZXx80HD/A3l2jLclYs061xNpWCfoB6WHJTjbH0mV\n35Q/lRXlC+W8cndbl9t2SfhU+Fb4UfhO+F74GWThknBZ+Em4InwjXIyd1ePnY/Psg3pb1TJNu15T\nMKWMtFt6ScpKL0ivSMXIn9QtDUlj0h7U7N48t3i8eC0GnMC91dX2sTivgloDTgUVeEGHLTizbf5D\na9JLhkhh29QOs1luMcScmBXTIIt7xRFxSBxnuJWfuAd1I7jntkyd/pgKaIwVr3MgmDo2q8x6IdB5\nQH162mcX7ajtnHGN2bov71OU1+U0fqqoXLD0wX5ZM005UHmySz3qLtDqILDvIL+iH6jB9y2x83ok\n898GOPQX3lk3Itl0A+BrD6D7tUjWh3fis58BXDigN9yF8M5PJH4B8Gr79/F/XRm8m241mw/wvur4\nBGDj42bzn+Vmc+NL9L8GcMn8F1kAcXgSteGGAABAAElEQVR4Ae2dB5wURdbAHyo5SlAyEhSVoKic\nwnkocigeJg4jqCiiiB535CiggGREQDwR/Tj1BBXlFAETCuopHnKAoARBonigSM5pvvcad5jdSb3s\n7EzPzL9+v97prqqu8O/eeVNVr97L49MgBAhAAAIQgAAE4krgtLjWRmUQgAAEIAABCDgEEMC8CBCA\nAAQgAIEEEEAAJwA6VUIAAhCAAAQQwLwDEIAABCAAgQQQQAAnADpVQgACEIAABBDAvAMQgAAEIACB\nBBBAACcAOlVCAAIQgAAEEMC8AxCAAAQgAIEEEEAAJwA6VUIAAhCAAAQQwLwDEIAABCAAgQQQQAAn\nADpVQgACEIAABBDAvAMQgAAEIACBBBBAACcAOlVCAAIQgAAEEMC8AxCAAAQgAIEEEEAAJwA6VUIA\nAhCAAAQQwLwDEIAABCAAgQQQQAAnADpVQgACEIAABBDAvAMQgAAEIACBBBBAACcAOlVCAAIQgAAE\nEMC8AxCAAAQgAIEEEEAAJwA6VUIAAhCAAAQQwLwDEIAABCAAgQQQQAAnADpVQgACEIAABM4AAQQg\nkB4EDh06JB988IEULFhQmjZt6nR637598v7778vmzZuldu3acs011zjxL7zwglSqVEmuu+46WbJk\niWzcuNEPKW/evFK0aFG5+OKLpUiRIk78ggUL5PPPP5cuXbpInjx5/Hk5gQAEIhDwESAAgbQgsGXL\nFp9+FfgqV67s7++1117rxFl848aNnXgVpk7cjBkznOv777/fn8fyZRwqoH2ffPKJk+fbb7914qdP\nn+4vmxMIQCAygTyWrP9QBAhAIMUJbN26VcqWLSsqgGXDhg1y9OhRZzRs3V65cqWULFlSzjzzTGd0\nbKPen376SWy027ZtW5k8ebLceOONokJabNT8zjvvyMKFC6Vq1aqydu1ah9wll1wiR44ckW+++UZO\nO43VrRR/neheDAgwBR0DiBQBgdwgcODAAXnxxRdl+fLljqBr3bq1lC9fPlNV//nPf2TmzJlieW36\n+Oqrr5ZChQr587z11lsyb948Zzr55ptv9sf/+uuvMn78eEcI2zTym2++KX/6058cYTpnzhzp2LGj\nI3z9N+jJlVdeKZ07d3aiOnXqJGeddZasW7dOVqxYIRdccIHce++9TrrVedtttwXeyjkEIBCCAAI4\nBBSiIJBoAnv27JErrrjCEb62ZmsCdvjw4aJTvlK3bl2neSZATRAeP35czjjjDBk9erQjhD/++GMn\nvVu3bk6cXRQoUEDGjh3rxNsfGw0/8cQTzvXevXulV69eUqZMGbF1Ygs6Ne18hvuzc+dOZ7Rr6fny\n5XOyZawrW/0I4HDkiIfASQLME51kwRkEPENgxIgRjvC95557ZPv27TJ16lSxUWvfvn2dNpoA7dmz\np6PwZCPcXbt2+RWm7NqmmJ9++mlnKtiuf/75Z2cEm9HBmjVryvfff+9cmuC18lq1auVMH1ukKWBl\nDa+99prcfffd0qJFC6lXr54zejbFrerVqztZbWrbgq4HO5/8gQAEIhNgBByZD6kQSAiBzz77zKnX\nhOPrr7/ujHLz588vc+fOdeIXLVrkjIrr168vV111lROnSlP+0ait0R47dswZRWek33ffffLGG284\neU8//XQpXbq0c27rtTadbGHp0qXOZ4YwdS5++7N48WKxw4LdYyPe559//rdUcTSjixcvLt99950/\njhMIQCA8AQRweDakQCBhBGxa2MLbb78t7733nnNerVo159OUoGxUbMEUpzJCxlSwXe/evduJNqWr\njFClSpWM07CfNlK2bUQmSLOGPn36SIcOHZytR5YearuRtWfTpk2O8DchT4AABMITYAo6PBtSIJAw\nAn/84x+duh977DFnKtpGvCNHjpTZs2dL4cKFpUaNGk76/PnzndGxXUyYMMGZhp41a5ZkCGu7z0bC\nFjJGz85FmD82NW0bI2y6O2swoVuxYkUpUaJESOFr+W3t2jSjEb5Z6XENgWACCOBgJsRAIOEETCPZ\nwoABA2TMmDGie3HlhhtukGHDhjnxNvV82WWXOSPdyy+/XEzhql+/fvLll1+KbQcyBa5atWo5BjTM\nYEb37t3968dOAWH+XHjhhU6KbUHKbrBR97Zt2+S8887L7q3kh0BaEkAAp+Vjp9NeJ2Drtn//+9/F\npqLNupSNfG+//XZHGFvbbQ122rRpzr7c//73v068jV5t21K5cuWcEailmxA2pSgry7SkMyxXheu/\nCWsLX331VbgsYeNt77CFBg0ahM1DAgQgcJIAhjhOsuAMAp4kYGYizz77bGerUagG2rTv4cOHpVSp\nUqGSHYMapmRlW5WiBdvuZFPIV199tZjWc3aCTZeb5rVpYIdrS3bKIy8EUp0AAjjVnzD9g0A2CYwa\nNUoGDhzo2Ic2m89ugq0b2wjcps5NCBMgAIHoBJiCjs6IHBBIKwIPP/yws53p5Zdfdt1vc+iwfv16\n6dq1q+t7yAiBdCfACDjd3wD6D4EQBMyms20zcjNtbbebprUdgVuhQhRLFAQgEEAAARwAg1MIQAAC\nEIBAvAgwBR0v0tQDAQhAAAIQCCCAAA6AwSkEIAABCEAgXgQQwPEiTT0QgAAEIACBAAII4AAYnEIA\nAhCAAATiRQABHC/S1AMBCEAAAhAIIIAADoDBKQQgAAEIQCBeBBDA8SJNPRCAAAQgAIEAAgjgABic\nQgACEIAABOJFAAEcL9LUAwEIQAACEAgggAAOgMEpBCAAAQhAIF4EEMDxIk09EIAABCAAgQACCOAA\nGJxCAAIQgAAE4kUAARwv0tQDAQhAAAIQCCCAAA6AwSkEIAABCEAgXgQQwPEiTT0QgAAEIACBAAII\n4AAYnEIAAhCAAATiRQABHC/S1AMBCEAAAhAIIIAADoDBKQQgAAEIQCBeBM6IV0VeqOfNN9+Uo0eP\neqEptAECEIAABDxA4KyzzpJrrrkmIS3J49OQkJrjXOlbb70lo0ePljZt2sS5ZqqDAAQgAAG3BFau\nXClffvmltGzZUooXL+72tlPON27cOHn11Vfl4osvPuUyTvXGtBkB28j33nvvlfbt258qK+6DAAQg\nAIFcJPDpp5/Khx9+KNOmTZPKlSvnYk0ni/7+++/l+PHjJyPieJY2AjiOTKkKAhCAAASySeCLL76Q\nsWPHOjOV8RK+2WxizLMjgGOOlAIhAAEIQCA7BBYsWCAjR46UUaNGSdWqVbNza1LnRQs6qR8fjYcA\nBCCQ3AQWLVokQ4YMkWHDhkmNGjWSuzPZbL3nBLCt1e7YsSOb3SA7BCAAAQgkG4Fly5bJwIEDZdCg\nQXL++ecnW/Nz3F5PCODDhw9Lnz59pFKlSpIvXz4pWbKkFC5cWGrXri2TJ0/OcScpAAIQgAAEvEVg\n+fLl0r9/fxkwYIDUqVPHW42LU2s8sQbcsWNH2bJli8yaNUuqVavmCN/du3eLPaBOnTrJwYMHpUOH\nDnFCQjUQgAAEIJCbBFavXu0Muvr27Sv16tXLzao8XbYnRsCmdj5x4kSpW7euFClSRPLkyePs/2rQ\noIGjFff22297GiKNgwAEIAABdwTWrVsnPXv2lB49ekj9+vXd3ZSiuTwhgG2qee7cuSERz5w5U8qU\nKRMyjUgIQAACEPA2AbWrIfPni2zbJrJx40bp2rWrM7PZsGFDbzc8Dq3zxBS0LcK3atVKxowZI9Wr\nV5dixYrJrl27ZMWKFY7pyNmzZ8cBBVVAAAIQgEAsCeggV41qiBw5IvLzz5vld7/rIk888ag0atQo\nltUkbVmeEMC2BrB48WL9lTRf1q9f76wH26jX1n3tQdmUNAECEIAABJKHgG7p1UGV6CDK2rxFj86y\ndGk7qVixSfJ0Ipdb6gkBbH0sUKCANG7cOKi7q1atkv3797taqP/ss8/kq6++CirDIj7//HPH1qeN\nqjNCzZo15bTTTs7CW12BJslIhw/vB/8fGd8XfD+4/348dkykV6+acuyYvT+/6NFZLrroailXroq8\n+OIKqVJFZN++mhp/8v0qXHiVXp80CRmv9JUrE6cE5hkBrORDBrMJumHDBpk0aVLI9MDIsmXLhlVn\nN+PeW7dulaJFiwbekuncFMAi+aYgHT68H75M/zOBF/x/8P+R8f/RooWKUkeWbtdXpLMeLXUms4FO\nQ/tUKIuOgkXy5hWNO/kGHT5cRC9Ovl/xSl+y5NeTjYjzWdp4Q+rcubMztT116tQ4I6Y6CEAAAulD\noHVrkSlTrL+79Oiox/V63KXHifD66yK3355xlfhPUwprrY2+5JJL4t6Yk+P/uFdNhRCAAAQgkEoE\n1KbSb8J3j3arix7X6HFS+Fpf0b8yCicCAjiDBJ8QgAAEIJAjAja1rKu7enTV4wo97tPjZDB7SrpS\nSPiNgCfWgEePHq1q6qqnHiaYjdBbbrklTCrREIAABCCQaAIvvqjqVr8c0GZ008Oc2z+YqUmmePXs\ns5mi0v7CEwLYth4988wz0qZNG8cMZdangiGOrES4hgAEIOAdAuY/p127Q9og3fgrpt38SFDjfvgh\nKCrtIzwhgMePH+9s/7EtQBMmTEj7hwIACEAAAslCYJ/OOJcsqYu/0luPynp0Cmp6F10OPv30oOi0\nj/DMGvDw4cPFHDDs3bs37R8KACAAAQgkC4GiRc3SxmN6lNLDpp+Dg64yEkIQ8MQI2NplewhfffXV\nEE0kCgIQgAAEvEigZ89jajthgDatkB42Ag7Y2Ptbg9XPDiEMAc8I4DDtIxoCEIAABDxGQD3Fqiej\n42qlcJC2zIRuPz2CJ1SbNhV56CGPNd5DzQkm5qHG0RQIQAACEPAWAbXqK7Vq+VT4DtWGmdazjYCD\nF3h1VVHU0ywhAgEEcAQ4JEEAAhCAwEkC776bYUhDPS2ImZm0EbDalMwSSpRQfWhTiCZEJMAUdEQ8\nJEIAAhCAgBH49luRm26yM3VxJD/qoUNcyadHcFC3vwQXBBDALiCRBQIQgEA6E1i9WtTRjRGwbaJ6\nIabWXECP4DB3rqjTm+B4YoIJMAUdzIQYCEAAAhD4jcAeNet83nl2MUmPpXqM0KOgHpmD2Xhes0bk\n6qszx3MVngAj4PBsSIEABCCQ1gQOqI5VjRqG4GU95uvxtB7mNjBzmDlTpHnzzHFcRSeAAI7OiBwQ\ngAAE0o7A11+L/O531u3X9Jijxzg9iumROXz/vci552aO48odAQSwO07kggAEIJA2BE4K3+naZ1V9\ndoSvqjZnCZ9+ivDNgiRblwjgbOEiMwQgAIHUJmDazidGvjqvLG/oMV6PUkGdrlUrY0tSUBIRLgkg\ngF2CIhsEIACBdCBwQtv5A+2qrfuO1aNMyG4vWhQymshsEEALOhuwyAoBCEAglQmc2D70iXbxeT1s\nv2+5oO6aX9+dO3UHcOgtwEH5iQhPAAEcng0pEIAABNKGwNOq4Lx3r9qZlGf0eEqPCkF9f+45EXXf\nLsWLByURcQoEmII+BWjcAgEIQCDVCHTv/pV2yQSvmZnUYW6W8PrrIrffniWSyxwRYAScI3zcDAEI\nQCD5CSxcuFCOHjXTksP0qB7UoTE6G43wDcKS4whGwDlGSAEQgAAEkpfAkiVL5L77BmsHhuhRM2RH\nOnUKGU1kDgkggHMIkNshAAEIJCuB7777Tm699XH54YcntAsXhuzG738fMprIGBBgCjoGECkCAhCA\nQLIRWLVqlTz6aF8Vvv206ReFbf7zphBNyBUCCOBcwUqhEIAABLxLYI16Tejdu7d89VUvbeSlYRt6\nmkqIC0MPjMPeQ4J7Aghg96zICQEIQCDpCazXfUQ9evSQBx/sIocOXRGxP5s3R0wmMYcEWAPOIUBu\nhwAEIJAsBH788Ufp2rWrPPJIR2na9MqIzdblYSlbNmIWEnNIgBFwDgFyOwQgAIFkIPC///1POnfu\nLO3bt5e//a1xxCaPUJe/TD1HRBSTREbAMcFIIRCAAAS8S+Dnn392hO9tt90nbdteK9Gmlrt3925f\nUqllCOBUepr0BQIQgEAWAtu2bXOEb/Pmd8gddzTPkhp8+c9/BscRkzsEEMC5w5VSIQABCCScwI4d\nO6Rduy6yadNNMmVKi6jtuflmkdato2YjQ4wIIIBjBJJiIAABCHiJwK5du6R+/S66z7epNusOV017\n+21X2cgUIwII4BiBpBgIQAACXiGwd+9eueSSrrJ27R+0Sfe4atamTa6ykSmGBNCCjiFMioIABCCQ\naAL79+/XaWQTvvW1KW2jNqdqVZF160QqVoyalQwxJoAAjjFQioMABCCQKAIHDx6U7qrCPGNGbW1C\n+6jNGDhQVFCLnHNO1KxkyAUCTEHnAlSKhAAEIBBvAocPH5ZevXrJ//2fuRPsGLH6M88UmTtXLUCH\nNwEd8X4SY0MAARwbjpQCAQhAIGEEjhw5Ig880EdHvmXl4MEuEdvx8cci11wTMQuJcSKAAI4TaKqB\nAAQgkBsEjh07Jo0a9VfHCsW1+B4Rq1i9WqRGjYhZSIwjAQRwHGFTFQQgAIFYEjDhe8UVT8jChXm1\n2D56hFfr6dYN4RtL9rEoK/zTikXplAEBCEAAArlCYM4cn1Sp8qQK3yNafn89Tg9bjylZjRwZNpmE\nBBFgBJwg8FQLAQhA4FQJmMWqGTOG6+279RiqR/iv8gIFTmwzOtW6uC/3CDACzj22lAwBCEAg5gQ6\ndTLhO0rL3aLHED1s+jl0qK27kdQmB8GjBBDAHn0wNAsCEIBAVgK//ioyduw4jd6gh41882XN4r/+\nz39Eli3TienwM9P+vJwkhkD4eYvEtIdaIQABCEAgDIFy5Z7TlOV6PKVHwTC5RA4cELGpZ4K3CTAC\n9vbzoXUQgAAEHAIFC06WI0cW6rlNPxcKSaVIEYRvSDAejUQAe/TB0CwIQAACGQSqVZuiBjbm6eVo\nPVTKhggvvSSi3gcZ+YZg49UopqC9+mRoFwQgAAEl8Mor09RZwnt6NlYPM7YRHP70J5F77w2OJ8bb\nBBDA3n4+tA4CEEhjAu+884507fovJWDCt2RIEmedJTJrVsgkIj1OgClojz8gmgcBCKQngdmzZ8uU\nKVPkl1/GKIAyISGUKyeydWvIJCKTgAACOAkeEk2EAATSi8BHH30kkydPlgoVTPieHbbzP/0UNomE\nJCDAFHQSPCSaCAEIpA+BTz/9VJ577jlp02aM3Hxz+bAdb9gwbBIJSUIAAZwkD4pmQgACqU/giy++\nUEMbY+WJJ0ZL3bqVI3Z4xoyIySQmAQEEcBI8JJoIAQikPoEFCxbIqFGj1GnCSGncuGrEDp95pkip\nUhGzkJgEBDy7Bnz8+HG1YYoR0yR4h2giBCCQQwKLFi2SIUOGyNChQ2Xnzhry44+RC9xiZqAJSU/A\nEwJ49+7dzi+/G2+8UebOnauGxmfI2WefLeXLl5eHHnoIQZz0rxkdgAAEwhFYpgabBw4cKIMGDZLz\nzz9f6tcPl/NEvA6UJV94E9CRbybVUwQ8IYCHDRsmNv3yJ91N/re//U169erlCOH169fL4cOHZdq0\naZ6CRmMgAAEIxILA8uXLpX///jJgwACpU6eOXH115FI7dpSoAjpyCaR6iYAn1oBts7kJ4MKFC+ue\ntq2ybds2adCggcOpd+/e0qVLF7n//vu9xI22QAACEMgRgdWrV0ufPn2kb9++Uq9ePVHlZ+eIVOg4\nc4RESBkCnhgBX3DBBWL73nbt2iWfffaZ/Pe///UDXrp0qVxyySX+a04gAAEIJDuBdevWSY8ePaRn\nz546oq2vgw6JOvr93/+Svde0PysBT4yAu3btKm3btpW1a9fKX//6V9mzZ4+YUL7ooovk3//+t8yb\nNy9ru7mGAAQgkJQENm7cqOYlu0rnzp2dmT7VN5UyoQ1d+ft31VUiZcv6LzlJEQKeEMA23WxrIdu3\nb1fV+lJy6NAh+eCDD1QbcKdjDaZgwfB+L1PkOdANCEAgDQhs3rzZWVJ79NFHpVGjRk6PbV03WlDd\nVEIKEvCEADauefLkcYSvnefPn19uuukmO5VVq1bJ/v37nTUSJyLCH9OmttFzqGBbmo4ePRoqiTgI\nQAACuU5gi+4dMn2Wdu3aSZMmTfz1Pfus/zTkia7M6fdjyCQik5yAZwRwOI6mAb1hwwaZNGlSuCz+\n+Pfee0/effdd/3XgycKFC6WcWS4nQAACEIgzgV9++cWZcm7durU0a9bMX3s0pSqdCJRixfzZOUkx\nAnl8GlKsTyG7Y+st9gt06tSpIdOJhAAEIJAbBGxpzXRbbrnlFrn11lv9Vdh+3ssv918GnZilK1PO\nIuQuAVuPtx9GiVD29YQWdCBemybesWNHYBTnEIAABJKSgO3s6NSpkzRv3jyT8NWl4IjC1zqrZqEJ\nKU7AEwLYjG3YfrhKlSqphZd8UrJkSWdPcO3atR0lrBR/BnQPAhBIQQKmj2Izb7bee9ddd/l7aP57\nK1b0X4Y8KVFCpGbNkElEphABT6wBd1Q1QJsenjVrllSrVs0RvqZQZZrR9uvx4MGD0qFDhxTCTlcg\nAIFUJrBv3z5H4aqh+gxs06aNv6uffy6q/ey/DHvy669hk0hIIQKeGAF/+OGHMnHiRHW/VVeKFCni\naEQXL17c2SNnrrnefvvtFEJOVyAAgVQmcODAAenWrZuzc8M0njOC6oi6Er6qcyqneeKbOaPlfOYW\nAU88ZptqNicMocLMmTN1k3qUXeqhbiQOAhCAQJwJmA0Ds3BVU+ePH3nkEX/tqoeltu79l2FPvvlG\npHJkN8Bh7yUh+Qh4YgraPIG0atVKxowZI9WrV1e1+2KOWcoVK1Y4e3dnz56dfGRpMQQgkFYETJfF\nbNdXVglqS2eB4ZxzAq9Cn6tJfJ0FDJ1GbGoS8IQANkPkixcvlvnz54t5QLL1YBv12rqvWYsxIx0E\nCEAAAl4lYLs3HnvsMceYkE0/W1B5rNrPoqZ0RQcSTlTYP7Y78jfbQ2HzkJB6BDwhgA1rgQIFpHHj\nxqlHmB5BAAIpTeDYsWOOO8FChQo5I+CMAYOqs8iRI9G7rmMPufji6PnIkXoEPLEGnHpY6REEIJAO\nBI6rJ4VBgwY5s3T9+vVT5akTX6m6kuZK+KruKcI3HV6UMH30zAg4TPuIhgAEIOBJAmZEcOjQoWJa\nz4MHD5bTTz/daeecOaKe3aI3OW9ekYceip6PHKlLAAGcus+WnkEAArlIYOTIkY4HNxPCeU2aalC9\nUWna1F2l6h+GkOYEmIJO8xeA7kMAAtknYDs2zLXgk08+6VjvsxLUjpBceGHksmyQrN5XHQUtNfpH\nSHMCjIDT/AWg+xCAQPYITJgwQVavXi2jR492lEft7q+/FrnhhujlRNOGjl4COVKJAAI4lZ4mfYEA\nBHKVgLlFXbp0qSN8CxYs6K/rd7/zn4Y9MQtXBAgEEkAAB9LgHAIQgEAYAi+//LJjq+Dpp592TOZm\nZPv97zPOwn+qxzssXIXHk7YpCOC0ffR0HAIQcEvgtddekzmq3jxu3DjHUl/GfWY848svM65Cf+pM\ntTpmCJ1GbHoTQACn9/On9xCAQBQC06dPl3fffdcRviXMT+BvYcYM0fiMq9Cf5ve3fPnQacRCAAHM\nOwABCEAgDAFzBvPGG2/I+PHjHTOTgdlatgy8Cj6/7DKEbzAVYgIJIIADaXAOAQhA4DcCH3zwgbz0\n0kuO8M3qkU3N10e072zm600zmgCBSAQQwJHokAYBCKQlgU8++USef/55Z9q5bNmyfgZ79ojY5f79\n/qiQJ+qVkACBqAQQwFERkQECEEgnAp9//rk888wzjnvUChUq+Ltue3jVU2rU8OqropaxomYjAwQE\nS1i8BBCAAAR+I2AuUc3KlZmZrFKlSiYuxYtnugx5UamSqG/zkElEQiCIACPgICREQAAC6Uhg4cKF\nMmLECBk2bJhUN3dGAeHRR6NPO5uZyY0bA27iFAJRCDACjgKIZAhAIPUJLFmyxPFoZLada9asmanD\nugVYnn02U1TQRcmSIvv2BUUTAYGIBBDAEfGQCAEIpDqB7777Th5//HF54okn1JlCZm8K5rHorrsi\nE7jzTpFffxXJnz9yPlIhkJUAAjgrEa4hAIG0IbBq1Srp27ev9OvXTy666KJM/T52TKR06UxRQRem\nozV1alA0ERBwRQAB7AoTmSAAgVQjsGbNGunVq5f07t1bLr300kzdM01mcxcYbTvRjz9muo0LCGSL\nAEpY2cJFZghAIBUIrF+/Xnr06CFd1UvC5Zdf7u+SjXovuUTU45E/KuzJwYNhk0iAgCsCCGBXmMgE\nAQikCoEfddhqgrdjx45y5ZVX+ru1YoVI3bqRLVxlZB4yhDXfDBZ8njoBpqBPnR13QgACSUbgp59+\nUs9EXaR9+/bSuHFjf+v/9jdRBSx3wtf2+uqsNQECOSbACDjHCCkAAhBIBgI///yzI3zbtGkj1157\nrb/Jt98uMm2a/zLiSY0aIqtXR8xCIgRcE2AE7BoVGSEAgWQlsG3bNuncubPccccd0rx5c383HnzQ\nnfA15wpjxyJ8/eA4iQkBRsAxwUghEICAVwns2LHDGfnedNNN0qJFC6eZu3eLKl+JrFwZvdVml+Pb\nb0XO4NsyOixyZIsAI+Bs4SIzBCCQTAR27drlCF+bcrbRr4VFi0TMrrMb4TtgwIl8CN9keurJ01Z+\n0yXPs6KlEIBANgjsVTNWpu38hz/8Qe6++27nzq1bRff8Ri/EppzVQJZccEH0vOSAwKkSQACfKjnu\ngwAEPEtgvzrsNeFbv359adu2rb+dAa59/XGhTo4fDxVLHARiS4Ap6NjypDQIQCDBBA6qhYzu3btL\n7dq1ne1G1hxzlGDeitwEs/9MgEA8CCCA40GZOiAAgbgQOHz4sGNe0twJmqENC2bdqkgRkWijWlO2\nMqcKhQvHpalUAgFBAPMSQAACKUHgyJEj0qdPHymr88xmbCMjlC+fcRb+01wOmlKWuRUkQCBeBBDA\n8SJNPRCAQK4ROKbD3P79+6t2c3HHxnNGRePGiaj9jYhBZ6tVQzpiFhIhkCsEEMC5gpVCIQCBeBEw\n4Wu+fPPmzeuMgE877cTX2ltviZiJyUjBTEGPGBEpB2kQyD0CCODcY0vJEIBALhPw+Xzy5JNPik0/\n2wj49N80rT75ROTWWyNXfv31Ip9/HjkPqRDITQJsQ8pNupQNAQjkKoFhw4bJnj17ZIi6JzrjN2sZ\n8+eLNGkSudoCBURmz46ch1QI5DYBBHBuE6Z8CEAgVwiMGjVKtqpljRE6h2zTzxbWrhVp2DB6ddu3\nR89DDgjkNgEEcG4TpnwIQCDmBMapdtWGDRsc4ZsvXz6nfNtmpLuPogazcFWwYNRsZIBArhNAAOc6\nYiqAAARiSeC5556T5cuXy1NPPaWC9KQkLVMmei3vvnvC72/0nOSAQO4TQADnPmNqgAAEYkRg8uTJ\nsnDhQnn66aelUKFC/lLvuUck2rTy3/8ucsMN/ls4gUDCCSCAE/4IaAAEIOCGwJQpU2TevHli089F\nzLTVb0GVoOWf/8y4Cv350Ucif/xj6DRiIZAoAgjgRJGnXghAwDWBadOmyXvvvSdjx451jG1k3Ki7\nj+SxxzKuQn/aPl+Eb2g2xCaWAAI4sfypHQIQiELgnXfekX/961+O8C2ZxVak+luIGGygbJauCBDw\nIgEEsBefCm2CAAQcArN1s65NPdu0c5ksWlaLFol8/31kUOZcgQABrxLAEpZXnwztgkCaE/hIF25N\n6WrMmDFy9tln+2m8+KJIpUoil17qjwp5snOnyG87lEKmEwmBRBNgBJzoJ0D9EIBAEIFPP/1UbLuR\nCd/yAe6Mzj9fZNWqoOxBEe3aia4VB0UTAQFPEUAAe+px0BgIQOCLL75w1ntHjx4tlStX9gOxEa8b\n4Ws3TJzov40TCHiWgGcF8MGDBx3D6hkm5jxLkIZBAAIxI7BgwQIxE5MjR46UqlWr+sutVUvU+Ib/\nMuKJDpzlN4dIEfORCIFEE/DEGvDGjRvl3nvvdTbY//LLL/LAAw84TrVLlCghbdu2lcOHDyeaE/VD\nAAK5TGCRalWZU4WhQ4dKjRo1/LVddJF74dujh0j79v5bOYGApwl4QgCbGzGbaqqlP3PHjx8vR48e\nlW+//VaWLl3qeDoZNGiQpyHSOAhAIGcEli1bJgMHDhT7Xz/fFnp/C4MHi34PZFyF/zznHNEf8CLD\nh4fPQwoEvEbAE1PQn332maxcuVI1FvM5+/3efvttqVixosPK/iEffvhhr3GjPRCAQIwImF1n+xE+\nYMAAqVOnjr9UE7z9+vkvw57Y1PQFF4RNJgECniXgiRHweeedJy+//LID6eqrr1Y/nScddc6cOVPO\nPfdczwKkYRCAwKkTWL16tfTp08c56tWr5y/I3Ara1HOkYO5/jx1D+EZiRJq3CXhiBDxhwgQ1kn6D\nvKgb/Gztp1u3bvJ///d/qkhxmuzevVtshEyAAARSi8C6deukhy7a9uzZU+rXr+/vnK5AuXIraGYo\nCRBIZgKeEMDV1YmnTUPZxvtVus/A1oPPPPNMZ+TbvHlzOcN+6hIgAIGUIWCKl127dpXOnTtLgwYN\nMvWrQIFMlyEvPvggZDSREEgqAp6RbHny5JFrr73WOZKKII2FAASyRWDz5s2O8H300UelUaNG/ntt\n5Js/v8jx4/6okCf33y/6PREyiUgIJBUBzwjgcNRsRLx//34JXB8Kl/f555937MaGSl+zZk2mfYWh\n8hAHAQjkLoEtW7ZIly5dnK2GTZo0yVSZTnpFFb66MiUmgAkQSAUCnhfA5oZsw4YNMmnSpKi8H3ro\nIbEjVLCpLvvnJ0AAAokhYHv87f+wdevW0qxZM38j9uwR+f3vRfbu9UeFPDHrVgjfkGiITFICnhPA\ntgd4j/5H2hqwhceiOftMUvA0GwLpRGD79u2O8G3ZsqXcdNNN/q5v3SpqdMd/GfbEBsthfluHvYcE\nCHidgCe2IZmlK9uKUEldnNheYPP5WbhwYamtzj7NGwoBAhBIXgK7du2STp06iSlU3nrrrU5HZs0S\n/f92J3xbtBCZMyd5+0/LIRCOgCdGwB07dnSmh2fpf2W1atUc4Wvbj0wz2v5xzS50hw4dwvWBeAhA\nwKMEbDbLpp1tvfeuu+4Sn0+kXDkRG/m6CaVKiUyf7iYneSCQfAQ8MQL+8MMP1XvJRKlbt64UKVJE\nTCO6uPoSs+0JY8eOFbOMRYAABJKLwL59+xyFq4YNG0qbNm1ETT07ThLcCl/rrS4bEyCQsgQ8IYBt\nqnnu3LkhIZslrDJlyoRMIxICEPAmgQMHDjgGdWz3Qjt1zvvxxyLmTtBtOOssUScsoj/G3d5BPggk\nHwFPTEGbEfZWrVo5zrfNKEexYsXE1o1WrFjhOGYINE2ZfIhpMQTSi8ChQ4ccC1c1a9aURx55RLcR\nivzxj+4ZmO8VdC/d8yJn8hLwhAC2X8mLFy+W+fPny/r16531YBv12rqvbdS3KWkCBCDgfQKmUNm7\nd2/Hmp3pb1iwdVw3wTY+2D7fW25xk5s8EEh+Ap4QwIaxgNqfa9y4cfITpQcQSFMCtoXQtg2WUolr\n9twtqHtfVaKMDOSyy0SmThW1Ax85H6kQSDUCnhHAqQaW/kAgnQgcU7dE5k6wUKFCzgjYZq1M47lv\n38gU3npL5M9/jpyHVAikKgFPKGGlKlz6BYF0IHBcjTeb324Tuv3Uga95MbMQbSp5yhSEbzq8H/Qx\nPAFGwOHZkAIBCEQh4NNh7tChQ8W0np988kk5/fTTnTvUsZnMmBH+ZrN+pduCCRBIawII4LR+/HQe\nAjkjMHLkSDEzkyaEM9yGqt+TqN6KvvsuZ/VyNwRSgQACOBWeIn2AQAIIjBkzRsy14PDhwx0TstYE\nXQpWP96RG9O0qai52ch5SIVAOhBAAKfDU6aPEIgxgQkTJsjq1atl9OjRzg6GjOLdbDl6//2M3HxC\nIL0JIIDT+/nTewhkm4C5Bl26dKkjfAsWLOi/v359UQM6/suQJzpgdsxRhkwkEgJpRgABnGYPnO5C\nICcEXn75ZcdgztNPP+3Ybc8o69prRRYuzLgK/fnJJyLly4dOIxYC6UiAbUjp+NTpMwROgcBrr72m\nbgHnyFNPPeWYi7UitmwRUV8LYlrPkYI6PFNDO5FykAaB9CMQdgRsJuXMPOR3qq64du1atVJTQ42p\nXyoXXHCBX9sx/XDRYwikJ4Hp6hPw3XfflXHjxkmJEiUcCHfcIfLGG9F5VK0qel/0fOSAQLoRCBoB\n276+qWoXrlatWo4fz0903sj2+JlDhNatW0vFihWdtR8T0AQIQCD1CZhHsjdU0tq0s5mZtGCazG6E\nr+XV3+8ECEAgBIFMI2BzfN+yZUv952rqrPOULl066JYdO3bIs88+K9dcc43zT1meRZ0gRkRAIFUI\nfPDBB/LSSy/J+PHj/W5BV64UnYp210P9uiBAAAJhCGQSwLaR3tZ5ihYtGia7yJnqsqSvGnjt0qWL\nmAk6AgQgkJoEbPbr+eefd6ady5rpqt+CrkK5Cvv2idqGdpWVTBBISwKZpqBNAGcIX1vv2bZtWyYo\nP/zwgzP9bJG2/aBw4cKZ0rmAAARSg8Bnn30mttfXFK4qVKjg79TZZ/tPw56YkY2dOxG+YQGRAIHf\nCGQSwIFUtm7dKuand968eU70iy++KL/73e/8AjowL+cQgEDqEDC/3LbeO2LECKlSpYq/Y336iPz8\ns/8y5IkuF8uvv4oULx4ymUgIQCCAQKYp6IB4adeunfPPd++990q5cuXUxNwxRxjXqVMnMBvnEIBA\nChFYqJt5TfAOGzZMqlev7u+Z/Q5Xc88Rw7JlIrVrR8xCIgQgEEAg7Ag4I0/evHnVofZBx9VYhpux\njDQ+IQCB1CGwZMkSGTx4sOPVqGbNmv6OaXTUPbzduyN8/cA4gYBLAmEFsK3/3HPPPTJ27Fj55ptv\n5OGHH5arrrpKXnjhBZdFkw0CEEgWArbf//HHH5cnnnhCLrzwQn+zTc9SV6IihmLFREfNEbOQCAEI\nhCAQVgBXrlzZsfd6ww03OLc98MAD8p///EdtvUYx9hqiEqIgAAHvEli1apWzs6Ffv35y0UUXZWqo\nm7VcvhIyIeMCAq4JZBLAZnBj7ty5zs033nijnHXWWZkKsjWhrl27OnGff/65/O9//8uUzgUEIJBc\nBNao895evXpJ7969HUt3ga2339579wbGBJ//9FNwHDEQgIA7ApkEcL58+RxFqxYtWjjWsLKOdn/8\n8Ud56623pFGjRk56xpYld1WRCwIQ8BKB9evXS48ePZwf1ZdffnmmpqntDZk1K1NU0IWuTKmCZlA0\nERCAgEsCmbSgTz/9dGcN6Cf9WTtgwADpqBbUzdjG2br5z4Sv2YBt0KCBsz8QbWiXhMkGAQ8S2LRp\nk3Tr1s35H7/yyiuDWnjffUFRmSJGjhSpWzdTFBcQgEA2CWQSwBn3mnlJ8/lpxxZ1d5LhjCHrlHRG\nfj4hAIHkIWA/sG0p6aGHHlLt5swuilasiK50VamSqPBOnv7SUgh4lUBIARzYWDNBF2iGLjCNcwhA\nILkI/KyWNMyMbJs2beRac+IbEGbMELn55oCIEKd58ohs3BgigSgIQCDbBMIK4J1qS+6RRx6RZbq7\nPtDz0fXXX+9Yycl2TdwAAQgklICZlu3cubPcoX4Emzdvnqktpk8ZTfjaDboSRYAABGJEIKwANms4\npoRl/j+LFCnir66kGXolQAACSUXAvJjZyPemm24SU7LMGqpVyxoTfH3nnSI4PwvmQgwETpVAWAG8\nefNmZwScdY3oVCviPghAIDEE7Ie0jXyvu+46Z/SbtRW2nqvG7iIG84CkbsIJEIBADAlk2oYUWO6f\n//xneeWVV9T4ehTr64E3cQ4BCHiKwF7dyGsKV7Z1sHXr1pnaZv/a+m+uHs4yRQdd6HKxLF8eFE0E\nBCCQQwJhBbBpSs6ePdtxxHDuuefK+eef7xydOnXKYZXcDgEIxIPA/v37HeFbv359adu2baYqbRuR\nuRb8178yRQdd2Kj3H/8IiiYCAhCIAYGwU9BmgvKyyy4LqoI14CAkREDAcwTMgUp39ZBQW90TtW/f\nPlP7bNvvF19kigp5UaaMiK37EiAAgdwhECSAmzVrJlOmTJFKutmvdOnSYpqTdk6AAASSg4DtWjDz\nkmY61ozpBAYdDIt6HHQV1CU4AQIQyEUCQVPQ5g/0yJEjTpULFiyQVq1a5WL1FA0BCMSSgP3v9unT\nx9m7b1rPgWHQIPfC1xws2J5fAgQgkHsEggRw7lVFyRCAQG4SOHbsmPTv31+Kqwsjs/EcGMzClSZF\nDeb9yAxtmItBAgQgkLsEEMC5y5fSIRAXAiZ8zZdv3rx5nRHwaaed/Ne2Ca0AF79h22OWsNT+ji45\nhc1CAgQgEEMCQWvAVrY5XjAlDrMDfejQIdmwYYO/ysKFCztrw/4ITiAAgYQS8Pl88uSTTzpLR4MH\nDxZzqhIYrrgi8Cr43GS1TTkH2NsJzkQMBCAQcwIhBXBW7edzzjnHX/Ftt90mb7zxhv+aEwhAILEE\nhg0bJnv27JEhQ4YECV/Vx5JFiyK3T2+VQoUi5yEVAhCIPYEgAbw1iupjHjQzYv8UKBECp0hg1KhR\nYv+zZjrWpp8Dw9GjIvnzB8YEn9vaMMI3mAsxEIgHgSABnHX6Kh6NoA4IQCD7BMxOuy0PmfDNly9f\npgJ0SVjjMkUFXdgGB7WvQ4AABBJEIEgAJ6gdVAsBCGSDwHPPPafmIZfLU089JQULFgy607SZdWk4\nYnj11YjJJEIAArlM4KSqZC5XRPEQgEBsCEyePFmNaSwUm34ulGX+2NZzK1QQ2bcvcl1ffhk5nVQI\nQCD3CSCAc58xNUAgZgTMSt28efPUgcLoTG5CrQKbUrb9u2rGPWIYPlykQYOIWUiEAATiQIAp6DhA\npgoIxILAtGnT5L333pOxY8c6xjYCy1Srk7J2bWBM6PNbbxU10hE6jVgIQCC+BBDA8eVNbRA4JQLv\nvPOOTJ8+XcaPHy9ZHaKMGeNO+Kppd1EZToAABDxCAAHskQdBMyAQjoC5BbWpZ9N6NgcpWYM6PYoa\nztD/9Cg7DKOWQQYIQCC2BFgDji1PSoNATAl89NFHYkpXY3SYe7Y58M0SzNCGbTmKFCpXPqGUFWCd\nMlJ20iAAgTgRQADHCTTVQCC7BEzZauLEiY7CVfny5UPeHs3MpMpuWb8++p7gkIUTCQEI5CoBpqBz\nFS+FQ+DUCHzxxRfOlLNpO1e2IWyYsHhxmASNnjpV5M47w6eTAgEIJJaAZ0fA5gxi9+7diaVD7RBI\nAAHzw217fM3CVdWqVcO2IIu736B8CN8gJERAwFMEPCuA33rrLcnqUNxT5GgMBHKBwCL1nGBOFYYO\nHSo1atQIW4MOkHVdOGyy3HRT+DRSIAABbxDwxBT0ueeeK9u2bctE5LBqlxxVa/ImiG+55RZHESVT\nBi4gkGIEli1bJgMHDpRBgwapjebwRprNO+iVV0bu/GuvRU4nFQIQSDwBTwhg0/Js27at3H333dKm\nTRuHyttvvy3z58+X4Wq2x3wQEyCQygTMrnO/fv3k8ccflzp16kTsau3aEZOlUiVR+9CR85AKAQgk\nnoAnpqCv1J/zZtt2zZo1zrSzCVzb71hEPYRXqVIl5N7HxKOjBRCIDYHVq1dLnz59pG/fvnLxxRdH\nLHTLFpG9eyNmkU8+iZxOKgQg4A0CnhgBG4piasT25ZdfljfeeEMaNWokl19+eZBzcW8goxUQiB2B\ndevWqWnIHtKzZ0+pX79+1ILvvTdylgcfFF07jpyHVAhAwBsEPDECDkRx++23y4cffuisCZctWzYw\niXMIpBSBjRs3SteuXaVz587qHMGddwS1yxE2mHB+/vmwySRAAAIeI+CZEXAgl4oVK8q7777rRK1a\ntUr2798v9erVC8wS8vzrr7+WJUuWhEwzBZf8+fOHTCMSAvEmsHnzZkf4Pvroo86Mj5v6dUtwxPDS\nSxGTSYQABDxGwJMCOJCReYDZoGqfkyZNCowOeZ43b96wCluWlidPnpD3EQmBeBLYogu5tsXugQce\nkCZNmkSt+tAhEZXT8uKL4bPWrBk+jRQIQMCbBDwvgB977DHX5EyBJZwSi42O7YuPAIFEEvjll1+c\nKefWrVtLs2bNojbFHChUqBDd3vPSpVGLIgMEIOAxAp5bA7a9vzt27PAYJpoDgZwT2L59uyN8W7Zs\nqYYy3FnKMDWIaM4W8uXD1nPOnw4lQCD+BDwhgM3ohm3DqKQbGPPpt4n5O7WtSLV1w6PtESZAINkJ\n7Nq1Szp16iQ33HCD3Hrrra6607y5q2zyzTfu8pELAhDwFgFPTEF37NjRmR6eNWuWVKtWzRG+Zgfa\njBPYl5bZhe7QoYO3yNEaCLgksGfPHmfka+u9d2bDQLO6AY4a7PdpBKNZUe8nAwQgkDgCnhgB27Yj\nc7tWt25dx/iGKUsVL17c2ZoxduxYMatYBAgkI4F9+/Y5ClcNGzb0W3lz049oW4L130NWrhS57z43\npZEHAhDwIgFPCGCbap47d25IPjNnzpQyZcqETCMSAl4mcODAAenevbuzha5du3aum6p2OdQyXPjs\n1auL7NwpguZzeEakQCAZCHhiCtoM0Ldq1Uq9u4yR6vrtYlaxbM1sxYoVjkOG2W7m4pKBNm1MGwKH\ndO+QWbg677zz5JFHHnHdb7VKKSNHRs6uJtIJEIBAChDwhAA2IxuL1bO4OV9Yv369sx5so15b9zWz\nlOzfTYE3LY26YEqFvXv3lsqVKzs6DG67PmGCyF/+Ejn3VVeJzghFzkMqBCCQHAQ8IYANVYECBaRx\n48bJQY1WQiAMAdtGZ3vXS5UqJd26dQuTKzhafY6IWqaMGubNi5qFDBCAQJIQ8IwAThJeNBMCYQkc\n0w27AwYMkEKFCjkjYLczN+r0S1RXK2pge3xURGSAQFIR8IQSVlIRo7EQCEHg+PHjYroMJnTNr+9p\np7n71xoyxJ3wnTJFpESJEBUTBQEIJC0BRsBJ++houFcI+Hw+GTp0qLNf/cknn8yWG011ARw1PP64\nyF13Rc1GBghAIMkIIICT7IHRXO8RGKlqy2Zm0oTwGWdE/5fSZWI1MiO6PzhyX8x3iNl41l16BAhA\nIAUJRP+2SMFO0yUIxIqAbZ0z14LDhw93zKhGK1cde8lll4n6u46cU3USRY3BiTrxIkAAAilKAAGc\nog+WbuU+gQm6b2i1btwdrY56TYs/WlCLlHLOOdFynUhX510IX3eoyAWBpCWAAE7aR0fDE0nA/FMv\n1flhE74FCxaM2hRdJlYnI1GzORl0BxPTzu5QkQsCSU0AAZzUj4/GJ4LAyy+/7BiNefrppx3b5dHa\nYMK3aVNRq27Rcp5IX7vWXT5yQQACyU0AAZzcz4/Wx5nAa6+9JnPmzJFx48Y5JlOjVb9unchFF4nY\n9HO0YH59Fy0SLTdaTtIhAIFUIOBus2Iq9JQ+QCCHBKZPny7vvvuu2Mi3hItNuWbZSr1rRhW+ZgVL\n5bqjdFWrVg4bye0QgEDSEGAEnDSPioYmkoB55XrjjTdk/PjxupYbfTH3yBERE6zRwplnito/j5aL\ndAhAIBUJIIBT8anSp5gS+OCDD+Sll15yhK9b15gPPOCuCW7sP7sriVwQgECyEUAAJ9sTo71xJfDJ\nJ5+IaTyPHTtWypYt67ruV16JnvWZZ0SVuKLnIwcEIJCaBFgDTs3nSq9iQOCzzz4T2+trW40qVKjg\nusRo7gIvvVRUkUvk0UddF0lGCEAgBQkggFPwodKlnBMw39SmbDVixAhdy3WxmKtVvvCCSOHCka1c\nlS4tsnChSJMmOW8jJUAAAslNgCno5H5+tD4XCCxUCWmCd9iwYVK9enVXNTRsKLo3OHrWxYuj5yEH\nBCCQHgQQwOnxnOmlSwJLliyRwYMHyxD1E1izZk1Xd115pTvh26aNSMWKrookEwQgkAYEEMBp8JDp\nojsC3333nQwYMEAGDRokF154oaubTKh+8UX0rOXKifzjH9HzkQMCEEgfAqwBp8+zpqcRCKxatUr6\nqnPe/v37S926dSPkPJn01VciapUyajCPRj/9FDUbGSAAgTQjgABOswdOd4MJrFmzRnr16iW9e/eW\nS01F2WX4wx+iZzQzlPv3R89HDghAIP0IIIDT75nT4wAC69RYc48ePaRr165y+eWXB6SEP508WaR+\n/cjOFWzUu2mTiC4pyxks9ISHSQoE0pgAXw1p/PDTveubVEJ2795dOnbsKFeaJpWL0Lq1yJQp0TMe\nOiSSJ0/0fOSAAATSlwAj4PR99mnd8590UdZGvQ899JA0btzYFQu1y+FK+OoyMsLXFVEyQSC9CSCA\n0/v5p2Xvf/75Z+nSpYu0URXma6+91hWDY8dErrrKVVbp3NldPnJBAALpTQABnN7PP+16v23bNhWQ\nneWOO+6Q5s2bu+q/bTNyu45r68MuPBW6qpdMEIBAahNgDTi1ny+9CyCwY8cOR/jefPPN0qJFi4CU\n8Ke2fcjN8nCBAiJTp4rcckv4skiBAAQgEEgAARxIg/OUJbBr1y5H+DZr1kxuv/121/10Y4nyn/8U\nufNOkdNPd10sGSEAAQgIApiXIOUJ7N2711G4atSokbQ2NWaXwaaTDx6MnFkVqLXMyHlIhQAEIBCK\nAGvAoagQlzIE9qsVDNN2rq8bd9u2beu6Xx99JJo/cvYbbxQZNy5yHlIhAAEIhCOAAA5HhvikJ3BQ\nh6+2z7d27drSvn171/3ZvVtUOzpy9gceEJkxI3IeUiEAAQhEIoAAjkSHtKQlcPjwYce8pLkTNEMb\nbsPx4yLFi0fObRrR5vuXAAEIQCAnBBDAOaHHvZ4kcOTIEenTp4+ULVvW2e+bnUaWKhU996+/Rs9D\nDghAAALRCCCAoxEiPakIHFOLGebRqLgOY83Gs9vw97+f2Ou7c2fkOxYvFilWLHIeUiEAAQi4IYAW\ntBtK5EkKAiZ8H3/8ccmXL58zAj7ttOi/L48eFalaVeTHH6N3Ud0Ey8UXR89HDghAAAJuCET/hnJT\nCnkgkGACPp9PnnzySTmqEtVGwKe72JQ7Z46IeS1yI3xbthR57LEEd5LqIQCBlCKAAE6px5m+nRk2\nbJjs2bNHBg4c6Er4vvaaSNOm7niVKSPy5pvu8pILAhCAgFsCTEG7JUU+zxIYNWqUbN26VUaMGKEj\nWh3SRgmrVoncdVeUTL8lFywoWra7vOSCAAQgkB0CjICzQ4u8niMwTi1hbNiwQYYOHeqs/UZr4Lx5\nIuefHy3XCXeCDz4oOqrGtWB0WuSAAAROhQAj4FOhxj2eIPDcc8/J8uXL5amnnpKCNlSNENQgltSp\nI7J2bYRMvyWZ28GPP8a2c3RS5IAABHJCgBFwTuhxb8IITFZDzQsXLhSbfi5UqFDYdqgZaBk7VqRw\nYXfC19aGbZTsQocrbJ0kQAACEHBDgBGwG0rk8RSBKVOmqJCcp3aYx0mRIkWC2nbokIgJ0gkTRIW0\niCpIuwrffCNSt66rrGSCAAQgkGMCCOAcI6SAeBKYNm2azJ49W8aPH+8Y28hatwnRyy8XMSGcnXDP\nPQjf7PAiLwQgkHMCTEHnnCElxInAO++8I9OnT5enn35azjzzzKBaf/nlhKGM7ApfdZYkL78cVBwR\nEIAABHKVACPgXMVL4bEiYKNem3q2aefSpUuHLFb9LmQrmFWrTz4RFebZuo3MEIAABGJCAAEcE4wU\nkpsEPlLnvKZ0NVa1qc4+++yQVekWYGfLUMjEEJFmBatJkxAJREEAAhCIEwGmoOMEmmpOjYApW02c\nOFFGjx4t5cuXD1nIn/8s0rNnyKRMkaYs/de/itiWJIRvJjRcQAACCSDg+RGwGdg3+7758+dPAB6q\nTCSBL774wplyNuFbuXLlTE3R10JuvllUISu6lnPz5iKDB4tceKGosY5MxXABAQhAIGEEPDEC3rRp\nk9x7HiJEZQAAHQ1JREFU773OlpKmaqB3zZo1fiCm9XqPqagS0orAggULnD2+Zl6yqrkryhJsNDtr\nVnTha/t/Z848oZyF8M0CkUsIQCChBDwhgMeMGSPlypVzDCs0aNBAGjVqJN9//31CwVB54ggsWrRI\nhgwZ4piXrFGjRlBDbDLk8OGg6JARO3aEjCYSAhCAQMIJeGIK2jRcF6unczMnaN5sLtS5wuuuu07+\n/e9/JxwQDYgvgWXLlskTTzzhuBY8P4vRZpt2LlrUvfDdtu2Eu8H49oDaIAABCLgj4IkRsAlcMyuY\nEe68807p2LGjXH/99fLrr79mRPOZ4gTMrnO/fv0cAVy7du1Mvd29+8T67YEDmaJDXpQocULRqlSp\nkMlEQgACEPAEAU8I4Icfflhuu+02GT58uB9Kly5dpKV6Qe/cubM/jpPUJbB69Wrp06eP9O3bVy62\nDboBQWek1eqVyPHjAZEhTs/Q+RzTcrZp5yi+GULcTRQEIACB+BLwxBT0tddeKz/88IN6qsnsqmbA\ngAFylbqmsTRC6hKw595T9xHZUb9+/Uwd3b5d5NJLM0WFvNAVC3n//ZBJREIAAhDwJAFPCGAjU1jV\nVeuYv7gswZSzitvwx0U4cuSI2BEqWPzxaEOoUDcSl6sENm7cKN26dZNOnTqJKeBlDfr4o4aLLkL4\nRoVEBghAwHMEPCOAw5GxbUjmcH3SpEnhsvjjX331VbH8oYKtL2bdSxoqH3HxI7B582bpqoaYH330\nUUfzPWvNNWtGV7i68kqRzz/PeifXEIAABLxPwPMC+LHHHnNN8b777hM7QgVbS96yZUuoJOISQMCe\nha3zP/DAA2qVKrNNyJ07RUywRtuJZoY43n47AY2nSghAAAIxIOAJJazAfpjVqx1s3gxEknLnv6jb\nIvtB1Lp1a2nWrJm/f0uXil6fcI7w3Xf+6JAnKrcRviHJEAkBCCQLAU8I4MNqVcE0YCtVqqSmAvNJ\nyZIlnTVh24piRvgJqUNgu2pV2XrvrbfeKjfddJO/YzbitbXcDz7wR4U9sXXhF14Im0wCBCAAgaQg\n4IkpaNvza1OSs9S2YLVq1Rzhu1s3ftq6rX1ZHzx4UDp06JAUQGlkeAK7du1ynueNN97obDHLyFmr\nluizzriK/Hma/mRUy6UECEAAAklPwBMj4A8//NDxeFO3bl3HHnSePHkczWfTijUXdG+z0Jf0L9qe\nPXucaWdb7zVDKxlh6lT3wtfu2btX5PTTM+7mEwIQgEDyEvCEALap5rlz54akOFMt6ZcpUyZkGpHJ\nQWDfvn2OwlXDhg2lTZs2/kbrcr+0auW/jHhihrHMChYGNiJiIhECEEgiAp6Ygjb7z630m9icMlSv\nXl2KFSsmNl25YsUKxxWh2YomJCeBAyo1u3fvLvXq1ZN27dpl6sQll2S6DHlh5iQ//vjE+nDIDERC\nAAIQSFICnhDA9uVszhjmz58v69evd9aDbdRr677mGcmmpAnJR+DQoUPSo0cPOe+88+SRRx5xOuDz\niZhpSZPF6nchbLBHru6A1ThH2CwkQAACEEhqAp4QwEawQIEC0rhx46SGSeNPEjDN9t69ezvGT0yR\n7r77RGbMENHZaFeuBH/6SaRs2ZPlcQYBCEAg1Qh4RgCnGth07o/t5TYDKqV0/rhJk26O0lR2rIA2\nbYrwTef3h75DIF0IIIDT5UnHqZ/H1Glv//79na1kLVr0llq1sr98gFOFOD0sqoEABBJKAAGcUPyp\nVbk5uzCFutN0s27v3o9J/vzZV7Lftk30/tTiQm8gAAEIhCLAV10oKsRlm4BPtauGDh3qGE2pUuVx\nFb7Z26xbpYqI2l7RaetsV80NEIAABJKSAAI4KR+b9xo9cuRIMTOTn38+SB56KPrEiloc1e1mIn/4\ng6j2u6j2u0jRot7rFy2CAAQgkFsEon9T5lbNlJsyBGz/trkW3LdvuApTlaxRgnmWzLIlOModJEMA\nAhBIPQII4NR7pnHt0YQJE2TlytWybdtomT69QNS6dUcSwjcqJTJAAALpQAABnA5POZf6OEmHsv/9\n71J56aXRcvx4wai12HSzDpYJEIAABCCgBBDAvAanRGDEiJdlxIj58uuvT+v9RaKWYe4DzYcvAQIQ\ngAAEThBAAPMmZIuAbvNVi2WvqbLVHL1vnB6qSRUlmPWrQoWiZCIZAhCAQJoRQACn2QPPaXeLFJmu\nW43e1WLG61EianFm7xnhGxUTGSAAgTQkgABOw4d+Kl02JwqFCs1U4fuG3m7Ct2TEYooXP+FsoVKl\niNlIhAAEIJC2BNgHnLaPPnsdP+us91X4vqw32ZpvZP/MtsVo504RhG/2GJMbAhBILwII4PR63qfU\n2x49PtFtRqpFJabCHN5FkU01v/SSiO3zJUAAAhCAQGQCTEFH5pP2qZMmfSYjR05QDk/pUSEkj4K6\nA+nrr0Vq1lS1et6okIyIhAAEIJCVAF+XWYlw7Scwc+Z8NStpU84j9VBjzSFC1aoiP/wgkif7To9C\nlEYUBCAAgfQhgABOn2edrZ4uXLhQbrxxhN4zTI/qIe8toIav1q4NmUQkBCAAAQhEIcAacBRA6Zjc\nps0SqV9/kHb9ST10XjlMsP29BAhAAAIQODUCjIBPjVvK3lW69Hdq3WqA9s8E8IVh+/njj/jtDQuH\nBAhAAAIuCDACdgEpXbJUrLhKhW9f7W5/PeqG7ba5D6wQWh8r7D0kQAACEIBAZgII4Mw80vbqvvvW\nqEvB3tp/Oy4Ny6FHD5ErrgibTAIEIAABCLgkgAB2CSpVs731lkjFiut0/65KVumix+Vhu9qqlcjw\n4WGTSYAABCAAgWwQYA04G7BSKeuRIycsVW3dukm71V2PjnpcGbKLefOKvPqqyG23hUwmEgIQgAAE\nToEAAvgUoCX7LXv3qg8jdWLk8/2kXemqx0N6NA7ZrSLqaXDXLhSuQsIhEgIQgEAOCDAFnQN4yXir\neScqWtSE78/afJtybqPHtWG7sns3wjcsHBIgAAEI5IAAAjgH8JLp1qNHRS66SHWbHeXmbdr0Tnrc\nqUfzsN2wkTIWrsLiIQECEIBAjggwBZ0jfMlzc4kSIicMZ+zQRnfW45bfjuA+VK4ssnKliNl4JkAA\nAhCAQO4QYAScO1w9U+qzz4qYQD0hfHUx1xG+zfTz9qA2/v73J5wqbNiA8A2CQwQEIACBGBNgBBxj\noF4pbvt2214kcuBARot0PtlRuGqkn60zIv2fI9Tsc3dThiZAAAIQgEBcCCCA44I5vpUcPixSqlRg\nnfv1wrSd6+vRNjDBOZ8+XaRFi6BoIiAAAQhAIBcJIIBzEW4iirap5szC96A2w4a2tfVon6lJhQqd\nWOutVClTNBcQgAAEIBAHAgjgOECOVxWm6Wz7dk8GHQpLLz3MnaAZ2sgcdu4UMSMbBAhAAAIQiD8B\nlLDizzzXarwwk/MiNXUlffQoq4ft9z0ZihcX2bMH4XuSCGcQgAAE4k8AARx/5rlS47FjIqtXZxSt\nF9JPD5W00jMjUs49V2TuXBEb+WYeKfuzcAIBCEAAAnEigACOE+jcrGbUKJH8+TNqMOH7uB4WYSPg\nPHqcMD35/fciV1/tXPIHAhCAAAQSTIA14AQ/gJxWf/PNIjNmZJTi05Mn9dDFYBmsx+l6iJx3nsiK\nFc4pfyAAAQhAwCMEEMAeeRCn0gzbt3tS+FoJw/TQxV0ZoscJ4asnsmqV/SVAAAIQgICXCCCAvfQ0\nstEWs3BlU88ng11s1UMtashJ1WZb9yVAAAIQgID3CCCAvfdMorbI1nIffTQw2zi9UPuRjvDN508w\nRwpMPftxcAIBCEDAUwRQwvLU44jcGJ8u8f7rXyI1awbme04vlusxXI+T3hOaqblnM8px+smZ6MCb\nOIcABCAAgQQTYASc4AfgtnoTvlWqiGzaFHjHZL1YqMfTeqhZq9/CEd0CfAZPNgMHnxCAAAQ8SYCv\naU8+luBGmbnIzZsD46foxTw9bPq5iB4nghnYQPhm0OATAhCAgHcJIIC9+2z8LXvnnazCd5qmzdZj\nvB5mbONEeP99DGxksOATAhCAgNcJsAbs8SfUS00533JLYCNVGou6L3Kmnc/0J1xzjch11/kvOYEA\nBCAAAY8TYATs4QfUtavIU08FNnCWXkzVY6wepf0JTZqIzJnjv+QEAhCAAASSgAAC2KMPqVYt1W02\n5WZ/+EjP/qGHCd+z9TgRli1TR4O1M674hAAEIACBZCGAAPbYkzKXgmbX+fjxwIbN04uJethwuLwe\nJwLCN4MEnxCAAASSj4Dn1oCPqgTasWNH8pGMUYuDhe8XWrJpOo/Uo7K/lldeYeTrh8EJBCAAgSQk\n4AkBfPjwYenTp49U0r02+fLlk5IlS0rhwoV1arW2TJ5se13TI9hWo8wj3wXacTMxaeYlq/ohLNSt\nv3ff7b/kBAIQgAAEkpCAJ6agO3bsKFu2bJFZs2ZJtWrVHOG7e/duXQNdLp06dZKDBw9Khw4dkhCv\n+yabpvOPPwbmX6QX5lTBHCzUcBLMtOQGtThpgpoAAQhAAALJTcATI+APP/xQJk6cKHXr1lVH8UUk\nj0qa4sWLS4MGDWTs2LHy9ttvJzflKK3fv1/E9vqeDKpZJU/oMViP853oiy8WMQtXCF8HB38gAAEI\nJD0BTwhgm2qeO3duSJgzZ86UMmXKhExLhUizbqWz7QHBVJ/76WEC+IR68xVXiCxejF3nAEicQgAC\nEEh6Ap6Ygh44cKC0atVKxowZI9WrV5dixYrJrl271JPPCjGlrNmzzepT6oW9e0UqVgzsl7o5kr6/\nHTrk1VBITTzPn++c8gcCEIAABFKIgCcEcL169XSEt1gFzXxZv369sx5so15b923UqJEzJZ1CzP1d\nGTDAf6ona/VQs1fSQ4/6epwI+juEAAEIQAACKUjAEwLYuBYoUEAaN24chHjVqlWyXxdJTUhHC6/o\n3pzp081MY3BYunSpNGkySgYP/sqfuGhRffH5Tvrrq1dvgZx22skNuLmd/uabJmit/o16dFNTktfp\nGm8ePT/Rxp0766tjhZPtW7BggWpJn2xf/fr11d0g6QrMCfDh/eD/g++HjO+D7H4/ZtwXz888Pg3x\nrDC7dQ0ePFg1fzfIpEmTot5q2tJ2hArTpk3T/cVHpXXrk/t3Chcuminr/v17VSCfxJHb6a+8UlT6\n9t2sQrWTtuNhXeu+XH8AnKhf9c/U92/m9u3VOevA9hUtSnrgA4QP7wf/Hye/v/h+cPf92FVt/rZu\n3VouueSSwK+TuJx7XgDHisLrr7/uGPh4+OGHY1VkjsvZtGmL7nX+m+zefb+W1cwpz7YamZENfR8I\nEIAABCCQywQSKYA9MwWdwdiUrvaoU9szzzwzIyolP3/55Rfp3r2TCtu7ddTezNFyLldOpH37E6Yo\nU7LTdAoCEIAABPwEPLENKd0sYW3fvt0xMHLbbbfJTTfdKLffLjJ0qMhf/4rw9b+ZnEAAAhBIcQKe\nGAGnkyUs215l1r1uvPFGadmyZYq/XnQPAhCAAATCEfDECDhdLGHZ1Hrnzp1VG7uJ3HnnneGeCfEQ\ngAAEIJAGBDwhgNPBEta+ffukS5cu0rBhQ2nTpk0avFp0EQIQgAAEIhHwxBR0qlvCOnDggCpcdXf2\nMrdr1y7S8yANAhCAAATShIAnBHAqW8I6dOiQ9OjRQ8477zx55JFH0uS1opsQgAAEIBCNgCcEsDUy\nnCWsaB3wcrppd/fu3VsqV67sKF55ua20DQIQgAAE4kvAE2vA8e1yfGqz/cyPPfaYlC5dWrp16xaf\nSqkFAhCAAASShgACOBce1bFjx6R///7qZrCw9OrVK2WdSeQCOoqEAAQgkDYE0sYU5ZIlS6R58+au\nnDrk5OmbLdpt27Y5NpvNo1Mesy3pInz88cfqelB9D6ZZMNvdefPmzeRUIh0Q2I80W6IoWLBgOnQ3\nUx9tR4D9OE23YM/bvg/sfU+nYN+J9r7/4Q9/8GS3165dKx999JFUqFAh7u1LGwEcd7LZrPDqq6+W\nefPmZfOu5M9u+6JtW9bFF5/wf5z8PXLXg2XLljkORsaNG+fuhhTKla7v+rPPPqsOV8qIWcBLp2C7\nQMzoUKr6dc/Js2QKOif0uBcCEIAABCBwigQQwKcIjtsgAAEIQAACOSGAAM4JPe6FAAQgAAEInCIB\nBPApguM2CEAAAhCAQE4IIIBzQo97IQABCEAAAqdIAAF8iuC4DQIQgAAEIJATAmxDygm9GN77v//9\nT8qVKxfDEpOjqF9//VWKFi0q+fLlS44Gx6iVtifU3FOWKlUqRiUmTzHp+q7v3r3b2e+ebnugbR/w\n1q1bpWzZssnzksappQjgOIGmGghAAAIQgEAgAaagA2lwDgEIQAACEIgTAQRwnEBTDQQgAAEIQCCQ\nAAI4kAbnEIAABCAAgTgRQADHCTTVQAACEIAABAIJIIADaXAOAQhAAAIQiBMBBHCcQFMNBCAAAQhA\nIJAAAjiQBucQgAAEIACBOBFAAMcJNNVAAAIQSAcCR48eFTO+QYhOAAEcnVHMcsybN0+uvPJKqVq1\nqrRo0UJ27NgRsmy3+ULe7MFI6+ftt98u5557rtSpU0e+/PLLkK1cvny53HXXXXLRRRdJkyZN5PXX\nXw+ZL5kihw4dKnXr1nWeuZ1HCw899JC0b98+WjbPp7t9h//zn//IZZddJrVq1ZIbbrhBVqxY4fm+\nRWqg23f9xx9/lHvuuUcuvvhiuf766+XTTz+NVGzSpG3atEmqVKkia9euDdtmt+9G2AJSKUF/qRDi\nQOCXX37xqalJ3zfffONTM4S+zp07++6///6gmt3mC7rRwxG33Xabb9CgQb7jx4/75s6d6zv77LN9\n+/fvD2px06ZNfS+99JITv3nzZt9ZZ53l27JlS1C+ZIl44403fL///e99O3fu9Kn5RZ/+sPDNnj07\nbPNnzpzpK1mypE+FcNg8yZDg9h0+ePCgr1q1ar758+c73Zo6daqvZcuWydDFsG10+663a9fON2TI\nEKecr7/+2uFw5MiRsOUmQ8ILL7zgq169ui9v3ry+NWvWhGyy23cj5M0pGGlTBYQ4EHjvvfd811xz\njb8m/YXoK168uP8648Rtvoz8yfCptp59avPZ39RLL73U9+GHH/qv7eTYsWO+f/3rX86Pk4wE+2eO\nJLAy8nn1s23btr6///3v/uYNGzbM9+CDD/qvA0+2bdvmu+KKK3wDBgxIegHs9h2ePn2679prr3Uw\n2I+UVAhu3nXr58033+wbNWqU0+WVK1f6ChUq5LMfJMkaDh065DzLVatW+cqUKRNWALt9N5KVQ3bb\nzRR0nKYzNm7cmMnZgo4CZdeuXaIvbqYWuM2X6SYPX9iUnPVRR3b+VppR9p9//tl/bSennXaa3HLL\nLaK/np34jz/+2Jmib9CgQaZ8yXSR9Vlav80ofajQoUMHefzxx6VIkSKhkpMqLmu/w73rGzZscN6L\nRo0aiX5pi/7gku+++y6p+hrYWLfvut0zePBgmTRpktx6662iMz/y7LPPSv78+QOLS6pzc6bywQcf\nyHnnnRex3W7fjYiFpFAiAjhOD9O8/gR6QSlYsKBTs07FZmqB23yZbvLwRdb+WFOt73v37g3b6u+/\n/95ZH3vmmWekRIkSYfN5PSFr33WUI/v27Qtq9pQpUxwm1113XVBaMkZk7Xe4d12nI0Wn6eXhhx8W\nu6dZs2YyfPjwZOyy0+as/bbIcO/6F1984SgqXXDBBVKhQgWxdVFTXkr1kJVRuHcj1Tlk9A8BnEEi\nlz9Lly4t5o4sI5grugIFCsiZZ56ZEeV8us2X6SYPX2TtjzXVOJQvXz5kq3U6Tq6++mrp37+/o5AV\nMlOSRGbte6h+2xdSp06dRJcnRNeAHSUkGxnqumiS9DK4mVn7He5dtx9XpnzVqlUrxyVlz549ZcaM\nGWKuGpMxZO239SHUMzdB26NHD3nttddEdSPEhLEuyci///3vZOx2ttqclVG4dyNbhSZx5jOSuO1J\n1fSKFSvK+vXr/W2280qVKvmvM07c5svI7/VP+5K1X7mm9Wl9s2B9r1y5clDTTXPyj3/8o/Tt29cZ\nFQVlSLII668J04wQ6pnr2qfUqFFDJk6c6GT76aefRNcC5ZVXXpFknX53+w5bPvMFnRFs+eHAgQOi\nynoZUUn16fZdt6lq1XlwtP2tg7b8onoRsm7dOufHZ1J1OpuNdftuZLPY5M2e3UVj8p8aAVOwMK3e\nOXPmOMoW9957r69Xr15OYdu3b/d99NFHznmkfKdWc+LvMmWkjh07+kzL88033/Sdf/75fmUr04o2\nzUgLDRs29OnIwFHY0pGh82nKHckaTOFEtyD5TKNbv1x9Kmh9pvFqQbdc+ZYtWxbUtZEjRya9Elak\ndzjwXdfRoa9UqVI+3YrkcNB1UZ+uBwcxSaYIt++6/tD06QjY6ZopZOpMmM/YpELIqoQV+MwjvRup\n0Pfs9gEt6OwSy0F+25aiSjY+XfPxNW7c2KfTL05pugfQ0YLMKDpcvoz0ZPs04VO7dm2fTjs72xRM\n6GYE25I0a9Ys34IFC2znftDxj3/8IyNr0n3ativbaqYjI58qYDkazhmd+Mtf/uK77777Mi79n6kg\ngK0z4d7hrO/6O++847AxjfeaNWv6TBglc3Dzrlv/7H3/05/+5PxAsx9ptoUnVUJWAZz1mYd7N1Kl\n/9npRx7LnLzj9+Rrua3/2LpH1rXfrD1xmy/rfV6+NqUb03ZNt2DrgKbhmsxarqfyzNy+w/YVpNuw\nUurdcPuu27tRrFixU8Gb1Pe4fTeSupMuGo8AdgGJLBCAAAQgAIFYE0ALOtZEKQ8CEIAABCDgggAC\n2AUkskAAAhCAAARiTQABHGuilAcBCEAAAhBwQQAB7AISWSAAAQhAAAKxJoAAjjVRyoMABCAAAQi4\nIIAAdgGJLBCAAAQgAIFYE0AAx5oo5UEAAhCAAARcEEAAu4BEFghAAAIQgECsCSCAY02U8iAAAQhA\nAAIuCCCAXUAiCwQgAAEIQCDWBBDAsSZKeRCAAAQgAAEXBBDALiCRBQIQgAAEIBBrAgjgWBOlPAhA\nAAIQgIALAghgF5DIAgEIQAACEIg1AQRwrIlSHgQgAAEIQMAFAQSwC0hkgQAEIAABCMSaAAI41kQp\nDwIQgAAEIOCCAALYBSSyQAACEIAABGJNAAEca6KUBwEIQAACEHBBAAHsAhJZIAABCEAAArEmgACO\nNVHKg4BHCaxbt07q1Kkja9ascVo4efJkue2228Tn83m0xTQLAqlNII/+8/Hfl9rPmN5BwE+gS5cu\nsnr1apk4caLUrVtX3n//fbnsssv86ZxAAALxI4AAjh9raoJAwgns27dPatWqJcWKFZPmzZvL0KFD\nE94mGgCBdCXAFHS6Pnn6nZYEChcuLB06dJBvv/1W/vKXv6QlAzoNAa8QYATslSdBOyAQBwI7d+6U\nCy+80DnKlSsnr7zyShxqpQoIQCAUAUbAoagQB4EUJdC1a1dp1qyZvPXWWzJnzhxnDThFu0q3IOB5\nAmd4voU0EAIQiAmBTz75RN555x1ZuXKlFC9eXEaPHi0PP/ywMx1dpEiRmNRBIRCAgHsCTEG7Z0VO\nCEAAAhCAQMwIMAUdM5QUBAEIQAACEHBPAAHsnhU5IQABCEAAAjEjgACOGUoKggAEIAABCLgngAB2\nz4qcEIAABCAAgZgRQADHDCUFQQACEIAABNwTQAC7Z0VOCEAAAhCAQMwIIIBjhpKCIAABCEAAAu4J\nIIDdsyInBCAAAQhAIGYEEMAxQ0lBEIAABCAAAfcEEMDuWZETAhCAAAQgEDMCCOCYoaQgCEAAAhCA\ngHsCCGD3rMgJAQhAAAIQiBkBBHDMUFIQBCAAAQhAwD0BBLB7VuSEAAQgAAEIxIwAAjhmKCkIAhCA\nAAQg4J4AAtg9K3JCAAIQgAAEYkYAARwzlBQEAQhAAAIQcE/g/wFDa2VW/oFOegAAAABJRU5ErkJg\ngg==\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -i P\n",
    "plot(ecdf(P), col='blue') \n",
    "abline(0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the kind of distribution function whose quantiles we find..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10c622710>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGQdJREFUeJzt3XuUnHV9x/H3J5sASRFRIxdDhFONgAgoSghgYBGoASvB\nS0uDIAoe8RKMKJoETyU9VsUbIKKQQsRLlSiCNEWEKmXLpYpQCBdJMIGiSRAoIHKLQODbP37P6mTZ\n3ZndnZnfPM98XufM2X12np39TEK+/PZ3VURgZmblNS53ADMzGxsXcjOzknMhNzMrORdyM7OScyE3\nMys5F3Izs5KrW8glfVPS/ZJuG+aeMyWtknSLpNc1N6KZmQ2nkRb5+cCsoZ6UdCjwyoiYBrwfOLtJ\n2czMrAF1C3lEXAP8YZhbDgO+Xdx7PbClpK2bE8/MzOppRh/5FGBNzfVaYLsmvK6ZmTWgWYOdGnDt\ndf9mZm0yvgmvsQ6YWnO9XfG1jUhycTczG4WIGNhY3kgzCvkyYC6wVNIM4JGIuH80YcpM0qKIWJQ7\nR6tU+f1V+b2B318nkZgEHAK8HdiR1Ah+EfAYQ/Zk6CX1XrduIZd0AbA/MFnSGuAUYAJARCyOiMsk\nHSppNfAE8N4G3o+ZWaVJvAzYmdRLsR2wK2kG4I3AhcBXSeOLD0Tw7DCvU7c3o24hj4g5Ddwzt949\nZmZVJrE9cDiwLzADmATcTirWa4D/AD4SwQPN/tnN6FqxpC93gBbryx2ghfpyB2ixvtwBWqwv1w+W\nmAy8D3gnsD2pq3kZcDJwV0R7Jn6oXQdLSIoq95GbWfeQ2Bw4EZgHXAJ8H7g6gg3N/1n1a6db5GZm\nDZDYBNgLOIi0iv0/gb0iuCtrMFzIzcyGJbE7sAg4EFgFXAnMiuCWnLlquZCbmQ1CYmvgn0nbkPwz\ncFwED+dNNThvY2tmVkNiosRC4NfAo8COEXytU4s4uEVuZgaAxDhgDvA50lzvGRGszpuqMS7kZtbV\nJF4EHE0awFwPHBXBNXlTjYwLuZl1FYm/AnYH9gD2AQ4FLiNtNfJf7Zr73UyeR25mXaFYMv954O9I\n/d83kbpQfhzBgzmzDcfzyM2s60lsBnyseJwLbBPBo3lTNZcLuZlVlsSrgB8DvwGmR3B35kgt4emH\nZlZJEm8FrgW+GsHbqlrEwS1yM6uYoivlZOBYYHYEv8gcqeXcIjezSpDYVOLDwGrSrJQ3dEMRBxdy\nMys5iXES7yHtg3IIcHgEsyO4L2+y9nHXipmVlsQM4EzgOeDvI/hl5khZuJCbWelIjAe+RtrQagHw\nvQiey5sqHxdyMysViYnAD4BNgJ0ieCxzpOzcR25mpSGxJXAF8DhwmIt44kJuZqUgsQXpfM7lpI2t\nns6bqHN4rxUz63jFFrMXAfcDHyzjxlaj5b1WzKwq5gPbAv/QTUW8US7kZtbRJA4GTiDtlfJU7jyd\nyIXczDqWxA7Ad0kt8bWZ43QsD3aaWUcqDj/+D+CzEfRljtPRXMjNrOMU0wwvBy6I4Gu583Q6z1ox\ns44iMYk0V/xmYF63D242UjvdIjezjiHxUuBS4G7go91exBvlQm5mHUFiX9I5mtcDx3Xz3ikj5Vkr\nZpaVhEjnaX6SVMAvzRypdFzIzSwbiZcA3wK2AvaK4J6sgUrKXStmloXEG0kDmncCM13ER88tcjNr\nO4kjgdNIXSk/yZ2n7Dz90MzaSuKvSQOaB0Zwa+48nc7TD82so0j0AN8BTnURb566hVzSLEkrJa2S\nNH+Q5ydLulzSckm3S3pPS5KaWRXMB54GTs8dpEqG7VqR1EMaiDgIWAfcAMyJiBU19ywCNo2IhZIm\nF/dvHREbBryWu1bMupjE64GfAq+PYE3uPGXRjK6V6cDqiLgnIp4BlgKzB9zze2CL4vMtgIcGFnEz\n624SO5MOhpjnIt589Qr5FNjoD31t8bVa5wK7SLoXuAWY17x4ZlZ2EvuTjmg7JYILMseppHrTDxuZ\n0nIysDwieiW9AviZpN0j4nmHohbdMP36IqKv4aRmVjrFNMMzgDkRXJk7TxlI6gV6R/I99Qr5OmBq\nzfVUeN7m7vsAnwWIiLsk/S+wI3DjwBeLiEUjCWdm5VW0xL9MmmZ4W+48ZVE0cPv6ryWdUu976nWt\n3AhMk7SDpE2AI4BlA+5ZSRoMRdLWpCJ+d8OpzaxyJCYA3wBOcBFvvWFb5BGxQdJc0t7APcCSiFgh\n6fji+cXA54DzJd1C+h/DJyPi4RbnNrPONo80vnZx7iDdwCs7zaypJLYDlgN7R7Aqd56y88pOM8vh\nNOAbLuLt402zzKxpJA4G3gAckztLN3GL3Myaojhr8xxgbgTrc+fpJu4jN7OmkPgysG0E78qdpUoa\nqZ3uWjGzMZN4A3AUsGvuLN3IXStmNibFnPHzgJMi+L/cebqRC7mZjdVJwH3A93IH6VbuWjGzUZPY\nD/goMD2iob2ZrAXcIjezUZGYStra+t0R/DZ3nm7mQm5mIyYxEfgxcHoEV+TO0+08/dDMRkRCwLeA\nTYAj3aXSWp5+aGatcCSwBzDDRbwzuEVuZg2TmAzcDrw1ghty5+kGjdROF3Iza5jEd4CHIjgxd5Zu\n4a4VM2saib8B9gNekzuLbcyzVsysLom/Im2I9YEIHs+dxzbmrhUzq0viNGBrb4jVfu5aMbMxkziA\ndF7vbrmz2ODctWJmQ5LYAjgfeF8ED+XOY4Nz14qZDUnifOCpCD6QO0u3cteKmY2axOHATOC1ubPY\n8NwiN7PnkdgUuBs4IoJrc+fpZo3UTveRm9lgjgJudREvB3etmNlGJMYBnwA+mDuLNcYtcjMb6K3A\n40Bf5hzWIBdyMxvok8AXvbNhebiQm9mfSewLbANcnDuLNc6F3MxqfQI4LYINuYNY4zz90MwAkNgN\n+DmwQwRP5s5jiacfmllDipkqZwOfdhEvHxdyMwM4jlQP/iV3EBs5d62YdTmJrUjHtx0Uwa2589jG\nfNSbmdVVHN92fwSfyJ3Fns+bZpnZsCTeRDq+bZfcWWz03Edu1qWK49vOBeZG8ETuPDZ6dQu5pFmS\nVkpaJWn+EPf0SrpZ0u2S+pqe0sxa4XPAf0dwae4gNjbD9pFL6gHuBA4C1gE3AHMiYkXNPVsC1wFv\njoi1kiZHxIODvJb7yM06hMRMYCmwawQP585jQ2vGPPLpwOqIuCciniH9xc8ecM+RwEURsRZgsCJu\nZp1DYhLwTeDDLuLVUK+QTwHW1FyvLb5WaxrwYklXSbpR0tHNDGhmTfcZ4IYILskdxJqj3qyVRuYm\nTgD2AA4EJgG/kPTLiFg11nBm1lwSuwBHA6/OncWap14hXwdMrbmeSmqV11oDPBgR64H1kq4Gdgee\nV8glLaq57IuIvpEGNrMx+SLwuQjcBdqhJPUCvSP6njqDneNJg50HAvcCv+L5g507AWcBbwY2Ba4H\njoiIOwa8lgc7zTKSOAg4B3h1BE/nzmONGfOCoIjYIGkucAXQAyyJiBWSji+eXxwRKyVdDtwKPAec\nO7CIm1leEj3AV4D5LuLV4yX6Zl1A4r2kjbFm+uSfcvFeK2bWv4LzTuAdEVyfO4+NjPcjNzOAU4A+\nF/Hq8qZZZhUm8TrgGGDX3FmsddwiN6soifGkTbEWRPBA7jzWOi7kZtV1AvAo8K3MOazFPNhpVkES\nOwA3AntHPH9xnpWHBzvNuteXgNNdxLuDW+RmFVPsp/KfwF/7wIjyc4vcrDudDJzhIt493CI3qxCJ\nacB/A6+I4NHceWzs3CI36z4nA2e5iHcXLwgyq4hipsphwCszR7E2c4vcrDoWAOdE8IfcQay93CI3\nqwCJVwHvBHbKncXazy1ys2o4HTjVJ/90J7fIzUpO4i2kfvG35c5iebiQm5WYxKak1vg8n/zTvdy1\nYlZu84A7I/hp7iCWjxcEmZWUxLbAbXhjrErzgiCzavsicJ6LuLmP3KyEJPYD9gdenTuL5ecWuVnJ\nSEwAvg58LILHc+ex/FzIzcrnw8B9wEW5g1hn8GCnWYkUA5y3AjMjWJk7j7VeI7XThdysRCSWAA9G\nMD93FmuPRmqnBzvNSkLi5aTVm6/IncU6i/vIzcrj48AS725oA7lrxawEJF4K3Am8JoJ7c+ex9vGC\nILPq+AhwoYu4DcYtcrMOJ7EFcDewVwR35c5j7eUWuVk1HA/8zEXchuIWuVkHk5gI3AUcEsEtufNY\n+7lFblZ+HwJ+6SJuw3GL3KxDFX3jq4A3RfDr3HksD7fIzcrtROAKF3Grxy1ysw4kMRlYCUyP4O7c\neSyfprTIJc2StFLSKklD7u8gaU9JGyS9fTRhzWwj84EfuohbI4bda0VSD3AWcBCwDrhB0rKIWDHI\nfV8ALgfc6jYbA4kpwLHArrmzWDnUa5FPB1ZHxD0R8QywFJg9yH0nAD8C/q/J+cy60ReAxV7FaY2q\nt/vhFGBNzfVaYK/aGyRNIRX3NwF7Au3pdDerIIleYCY+ws1GoF6LvJGifAawINKoqXDXitmo1Bzh\ndmIET+TOY+VRr0W+Dphacz2V1Cqv9XpgqSSAycAhkp6JiGUDX0zSoprLvojoG2lgswr7KPA74Me5\ng1g+knqB3hF9z3DTDyWNJ22deSBwL/ArYM7Awc6a+88H/j0iLh7kOU8/NBuCxHbAcmBGBKtz57HO\nMeYTgiJig6S5wBVAD7AkIlZIOr54fnHT0pp1t68AX3cRt9HwgiCzzCT2AX4A7BjBk7nzWGfxEn2z\nDich4DTgZBdxGy0XcrO8jiB1cX4vdxArL3etmGVS7DW+Anh3BFfnzmOdyV0rZp1tHnCTi7iNlVvk\nZhlIbAXcAewdwarceaxzNVI7XcjNMpD4BvBUBCfmzmKdbczzyM2s+SR2Bt4J7JQ7i1WD+8jN2u9L\nwOcjeDh3EKsGt8jN2kjiQGBn4B25s1h1uEVu1iYSPcCXgfkRPJU7j1WHC7lZ+xwDPAlclDuIVYtn\nrZi1gcSLSdMND43gptx5rDw8/dCsQ0icA2yIYG7uLFYunn5o1gEk9iQdh7hz7ixWTe4jN2uhYoDz\nbGBBBI/kzmPV5EJu1lrvB9YD38kdxKrLfeRmLSKxE3ANcEAEt+fOY+Xk3Q/NMpGYBFwILHQRt1Zz\ni9ysBSS+CUwg7TXenn9kVkmetWKWgcQxwN7Ani7i1g5ukZs1kcSrgOtwv7g1ifvIzdqomGr4beCf\nXMStnVzIzZrnJNJUw2/kDmLdxV0rZk0g8RrgKlK/+D2Z41iFuGvFrA0kJpAW/Cx0EbccXMjNxm4B\ncD+wJHcQ607uWjEbg+L8zauBPSJYkzuPVY+7VsxaSGIccC6wyEXccnIhNxu9D5D+DZ2dO4h1N3et\nmI2CxFTgJmD/CO7Inceqy10rZi0gMR44D/iai7h1AhdysxGQEPD14vLzObOY9fOmWWYjsxCYDuwX\nwTO5w5iBC7lZwyTeTTrxZ+8IHsudx6yfBzvNGiCxN/BvpMHNFbnzWPfwYKdZE0hsCXwfeL+LuHWi\nhgq5pFmSVkpaJWn+IM+/S9Itkm6VdJ2k3Zof1az9isHNc4DLIrgkdx6zwdTtI5fUA5wFHASsA26Q\ntCwialsmdwP7RcQfJc0C/gWY0YrAZm32XmAX0gCnWUdqZLBzOrA6Iu4BkLQUmA1/+RUzIn5Rc//1\nwHZNzGiWhcROwBeA3gjW585jNpRGulamwEb7SKwtvjaU44DLxhLKLDeJF5EGNxdE8OvcecyG00iL\nvOFpLZIOAI4F9h3i+UU1l30R0dfoa5u1S7G/+IWkfnFvTWttJakX6B3J9zRSyNcBU2uup5Ja5QN/\n+G6kneBmRcQfBnuhiFg0knBm7VYMbn4VeIp0dJtZWxUN3L7+a0mn1PueRgr5jcA0STsA9wJHAHNq\nb5D0cuBi4KiIWN1oYLMO9GFgf9Kin2dzhzFrRN1CHhEbJM0FrgB6gCURsULS8cXzi4FPAy8CzpYE\n8ExEeJTfSkXiHaQl+DMjeDR3HrNGeWWnGSBxAPAD4M0R3Jw7j1k/r+w0a4DE60hF/AgXcSsjF3Lr\nahLTgJ8AH4jgqtx5zEbDhdy6lsT2wM+Bf4zg4tx5zEbLhdy6ksS2wJXAVzxX3MrOhdy6jsRkUkt8\nSQRn5s5jNlaetWJdRWIicBVwVQQLc+cxq6eR2ulCbl1DYhzwQ+BPwNERjW8/YZZLI7XTR71ZNzkV\n2Ao42EXcqsSF3LqCxPGk7Zf3ieCp3HnMmsmF3Cqt2ATr48A84IAIHsocyazpXMitsiTGk3YynElq\nia+p8y1mpeRCbpUk8QLge8BmpE2w/pg5klnLeB65VY7EnsDNpL303+IiblXnFrlVRjG98KTi8aEI\nfpQ5kllbuJBbJUi8FPgusDmwZwS/zRzJrG3ctWKlJ7EvcBOwnHTivYu4dRW3yK20iqmFJwLzgeMi\nuDRzJLMsXMitlCQmAecBOwJ7RXBP3kRm+bhrxUqn2Ef8OuBZ4I0u4tbtXMitVCTeCvwS+Dbw7gjW\nZ45klp27VqwUJF5MWqW5D+lszaszRzLrGG6RW0eT6JE4ErgNeAjYzUXcbGNukVtHktgUOJo0I+V+\nUiv82rypzDqTC7l1FInNgOOABcDtwLERXJM3lVlncyG3jjCggC8H3hbBjXlTmZWDC7llVSyt/xDw\nQeBXuICbjZgHOy0LiZ0lFgO/AaaQDn04zEXcbOTcIre2KZbUHwx8DHgtcDawUwT3Zw1mVnIu5NZy\nxUk97yT1f/cApwGHR/CnrMHMKsKF3FpGYmvgXaQ+8N8DnwIu8wn2Zs3lQm5NJTEROAQ4BtgfuAQ4\nJoLrsgYzqzAXchuzYurgm4G/B95C2hv8X4GjIngsZzazbqCI9vyWKykiQm35YdZyEpsABwJHAIcB\ntwI/AC724KVZ8zRSO13IrWFFn/cs4G9Js0/uIBXvCyO4N2c2s6pyIbdRk+gBpgF7AjOB/YBtgJ8D\nlwI/dcvbrPWaUsglzQLOIE0bOy8ivjDIPWeSBrieBN4TETePJozlIbE5sAuwe/F4LbAbabOqm4Br\nisdtETybK6dZN2qkdg472CmpBzgLOAhYB9wgaVlErKi551DglRExTdJepEUeM8acvmQk9UZEX+4c\ngyn6s6cA2wMvLz5uT2pxvxJ4MbACuKV4XAgsj+CRv7yGeiOikkW8k//umsHvr/rqzVqZDqyOiHsA\nJC0FZpP+0fc7jHRaCxFxvaQtJW0dEd32a3cv0NfOH1gU6BcCWwHbFo+XkYp2/2Mq8FLgPuB3wG+L\nxw3ABcAqYG0Ez9X5cb20+f21US/VfW/g91d59Qr5FGBNzfVaYK8G7tkOurv/tFiO3gNMKB6b1Dwm\nApNqHhOBzYqPm9c8XgBsUXzs//yFxccti9d6BHgQuJe06Ob3wN2krpB7ScX79xFsaPFbNrNM6hXy\nRkdCB/bfDPp9EpcOcf9Qr6NBPq/9OPBrQ73mwO/REF8fN+C5cTWPngHX/V/rAcbDwhdKnMBfivUE\n0p/DM8Xj6ZrHetJ4whPF5+uBPxUfH695PAA8CjxWPP5Yc/0H4AmvkjSzYQc7Jc0AFkXErOJ6IfBc\n7YCnpHOAvohYWlyvBPYf2LUiyQXHzGwUxjTYCdwITJO0A+nX9COAOQPuWQbMBZYWhf+RwfrHPWPF\nzKw1hi3kEbFB0lzgClIXwpKIWCHp+OL5xRFxmaRDJa0mdRW8t+Wpzczsz9q2IMjMzFqjrScESfqM\npFskLZd0paSp7fz5rSTpS5JWFO/vYkkvzJ2pmST9naRfS3pW0h658zSLpFmSVkpaJWl+7jzNJOmb\nku6XdFvuLK0gaaqkq4r/Lm+X9JHcmZpF0maSri9q5R2SPj/s/e1skUt6QUQ8Vnx+ArB7RLyvbQFa\nSNLBwJUR8ZykUwEiYkHmWE0jaSfgOWAx8PGIuClzpDErFrzdSc2CN2BO7YK3MpM0kzT76TsRsWvu\nPM0maRtgm4hYLmlz4H+Awyv09zcpIp6UNB64FjgpIq4d7N62tsj7i3hhc9L850qIiJ9FRP+imutJ\nc+krIyJWRsRvcudosj8veIuIZ4D+BW+VEBHXkKapVlJE3BcRy4vPHyctVHxZ3lTNExFPFp9uQhqj\nfHioe9t++LKkz0r6HenggVPb/fPb5FjgstwhrK7BFrNNyZTFxqCYWfc6UiOqEiSNk7SctLjyqoi4\nY6h7m36whKSfkXbJG+jkiPj3iPgU8ClJC4DTKdEsl3rvrbjnU8DTEfH9toZrgkbeX8V4pL8Cim6V\nHwHzipZ5JRS/4b+2GG+7Yrg9ZZpeyCPi4AZv/T4la7XWe2+S3gMcSjpwoXRG8HdXFetIe9H0m0pq\nlVtJSJoAXAT8a0RckjtPK0TEHyX9BHgDQ+wp0+5ZK9NqLmcDz9vutqyK7X4/AcyOiKqfDl+VxV1/\nXvAmaRPSgrdlmTNZgyQJWALcERFn5M7TTJImS9qy+Hwi6SCXIetlu2et/AjYEXgWuAv4YEQ80LYA\nLSRpFWlQon9A4hcR8aGMkZpK0tuAM4HJpD1fbo6IQ/KmGjtJh/CX/faXRMSw07zKRNIFpAOwX0La\nt+fTEXF+3lTNI+mNwNWkYwb7C9nCiLg8X6rmkLQraVfZ/n2dvhsRXxryfi8IMjMrt7bPWjEzs+Zy\nITczKzkXcjOzknMhNzMrORdyM7OScyE3Mys5F3Izs5JzITczK7n/B54aFMnpU3y9AAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c2acd10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "V = np.linspace(-3,3,101)\n",
    "Z = np.random.standard_normal(5000)\n",
    "plt.plot(V, [pivot(v, truth, Z) for v in V])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}