CamDavidsonPilon/findN.ipynb

## findN.ipynb
{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Suppose we have two groups, A,B, in our experiment, of approximately equal size, $N$. The groups have true conversion rates $p_A, p_B$ respectively. We use a beta-binomial model to find the posterior of $p$,  e.g. $p_A \\sim Beta(\\alpha=1 + c_A, \\beta= 1 + N - c_A)$, where $c_A$ is the number of conversions observed. For large $N$, this posterior is approximately (actually really close to being) Normally distributed, e.g. \n",
      "\n",
      "$$p_A \\sim Nor \\left( \\mu_A = \\frac{\\alpha}{\\alpha+\\beta}, \\sigma_A = \\frac{\\frac{\\alpha}{\\alpha+\\beta}\\frac{\\beta}{\\alpha + \\beta}}{(\\alpha+\\beta+1)}  \\right)$$\n",
      "\n",
      "Ultimately, we are interested in when $Pr( p_A > p_B \\;| \\;c_A, c_B, n ) \\ge 0.95 \\Rightarrow Pr( p_A - p_B > 0  \\;| \\;c_A, c_B, n ) \\ge 0.95$. As both $p_B$ and $p_A$ are normal, denoting $D = p_A - p_B$, then $D$ is Normal, $D \\;| \\;c_A, c_B, n  \\sim Nor\\left( \\mu = \\mu_A - \\mu_B, \\sigma^2 = \\sigma_A^2 + \\sigma_B^2  \\right)$. Suppose $\\mu_A > \\mu_B$ so that $\\mu > 0$. \n",
      "\n",
      "We'd like to find a $\\sigma^2 = \\sigma^2(N)$ (as $\\sigma_{A,B}$ are functions of the sample size $N$) s.t. $Pr(D > 0 \\;| \\;c_A, c_B, n ) \\ge 0.95 \\Rightarrow Pr(D < 0 \\;| \\;c_A, c_B, n ) \\le 0.05$.\n",
      "\n",
      "$$ Pr(\\frac{D - \\mu}{\\sigma} < \\frac{-\\mu}{\\sigma}) \\le 0.05 $$\n",
      "\n",
      "Inverting the normal CDF:\n",
      "\n",
      "$$\\frac{-\\mu}{\\sigma} \\le -1.65$$\n",
      "\n",
      "$$\\frac{\\mu^2}{1.65^2} \\ge \\sigma^2 = \\sigma_A^2 + \\sigma_B^2 $$\n",
      "\n",
      "\n",
      "$\\sigma_A^2$ can be approximated by $\\frac{\\hat{p}_A(1-\\hat{p}_A)}{N}$, so:\n",
      "\n",
      "$$ N \\ge \\frac{(1.65^2)\\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right)}{\\mu^2} $$\n",
      "\n",
      "\n",
      "where $\\mu = \\hat{p}_A - \\hat{p}_B$. Denote\n",
      "\n",
      "\n",
      "$$ N^* = \\frac{(1.65^2)\\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right)}{\\mu^2} $$\n",
      "\n",
      "\n",
      "It can be empirically shown that using the above formula that \n",
      "\n",
      "$$E[ 1_{Pr( D> 0 \\;| \\;C_A, C_B, N^*) \\ge 0.95 )} ] = 0.5$$\n",
      "\n",
      "i.e. there is a 50% chance that using $N^*$ will provide significance, a desirable property when we talk about *expected* sample size. Practically though, choosing $N^*$ should be a lower bound, and anything above it should be chosen. \n",
      "\n",
      "Thus the \"power\" (defined at the probabilty of correctly rejecting insignificance) of the test is 50%. Can we tweak the above formula to add more power?\n",
      "\n",
      "Of course by increasing $N$ we do this, but by how much? What is a good pre-determined power? Most pracitioners choose 80%, so let's try that first. \n",
      "\n",
      "Empirically, it seems like multiplying the $N^*$ by 2.25 achieves 80% power. \n",
      "\n",
      "### Unequal groups\n",
      "\n",
      "What about the case where $N_A \\ne N_B$? We get as far as:\n",
      "\n",
      "\n",
      "$$\\frac{\\mu^2}{1.65^2} \\ge \\left( \\frac{\\hat{p}_A(1-\\hat{p}_A)}{N_A}+ \\frac{\\hat{p}_B(1-\\hat{p}_B)}{N_B}\\right) $$\n",
      "\n",
      "$$\\frac{\\mu^2}{1.65^2} \\ge \\frac{ \\left( \\frac{\\hat{p}_A(1-\\hat{p}_A)}{N_A}+ \\frac{\\hat{p}_B(1-\\hat{p}_B)}{N_B}\\right) }\n",
      "                                 {\\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right)} \\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right) $$\n",
      "\n",
      "\n",
      "We can call the term:\n",
      "\n",
      "$$ \\frac{ \\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right) }\n",
      "                                 {\\left( \\frac{\\hat{p}_A(1-\\hat{p}_A)}{N_A}+ \\frac{\\hat{p}_B(1-\\hat{p}_B)}{N_B}\\right)} $$\n",
      "\n",
      "\n",
      "\n",
      "The *effective sample size*. So if we assign $N_A = N, N_B = rN_A$, where $0 < r \\le 1$, (so in the equal case $r=1$), then: \n",
      "\n",
      "$$ \\frac{ \\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right) }\n",
      "                                 {\\left( \\frac{\\hat{p}_A(1-\\hat{p}_A)}{N}+ \\frac{\\hat{p}_B(1-\\hat{p}_B)}{rN}\\right)} $$\n",
      "                                 \n",
      "$$ N \\frac{ \\left(\\hat{p}_A(1-\\hat{p}_A)+ \\hat{p}_B(1-\\hat{p}_B)\\right) }\n",
      "                                 {\\left( \\hat{p}_A(1-\\hat{p}_A) + \\frac{\\hat{p}_B(1-\\hat{p}_B)}{r}\\right)} $$\n",
      "\n",
      "Notice that by symmetry, $r<1$ (consider if $r>1$, then we could rewrite the above $N_B = N, N_A = rN_B$, in which case $r < 1$ again). It is easy to see that for every $r < 1$, the effective sample size is smaller than $N$. \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#tests\n",
      "\n",
      "needed = lambda p_A, p_B: (1.65)**2*( p_A*(1-p_A) + p_B*(1-p_B) )/(p_A - p_B)**2\n",
      "Var = lambda p,N: p*(1-p)/N"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import norm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figsize(12,6)\n",
      "x = np.linspace(-0.5,.5, 1000)\n",
      "\n",
      "\n",
      "p_A, p_B = 0.15,0.07\n",
      "n = N(p_A, p_B)\n",
      "print n\n",
      "mu = (p_A - p_B)\n",
      "var = np.sqrt( Var(p_A, n) + Var(p_B,n))\n",
      "plot(x,norm.pdf(x, loc=mu, scale= var ))\n",
      "print norm.cdf(0, loc=mu, scale= var )\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "81.930234375\n",
        "0.0494714680336"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAFwCAYAAABejXgsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4XfV97/vP2oPmybIt2dbg2WDZwkNsHMAJEAdSIEBJ\nfO5NoCUBTJ/WyU3J6U3Tcofk9NwmUNLTmuQ2eU5I+2SEnkvSMAScFh9MEIRRNhajjbGswfIgW/O0\np3X/kNa28bgl7b3X+u31fj1PnnbZW1tf/PFP/uqn7/oty7ZtWwAAAECOC7hdAAAAAJANNL4AAADw\nBRpfAAAA+AKNLwAAAHyBxhcAAAC+QOMLAAAAX5hU47tt2zY1NjZq5cqV2rZtW6ZqAgAAANIu5cb3\nzTff1EMPPaRXX31Vb7zxhp588knt378/k7UBAAAAaZNy4/vuu+9qw4YNKigoUDAY1JVXXqlf/epX\nmawNAAAASJuUG9+VK1fq+eef14kTJzQ8PKzf/OY36ujoyGRtAAAAQNqEUn3hxRdfrK9//eu69tpr\nVVxcrDVr1igQ4N44AAAAmMGybdueygfee++9qq+v15/+6Z8mf+0Xv/iFqqur01YcAAAAcLrBwUHd\nfPPNk/64lHd8Jeno0aOqqqpSW1ub/u3f/k0vv/zyh36/urpaa9eunXQRcN99992nv/qrv3K7DEwR\n+ZmL7MxGfmYjP3M1NzdP6eMm1fhu3rxZx48fVzgc1j/90z+prKxsSp8U3tPW1uZ2CZgG8jMX2ZmN\n/MxGfv4zqcb3d7/7XabqAAAAADKKu9MgSbr11lvdLgHTQH7mIjuzkZ/ZyM9/pnxz29ns2LGDGV8A\nAABkVHNzszZt2jTpj2PHF5KkpqYmt0vANJCfucjObORnNvLzHxpfAAAA+AKjDgAAADAKow4AAADA\nedD4QhJzTqYjP3ORndnIz2zk5z80vgAAAPAFZnwBAABgFGZ8AQAAgPOg8YUk5pxMR37mIjuzkZ/Z\nyM9/aHwBAADgC8z4AgAAwCjM+AIAAADnQeMLScw5mY78zEV2ZiM/s5Gf/9D4AgAAwBeY8QUAAIBR\nmPEFAAAAzoPGF5KYczId+ZmL7MxGfmYjP/+h8QUAAIAvMOMLAAAAozDjCwAAAJwHjS8kMedkOvIz\nF9mZjfzMRn7+E3K7AACAGSKxhPZ2DytgWVo6q1DhIHsnAMzCjC8A4IJ+d6BHP/h9p7qHo5KkWUVh\nfe3K+VpTU+pyZQD8iBlfAEBGbH/vuP6fHa3qHo6qtjxfteX56h6O6q+3v6/nD/S6XR4ApIzGF5KY\nczId+ZnL69m9e3RI/9jUJkm6Y91cPbR5uR7avFybG6uUsKUHnjuotp5Rl6t0j9fzw/mRn//Q+AIA\nzioST+iB5w4qYUu3rJytz6+eo4BlKWBZuvvSebp68QyNxhL6b8+3KY1TcwCQMZNqfL/97W9rxYoV\namxs1K233qqxsbFM1YUs27hxo9slYBrIz1xezu6pd4+rvW9MteX5unPdvA/9nmVZ+soVdZpRGNLb\nR4e084Mel6p0l5fzw4WRn/+k3Pi2trbqhz/8oZqbm9XS0qJ4PK5HHnkkk7UBAFwyGkvo4d2HJUl3\nrZ+n/NCZ/1wU5wX1xYmG+MevH1Y8wa4vAG9LufEtKytTOBzW8PCwYrGYhoeHVVNTk8nakEXMOZmN\n/Mzl1ez+Y+9x9YzEtHRWoS6fX37O1127tFLVJXk61D+m37f1ZbFCb/BqfkgN+flPyo1vZWWl/uIv\n/kL19fWaN2+eKioq9MlPfjKTtQEAXGDbth5/u1uS9J8aq2VZ1jlfGwxY+mxjlSTply1Hs1IfAExV\nyo3v/v379Y//+I9qbW3VoUOHNDg4qJ///OeZrA1ZxJyT2cjPXF7Mbk/XoA72jqqyKKQrFpx7t9fx\nqWWVKgoH9NaRId+d8ODF/JA68vOflJ/c9tprr+nyyy/XzJkzJUmf+cxn9OKLL+q222770Ou2bt2q\n+vp6SVJ5ebkaGxuTf7GcHylwzTXXXHPt3etn3j+h/v27demSGQoHGy/4+sJwUPOH3tfL7f36931V\n2nJpjaf+e7jmmmvzr1taWtTXNz5O1dbWpi1btmgqUn5y2xtvvKHbbrtNr776qgoKCvTFL35Rl156\nqb70pS8lX8OT28zV1NSU/AsG85CfubyWXSSe0Od+/qYGI3E99Nnlqp9RkNLHvX1kSPc8sVeVhSH9\n/PMrFQycezwil3gtP0wO+Zkr409uW7VqlW6//XatW7dOl1xyiSTpT/7kTyb9CQEA3vV6x4AGI3Et\nqixMuemVpOVVRZpbmqcTIzG9c3QogxUCwNSlvOObCnZ8AcBs3362Vc/u79Gd6+fqc6vmTOpj//vL\nnXq05ag+s3K2/vSjtRmqEACysOMLAMhtY7GEfn9wfIbuykUzJv3xzo1wL7T28SQ3AJ5E4wtJJwfJ\nYSbyM5eXsnuja0CjsYSWzCzU3NL8SX/88qpiVRaGdGQwog9OjGSgQu/xUn6YPPLzHxpfAIAk6bWO\nAUnSpXVlU/r4gGXpsvknd30BwGtofCFJ3NVqOPIzl5eye7W9X5K0foqNr6Rk4/tqR39aavI6L+WH\nySM//6HxBQDoUP+YOvvHVJIX1MWzi6f8Po1zShQOWNp7bFj9o7E0VggA00fjC0nMOZmO/Mzlleyc\n3d6P1JRO6wzewnBQDdXFsiW90TWYpuq8yyv5YWrIz39ofAEAeq1j+mMOjjXzSiVJuzoHpv1eAJBO\nNL6QxJyT6cjPXF7ILp6w1XJ4fHd2TU3ptN9v7cR7NB/K/TlfL+SHqSM//6HxBQCf239iRMPRhOaW\n5ml2cd6032/prCIV5wV1qD+iroGxNFQIAOlB4wtJzDmZjvzM5YXs9kzM4l4ytyQt7xcMWFo98V5v\nHMrtOV8v5IepIz//ofEFAJ9zxhwa56Sn8ZWklRPv9daR3G58AZiFxheSmHMyHfmZy+3sEratN53G\nN007vpK0cs74kWhvHRlK23t6kdv5YXrIz39ofAHAxw72jGpgLK7ZxWHNKZn+fK9j8cwi5YcC6ugb\nU+9ING3vCwDTQeMLScw5mY78zOV2dqeOOVjW1M/vPV0oYOni2UWScnvX1+38MD3k5z80vgDgY+8c\nHW9KV1RP/Wlt5+K8Zy43vgDMQuMLScw5mY78zOV2du8cHZYkXVyV/sbXDze4uZ0fpof8/IfGFwB8\nqn80pkP9Y8oLWlpYWZj2919eVSxL0r7uEUViibS/PwBMFo0vJDHnZDryM5eb2b17bHwEYdmsIoUC\n6ZvvdRTnBVU/o0CxhK39J0bS/v5ewNozG/n5D40vAPjUuxkcc3BcNGv8Brf3jg1n7HMAQKpofCGJ\nOSfTkZ+53MzO2fG9uKooY59j2cTJDnuP5eYNbqw9s5Gf/9D4AoAPJWz75I7v7Azu+M5mxxeAd9D4\nQhJzTqYjP3O5lV1X/5gGI3FVFoU0uzicsc+zsLJQoYCljr4xDUXiGfs8bmHtmY38/IfGFwB8aF/3\n+M1mS2cWpfXBFafLCwa0qLJQtqR93ez6AnAXjS8kMedkOvIzl1vZvX98vAldMitz872OZTk87sDa\nMxv5+Q+NLwD40PvHx3d8l8xM//m9p7s4hxtfAGah8YUk5pxMR37mciM727b1/sTYwZKZmd/xdT7H\nBydyr/Fl7ZmN/PyHxhcAfObYUFT9Y3GV5gdVVZK5G9scdRX5CgcsHeqP5OQNbgDMQeMLScw5mY78\nzOVGdsn53pmFGb2xzREOBjR/RoEk6UCOPcGNtWc28vMfGl8A8Jn3u5353syPOTgWT8wS7z+eW40v\nALOk3Pi+9957WrNmTfJ/5eXlevDBBzNZG7KIOSezkZ+53Mgumyc6OBZV5mbjy9ozG/n5TyjVF150\n0UXatWuXJCmRSKimpka33HJLxgoDAGRGNk90cCye2F3en4M3uAEwx5RGHZ555hktXrxYdXV16a4H\nLmHOyWzkZ65sZ9czElX3UFSF4YBqyvOz9nmdUYfWnlHFEnbWPm+msfbMRn7+M6XG95FHHtGtt96a\n7loAABnmjBosqixUIAs3tjmK84KaU5qnaNxWe+9o1j4vAJwq5VEHRyQS0RNPPKH777//rL+/detW\n1dfXS5LKy8vV2NiY/I7KmaXh2nvXp845eaEersnPL9fOr2Xr8x0qWypJCnS+qaamo1n9783vOiSV\nLNX+4yPqfPv1rPz35lp+XJOfX69bWlrU19cnSWpra9OWLVs0FZZt25P6mdNjjz2m73//+9q+ffsZ\nv7djxw6tXbt2SoXAXU1NTcm/YDAP+Zkr29n93c5WPfN+j75yRZ0+vXxW1j6vJP2suUs/aT6szY1V\n+pMNNVn93JnC2jMb+ZmrublZmzZtmvTHTXrU4eGHH9bnP//5SX8ieBsL32zkZ65sZ3egZ3zMYGFl\nQVY/r3TKDW7Hc+cGN9ae2cjPfybV+A4NDemZZ57RZz7zmUzVAwDIkHjCVtvEfO2CGdk70cFx6lm+\nk/xhIwCkxaQa3+LiYnV3d6u0tDRT9cAlp847wTzkZ65sZtfZP6Zo3FZVSVjFecGsfV7H7OKwisIB\n9Y/F1TsSy/rnzwTWntnIz394chsA+ETrxOOCF7qw2ytJlmUld5pbezjZAUD20fhCEnNOpiM/c2Uz\nO2e+d0GlO43v+Oceny1u7cmNJ7ix9sxGfv5D4wsAPnEgueOb/RvbHOz4AnATjS8kMedkOvIzVzaz\na02e6ODeju/8iab7YI40vqw9s5Gf/9D4AoAPjETj6uofU9CSarP4qOLTLZhxctSBkx0AZBuNLyQx\n52Q68jNXtrJr6x2VLam2okDhoHtf+mcUhlVeENJwNKFjQ1HX6kgX1p7ZyM9/aHwBwAcOnJgYc3Bx\nvtdx6q4vAGQTjS8kMedkOvIzV7ayc5pMN+d7HScbX/PnfFl7ZiM//6HxBQAfcJ7YNt8DO77zOdkB\ngEtofCGJOSfTkZ+5spWdc4pCfYX7je/C5MkO5o86sPbMRn7+Q+MLADluJBrXsaGowgFLc0vdO9HB\n4ew6t/WMKp7gZAcA2UPjC0nMOZmO/MyVjezae8ckSTXl+QoGrIx/vgspyQ9pVlFYY3Fbhwcibpcz\nLaw9s5Gf/9D4AkCOO9g7PlLghTEHR72z69vLnC+A7KHxhSTmnExHfubKRnZtEzu+Xmp868rHa2nv\nM7vxZe2Zjfz8h8YXAHKcs6vqpca3vmJ81ridHV8AWUTjC0nMOZmO/MyVjezaPHSig8OpxfRRB9ae\n2cjPf2h8ASCHReIJdQ2MKWBJteXun+jgqJtofNt7x2TbnOwAIDtofCGJOSfTkZ+5Mp1dZ9+YErY0\npzRfeSHvfMmfURhSSV5Qg5G4ekZibpczZaw9s5Gf/3jnqyAAIO1Ozvd6Z7dXkizLyplxBwDmoPGF\nJOacTEd+5sp0dslHFXtovtdRlwM3uLH2zEZ+/kPjCwA5zLmxrc6Tja+z4zvmciUA/ILGF5KYczId\n+Zkr09kld3xneLDxzYGzfFl7ZiM//6HxBYAcFU/Y6ugf3011mkwvYcYXQLbR+EISc06mIz9zZTK7\nwwNjisZtzSoOqygvmLHPM1VzSvMUDljqHopqOBJ3u5wpYe2Zjfz8h8YXAHKUMzvrxRvbJCkYsFQz\ncbZwRx9zvgAyj8YXkphzMh35mSuT2Tmzs7UeHHNwmD7uwNozG/n5D40vAOSozold1DqPneF7KtMb\nXwBmofGFJOacTEd+5spkds74QE2Zdxtf08/yZe2Zjfz8Z1KNb29vrzZv3qzly5eroaFBL730Uqbq\nAgBMUyejDgDwIaHJvPjP//zPdf311+vRRx9VLBbT0NBQpupCljHnZDbyM1emshuKxHViJKZw0NLs\nknBGPkc61Ew05Yf6xxRP2AoGLJcrmhzWntnIz39S3vHt6+vT888/rzvvvFOSFAqFVF5enrHCAABT\n19l/cswhYHm3mSwIBVRVElbcHj9+DQAyKeXG98CBA5o9e7buuOMOrV27VnfffbeGh4czWRuyiDkn\ns5GfuTKV3ckxB+/O9zpqysZ3fU080oy1Zzby85+UG99YLKbm5mZt3bpVzc3NKi4u1n333ZfJ2gAA\nU5S8sc3D870OpzlvN7DxBWCWlGd8a2trVVtbq/Xr10uSNm/efNbGd+vWraqvr5cklZeXq7GxMTlD\n43xnxbX3rjdu3OipergmP66nd/3iCy+o/9CAaj9+kyfqOd91bXm++vfv1gvRD7S58bOu18M111x7\n77qlpUV9fX2SpLa2Nm3ZskVTYdm2baf64o9//ON66KGHtGzZMn3zm9/UyMiI7r///uTv79ixQ2vX\nrp1SIQCA9Pnyr9/T3u5h/cOnl2rFnBK3yzmvV9v79X/8dr9WzS3RAzcsdbscAAZobm7Wpk2bJv1x\nkzrO7Lvf/a5uu+02rVq1Snv27NG999476U8Ib3K+u4KZyM9cmcjOtm11TMz41hgw41tr8GOLWXtm\nIz//CU3mxatWrdKrr76aqVoAAGnQOxLTcDShkrygygsm9WXeFVUleQoHLB0fjmokGldhOOh2SQBy\nFE9ug6STczQwE/mZKxPZJY8yK8+X5eGjzBzBgKV5E0+X6zRs15e1Zzby8x8aXwDIMc7IgAlHmTlq\nONkBQBbQ+EISc06mIz9zZSK7zuR8r/ePMnPUlTs7vmY9upi1Zzby8x8aXwDIMckd3zKTdnzNfYgF\nAHPQ+EISc06mIz9zZSK7jn7zRh1MPdmBtWc28vMfGl8AyCHxhK1Dp9zcZoqTje+oJnG8PABMCo0v\nJDHnZDryM1e6szs2FFE0bquyKGTUsWDlBSGV5AU1HE2odyTmdjkpY+2Zjfz8h8YXAHLIyflec25s\nkyTLsjjZAUDG0fhCEnNOpiM/c6U7O+ccXJPGHBy1Bp7swNozG/n5D40vAOQQE8/wddRysgOADKPx\nhSTmnExHfuZKd3ad/eO7pbUGneHrMPFkB9ae2cjPf2h8ASCHdOTAqEOHQaMOAMxC4wtJzDmZjvzM\nlc7sIvGEjgxEFLCkuaV5aXvfbJk38cCNroGI4gkzjjRj7ZmN/PyHxhcAckRX/5hsSXNK8xQOmvfl\nvTAc1KzisGIJW4cHIm6XAyAHmfeVERnBnJPZyM9c6cwuOeZg2FFmp0qe7NBvxrgDa89s5Oc/NL4A\nkCM6DT7RwcHJDgAyicYXkphzMh35mSud2Zl8Y5sjeYNbrxmNL2vPbOTnPzS+AJAjOvtzYcd3ovE1\nZNQBgFlofCGJOSfTkZ+50pmd88QzE8/wdSRHHQzZ8WXtmY38/IfGFwBywFAkrhMjMeUFLc0qDrtd\nzpRVl+QpFLDUPRzVSDTudjkAcgyNLyQx52Q68jNXurJzxhxqyvIVsKy0vKcbggFLcybOID7U7/1d\nX9ae2cjPf2h8ASAHOGMONQaPOTiSR5pxsgOANKPxhSTmnExHfuZKV3YdOXCUmcOkI81Ye2YjP/+h\n8QWAHJBLjW9N8iEW3m98AZiFxheSmHMyHfmZK20zvjlwhq+jtsycUQfWntnIz39ofAHAcLZtq8OZ\n8S3LgcY3OerAWb4A0ovGF5KYczId+ZkrHdn1jsQ0HE2oJC+o8oJQGqpyV2VRSAWhgPrH4uofjbld\nznmx9sxGfv5D4wsAhuvoPznmYBl8lJnDsqyTJzsw5wsgjSa1NbBgwQKVlZUpGAwqHA7rlVdeyVRd\nyDLmnMxGfuZKR3a5dGObo6Y8X+8fH1FH36iWVxW7Xc45sfbMRn7+M6nG17Is7dy5U5WVlZmqBwAw\nSbl0hq/DmVU24UgzAOaY9KiDbduZqAMuY87JbORnrnRkl9zxzYEb2xzODW5eP9mBtWc28vOfSTW+\nlmXpk5/8pNatW6cf/vCHmaoJADAJnTk66iAx4wsgvSY16vDCCy9o7ty5OnbsmK655hpdfPHF+tjH\nPpap2pBFzDmZjfzMNd3s4glbh/pz5wxfx6mjDrZte/amPdae2cjPfybV+M6dO1eSNHv2bN1yyy16\n5ZVXzmh8t27dqvr6eklSeXm5Ghsbk3+xnB8pcM0111xznZ7rJavWK5qwZXW8qddfHnK9nnRd73nt\nJSXaP9BYXaOOD0f17q5XPFUf11xznd3rlpYW9fX1SZLa2tq0ZcsWTYVlpzi0Ozw8rHg8rtLSUg0N\nDenaa6/VN77xDV177bXJ1+zYsUNr166dUiFwV1NTU/IvGMxDfuaabnavdfTr3u37tWpuiR64YWka\nK3PfPY/v1dtHh/R31y/R6nmlbpdzVqw9s5GfuZqbm7Vp06ZJf1wo1RceOXJEt9xyiyQpFovptttu\n+1DTCwDIvlx6VPHpasvz9fbRIXX0jXm28QVglpQb34ULF2r37t2ZrAUu4jtes5GfuaabXS6e6OBI\n3uDm4UcXs/bMRn7+w5PbAMBgnf25d4avwznSjLN8AaQLjS8knRwkh5nIz1zTzS4Xn9rmcE528PKR\nZqw9s5Gf/9D4AoChIvGEjgxEFLCkOaV5bpeTdvMmmvmu/jHFEzw8CcD00fhCEnNOpiM/c00nu67+\nMdmS5pTmKxzMvS/nBaGAZheHFbelwwMRt8s5K9ae2cjPf3LvKyUA+EQujzk4apNPcPPuDW4AzEHj\nC0nMOZmO/Mw1nexy+SgzR43Hb3Bj7ZmN/PyHxhcADJXLR5k5nB1frza+AMxC4wtJzDmZjvzMNZ3s\nOpJHmeV+4+vVs3xZe2YjP/+h8QUAQ3UmZ3xz7wxfR02Zt0cdAJiFxheSmHMyHfmZa6rZDUXi6hmJ\nKS9oaVZxOM1Vecec0jwFLenYUFRjsYTb5ZyBtWc28vMfGl8AMFDyxrayfAUsy+VqMicYsDR3Yob5\nkIcfZAHADDS+kMSck+nIz1xTza6jL/fnex3OE9y8OO7A2jMb+fkPjS8AGKjDB/O9jpMnO3jzBjcA\n5qDxhSTmnExHfuaaanZOE5jLD69wOGf5dnpwx5e1Zzby8x8aXwAwkD93fL3X+AIwC40vJDHnZDry\nM9dUsrNtW539uf+4YsfJxxZ7r/Fl7ZmN/PyHxhcADHNiOKaRaEJl+UGVFYTcLifjZhaFlR8KqG80\npoGxmNvlADAYjS8kMedkOvIz11SyOznfm/tjDpJkWdYpT3Dz1q4va89s5Oc/NL4AYJgOH405OLx8\npBkAc9D4QhJzTqYjP3NNJbuOXv+c4euo8eicL2vPbOTnPzS+AGAYP53o4OAsXwDpQOMLScw5mY78\nzDW1GV//jTrUevQsX9ae2cjPf2h8AcAgsYStroExWZLmlfmn8T11xte2bZerAWAqGl9IYs7JdORn\nrslmd3hgTAlbqirJU37IP1/CywpCKssPajSW0Ilh7xxpxtozG/n5j3++agJADmjv9d+Yg8MZd2DO\nF8BU0fhCEnNOpiM/c002u87kGb7+a3y9eLIDa89s5Oc/NL4AYBDnDN8aH53o4Dh5soN3Gl8AZqHx\nhSTmnExHfuaabHYdPh51cG5w89LJDqw9s5Gf/0yq8Y3H41qzZo1uvPHGTNUDADiPjn5GHZjxBTBV\nk2p8t23bpoaGBlmWlal64BLmnMxGfuaaTHbDkbhODMcUDlqaXZyXwaq8yTm+rWsgonjCG0easfbM\nRn7+k3Lj29HRoaeeekpbtmzhDEUAcEFyvrcsX8GA/zYgCsNBzSoOK5awdWQw4nY5AAyUcuP71a9+\nVQ888IACAcaCcxFzTmYjP3NNJjs/n+jg8Nqji1l7ZiM//0mpi33yySdVVVWlNWvWsNsLAC45+ahi\n/53o4Kgt8+ajiwGYIZTKi1588UU9/vjjeuqppzQ6Oqr+/n7dfvvt+slPfnLGa7du3ar6+npJUnl5\nuRobG5PfUTmzNFx77/rUOScv1MM1+fnl2vm1VF7/4q7DUuFi1Zbne6b+bF/XlC9LXs/urXK9HufX\nvPLnwzX55ep1S0uL+vr6JEltbW3asmWLpsKyJ7mF+9xzz+k73/mOnnjiiTN+b8eOHVq7du2UCoG7\nmpqakn/BYB7yM9dksvvSr9/Vvu4R/cONS7WiuiTDlXnTy219+r/+/QOtrSnVfdctcbsc1p7hyM9c\nzc3N2rRp06Q/bkoDu5zqkHtY+GYjP3Olmp1t28lRhzo/jzqUe+ssX9ae2cjPf0KT/YArr7xSV155\nZSZqAQCcw4mRmEaiCZXmB1VWMOkv3TmjujRfAUs6OhjRWCyh/BA3XANIHV8xIOnD804wD/mZK9Xs\nONFhXChgaW5pvmxJh/rd3/Vl7ZmN/PyHxhcADNDOiQ5JXht3AGAOGl9IYs7JdORnrlSz60w2vv7e\n8ZVOeXRxv/tn+bL2zEZ+/kPjCwAGcB7YUEPjm9z1ZscXwGTR+EISc06mIz9zpZpdey8nOjiSO74e\naHxZe2YjP/+h8QUAj4vEE+oaGFPAkmrK2PGtm2h823vdH3UAYBYaX0hizsl05GeuVLI71D+mhC3N\nKc1THsd3aWZRWEXhgPrH4uodibpaC2vPbOTnP3wFBQCPa5vY2WTMYZxlWaqrGP+zaOt1f9wBgDlo\nfCGJOSfTkZ+5UsnOae7qK2h8HfXJxtfdcQfWntnIz39ofAHA45xZ1voZNL6OugrmfAFMHo0vJDHn\nZDryM1cq2bUz6nAGZ8e3vc/dxpe1Zzby8x8aXwDwsIRtn9zxreBEB4dXRh0AmIXGF5KYczId+Znr\nQtkdHYxoLG6rsjCkkvxQlqryvrml+QoHLB0djGokGnetDtae2cjPf2h8AcDDkic6cGPbhwQDluY5\n5/l64EEWAMxA4wtJzDmZjvzMdaHskk9so/E9Q3Lcoce9cQfWntnIz39ofAHAw9qS8700vqdL3uDG\nnC+AFNH4QhJzTqYjP3NdKDtubDs358/EzRvcWHtmIz//ofEFAA9jxvfcnOPdmPEFkCoaX0hizsl0\n5Geu82XXNxpT/1hcReGAZhWFs1iVGWorCmRJ6uwbVSxhu1IDa89s5Oc/NL4A4FGn7vZaluVyNd5T\nEAqoqiRPcVs61M+uL4ALo/GFJOacTEd+5jpfdow5XJjbD7Jg7ZmN/PyHxhcAPOrko4q5se1cnBvc\nONkBQCqU9NuCAAAe6UlEQVRofCGJOSfTkZ+5zpcdR5ldmNs7vqw9s5Gf/9D4AoBHOQ+voPE9N7cb\nXwBmofGFJOacTEd+5jpXdiPRuI4MRhQKWJpbxqjDudQlH2IxJtvO/skOrD2zkZ//0PgCgAd1TJxN\nO68sX6EAJzqcS1lBSOUFIY3GEjo2FHW7HAAeR+MLScw5mY78zHWu7Fp7RiRJC2cw5nAhCyb+jA72\nZH/cgbVnNvLzHxpfAPCg1hPjTdx8Gt8Lchpf55sFADiXlBvf0dFRbdiwQatXr1ZDQ4P++q//OpN1\nIcuYczIb+ZnrXNm1TuxeLphRmM1yjDR/4s+o1YUdX9ae2cjPf0KpvrCgoEDPPvusioqKFIvFtHHj\nRjU1NfFjAgDIAGf3ckElO74XspAdXwApmtSoQ1FRkSQpEokoHo+rsrIyI0Uh+/gGxmzkZ66zZTcU\nievYUFThoKW5pZzocCHOOEhbz6jiieye7MDaMxv5+c+kGt9EIqHVq1erurpaV199tRoaGjJVFwD4\nlnOT1vyKAgU50eGCSvJDmlUU1ljc1uGBiNvlAPCwSTW+gUBAu3fvVkdHh373u99p586dGSoL2cac\nk9nIz1xny875kT03tqXOGQnJ9rgDa89s5Oc/Kc/4nqq8vFw33HCDXnvtNV111VUf+r2tW7eqvr4+\n+brGxsbkjxKcv2Bcc80111x/+B/cU3+/tWdU/ft3KxKeJWmBp+r16nW8/U31f9Cjgx+ZqysWuJsf\n1+ZcO7xSD9fnvm5paVFfX58kqa2tTVu2bNFUWHaKj7rp7u5WKBRSRUWFRkZG9KlPfUrf+MY3tGnT\npuRrduzYobVr106pEADAuL98ap92HxrUf712kTbUl7tdjhH+fe9xfed3bbpqUYXu/cRCt8sBkGHN\nzc0f6kFTFUr1hV1dXfrCF76gRCKhRCKhP/7jP57SJwQAnN9BjjKbtAUuHmkGwBwpz/g2NjaqublZ\nu3fv1p49e/S1r30tk3Uhy07/sQ/MQn7mOj273pGoekZiKgwHVFUSdqkq89RV5MuS1N47qmg8kbXP\ny9ozG/n5D09uAwAPObnbWyDL4kSHVBWGg5pTmqe4LXX2j7ldDgCPovGFpJMD5DAT+Znr9Oxak0eZ\nMeYwWclxhxPZG3dg7ZmN/PyHxhcAPIQntk3dAp7gBuACaHwhiTkn05GfuU7PrvWUUQdMzsmzfLO3\n48vaMxv5+Q+NLwB4hG3bJ5/axokOk8bJDgAuhMYXkphzMh35mevU7I4ORjUYiau8IKTKwpRPm8SE\n2vJ8hQKWuvrHNBKNZ+VzsvbMRn7+Q+MLAB6x/8SwJGnxzEJOdJiCcDCg+op82WLXF8DZ0fhCEnNO\npiM/c52a3QfHx2/KWlTJmMNULZpZJEnafzw7N7ix9sxGfv5D4wsAHuE0a4tn0vhO1eKJbxr2Hx92\nuRIAXkTjC0nMOZmO/Mx1anYfnGDHd7oWTXzT4PxZZhprz2zk5z80vgDgAUORuLoGIgoHLdVVcJTZ\nVDk7vh+cGFU8YbtcDQCvofGFJOacTEd+5nKyOzCxQzm/okChADe2TVVZQUizi8MaiyXUNZD5Rxez\n9sxGfv5D4wsAHsB8b/osSs758gQ3AB9G4wtJzDmZjvzM5WTHfG/6ON88ZKPxZe2Zjfz8h8YXADyA\nHd/0WTxxpFm2bnADYA4aX0hizsl05GeupqYmxRO2WnvY8U2XbI46sPbMRn7+Q+MLAC7r6BtVJG6r\nuiRPJfk8qni65pblqTAc0PHhqHpHom6XA8BDaHwhiTkn05GfuTZu3JjcmVzEmENaBCwrueub6XEH\n1p7ZyM9/aHwBwGVOc7aYMYe04WQHAGdD4wtJzDmZjvzM1dTUdHLHl8Y3bbJ1sgNrz2zk5z80vgDg\nItu2ta97WJK0dFaRy9XkjiUTJzu8z44vgFPQ+EISc06mIz9zLV19qfrH4iovCKmqJOx2OTljQeX4\nE/Dae0c1HIln7POw9sxGfv5D4wsALtqb3O0tlGXxqOJ0yQsGtKiyULaU3FEHABpfSGLOyXTkZ67t\nzzwnSVrGmEPaOX+mezPY+LL2zEZ+/kPjCwAuau8blSQtm03jm27On+neY+z4AhhH4wtJzDmZjvzM\nZNu2hqoaJLHjmwnZ2PFl7ZmN/PyHxhcAXNI1ENFgJK7KwpBmFnFjW7rNn1Gg/KClroGI+kdjbpcD\nwANofCGJOSfTkZ+Z9h4bVv/+3Vo6q4gb2zIgGLC0JMO7vqw9s5Gf/9D4AoBL9nJ+b8Ylxx2Y8wWg\nSTS+7e3tuvrqq7VixQqtXLlSDz74YCbrQpYx52Q28jPTvu5hlS1ezY1tGZS8wS1DO76sPbORn/+E\nUn1hOBzWP/zDP2j16tUaHBzURz7yEV1zzTVavnx5JusDgJyU4IltWcGOL4BTpbzjO2fOHK1evVqS\nVFJSouXLl+vQoUMZKwzZxZyT2cjPPAd7RjUcTSh46E1ubMugmvJ8FYUD6h6O6vhwNO3vz9ozG/n5\nz5RmfFtbW7Vr1y5t2LAh3fUAgC+8e3RIkjS/osDlSnJbwLKSO+rs+gJIedTBMTg4qM2bN2vbtm0q\nKSk54/e3bt2q+vp6SVJ5ebkaGxuTMzTOd1Zce+9648aNnqqHa/LL9eundzyn/o5+XXvbDZ6oJ5ev\nL5pdpOebmvQb66Auu+sPXa+Ha665nvx1S0uL+vr6JEltbW3asmWLpsKybdtO9cXRaFSf/vSndd11\n1+mee+454/d37NihtWvXTqkQAPCTux99Rwd7R7XtpmVaXlXsdjk5ram1V3/zzAGtmVei+69f6nY5\nANKgublZmzZtmvTHpTzqYNu27rrrLjU0NJy16YXZnO+uYCbyM8vgWEwHe0cVDlg6/O7rbpeT8xom\nvrF499iw4omU93pSwtozG/n5T8qN7wsvvKCf/exnevbZZ7VmzRqtWbNG27dvz2RtAJCT3p2YNV0y\nq1DhAMepZ1plUVhzS/M0Ek2otWfE7XIAuCiU6gs3btyoRCKRyVrgImeOBmYiP7O8M3Fj2/KqYm38\n6EUuV+MPDdXF6hqI6K0jQ1o8M33Hx7H2zEZ+/sNWAwBkmdP4NjDbmzXOn/XbR4ZcrgSAm2h8IYk5\nJ9ORnzkStq13j46POlxcVUx2WdJQPdH4Hk1v40t+ZiM//6HxBYAs6ugd02AkrllFYVWV5Lldjm8s\nmFGoonBAhwciGXmQBQAz0PhCEnNOpiM/c+w5PCjp5A4k2WVHMGDp4gyMO5Cf2cjPf2h8ASCLWiYa\n30vmnvkAIGTWyTnfQZcrAeAWGl9IYs7JdORnBtu21dI13nQ1zhlvfMkue5xd9jfTuONLfmYjP/+h\n8QWALDk8GFH3cFSl+UHNn1Hgdjm+01BVrIAl7ese1nAk7nY5AFxA4wtJzDmZjvzM4Oz2rpxTooBl\nSSK7bCrKC2rZrCIlbOmtNO36kp/ZyM9/aHwBIEuc+V5nzAHZt2pitnpP14DLlQBwA40vJDHnZDry\nM8PZbmwju+xaNa9UkrS7Kz03uJGf2cjPf2h8ASALuociOtQfUVE4oMWVhW6X41srqk/O+Q4x5wv4\nDo0vJDHnZDry876WU87vDQas5K+TXXYVhoO6aLYz5zv9XV/yMxv5+Q+NLwBkQXPn+EzpqrmlLleC\nSyYy2JOmcQcA5qDxhSTmnExHft5m23ay8V1b8+HGl+yyz7nB7Y00NL7kZzby8x8aXwDIsI6+MR0b\niqq8IKTFM5nvdduK6mIFJ+Z8B8dibpcDIItofCGJOSfTkZ+3Obu9a+adPL/XQXbZVxgOqqG6RAlb\n2nVoeru+5Gc28vMfGl8AyLDXO/slSWtrylyuBI51teMjJ6919LtcCYBsovGFJOacTEd+3hVL2Mmb\nqE6f75XIzi3rase/CXmto1+2bU/5fcjPbOTnPzS+AJBB7x4d0nA0obryfFWV5LldDiYsnlmo8oKQ\njg1F1dY76nY5ALKExheSmHMyHfl51+vJ0xzOPuZAdu4IWNYp4w5Tf3wx+ZmN/PyHxhcAMujV9vEZ\n0o/Ucn6v13yk5uS4AwB/oPGFJOacTEd+3nR8KKq93cPKD1paPe/sjS/Zucf5ZqTl8KDGYokpvQf5\nmY38/IfGFwAy5KX2PknjYw4FIb7ces2MwrCWzipUJG5r96GpjzsAMAdfiSGJOSfTkZ83vXRwvPH9\naP25jzEjO3ddVl8uSXpxIqvJIj+zkZ//0PgCQAaMRONqnthF3DDRXMF7Lp9fIUn6/cE+xRNTP9YM\ngBlofCGJOSfTkZ/37Do0oGjc1sWzi1RZFD7n68jOXQsrCzS3NE+9ozG9c3Ro0h9PfmYjP/+h8QWA\nDPh9csyB3V4vsyxLl8+f3rgDAHPQ+EISc06mIz9viSdsvdQ2fkTWhRpfsnPf5QvGxx1eaO2d9FPc\nyM9s5Oc/KTe+d955p6qrq9XY2JjJegDAeG90DahvNKba8nwtrCxwuxxcQENVscoLQuoaiKi1h6e4\nAbks5cb3jjvu0Pbt2zNZC1zEnJPZyM9bdu7vlSRduWiGLMs672vJzn3BgJU83eH5A72T+ljyMxv5\n+U/Kje/HPvYxzZgxI5O1AIDxovGEmlrHm6erFlW4XA1SdeVEVv9zf8+kxx0AmIMZX0hizsl05Ocd\nzZ0DGozEtXBGgebPKLzg68nOG1bPK1VlYUiH+sf03rHhlD+O/MxGfv5D4wsAabTzgx5J42MOMEcw\nYOnKxeOZ/c/9PS5XAyBTQul+w61bt6q+vl6SVF5ersbGxuR3VM4sDdfeuz51zskL9XBNfiZeR+MJ\nvXhwfFa06Og7amp6/4If7/yaF+r3+3VF76ikmdq5v0crowcUCFjkl+PXzq95pR6uz33d0tKivr7x\nIwfb2tq0ZcsWTYVlT2KYqbW1VTfeeKNaWlrO+vs7duzQ2rVrp1QI3NXU1JT8CwbzkJ83PLPvhP7u\nuYO6aHaRvnvzRSl9DNl5h23buuvRd9TRN6Zv/cFiras996OmHeRnNvIzV3NzszZt2jTpj0t51OHz\nn/+8Lr/8cu3du1d1dXX6l3/5l0l/MngXC99s5OcN2987Lkn61LKZKX8M2XmHZVn6xMS4wzP7TqT0\nMeRnNvLzn1CqL3z44YczWQcAGK2jb1R7Dg8qPxTQ1YuZ7zXVJ5dW6qfNh/V8a6+2jsZUVpDyP5MA\nDMDNbZD04XknmIf83Pfbid3eqxZVqDgvmPLHkZ23zCnN17raMkXjtv49hV1f8jMb+fkPjS8ATFMs\ncbJJ+oOLUh9zgDd9evksSdJv3unmTF8gx9D4QhJzTqYjP3e9dLBPPSMx1VcUqKGqeFIfS3bec2ld\nmWYVh9XZP6bdXYPnfS35mY38/IfGFwCm6d/eOiZJuuHimRd8RDG8LxiwdP3Ezv2T73S7XA2AdKLx\nhSTmnExHfu7Z1z2slsODKgoHJnWag4PsvOkPLpqpgCW90NqrwwNj53wd+ZmN/PyHxhcApuFXbx6V\nJF130UwVTeKmNnjbrOI8Xb14hhK29MuWo26XAyBNaHwhiTkn05GfO44ORvTcB70KWNLNK2ZP6T3I\nzrv+l0uqJUlPv3dcPSPRs76G/MxGfv5D4wsAU/TIG0cUS9i6ctEMzSnNd7scpNnCykJtqCtTJG7r\nsYk5bgBmo/GFJOacTEd+2XdsKKLfvndclqRbV1dP+X3Izts+t2o828ff7tZQJH7G75Of2cjPf2h8\nAWAK/scbRxRN2Pr4wgrNn1HodjnIkBVzStQ4p0SDkbj+vz1H3C4HwDTR+EISc06mI7/s6uwb1ZPv\ndI/v9q6ZM633Ijvvu3P9XEnjN7kdH/rwrC/5mY38/IfGFwAm6UevdiluS9cuq9TCSnZ7c92K6hJt\nXFCusbitnzR3uV0OgGmg8YUk5pxMR37Z89bhQTW19io/aOkLH5k77fcjOzPcsW6eApb0273H1doz\nkvx18jMb+fkPjS8ApCiWsLXthXZJ0uZLqjWrOM/lipAtdRUFuv7iWUrY0oNN7UrYttslAZgCGl9I\nYs7JdOSXHb9sOarWnlHNLc1L3u0/XWRnjjvWzdWMwpDePDKkp949Lon8TEd+/kPjCwApaO8d1c8m\n5jv/tyvqlB/iy6fflOaH9KXLaiVJD73SecaNbgC8j6/ckMSck+nIL7NiCVv37zyosbitTy6ZoXW1\nZWl7b7Izy8cWVuij9WUajib0988f1O+ef97tkjANrD//ofEFgAv4yetd2ts9rKqSsL50eZ3b5cBF\nlmXpK1fUqSw/qNc6BvTs/h63SwIwCTS+kMSck+nIL3NeaO3VI28cUcCS/vLK+SrOC6b1/cnOPLOK\n8/SXV82XJP0+Xqc9XQMuV4SpYv35D40vAJxDa8+IHnjuoCTpzvXzdMncUpcrgldcWleu/3VVtRK2\n9F93tOpQ/5jbJQFIAY0vJDHnZDryS7+jgxHd+/R+DUcTunJhhf5TY1VGPg/ZmeuLH5mrOX171Tca\n073b31fPCDe7mYb15z80vgBwmu6hiP7q6ffVPRzVyupi/e9XzpdlWW6XBY8JBizdvnaulsws1KH+\niP7P3+5X/2jM7bIAnIdl2+k7hXvHjh1au3Ztut4OALKuq39MX3/6fR0eiGhRZYG+c8NSleSH3C4L\nHtYzHNU9T+xV10BEC2YU6NvXLdHMorDbZQE5rbm5WZs2bZr0x7HjCwAT3js2pK8+uVeHByK6aHaR\n/u56ml5c2IyisL7z6aWqK89Xa8+o/vMTe9XWM+p2WQDOgsYXkphzMh35TY9t23rynW795yf26cRw\nTKvmluj+65aorCDzTS/Zmc3Jb3Zxnv7bjcu0bFaRugYi+vJj7+nZ/Sdcrg4XwvrzHxpfAL7WPRTR\nf3nmgB58oV3RhK2bG2bpW3+wWEVpPrYMua+8IKQHbliiTyyeodFYQt9+9qD+7rmD6mPuF/AMZnwB\n+FIkltAT73Trp81dGo4mVBgO6CtX1GnTkkq3S4PhnJ8g/ODlTkXjtsryg7pz/Txdu2ymQgFukgTS\nYaozvgyvAfCVkWhcz+w7oYffOKLuofHjpy6rL9eXr6jV7OI8l6tDLrAsSzc2zNbamlJte6Fduw8N\n6h+b2vU/9hzRravn6KpFM5QX4geugBsmtfK2b9+uiy++WEuXLtX999+fqZrgAuaczEZ+55ewbb11\nZFD/74vt+vwv3tR3X+xQ91BUiyoL9DfXLtI3r1noWtNLdmY7X3415QW6/7oluvfqBaotz9eh/oi+\n87s23frwm/rvL3dq77FhpfGHrpgC1p//pLzjG4/H9eUvf1nPPPOMampqtH79et10001avnx5JutD\nlrS0tPDoRoOR35mOD0e1p2tQLV2D+n1bn44Pn3y4wIrqYt2yYrY2LqxQwOXzecnObBfKz7IsXbV4\nhj62sEI73j+hX791TO8fH9GjLUf1aMtRVZfk6aP15Vo1t0Qr5xSropBj0LKJ9ec/KTe+r7zyipYs\nWaIFCxZIkj73uc/pscceo/HNEX19fW6XgGnwa36xhK3uoYiODkZ1dDCig72jOnBiRAdOjOjY0Ief\nolVdkqePLazQpiUztHhmkUsVn8mv2eWKVPMLBixdu2ymrllaqXePDWvH+yfU1NqrI4MRPfb2MT32\n9jFJ0pzSPC2cUagFlQWaX1Gg6pI8zS7J08yisILMB6cd689/Um58Ozs7VVdXl7yura3Vyy+/nJGi\nvKazb0xDkfg5f9/W+X9UdaGfZE33B10Xfv8LF3BsMKK3Dg+e67cv9OHn//0L/gdm9s9vmp/eiHw7\n+8f0SvvUvoBP+893mv+BtmzFE1IknlA0bisaTygStxVNJBSJ2YombI3FEhoci2lgLK6BsbgGI3EN\njMXUOxI7Z32F4YBWVBercU6JPlJTpqWzCnn6GlxnWZaWVxVreVWxtl5Wq3eODmlX54D2HB7UO0eG\ndHggosMDEf2+7cPrOWBJFYUhleaHVJofVGl+SGX5QRXlBZUXDCg/aCkvGFBeKKC8if8/FLAUCEgB\nTfxfy1LAOvl/LctScOLakqTTlod1+i9ISnUJnf6ys33c2d7/bL+USceHo9rbPZzdT5pFc0vzVMpZ\n5B+S8p+Gn//BeOiVTr1wMLe/K/zg5bfVsmCf22Vgij549R29t+gDt8vIOkvSrKKwqkryNLs4rNqK\nAi2sLNDCGYWaV5ZvxA5ZW1ub2yVgGqaTX8CytKK6RCuqSySN/wSjo29UB06MqvXEiDr6x3RsMKKj\nQxGdGI4l/4f0+eD3b+mN+vfcLiNj/u9NC7VxYYXbZXhKyseZvfTSS/rmN7+p7du3S5K+/e1vKxAI\n6Otf/3ryNY899phKSkoyUykAAAAgaXBwUDfffPOkPy7lxjcWi+miiy7Sjh07NG/ePF166aV6+OGH\nmfEFAACAEVIedQiFQvre976nT33qU4rH47rrrrtoegEAAGCMtD65DQAAAPCqKT865sSJE7rmmmu0\nbNkyXXvttert7T3r63p7e7V582YtX75cDQ0Neumll6ZcLNIn1fyk8TOc16xZoxtvvDGLFeJ8Usmv\nvb1dV199tVasWKGVK1fqwQcfdKFSOFJ5ANBXvvIVLV26VKtWrdKuXbuyXCHO50L5/fznP9eqVat0\nySWX6IorrtCePXtcqBJnk+rDt1599VWFQiH96le/ymJ1uJBU8tu5c6fWrFmjlStX6qqrrjr/G9pT\n9LWvfc2+//77bdu27fvuu8/++te/ftbX3X777faPfvQj27ZtOxqN2r29vVP9lEijVPOzbdv++7//\ne/vWW2+1b7zxxmyVhwtIJb+uri57165dtm3b9sDAgL1s2TL77bffzmqdGBeLxezFixfbBw4csCOR\niL1q1aozsvjNb35jX3fddbZt2/ZLL71kb9iwwY1ScRap5Pfiiy8m/317+umnyc8jUsnOed3VV19t\n33DDDfajjz7qQqU4m1Ty6+npsRsaGuz29nbbtm372LFj533PKe/4Pv744/rCF74gSfrCF76gX//6\n12e8pq+vT88//7zuvPNOSeNzwuXl5VP9lEijVPKTpI6ODj311FPasmULj9b0kFTymzNnjlavXi1J\nKikp0fLly3Xo0KGs1olxpz4AKBwOJx8AdKpTM92wYYN6e3t15MgRN8rFaVLJ77LLLkv++7ZhwwZ1\ndHS4USpOk0p2kvTd735Xmzdv1uzZs12oEueSSn6/+MUv9NnPfla1tbWSpFmzZp33Pafc+B45ckTV\n1dWSpOrq6rN+gT5w4IBmz56tO+64Q2vXrtXdd9+t4eHcPSjaJKnkJ0lf/epX9cADDygQmPJfFWRA\nqvk5WltbtWvXLm3YsCEb5eE0Z3sAUGdn5wVfQ/PkDankd6of/ehHuv7667NRGi4g1bX32GOP6c/+\n7M8k+fu5BV6TSn779u3TiRMndPXVV2vdunX66U9/et73PO+pDtdcc40OHz58xq//7d/+7YeuLcs6\n61+UWCym5uZmfe9739P69et1zz336L777tPf/M3fnLcopMd083vyySdVVVWlNWvWaOfOnZkqE+cw\n3fwcg4OD2rx5s7Zt28Y52y5J9R/S03+qwj/A3jCZHJ599ln98z//s1544YUMVoRUpZKd05tYliXb\ntvnppoekkl80GlVzc7N27Nih4eFhXXbZZfroRz+qpUuXnvX15218/+M//uOcv1ddXa3Dhw9rzpw5\n6urqUlVV1Rmvqa2tVW1trdavXy9J2rx5s+67774L/kcgPaab34svvqjHH39cTz31lEZHR9Xf36/b\nb79dP/nJTzJZNiZMNz9p/AvCZz/7Wf3RH/2R/vAP/zBTpeICampq1N7enrxub29P/ljuXK/p6OhQ\nTU1N1mrEuaWSnyTt2bNHd999t7Zv364ZM2Zks0ScQyrZvf766/rc5z4nSeru7tbTTz+tcDism266\nKau14kyp5FdXV6dZs2apsLBQhYWF+vjHP6433njjnI3vlH9+fdNNN+nHP/6xJOnHP/7xWf9RnTNn\njurq6rR3715J0jPPPKMVK1ZM9VMijVLJ71vf+pba29t14MABPfLII/rEJz5B0+sRqeRn27buuusu\nNTQ06J577sl2iTjFunXrtG/fPrW2tioSiehf//Vfz/hH9aabbkqur5deekkVFRXJcRa4K5X82tra\n9JnPfEY/+9nPtGTJEpcqxelSye6DDz7QgQMHdODAAW3evFnf//73aXo9IpX8br75ZjU1NSkej2t4\neFgvv/yyGhoazv2mU73T7vjx4/amTZvspUuX2tdcc43d09Nj27Ztd3Z22tdff33ydbt377bXrVtn\nX3LJJfYtt9zCqQ4ekWp+jp07d3Kqg4ekkt/zzz9vW5Zlr1q1yl69erW9evVq++mnn3azbF976qmn\n7GXLltmLFy+2v/Wtb9m2bds/+MEP7B/84AfJ13zpS1+yFy9ebF9yySX266+/7lapOIsL5XfXXXfZ\nlZWVybW2fv16N8vFKVJZe44vfvGL9i9/+ctsl4jzSCW/Bx54wG5oaLBXrlxpb9u27bzvxwMsAAAA\n4Avcqg8AAABfoPEFAACAL9D4AgAAwBdofAEAAOALNL4AAADwBRpfAAAA+AKNLwAAAHyBxhcAAAC+\n8P8DauRUUNCNOBoAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10d2b7310>"
       ]
      }
     ],
     "prompt_number": 98
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 4000\n",
      "p_A, p_B = 0.11,0.1\n",
      "N = needed(p_A, p_B)*2.25\n",
      "print N"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "11510.049375\n"
       ]
      }
     ],
     "prompt_number": 85
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import beta\n",
      "s = 0\n",
      "for i in range(n):\n",
      "    c_A = np.random.binomial(1, p_A, size=N)\n",
      "    c_B = np.random.binomial(1, p_B, size=N)\n",
      "    s+= (((beta.rvs(1 + c_B.sum(),1 + N - c_B.sum(), size=5000 ) - beta.rvs(1 + c_A.sum(),1 + N - c_A.sum(), size=5000 )\n",
      "    ) > 0 ).mean() > 0.95)\n",
      "    \n",
      "print 1.0*s/n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.8095\n"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(10):\n",
      "    v = np.random.random(size=2)\n",
      "    p_A = v.min()\n",
      "    p_B = v.max()\n",
      "    N = needed(p_A, p_B)*2.25\n",
      "    s = 0\n",
      "    for i in range(n):\n",
      "        c_A = np.random.binomial(1, p_A, size=N)\n",
      "        c_B = np.random.binomial(1, p_B, size=N)\n",
      "        s+= (((beta.rvs(1 + c_B.sum(),1 + N - c_B.sum(), size=5000 ) - beta.rvs(1 + c_A.sum(),1 + N - c_A.sum(), size=5000 )\n",
      "        ) > 0 ).mean() > 0.95)\n",
      "        \n",
      "    print \"Power: %.2f\"%(1.0*s/n), \"p_A: %.2f, p_B %.2f\"%(p_A,p_B), \"N: %d\"%int(N)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Power: 0.76 p_A: 0.10, p_B 0.44 N: 18\n",
        "Power: 0.04"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.16, p_B 0.92 N: 2\n",
        "Power: 0.78"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.27, p_B 0.52 N: 42\n",
        "Power: 0.74"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.42, p_B 0.78 N: 19\n",
        "Power: 0.70"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.18, p_B 0.58 N: 14\n",
        "Power: 0.00"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.02, p_B 0.72 N: 2\n",
        "Power: 0.75"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.18, p_B 0.55 N: 17\n",
        "Power: 0.77"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.22, p_B 0.47 N: 41\n",
        "Power: 0.72"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.37, p_B 0.74 N: 19\n",
        "Power: 0.79"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " p_A: 0.89, p_B 0.97 N: 135\n"
       ]
      }
     ],
     "prompt_number": 86
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#effective sample size.\n",
      "r = 0.5 \n",
      "N =\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}