{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2b23f073",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Program to find measure of expected significance as a function\n",
    "# of a cut value x_cut applied to measured variable x\n",
    "# by Monte Carlo simulation of the likelihood ratio statistic.\n",
    "# G. Cowan / RHUL Physics / December 2022\n",
    "\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b565a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Define pdfs and likelihood-ratio statisic\n",
    "s_tot = 10.\n",
    "b_tot = 100.\n",
    "ps = s_tot/(s_tot+b_tot)\n",
    "def f_s(x):\n",
    "    return 3.*(1.-x)**2\n",
    "def f_b(x):\n",
    "    return 3.*x**2\n",
    "def q(x):\n",
    "    return -2.*np.log(1. + (s_tot/b_tot)*f_s(x)/f_b(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5880ac4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "....................................................................................................\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Generate data under b and s+b hypotheses\n",
    "qb = []\n",
    "qsb = []\n",
    "numExp = 1000000\n",
    "np.random.seed(seed=1234567)        # fix random seed\n",
    "for i in range(numExp):\n",
    "    n = np.random.poisson(b_tot)       # first b only\n",
    "    r = np.random.uniform(0., 1., n)\n",
    "    xb = r**(1./3.)\n",
    "    qb.append(np.sum(q(xb)))\n",
    "    n = np.random.poisson(s_tot+b_tot) # then s+b\n",
    "    r1 = np.random.uniform(0., 1., n)\n",
    "    r2 = np.random.uniform(0., 1., n)\n",
    "    xsb = np.where(r2 < ps, 1. - r1**(1./3.), r1**(1./3.))\n",
    "    qsb.append(np.sum(q(xsb)))\n",
    "    if len(qsb)%(numExp/100) == 0:\n",
    "            print(\".\", end=\"\", flush=True)\n",
    "print(\"\\n\")\n",
    "qb = np.array(qb)\n",
    "qsb = np.array(qsb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bce0e84a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median[q|s+b]   = -42.664\n"
     ]
    }
   ],
   "source": [
    "# Make and analyse histograms of q for b and s+b hypotheses\n",
    "nBins = 400\n",
    "qMin = -200.\n",
    "qMax = 0.\n",
    "qbHist, bin_edges = np.histogram(qb, bins=nBins, range=(qMin,qMax), density=True)\n",
    "qsbHist, bin_edges = np.histogram(qsb, bins=nBins, range=(qMin,qMax), density=True)\n",
    "med_q_sb = np.median(qsb)\n",
    "print(\"median[q|s+b]   = {:.3f}\".format(med_q_sb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "68563526",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAECCAYAAADgnZClAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAk90lEQVR4nO3dC5QU5Zk38P8jgtcskiBLIuCgA0SQiGGCxKigjIoKXlluxssxYUg8msNZPYl4+aJHhXW/mBBPvJGEJV64eAdc/DAkwRsaEcUAKsriBXSzIxhdCURFn+88NVVD0XT39MxUdb1V9f+dU6e7q3u6367uqaff6yOqCiIiolL2KHkPERERAwUREbWEgYKIiMpioCAiorIYKIiIqCwGCiIiKmtP5JCIjAYw+ktf+tKkvn37Jl0cIorIunXrvMt+/folXZTUWbly5WZVPbDYfZLneRR1dXX6wgsvJF0MIorI8OHDvctly5YlXZTUEZGVqlpX7D42PRERUVm5bHoiomyaNm1a0kXIJAYKIsqMo48+OukiZBKbnogoM5YvX+5tFC3WKIgoM6688krvkp3Z0WKNgoiI8lGjEJFDAFwFoLOqjkm6PERErTVnNbCgaSoIzugHTBwIJzhRoxCRWSLSKCJrCvaPFJF1IrJeRK4o9xyqukFVvxd7YYmIYgoSU/8IPPdu02bXbZ8LnAgUAGYDGBneISIdANwK4BQA/QFMEJH+IjJQRB4t2LolV3QiomiChJl+QtNmgtpF0pwIFKr6JIAPCnYPAbDeryl8CmCe1cZUdbWqjirYGpFBl112Gfr3749JkyZh2LBh+Pzzz1uclfrWW28137YOvfPOO2+Xx3z66ac47rjjsGPHjtjKTZSUGTNmeFvaLPADggUIa26ybehBTTULF2oVTgSKEg4CsDF0e5O/rygR+YqI3AHgSBGZWuZxDSLygm3vv/8+XLVhwwY888wzeOWVVzBo0CCcffbZ6NDBKlmVW7VqFY488shd9nXq1AkjRozA/PnzIy4xUfLsf8W2NJmzuikgWGAI90lYH4UrtQqXA4UU2VdyYSpV3aKqP1DVQ1V1epnHzbT1TGw78MCi6185sbCZ1SDefvtt70T/m9/8BmeccUbz/W+88YZXe6irq8OPf/xj1NbWFn2el19+Ge+++y6OOuooHHLIIc1DBs8880zce++9VXs/RNWydOlSb0uTBaHO67CgVuECl0c9WQ2iZ+h2DwDvRbl6bKkTbOC6J4BXIq509D8Q+Omw8o+xlS8vuOAC1NTU4Pzzz0evXr2868aan2zfrbfeim9+85u49NJLMWDAgJI1Cgswf/7zn/H444/jmmuuwVNPPYXDDz8cK1asiPaNETnghhtu8C7r6+uR5tpEWND8lOQIKJdrFHYm6yMivUWkE4DxABZG8cSqukhVGzp37gxXrV69GkcccQQ2b96MAw44oHn/I4884vVbWJAwhx12GL7xjW/s9vfWB7Fly5bmCUhWHbfnMtaEZU1QH3/8cdXeDxFVXptwrfnJiRqFiMy1vlgAXUXEahI/VdXfisglAJbYuQ3ALFVdW80aRUu//OO0du1ar6bwySef4B//+Efz/pdeemmXNlhrXjrxxBN3+3vr27D3ZwHBvPjii17gCdjz7r333rG/DyJqe23C9icdJJypUajqBFX9qqp2VNUeFiT8/YtVta/f73BjhK/ndI3Cful37NgR++67L7p06eI1NwXB4itf+Qpee+0177o1Kd11111FaxQWQN58800vIGzduhXXXXcdpkyZ4t1nNQ3rn7HXICI3axNhSY9+ciJQVJvVKERk5kcffQQXrVmzxutHCJx00kl4+umnves23NWSLQ0cOBAPPfSQFziK1YwsUJx77rneappDhgzBj370IwwdOtS7709/+hNOPfXUKr4jIiqmXG3CpeYnZrhLQYY7a276+c9/jrvvvnuX/Rs3bsSYMWO8moWxkVCzZ89u7vguxYbaTp8+nekiKXPSlAp13ANNl/PHRPvYODLcOdFHQeXZENnjjz/ea4IKz6WwWkOxZqdybMKdDY9Nwz8SUWul5Xs9J9Q/kQa5DBSVdma75KKLLtpt36hRo7wtcOGFF+4yQqoY69y24bVEWbRo0SLvcvRo+xfPRv+EC3LZR+F6Z3ZbVRIoiLLs5ptv9ra0j3ZyTS4DBRFRUhakrDaR20Dh+qgnIsq2oSmqTeQ2UGS16YmIkIpmp7TJZaAgIkpbs9NzCU66y+WoJyLKpsK5RllpdjqjX1OgsECTRJNVLgNFGofHElHLevYMLzidHRMTXvMpl01P7KMgyiZLyMWkXNHLZaDISyrUKDCdKqXJ7bff7m0umpPSjmzDQJHhVKiVBgKbqFcK06kS5Xf+RICBIsOpUH/3u99h8ODB3npQxx57bJvKwnSqRPmcjR3GzuxyVk4B/rYq2hfvMggYPCP2VKiW0+Kmm27yagRWA/jwww/bVFymUyVypzaRVFrUXNYo8pAK1Zqptm/f7vVz2FLqhWtAWQ3BmrS+//3vY+HChd5125YssYSCTZhOlSg67a1NJJmXIpc1ioq18Mvf5VSolh3PEiDZapoNDQ1eQLj44oub7w9yWFhTkuWwsK0Q06lS2jzwgJ+4Yf1M4K05wI6twJ777/7A8P6aiUBtg/NLik9McIgsA4XjqVBtC1Kh2km5WCpUqzUUsn6MPn36YPz48d4JPxxsKhVOp/rZZ5956VR/8YtfePcxnSo5IwgKALoG+xqf2PUx3YbtvF54n932/z6OwLEgxZ3YAQaKFKVCra+v94aqnnbaaV4qVEtnWioV6o033ohnn30W++23n1cz+fWvf93qcoTTqVozlvVPMJ0quRQYdjnxdxuG2Uv+6l298ORhO2sNhSf+4O9tvyl8riBwRBgwhqa0EzvAVKg5TIUaBaZTpcSCQygwNPNP6vY/YILReW1+nfBr1LQvYESZxjTOlKhMhZpyUaZCjQLTqVLVFDtxR3DyLsme07bw6zaGmqZi7s9wVS4DRZ5ToUaB6VQpdjH8sm9XwDDBUHkGivwMj7X0unV1dZOQIeVmWBOlQtIBolTAMEuHN5XLylhheaIa8ZS0XAYKIkpB/4NrzTw1E5vK9/zkiju74xjxlMSkOwYKIko2SNiJN6L+h8WLFyM2tX6Zwn0X4f1VGPGUVF6KXM7MJiLHgsSQO4H6ZU1bO2oRwdyj2NQ2NJXRymus/NYkZe+lCqvFWnBIohmLNQoiSrYfwk66ETUx3Xbbbd5leBWC2GsXfyveyZ2FiXYB1iiIqHoBwn59269wCxLWzBRhkDD33Xeft1VFrV+7sIU+g07ujE20C7BGQUTV74twraM6qk5uU9uQmdFOmQwUInImgNPsqwjgVlV9POkyEeVeYV9EVgJEIHg/9h7997lgXUNmmp2canoSkVki0igiawr2jxSRdSKyXkSuKPccqvqIqtrcCJtQMC72QqfAtddei5/97GfNt23dprawFKv77LPPLivXJp2S1dafsvLYBMBg+XNytKkpq0EiYO8r6OD2h/pmpdnJqUABwNa5HhneISK2XsWtAE4B0B/ABBHpLyIDReTRgs1qEYGr/b+jAsuXL2/z3x566KFeIqM4tCUlqwUu2/+1r30tljKRu30RTqptaHqvjU9gxOe791ekmTOBQlWfBPBBwe4h9rVT1Q2q+imAeVabU9XVqjqqYLPaiLkJwGOq+iJSyn6Nf/3rX/dySNgqsraC69KlS/Gd73zHWzr8+eef9x53zz33YMiQId6v6smTJ3trQQUrx9o6TLbarKVVDdt//51r89t6TZYq1VaXnTlzZvNrWzKkSZMmeftt5Vr75V6MCylZydFmpnCAaOeQ19aw70Wi342aphVpGz6fnKlg4XofhXUFbQzd3mTJ2co8/lIA9QA6i0itqt5R+AARsW+s9621FKMtCVajDBs7dqw3/G7btm1Fl9q2X8a2WXOIre4aVumXeP369bj//vu9E/i3vvUtzJkzx1tq3LLRTZs2zVu5df78+XjmmWe8nBBWHsthbSf3efPmeSvOWoY6y4RnJ+liZs2ahS9/+cteILDXOOecc5oDwNy5c72lye29PvjggzjmmGN2+VtXUrKSQ7LeF9HK/orvfGFNUNk4Bq4HCimyr+S66Kp6CwDbSlJVC/Mzg2XG4ajevXt7OSeMnYBHjBhhQc7bZ7/6//CHP2DlypXeCd7Yyb5bt2744IMPcNZZZzVPOjr99NNLvsYtt9yChx9+uHnJcgsQ3bt391476IuwIGOvVxgoiqVkLbYgYTglq+UBt9pHmNUQLDHS1q1bvbIHr2vB5eSTT24xJStle15EawX9cZdffjmSMmd7A3p/Ogf/tBcyw/VAYTWInqHbPQC8V83VY8vVAOxkXO7+rl27trkavNdeO79le+yxR/Ntu24nTssjYideq1mEzZgxwwsoLbFyWXOWJTey92E1pyALXvi1gxN9IZdSslKCHBv2+uijjyYeKBasA6bYDzxt3QKCLnOmj6KEFQD6iEhvEbGzxHgAC6NYPVZVGzp37hxNKRNgNQzLD9zY2Ojdtl/jb7/9No477jivlmAnd2v2sRN0MR999BG6dOnincgttepzzz3XqtcvlpK1WG4Mq6VYlj1LyWpLorc3JavVPCwl65Qp9q9Izi3BkYGTYhRW7+Nnz7PjU2QiXto4U6MQkbnWJWA/xEXEahI/VdXfisglAJbYj1trVlfVtXnMR1HImn1uuOEGr7P5iy++8PoprL/AUpWOGzfO+7V/8MEHl+w8HjlyJO644w7v5G4d30GK00q5kpKVEsL+iKKa13c6qAENg3edW5HqY2RNGHndBg8erFSZN998UwcMGFD0vnfeeUeHDBnSfHvYsGHe46vl4IMP1vfff79qr5d7b9ypei+aNrvuEPvu2ZaUsfer9pqheu9fCo7V74dF+hq2RQ3AC6XOla43PcXCahQiMtOaX6gy1ldhx6vYhLukUrIGE+4+++wzr++GYpaCCXQ2t8a2JA0NT7QL5laknDNNT9WU1Qx3cerZs6c3MqqYpFKyBhPuKH+d1qU89thjcFJjuju2cxkoKF5MyZohDg19Tf2igW/NSe2xy2V9nU1PRBUK8i2kZBmO66+/3ttcXdoDKR0BlctAkYXhsUSxs5Oandws30JKhr7aRFTbnFOT7uGyuQwURNSKTuvgJEdllU19WmR12TTJZaBg0xNRdpqbXLGgpdSnKR4BlctAwaYnouw0N7lkaCU5KFLYV5HLQEFELQyBTWlzk60QYJuzaiZG0vxkTVzW1FUtHB5LlGfB8FeTgSGwtiR+EuZUmiPbjms7g4Q1bdlrWVNXtTLo5bJGwT4KooIkQ4Z9EvH1T0TY/GTBocWAFLFcBgr2URAVNH9Yf0QG+iSmTp3qbUkYWmmO7Iian6opl4GCKPeCTmvT3ZJCZoOtVGybM8NiMzIBj4GCKO9LhJ/w+6RLlK9mpxTWKtiZTZQXXLcp+WanCDu1qymXNQp2ZlMucSIdtVEuAwU7syl3/nhiLibS9ejRw9tSozEd/RRseiLKelPTjq3AByuRB/fccw9SoyY9y4/nskZBlAtBf0Q4SKR0xnUm1aZn9BMDBVHWBX0SE9X5X67tNWXKFG9zdmhsSkc/semJKOtBwvokcqLaqXEXtGVobApHP7FGQZTVfBI2woncGxqbQrkMFBweS7kYBmsjnNgnkQ6NbvdT5DJQcHgsZRbzSaRPjfv9FOyjIMrirOuc1iT69u3r3tLiGeinYKAgytLaTdZ5bUEipzWJmTNnpqcjO0Vy2fRElNkF/tjclJuO7OeqmOWOgYIozctyhIMEAwQaGhq8zfn5E+0U1GKCWk3c2PRElFZ/XbqzuYlBwvP666/notlp4sDqBQnDQEGU5jzXOe64TtrQHMyfyFzTk4gcJiJ3iMgDIvLDpMtDFHueay4Xni2N7s6lcCJQiMgsEWkUkTUF+0eKyDoRWS8iV5R7DlV9VVV/AGAsgLrYC02UhKAmwY7rbKlxey6FE4ECwGwAI8M7RKQDgFsBnAKgP4AJItJfRAaKyKMFWzf/b04H8DSAPyT2Tojixj6JkgYNGuRtqVPr9kqyTvRRqOqTIlJTsHuIVbRVdYPdEJF51nekqtMBjCrxPAsBLBSR/7SBCVUpPFE1+yWCpTmoqBkzZiC1atzNT+FEoCjB5jtuDN3eBOCoUg8WkeEAzgawF4DFZR5nn4D3KfTq1SvqMhPFN+O6Y2eu3+SAOVHNyE7RDG2XA4UU2aelHqyqtpZyi+spq6rV67y6XV1dXcnnI3JqMp0J1m+ikr773e/GnuluQY5mZKchUFgNomfotiXCfS+q1WMBjK6trY3i6YjiEfy6/PJgYM/9WZOowKZNdtqI39AcDY11PVCsANBHRHoDsDmQ422eSVSrxwJYVFdXNymK5yOKTc4SD5GbnBj1JCJzATwLoJ+IbBKR76nqDgCXAFgC4FUA96nq2ohej/koyF1MPJRvje6NfHIiUKjqBFX9qqp2VNUeqvpbf/9iVe2rqoeq6o0Rvh7zUZD7E+rYce2cOXGv8eTofIpWNz2JyH4A/qGqnyOl2EdBqZhQ59gQyTT49re/HevzL4i7I9vRkU8tBgoR2cPvHzgXwLcAfGJDUEXkfX8Y6kxVfQMpwj4KcnoYLCfUtdn06TbNKv5hsRNz1JFdadPTnwAcCmAqgO6q2lNVbSb0sbYkOoB/E5GmMWlE1Dbh7HTkpAU5HBbbmqanelX9rHCnqn4A4EHbRKQjUoRNT+SsIEMdtck555zjXT74oJ2aoje0WrWJoEPbkZplizWKIEiIyBsi8rCIXCciZ4tI81m2WCBxGTuzyRl2MpgjO5ucuNBfu2zZssXbUq3GvQ7t1ox6eshfUuOvAE4C8BcReUdEnhWRO2MsI1E+Zl6zJkHhBQJTOurpeFW1hfo8IvIfAM4C8CsARyBF2PRETqQxDTLUme71rEmQs1pTo/i7iDQHBFX9sy0BrqqbVNVWa00NNj2RU2lMJypwwu+TLhFRJDUKG0p6l4jY7GibMnoYgO2t+HsioliNGDEiXSvGZi1QqKplmTsGwJkAjrQWVgA/jbd4RBmaI2F9ENa8FF6egf0SkbrmmmuyMzS20Z2RTy02PYlI83LfqvqFqj6kqteo6gxV3VL4mDTgWk+UyByJYBQLZ1+n0tBqTrRzbORTJX0UfxSRS0Vklyw/ItJJRE4Qkd8BuAApwj4KSoQFCxsKa4v9cfZ1LE455RRvy8vIp+febWoWcyFQ2PIctq6TzaF4T0ReEZE3/f0TAPxCVS3nNRFVgov9xWb79u3elqqFANsoaAYLmsWS7qM42n59i8j3LXsogAPt81DVD+MvHlFG2K/DPfwFDDjCKVUWOLp0hzWDVSNIVBooltikOgD/DOB8AC8DiCQvBFGuMECkTp4XAmxVoFDVy0TkED8ftWWbOx3AABH5FMAaVR2HlOGEO6raSCcTLM9BqeNqbcLJ4bGqukFEbHHA14N9IrI/gMORQlxmnKq6NAdVzahRoyJ/zqE5r020dh7F6wW3t/rLjBNRWHj4a3CbnddVcfnllyddhExqdYY7Iqpg/abw8FcOg6WUT7pzImc2USaEF/ljDSIRw4cP97ZMqHFn0h0DBVHUOJmOMrbcOJueiKIc4WRYm6CMyWWNgms9UWz5rlmboAzKZaDgWk8UaQd2ECSIMopNT0RRdWBT4saOHZt0ETKJgYIoKtbsxP6JRF188cWRPVfekxWFMVAQtaXjesdWYE9bnMDH3BJO2LZtm3e57777tvu5uHzHTgwURG2ZUBfGDmxnnHrqqd7lsmW2NF37cfmOHHdmE7VJqSDB5iaqxuzsBDFQELVWeBJU/TLWJijzs7MZKIiIXFXrxuzsTAUKEdlPRFaKSPRrDROFda9v2ohy0PzkRKAQkVki0igiawr2jxSRdSKyXkSuqOCpfgLgvvhKSuQ3B1i2Omasc86FF17obZlSk3zzkyujnmYD+BWAu4IdItIBwK0ATgSwCcAKEVkIwPZPL/j7iwB8A8ArAPaufvEpk0Ngw/+gdj1oAmCfhLMyFySC71vCfRROBApVfVJEagp2D7F/WcuuZzdEZJ4NaVZVCxK7NS2JyPEA9gPQH8B2EVmsql8UeZz9l3v/6b169YrxXVEm1m4ydt22jp2BLoOSLBm1YPPmzd5l165dky5KpjgRKEqw+ZAbQ7etVnFUqQer6lV2KSL2k2JzsSDhP84a+rzGvrq6Oo2h3JQV4WAR1CY4FNZpY8aMiXQeBbkfKKTIvhZP7Ko6u8UnFhkNYHRtbW2bC0c5YkHChsFSLtjSHVP/2HSdy3c41JldgtUgeoZu9wDwXhRPzNVjqaLg4MCwRKq+YOkOw+U73A8UKwD0EZHeItIJwHgA1pndbsxHQUSV4PIdDgUKEZkL4FkA/URkk4h8T1V3ALgEwBIAr9qwV1VdG8XrsUZBRJSyPgpVnVBi/2IAtkWKfRRUdmisdWKH13BiB3Zq/PCHP0y6CJnkRKCoNqtRAFhUV1c3KemykGNzJ4KRThYcbPw650ykyrhx4yLJQWGOCfeQ5lwuAwXRbixI/G3VzpoEA0QqbdzYNKK+Z8+e7erInn4C+yeQ90DBpicqyibTcRhsqp133nntnkfBHBSOdmZXGzuzaZdkRLZRrlmT08G/3NnsRLvKZY2CqGgyIs6byC3OnSgvlzUKzqMgT+GyzRzdlHvWgc1mp93lMlCw6SnHgWHp8J0BIrwiJ/Ne5571Tdx7dtKlcBObnig/wsNfGRQy6bLLLku6CJnEQEH5U5gtLDy5jlJt9Ggb0EhRy2Wg4PBYwvOTm/JLcGXYTFm3rqlXul8/9khHiX0UlC/hkU02b4I1iUyZPHmyt2VSY3J5s3NZo6Cc617fdMmc15QWNRObAoX1syXQv5bLGgXleLG/IEAwSFCR9Z2cVduQ6Dwf1igof4v9ERWZbMeJdqXlskbBCXc5WoojHCQ4V4JK1Ca4vlN5uaxRcJnxjC/FQbl19dVXt+rxrE1UJpeBgjIsqFUU64MIahaUWfX1/kCFVmBtomUMFJTNmoUt1cHAkDurVq3yLgcNGpR0UTKFgYKyyZIQBYbcCbxzf5KloSqZMmVKu/NRpGIuRW11+9py2ZlNOWCT6QL2T8UhsZR2NRN3X8yyShgoiIjSoDa5uRS5DBQcHktEVLlcBgqu9ZQRTGNKVBXszKb04twJKjBt2rSki5BJDBSUfgmtqEnuOfroo9s0K5vKy2XTE2VMAqNAyE3Lly/3tkpwVnblWKOg9CzsZ8MDi40f58Q68l155ZWtmkfBWdmVYaCg7Oa6DvJOEFG7MFBQetmY8nK1CU6wI4oEAwWltxmKiKoiM53ZIjJcRJ4SkTvsetLloZibocId2IWzVe02AwlRtmoUIjILwCjrllTVw0P7RwL4JYAOAH6jqv9W5mkUwFYAewPYVJ2SU+IsKNQv27labHCbcmnGjBlJFyGTXKlRzAZgQaGZiFhwuBXAKQD6A5ggIv1FZKCIPFqwdQPwlKraY38C4Lrk3goVbS6yE3nhfIdS+1u6j6gEW168kiXGU5En2yFO1ChU9UkRqSnYPcROF6q6wW6IyDwb8qyq0/3aRyl/A7BXqTtFxIbNeENnevXqFdVboLaMWio3mil8H1GFli5dWlECI86hSGGgKMHmS24M3bbmpKNKPVhEzgZwMoADAPyq1ONU1X6iej9T6+rqrLmKiDLihhtuqDjTHedQZCNQSJF9JU/sqvoQgIcqemKR0QBG19bWtquARER54EofRTFWg+gZut0DwHtRPDFXj02INSWV63doqV/CstaFM9cREfJeo1gBoI+I9AZg3U7jAUQy5pE1igSVm2HdUr9EkLWu2GOC4bAcFkuUzRqFiMwF8CyAfiKySUS+p6o7AFwCYAmAVwHcp6pro3g91ihSKBj2aluxLF8WeOy+KucSJsoDJ2oUqjqhxP7FAGyLFGsURNl05513Jl2ETHKiRlFtrFEQZVO/fv28jTJYo6Acsg7rUv0RlcyfCB6TULJ5ctOiRYu8y9GjrdGgOCYsar1c1iis6UlEZn700UdJFyW/2pNsiB3WVMLNN9/sbeVwsl3r5TJQsOnJEVYbaEuNwDqsWZOgduBku9bJZaAgIqLK5TJQsOmJiKhyuQwUbHoiIqocRz0RUWbcfffdSRchkxgoKN6UpZUOdf1/dcAHK3e/z9Z2+oxNhFSZnj3Dy8MB5/rLhJ7Wp2m0UyZGOjU+0fT/5a9CYEN9bchvnJ3zuWx6Yh9FFRQLEsVGKgVDXQuDRPf6Xdd3KvZ3THlKBebPn+9tgac3Nm0WJOyEGgyNTa2aibsMLw8CX9zvK5eBgn0UCa7VVBgswkNdw/ed8HtgopZOa8q1naiI22+/3dsyq3bXoeFWi6jGxMFcBgoiIqocAwUREZXFQEFERGXlMlCwM5uIqHJ75rUz2xaarKurm5R0WYgoOg888EDSRcikXAYKIsqmrl27Jl2ETMpl0xMRZdPs2bO9jaLFQEFEmcFAEQ8GCiIiKouBIkts/Zelw5su2/P41j5Psectt8aTrd9kG1HEGv8OvPI+MGouMC7Ur23Ldxi7z7bMrPdUJXvmdXispdWtra1FZtdXqmRpi1KPb+3zFHveUsJrMxVbpylY44moDTZvA/73E2B1Y/H7+x/YdJnqxQFrJjb9f9r/WZWWsMlloODw2ATZF7vcl9vWeCKKyfwxSL/ahvblnG+DXAYKIsqmY69cjOffS7oU2cNAQUSZsede+2KPTkmXInvYmU1EmbF+yW34+Onbki5G5jBQEFFmbFx+H7a9dF/SxcgcBgoiIiqLgYKIiMpioCAiorIYKIiIqCxRVeSViNhk/rfb+Oe2nvFmuMO18hiWqTIsUzrL1NWx8rS3TAerqj93fVe5DhTtISIvqGodHOFaeQzLVBmWKZ1lEsfKE2eZ2PRERERlMVAQEVFZDBRtV701ftNZHsMyVYZlSmeZZsI9sZSJfRRERFQWaxRERFQWA0ULROT/ishrIvIXEXlYRA4I3TdVRNaLyDoROTm0f7CIrPbvu0VEJOIy/YuIrBWRL0SkeYSDiJwrIqtCm90/yL9vmV/O4L5uVSpTjYhsD73uHdU4TmXKc6KIrPRf1y5PCN2XyDFK8rtUUIb5off+ll229BnGTUSuFZF3Q699akvHLKlzgiR4nPzXH+kfCzsmV0T65Nb0xK30BuAkW73Yv36Tbf71/gBeBrAXgN4A/gtAB/++5wF82z47AI8BOCXiMh0GwHJ0LQNQV+IxAwFsCN0u+dg4y2T5uACsKfE3sR2nMuU5EsDX/OuHA3jXgWOU2HepTFlvBvB/WvoMq1COawFcXmR/yWOW4DmhJsHj1ME/BocA6OQfm/5RPT9rFC1Q1cdVdYd/8zkAPfzrZwCYp6qfqOqbtsIxgCEi8lUA/6Sqz2rTJ3gXgDMjLtOrqrquhYdNADAXVVJhmZrFfZxKlUdVX1LVILXNWgB7i4idbJI8Rol9l4rxay1jq/n9aYOixyzhc0KS7L2vV1X7cfipHRv/GEWCgaJ1LvJ/1ZmDbFXj0H2b/H0H+dcL91fbuCL/6P/hV4mvibMJo4jeIvKSiDwhIsf6+1w4TucAsMDxScLHyLXvkn1G/6Oqb7TwGVbLJX4zzywR6dLCMUvynJDkcYr1eDDDXdMvqKUAuhe56ypVXeA/5ioA9ivi3uDPijxey+yPvExl/vYoANtUdU1o97mqam29XwLwIIDz/F+ocZfpvwH0UtUt1t4O4BERGRDFcWrnMRrgNxtYM0LSxyjW71IbyldYGy36Garq/7anLJWUCcDtAK733/f1fpPYRXEcmwjOCf8d53FqqchxHg8GiqaqZH25+0XkAgCjAIzwmwCCiN0z9DCrfr7n7+9RZH+kZWrB+MLahJ0A/cuPRWSOX1W9K+4y+b/WvV/sqmqdx9aO2jeK49TWYyQi9loPAzhfVf8r6WMU93epNeUTETsnnA1gcAWf4QvtKUulZQqV7dcAHm3hmCV2TtCYj1MLYj0ebHqqYCQBgJ8AOF1Vt4XuWmgnZGvfFhHrTOtjHY+qar8q7EQz1G+6OB/AgiqW1z7Tf/HbKIN9e4pIV/96R/8LvqZK5TlQRDr41w/xj9OGpI6TP0LlPwFMVdVnXDhGjn2X7AT5mqpuaukzjLkcwWtbP03grNBnUvSYJXlOkASPE4AV9np2LESkk/9j0Y5RNJLooU/T5neSWdvfKn+7I3SfVT3tV8O68GgUG9Hif6Htvl8FExsjLJP9w2zyf738D4AlofuGWwdbweP3A7ASwF/8DtxfRj1CpFSZ/H6Atf4ojBcBjK7GcSpTnqsB/D30edrWLcljlOR3qUgZZwP4QcG+kp9h3BuAuwGs9j8XO/F9taVjltQ5AQkeJ//1bejw6/4xuSrK5+bMbCIiKotNT0REVBYDBRERlcVAQUREZTFQEBFRWQwURERUFgMFERGVxUBBRERlMVAQEVFZDBRECRGRPn6ypBdE5N8t4UzSZSIqhoGCKAH+mkC24OC/qqot07GPv/wDkXO4eixRMiwB0SuqamsCmVcBfJhwmYiKYo2CKBmWktXLSe07wl9Mjsg5DBREydgC4OuhRFPn+yukEjmHq8cSJcDPfWF5MfYFsNiy6/nZ0b5IumxEhRgoiBImIpaZ7AFVtZoFkXPY9ESUPOufYLMTOYs1CiIiKos1CiIiKouBgoiIymKgICKishgoiIioLAYKIiIqi4GCiIjKYqAgIqKyGCiIiAjl/H++EB+tYnVIswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot histograms of q\n",
    "binLo, binHi = bin_edges[:-1], bin_edges[1:]\n",
    "xPlot = np.array([binLo, binHi]).T.flatten()\n",
    "ybPlot = np.array([qbHist, qbHist]).T.flatten()\n",
    "ysbPlot = np.array([qsbHist, qsbHist]).T.flatten()\n",
    "fig, ax = plt.subplots(1,1)\n",
    "plt.yscale(\"log\")\n",
    "plt.gcf().subplots_adjust(bottom=0.15)\n",
    "plt.gcf().subplots_adjust(left=0.15)\n",
    "plt.xlabel(r'$q$', labelpad=5)\n",
    "plt.ylabel(r'$f(q)$', labelpad=10)\n",
    "plt.plot(xPlot, ybPlot, label=r'$f(q|b)$', color='dodgerblue')\n",
    "plt.plot(xPlot, ysbPlot, label=r'$f(q|s+b)$', color='orange')\n",
    "ax.axvline(med_q_sb, color=\"black\", linestyle=\"dashed\", label = r'median$[q|s+b]$')\n",
    "ax.legend(loc='upper left', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fb2a1953",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add code to calculate the p-value of the b-only hypothesis\n",
    "# for the median q from data generated according to s+b.\n",
    "# Find the corresponding significance Z."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edb93068",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}