fonnesbeck/Bayesian Power Analysis.ipynb

## Bayesian Power Analysis.ipynb
{
  "cells": [
    {
      "metadata": {
        "editable": true,
        "deletable": true
      },
      "cell_type": "markdown",
      "source": "We are looking for a sample size such that the probability that the effect size is greater than zero in 95%."
    },
    {
      "metadata": {
        "collapsed": true,
        "trusted": true,
        "editable": true,
        "deletable": true
      },
      "cell_type": "code",
      "source": "%matplotlib inline\nimport numpy as np\nimport pymc3 as pm\nimport matplotlib.pylab as plt",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "titers = np.arange(10, 320)",
      "execution_count": 185,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": true
      },
      "cell_type": "code",
      "source": "invlogit = lambda x: 1/(1+np.exp(-x))\nlogit = lambda p: np.log(p/(1-p))",
      "execution_count": 186,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Find parameters that result in neutralization of 80"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "β_1_true = np.array([3, -0.0323])\np1 = invlogit(β_1_true[0] + β_1_true[1]*titers)\n(logit(0.6) - β_1_true[0]) / β_1_true[1]",
      "execution_count": 187,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": "80.326157643710076"
          },
          "metadata": {},
          "execution_count": 187
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Find parameters that result in neutralization of 160"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "β_2_true = np.array([3, -0.0162])\np2 = invlogit(β_2_true[0] + β_2_true[1]*titers)\n(logit(0.6) - β_2_true[0]) / β_2_true[1]",
      "execution_count": 188,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": "160.15647480813803"
          },
          "metadata": {},
          "execution_count": 188
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "This is what the \"true\" curves look like for both groups"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "plt.plot(titers, p1, label='80')\nplt.plot(titers, p2, label='160')\nplt.ylim(0, 1)\nplt.xlabel('dilution')\nplt.ylabel('% neutralization')\nplt.legend()",
      "execution_count": 190,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": "<matplotlib.legend.Legend at 0x121af3ba8>"
          },
          "metadata": {},
          "execution_count": 190
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<matplotlib.figure.Figure at 0x1217bbb70>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FGXXwOHfSU8gISShhxJ6r6GIioAgiEoTaVIUxILo\na8Peu6Kf7VWUJlaKKIJIkSqotCA1QOhC6IROSOX5/piFN0JIlpDN7G7OfV1z7e7s7O4ZNuFknnIe\nMcaglFJKAfjYHYBSSin3oUlBKaXUBZoUlFJKXaBJQSml1AWaFJRSSl2gSUEppdQFLksKIjJORA6J\nyIbLPC8i8rGIbBORdSLS2FWxKKWUco4rrxTGAx1zeP5moJpjuxcY6cJYlFJKOcFlScEYsxg4msMh\nXYCvjWUZEC4iZVwVj1JKqdz52fjZ5YA9WR4nOvbtv/hAEbkX62qCIkWKNKlZs2aBBKiUUt5i1apV\nR4wxJXI7zs6kINnsy7bmhjFmFDAKIDY21sTFxbkyLqWU8joi8o8zx9k5+igRKJ/lcTSwz6ZYlFJK\nYW9SmA4McIxCagGcMMZc0nSklFKq4Lis+UhEJgCtgSgRSQReAvwBjDGfAzOBTsA2IBm421WxKKWU\nco7LkoIxpk8uzxvgQVd9vlJK5SQ9PZ3ExERSUlLsDiVfBQUFER0djb+/f55eb2dHs1JK2SYxMZHQ\n0FAqVaqESHbjXjyPMYakpCQSExOJiYnJ03tomQulVKGUkpJCZGSk1yQEABEhMjLyqq5+NCkopQot\nb0oI513tOWlSUEopdYEmBaWUstEHH3xAnTp1qFu3Ln369CElJYWdO3fSvHlzqlWrRq9evUhLSyuw\neApPUji0GTb8CEnbwWQ7cVoppQrU3r17+fjjj4mLi2PDhg1kZmYyceJEnnrqKR599FG2bt1K8eLF\nGTt2bIHFVHhGH238GRa9Zd0PLAZl6kPZhlCmIZRtBMVjwKfw5EillHvIyMjg7Nmz+Pv7k5ycTJky\nZViwYAHff/89AAMHDuTll1/mgQceKJB4Ck1S2FnrPrZmNqZN2F78D66DfWtg+SjITLUOCAyDMg2s\nRFG+ubUVLWlv0EqpAvHKL/Fs3HcyX9+zdtkwXrqtTo7HlCtXjieeeIIKFSoQHBzMTTfdRJMmTQgP\nD8fPz/rvOTo6mr179+ZrbDkpNEnh53WH+Wh+OiVCy3NXy+vp168ixQKBQ5tg/xorSexfA8u/gL8+\nsV5UPMaRIJpZtyVrgY+vreehlPIex44dY9q0aezcuZPw8HDuuOMOZs2adclxBTlKqtAkhUfaVaNp\npQi+WLydEXMS+HThNno3rcD9N1SnZOP60HiAdWBGKuxfC3uWW9v2BbBuovVcQChEx0KFayDmeigX\nC34B9p2UUipf5PYXvavMmzePmJgYSpSwKlp3796dv/76i+PHj5ORkYGfnx+JiYmULVu2wGIqNElB\nRLiuWhTXVYti0/6TjF68g6+W7uK75f/Qr0VF7r+hCiVCA8Ev0HFl0Ax4yOqUPv4P7FlhJYndy62+\niUVvgl8wVGgOMa2gUiurb8K30PyTKqWuUoUKFVi2bBnJyckEBwczf/58YmNjadOmDVOmTKF37958\n9dVXdOnSpcBiEuNhI3Hycz2Ff5LO8MmCbfz0dyIBfj4MbFmJB9tUJSwol5ohyUfhn79g1xLYuQQO\nxVv7A4paVxGVW0O19hBVHbxwcoxS3mDTpk3UqlXL7jB46aWXmDRpEn5+fjRq1IgxY8awd+9eevfu\nzdGjR2nUqBHffvstgYGBTr9nducmIquMMbG5vbZQJ4Xzdhw+zcfztzJt7T6KhwTwaPvq9GlaHj9f\nJ0cjnTnyvwSxczEkbbX2F6sAVW+0EkRMKwgMzde4lVJ55y5JwRU0KeSTDXtP8OqMjazYeZTqpYry\nwq21ub5arqvXXer4btg2D7bOg52/Q9pp8PGHCi2sBFH9ZihRPf9PQCnlNE0K2dOkcBFjDHPiD/LW\nrE38k5RM5wZleeHW2lZ/Q15kpMGeZbB1Lmyb/7+mpshqUOtWqHmb1RehcySUKlCaFLKnSeEyUjMy\nGbloO58t3E6Qvw/PdKpFr9jy+PhcZR/BiURImAWbfoFdf4DJhNCyUPMWa6t0HfjmrQ66Usp5mhSy\np0khF9sPn+a5qetZtuMozWIieP+OBpSPCMmfN08+ClvmwOYZ1lVExlkICreuIOr2sPohdF6EUi6h\nSSF7mhScYIzhh7hEXp2xEYCXbqtNjybR+TuhJC3ZmhOxaTps/tXqhyhSEup0tRJE+WY6kkmpfKRJ\nIXs6qN4JIkLPpuW5pkokj/+wluFT1jFv00He7FaPyKJ57Gu4WECIdYVQ61ZIP2tdQWz4EVZ9BStG\nWSOZ6t0ODfpAiRr585lKKXUR7d28AuUjQpgwpAXPdqrJws2H6fTxEuJ2Hc3/D/IPtq4Qen0Dw7dB\n18+t0Up/fgyfNoPRN0LcOEg5kf+frZQqMIMGDaJkyZLUrVv3X/s/+eQTatSoQZ06dXjyyScv7H/r\nrbeoWrUqNWrUYM6cOS6JSZPCFfL1Ee5tVYWpD7YkyN+X3qOWMWbJDlzWDBcUBg37QL8f4fHNcNMb\nkHYGZjwK71WHH++B7Qvh3DnXfL5SymXuuusuZs+e/a99CxcuZNq0aaxbt474+HieeOIJADZu3MjE\niROJj49n9uzZDB06lMzMzHyPSZNCHtUpW4zpw66jbc2SvP7rJoZ+9zenUtJd+6FFS0LLYTB0KQxZ\nAA3vhK2/wTdd4aP6sHgEnD7k2hiUUvmmVatWRERE/GvfyJEjefrppy/MYC5Z0qrWPG3aNHr37k1g\nYCAxMTFUrVqVFStW5HtM2qdwFYoF+/NF/yaMXrKDd2YnsPWzvxg3sCkVIvNpdNLliEC5JtbW4U1I\n+NXqe1jwOix6B2p3htjBULGldk4r5YxZT8OB9fn7nqXrwc1vX/HLtmzZwpIlS3juuecICgrivffe\no2nTpuzdu5cWLVpcOM5VJbX1SuEqiVjNSd8MbsbhU6l0+fQPVrqin+Fy/IOg7u0wcDoMWwXN7rVm\nU4/vBJ9dAytGQ0r+1olXSrlORkYGx44dY9myZYwYMYKePXtijMm2idoVJbX1SiGftKwSxc8PXsug\n8Su5c/Ry3upej9ubRBdsEFFVoeOb0PZ5a+TSyjEw8wmY9zI06g8t7ofilQo2JqU8QR7+oneV6Oho\nunfvjojQrFkzfHx8OHLkCNHR0ezZs+fCca4qqa1XCvkoJqoIU4e2pEnF4jz+w1o+mLvFdR3QOQkI\ngcb94b7frb6HmrfAytHwcSOYPBD2rCz4mJRSTunatSsLFiwArKaktLQ0oqKi6Ny5MxMnTiQ1NZWd\nO3eydetWmjVrlu+fr0khn4WHBPDVoGb0aBLNR/O38tL0eM6ds3GCYLkm0H0UPLIeWj4MOxbC2HYw\n9ibYOA3O5f/oBaWUc/r06cM111xDQkIC0dHRjB07lkGDBrFjxw7q1q17YT0FEaFOnTr07NmT2rVr\n07FjRz799FN8ffO/4oHOaHYRYwxvzdrMqMU76NygLO/3bIC/s6W4XSn1NKz5DpZ9Bsd2Wc1J1z0K\nDfrqKnKqUNEZzdlzg/+lvJOI8GynWjzVsSbT1+5jyNdxnE1zg7/KA4tC8/vgob+h59cQXBx++Q98\n3BCWj7JmUyulCi1NCi72QOsqvNmtHr9vOczd41e4R2IAq9Be7S4wZCH0+wnCK8Cs4fBhfWvmdOpp\nuyNUStlAk0IB6Nu8Ah/2asiKnUcZ/NVK90kMYM1jqHojDJoNd82EUnVg7gvwYV1rMlzqKbsjVMpl\nPK353BlXe06aFApIl4bleL9nA5buSGLI13GkpLtRYjiv0rUw4Ge4Zz6Ub25NhvuoASz9FNJT7I5O\nqXwVFBREUlKSVyUGYwxJSUkEBQXl+T20o7mATVmVyPApa7muahSjB8QS5O/G6yUkxsGC12DHIggr\nB62GQ6N+ugiQ8grp6ekkJiaSkuJdf/AEBQURHR2Nv/+/f091PQU3NjluD0/9uI4ba5bk835N8HOH\nUUk52bkY5r8GiSugeAy0ec6aRa1LiCrlMXT0kRvrGVueVzvXYd6mQzw7db37X77GtILBv0HfyRBQ\nFH66B7643loUSCnlVVyaFESko4gkiMg2EXk6m+criMhCEVktIutEpJMr43En/a+pxMM3VmNyXCLv\nzkmwO5zciUD1DnDfYugxzuqA/qYbfHcHHNpsd3RKqXzisqQgIr7Ap8DNQG2gj4jUvuiw54HJxphG\nQG/gM1fF444ebVeNvs0rMHLRdsYs2WF3OM7x8bGajoathPavwe7lMLIlzHgMTh+2Ozql1FVy5ZVC\nM2CbMWaHMSYNmAh0uegYA4Q57hcD9rkwHrcjIrzWpS431y3N679uYvpaDzp9v0C49mF4eDU0HQyr\nxsMnjeGPD3SkklIezJVJoRywJ8vjRMe+rF4G+olIIjATeCi7NxKRe0UkTkTiDh/2rr9GfX2ED3o1\npFlMBE/8sJa/dx+zO6QrUyQSOo2Aocug4rVWRdZPm0HC7FxfqpRyP65MCtkV+r64R7UPMN4YEw10\nAr4RkUtiMsaMMsbEGmNiS5Qo4YJQ7RXk78vn/ZpQplgQ934dR+KxZLtDunIlqkPfiTBgGvgFwYRe\n8H0vOLrT7siUUlfAlUkhESif5XE0lzYPDQYmAxhjlgJBQJQLY3JbEUUCGDuwKakZ5xg8Ps71S3u6\nSuXWcP8fVn/DziXwaXNY+JbWVFLKQ7gyKawEqolIjIgEYHUkT7/omN3AjQAiUgsrKXhX+9AVqFqy\nKJ/d2Zhth0/zn4lryLSz5PbV8Auw+hseioNat8Lvb1vJIWGW3ZEppXLhsqRgjMkAhgFzgE1Yo4zi\nReRVEensOOxxYIiIrAUmAHcZtx+071rXVyvBK53rsGDzId6d4+FDPcPKWsNXB/4C/sEwoTdM6AMn\n8n9dWaVU/tAZzW7quanr+W75bkbe2Zib65WxO5yrl5lureGw8C3w8YP2L0OTQTorWqkCojOaPdyL\nt9WmYflwnvhhLdsOeUGlUl9/uPY/MHQpRMfCr4/DlzfrxDel3IwmBTcV6OfLyH6NCQ7w5b5vVnE6\nNcPukPJHRAz0nwpdP4cjCfD5ddbVQ0aq3ZEppdCk4NbKFAvmkz6N2ZWUzPAf1rp/jSRniUDDPvDg\nSqjT1eqI/vx62LPS7siUKvQ0Kbi5a6pE8nTHmszacICxf3jZmP+iJeD2MXDnFEhPhnE3wdyX9KpB\nKRtpUvAA91wfQ/vapXhn9mbWJR63O5z8V609PPCXtVbDnx/CFzfAvtV2R6VUoaRJwQOICCN61KdE\n0UAemrDacye25SQoDDp/Yl01pByH0TdafQ2ZXniuSrkxTQoeIjwkgI/6NGLP0WSem7rBe/oXLlat\nvTVCqV4Pq69hdFs4GG93VEoVGpoUPEjTShE82q4609fu44dViXaH4zrBxaH7KOj1HZzabzUn/fkR\nnDtnd2RKeT1NCh5maJuqXFM5kpemxbP98Gm7w3GtWrda1VdrdIS5L8I3XeGkB5UXV8oDaVLwML4+\nwoe9GxLo78Njk9aQnunlfz0XiYKe38BtH0PiSmtBn00z7I5KKa+lScEDlQoL4s1u9VibeIJPF26z\nOxzXE4EmA+G+JRBeESbdCb/8B9LO2B2ZUl5Hk4KH6lSvDN0aleOTBdtYs8cLh6lmJ6oqDJ4L1z4C\nq75yDF1dY3dUSnkVTQoe7OXOdSgZGshjk9ZwNi3T7nAKhl8AtH8FBk63rhTGtIOln4G3jsZSqoBp\nUvBgxYL9ef+OBuw4coa3Z22yO5yCFdMKHvgTqt0Ec56BSf3grIctZaqUG9Kk4OFaVo1i0LUxfLX0\nH/7adsTucApWSAT0/g46vAlbZsMXrWDvKrujUsqjaVLwAk92rEFMVBGe+mkdyWleUk3VWSJwzYMw\naI7VhDS2Ayz7XJuTlMojTQpeIMjfl3dur8+eo2cZMSfB7nDsER0L9y2Gqu1g9lOO5qRC0gGvVD7S\npOAlmsVEMOCaioz/axer/jlqdzj2CImAPhPgpje0OUmpPNKk4EWe7FiTssWCeXLKOlLSC8lopIuJ\nQMthcPcsOJdpNSfFjdPmJKWcpEnBixQN9OOt7vXYfvgMnyzYanc49irfDO5fYo1SmvEoTB8G6Sl2\nR6WU29Ok4GVaVS9BjybRfP77DjbsPWF3OPYKiYA7f4BWw2H1t/BlRzi+x+6olHJruSYFEakuIqNF\n5DcRWXB+K4jgVN68cEttIooE8OSUdd5fGyk3Pr7Q9nno/T0kbYdRN8CORXZHpZTbcuZK4Qfgb+B5\nYHiWTbmpYiH+vNalDhv3n+TLP71sCc+8qnkLDFkARUrAN93gjw+1n0GpbDiTFDKMMSONMSuMMavO\nby6PTF2VDnVK065WST6Yu5W9x8/aHY57iKoG98yDWrfBvJdg8gBIPWV3VEq5FWeSwi8iMlREyohI\nxPnN5ZGpqyIivNy5DgAvT9eVyy4IDIU7voL2r8HmGdayn0nb7Y5KKbfhTFIYiNVc9BewyrHFuTIo\nlT+ii4fwSLtqzN14kN/iD9gdjvsQgWsfhv5T4cwha8lP7WdQCnAiKRhjYrLZKhdEcOrqDbouhhql\nQnl5ejxnUgtZCYzcVG4NQxZCaBn4pjssH6X9DKrQc2b0kb+IPCwiUxzbMBHxL4jg1NXz9/XhjW51\n2XcihY/mF/K5C9mJiIHBv1nVVmcNhxmPQEaa3VEpZRtnmo9GAk2AzxxbE8c+5SFiK0XQu2l5xv6x\nk037T9odjvsJCrOGrF73GKwab60FfaaQVZxVysGZpNDUGDPQGLPAsd0NNHV1YCp/PdWxJsWC/Xlu\n6nrOndMmkkv4+EC7l6D7GKte0ug2cGCD3VEpVeCcSQqZIlLl/AMRqQwU0sI6nqt4kQCe7VSLv3cf\n58e/E+0Ox33VvwPungmZ6TD2Jtg0w+6IlCpQziSF4cBCEVkkIr8DC4DHXRuWcoXujcrRuEI478ze\nzMmUdLvDcV/lmlgd0CVrwqQ7YfEI7YBWhYYzo4/mA9WAhx1bDWPMQlcHpvKfj4/wape6JJ1J46N5\n2umco7AycNevUK8nLHgdfn4AMlLtjkopl7tsUhCRto7b7sAtQFWgCnCLY5/yQHXLFaN30wqM/2sX\nWw7qbN4c+QdD91HQ5jlYO8Eqj5FcSNeqUIVGTlcKNzhub8tmu9WZNxeRjiKSICLbROTpyxzTU0Q2\niki8iHx/BbGrPBreoQZFA/14eXo8RptFciYCNzwJt4+FxJUwpp3OgFZeTXL7T0FEYowxO3Pbl83r\nfIEtQHsgEVgJ9DHGbMxyTDVgMtDWGHNMREoaYw7l9L6xsbEmLk4nVF+tr5fu4sVp8Xx2Z2M61Stj\ndzieYfcymNgXzDlrCGvFlnZHpJTTRGSVMSY2t+Oc6Wj+MZt9U5x4XTNgmzFmhzEmDZgIdLnomCHA\np8aYYwC5JQSVf/o2q0DN0qG88esmzqbpYDKnVGhhFdQLiYKvOsPaSXZHpFS+y6lPoaaI3A4UE5Hu\nWba7gCAn3rsckHVFk0THvqyqA9VF5E8RWSYiHS8Ty70iEicicYcPH3bio1Vu/Hx9eKVzHfYeP8vI\nRdvsDsdzRFSGe+ZaCWLqvbDwTR2ZpLxKTlcKNbD6DsL5d39CY6y/8HMj2ey7+LfHD2tkU2ugDzBG\nRMIveZExo4wxscaY2BIlSjjx0coZzStH0rlBWT5fvIPdScl2h+M5gotDv5+gYT/4/R34aYgu9am8\nht/lnjDGTAOmicg1xpileXjvRKB8lsfRwL5sjllmjEkHdopIAlaSWJmHz1N58GynWszbdJDXft3I\n6AG5Njeq8/wCoMt/IbIKzH/FWuaz93dQJMruyJS6Ks70KawWkQdF5DMRGXd+c+J1K4FqIhIjIgFA\nb2D6Rcf8DLQBEJEorOakHVcQv7pKpYsFMaxtVeZuPMgfW7XezxURgesfgzvGw/41OjJJeQVnksI3\nQGmgA/A71l/8uQ5wN8ZkAMOAOcAmYLIxJl5EXhWRzo7D5gBJIrIRWAgMN8YkXflpqKsx6NoYoosH\n8/qvG8nUukhXrk43GDgDUk/C2PawRy90ledyZkjqamNMIxFZZ4yp7yibvcQY06JgQvw3HZLqGjPX\n72fod3/zZrd69G1ewe5wPFPSdvj2djh1AHqMtdaFVspN5OeQ1PNFco6LSF2gGFDyaoJT7ufmuqVp\nVimC939L0LpIeRVZBQbPhVK1YVI/WDHa7oiUumLOJIVRIlIceB6rT2Aj8K5Lo1IFTkR44dbaHE1O\n49OFOkQ1z4qWgIG/QLUOMPMJmPsinDtnd1RKOc2ZgnhjjDHHjDGLjTGVjTEljTGfF0RwqmDViy7G\n7Y2j+fKPXfyTdMbucDxXQBFrJFLsYPjzI/jpHi2mpzyGM8txZorI2yIiWfb97dqwlF2Gd6iBn6/w\n1szNdofi2Xx84Zb3od3LsOFHaw3os8ftjkqpXDnTfBTvOO43EYlw7MtuYpryAqXCgnjghirMjj/A\nsh06EOyqiMB1j1qrue1ZDuM6WPMZlHJjziSFDGPMk8BoYImINOHSmcnKiwxpVZmyxYJ4bYYOUc0X\n9e+A/j/ByX3WXIb96+yOSKnLciYpCIAxZjLQE/gSqOzKoJS9gvx9eermmsTvO8mPq3TpznwR0woG\nzbGalb68GbbNtzsipbLlTFK45/wdY0w8cB3WCmzKi3VuUJZGFcIZ8VsCp1Mz7A7HO5SqbVVZLV4J\nvu8Jq7+zOyKlLpHrymtAxaxVUoF2wOkCiU7ZRkR48dbaHD6VqlVU81NYWbh7FlS6DqYNhUXvaJVV\n5VZcuvKa8myNKhSna8OyjF6yk8RjWkU13wSFQd8foH5vWPQm/PIfyNSrMeUecqqS+pLj9u6CC0e5\nmyc71mR2/AHenrWZ//ZtbHc43sMvALp9DsXKwZL3rdIYd3xpzXFQykaXTQoi8lhOLzTG/F/+h6Pc\nTdnwYO5tVYWP52/l7muP0qRiRO4vUs4RgRtfhLBy1uzn8bdC38nWrGilbJJT81FoLpsqJO6/oTKl\nwgJ5dcYmzukQ1fzXdDD0+g4ObbKqrGr5bWWjXKukuhutkmqPKasSeeKHtXzYqyFdG128qqrKF3tW\nWqOSRKwrhmhd9Ejln3yrkioiQXlcZEd5ke6NylG3XBjvzN7M2bRMu8PxTuWbWlVWA4paTUkJs+yO\nSBVCLltkR3kXHx/hhVtqs/9ECmOW6OJ4LhNV1ZrLULImTOwLcfr3lypYziSFqsaYF4AzxpivgFuA\n5q4NS7mj5pUj6VinNCN/387Bk7pQvcsULWmt5Fa1Hcx4FBa8rnMZVIHRRXbUFXmmU00yMg3vzUmw\nOxTvFlgUek+ARv1h8Qj4eShk6uJHyvXyusjOOy6NSrmtipFFuOvaSkz5O5ENe0/YHY538/WDzp9A\n62dg7fdWJ3Sqttwq18oxKYiID3Aym0V2viig+JQbGta2KsVDAnhtxkY8bfSaxxGB1k9D5//Cjt/h\ny07WRDelXCTHpGCMOQc8WUCxKA8RFuTPo+2rs3znUebEH7Q7nMKhcX/oOwmStsGY9nB4i90RKS/l\nTPPRPBF5QkTKi0jE+c3lkSm31qdpeaqVLMpbszaRlqFrEBeIau3hrl8h4yyMuwl2L7M7IuWFnEkK\nvYAHgcXAKsems8cKOT9fH567pRb/JCXz9dJddodTeJRrbM1lCI6Ar7vApl/sjkh5GWeSQi1jTEzW\nDajt6sCU+2tdoyQ3VC/BR/O3cvRMmt3hFB4RMVZiKF0PJvWH5aPsjkh5EWeSwl9O7lOF0PO31CI5\nLZMP52kbd4EqEgkDpkONTjBrOMx9Ec5pM566ejktslPasR5zsIg0EpHGjq01EFJgESq3Vq1UKH2b\nVeC75bvZelCHSxaogBDo9Q3EDoY/P4Kp90JGqt1RKQ932dLZWGUt7sIqa5G1TPYp4FkXxqQ8zKPt\nq/Pzmr28OXMTX97dzO5wChcfX7jlfSgWDfNfsYar9v4OgorZHZnyUJe9UjDGfGWMaQPcZYxpk2Xr\nbIz5qQBjVG4uokgAD7etxsKEwyzectjucAofEbj+Mej2BexeCuNuhpP77I5KeahcS2eLyEvAJQcZ\nY151VVA50dLZ7ik1I5ObPlhMoJ8PMx++Hj9fZ7qrVL7bvgAmDbCW/Oz3I5SsZXdEyk3kW+ls4DRw\nxrFlAjcDla4qOuV1Av18eebmmmw5eJqJK/fYHU7hVaUt3D0TzmXC2A6w6w+7I1IeJtekYIx5P8v2\nBtAaqOLyyJTH6VCnNM1jIvhg7hZOpmjxNtuUqQ/3zIXQ0vBNN9jwo90RKQ+Sl2v8EECX3lKXEBFe\nuLU2R5PT+HTBNrvDKdzCK8Cg2VAuFqYMgr/+a3dEykM4s/LaehFZ59jigQTgQ9eHpjxR3XLFuL1x\nNF/+uYvdScl2h1O4hURA/6lQuwv89hzMfkbnMqhcOXOlcCtwm2O7CShrjNE/O9RlDe9QAz9f4e3Z\nm+wORfkHQY/x0PwBWPYZTLkb0nWBJHV5zvQp/AOUB9oaY/YC4SIS4/LIlMcqFRbE/TdUYeb6A6zY\nedTucJSPD9z8Ntz0Bmz8Gb7tDmeP2R2VclPONB+9BDwFPOPYFQB868ybi0hHEUkQkW0i8nQOx/UQ\nESMiuQ6XUp5hyPWVKVMsiFd+iSfznK654BZaDoMe4yBxpTUy6dguuyNSbsiZ5qNuQGesIakYY/YB\nobm9SER8gU+xhrDWBvqIyCWF9EQkFHgYWO582MrdBQf48mynWsTvO8nElbvtDkedV/d2q5/h9EEY\nfSPsWWF3RMrNOJMU0ow1w80AiEgRJ9+7GbDNGLPDGJMGTAS6ZHPca8C7gDZ0eplb65ehReUIRsxJ\n4JhWUXUfla6De+ZBYCiMvxXWT7E7IuVGnEkKk0XkC6y+hCHAPGC0E68rB2SdxZTIRUNZRaQRUN4Y\nMyOnNxKRe0UkTkTiDh/WMgqeQkR4uXMdTqVk8P7cBLvDUVlFVYN75kO5JvDjYFj0DujSqgrnOprf\nA6YAPwJmHQRrAAAZrUlEQVQ1gBeNMZ848d6S3dtdeNJa//kD4HEnYhhljIk1xsSWKFHCiY9W7qJm\n6TD6t6jI98t3s2HvCbvDUVkViYQBP0ODPrDoTZh6n1ZZVc5NXjPGzDXGDDfGPGGMmevkeydijVo6\nLxrIWqUrFKgLLBKRXUALYLp2NnufR9tXp3hIAC9Pjye3WluqgPkFQteR0PZ5WDfJWs3tTJLdUSkb\nOTP6qLuIbBWREyJyUkROichJJ957JVBNRGJEJADoDUw//6Qx5oQxJsoYU8kYUwlYBnQ2xmi1Oy9T\nLNifpzrWJO6fY/y8Zq/d4aiLiUCr4dDjS9i3Gsa0hcO6aFJh5cyVwrtY/1kXM8aEGWNCjTFhub3I\nGJMBDAPmAJuAycaYeBF5VUQ6X13YytP0aBJNg/LhvDlzM6e0LpJ7qtsd7voV0s7AmHawY5HdESkb\nOJMUDhpj8jQ11Rgz0xhT3RhTxVFMD2PMi8aY6dkc21qvEryXj4/wSuc6HD6VyidaF8l9RcdaHdBh\nZeHb22HVeLsjUgXMmaQQJyKTRKSPoympu4h0d3lkyus0LB9Or9jyjPtjJ9sOnbY7HHU5xSvC4N+g\ncmv45T8w80nI1Ku7wsKZpBAGJGPVPTpfA+lWVwalvNfwjjUICfDlxWkbtNPZnQWFQZ9JcM0wWPGF\nVRojWUuWFAY5rdEMgDHm7oIIRBUOUUUDebJjTZ7/eQPT1uyjayOtwu62fP2gwxtQqo51xTC6DfSe\nAKUuKUygvIiumagKXN9mFWhYPpzXf93IiWRtlnB7DfvCXTMh/SyMbQ+bf7U7IuVCmhRUgfPxEd7s\nVo9jyem8PXuz3eEoZ5RvCvcugqjqMLEv/D5CZ0B7KU0Kyha1y4Zxd8tKTFixm1X/aFu1Rwgra63/\nXL8XLHwdfrjLGr6qvIrTSUFEWojIbBFZJCJdXRmUKhwebV+dMsWCeG7qBtIzdUUwj+AfDN2+gPav\nwabpMK4DHNcquN7ksklBREpftOsxrDLanbAqmyp1VYoE+vFy5zpsPnCKcX/stDsc5SwRuPZh6DsZ\nju2GUa11opsXyelK4XMReUFEghyPjwM9sBKDM2UulMpVhzqlaVerFB/O20riMV3T2aNUaw9DFkCR\nkvBNN1jyf9rP4AUumxSMMV2BNcAMEekPPAIEApGANh+pfPNKlzqIwHNTde6Cx4mqaq3NUKcbzH8F\nJvWDFK2G68ly7FMwxvwCdADCgZ+ABGPMx8YYXdRA5Zty4cE82aEGv285zE9/a8E8jxNYFG4fCx3f\nhi2zYVQbOLjR7qhUHuXUp9BZRBYCs4ENWFVOu4nIBBGpUlABqsJhwDWViK1YnFdnbOTQKV2Ez+OI\nQIsHYOAvkHYaxtyoK7p5qJyuFF7HWl+5J/COMea4MeYx4EXgjYIIThUePj7COz3qczY9kxd+1mYk\nj1WxJdy3GMo0sFZ0m/W01k3yMDklhRNAd+B24ND5ncaYrcaY3q4OTBU+VUoU5dF21ZkTf5CZ6w/Y\nHY7Kq9DS1hVDi6GwfCR82QmO78n9dcot5JQUumF1KvsBfQsmHFXYDbk+hnrlivHS9A0cPZNmdzgq\nr3z9oeNb1sI9hzbB59fB5pl2R6WckNPooyPGmE+MMZ8bY3QIqioQfr4+jLijPifOpvPqL/F2h6Ou\nVt3ucN/vVjnuiX1g9jOQocnenWmZC+V2apYOY2jrqvy8Zh9zNx60Oxx1tSKrwOC50Ow+WPYZjLsJ\njupkRXelSUG5pQfbVKVWmTCe+WkdR06n2h2Oulp+gdDpXej1LRzdAV+0gvif7Y5KZUOTgnJLAX4+\nfNirISfPZvDMT+t1NJK3qHUb3LcEoqrBDwNhxmOQrkOQ3YkmBeW2apQO5cmONZi78SA/rEq0OxyV\nX4pXhLtnQ8uHIG4sjG4LB7X/yF1oUlBubdC1MbSoHMEr0+PZc1RrI3kNvwC46XW4cwqcOWwV1Vv6\nKZzTarl206Sg3JqPj/DeHQ3wEeGxyWvIPKfNSF6lWnsYuhSqtoM5z8K33eDkPrujKtQ0KSi3F108\nhFe61GHlrmN8tnCb3eGo/FYkCnp/D7d9BHtWwMiWsHGa3VEVWpoUlEfo1qgcnRuU5YN5W1ixU1dq\n8zoi0OQuqxO6eAxMHgA/D4UUnSJV0DQpKI8gIrzRrS7lI0L4z8TVHNPZzt4pqioM/g1aPQlrJ1hX\nDdsX2h1VoaJJQXmM0CB//tunMUdOpzJ8yjodpuqtfP2h7XMw6DfwC4JvusIvj0DqKbsjKxQ0KSiP\nUi+6GM/cXIt5mw7y5Z+77A5HuVL5pnD/Emvo6qrx8Nk1etVQADQpKI9z97WVaFerJG/N2sTaPcft\nDke5kn+wNXR18G/WrOhvusKMR/WqwYU0KSiPIyKM6NGAkqFBPPDtKpK0DIb3K98M7v8DrhkGcV/C\nZy1h2zy7o/JKmhSURypeJIAv+jch6UwaD01YTUamTnryev7B0OENGDTHumr49nb48R44rasD5ydN\nCspj1S1XjNe71uWv7UmMmJNgdziqoFRoDg/8CTc8bc1n+G8s/P016MCDfKFJQXm0O2LL069FBb5Y\nvINf1+23OxxVUPwCoc0zcP+fUKoOTH8Ixt8Ch7fYHZnH06SgPN6Lt9ahUYVwhk9Zy+YDOtmpUClR\nHQbOgM6fwMEN8Pm1sPAtrbx6FTQpKI8X4OfD5/2aEBrkx+DxcRw+pR3PhYqPDzQeAMPioFZn+P1t\n+Ky5tfynNildMZcmBRHpKCIJIrJNRJ7O5vnHRGSjiKwTkfkiUtGV8SjvVSosiDEDmnL0TBpDvo4j\nJT3T7pBUQStaEnqMhQHTwDfQWv7zux5wROtlXQmXJQUR8QU+BW4GagN9RKT2RYetBmKNMfWBKcC7\nropHeb960cX4oFdD1iYe54kf1uqM58KqcmurI7rDm1aBvc9awNwXdW6Dk1x5pdAM2GaM2WGMSQMm\nAl2yHmCMWWiMOV8kfxkQ7cJ4VCHQsW5pnupYkxnr9vPBvK12h6Ps4usP1zxoNSnV7wl/fgT/bQrr\nJmuTUi5cmRTKAXuyPE507LucwcCs7J4QkXtFJE5E4g4f1jHJKmf3tapMz9hoPp6/lYkrdtsdjrJT\naCno+hkMngehpeGnIdZKb7v+tDsyt+XKpCDZ7Ms2RYtIPyAWGJHd88aYUcaYWGNMbIkSJfIxROWN\nRITXu9bjhuoleHbqeubEH7A7JGW38k3hngXQdSScOgDjO8GEPjqENRuuTAqJQPksj6OBS5ZUEpF2\nwHNAZ2OMDhtR+SLAz4eR/RpTPzqchyasZun2JLtDUnbz8YGGfeGhVdD2Bdi5xOpvmPEYnD5kd3Ru\nw5VJYSVQTURiRCQA6A1Mz3qAiDQCvsBKCPqtqHwVEuDHl3c1pUJECEO+jmPD3hN2h6TcQUAItHoC\nHl4NsYOsCqwfN4Lf39XOaFyYFIwxGcAwYA6wCZhsjIkXkVdFpLPjsBFAUeAHEVkjItMv83ZK5Unx\nIgF8M7gZYUF+3PXlCrYd0l965VC0BNzyHjy43BqxtPAN+KgB/PkxpCXn9mqvJZ42bC82NtbExcXZ\nHYbyMNsPn6bXF8sAmHhvC6qWLGpzRMrt7FlpJYYdC6FoKbj+cWuJUL9AuyPLFyKyyhgTm9txOqNZ\nFQpVShRl4r3NAUOf0cvYfvi03SEpd1O+KQz4Ge6eBZFVYdaTVrNS3JeQmW53dAVGk4IqNKqWDGXC\nkBYYY+gzShODuoyKLeGuX6H/zxBWFmY8Ap80hpVjCkVNJU0KqlCpViqU74e0IPOclRgSDmgfg8qG\nCFRpA4PnQt8frOakXx+Hj+pbE+G8uENak4IqdKqXCmXCvS0A6PnFUlb9c8zmiJTbEoHqN1nJYeAM\nKFnbKpnxQR1Y8Aac8b6hzpoUVKFUvVQoPz7QkvAQf/qNWc7vW3SmvMqBCMRcb/U5DFkAla6Hxe/C\nh3Vh1tNwdKfdEeYbTQqq0CofEcKU+1tSKaoI93y1kulrL5lbqdSlyjWB3t/B0OVQuwusHG11SE+8\n0yqf4WEjOi+mQ1JVoXcyJZ17vopjxc6jDO9Qg6GtqyCSXZUWpbJxcr+VGOLGwdljULq+VYyvTnfw\nC7A7ugucHZKqSUEpICU9k6d+XMe0Nfvo3rgcb3WvR6Cfr91hKU+SlgzrJsGykXAkweqcjh1kLQAU\nVtbu6DQpKHWljDH8d8E23p+7hdiKxfmifxMii3rHxCVVgIyB7fOt5LBtHogv1LgZmtwNVdpaNZhs\noElBqTz6dd1+Hpu8hqiigXx2Z2MalA+3OyTlqY7ugFVfwepvIfkIhFeEJgOhUX9rpbgCpElBqauw\nLvE4D3z7N4dPpfL8rbXo36Ki9jOovMtIhc0zrNnRu5aAj5919dCgL1Rrby0K5GKaFJS6SsfOpPHY\n5DUsTDjMbQ3K8lb3ehQN9LM7LOXpjmy1KrOumwRnDkNIFNS7wyrrXaa+yz5Wk4JS+eDcOcPI37fz\n/m8JVIwswge9GtJQm5NUfshMt/oc1nwPW2ZDZhqUqgsN+kC9HtZKcflIk4JS+Wjp9iQen7yGg6dS\nebBNVR5qWxV/X53mo/JJ8lHY8COsnQB7VwECla6Dut2hVhcoEnnVH6FJQal8duJsOq9Mj+en1Xup\nH12M/+vZgKolQ+0OS3mbw1usBLHhR0jaao1eqtwa6t4ONW+B4LxdqWpSUMpFZq7fz7NT15OcmsnQ\nNlV4oHUVndOg8p8xcHCDI0H8BMf/gQ5vWhPj8kCTglIudOhUCq/P2MT0tfuoHFWE17vVpWWVKLvD\nUt7KGNj7NxSvCEXy9nOmi+wo5UIlQ4P4uE8jvh7UjIxzhr6jl/PIxNXsO37W7tCUNxKB6CZ5TghX\nQpOCUlehVfUS/PZoK4a1qcrMDQdo894i3v8tgdOpGXaHplSeaFJQ6ioF+fvyRIcaLHj8BjrUKc0n\nC7bResQivl32D2kZ5+wOT6kroklBqXwSXTyEj/s0YurQllSKDOH5nzfQ5r1FTFixm/RMTQ7KM2hH\ns1IuYIxh8dYjfDB3C2v2HCe6eDD331CFHk2iCfLXkUqq4OnoI6XcgDGGRVsO8+G8razdc5yIIgH0\nb1GR/tdUJEorsKoCpElBKTdijGHFzqOMXrKDeZsOEeDnQ+cGZenTrAKNK4RrsT3lcs4mBa3upVQB\nEBGaV46keeVIth06zbg/dzJt9V6mrEqkeqmi9GlWge6NoikW4vpqmUrlRK8UlLLJ6dQMflm7j4kr\ndrM28QSBfj60q12K2+qXpXWNEtr3oPKVNh8p5UHi951g0so9/LpuP0ln0ggN9OOmOqW5rUEZWlaJ\nIsBPBwqqq6NJQSkPlJF5jqU7kvhl7T5mbTjAqZQMigb60ap6FG1rlqJ1jRLaQa3yRJOCUh4uNSOT\nP7YeYd6mQyzYfJCDJ1MRgQbR4VxXNYprqkTSpGJxbWZSTtGkoJQXMcYQv+8k8zcdYmHCIdbvPUHm\nOUOArw8Ny4fTokokLSpH0CA6nCK6OpzKhiYFpbzYqZR04nYdY+mOJJZuTyJ+3wnOGfARqF4qlAbR\n4TQoH06D8sWoXipUFwRSmhSUKkxOnE3n73+OsWbPcdYmHmftnuMcS04HIMDPh2oli1KjVCg1SodS\nvXQoNUuHUjosSOdHFCKaFJQqxIwx7Dl6ljWJx1mfeJyEg6fZcuAUB06mXDgmNMiP6qVCqRgZQqXI\nIlSMDKFiZBEqRYYQHhJgY/TKFTQpKKUucTw5jS0HT5Nw8BQJB06y9eBpdh9NZv+JlH8dFxbkR4XI\nEMoUC6ZMsSBKFwuibLFgShcLokyxIEqFBWkHt4fRGc1KqUuEhwTQLCaCZjER/9qfkp7JnqPJ7EpK\n5p+kM/yTlMzuo8nsOZrM8h1JnEy5dH2I8BB/IosEEFk00HEbQGSRwH/dhof4ExbkT1iwP0UCfLW5\nygO4NCmISEfgI8AXGGOMefui5wOBr4EmQBLQyxizy5UxKaUuFeTvS7VSoVQrFZrt82dSMzhwMoUD\nJ1LYfyKF/cfPcuhUKkfPpHHkdCpbD51m2Y5Ujp9N53KNDz4CYcFWkggN8nMkCz/HY3+KBvoSFOBL\niL8vIQF+BAf4EhLgS7C/r+O+HyEBvgT5/2+/j48mmfzmsqQgIr7Ap0B7IBFYKSLTjTEbsxw2GDhm\njKkqIr2Bd4BeropJKZU3RQL9qFKiKFVKFM3xuIzMcxxLTifpTCpJp9M4cTadk2fTOZmSzsmzGY7b\ndE6mZHAqJZ1dR5Iv7DuTlnnFcfn7CgG+PgT4Zdl8ffD39SHwon3WfV/8fcV6ztcHXx8f/HwFHxH8\nfARfH8etr3V7Yb+vz4XnfUXw881yrI/1nI/jOR+xal35CPj4WLfgeCzWe4rjvlzYl+U1WY6RLI99\nBEKD/AkOcG2znSuvFJoB24wxOwBEZCLQBciaFLoALzvuTwH+KyJiPK2jQykFgJ+vDyVCAykReuWz\nro0xpKSfIzktg+S0TFLSM0lOs7az6RmcTbOeO5tlf3rmOdIysmyOx6kZ5y48l5J+jpNnM/71/Pnb\nc+cMGecMmecMmca6dWevda1L/xYVXfoZrkwK5YA9WR4nAs0vd4wxJkNETgCRwJGsB4nIvcC9joen\nRSThoveJuvg1HsobzkPPwT14wzmAd5xHvp3DgHdgQN5f7lQ2cWVSyK6x7+I07MwxGGNGAaMu+0Ei\ncc70qrs7bzgPPQf34A3nAN5xHp52Dq6c5pgIlM/yOBrYd7ljRMQPKAYcdWFMSimlcuDKpLASqCYi\nMSISAPQGpl90zHRgoON+D2CB9icopZR9XNZ85OgjGAbMwRqSOs4YEy8irwJxxpjpwFjgGxHZhnWF\n0DuPH3fZpiUP4w3noefgHrzhHMA7zsOjzsHjZjQrpZRyHS2dqJRS6gJNCkoppS7w+KQgIh1FJEFE\ntonI03bH4ywR2SUi60VkjYjEOfZFiMhcEdnquC1ud5wXE5FxInJIRDZk2Zdt3GL52PHdrBORxvZF\n/j+XOYeXRWSv4/tYIyKdsjz3jOMcEkSkgz1R/5uIlBeRhSKySUTiReQ/jv0e813kcA4e812ISJCI\nrBCRtY5zeMWxP0ZElju+h0mOwTaISKDj8TbH85XsjD9bxhiP3bA6sLcDlYEAYC1Q2+64nIx9FxB1\n0b53gacd958G3rE7zmzibgU0BjbkFjfQCZiFNR+lBbDc7vhzOIeXgSeyOba24+cqEIhx/Lz5usE5\nlAEaO+6HAlscsXrMd5HDOXjMd+H49yzquO8PLHf8+04Gejv2fw484Lg/FPjccb83MMnu7+HizdOv\nFC6U0jDGpAHnS2l4qi7AV477XwFdbYwlW8aYxVw6l+RycXcBvjaWZUC4iJQpmEgv7zLncDldgInG\nmFRjzE5gG9bPna2MMfuNMX877p8CNmFVCPCY7yKHc7gct/suHP+epx0P/R2bAdpile6BS7+H89/P\nFOBGcbPSsZ6eFLIrpZHTD5U7McBvIrLKUcYDoJQxZj9YvzBASduiuzKXi9vTvp9hjqaVcVma7tz+\nHBxNEI2w/kr1yO/ionMAD/ouRMRXRNYAh4C5WFcwx40x5+uNZ43zX6V9gPOlfdyGpycFp8pkuKlr\njTGNgZuBB0Wkld0BuYAnfT8jgSpAQ2A/8L5jv1ufg4gUBX4EHjHGnMzp0Gz2ucV5ZHMOHvVdGGMy\njTENsao2NANqZXeY49YtzyErT08KzpTScEvGmH2O20PAVKwfpoPnL+kdt4fsi/CKXC5uj/l+jDEH\nHb/c54DR/K9Zwm3PQUT8sf4z/c4Y85Njt0d9F9mdgyd+FwDGmOPAIqw+hXBH6R74d5xuX9rH05OC\nM6U03I6IFBGR0PP3gZuADfy77MdAYJo9EV6xy8U9HRjgGPnSAjhxvmnD3VzUvt4N6/sA6xx6O0aN\nxADVgBUFHd/FHO3QY4FNxpj/y/KUx3wXlzsHT/ouRKSEiIQ77gcD7bD6RhZile6BS78H9y7tY3dP\n99VuWKMqtmC14z1ndzxOxlwZaxTFWiD+fNxYbYvzga2O2wi7Y80m9glYl/TpWH/1DL5c3FiXyp86\nvpv1QKzd8edwDt84YlyH9YtbJsvxzznOIQG42e74HTFdh9XssA5Y49g6edJ3kcM5eMx3AdQHVjti\n3QC86NhfGSthbQN+AAId+4Mcj7c5nq9s9zlcvGmZC6WUUhd4evORUkqpfKRJQSml1AWaFJRSSl2g\nSUEppdQFmhSUUkpdoElBqYs4qnQ+ISKvikg7x75FIpLj4usi0lVEamd5fOH1SnkKly3HqZSnM8a8\neIUv6QrMADbm8fVK2U6vFJQCROQ5EdkiIn8ANRz7xotIj2yOPZ3lfg/HcS2BzsAIxxoAVbK+XkRu\nFJHVYq2hMU5EAh37d4nIKyLyt+O5mgVywkpdhiYFVeiJSBOsEikNsWbUNr3S9zDG/IU1+3a4Maah\nMWZ7lvcPAsYDvYwx9bCu0B/I8vIjxiqOOBJ4Iq/noVR+0KSgFFwPTDXGJBurSmd+18+qAew0xmxx\nPP4Ka6Gf884Xs1sFVMrnz1bqimhSUMpyJfVesh4b5MTxuS2ikuq4zUT7+ZTNNCkoBYuBbiIS7Khe\ne1suxx8UkVoi4oNVxfO8U1jLSl5sM1BJRKo6HvcHfr/aoJVyBU0KqtAz1pKQk7Cq1s7CKsmek6ex\nRhn9hVVt9byJwHBHh3KVLO+fAtwN/CAi64FzWOv2KuV2tEqqUkqpC/RKQSml1AWaFJRSSl2gSUEp\npdQFmhSUUkpdoElBKaXUBZoUlFJKXaBJQSml1AX/D2QC4xtuXlAsAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "These parameters can then be used to simulate data:"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "dilutions = np.array([10, 20, 40, 80, 160, 320, 640])\n\nn = 30\nx1 = np.random.binomial(n, invlogit(β_1_true[0] + β_1_true[1]*dilutions))\nx2 = np.random.binomial(n, invlogit(β_2_true[0] + β_2_true[1]*dilutions))",
      "execution_count": 203,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Here is one simulated dataset:"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "plt.plot(dilutions, x1/n, 'bo', label='group 1')\nplt.plot(dilutions, x2/n, 'ro', label='group 2')\nplt.ylim(0, 1)\nplt.xlabel('dilution')\nplt.ylabel('% neutralization')\nplt.legend()",
      "execution_count": 204,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": "<matplotlib.legend.Legend at 0x12165ecc0>"
          },
          "metadata": {},
          "execution_count": 204
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<matplotlib.figure.Figure at 0x121074240>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGyhJREFUeJzt3XuUVOWZ7/Hvj4tBEDVe5sRlC41KEEeRS8toJAxHkglq\nBE1co9jjLWovRzkZl8Y5ephjIhnWyiQeczRxTDqcREc7olEnoktD4gVzODHRJvEKGlFp7DFRZLwh\nIYA854+9e1s01d3Vl93VVf37rFWrar/11q7nbTb99N773c9WRGBmZgYwpNwBmJnZwOGkYGZmGScF\nMzPLOCmYmVnGScHMzDJOCmZmlsktKUj6oaQ3JT3XwfuSdIOktZKekTQ1r1jMzKw0ee4p3AzM6eT9\nE4Dx6aMBuCnHWMzMrAS5JYWI+CXwn510mQf8WyR+Dewt6YC84jEzs64NK+N3Hwi8VrDcmrb9oX1H\nSQ0kexOMGjVq2mGHHda7b161quP3pk3r3brNzAagVatWvRUR+3fVr5xJQUXaitbciIhGoBGgrq4u\nmpube/fNtbXQ0rJr+9ix0Nt1m5kNQJKK/NLbVTlnH7UCBxUs1wCv98s3L14MI0fu3DZyZNJuZjaI\nlTMpLAPOTmchHQO8GxG7HDrKRX09NDYmewZS8tzYmLSbmQ1iuR0+knQ7MAvYT1Ir8FVgOEBEfA94\nADgRWAtsBs7LK5ai6uudBMzM2sktKUTE/C7eD+CSvL6/UFMTLFwI69fDmDHJUSLnA7PKsm3bNlpb\nW9myZUu5QxnQRowYQU1NDcOHD+/R58t5orlfNDVBQwNs3pwst7Qky+DEYFZJWltbGT16NLW1tUjF\n5qlYRLBx40ZaW1sZN25cj9ZR9WUuFi78KCG02bw5aTezyrFlyxb23XdfJ4ROSGLfffft1d5U1SeF\n9eu7125mA5cTQtd6+zOq+qQwZkz32s3MBrOqTwq+JMHMqsV3v/tdDj30UCTx1ltv5fIdVZ8UfEmC\n2eDU1JQULxgyJHluauqf792+fXtu6z7uuON46KGHGDt2bG7fUfVJAaCeJtZRyw6GsI5a6unG1lGu\nLcvMeqxt1mFLC0R8NOuwt/99v/71rzNhwgRmzJjB/PnzufbaawGYNWsWl156KXV1dVx//fW0tLQw\ne/ZsJk2axOzZs1mfnsQ899xzueuuu7L17bHHHgCsWLGCmTNnctJJJzFhwgQuuugiduzYscv3T5ky\nhdra2t4NogvVnxR6s3XktWWZWa7ymHXY3NzM3XffzVNPPcWDDz5I+xpsW7dupbm5mcsvv5wFCxZw\n9tln88wzz1BfX8+Xv/zlLtf/xBNP8J3vfIfVq1fz8ssvc8899/Q82F6o/qTQm63D81nNKlIesw5X\nrlzJvHnz2H333Rk9ejQnn3zyTu+ffvrp2evHH3+cM888E4CzzjqLlStXdrn+6dOnc/DBBzN06FDm\nz59f0mfyUP1JoTdbh+ezmlWkPGYdJkUYOjZq1KgO32ubJjps2LDssFBEsHXr1l36dLTcX6o/KfRm\n6/B8VrOKlMeswxkzZnDfffexZcsWNm3axP33399h30996lMsXboUgKamJmbMmAFAbW0tq9L7udx7\n771s27Yt+8wTTzzBq6++yo4dO7jjjjuyz/S36k8Kvdk6PJ/VrCLlMevw6KOPZu7cuUyaNIkTTjiB\nI488kr322qto3xtuuIEf/ehHTJo0iVtvvZXrr78egAsvvJDHHnuMo446iscff3ynvYujjz6aBQsW\nMHHiRMaNG8epp55adL01NTW0trYyadIkLrjggp4PqCMRUVGPadOmRbfddlvE2LERUvJ8223981kz\n6zOrV68udwjx/vvvR0TEBx98ENOmTYtVq1b1yXofffTROOmkk/pkXRHFf1ZAc5TwO7bqC+IBvSqT\n3UQ9C6lnPTAGWAz4EgezwamhoYHVq1ezZcsWzjnnHKZOnVrukPrc4EgKPeQKq2ZW6Mc//nEu6501\naxazZs3KZd3dVf3nFHrBM1LNbLBxUuiEZ6Sa2WDjpNAJz0g1s8HGSaETnpFqZoONk0InXGHVzAaS\n+vp6JkyYwBFHHMGXvvSlnS5+6ytOCl2or4d162DHjuTZCcGsQpSpwnGepbPr6+t54YUXePbZZ/nT\nn/7EkiVL+vw7nBTMrPrkVOG43KWzTzzxRCQhienTp9Pa2tqr8RTjpGBm1SeH+eQDqXT2tm3buPXW\nW5kzZ06Px9MRJwUzqz45zCcfSKWzL774YmbOnMmnP/3pHo6mY04KZlZ9cphPHgOkdPY111zDhg0b\nuO6660qKu7ucFMys+uQwn3wglM5esmQJy5cv5/bbb2fIkHx+fTspmFn1yWE++UAonX3RRRfxxhtv\ncOyxxzJ58mQWLVrU4/F0RF3tEg00dXV10f4Ej5lVvzVr1jBx4sSyxrBp0yb22GMPNm/ezMyZM2ls\nbOyTSqkrVqzg2muv7XTvozuK/awkrYqIuq4+O2j3FMoyhblM86bNrG80NDQwefJkpk6dyhe/+EWX\nzq4WZSmJ7TrcZhXPpbOrVFlKYrsOt1mvVdrh7nLo7c9oUCaFspTEdh1us14ZMWIEGzdudGLoRESw\nceNGRowY0eN1DMrDR2PGJEdvirVX15eaVY+2G9Zv2LCh3KEMaCNGjKCmpqbHnx+USWHx4p0P70M/\nlMQuy5eaVY/hw4czbty4codR9Qbl4aOylMR2HW4zqwC5XqcgaQ5wPTAUWBIR32j3/hjgFmDvtM+V\nEfFAZ+v0dQpmZt1X9usUJA0FbgROAA4H5ks6vF23fwLujIgpwBnAv+YVj5mZdS3Pw0fTgbUR8UpE\nbAWWAvPa9Qlgz/T1XsDrOcZjZmZdyDMpHAi8VrDcmrYV+hrwd5JagQeA/1ZsRZIaJDVLavbMAzOz\n/OSZFIrVfW1/AmM+cHNE1AAnArdK2iWmiGiMiLqIqNt///1zCNXMzCDfpNAKHFSwXMOuh4fOB+4E\niIjHgRHAfjnGZGZmncgzKTwJjJc0TtJuJCeSl7Xrsx6YDSBpIklS8PEhM7MyyS0pRMR2YAGwHFhD\nMsvoeUmLJM1Nu10OXCjpaeB24NzwNexmZmWT6xXN6TUHD7Rru7rg9WrguDxjMDOz0g3KK5rNzKw4\nJwUzM8s4KZiZWcZJwczMMk4KZmaWcVIwM7OMk4KZmWWcFMzMLOOkYGZmGScFMzPLOCmYmVnGScHM\nzDJOCmZmlnFSMDOzjJOCmZllnBTMzCzjpGBmZhknBTMzyzgpmJlZxknBzMwyTgpmZpYZ1lUHSZ8E\nrgDGFvaPiONzjMvMzMqgy6QA/AT4HvAD4MN8wzEzs3IqJSlsj4ibco/EzMzKrpRzCvdJuljSAZL2\naXvkHpmZmfW7UvYUzkmfryhoC+Dgvg/HzMzKqcukEBHj+iMQMzMrv1JmHw0H/h6YmTatAL4fEdty\njMvMzMqglMNHNwHDgX9Nl89K2y7IKygzMyuPUpLC0RFxVMHyI5KezisgMzMrn1JmH30o6ZC2BUkH\n4+sVzMyqUil7ClcAj0p6BRDJlc3n5RqVmZmVRSmzjx6WNB6YQJIUXoiIP+cemZmZ9bsOk4Kk4yPi\nEUlfaPfWIZKIiHtyjs3MzPpZZ3sKfw08Apxc5L0AukwKkuYA1wNDgSUR8Y0iff4W+Fq6zqcj4syu\nwzYzszx0mBQi4qvpy0UR8Wrhe5K6vKBN0lDgRuCzQCvwpKRlEbG6oM944CrguIh4W9Jf9GAMZmbW\nR0qZfXR3kba7SvjcdGBtRLwSEVuBpcC8dn0uBG6MiLcBIuLNEtZrZmY56eycwmHAXwJ7tTuvsCcw\nooR1Hwi8VrDcCvxVuz6fTL/r/5EcYvpaRPysSCwNQAPAmDFjSvhqMzPric7OKUwAPg/szc7nFd4n\n+Qu/KyrSFkW+fzwwC6gB/q+kIyLinZ0+FNEINALU1dW1X4eZmfWRzs4p3AvcK+nYiHi8B+tuBQ4q\nWK4BXi/S59dpHaVXJb1IkiSe7MH3mZlZL5Vy8drvJF1CcigpO2wUEV/q4nNPAuPTk9L/AZwBtJ9Z\n9FNgPnCzpP1IDie9UmLsZmbWx0o50Xwr8Angc8BjJH/xv9/VhyJiO7AAWA6sAe6MiOclLZI0N+22\nHNgoaTXwKHBFRGzs/jDMzKwvlJIUDo2I/wl8EBG3ACex6wnjoiLigYj4ZEQcEhGL07arI2JZ+joi\n4rKIODwijoyIpT0diPVSUxPU1sKQIclzU1O5IzKzMijl8FHbfRPekXQE8EfA1xNUk6YmaGiAzZuT\n5ZaWZBmgvr58cZlZvytlT6FR0seBfwKWAauBb+YalfWvhQs/SghtNm9O2s1sUCmlIN6S9OUv8X2Z\nq9P69d1rN7Oq1eWegqQPJX1DkgrafptvWNavOrog0BcKmg06pRw+ej7t93NJ+6RtxS5Ms0q1eDGM\nHLlz28iRSbuZDSqlJIXtEfGPwA9Irjiexq5XJlslq6+HxkYYOxak5Lmx0SeZzQahUmYfCSAi7pT0\nPHA74OMK1aa+3knAzEpKChe0vUgvPpsBnJJfSGZmVi5d3nkNGCtpbLu3N+UblpmZlUOud14zM7PK\n0uWd1yLivP4Lx8zMyqmzw0eXdfbBiLiu78MxM7Ny6uzw0eh+i8LMzAaEzg4fXdOfgQwGTU1JOaH1\n65OLhRcv9ixQMxtYupySKmkEcD7dv8mOFXAhUjOrBLndZMd25kKkZlYJcr3Jjn3EhUjNrBKUkhTa\n32RnL3yTnW5zIVIzqwQ9vcnOv+QaVRVyIVIzqwSdnmiWNAR4LyLexjfZ6ZW2k8mefWRmA5kiOq+C\nLak5Iur6KZ4u1dXVRXNzc7nDMDOrKJJWlfK7vJTDRw9J+oqkgyTt0/bogxjNzGyAKaV09unp8yUF\nbYEPJZmZVZ1SksLEiNhS2JBe0GZmZlWmlMNHvyqxzczMKlxnVVI/ARwI7C5pCultOYE9gZEdfc7M\nzCpXZ4ePPgecS1LWorBM9vvA/8gxJjMzK5POqqTeAtwi6YsRcXc/xmRmZmVSyonmIyT9ZfvGiFiU\nQzxmZlZGpSSFTQWvRwCfB9bkE46ZmZVTl0khIv5X4bKka4Gf5xaRmZmVTSlTUtsbSTIryczMqkwp\nd157luQKZoChwP6AzyeYmVWhUs4pfL7g9XbgjYjYnlM8ZmZWRl0ePoqIFuAg4PiI+A9gb0njco/M\nzMz6XZdJQdJXgf8OXJU27QbcVsrKJc2R9KKktZKu7KTfaZJC0oAp0W1mNhiVcqL5VGAu8AFARLwO\njO7qQ5KGAjcCJwCHA/MlHV6k32jgy8BvSg/bzMzyUEpS2BrJnXgCQNKoEtc9HVgbEa9ExFZgKTCv\nSL+vA98EthR5z8zM+lEpSeFOSd8nOZdwIfAQ8IMSPncg8FrBcivtprKmhfYOioj7O1uRpAZJzZKa\nN2zYUMJXm5lZT5Ry8dq1kj4LvAdMAK6OiF+UsG4Vacvu/Zne//nbJEX3uoqhEWiE5HacJXy3mZn1\nQClTUkmTQCmJoFAryaylNjXA6wXLo4EjgBWSAD4BLJM0NyJ8E2YzszIoZfbRFyS9JOldSe9Jel/S\neyWs+0lgvKRxknYDzgCWtb0ZEe9GxH4RURsRtcCvAScEM7MyKuWcwjdJflnvFRF7RsToiNizqw+l\nF7gtAJaTFNC7MyKel7RI0tzehW1mZnko5fDRGxHRo6qoEfEA8EC7tqs76DurJ99hZmZ9p5Sk0Czp\nDuCnwJ/bGiPintyiMjOzsiglKewJbAb+pqAtACcFM7MqU8qU1PP6IxAzMyu/ntxPwczMqpSTgpmZ\nZZwUzMwsU3JSkHSMpJ9JWiHplDyDMjOz8ujwRLOkT0TEHwuaLiMpoy2SMtc/zTk2MzPrZ53NPvqe\npFXAtyJiC/AOcBqwg6Q4npmZVZkODx9FxCnAU8D9ks4CLgU+BuwL+PCRmVkV6vScQkTcB3wO2Jvk\nYrUXI+KGiPBNDczMqlCHSUHSXEmPAj8DniOpcnqqpNslHdJfAZqZWf/p7JzCP5PcUnN3YHlETAcu\nkzQeWEySJMzMrIp0lhTeBb4AjATebGuMiJdwQjAzq0qdnVM4leSk8jDgzP4Jx8zMyqnDPYWIeAv4\nTj/GYmZmZeYyF2ZmlnFSMDOzjJOCmZllnBTMzCzjpGBmZhknBTMzyzgpmJlZxknBzMwyTgpmZpZx\nUjAzs4yTgpmZZZwUzMws46RgZmYZJwUzM8s4KZiZWcZJwczMMk4KZmaWcVIwM7NMrklB0hxJL0pa\nK+nKIu9fJmm1pGckPSxpbJ7xmJlZ53JLCpKGAjcCJwCHA/MlHd6u2++AuoiYBNwFfDOveMzMrGt5\n7ilMB9ZGxCsRsRVYCswr7BARj0bE5nTx10BNjvGYmVkX8kwKBwKvFSy3pm0dOR94sNgbkhokNUtq\n3rBhQx+GaGZmhfJMCirSFkU7Sn8H1AHfKvZ+RDRGRF1E1O2///59GKKZmRUaluO6W4GDCpZrgNfb\nd5L0GWAh8NcR8ecc4zEzsy7kuafwJDBe0jhJuwFnAMsKO0iaAnwfmBsRb+YYi5mZlSC3pBAR24EF\nwHJgDXBnRDwvaZGkuWm3bwF7AD+R9JSkZR2szszM+kGeh4+IiAeAB9q1XV3w+jN5fr+ZmXWPr2g2\nM7OMk4KZmWWcFMzMLOOkYGZmGScFMzPLOCmYmVnGScHMzDJOCmZmlnFSMDOzjJOCmZllnBTMzCzj\npGBmZhknBTMzyzgpGABNTVBbC0OGJM9NTeWOyMzKIdfS2VYZmpqgoQE2b06WW1qSZYD6+vLFZWb9\nz3sKxsKFHyWENps3J+1mNrg4KRjr13ev3cyql5OCMWZM99rNrHo5KRiLF8PIkTu3jRyZtJvZ4OKk\nYNTXQ2MjjB0LUvLc2OiTzGaDkWcfGZAkACcBM/OegpmZZZwUzMws46RgZmYZJwUzM8s4KZiZWcZJ\nwczMMk4KZmaWcVIw6ynXG7cq5IvXzHrC9catSnlPwawnXG/cqpSTgllPuN64VSknBbOecL1xq1JO\nCmY94XrjVqWcFMx6wvXGrUrlmhQkzZH0oqS1kq4s8v7HJN2Rvv8bSbV5xmPWl5qop5Z1DGEHtayj\nCScEy8fKi5toHVbLDg2hdVgtKy/Ob/pzbklB0lDgRuAE4HBgvqTD23U7H3g7Ig4Fvg38S17xmPWl\nthmpLS0Q8dGMVF+qYH1t5cVNTLmpgZoPWxhCUPNhC1NuasgtMeS5pzAdWBsRr0TEVmApMK9dn3nA\nLenru4DZkpRjTGZ9wjNSrb/UNi5kFDtvbKPYTG1jPhtbnhevHQi8VrDcCvxVR30iYrukd4F9gbcK\nO0lqANIrg9gk6cUSvn+/9uupQJU+hiqOf9q0Yq0tLSCtWpVjTN1Vxf8GFaHX8U+DotsaH7awSurO\ntja2lE55JoVif/FHD/oQEY1AY7e+XGqOiLrufGagqfQxOP7yq/QxOP7+l+fho1bgoILlGuD1jvpI\nGgbsBfxnjjGZmVkn8kwKTwLjJY2TtBtwBrCsXZ9lwDnp69OARyJilz0FMzPrH7kdPkrPESwAlgND\ngR9GxPOSFgHNEbEM+D/ArZLWkuwhnNGHIXTrcNMAVeljcPzlV+ljcPz9TP7D3MzM2viKZjMzyzgp\nmJlZpiqTQlflNQYCST+U9Kak5wra9pH0C0kvpc8fT9sl6YZ0PM9Imlq+yLNYD5L0qKQ1kp6X9A9p\neyWNYYSkJyQ9nY7hmrR9XFp25aW0DMtuafuALMsiaaik30m6P12umPglrZP0rKSnJDWnbRWzDQFI\n2lvSXZJeSP8/HFtpYyhUdUmhxPIaA8HNwJx2bVcCD0fEeODhdBmSsYxPHw3ATf0UY2e2A5dHxETg\nGOCS9OdcSWP4M3B8RBwFTAbmSDqGpNzKt9MxvE1SjgUGblmWfwDWFCxXWvz/NSImF8znr6RtCOB6\n4GcRcRhwFMm/RaWN4SMRUVUP4FhgecHyVcBV5Y6rg1hrgecKll8EDkhfHwC8mL7+PjC/WL+B8gDu\nBT5bqWMARgK/Jbnq/i1gWPvtiWQm3bHp62FpP5U57hqSXzrHA/eTXBBaSfGvA/Zr11Yx2xCwJ/Bq\n+59jJY2h/aPq9hQoXl7jwDLF0l3/JSL+AJA+/0XaPqDHlB6GmAL8hgobQ3ro5SngTeAXwMvAOxGx\nPe1SGOdOZVmAtrIs5fS/gX8EdqTL+1JZ8Qfwc0mr0nI2UFnb0MHABuBH6SG8JZJGUVlj2Ek1JoWS\nSmdUmAE7Jkl7AHcDl0bEe511LdJW9jFExIcRMZnkL+7pwMRi3dLnATUGSZ8H3oyIwvo3ncU4oOJP\nHRcRU0kOq1wiaWYnfQdi/MOAqcBNETEF+ICPDhUVMxDHsJNqTAqllNcYqN6QdABA+vxm2j4gxyRp\nOElCaIqIe9LmihpDm4h4B1hBcn5kbyVlV2DnOAdaWZbjgLmS1pFUIT6eZM+hUuInIl5Pn98E/p0k\nMVfSNtQKtEbEb9Llu0iSRCWNYSfVmBRKKa8xUBWW/TiH5Dh9W/vZ6cyFY4B323ZNy0WSSK5IXxMR\n1xW8VUlj2F/S3unr3YHPkJwkfJSk7ArsOoYBU5YlIq6KiJqIqCXZzh+JiHoqJH5JoySNbnsN/A3w\nHBW0DUXEH4HXJE1Im2YDq6mgMeyi3Cc18ngAJwK/Jzk+vLDc8XQQ4+3AH4BtJH89nE9yfPdh4KX0\neZ+0r0hmVL0MPAvUDYD4Z5Ds9j4DPJU+TqywMUwCfpeO4Tng6rT9YOAJYC3wE+BjafuIdHlt+v7B\n5R5DwVhmAfdXUvxpnE+nj+fb/q9W0jaUxjUZaE63o58CH6+0MRQ+XObCzMwy1Xj4yMzMeshJwczM\nMk4KZmaWcVIwM7OMk4KZmWWcFMzakfQ1SV+RtEjSZ9K2FZI6vQG7pFMKiy8Wft6sUuR2O06zShcR\nV3fzI6eQFKVb3cPPm5Wd9xTMAEkLJf1e0kpgQtp2s6TTivTdVPD6tLTfp4C5wLfSewMcUvh5SbPT\ngmnPKrmXxsfS9nWSrpH02/S9w/plwGYdcFKwQU/SNJIyEZNJrso+urvriIhfkZQwuCKSewO8XLD+\nEST3zzg9Io4k2UP/+4KPvxVJUbibgK/0dBxmfcFJwQw+Dfx7RGyOpNJrX9fKmgC8GhG/T5dvAQqr\ngbYVE1xFco8Ns7JxUjBLdKfeS2HfESX0L1YuudCf0+cP8Xk+KzMnBTP4JXCqpN3Tqp0nd9H/DUkT\nJQ0BTi1ofx8YXaT/C0CtpEPT5bOAx3obtFkenBRs0IuI3wJ3kFTrfJCk/HpnriSZZfQrkkq3bZYC\nV6QnlA8pWP8W4DzgJ5KeJblL2vf6bgRmfcdVUs3MLOM9BTMzyzgpmJlZxknBzMwyTgpmZpZxUjAz\ns4yTgpmZZZwUzMws8/8BnDtZWYn+XcEAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Now, lets simulate 100 datasets, and see how often we are able to detect a difference, where \"detect a difference\" means that the probability of group 2 being greater than group 1 is greater than 95%."
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "n_sim = 100\n\nprobs = np.empty(n_sim)\nfor sim in range(n_sim):\n\n    # Simulate data\n    x1 = np.random.binomial(n, invlogit(β_1_true[0] + β_1_true[1]*dilutions))\n    x2 = np.random.binomial(n, invlogit(β_2_true[0] + β_2_true[1]*dilutions))\n\n    # Estimate model\n    with pm.Model() as mod:\n\n        β1 = pm.Normal('β1', 0, sd=10, shape=2)\n        β2 = pm.Normal('β2', 0, sd=10, shape=2)\n\n        # Number of times β1 is greater than β2 (i.e. false positive)\n        p = pm.Deterministic('p',  β1[1]>β2[1])\n\n        pm.Binomial('x1', n=30, p=invlogit(β1[0] + β1[1]*dilutions), observed=x1)\n        pm.Binomial('x2', n=30, p=invlogit(β2[0] + β2[1]*dilutions), observed=x2)\n        \n        # Run model\n        tr = pm.sample(2000, init=None, step=pm.NUTS(), progressbar=None)\n        \n        # Tally count\n        probs[sim] = tr['p', 1000:].mean()",
      "execution_count": 217,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Here is the estimated power, based on these simulations:"
    },
    {
      "metadata": {
        "trusted": true,
        "collapsed": false
      },
      "cell_type": "code",
      "source": "(probs<0.05).mean()",
      "execution_count": 219,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": "0.95999999999999996"
          },
          "metadata": {},
          "execution_count": 219
        }
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.0",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Bayesian Power Analysis.ipynb",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}