{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "import scipy.optimize\n",
    "from ballot_comparison import ballot_comparison_pvalue\n",
    "from hypergeometric import trihypergeometric_optim\n",
    "from fishers_combination import fisher_combined_pvalue, maximize_fisher_combined_pvalue, \\\n",
    "    plot_fisher_pvalues\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N1 = 10000\n",
    "N2 = 1000\n",
    "N_w1 = 4550\n",
    "N_l1 = 4950\n",
    "N_w2 = 750\n",
    "N_l2= 150\n",
    "margin = (N_w1 + N_w2 - N_l1 - N_l2)\n",
    "\n",
    "n1 = 500\n",
    "n2 = 250\n",
    "\n",
    "cvr_pvalue = lambda alloc: ballot_comparison_pvalue(n=n1, gamma=1.03905, o1=0, u1=0, o2=0, u2=0, \n",
    "                                                    reported_margin=margin, N=N1, null_lambda=alloc)\n",
    "nocvr_pvalue = lambda alloc: trihypergeometric_optim(sample= np.array([0]*int(n2*N_l2/N2)+[1]*int(n2*N_w2/N2)+[np.nan]*int(n2*(N2-N_l2-N_w2)/N2)), \n",
    "                                                     popsize=N2, \n",
    "                                                     null_margin=(N_w2-N_l2) - alloc*margin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.23557770396261943, 0.001425699715253263]\n",
      "0.00302236691268\n"
     ]
    }
   ],
   "source": [
    "# This is one possible allocation lambda=0.3\n",
    "pvalues = [cvr_pvalue(0.3), nocvr_pvalue(0.7)]\n",
    "print(pvalues)\n",
    "fisher_pvalue = fisher_combined_pvalue(pvalues)\n",
    "print(fisher_pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-9 11\n",
      "-9.0 0\n",
      "-8.5 0\n",
      "-8.0 0\n",
      "-7.5 0\n",
      "-7.0 0\n",
      "-6.5 0\n",
      "-6.0 0\n",
      "-5.5 0\n",
      "-5.0 0.0\n",
      "-4.5 0.0\n",
      "-4.0 0.0\n",
      "-3.5 0.0\n",
      "-3.0 0.0\n",
      "-2.5 0.0\n",
      "-2.0 0.0\n",
      "-1.5 0.0\n",
      "-1.0 1.07247544179e-13\n",
      "-0.5 2.09982869981e-08\n",
      "0.0 0.000201206808265\n",
      "0.5 0.010103481892\n",
      "1.0 0.0259274436535\n",
      "1.5 0.00575285367914\n",
      "2.0 0.000645115294631\n",
      "2.5 6.78091214908e-05\n",
      "3.0 0\n",
      "3.5 0\n",
      "4.0 0\n",
      "4.5 0\n",
      "5.0 0\n",
      "5.5 0\n",
      "6.0 0\n",
      "6.5 0\n",
      "7.0 0\n",
      "7.5 0\n",
      "8.0 0\n",
      "8.5 0\n",
      "9.0 0\n",
      "9.5 0\n",
      "10.0 0\n",
      "10.5 0\n",
      "11.0 0\n",
      "11.5 0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEwNJREFUeJzt3G2MXNV9x/Hfzw8bbxJsFby4rQFvAqGhKClFqusmKBpK\ngSWqBE2rYFCjQp3UUuq2al/U0DbyEqUq9EXapIgSky3Ki65RGvIASR2cB6aV29DdipiHeo0NyA5P\n8S4tIOFY8uL8++KO7dn1zHpn7t07M3u+H2k0c8/9zz3n+np+c+fszHVECACQliWdHgAAoHyEPwAk\niPAHgAQR/gCQIMIfABJE+ANAggoJf9sjtg/bfrLJ+pttP1G77bb9viL6BQC0p6gz//slXTvH+ucl\nfSgifknSZyTdV1C/AIA2LCtiIxGx2/a6OdY/Vrf4mKS1RfQLAGhPJ+b8Py5pZwf6BQDUFHLmP1+2\nr5R0q6QryuwXADBTaeFv+/2StksaiojX5qjjYkMA0KKIcCv1RU77uHY7fYV9gaQHJX0sIp4704Yi\nYlHetm3b1vExsH/sH/u3+G7tKOTM3/aopIqkc2z/SNI2SX1Zjsd2SZ+SdLake2xb0nRErC+ibwBA\n64r6ts/NZ1j/CUmfKKIvAEB+/MK3RJVKpdNDWFDsX29j/9LidueLFort6LYxAUA3s63o4B98AQA9\ngvAHgAQR/gCQIMIfABJE+ANAggh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIII\nfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCgl/2yO2D9t+co6az9s+YHuP7cuK\n6BcA0J6izvzvl3Rts5W2r5N0YUS8R9JmSfcW1O/iNTUljY9n93O1lV3L2NSukrsD5hYRhdwkrZP0\nZJN190q6sW55QtKaJrWRvNHRiP7+iFWrsvvR0cZtZdcytpntJR5SYC613Gwts1t9QtMNzR3+D0v6\nQN3ydyVd3qR2gf55esTkZPaKl07dVqw4va2/P2Lv3vJqy+6vm8fW358dp5IOaYvdIUHthP+y8j5j\nzN/w8PDJx5VKRZVKpWNjKd3Bg1Jfn3T06Km2pUtPr1u+XBobK6+27P66eWzLl2fHaWDg9Oc0kPeQ\nttgdElCtVlWtVnNto6zwf0nS+XXL59XaGqoP/+QMDkrHjs1sO35csme2TU9L69eXV1t2f908tunp\n7DjNU95D2mJ3SMDsk+I77rij9Y20+lGh2U3SoKSnmqz7sKRv1R5vkPTYHNtZoA9GPeTEpO/KladP\nENe3lV3L2HLP+ZfUHRKjNqZ9nD0vH9ujkiqSzpF0WNI2SX21AW2v1dwtaUjSEUm3RsTjTbYVRYyp\n501NZZ/1BwdPfd5v1FZ2LWNre/6l5O6QENuKCJ+5su453Ra0hD8AtKad8OcXvgCQIMIfABJE+ANA\nggh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYQ/ACSI\n8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJIvwBIEGEPwAkqJDwtz1ke5/t/ba3Nlh/ju2dtvfY\nfsr2LUX0CwBojyMi3wbsJZL2S7pK0suSxiVtjIh9dTXbJK2IiNttr5b0jKQ1EfFWg+1F3jEBQEps\nKyLcynOKOPNfL+lARByKiGlJD0i6flbNjyWdVXt8lqT/bRT8AIByLCtgG2slvVC3/KKyN4R690n6\nnu2XJb1T0o0F9AsAaFMR4T8ft0t6IiKutH2hpO/Yfn9EvNmoeHh4+OTjSqWiSqVSyiABoBdUq1VV\nq9Vc2yhizn+DpOGIGKot3yYpIuKuupp/lfTXEfEfteXvSdoaEf/dYHvM+QNACzo15z8u6SLb62z3\nSdoo6aFZNROSfqM2yDWSLpb0fAF9AwDakHvaJyKO294iaZeyN5ORiJiwvTlbHdsl/Y2k+20/IcmS\n/jwi/i9v3wCA9uSe9ika0z4A0JpOTfsAAHoM4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgD\nQIIIfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJIvwBIEGEPwAk\niPAHgAQR/gCQIMIfABJUSPjbHrK9z/Z+21ub1FRs/9D207YfLaJfAEB7HBH5NmAvkbRf0lWSXpY0\nLmljROyrq1kl6T8lXRMRL9leHRGvNtle5B0TAKTEtiLCrTyniDP/9ZIORMShiJiW9ICk62fV3Czp\nwYh4SZKaBT8AoBxFhP9aSS/ULb9Ya6t3saSzbT9qe9z2xwroFwDQpmUl9nO5pF+X9A5JP7D9g4h4\ntlHx8PDwyceVSkWVSqWEIQJAb6hWq6pWq7m2UcSc/wZJwxExVFu+TVJExF11NVslrYiIO2rLX5S0\nMyIebLA95vwBoAWdmvMfl3SR7XW2+yRtlPTQrJpvSLrC9lLbb5f0q5ImCugbANCG3NM+EXHc9hZJ\nu5S9mYxExITtzdnq2B4R+2w/IulJScclbY+IvXn7BgC0J/e0T9GY9gGA1nRq2gcA0GMIfwBIEOEP\nAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4AwWbmpLGx7N7oFsR/kCBduyQ1q2Trr46u9+xo9Mj\nAhrj8g5AQaamssA/evRUW3+/dOiQNDDQuXFh8ePyDkAHHTwo9fXNbFu+PGsHug3hDxRkcFA6dmxm\n2/R01g50G8IfKMjAgDQykk31rFyZ3Y+MMOWD7sScP1CwqalsqmdwkOBHOdqZ8yf8AaDH8QdfAMC8\nEP4AkCDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYWEv+0h2/ts77e9dY66X7E9bfsjRfQLAGhP7vC3\nvUTS3ZKulXSppJtsv7dJ3Z2SHsnbJwAgnyLO/NdLOhARhyJiWtIDkq5vUPdHkr4iabKAPgEAORQR\n/mslvVC3/GKt7STbPy/phoj4R0kt/QQZAFC8ZSX18/eS6v8WMOcbwPDw8MnHlUpFlUplQQYFAL2o\nWq2qWq3m2kbuC7vZ3iBpOCKGasu3SYqIuKuu5vkTDyWtlnRE0h9ExEMNtseF3QCgBR25qqftpZKe\nkXSVpFckjUm6KSImmtTfL+nhiPhqk/WEPwC0oJ3wzz3tExHHbW+RtEvZ3xBGImLC9uZsdWyf/ZS8\nfQIA8uF6/gDQ47iePwBgXgh/AEgQ4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBI\nEOEPAAki/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJIvwBIEGEPwAkiPAHgAQR\n/gCQoELC3/aQ7X2299ve2mD9zbafqN12235fEf0CANrjiMi3AXuJpP2SrpL0sqRxSRsjYl9dzQZJ\nExHxhu0hScMRsaHJ9iLvmAAgJbYVEW7lOUWc+a+XdCAiDkXEtKQHJF1fXxARj0XEG7XFxyStLaBf\nAECbigj/tZJeqFt+UXOH+8cl7SygXwBAm5aV2ZntKyXdKumKueqGh4dPPq5UKqpUKgs6LgDoJdVq\nVdVqNdc2ipjz36BsDn+otnybpIiIu2bVvV/Sg5KGIuK5ObbHnD8AtKBTc/7jki6yvc52n6SNkh6a\nNbALlAX/x+YKfgBAOXJP+0TEcdtbJO1S9mYyEhETtjdnq2O7pE9JOlvSPbYtaToi1uftGwDQntzT\nPkVj2gcAWtOpaR8AQI8h/AEgQYQ/ACSI8AeABBH+AJAgwh8AEkT4A0CCCH8ASBDhDwAJIvwBIEGE\nP5DH1JQ0Pp7dAz2E8AfatWOHtG6ddPXV2f2OHZ0eETBvXNgNaMfUVBb4R4+eauvvlw4dkgYGOjcu\nJIkLuwFlOXhQ6uub2bZ8edYO9ADCH2jH4KB07NjMtunprB3oAYQ/0I6BAWlkJJvqWbkyux8ZYcoH\nPYM5fyCPqalsqmdwkOBHx7Qz50/4A0CP4w++AIB5IfwBIEGEPwAkiPAHgAQR/kAOXNoHvaqQ8Lc9\nZHuf7f22tzap+bztA7b32L6siH6BTmp6aR/eEdADcoe/7SWS7pZ0raRLJd1k+72zaq6TdGFEvEfS\nZkn35u13sWuUH80ypcxaxnZqedOm7NI+b7yR3W/aJE194avNL/aWc3BTE69q/Et7NTXx6szSFtoX\nqrbs/rplbD0tInLdJG2QtLNu+TZJW2fV3CvpxrrlCUlrmmwvUjc6GtHfH7FqVXY/Otq4rexaxnaq\nfWwsW5ZO3VaedTzG3nbFzMb+/ojJydwdjm7ZHf06Eqv0evTrSIxu2Z2VttC+ULVl99ctY+smtdxs\nLbtbfcJpG5B+W9L2uuXflfT5WTUPS/pA3fJ3JV3eZHsL9g/UCyYns9d7fX6sWHF6W39/xN695dWW\n3V83j61p+9veismz3j3rHWFlxCOP5Opw8m3nRb+OzCzVkdj78LPzbl+hI9GvnxReW3Z/3TK2yb1T\nnY6KGdoJ/2VlfLpo1fDw8MnHlUpFlUqlY2Mp24mLRdZfKXjp0tPrli+XxsbKqy27v24e2/Ll0ptv\nZpfy2bQpW56elkb+7ogG/vSVmRuZns7uc3R40O9Sn6Z1tL5U0xr7+svq0+p5tS/VT0/vroDasvvr\nlrEdHJvUwCWrT1tXlmq1qmq1mm8jrb5bzL4pm/b5dt3yfKZ99olpn4Y48+/+sZ2YyTlxvMbGTi2f\nnLJZufLUVE7Og8qZf/eNbTGc+bdU3HAD0lJJz0paJ6lP0h5Jl8yq+bCkb9Ueb5D02BzbW8B/ot7Q\nKD8atZVdy9hmtjd12jtC/g5PzDuvbDJHPZ/2haotu79uGVs3aSf8C7mwm+0hSZ9T9u2hkYi40/bm\n2oC212ruljQk6YikWyPi8SbbiiLG1OsaXSyy2QUky6xlbDku3pmzw6mJV3VwbFKD68+dMeXQSvtC\n1ZbdX7eMrVtwVU8ASBBX9QQAzAvhDwAJIvwBIEGEPwAkiPAHgAQR/gCQIMIfABJE+ANAggh/AEgQ\n4Q8ACSL8ASBBhD8AJIjwB4AEEf4AkCDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYQ/ACSI8AeABBH+\nAJCgXOFv+2ds77L9jO1HbK9qUHOe7e/b/h/bT9n+4zx9AgDyy3vmf5uk70bEL0j6vqTbG9S8JenP\nIuJSSb8m6Q9tvzdnvz2pWq12eggLiv3rbexfWvKG//WSvlR7/CVJN8wuiIgfR8Se2uM3JU1IWpuz\n35602P/zsX+9jf1LS97wPzciDktZyEs6d65i24OSLpP0Xzn7BQDksOxMBba/I2lNfZOkkPRXDcpj\nju28U9JXJP1J7RMAAKBDHNE0r8/8ZHtCUiUiDtv+WUmPRsQlDeqWSfqmpJ0R8bkzbLP9AQFAoiLC\nrdSf8cz/DB6SdIukuyT9nqRvNKn7J0l7zxT8Uus7AABoXd4z/7MlfVnS+ZIOSfpoRLxu++ck3RcR\nv2n7g5L+XdJTyqaFQtJfRMS3c48eANCWXOEPAOhNXfELX9u/Y/tp28dtX17Xvs72T2w/Xrvd08lx\ntqvZ/tXW3W77gO0J29d0aoxFsb3N9ot1x2yo02PKy/aQ7X2299ve2unxFM32QdtP2P6h7bFOjycv\n2yO2D9t+sq7tjD9I7RVN9q/l111XhL+yKaHfkvRvDdY9GxGX126fLHlcRWm4f7YvkfRRSZdIuk7S\nPbYXw988Plt3zHp6es/2Ekl3S7pW0qWSblqEP1L8qbIvbvxyRKzv9GAKcL+y41VvPj9I7RWN9k9q\n8XXXFeEfEc9ExAFlXyOdrefDcI79u17SAxHxVkQclHRA0mJ48fX8MauzXtKBiDgUEdOSHlB23BYT\nq0uyoAgRsVvSa7Oaz/iD1F7RZP+kFl93vXDAB2sfYx61fUWnB1OwtZJeqFt+SYvj189bbO+x/cVe\n/nhdM/sYvajFcYzqhaTv2B63/YlOD2aBtPSD1B7V0usu71c9522OH4v9ZUQ83ORpL0u6ICJeq82V\nf932L3bjj8Ta3L+eNNe+SrpH0qcjImx/RtJnJW0qf5RowQcj4hXbA8reBCZqZ5eL2WL7pkvLr7vS\nwj8irm7jOdOqfbyJiMdtPyfpYkmPFzy83NrZP2Vn+ufXLZ9Xa+tqLezrfZJ6/Y3vJUkX1C33xDFq\nRUS8Urufsv01ZVNdiy38D9teU/eD1MlOD6hIETFVtziv1103TvucnLeyvbr2BzfZfrekiyQ936mB\nFaR+Xu4hSRtt99l+l7L96+lvW9ReWCd8RNLTnRpLQcYlXVT75lmfpI3KjtuiYPvttUuvyPY7JF2j\n3j9mUvY6m/1au6X2eK4fpPaKGfvXzuuutDP/udi+QdI/SFot6Zu290TEdZI+JOnTto8p+0bC5oh4\nvYNDbUuz/YuIvba/LGmvpGlJn4ze/+HF39q+TNnxOihpc2eHk09EHLe9RdIuZSdLIxEx0eFhFWmN\npK/VLquyTNI/R8SuDo8pF9ujkiqSzrH9I0nbJN0p6V9s/75qP0jt3AjzabJ/V7b6uuNHXgCQoG6c\n9gEALDDCHwASRPgDQIIIfwBIEOEPAAki/AEgQYQ/ACSI8AeABP0/ut8+1O6BH4oAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1103ed5f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEyFJREFUeJzt3W+MXOd5nvHrpkjG26Yk3JhmDMnabSTVcgkEslAIbF3Y\nmwSMKKUIndZIRbRw7RApgURogKKBZCWAtnFRWIGR1irhqko2hlSUVJIWaWnHf2TXWgH5YJGozMS2\nSIlqsazEyuSitVTEWUQb+emHHTLzrnZX3J3hzszy+gEDznnnOec8Z4Yz954zc2ZSVUiSdMmWQTcg\nSRouBoMkqWEwSJIaBoMkqWEwSJIaBoMkqdGXYEiyP8mZJC8kuW+FmoeTnE1yKsltS27bkuTZJMf7\n0Y8kaf16DoYkW4AjwJ3AHuBgkluX1NwF3FRVtwCHgUeWLOaXgOd67UWS1Lt+7DHcAZytqnNVtQA8\nARxYUnMAeBygqp4BdibZDZDkBuBu4Lf60IskqUf9CIbrgZe6pl/ujK1Wc76r5l8Dvwx4CrYkDYGB\nvvmc5KeAC1V1CkjnIkkaoK19WMZ54Mau6Rs6Y0tr3r1MzYeBn05yNzAG/JUkj1fVR5auJIl7FJK0\nDlW1pj+6+7HHcBK4Ocl4ku3APcDSTxcdBz4CkGQv8GpVXaiqB6rqxqr6kc58X1suFC6pqk17efDB\nBwfeg9vntrl9m++yHj3vMVTVG0nuBZ5kMWimq+p0ksOLN9ejVfWFJHcneRH4HvCxXtcrSbo6+nEo\niar6EvCeJWP/fsn0vW+xjKeBp/vRjyRp/TzzeUhMTk4OuoWrajNv32beNnD7rkVZ7zGojZakRqVX\nSRoWSagBvPksSdpEDAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1\nDAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1DAZJ\nUsNgkCQ1DAZJUqMvwZBkf5IzSV5Ict8KNQ8nOZvkVJLbOmM/kOSZJN9I8u0k/6of/UiS1q/nYEiy\nBTgC3AnsAQ4muXVJzV3ATVV1C3AYeASgqv4M+LGqeh/wo8CPJ3l/rz1JktavH3sMdwBnq+pcVS0A\nTwAHltQcAB4HqKpngJ1Jdnem/7RT8wOdfr7bh54kSevUj2C4Hnipa/rlzthqNecv1STZkuQbwHeA\nmap6rg89SZLWaeBvPlfV9zuHkm4APpDkg4PuSZKuZVv7sIzzwI1d0zd0xpbWvHu1mqr6f0n+APib\nwNPLrWhqaury9cnJSSYnJ9fbsyRtSjMzM8zMzPS0jFRVbwtIrgOeB34CeAU4ARysqtNdNXcDv1hV\nP5VkL/BvqmpvkncAC1X1WpIx4MvAv6iq/7bMeqrXXiXpWpOEqspa5ul5j6Gq3khyL/Aki4empqvq\ndJLDizfXo1X1hSR3J3kR+B7wsc7s7wIeS5LOvP9huVCQJG2cnvcYNop7DJK0duvZYxj4m8+SpOFi\nMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDtIHm5uY4efIkc3Nzg25FWpHBIG2QY8eO\nMT4+zr59+xgfH+fYsWODbklall+JIW2Aubk5xsfHmZ+fvzw2NjbGuXPn2LVr1wA702bnV2JIQ2p2\ndpbt27c3Y9u2bWN2dnYwDUmrMBikDTAxMcHrr7/ejC0sLDAxMTGYhqRVGAzSBti1axfT09OMjY2x\nY8cOxsbGmJ6e9jCShpLvMUgbaG5ujtnZWSYmJgwFbYj1vMdgMEjSJuabz5KknhkMkqSGwSBJahgM\nkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJavQlGJLsT3ImyQtJ7luh\n5uEkZ5OcSnJbZ+yGJF9L8u0k30zyT/vRjyRp/XoOhiRbgCPAncAe4GCSW5fU3AXcVFW3AIeBRzo3\n/Tnwz6pqD/C3gF9cOq8kaWP1Y4/hDuBsVZ2rqgXgCeDAkpoDwOMAVfUMsDPJ7qr6TlWd6oz/CXAa\nuL4PPUmS1qkfwXA98FLX9Mu8+cV9ac35pTVJJoDbgGf60JMkaZ22DroBgCQ/CPwn4Jc6ew7Lmpqa\nunx9cnKSycnJq96bJI2SmZkZZmZmelpGzz/tmWQvMFVV+zvT9wNVVQ911TwCPFVVv9OZPgN8sKou\nJNkKfB74YlV9epX1+NOekrRGg/ppz5PAzUnGk2wH7gGOL6k5Dnyk0+Re4NWqutC57beB51YLBUnS\nxun5UFJVvZHkXuBJFoNmuqpOJzm8eHM9WlVfSHJ3kheB7wEfBUjyfuAfAt9M8g2ggAeq6ku99iVJ\nWp+eDyVtFA8lSdLaDepQkiRpEzEYJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAY\nJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEkN\ng0GS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEmNvgRDkv1JziR5Icl9K9Q8nORsklNJ3tc1Pp3kQpI/\n7kcvkqTe9BwMSbYAR4A7gT3AwSS3Lqm5C7ipqm4BDgP/ruvmz3bmlTaNubk5Tp48ydzc3KBbkdas\nH3sMdwBnq+pcVS0ATwAHltQcAB4HqKpngJ1Jdnem/xD4bh/6kIbCsWPHGB8fZ9++fYyPj3Ps2LFB\ntyStST+C4Xrgpa7plztjq9WcX6ZGGnlzc3McOnSI+fl5XnvtNebn5zl06JB7DhopWwfdwFpMTU1d\nvj45Ocnk5OTAepGWMzs7y/bt25mfn788tm3bNmZnZ9m1a9cAO9O1YmZmhpmZmZ6WkarqbQHJXmCq\nqvZ3pu8Hqqoe6qp5BHiqqn6nM30G+GBVXehMjwOfq6ofXWU91Wuv0tU2NzfH+Ph4EwxjY2OcO3fO\nYNBAJKGqspZ5+nEo6SRwc5LxJNuBe4DjS2qOAx/pNLkXePVSKHSkc5FG2q5du5ienmZsbIwdO3Yw\nNjbG9PS0oaCR0vMeAyx+XBX4NItBM11Vn0xymMU9h0c7NUeA/cD3gI9V1bOd8aPAJPBDwAXgwar6\n7DLrcI9BI2Nubo7Z2VkmJiYMBQ3UevYY+hIMG8FgkKS1G9ShJEnSJmIwSJIaBoMkqWEwSJIaBoMk\nqWEwSJIaBoMkqWEwSJIaBoM0YP52g4aNwSANkL/doGHkV2JIA+I3sWoj+JUY0gi59NsN3S79doM0\nSAaDNCATExO8/vrrzdjCwgITExODaUjqMBikAfG3GzSsfI9BGjB/u0FXk7/HIElq+OazJKlnBoMk\nqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoM0pPzJTw1KX4Ihyf4k\nZ5K8kOS+FWoeTnI2yakkt61lXula409+apB6DoYkW4AjwJ3AHuBgkluX1NwF3FRVtwCHgUeudF61\nlvsrcqW/LNcyfrVqN3p9w9JbL+bm5jh06BDz8/O89tprzM/Pc+jQIebm5kbyvhi13vq1jJFWVT1d\ngL3AF7um7wfuW1LzCPAPuqZPA7uvZN6u2+pad/To0RobG6udO3fW2NhYHT16dNmxlWr7sYxhXt+w\n9NarEydO1M6dOwu4fNmxY0d94hOfGLn7YtR669cyhknntXNtr+trneFNC4C/DzzaNf2PgIeX1HwO\n+Ntd018Bbr+Sebtuu0p322i4ePFijY2NNS8Wb3vb2940NjY2Vs8999wVj69lGcO8vmHp7eLFiwN/\nrIflvhi13vq1jH78H+gn1hEMWxmMNf1oxCVTU1OXr09OTjI5OdmndobfpR+On5+fvzx23XXXvalu\n27ZtnDhx4k21K42vZRnDvL5h6W12drbnX2G79JOfhw4dYtu2bSwsLPDAAw/wqU99qu/bca0+Tldr\nff36P9CLmZkZZmZmelvIWpNk6YXFw0Ff6pq+kkNJZ/iLQ0mrztt121XK09HgHsNo9NbPvxYvXrxY\nJ06cqIsXLy77+A/7fTFqvbnH0N9DSdcBLwLjwHbgFPDeJTV3A39QfxEkX7/SebuWcRXvutFw6Xjm\njh07Lh/PXG5spdp+LGOY1zcsvW3k4z/s98Wo9davZQyTgQTD4nrZDzwPnAXu74wdBv5JV82RTgj8\nEXD7avOusI6rdseNku6/IlcbW+v41ard6PUNS29XyyjeF6PWW7+WMSzWEwxZnG/4JalR6VWShkUS\nqmpN7+t65rMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEw\nSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIa\nBoMkqWEwSJIaBoMkqWEwSJIaPQVDkrcneTLJ80m+nGTnCnX7k5xJ8kKS+7rGP5zkW0neSHJ7L71I\nkvqj1z2G+4GvVtV7gK8BH19akGQLcAS4E9gDHExya+fmbwI/AzzdYx+SpD7pNRgOAI91rj8GfGiZ\nmjuAs1V1rqoWgCc681FVz1fVWSA99iFJ6pNeg+GdVXUBoKq+A7xzmZrrgZe6pl/ujEmShtDWtypI\n8hVgd/cQUMCvLlNefepLkjQgbxkMVbVvpduSXEiyu6ouJPlh4OIyZeeBG7umb+iMrdnU1NTl65OT\nk0xOTq5nMZK0ac3MzDAzM9PTMlK1/j/ykzwE/N+qeqjzaaO3V9X9S2quA54HfgJ4BTgBHKyq0101\nTwH/vKr++yrrql56laRrURKqak3v4/b6HsNDwL4kl174P9lp5F1JPg9QVW8A9wJPAt8GnrgUCkk+\nlOQlYC/w+SRf7LEfSVKPetpj2EjuMUjS2g1ij0GStMkYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEg\nSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoY\nDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkRk/BkOTtSZ5M8nySLyfZuULd\n/iRnkryQ5L6u8V9PcjrJqST/OcmOXvqRJPWu1z2G+4GvVtV7gK8BH19akGQLcAS4E9gDHExya+fm\nJ4E9VXUbcHa5+a8VMzMzg27hqtrM27eZtw3cvmtRr8FwAHisc/0x4EPL1NwBnK2qc1W1ADzRmY+q\n+mpVfb9T93Xghh77GVmb/T/nZt6+zbxt4PZdi3oNhndW1QWAqvoO8M5laq4HXuqafrkzttTPAV/s\nsR9JUo+2vlVBkq8Au7uHgAJ+dZnyWk8TSX4FWKiqo+uZX5LUP6la12v54szJaWCyqi4k+WHgqap6\n75KavcBUVe3vTN8PVFU91Jn+KPDzwI9X1Z+tsq71NypJ17Cqylrq33KP4S0cBz4KPAT8Y+C/LlNz\nErg5yTjwCnAPcBAWP60E/DLwgdVCAda+YZKk9el1j+GvAr8LvBs4B/xsVb2a5F3Ab1bV3+3U7Qc+\nzeJ7GtNV9cnO+FlgO/B/Oov8elX9wrobkiT1rKdgkCRtPkN95nOSDyf5VpI3ktzeNT6e5E+TPNu5\nfGaQfa7XStvXue3jSc52TgD8yUH12C9JHkzyctdjtn/QPfXDSidvbhZJZpP8UZJvJDkx6H56lWQ6\nyYUkf9w1dkUn6o6CFbZvzc+9oQ4G4JvAzwBPL3Pbi1V1e+cyqoeflt2+JO8FfhZ4L3AX8Jkkm+E9\nlt/oesy+NOhmevUWJ29uFt9n8QMm76uqOwbdTB98lsXHq9tbnqg7QpbbPljjc2+og6Gqnq+qsyx+\nRHapkX+hXGX7DgBPVNWfV9Usi2eFb4Yn5cg/ZkusePLmJhKG/HViLarqD4HvLhm+khN1R8IK2wdr\nfO6N8gM+0dkteirJ3xl0M3229KTA8yx/UuCoubfzvVi/Ncq7612u9OTNUVbAV5KcTPLzg27mKrmS\nE3VH3Zqee71+XLVnq5xA9ytV9bkVZvvfwI1V9d3Osfn/kuRvVNWfXOV212yd2zeSVttW4DPAr1VV\nJfmXwG8Ahza+S63R+6vqlSS7WAyI052/SjezzfaJnDU/9wYeDFW1bx3zLNDZXaqqZ5P8D+CvA8/2\nub2erWf7WNxDeHfX9A2dsaG2hm39TWAzhOJ54Mau6ZF4nNaiql7p/DuX5PdZPHy22YLhQpLdXSfq\nXhx0Q/1UVXNdk1f03BulQ0mXj5EleUfnjT+S/AhwM/A/B9VYn3QfAzwO3JNke5K/xuL2jfQnQjpP\nuEv+HvCtQfXSR5dP3kyyncWTN48PuKe+SfKXkvxg5/pfBn6SzfG4hTc/3z7aub7SibqjpNm+9Tz3\nBr7HsJokHwL+LfAO4PNJTlXVXcAHgF9L8jqLn5o4XFWvDrDVdVlp+6rquSS/CzwHLAC/UKN/wsmv\nJ7mNxcdrFjg82HZ6V1VvJLmXxa+Pv3Ty5ukBt9VPu4Hf73wdzVbgP1bVkwPuqSdJjgKTwA8l+V/A\ng8Angd9L8nN0TtQdXIe9WWH7fmytzz1PcJMkNUbpUJIkaQMYDJKkhsEgSWoYDJKkhsEgSWoYDJKk\nhsEgSWoYDJKkxv8HGHDC/qVfqnkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110886048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find range of possible lambda\n",
    "lambda_upper = int(np.min([2*N1/margin, 1+2*N2/margin]))+1\n",
    "lambda_lower = int(np.max([-2*N1/margin, 1-2*N2/margin]))\n",
    "\n",
    "print(lambda_lower, lambda_upper)\n",
    "\n",
    "fisher_pvalues = []\n",
    "cvr_pvalues = []\n",
    "nocvr_pvalues = []\n",
    "for lam in np.arange(lambda_lower, lambda_upper+1, 0.5):\n",
    "#for lam in np.arange(0, 1, 0.05):\n",
    "    cvr_pvalues.append(np.min([1, cvr_pvalue(lam)]))\n",
    "    nocvr_pvalues.append(nocvr_pvalue(1-lam))\n",
    "    fisher_pvalues.append(fisher_combined_pvalue([cvr_pvalues[-1], nocvr_pvalues[-1]]))\n",
    "    print(lam, fisher_pvalues[-1])    \n",
    "    \n",
    "plt.scatter(np.arange(lambda_lower, lambda_upper+1, 0.5), cvr_pvalues, color='r')\n",
    "plt.scatter(np.arange(lambda_lower, lambda_upper+1, 0.5), nocvr_pvalues, color='b')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(np.arange(lambda_lower, lambda_upper+1, 0.5), fisher_pvalues, color='black')\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'allocation lambda': 0.93248039808125949,\n",
       " 'max_pvalue': 0.027011101425504891,\n",
       " 'min_chisq': 10.960538840556488}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maximize_fisher_combined_pvalue(N=(N1, N2), overall_margin=margin, \n",
    "                          pvalue_funs=(cvr_pvalue, nocvr_pvalue))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}