{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c752a0ff",
   "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",
    "# Full version containing exercise solutions\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": "42cfb911",
   "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": "274ca242",
   "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": "da9f66cd",
   "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": "ffe5f88e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAECCAYAAADgnZClAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnnklEQVR4nO3df5QU9Znv8fcjgoLJGiJyTQQEM0hkREcZwZiYmQgaNBCMuv7AqJxEh40n2cO5enajiTfxaGTdu27YnBh1TFyiEcX4I4Ix1wTNmCj4A+IQkEiYKAZwdxEwLKwTFHzuH9011DTdPT0z1V3VXZ/XOX2Yrqquerqm6We+v83dERERKeSAuAMQEZFkU6IQEZGilChERKQoJQoRESlKiUJERIpSohARkaIOjDuAOJjZDGDGBz/4wSuPOeaYuMMRkYisW7cOgHHjxsUcSfVZuXLlVnc/PN8+S/M4isbGRl+xYkXcYYhIRJqbmwFoa2uLNY5qZGYr3b0x3z5VPYmISFGprHoSkdp08803xx1CTVKiEJGaceqpp8YdQk1S1ZOI1Ixly5axbNmyuMOoOSpRiEjNuO666wA1ZkdNJQoRESmqZkoUZnY08A3gUHc/P+54RER6a+FqeCwzFISZ42DWhHjjCSSiRGFmd5vZFjNbk7N9mpmtM7MOM/t6sXO4+2vu/uXyRioiUh4LV8O1T8PzmzOPa5/ObEuCRCQKYAEwLbzBzAYAtwFnAeOBi81svJlNMLPHcx7DKx+yiEg0giQBMO/0zAP2lS7ilohE4e6/AbbnbJ4EdGRLCu8CDwAz3X21u0/PeWypeNAVcPXVVzN+/HiuvPJKmpqa2Lt3b9Hjm5ub2bBhQ9fztrY2Lr300m7HvPvuu3z6059mz5495QhZJFbz589n/vz5cYfRa0FCmHd6prpp1gQ45chMySIJpYpEJIoCjgQ2hp5vym7Ly8wOM7M7gBPN7Noix7WY2QozW/HWW29FF23EXnvtNZ577jnWrl1LQ0MD5557LgMGDOjVOdrb2znxxBO7bRs0aBBTpkxh0aJFUYYrkggNDQ00NDTEHUavLFydSQinHNm9TWJmdrqqJJQqkpwoLM+2ghNTufs2d/87d/+Yu88rclyruze6e+Phh+ed/yp269ato6mpiTfeeIMTTzyRH/7wh8ycObNr//r162lubqaxsZF/+Id/oK6uLu95Vq1axebNm5k8eTJHH310V5fBc845h/vuu68Sb0WkopYuXcrSpUvjDqNXwo3XYUGpIgmS3OtpEzAy9HwE8GYUJw5mjy30BRu44RlYG3GhY/zh8K2m4seMGzeOyy+/nNGjR3PZZZcxatQoRo8eDcDevXu57LLLuO222zjppJP42te+Rn19fd7ztLe3M3PmTF544QV++ctfcv311/Pb3/6W4447jpdeeinaNyaSADfddBMAU6dOjTmS0hQqTYQF1U9x9oBKconiJWCsmY0xs0HARcDiKE7s7kvcveXQQw+N4nRlsXr1ak444QS2bt3Khz70oa7tP/vZzxg/fjwnnXQSAMceeyzHH3/8fq/fs2cP27Zt6xqA1NDQwNatWwEYMGAAgwYNYufOneV/IyJSUKHSRCAp1U+JKFGY2f1AMzDMzDYB33L3H5nZV4EngQHA3e7+SkTXK6lE0dNf/uX0yiuvUF9fz+7du/nrX//atf3ll1/uVge7atUqzjjjjP1ev3btWurq6hg0aBAAv/vd7zjhhBO69u/evZuDDz64fG9ARIoqpTQxa0L8SQISUqJw94vd/SPuPtDdR7j7j7Lbn3D3Y7LtDt+J8HqJLlHs3LmTgQMHMmTIEIYOHcrevXu7ksVhhx3Gq6++CsALL7zAPffck7dEsWrVKl5//XV2797Nrl27uOGGG5g7dy4A27Zt4/DDD2fgwIEVe08i0l1PpYmwuHs/JSJRVJqZzTCz1h07dsQdSl5r1qzhuOOO63p+5pln8uyzzwJw6aWXsmLFCiZMmMAjjzzCYYcdlrcxe9WqVVxyySWceuqpTJo0ib//+7/nlFNOAeDXv/41Z599dmXejIgUVKw0EUhC9ZNWuKuCFe5efvll/vVf/5V777232/aNGzdy/vnn88ILLwCZcRQLFizoavgu5Nxzz2XevHlaLlJqTjUthXrhQ5l/F5Uw4VBvju2rYivcJaKNQoo78cQT+cxnPsPevXu7jaVYtWpV3mqnYt59913OOeecqviPJNJb1fK5DrdPVINUJopSG7OT5Etf+tJ+26ZPn8706dO7ns+ePbtbD6l8Bg0axGWXXRZ1eCKJsGTJEgBmzJgRcyTF9aZ9IglS2UaR9MbsviolUYjUsltvvZVbb7017jCKKqW3U9KkMlGIiMSl2koTkNJEkfReTyJS26qpNAEpTRS1WvUkIskWVDtVm1QmChGROPSn2inOQXep7PUkIrUpd6xREvWl2mnmuEyieGxdPFVWqUwU1dg9VkR6NnLkyJ4PqkJxz/mUyqontVGI1KZFixZpUa4ySGWiqBb9XQo1ClpOVarJ7bffzu233x53GHlVa0M2KFEkVhRLoZaira2N2bNnF9yv5VRFolGN4ycCShQJFNVSqD/+8Y+ZOHEixx9/PKeddlqfYtFyqiL9V42jscPUmF3Myrnwdnu0Fx/aABPnFz0kiqVQd+7cyS233EJ7ezuDBg3iL3/5S5/C1XKqIv0XVWkirmVRU1miqIbG7P4uhTpgwAA6Ozu5+uqrWbFixX5zQE2ePJmGhgauuOIKFi9eTENDAw0NDTz55JNdx2g5VZHo9Lc0Eee6FKksUZSsh7/8y6m/S6EOGTKENWvWsGTJElpaWrjiiiu46qqruvYHa1i0tbWxYMECFixYsN85tJyqVJuHHsou3NDRChsWwp5dcOAH9j8wvH30LKhrKVtMUU0pHmcXWSWKBAovhTpkyJCupVAPPvjgvEuhXn311fudY/369YwdO5aLLrqItWvXdks2pQovp/ree+9xww038N3vfhfQcqqSIEFSAIYF27Y80/2Y4U37fs7dt+WZrtcDkSeOam7EDihRJFChpVCnTp3KpZdeyuc+9zkmTJjA2WefXXAp1O985zssX76cQw45hPr6eu66665exxFeTrWzs5Prr79ey6lK/EKJAdj3xT+8iQVP/icAsz/btK/UkPvFH7x+9KzM89xzBYkjwoRRrY3YAS2FmsKlUKOg5VSlosLJIZQYumS/1JubmwG6euf1+Trha/QzYUS5jGk5l0TVUqhVLsqlUKOg5VSlYvJ9cUfw5V1QXUvmEb5uuGqqzO0ZSZXKRFGNcz1FtRRqFLScqpRdGf6y75XchAH7usorUaSDuy8BljQ2Nl4ZdyxRKjbCWqQqxJ0gcgUJA2BpcyaujtaS44mqx1PcUpkoRCRBCrU/JK2aZ/SsTHwvzim5sbscPZ7iGHSnRCEi8elozXzxQiTtD0888USEweUIYgq3XYS3FxBlj6e41qVI5chsEUmAcJKYdCdMbcs8+lGKCMYelU1dSybGSXdmnr84J1Ml1dG636HlmC121oR4qrFUohCRyspth5h0Z2RVTD/4wQ8Aus1CUBbh0kWBRu5aGGgXUIlCRCqjozXz1/eLczJJYnhTpEkC4MEHH+TBBx+M7HxFBaWLoQ37GrlzVPtAu4BKFCJSfrltEUlrqO6PcCM3QF1LzfR2CtRUojCzc4DPAcOB29z9l/FGJCL7tUXUSoIIBO/nxTld7/OxdZlttVDtBAmqejKzu81si5mtydk+zczWmVmHmX292Dnc/WfufiUwG7iwjOFWjW9/+9v8y7/8S9fzU089tU/n2bBhA4MHD+42c20hlVqStbOzk4aGBgYNGtQ1/bkkSLiqCWozSQTqWvY1cGe7+tZKtRMkKFEAC4Bp4Q1mNgC4DTgLGA9cbGbjzWyCmT2e8xgeeuk3s6+THMuWLevzaz/2sY/R3t4eXTAhfVmSdfDgwbS3t/PRj360LDFJH1WgLSKR6loy73XLM0zZu397RTVLTKJw998A23M2TwI63P01d38XeACY6e6r3X16zmOLZdwC/MLdf1fp9xCVDRs28PGPf5wrrriC4447jksuuYSlS5fyyU9+krFjx/Liiy8C8JOf/IRJkybR0NDAnDlz2Lt3L5CZOXbcuHFMnTqVdeu6T2D/gQ/sm5v/nHPOYeLEidTX19Pa2tp17WOPPZYrr7yS+vp6zjzzTDo7O/PGmYQlWSVhgmqmcILoZ5fX3mhra4v3s5GdkbZl75yaShZJb6M4EtgYer4JmFzk+K8BU4FDzazO3e/IPcDMWoAWgFGjRvUYQDAbZdgFF1zAVVddxTvvvJN3qu3Zs2cze/Zstm7dyvnnd5/msdQPcUdHBz/96U9pbW3l5JNPZuHChTz77LMsXryYm2++mXnz5rFo0SKee+45Bg4cyFVXXcV9991HfX09DzzwAC+//DJ79uzhpJNOYuLEiXmvcffdd/PhD3+Yzs5OTj75ZM477zwgkwDuv/9+7rrrLi644AIefvhhPvWpT3V7bVKWZJUEqfW2iFKE2is++f5Csl81VS/picLybCs4L7q7fw/4XrETunsr0AqZacb7FV0ZjRkzhgkTMhWc9fX1TJkyBTNjwoQJbNiwgaeeeoqVK1dy8sknA5n6+uHDh7N9+3a+8IUvdA06+vznP1/wGt/73vd49NFHgcyU5evXr+eII45gzJgxXW0REydOZMOGDfslinxLsuabkDC8JOvll19OY2P3WYwnT57M7t272bVrF9u3b++67i233MJnP/tZoPiSrJIAZRwX0VtBe9w111wTy/UBFna2MObdhfzNQbGFELmkJ4pNwMjQ8xHAm/09aW9mjy1WAhgyZEjR/cOGDetzMfigg/Z9yg444ICu5wcccAB79uzB3bn88suZN29et9fNnz8fs3z5tbu2tjaWLl3K8uXLGTJkCM3NzV2r4IWvHXzR50rSkqwSo4R1e3388ceBeBPFY+tgLlDvvZtAMMkS00ZRwEvAWDMbY2aDgIuAxf09qbsvcfeWQw89tN8BxmXKlCk89NBDbNmyBYDt27fzxhtv8OlPf5pHH32Uzs5Odu7cyZIlS/K+fseOHQwdOpQhQ4bw6quv8vzzz/fq+vmWZM23Nsb69es55JBDuOiii5g+fXq/l2TdtWsXN9xwA3Pnzu31eSRi+abgqIEvxSisHpxdPe/FOXkH4lWbxJQozOx+oBkYZmabgG+5+4/M7KvAk8AA4G53fyWCa1XdehS5xo8fz0033cSZZ57J+++/z8CBA7nttts45ZRTuPDCC2loaOCoo44q2Hg8bdo07rjjDo4//njGjRvXtcRpqZKyJKvERO0ReXXN73RkCy0T6Ta2oqrvkbun9jFx4kSX0rz++uteX1+fd9+f//xnnzRpUtfzpqYmf/311ysUmftRRx3lb731VsWul3rr73S/j8xj/Z1xR9NNU1OTNzU1xXb9C37qPmq++32/z24I7tWvoovpgp9mHlEDVniB78qkVz2VhZnNMLPWHTt2xB1K1RgwYAA7duzIO+AuriVZgwF37733HgcckMqPcmVVwQC6wYMHM3jw4Fhj6DbQLhhbUeUSU/VUSV6jK9yV08iRI9m4cWPefXEtyRoMuJMKSFijdSG/+MUv4g4hv16ujJc0qUwUUl5akrWGJKjra9UKJg3csLBq710qy+uqehIpUbDeQpVMw3HjjTdy4403xh1Gd6GpPaq1B1QqE4XXQPdYkbLraM18uQ1tqJqur0899RRPPfVU3GHsb3R1d5dNZaIQkSJyG62DLzkpqujSp3lml60mqUwUqnoSKaLKqpuSoselT6u4B1QqE4WqnkQKqMLqpiQpaQ2KKmyrSGWiEJE8wl1gq7S66bDDDuOwww6LO4zCgvvaz+qn5zdnqroqRd1jRdIs6P4KNdEF9uGHH47luiWvkV3X0u8kMXNc5lqPravcCnqpLFGojUKE7osMgdok+qHH9olc/ah+mjWhhIQUsVQmCrVRiND9L9upbTXRJnHttddy7bXXxnLtktfIjqj6qZJSmShEUi9otAY4Ymq8sURo+fLlLF++vKLXLNotNp8qHICnRCGSNrlThJ/+q3jjqXK9rnaCqitVqDFbJC00b1PZlFztFIigUbuSUlmiUGO2pJIG0kkfpTJRqDFbUufpM1IxkG7EiBGMGDEi7jBKVyXtFKp6EqlVQVXTnl2wfWXc0VTET37yk7hDKF0VTT+eyhKFSCoE7RHhJFGlI65rUhX1flKiEKl1QZvELE/8X679NXfuXObOnVux6/W6a2yuKun9pKonkVo2vCnTJpESlV4at09dY8OqpPeTShQitSZYT+Lt9rgjSYVed42tQqlMFOoeKzUt6AY7tEFtEtUi4e0UqUwU6h4rNUvrSVSfKminUBuFSC3IHXWd0pLEMcccU7FrlTy1eE+qoJ1CiUKk2oXnbhrelEkSKS1JtLZWrvqm3w3ZVSSVVU8iNSN3gj9VN1VUnA3ZlVzlTolCpFo9fUb3JKEEQUtLCy0t5b8P/R4/0U9BKSYo1ZSbqp5EqtV/Ls38O7xJSSLrj3/8Y0WuE3e106wJlUsSoEQhUn3C61xDahuu45aG8ROBmql6MrNjzewOM3vIzL4SdzwiZRFe51rThdeWBI+lSESiMLO7zWyLma3J2T7NzNaZWYeZfb3YOdz9D+7+d8AFQGM54xWJTVCSUMN1bUn4WIpEJApgATAtvMHMBgC3AWcB44GLzWy8mU0ws8dzHsOzr/k88CzwVGXDF6kgtUkU1NDQQENDQ9xh9F7CZ5JNRBuFu//GzEbnbJ4EdLj7awBm9gAw093nAdMLnGcxsNjMfg4kMzWL9EXQLhFMzSF5zZ8/P+4Q+i7B61MkIlEUcCSwMfR8EzC50MFm1gycCxwEPFHkuBagBWDUqFERhClSRuER1wMP1fxNCRDZiOxcCR6hneREYXm2eaGD3b0NaOvppO7eCrQCNDY2FjyfSOzCg+lg3/xNUtAXv/hFoLwr3cXdNTYOSU4Um4CRoecjgDejOLGZzQBm1NXVRXE6kfII/rr88EQ48AMqSZRg06ZNFblOmrrGQrITxUvAWDMbA2wGLgIi+Z/i7kuAJY2NjVdGcT6RsknZwkOSTIno9WRm9wPLgXFmtsnMvuzue4CvAk8CfwAedPdXIrqe1qOQ5NLCQ+mWwJ5PiUgU7n6xu3/E3Qe6+wh3/1F2+xPufoy7f8zdvxPh9bQehSRTeECdGq4Tp+xzPCV0PEWvq57M7BDgr+6+twzxVITaKCSxwgPqEtZFshp84hOfKOv5y96QndCeTz0mCjM7gEz7wCXAycBu4CAze4tMN9RWd19f1igjpjYKSZxwN1gNqOuzefPmle3c4W6xaWrIhtKqnn4NfAy4FjjC3Ue6+3DgNOB54J/M7ItljFGk9oVXp5NESmO32EApVU9T3f293I3uvh14GHjYzAZGHlkZqepJEitYoU765LzzzgPg4YcfLsv5K1aaCBq0E1Ky7LFEESQJM1tvZo+a2Q1mdq6Z1eUeUy3UmC2J0dEKC21flZMm+uuXbdu2sW3btrjD6J8ENmj3ptfTI2Sm1PhP4Ezg92b2ZzNbbmZ3liU6kVqWO/JaJQmBfRMEJkhvej19xt0nBU/M7N+BLwDfB06IOrByUtWTxO7pM/atUAdwxFSVJCSxelOi+B8z60oI7v4CcJa7b3L3n0cfWvmo6kliF17GdJbD6b+KNx6RInpTorgSuMfMXgHagWOBznIEJSLSF1OmTCnLecs2Y2yVKDlRuHuHmX0KOAc4EegAvlWmuERqRzBGYvSsTPVSeHoGtUtE6vrrry/LeWPpGpugnk89Vj2ZWdd03+7+vrs/4u7Xu/t8d9+We0w10FxPUlHBGImgF4tGX1elig60S1jPp1LaKJ42s6+ZWbdVfsxskJmdbmY/Bi4vT3jloTYKicWWZzJdYd9u1+jrMjnrrLM466yz4g6j/0rs+fT85ky1WLmVkijWA3uBR83sTTNba2avZ7dfDHzX3ReUMUaR2qLJ/sqms7OTzs5om07LPhFgHwXVYEG1WDmV0kZxqru3mNkVwCjgcKDT3f9S1shEasnwJjggO4GBejhVlaRO3TFrQmWSBJSWKJ40s+XA/wIuA1YBkawLIZIqShBVJ80TAYb1mCjc/WozO5rMetRjgM8D9Wb2LrDG3S8sb4jR04A7KbugpxPsm55Dqk5SSxOVVlL3WHd/zcymuvsfg21m9gHguLJFVkaaZlzKKndqDqmY6dOnR37OtJcmoHfjKP6Y83wXmWnGRSQs3P01eK7G64q45ppr4g6hJvV6hTsRKSA8f1O4+6u6wUpfJWTQXSLWzBapCeFJ/lSCiEVzczPNzc1xhxGNBA26U6IQiZoG00kUEjTduKqeRPor3MMJVJqQmpPKEoXmepJIhde7VmlCalAqE4XmepLIPH3GviQhUqNU9STSH+EGbIndBRdcEHcINUmJQiQqw5vUPhGzq666KrJzpX2xojAlCpHeCBqu9+yCAz+wb7vWlkiEd955B4AhQ4b0+1yavmMfJQqRUoUH1IWpATsxzj77bADa2toiOZ+m78hIZWO2SJ8UShKqbpJyCkZnx0iJQqS3woOgprapNCHlk5DR2UoUIiJJlZDR2TWVKMzsEDNbaWbRzzUsEnbE1MxDpBJirn5KRKIws7vNbIuZrcnZPs3M1plZh5l9vYRT/SPwYHmiFMkaPSuzWp1WrEuc2bNnM3v27LjDiFYCqp+S0utpAfB94J5gg5kNAG4DzgA2AS+Z2WJgADAv5/VfAo4H1gIHVyBeqWVBF9jwf9DRs/ZVAahNIrFqLklA5vMWcxtFIhKFu//GzEbnbJ4EdLj7awBm9gAw093nAftVLZnZZ4BDgPFAp5k94e7v5zmuBWgBGDVqVKTvQ2pEeO4myPy85RkYeCgMbYgtLOnZ1q1bARg2bFjMkdSWRCSKAo4ENoaebwImFzrY3b8BYGazga35kkT2uFagFaCxsdGjClZqUDhZBKUJdYVNtPPPPx+IbhyFZCQ5UViebT1+sbv7gh5PbDYDmFFXV9eHsCR1hjdlusFKKixcDdc+nflZ03dkJKIxu4BNwMjQ8xHAm1GcWLPHSo+GNyWiW6JUXjB1B2j6jkCSE8VLwFgzG2Nmg4CLgMVRnFjrUYhIKTR9R0YiEoWZ3Q8sB8aZ2SYz+7K77wG+CjwJ/AF40N1fieJ6KlGIiJQuEW0U7n5xge1PAE9EfT21UUhBHa2ZRuzwHE5qwK4aX/nKV+IOoSYlIlFUmrsvAZY0NjZeGXcskhDB2Imgp9PoWZn+6xozUVUuvPDCfr0+WIMC4FMjix+bJqlMFCL72bAQ3m7fV5JQgqhKGzdmetSPHNm3b/mgIXve6WqfCEtlolDVk+Q1tEHdYKvcpZdeCvRvHIXWoNhfIhqzK02N2dLl6TMyD0m1havhqH/bV+0k3aWyRCHSJbwYkcZNpJbGThSXyhKFxlEIsP+0zerdlHqfGqlqp3xSmShU9ZRSHa2wtHlfggjPyKl1r1PvlCPhvnPjjiKZVPUk6RHu/qqkUJOuvvrquEOoSUoUkj65q4WFB9dJVZsxY0bcIdSkVCYKdY8VXpyTWV9CM8PWlHXrMq3S48apRTpKaqOQdAn3bBraoJJEjZkzZw5z5syJO4zyiHHd7FSWKCTljpia+VdrXku1GD0rkyg2LIylfS2VJQpJoWCyP8gkCCUJyQrP75RYdS2xjvNRiUJqW77J/kRCgsF2GmhXWCpLFBpwV4MKTcURThIaKyE5gtKE5ncqLpUlCk0zXoPCU3FIan3zm9/s1fEqTZQmlYlCalhQqsjXBhGULKRmTZ06tdevUWmiZ0oUUluCksXSZiWGFGpvbwegoaEh1jhqjRKF1Ka32/f9POlO+PNPYwtFKmfu3LlA/9ajSLRgLEWF29pS2ZgtKTC0Yd/PdS3qEivVL+ixF57MskKUKEREqkGMYylSmSjUPVZEpHSpTBSa66lGaBlTkYpQY7ZUL42dkBw333xz3CHUJCUKqX4xzagpyXPqqaeWfGx4VLYUl8qqJ6kxMfQCkWRatmwZy5YtK+lYjcounUoUknzBxH6jZ+XvP66BdZJ13XXXAaWPo9Co7NIoUUjy9XWt6yN6P52DiOxPiUKq1/Cm4qUJDbATiYQShVSXcDWUiFREzTRmm1mzmf3WzO4ws+a445EyCaqhwg3YuaNVhzcpkYhEKBElCjO7G5gObHH340LbpwH/BgwAfuju/1TkNA7sAg4GNpUxXEmS4U0wtW3fbLHBc0ml+fPnxx1CTUpKiWIBMC28wcwGALcBZwHjgYvNbLyZTTCzx3Mew4HfuvtZwD8CN1Q4fimmozXzRZ473qHQ9p72iRTQ0NBQ0hTjVbFOdoIkokTh7r8xs9E5mycBHe7+GoCZPQDMdPd5ZEofhbwNHFRop5m1AC0Ao0aN6k/YUqpCvZaK9WYK7xMp0dKlmdH6PS1gpDEUvZOIRFHAkcDG0PNNwORCB5vZucBngQ8B3y90nLu3Aq0AjY2NHkWgIpIMN910E1DaSncaQ1G6JCcKy7Ot4Be7uz8CPFLSic1mADPq6ur6GJqISHokpY0in03AyNDzEcCbUZxYs8fGZMszxdsdemqXeLu9+8p1IlIRSS5RvASMNbMxwGbgIiCSPo8qUcSo2AjrntolglXr8h0TdIdVt1iRyCWiRGFm9wPLgXFmtsnMvuzue4CvAk8CfwAedPdXorieShRVKOj2OrUt/ypfdS2ZfRVeS1gkDRJRonD3iwtsfwJ4IurrqUQhUpvuvPPOuEOoSYkoUVSaShQitWncuHGMG6c+r1FLRIlCUqijtXB7RCnjJ4JjYlpsXpJpyZIlAMyYMaPgMVqwqPdSWaIwsxlm1rpjx464Q0mv/iw2pAZrKeDWW2/l1ltvLXqMBtv1XioThaqeEmJ4U99KBHUtKklIv2iwXe+kMlGIiEjpUpkoVPUkIlK6VCYKVT2JiJROvZ5EpGbce++9cYdQk5QopDyCJUtL7er6/xph+8r9973dDu+pilBKM3LkyG7PL8lOE/q5sZneTjXR02nLM5n/X9lZCJ7fnOnyW87G+VRWPamNogLyJYl8PZWCrq65SeKI7DTRwfxO+V6nJU8lx6JFi1i0aFHX82c3Zh6Prct8oQZdY6tW8HnPdi8PEl+531cqE4XaKGIQzNWUmyzCXV3D+07/Fczywsuaam4nyeP222/n9ttvjzuM8snpGj5rQmUGDqYyUYiISOmUKEREpCglChERKSqViUKN2SIipUtl91h3XwIsaWxsvDLuWEQkOg899FDcIdSkVCYKEalNw4YNizuEmpTKqicRqU0LFixgwYIFcYdRc5QoRKRmKFGUhxKFiIgUpURRSzpaYWlz5t/+HN/b8+Q7b7E5nt5uzzxEIrblf2DtWzD9frgw1K79/ObMv2vfyjyqXjDfU4WksjHbzGYAM+rq6uIOJVrh+ZVKmdqi0PG9PU++8xYSnpsp3zxNwRxPIn2w9R34792wekv+/eMPz/xb1ZMDjp6V+f+5YWHFprBJZaJQ99gY1bUU/3Cf/qvKxSKps+j8uCOIQF1L/9ac74NUJgoRqU2nXfcEL74ZdxS1R4lCRGrGgQcN4YBBcUdRe9SYLSI1o+PJH7Dz2R/EHUbNUaIQkZqxcdmDvPPyg3GHUXOUKEREpCglChERKUqJQkREilKiEBGRoszd444hNmb2FvBGH18+DNgaYTj9lbR4QDGVSjGVJmkxJS0e6F9MR7n74fl2pDpR9IeZrXD3xrjjCCQtHlBMpVJMpUlaTEmLB8oXk6qeRESkKCUKEREpSomi7yo3x29pkhYPKKZSKabSJC2mpMUDZYpJbRQiIlKUShQiIlKUEkUPzOz/mtmrZvZ7M3vUzD4U2netmXWY2Toz+2xo+0QzW53d9z0zs4hj+lsze8XM3jezxtD2S8ysPfR438wasvvasnEG+4ZXKKbRZtYZuu4doX1lu09F4jnDzFZmr7vSzE4P7YvlHmX3xfJZyolhUei9bzCz9uz2gr/DcjOzb5vZ5tC1zw7ty3vPKhBT3u+EOO9T9vrTsveiw8y+HunJ3V2PIg/gTODA7M+3ALdkfx4PrAIOAsYAfwIGZPe9CHwCMOAXwFkRx3QsMA5oAxoLHDMBeC30vOCx5YwJGA2sKfCast2nIvGcCHw0+/NxwOYE3KPYPktFYr0V+D89/Q4rEMe3gWvybC94zyoQU6HvhDjv04DsPTgaGJS9N+OjOr9KFD1w91+6+57s0+eBEdmfZwIPuPtud38d6AAmmdlHgL9x9+We+Q3eA5wTcUx/cPd1PRx2MXB/lNctpsSYupT7PhWKx91fdvdgaZtXgIPN7KCortuXmIjxs5RPttRyARX8/PRB3ntWiQsX+U6I0ySgw91fc/d3gQfI3KNIKFH0zpfI/FUHcCSwMbRvU3bbkdmfc7dX2oXs/x/937NF4uvLWYWRxxgze9nMnjGz07LbknCfzgNedvfdoW1x3KOkfZZOA/7L3deHtuX7HVbKV7PVPHeb2dDstkL3rNLC3wkQ330q6/3QCneAmS0Fjsiz6xvu/lj2mG8Ae4D7gpflOd6LbI88piKvnQy84+5rQpsvcffNZvZB4GHgUjJ/oZY7pv8ARrn7NjObCPzMzOqJ4D718x7Vk6k2ODO0Oa57VNbPUrcLlRZfbmk07+/Q3f+7P7GUEhNwO3Ajmfd9I5kqsS9RhntTakxFvhPKep96CjnPtsjuhxIF4O5Ti+03s8uB6cCUbBUAZDL2yNBhI4A3s9tH5NkeaUw9uIic0oS7b87+u9PMFpIpqvbqS7AvMWX/Wt+d/Xmlmf0JOIYI7lNf75GZjQAeBS5z9z+FzhfLPaLMn6WwEj7rBwLnAhNDryn0O1zRn1hKjSkU213A49mnhe5ZJPrynVDu+9SDst4PVT31wMymAf8IfN7d3wntWgxcZGYHmdkYYCzworv/B7DTzE7JVl1cBhT96zbieA8A/pZMHWWw7UAzG5b9eSCZD/ia/GeIPJ7DzWxA9uejydyn1+K6T9keKj8HrnX350LbY7tHJOuzNBV41d27qrwK/Q7LHEdw7Y+Enn6Bfb+TvPesQjHl/U6I8z4BLwFjzWyMmQ0i88fi4sjOHkcLfTU9yDSSbQTas487Qvu+QaanwTpCvVGARjIf6D8B3yc7sDHCmL5A5i+I3cB/AU+G9jUDz+ccfwiwEvg9mQbcfyPiHiKFYiLTDvAKmV4YvwNmVOI+FYnnm8D/hH6f7cDwOO9RnJ+lPDEuAP4uZ1vB32G5H8C9wOrs72Ux8JGe7lkFYsr7nRDnfcpe/2zgj9l78o0oz62R2SIiUpSqnkREpCglChERKUqJQkREilKiEBGRopQoRESkKCUKEREpSolCRESKUqIQEZGilChEYmJmYy2zWNIKM/tnM+uIOyaRfJQoRGKQnRPoHuB/u3sjMJjM9A8iiaPZY0XicQ6w1t1/l33+B+AvsUUjUoRKFCLxOJHMhHKBE8hMJieSOEoUIvHYBnwcuhaauozMDKkiiaPZY0VikF374ufAEOAJ4BIyq6O9H2tgInkoUYjEzMxGAg+5++S4YxHJR1VPIvE7AVU7SYKpRCEiIkWpRCEiIkUpUYiISFFKFCIiUpQShYiIFKVEISIiRSlRiIhIUUoUIiJSlBKFiIgU9f8BvhAfrQqao/QAAAAASUVORK5CYII=\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": "dc94e005",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median[p_b|s+b] = 1.000e-06\n",
      "median[Z_b|s+b] = 4.753\n"
     ]
    }
   ],
   "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.\n",
    "p_val_b = len(qb[qb<=med_q_sb])/numExp\n",
    "Z_b = stats.norm.ppf(1.-p_val_b)\n",
    "print(\"median[p_b|s+b] = {:.3e}\".format(p_val_b))\n",
    "print(\"median[Z_b|s+b] = {:.3f}\".format(Z_b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab10f62a",
   "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
}