fperez/Randomization Tests.ipynb

## Randomization Tests.ipynb
{
 "metadata": {
  "name": "",
  "signature": "sha256:05ad055d743b0d9da983a55e70f961ce2b53de55b912d63f7b4ef698c5d5a4a4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Randomization tests"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Author: Douglas Bates**\n",
      "\n",
      "**Date: 04/03/2014**\n",
      "\n",
      "## Note by F. Perez\n",
      "\n",
      "This is a port as an IPython notebook of [a post by Douglas Bates on Rpubs](http://rpubs.com/dmbates/15250). Aside from this note, the entire code and text are Doug's authorship, I just turned his Rmd sources into and IPython notebook.  To run this notebook, you will need the new [R kernel for IPython](https://github.com/takluyver/IRkernel)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Purpose\n",
      "\n",
      "In the interests of full disclosure, the purpose of this contribution to [Rpubs](http://rpubs.com) is to present code for a couple of simple tests as implemented in R so they can be compared to [Julia](http://julialang.org) code that I will publish in a gist.\n",
      "\n",
      "# A randomization test in paired designs\n",
      "\n",
      "One of the easiest types of statistical tests to explain is a [randomization test](http://en.wikipedia.org/wiki/Randomization_test) comparing two samples of responses.  In a paired design we reduce the two, equally-sized sets of responses to their differences.  That is, if we write the two samples as $(x_i,y_i), i=1,\\dots,n$ then we consider the differences $d_i = y_i - x_i, i=1,\\dots,n$.\n",
      "\n",
      "Under the hypothesis $H_0:\\mu_X=\\mu_Y$ that the means of the two populations are equal the signs of the individual differences are arbitrary.  Under the alternative, say $H_a:\\mu_X < \\mu_Y$ we expect the differences to be, on average, positive.\n",
      "\n",
      "We can compare the sum (or, equivalently, the average) of the observed differences to the set of all possible differences that we could generate under $H_0$ by flipping the signs of some of the differences.  If we have $k$ differences there will be $2^k$ possible sums corresponding to changes of sign on some subset of the differences.\n",
      "\n",
      "## Sample data from a paired design\n",
      "\n",
      "The [Secchi Disk](en.wikipedia.org/wiki/Secchi_disk) is used to measure water transparency in lakes or oceans.  The measurements are conducted by lowering the disk into the water and recording the depth at which the disk is no longer visible.  Such measurements were conducted in 1980 on 22 Wisconsin lakes and in 1990 on the same lakes.  If the water clarity (or, conversely, the turbidity) remained constant than the differences in the Secchi depths would be neither systematically positive nor negative."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Secchi <- within(\n",
      "    data.frame(d80=c(2.11,1.79,2.71,1.89,1.69,1.71,2.01,1.36,2.08,1.10,1.29,\n",
      "                     2.11,2.47,1.67,1.78,1.68,1.47,1.67,2.31,1.76,1.58,2.55),\n",
      "               d90=c(3.67,1.72,3.46,2.60,2.03,2.10,3.01,1.82,2.64,2.23,1.39,\n",
      "                     2.08,2.92,1.90,2.44,2.23,2.43,1.91,3.06,2.26,1.48,2.35)),\n",
      "    diff <- d90 - d80)\n",
      "str(Secchi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "'data.frame':\t22 obs. of  3 variables:\n",
        " $ d80 : num  2.11 1.79 2.71 1.89 1.69 1.71 2.01 1.36 2.08 1.1 ...\n",
        " $ d90 : num  3.67 1.72 3.46 2.6 2.03 2.1 3.01 1.82 2.64 2.23 ...\n",
        " $ diff: num  1.56 -0.07 0.75 0.71 0.34 0.39 1 0.46 0.56 1.13 ...\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(Secchi$diff)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "[1] 10.94"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As we are primarily interested in the differences we check the distribution of the differences with an empirical density plot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "library(lattice)\n",
      "densityplot(~diff, Secchi, \n",
      "       xlab=\"Difference between 1990 and 1980 Secchi depths on 22 Wisconsin lakes\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": []
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3deVxU5eIG8GfYF9nBhUU2FbcUlxIMUXBN1FDUMsu0rFtXy26ZbXrrprdb11v3l5W3xczENnfFPQUXXEE0FUEUBEVUFkFQFBk4vz/GJoIZ1pl5zzDP9w8/8DpzzjOH4fH4zlkUkiSBiIjkx0x0ACIi0owFTUQkUyxoIiKZYkETEckUC5qISKZY0EREMsWCJiKSKRY0EZFMsaCJiGSKBU1EJFMsaCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZIoFTUQkUyxoIiKZYkG3iFKpFB2hOYw0dlVVlSRJolM0h5FucMYWjgXdfJWVlZMmTRKdojliYmKqqqpEp2iyL7/8cseOHaJTNFl2dvbf/vY30SmaIzo6WnSE5njnnXdSU1NFp9ANFnTzVVdX37lzR3SK5igvLzfGXdF79+7du3dPdIomM9LYAG7fvi06QnNUVlYa6QaviwVNRCRTLGgiIpliQRMRyZSF6ACaXblyZevWraJTNECpVObm5n799deigzTZlStXli1bZmZmZP88Hz58+MKFC9evXxcdpGmuX7+elpZmjO+TvLw8Y4x9+vTp9evXJyUlGXKlUVFRXl5eOl+sQp4fFq1cufLIkSMRERGig9RHkqSsrKzAwEDRQZosMzPTGGMXFBRYW1s7OjqKDtI0lZWVeXl5vr6+ooM0mZG+T3Jzcz08PKytrQ22xoSEhJCQkGnTpul8yTLdgwbQv39/Iz2IjYhMSllZmZ6WbGT/ySUiMh0saCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZIoFTUQkUyxoIiKZYkETEckUC5qISKZY0EREMsWCJuNTLSGvDDcrROcg0jML0QGIGpBzExvTcfkmrpQhrwzKaliZw9MB125h55Ow4D4GtV4saJI1ZTWmb8RTvTA2CB2d4O0Iy98beWkS5sfjw2FC8xHpE3c/SNb+lYgnHsAzfTDYF/7Of7QzgL8+iLwybDonLhyRnrGgSb7O5OPYFczso/UBS6PwYSJybhowE5EBsaBJpqokzNmBzx6BQqH1MW2s8NUYTN+IymoDJiMyFBY0ydR/D+PRIPg5N/CwXu0woRv+sdcQkYgMjAVNcnT+Bn7NwuyHGvXglx5CeiF+zdJzJiKDY0GT7FRLeGELPnsEZtonN2pZNg4L4nH1lj5jERkcC5pk53/JGBGILm5NeIqzDZY8gukbUSXpLRaRwbGgSV5yS7E6Fa+GNvmJD3khrCM+P6aHTESCsKBJXmZvw39H/ul458Z7Kww/n0E1d6KptWBBk4z8cBrdPdC3QzOfbmGGQR2RkK3LSEQCsaBJRv6XhPnhLVrC9GCsOKmjNESisaBJLo5eQe/2sLNs0UK6eyCvDMV3dZSJSCgWNMnFipN4urcOljOlJ34+o4PlEAnHgiZZKK/EmXw85KWDRU3ugTWpOlgOkXAsaJKFTecwLkg3i3K0hpcjTufrZmlEArGgSRZ+OIVpupjfUJkRjJW/6WxpRKKwoEm8y6UwN0M7e50tMMIPBy/xEndk9FjQJF7sb7rcfQagUGBEIHZc0OUyiQyPBU2CSRK2ZGBsFx0vdjpnOcj4saBJsEO56O8JK3MdL9bPGaUVKCjX8WKJDIkFTYKtOInpwXpZ8tQH8ONpvSyZyDBY0CRSeSXOFTb/4hv1i+mO1TwgmowZC5pEWp+G8d30tXB7S3Rzx4lr+lo+kb6xoEmkH05j6gN6XP6MPrx2EhkxFjQJk1UMWwu01d3hz3U97IPkPFRU6XEVRPrDgiZhVp3C0/r5eLCmqM7YkqH3tRDpAwuaxJAk7LiARzrpfUVTe+GHU3pfC5E+sKBJjP2XEOKt+8Of6/J1QmkFyiv1viIinWNBkxixv+EpnZ7eXY9BvjhwyUDrItIhFjQJoKzGuSL0aW+g1Y3qhJ28LgcZIRY0CbAvB4N9Dbe6Bz2RnGe41RHpCguaBIg7h7E6ujx/Y5gp4O2IrGLDrZFIJ1jQJMDRK3jQ06BrHNkJv2YZdI1ELceCJkM7dR0PtIWZwqArHRmIXZkGXSNRy7GgydA2G3Z+Q6V9GxTc5j1WyMiwoMnQ9lzEsAAB6x3ogyO5AtZL1GwsaDKoK2VwtoGthYBVj+TBdmRsWNBkUFszMEbXd7dqpId9cPCymFUTNY/uC1qSpFmzZkVGRkZFReXn56vHKysrp06dGhISEhYWdvHiRZ2vl4zC1vO6v/1gI1mZw8ka+bfFrJ2oGXRf0AkJCQUFBfHx8TExMZ988ol6fMuWLZaWlkeOHHnhhRc++ugjna+X5O92JUor9Ht90foND+TBdmRMdD8XmJiYGBoaCiAkJGTFihXqcUdHx9LS0qqqqpKSEkdHx1rPOnDgQEVFhfrb1NRUHx8fnWcjsX7NxHARHw+qjeqEf+zV7y0CyASVl5dfvnx59+7d6hFra+tBgwa1fMm6L+jCwsKePXsC8PX1LSwsVI8PGTJk/vz5Xbp0KSgoSEtLq/kUSZKOHj1aVfXHZdWzs7OdnJx0no3E2nwOr4aKDBDogqxiSBIUhj0Km1q3mzdvZmdnHz9+XD1ibm4eFhamaPH7TPcF7eLikpOTAyAnJ8fV1VU9/p///GfEiBHvvvvu4cOHp0+f/uuvv6r/SqFQzJ07t+ZCVq5cqVQqdZ6NBKqScK4IPdsKjhHcHieu6es2tWSaOnTo4OPjM23aNJ0vWfdz0OHh4ceOHQOQnJwcFhamHi8qKnJ3dzczM3NzcysoKND5eknmjuTiIS/RIVQH2/GUQjISui/oiIgIDw+PqKiodevWzZ07NzU1tW/fvgBef/31bdu2hYSEPPnkk1988YXO10syF3dO2PEbNUX6Y2+26BBEjaP7KQ4zM7MlS5aov3V3d09JSQHg4eGxfft2na+OjMWBS1gUKToEYG8JADcr4GQtOgpRQ3iiChlCRhH8nGEhj7fbUO5Ek5GQx28MtXZxGbKY31DhOd9kLFjQZAg7LuCRzqJD/K5XW/x2XXQIokZgQZPeFZbDTCGjOV+FAp1dkVEkOgdRQ1jQpHfbL2C0bHafVXiwHRkFFjTp3eZziO4qOsSfjQjEryxokj0WNOlXZTWulMJXZuftu9niZgVvsEJyx4Im/Uq4iCF+okNo8pAXkq6IDkFULxY06dfW88Ku0F+/wb7YlyM6BFG9WNCkX8euIMRbdAhNBvki8ZLoEET1YkGTHqUWoIsbzGR5bU8na5RVQMlpaJIxFjTp0RZxdyBsjH6eOHFNdAgi7VjQpEe7MjEyUHQI7cJ9sS9bdAgi7VjQpC+F5bAwg6NsTiCsa7AvDnAammSMBU36suMCHukkOkS9XG1RVI4qSXQOIi1Y0KQvWzIwLkh0iIYEt8cpXjiJ5IoFTXpRWY3sEgS4iM7RkMF+2M+joUmuWNCkF4mXEO4rOkQjDPZlQZN8saBJL2R+gJ1aW3tcv4VqTkOTLLGgSS8OX8bDHUWHaJyebXGWd5knWWJBk+5lFMHfBeayPIGwLk5Dk2yxoEn34oxkfkNliB+vmkQyxYIm3dtxAaPkfQR0TR3aIK8MEqehSX5Y0KRjxXcBwMVGdI6mCHJDxg3RIYjqYEGTju28IOvrb2jEi3KQPLGgSceM5QC7mvg5IckTC5p0SVmNCzfQ1V10jibydUJ2iegQRHWwoEmXDudioI/oEM0S6IqsYtEhiP6MBU26tCUDUcY2v6HCWxSSDLGgSZeM5RIcdXEammSIBU06k1UMLwdYGud7KtAFmTzSjmTGOH+ZSJZ+PoPHe4oO0QI+Trh0U3QIohpY0KQzOy5gdGfRIVqAlx4luWFBk26cLUBHJ9hYiM7RApyGJrlhQZNurDmLx4x5fgNAkBvOcxqa5IQFTbqxKxPDA0SHaDF3O1y/LToE0e9Y0KQDp/PRxc245zdUBvOiHCQnLGjSgV/O4LEeokPoQrgvDlwSHYLodyxo0oFfszDU+Oc3ADzQFmfyRYcg+h0Lmloq5Sp6tzPW81NqUSjgaI2iO6JzEAFgQVPL/ZJq9Mdv1BTWEQc5y0HywIKmFpEkJFzEED/ROXSH09AkHyxoapFjeXjQy2hu4N0Y/Tog5aroEEQAWNDUQq3m+A01CzNYmaPsnugcRCxoaolqCYcuI6yj6By6FuqNw5dFhyBiQVNLHLqMEG+YtaL5DRVOQ5NMsKCp+VrZ8RtqId44mis6BBELmpqtSkLSFYR4ic6hBzYWqJZwVyk6B5k8FjQ10/4cDPKFotXNb6g85IWkPNEhyOSxoKmZWt/xGzUN4sX7SQZY0NQcymr8dh39PUXn0JuHfXCIB3KQaCxoao49F4317t2N5GiN2/dQJYnOQaaNBU3N8dlR/PVB0SH0rA9PKSTRWNDUZGmFsLGAr5PoHHo2qCOnoUkwFjQ12ZKjeHmA6BD6F+6LRJ6uQkKxoKlpiu4graCVT0CruNuh4DaqOQ1N4rCgqWmWpeCZPqJDGMoD7ZBaIDoEmTAWNDVBZTXWpOLx1nh6t0achiaxWNDUBBvSMDYIVuaicxhKuC8OsKBJHBY0NcFXx1v/0XU1eTvi0k3RIciEsaCpsY5ega8TPOxE5zCsLm44f0N0CDJVLGhqrCVH8Wqo6BAGN4izHCQOC5oaJbcU+bfRs63oHAbHi/eTQCxoapQvkzHLlGaf1Tq74nyR6BBkqljQ1LDySuzKxLgg0TkE8XLE5VLRIcgksaCpYb+kYnKPVnjvwUYK68hzvkkMFjQ1QJLwzXH8pb/oHOIM9cfuLNEhyCTpvqAlSZo1a1ZkZGRUVFR+fn7Nv/r4448jIiIGDBiQlcX3u9GIz0bfDnCwEp1DnG7uSC8UHYJMku4LOiEhoaCgID4+PiYm5pNPPlGPnzhxIi4ubs+ePW+//faiRYt0vl7ShyoJ7+3Fm2GicwilUMDHETk8Y4UMTvcFnZiYGBoaCiAkJOTQoUPq8W3btk2ePNnMzGzs2LGLFy/W+XpJH5YmYVwQvB1F5xBtWABnOUgAC50vsbCwsGfPngB8fX0LC//4n+H169czMzOHDx+uUCj+/e9/u7m5qf9KkqRPP/20oqJCPZKSktKnj8lcM02uCsrx42nsmy46hwxE+OPdBDzLtyRpkp+ff+LEiatX/7gBj7W19Zw5cxQtvum97gvaxcUlJycHQE5Ojqurq3rcwcHh7t2727dvT05OnjlzZnJysvqvFApF7969q6qq1CM3btywtbXVeTZqkvnxeD/ChC6NVI9AF2QViw5BcmVra+vn59evXz/1iLm5ecvbGfoo6PDw8GXLlgFITk4OC/tj8nLgwIF79uyxsLBwdXWtrq6u9ayIiIia3+bl5SmVSp1no8ZLysP1WxgeIDqHbHR1R3ohurqLzkHy4+Dg0KNHj2HDhul8ybqfg46IiPDw8IiKilq3bt3cuXNTU1P79u0L4JFHHlEqlaGhoVOnTv388891vl7SIUnCG7/ik5Gic8hJhD/iL4oOQSZG93vQZmZmS5YsUX/r7u6ekpJSd5zk7IfTCPVBgIvoHHIy1B9zdpjW1VZJON0XNBm70gp8dgwJT4vOITOeDsgtRbVkumdUkuHxTEKq7YMDeCUEdpaic8hPcHuczm/4YUS6woKmPzlXhJSrmGIydx1skgg/TkOTQbGg6U/m/Yp/DxcdQq4i/LE3W3QIMiUsaPrDrky0tUdwe9E55MrDDkXlqJJE5yCTwYKm+25X4v19+NdQ0TnkrW8HnLja8MOIdIIFTff9bQfeHQJ3E7snbFMN5UU5yIBY0AQA286jWuJ5gw0b4of9vIcsGQqPgybcuIOF+7HzSdE5jIGTNW5XoqIK1rxECekf96AJc3bgn5FwtBadw0iEeOPYFdEhyDSwoE3dmrNwsEKkv+gcxiOCR0OTobCgTdr12/jkMBaPEJ3DqIT74iDvIUsGwYI2aS9uwX9GwJ5ndTeFnSWU1bhdKToHmQAWtOmKPQVfZzzsIzqHEQrriMOXRYcgE8CCNlGXS/FlMk9LaaYIfyRkiw5BJoAFbaJe34UPh8GGh1k2S4g396DJEFjQpmhXJjzsMaij6BxGy9YCVuYorWj4kUQtwYI2OcpqLNqP9yMafiTVY5AvTykkvWNBm5wVJxHdFS42onMYuUhOQ5P+saBNS9EdxJ7CywNE5zB+A7yQxPMJSc9Y0KZl0X78fTAs+GNvMTMFvB2RVSw6B7Vq/E01IakFuFiMoTyrW0dGdcKuTNEhqFVjQZuQ+fH4eKToEK3ISBY06RkL2lRsv4AAFwS6iM7RirSzx407uFclOge1Xixok1BZjY8S8e5g0TlanYc74hDPWCG9YUGbhKVJmPIAr/iseyMDseOC6BDUerGgW7/CcmxIw8y+onO0RgN9cDhXdAhqvVjQrd/C/Xh7EMwVonO0RhZmcLPFtVuic1ArxYJu5a7dwtkCjAgUnaP14rEcpD8s6Fbuq+P4Sz/RIVq1qM7Ydl50CGqltBb0yy+/vH///qoqHkNkxCqrsSUD0V1F52jVvB1xpQzKatE5qDXSWtAuLi4vvfSSl5fXX//61/j4eKVSachYpBNbMjC6M0/s1rsBXrzPN+mF1t/df/zjH7/99tuhQ4c6der03nvveXt7P//887t27aqs5L3YjMY3x/EcD97Qv1GdsJPT0KQHDexcubq6+vj4BAYG3rt379ChQ++9956/v/+mTZsME45aIqMINhbwdhSdwwSE+yKR9/kmPdBa0IsXLx4yZIi3t/eyZcv69u17/PjxM2fOHDp0aNWqVS+88IIhI1LzfH0cL/QXHcI0WJnD3hL5t0XnoFZH6z3pkpOTX3755eHDhzs4OKhGbt++bW9v/+CDDy5dutRQ8aiZ7ipx6DIWDxedw2SMCMTuLDzxgOgc1Lpo2INWKpVKpfLIkSPjxo2ztbVVfVtSUtKhQwcA9vb248ePN3hOaprVqYjpDgVPTjGUUZ14zjfpnoY9aBsbGwBVVVWqL9QmTpxooFDUYt//hnWTRYcwJZ1cceEGqiWY8R9F0h0NBa06om7EiBG7du0yeB7SgZSr8HGEM+86aFj9PHHiGvp1EJ2DWhGtHxKynY3XV8fx4oOiQ5geznKQzmnYg+7fv//777//97//ve5fJScn6z8StcjNCpwrxAAv0TlMT4Qf/u8I3hkkOge1IhoK+ssvv/T39//yyy8Nn4ZabtUpPNlLdAiTZGcJSzOU3OXkEumMhimO/v37u7m59e/f39nZuVevXj169Dhy5MipU6d69+5t+HzUJJKEn07zYC9hhgZgz0XRIagV0ToH/f777/fs2bO0tPQ///lPbGzsp59+OmfOHEMmo2Y4nIvg9rCzFJ3DVI0MxE5OQ5PuaD1R5dNPPz1y5Iibm9vSpUuPHj2qVCofeughnqIic18kYUG46BAmrGdbnC3gwXakM1r3oKuqqpydnZOSktq1a9exY0c7O7t79+4ZMhk1VfFdXLuFru6ic5i2fp5IzhMdgloLrXvQjz/++KhRoyorK995552LFy9OmTJl+HCeOCxrWzIwtovoECYvuis2puMhHkVDuqC1oD///PMNGzYAmDBhQlZW1qRJk/7yl78YMBg1Wdw5fDhMdAiTN6gj/p4gOgS1FloL2sLCYtKkSaqvO3fu/NprrxkqEjXHvSpcKUOAi+gcJs/CDL5OOH8DnV1FRyHjp7Wg9+zZs2DBghs3btQcTE9P138kao692RjiJzoEAQDGBmFLBv4WIjoHGT+tBf3MM89MmTLlySeftLDQ+hiSj7gMTOXhz/LwSCdMXsOCJh3QWr6VlZXvvvuura2tIdNQ80gSkq7g01GicxAAwNEa5mYovgsXnlJILaP1MLtXX311yZIlvKu3UTiVjwfa8dhbGXmkE7adFx2CjJ/Wgt64cePChQtdXV2DgoK6/s6Qyajx4s7xADt5GReEuHOiQ5Dx0zrFsWzZMkPmoJbYcxGvDRQdgmro6ISrt3CvClbmoqOQMdNa0Kr95aqqqvz8/Pbt2yt49yS5ulIGR2vY8qNcmRnih73ZGBEoOgcZM61THFeuXImIiHB0dOzevfvx48cHDRp08SKv0yVH285jDOc35GdsF8RliA5BRk5rQc+YMaNnz55FRUVOTk7BwcEhISHPPfecIZNRI23NwLgg0SGojn4dkHIVkiQ6Bxkzrf8xTkxMXL16teq+sRYWFm+88Yavr68Bg1Gj3K7EzQq0sxedg+pQKNCrHX67juD2oqOQ0dK6B925c+fExET1t0ePHg0ICDBIJGqCPVkYxh+LXI0LwmYey0EtoLWglyxZMn369IkTJ964cSMmJmbGjBkff/yxIZNRY8RxfkPGIv2xN1t0CDJmWqc4Bg8efO7cubi4uODg4A4dOnzxxRft2/O/avJSLeFMPh5oKzoHaWFtDjc75JbC21F0FDJO9R2c5ebmNn36dEMloSZLyuPdu+VuTBdsycAL/UXnIOOkeYojOTl54sSJAQEBNjY2gYGBkydPTklJMXAyalDcOYzl/Ia8jenCc76p+TQUdHx8/JAhQ7p06bJq1aozZ87ExsYGBgaGh4fv27fP8PmoHvtzEM4ja+TNzRZ3lCjj3eKoWTRMcbz99tsfffTRrFmzVN926tRp4MCBnp6eb7311qFDhwwbj7S6WAIvR1hq/ZSX5GJEIHZlIqab6BxkhDT8fp88eXLs2LG1BseNG8dZDlnZdh5RnUWHoEYY0wVbeUohNYuGgq6oqHB0rP2ps5OTU0VFhUEiUaNsO4/RLGhj0M0d52+giqcUUtNpPorj9OnTDg4ONUfKysoauURJkmbPnp2WlmZra/vdd9+1bfuno8Bu3rz5wAMPXLp0qXlxSaXkLqqq4cq7KRiJwb7Yl41If9E5yNhoKGgnJ6e6Uxyq8cYsMSEhoaCgID4+fvny5Z988smHH35Y828XLFhQVFTUvKykti8HEfxtNx7jgvDjaRY0NZmGgi4pKWnJEhMTE0NDQwGEhISsWLGi5l8lJSWVlZX5+Pi0ZPkEYE8WpvUWHYIa7UFPvLIDkgRetZeaRPdXES4sLOzZsycAX1/fwsJC9bhSqXz77bd/+OGH8PDwWk+RJCk6Orq8vFw9cvXq1ejoaJ1nazVOXsN/eQdC46FQ4CEvJOXhIZ5Y1BplZWVt3LgxNjZWPWJnZ7dx48aWX0Zf9wXt4uKSk5MDICcnx9XVVT3++eefT548udaUtIpCodi0aVPNkZUrVyqVSp1nax3yb8PZBubcFzMq47thQzoLunUKCAiYN2/etGnTdL5k3R9GGx4efuzYMQDJyclhYWHq8ZSUlDVr1owaNery5cujR4/W+XpNx/4cDPETHYKaKKwj9ueIDkHGRvd70BEREZs2bYqKirKwsPj2229TU1OfeuqplJSUlStXqh7QtWvXbdu26Xy9piMhGzP7ig5BTWSuQDd3pBagh4foKGQ8dF/QZmZmS5YsUX/r7u5e6wyX9PR0na/UpJy+jt7tRIegphvfDevTWNDUBDxT2MhcuwU3O5hxAtoIDQvAnizRIciosKCNzD5OQBsta3O0b4OsYtE5yHiwoI3M3mxE+IkOQc01vhs28SZY1GgsaCNz+jpvoWLEojpjOy8PTY3GgjYmV8rQrg3PRjNibaxgY4Grt0TnICPBgjYm+zi/Yfwe7Yo4znJQ47CgjcnebH5CaPQeDUIcLw9NjcOCNiZneZqD8XO3w10liu+KzkHGgAVtNC7dhKcDJ6BbA95JlhqJBW00eAR0qxHTDRt5Oi01AgvaaOzN5kX6WwlvR1y/hfJK0TlI9ljQRiO9EF3dRIcgHRnZCbsyRYcg2WNBG4ecm+joxAno1mN8V2zgLAc1hAVtHPbnINxXdAjSne4eOF8EZbXoHCRvLGjjsCcLQzkB3bpE+CP+ougQJG8saOOQUYQunIBuXSZ0w/o00SFI3ljQRuBiCQJcRIcgXevXAaeuo0oSnYNkjAVtBBIu8gjo1inUBwcviQ5BMsaCNgIJPAK6lYrphnWc5SDtWNBGIKsYgZziaI1CvHEkFxJnOUgLFrTcXbjBdm61zBR40BPH8kTnILliQcsdj4Bu3R7tyutykFYsaLnjNZJat0geDU3asaDlLvMGOrmKDkF6Y67AA23x23XROUiWWNCydrkUng6iQ5CexXTHurOiQ5AssaBl7eAlDOIEdGs3LICzHKQZC1rWDlzCoI6iQ5CeWZrBzxnphaJzkPywoGXt9HX0aic6BOlfTHdel4M0YEHL1407cLSGGa8BbQJGdcJOXr+f6mBBy9ehy3iY8xumwdYCbe2RVSw6B8kMC1q+OAFtUmK68R4rVBsLWr6S8/Cgl+gQZChjumDbedEhSGZY0DJVXgkFYG0uOgcZShsrtLFCXpnoHCQnLGiZOnYFA7xFhyDD4p1kqRYWtEwduIQwTkCbmOiuiDsnOgTJCQtapo7kIpR70CbG2QYAiu6IzkGywYKWI2U1Sivu/7qSSYnpjjWpokOQbLCg5ejkNfRpLzoEiTCBB9tRDSxoOTrAaySZKjdbmCt4LAfdx4KWo0SeomLCJnbnTjTdx4KWHUnC1TK0byM6BwkyoRs2saAJAAtahjJuoKu76BAkjrMNbCxwhbMcxIKWoQM5nIA2dRN5jxUCwIKWIV4jiaK7YjPPWCEWtAzxLrHkaI02Vsi5KToHicaClhfeJZZUJvXgLAexoGWGd4kllXFBiMsQHYJEY0HLCyegScXBCq62yC4RnYOEYkHLC+8SS2qTe2ANZzlMGwtaRniXWKppbBds5SyHaWNBywjvEks12VmirT0u3BCdg8RhQcvIwct42Ed0CJKTGF6Xw7SxoGXkSC5CeGV8xFUAACAASURBVJF+qmFsF95jxaSxoOWi5C6szWHFu8RSDXaW8HJERpHoHCQIC1ouEi8hnEdAUx2Te2A177FiqljQcpGQjQh/0SFIfh7phJ2ZokOQICxouUjOQ78OokOQ/NhYwMcR5zjLYZJY0LJQdAdtrDgBTZpND8aKk6JDkAgsaFk4kIPBnIAmLYYFIP4iKqtF5yCDY0HLAiegqR5mCozqhB0XROcgg2NBy8KJq+jLCWjSbgZnOUwSC1q8gnI428Ccl+Ag7fyccfsert0SnYMMiwUt3r5sDPETHYJk76ne+OG06BBkWCxo8TgBTY0R0w1refVRE8OCFu/UdfTmNaCpITYW6NMeR3JF5yADYkELdu0WPOx4DWhqlGf64Dt+VGhKWNCC7eUENDVaf0+kFaC8UnQOMhQWtGCcgKYmmcCZaFOi+4KWJGnWrFmRkZFRUVH5+fnq8YqKihkzZkRGRvbt2/fYsWM6X6+ROpOPnh6iQ5Dx4LEcJkX3BZ2QkFBQUBAfHx8TE/PJJ5+ox3ft2tWmTZv4+Phvvvlmzpw5Ol+vMcothacDFJyApkZzs4WTNe+DZSp0X9CJiYmhoaEAQkJCDh06pB739vaePXs2ADc3NwU7CQCwNxsRfqJDkLGZ0YdnFZoKC50vsbCwsGfPngB8fX0LCwvV43369AGQlJT0wgsvLFq0qOZTJEmKjo4uLy9Xj1y9ejU6Olrn2eQmIRtzB4oOQcZmZCAW7kOVxLNP5SIrK2vjxo2xsbHqETs7u40bN7Z8T1T3Be3i4pKTkwMgJyfH1dVVPS5J0vz58w8cOLB8+fLevXvXfIpCodi0aVPNkZUrVyqVSp1nk5v0QnR1Ex2CjI2ZApH+2JWJRzqJjkIAgICAgHnz5k2bNk3nS9b9FEd4eLjqM8Dk5OSwsDD1+Jo1azIzM+Pj42u1s8nKLkFHJ05AU3M825ezHCZB9wUdERHh4eERFRW1bt26uXPnpqam9u3bF8CuXbuOHDnSv3//4ODgsWPH6ny9RocT0NRs/s4ovoOC8oYfSUZN91McZmZmS5YsUX/r7u6ekpICYNmyZTpfl1Hbl4M3wxp+GJFGT/bCipN4nZ9htGo8UUWYc4UI4gQ0NdfUXtiYjrJ7onOQPrGgxcgsRqBrww8j0sZcgb/0w9Ik0TlIn1jQYiRc5AQ0tdTUXth8jjvRrRkLWoy92RjsJzoEGTlzBWYE45vjonOQ3rCgBVBWI7sEgS6ic5Dxm9EH69JwizvRrRQLWoCDlxHWUXQIahXMFZjZF/9LFp2D9IMFLUDcOYwNEh2CWosne2HtWe5Et04saAEO5yLUW3QIai0szfB8P3zFmejWiAVtaGmF6OLGe1yRLk3rjXVneaeVVogFbWibz2FsF9EhqHWxNMMzffAlZ6JbHRa0of2aiZG8CBnp2tPBWMud6FaHBW1Q12/DxgL2lqJzUKtjaYYZffA1Z6JbFxa0QW07jyjOb5B+TA/GL6nciW5VWNAGte08ojqLDkGtlKUZnuyF73id6FaEBW04FVW4fgsdnUTnoNbrub7YlI70woYfSUaBBW04e7IQ6S86BLVqVuaInYAXtqDkrugopAssaMOJy8A4nkBIetbOHv8ahmc3o1oSHYVajAVtIJKElKvo0150DjIBod4I64h/JYrOQS3GgjaQ41fRrwNvEUsG8rcQnMnH9guic1DLsKANJC6DF0gig/pmLD44gKxi0TmoBVjQBpJwkZ8QkkG1scK34zB9I4+MNmIsaEPILoGHPazNRecgE9PFDa8NxAtbROeg5mJBG8L2Czw/hcR4NAhONvj2hOgc1CwsaEPYmsEzvEmYj0dgbzambUBBuego1EQsaL0rrUB5JdrZi85BpsrKHLHj8Xw/jP0RHx1EFY+PNh4saL3byeuLkgyEdUTiM7CxQOT3SMoTnYYahwWtd3HneAIhyYKFGeYMwI8x+F8SXtmBMt7GUPYsRAdo5e4qkVWMbu6icxD9zssByx/F6lSMWgV/Z/T3RH9P9OnAy5TLEQtav9aeRUx30SGI6pjcA5N74NJNJOdh63ks2o/ySnRyRX9PdHZDRyf4OsGG9SAafwL69f1v+GWi6BBEWnR0QkcnTOgGAJKEC8VIuYrkPKxPQ04JKqogSXC3Q0cnDPTBqE6w4162YbGg9eh0Prwc4GorOgdRIygU6OyKzq61xwvLkVuKQ5fx+FpYW2BMF4zpAje+qw2CBa1HXyXjL/1FhyBqGXc7uNshuD3++iCK72JLBl7cgrtKjAjEpB48flS/eBSHvty6h9P5CPUWnYNId1xs8FQvrJ6EnyeifRuM/RFfH4fEA6v1hgWtL6tO4cleokMQ6YedJSZ2x8FnUXwXA5cjmQdW6wcLWl9+PI0nHhAdgkifLM3wxsNYNQELEjBnB27xwGpdY0HrxeFc9GrHA0vJJAS6YNsT6NcBQ1diZ6boNK0LC1ovlqVgZl/RIYgMRaHAtN5Y/xi+Po43d3NWWmdY0Lp34w5yShDM2w+SifFywLrJ6OiESWtwVyk6TavAgta973/D9GDRIYgE+euDmNoLMatRWiE6ivFjQeuYJGFDGib1EJ2DSJzxXTF3IMb/witQtxQLWscSsjHAm3e3IlMX4Yd/D8e4n5DJu9a2AAtax74+jhd49iAR0K8DfpiA6Rvx23XRUYwWC1qXCspRcheBLqJzEMlDgAtWTcDsbbxFQDOxoHXpuxOY0Ud0CCI58XXC+sfwtx04ekV0FCPEgtYZZTU2pmN8V9E5iGTGww7rH8OrO3H8qugoxoYFrTNfJGFid1jx40GiOtraY+1kvLSN89FNw4LWjStlWJ+GOSGicxDJVYc22Pg4/rqVV1ZqAha0bry+C/8ZAXOF6BxEMtbWHusfw8vbkcK5jsZhQevAjguwt8KDnqJzEMleO3usnoS/bkVqgegoxoAF3VIVVVi4Hx8NE52DyEh4O2L1JDy7CWfZ0Q1hQbfU58fwZC/eeJCoCTo6YUU0nt2MLJ5nWC8WdItkl2D7efyln+gcRMamqzt+mYhpG3DymugoMsaCbpE3duO/o2DGzwaJmq6jE9Y9hjk7cDhXdBS5YkE337bzaGePB9qKzkFktNrZY+PjmB+P3Vmio8iShegAxuquEh8dRNwU0TmIjJyLDdY/hphfcLsSjwaJTiMz3INupsWH8Hw/OFqLzkFk/JyssXkKlqVgzVnRUWSGBd0cJ67h4CU80VN0DqLWws4Sv0xE7G/4+YzoKHLCgm6yk9cwexu+i4aCnw0S6Y6dJdZOxoZ0/CuRt529jwXdNMeuYPY2bHwcHdqIjkLU6liZ45eJcLPFqB9w7ZboNDLAgm6CY1fw6k5seBwedqKjELVez/fDPyMxcTUSL4mOIhoLurFU7bz+MbYzkd7198SmKfjoIN7bi2oTnu5gQTdKUt79dm5rLzoKkWlws8Xmx2FriZjVKL4rOo0gLOiGncnHnO34eSLbmcigFAq88TBm9sWYH7Hjgug0IrCgG3DgEp6Pwy+T4O0oOgqRSYrqjA2PYdt5jFxlcieF80xCrTKK8PYe2Fth/WNob+THbEgSrt6CmaIJL6S8EuZmsBZ0B6+rt3A0F/ZWiPTnbRAIbe2x5BFkl+DdvVh8EAsj0cNDdCaDYEFrkF2Ceb/CTIF/DkWQm+g0TXHwMs4XYaAPutSIXXwXk1bD3Q73qlAl4eeJsK33x37jDqauh5kCZRXo0RZfjG7U1aBiT+GrZFiYIbor5gzQepD4zQqU3EWHNvXdvHFvNt7eg0c6Y382Zm/D5scR5N5wgGYrvouF+3AmHz5OeHcwOjrpcV3UEn7O+D4a2SX4MBEld/F+xJ/e560Spzj+pLwS/zyAJ9fjmT74eaKM2vl2JX48jRUncemm1sc8uxmbz+FeFf6egGUpf4y/txevDcTPE7H+McR0w8eHGljXm7vxaii2PoH9M+DjiOUnGo636RwSL2H3NCQ8jRt3sDRZ88MW7ce4n/CPvRj0HY5d0bq0hfvx8QjsysRz/TDIF+N+RspVnLiGeb/inXhcLPnTg0vu4tJNVFY3HFKbJ9cjuit2PYW5AzFtA8orm78oMgA/Z3w5Bi8PwKxtmL4RWzJQUSU6k96woO87V4T/O4KI7+HpgP0zMKqT6EA1FJRjZCyK7gDA0xtx6LKGxxzJhYMVPhqG5/vhxxh8/xuUv3fWmXyMDLz/9dggHG/odnDnb2B4wP2vH++peXW1rDuLBeGwsYBCgfnh2HxOw2OOXUF6IfZNx/JHsXkKXtuldWnVEr47iS9GY3IPPNsHQ/3x2i68txdTe2FsFzxT404cCxIw/he8uxdhy3GiWZcVvlIGV1uE+wJAN3cMC8CRls1yJmRj2EoMXYlpG1BQ3qJFUT0G+uDXp/C3UBy7giErMHU91qe1wn9cdT/FIUnS7Nmz09LSbG1tv/vuu7Zt29Y/LlBhOfZcxK+ZOJOPLm4YHoiEp2FnKTpWHf9LwpthGNMFAMZ0wbQN2Da19mOyS9Dt91k5MwV8HJF/G54OAODrhNSC+5dFTbna8H8LHK1x7db92er0Qvg6N5zQ1hKlFfe/Lq/UPH1x4tr9lwCgnT0crXG7EvaatraLDS6W3J9q2HwOwwKwNg0XXrp/aar/G4Uvk7HkEezPQV4ZEp4GgLwyPLEOe6c3HLUWcwUqlMgtxbkidHRCRRXMW7DTkl2Cf+7H2slwtsGhy5i5GZseb/7SqEG926F3O7wfgbRCrD2Ljw/D0wGjO6NXO/TwgI3xz+Dq/hUkJCQUFBTEx8cvX778k08++fDDD+sfN5ibFcgpQc5NZJfgYjGS8+BojaEBmBMi92s655Xh8d8vzORuhztKDY/p74lXd+K5vjBToPgusktgpsCVMni2wYLBmLIWMd1xrwpbM7CpoUukLghHzGrM7IuyCvx0BlueaDjhI50xaQ2eeADDAvBRouZbzHR2xc7M+y+kogo372puZwCLR2B4LIbHwtUWXdxQUA4nK9hb3f9bV1uU3AWAY3nwdMDmcwj3hacDrMxRWQ3LJtZr+za4cANjf8IjnXHwEnJK8M6gpi2hpoRsTOsNZxsAGOiDe1W4o4S1Oe4q5fgPf2vSzR0LwrEgHJnF2JeNH07hbAHKK9G+DR5oh55t0c0d7drAxUZ00CbSfUEnJiaGhoYCCAkJWbFiRYPjKufOnbt1649z77Ozs9u3b6+TPG/uvj/d6WQDXyf4OcPXGQN98K9hwg5RaKpwX/x8Bn8fDADxF9FN0ydmnVwxoRsGLIOPI66UwdYCf90Kawvk38aaSUiYjgM5sDDDvIcbrrD+ntj4OHZnwcUGCU83vBuSkI1PDmFKT2zNwGfH8N2jeETTBFGEH749gdnb0M0DG9Pxt1CtCwx0QfoszP0Vv11DeiGszTGjDz49gr+FoErC4oMYG4T821hxAu3awMka/z6IpWNwV6n5pVVU4evjSM1HgAtmPQR7S+Tfxmu7kFcGBfBGGKwsENUJZwoQ4g0XW+TcbP5nD07WuFz6x7d3lfgqGT+chr0l2ljhm3GNuoTLvSpUVmv910usI7nYdh7WFni6t0wPPA10QaDLH9/m38ap6zidj91ZuH7r/j/tAMwUaGuP/p54JUQHK1Uqlbm5ucePH1ePtGnTJihIBxe31n1BFxYW9uzZE4Cvr29hYWGD4wAkSYqNjVUq/9gzPHPmzIABA3SS50Pjv9/2Ew9gfjyGx8LWAvZW+GqM5odND8bUXrh+C7G/oV0bPNMHAH7NwoIEfDEaIwI1P0sjDztMafTFVD9KxOYpcLbB/HCsOoVzhZoLWqHADxNw8hpybmLVBLSr96wfK3MsGfXHt1US3tuLiO9RLeGxnpjUHW/sxicj8U0KcksR6Y8Rsfhay2aZshaPdMaroUjOw4RfsPUJvLgVsx9ChB8KyjH2R3R2w4LB9x/8yWFkFDW/oEd1wpgf4euEbh74/iR8nXE6H8dmQqHAiWt4aRvWTq7v6ZKEub8i6QqcbFChxMrx8jq+My4D/0vCW4NQfAcxq7Fuskw7uqa29hgWgGEBtcerJBTc1tlarl+/fuzYsfLyPz5zsLCwWLhwoaLFV7zUfUG7uLjk5OQAyMnJcXV1bXAcgEKhWLRoUc2RlStX1uxrWhQJZTWU1Q3sz1qawdsRJ67hi6j7I8P88cEB/WarqLr/n3oA3T3ww+n6HhzcHsE1/muUW4p/7ENqPh70wt8Hw03LzdHNFVgY8aeRnBIEt8eaSTh+Fbk30ckF4zTtr+SWwtoCz/UFgC5uOHYFx/JwV4kIPwDwsMNz/fBRIiqq7v93al8OYro34jVrYWeJTVPw9XHsuYjhAfjtOh72uX/EYZ/2KGzoM8M1Z2Fvif0zACDlKl7diR9jmh9G5z4/hg2P3Z+rsbPEshS8N0RwpGYzb8o5AQ3y8vKaPHnytGnTdLbE3+n+KI7w8PBjx44BSE5ODgsLa3CcGsnCrLEfenR2w4nfD9VIK9T7gb0+jjj6+zFzP53BYN/GPvGuEk+uxyshOPQsJvfA1HVNuApwiDfiMgCgXwfYW+FBL80PK6+EU4273jjZ4J4Sd2p81n/rHib3xOgfMHcXRq3CmC7wbdnmcrDCa6H47BGMC4KnA7J/PyiwoqrhC4gfycX4bve/7tsBV2V2vc17VX/MpHs68BgVQ9D9HnRERMSmTZuioqIsLCy+/fbb1NTUp556KiUlpda4ztdLaq+GYtxPSMqDtTnWnsUPet4L+/dwTF0PO0vcuIOwjpr3ZDVKzkO47/1Twh72gZcjsm/CvxEHjQCY/RCej8PWDFhb4F4VVkRrflgnV5y/gTP56NkWF0uwLxtvPIx+npgfj6d642wB1qRi11NQKHDhBnyddHwPsyk98ejPuKuEhz1+OIU5DU3a+TojvRB92gNA8V3ZHYTwoCdWncKTvaCsxv8d+ePfEtIfhSTLWxeopjieeeYZ0UGMVWU1DuRAWY2wjgY6fuB2JSzN6js/sK4judiQjo9+/5DgyfX4YGjT9vdvVkBZrXViROVKGebH40op3OywMAKdXFEt4fvfkHQFbe0x+yG46/P6sfeqsC4NBbcxslPDU9u3KzH+Z4zsBA87fP8b/hmJEG89Zmuq8kq8tB1ZxahQYnIP3Xy81josX77cwsJCH1McMvs3mnTE0gyR/gZdYzOOOujbAW/twc5M9O2AhIu4WdHk2RinRuzwejngu0f/NGKmwIxgzAhu2rqax8q8CR+32lti21Rsv4DSCvwY08DnqIZnZ4lvx4kOYWJY0CSMlTnWTMJHB7HkKB70xE9y+kBMFAszjO3S8MPIRLCgSSR3OyweLjoEkVzxWhxERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZIoFTUQkUyxoIiKZYkETEcmUfK8HvW/fPpnf2FuSpIsXLwYE1Lmlu+xlZWUZY+yCggIbGxsHBwfRQZqmsrLy6tWrHTt2FB2kyYz0fZKbm+vh4WFtrdP7S9br4MGDQ4cO1ceSZVrQw4YNs7Wt905zMlBZWbljx4533nlHdJAm27Zt29tvv21u3pQbCMrAoUOH2rZta3RNd/Xq1aNHj/bu3Vt0kCbbsmXLu+++KzpFk23evDk8PLx9+/YGW+OYMWPCwsL0sWSZFrSnp+ekSZNEp2hARUXFypUr5Z+zrq+//nrixIkWFjL96WuTm5sbEBDw6KOPNvxQOcnIyEhNTTXG98nSpUuNMfbhw4eHDx/ep08f0UF0gHPQREQyxYImIpIpFjQRkUyxoJvPzMxM/p9kamRnZ6dQKESnaDIrKysrKyvRKZrMSGMDsLe3Fx2hOSwtLY10g9elkCRJdAYiItKAe9BERDLFgiYikikWNBGRTLGgiYhkigXdNJIkzZo1KzIyMioqKj8/Xz1eUVHh7OwcHBwcHBy8ePFigQnr0pZZ27hMGOOmrmnRokVr166tOSLzDa5SN7bMN3hFRcWMGTMiIyP79u177Ngx9bhRbO2GSdQUe/bsmTRpkiRJ33777RtvvKEeP3v27LRp08Tlqo+2zNrGZcIYN7WKUqkMDw83Nzdfs2ZNzXGZb3BtsWW+wTdv3jx79mxJkpKTk0NCQtTjMt/ajcQ96KZJTEwMDQ0FEBIScujQIfX4+fPn09PTo6OjJ0+efPnyZXEBNdCWWdu4TBjjplYxMzPbs2fPvHnzao3LfINriy3zDe7t7T179mwAbm5uNY/ul/nWbiQWdNMUFhb6+voC8PX1LSwsVI97eHjMnTt348aNEyZMUL1d5ENbZm3jMmGMm1pFoVBYWFiYmdX+5ZL5BtcWW+YbvE+fPkFBQUlJSTExMQsWLFCPy3xrN5KRXc9MlO+++27r1q3BwcEuLi45OTkAcnJyXF1d1Q9Q/VsN4NFHH3377bfFpNRCW2Zt4zJhjJu6fjLf4NrIfINLkjR//vwDBw4sX7685jVdjXRr18I96EaZMWPG2rVr58+fHx4ervogIjk5ueYVYD/66KOlS5cCOHLkSM+ePYUF1URbZm3jMmGMm7p+Mt/g2sh8g69ZsyYzMzM+Pr7WFbeNdGvXJnoS3MhUVVW99NJLo0ePHjduXEFBwZkzZ/r06SNJUlFRUXR09KBBg4YPH37hwgXRMf9EW+Za46Jj1maMm7qmd955R/Vpm7FscJW6sWW+wZ999llfX9/evXv37t17zJgxxrW1G8RrcRARyRSnOIiIZIoFTUQkUyxoIiKZYkETEckUC5qISKZY0EREMsWCJiKSKRY0EZFMsaCJiGSKBU1EJFMsaCIimTKCgnZ2dlYoFAqFwsbGJjQ0dO/evarx5OTk/v37q75+/fXXXVxc8vPz1V+ISrtq1aonn3yykQ+2sLBQKpV6zVPL6NGj09PT1d9u2LCha9eurq6u06ZNu337dj2DycnJffv2dXFxmTFjxp07d1qepKSkxNnZudkhv/rqKz8/Pzs7uyFDhqgfXH/IU6dOjRgxwsnJyc3Nbdy4cefPn29J/vT09K5du9YarPmerKuel1wP1ZtE4+p0bvfu3cHBwfb29mFhYampqfUMqsycOfPpp59WfX3z5k0LC4tPP/1U9W1SUpKVldW+ffvq2SA6VP+WVzt58mQ9F+TbsmVLdHS0TnO1mOirNTXMyclp//79xcXF2dnZn332mb29fXJysiRJhYWFW7duVT3G2dn56tWrNb8QJTY2durUqY18sLm5eWVlpV7zqO3evXvmzJkA0tLSVCNpaWnOzs4JCQkFBQXjx4+fM2eOtsHKykpfX99vvvkmNzd36NCh//znP1uep7i42MnJqXkhz58/b2lpuXv37qtXr86aNSsiIqLBkEql0tvb+5133snMzMzLy5s7d26PHj2qq6ubnT8tLS0oKKjWYM33ZCNfcoNUbxKNq9OtvLy8Nm3arF69uqSkZP78+d27d9c2qBYbG9u5c2fV15s2bbK1tR09erTq2yVLlgwcOLD+DaJDjVzRiRMnevTooe1v4+LiHn30UZ3mainjKOgTJ06ov33rrbcmTpwoSVJSUlK/fv0kSYqOjlYoFD4+Pg8//LDqi/z8/P379wcHB9vZ2Y0cOTI3N1eSpNOnTw8ePHjhwoUPPPCAJEl1H5CWlvbwww8vXrzY09PTz89vz549qjX+8ssvnTp1cnV1feGFF+7evavxuWqxsbETJkx47LHHHB0dBwwYcPr0adV43acMHz4cQMeOHbt166Z6b/373/+2tLS8c+eOJEkRERHff/+9tnU1Prza4sWLZ82aZWdnp+6+xYsXT58+XfV1enq6m5ubtsHdu3d37dpVNZiQkKD+nVT7+uuv/fz8bGxsBgwYkJ6eXk+eJUuWeHt7e3t7f/zxx3XbqpEh8/LyHBwcjhw5UlpaOnfu3JiYmAZDqq7dXlpaqvpWqVSOGTOmuLhY2xau+0OvNZKWltalS5eFCxe2bdvW19dX9QLV78maNL7kuiuNjY2dOXPmU0895eTkNHDgQNVmVL9JkpOT666usrLyhRdecHZ2dnNze//992utd926dV26dHF0dJwwYUJ+fn49PxSVn3/+WX1Pv4qKCoVCcePGDY2D6qdcunQJgOpini+99NIbb7zRpk0b1eaaMmXKW2+9pd4gGqPW3ciNzFx3aeoV1f8aaxZ03TetuqBTU1O9vLwOHjyo8SdV/2bXLeMr6IMHD/r7+0t//mVwcnIqKytTf1FYWOjm5rZ58+YbN27MmjVryJAhkiSdPn3ayclp+vTpZ86c0fiAtLQ0e3v7f/3rX7dv3543b15oaKgkSefOnXNzczt06FBmZma/fv2WLVum8blqsbGxAL766qv8/Px58+Z169ZNqVRqe4pq52jevHlvvfWWJEkTJkxo06bNwYMH7927Z29vn5ubq/GJjQ9fl5eXl7r7Fi1aNHPmTNXXWVlZAG7duqVxcNmyZar7b6rWbmVlVXPf89KlS6r/zBYUFEyfPv3555/Xlmf//v0uLi779u3Lzc2NjIzUtjvZYEhJklSXkFcoFG5uboWFhZIk1R+ysrKyZ8+eI0eO3L59e3l5uXpc48as+0OvO5KWlqZQKD744IPy8vKFCxeGhYVJmgpa40vWuNLY2FgLC4ulS5eq3jm9e/dW5VfvQddd3erVq4OCgi5evJiSkmJtbV3zYs1ZWVlOTk67du0qKiqaPn36Y4891uCbpLS09Pr166qv9+3b5+/vX11drXGw5rMCAwM3b94sSVK3bt1OnDjRt2/f+Ph4SZL8/f137Nih3iB1o9bdpI3PXHdpNQu6nteoLmiNb1pVQV++fNnPz0/1ojT+pOrZ7DpnfAV9/vx5a2trqd6CXrFihWqvSpKkO3fu2NnZKZXK06dPW1lZqf6h1viAtLQ0R0dH1ZzD6dOnVf+jXLhw4UsvvaR65MmTJ/ft26fxuep4sbGxffv2OkuaQQAACCZJREFUVX197949V1fXtLQ0bU9R/e7t2LFj0KBB1dXVHTt2nDVr1uLFi5OSkrp166YtZ+PD11Wz+44ePeri4nLs2LHCwsLJkycDyM/P1zj44YcfPvPMM+oXBeDmzZvqZd65cycnJ0eSpFu3bs2dO1f9e1U3z5w5c958803Vsw4ePNiYgtaYJy0trUOHDocOHSovL3/99dfHjh0rSVL9ISVJunv37tKlS0eNGuXm5jZy5Mhjx45p28J1f+h1R2q+wNTUVNWvfd2C1viSNa40Nja2V69e6vyurq4ZGRlSjYKuu7rVq1cHBAQcPXq0urq6oKCgoqJCvd7//ve/Tz/9tOrr/Px8S0vLRr5JqqurN27c6OXlpWqoegZVnn322bfeeis3N9fd3b2qqurNN9988803r127Zm5uXlpaWrOga0Wtu0kbn7nu0moWdD2vUV3QGt+0cXFx4eHhPXr0UN/FXONPqp7NrnPGd0/C/Px8T0/P+h9z+fLlXbt2+fn5qb61srJSfWzo4+NjbW1dzwPat29vYWEBQPUngNzc3M6dO6u+Vt1TZ//+/XWf26FDB/Xa/f39VV9YWlr6+fldu3ZN4+rUTwkLCzt58uSFCxc8PDxGjBixYsUKKyurYcOGacvZ+PD1e+ihhz744IPx48crlcqXX3553bp1rq6uHh4edQddXFxu3bqlelZpaamFhUWbNm3Uy7GwsFi2bNn27dudnJysra0dHBxU43XzXLt2TfW6AAQEBDQ75IoVK0aPHq26V97ChQudnJxu3rxZf8h79+5JkvTiiy+++OKLFRUVP/3006BBgw4cOKBxY9b9of/444+1RtLT09UvsO6NVtU0vmRtP8Fa75wrV66oV1pze6pXN2HChNLS0ueff/769euzZs167bXXaq5XvXwPDw8rK6uCggI09CYpKip67rnnLl26tHHjRvVnbhoH1YYMGbJs2bKuXbsOHTrUzMxsxIgRc+fODQkJ6devn/qdoDFq3Y28bdu2Rmau54U3+BpVtL1p9+/f/+qrr3799dcffPCBl5eXxp9U/WvXLSM4iqOWLVu29OvXr/7HdOjQYfjw4dnZ2dnZ2ZmZmb/++mv79u1R4wem7QE1b9uu0q5du9zcXNXXhw8fjo2N1fZctYsXL6q+qKyszMnJ8fHxqf8p9vb2wcHBn3/+eUhIyMCBAw8fPnzw4EHVb7XGJzY+fP3KysrGjh2bm5t77dq1iIiI7t27m5ubaxwMCAhQf3yflpbm5+dXs5LWrFmzdevWnTt37t69e8qUKerxunk8PT0zMzNrbaVmhKyqqqqqqlI9QJKkqqoqSZIaDBkVFaX62traevr06aGhoSdOnNC4Mev+0OuOaHyBdWl8ydp+guoHKJXKS5cu1fxXX+PqsrKyIiMjT548efTo0bi4uG+//Vb9V+3bt1dNuwNQ7eW5u7vXn7miomLEiBHdunU7evSouog1DtY0ePDgpKSkHTt2DB06FMDAgQMzMjLi4uIGDx5cf9S6m7Txmet54fW/RjVtb9oRI0Z8/PHHEydOnD9/PrT8pOpfu24ZR0GXlZWVlJRcunTpiy++WLJkyVtvvVX/46Oiog4cOLBt27bCwsI333zzlVdeqfUza/ABajExMbGxsUePHs3KynrllVcKCwsbfO7JkydVU9XvvPOOv79/QEBAPU9R7fRFRkZ+++23ISEh7u7uLi4uO3fuVL2/NT6x8eHrl5ubGxQUlJiYeOnSpTfffPPFF1/UNjhkyJDi4uK1a9feunVr8eLFtY4jLCoqatOmja2tbX5+/meffVbPQXgTJ0786quvDhw4kJeX9/e//70xsTXmGTt27IYNG3bv3l1UVDRv3rxBgwY5OzvXH3L48OEpKSnvvvvu+fPnz50797///e/48eMREREaN2bdH3rdkUZuZI0vWdtP8NSpU1999VVhYeGCBQs8PT3VO5jq/xnUsnnz5ilTply/fr2qqqqiosLW1lb9V2PHjl2/fv2ePXuKi4vnzp07fvz4Bv9ftXHjxqqqqueee+7y5cuqSqqqqtI4WPNZqv2PtWvXqgra2to6PDx85cqVtQq6btS6m7Txmet54Y2k7U2rWtSiRYvWrl2bkpKi8SfV8rU3gf5mT3TFyclJFdXa2nrAgAEJCQmq8XrmoCVJ2rlzZ48ePezs7CIiIjIzM6U6E1J1H1DzSKaaX3///ff+/v6Ojo4zZsxQzTfVfa5abGzs888/Hx0d7eDgMGjQIPUHCBqfMnnyZAcHh1u3bu3btw/A+fPnJUl67rnnan6yofGJjQ9fS83pXUmSPv/8cw8PDy8vr/fee0/94Y/GwaSkpN69e7u6uk6fPl01j69WUlIyfPhwV1fXgQMHxsXFtWvXbuXKldryfPbZZ97e3l5eXsuXL/fy8mp2yHXr1gUFBTk4OKj2rxsMKUlSRkZGVFRUu3btHBwcQkNDt2/fXs8WrvtDrzVS6wVqm4PW9pLrrjQ2Nnb06NETJ050cHAICQk5e/as6pGqN0lycnLd1ZWVlUVHR9vb26sOhFBN46itXbu2S5cuDg4O0dHRqg/66n+TzJs3r1Y5FBQUaBys9QKfffZZ1ef2Kv/3f/9nZmZWUlJSc4NojFp3Izcyc92l1ZyDruc1quegNb5pax5m9+677w4ePLi6urruT6r+za5bvGkskSysWrVqy5YtP//8s+ggJCPGMcVBRGSCWNBERDLFKQ4iIpniHjQRkUyxoImIZIoFTUQkUyxoIiKZYkETEckUC5qISKZY0EREMsWCJiKSqf8HR7fWQCuke+4AAAAASUVORK5CYII="
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It seems that the differences are systematically positive, which is good news because it means that the lakes were, on average, cleaner in 1990 than in 1980.\n",
      "\n",
      "## Reference distribution of signed sums\n",
      "\n",
      "To conduct the test we want the distribution of the possible sums resulting from arbitrary changes in the signs on the differences.  There are $2^{22}$ or `r 2^22` such sums.  Enumerating all of them is not easy in R.  Certainly it would be a mistake to try to do so by creating 22 nested loops!\n",
      "\n",
      "However, there are ways of generating all the possible combinations of -1 and +1, say by using `expand.grid`.  For 3 differences it would look like"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(gg <- expand.grid(c(-1,1),c(-1,1),c(-1,1)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "  Var1 Var2 Var3\n",
        "1   -1   -1   -1\n",
        "2    1   -1   -1\n",
        "3   -1    1   -1\n",
        "4    1    1   -1\n",
        "5   -1   -1    1\n",
        "6    1   -1    1\n",
        "7   -1    1    1\n",
        "8    1    1    1"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(mm <- Secchi$diff[1:3] * t(data.matrix(gg)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "      [,1]  [,2]  [,3]  [,4]  [,5] [,6]  [,7]  [,8]\n",
        "Var1 -1.56  1.56 -1.56  1.56 -1.56 1.56 -1.56  1.56\n",
        "Var2  0.07  0.07 -0.07 -0.07  0.07 0.07 -0.07 -0.07\n",
        "Var3 -0.75 -0.75 -0.75 -0.75  0.75 0.75  0.75  0.75"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "colSums(mm)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "[1] -2.24  0.88 -2.38  0.74 -0.74  2.38 -0.88  2.24"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For 22 differences, writing the vectors to expand will get rather tedious but we can use the `do.call` and `lapply` functions to generate the call to `expand.grid`.  The reference distribution then becomes"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "system.time(\n",
      "    refd <- colSums(Secchi$diff * \n",
      "                        t(data.matrix(\n",
      "                            do.call(expand.grid,\n",
      "                                    lapply(1:22, function(i) c(-1,1))\n",
      "                                    )))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "   user  system elapsed \n",
        "  6.803   2.654   9.731 "
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gc()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "          used (Mb) gc trigger   (Mb)  max used   (Mb)\n",
        "Ncells  373966 20.0     667722   35.7    549655   29.4\n",
        "Vcells 5114376 39.1  226413700 1727.4 281890290 2150.7"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The observed reference density is, not coincidentally, very like that of a \"normal\" or Gaussian distribution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "densityplot(~refd, plot.points=FALSE,\n",
      "            xlab=\"Reference distribution of signed sums\",ref=TRUE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": []
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3deXxU5aE+8GeSyb7vCdkTAggRI1sSCUvYpI1KkEX5uQFel4qKLYtWW7VKF6+KFau1vSgq9F6LKKLVVpQgWxYJECCYsCQhJIEsQyCELckk5/fHtGMkM1nPzHvOmef7B5/k5OScZ2aYh+GdM++rkyQJRESkPE6iAxARkWUsaCIihWJBExEpFAuaiEihWNBERArFgiYiUigWNBGRQrGgiYgUigVNRKRQLGgiIoViQRMRKRQLmohIoVjQREQKxYImIlIoFjQRkUKxoImIFIoFTUSkUCxoIiKFYkETESkUC5qISKFY0ERECsWCJiJSKBY0EZFCsaCJiBSKBU1EpFAsaCIihWJBExEpFAtaNkajUXQEm9Dk7eKNUhGt3q7eYEHLZvHixefOnROdQmY7dux47bXXRKeQX3Z2tugI8nvttdd27NghOoXMGhsbFy9eLDqFMCxo2bS0tLS1tYlOIbPW1tbW1lbRKeR36dIl0RHkp8kHy2g0trS0iE4hDAuaiEihWNBERArFgiYiUii96AC9UlNT88UXX4hO0YPy8vL169f7+PiIDiKnkpKSioqKv/71r6KDyOz06dPau1F79+41GAwVFRWig8ipubm5vLxc+Q9WVlZWZGSk7IfVSZIk+0Fl98EHH+Tn52dmZooO0p2KiorY2FgnJ039p+TixYsXL14MDw8XHURmZWVliYmJolPIrLa21tvb29vbW3QQOXV0dFRWVsbHx4sO0p3t27enpaXde++9sh9ZHa+gAYwZM2bevHmiUxARXau5udlGR9bUyz0iIi1hQRMRKRQLmohIoVjQREQKxYImIlIoFjQRkUKxoImIFIoFTUSkUCxoIiKFYkETESkUC5qISKFY0ERECsWCJiJSKBY0EZFCsaCJiBSKBU1EpFAsaCIihWJBExEpFAuaiEihWNBERArFgiYiUijVrOpNNEAXW3H0rIXtQ4Pg7Wr3NES9wIImh3C5DXM2YmQYnHU/2t4u4VAdNt8BTxdByYisY0GT9rW0Y+5GPJWBzDgLP91+EvM+wid3wM3ZzrmIesAxaNI4Ywfu+hgPjbHczgAy4/DAKNz1MYwddoxF1AssaNIyScJ/fYbpiZg1tLvdsodheiL+6zNIkr2SEfUCC5q07JkcJATgodE97/nQaMQH4FfbbZ+JqNdY0KRZr+TiihHPTurt/s9NwuU2vJJry0xEfcGCJm165wCK67F6Rt9+69UZ2HcG7xywTSaiPmJBkwZtP4m/F+PtW6DT9bxzZ046rJuFvxfj25M2CUbUJyxo0qBVO/HhXLj36yJSdz3+by5e3Cl3JqK+Y0GT1vzzBNKiEOjR/yMEeSA1Ev86IV8mon5hQZPW/DEfS1MHepAn0vDHfDnSEA0AC5o0ZdcpDAlCqNdAjxPqhcGB2H1KjkxE/cWCJk15NRcrx8tzqCczeMkdCcaCJu3YdwYBHoj2ledo0b7wd8f+M/IcjagfWNCkHa/k4kmZXj6bPMUX0SQUC5o04uhZGDswLFjOYw4LRms7jlmaRZrIDljQpBH/vUe20efOnszAf++R/7BEvcGCJi2obMKZZowdJP+Rxw5CTTNONcl/ZKIesaBJC17Lw/KbbHXw5TfhNV4TTSKwoEn16i+huB5T4m11/KnxOFyHhsu2Oj6RNSxoUr01BXh0nG1PsWQc1hTY9hREXbGgSd3OX8X2kz0smDJw2UORU4GmFtuehegaLGhSt3VFuP/GPk8r2lc6HRbfiHWcJ5rsiwVN6vbZUdw10h4nunskthy1x4mIzFjQpGLHGxHtCzdne5zLzRlRvjjRaI9zEZmwoEnFNn2PucPtd7q5w7Hpe/udjogFTSr2dRluHmy/080cjK/L7Xc6IhY0qdWJRgzysdP4hombM8K9UXbOfmckB8eCJrX6uATzRtj7pPOG42OOcpC9sKBJrbaWYaYdxzdMZg7GV2X2Pik5LBY0qVLZOYR723V8w8RdjzAvjnKQnbCgSZU+/h7z7Hj9RmfzRuCTEjGnJkfDgiZV+krE+IbJTwbjqxNiTk2OhgVN6lN2DmFecNeLObu7HqFeKOcoB9keC5rU5xMR1290Nnc4RznIHljQpD5by/ATQeMbJj9Jwr84ykG2x4ImlTl5HsGewsY3TDz0CPZEJdfBIhtjQZPKfFxi1/k3rJnLT6yQ7bGgSWX+eVzw+IbJT5PwT45ykI2xoElNKpsQ7AlPF9E5AE8XBHpwlINsiwVNavKxfecX7R6v5SBbY0GTmvzzBH6aJDrEf2Ql4cvjokOQprGgSTUqmxDooYjxDRNPFwS44xRHOchmWNCkGp8o4/qNzjjKQTbFgibV+OoEbk4UHeLHbh7MT6yQDbGgSR0utqKtA75uonP8mJ8bWttxqU10DtIoFjSpw65TmBgrOoQlE2Oxq1J0CNIoFjSpwzflmJ4gOoQl0xPxDVeSJdtgQZM67K3BuEjRISxJjcR3NaJDkEaxoEkFzlxEoAf0ivzbqndCoAdqL4rOQVqkyL/yRD+WU4FpihzfMJmagJwK0SFIi1jQpALflCu6oKclcBiabIIFTSpQasCwYNEhrLsuGCUG0SFIi1jQpHQlBlyn4HY2GRaMUnY0yY0FTUqn8PENE45ykC2woEnptqmkoLfxfUKSGwuaFK2tA2evINRLdI6ehHmh4RLaOkTnIG1hQZOifVeDVEV+PqWr1Cjs5SdWSFYsaFI0VQxAm3AYmmTHgiZF21Wp0DmSupoUi12nRIcgbWFBk3JdaIGTTkFLqHTPlPNCi+gcpCEsaFKuHZWYHCc6RF9MjsNOTj1K8mFBk3KpaADahMPQJC8WNCnXvtMYPUh0iL4YMwj7zogOQRrCgiaFqr6AcG8460Tn6AtnHUK9UH1BdA7SChY0KdQ2ZU8xas00Tj1K8pG/oCVJWrJkyZQpU7Kysurr66/56apVqzZt2tTjbkTflGO6wtbw7o3pHIYm+chf0Nu3b29oaMjJyZkzZ87q1avN29vb2ydNmvT88893vxsRAEnCiUYkBojO0XeDA3G8EZIkOgdpgvwFvXv37vT0dABpaWm5ubk/nMnJadu2bStXrux+NyIAxQ24PlR0iP5KDsWRBtEhSBP0sh/RYDAkJycDiI2NNRh+mCJXp9Pp9XonJ6fudwMgSdKKFSuam5vNW44dOzZt2jTZo5JifV2mygFok2kJ+Locyar9B4b66syZM998882ePXvMW3x8fF5++WWdbqDvcctf0AEBAZWVlQAqKysDAwP7sZtOp1u4cGFLyw8fyfryyy+DgoJkj0qKtaMSa28THaK/MuPw4Of4eZroHGQvgYGBmZmZWVlZ5i1eXl4Db2fYoqAnTpy4du1aAIWFhRkZGf3bzfTi2uzIkSNGo1H2qKRMHRIaryDEU3SO/gr1guEyOiQ4qeoaQeo3Nze3uLi40aNHy35k+cegMzMzQ0JCsrKyPv744+XLlx85cmTUqFE97iZ7DFKv4nrVjw9wGJpkIf8raCcnpzVr1pi/DQ4O3r9/v/nbVatWWdyNyGz3KYyPFh1iYMbHYPcpFb/PSQrBD6qQ4uypQkaM6BADkxGDPZx6lAaMBU2KU3kecf6iQwxMvD9OnhcdgtSPBU3KcqoJ0X6iQ8ghyhdVnJSDBoYFTcqigQFoE9MwNNFAsKBJWTQwAG3CYWgaOBY0KcvhOowMEx1CDiPDcJiTgNHAsKBJQc5fhY+bRj7f4ayDlwvOXxWdg9SMBU0KkluF9CjRIeSTHo28atEhSM1Y0KQgmhmANuEwNA0QC5oU5LsapGroFXRaFL6rER2C1IwFTUrR0g5jBzzkn31AGA892jrQ2i46B6kWC5qUYt9pjFHVGt69MTqC63xT/7GgSSl2ndLUALRJRgx2VYoOQarFgialyK3SyGcIO8uIQW6V6BCkWixoUoQOCYbLCFbtJP3WBHvCcJlryFI/saBJEUoMGBEiOoRtXBeCEkPPuxF1xYImRdh9CuM1NwBtMj6asyZRP7GgSRH2aPEdQpOMGOzhMDT1CwuaFKHsHBIDRIewjcGBKGsUHYLUiQVN4lVfQKSP6BC2FOGD082iQ5AKsaBJvD1Vmh2ANuEwNPUPC5rE0/AAtAmHoal/WNAk3sE6pISLDmFLKeEoqhUdglSIBU2CNbXAQw9nTUzSb43eCR56NLeKzkFqw4ImwfKrka65T3h3lRaFPI5yUB+xoEkwzQ9Am3AYmvqBBU2CFdRgXKToELaXGoUCLn9FfcSCJpGMHbjcBh9X0Tlsz8cVF1th7BCdg1SFBU0iHa7HyDDRIexlZBiK60WHIFVhQZNI+dVIdYDxDZPUKORzlIP6ggVNIhVUI01Dq8R2Ly0KBVxDlvqCBU0inWhEUqDoEPYyJBDHz4oOQarCgiZhzl6Bnzt0mv6ISmc6Hfzc0XhFdA5SDxY0CfNdjQMNQJuMi8R3HOWgXmNBkzAONQBtksb3CakvWNAkTEENUh2soFP5Cpr6ggVNYnRIaLoKPzfROezL3x3nr6KDi3xT77CgSYxSA67T6DLe3RsWjKO8loN6hwVNYuQ73gC0CYehqfdY0CRGQY3jFjRnTaJeYkGTGN83YIRDDnEkh+L7BtEhSCVY0CTAhRZ46OHkMB9R6cxJB3c9LnJ1FeoFFjQJUHjaIeaAtmZsJPaeFh2C1IAFTQLkVzvcFdCdpUbyfULqFRY0CeCw7xCa8H1C6iUWNAlQfwkhnqJDiBPqhbpLokOQGrCgyd7KzyExQHQI0RICUHFedAhSPBY02VueYw9Am6RGIo+LfFNPWNBkbwXVSHf4gk6P5uoq1DMWNNnbwTrcEC46hGgp4ThYKzoEKR4LmuzqihHOOrg4/N87FyfodLhiFJ2DlM3hnyhkX/vPYPQg0SGUYXQEDpwRHYKUjQVNdpVf7XDLXFmTymntqCcsaLIrB1zmypq0KL5PSD1gQZNdVV9AlK/oEMoQ7YuqJtEhSNlY0GQ/Nc2IZDt3MsgHNc2iQ5CCsaDJfjgAfY1UTspB3WJBk/047DJX1nD5K+oeC5rs58AZjIoQHUJJRkVgP6+0I+tY0GQnxg60tsPTRXQOJfFyQUs7jB2ic5BSsaDJTg7VYWSY6BDKMzIMh+tFhyClYkGTnRTUcBI7C1Ij+T4hWcWCJjsp4CUclqTy4ypkHQua7KTsHJICRYdQniGBOH5WdAhSKhY02cO5q/B1g04nOofy6HTwdcO5q6JzkCKxoMke9tZgLCexs2JcJApPiw5BisSCJntw8GW8u8fPE5I1LGiyh+9qMI7vEFqRGonv+D4hWcKCJntovIJAD9EhlCrQA2eviA5BisSCJpsrO4fEANEhlC0hAOXnRIcg5WFBk80VVPMjKj1IjeTV0GQBC5psrqCGH1HpAd8nJItY0GRzh+twQ7joEMqWEo5DdaJDkPKwoMm2WtoBwIV/0brl4gTpP/cVkRmfN2RbB2uRwpfPvcAX0dQVC5psK5/vEPZOaiRXV6FrsaDJtvgZwl5K4/uE1AULmmyr8jxi/USHUIM4f1Q2iQ5BCsOCJhsyXEawp+gQ6hHoAcNl0SFISawW9OOPP75z5872dr6vTP3HKTj6ZFwk9nJaO+rEakEHBAQ89thjkZGRjzzySE5OjtFotGcs0gYOQPcJh6HpGlYL+je/+c3Bgwdzc3MHDx78/PPPR0VFPfjgg1u3bm1ra7NnPlK1vTUYy1fQvcZX0HSNHsagAwMDo6OjExMTW1tbc3Nzn3/++fj4+C1bttgnHKlah4TmVvi4is6hHj6uuNACSRKdgxTDakG//PLLkydPjoqKWrt27ahRo/bt21dcXJybm7thw4aHH37YnhFJpY43YmiQ6BBqMyQIxxtFhyDF0Fv7QWFh4eOPPz59+nQfHx/TlkuXLnl5eY0dO/att96yVzxSMU5i1w+mae2G8B82AmDxFbTRaDQajfn5+bfddpuHh4fp2/Pnz0dERADw8vKaPXu23XOS+nASu37gtHbUmYVX0O7u7gDa29tNX5jNnTvXTqFIE47UY0So6BBqkxyKIw2iQ5BiWCho0xV1M2bM2Lp1q93zkEZcboOrM5x1onOojbMOeidcMcLD6ugjORCrbxKynWkgDtRiVIToEOo0KgIHzogOQcpg4Z/pMWPGvPDCC88++2zXHxUWFto+EmkB3yHsN9P7hDdFi85BCmChoN9+++34+Pi3337b/mlIMwpqsOB60SHUKT0aG4+IDkHKYPkVNICgoKATJ07ExMS0t7e/8847np6e99xzj93jkVqdbkaEt+gQ6hThjZpm0SFIGayOQb/wwgvJyckXLlx45ZVX1q9f//rrry9dutSeyUi9apoxyEd0CDWL8MaZi6JDkAJYLejXX389Pz8/KCjorbfe+uijjzZv3rxx40Z7JiP1KqjmHEkDkhbF1VUI6Kag29vb/f399+7dGxYWFhMT4+np2draas9kpF751UhnQQ9AejQLmoBuCvrOO++cOXPmggULHn/88YqKiuzs7OnTp9szGanXgVrcyGvsBmBUBPbzSjvqZi6OP/3pT5s3bwZw++23l5eXz5s376GHHrJjMFKr1na0d8DNWXQONXNzhrEDre1w5d3o2KwWtF6vnzdvnunrpKSkZcuW2SsSqdvBOqSEiw6hfinhOFSHMYNE5yChrBb0tm3bfv3rXzc2/mjqw9LSUttHInXL5zuEcjC9T8iCdnBWC3rx4sULFiy4++679XpOCkB9kFeFl/huxYClR+GX2/DoONE5SCir5dvW1vbcc895eHjYMw1pQPUFRPuKDqF+MX6oahIdgkSzehXHL37xizVr1nBVb+qT2osI4wcIZRLqhVp+XMWxWX0F/emnnxYVFf3ud78LDw/X6f49ayTHoKl7HICWUVoUCmowa6joHCSO1YJeu3atPXOQNuRX4zYWikzSo/GPYyxoh2a1oIcNGwagvb29vr6+84toom7sO4PnJ4sOoRWjI/Cbb0WHIKGsjkHX1NRkZmb6+voOHz583759EyZMqKiosGcyUp22DrS1w51X/cjEXY+2Dhg7ROcgcawW9KJFi5KTk8+ePevn55eSkpKWlvbAAw/YMxmpziF+REVuI8NwqE50CBLH6qud3bt3b9y40bRurF6vf/LJJ2NjY+0YjNQnvxrpXAdEVulRyK/m4mGOy+or6KSkpN27d5u/LSgoSEhIsEskUqu8Kl7CIbO0KORxWjsHZvUV9Jo1a+bMmTN58uTGxsY5c+bs2rVrw4YN9kxGqlN1AbF+okNoS5w/Tp4XHYLEsVrQkyZNOnr06Oeff56SkhIREfHmm2+Gh3N8kayqu4RgT9EhtCjUCw2XEcL71iF19457UFDQwoUL7ZWE1I2rqNhIaiTyq3HrENE5SATLY9CFhYVz585NSEhwd3dPTEycP3/+/v377ZyM1IWrqNgIV1dxZBYKOicnZ/LkyUOGDNmwYUNxcfH69esTExMnTpy4Y8cO++cjtdh3BqM5N6YNjBmEfadFhyBBLAxxPP300y+99NKSJUtM3w4ePPimm24aNGjQL3/5y9zcXPvGI3UwLf/hwY+o2ICHHleNMHZAb/WSK9IsC495UVHRrbfees3G2267jaMcZM3helwfKjqEdl0fhuJ60SFIBAsF3dLS4ut77YS+fn5+LS0tdolE6pNXxY+o2FA6r4Z2VJb/U3r48GEfH5/OW5qbm+2Sh1Qpvxq/yRQdQrvSovCbHfjZGNE5yO4sFLSfn1/XIQ7TdtvnIVU6eR7x/qJDaFdCACrOiQ5BIlgo6PPn+dEl6oOGy/yIis0FecLA+9nx8I1hGqiCaqTyCmgbS41EQY3oEGR3LGgaqDx+RMX20qORVyU6BNkdC5oGat9pjI0UHULrxg5CIT+u4nhY0DQg7RKuGvkRFZvzdMFVIzok0TnIvljQNCAHa3Ejp5O3i5RwHOTqKg6GBU0DsusUJsSIDuEYJsRiV6XoEGRfLGgakN2nMIFLodlFRgx2nxIdguyLBU39J0mov8S55O0kzAt1l0SHIPtiQVP/HW/EkCDRIRxJUiCON4oOQXbEgqb+4wC0nXEY2tHIX9CSJC1ZsmTKlClZWVn19fXWtre0tPj7+6ekpKSkpLz88suyxyA72FXJAWi7mhCDXRyGdiTyF/T27dsbGhpycnLmzJmzevVqa9vLy8tnzZpVVFRUVFS0YsUK2WOQHXCOJDtLCOAi345F/oLevXt3eno6gLS0tM4rsFyz/fjx46WlpdnZ2fPnz6+q4odY1af6AiJ8et6N5BXmhdOc+tdhyP8JMIPBkJycDCA2NtZgMFjbHhISsnz58nnz5n344YePPvroli1bzHtKkrR06dLO6wMcO3Zs8uTJskelgdhTxQFoASbEYvcpzB8hOgd1UlVV9e233+7Zs8e8xc3N7fXXX9fpdAM8svwFHRAQUFlZCaCysjIwMNDadtOraQCzZs16+umnOx9Bp9MtW7asvb3dvGXz5s1eXl6yR6WB2FWJhziFvN1NiMH/7GdBK0tYWNgtt9wye/Zs8xZnZ+eBtzNsUdATJ05cu3YtgMLCwoyMDGvbX3rpJR8fn0ceeSQ/P9/0yrqz2NgfvfcUEhJiNBplj0oDcaQBySGiQzie5FAc4fqECuPq6urp6ZmQkCD7keUv6MzMzC1btmRlZen1+nfeeefIkSP33HPP/v37r9nu5OR0//33f/jhh+7u7n/+859lj0E21XgFvm6Q4yUC9Y2TDt6uaLyCQA/RUcj25C9oJyenNWvWmL8NDg42LQd+zXYAmzdvlv3sZB+5VRjPVWIFGR+DvGpkJYnOQbbHD6pQf+ziFBziTIjhx1UcBQua+mP/GYwZJDqEoxobiX1nRIcgu2BBU59dboOzDi78uyOIixOcdLjcJjoH2R6fZNRn39VwlVjBUiOxlytgOQAWNPUZ50gSjrMmOQgWNPVZfjXSeQmHUOlRyK8WHYJsjwVNfWPswJU2eLmIzuHYvF1xqQ3GDtE5yMZY0NQ3RVwlVhlu5BqyDoAFTX3DAWiF4DC0I2BBU9/sPoUMFrQCTOAasg6ABU19IEkwXEYwV4lVgGBP1F+CJInOQbbEgqY+OHoWw4JFh6D/GBaMY1xDVtNY0NQHHIBWFA5Dax4LmvqAq8QqCteQ1TwWNPVBZRNi/USHoP+I80cl15DVNBY09VbVBUR4iw5BPxbmjeoLokOQzbCgqbe2lmGq/Gv60IBMjcfWMtEhyGZY0NRb28oxI1F0CPqxaQnYViE6BNkMC5p6RZJw8jwHoBUnIQDl53g1tGaxoKlXDtdjZJjoEGTJyDAUN4gOQbbBgqZe+aacA9AKNTUe35SLDkG2wYKmXtl+ElPiRYcgS6YmYDuHoTWKBU09a2nH+asI8hCdgywJ8sD5q2jj3NBaxIKmnhVU4yYuoaJgaVxgRaNY0NSzbRWYyvENBZuagG0chtYiFjT1bPcpTOQUHArGSTm0igVNPbjQAr0T3PWic5B1ni5w1uFCi+gcJDcWNPVgRyUm8eWz4k2Kw05OPao5LGjqwbZyTOMV0IrHz3xrEguaerDvDEYPEh2CejJmEPadFh2C5MaCpu7UNCPEE8460TmoJ846BHnizEXROUhWLGjqTk4FP+GtGlPjebGd1rCgqTscgFYRDkNrDwuaunP0LIYGiQ5BvTMsGKUG0SFIVixosqrUgOuCRYegvhgWjKNnRYcg+bCgySpOMao6nHpUY1jQZFVOBQegVWZ6InI4DK0hLGiyzNiBhssI8xKdg/oizAt1F9HOFbC0ggVNlhWexhh+PkWFxgxCIT+xohUsaLLsm3JOMapKUxM4DK0dLGiybGclJsWJDkF9N5mzJmkIC5osuGKEsQM+rqJzUN/5uKKtHVeMonOQHFjQZEFOBSZwilHVyojhMrIawYImCzaXYPYw0SGov2Zfh82lokOQHFjQdK12CYfrkRIuOgf1143hOFSHDl5sp34saLpWQTVSI0WHoIEZF4mCGtEhaMBY0HStLUcxi+MbKjdrKLZwlEP9WNB0rZ2VXMNb9bhEoTawoOlHjp5FQgBc+PdC5VycEOePY5zZTuX4RKQf+ewobhsqOgTJ4bah+Oyo6BA0MCxo+pF/HsdPBosOQXL4SRL+eUJ0CBoYFjT9oPYi3PTwdROdg+Tg5wYXJ9RdEp2DBoAFTT/48jhuHSI6BMnn1qH48rjoEDQALGj6wWe8wE5bsodxGFrdWND0b5facPYKIn1E5yD5RPrAcBmX20TnoP5iQdO/fVOOGYmiQ5DcpnN6aDVjQdO/bSnFLF5gpzmzhmELRzlUiwVNANAu4fsGjAwTnYPkdkMYjtRzlUK1YkETAORVIS1KdAiyjdQo5FeLDkH9woImgBMkaRonTlIvFjQBwC5OkKRdE2OxgxMnqRMLmlBiQFIQnHWic5Bt6J2QFIhSg+gc1HcsaMKWUk6QpHG3DeW1HKrEgiZsLeMESRr30yRsLRMdgvqOBe3o6i/BTQ9vV9E5yJa8XeHqjIbLonNQH7GgHd2m7zFvuOgQZHtzh2PT96JDUB+xoB3dR9/jjmTRIcj27kzGR0dEh6A+YkE7tMP1iPGDl4voHGR7Xi6I8kVxvegc1BcsaIe24RDuu0F0CLKX+1Kw4ZDoENQXLGjH1S5hZyUmx4nOQfaSGYcdlZyXQ01Y0I5rWzmmxsOJn09xGE46TIlHToXoHNRrLGjH9cFB3JciOgTZ18IUfHBQdAjqNRa0g2pqwZmLSAoUnYPsKykQZ5rR3Co6B/UOC9pBffw95o8QHYJE4AXRKsKCdlAffY87efmzQ7ojGX8vFh2CeocF7YjKz8HXDX5uonOQCAHu8HVD+TnROagXWNCOaMMh3MvLnx3YvTfgb4dFh6BeYEE7HEnC1+W4mQt4O7CZg7G1DBIviFY8FrTDya3GuEjo+cg7ML0Txg5CHhcqVDw+TR3OBwexiJc/O7z7eEG0GrCgHZCBV9IAABPiSURBVMsVI44akBwqOgeJdkMYSgy4YhSdg7rFgnYsnx1FNlfvJgBA9jB8znWwlI0F7Vj+7zD+3/WiQ5Ay3HU9/o8XRCsbC9qB1F8CgFAv0TlIGUK90CH9+28FKRML2oG8uZeXP9OP3HsD3torOgRZx4J2FE0t+KYcszkATZ3cPgxby9DUIjoHWcGCdhR/3YcHRkHH2Z+pE50OD4zG/+wTnYOsYEE7hKtGbPoed40UnYOU5+6R+KQELe2ic5AlLGiHsOEQ7kyGCx9t6sLFCXOHc61CheJTVvs6JKwrwoOjRecgpXpoDN49gA5OzaE8LGjt+6QE0xLg5SI6BymVlwumxmNzqegc1AULWvv+9B0eTxUdgpRtaRr+9J3oENQFC1rjtlVgZBiCPETnIGUL8kByKBf8VhwWtMa9locV40WHIDVYlo5X80SHoB9jQWvZvjMI8EC0r+gcpAZx/gj0wL4zonNQJyxoLXt5D1by5TP12srxeCVXdAjqhAWtWccbccWI6zn1M/Xa9aG41IrjjaJz0H+woDXr1Vwsv0l0CFKb5TdhNUeiFYMFrU31l3D0LCbEiM5BajMxFiUNnINUKVjQ2vSrHDyVIToEqdMvJ+DX20WHIAAsaE369iQuteHmRNE5SJ1uTkRzC749KToHsaC1p60Dv8rB6ptF5yA1W30zfpWDtg7RORweC1prXs/H/BEI47pWNADh3pg3AmsKROdweCxoTalswpajWDJOdA5Sv0fH4bOjqL4gOodjY0FrysqvsfpmOHPZFBowZx1emYHlW0XncGwsaO344jj83TF2kOgcpBVjB8HXDV8eF53DgbGgNeKqEat24vdTRecgbfnDNKzaiatG0TkcFQtaI/6wGw+MQiCnFSVZBXrg/lF4aY/oHI6KBa0FR88irxqLUkTnIC1anIIdJ3HsrOgcDokFrQXLvsKrM6Dje4NkAzod1vwEv/hKdA6HxIJWvb8fQXwAkjlrHdlMciji/LHxiOgcjocFrW77z2BNAd8bJJv7wzT8MR/7OZ2/fbGgVayyCQ//A5vmw9tVdBTSOm9XfHwHHv4HTjWJjuJIWNBq1dSC+R/hr7ciwlt0FHIMEd74y62Y9xGaWkRHcRgsaFVql3DXx3h2ElLCRUchR3JjOJ7KwF0fo10SHcUxsKBVaeXXmBCLrCTROcjxzB6GjBg8+bXoHI6BBa0+bxfiYiue5GqwJMhTGbjQgrcLRedwACxoldlWgc2l+NNPRecgx/ZmFjaXYluF6Bxax4JWk4rzeOobfDAbLnzcSCgXJ3wwG099g4rzoqNoGp/oqnGgFgs/xfvZnIyfFCHMC+9nY+GnKKoVHUW7WNDq8F4RfpWDT+/E8BDRUYj+Y3gINt+BZ3LwXpHoKBrFgla6dglL/4UDtdhyJwLcRach+rFAD2y5EwdqsfRfvPZOfixoRTt/FbM/xIgQvD4Tej5WpEh6J7w+EyNCMPtDnL8qOo228EmvXCUGZP0vVozHg6NFRyHqyYOjsWI8sv4XpQbRUTSEBa1QORVY9Cn+51ZMiBEdhah3JsTgr7di4afI4eV3MmFBK06JAXM24r0ifHEX3xIklRkRgi/uwroizN3Il9Iy0IsOQD8oO4entwHA76ZiaJDoNET9EuSB9bNx9Cye3Q4Av5uKxADRmVSLBa0ItRfx3Lc4eR4vZmJcpOg0RAM2NAh/n4vvavDIF4jzxwuZvH6/P1jQgp08j78dxhfHsOwm3D6My1aRpoyLxL/uwsclyP4QtwzBXdcjzl90JlVhQYtRfg6bvsfWMkT6Yv4IPDmeV9GRNul0mDsc2cPwVRme+xY1FzAjEXOHI4HjHr0gf0FLkvToo4+WlJR4eHisW7cuNDTU4vaQkBCLu2lYh4SjZ/FpKf51ApE+uCMZj6fCnf9EkgPQOyErCVlJuGrEv07gVzmoacbMwcgehqFBcOJ/HK2Qvx62b9/e0NCQk5Pz7rvvrl69+g9/+IPF7TNmzLC4m5a0tKO4HgfOoKgWJQa0dyAxED8ZjCfuhgd7mRySux7Zw5A9DFeM+OIYXslFWSOcnXBdMFLCcWMEkkPh5iw6pWLI3xO7d+9OT08HkJaW9t5771nb7unpaXE31WlqweU21FzA6WbUNONMM6ouoPYirhrh4oTkUKSE48HRuC6EU9AR/cBDj7nDMXc4ALR1oKQBRbX42yEU16OtA+56hHsj2hcRPoj0wSAfRPrC0wV+bqJz25f8BW0wGJKTkwHExsYaDAZr263tBkCSpKVLl7a0/LDw2bFjx2bOnNnUNNDlKq+2657Z4d7LGQOuGHUtRgCQgKaWH/4PJknQ6SBJAKDTwcdV8nRBmGfHIB8pwqtjfHhHZJIU5tnh+uNXAZebB5idSMti3REbh1lx//62tR11l51qmnV1l5xKap22HdfVXXa63IbmVp35qWd6Jpr5uUmm79z08ND36knurMNvJ111dx7oHCINDQ15eXl79uwxb3Fzc3v99dd1A37TX/6CDggIqKysBFBZWRkYGGhtu7XdAOh0umXLlrW3t5u3bN682dPT09t7oMujekl4auK/u7VHni5w0wOADvDvbpYi02PAl8dEcgr0w3URfdj//FWd6ZndYsTltl79ik6HIF+vgV865evrO2HChNmzZ5u3ODs7D7ydYYuCnjhx4tq1awEUFhZmZGRY225tN5PY2NjO34aEhBiNRmdnGYamBvMDIERaFCTuOms3NzcvL6+EhATZjyx/QWdmZm7ZsiUrK0uv17/zzjtHjhy555579u/ff832wMDAzt/KHoOISO10Ui//wy/UBx98YDQaFy9eLDoIEdG13n33Xb1ef++998p+ZI6cEhEpFAuaiEihWNBERArFgiYiUigWNBGRQrGgiYgUigVNRKRQLGgiIoViQRMRKRQLmohIoVjQREQKxYImIlIoFjQRkUKxoImIFIoFTUSkUCxoIiKFYkETESkUC5qISKFY0ERECiX/orE2smPHDqPRKDpFd06ePBkTE+PkpKl/8y5dunTx4sWwsDDRQWRWXl5uizWYxaqrq/P29vbyEre6tQ10dHScOnUqLi5OdJDu7NmzZ+rUqbY4sjoKetq0aR4eHqJT9OC9995bvHixn5+f6CByqqqqOnHixLBhw0QHkdk//vGP5557TnQKmW3fvn3w4MFRUVGig8ipqanp22+//fnPfy46SHduueWWjIwMWxxZHQU9aNCgefPmiU7Rg02bNt12222hoaGig8jJ39/f3d1d+Xd+X7311lvau1HHjh0bN27c9OnTRQeRU319/ddff629B6uXNPX/cSIiLWFBExEpFAuaiEihWNCycXNzc3FxEZ1CZq6urq6urqJTyE9jlzqYaPLB0uv1bm5uolMIo5MkSXQGIiKygK+giYgUigVNRKRQLGgiIoViQRMRKRQLWh6rVq3atGkTAEmSlixZMmXKlKysrPr6etG5BqqlpcXf3z8lJSUlJeXll18WHUcGGnuATLT3MEG7z6k+YUEPVHt7+6RJk55//nnTt9u3b29oaMjJyZkzZ87q1auFRpNBeXn5rFmzioqKioqKVqxYITqODDT2AJlo7GHS9nOqT1jQA+Xk5LRt27aVK1eavt29e3d6ejqAtLS03NxcodFkcPz48dLS0uzs7Pnz51dVVYmOIwONPUAmGnuYtP2c6hMW9EDpdDq9Xm+eZdRgMMTGxgKIjY01GAxCo8kgJCRk+fLln3766e233/7oo4+KjiMDjT1AJhp7mLT9nOoTFnQ/rVu3bu7cuatWrbpme0BAQGVlJYDKysrAwEAR0WRgvnXp6emmicRmzZp1+PBh0blkoI0H6Brae5g60+RD1kvqmG5UgRYtWrRo0aKu2ydOnLh27VoAhYWFNpoi1g7Mt+6ll17y8fF55JFH8vPzk5OTReeSgTYeoGto72HqTJMPWS+xoGWWmZm5ZcuWrKwsvV7/zjvviI4zUA888MD999//4Ycfuru7//nPfxYdRwYae4BMtPcwdabJh6yXOBcHEZFCcQyaiEihWNBERArFgiYiUigWNBGRQrGgiYgUigVNRKRQLGgiIoViQRMRKRQLmohIoVjQREQKxYImIlIoFrRGREVF6f7Dx8cnKyvr9OnT3ey/YsWKgIAAhSwgdP78eX9/fwCFhYVjxoyxtpterzcajZ23mPffsGHD3Xff3cvTmY7T/bn6p8d71RYnNd97pD0saO348ssvz50719jYuH///gsXLjzzzDPd7Lx27dqSkpLQ0FC7xeuN+Pj4F154wXb7y/W71vR4r9ripKRhLGjt8PHx8ff3DwgISEpKuvvuu8vLy03bd+3adeONN3p5ec2cObOmpgbA7Nmzm5qaxo0b19DQ0PWnxcXFkydPXrVq1ciRIy3+emlpaUZGxiuvvBIZGRkfH5+Tk2M60caNG5OSkoKCgn72s5+1tLRY/N3O3njjjejo6Ojo6Hfffde0paKi4tlnnwVgNBp/9rOfBQQEBAcHv/jiiwBmzJjR3t6emJhYUFBgjmfeH8CVK1fuvPNOPz+/tLS04uJiAPn5+Wlpaaafmr82H6e4uNj8u5988snQoUP9/PzmzJnT0NDQzW006/orne9V0z5db0XnwO+9915cXFxcXNz7778fFxfXzUkt3o1d7z2zrue1eFeUlpaOHz9++fLlwcHBGRkZeXl5Y8eO9fHxeeKJJywehASQSBMiIyN37dpl+rqmpiY7O/u3v/2tJEkGgyEoKOizzz5rbGxcsmTJ5MmTTfv4+fk1Nzdb/Onhw4f9/PwWLlxYXFxscYeSkhIvL6/f//73ly5dWrlyZXp6uiRJR48eDQoKys3NLSsrGz169Nq1a62d2mTnzp0BAQE7duyorq6eMmWKn5+fJEl79+4dPXq0JEkbN24cOnRoRUXF/v373dzcTpw4IUmSs7NzW1tb53jm/devXw/gL3/5S319/cqVK6+77jqj0ZiXl5eammo6XeevTccx/255ebmfn9/WrVvPnj27cOHCO+64w9ptNLP4K+Z71bxb11thPumhQ4eCg4MLCgpqamoyMjJiY2OtndTi3Wjx3uvmvBbvipKSEicnp7/97W9nz54dPXp0aGjoyZMn8/LyANTX11t8CMjOWNAaERkZ6eXl5efn5+vrCyAtLc1oNEqS9N57782ZM8e0z5UrVzw9PU3bTVVi8aeHDx92dXW9evWqtV8vKSnx9fVta2uTJOnw4cNDhw6VJOnFF1987LHHTHsWFRXt2LHD2qlNli5d+tRTT5m+3rNnT9eCTkhIKCgo6OjoaGhoaGlpkToVtDle54IeNWqU6Witra2BgYElJSW9LOjXXnvtvvvuM/2ovr7excXF2m00s/grkqWCvuZWmE/69NNPr1ixwrTbli1bzAXd9aQW70aL914357VW0JGRkaaNTz755EMPPWT6OjY29vjx4xYfArIzrqiiHevWrRs7diwAg8GwYMGCDRs23HfffVVVVVu3bjX9DxqAq6trfX19RESE6VuLPwUQHR3t5ubWzQ7h4eF6vR6A6U8A1dXVSUlJpq9vuOEGADt37uzm1LW1tdOmTTN9nZCQcM1tuf322y9cuPDggw/W1dUtWbJk2bJlnX9qjtdZfHy86QsXF5e4uLja2lp3d3fzTyXrC1PU1taaQ4aEhLi6uprGKLrexu5/JTw8vPe3oqamxvRgAYiJiTFv73pSiw/BQO69zneFt7e36Qu9Xm/Obzp19wch++AYtHZERESYxjTHjBkzZ86cAwcOmDZOnz795MmTJ0+eLCsr+/rrrzv3iLWfmtvB2g46ne6as4eFhVVXV5u+zsvLW79+ffenHjRoUFlZmenrioqKa45WXl4+ZcqUoqKigoKCzz///JqFjro2ZueDtLW1VVZWRkdHAzBf9WHO1lV4eLhpTVIAppeKwcHBFm9jj7/S+1sRERFx6tQp09dVVVXm7V1PavFu7Me915u7opfhyW5Y0NoUHh5uetpnZWXt2rXryy+/NBgMTz311BNPPNG5Arr/aW92MJszZ8769esLCgrKy8ufeOIJg8HQ/e/OnTv3L3/5y65du06fPv3ss89ec9jPPvtswYIFdXV17e3tLS0tHh4epu0XL160dpOLiopMA9/PPPNMfHx8QkKCn5/fwYMHi4qKzp49++abb3beufNxbr311k8++WTbtm3nzp1bvnz57NmzLf4D0Fkvf8XarTDd/HXr1hUWFp45c+bVV1/t5lwW78a+3nvd3BXWdBOe7Ef0GAvJo/ObhJIkffHFF6GhoU1NTZIkffXVVyNGjPD09MzMzCwrKzPtYB4t7frTa4Zcu+5QUlJi3qHz1++//358fLyvr++iRYtMQ5YWT232xhtvREVFRUZGvvvuu6bBUPMQbXNzc3Z2tpeXV2Bg4MMPP9za2ipJ0vz58318fPLz881n7DwG/eCDD2ZnZ/v4+EyYMMH0jlZHR8djjz3m7e19/fXXf/TRR+ZBWNNxvv32W9PvSpK0adOmIUOG+Pj4ZGdn19XVdXMbzbr+itRlDLrrrTAHliTp7bffjoiIGDJkyNtvvz1ixIhuTmrxbux673VzXot3RedTPPPMM88995zp68TExOPHj1t8CMjOuGgskQClpaV1dXWTJk0CsHXr1t///vfbt28XHYoUh0McRAKcO3duwYIF9fX1V65ceeONN37605+KTkRKxIImEiA9Pf3xxx+/8cYbk5KSIiIilixZIjoRKRGHOIiIFIqvoImIFIoFTUSkUCxoIiKFYkETESkUC5qISKFY0ERECsWCJiJSKBY0EZFC/X/IVrbXH+rmNwAAAABJRU5ErkJggg=="
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The _p-value_ for a test of $H_0:\\mu_{80}=\\mu_{90}$ versus $H_a:\\mu_{80}<\\mu_{90}$ is"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(refd >= sum(Secchi$diff))/length(refd)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "[1] 1.192093e-05"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Sampling from the reference distribution\n",
      "\n",
      "A more common approach to working with the reference distribution is to generate a reasonably large random sample from the distribution.  Generation of one instance of a sum with randomly allocated signs is often written using a random sample from a binomial distribution with size 1 and probability of success, 0.5"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "set.seed(1234321)\n",
      "rbinom(22,1,0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        " [1] 1 0 1 1 1 0 1 1 1 0 0 1 1 1 1 0 1 0 1 0 0 1"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "but it is somewhat cleaner to simply generate a sample from a uniform (0,1) distribution and compare them to 0.5"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "set.seed(1234321)\n",
      "c(-1,1)[1 + (runif(22) > 0.5)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        " [1]  1 -1  1  1  1 -1  1  1  1 -1 -1  1  1  1  1 -1  1 -1  1 -1 -1  1"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can vectorize this calculation to obtain, say, 100,000 samples from the distribution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ns <- 100000L\n",
      "system.time(refsamp <- \n",
      "            colSums(Secchi$diff * matrix(c(-1,1)[1 + (runif(22*ns)>0.5)],\n",
      "                                         nrow=22)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "   user  system elapsed \n",
        "  0.163   0.000   0.163 "
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "with the corresponding density"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "densityplot(~refsamp, plot.points=FALSE,\n",
      "            xlab=\"Sampled Reference distribution of signed sums\",ref=TRUE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": []
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3deVxU9eI+8GcAWYVhFRDZ3VJU3EFxAZcsKjWX8t4Wq291y7ZbWn3rlt6yW/3qdq/1bb2WlXYzs8Vb6dVSVBRRESlRcAFkVWFUXJBthvP7Y2pCmGGdmc85Z573H73gcOacZzjD0/iZcz5HI0kSiIhIfpxEByAiIvNY0EREMsWCJiKSKRY0EZFMsaCJiGSKBU1EJFMsaCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZIoFTUQkUyxoIiKZYkETEckUC5qISKZY0EREMsWCJiKSKRa0Fej1etERbEXFTw2qfnYGg0GSJNEpbEXFB64FFnR3GQyGuXPnik5hK6+//vquXbtEp7CVWbNmiY5gKx988MHGjRtFp7AVFR+4FljQ3dXU1HTlyhXRKWyloaGhoaFBdApbqampER3BVnjg1IEFTUQkUyxoIiKZYkETEcmUi+gAHVJeXv7DDz+ITmGewWAoKyv74IMPRAexiaysrOrq6hMnTogOYhMVFRVqPXAZGRn5+flVVVWig9iEDA9campqWFiY1TerUcS5OJ9++mlmZmZycrLoIOYVFBTExsaKTmETp06d8vHx8fLyEh3EJlR84KqqqlxdXbVareggNiG3A5eWlpaQkHDHHXdYfcvKeAcNYNSoUfPmzROdgoiopUuXLtloyxyDJiKSKRY0EZFMsaCJiGSKBU1EJFMsaCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZIoFTUQkUyxoIiKZYkGTQ6uuw+UG0SGILHARHYDI3ur0yCjFT4XIKIWbCzTA+vno6So6FlErLGhyFKcu49OfkVaEOj3GhWNKDJ6fBHcXpJdg/pdYPx+ePURHJLoaC5ocwpka3PwFHk/E/SPh637VjyZEoHEcbl2PL+bBg38QJCd8PZL6nanBzM/xwY0YGmx+hZRouDrj1vVYO5cdTTLCDwlJ5c7XYc4X+OcMi+1slBSBh8bgli9Rp7dXMqL2sKBJzS434OYv8EIyEvq0v/K0GCwag1vXo95g+2REHcCCJtWq02PuOvw5ASnRHX3ItbFYGI87vkFjky2TEXUMC5rUySBh4be4JQ43DejcA2cNxPRYLE2zTSyizmBBkwoZJCxYjxl9cVd8Vx5+z3D8fAYF560di6iTWNCkQv9vN4YEY2GX2tlo2WS8sMN6gYi6hAVNanO5Ad8fw18mdGsjo3tDkpBVYaVMRF3Cgia1WZWDO4ZBo+nudl5MwbPbrBGIqKtY0KQqBgn/PtStwQ2TSC2GBuOH41bYFFHXsKBJVb7Jw4y+cHO2ztaenYCX06HnKXckCAuaVOXt/Xh4jNW25uuOWQPxyc9W2yBRp7CgST12l+KaQPh7WHObj4zFxzmoabTmNok6iAVN6vGPPfhzopW36eqM+0fin5lW3ixRR7CgSSWOnwOAfv7W3/Ifh+DHApypsf6WidrGgiaV+Gem9d8+G2k0WDYZL+20ycaJ2sCCJjXQXcFRHcaH22r7k6NwsvrXN+lEdsOCJjVYmY37Rtp2F3+ZiNczbLsLohZY0KR4dXp8dwxzB9l2L2PCcFSH6jrb7oWoORY0Kd5nh/DHIXDq9rXd7bptKP59yOZ7ITJhQZOyNUn4OMc613a3a8EQfJ5rjx0RGbGgSdn+ewIp0fDsYY99efXAkF7YW26PfRGBBU1Ktyqni7Pyd83dw/Fhtv12Rw6OBU0K1mBA+UVE+dpvj6N64+hZnOdHhWQXLGhSsIxSJEXYe6f8qJDshgVNCrbpBK7rZ++dLojD5yxosgsWNCmYkHfQPV0R1wuZZfbeLzkgFjQpVckFBHmih4iX8L0j8cEBAfslR8OCJqXaXCBgfMNoZCiOn+NHhWRzLGhSqv+ewHV9he39dn5USLbHgiZFajDgzGX08REW4A9D+FEh2RwLmhRpt4iPB5vr6YrB/KiQbIwFTYq06biwAWiTe0fwo0KyLRY0KdKeMhtOz99Bo3rjGK8qJFtiQZPylF9CLy+4yODFu2AIvuD8dmQzMniNE3XSf09ghrjzN5r74xB8my86BKkXC5qU54djuKG/6BAAAF93ALhQLzoHqRQLmhSmwQDdFYT2FJ3jN1NjsK1IdAhSKRY0KYzwE+xauK4fNh0XHYJUigVNCiOHE+yaGxyEI1WQJNE5SI1Y0KQwe8qQ2Ed0iKsN7oXDVaJDkBqxoElJii8gtKcsTrBr7rq+2HRCdAhSI5m90onatFnEDP3tmhKDrYWiQ5AasaBJScTOYGeJtysknmxHNsCCJsWoN0B3BSGyOcGuuakxSOPJdmRtLGhSjN0lmBgpOoQFHIYmW2BBk2JsLsD0WNEhLBgchNxK0SFIdVjQpBh7SjFO9Ax2lmg0GBzEk+3IyljQpAxnauDrLrsT7JrjJYVkdTJ+vRM1s/0kUqJFh2jT1Bhs5eeEZFUsaFKGrYWYGiM6RJu8XdEk4VKD6BykIixoUoYjVRgcJDpEe6ZEc2Y7siYWNClAUTUitNBoROdoD4ehybpY0KQAWwsxRd7jG0ZDevFkO7ImFjQpwNYiuQ9Am1wThCM82Y6shAVNcidJKK5GpFZ0jo7hJYVkRSxokrvcKgwJFh2iw6ZyZjuyHhY0yd1PhZgi7zOgm/Nxg76JJ9uRdbCgSe7SipTxCaHJFM5sR1bCgiZZ0zfhfB0CPETn6IzpsdhcIDoEqQILmmQtqwKjeosO0Unxwcg5LToEqQILmmTtx0L5TjFqiUaDgYHI04nOQcrHgiZZSy+W7yT9bbiuLy8pJCuwfkFLkrRo0aKUlJTU1NTKypaXVS1fvnz9+vXtrkYEoKYRBglePUTn6LzpsfiJJ9tRt1m/oNPS0qqqqrZt2zZnzpw33njDtNxgMEyaNGnZsmVtr0ZksqsEEyJEh+gSHzc0NuEyT7aj7rF+Qe/atSsxMRFAQkJCRkbG73tyctq6deuTTz7Z9mpEJvKfYrQNKdFIOyk6BCmci9W3qNPp4uLiAERGRup0v39QotFoXFxcnJyc2l4NgCRJy5Ytq6///S72ubm5CQkJVo9KMre3HC9NER2iq67riw8O4Mb+onOQ7VVUVGRmZh45csS0xM3NbdmyZZpuT8Bo/YL28/MrLi4GUFxc7O/v34XVNBrN7NmzDQaDaYmnp6efn5/Vo5Kc6a7A2xU9FPsx9rBg/HxGdAiyCz8/vzFjxqSmppqWODs7d7+dYYuCnjhx4sqVKwFkZWUlJSV1bbX4+Pjm3x4+fFiv11s9KsnZ9pNIVs4V3q1pNBgQgHwdBgaKjkI25uHhERUVNXLkSKtv2frvT5KTk4OCglJTU7/66qvFixcfPnx4xIgR7a5m9RikdFuLlDQFh1nX9ePMdtQt1n8H7eTk9Oabb5q+DQwMzM7ONn27fPlys6sRtfDLGQxTziR2Zk2LwR++wp/56Ql1lWJH+EjVii+gj48C7nHVNl93NBh4sh11HQua5Gib8sc3jFKisf2k6BCkWCxokqO0ImV/QmhybV9s4cx21FUsaJIdScKJc+hn8RRNJRkRggOnRIcgxWJBk+z8fAbDQkSHsBKNBv0DcPSs6BykTCxokp0tBZim2Cu8W+PMdtRlLGiSnW1FSFHFALTRtFj8yJntqEtY0CQvNY1oMMDXXXQO6/FzR52eJ9tRV7CgSV52FmNSlOgQ1pYSjR3FokOQArGgSV5+LFDePa7axWFo6hoWNMnL/gqMCRMdwtqGh+AgbyNLnceCJhkpvYhgLzgr/Arv1jQa9PPHMZ5sR53EgiYZ+akQ01Q3vmHEme2oC1jQJCOKvsdV26bG8Day1GksaJKLJgmF5xGr0jvnBHjgXC30TaJzkKKwoEkusk9hZG/RIWxpWDAOVYoOQYrCgia5UNkV3q2Nj0BGqegQpCgsaJILpd+EsF3jwrG7RHQIUhQWNMnC5QZIgLer6By2FO2LomrRIUhRWNAkC9tPYlKk6BC218cHZRdFhyDlYEGTLGwpwLV9RYewvfHhHIamTmBBkywcPI2RoaJD2N44FjR1BguaxCuqRpg3nFR3hXdrw0M5KQd1AguaxFPxFd4t9HBCDyfODU0dxYIm8X4swLWOUdAAxvZBVoXoEKQQLGgSzCCh7CL6+IjOYS/jwrGbw9DUMSxoEixLjRNAt2FcOPawoKljWNAk2OYTKryFShv83FFdhyZJdA5SAhY0Cba7FBMd4BKV5gYG4ign76cOYEGTSPomNBrQU9VXeLc2IRLpvIcsdQALmkTKrURcL9Eh7I6Xq1AHsaBJpH3lGO1InxAa9fXDiXOiQ5ASsKBJpP0VGK3qSfrN0mjQywuVNaJzkOyxoEmkY2cxIEB0CBE4ykEdwYImYWoa4eECjQNMwdEaC5o6ggVNwmSfwggHmMHOrJG9ceCU6BAkeyxoEsYxPyE0cnMGgDq96BwkbyxoEmZ/uWNd5N3CKL6JpvawoEmY8ksI8xYdQhwOQ1O7WNAkRmUNAj1FhxCKt7+idrGgSQzHPAO6uUBPVNZA4qxJZBkLmsRw8AFoowEBOHFedAiSMRY0iZFV4bincJiMC8fuEtEhSMZY0CTGhXpo3USHEG18BIehqS0saBLgZDUitaJDyMCAAOTpRIcgGWNBkwCZZRjbR3QIGXDSwN+DsyaRRSxoEmC/g92HsA082Y7awIImAXJOIz5EdAh5SIrgTb7JIhY02Zu+CfqmXyejoJG9kVUhOgTJFQua7O1wlSPe5soSN2c4aXClUXQOkiUWNNnbvnJHv4awhbFh2M830WQOC5rsjdcQtjA+Art4uQqZw4Ime8vXYWCg6BByMj4ce/g5IZnDgia7qmmEmwucHPI2V5b4uqO6DgbOmkStsKDJrrJPYaSj3uaqDUOCcbhSdAiSHxY02ZUj3+aqDePDOQxNZrCgya74CaFZvFyFzGJBk12VXUS4j+gQ8hPli+Jq0SFIfljQZD+VNQhw7NtctSFci5ILokOQzLCgyX6yHP42V20YH85RDmqJBU32s48D0JaNj+DdVaglFjTZT1YFRvEdtAVDg3GIZ9rR1VjQZCeShHO18PcQnUOunDXw7IGL9aJzkJywoMlOjp9D/wDRIeRtHIeh6WosaLKTjFIkhosOIW9JHIamq7GgyU72lCGR9yFs05gw7CsXHYLkhAVNdnKE8/S3x6sHGgxoMIjOQbLBgiZ7uFAPzx6cxK59I3sj+5ToECQbLGiyh33lGMszoDuAsyZRcyxosoc9/ISwY5IikMETOeg3LGiyh73l/ISwQ3p54fRlSJy8nwCwoMkOmiScr4Wvu+gcCjEwEMfOiQ5B8sCCJpvL1+GaINEhlIP3kCUTFjTZHM+A7hRerkImLGiyub1lGMuC7rD+/jh2VnQIkgcWNNnckSrEcYijwzQaBHribK3oHCQDLGiyrXO18HGDhpeodMboMOznNd/EgiZb21fO8Y1OGxOG/RWiQ5AMsKDJtvaUYRwvUemk0b2RxYKmNgr6kUce2blzp8HAiVuoW/aW8SLvTvN1x9krokOQDFgsaD8/v4cffjgsLOzBBx/ctm2bXq+3ZyxSB4OESw3wcROdQ4GifHmTb7Jc0H/9619//vnnjIyMvn37Llu2rE+fPvfdd9+WLVsaGxvtmY8UjVOMdtlozg1N7Y5B+/v7h4eHx8bGNjQ0ZGRkLFu2LDo6esOGDfYJR0qXUcpLVLqInxMS2ijo1157bfLkyX369Fm5cuWIESMOHDiQm5ubkZGxZs2aP/3pT/aMSMrFSey6bHgIck6LDkGiuVj6QVZW1iOPPDJt2jRvb2/jkpqaGi8vr9GjR7/zzjv2ikfKduIc+vuLDqFM7i6o08MgwZmnkDswM++g9Xq9Xq/PzMy86aabPDw8jN9WV1eHhoYC8PLymj17tt1zkvLorsDPg5eodN01gcjXiQ5BQpl5B+3u7g7AYDAYvzCZO3eunUKRKmSWIYED0N1gvJ5wMK+Sd2AW30FPmzZNf7W1a9faPx8pFyex6yZ+TkgWPyTcsmWLPXOQ+uznRd7dMygIhytFhyChzAxxjBo16oUXXnj++edb/ygrK8v2kUgN9E240givHqJzKJmzBj2cUW+Am7PoKCSImYJ+7733oqOj33vvPfunIdU4VImhwaJDKF98CHJO81p5x2X+HTSAgICAEydOREREGAyGDz/80NPT8/bbb7d7PFIqngFtFaN7Y185C9pxWRyDfuGFF+Li4i5evPj666+vXr16xYoVjz76qD2TkaJxEjur4MTQDs7ihSorVqzIzMwMCAh455139u7dq9frx4wZw0tUqIMKziHGV3QI5YvS4mS16BAkjsV30AaDwdfXd//+/cHBwREREZ6eng0NDfZMRspVcQm9vXmJihVoNNC6o7pOdA4SxOI76FtvvXXGjBmNjY3PPvtsUVHRggULpk2bZs9kpFzpJZgQKTqEWowMxYFTmBItOgeJYLGg/+///u+bb74BcPPNNxcWFs6bN+/++++3YzBSsPRi3DNCdAi1GBOGfeUsaAdlsaBdXFzmzZtn/Lpfv35PPPGEvSKR4v1yBsN4jp2VjAnDhwdFhyBBLBb01q1bn3vuuXPnzjVfmJ+fb/tIpGxna6F1hxMHoK0k0BNVNaJDkCAWC/ruu+9esGDBbbfd5uJicR2i1naXIClCdAh16e396+eu5Ggslm9jY+PSpUs9PDzsmYZUIL0EN18jOoS6jA7D/grMHCA6B9mdxdPsHn/88TfffJN39abOOlCBUb1Fh1CX0b15uYqDsvgO+ttvv83Jyfnb3/4WEhKi+e2MVo5BU9suNcDVGT3audUldc6o3nhll+gQJILFgl65cqU9c5A6cAoOW/DsgcsNkCRe++NwLBb0wIEDARgMhsrKyuZvoonakF6CFJ6xawP9A3D8HPoHiM5B9mXx36Ll5eXJyck+Pj6DBg06cODAhAkTioqK7JmMlGgvb3NlG6PDsI/D0I7HYkHfddddcXFxZ8+e1Wq18fHxCQkJ9957rz2TkeIY70LtwdMybYC3v3JMFv+Ydu3atW7dOuN9Y11cXJ566qnISE6vQG3ZX8GZi20lrhdyefsrx2PxHXS/fv127fr9k+O9e/fGxMTYJRIp1e4SjOclKrZhPDGmnme9OhiL76DffPPNOXPmTJ48+dy5c3PmzElPT1+zZo09k5HipJfgwdGiQ6jX2DBklmES/x3rSCwW9KRJk44ePfrdd9/Fx8eHhoa+/fbbISEh9kxGyqJvQk0DfNxE51CvCZHYWcyCdixtfaATEBCwcOFCeyUhZcs5jRGhokOo2vhwrMgUHYLsy/wYdFZW1ty5c2NiYtzd3WNjY+fPn5+dnW3nZKQsO4s5Sb9t+bjhSiMaOAztSMwU9LZt2yZPnty/f/81a9bk5uauXr06NjZ24sSJO3bssH8+UopdJZjATwhtbHQYsk+JDkF2ZGaI45lnnnn11VcXLVpk/LZv377jxo3r3bv3//7v/2ZkZNg3HilDk4SqKwj0FJ1D7SZEYEcxLwVyIGbeQefk5Nx4440tFt50000c5SBLjlQhrpfoEA5gYiR2l4gOQXZkpqDr6+t9fHxaLNRqtfX19XaJRMqzs5jjG/bg74FztTBIonOQvZg/i+PQoUPe3lfdv+HSpUt2yUOKlF6C16eLDuEY4kPwyxkM5ymvjsFMQWu12tZDHMblts9DilR+EWG8IZNdTIzEjpMsaEdhpqCrq6vtn4OUq+A8Yv1Fh3AYEyOxaCMeSxCdg+yCt76g7uIAtD2F9MTpy5A4DO0YWNDUXenFmMhLVOxoUBCO6ESHILtgQVN3FZxHXw5x2FFSBHbxZDvHwIKmbqm4hNCeokM4mImR2FksOgTZBQuaumVrEaZwnnD7ivZFwTnRIcguWNDULWlFSI4SHcLx9A/ACXa0A2BBU7cUcgBaBI5yOAgWNHVdUTUiePWSCCxoB8GCpq5LK0JytOgQDolDHA6CBU1dl3aSA9DC9PFByQXRIcjGrF/QkiQtWrQoJSUlNTW1srLS0vL6+npfX9/4+Pj4+PjXXnvN6jHIDgrPI8pXdAhHNSES6TwbWu2sX9BpaWlVVVXbtm2bM2fOG2+8YWl5YWHhzJkzc3JycnJylixZYvUYZGvHz2FAgOgQDmwSh6EdgPULeteuXYmJiQASEhKa34GlxfLjx4/n5+fPmjVr/vz5paWlVo9BtsYBaLEGB+FwZfurkaK1dVfvrtHpdHFxcQAiIyN1Op2l5UFBQYsXL543b97atWsfeuihDRs2mNaUJGnJkiXNZ6A+duzY1KlTrR6VuiPtJF6bJjqEA9NoEOSFMzUI9hIdxeGdOnXqp59+2r17t2mJt7f3a6+9ptFourll6xe0n59fcXExgOLiYn9/f0vLje+mAcycOfOZZ55pvgWNRrNw4cLmN3DZuHFjQAD/OS0jkoTSC+jT8sY7ZFfGk+3mDRKdw+H5+/snJyenpqaalnh5eXW/nWGLgp44ceLKlSsBZGVlJSUlWVr+6quvent7P/jgg5mZmcZ31s21WHL48GG9Xm/1qNRleToMChIdwuFNjMQnOSxo8dzc3KKiokaOHGn1LVu/oJOTkzds2JCamuri4vLhhx8ePnz49ttvz87ObrHcycnpnnvuWbt2rbu7+7vvvmv1GGRTaSc5AC1efAgePyM6BNmS9QvaycnpzTffNH0bGBhovB14i+UAvvnmG6vvnexjZzHeuFZ0CIfnrEFPV5yvg5+76ChkG7xQhTqtSeJNCOUiJRpbC0WHIJthQVOn5VZiaLDoEAQAmBqDrUWiQ5DNsKCp0zgALR9DeyGXZ0OrFwuaOm37SUyOEh2CAAAaDcJ9UFQtOgfZBguaOscgQXcFQZ6ic9BvpsTgJw5DqxQLmjon5zSGh4gOQc1Mi+HnhKrFgqbO4RQcchOhRfEFNEmic5ANsKCpc3YUY1Kk6BB0teEh+JlXrKgRC5o6obEJF+rg7yE6B12Nw9BqxYKmTsg+hVG9RYegVqZEY/tJ0SHIBljQ1AnbOAAtS77uqGlAvUF0DrI2FjR1QnoxJnIAWpbGhSOD971QHRY0dVRjE2oaoXUTnYPMSY5GGq/5Vh0WNHXU3jKMDRMdgiyYEMF7yKoQC5o6avtJTIoSHYIscHeBuwvO14nOQVbFgqaO2lWCCRGiQ5Blk6N4LofasKCpQ+oNqNPDhwPQMjaV13yrDguaOmR/Ocb2ER2C2sTrCdWHBU0dsv0kr/CWOycNenuj9KLoHGQ9LGjqkN2lmMCClr0p0bzmW1VY0NS+egNqG+HtKjoHtYfD0CrDgqb2cQBaKWL8UHgeEqceVQsWNLWP97hSkCHByK0SHYKshAVN7dtdyjOgFYPD0GrCgqZ2GM+A7skBaIWYEsNJOdSDBU3t2F/OKTiUJMAD1XWcelQlWNDUjvQSnmCnMInh2FcuOgRZAwua2pHOmxAqzY398d1R0SHIGljQ1BYOQCtRYjgyy0SHIGtgQVNbeAa0EjlrEOaDomrROajbWNDUFk7BoVCp/bDxuOgQ1G0saGpLRimSeAa0As3oi80nRIegbmNBk0UcgFauQE9cbsCVRtE5qHtY0GQRB6AVLTkaaSdFh6DuYUGTRZyCQ9FS++GHY6JDUPewoMmijFKMDxcdgrpqeAgOnhYdgrqHBU3m1RtQb+AAtIJpNBgUhMOc2U7JWNBkHqfgUAGOcigdC5rM4wC0CkyLxY+celTJWNBk3p4yngGteN6ucNLgfJ3oHNRVLGgyo7EJ9Xp49hCdg7qNdylUNBY0mbGtCIk8f0MVruc130rGgiYzvjyMuYNEhyBrGByEfB2aeBtZZWJBU0v6JuRWYliw6BxkJaN6Y3+F6BDUJSxoamlHMc/fUJXrebKdYrGgqaWvjmAOxzdUZHIUtp8UHYK6hAVNVzFIOHAKo0JF5yDrcXeBrzsqLonOQZ3HgqarGOff0GhE5yCrur4f/svpoRWIBU1X4fiGKqX258l2isSCpt9JEvaUIZFzQKtOuA9OXUaDQXQO6iQWNP3uwCmMCIUTxzfUaEIE0ktEh6BOYkHT79YdxvzBokOQbcwfjC9yRYegTmJB0+92FGMi7+GtUiNCcagS+ibROagzWND0q4OnMSwYzhzfUK+UaGwrEh2COoMFTb/i+RuqN28QvjwiOgR1BguafrW1CCnRokOQLcWHIK8KjRzlUA4WNAHA4Sr0D0APvhzULpmjHIrCv0gCjOMb14gOQbY3bxC+PCw6BHUYC5oA4MdCXNtXdAiyvaHByNPxihXFYEETjp5FHx+4OYvOQXYxNYajHIrBgiZ8k8fxDQfCczkUhAVN+O8JXN9PdAiyl7heOMpRDoVgQTu60ovw8+ANvB1LcjSn8FcGFrSj+yQHdwwTHYLs69Y4rOW8HErAgnZokoTvjuGG/qJzkH0NDsKxsxzlUAAWtEPbVYrEPrw+xRFNi8VPhaJDUHv4p+nQPs7BwnjRIUiEW+N4LocCsKAd1+UGHDuL+BDROUiEAQEoOMdRDrljQTuur/Iwj9PXObBpsfiRoxzyxoJ2XJ8fwh+Hig5B4twymPNyyB0L2kEVnEdPVwR4iM5B4vQPwMlq1HOUQ8ZY0A7q34dwG98+O7wpMdh0XHQIsowF7YiaJGw8ztOfCf8zAu/sFx2CLGNBO6LtJ5EUARcefIcX2hOBnvjljOgcZAH/Rh3Rxzm4i6c/EwDg4bF4m2+i5YoF7XCq61B8AYOCROcgeUjsg3wdztWKzkHmsKAdzvojuDVOdAiSk7uHY1WO6BBkDgva4XxxGH8YIjoEycmtcVh3GAZJdA5qhQXtWPJ1CPSE1k10DpITN2dM5fl2ssSCdtI1YWAAABX9SURBVCyrf8GdnP2ZWlk0Gu8fEB2CWmFBOxCDhG1FmB4rOgfJT29veLggXyc6B12NBe1ANp/A5Cg4aUTnIFl6YDQ+4JtomWFBO5AVe/FYgugQJFfJUdhfgYv1onNQMyxoR5FZhjBvBHuJzkEytjAeq38RHYKaYUE7in9m4vFE0SFI3v4wBGtzIfF8O9lgQTuEgvOoaURcL9E5SN48XDAuHFuLROeg37CgHcKbHH2mjlk0Gu9miQ5Bv2FBq191HQ5UICVKdA5SgggtmiSUXhSdgwCwoB3BBwdw30hoeHYddcyi0Xhjj+gQBIAFrXp1eqw/wsk3qBOmxiCvCkXVonMQC1r1/n0It8Rxbn7qnOcn4W/pokMQC1rdmiT8Kxv3jRSdg5RmXDjOXMbRs6JzODwWtJr99wQmRMDbVXQOUqBlk/kmWjwWtJq9tQ+P8uw66pIRoahpwJEq0TkcGwtatbJPIcgTYd6ic5Bi/TUZL+wQHcKxsaBV6+97sGS86BCkZIOD4OyEg6dF53BgLGh1KruIs1cwhNd2U/c8O4Ej0SKxoNVp6XY8ybfP1G2DguDhggOnROdwVCxoFfruGJw0SIkWnYNU4cUULNsuOoSjYkGrzaUGvJyO16eLzkFqEalFHx9sPyk6h0NiQavNX7Zh8Tjet5us6bmJHIkWgwWtKlkVKLmAm68RnYPUpbc3+gfgx0LRORwPC1o9Gpvw+Ga8eZ3oHKRGf03G0jRcaRSdw8GwoNXjb+m4fRjCfUTnIDUK8MAT4/DMVtE5HAwLWiUOV2FfOf5nuOgcpF5zrkHZRWSUis7hSFjQamCQ8NBG/HMGZ+Un23rrejz9ExoMonM4DBa0Gry7HzP6op+/6BykdqE9cfswvLJLdA6HwYJWvPJL+OIwHk8UnYMcwz3DkV6CY5wq2i5Y0MomSfjT93h9OnrwSJJdOGnw5nV4eBMkSXQUB8A/a2V7eReSIjA2THQOciTXBCIlGm/vF53DAbCgFWx3KXaXYsk40TnI8TyRiK/yUHJBdA61Y0Er1blaPL4Zq2bCiWdukN25OOGfM/DwJtE51I4FrVSPbMJzE9HLS3QOclTDghHrh88Oic6haixoRVqVg0BP3NBfdA5ybMtTsDIbO4pF51AvFrTy/HwGq3/mhKIknmcPbLgVf9mGbM7obxssaIW50ogHvsfKm+DCQ0cy4OOGL+bi/u95ZrRN8K9cYf68GX9ORIyf6BxEv+ntjU9m4bavcaZGdBTVYUErydpc1DZi3iDROYiuNigI/28aFqznfKRWxoJWjP/bh/8cxQc3is5BZM7kKDwwGn/8Gvom0VFUhAWtAPom3PsddFfw2c1wdxGdhsiCeYOQ2g93fsurwK2GBS13dXrc9jUGBWHZZM4mSnL3PyMQqcVTP6GJHW0NLGhZO1eLGz/HzdfgzwmioxB1zEspCPTEzLXQXREdRflY0PKVr8P1n+GlFMwfLDoKUYdpNHhyPJZNxg3/xk+8z2z3sKBlKr0Ed3yDD2diDGeqIwUaGYoNC/B6BlbsFR1FyVjQspOvw9x1+NcBfHsrBgeJTkPUVcFe+OGPqKrBLetxsV50GmXiOQEyUn4Jy7bjzGUsT8HQYNFpiLrNWYPlKdh0AjPW4N0bMIyv6k5iQctCdR1e2YV95Vg6GZMiRachsqrr+mJwEBZvAYBHEzA+XHQg5eAQh2ClF/G3dFz/GUaHYesdbGdSpwgt1s3Dy1Ox7jCmfIq1ubyepUP4DloM3RWsP4Jv8+HngT8MwZPjOfkRqV+sH1bMwIV6/OsAJn+MmQNx7wj4uouOJWMsaLu61IAN+Vh3GI1NuGUw1s2Dj5voTET2pXXD4nF4LAFf5WHOOrg6Y3osZvTFNYGik8mP9QtakqSHHnooLy/Pw8Nj1apVvXr1Mrs8KCjI7Goqo2/CkSrsr8D+cuTp0MMJ1/XDuzcgzFt0MiKhXJxwy2DcMhhna7GlAK/swlEd4kMwoy+mxMDbVXQ+ebB+QaelpVVVVW3btu2jjz564403XnnlFbPLp0+fbnY1pbvUgCNVOHQGuZU4VIkmCYOCMLo3Hh6LgYFw5rXaRFcL8MCCOCyIgyTh4Gn89wTeP4A6PXq64ppADAzEoCAMDHTQkRDrF/SuXbsSExMBJCQkfPzxx5aWe3p6ml1N/hoMuFiPi/U4X/fbF7XI0yG3Elca4dUDg3shrhfuGIa4XnB1Fh2XSCE0GowIxYjQX7+91IB8HfKqsOEoXt2NC3UA4O+BPj6I9EUfH/TxQYQWXj3U3N3WL2idThcXFwcgMjJSp9NZWm5pNQCSJD366KP19b+f2n7s2LEZM2ZcuNDdm7zXGTTP7nA3tJrG5UqjpsEAAE0SLjZc9S5X34SaRo3WTbpYr2mSIAEeLvB1l3zdJNN/tW7SteGGR+Kbeva4atO1l1HbzcREDqy/F/p7YWbU70vO12kqLjuVXtSUnnXaW+xUfklTq9dcqNcY58/TaCBJV80p5uMqGW977+oMzx4t//KdNXhpUp27c3cndqqqqtqzZ8/u3btNS9zc3FasWKHp9vRm1i9oPz+/4uJiAMXFxf7+/paWW1oNgEajeeKJJwwGg2nJN9984+np2bNnz25m85Lw9ESLcyH2cEZPVwDwdWs9b5ym1bccrSCyt549ER6IsR1YU5JQXQ9Ac7kBjQbz62g0CPDx6v4kkT4+PhMmTJg9e7ZpibOzc/fbGbYo6IkTJ65cuRJAVlZWUlKSpeWWVjOKjLzqfOCgoCC9Xu/sbIXxgr4B3d8GESlAoAsABHrZfEdubm5eXl4xMTFW37L1Czo5OXnDhg2pqakuLi4ffvjh4cOHb7/99uzs7BbL/f39m39r9RhEREqnkZRw84NPP/1Ur9fffffdooMQEbX00Ucfubi43HHHHVbfMi9fIyKSKRY0EZFMsaCJiGSKBU1EJFMsaCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQREQyxYImIpIpFjQRkUyxoImIZMr6N421kR07duj1etEpzCssLLTFDX3l4NSpU1qt1tPTU3QQm1DxgauqqnJzc/Px8REdxCbkduB27949ZcoUW2xZGQU9depUDw8P0SnMMxgMmzZt+stf/iI6iE1s3bp14MCBYWFhooPYxPfff7906VLRKWxiz549gYGBkZGRooPYhNwO3A033JCUlGSLLSvjrt5y1tjYmJqaumXLFtFBbOKFF15ISkpKSUkRHcQmkpOT09LSRKewiRUrVkRERMyePVt0EJtQ8YFrgWPQREQyxYImIpIpFjQRkUyxoLvLyclJrSc5AHB1dXV1dRWdwla8vLxER7AVHjh14IeEREQyxXfQREQyxYImIpIpFjQRkUyxoImIZIoF3V3Lly9fv349AEmSFi1alJKSkpqaWllZKTqXddTX1/v6+sbHx8fHx7/22mui41iNKg+WkVoPGdT+t2YWC7rrDAbDpEmTli1bZvw2LS2tqqpq27Ztc+bMeeONN4RGs5rCwsKZM2fm5OTk5OQsWbJEdByrUeXBMlLlIXOEvzWzWNBd5+TktHXr1ieffNL47a5duxITEwEkJCRkZGQIjWY1x48fz8/PnzVr1vz580tLS0XHsRpVHiwjVR4yR/hbM4sF3XUajcbFxcXJ6dffoU6nM04eFhkZqdPphEazmqCgoMWLF3/77bc333zzQw89JDqO1ajyYBmp8pA5wt+aWSzoTlu1atXcuXOXL1/eYrmfn19xcTGA4uJif39/EdGsxvQcExMT582bB2DmzJmHDh0Snctq1HSwWlDrIWtOxYevBWXMBy0rd91111133dV6+cSJE1euXAkgKyvLRpPD2o3pOb766qve3t4PPvhgZmZmXFyc6FxWo6aD1YJaD1lzKj58LbCgrSY5OXnDhg2pqakuLi4ffvih6DjWce+9995zzz1r1651d3d/9913RcexGlUeLCO1HrLmVHz4WuBcHEREMsUxaCIimWJBExHJFAuaiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRyRQLmohIpljQwvzyyy/Tp0/XarUBAQE33XTT8ePHrbjxnJycDs5kVl1d7evr22Jhnz59NL/x9vZOTU2tqKhoYyNLlizx8/OTyc2HTM8oKytr1KhRllZzcXHR6/XNl5jWX7NmzW233dbB3Rm30/a+uqbd36otdmr29UCisKDFMBgMqampY8aMOXjwYG5u7oABA2bPni2rias2btx4/vz5c+fOZWdnX7x48dlnn21j5ZUrV+bl5fXq1ctu8ToiOjr6hRdesN361nqsJe3+Vm2xU5IVFrQY5eXlZWVlTz31VExMTGho6CuvvBIdHX3hwgUA//rXv6Kjoz08PBISEo4ePQogPz9//PjxixcvDgwMTEpK2rNnz+jRo729vR977DEAa9asuffee++44w5fX9/x48cbH9Jcenr68OHDvby8ZsyYUV5eblz41ltvhYeHh4eHf/TRR2YTent7+/r6+vn59evX77bbbissLLS0tdmzZ1+4cGHMmDFVVVWtf5qbmzt58uTly5cPHTrU7MPz8/OTkpJef/31sLCw6Ojobdu2GXe0bt26fv36BQQEPPDAA/X19ZaeiEnrZ1RUVPT8888D0Ov1DzzwgJ+fX2Bg4Isvvghg+vTpBoMhNjZ27969pnim9QHU1tbeeuutWq02ISEhNzcXQGZmZkJCgvGnpq9N28nNzTU99uuvvx4wYIBWq50zZ05VVVUbz9Gk9UOa/1aN67R+Fs0Df/zxx1FRUVFRUZ988klUVFQbO+3s66H1fs3+Ktp+lbbeCHWIRCI0NjbGxcVde+21mzZtunLliml5SUmJq6vrjh07qqqqFi5ceN9990mSlJeX5+Tk9Nlnn509e3bkyJG9evU6efLknj17AFRWVq5evdrFxeWdd96prKx88sknhw0b1tTUdPDgwcGDB0uSpNPpAgIC/vOf/5w7d27RokWTJ0+WJGnnzp1+fn47duwoKytLSUnRarUt4oWFhaWnpxu/Li8vnzVr1ksvvWRpa5IkabXaS5cumf3poUOHtFrtwoULc3Nzza6Ql5fn5eX18ssv19TUPPnkk4mJiZIkHT16NCAgICMjo6CgYOTIkStXrrS0ayOzz2j//v0jR46UJGndunUDBgwoKirKzs52c3M7ceKEJEnOzs6NjY3N45nWX716NYD333/f+Cu95ppr9Hr9nj17xo4da9xd86+N2zE9trCwUKvVbtmy5ezZswsXLrzlllssPUcTsw8x/VZNq7V+Fqad/vLLL4GBgXv37i0vL09KSoqMjLS00y68Hlrv1+yvou1XqdlDQO1iQQtTV1f3zjvvzJgxIyAg4Nprr923b58kSbW1tcXFxZIkXb58efHixaY/77CwMOOjnnrqqfvvv9/4dWRk5PHjx1evXj106FDjkoaGBn9//2PHjpkK+uOPP54zZ47xp7W1tZ6ennq9/tFHH3366aeNC3fv3m22oL28vLRarY+PD4CEhAS9Xm9pa9JvVWL2p4cOHXJ1da2rq7P08Ly8PB8fn8bGRkmSDh06NGDAAEmSXnzxxYcffti4Zk5Ozo4dOyzt2sjsM2pe0DExMXv37m1qaqqqqqqvr5eaFbQpXvOCHjFiRPNfaV5eXgcL+h//+Medd95p/FFlZWWPHj0sPUcTsw+RzBV0i2dh2ukzzzyzZMkS42obNmwwFXTrnXbh9dB6v5YKuo1XqdlDQO3iHVXEaGhokCTpgQceMP77/fPPP58wYYLx354rV67ctGmTVqt1c3Pz9vY2rt+zZ0/jFy4uLiEhIaavjV9ER0cbv+jRo0dUVFR5ebnpc57S0tItW7YY/80LwNXVtbKy8vTp01OnTjUuiYmJMZtw1apVo0ePBqDT6RYsWLBmzZo777zT7NZCQ0Pb2BeA8PBwNze3NlYICQkxPhfTMyorK+vXr5/x62HDhgHYuXNnG7tu+xndfPPNFy9evO+++86cObNo0aInnnii+U9N8Zpr8Ss9ffq0u7u76aeS5U8LTp8+bQoZFBTk6upqHKNo/RzbfojpKHfkWZSXlxsPFoCIiAjT8tY77cLroe3fXvNfRRuv0rY3QpZwDFqML7/8MjU11fi1m5vbwoULExMTDx48+OWXX/7www+bN2/+6aefFixY0MGtFRUVGb/Q6/UlJSWm2gIQGho6bdq0kydPnjx5sqCg4McffwwJCendu3dBQUGLx7YQGhpqHNMcNWrUnDlzDh48aGlrbe8LzdrB0goajabF3oODg8vKyoxf79mzZ/Xq1W3vuu1nVFhYmJKSkpOTs3fv3u+++67FTZJaN2bzjTQ2NhYXF4eHhxt/vcaFpmythYSEGO9nCsD4VjEwMNDsc2z3IR1/FqGhoSUlJcavS0tLTctb77QLrwez++3Ir6KD4akNLGgxpk2blp2dvXTp0uPHjx89evTdd989cOBAcnLy2bNne/bs6eHhUVlZ+dZbb9XW1nZka7/88sv777+v0+mee+653r17m957AkhNTU1PT9+4caNOp3v66acfe+wxjUYzd+7c999/Pz09vaKi4vnnn2+jO4xCQkKMf/Zmt9b2vppvp90VTObMmbN69eq9e/cWFhY+9thjOp2u7ce2/Yz+85//LFiw4MyZMwaDob6+3sPDw7j88uXLlp5yTk6OceD72WefjY6OjomJ0Wq1P//8c05OztmzZ99+++3mKzffzo033vj1119v3br1/Pnzixcvnj17ttn/ATTXwYdYehbGp79q1aqsrKxTp079/e9/b2NfXXg9tN5vG78KS9oIT20RPMTiwI4dO5aamhocHOzt7Z2YmLhp0yZJkqqrq6dNm+bv7z9u3LjvvvsuODj4008/zcvLM41aPvvss0uXLjV+HRsbaxyDvv766+fOnevt7Z2QkHDkyBFJkkxj0JIkbd68efDgwZ6ensnJyQUFBcaFb731Vp8+fcLCwj766CPT0KFJ8w8JJUn64YcfevXqdeHCBUtbM42Wtv5piyHX1is0f3bNv/7kk0+io6N9fHzuuusu45Cl2V2btH5GpiHaS5cuzZo1y8vLy9/f/09/+pNxfGn+/Pne3t6ZmZmmPTYfg77vvvtmzZrl7e09YcIE4ydaTU1NDz/8cM+ePYcMGfLll1+aBmGN29m+fbvxsZIkrV+/vn///t7e3rNmzTpz5kwbz9Gk9UOkVmPQrZ+FKbAkSe+9915oaGj//v3fe+8946G3tNPOvh5a79fsr6LtV6nZQ0Dt4k1jFW/NmjXff//92rVrRQchYfLz88+cOTNp0iQAW7Zsefnll9PS0kSHIivgEAeR4p0/f37BggWVlZW1tbVvvfXW9ddfLzoRWQcLmkjxEhMTH3nkkeHDh/fr1y80NHTRokWiE5F1cIiDiEim+A6aiEimWNBERDLFgiYikikWNBGRTLGgiYhkigVNRCRTLGgiIpliQRMRydT/B3K3OpaxgYOjAAAAAElFTkSuQmCC"
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "and _p-value_"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(refsamp >= sum(Secchi$diff))/length(refsamp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "[1] 0"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that a sample of size 100,000 may not be large enough to evaluate this very small _p-value_ accurately.\n",
      "\n",
      "# Randomization tests for independent samples\n",
      "\n",
      "If we have collected responses under two different conditions and the allocation of conditions has been randomized then a test of, say, $H_0:\\mu_1=\\mu_2$ versus $H_a:\\mu_1<\\mu_2$ can be based upon the set of possible combinations of response values.  Under the null hypothesis, the sample mean of the particular combination of responses that we saw for the first condition should be similar to the other possible sample means of combinations of responses of this size.  Under the alternative hypothesis the sample mean from the combination of responses we saw should be systematically smaller than the other possible sample means.\n",
      "\n",
      "In the first chapter of Bob Wardrop's [course notes for Statistics 371](http://www.stat.wisc.edu/~wardrop/courses/371chapter1fall2013a.pdf) he describes an experiment performed by a student on the consumption of treats by her cat according to whether the treats are chicken or tuna flavored.  The experiment lasted for 20 days. The assignment of chicken or tuna on each day was randomized, subject to the constraint that each flavor is provided exactly 10 times.\n",
      "\n",
      "Each day 10 treats of the indicated flavor were provided and the number consumed was recorded"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "treats <- c(4,3,5,0,5,4,5,6,1,7,6,3,7,1,3,6,3,5,1,2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The tuna treats were provided on days"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tuna <- c(2,3,4,6,10,12,14,17,19,20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus the number of tuna treats consumed (out of a possible 100) was"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(treats[tuna])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "[1] 29"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The number of chicken treats consumed, also out of a possible 100, was"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(treats) - sum(treats[tuna])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "[1] 48"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The test\n",
      "We want to compare the observed sum, 29, to the possible sums of any 10 of the 20 days.  It is not easy to enumerate the possible sums.  I don't know of a good way in R of stepping through all the"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "choose(20L,10L)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "[1] 184756"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "possible selections.\n",
      "\n",
      "The alternative is to produce a random sample from the reference distribution.  About the best way I can think of doing this is"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "system.time(tunaref <- replicate(100000,sum(sample(treats,10L))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "   user  system elapsed \n",
        "  0.697   0.001   0.698 "
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "with the histogram"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "histogram(~tunaref,breaks=seq(min(tunaref)-0.5, max(tunaref)+0.5,by=1.0),\n",
      "          xlab=\"Sample from the reference distribution of sums\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": []
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHgCAIAAADytinCAAAgAElEQVR4nO3de1wU973/8e8icscF8ZqIKMajUR+Kl1hJ0AhBG0OSRjGm5uLBaIxWezRt1DysxjQxF6tHf02MNvWSPMQkPvSYJo2XRhFjRcULiIqCIaAgoHIRlZvALvP7Y3v2UJxZEJfZ77qv51+785n57me/C2+G2d0Zg6IoAgAgHzdHNwAAUEdAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkJQTB7TZbFYUxdFdOBmTyeToFpwMM3a3mDE7cuKAXr9+/a5duxzdhTO5cOHCwoULHd2Fk3nuuecc3YKTiYuLu3HjhqO7uE84cUDX1tbW1tY6ugtnUltbW1dX5+gunExlZaWjW3AyNTU17ETbixMHNADc3whoAJAUAQ0AknJ37MMXFBS0+I2+I0eOBAYGlpSU2Lel+1hBQUF6evpf//pXRzfiTAoLC5mxu5KTk7N582Y/Pz9HN6KrmJiYBx980O7DGhz7SbXNmzcnJydHRka2YNvi4mIPDw+j0Wj3ru5XNTU1RUVFwcHBjm7EmWRnZ/fq1cvRXTiTixcvhoSEuLm50H/nBw4cGDFixJQpU+w+soP3oIUQw4YNe/755x3dBQC0UHl5eSuN7EJ/5QDAuRDQACApAhoAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQlOOv6g3YV1lZWU5Ojla1vr7e3d3daDRqrRAYGBgYGNg6rQF3h4DG/Wbt2rWL09NFSIh6+cQJceuWGDJEvWoyxVy7tnPnztZrD2g+Ahr3m/r6ejF9unjiCfXyypWiulosWaJera6uj41tvd6Au8IxaACQFAENAJIioAFAUgQ0AEjK/gFdV1f30ksvjRgxIiIi4uLFi3YfHwBchP0DeufOnW3btk1OTp45c+by5cvtPj4AuAj7B3S7du1u3bplNptv3LjRrl07u48PAC7C/p+DHj169OLFi//jP/6juLg4IyOjYUlRlD//+c81NTXWJampqYMHD7Z7DwCgm6KiolOnTl25csW6xNPTc+7cuQaD4R5Htn9Ar1y5cuzYsUuXLj169GhcXNy+ffusJYPBMGjQILPZbF1y/fp1b29vu/cAALrx9vbu0aPH0KFDrUvatGlz7+ksWiOgS0tLg4OD3dzcgoKCiouLG1UjIyMb3i0sLDSZTHbvAfe3lJQUG2fbOHfunHj0UT37gYvz9/fv379/dHS03Ue2f0DPnz9/ypQpX375pclk+vTTT+0+PvD222/vHjdOeHiol0+f1rcdoLXYP6A7duy4Z88euw8L/JtXXxU+Puqlbdv0bQVoLXxRBQAkRUADgKQIaACQFAENAJLihP1AA3V1hw8fHjNmjFbdZDIdOHBAz47gyghooIG6ulsREQm7dmmu8O8f5AdaFYc4AEBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCS4oT9kNEbb7xRUFCgVT1z5oyezQCOQkBDRv/v2DGxZYtmefhwHXsBHIaAhpQ8PUVoqGbVjUNzcAn8oAOApAhoAJAUAQ0AkiKgAUBSBDQASIpPcQB34+efx4wZo1V0c3PbvHlz586d9ewI9zECGrgbFRUJ+/ZpVmfPLi0tJaBhLxziAABJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQFAENAJJyd3QDcEUVFRXLly+vq6tzdCP2VlCwatWqDh06aNWnTZvWu3dvPTuCUyOg4QB5eXnLTp8W//VfmmscO6ZjO/aTnb3xmWdESIh6df/+hw4eJKDRfAQ0HCQ4WERHa1bff1/HVuwqPFz066deunpV3L6tbzdwbhyDBgBJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJNUqAf3f//3fkZGRv/jFL3JyclpjfABwBfYP6FOnTn3//ff79+9ftGjRsmXL7D4+ALgI+5/Nbvfu3ZMmTXJzc3vmmWciIiLsPj4AuAj7B/S1a9eys7PHjBljMBj+9Kc/BQUFWUuKonzwwQfl5eXWJenp6YQ4AKdWVFSUlJR0/vx56xJ/f/9FixYZDIZ7HNn+Ae3v73/79u09e/acPHly+vTpJ0+etJYMBsOoUaNqamoarW/3HgBAN35+fv37949ucH5zX1/fe09n0RoB/eijj+7fv9/d3b19+/b19fWNqiNHjmx4t7Cw0GQy2b0HANCNj49Po4C2F/sH9Lhx43744Yfw8HCTybRmzRq7jw8ALsL+Ae3m5vbxxx/bfVgAcDV8UQUAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkZf8rqgBQV1ubm5ubkpKiVe/WrVvnzp317AiSI6ABvZw8uezSpWV1derVGzemm83r16/XtydIjYAG9KIoIjZWTJumXs3OVj78UN+GIDuOQQOApAhoAJAUAQ0AkuIYNFpLWVmZVunWrVt6dgI4KQIareLrr79+8f33Rdeu6uWSEvHoo/p2BDgfAhqtoqqqSixcKF55Rb28c6fYs0ffjgDnwzFoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKRUvqgyYMAArbXT09NbsxkAwP9RCegtW7bo3wcAoBGVgA4LC1NdNT4+XqsEALA7zXNxZGZmrl69+ubNm5a71dXVx44de0Xr1AoAAHvTfJNwypQptbW1wcHB5eXlTz/99NWrVzdu3KhnZwDg4jT3oM+cObNnzx4vL6/nnnvu5ZdfHjNmzMSJE2NiYvRsDgBcmeYedJcuXc6fP+/r61tWVlZWVtauXbvz58/r2RkAuDjNPehFixZFR0dnZWWNGzcuOjraz89v2LBhenYGAC5OM6BfeumlcePGdezYcenSpQ8//PD169djY2P17AwAXJxKQJtMJiFEv379srOzLUsmTZpUUVHRvXt3LiUHALpRCWgvLy8hhNlsttywmjhxok5NAQBs7EGPHTt27969uvcDAPgXzU9xWNLZbDZfvXrVEtkAAD1pBnRxcfELL7zg5eX10EMPeXt7v/DCC8XFxXp2BgAuTjOgZ8yY4ePjU1hYWFFRUVhY6OXlNXPmTD07AwAXp/kxu/379+fm5gYGBgohOnbsuGrVqtDQUB0bAwBXZ+ubhCkpKda7p06d6tq1qy4tAQCEsLEH/eGHH06cOPG5557r0aPHpUuXvv32202bNunZGQC4OJU96Pz8/Pr6+tjY2NTU1LCwsNu3b4eFhaWmpk6YMEH//gDAZansQQcHB5eVlQUEBISGhs6bN0//ngAAgovGAoC01I9BJyUl+fn53bl89OjRrdsOAOB/qQf0b37zGzc3lZ3rS5cutW47AID/pR7QZ86cCQgI0LkVAEBDHIMGAEmpBPTjjz/u7q75+WgAgD5UgvjHH3/UvQ0AQGMc4gAASansQRcVFXXq1Km4uLhjx476NwRn8e677yYmJrZt21a1mpOTI955R9+OgPuNSkD36dMnIyNj+PDh58+fb1RS/XA0XFNeXt7BjRtFr17q5bg4Xbu5D5jNV65cSUhI0Kp36NAhLCxMz47gcCoBPX369IcffvjmzZvdunVrVLpx44YuXQGup7Bwd3Hx7gankGzkkR07jh8/rmdHcDiVgF6xYsWKFSvGjx//t7/9Tf+GABelKGLgQLFwoVbd9x//0LMdyEDzTUJLOnNNQgBwFK5JCACS4pqEACAprkkIAJLimoQAICmuSQgAktLcg+aahADgWLbOWsc1CQHAgThZEgBISjOgKysrm1wCAGg9KgFtMplMJlO/fv1MDdy4ceOuPsVx8+bN7t27269PAHA5Ksegvby8hBBms9lyw2rixInNH3fJkiWlpaX32BwAuDKVgLaceWPs2LF79+5t2aAnTpwoLy8PDg6+p9YAwLVpfoqjxelsMpkWLVr05Zdfjho1qlFJUZR33nmnpqbGuiQ9PT08PLxlDwQAMigsLExOTm54An1PT8933nnHYDDc48i2vuq9ZMmS69evN1yYmZnZ5Ihr1qyZNGlSp06d7iwZDIbx48ebzWbrEh8fn4CAgLtpGADkEhgYOHz48JiYGOuSNm3a3Hs6CxsB/eqrr06ePPnll1++2yt8p6amXr16dceOHZcvX37qqad2797dsNrokhDnzp3jXKYAnJq3t3ePHj2GDh1q95E1w7eurm7p0qXe3t53O+LmzZstN/r27dsonQEAzaf5Oejf/e53H3/8ccPDEXerOcdDAABaNPegv/3227S0tA8++KBLly7WgylkLgDoRjOgN2zYoGcfAIBGNAO6b9++Qgiz2VxUVNRwJxoAoA/NY9AFBQWRkZHt2rXr169fSkrKyJEjL168qGdnAODiNAN66tSpAwYMKC0tNRqNYWFhI0aMeO211/TsDABcnOYhjqSkpG3btllOx+Hu7r5w4cKQkBAdGwMAV6e5B927d++kpCTr3WPHjnHRWADQk+Ye9McffxwbGzt69Ojr16/HxsYeOnRoy5YtenYGAC5OM6Aff/zxCxcufP/992FhYV27dv3000+7dOmiZ2cA4OI0D3HU1NR8++23oaGhixcv9vT03Lp1a21trZ6dAYCL0wzouXPnrlu3zmg0CiF69er11VdfzZo1S8fGAMDVaQb09u3bt23bNmjQICFEeHj41q1bd+zYoWNjAODqNAM6MDCwuLjYevfatWsdOnTQpSUAgBA23iR8//33Y2JiXnzxxZCQkPz8/Pj4+NWrV+vZGQC4OM096BdeeOHw4cOdOnXKysoyGo2JiYmvvPKKnp0BgIvT3IMeOHDg119/vXjxYj27AQBYae5BT5o0aeXKlQ0v8AoA0JPmHnRCQkJaWtpXX30VHBxsvSwhJ+x3KVu3bs3NzdWqnj17Vs9mABekGdBr1qzRsw9I6JNPPjnyxz9qlvkJAVqZZkAPGDBAcMJ+1+bh4SGiozXLPj469gK4Ik7YDwCS0tyDtpywf8+ePX379rWesD8hIUHP5gBYXbhwYcyYMVpVHx+fbdu2eXp66tkSWhsn7Aecw5U+fa7s26dZnjChqqqKgL7PcMJ+AJAUJ+wHAElxwn4AkJR6QJ89ezYtLa1Xr15xcXH69gMA+BeVY9Dr1q0bOnToypUrn3jiiYULF+rfEwBAqAb0Rx99tGXLltOnT//4448rVqyorKzUvy0AgEpA5+XlRUVFCSGGDx/u7e3d8LT9AADdqH/MznJ2JIPB0LZtW337AQD8i/qbhKdOnfL39xdCmEymM2fOlJSUWJYPGzZMv9YAwLWpBHTnzp0nT55sue3n5zdjxgxr6erVqzr1BQAuTyWgSWEAkIHmV70BAI5FQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEm5O7oBOExtbe3UqVPr6uq0VsjKytKzHwCNENCuq6qq6qsbN8SWLZpr9OypYzsAGiOgXZuHhwgM1KwaDDq2gntz4UJUVFSbNm1UiwaD4fPPPx8wYIDOTeEeEdDAfaGsLO38eREQoF7905/y8vIIaKfDm4QAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIyv4BXVNTM3Xq1KioqCFDhhw/ftzu4wOAi7B/QO/du9fPzy8xMXH9+vVz5861+/gA4CLsH9DdunWbM2eOECIoKMjANZMAoKXsf8mrwYMHCyFOnDgxc+bMZcuWNSwpijJ//vzy8nLrkp9++ik6OtruPQCAbq5cuZKQkHD48GHrEn9//xUrVtz7Hqr9A1pRlMWLFx86dGjTpk2DBg1qWDIYDHFxcTU1NdYlu3fvDgoKsnsPAKCb9u3bR0ZGxsTEWJf4+vra5fiB/QN6+/bt2dnZiYmJ7u4qgze6bOW5c+dMJpPdewAA3Xh6evbo0WPo0KF2H9n+Ab13797k5ORhw4YJIYKDg7///nu7PwQAuAL7B/SGDRvsPiYAuCC+qAIAkiKgAUBSBDQASMr+x6ABSKek5B//+Ed+fr5W/bHHHuvfv7+eHaE5CGjABZw798mgQSIwUL2ambl469b33ntP357QNAL6flZVVXXkyBGtakVFhZ7NwMEiIsRTT6mXkpLEDz/o2w2ahYC+nx08ePCpTz4Rjz+uXi4r07cdAHeHgL6fKYoiRo8WCxaol3Nzxbx5+nYE4C7wKQ4AkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqABgBJEdAAICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASXFVb+e2evXq323bJvz81Mt5eWLaNH07AmA3BLRzu379ulixQkREqJf/+Ed92wFgTxziAABJsQcNuLzKytTU1L/+9a9a9YceeigqKkrPjmBBQAMuLzt7t5fX7sBA9aqiRL3/PgHtEAQ0ACH69RPPP69eqq8Xn32mbzf4F45BA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIqz2cnu8uXLRUVFWtUrV67o2QwAPRHQsnv11VcThgwRBoN6+YcfRFycrg3B9eTn5y9fvlyrGhgYOGPGDD37cR0EtOwURREffCDatFEvHzmibztwPWbzT76+bw0dqrnCggUEdCshoAE0pX17ER2tWTUadWzFtfAmIQBIioAGAEkR0AAgKQIaACRFQAOApAhoAJAUAQ0AkiKgAUBSBDQASIpvEjpeTExMbW2tVvXUqVN6NgPctfr6srIyraLBYAgICNCznfsJAe14u6uqxIEDmuWgIB17Ae7emTPtX39ds3r6dPo33/Tv31/Hhu4fBDSAe7Ztm2bpd7+z8Q8ibOMYNABIioAGAEkR0AAgKQIaACRFQAOApPgUB4DWdObM66+/btS46kptbe38+fOffvppnZtyFgR0q8vMzHz4ySdF796ObgRwhLKyExs2iMGD1avbt0/KzdW3IWdCQLe66upqMWGCWLVKc43ISB3bAeA0OAYNAJIioAFAUhzisIOampqCggKtqo0S4OpKSvbt21dRUaFVHzZs2BNPPKFnR1IhoO3giy++mPk//yNCQ9XL2dli4EB9OwKcRFbWdyEh3w0dql4tKXn5iy8IaHtSFGXOnDkZGRne3t6ff/55p06d7P4QsjGbzeK118SkSerlrVvF8eP6dgQ4jz59RHS0eik/X+zapW83crF/QB84cKC4uDgxMXHTpk2rVq366KOP7P4Q+luyZElSUpK7u/p0/fzzz2L5cp1bAu5/ilJWVpaSkqJVNxqNDz30kJ4d6cz+AZ2UlBQeHi6EGDFixBdffGH38VvJp59+evnyZa3qrl270vfsEd26qZd//evWagtwZUVFuy5d2rV9u1a925dfvvTSS1pVd3f3t99+28PDo3Wa04P9A7qkpGTAgAFCiJCQkJKSkoYlRVHmzp1bU1NjXfLTTz+NHj26xY+1cOHCGzduWG7fvHnzxIkTWmvW19cXFRVp7QILIaqrq+tsHCnOyhJvvSV8fdWrKSmitFTs369ezc4WRUXCxhnNL12yVa2sFLNmCYNBvfrzz2LFChEfr15NTRVt2ojsbPVqRYVIS7P10Ldv26peuyYWLxbt2qlXk5NFXp5ISlKv5uWJvDxbg+fn26pWVIjf/lZovZqZmWL1as0zFJ85I0wmkZ+vXq2pEenpLZ+TwkLxzjsiMFC9mpQkMjI0j3cVFoqsLFuDX71qq3rzpnjjDaEVRunpYs0a8d136tXz50VlpSgqUq/W14vMzJbPSV6eWLZMdOigXk1OFidOiNOn1avFxaK2VmhfriXf23t5QoLmQ58//5e//MWg9bsjREREhPUYbEBAwPKW/h98+fLlH3/88fDhw9Ylnp6ef/7zn208dDMZFEW5xyEaWbp0aUBAwBtvvHH+/PkZM2Yk/fuvaG5urtlstt7929/+5uvrO3PmzJY91uXLl+vq6iy3a2trr127ZmPlW7dutdNKEyHKy8v9/f1tbOvv76813bdv33Zzc9P6Q11fX19VVeXn59cajVVWVnp5ebVp00a1WldXZzKZvL29WzZ4k1U/Pz+tObH8Gfb09FSt1tfXV1ZW2p7wFs9JVVWVh4eH1h9jk8lUW1vr4+PTssFtN1ZRUeHj4+Pmpv7p1draWkVRtOZEUZSKiopWmpPq6mp3d/e2bduqVs1m8+3bt3219j/ubU6a/BE1m81eXl6qVUVRysvLbQxu+6FtV4UQnTt3tv7atm3bNjg42MbKNvzlL3+prKwcP368dUmbNm1CQkJaNlpD9t+DHjVq1IYNG4QQJ0+ejIiIaFRt1HTHjh1NJlOLH6vRhPbt27fFQwFAy3h4ePj4+IRqfY7rHtg/oCMjI7/77ruYmBh3d/eNGzfafXwAcBH2D2g3N7ePP/7Y7sMCgKvhq94AICkCGgAkRUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCSIqABQFIENABIioAGAEkR0AAgKQIaACRFQAOApOx/TcK7dfDgwZZd2LukpMTDw8P2ZdXRUG1t7bVr11p8bXnXlJOT0xpXa76PXbp0qXv37m5uLrTzd/jw4SeeeKI1RnZwQEdHR3t7e7ds2+Tk5Pbt24eEhNi3pftYXl5eSkrKwIEDHd2IM9m5c+fSpUsd3YUz+fzzz1977TV/f39HN6Kfp59+OiIiojVGdnBAP/DAA88//3zLtr1y5Uq3bt0mTJhg35buY2fPns3NzW3xhLumtWvXMmN3Zdu2bc8++2yHDh0c3cj9wIX+DQEA50JAA4CkCGgAkJQTB7SHh4eHh4eju3AmHh4ebdu2dXQXTsbX19fRLTgZT09Pd3fHfzzs/mBQFMXRPQAAVDjxHjQA3N8IaACQFAENAJIioAFAUk4W0DU1NVOnTo2KihoyZMjx48cVRZk9e3ZUVFRMTExRUZGju5NRRUXFs88++/jjjz/22GMXL15kxprp5s2b3bt3F0IwY81RU1MTEBAQFhYWFha2YsUKJs1enCyg9+7d6+fnl5iYuH79+rlz5x44cKC4uDgxMTE2NnbVqlWO7k5G8fHxjzzyyMGDB6dPn7569WpmrJmWLFlSWloqhGDGmiMnJ+dXv/pVWlpaWlra/PnzmTR7cbKA7tat25w5c4QQQUFBBoMhKSkpPDxcCDFixIgjR444ujsZjRo1atq0aUIINzc3o9HIjDXHiRMnysvLLaf9Y8aaIysrKzMz87nnnps0adLly5eZNHtxsoAePHhwnz59Tpw4ERsbu2TJkpKSEsvZ7EJCQkpKShzdnYz69+//wAMPxMbGvvHGG1OnTmXGmmQymRYtWrR8+XLLXWasOTp27Pjmm29+++23EyZMmDNnDpNmL04W0Iqi/OEPf/j973+/adOmcePGBQYG5ubmCiFyc3Pbt2/v6O5kdOvWLZPJtGPHjh07dsyaNYsZa9KaNWsmTZrUqVMny11mrDnCw8Mt5/z71a9+dfbsWSbNXpwsoLdv356dnZ2YmDho0CAhxKhRo44fPy6EOHnyZCudj9XZLVu2LD4+Xgjh5eVVW1vLjDUpNTV1+/btTz755OXLl5966ilmrDmWL1++du1aIURycvKAAQOYNHtxsq96T58+PSEhISAgQAgRHBz83XffzZs3Lzs7293dfePGjZyC9k5Xrlx55ZVXqqurTSbTunXrwsLCmLFm6tu3b2ZmZn19PTPWpOvXr0+bNq20tNTLy2vdunU9e/Zk0uzCyQIaAFyHkx3iAADXQUADgKQIaACQFAENAJIioAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBraszZ86MHTvWaDQGBQU9++yzWVlZdhw8LS1twIABzVlz/vz5gYGBdrwWkbu7u8lkyszM7Nu3r73GtLJ7t/fixo0blnN1nTx5ctiwYVqrWSak4RLr+lu2bHn55Zeb+XCWcWw/VstINavQQkDrx2w2x8TEDB8+/NSpU+np6X369Bk/frxDTla1YcOGjIwM6ymPJSdntz179nz33Xdbb317batFzllFIwS0fgoKCvLz8xcuXBgaGtq1a9ePPvqoZ8+eN2/eFEKsX7++Z8+e3t7eI0aMuHDhghAiMzPzsccee/PNNzt06BAREXH06NFHHnnE399/3rx5QogtW7a89tprU6ZMCQgIeOyxxyybNHTo0KHBgwf7+vo++eSTBQUFDUvjx4+/efPm8OHDDxw4MHr06GXLlg0cOFAI8c033/Tp08doNMbGxhYXFzfZg9XYsWPNZnOvXr0qKysVRVm2bFnnzp179OiRmJjYZDPp6ekNe7hzTWu3xcXFd1ab3DwzMzMiImLlypUPPvhgz3p9m8QAAAcrSURBVJ49rS1t27atd+/eQUFBs2bNqqmpsd2kEOKTTz4JDg4ODg7etGmTZcnFixfffvttIYTJZLJcCaFDhw7vvfdewwk5duyYtT3r+kKI6urqX//610ajccSIEenp6UKI5OTkESNGWKrW29Zx0tPTrduqvkyqz9Hqzk0azqplnTufhWpLtn8k7hwE90qBXurq6gYMGPDLX/5yz549VVVV1uV5eXkeHh4HDx4sLi6Oi4ubMWOGoigZGRlubm5ffvllaWnp0KFDO3XqdOnSpaNHjwohioqK4uPj3d3d165dW1RUtGDBgkGDBtXX1586dap///6KopSUlAQFBf3973+/fv367NmzR48e3agTo9FYXl5+9uxZo9EYFxeXnp6ek5NjNBr37t1bWloaFxf3wgsvNNlDwwHbtGlTV1eXkZFhMBg++OCDqqqq9957LyIioslmGvagtaalW9Vqk5tnZGT4+vp++OGHlZWVCxYsCA8PVxTlwoULQUFBR44cyc7OHjp06IYNG2w3+c9//jMwMPDgwYP5+flRUVFGo1FRlBMnTgwdOlRRlG3btvXp0+fixYupqamenp4///yzdUIatmdd33L9hM8++8zy2j388MMmk+no0aO/+MUvLA/X8LZlHOu2Wi/Tnc/RSnUT66xaV7vzWai2ZPtHQnUqcC8IaF3dvn177dq1Tz75ZFBQ0C9/+cvjx48rilJdXZ2bm6soSkVFxZtvvmn9rXvwwQctWy1cuPD111+33A4JCcnKyoqPjx84cKBlSW1tbfv27X/66SdrQH/xxRexsbGWanV1tY+Pj8lkatiGNaA9PDxu376tKMrq1av/8z//01ItKipq27atyWSy3UPDAa0B3a5du7q6OkVRzp0715xmGvagtaalW9Vqk5s3bOns2bN9+vRRFOW999777W9/a1kzLS3t4MGDtpucO3fuW2+9Zbl9+PDhOwM6NDT02LFj9fX1xcXFNTU1SoOAtrbXMKCHDBnS8LXLyMhoZkBrvUx3Pkcr1U0UtYBu9Cy0AtrGj4TqVOBeuDt2/92l1NbWKooya9Ysy7/VX3/99ciRIy3/WW/YsGHPnj1Go9HT09Pf39+yvp+fn+WGu7t7ly5drLctN3r27Gm50bZt2x49ehQUFFjevBJCXL58ee/evT169LDc9fDwKCoq6tq1650tBQcHe3p6CiGuXr1qXb9jx44eHh6Wf35t93CnLl26WKpubm7NbMbag+01VatNbt6wJWvb+fn5vXv3tty2XDvtn//8p42Hvnr1anR0tOV2aGhoo6c8YcKEW7duzZgx49q1a7Nnz/7973+vOsMNNXrtrl696uXlZa0q2m9LaL1Mdz5H25tYX8pmPouGLdn4kbA9CFqAY9D62b59e0xMjOW2p6dnXFxceHj4qVOntm/fvmvXrh9++CEhIWHy5MnNHO3ixYuWGyaTKS8vr2Hkde3adcyYMZcuXbp06VJ2dva+ffvu/IW0sP4+d+nSxXKVTyGEZd+nZZcpMhgMjZY02Yy1B9tralWb3PzOljp37pyfn2+5ffTo0fj4eNsP/cADD2RnZ1tuW6fdKicnJyoqKi0t7dixY99///3GjRtVn11D1kHq6upyc3ODg4OFENZPfVh7u5PWy3Tnc2xyk+Y8i+a01OQguBcEtH7GjBmTmpq6dOnSrKysCxcurFu3LiUlJTIysrS01M/Pz9vbu6io6JNPPqmurm7OaGfOnPnss89KSkqWLFnywAMPWHcJhRAxMTGHDh3avXt3SUnJW2+9NW/ePBu/wBbPPPPMN998s3///rKysjfffHP8+PE2dpNVVVRUqC5vfjO212xynOY/UGxsbHx8/LFjx3JycubNm1dSUmJ724kTJ3722WeHDh0qLCx8++23Gw3797//ffLkydeuXTObzTU1Nd7e3rYnRAiRlpZmOfD9hz/8oWfPnqGhoUaj8fTp02lpaaWlpZ9++qnWxLbgZWrmJnc+CxstadGaCrQYAa2fTp06HT9+PCUlZeTIkY888kh8fLzlswSvvPKKp6dnt27dxo8fv2TJkmPHjlneR7LtqaeeSkhICA0N/fHHH7du3Wo9pCCE6NKly5YtWxYsWBASEpKSkrJ58+YmR+vVq9emTZt+85vfhISE3Lp1a82aNXf11GJjY7t3715ZWXlnqfnN2F6zyXGa/0ADBw5ctWrV5MmTBw8e3L9//9mzZ9ve9tFHH3333XdffPHF4cOHv/jii76+vg2rr7/+eteuXXv16jVs2LDw8PApU6bYnhAhxPTp03ft2hUaGpqcnLx161aDwdC3b99Zs2aNHDkyMjJyzpw51jUbjdOCl6mZm9z5LLRaskF1KnAvuGisU9qyZcvOnTu3bt3q6EYAtCL2oAFAUgQ0AEiKQxwAICn2oAFAUgQ0AEiKgAYASRHQACApAhoAJEVAA4CkCGgAkBQBDQCS+v/K5aTwNTRJmgAAAABJRU5ErkJggg=="
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "and the _p-value_ of"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(tunaref <= 29)/length(tunaref)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "[1] 0.02776"
       ]
      }
     ],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}