{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Guided Sequential Search\n",
    "Kyle Cranmer, June 2-15\n",
    "\n",
    "##The Scenario\n",
    "\n",
    "We imagine that we see an excess in some initial dataset, which we refer to as \"run 1\", and we want to know how to best use a subsequent \"run 2\" dataset for discovery. We are imagining that either:\n",
    "\n",
    " * several signal regions were defined in Run 1, some of which had a significant upward fluctuation, or \n",
    " * after investigating an unexpected excess that the events can be subdivided into multiple (non-overlapping / disjoint) categories.\n",
    " \n",
    "In both cases, the difficulty is that one does not have a signal model that allows one to relate the excess in the different categories or \"bins\". This leads to a few strategies:\n",
    "\n",
    " 1. the inclusive strategy: merge the bins with a significant excess into one inclusive bin.\n",
    " 1. the ideal strategy: one performs a combined likelihood fit with the correct signal model. This is not possible usually, but we use it for comparison.\n",
    " 1. the nested strategy: where one starts with the bin with the highest significance excess and defines additional bins by merging with the bin with the next highest significance. This gives a nested set of selection regions. The problem here is that one has multiple tests and a look-elsewhere effect / trials factor to correct for, and the events are not statistically independent, so the analysis is tricky.\n",
    " 1. the data guided sequential search strategy (DGSS): where one uses the excess in run 1 to define an emperical signal model to be used in run 2.\n",
    " \n",
    "The basic idea of DGSS is that if the excess in run 1 is from the presence of signal and representative of the true signal, that it is a nearly optimal way to combine (ie. approximate the ideal strategy). While there are statisitcal fluctuations, the emperical signal model is not biased (ie. expectation = true signal) when there is no selection on run 1. When one only uses this procedure with a significant excess in run 1 it will introduce some bias on signal estimate. Thus we resort to some numerical simulations to check the power of this method in comparison to the ideal strategy and the inclusive strategy.\n",
    "\n",
    "While it has not yet been formalized, we anticipate a recursive argument can be made in going from $N\\to N+1$ bins. Therefore, we focus on the base of the recursion $1\\to 2$ bins.\n",
    "\n",
    "\n",
    "##The Model\n",
    "\n",
    "We start with the model for run 1\n",
    "\n",
    "\\begin{equation}\n",
    "Pois(n_1^1 | s_1 + b_1) Pois(n_2^1 | s_2 + b_2)\n",
    "\\end{equation}\n",
    "\n",
    "From this run, we estimate $\\hat{s}_i = n_i = b_i$ and then use those signal estimates in run 2, where the model is:\n",
    "\\begin{equation}\n",
    "Pois(n_1^2 | \\mu \\hat{s}_1 + b_1) Pois(n_2^2 | \\mu \\hat{s}_2 + b_2)\n",
    "\\end{equation}\n",
    "\n",
    "Based on run 2 we estimate the MLE $\\hat{\\mu}$, which is solved below using sympy.\n",
    "\n",
    "Then, from the estimate of $\\hat{\\mu}$ we construct the liklelihood ratio test of the background only based on:\n",
    "\n",
    "\\begin{equation}\n",
    "q_0 = -2 \\ln \\lambda(0) = -2 \\ln \\frac{L(\\mu=0)}{L(\\hat{\\mu})}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline --no-import-all\n",
    "from sympy import Symbol, symbols, solve, lambdify\n",
    "#import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "#plt.rc('text', usetex=True)\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#observables\n",
    "(n11, n12, n21, n22) =  symbols('n11 n12 n21 n22')\n",
    "#parameters\n",
    "(s1, s2, b1, b2, mu) =  symbols('s1 s2 b1 b2 mu')\n",
    "#helpers\n",
    "(n1, n2, nu11, nu12, n21, nu22) =  symbols('n1 n2 nu11 nu12 nu21 n22')\n",
    "#estimators\n",
    "(s1hat, s2hat, muhat) = symbols('s1hat s2hat muhat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#first run\n",
    "nu11 = s1hat + b1\n",
    "nu12 = s2hat + b2\n",
    "#second run\n",
    "nu21 = muhat*s1hat+b1\n",
    "nu22 = muhat*s2hat+b2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 5)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#first run\n",
    "shatSubs = [(key, val) for key,val in solve([n11-nu11, n12-nu12], s1hat, s2hat).iteritems()]\n",
    "moreSubs = shatSubs+[(n11,10),(n12,10),(n21,10),(n22,10),(b1,5),(b2,5)]\n",
    "s1hat.subs(moreSubs), s2hat.subs(moreSubs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#second run\n",
    "[expr1, expr2] = solve( s1hat+s2hat - n21*s1hat/nu21 - n22*s2hat/nu22, muhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1, 1)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr1.subs(moreSubs),expr2.subs(moreSubs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "muhatFunc = lambdify((n11,n12,n21,n22,b1,b2), expr2.subs(shatSubs), \"numpy\")\n",
    "muhatFunc(10,10,10,10,5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "muhatFuncTruth = lambdify((n11,n12,n21,n22,b1,b2,s1hat,s2hat), expr2, \"numpy\")\n",
    "muhatFuncTruth(10,10,10,10,5,5,5,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Pois(n|\\nu) = e^{-\\nu} \\nu^n/n!$\n",
    "\n",
    "$-2 \\log L(\\nu) = 2*(\\nu - n \\ln \\nu)$\n",
    "\n",
    "$q(0) = -2 \\log L(\\nu)/L(\\hat\\nu) = 2*(\\nu(0)-\\hat{\\nu}) - n \\ln \\nu(0) / \\hat\\nu )$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def inclusive_muhatFunc(n11,n12,n21,n22,b1,b2):\n",
    "    shat = 1.*(n11+n12-b1-b2)+.01\n",
    "    return np.maximum(-7,np.minimum(7,np.divide(n21+n22-b1-b2,shat)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def q0inclusive(n11,n12,n21,n22,b1,b2):\n",
    "    #shat = 1.*(n11+n12-b1-b2)+.01\n",
    "    #muhat = np.divide(n21+n22-b1-b2,shat)\n",
    "    return 2*(-(n21+n22-b1-b2) - (n21+n22)*np.log((b1+b2)/(n21+n22)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def q0DGSS(n11,n12,n21,n22,b1,b2):\n",
    "    muhat = muhatFunc(n11,n12,n21,n22,b1,b2)\n",
    "    s1hat = 1.*(n11-b1)\n",
    "    s2hat = 1.*(n12-b2)\n",
    "    q_1 = 2*( -muhat*s1hat - (n21)*np.log((b1)/(b1+s1hat*muhat)))\n",
    "    q_2 = 2*( -muhat*s2hat - (n22)*np.log((b2)/(b2+s2hat*muhat)))\n",
    "    return q_1+q_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def q0Truth(n11,n12,n21,n22,b1,b2,s1hat,s2hat):\n",
    "    muhat = muhatFuncTruth(n11,n12,n21,n22,b1,b2,s1hat,s2hat)\n",
    "    #s1hat = s1\n",
    "    #s2hat = s2\n",
    "    q_1 = 2*( -muhat*s1hat - (n21)*np.log((b1)/(b1+s1hat*muhat)))\n",
    "    q_2 = 2*( -muhat*s2hat - (n22)*np.log((b2)/(b2+s2hat*muhat)))\n",
    "    return q_1+q_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate toy experiments *with*  Run-1 filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64493 b experiments\n"
     ]
    }
   ],
   "source": [
    "#b experiments, with filtering\n",
    "nExp = 3000000\n",
    "btruth=[(b1,100),(b2,100)]\n",
    "struth=[(s1hat,0),(s2hat,0),(muhat,0)]\n",
    "truth=btruth+struth\n",
    "n11_expts = np.random.poisson(nu11.subs(truth),nExp)\n",
    "n12_expts = np.random.poisson(nu12.subs(truth),nExp)\n",
    "#n21_expts = np.random.poisson(nu21.subs(truth),nExp)\n",
    "#n22_expts = np.random.poisson(nu22.subs(truth),nExp)\n",
    "temp = np.zeros(nExp*2).reshape(nExp,2)\n",
    "temp[:,0] = n11_expts\n",
    "temp[:,1] = n12_expts\n",
    "#temp[:,2] = n21_expts\n",
    "#temp[:,3] = n22_expts\n",
    "temp = np.array(filter(lambda x: x[0]>b1.subs(truth)*1.1 and x[1]>b2.subs(truth)*1.1, temp))\n",
    "n11_expts = temp[:,0]\n",
    "n12_expts = temp[:,1]\n",
    "#n21_expts = temp[:,2]\n",
    "#n22_expts = temp[:,3]\n",
    "n21_expts = np.random.poisson(nu21.subs(truth),n11_expts.size)\n",
    "n22_expts = np.random.poisson(nu22.subs(truth),n11_expts.size)\n",
    "muhat_b_expts = muhatFunc(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "muhat_b_expts = filter(lambda x: x<1e6, muhat_b_expts)\n",
    "muhat_b_expts =  np.maximum(-10,np.minimum(10,muhat_b_expts))\n",
    "inclusive_muhat_b_expts = inclusive_muhatFunc(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "inclusive_q0_b_expts = q0inclusive(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "q0_b_expts = q0DGSS(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "nu11.subs(truth),nu12.subs(truth),nu21.subs(truth),nu22.subs(truth)\n",
    "print len(n22_expts), 'b experiments'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#b truth experiments, need to know true signal rate\n",
    "#experiments are still the b-only ones\n",
    "struth=[(s1hat,50),(s2hat,1),(muhat,1)]\n",
    "truth=btruth+struth\n",
    "truth_q0_b_expts = q0Truth(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)),float(s1hat.subs(truth)),float(s2hat.subs(truth)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.629604480626714"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(truth_q0_b_expts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "68538"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#s experiments, with filtering\n",
    "nExp = 400000\n",
    "#ok, now make s+b experimentes\n",
    "n11_expts = np.random.poisson(nu11.subs(truth),nExp)\n",
    "n12_expts = np.random.poisson(nu12.subs(truth),nExp)\n",
    "#n21_expts = np.random.poisson(nu21.subs(truth),nExp)\n",
    "#n22_expts = np.random.poisson(nu22.subs(truth),nExp)\n",
    "temp = np.zeros(nExp*2).reshape(nExp,2)\n",
    "temp[:,0] = n11_expts\n",
    "temp[:,1] = n12_expts\n",
    "#temp[:,2] = n21_expts\n",
    "#temp[:,3] = n22_expts\n",
    "temp = np.array(filter(lambda x: x[0]>b1.subs(truth)*1.1 and x[1]>b2.subs(truth)*1.1, temp))\n",
    "n11_expts = temp[:,0]\n",
    "n12_expts = temp[:,1]\n",
    "#n21_expts = temp[:,2]\n",
    "#n22_expts = temp[:,3]\n",
    "n21_expts = np.random.poisson(nu21.subs(truth),n11_expts.size)\n",
    "n22_expts = np.random.poisson(nu22.subs(truth),n11_expts.size)\n",
    "muhat_s_expts = muhatFunc(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "muhat_s_expts = filter(lambda x: x<1e6, muhat_s_expts)\n",
    "muhat_s_expts =  np.maximum(-10,np.minimum(10,muhat_s_expts))\n",
    "inclusive_muhat_s_expts = inclusive_muhatFunc(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "inclusive_q0_s_expts = q0inclusive(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "q0_s_expts = q0DGSS(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)))\n",
    "truth_q0_s_expts = q0Truth(n11_expts,n12_expts,n21_expts,n22_expts,float(b1.subs(truth)),float(b2.subs(truth)),float(s1hat.subs(truth)),float(s2hat.subs(truth)))\n",
    "nu11.subs(truth),nu12.subs(truth),nu21.subs(truth),nu22.subs(truth)\n",
    "len(muhat_s_expts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots for the signal strength estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD/CAYAAAAZg9YLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAFI9JREFUeJzt3cFvHNd9wPEvaYqktaRFcVV5I8eyRcVBGbhBuFEL9BBg\n",
       "BUmofIjRQIRy9cGUgCDHCkjQFN6DDxbq/gOWDzkECBxItyIXWcQiRdGgcGAdiqqwZQapUcZVRYlR\n",
       "JZWRRLGHmSWX613u7uzuzHDn+wEWO/PjivPT8vHH2Tfz3gNJkiRJkiRJkiRJkiRJu8NQG69ZCJ+P\n",
       "Aj8Kt88Aq8AMcKkHMUlSjEZafP0E8CHwW+AX4f6d8GvXCAr4HFt/RKLEPu72PyFJ6sxwi6/PACfD\n",
       "7aVw//sEZ+7V2EngbBcxSVLMWp3513bLFIEPgG8DKzXxPDDF1ieCTmOSpJi1OvOvKgK/YauLpp1r\n",
       "BZKklGp15l91AvhxuL0KTIfbU2x9Cug0tp/tnyCqbhJcXJb64TPgawkd27atfup52z5Xs32C4CJt\n",
       "9Q6gC8C3uozV2+hV4j1WTjqBJspJJ9BEOekEmkiyfdm2O1NOOoEmykkn0ERH7atVt89J4B2CM5Y7\n",
       "4Tevdv2cIPgUcL3LmCQpZq26fT5kq5umVvVC8LUexdRzs/PB843LyeYhKY3a7fMXVJJOoIlK43Ah\n",
       "vJPqRmyJ1KkkdWB1rJJ0Ak1Ukk6giUrSCQyqtPaL7jLHzwcP1bHPX4Oqp33+kqQBZPGXpAyyz1/S\n",
       "rjEL8wXIfwErN8CbGbrgmb+kXaMA+UVYLjg1TNcs/pKUQRZ/Scogi7+kXWEW5teCSSbVA17wlbQr\n",
       "FCC/Abmk8xgUnvlLUgZZ/CUpgyz+kpRBFn9JyiCLvyRlkMVfkjKo3eJ/scn+Qk3sDMEKXVFikqQY\n",
       "tVP8zxEU7FoLwKcECwbD1sCL6upccx3EJEkxa6f4vwcs1cUWgFeAxXD/LHA33F4iWPv3LME6va1i\n",
       "kqSYRe3znybourkQ7k8RLPBele8gpr65VwxW86qu5ytJgajF/xJB102e4I8AwFBPMlIPTeZgcXlr\n",
       "PV9JCkSZ22eB4Oz9CrACzBB05UyHX58K47SI7a+JSZJiFKX4LwEfhdt54Gq4f4zg08BMGBtqETsS\n",
       "xhop12xXwocURSl8pEW5ZruCbVvRleiibbdT/OcJCvabwPsEhbt6989t4Hq4fYygC2g1Qqxeud3/\n",
       "gNRChe0F9q1k0thUTvj4GhwVumjb7RT/y3x5rcwrDV53KXy+FiEmSYqRI3wlKYMs/pKUQRZ/Scog\n",
       "i78kZZDFX5IyyAXcB87sfDCid60IY0knIymlPPMfOIV8MKXDWC7pTCSll8VfkjLI4i9JGWTxl6QM\n",
       "svhLUgZZ/CUpgyz+knade1CcDWYcVkTe57+bHWCeCfLcZ4XbX5p5VRpYk5B7DvI3kk5kF/PMfzeb\n",
       "IM8bLDPhWsgaLLOHmD/+Dc7PHvLsvl8s/pJSpzBFfvFvWS5MeWLTL+0W/4t1+2cIVuNa6FFMkhSj\n",
       "dor/ObaWbQQohs/VlbjmuoxJkmLWTvF/j2DR9qqzwN1wewk4GcZWI8YkSTGL0uc/Bdyp2c93GVMU\n",
       "B5jnyeYnKWlgzB5ifu2Rbbvfot7qOdTTLNS54A4fZ+7UwClMkd/YsG33W5Qz/1VgOtyeAlYixvaH\n",
       "MXXrjxR5mfMc8LY4DZZ7/0fR2z37I8qZ/wfAMYKLtjPAVYJPAp3GjoSxRso125XwoWbGyPEGy/yU\n",
       "Q9xOOpnUKYWPtCjXbFewbe9ocpzcc8+Sv7GcdCapVKKLtt1O8Z8nKNhvAu8DH4f7JwjO5K+Hr+sm\n",
       "Vq/cwf9B2kmF7QX2rWTS2FRO+PgaHBW6aNvtFP/L4aPWpfD5Wo9ikqQYOcJXkjLI4i9JGWTxl6QM\n",
       "svhnwr0izHq7nKRNFv9MmMxBwdHUkjZZ/CWl2to9isfh/BpO+dBLFn9JqTb2lNwiLI85nUlPWfwl\n",
       "KYMs/pKUQRZ/ScqgqFM6KykHmGeCPE8ofumn90eK5G7CA5wGS7vO7CHmC1Pk1x5RHNuzFX+wwYs/\n",
       "+SqnH9zlxdwDPk8uw8Himf9uM0GeN1hmpMHFrzFyPPvIi2LalaqLto/t2d62c+OMvX2KldxexpLK\n",
       "bRBZ/CUpgyz+kpRBFn9JyiCLvyRlkMVfkjIoavG/GD4v1MTOECzPGCUmSYpR1OK/AHwKfBbuVydc\n",
       "qi7NONdBTJI6do9gwrfZYJ1xdaib4v8KsBjunwXuhttLwMkwttpGTJI6Nkkw4VsBnK48gqjFf5qg\n",
       "6+ZCuD8F3Kn5er6DmCQpZlGnd7gUPp8i+CMAMNR9OpKkOEQp/gsEZ+9XgBVghqArZzr8+lQYp0Vs\n",
       "f02sXrlmuxI+pChK4SMtyjXbFWzbiq5EF207SvFfAj4Kt/PA1XD/GMGF3JkwNtQidiSMNVKOkJfU\n",
       "SIXtBfatZNLYVE74+BocFbpo21GK/zWC2zUBbgPXw+1jBF1AqxFikqQYRe3zv9IgVr0OcC1CTFJG\n",
       "1U7lDE5HHhdH+O4mB5jniYtYa7A0m8q53sMhph/keDGuvAadxX83mSDfcB7/bR68yMHXTzP+ed0v\n",
       "yb0iHD8Psw6I0a60d5wR5/TvHYv/IHm85zDPDH2FP39hhfHHdb8kkzlYXIaCYyskWfwHysboKAy7\n",
       "NKekliz+kpRBFn9JyiCLvyRlkMVfkjLI4i9JGWTxl7Sr3YOiC7p0zuK/Gxxgnpc539Ho3qH7018e\n",
       "6CWly+wh5sNpHdrycIjpn3yV07UjfSch54IunbP47wYT5HmD5daje2uMro98eaCXlC6FKfKtpnWo\n",
       "tXeckbdPseJI3+5Z/CUpgyz+mXOv6Pw+2i3WH3H46dPNBaDUQxb/zJnMOb+PdotnnjI6vBF56nnt\n",
       "wOIvSRnkX9RBNnR/moOvn+be5y+y9uLnSacjVdUu4DK2p/N//3CI6Y0cL/Kg97llRRJn/mcIlnFc\n",
       "SODYu083C7iMro80nt5ZSla7C7g0Uz+3/z0oHofz3u/fvriLf7WIVZdwnIv5+N0oJXLUVgu4fMir\n",
       "nX/TWBZ2KfXxe6u3Skkn0MiPP2i/bU9CbhGWY7rfvxTDMfou7uJ/Frgbbi8BJ2M+fjdKSSfQ0O/4\n",
       "s5av+dKAr1gWdin18Xurt0pJJ9DIr/5j57bdaMBXTEoxH68v4i7+U8Cdmn3vOmmmk1G9j/cc5lcv\n",
       "nWZ9uPEtcU0HfHnbp+LX6ajeZhoN+HKqh/YlccF3qI3XfAf4p34nkloHmGecU7zBP/JTvtPy9Ruj\n",
       "ozz72gr8vPnPs3rx9+G/zfFk323WIPgEwCk4nocvVuDG5Z79H9TMV4DfJ51EEuou8u7Y17/+iMND\n",
       "Txh/ClOtzlAfDjH9aIq5n0zwh+G7/OneB3Ac8l/AK8Cvb4DtuoF2CnEvvQNcJejznweOAH9f95qb\n",
       "wNGY81J2fAZ8LaFj27bVT0m27Zbm2LrL5wLwrQRzkaTMeibm430BfBvYDzwL/DLm40uSJEmSJEmS\n",
       "JEmSJEmSJEmSJEmSJEmSJEmSJElSz7UzpXN1Fs6jwI/C7TPAKjADXOpBTJIUo1aLuZwAPgR+C/wi\n",
       "3K+uxHWNoIDPsfVHJErs427/E5KkzrRaJGeGrXV2l8L97xOcuVdjJwnW5o0akyTFrNWZf223TBH4\n",
       "gGA+/pWaeJ7Ga/O2G5MkxazdBdyLwG/Y6qKJe/lHSVIPtbuA+wngx+H2KjAdbk+x9Smg09h+tn+C\n",
       "qHKdU/WTa/hqUPW8bZ+r2T5B43V4u4nV2+hV4j1WTjqBJspJJ9BEOekEmkiyfdm2O1NOOoEmykkn\n",
       "0ERH7atVt89J4B2CM5Y74Tevdv2cIPgUcL3LmCQpZq26fT5kq5umVvVC8LUexdRzs/PB843LyeYh\n",
       "KY3a7fMXVJJOoIlK43AhvJPqRmyJ1KkkdWB1rJJ0Ak1Ukk6giUrSCQyqtPaL7jLHzwcP1bHPX4Oq\n",
       "p33+kqQBZLfPgDrIzSLAraQTkZRKFv8BNc2jHFj8JTVmt48kZZDFX5IyyOIvSRlkn/+AmYX5AuT/\n",
       "h4eHn7D3P5POR1I6eeY/YAqQX4TlcZ6OJp2LpPSy+EtSBln8JSmDLP6SlEEWf0nKIIu/JGVQu8X/\n",
       "YpP9hZrYGYJFWqLEJEkxaqf4nyMo2LUWgE8J1oyEYIF32FqgZa6DmCQpZu0U//eApbrYAvAKsBju\n",
       "nwXuhttLBMs/niVYqrFVTJIUs6h9/tMEXTcXwv0pgjV+q/IdxCRJMYta/C8RdN3kCf4IAAz1JCP1\n",
       "zCM4HKzmVV3PV5ICUYr/AlvXAFaAGYKunOpC71NhvFVsfxhTn2wwNgqLy1vr+UpSIMrEbkvAR+F2\n",
       "Hrga7h8j+DQwE8aGWsSOhLFGyjXbFVwwWdGVwkdalGu2K9i2FV2JLtp2O8V/nqBgvwm8T1C4q2f+\n",
       "t4Hr4fYxgi6g1QixeuV2/wNSCxW2F9i3kkljUznh42twVEhX2+5aRyvQa7vjcH4DvnuM3K+PcPjX\n",
       "sPHdoN9foSTbl21b/dRR+3KEryRlkMVfkjLI4i9JGWTxl6QMsvhLUgZZ/CUpg6IM8lKK3WKi+Db5\n",
       "bz7hzjTb51GSpE2e+Q+YdXK5PK+tbDDsH3ZJTVn8JSmDLP6SlEEWf0nKIIu/JGWQxV+SMsg7QgbY\n",
       "MPenX+X107e4OXIr6WQkpYpn/gNsgvWRH/DCyjSPcknnInXkAPO8zHkO4BKkfWLxl5Q6E3s59dJf\n",
       "8M2JvZxKOpdBZfGXlDq5PeRe+0tWcnvwU2uftFv8L9btnyFYinGhRzFJUozaKf7n2FqzF6AYPl8L\n",
       "n+e6jEmSYtZO8X8PWKrZPwvcDbeXgJNhbDViTJIUsyh9/lNsny0y32VMkjbtO8i7a483ewnUJ1Ev\n",
       "+A71NAtJCo1PUuAZxpPOY9BFKf6rwHS4PQWsRIztD2OS1NDDNYr7DvJu0nkMoigjfD8AjhFctJ0B\n",
       "rhJ8Eug0diSMNVKu2a6EDymKUvhIi3LNdgXb9o6GRxkfH6fwB4eoN1Kii7bdTvGfJyjYbwLvAx+H\n",
       "+ycIzuSvh6/rJlav3MH/QdpJhe0F9q1k0thUTvj4GhwVumjb7RT/y+Gj1qXw+VqPYpKkGDnCV5Iy\n",
       "yOIvSRlk8ZekDLL4Z8AjOAyzTo0raZPFPwM2GBuFgqOpJW2y+A+QWZgf4eHhpPOQlH4W/wFSgPw4\n",
       "T0eTzkNS+ln8JSmDLP6SlEFR5vaRpJ7bd5B3xycprD2myDNJZzP4PPOXlArjkxS+90OW6qdzfrhG\n",
       "8fmj/MzZPXvL4i8p1YZHGf/eD1kan6SQdC6DxOIvSRlk8ZekDLL4S1IGWfwlKYMs/pKUQVGL/8Xw\n",
       "eaEmdoZgecYoMUlSjKIW/wXgU+CzcL8YPleXZpzrIKY+G+b+9EFuFlu/UlJWdFP8XwEWw/2zwN1w\n",
       "ewk4GcZW24ipzyZYH5nmUS7pPCSlR9TiP03QdXMh3J8C7tR8Pd9BTJIUs6hz+1wKn08R/BEAGOo+\n",
       "HUlSHKIU/wWCs/crwAowQ9CVMx1+fSqM0yK2vyZWr1yzXQkfUhSl8JEW5ZrtCrZtRVeii7Ydpfgv\n",
       "AR+F23ngarh/jOBC7kwYG2oROxLGGilHyEtqpML2AvtWMmlsKid8fA2OCl207SjF/xrB7ZoAt4Hr\n",
       "4fYxgi6g1QgxSVKMovb5X2kQq14HuBYhJimjts3jH/QsKAaO8JWUqGbz+Nd7uEbROf17x+KfEY/g\n",
       "MBw/D7PzSeciRTE8yrhz+veOxX+A3GKi+ITh6UZf22BsFBaXoeDYCkkW/0GyTi63wbDrMktqyeIv\n",
       "SRlk8ZekDLL4S1IGWfwlKYO8OJgRw9yffpXXT9/i5sitpJORlDjP/DNigvWRH/DCivP6a1db5/Dx\n",
       "b3B+9hCOV+mSxV9SYvYd5N1wWoe2PF6nMHGYb94f4VQ/88oCi7+kxIxPUmg1rUOt9RFGXvgrVh7t\n",
       "wU+wXbL4S1IGWfwHwCzMH4fzIzw83Oq1wRw/zu8jZZ3FfwAUIL8Iy+M8HW312mCOH+f3kbLO4i9J\n",
       "GZRE8T9DsJLXQgLHlpQC+w7y7vNH+Vknd/rUcm7/7sVd/Ks/6OoqXnMxHz/zhrk/fZCbkX7hpF5p\n",
       "dwGXZpzbv3txF/+zwN1wewk4GfPxu1FKOoEmSp28eIL1kQkezcawsEupj99bvVVKOoFG/v1feCnp\n",
       "HJooJZ1AL8Rd/KeAOzX7u+nCYynpBJoodfoPYlrYpdTH763eKsV5sHYHdv3XJ7y809cfrlF8/ig/\n",
       "S6D7pxTz8foiiT7/oQSOOdBuMVF8m5dON1vFq151np9pPjnjbZ+KW6cDu5oZHmX8ez9kye6faOKe\n",
       "2G0VqBao/cBKk9edAa7EktEuNgvz/w3fHmJoX57XKhv8vK2f5wTrI+d5YeUf+OdTz3L/7yZ5/q//\n",
       "l/tfX2Pv71d49ZfwxQrcuNzv/DNqhqDLMxP2HeTd8UkK9+/xdYCJ5/hk7TFFnundMR6uUcz9Cf86\n",
       "8Ryf3L/H10eG+NUfbvE3vTvCYIr7LHwOOAZcAi4AV4Hrda+5CRyNOS9lx2fA1xI6tm1b/ZRk227L\n",
       "At7qKUmS+uxi3b5jTTQoIrfttI7w9Ze1PWl9X6o/vzTkdY7gfapKeqyJbbu1tL4naWrX0GXbTmPx\n",
       "T9svayNpaARpfF+qFoBPCfogk/Ye2y+wJjnWJO1t23a9szS1a+iybaex+Kfpl7WZNDSCNL4vVQvA\n",
       "K8Bi0ok0kORYk7S3bdv1ztLcrqHDtp3G4l8vjQPD0tAI0vi+VE0TfGy/kHQiTaRlrEnafoa2652l\n",
       "vV1DB217NxR/SM8va1VaGkHa3peqSwQf2/ME71OatDvWJC5p+hnarneW5nYNHbbtuAd5QeP+xDs0\n",
       "H9SVtl9WCBoBwCmCRnBth9f2SxrfFwh+vtWf5wrBoKYk3p9mPiAYa3INOEIw1qRXdnvbtl03l/Z2\n",
       "DR227SSK/6XWL9mmn7+szez0S3yO4IefdCNI4n1pxxLwUbidJ/m85gnepzeB94GPw/0TBIWmfpBh\n",
       "N9Letm3X0aWtXUO8bTsW8wQN8s2aWJoGhp0A9oXb7wDfSjCXNL0vtc6ED4fYb5fmtm27bm2g2nVa\n",
       "+9bSrnq73hFwQQkNDNu1JEmSJEmSJEmSJEmSJEmSJElSAv4f14+e6thfFfUAAAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fc3e190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#figsize(600,600)\n",
    "edges = np.linspace(-10,10,100)\n",
    "fig, ax = plt.subplots(2,2,sharex=True, sharey=True)\n",
    "#plt.figure(None,figsize=(1800,1800))\n",
    "\n",
    "# DGSS s+b vs. b experiments\n",
    "contents, edges, patches = ax[0][0].hist(muhat_s_expts,bins=edges,alpha=0.3,color='b')\n",
    "contents, edges, patches = ax[0][0].hist(muhat_b_expts,bins=edges,alpha=0.3,color='g')\n",
    "\n",
    "# inclusive s+b vs. b experiments\n",
    "contents, edges, patches = ax[0][1].hist(inclusive_muhat_s_expts,bins=edges,alpha=0.3,color='r')\n",
    "contents, edges, patches = ax[0][1].hist(inclusive_muhat_b_expts,bins=edges,alpha=0.3,color='orange')\n",
    "\n",
    "# comparison for s+b experiments\n",
    "contents, edges, patches = ax[1][0].hist(muhat_s_expts,bins=edges,alpha=0.3,color='b')\n",
    "contents, edges, patches = ax[1][0].hist(inclusive_muhat_s_expts,bins=edges,alpha=0.3,color='r')\n",
    "\n",
    "# comparison for b experiments\n",
    "contents, edges, patches = ax[1][1].hist(inclusive_muhat_b_expts,bins=edges,alpha=0.3,color='orange')\n",
    "contents, edges, patches = ax[1][1].hist(muhat_b_expts,bins=edges,alpha=0.3,color='g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes to self:\n",
    "\n",
    "current filtering is s1hat >0 and s2hat > 0. Instead we might want cumulative n1+n2-b1-b2>0, but then need to make this method ignore bin if sihat<0.\n",
    "\n",
    "Might want to cut harder than 0\n",
    "\n",
    "muhat not best test statistic. If shat small and n2-b2 also small, can get muhat large. Should move to lambda(0).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plots of the q_0 test statistic for inclusive, ideal, and DGSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "testSize=0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13.0199453786 0.205169274628\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEWCAYAAACEz/viAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAGZZJREFUeJzt3W+MHOV9wPGvExMIbvDZruBwQorPuBFVK+UOU1WVULHP\n",
       "po2USFVuMYpU5RVnk/cxhIAqXpCAQ14igW36ApRIMT5HeRG1iu27OK36pgHMi6QWMpxJlbguwuez\n",
       "E1JDY19fPM9459a7d7O3Ozszu9+PtNqZ3Z3d383tPr95/swzIEmSJEmSJEmSpDIZA17rcPu3C/x8\n",
       "SRW2uugAdM07wGMD/Pmd2ge8CFwAxoEj8fFafGwImAemm2z7KLADeCDnGMeAe4CDqcdaxZcl7l5q\n",
       "Fntar/Zhq1iqsh+lZaWPzEeA14G9wBzhiH809doaofCeIxSAyfZvp5bTR/k7gKOp9f1x27n4Gck2\n",
       "r8fXTaZemxSyO+JnXo2vWZt6zTHg2ax/aE6OEv7+F1KPjVDfP8lrmpkAvp5TXIlx4FXq+xtax5c1\n",
       "7l5pFnujXuzDVrFUZT+W2seKDkAtjQILwHrgOKFQhvAFP0D48d1DKKQn2njfGuEHdWe8fxa4JT63\n",
       "QPih7Uy9fhz4l/j4w4QjrFkWH5U9C/ygjRjysB+4C/ha6rEdhKPBxDyLk2riAfI/YpwmJM20VvFl\n",
       "jbtXmsXeqBf7sFUsVdmPpWYyKK954LtxeT8hKUAozPcDbwJngAcJhXNWC/F+M3Ayvu+l1POH42dA\n",
       "KPhHgFsJCemn8bWPEH5oiekYTyuThB/heBtxtmt9/IwJ6slxLXA+9Zo5wt/TaBxYF7fLM8ZGreLL\n",
       "Gjf0Zt9mUdQ+hO7sx4Fnn0F5zaWWV6WWkyakxMl4P5bxfY8QCs7D8f4Z4LnU8xeBNwg/6M2EGsFm\n",
       "QoJIx7RANhOEpHKSULWfJvxIt8aYn2t4/Vpg1xLvd5yQBBslNZWThCay4y22b4w7qenMxPWjcftW\n",
       "8XUSYyea7e9m+3aIUGMcI/wfG4/W84h9uX3YqziyyPq9HTgmg+qZJxTOiTFgE0vXDoZSyyOEH9rB\n",
       "uN2xuG36h3eIUO0fJXQM3gtMsfjHm7W6nTRxvU29g/pi/MydTV5/kdadlK3UCH9LUnAnR4DzLP7b\n",
       "13P9ftpBSHjJ8oVl4ltpjM20im99i8cbpffto/GxBwn/02lCodxYCHcr9rTGfThP+K70Ko5O96Mw\n",
       "GVTRfkLN4BDhB3WYcHSfrj3MU08SF4DHgffjczXCD3U8tU3j0dIU4WjuKqH5511CoTMeP3tffO9k\n",
       "5MhYfI+TDe+ziVAoP0I4Enyc+iifVlZyxPhOvCXWx1hmqfe1QCgYGpuz1lEvIGqEfbmclR7VrmpY\n",
       "f7VJfFnjbrZvf0i9gB1j8T7pduxpjfvw26l4ux1Hs1g62Y9S6YwBP4/LI8DpFs9BaCeeIxTWL6Re\n",
       "k97mxfj8aUJzwk9Szx2Nz81RL/waP6NxZM54fOxqfK9bUs8dpXkhuilul7Tl39LwXDdHICV9BXuB\n",
       "7anHx1O37U22W0vYn+PAl3OMLxkF8xMWt6m3im+5uJfat7B4YECnWsWeaLUPux3HUrGsdD+qREao\n",
       "j1TR4Bih+OGoSyl7fEupEQroTcahMtnXsJ6MNkjGsm9i8Zh1DYZJwtFdWQuKssfXSo1Qg2tVWxu0\n",
       "OFQSu1k8RcIY9WF/yZA4CD+4SUwKklSIvM8zOMDi3vtdhA5N4uPJCVNnCEPRduccjySpiV6fdDbE\n",
       "4rHqG6iPad9KGBkjSeqxIoaWNg4LOxNvTiAlSQXpdc1gnvq0CutYfKq4JKkgva4ZHCI0B01TP/t1\n",
       "Ob8BNuYZlCT1obPAp4sOIpHMZ5M+hyA5OWWy6RbXq/pcIk8VHUCHnio6gA49VXQAHXiq6AA69FTR\n",
       "AXToqaID6FBbZWfeNYOpeEtLTpe3j0CSSsIprCVJJoMeOFF0AB06UXQAHTpRdAAdOFF0AB06UXQA\n",
       "HTpRdABarOp9BpJUhFL1GUhSu+YIQ8+VzQXqQ/b7mjUDabD4m29Pq/3V1n60z0CSZDKQJJkMJEnY\n",
       "gSypAu6G2nCY5TgX5+D8qetPkO2W5LrjS13zuXAmA0mlNwwbZsJcO7nYDhtP5fXmFekQt5lIktqz\n",
       "nzD8dQ7Y2+T5sRbbvQpcJVzIq3SXUzUZSFJ2NcJEm3fG+2eBWxpe82rD+qq43VHCBb5eo4QX8qpK\n",
       "M9Gt8X4e+KjIQCQNtKTJZzNwknCy16UM270OvBSXHyHUEG7JuG1PVKJm8OgX+fs943yFkI0lqShH\n",
       "gH2EI/s56tdtT6brnwNGUssThARypuF9ZuPrSqMSyWDfV/jvv97CZa6/ZKYk9dIIcBy4C7gH2EMo\n",
       "8KcItYT1hII+WT5CKLcaC/6R+LrSqEQykKSSqBFqBWupH5xmGS00RkgaQ4QO6GOUqIkIqtNnIGmA\n",
       "nYPz23O8/O257Ndj/w6wgzA53DyhYP9hw2sak8MCcIBQizhMSAQPrjjYnJgMJJXeKZjK8TyAdj2w\n",
       "zPNbGtaPxFup2UwkSTIZSJJMBpIkTAaSJEwGkiRMBpIkTAaSJDzPQFIF3L2R2vBQjhe3mef8qbOZ\n",
       "Lm7TiwvVjBFmPr0rx8+4jslAUukND7Fh5okcL27zLTaeyvbulbhQzUqUqZnoxaIDkKSMsl6opga8\n",
       "E1/7KmFOIwgT1b1OuDjOHPA2MNpk+/0svoDOPsI1FLquF8lgX8P6BOGiEJOpx1pdGUiSyqSdC9WM\n",
       "EBLAJLAuPpYuD0cJNY31hJlQG8tK4ufsTK1PAD9YYexLyjsZ7CYEn0gK/el4P0rIqsmkT5JUdsmF\n",
       "ai4RLlQzxvVXO4OQNPYDM8BF4DEW9zXMA9+Ny/sJSaHREcLEeBCSy3rgzc7Cby7vZHCAxXN27yIU\n",
       "/MTHdxD+wKF4X7rrgkpSg2YXqtlMOPidI8yAOkH92gbp7YZS63Op5aWu1XKc0JqSJJdc9LoDeYjF\n",
       "O2ADoZYwFG992zkjqW80u1DNO4TLYB5oeHxzan2IlbWAHCMcSN8DPLyC7TMpYjRRsww4z/LTwkpS\n",
       "GSQXqpkmtPO3ulDNFCFJHCY0LR0EDq3g86aANwid0Lk0EUHvk8E89XaxdWS/oAQAf/LH/O3IrdwP\n",
       "bY0LlqRuaedCNWfic/sJtYTDhH6D9Hs1vnez5TOEsvLYiqPOoBfXFD5K/ah/FNhKyJB7CX/ccplu\n",
       "4dZb+MUnb+C3ly7z4Z5xjj7zEL/Y/i02/vQ/82s/k1SYBRrKphKddFaU1whNRM3Ky2R/3R9vd8bb\n",
       "31Ci68bXCH0E6XauSa4fWrqUhYXv86WXH2HPvSP848L3+dLC9/nStj9jT7eDlVQK9h0utpOlW1Fa\n",
       "7a+29mPezURT8ZZ2MN5PI0laSo3QLFXL+4PKdAayJGmxKUI/60zeH2QykCSZDCRJJgNJEk5hLal8\n",
       "LuCIonZcWP4ly6tsMrj0v4wlw0srMEZYUnbNJmxTzirbTPSpm1gz8wRnZ57gbJ4no0jSIKhsMpAk\n",
       "dY/JQJJkMpAkmQwkSZgMJElUeGhpmsNMJakzfVEzcJipJHWmL2oGZXU31IbDdZ45B+dPNUznvdzz\n",
       "WV8jSZ0yGXRBqwJ7GDbMwFmA7bDxVMN2yz2f9TWS1Km+SwZF9B9YYEuqur5LBkn/AcD2b7Hx1Nnu\n",
       "vXfVmmyqFq+k4vRdMui2dIF6GcZm4MdQjRqANRZJWfXFaKI8JQXqDJy9EdYUHY8k5cGaQRONtQHi\n",
       "0bUk9StrBk1YG5A0aKwZRNYGJA0yawaRtQFJg2ygawZlqg1cgrFtxPMjHAYqqccGOhmkh15ug/vy\n",
       "/KyksD8HW4bhNCxOQJ+CNUksW+GL25okqYaEce19TB6SOlWGZqK1wDiwNy5XWlJgb4M9sSAH6oX9\n",
       "bTC8XHNU8trG16QfT79PUruRpJXqRc1gH/BYan0CmAdGgIPAVuC1uD4CnMwzmG41DaWP0lsd4edd\n",
       "22gWS6uahyQtJe+awW5C4Z9IjpSn4/1oanmInBMBdK+juNXRexFa1RiKjktSdeRdMzgA1FLru4Cj\n",
       "cXkW2EGoGRwE3iA0FT3XrQ9fNGndb9gyfJHTg3q0nOqzsH9B0nV63YE8BMyl1jcAxwh9BiPA4U4/\n",
       "4MpHfPZ7P+LvAG7+GHfPPBHmEtq2l/tmLvKzXjXdlE1Se0h3TtsJLSlRxGiiVQ3rb8b76cYXrsTH\n",
       "r/KJf/gjzgP801U+0Y337CeNfRoz8DNYPILJxCANnl4ng3lgfVxeB6HQVvFaDW219iANhl4ng0OE\n",
       "PoJpYBOhiWhZt32Nb3/yBn576TIfPn6I/3vmIX6RZbs/XGV90mT0uyuMPfkZLn5wgTv4YKXhD4ZW\n",
       "tQenwZZK7f54uzPe2pJ3MqgRCv+HgZcIo4W2EvoI5qk3ES3pf17gm6/8GxufP8ptWRMBwOoFVl9r\n",
       "MrqRm54e5/y//5gbTQaS+tCJeEsstLNx3slgiuubFQ7G+670EcDiTuOrV681Q6nLPANa6l99MR3F\n",
       "ok7jhf74m8rI5iOpfw1cwfn7Vax/8jOhFmH/QXc4yZ5UfZVNBittGrr5JlY/PR5qEfYfdMcSI5FM\n",
       "DFJFVDYZ2DRUTunEYPORVB0WospNsw5nawtSOZkMlJtmHc42I0nlNNDJwM7k3rMZSSqngU4GdiYX\n",
       "y1FIUnkMdDJQsRyFJJWHyUClYPORVCyTgUrHaS+k3sv7spdS21pdxjO5drWk7qtEzeDwUf5q9hzr\n",
       "P/oDvys6FknqR5VIBtuvsPrShwwtfMhnv/cjbgZnJ5WkbqpEMthwAx+uXsWVTyxwk1NQDC6Hokr5\n",
       "sUBt4oM13PHkOk9GKxuHokr5MRlE6bORb7jC7U/vDHP1ezJaOZkYpO4yGUTps5G3TbtfqsRzFKTO\n",
       "Weipr9ivIK2MyUB9xVqCtDKedCZJMhlIkmwm0oC4G2rDsMG5jqTmTAbL8AI4/WEYNszA2eSKa2Cf\n",
       "gpRmMliGF8DpX448kuqWSwaTwD3AemAudX8MOJJvaFK+HHkk1bVKBuPAEPAqcLDJ85uACWAWOJlP\n",
       "aOVjk5GkftUqGbwGXFxiuzPxtqkLMWwCRoAxQjX9TBfeMxc2GUnqV62Gli6VCNKyFNz7GtYnCDWP\n",
       "ybg+Rkg+x4Faxs+VuirpP9gGe+72e6gBlLUDeS2wG9gMvAMcIFvC2E0o/B+L62PxfppQGxil3vew\n",
       "AzicMR6pq5z4ToMu60lnOwgJ4BFCH8KOjNsdIPQrJHYBF+LybOp9dhBqBu9mfF8pN+nLbnqpTQ2K\n",
       "rDWDeeo1gYUOPm+IMBopsYHQZPQoITk4Skml4vBTDYosyWACeB94llCYXyAc8a/Uqob16XiTSqfV\n",
       "8NPkjGYwSag/ZK0ZfIHQX7BAqCUk7f0/bPPz5uHatYvXQRiZI1VNckYzmCTUH7Ikg3ngG6n1vYRk\n",
       "sIf2k8EhYCuhJrCJ0Cy0rC1TfPWm1Vx57zKrDrzFr3d/jl+1+blST7RKElIP3B9vd8ZbW7Ikg3XA\n",
       "i8BRwglm84RO5CxNOzVC4f8w8FLcfiuhn2AeeDNLkKdrvPLyu4w8/0tuNxFIUlMn4i3RVv9ulmQw\n",
       "BbxBKNj/EtgfH59tucXibRuryckZzX3TT/DBGu54cp1nJg+qdCfz5TB8+mzBIUlty9pnMAt8p8nj\n",
       "a8l+glrfWnMzNz690zOTB1W6k3kb3Fd0PNJKtDrPIDlLeCkThCYfSVLFtaoZHCF08O4ljCJKmyec\n",
       "hfwq1gokqS8s1Ux0BniuV4FUTXoG04+uXBsuK0mVlKXPYBJ4MC6/SPvDSftSegbTbdNeJEhStWWZ\n",
       "m2gWeIAwr9DHqM82KmkJySgjZ0FVFWQ5oh0CPk84J2CK5TuWJVEfZeQZyqqCLMngXsIZx98kJIbE\n",
       "WmwyktriGcoqqyzJ4BBhBFHSmTxCmHJ6NyYDacUaZkTdMgyn47I1BvVclj6Dkyy+otksYdbSrNc0\n",
       "kNRE+roJt8Gw11BQkbJe3KaZ+a5FIUkqlEMipZLxgjoqgslAKplWF9SR8mQyyJGzmaqbHJaqPJkM\n",
       "cuRspuomh6UqT510IEuS+oTJQJJkMpAkmQwkSZgMJEk4mqjrvOiNpCoyGXSZF72RVEU2E0mSPHLt\n",
       "laT5yDOR1Q3OX6RuMxn0SNJ8lD4T2ekqtFLp+Yu2whe3OU2FOmQyKJDTVagb2pnYzvmN1EpZ+gxG\n",
       "8NrKUu6S+Y28iI4a9aJmsA94LLU+QbgwzghwMPXYbA9iKZxDT9Vr6drAZRgj1iKktLxrBrsJBX1i\n",
       "LN5Px/vReH885zhK4+abWP30Ts4/vZPzn7jBZjrlL10buBHWFB2PyinvZHCAxUf8u4ALcXmW+nWU\n",
       "V+UchyRpCb3uMxgC5lLrSZvlOLAVWNvjeCRJFDOaqFkt4LmeRyFJuqbXyWAernWaroMwrFJS9yUd\n",
       "x3YaK4teJ4NDhOagaWATcCzLRlum+OpNq7ny3mVWHXiLX+/+HL/KM8iieTKauiHpON4G9xUdi3ri\n",
       "/ni7M97akncyqBEK/4eBl4CTcX2cUEt4M8ubnK7xysvvMvL8L7m93xMBeDKapBU5EW+JhXY2zjsZ\n",
       "THH9GY7JuQXTSJJKoSxnIEuSCuRJTyWXPmPZ/gP1Wqu5jJzjqP+YDEoufbEc+w/Ua0knNCyeBK/V\n",
       "46ouk0FJOGeRpCKZDErCy2VKKpKFjqTGPoAtw3AanOV0kDiaSNKimU1vg2FnOR08JgNJks1EkvKR\n",
       "ND2lm50cnlpe1gwk5SJpeko3O6UvteklOMvFZCBJMhlIkuwzqCynuZbUTSaDinKaa0ndZDORJMlk\n",
       "IEmymagvOM21pE6ZDPqA01xL6pTNRJIkk4EkyWQgScJkIEnCZCBJwtFEfcdhpqqiSzC2DfaA01kX\n",
       "xWTQZxxmqir6FKyZiZfX3A4bTxUd0ACymUiSZDKQJJWjmWgtsBUYAo4DF4sNR5IGTy9qBvsa1ieA\n",
       "cWAyrj8OTBMSwe4exCNJapB3MthNKPwTY/F+Ot6PEmoEEGoEm3OOR5LURN7J4AAwm1rfBVyIy7PA\n",
       "jri+lpAU3sk5HklSE73uMxgC5lLrG4D9hKQAIXlIknqsiA7kVQ3rZ+JNklSQXg8tnQfWx+V1EE6O\n",
       "kiQVq9c1g0OEYaTTwCbgWJaNtkzx1ZtWc+W9y6w68Ba/3v05fpVnkGWVnmrioyvXkmpbPljDHU+u\n",
       "c7oKdU+vppK4G2rDoWmZc7BlGE43fmbDa5aMJctr23m/Erg/3u6Mt7bknQxqhML/YeAl4GRcHyfU\n",
       "Et7M8iana7zy8ruMPP9Lbh/URACLp5rYNr2y/92am7nx6Z1OV6Hu6dVUEsOwIfmcbXDfDPys8TPT\n",
       "r1kuliyvbef9SuBEvCUW2tk472QwxfWZ9GC8n0aSVApORyFJMhlIkkwGkiRMBpIkTAaSJEwGkiRM\n",
       "BpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwG\n",
       "kiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJEwGkiRMBpIkTAaSJGB10QEAE8A8MAIc\n",
       "LDgWSRpIRdcMxuL9dLwfLSqQvBx4iz8pOoZOvAe3FR1Dh+4vOoCVehz+vOgYOvFf8KdFx9Ch+4sO\n",
       "oJeKTga7gAtxeRbYUWAsuTh+ljuLjqETcyaDwvwr/EXRMXTigsmgUopOBkPAXGp9Q1GBSNIgKzoZ\n",
       "AKxa7gWzv+Pm3/+BG3oRjCQNomUL4pw9Cxwj9BnUgE3Acw2v+Q2wscdxSVLVnQU+XXQQWY0Ck3F5\n",
       "L/D5AmORpIH18YI//xxwD7AO+CTwz8WGI0lS9+xrWJ8AxqnXgqSl7E0t+91RT5ShA3kpVfwh7CbE\n",
       "najauRST8fZs6rEq/R9qhFhfTD1Wpfh3ADvjctW+O8lBUHo/V2nfjxHirWL8Y8BV4O14eyE+njn+\n",
       "MieDqv0QEgcI50wkqnQuxThwnHAm+EhcT/Z7Ff4P4/E2TYh/lOp9jxZSyw9Rne8OhALnNPBOXK/a\n",
       "vv8GcIQw5L1q3511hPL8LuBB4Du0GX+Zk0GVCtGlVOlcihHq+3k2rj9EmC4keazM/4dp4GtxeT1w\n",
       "kmoVqKPUf7gAa6nOdwdCMtgCzMT1Kv2Ga8DP4/JzhO9OleJPf2+2AmdoM/4yJ4MqFaLLKXoIb1YH\n",
       "qc8PNQa8Rvg/nE+9puz/h7WENvdnUutV+R6tb/JYVb47EOIfp97nUaXf8FZCfKNUM/7EOHAoLrcV\n",
       "f5mTAVTrh9DKPPUf+ToWF6xlNQa8Tjg6gmr9Hy4Sjuz2EM5bgWrE31grgOp9dw4S/oYNhEIJqrHv\n",
       "E+9T/84n/X5Vih9Cf9Ol1Hrm+MucDKr2Q2jlEKG5BULhdKzAWLIaBx6Py1X6P4xRbxd9g1D1r0r8\n",
       "I4QCaDch3lGq9d2ZpF6AnifEXZV9DyG2M3F5HriXasWfGEsttxV/mZNBlX4IaTVClfPhuJ4caYwT\n",
       "/jlvFhFUG3ZTPws8qXJW5f8wTv3LP0ToyKxK/EfibYHQtLVAtb47s4TBBxBqBj+nOvseYIp6rEPA\n",
       "f1Ct+KEea6Jq8S9pkmoM6+oXOwhtjG/H++3x8ar8H9ZSHxr7TOrxqsRfdRPx9vXUY1Xa90ntpqrf\n",
       "nU3Uh5QmqhS/JEmSJEmSJEmSJEmSJEmSJEkqiarNuyEVrUaYCXIorh8pMBapa8o8HYVUNskU39OE\n",
       "uWveJ5zdOUE4+1mqLJOBlF0NOBqXNwNfICSG44Q5naTKMhlI2SUXzIEwkVzSVHSRkBykylpddABS\n",
       "hewnNBMNEWbpTGYYXUX9Uo9SJdmBLLVvkpAMZqnPH3+cUEOQJA2AEcLlQL9cdCCSJEmSJEmSJEmS\n",
       "JEmSJEmSJElSHv4fLRE2FHv+TYwAAAAASUVORK5CYII=\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e0a4850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.yscale('symlog',linthreshy=1) #works for histograms\n",
    "i_s_contents, edges, patches = plt.hist(inclusive_q0_s_expts,normed=False, bins=100, alpha=0.3, color='r')\n",
    "i_b_contents, edges, patches = plt.hist(inclusive_q0_b_expts,normed=False, bins=edges, alpha=0.3, color='orange')\n",
    "i_s_cumsum = np.cumsum(i_s_contents)/np.sum(i_s_contents)\n",
    "i_b_cumsum = np.cumsum(i_b_contents)/np.sum(i_b_contents)\n",
    "i_index = np.searchsorted(i_b_cumsum,1.-testSize)\n",
    "print np.mean(inclusive_q0_s_expts), np.var(inclusive_q0_s_expts)/np.sqrt(inclusive_q0_s_expts.size)\n",
    "plt.xlabel('$q_0$')\n",
    "plt.ylabel('$p( q_0 )$')\n",
    "plt.legend(('s+b', 'b-only'),loc='upper right')\n",
    "plt.title('inclusive: $s_1=%d$ $b_1=%d$ $s_2=%d$ $b_2=%d$' %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))\n",
    "#plt.ylim(1e-4,1)\n",
    "plt.savefig(\"q0_dist_incl_s1_%d_b1_%d_s2_%d_b2_%d.pdf\" %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.444914618 0.348139477587\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEWCAYAAACEz/viAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAGVJJREFUeJzt3U1sXWeZwPF/SmhS0qkdB0HUqpC47UggRqrddDYjNE3t\n",
       "MiAhFtgNYoNmESewJ5Sy6gqadvY0SReITUmasEIjTRKHgMSGlqZIFV20cYCKTjWqHacqTD8mzSye\n",
       "93CPb+61z7V9Pu69/590dM+5vh+Pz7Xf575f5wVJkiRJkiRJkiT1g4/atheBiQ6POwRcTo9ZAp4B\n",
       "RnI/nwTO5X5+qsefS5Jq9BFwP3AHsAc4ku7bm3vMUaIA/3p63ARwlkgcmavAj9JrjBDJopef1+0o\n",
       "8TuPAjO5+2eBqXTfVJfnfo84H1WYBOba7usWY5HYq9Qp9ryqzmM/n0OpNB8RBXTeM2mDKBw7PQbg\n",
       "ZLp/PD2m3X+l27V+DlFreLJAvGU5C7wO/Dh33zit85A9ppMZ4LslxZU3RdSojuTu6xZj0dir0in2\n",
       "dlWcx34+h33tlroD0Lo8D+xL+/uAl4A/dnjcN9L9C2k7y8pvT/+Wbtf6OUQi+NnGwt6QY8C9wHdy\n",
       "900Dy7njZTo3oX0JmC8vtL+bJ5JmXrcYi8ZelU6xt6viPPbzOexrJoP+dJX4VgRRpV7I/Sz7lp9t\n",
       "WXX7XiKJPEY0KZ1lZVPTWj+fB17uEs8c8U9YZjV9LL3HDK1mohFgMfeYJVrnJW8K2Ek9TQndYiwa\n",
       "O1Rzfouo6zxuxjnUGrbWHYDWZYxWArgMHMj9bIHWP8RTrOwEPpE2iALmMnAPcKXgzzuZIZqqLhFV\n",
       "+/n0nvuIRPV0h+eMtMXc7nyH98ziukT0ZZzv8twbbcejxDm5kI6zvpTV4ltvjBvVHjt0Pr+jwANE\n",
       "/C+x8tt6WXGvdR7b4ygzltV0OocqwGTQnx4FXkj788Q3+gmiwIBWk9EDwG+JAuXbwCO51zgBHE7P\n",
       "m1zj56v9w54nOndfJ2oVANeIguORLs+5RqtwL2KWqKVkBXf2DXCZKKQy+SSZmSbaoLP9qwXiW0+M\n",
       "3XSLcazL/e3y5/d76b5HiaaUeaJQzhfCmxV3u/bzuEwU9N3i2MxYNnoOVYDJoD/sJP4hxoiC8SDx\n",
       "jZ10/2PEP+Jcuh0nCpDsW9J5ImEcAc6k+2aJgv480Vy42s8hEsYNWgkHooAeJxLJCPB47vmr6fUb\n",
       "4+W0ZcZSHAvE75kZ5eamrJ20CohZYsRUEev9Vrul7fhUhxiLxt7p/P6cVgE7ycrzspG4O8We134e\n",
       "f5iLt1McG4llM8+hNDDa5xm8QAw1bTdHVNuzx3ydlaM/9hLf3pZyj3k49/y1fn6WmwvSvUTbcdaW\n",
       "f0fbzzZz9FHWV3CkLa6p3PZwh+eNEOdmijgnZcWXxXKKGIU11XZ/pxjXin218wsRf/t969Ut9ky3\n",
       "81hVHOs9h+oj48SHf7DuQLSpxql3KOpamh7fWmaJAnrvWg8ckjjUB462HWcjEbJRLntxpusgmiO+\n",
       "3TW1kGh6fKuZJfoQOtXWhjEO9YFDxB9LZpLWsMBsuBzEP+QcJgVJqkXZ8wyOs7J3/wAxmoN0/zSR\n",
       "HK4Q7d2HSo5HktRB1ZPORokOyswuYnzyFDFe+fmK45EkUc/Q0vZhY1fSVsXlAiRJHVRdM8jGykOM\n",
       "W15c5bGSpIpUXTM4STQHzROdxmtdGAvgL8CdZQYlSQPoTeCuuoPIzBJ9BPk5BNnEldWum57XxGuN\n",
       "PFF3AB08UXcAXTxRdwAdPFF3AB08UXcAHTxRdwAdPFF3AF08UXcAHfRUdpZdMzidtrxsKr19BJLU\n",
       "EF7CWpJkMlini3UH0MHFugPo4mLdAXRwse4AOrhYdwAdXKw7gA4u1h1AFxfrDmAYNLHPQJKarlF9\n",
       "BpLUqyVi6LmKuUpryP5As2YgDRf/53vT7Xz1dB7tM5AkmQwkSSYDSRJ2IEvqC5+bhd27ynv9txbh\n",
       "1fYJsptllrh8/2rrQdfOZCCpD+zeBRfeLO/1H74TXi3rxfuiQ9xmIknqzTFi+OsScKTDzye7PO8U\n",
       "8BGxkFfjlls1GUhScbPEhTb3pNsngTvaHnOq7XhLet5ZYoGvF2ngQl790Uy0Jy2H+TYv8C6Xao5G\n",
       "0vDKmnzuAS4Rk73eKfC83wHPpv1vEzWEOwo+txL9UTP4Fv/Nv/A+29lRdyiShtoZ4CjxzX6J1rrt\n",
       "2eX6l4Dx3P4MkUCutL3OQnpcY/RHMriFfolU0mAbB84D9wIPAIeJAv80UUsYIwr6bP8M0UzUXvCP\n",
       "p8c1hkWsJBU3S9QKRmit515ktNAkkTRGiQ7oczSoiQj6pc9A0pB7azGGf5b5+oU8BUwTF4dbJgr2\n",
       "n7c9pj053ACOE7WI54lE8Oi6Qy2JyUBSH3j1dInzAHr1pTV+fl/b8Zm0NZrNRJIkk4EkyWQgScJk\n",
       "IEnCZCBJwmQgScJkIEnCeQaS+sEnmeV2ylvc5l0WeZsii9tUsVDNJHHl03tLfI+bmAwkNd/t7OLf\n",
       "KW9xm59wJ28XemRfLFSzHk1qJnqm7gAkqaCiC9XMApfTY08R1zSCuFDd74jFcZaA14GJDs8/xsoF\n",
       "dI4SayhsuiqSwdG24xliUYi53H3dVgaSpCbpZaGacSIBzAE703358nCCqGmMEVdCbS8rSe/zSO54\n",
       "BvjZOmNfVdnJ4BARfCYr9OfT7QSRVbOLPklS02UL1bxDLFQzyc2rnUEkjWPABeAa8Bgr+xqWgf9I\n",
       "+8eIpNDuDHFhPIjkMga8vLHwOys7GRxn5TW7DxAFP+n+aeIXHE23jVsXVJLadFqo5h7iy+8SsEh8\n",
       "Cc7WNsg/bzR3vJTb30J354nWlCy5lKLqDuRRVp6AXUQtYTRtA9s5I2lgdFqo5jKxDObxtvvvyR2P\n",
       "sr4WkHPEF+kHgIPreH4hdYwm6pQBl1n7srCS1ATZQjXzRDt/t4VqThNJ4nmiaekEcHId73caeIno\n",
       "hC6liQiqH020TKtdbCdRnerFQ+zhMHs4zCeZ3dzQJGlN+YVqloA9dF+o5kr62bH02I+IfoP8a7W/\n",
       "dqf9K0RZeWq9QRexWjvVZjlL61v/BLCPyJBHiIy6Vqa7wQ5e4eMs8wG3MMkvmOYVfsKd/LG89jNJ\n",
       "tblBe9nUnElndXmRaCLqVF5m5+uhtO1J27/SQxlfdjPRLFH4HyR63y+l4ymillCsynOEH3CZMS7y\n",
       "BaZ5paRYJTXV25wuOClsED1CDK5Zq7y8mLZMT32wZSeD02nLO5Fu55EkrWaWaJYqvVm8STOQJUkr\n",
       "nSb6WS+U/UYmA0mSyUCSZDKQJOElrCU1z1W8GkEvrq79kLWZDCQ1TacLtqlkNhNJkvq4ZvA+k+zh\n",
       "MNAPswclqdH6NxlsY8ffl8ErvmSdJKkDm4kkSSYDSZLJQJKEyUCShMlAkoTJQJJEPw8tzXPOgSRt\n",
       "yGAkA+ccSNKG2EwkSTIZSJIGpZloIHxuFnbviv237oPdr/Wwvwiv2k8iad1MBpUoUtC/NwkXfhH7\n",
       "+78IF35VfH/fV2F/p9c3SUgqxGRQmnwCKFrQr9c/7IALb978miYJScWYDDbVagmgDkWSBJgcJJkM\n",
       "NqxpCaCIfJKAthqEiUEaQoOXDCqZgNaPCWA1+eRgYpCG0eAlg0omoO3etbL5ZZCYGKRh1IRkMALs\n",
       "AyaB48C1esPppr02wJurPnwgmBikYVFFMjgKPJY7ngGWgXHgBJEIXkzH48ClCmJah0GuDRRhYpAG\n",
       "WdnJ4BBR+GfJYDLdzhMF/0TaHwFGaVwiGMbaQBH5xPDwnfBqveFI2rCyk8FxYDZ3fAA4m/YXgGmi\n",
       "ZnACeAk4Ajxdckw9GPbaQBHvTML+6LB3LoPUt6ruMxgFlnLHu4BzwBRRU3i+4ng6sDbQm25zGawx\n",
       "SP2kjg7kLW3HL6fb+aoD6czawOZYUWOwliA1XNXJYBkYS/s7gcVS381Fb2pkv4LUT6pOBieJPoJ5\n",
       "YC/RRLS2p/khH2eZD7iF8ywxzSuFnld4zoFNQ+WyliBV4KG07UlbT8pOBrNE4X8QeJYYLbSP6CNY\n",
       "ptVEtLoj/IDLjHGRLxROBD2xaahcDkuVKnAxbZkbvTy57GRwOm15J9JtjX0E+ZoAWBuokolBaqIm\n",
       "zECuQb4mANYG6mJikJpiSJOBmscOZ6lOJgM1kB3OUtWGJxn89fZJeDAVMPYRNJu1BKlqw5MMbuzY\n",
       "4YihfmQtQarC8CQD/no3n/ralwF45427ea/mcFSQtQSpCsOTDG69sY2Ju2LG829+v81k0I+sJUhl\n",
       "GZ5koAHgUFSpLIOdDD78+Gf49Z3RNHR9aWyNR6uv2HwkbabBTgY3br2V276SLob33GD/rkPN5iNp\n",
       "oywgNQCsJUgbNZzJYMu7Y38fWQSOLhoo1hKk9Ri8ZFCkn+DW61vZd1drLQVHFw0QO5ml9Ri8ZGA/\n",
       "gf7O5iOpqFvqDkCSVD+TgSRpAJuJpI7yHctgH4K0kskAVo4ucmTRgMr3H8DKPoT8yncmCQ2nwUgG\n",
       "G51pnB9d5MiiIZGvKbw3CRd+Eft2NGs4DUafQTaC6LavLMItg5HgVLKspnDhTdi2o+5opLoNRjKQ\n",
       "JG2I36KlFZzBrOFkMpBWcAazhpPJQOrKGcwaHvYZSJLWrBnMAQ8AY8BS7vYccKbc0CRJVemWDKaA\n",
       "UeAUcKLDz/cCM8ACcKmc0KQmWdGxfB/sfi3t25eggdAtGbwIXFvleVfStncTYtgLjAOTwOn0umtz\n",
       "SUtVKt9/sP+LcOFXsW8nswZDtz6D1RJBXpGC+2jb8QxR85hLx5NE8jkPzBZ83/ImmmWXpvjU177M\n",
       "9jfu3rTX1YDKT17LLmkh9Z+ihegIcAi4B7gMHKdYwjhEFP6PpePJdDtP1AYmaPU9TAPPF4ynPF6a\n",
       "QuvmHAX1r6KjiaaJBPBtog9huuDzjhP9CpkDwNW0v5B7nWmiZvDHgq8rNZC1BPWvoslgmVZN4MYG\n",
       "3m+UGI2U2UU0GX0POEzUIiRJFSvSTDQDvA08SRTmV4lv/Ou1pe14Pm2SpJoU7TP4CtFfcIOoJWTt\n",
       "/T/v8f2WibkKADuBxVUeWz/XOZA0JIokg2Xg+7njI0QyOEzvyeAksI+oCewlJq+t7Wl+yMdZ5gNu\n",
       "4TxLTPNKj++7PnYma1O4eI4q8VDa9qStJ0WSwU7gGeAsMcFsmehELtK0M0sU/geBZ9Pz9xH9BMvA\n",
       "y4WiPMIPuMwYF/lCZYlA2hAXz1HlLqYt01P/bpFkcBp4iSjY/xk4lu5f6PqMlc9t/xaUzWi2n0AD\n",
       "rH2SmtRsRfsMFoCnOtw/QvEJapKkhuqWDGaIZpzVvr0Xeczmen/bffx69ycAL0GhPtVtYpr9CqpX\n",
       "t2RwhujgPUKMIspbJmYhn6LqWsGN7Z+Iy08APOdaDOpD3dZI2L3LtRNUp9UK1CvA01UFIkmqT5EZ\n",
       "yHPESKKzwNfLDUeSVIciTS0LwJeI2cfTRHLotMaBpJ60Dz/lzVrD0VArkgxGgfuJOQGniTkC1XqP\n",
       "rXzIxyp/3zxnI2vTOfxUzVEkGTxIzDj+AZEYMiP0PgN5fX77j1O8+/4OPrp2RyXv14mzkSUNsCLJ\n",
       "4CQxgijrTB4nmosOUVUyuO3hJf62cCu8UPQqq5KkHhRJBu1rHC8QVy09tfnhSJLqsJGx+subFoWk\n",
       "nEIT0+6D3a/d/BhpfWx2kRqn24pp2cS0C2/Cp3e7qpo2k7N4pUZz+KmqYTJYD4eZqjIOP1U1bCZa\n",
       "j1uvb+XBuxZ58K5Ftn+4re5wJGmjrBlIfa9bh7NUnMlA6nvdroQqFWczkSTJmsGG2ZksaQBYM9go\n",
       "O5MlDQCTgSTJZCBJss9A0sprHjk0dUiZDKSBsp45B9k1j8ChqcPLZFCW7W/czR2OMlLVnHOg9bHP\n",
       "oCzbP9zmKCNJ/aIpNYNxYC8wX3cg0uDINxmB/QFaTRXJ4CjwWO54hlgYZxw4kbtvoYJYypWfgHb9\n",
       "3bGao9HQyzcZgc1GWk3ZzUSHiII+M5lusxrARLo9X3Ic1chPQNt6vSm1LqkHWW1i/+EYZaRhUXYy\n",
       "OM7Kb/wHgKtpfwGYTvtbSo5DUqGCvtsqaxp0VX97HQWWcsfZH9sUMEbUEK5VHJM0JPLNRvu+CvvT\n",
       "/58rqKmeDuROtYCnK49CGmquoKaVqk4Gy0QNAGAnsFjx+0vaMGcsD6Kqk8FJYB/RgbwXOFfoWb98\n",
       "7lts/dj7fPC32/nzHz7LZz7/pxJjlLQqZyw31ENp25O2npSdDGaJwv8g8CxwKR1PEbWElwu9yv5v\n",
       "/pTFhd289sL9JgJJ6uhi2jI3enly2cngdNrysrkFwzPBzAVwJDWcl6OoggvgSGo4k4EkyWQgSWrO\n",
       "heok9T2HnPYzk4GkTeKQ035mM5EkyWQgSTIZSJKwz0BSV/mV0t66D3a/Fvte5XQQWTOQ1EV+bYNP\n",
       "727tb9tRd2TafCYDSZLJQJJkn0G9tr9xN3d4ATv1s3y/Qr4voVt/Q6+T0ZzIVhVrBnXa/uE2L2Cn\n",
       "/pbvV8j3JXTrb+h1XeVsIptrMpfNZCBJMhlIkuwzqF5+oZvr746t8WhpwGxWX4I2m8mgarde38q+\n",
       "uxYBePG6519DJutLANj/Rbjwq9j3wnZ1s5lIkmQykCSZDCRJmAwkSZgMJEk4mqiZvEyF1MGKYakO\n",
       "Rd1k1gyayMtUSB3kL3HhpSk2m8lAkmQykCQ1o89gBNgHjALngWv1hiNJw6eKmsHRtuMZYAqYS8eP\n",
       "A/NEIjhUQTySpDZlJ4NDROGfmUy38+l2gqgRQNQI7ik5HklSB2Ung+PAQu74AHA17S8A0+l4hEgK\n",
       "l0uOR5LUQdV9BqPAUu54F3CMSAoQyUOSVLE6OpC3tB1fSZskqSZVDy1dBrIFXXYCixW/vySpg6pr\n",
       "BieJYaTzwF7gXKFn/fK5b7H1Y+/zwd9u589/+Cyf+fyfSoyxWfIro3lpCg2sXi81sZkrpn1utjWj\n",
       "udfnd3vuRl5z3R5K25609aTsZDBLFP4HgWeBS+l4iqglvFzoVfZ/86csLuzmtRfuH6pEACtXRvvN\n",
       "77eZDDSY8iugFVn1bDNXTNu9q7f3LvLcjbzmul1MW+ZGL08uOxmcTlveiXQ7jySpEbwchSTJZCBJ\n",
       "MhlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIw\n",
       "GUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZSJIwGUiSMBlIkjAZ\n",
       "SJIwGUiSgK11BwDMAMvAOHCi5lgkaSjVXTOYTLfz6XairkB68v6H2+sO4WaPf6HuCDprYlxNjOl/\n",
       "Pl13BDdrYkxN/OyguXEVV3cyOABcTfsLwHSNsRT3wf/dVncIN/v1P9UdQWdNjKuJMS01sOBtYkxN\n",
       "/OyguXEVV3cyGAWWcse76gpEkoZZ3ckAYMuaj/jfq5/g+vu3VhCLJA2ltQvicj0JnCP6DGaBvcDT\n",
       "bY/5C3BnxXFJUr97E7ir7iCKmgDm0v4R4P4aY5GkofWxmt//LeABYCdwG/Cf9YYjSdLqjrYdzwBT\n",
       "tGo2ar4juX0/PymnCR3Iq2nKP+yhFEumKfMj5tL2ZO6+us/ZbHr/Z3L31R0TxLDlR9J+Ez6/7MtF\n",
       "/pzUfZ4mUwxNi+kj4PW0/bghcXV6/7pjOsIGPr8mJ4Mm/MNmjhPzIDJNmB8xBZwnZm2Pp+PsHNV1\n",
       "zqbSNp9imqA5n+ON3P43qP/zmwNeAy6n4yacp+8DZ4gh30357HYS5dS9wKPAUw2Ia4L4u5lPtxPU\n",
       "/7+X/Q2fAe4hBuP0dJ6anAyaUOB204T5EeO0zslCOv4GcWmP7L6qz9k88J20PwZcohkF7wStfwiA\n",
       "Eer//OaA+4AL6bjuv/dZ4IW0/zTx2dUdE6z83PYBV2hGXFnNbpxm/J1P0/picTkdH6CH8qDJyaAJ\n",
       "Be5q6h6We4LWtZwmgReJc7aYe0wd52yEqK7+KHdc9+c41uG+uj+/MaIWlfVj1P33vi+950SDYsqb\n",
       "Ak6m/brjukQkpaVcHHXHtJh7z51E7aCnmJqcDKD+f9hulmkVMDtZWQBXbRL4HfEHCvWfs2vEN8vD\n",
       "RFUV6o2pvVYAzfj8ThBx7SIKOqj/s3ub1t9R1kdWd0yZR4B3csd1xjVK9F/MEZ9jE/7OTxMJAKK2\n",
       "0vPfdBOuWtpNE/5huzlJfJOaJ/4QztUYyxTweNqv+5xNEm3zl4CXiKaHumMaT9uuFMcE9X9+c8Q3\n",
       "tjPE+Rin/vO0SHzbJcXyYANiypvM7dcd1xxwjEhOyzTj7/wK8Xc9kWJZoPU3XyimJtcMThL/JFB/\n",
       "gTtLFB4H03H27WmKOPEv1xEUMcopm7GdVaPrPGdTtP74Rom2y7pjOpO2G0STVZasoL7Pb4Ho/If4\n",
       "h32B+s/T6dz7jwK/bUBMmfG24ybEldVS5om/obpjmiDKqEvE53emATFtqjnqH5LYVNPEt8vX0+3D\n",
       "6f46z9kIreGuP8rd7+d4s5m0fTd3X93naY6IqWmf3V5aQ0ozdcfVaRhn3TFlf1P5KznUHZMkSZIk\n",
       "SZIkSZIkSZIkSZIkSepDTbnuiNQvZomrU46m4zM1xiJtmiZfjkJqmuyy4fPEtXveJmZ3zhCzr6W+\n",
       "ZTKQipsFzqb9e4CvEInhPHGdKKlvmQyk4rIFeyAueJc1FV2jdflgqS81+RLWUtMcI5qJRokrj2ZX\n",
       "Qt1Ca5UpqS/ZgSz1bo5IBgu0rrN/nqghSJKGwDixxOjX6w5EkiRJkiRJkiRJkiRJkiRJkiRJKsP/\n",
       "AwVCJS+piiz4AAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dfd70d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.yscale('symlog',linthreshy=1) #works for histograms\n",
    "s_contents, edges, patches = plt.hist(q0_s_expts, bins=100, alpha=0.3, color='b')\n",
    "b_contents, edges, patches = plt.hist(q0_b_expts, bins=edges, alpha=0.3, color='g')\n",
    "s_cumsum = np.cumsum(s_contents)/np.sum(s_contents)\n",
    "b_cumsum = np.cumsum(b_contents)/np.sum(b_contents)\n",
    "index = np.searchsorted(b_cumsum,1.-testSize)\n",
    "print np.mean(q0_s_expts), np.var(q0_s_expts)/np.sqrt(q0_s_expts.size)\n",
    "plt.xlabel('$q_0$')\n",
    "plt.ylabel('$p( q_0 )$')\n",
    "plt.legend(('s+b', 'b-only'),loc='upper right')\n",
    "plt.title('DGSS: $s_1=%d$ $b_1=%d$ $s_2=%d$ $b_2=%d$' %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))\n",
    "plt.savefig(\"q0_dist_DGSS_s1_%d_b1_%d_s2_%d_b2_%d.pdf\" %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22.629480736 0.381990228924\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAGLFJREFUeJzt3W2MHPV9wPHvEcKDTeLzughEFMMdkCclET6cpIpAAp+h\n",
       "SVVeJHdxIkVKpYozad4R1RCQUIiIAIdKl7xq7ONVpEjY2HlRpYXaPhck+qYBTCMk0gTOwS0OIO58\n",
       "JhYNjuH64j8zOzfePc/e3ezM7nw/0mp31ju7vx3fzG//zyBJkiRJkiRJkqR+9Aow1uL5QWBuhe89\n",
       "Ary8wveQ1MPOLzsAdWwhutXBTuCnwAlgFNgfPT8ePTcIzAPTLfa9C9gK3Fp8mIwA1wNTqefaxZgn\n",
       "9m5qFXtat45jLx9DqVQvA1+NHt9FKB3MpR7HthJKFe8DB4B1qX/bHr32feBZYCh6PltSOAg8vLrh\n",
       "d+QAIZ5/Sj03TEgU6de0Mgb8Q0FxpY0Ce4EdqefaxZg39m5pFXtWN45jLx/DvnNe2QFo2UaA7wFb\n",
       "CBf1r9MsQQwSTrLbo8czNH+BDRJOqpuB9dG/3dHmMx4GHisg9rx2AdcAf596bivh12FsHtjUYt9b\n",
       "6c4vyGlC8kxrF2Pe2LulVexZ3TiOvXwM+45JoTcNEC7ku4AXgJOEksJA9O/bgEPAvwNvA98mnEwQ\n",
       "TqL1wH9F+8XF8Famo/dvZYJwMo6u4HucSyP6jDGa7SjrgNnUa+YIvx6zRgnfc6zgGFtpF2Pe2KE7\n",
       "xzePso7jahxDLYNtCr1rPaHqJ3Y09fhqQr1rujop3Q5xL+EEj/99psPPHiMkkiOEIv804WTdTCjB\n",
       "PNJin3WEZNXOIRZ/B2iWbo4QvuuhNvtm21ji0tHhaPtAtP9qx7caWrUPtTq+g4Q69xHgeRb/ei8q\n",
       "9nMdx2wcRcaylLq0sXWFSaF3zRAu/rH0r6WXgX0sPjnjIvY4ISFsIZQithNO8E4cIjQCvwzcHT13\n",
       "Morpljb7nKR9Y2Yr44RqsfgCHv8inGdxyabB2UltK6H6LH58ooD4ltIuxkab57PSx/eu6LmvEapY\n",
       "pgkX5/TFeDVjT8sex3nC31S7OFYzlpUeQy2TSaE3LQB7CCfkHsIvr52pf3882h4FnoseXwX8FaGE\n",
       "MUdICIOEpNDupBqJPutI6rkhwsX524RfhffQ7BW0lE5/Qb4S3WKNKI4ZFn/XQc6u4orbSiAkl4cK\n",
       "iC9tILO9t0WMeWNvdXx/QfNCO8Li47Lasadlj+ODqXhbxbGSWFbzGEq1ku59NEG4wM8SqhzSda2j\n",
       "0WvfB/4N+HDq3w5E+/2KZjXSFkJp4neZ12UvqEPRPnFd/4cz/7aavZXitoQdUXyx0dRtS4v91hGO\n",
       "zSjNY1VEfHEsewnHeDTzfKsYzxX7UscXQvzZ55arXeyxdsexW3Es9xiqxw3T7Cmj3jZMuV1Yz6Xq\n",
       "8Z3LOOFCPXSuF9YkDvWonZntuBfDRLQ9xOI+9OpdE4Rfe1W9WFQ9vqWME0p+rUpvdYxDPWo7iwdD\n",
       "jdDsWhh3uYNwkk5gcpCkUhU9TmE3ixsxtxF6ghA9v5WQJI4SurptLzgeSdISuj14LTtp2wZCX+dR\n",
       "Qt/nx7scjyQppYwuqdmuZ0ejm5NaSVLJul1SmCf0N4fQB3p2iddKkrqs2yWFPYRqomlC4/K5JuMC\n",
       "eA24osigJKkPHQc+UnYQWfH8O+kxCPFgmImWe5zNeU2a7i87gAq5v+wAKuT+sgOokPvLDqBClnXt\n",
       "LLqksC+6pcXD9W1DkKSKcepsSVLCpNBbnio7gAp5quwAKuSpsgOokKfKDkDFs01BkjpXyTYFSerU\n",
       "HKHLuvI5QbOrfy1YUpDqxXO+M+2O17KOo20KkqSESUGSlDApSJISNjRLqryNGxlvNNhQ1PvPzTF7\n",
       "7NhZA21Xyzhh2YCl1q6uDJOCpMprNNgwOcnxot7/zju54tixot69txrOrT6SpM7sInSbnQN2tPj3\n",
       "kTb77QXeJywoVtklYU0KkpTfOGFCz6ui+4eBD2deszezPRDtd4Cw0NizVHhBsd6oPhoa+gIAf/jD\n",
       "Mf70pz+UHI2k+oqrgq4GjhAGjb2dY7/ngEejx98mlBg+nHPfruqNksJPfvIpvvvdz7NxY2WLXJJq\n",
       "YT+wk/BLf47muvLxMgFzwHDq8RghkRzNvM9M9LrK6Y2kcNttb3HNNZXLqJJqZxg4BFwDXA/cQbjw\n",
       "7yOUGhqEC378eD+h+iibAIaj11VObyQFSaqGcUIpYR3N9ebz9C4aISSPQUJD9UEqWHUEvdKmIKnW\n",
       "5uaYvfPO4pblnZvLvV78j4CthEno5gkX+F9kXpNNEgvAbkKp4nFCQvjasoMVCyws3MYTT/wtH/vY\n",
       "F8sORlLheqpffwU4IZ4kqRgmBUlSwqQgSUqYFCRJCZOCJClhUpAkJUwKkqSEg9ckVd8nPznO5ZcX\n",
       "tsgOr78+y0sv5VlkpxsL5owQZlq9psDPaMukIKn6Lr98A4cPF7bIDlu2XMFLL+V5Zd8PrKtS9dFP\n",
       "yw5AknLKu2DOOPBK9Nq9hDmTIEyI9xxhkZ454GVgU4v9d7F4IZ+dhDUcCtONpLAzsz1GWJxiIvVc\n",
       "u5WKJKlKOlkwZ5iQCCaA9dFz6evhJkLJo0GYeTV7rST6nFtS22PAY8uMPZeik8J2wpeIxRf/6eh+\n",
       "EyHLxpNLSVLVxQvmvE1YMGeEs1dfg5A8dgGHgZPA3Sxui5gH/jF6vIuQHLL2Eybgg5BkGsALKwt/\n",
       "aUUnhd0snjN8GyEBED2/lfBFB6N7F9GRVHWtFsy5mvAjeA6YJfwYjtdWSO83mNqeSz0eoL1DhNqV\n",
       "OMkUqtsNzYMsPhAbCKWGwejW9404knpeqwVzXiEsz7k78/zVqe1BllcjcpDwg/p64PZl7N+RMnof\n",
       "tcqI88Ct3Q5EkpYhXjBnmtAO0G7BnH2EZPE4ocppCtizjM/bBzxPaKwutOoIup8U5mnWm62H3Atb\n",
       "BAMDN3PzzZ8BOulXLEmrpZMFc45G/7aLUGp4nNCukH6v7Hu3enyUcK08uOyoO9DtpLAH2EzIsEPk\n",
       "/ZKXXfYga9bMc+rUAF/5yr/w0EMvdtCvWFKve/31WbZsKWzlNV5/Pe8P1P3RLa92r58Brk1tPw98\n",
       "LvX42szr41XelnJTdLsqulXSOKENIV0PNsHZXVKX0lx57Ytf3MnCwm0sLNzGzTffsdrBSqoE2xYX\n",
       "u4Wla1VWdeW1oksK+6Jb2lR0P40kaSnjhOqq8W59YJVGNEuSFttHaIc93K0PNClIkhImBUlSwqQg\n",
       "SUo4dbakqjmBPZA6ceLcL8nPpCCpalpNDKcusfpIkpQwKUiSEr1bffT22yPJqGbnQZKkVdG7JYUP\n",
       "fWgthw8f5/Dh44Uu6C1JNdK7SUGStOpMCpKkhElBkpQwKUiSEiYFSVLCpCBJSpgUJEmJ3h28ltaD\n",
       "A9k2bmS80WADwNwcs8eOnbVCXa59l7O/JLXTH0khHsgGsGXLFbz0UskBnVujwYbJSY4D3HEHf3Pd\n",
       "dfkTRHrf5ewvSe30R1LocWvWsDa+yN95J1ccOxaez5Qmrm00+B3A6dOMQDMppPc3QUhaCZNCAdpV\n",
       "DaWfz17YY++8w8h113FH/JrJSX4JcOed3Dg5ydPx43afbYKQtBImhQKkq3fSv/wzz7e8sGdKDW0v\n",
       "/nmYICR1yqSwStqVArK//GlROuiGJRJEUi1lspBkUuhQu6qhdqWA1fzlv1qyMcXVUpYmJJkUOtSu\n",
       "aqgfWN0kycFrailOEJOTHE+PiZDU36pQUlgHbAZGgN3AyXLDya8q7QVFS39PSw1Sf+tGUtgJ3J3a\n",
       "HgPmgWFgipAQno22h4EjK/q0Lo5urmJ7QRGsVpLqo+iksJ2QBOKkMBLdTxMSwKbo8TpgkJUmBChk\n",
       "dHOe8QV1YYKQ+lvRSWE3MJ7a3gYciB7PAFsJJYUp4HlgB/BIwTHlkk0E6UFk5UZWHe1GYkvqXd1u\n",
       "UxgE5lLbG4CDwCih5PB4l+NZxESwfJl2B8c+SD2qjIbmgcz2C9H9dLcDycoz4littRv7YAlC6i3d\n",
       "7pI6DzSix+uB2S5/viRpCd0uKewhtCFMA0OEqqNzu+yyB1mzZp5Tpwa45555HnroxQJjVEFWsoaE\n",
       "pHO6KbpdFd2WpeikME5IArcDjxJ6F20mtCHM06w6Wtobb9zLk082eOCBT5kQessSs75arSStrqei\n",
       "W2xhOW9SdFLYF93SpqL74tsQenBFtn7TbiyHA+KkaqrCiObi5Biz4BiEcjjeQaqm/k4KOdjjqHwm\n",
       "CKk6nBBPleJEfFK5allSsMqoN9juIHVfLUsKcZXR5CTHL7iAtWXHo9YsNUjdV8uSgnqPpQapO2qT\n",
       "FC4+9ebIx2uw9kG/cvI9qTtqkxR478wnLv7Cp88AvPkfr30CTvyy7JAkqWpqkxQW1n7owkvv+7tZ\n",
       "gFfGH7kQTpQdkiRVTm2SgvqH7QtScfo6KXzgz+9s/IupH34J4O33zjTO9Xr1Bge7ScXp66TwwfPe\n",
       "v+ALf71hFmD6Zwt9/V3rygZoaXV5oVTfsFpJWjmTgvrGEtVKLg8q5VTLEc3qf5nR0Jc7MlrKp+9K\n",
       "CjYuS9Ly9V1SyNO4/M675zXufezTXwJ488xrH3XMgiQFfZcU8lhYe8n5l953hwPZasjGaGlptimo\n",
       "Vpx5VVqaSUGSlDApSJISJgVJUqKWDc0S2OgstWJJQbVlo7N0ttqXFByzIDir1OC0GKqtcyWFCeB6\n",
       "oAHMpe4PAvuLDS2/lYxidsyC4KzZVm+cnORpcGpu1U+7pDAKDAJ7gakW/z4EjAEzwJFiQsvPKbJV\n",
       "FKfmVt20u4A+C5xcYr+j0W1oFWIYAoaBEWBf9L6SpBK0a2heKiGk5bmA78xsjxFKIhPR9gghCR0C\n",
       "xnN+riSpAHmrWtYB24GrgVeA3eRLHNsJSeDuaHskup8mlA420Wyb2Ao8njMeqevswqo6yJsUttJM\n",
       "BIPRdp6G5t0s/vW/DTgQPZ6J3udIdH8I+H3OeAphTyQtxfYF1UHepDBPs2SwsILPGyT0XoptIFQl\n",
       "3UVIEqX2arInkqS6y5MUxoC3gIcJF/UThBLAcg1ktqejmySpZHlLCl8mtCcsEEoNcXvALzr8vHnC\n",
       "WAeA9cBsh/tLleBgN/WrPElhHvheansHISncQedJYQ+wmVAyGCJUF53bZZc9yJo185w6NcA998zz\n",
       "0EMvdvi50qpqN9jNtgaV6KbodlV0W5Y8SWE98FNCA/ERQpKYIl+VzzghCdwOPBrtv5nQjjAPvJAr\n",
       "yjfeuJcnn2zwwAOfihOCazGr6jZuZDyeU8kShLrgqegWW1b7b56ksA94nnCB/zywK3p+Jue+2RMh\n",
       "HiG9onYERzGritLVSqdPMzI5yS/BEoR6R96L6QzwoxbPryP/QDep72WrlcqOR+pUuxHN8ajjpYwR\n",
       "qoIkSX2iXUlhP6EheAeh11HaPGFU814sJUhSX1mq+ugo8Ei3ApEklS/PymsThJ5HB4CvFhuOJKlM\n",
       "eRqaZ4Bbac55NEHrNRb6ivMgSaqjPCWFQeA6QlvCPvJ1Re15YR6k78xeet93Zs9cMnhh2fFIUjfk\n",
       "KSl8jjCC+V5Cgoito/MRzVItOS2GekWepLCHUEqIG52HCdVI2zEpSLk4LYZ6RZ6kkF2DeYYwS+re\n",
       "1Q+nmmxfUFEsQahqVjI9xPyqRVFxrrOgoliCUNX01JxBF7z7x4+vn/rhJeAkeJJUhJ5KCh88/72L\n",
       "nARPkoqTp0uqJKkmTAqSpIRVMFIFZXol2RNJXWNSkCoo0ysp6Ynkam4qmkmhQ45ZUJkaDTa0ShbS\n",
       "ajEpdMgxC5L6mUlB6jNWMWkl7H0k9Zm4imlykuNxcpDyMilIkhImBUlSoifaFM5/4p+HFn7z3+t4\n",
       "f6FSi93YE0llciyDitATSeEvL3nxmvmL/rfxe96/uOxY0uyJpG5IX/xPn2YEQpfUdmMZpJXoieqj\n",
       "xvD6+YvWXfhO2XFIZYgv/pOTHL/gAtaWHY/6W1WSwjAwWnYQklR33UgKOzPbY4QEMJF5bhBJUqmK\n",
       "TgrbCRf82Eh0Px3db4ruDxUchyQph6KTwm7Cms6xbTRbY2eArdHjgYLjkCTl0O3eR4PAXGo7Hm05\n",
       "CjQIJYaTXY5J6nnteijl4bQYSiujS2qrUsEjXY9C6iOZ7qk3drKvM68qrdtJYZ5QIgBYD8x2+fOl\n",
       "2rJEoDy6nRT2AJsJDc1DwME8O03f8P1vnbf2onfffevU2ld//syVV37zhleLDFLqR5YI+t5N0e2q\n",
       "6LYsRSeFcUISuB14FDgSbY8SSg0v5HmT0Wd+8LPjT//28t/8+InPmhAkqaWnoltsYTlvUnRS2Bfd\n",
       "0qai+2n6iPMgSeoHPTH3US9wHiRJ/cCkIPWxlXRVVT1VZe4jSQVwMj11yqQgSUqYFCRJCdsUCmBP\n",
       "JPUDB7vVk0mhAPZEUj9wsFs9WX0kSUpYUpBqqNOuqunXW5XU3ywpSDXUaVfV9Ovjdgb1J5OCJClh\n",
       "UpAkJWxTkNQR2xf6myUFSR2xfaG/mRQkSQmTgiQpYZtCwZzyQv1sJe0LTqNRTSaFgjnlhfpZ3L4A\n",
       "nU+F4TQa1WT1kSQpYVKQJCVMCpKkhElBkpQwKUiSEvY+qoA3z6z/6L2PfeRLAEdPz22697HGyfC8\n",
       "XVjVOzLdU69tNPhd9Njupj3EpFABZy4ZvPDS+74zC3DmG7sutgurelGme+qNk5M8HT22u2kPsfpI\n",
       "kpQwKUiSElWoPloHbAYGgUPAyXLDkaT66kZJYWdmewwYBSai7XuAaUJC2N6FeCRJbRSdFLYTkkBs\n",
       "JLqfju43EUoIEEoIVxccjyRpCUUnhd3ATGp7G83uNDPA1mh7HSE5vFJwPJKkJXS7TWEQmEttbwB2\n",
       "EZIDhCQiSSpJGQ3NA5nto9FNklSybndJnQca0eP1wGyXP1+StIRulxT2ELqfTgNDwME8O03f8P1v\n",
       "nbf2onfffevU2ld//syVV37zhleLDLIo6VXY0tNZ/N97f24svefS0tNkODWGViI9VcXp04xAGKG8\n",
       "WtKrraXfP88Kbul9s6/LrOJWyhQbFVhJ7qbodlV0W5aik8I4IQncDjwKHIm2RwmlhhfyvMnoMz/4\n",
       "2fGnf3v5b378xGd7NSHA4lXY0tNZvP+NXSv6f0hPk+HUGFqJ7FQVq/3+mdXWkvfPs4Jbet/s67Lv\n",
       "W8YUGxVYSe6p6BZbWM6bFJ0U9kW3tKnofhpJUqU4zYUkKWFSkCQlTAqSpIRJQZKUMClIkhImBUlS\n",
       "wqQgSUqYFCRJCZOCJClhUpAkJUwKkqSESUGSlDApSJISJgVJUsKkIElKmBQkSQmTgiQpYVKQJCVM\n",
       "CpKkhElBkpQwKUiSEiYFSVLCpCBJSpgUJEkJk4IkKWFSkCQlTAqSpIRJQZKUMClIkhImBUlSwqQg\n",
       "SUqYFCRJifPLDgAYA+aBYWCq5FgkqdbKLimMRPfT0f2msgLpBVNTfLrsGKrCY9HksWjyWKxc2Ulh\n",
       "G3AiejwDbC0xlsr79a/5TNkxVIXHoslj0eSxWLmyk8IgMJfa3lBWIJKk8pMCwMC5XvDH/5lbc/qP\n",
       "f7qgG8FIUp2d84JcsIeBg4Q2hXFgCHgk85rXgCu6HJck9brjwEfKDqJTm4CJ6PEO4LoSY5Gk2vtA\n",
       "yZ//OnA9sB64GPjXcsORJEm9ZEfq8RgwSrO0KdXRzsx2q/Mi97lShYbmpdT9pJ+Ibg+nnqvzMdkK\n",
       "3BI9rvMYlxHC38GyTvo+s6ILYB/YTvi+sVbnRUfnSpWTQp1Pegh/1IcIo7yHo+34GNT1mCykHn+d\n",
       "+o5x+R6wn9Clu+OTvo9sIvzfT0f3m6jfObKb8N1jrcZ+bSPMGpF+rq0qJ4W6D2wbpvmdZ6Ltr9PB\n",
       "f26f2UTzRAdYRz3HuIwDv4oePwIcod7nSlx1Mkw4FnX+sQCtx351NB6sykmh7gPbpmjOBTUCPEs4\n",
       "JrOp19TpmDRaPFd2l+oybCb8v2+i2b5S13PlCHCU8N3j71/XY5G2ovOiykkB6nnSZ40AzxFOAKjn\n",
       "McmWEiCUmOJEsZ7FybLfvUXz7yGuT67j38Ug8DKh7WCKMM4J6nksYunzIv4R2dG5UoVZUtup80mf\n",
       "NgrcEz2u6zEZjm4bCN9/E7CH8Kt5mnAxOFhadN01S/h1DOHv4XPU9+9iAtgFvE04BuPU91jE0ufF\n",
       "MOG8GKCDc6XKJYU9hC8F9Trp07bTHOE9Sn2Pyf7otkBoS1ig+Ut5lHAheKGc0LpuH82/gUHgP6nv\n",
       "3wWEhADhgjdP/Y7FOOGCf3u03eq86KtzZYL6dC3L2kqoG305ut8SPV/nY6JgglBt9FDmuTr+Xezg\n",
       "7O65dT0WkiRJkiRJkiRJkiRJkiRJkiRJlVLnOUKk5RgnzMI5GG3vLzEWadVVeZoLqWri6cynCXMO\n",
       "vUUYOTtGmH5D6nkmBSm/ceBA9Phq4MuEBHGIME+V1PNMClJ+DZqTiy3QrEI6SUgSUs+r8tTZUtXs\n",
       "IlQfDRJW9YpnbR0AXikxLmnV2NAsdW6CkBRmaK6PfIhQYpAk1cgwYWnUr5YdiCRJkiRJkiRJkiRJ\n",
       "kiRJkiRJklRl/w9ktVfppUT7WgAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f6ce1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.yscale('symlog',linthreshy=1) #works for histograms\n",
    "truth_s_contents, edges, patches = plt.hist(truth_q0_s_expts, bins=100, alpha=0.3, color='y')\n",
    "truth_b_contents, edges, patches = plt.hist(truth_q0_b_expts, bins=edges, alpha=0.3, color='cyan')\n",
    "truth_s_cumsum = np.cumsum(truth_s_contents)/np.sum(truth_s_contents)\n",
    "truth_b_cumsum = np.cumsum(truth_b_contents)/np.sum(truth_b_contents)\n",
    "truth_index = np.searchsorted(truth_b_cumsum,1.-testSize)\n",
    "print np.mean(truth_q0_s_expts), np.var(truth_q0_s_expts)/np.sqrt(truth_q0_s_expts.size)\n",
    "plt.xlabel('$q_0$')\n",
    "plt.ylabel('$p( q_0 )$')\n",
    "plt.legend(('s+b', 'b-only','other'),loc='upper right')\n",
    "plt.title('Ideal: $s_1=%d$ $b_1=%d$ $s_2=%d$ $b_2=%d$' %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))\n",
    "plt.savefig(\"tq0_dist_truth_s1_%d_b1_%d_s2_%d_b2_%d.pdf\" %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(18, 13, 12)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_index, index, truth_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcHVWZ97+dsCPpTgIIopJuQAEXSEdkEJWYBFw+6igJ\n",
       "uM84SgfEZVQggMxo5FUh4PKqM5IA7msScGBGRkmIoqKvkEDAXekERAmIpJOII+JI+v3jqUPVra67\n",
       "VN+qe+qe+/t+PvXpe6tO1XnO6brPc9bnASGEEEIIIYQQQgghhBBCCCGEEKIOfb4FEEGzDFgObAPm\n",
       "A1dH5xdF5waA7cC6jHuXAAuAk0qWcRiYA1yROFdPvlbk7iRZsifpVB3Wk6Vb6lEI0QHWAKPAZYlz\n",
       "Q5iRSKbJYiFwdklyOeYDq4BzEufqydeq3J0iS/Y0najDerJ0Sz2KBFN8CyCCZgVwKPCWxLkFWKvQ\n",
       "sR2YnXHvSZTfclwHrE2dqydfq3J3iizZ03SiDuvJ0i31KBLIIPQuI9gPcX6JecyI8lgYHQD9wNZE\n",
       "mjGs1ZhmPjA9uq9MGdPUk69VuaEzddsKvuoQiqlH0WF28S2A8MJCbPx2I9bNX4f9UJ+FjQVfmnFP\n",
       "P3Bqg2feANyVOufGkzcCG6I0WYynvg8Am4FvR9/XRPcXLV+7pOWG7LodwMbXh4HbmNhqL0P2ZnWY\n",
       "JUdZsjQjqx6FB2QQepMbsAnfUeDc6NwOTIGcWOeeHdSfvMxiETBIrLxdS3A7pqwcM6J8kyzAxqTd\n",
       "520lyFePevLNqHM+TbJul0TnTsGGVNZhijmtiIuSPUm6Drdjir6RHEXK0m49Cg/IIPQeg5hiPgNr\n",
       "DZ5PvPqnEXlbjpuiwzEDazVvxhSmYwC4PfWs6cRKYhFwUQnyOdIr7VZlyNeq3Fl1+3ViBTtMbZ0U\n",
       "LXuSdB1+KCFvPTmKlKWdehSe0LLT3sMpLddiXwv8MXHtdOC8gvJy8wZDwK3EwxfJ8ezxxHmHU0qb\n",
       "o89fL0m++dHz+oFLiFvM9eRrJnejugW4GFPMf6R96snuqFeHRcvRSJbJ1qPoUfqxl2Nh9Fn4ZQhT\n",
       "FlWl6vI1YhH2jg9KDtGrLEt9d6sdRqLv7sfdT+P11KIzjADXU11lUXX56rEIm1NYQ2vDX6HLIXqQ\n",
       "xdjL5xgmHkJwy/KSG1SSn4UQQnSYMvchXE7t6oFTsdUiROfd6pF+bGKp3iSXEEKIDtDJVUYD2GSb\n",
       "Yya2k3VB9P3yDsoihBAiRaeXnaZXNd1F8ZtchBBCTIJOuq7Yjq1FB1sjvbVBWiGEEB2mkz2Eldi2\n",
       "+XXYKpFmjrnAJqUPKVMoIYQIkE2YY8nKsAibMzgtcW6E2mWnzeiUj5OlHbivWdp61/OcT59r9r0M\n",
       "2smj1XtbSVcvTZ7z6XPNvpdBO3m0em+9dH3Yoo9PAHvVSbMe23R4J3AvtlDkUeB5wBHAXOBk4D8x\n",
       "Vx63A/cAvwMeAh7ENsj9Ffu9F3nsBP4EPAD8BvgFsAW4CWuYXg9cB1wLXAV8DfgS8DngSsxt+yeB\n",
       "jwE/wJbqXgi8F3hPVJ53A+8Azow+vxz4Cubh9wPASzLqbGmduqx3rZVz6e+T0p1l9hCuio4kbgt/\n",
       "1QJi3NiB+5qlrXc9z/n0uWZ5lkE7ebZ6byvp6qXJcz59rpV8i6adPLPunQocgCnuP2ek+zjwfGB/\n",
       "YF/gEUxxfxZzPTE1unZQdIwCf8EWjewbXdsd+F5G3i/LOPe4xOe/YY3IsUi+sdT3hyKZ/ww83MJn\n",
       "Z2SSzCV/ne4OvBIzYk8kLnvy84GYYbs3kuEhYgOZplH+WddaOdfomS1TddcV41Rfxm5iKZ1p1fYK\n",
       "S6l+fb4FU4JPwhTYAZiCXYi1epP0AS+I0u4DTAMeT6z0DgKegBmFZjyMKcT7sPnCsSbH6cAF+PV8\n",
       "2gc8BTg+OoaxOuvHynEvZhTuzfi8BTOMVWFSulPO7XqLG30LEBg3esjzSOBpxAreHRcC38pIfw/w\n",
       "H5jy+h2muGZi48v/FP1NHtNakOFBspVi8vt28in3NTnTF8Hu2LymMwDPAf4HM5Q/wJbF34MNOe3s\n",
       "sGxeqHrrWz0E0Uv0YUMuh2KKO2tJ9nmYEnMK/rfR358RRyKbgrXm08reHfXmA8CGOu7ENo865V7l\n",
       "lnAe9iVW/scDRwO/JDYAP8DKGgKT0p1VV7YyCCJ0Xga8ATgMU9aPYAp5GXBNC/fvCjwDOBZ4NmYs\n",
       "DsNav/XYio39b4r+Jo8HCSdgzZOxja/OABwI/IhY+d+MTTqHiAyCEBViCqaYh4GjsNU1X8tIdyzm\n",
       "xfVOTEFnTUI6+oCDo3ucAZgD7JGR9vdMVPbOCDTKo5uZitXJS6PjIGx5+02YAfgJtgKqF5BBEKIC\n",
       "PB9bbng01tq+DTMG64D/l/NZA8Ax1BqA/TPS3QncgrV4bwF+jg399AL9wAsxA/BibPL3G9FxM71j\n",
       "ANLIIAjRAXbFJnWnA9/JuP4kbKXKRmp9dzWjD5swfj7wd5gBeGpGuq2YonPK/5ac+YTAU4h7AccA\n",
       "38cMwHXYfgMhgyBEKUwDXo0N/czBjMHdwH8Rx6OeDFOBZ2IG4ARsI9e+qTSPYIbFKf+bscneUMb4\n",
       "W2U34LnERmBvYgOwDlsZJGrRslMhSmAXbDnircAXgTuY3ETkrphRcQbguUyMErgF+C7wQ0z534Ft\n",
       "rupFpmB19Dpsz8RmzAi/ChuC6zWjKNA/XZSHm6B9HeaiYCONl2PmZXes1X8Btsb+T0x0rXAX5ibh\n",
       "TZjPrl7vDfdhvaZl2Pr/H2PLbA/2KVSXEqTuDLJQwjv/hq3fvx+4GngXNhbdyg7ceuyB+el6P7Zh\n",
       "7S9MNAC/wuJ+vB5bEimMWcD5wE+xOYCLsKW0YvIEqTuDLJToCHthY81ZzMPW/LfTIu/DnLe9C9sh\n",
       "/DATDcBPgH/HogUe0EZeIbIv5lbjJmw11mVYj6qTLvlDJkjdGWShRGkchPnEuQ5zNPbKgp8/HfPi\n",
       "ewU2pJE2ALdjnjFfwcQJYmEG+jXYhPAO4KvYxrzdfAoVKEHqziALJQrnJMwN81bMffGpTJywnQxT\n",
       "sSWg78Mmeh+l1gA8EOX3BtQDqMeu2P6AL2GuNb6J1dc+PoXqAYLUnUEWShTOUzGPnrsW8KwnAm8G\n",
       "VmHr+5MG4H+x+YHzsRVDGt6ozx7AGdgS3ZuBt2OeU0VnCFJ3BlkokZsZ2GqgD5bw7D2wHsZHMAdx\n",
       "6WGgUWwe4OWoVdsKewHvxBzhXYct2RWdJ0jdGWShREscikWguhGbD7gWW55ZxNLMPmwo6PLo2UkD\n",
       "8BAW3eutVCwEYcWZhm3Uux/4OraJT/gjSN0ZZKFEU/qwYYbLsZ2pexb03AOAczBfP0kjsBFb6jgX\n",
       "TXDmZTo2x/IHLHTk0/2KIyKC1J1BFkp0lF2xVT/XYiEanRG4H7gU8x8k8rMfZkS3Ap/B/AuJ6hCk\n",
       "7gyyUAIwJ3CXYEMzZfA0bF7g99ROCl+DzQcUMQHdi/QDH8Im3D+FbSoT1SNI3RlkoXqcYeDLmEL5\n",
       "GDBY4LMHsJUtN1M7JPQz4Cy0yqUddgf+GTOwn8EMuqguQerOIAvVozwOcxf9W+BsTHkXwRRs5/GX\n",
       "qN0tvANYjsUQ6HUfQe0wBfP2uhlbNSSXEt1BkLozyEL1MC+muKGafbC17b+mtjewDluiWqSjul5l\n",
       "HrAB2/Q3168oIidB6s4gCyXaYggbatpBbATuAZZS7PBTL/MM4L+xcJuvQhvwupEgdWeQhQqUPmx+\n",
       "4EqK30DWh7VWrwV2EhuC72G+8hXXoxieBHwWmyd4B1qC280EqTuDLFRg7ItNNt6O+ff/F4pz7LYn\n",
       "5kbix8RG4BEshsDsgvIQtnLoYmwJ6Qcpxg+U8EuQujPIQgXEdGKHcvMobmjhiZhiepDYENwHvBet\n",
       "FCqSPmy+ZQu2cuggv+KIAglSdwZZqMAocvL274CvYfsFnCFYjwWU0fBFsTwNcwtyG1bvIiyC1J1B\n",
       "FqrL2AfzIVTWEM1uwGup3TvwN8zb6PFoyWjR7AN8GHM18VbaixInqkuQujPIQnUBfViA888A27DJ\n",
       "3KKdle0CvA0brnCGYAwby9amp+Lpw1YM/Q6bON7frziiZILUnUEWquLMwdb2/xzbQFZG4Jd5WHjJ\n",
       "5E7ixWjvQFkcDtwA3IH1ukT4BKk7x7H15XP9itFTTAeOpZyhmoOB1cSGYDMW5lLDQuWwN9bj+gO2\n",
       "jFTLc8NnLqYzgzUIovvZC3tJnWuJ/wEuwILTiOLpw/Zn3AN8EYX37EWC1J1BFqoCHIjND5xScj59\n",
       "WFD63xD3Cr6CLSsV5TAIfAv4KfB8z7IIfwSpO4MslEf2xFrmW4FllLsB6RmYMztnCG4Hnldifr1O\n",
       "H/BGbHjoHOTeu9cJUncGWaiCOauFNH2Yx8q7gauBQ0qUZwbwSeBR7P/3IHA6Wt5YJjOx/+uPkTdS\n",
       "YQSpO4MsVMFsaSHNrsBK4IQS5ZiKxSJwu4sfxQzDjBLzFPBCbCnppVjMAiEgUN0ZZKEKphWDUDbP\n",
       "w4aE3PDQt1FLtWz2xAzuPcALPMsiqkeQujPIQhVM2iB0cgnnXphScobgN9gkspaRlsswtk/kK9gy\n",
       "YSHSBKk7gyxUwTiDMAXz+XMznRk6GAZ+QRyr+P1oY1nZTAXOBx4AXuNZFlFtgtSdQRaqYLYAx2GG\n",
       "4BbK34k6FXgPsQO6n2PGQZTLIPB9bOXWkz3LIqpPkLozyEIVyKHAn7FJxTdQfmSrQeAm4iGiT2Bj\n",
       "2aI8+oB/xJaTnoWil4nWCFJ3BlmoSXAW1hNIHw8Cf8RiBWRdb2VJaiu4Ne5/xP4nW7DVLaJcZgJX\n",
       "YX6fnulZFtFdBKk7gyxUwZS9ymhfbI276xWsxhSVKJeZwC+x+NFy8SHyEqTuDLJQOTkS+BT1I4WV\n",
       "aRBehPU+xrHewT+gFUSdYE9saG6Zb0FE1xKk7gyyUC0wFXgZsBZTyO+j/gavMgxCejnp94FZJeQj\n",
       "JjIF64V9Fc0XiMkTpO4MslBNOBHYRBw6stkS0qINQnI56V+B85DbiU7yUSy0pXYdi3YIUncGWagm\n",
       "HIrFuG11aKYog5C1nLSssJkim3/GggVps5lolyB1Z5CFiihqLL4Ig3Aw8D3iIaKPo+WknWYhtnz4\n",
       "YN+CiCAIUneGWKhpWEvw1xTjdbTdpaWvBrYTLyc9qW2JRF6eg+0+Vo9MFEWIujOoQu2BhTMcwyYM\n",
       "n4PfFTv7AJ8n7hVciy0xFZ3lKcD9aF+HKJaQdOdjhFKoQ7AA56uBgzzLAjZHsQmr3z9j8Qq0nLTz\n",
       "7I/9H97kWxARHKHozhpCKdRMOreG3+1SzmIq8K/A37C63Qgc0QGZxET2xnxPvd+3ICJIQtGdNQRZ\n",
       "qJKpZxAOxvYTuCGiD6Oljb6Yig3RfQ71zEQ5BKk7gyxUyWQZhPTE8YJOCyUeY1ds7mYtsJtnWUS4\n",
       "BKk7u61QuwHvxO8PPWkQpgFfIO4VXIMmjn3yOOCbwHXYkJEQZdFturMluqlQTwVuw4YC+j3K4QzC\n",
       "ccBm4onjxWh4wiePBzYAVwK7eJZFhE836c6W6YZC9QGnYa6oz8C/0t2COaJzE8e3AYd7lUgchq0m\n",
       "Wor/90P0Bt2gO3NT9ULthbmGvgPzStpJsuIf3Ac8QjxE9KeMNKKzHIv9X07zLYjoKaquOydF1QvV\n",
       "hw3F+PBXn1b09xP7IXoU67FkGQ3ROV6KRTp7qW9BRM9Rdd05KYIsVAkchAVTcQHv7/crjgBGsJ7B\n",
       "sb4FET1JkLozyEIVzCziXcc/xoyBegL+6MM2m41icwdC+CBI3VmlQr2c6gWJOQz4LVZP67Ed0Roa\n",
       "8scu2Cqi9dSPcCdEJ6iS7iyMKhRqCtbiu4dqeaN8GnF4y5uIl7rKIPhhb2x/wX9j+w2E8EkVdGfh\n",
       "+C7UNGxfwfepVotvGJs0HgfWUbvJSQah8+yP+SX6DLYTWQjf+NadpeCzUE/FQkl+imq5GDiO2A3F\n",
       "dUwMZCOD0FkOBe4ELkR7DER1kEEomDOp3trxudjegnHgKqplqHqRYzDje7pvQYRIIYMQOC8EHsbq\n",
       "5EvI/YFv5mF7DF7uWxAhMghSdwZZqEnw98Q7kK/A3CcLfxyOhbx8gW9BhKhDkLpzHPP/MrfkfKo8\n",
       "EfhqYr9EH0fj1L6Zic0ZKMqZqCJzMZ0ZrEEomxcBv6J67oj3xYLY7MTq4UPIGPhmN+A7wKW+BRGi\n",
       "CTIIOekDzsUmBZ9XYj552Qd4L+ax1Dmpu8CrRALsfbkc+C80ZCeqjwxCDvYGVmI7Sp9YUh552QN4\n",
       "FzZR6QzBt4A5PoUSj/FOzKvtPr4FEaIFZBBaZCpmCD6HHy+laXYB3ozthHaG4IfACT6FEjW8BOtJ\n",
       "HuxbECFaRAYhB0/H/3j8FOBUbP7CGYI7MFfJvmUTMU/HVhQd51sQIXIgg9Al9AEvxiKZOUMwCrwG\n",
       "MxKiOuwP3AW8zrcgQuQkRN0ZXKGOB75HbAjuxXa5VnnZa6+yO+Y08AO+BRFiEoSmO4H2C9VHNeIJ\n",
       "H435HXKGYCtwNhP9EIlq0Ad8HnMPol6b6EZkEFL0YRu5voe/MfnDgK9SG+P4QmJX1aKanAvcSvX2\n",
       "pgjRKjIICaZgXkp/BAwUJ07LPBFbs+52GD8CfAwbkxbV5hXA77CwpEJ0KzIIEVMwfz8/wOIZdBK3\n",
       "u/gvxMHurwSe3GE5xOSYje0DeZZvQYRoExmEiMuA79LZqFVZu4tXYTEVRHdwIPAbYJFvQYQoABmE\n",
       "iOPp3NivdheHwZ7AzcC/+BZEiIKQQegwCzCvl84Q/ADtLu5GpgCrgS+jDYEiHKqsOydNFQu1P/BF\n",
       "YkPwM7S7uJu5GIuZvbtvQYQokCrqzrapUqGmYD6HxjC5HgbOR2Esu5kR4NdYjAMhQqJKurMwGhVq\n",
       "D+DrwLM7IMeR1O4wvh44pAP5ivI4Ebgf2ysiRGj0lEHYE1PKqyjX7cMewP8B/hrJ8nvM55CGh7ob\n",
       "57CuSnEwhCiSnjEIfZgx+DLlBppPTxqvAKaXmJ/oDAcAdyOHdSJsesYgPBf4BeVFrdqP2knjn2JL\n",
       "WUX3sxdwC7ZnRIiQ6RmDcAVwTgl5uUnjrWjSOESmYHNOX0BDfiJ8esYgzKR453BHoEnj0LkUuBEt\n",
       "LxW9Qc8YhCLZA/M+qknjsDkDi0w3w7cgQnQIGYSczEeTxr3Ai4D7UI9P9BYyCC2yHzaOrEnj8Hkm\n",
       "5mdK/1/Ra5RqEGaV+fAGFFmoPuBNaNK4V3gC5r30Nb4FEcIDpRqE0TIf3gBXqClY13+yY/tHYC6x\n",
       "NWncG+yNRTy7wLcgQniiVIOwAtiAOQK7GLiozMwSuEKdAPyE/AahD3NprEnj3mEqcC3wWfR/Fr3L\n",
       "pAxCqzt910ZHW5m1wT9hP/C8+Z6HuZ4AC2l5HrCtQLlE9fgwFrDoFAKdWBOiVxnHIp9tBx6f895X\n",
       "R/fvRFGweoW3YbvYtVpM9DqlNoZmY/MI67FdwieXmVmCcax3cG3O+44njmv87qKFEpXkpdjy0iHf\n",
       "gghRAUo1CBuAAcy7qPveCcaxyeBX5rjnUODB6N5/R+PIvcAJmPfSY30LIkRFKN0gQGwQ1pSZWYJx\n",
       "4LW0vjR0JhbwZBy4jnK9oYpqMAczBvN9CyJEhSh9ldFyzBBcTGwYyiZPoXYn9ke0EZtYFGFzBDZM\n",
       "9ArfgghRMUpfUDGCGYWRsjNK0Gqh+rD4COPA74CDSpNIVIVZwD3AP3iWQ4gqUvrGtOV0vlveaqEu\n",
       "jNI+BBxVnjiiIhyA+aF6u29BhKgopfcQ5mPDReuBlWVnFtFKod4YpXsUeEmp0ogqMB34MfCvvgUR\n",
       "osKU7svoHGwO4Vast9AJmhVqHvC/Ubq3lC+O8MzewA+Bj6LVY0I0olSDsBNbaVSlIaMjsA1r49ju\n",
       "VBE2u2MNkk8TrjEYI/a3pUNHK8cY2YzXOV8IA9hu3+VUY8jo8cBd0fWvY87vRLjsAlwNXEV5sbSr\n",
       "QKk/YhEk9d6ZSb1Lra7TH8fGbmdiytenP6C9gP/EhrFuAV6P9WBEuHwCc2HycmyuSAjhkQ3A2XQ+\n",
       "LkLayk3BWorjwN3k928kuo9/BH4JTPMtSAdQD0HkpdAeQqsswHoFY1iAmXllZpZgHFgKzAX2BP4t\n",
       "OrcdOLJDMgh/HIVFPHuab0E6hAyCyEv6nZmL6cxS36UNQH/0eYjO+jLaFTgd23A2jq0qkpuC8BnA\n",
       "9hq81rcgHaTKBmFn6tiAOb1MsxjYFKUZw+Yd+xPXhzFX+u76qpzXRS1eeghp30Wd9GW0iXhGfSOd\n",
       "650If/QB1wCf9C1Ih6m6QTgaG7qbhS1D3wkMJtIsw5T4yVG62ZiuSDYgt2EBtmZhin55zuuiFi8G\n",
       "YTU2hzAfexE66ctoHBtDPhWtJuoVlgA/ovfiXVfdIMxKnVtOvCdpoE4asFWJs7DRhawFINdHf5td\n",
       "FxPxYhDAfqTLMYPQKcaxeAjyWto7zMUc1j3Jsxw+aPYjLmrt+mTIUvbziVvvC2itJT+K9RrqDfs2\n",
       "uy5q8WYQoLZ72Amq3GISxfMEYAumXHqRbjMIw8Qbo5ZQO3LgWvvuSDrFHMGU/lj0N61Xml0XMV4M\n",
       "wkJiyz1KZyOmid5gV+Am4ALfgnikyu97lkFI9goWMrGHMCs6VmFDzlmMMHEuIs/1XseLQUj/ozu5\n",
       "ykj0Bh8FvkFvzxNV+X3PMggrgMuiz24OIWvl0SbMICzEVhCl2YA1MptdFxPpuVVGInxOATZju+F7\n",
       "mSq/707ZD2DDQUuwXeOzEmnOwYZ5Fkbp3BLSUcwg9EfPOSd6RvI504iNSr3rYiJeDMIKbKbfeTxd\n",
       "g3XlTiszU6r9AxHFcDi2+WzYtyAVoMrve3ofwnpsGWqaEaxF79K4lr8bMhoknh9waZJLyZtdF7UU\n",
       "ahBa9Rq5JJFJXyqzSyeTcYu4/ESYPA64GRsu+rRnWaqA3neRl3rvTJDvUpVbTKI9+oCvIkOQRO+7\n",
       "yIvXZaedJshCCcDCX27EfFQJQ++7yIsMguh6jgMewCYNRYzed5EXGQTR1ewP/BZ4mW9BKojed5EX\n",
       "GQTRtUwF1gEf9C1IRdH7LvIigyC6lg8CNxB2GMx20Psu8iKDILqSlwH3YENGIhu97yIvMgii6xjC\n",
       "JpGP8y1Ixany+z5Mey5rhrEdy77yDxUZBNFV7IktL327b0G6gCq/7/2055K6XYPQbv6hIoMguorz\n",
       "sehnwe2aLIEqv+/JFvoQcCux76JRap3aLcIc2rkQmu7+0cTnZGt/AbX+0VZE944Rx18ZjvJ0bnMc\n",
       "y6I8FhCH7lxD74TdLNQg9LJnSVE+A8C7gfOotrIT+ZmN/U9nYAsFlkXnh4DLMf9FczBFvTDHcxdh\n",
       "PYFZ0d+LiR3bjWOutE9MpJ8PfDM6fxr2zm0GrshZHoEikYlyOQtzaf1L34IERBGGtYje2nbgw9Hn\n",
       "FcQKeFH0/fbo+yk583TlOwQbapwB/DFxfTVmcCD2vLo/ZpS+E50/gzhwj8iBDIIoi/2AM7FWoiiO\n",
       "qgy9JRVuUiY3nOTYGP1t1Zvt1ZgRWB39vYhaB5o7gNuwnsEhWM/gEMwQJWVSj3QSyCCIsjgPc153\n",
       "t2c5RGfZjiloxzDm0npzg3sGEp+HsNb+FdF9a6N770qkWQmchA1bLQGOAa4CTk2kyQrUI7ocWfnu\n",
       "5CBgK3Cgb0G6jCq/7+lJ5dE61waxlvrsKN0mbGw/Oans4i0PYsbgVizeCpiC34BNCrt8Ts7IYxv2\n",
       "jhE9YwzrNQxgQ1adCuLlG60yEpXnMuAS30J0IVV+34exYDVgivrOOtfAVgG5ADeXJdIk71keXb8T\n",
       "m3S+PnFtTXRtDBsyyspjNPFsMGMwGt13Pb0TYU0GQVSaIeBBYKZvQboQve8iLzIIotJ8HljqW4gu\n",
       "Re+7yIsMgqgsR2AuKnqlu140et9FXmQQRGVZTbyzVORH77vIiwyCqCSzgXuBvXwL0sXofRd5kUEQ\n",
       "leQ/gHf4FqLL0fsu8iKDICrHM4H7MM+mYvLofRd5kUEQleNrwNm+hQgAve8iLzIIolIcjq0sepxv\n",
       "QQKgyu97owA1bqdwu89vJ15CryL316JSnA98HPiTb0FEqWwCzvUthOhtqtxiEvGu5F4JRlI2VX7f\n",
       "0z2EJcRBbNxnR6NgNYuJ3VpswPwSueerh5AfDRmJynA5cKFvIQKiyu970iAMY0r9aEzZ38pER3Mv\n",
       "wDYoLsdcVLtrO4GjovtWYQFw3DNlEPIjgyAqwZMxJSCfRcXRyvu+NEqXPpa2mL5eumYkDcIKYqdz\n",
       "YI7lXA9hMbEBcCR7D8newgqyQ2yK1inUICgegpgsS4AriVuGojMsJZ9Sz5u+FaZT63k0GaugWbCa\n",
       "91BrQBrFSRAdRgZBTIYDgddivotE77EZODTxfSjxeZT6wWpcvOR5WFjMxbQeSU0IDRlVlI8A/9e3\n",
       "EAFS5fc9OWQ0mzgIjgtw41r8WcFqXKyDEeLANQPR89zwkoaMJofmEIRX9sN+8Af5FiRAqvy+1wuC\n",
       "sxULcJMcOmwUrGZNdN964qGjeZhxSQbQEa0hgyC88iFqI1WJ4tD7LvIigyC8MR1rCc7yLEeo6H0X\n",
       "edFOZeGNdwDXAnd7lkMI0YOoxVQdpgF/AA7zLUjA6H0XeVEPQXjhTGxCUBN/QggvqMVUDfYGfg8c\n",
       "6VuQwNH7LvLSc5PKS4G5fsXoed6FbTYS5RLkj1iUSvqdmUvsriQ4gixUl7EHFiv5aN+C9ADtvu9b\n",
       "okP0DppDEB3lTcBtwO2+BREixc7UsYHYTUaSxcTuuMcwh3pJJ3vDwNrE9VU5r6fz2kbs+nswcW0B\n",
       "tqs765poAfUQ/LIb8BvgWN+C9AjqIeRjJ9ZznYbtjTknOpdUtMswJX5ylG42poyTsR22Yd5bZ2GK\n",
       "fnnO646hKP9Grr+dHBfXeUZeem4OQfjjzcS+Z0T5VNkgLML8Eq1mYnAbd921wlcR+ziaH11fEl1z\n",
       "bMJcVkD9gDpDWMt8GdnKcycTN0kuJ3apPVAnDcDK6LxT4mmc/6Vm15O4/FwvZQXxrv5F1Lr+cGnb\n",
       "RQZBdIRdMH80z/MtSA9RdYOwEziNuPXrlLRTmvOIA98sx1rBLgDOWuBRYuXslGGjgDruuZeRrdSz\n",
       "lP38hFwEH6bVAAAMgklEQVQLaK0VPooZovmTvJ5kMfEQVto1fHKYaQHFLOGWQRAd4fXAd30L0WNU\n",
       "3SCsT53biSm5JdT6txok9njqFPK2KM0ItYq6UUCdeq3zZP6zUudcNDciuZLPds9zx0jimvPEOkb2\n",
       "+H6z6+75LpIc1Bq3JC5exLyMa3lRgBxROlOAC4C3+xZEZNJM6de7/oQ2870r9X0zpgRnUBvo5i6s\n",
       "5b8OU9CD2JDQDcCrMMOyNkrbLKBO3gA6SVk2URuXwckLcAm1LfYrogNM+W+KZLurxetE5VhJvADj\n",
       "DGoN2gA25DaAGYPKLdTQKiORxQjwIPaDFsIxlPF9EzY0ckji/ACwPfp8G3AuZgBuAOYAJ2KKE+KA\n",
       "OjMSx4I2ZDyFuCfjDFJy5dHd0TEn+r6Q2Dg5rojknt3C9STjQF/i+0Dq+jqsZ3QMFTQG3YCGjDrP\n",
       "gcADwNN9C9KDVH3IyA2zuMA3TvG6oRgXFGc18RDSMuLVNWBDR8mx9X4mBtRxCxmGaBw0x03gDkRp\n",
       "l1A7TwG28mgMU+wDxEtIR4Gzo/x3RumGUs+ZRjz5W+96kuRQmSuLM3yLiCfihxJHu2gOQZTKSizm\n",
       "geg8VTYICzFFvQpTkOupVbwLiYPirCRWlvOpVZ6rmLhCp15AnSEaT7ym9yGsJ3sD5QimjF2akyN5\n",
       "z46uDxLPD7g0yfH9ZteT1KuHizPkfbRB2VpFBkGUxkuwl3lP34L0KFU2CIvIniAVftGksiiFvYFP\n",
       "YcsKH/YsixDCA5pUFo6lwPexiT8h0oyjHrvwjF7AznA05t56f9+C9DhVHjIS1URzCKJQpmKTZG/y\n",
       "LYjQ+y5yI4MgCuUdwI3Urp8WftD7LvIigyAK40nYBrTDfQsiAL3vIj8yCKIwrsEmk0U1KOJ9P6uA\n",
       "Z4juQQZBFMIrgF8Cu/sWRDxGEe97r0wqD9NePIFhGu+CLjNvhwuW446ViWsujkOzYDoyCKJtpgG/\n",
       "BU7wLYioQQahdfppzR11PdoxCO3m7XCeUWdFh9vVvACTbRaxO/F6mwJlEETbfAL4tG8hxASqbBCq\n",
       "FiAn2UofivJyfotGqXU852Rz4TPd/S5d8tkLqA0KtSK6byx6fjrvtdS60V4GXEfsviJZnjT1XHsP\n",
       "ErvgGIieeVmdtDIIoi3mAPcDM30LIiZQdYNQpQA5aYOwk9g30XJqHeS5lvggpqgX0ppBWBSlcaE3\n",
       "nU+mZN4j1Lbe7wAeiuoiXZ4kTq41UR5ZcZpdnW/NuOaQQRCTZgpwM9pzUFWqbhCqFCAnbRCS8RSS\n",
       "Sn4JFg85eW02rRmEhdT2NpxSTuadDIU5APyZ2rmAZHmSzCZ2kjeI1UG9cLXLG1yTLyMxad6MtXI+\n",
       "51kO0Z1UOUBO8v7knho3nOTYGP0dbuGZV2NlWx39vQi4NJVmOxYbYT5Wll9jMRlOTKTJUs4bsbgI\n",
       "jhHMaDqZtwI7ou/nMrHuS0EGoXeYCXwAeCHFBPcW/jiLxstL6/USPhIdk6WIADmXRNeXRNddgJxk\n",
       "ZLN04Jl22J6SzRmoeoYmGdRmCJP5iuietdF9aeW8EjgJk/sbmFFoVp7ZmOG6Lfq+I3FtEWaAzou+\n",
       "a3g3QkNGxbECm0wW1aXqQ0ZVCpCTHjIarXPNDV/NJjZgpxEPGQ1G8g0ST4S7eA1Louf0J/I4mYnL\n",
       "TgcT5WpUnrT8Tq50MJ3ZqWvJ+kyjOQSRm2djiiId0k9UiyobhKoFyBmm1iDdWecamBFza/ovS6Rx\n",
       "9yyPrt0ZlSMp35ro2hjxXET6+UTyu2fXK0+aEbLrLHltLONaEhkEkYupWGvmDb4FEU2pskFQgJxq\n",
       "UqhBUDyE8FmMBbz5km9BhBDVRgYhTNyE437A+4EzUW9LtIcC5Ajv6AWcHG7Y4NPAR30KInJR5SEj\n",
       "UU00hyCasgV4DnAv9SejRPWQt1ORFxkE0ZQt2MaX1/gWRORC77vIiwyCaMoO4NsoClq3ofdd5EUG\n",
       "QTTk8di67yN9CyJyoyEjkRcZBNGQzwN/8i2EmBSaVG6dRkFqnAfVdp6dtUN6MbYjOR20xrn2Th9H\n",
       "ZzyjaGQQBGAtwS2p4wGsd3BfxjV3qAVZXWQQWqdRkJoyDIJzz5HlprufOMjNrEiuejusi0YGQdTl\n",
       "SuB99I5SCI0qG4QqB8hxz3eBbNxnR73nLyZ2aZEsT5ZBcG6unaO6FdT3L7SazvQOQAZB1OEJ2Mu9\n",
       "LzII3UrVDUJVA+Q4R3FHR/nfSuxAr97z+6NnH5WQ+eLE8+oNGbnhoK0Z19296XgIZSKDIDJZRuzN\n",
       "VAahO6m6QahqgJwV1AbBmZ94RqPnJ6OQrWBieM0kychrUD8S2gYmGq8yUYAcMYF+rOU2x7cgImiq\n",
       "GiBnOrXGKilno+e/h1rj0SivRVjL//bo+xlMNFYLsJ7D3S3IXEnkyygMTge+RRe/iCIXWQsK3EGD\n",
       "a+0uKCgiQM4cLJqYG1ZxAXJmJI4FOeXaDBxaR856z1+EGYN5WHCbq5rkMU7tvp4sV/KnYD0NURIa\n",
       "MmrO7tiP/ajEOQ0ZdSdVHzKqaoCcdECZW4lb/W4OIf38kUQ+A9Gz3BBQ1pBRchgsHdDGMUbnXcVo\n",
       "DkHUcBrwzdQ5LS3tTqpsEKocIAfiIDhbI1mSRqfe89dE96wnHjqahxmWrHzrldHJU2+iuUxkEMRj\n",
       "TAF+ha2gEN1PlQ2CAuRUEwXIEY/x95jfohs9yyGECAAZhO6lD5usuwT1pET5KECO8I5ewPo8H/g1\n",
       "FjNZhEGVh4xENdEcggDgOmzTjQgHeTsVeZFBEDwDc2C3h29BRKHofRd5kUEQfAE437cQonD0vou8\n",
       "FGoQqh5RK707UMCTse3zQ8S7QUUYjGFuGIRolW3Y7us0QepOtZgm8jHgw76FEEJUmiB1Z5CFaoMZ\n",
       "WCvyoEneP7c4UQSqzyKZ61uAwNDGtB7gTOAa4N5J3j+3OFEEqs8imetbACGD4Jjbgfuapa133Z3f\n",
       "E3gbcGmD9OlzzfIsg3bybPXeVtLVS5PnfPpcK/kWTTt5tnpvs3SNrmdda+dc2Uw2zzz3NUvb6HrW\n",
       "tVbONcuzJWQQjLkduK9Z2nrX3fk3AjcDv2iQPn2uWZ5l0E6erd7bSrp6afKcT59rJd+iaSfPVu9t\n",
       "lq7R9axr7Zwrm8nmmee+ZmkbXc+61sq5Znm2RNVnoUep9bMuhBCiOZuojREhhBBC9A7LmycRLdCP\n",
       "+YM/h9o4s2JyDBLX56BnWUJgCKtPMTnc73shTX7fvucQlqW+L8QEH2nh3uHixel6Jlufz8IiRm1n\n",
       "YpjEXmay9emied2AxREQ7f3WF5IdslIYzer2fCy+9Q008X/m0yAsxgR3OAW/Lvo7u8G9g9gOPe3U\n",
       "jWmnPl2aAWBjwXJ1K+3U59VYnIoFWDjJXqedugRTZCKbVurWGdMdNJmT9WkQLseCYztOxZQ80XkX\n",
       "aHth6ujHWrED0V91yY126nMEe1luw4Y5RHv1SXT9BuDusgXtAtqty6ovfvFJK3W7DavLAWyyuS67\n",
       "lCDgZHHBsB0zo79XZ6RdR1xA7WbOJk99bsC6mEOoRVuPPPU5H1iC/SDX1knTy+SpS7D6nIEZ2B0l\n",
       "yhUCWXW7gtjoXt7o5ioZBMjXEtgBnFSWIIHQan26YaJ1DVOJVutzHarLZuT5rV/aPIlIkK7bu6Kj\n",
       "Kb4nlZNsJ/baNx3Y6lGWEFB9FovqszhUl+XRVt1WySCsJF7hMoh1tcXkUX0Wi+qzOFSX5dFW3fo0\n",
       "CIuw5Y6nRd/dsMV8zMrd7kOoLkb1WSyqz+JQXZaH6lYIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC\n",
       "CCGEEEKIgumn1jVwkews6blCeKNKriuEKJqZwKt8CyGEEMI/KzBXwAuxcKsuDOMi4GIs9sMKYA3m\n",
       "Anxh6t4N0ZEVwCWrh5B1z6Lo/Cjw1sTnQczVuMt7fkb6Wa0XVQghRCMGgVXR5/nEMbjXYMp2CXB9\n",
       "Iv1o9HdxIu1A4nyStEGod88iTOGnPy8Bzo4+9xP7sE+mEaKjVC0eghBFkvQLvw5refdj7oHvxoIr\n",
       "JcMzbsaMyDDmMdIZk200p9E9azM+J4MRpYO+KGSk8ILmEEQvcQOmhF1Lvg84MXF9CAskcisWTvTU\n",
       "6FhFcxrd05fxeRNx/NshFPlPVAAZBBEyY5jSPTn6fjkWSvDK6LtTwm4cf0n0/QpMSbvzWXFo0wq8\n",
       "3j3jibTJz5cCx0TpVwGnZKQRQghREguAyxLfzwFGPMkiROXQHILoFRYB50V/k6g1LoQQQgghhBBC\n",
       "CCGEEEIIIYQQQgghhBBCCCGEEK3z/wHIdLp+tba7OwAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f391150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(1.-b_cumsum,1.-s_cumsum,c='black',lw=2)\n",
    "plt.plot(1.-i_b_cumsum,1.-i_s_cumsum,c='black',lw=1)\n",
    "plt.plot(1.-truth_b_cumsum,1.-truth_s_cumsum,c='black',lw=1,ls='--')\n",
    "plt.scatter(1.-b_cumsum[index],1.-s_cumsum[index],alpha=0.9,marker='+',s=200,lw=2, c='black')\n",
    "plt.scatter(1.-i_b_cumsum[i_index],1.-i_s_cumsum[i_index],alpha=0.9,marker='+',s=200,c='black')\n",
    "plt.scatter(1.-truth_b_cumsum[truth_index],1.-truth_s_cumsum[truth_index],alpha=0.9,marker='+',s=200,c='black')\n",
    "power_DGSS=100*(1.-s_cumsum[index])\n",
    "power_inclusive=100*(1.-i_s_cumsum[i_index])\n",
    "power_truth=100*(1.-truth_s_cumsum[truth_index])\n",
    "lab_power_DGSS='power DGSS %.0f %%' %(power_DGSS)\n",
    "lab_power_inclusive='power inclusive %.0f%%' %(power_inclusive)\n",
    "lab_power_truth='power ideal %.0f %%' %(power_truth)\n",
    "plt.xlabel('type I error')\n",
    "plt.ylabel('power')\n",
    "plt.title('$s_1=%d$ $b_1=%d$ $s_2=%d$ $b_2=%d$' %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))\n",
    "plt.legend(('DGSS', 'inclusive','ideal',lab_power_DGSS,lab_power_inclusive,lab_power_truth),loc='lower right', scatterpoints=1 )\n",
    "plt.ylim(5e-1,1)\n",
    "plt.xlim(1e-4,1)\n",
    "plt.loglog()\n",
    "plt.savefig(\"q0_power_s1_%d_b1_%d_s2_%d_b2_%d.pdf\" %(s1hat.subs(struth),b1.subs(btruth), s2hat.subs(struth),b2.subs(btruth),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.4590446175799521e-05, 1.05)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEMCAYAAADeYiHoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xe4FNX5wPHvQVCx0VSwJCggqDFGEHvJKhhbjCZgvRpj\n",
       "IhoTMb9osESjR00MRmNUYiHGGMuxo7EQo3ANFtSAChqjIoiCSO92RN7fH2eWu1x27929d2ZnZuf9\n",
       "PM88d3fvzJl325w9HZRSSimllFJKKaWUUkoppcK3TtwBNKEDsC/QD5gNfBFvOEopparlykb3BwED\n",
       "gCHB/eHB3w7AsGoFpZRSymsT03lPw2cIef2Cv/XB375Ax+D2MqBnleJSSikViCuD+AswveD+McCS\n",
       "4PZ0YGBwvwM+o3i3qtEppZSibdwBBDoCiwvudwFG4jMK8BmKUkqpKkpKBgFgGt1/L9iUUkrFICkZ\n",
       "xFKgc3C7E7CouQM2ZOM5n/BRt0ijUkqp2vIasEu5Oyclg7gP6I9vpN4WGNPcAZ/wUbcz6Ph5W75q\n",
       "046VbQ5iZZuD+Mp8xMYrl9Lxy8V0XrmILqtms+WqGXSXGXRvM4Pu685my/XnsMWqZXT4WGizFJ85\n",
       "LcVnSguAhcB8YG7B9iGYz8p4HjbYksSiMZXDkryYIJlxWTSmcljijykXbHmXVHJwXBnEYHyGcCrw\n",
       "V2BScH8A/mI9uZxEbmJp+/ztawCBdTuyrFNHlnXZhhmbApsBWwJbA18DthTYAr+1X8G6Sz9hwy8W\n",
       "suln0+mxYjK7tH2efbs8y/5dP2KTA4BuwfFbgiwHZgHv4xvN3wWmAW/hMxBp5WuilFJhGxdseanI\n",
       "IB4MtkK3BH/rKZ+l4AUwsAKYF2xNEthoPVZ8fT1WdO/Mkm16M7XnITzZG67sDXQHZgBvAE9+SdtX\n",
       "r2LY9Au5oj2+hNML3zX3WGB7YAOQt+GUtnDbIuB1YDKYZRU8F6WUikqONUsSNS+yX+wC6wrsJHC8\n",
       "wHCBsQJLBd4XcAKnC2wvqxvWpRPIPnD1NSA3gbwA8jHINJD7Qc4D+TbIhlHF3IRcDOdsTi7uAIrI\n",
       "xR1ACbm4AygiF3cAReTiDqCIXNwBFFHRdbNxz6E0EaoYv/gxI9vhp//4drC1A54E/gWMNWs0rss6\n",
       "QG98SWM3YE/gm8A7wPPB9hyY2dV6DkqpzKvqdTNOgq9iysUYQC+Bnws8JrBMoF7gpwJdSxyxPshe\n",
       "IL8CeQRkIcg7IDeAHAXSobrPQCmVETn89VJLEHEQ2AA4BN8Afxi+Qft1fLey8cB4s9abI22AnYGD\n",
       "gm0vfAN9vlTyKphV1XkGSqkMSNR1M0qxlyBKEVhPoK/AyQLXCLwp8LbAOQKbNnHkBiCHgPwJ5C2Q\n",
       "uSC3gnw/pvYLpVRtyNGCEkSapeaJiu8Du4/A7QLLBV4RuEHgJPHTipQ6sifIWSBjQJaBjAI5DmSj\n",
       "6kWvlKohqblutlYqn6hAe4G9BX4pMEpgscAfSrdbrD6yM8iPQZ4IMot7Qb4Hsl51IldK1YBUXjdb\n",
       "IrFVTJUQ+LrAiCCjuEVgcNOlCgDpAnI6yDMgi0BuBtkbJBN1i0qpiuXQKqb0EthC4GyB0UGPqElB\n",
       "yWKgwPpNHPl1kAv8QD2ZCnIhyFbVi1wplSI1dd1sSs0+UYF2QZuFFXhB4COBkeKnDil1lAHZPRio\n",
       "txhkdNB1NinzbSml4lez183GsvREuwhcK7BAYKg0O0WKbADyQ5DxIDNAzgdpoveUUiojsnTdTH8b\n",
       "RCWC6T+eFpgqcLnAztJsn2bpB3IbyBKQv4J8szrRKqUSJIe2QdS+oMvsHgJXC8wQmCJwhUD/pjML\n",
       "2QzkIpDZIGNBDg8G6imlsiMz183MPNFSgsxi92BCwanBZIJnNZNRrAtyEsgkkP8FXWe1q6xS2ZCZ\n",
       "62Zmnmg5gsxiF4GXgp5QTTRoQ9CoPSAYVzEH5NcgnaoTrVIqJpm5bmauDaIcQQ+o4QKzBA5pvo0C\n",
       "fLuE3B70froGZOvoI1VKVVGOrLVBiGPnuINIKoGDgvaJtwTOFb+KXnNHfQ3kj0FGcaMfY6GUqiGZ\n",
       "yiCuijuIJCuYA+pWgSXB1B4HltHzaXOQK4OM4mafcSilakCmMojZ4lgn7kDSQGBjgTME3ghKFWcJ\n",
       "NNPmIJuCDA+m8/izjtBWKvWyk0GsuItXxfGduANJk6BUsb/A3UGp4jaBPZo5anOQq4MSxdU66E6p\n",
       "1Koog0h1P/jHP+F14KS440gT47svPWvgBKAP8BZwb9D76Tjxy6g2Pmo+mF8BOwHtgSkgFmTjKoau\n",
       "lFJlk43687+xF/CxOHR9hFYQWEfgKIFxQe+nZkpl0gPkDvyCRr/QcRRKJV6OrPViwjJv1V08Jo6T\n",
       "4w6mVggcIDBPoK6MvXcGeQzkfZATdWS2UomXqQzirUdHco44no87mFoSzPn0gcA5ZR6xH8gEkJdA\n",
       "9ow2OqVUK2SnDQL413Fz2RDYVhw7xR1MrTDwBrA38BOBOwR6NHPEc8CewI3AKJA7dbCdUumX9gzi\n",
       "iU+Fg4FbgdPjDqaWGPgAn0nMACYEGcX2TRyxCswd+IbvGcBrIJeh62crlVppzyCeBb55z0c8AJwg\n",
       "jg3jDqiWGFhq4DdAL2AK8Kz4Hk/faOKoj8FcBPTFlzymBGtT6HKoSqmq8XVplsexHCuOx8Tx45hj\n",
       "qmnBYLtzBeYK3C+UU60ne4C8DPKsrkWhVOwy1QYB8DhwBHAz8NOYY6lpBj4y8AegJzABqBdwAr2b\n",
       "OOo/+IF49wD1wWSAHaoRr1KqddKeQVhuZT5w+DFzqAe6iqNf3EHVOgOfGLgaX/X0Jn7d7L8JdC9x\n",
       "xFdgbsKXODYG3g7WoUj750+ptMjhx0FkRkNRyTIeyyHiuEgcI2OMKZMEOgn8VmCRwHUCXZs5YleQ\n",
       "F4Kusd+qTpRKKTJYxQTwEPB9fG+mo8WhVRhVZGCJgYuAHfEfwDcFLhZKdRowrwD7ADcBY0Au9Svd\n",
       "KaWSpFYyiIeBo8xU5gNPAqfEHE8mGZhn4P+A/vjMYorAj6To58wImNuAXfA9nl4G6V/NeJVStWvN\n",
       "opJlMpb9xLG3OKaKq5nML7UE9hQYL34cxW5N7GlATgjmdroaZIPqRalUpmSyigkaqpleBJYDh8Qb\n",
       "jjLwErAv8GfgUYGRAl2K7Clg7ga+iV/57r8gA6sZq1JqbbWUQTwM/MBMBWAEMDTWaBSwenrxO4Ad\n",
       "gBU0LFZUbFrxBWDqgLOAW0H+DlIkQ1FKVUOSM4gewIAK9n8D+BJfn30vsKu4pvrnq2oKRmUPBQ4A\n",
       "DsOvbHdE8eVPzWh8l9hlwP9A6nQktlK158pG9wfhL/pDyjh2WLB/KWvXpVmuxHI5gDiuEMd1Zcap\n",
       "qkzg0GDp03+KH09Ras/dQV4HeRSkW/UiVKomJaYN4jTWvMDnB7DVB3/7NnP82Bac8+GCc94EnCQO\n",
       "XfUsgQw8AXwL+Dd+NbvLBYo0TpsJ+F5Rr+MnADxOSxNKVUeUGcRfgOkF948BlgS3pwP5RshBjbb8\n",
       "GIaWXAQmABtj2dHU8QHwHE2XQlSMDKwwcBW+q2tv/PiJQWtXO5kVwQSA38VPHniftk0oFb1qtkF0\n",
       "BBYX3M9/wUc12pYFjw/A/3Isf9CbZRXwIDA4eOROdM3qxDMwy8Cx+PErFhgjfhxF4z0nArsCs/Cl\n",
       "iYOrGKZSmVPtRupKSgVXARfQkGGUqzCDeBzYRRy6eE0KGF/d1Bd4FHhG4Cqh8Xrj5nMwZwMnA7eA\n",
       "jABpX/VglcqAamYQS4HOwe1OwKKIzvMisCmWPqaOz/GlkjLWV1ZJYGClgevxvZg2xzdkH12k2qke\n",
       "34axGTBRpxJXKnzVzCDuo2Hpym2BMSGkaQu2XPDIKnymsEY1k7gWtWmomATTdpwMnABcDDy59rTi\n",
       "ZglwPL60+TTIUG3AVmoNOda8TibGYHybw6kFjw2h/G6uzSndXcuyP5ZJAOJoI473xbFLCOdUMRBo\n",
       "J3C2wMJg1tgivZ2kVzA77GiQZmaTVSqzKurmmmZCYcmhkGUdLHOxvn+9OH4rjj9WNToVOoGtBO4T\n",
       "mCbw7SJ7tAP5LcgckMOrH6FSiZXDXy8zlUGUZrkRy/kA4ugjjjniaFuVyFSkBI4U+FDgBqHYOBfZ\n",
       "H2SGNmArtZbEDJSL2wME7RCmjinAB1Q2dYdKKAOP4Bux2wOvi58QsHCPZ/EN2JsDE0DKWDtbKVVL\n",
       "SlcxAVjaYpmH9Q3j4jhNHGO0sbq2BPM5zQ1GYjeaAFAMyI9AFoCcqQ3YKsNyaBVTI5abCqqZ2onj\n",
       "LXF8N/LIVFUJdBP4l8B/BLYrskcvkIkgj+gIbJVxmapispQqQXj34UfoYur4EjgH+KM4dHnLGmJg\n",
       "Ln6GWAe8IHD6muMmzDT8EqfvAJNBcjGEqVScciSsm2vUyilBrINlNpY+AOIw4nhSHL+IPDoVC4Ed\n",
       "BF4WGC1QZPZXOSTo5XQJyDrVj1CpWGkV0xos12G5ePVBjp3EMV/c6lHdqsYE4yYuF5gjUGS+JtkC\n",
       "ZBzIGB0zoTImUxmEpekqJrDsjeVNbEOVgzhuEsf1kUanYieQE/hAYHiRBuy2wZiJWSBFxlQoVVNy\n",
       "aCN1EZY2WGZgWT1Xjzg2E8cicWwTVXAqGQQ2C6qbXhToXmSPQ0DmgvwaJO1tcko1J1ON1M3zczPd\n",
       "T9BYDWDqWIBfUOiiuMJS1WFgAXAEfpbfCQJHNdrjX/hp5Q8DRoNsWu0YlVLhKz8ntPTHMrVRNVNn\n",
       "cSwUR89IolOJI7CnwPsC1wmNe7JJO5ArQWaC7BlPhEpFLlMlCEtzbRDeK/jnml/2FFPHYmAEfoUy\n",
       "lQEGXsKvN7EN8LT4kdb5/34J5jxgKPAoyM91YJ2qITm0m2sTLL/DctUaCTg6iGOBON8NVmWDQBuB\n",
       "ywRmCMVm+ZVeIK+B3AmyYfUjVCoymSpBVOJe4Fhsw3M2dSwDroWGbrCq9hlYZfx7Pgy/vOnxjfaY\n",
       "BuyF/zKNB9m26kEqpVql8u5aljew7LNGIo6NxTFPHDuEFplKDYFdBN4RuFWgUWlBDMhZQS+ngfFE\n",
       "qFSoMlWCsJTXBpF3L41+LZo6PsIvcXluaFGp1DAwGd821QZ4VQraqcAImOuB44A7Qc7TrrAqpXJo\n",
       "G0QzLL2CGV7XWBdCHJ3EsVgcXwstOpU6AscJzBf4cZH/fh1kPMgTIJuv/X+lUkEHyjXJMhHLWtUF\n",
       "4rhaHH9qdVQq1QS2F3hb4MYSXWGvCEZfHxhPhEq1SqaqmFriHtZqlATgT8DJ4tDpoDPMwNvA7sBW\n",
       "wL/XnPDPfAnm18ApwF0gl/kpO5SqTVnMIO4HjsKyXuGDpo4PgYeBM2OJSiWGgeXA94Ex+DUmvtVo\n",
       "jzH4toq9gXqQraodo1KqaS2vS7M8g+V7ayXo2D6Y6VX7visABI4VWCCs/Xnx04XLhUEvJ13OVqVB\n",
       "ptogLJX1YvIsZ2C5u2iijlG6XoQqJLC7wIcCw9ZciGj1HjmQeSBFGreVSoQcOptrmSybYlmKZZO1\n",
       "EnXsIY73xaF1y2o1ga8JvBY0Xhf5bEgfkHdBfqddYVWCaSN1sywLgXp8//Y1mDr+A8wEBlc7LJVc\n",
       "Bj4A9gN6AQ+vPajOTAH2BA4A7gXZuNoxKhW2bGYQ3q0U7e8OwNXAMHHFqhNUVgWN14fjpxCfFKwz\n",
       "MVrgcYEjwSwADgSWARNAvhFnvEplWevq0ixtsXyIZa0vsTjaiOMtcRzQqnOomiR+iHVO4PBgOy6Y\n",
       "+G94Q/WTnAKyAOSEeKNVag3aBlE2yxVY/lg0ccep4hjd6nOoTBDYVOAp8dOIB+tcy7e0XUIljLZB\n",
       "VOA24ERs4xGzANwF9BO3dglDqcYMLAQOBcbjx058A8xrNLRL3A2yfpwxKlWpbGcQlqn4kbPfbfwv\n",
       "U8fnwJ+Bc6odlkonA18ZvwDVhfiSxIEF7RLgB9XpPE4qNdKeQVhaMg5iTX+jdGP1TcAR4tipledQ\n",
       "GWLA4ddAv0fgZDCfAycATwOvgOzTZAJKhS+HzubaApYNsSzBUnS6BHGcIY7nxaU+M1VVJrCDwHsC\n",
       "FzUMsJPvBoPqztElTVUMtA2iIpZP8HMwrTUmIjAS3zPllKrFpGqCgbfw8zUNAm7yPZzM48Ae+BLG\n",
       "AyDt44xRqaZoBuE5fBXAWkwdq4CfAleIY7OqRqVSz8Ac4NtAT2CUwAZg3scPuvsC3y6xaYwhKlWS\n",
       "ZhDeOGBLLNsX+6epYzK+V9PV1QxK1YaCAXbLgGcEtgDzBXAS8AzwAkjPOGNUqtaEO+DD8icsl5Y8\n",
       "mWMjccwUx16hnldlRjDA7kKBmQJ9C/5zBshsP25CqUjpQLkWseyGZSq29PQa4viFuOKzwCpVLoHB\n",
       "wRTi3y94dHDQeN0/vshUBtREI/W2wABgWHC7Gl4O/jb1Bb0dOFQc2pddtZiBB/GD6kYInO17OJkH\n",
       "gdOAf4LsGW+ESnlRZxBXNro/CH/hH9LMcf3wF+yxVGtWVYsAd1OisRrA1LEUeIjS4yaUKovxn++9\n",
       "8b3j/hz0cHoEOBl4FGT/WANUKmKnAdMK7vfDZxDgM4i+ax2xtmHANiX+F35dmmV7LHOwrFNqF3H0\n",
       "D9aLKLmPUuUS6BDM4fS4wEbBowOCif6+E290qgYlporpL8D0gvvHAEuC29OBgcHtQY22DsHjA/El\n",
       "iPcjjHFNlreB2TQxOtvU8TIwHzikSlGpGmZ8z6bDgXnAvwU2B1OPb5+4C+TIWANUmVbNVdM6AosL\n",
       "7ncJ/o4qsu8A4Fx8RjKmxD5RuQc4Hr+gUCk3AWeAzvaqWs/AlwKnApcB4wUOMZjnQQ4DHgfZCIyL\n",
       "OUyVQdVupC53aoF64Dv4AWrVzBwAHgCOwtKuiX3uA/YUV7L6S6mKGD/vxm+APwLPCVwmmO9OofdD\n",
       "p/C3a4Mpw7VaU1VVNTOIpUDn4HYnYFEVz10+ywzgXSi9WJCp41PgTnwGplRoDNwM/BBYBdCbqfNH\n",
       "cvqS6zirzrDqMZCO8UaosqSaVUz34buQ1uO7ro4JIU1bcHtcsIXhAXzvqaea2GcEMEEcvzd1LAvp\n",
       "vEphfNvb2Px9YeVNZ/LnsV2Zt9nx3DNBkEPBvBtjiCo9crR+xutIDMa3OZxa8NgQyuvmWo7oRgRa\n",
       "tsWyANt0BiqOO8Xx68jiUCoQrFg36WlyYwxffQBSdFoYpZqRqZHUlqhyR8tE7OqeVsUDcOwojnni\n",
       "2DCSGJQqINBZYMJEdq1fhy9ng+wcd0wqNXL462UkGURf/JiGifixCT+I4iQVijYntJyL5eZmg3CM\n",
       "Esf/RRqLUgGBTQT+/Ta9n1+Xz+eC7B53TCpVIrluvozvpnp/wf24RV2C6IFlXhnVTLuKY5Y41osk\n",
       "DqUaEWgv8OgHbDVxAz6e7xchUqpJOSIsQeQzhHwG0VTjbbVEX5dmeQW7ej3h0oE4/iWO0yKPR6mA\n",
       "QDuBu+fQ9d+Gr2aDaI86VY5IRlK/gu9+1xEYju+ymgX53kzN+R1wnk6/oarFwJfAKd2Yt9HbbH8b\n",
       "cDbIlTpWQsXlNPzym2H0QApDtFVMAJaeQTVTU4PmfDCOieI4NLJYlCpCYBuBeRPZ9XCQepCnQbrF\n",
       "HZdKnBwRVjElUXWeqOU5LM3OhyOOIeJ4qBohKVVI4DCBD2aydTeQy0A+BPl23HGpRIqsDaJwmxjF\n",
       "SSpUrQziR1gea243cWwsjiXi2KIaYSlVSOC3As8KbAxyMMhckO83f6TKmMivmwPx7RBxi76KCcCy\n",
       "IZbFWLZqNiDHSB04p+IgsI7ASIFXBLqC9AOZD7JH3LGpRMhRxSqmbPRiyrOMxDZ/4Q/WinhPXGJX\n",
       "6lM1LFjz+hKBaQK9fPdXmQ3SI+7YVGJUdN0sdy6mYQW3u5Tcq3bdCtyNZTjWT6JWwiv4NS8GkoxM\n",
       "VGWI8V/+SwXmAuMFM3kaPZd+yFav96HrJVsw749xx6jSpdxfutMLtonA0ZFFlEwTgU+BJhv+TB2C\n",
       "XygpKT29VAYZ39twIHBNL9795T846snPWf+Kd+j1+7hjU+nS3PoMpdoaBLgg5FgqJcClhDuLa2mW\n",
       "XwC7YTmxyaAcmwAzgO1NHfMij0upZonZm/EX3MlJl82g+20HMk5/wGRPLtguofx1eZrdcTBr11mZ\n",
       "4LFqL+TTmFDBE201Sxf8OhHbYlcvnVqUOEYA65o6Tq9KbEqV4Rju/cFwzr//efZ96ofcebhvslAZ\n",
       "U5XrZt+oT1CG6n+4LfdhOaO53cSxiThmiOOgaoSlVLmu4Pz+M9n68ysZ9jrIJnHHo6oukuvmIHyj\n",
       "65PB36lRnKRCcWQQ38XyfDm7iuPgIJPQL6FKlP/yjV7z2Gz5uQyfC7JD3PGoqorkuvkU0A8/H1O2\n",
       "xkEUsrTDMh9LWd0GxXGLOEZGHZZSlRLYZhkbLxzKdR+BlFxeV9WMHBGOg8h32cyvj3B/qR2rKJ76\n",
       "U8sILL8pZ9eCqqbvRB2WUpUS6PEp6887HrdM15XIjEhmczX4aibwXTizvHD6XcCJ2OYbekwdy4Gz\n",
       "AO1eqBLHwPT2fH7obZzCTvx3NMhOccekkqXcDOIg4FXgfKAncF5kESXfBHyGuVuZ+z8OdBPHN6IL\n",
       "SamWMfDqeqz45XPst2ITlj0F0ivumFRylJtB3A90wK8DcT4wKbKIks4i5EsRZTB1fAU44KQow1Kq\n",
       "pQz8rSPLRr/A3nMMq0aDZLmGQBUoN4P4C3AhfibXX0Hme+Y44Nhy1okI3AHU6YJCKsHO+gZvrrqB\n",
       "n88H7taFh1RLjQS+ijsI4ujFVMgyHsth5e4ujlfFMSDKkJRqDYHuq2DukTz8KkgSeiqq8OSIsBdT\n",
       "X3zX1mn4nkxJuNDFOwrU8lNs+b25xPF/4rg9ypCUai2BgStpM7c7780AOT7ueFToIuvmOqjZvaor\n",
       "7gxik2CdiLIWCBJHV3EsFcdGUYemVGsInL+cjV5fj8/mg+wbdzwqVJFeN5PUXTP+eWQsN2G5pNzd\n",
       "xfG4OG2sVskWrCvx0Pt8/ZkzuX75RHa9VOBU8R1VVLpFMg4ir2eF+9e6G4DTKmisvhPtzaQSLlhX\n",
       "4kfdmfnKj/nbxDfY6ZxPaX8KcFPcsalkS1I1U/wlCADLM9jy1scQR3txzBbHsVGHpVR45Jcb8PHb\n",
       "K2g7TZJ1DVCVi/S6uW2UiVcoKRnEMdjy16MQx87imCMuc4suqVST3+zHM3O+oN1Cgc3jjka1WGSz\n",
       "uU7DN1ZPA34QxUkqlJQMYl0ss7GUPU2BOL4ljrniGBxlaEqFS35wNWd/+hZ9XpZqrsWiwhTJdfPl\n",
       "Zu7HId5xEIUsl2K5sZJDCjIJnUlTpcYO/G+XKWz3xcXYepD14o5HlS1HFWZzLXU/DskoQQBYtgq6\n",
       "vFbUy0Mcx4jjRXH6a0ylxwT677OYjp+fxs3vgmwTdzyqIpH0YnoPv1jQMBoyhyHAqZWcrGZZPsS/\n",
       "Lj+q8MhRQCeSUApSqky78/L4jfi4/1UM63IW170OckTcMalolJtBvAvUB7fHAmPwU353iiKolLoe\n",
       "GIotf76lYCK/4fh5rpRKjXVZ+cYmfNRvOOcvv5SL7wAZAdI+7riUyktOFROAxWCZiKWiX1PiaBcs\n",
       "KrRHVKEpFRWBLb5knbev58zJIG+C9Is7JtWkZF03I5S8J2o5EcuYSg8Tx5nieCSKkJSKmsCWAtNG\n",
       "8f3bQeaDnAtS6SBcVR3Ju25GJHlP1Hd5nYOtbHGgYADdHHF8M6rQlIqSQHeBGW+y/bkgz4M8CdIt\n",
       "7rjUWpJ33WyBDvgZY4dRev6XZD5RyyXY1Wt3l00c54njwShCUqoaBLYT+HAVvDaLLef+l298+R7d\n",
       "6wW6xB2bWi3SuZgqdWWj+4PwF/4hzRzXHz/WYinQI4K4onQzcAyWzhUeNwL4po6wVmllYCqws4GT\n",
       "t2L2IXdx4s8f44i9VrLO6wI6K6xaw2n4Udd5/WiYx2UIfo2JpnTAlyBKSWYJAsByB5ZfVXqYOPYQ\n",
       "xzxx5U0hrlTyyYUncfukVTBP4AIdgR27RF03CwfUDQcODG7nq4/AZxqFWwcaShiF+zWWqCe6Bsse\n",
       "WKZhKy+hieMycfxTB8+p2iBtQV7ai/EXCkwWODfuiDIuUVVMhToCiwvu5+slRzXaluGrlwbgq5ce\n",
       "qGKMYZmAfx7facGxlwNd8SUwpVLOrAROfpG9f/krrhoKnCVwVNxRqWQoLEHcTEO10gB8iaI1kluC\n",
       "ALD8BMujLTlUHDuKY6E4eoUdllLxkKEgLy6k8x4CC6T5KmYVjcSWIJbC6obbTsCiKp47DvcA+2DZ\n",
       "ptIDTR1vAr8D/i6u/JHZSiXYDcCSTVl0NHAG8KigbW1JV80M4j4aeiRtC5UPKCvCFmy5ENILj+VT\n",
       "4A5aXlV0HbASOCe0mJSKjVmFX01xkPE/Yu/AT0+jopVjzetkYgzGtzkUTug3hPK6uZYj2VVMAJY+\n",
       "WOZhadG0yOLYRhwLdACdqh3SH2R+X17ZWeB9Qae7r7LkXzdDkpz1IJpiGYPlhJYeLo6fiGOyONYN\n",
       "Myyl4iNDQP43l81PFHhdoG3cEWVAjgjXg0iidDxRy1HBJH4t6rYqDiOO0eIqH1ehVDKJARlh+Grm\n",
       "TLZ+awVtz4w7ogxJbCN1FCxJL0HAo8D6wCEtOdjUIcClwM+0wVrVBiNghgptjj+aBz5dzibX7s5L\n",
       "x/uMQ0UkRwvaINL8hghpid9yDHA2sBe28pJPMGhuInCxqeOfYYenVHzEvMbOj61g3QNnsfWivXjx\n",
       "xa7Mnw9cZWBG3NHVoIqum2kvQaTFKGAT4KCWHByUIm4EfhZmUErFz8jO/Pe4HXlz6Gy2fNFiD32Q\n",
       "QX0FJggcGXd0Kr3S0UidZzkBy/hWtEVsEAye2zbs0JRKDtkC5IWfcMufgl5O1wot6wWo1pBDG6kT\n",
       "zLIOlinY1fNRVUwc14hr9Qh0pRJOtgNZ+EP+vovAIwLPCWwad1Q1Il3XzVZI3xO1/BDLMy09XBy9\n",
       "xTFfHOuHGZZSySPngTx1P4PbCPxeYJpAn7ijqgHpu262ULqqmAAsbYOBc9u0NAlxPCWOE0OMSqkE\n",
       "krYgr4KcDCDwY/FThuvAupbJoVVMKWC5F8spLT1cHEeK43UtRajaJ32DNa67AQgcGEz0t13ckaVY\n",
       "psZBpNEbxDRqAAAQPklEQVS/aV2p51HgHeDaUKJRKrHMJPycZC+BHGDgafyv4Hu14Vo1J60liD5Y\n",
       "Zra0NxOAODYRx9viODnM0JRKJjkMZBbItV2Zs4HAwwLXxB1VSmWqBGFJUxuE9w5+7pkWd1c1dSzH\n",
       "r753tTh2CSswpZLJ/BPYGdh8Ht2eG81hPwMGCRwWc2BpkiNhs7lGLZ0lCADL3Vh+0tpkxHGcON4V\n",
       "x0ZhhKVUsokBeQjkaoH9gkbrBwu2Fk1nkzGZKkGkVWvbIQAwddwLvAEc39q0lEo+I/j1VY43yLr4\n",
       "kdb3BtvTwK2C/lhSXppLENthmdWadog8cRwqjolhhKVUOsjBIDNBOq3xKNwtfmJLVVqmShCW9LVB\n",
       "AEwL/vYMIa2ngM3F6Rq/KivMk8DDwM2NZoC9ADhTYOt44kq0HNoGkSKWu7ChrKyHOC4Wx41hpKVU\n",
       "Okh7kMkglxdmEgJXCPw9xsCSLlMliDQLpR0i8DfgOHFsGFJ6SiWc+Qw/O/JRQGEmMRw4WGDX2EKr\n",
       "IZpBxGcccEAY7RCmjlnA88CxrU1LqfQwC4ADge8BvwUxBpYDFwN3CgwV2EHSsm5MAmkGEZ/pwErC\n",
       "mzbgFginykqp9FidSXwXuCh48FbgEvzYiSeA2QJvFGy6fG8GpLsNAsByB5afhpGUONqKY5Y4dg4j\n",
       "PaXSRbqCvA+yRilafN/Y7gI7BdsewfiJfjEFGrdMtUFY0tmLKa8eWr4+RCFTx0p849wPw0hPqXQx\n",
       "8/DjIkaA7Lb6Ud84McPAG8H2H2AY8DeBdnFFG4Mc2ospZSxbY1mADSejFkc/cUwN1rBWKoPke8G8\n",
       "TSW7ugaliicELqxmZAmRqRJEullmAYshtGqhScD6wPYhpadUyphHgeuBx0A6FN3DXyRPB34psGM1\n",
       "o0sbzSDiVw8MCCMhU4fgpwPXxd5Vll0FjAceASm6boqBmfjeTrdoL6fSNIOIX2gZROARNINQmWYE\n",
       "OAuYA9zrV6cr6magA/CdakWWNppBxG8csC+WdUNMbwdxdAspPaVSyKwCTsZXuf6l0ZQcfg9YBVxB\n",
       "NtsiyqIZRNwsi/BzM+0eRnKmjhXAk/h+4UplmFmBXzdlJ0qPfbgf2FJg/6qFlSKaQSRDFNVM3wsx\n",
       "PaVSynwCHA38CmTvtf7rB6sOR0sRRaU9g7CkexxEXtgZxBNATudmUgrAzABOxbdHbFpkhzvwU3Ls\n",
       "VuR/tSKHjoNIKcuGWD7GhndBF0e9OI4KKz2l0k+uAhkNstYP42DepofjiKrKdBxE6lg+AV4F9gsx\n",
       "Ve3NpNSafg10xDdad230v78CewqMFRgTbCdWP8Rk0QwiOcKuZhoFHCyO28WxVYjpKpVS5kt829wn\n",
       "wJsgV+RXpTPwGbA3cCXwB/wypheVSikrNINIjtHAYCyl+mxXxNTxIdAH+BB4PVhUaJ0w0lYqvcwi\n",
       "ML8A+gKbAfNBBEQMMsUgXxgYg19jZSPRWQlSq3baIPIsz2CpCztZcWwjjtfE8f2w01Yq3cQUbENB\n",
       "3Or/wA0C58cZXQS0DSLFhgPnh7GIUCFTx/v4UaODw0xXqfQz0rBxD3A4yCbBP/9Bxtvxkp5B3Bx3\n",
       "AFX2L+Ar4LAI0n4YOFwcReemUUqZhfilgI8OHngG6COwRXwxxSvqDOLKRvcH4Rtiy1n5LHsLelgE\n",
       "X4q4IOykTR1zgdfw6/gqpYr7O36KDgyswI8pyuyg0ygziNPwGUJe/oJfH/zt28Sx2wJLgKURxJV0\n",
       "DwLdsOwbUdpazaRUaU8A24P0DO7/A7I7nijKDOIv+HWX847BX/QJHh8Y3B7UaOsA9MD3V+6Bzyyy\n",
       "w7ISP11x6KUIfNfXI8SxXgRpK1UDzAp8W0R+ZcZ/AfsIbFL6mNpVzTaIjvjFcfK6BH9HNdqW4UsZ\n",
       "04Njaq+3UvNuB/bAhjt+wdQxG3iTcMdbKFVrbgd+CNLGwEfA88AhMccUi2o3UlfSO2cZfp7296MJ\n",
       "JcEsnwMTgf4RpK7VTEo1bRLwMQ0zG2S2mqmaGcRSoHNwuxOwqIrnTqNXgF0jSPch4EhxmVqwXakK\n",
       "GAFuwy9LCn6VxkPEV39nSiijdst0H/4XcT2+XWFMCGnagtvjgq1WvIKfgTJUpo6Z4pgKHAA8FXb6\n",
       "StWIvwIXgPQxmCkCj+FXqbs85rgqlSOhM14Pxrc5FF7khlB+N9fm1HbbhOXrWOZEkbQ4zhXHDVGk\n",
       "rVTtkItA7gQQ2E5gQQ2UImr7ullAqJ31INZmMVgWYtky7KTFsZs4/ht2ukrVFukAsgCkN4DAHQK/\n",
       "iTuqFsrhr5eZyiBqm+VJLEeEnaw42opjubjVPcmUUkXJRSB3AAj0roFSRKbmYrLUagnCi6Sh2tSx\n",
       "EngRIhmMp1QtGQEcCtLbwDv4gXRDY46pJXK0YEW5ajZSR8HGHUDEXiEY9h+BZ/ELtT8SUfpK1QCz\n",
       "DOQ6fNXSScBvgfECE/DzpjW21PjvbdKMC7ZLKjko7RlErXsVuD6itJ8FrokobaVqyQj8AkPfMZin\n",
       "xH9vziux704CZxjfnVzFqLYbqSHfUL0IS7ewkxbH+uL4WBwbh522UrVHciBzQJqc3UBgN4H54qcJ\n",
       "SpIc2khdgyxjsBweRdLieEYcB0eRtlK1Ry4CeRakyZoXgbMEXhYSOedZphqpsyCqEdXg57vfP6K0\n",
       "lao1V+DXrr6smf1GADOAqyOPKGJpzyAstVzF5L1KdGtj5BuqlVLNMquAE4GTQA4ouZf/lf4T4DBJ\n",
       "zloSOWq/U88aslLF1BPLB1EkLY4Ng3aI9lGkr1RtkqNA3iijqmmgwHQhUas4ahVTjZkObIRl87AT\n",
       "NnV8ArwB7B522krVsEeAOcAZTe1kYCwwGTinGkFFoZLpt5NGSHf85bPU4yc7fCl4ZD6WuWEkLY4/\n",
       "AB+ZutRNQqZUjGRHfBvejmAWlNzLT0z6MvAtA7OqFV0TKrpupr0EYan9NgjwCymdCdwFPIAfzRkW\n",
       "bYdQqmLmTcDRzOyuBt4DbgT+UI2ompCjBW0Qaf4Fnp0SRCHLRsB8fLXTqtYmJ45NgXeBjqYuI+06\n",
       "SoVCOgJv45dKfrPUXqcxcoMb+dl/ptFr6AMcPQlAMJ9czG9LljwilJnrZnYvZpZZWLYJKzlxLBQX\n",
       "/mA8pWqfnAiyGGRJU9tg7v94EZ1WLaHDqiV0WLWMjWU45+7XfPrhB1zJzjrVRjpNAfoQ3nKs+fRC\n",
       "addQKjvMXfiq3yY9GGx5D3Pky9sx9Q/AXlFFFgZtg0in/AU9LG+HnJ5Sqgnv0Pv8/Xhuj0v5TbWm\n",
       "Ds/RgjaIWsggxsUcQxzCziDCTk8p1YTzuGrsdHos3pLZv6/SKceRwQwiq6YA24ecnmYQSlXRJPre\n",
       "0p+XT4w7jqZoBpFOYVcJaQahVJXNo6vdnPkbXM3Zx8UdSymaQaTTTKBz0OU1DNOBr4lL5OyTStUk\n",
       "y2VfPMd+/96RNytaxKeaNINIIz/+YRrQO4zkTB0r8JlOzzDSU0qV5116nr0XL/a5nAu7xx1LMWnP\n",
       "ICzZ7MUE2lCtVOpdyO//+xJ7Tj2PK9//nHWlue1n/PlLkC+CbVwFp8qha1JnimYQStWAl9hz+wns\n",
       "3mx18S/505HXc9bxN/LzHwDtgDl+NLdZWsZpxqFrUmfKFOCwENN7G9gnxPSUUmW4lEsF+Ki5/S6G\n",
       "McAIwaw08AXIf/Df2dFRxZb2KqYs0xKEUhliYB5+2yl4KPKJNjWDSC9/QbehvYeaQSiVfIWZgmYQ\n",
       "qgTLMmA5sFVIKS4A2gSzuyqlkqkwU/gP8E2QDaM6mWYQ6Rbar/5gqu+wR2grpcL1LLC/gAHzGX7F\n",
       "uj2jOplmEOmm7RBKZYjx45U+o2EMVKTVTGnPICzZHQcBmkEolUUtaYfIkcHJ+izZnM01TzMIpbKn\n",
       "MFN4AdgNpLlpcsaRwQwiaXJVPl85F/RcBelVa12IXBXOUalc3AGUkIs7gCJycQdQRC7uAIrIRZRu\n",
       "QQZhluO/t/2jOJFmEOHKVfl87wNdsWzQxD65CtKbBmwjjnatCaoMuYjTb4lc3AGUkIs7gCJycQdQ\n",
       "RC7uAIrIRZTuVGA9gfz8TZG1Q+hI6jSzrMTyENAZ+LS1yZk6vhDHh8DPxLGwyC4rTR33tfY8SqmW\n",
       "MyDSkCncib99OhD64kOaQaSdpS7kFK+ldLe5z0AzCKUS4FngUIEXRzHowwv4/b4/4oKBbVm5CuBL\n",
       "2s28iCumtfYkptVhxmcy8K24g1BKqRR5Ddgl7iCUUkoppZRSShVzZaP7g4ABwJAYYlGVG1ZwW987\n",
       "pQqksZtrkr7Ep+HjyesX/K0P/vatbjiAf12GAMMLHov7NRscnP/mgsfijglgIHBQcDsJ7x00/OAo\n",
       "fF3ifq36BTEkJaZ+sHrZ3WnATQmIqdT5445pGK1479KWQSTlS5z3F2B6wf1jgCXB7en4C1A1DQDG\n",
       "ArcAPYL7+dcortdsQLDVBzH1JTnvoxTcPpZ437u8Ifh+7u8G95PwWp0PjAI6koz3rxP+2tULOBr4\n",
       "QwJi6ov/3NQHf/sS/3cv/xkehV9vflsqfJ3SlkHEfQFuTkdgccH9LlU+fw8aXpPpwf1jgaUFj1X7\n",
       "NasHzghudwYmkYyLcV8aviQAHYj3vcsbAmwHPB3cj/szPxiYGNy+Cv/+xR1T4fvWH3gvATFBQ+mv\n",
       "B8n4nA+k4YfGu8H9Y6jgepC2DCLuC3A54uw6fEuwgf+l8DL+NVtUsE8cr1kHfFH39wX3434fOxd5\n",
       "LAndvjvjS1z5tpG4P/P9g3P2TVBMeQNoGJcTd0yT8BnV4oI44o5pUcE5O+FLERXFlLYMApLxJS5l\n",
       "KQ0Xnk6seWGupn7AK/gPLcT/mi3D//o8HV/MhXhjalx6gOS8d7fgY+uCvwBC/O/fQho+S/k2t7hj\n",
       "At9+tLzgfpwxdcS3hwzBv4dJ+Jw/iM8UwJdqKv5Mp20kdVK+xKXch//FVY//gIyJKY4BwAXB7bhf\n",
       "s374uv5JwKv4Kou4Y+oRbF2COPqSjPduCP7X3Sj8a9KD+F+rRfhfxgSx7JaAmPL6FdyOO6YhwEh8\n",
       "hrWUZHzO38N/rvsGsUyn4TNfVkxpK0Hch//SQLwX4LzB+IvKqcH9/K+sAfg3ZHIMMZ2G/7WejyPu\n",
       "12wADR/Ijvi60LhjGhVsgq/uymdgEO97Nx3fyQD8F3ki8b9WDxacvyMwIQExUXD+vCTElC/N1OM/\n",
       "Q3HH1Bd/fZqEf+9GJSCmyA0h/u6RSTUQ/wt0WvD3wODxOF+zDjR0vS2cTEzfx+IGBduvCh6L+7Ua\n",
       "go8pSe/ftjR0b82LO6ZiXUrjjin/eSqcXiPumJRSSimllFJKKaWUUkoppZRSSimllFJKKaWUUkop\n",
       "pUKXhPlUlKoFg/Ezd3YM7o+KMRalQpG2qTaUSqL8NOv1+LmKFuJHqg7CjyRXKpU0g1Cq9QYDTwW3\n",
       "ewKH4jOLsfi5sZRKJc0glGq9/EJI4Cf+y1czLaNhumWlUidt030rlUQj8VVMHfEzsuZniTU0rOil\n",
       "VOpoI7VS4RmCzyCm07BWwVh8SUIppVRG9cAv8fqDuANRSimllFJKKaWUUkoppZRSSimllFJKKaWU\n",
       "UkoppZRSKmv+HxOzernk7vB2AAAAAElFTkSuQmCC\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10dfd6490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(edges[:-1],1.-s_cumsum,c='b')\n",
    "plt.plot(edges[:-1],1.-b_cumsum,c='g')\n",
    "plt.plot(edges[:-1],1.-i_s_cumsum,c='r')\n",
    "plt.plot(edges[:-1],1.-i_b_cumsum,c='orange')\n",
    "plt.xlabel('$q_0$')\n",
    "plt.ylabel('p-value')\n",
    "plt.semilogy()\n",
    "plt.ylim(0,1.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes:\n",
    "\n",
    "For bi=100, s1=25, s2=0, the b-only q0 distributions aren't the same, one of them isn't a chi^2. Probably DGSS distribution - numerator is invariant, but denominator of both depends on shat... but seems like inclusive approach probaby just soem rescailing while DGSS depends on n11/n12.  Didn't see that for s1=50, why would b-only depend?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#np.searchsorted?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#plt.legend?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}