{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# **Statistics lecture 5 Hands-on session : solutions notebook**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the companion notebook to lecture 4 in the statistical course series, covering the following topics:\n", "1. Parameter estimation using simple tools\n", "2. Parameter estimation using the profile likelihood\n", "3. Limit-setting, using simple tools\n", "4. Limit-setting using the profile likelihood\n", "\n", "First perform the usual imports:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.stats\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Profiling basics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous lectures, we have discussed 3 types of statistical results:\n", "* Testing for discovery\n", "* Computing confidence intervals for parameters (central value +/- uncertainty)\n", "* Setting upper limits\n", "\n", "For simplicity, all this was done using models with only 1 parameter (the signal $s$), which was the *parameter of interest*. Models usually also have additional *nuisance parameters* which aren't interesting in themselves, but need to be included to describe the data.\n", "\n", "These parameters need to be included in the treatment, but fortunately this is quite easy in the likelihood ratio framework. The basic quantity we used before was\n", "$$\n", "t(s) = -2 \\log \\frac{L(s)}{L(\\hat{s})}\n", "$$\n", "\n", "This allows to compare 2 hypotheses:\n", "* The free hypothesis $\\hat{s}$ which corresponds to the best-fit $s$ in the data.\n", "* The fixed hypothesis $s$ (for instance $s=0$ for discoveries)\n", "\n", "If we have an additional nuisance parameter $\\theta$, the likelihood is now a function of both parameters: $L(s, \\theta)$. How do we set $\\theta$ ? The interesting answer is that in both cases, we use the best-fit value of $\\theta$ :\n", "* For the free hypothesis $\\hat{s}$, we use $\\hat{\\theta}$, the best-fit value of $\\theta$.\n", "* For the fixed hypothesis $s$, we use $\\hat{\\hat{\\theta}}(s)$, the best-fit value of $\\theta$ at a fixed $s$.\n", "\n", "So in the end we use\n", "$$\n", "t(s) = -2 \\log \\frac{L(s, \\hat{\\hat{\\theta}}(s))}{L(\\hat{s},\\hat{\\theta})}\n", "$$\n", "with the appropriate best-fit values. This is now called the *profile likelihood*.\n", "\n", "A very nice feature of this is that $t(s)$ is still a function only of $s$, so we can perform the rest of the computation in the same way as before. The only difference is that we need to compute best-fit values for $\\theta$ in order to compute the appropriate ilikelihood values.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Profiling on a single-bin example\n", "\n", "We reuse the same examples as in the previous sections: the full binned example defined below:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdpElEQVR4nO3deZwV5Z3v8c+XRZGoURYZtMF2wS0KGjvGiYkyISbxFQeXa7gaSBiXIZNxErPoiOOoYybM4Otm0TuZOEFjbBPcrkswiWMG0dbojRowiInoYBSwlcjiDlcR+N0/qrptml7qdHedok9936/XeZ16nnOqnt9h+XX1c6p+jyICMzMrjwFFB2BmZtXlxG9mVjJO/GZmJePEb2ZWMk78ZmYlM6joALIYMWJE1NfXFx2GmVm/smjRorURMbJ9f79I/PX19SxcuLDoMMzM+hVJKzrq91SPmVnJOPGbmZWME7+ZWcn0izl+M7OOvPvuuzQ3N/P2228XHUqhhgwZQl1dHYMHD870fid+M+u3mpub2WWXXaivr0dS0eEUIiJYt24dzc3N7LPPPpn28VSPmfVbb7/9NsOHDy9t0geQxPDhwyv6rceJ38z6tTIn/RaV/hk48ZuZlYzn+M2sdtzYx2f/n8tnvZJ58+ZxySWXMGDAAAYNGsSVV17JRz/6UQDuuecezjvvPDZv3sw555zDzJkz+3x8J/52Jk6cCEBTU1OhcZhZ7Zo0aRKTJ09GEkuWLGHKlCk8/fTTbN68mXPPPZf58+dTV1fHhz70ISZPnswhhxzSp+N7qsfMrBfWr1/PZz7zGSZMmMChhx7KLbfc0u0+O++8c+u8/Pr161u3H3vsMfbff3/23XdfdthhB04//XTmzZvX5zH7jN/MrBfuuece9txzT375y18C8Prrr/O1r32N+++/f5v3nn766a1TN3feeScXXXQRq1evbt33xRdfZMyYMa3vr6ur49FHH+3zmJ34zcx64bDDDuP888/nwgsv5MQTT+RjH/sY3/ve97rd75RTTuGUU07hwQcf5JJLLuHee++lozXQ87hqKbfEL+lAoO3vPPsClwI3pP31wHJgSkS8mlccZmZ5OuCAA1i0aBF33303F110EZ/85Cd59dVXuz3jb3Hsscfyxz/+kbVr11JXV8cLL7zQ+lpzczN77rlnn8ecW+KPiGeAwwEkDQReBO4EZgILImK2pJlp+8K84jAzy9NLL73EsGHDmDZtGjvvvDPXX389P/vZz7rc59lnn2W//fZDEo8//jgbN25k+PDh7Lbbbixbtoznn3+evfbai5tvvpkbb7yxz2Ou1lTPJOCPEbFC0knAxLS/EWjCid/M+kJOl1925cknn+SCCy5gwIABDB48mKuvvrrbfW6//XZuuOEGBg8ezE477cQtt9yCJAYNGsT3v/99PvWpT7F582bOOussPvCBD/R5zOpoTqnPB5GuAx6PiO9Lei0idmvz2qsRsXsH+8wAZgCMHTv2yBUrOlxPoM/5ck6z/mPp0qUcfPDBRYexXejoz0LSoohoaP/e3C/nlLQDMBn4P5XsFxFzIqIhIhpGjtxm5TAzM+uhalzHfwLJ2f7LaftlSaMB0ufVVYjBzMxS1Uj8ZwA3tWnfBUxPt6cDfX93gpmZdSrXxC9pKHA8cEeb7tnA8ZKWpa/NzjMGMzPbWq5X9UTEBmB4u751JFf5mJlZAVyrx8ysZJz4zaxmSH376M7y5cs59NBD+/xzNDU1ceKJJ/b5cVs48ZuZlYwTv5lZL2zatInp06czfvx4TjvtNDZs2LDNexYvXszRRx/N+PHjOeWUU3j11aQ82cSJE7nwwgs56qijOOCAA/j1r3+91X5btmxh3LhxrFmzprW9//77s3bt2l7F7MRvZtYLzzzzDDNmzGDJkiXsuuuu/OAHP9jmPV/4whe44oorWLJkCYcddhiXX35562ubNm3iscce48orr9yqH2DAgAFMmzaNuXPnAnDvvfcyYcIERowY0auYnfjNzHphzJgxHHPMMQBMmzaNhx56aKvXX3/9dV577TWOO+44AKZPn86DDz7Y+vqpp54KwJFHHsny5cu3Of5ZZ53FDTfcAMB1113HmWee2euYnfi3IxMnTmytFWRm/UP7evmV1s/fcccdARg4cCCbNm3a5vUxY8YwatQo7rvvPh599FFOOOGEngebcuI3M+uFlStX8pvf/AaAm266qXXR9Bbvf//72X333Vvn73/yk5+0nv1ndc455zBt2jSmTJnCwIEDex2zE7+Z1YyIvn1kcfDBB9PY2Mj48eN55ZVX+NKXvrTNexobG7ngggsYP348ixcv5tJLL63oc02ePJm33nqrT6Z5wEsvmpn1WH19PU899VS37zv88MN55JFHtulvW/59xIgRrXP87ad9n3jiCSZMmMBBBx3U25CBEiT+ni5XWcl+VVjSwMxKavbs2Vx99dWtV/b0BU/1mJltx2bOnMmKFSu2+e6gN5z4zaxfq8Yqgtu7Sv8MnPjNrN8aMmQI69atK3XyjwjWrVvHkCFDMu9T83P8lZkLPAK8A9QDs4CpRQZkZl2oq6ujubm5taRBWQ0ZMoS6urrM73fibzWXZG33d9L2irQNTv5m26fBgwezzz77FB1Gv+OpnlYXA+2LK21I+83MaocTf6uVFfabmfVPTvytxlbYb2bWP+W92Ppukm6T9LSkpZL+XNIwSfMlLUufd88zhuxmAUPb9Q1N+83MakfeZ/xXAfdExEHABGApMBNYEBHjgAVpezswFZgD7Ji2907b/mLXzGpLblf1SNoVOBb4K4CI2AhslHQSMDF9WyPQBFyYVxyVmQpck243FRiHmVl+8jzj3xdYA/xY0u8kXSvpfcCoiFgFkD7vkWMMZmbWTp6JfxDwQeDqiDgCWE8F0zqSZkhaKGlh2W/OMDPrS3km/magOSIeTdu3kfwgeFnSaID0eXVHO0fEnIhoiIiGkSNH5himgVf/MiuT3BJ/RPwJeEHSgWnXJOAp4C5geto3HZiXVwxmZratvEs2fBmYK2kH4DngTJIfNrdKOpvk7qjP5hyDmZm1kWvij4jFQEMHL03Kc1wzM+uci7TlyKt/mdn2yCUbzMxKxonfzKxknPjNzErGid/MrGT85e42mooOwMwsVz7jNzMrGSd+M7OSceI3MysZJ34zs5Lp9MtdST8HOr0vNCIm5xKRmZnlqqurer6dPp8K/Bnw07R9BrA8x5jMzCxHnSb+iHgAQNI/R8SxbV76uaQHc4/MzMxykWWOf6SkfVsakvYBvDKK9RkvAmNWXVlu4Poa0CTpubRdD3wxt4jMeuPGHpZErcTnXBLV+rduE39E3CNpHHBQ2vV0RLyTb1hmZpaXrq7qObWTl/aTRETckVNMZmaWo67O+P+yi9cCcOLfjnkRGDPrTFdX9ZxZzUDMzKw6ur2qR9IoST+S9J9p+5B0oXQzM+uHslzOeT3wK2DPtP3fwFezHFzScklPSlosaWHaN0zSfEnL0ufdexC3mZn1UJbEPyIibgW2AETEJmBzBWP8RUQcHhENaXsmsCAixgEL0rYByVoATQXHYGa1LkviXy9pOGndHklHA6/3YsyTgMZ0uxE4uRfHMjOzCmVJ/F8H7iK5jPNh4AbgyxmPH8B/SVokaUbaNyoiVgGkz3t0tKOkGZIWSlq4Zs2ajMP1Z3NJ7o0bkD7PLTIYM6thWW7gelzSccCBgIBnIuLdjMc/JiJekrQHMF/S01kDi4g5wByAhoaGGr9wcC4wA9iQtlekbYCphUTUK5XePbu6B/v57lmzHuv0jF/Sx9PnU4HJJIn/AOAvu7i5aysR8VL6vBq4EzgKeFnS6PTYo3nvv32JXcx7Sb/FhrTfzKxvdXXGfxxwHx3fyNXtDVyS3gcMiIg30+1PAt8kmTaaDsxOn+f1IO4as7LCfjOznuvqBq7L0s1zIqKSq3hajALuVHIr6CDgxrTuz2+BW9N7AVYCn+3BsWvMWJLpnY76zcz6VpbqnM9Kug34cUQ8lfXAEfEcMKGD/nXApOwhlsEstp7jBxia9lu/4cqg1k9kuapnPMlNW9dKeiS92mbXnOMqmakk32PvTfL9+d5pux9+sWtm271uE39EvBkR10TER4C/By4DVklqlLR/7hGWxlSSFS23pM9O+tUw8VvJw6xMstTqGShpsqQ7gauA7wD7Aj8H7s45PjMz62NZ5viXAfcD/ysi/m+b/tskHdvJPmZmtp3KkvjHR8RbHb0QEV/p43isBmhqpV9ATkz3a8q8R3yuwiHMrFWWOf4Ok77VkrnAI8ADVL9cRJFjm5VTlqt6rKa1lItoWUa5pVxENRJwkWOblZcTf+kVWS7CpSrMitDVYutf72rHiPhu34dj1VdkuYhiS1XMfRgeWQbvbIL682DWFJh6TFWGNitUV2f8u6SPBuBLwF7p42+AQ/IPzaqjs7IQ1SgXUdzYcx+GGdcmSR9gxdqkPffh3Ic2K1yniT8iLo+Iy4ERwAcj4hsR8Q3gSKCuWgFa3maRlIdoq1rlIoob++JbYcPGrfs2bEz6zWpdljn+sUDb/yIbSS6/sJpQZLmIlrF3TNvVG3vl2sr6zWpJluv4fwI8lt65G8ApJKtwWc2YSnElIqYC16TbTVUbdeyIZHqno36zWpflOv5ZwJnAq8BrwJkR8S85x2WWq1lTYOgOW/cN3SHpN6t1WS/nHAq8ERFXAc2S9skxJrPcTT0G5pwDO6a/8+49Imn7qh4rg26neiRdRnJlz4HAj4HBwE8B/xexfm3qMXDN/cl20z8WG4tZNWU54z+FZM3d9dC6ju4ueQZlVgYuCW1FyZL4N0ZEkHyx27KWrpmZ9VNZEv+tkn4I7Cbpr4F7gWvzDcvMzPLS7Rx/RHxb0vHAGyTz/JdGxPysA0gaCCwEXoyIEyUNA24huRdgOTAlIl7tQexm23BJaLPuZVmB64qImB8RF0TE+RExX9IVFYxxHrC0TXsmsCAixgEL0raZmVVJlqme4zvoOyHLwSXVAZ9h66mhk4DGdLsRODnLsczMrG90mvglfUnSk8BBkpa0eTwPPJnx+FeSLNC+pU3fqIhYBZA+79Gz0M36r5bKoA8sTSqDujicVVNXc/w3Av8J/CtbT8e8GRGvdHdgSScCqyNikaSJlQYmaQbJqhyMHVuNSpFWnKaiA6iqziqDgm8gs+roqjrn6xGxHLgKeCUiVkTECuBdSR/OcOxjgMmSlgM3Ax+X9FPgZUmjAdLn1Z2MPyciGiKiYeTIkRV9KLPtmSuDWtGyzPFfDbRdd3d92teliLgoIuoioh44HbgvIqYBdwHT07dNB+ZVFLFZP+fKoFa0LIlf6Q1cAETEFrJV9ezMbOB4SctIvjie3YtjmfU7nVUAdWVQq5Ysif85SV+RNDh9nAc8V8kgEdEUESem2+siYlJEjEufu/2+wKyWuDKoFS1L4v8b4CPAi0Az8GHSL13NrHKuDGpFy3Ln7mqSOXoz6yOuDGpFynLn7gGSFkj6fdoeL8n/VM3M+qksUz3XABcB7wJExBL8G4BZvzZx4kQmTpxYdBhWkCyJf2hEPNaub1MewZiZWf6yXJa5VtJ+vFeP/zRgVa5RmfVDrgxq/UWWxH8uMIekZs+LwPPA1FyjMjOz3GS5quc54BPpylsDIuLN/MMyM7O8ZLmqZ7ik/w38GmiSdJWk4fmHZmZmecjy5e7NwBrgfwCnpdu35BmUWfU0UbbqoGZZ5viHRcQ/t2l/S9LJOcVjZj11o7K/t6UmbiX7fK7SL69te5XljP9+SadLGpA+pgC/zDswMzPLR5bE/0WSRVneSR83A1+X9KakN/IMzszM+l6Wq3p2qUYgZmZWHVmu6jm7XXugpMvyC8nMzPKUZapnkqS7JY2WdBjwCODfAszM+qksUz2fk/Q/gSeBDcAZEfFw7pGZWU1qKQ7X1NRUaBxllmWqZxxwHnA7sBz4vKShOcdlZmY5yTLV83Pg0oj4InAcsAz4ba5RmZlZbrLcwHVURLwBkC66/h1Jd+UblplVqrLqoBPTfZoy7+HKoLUjyxn/TpJ+JOkeAEmHAMd2t5OkIZIek/SEpD9IujztHyZpvqRl6fPuvfsIZmZWiSyJ/3rgV8DotP3fwFcz7PcO8PGImAAcDnxa0tHATGBBRIwDFqRtMzOrkiyJf0RE3ApsAYiITcDm7naKxFtpc3D6COAkoDHtbwROrjBmsxrRhAvEWRGyJP71aRnmlhW4jgZez3Lw9GavxSQloeZHxKPAqIhYBZA+79HJvjMkLZS0cM2aNVmGMzOzDLIk/q8DdwH7SXoYuAH4cpaDR8TmiDgcqAOOknRo1sAiYk5ENEREw8iRI7PuZmZm3chyA9fjko4DDgQEPBMR71YySES8JqkJ+DTwsqTREbFK0mjeKxBrZmZVkOWMn4jYFBF/iIjfZ036kkZK2i3d3gn4BPA0yW8P09O3TQfmVRy1mZn1WJbr+HtqNNAoaSDJD5hbI+IXkn4D3JoWf1sJfDbHGMzMrJ3cEn9ELAGO6KB/HTApr3HNzKxrnSZ+SR/saseIeLzvwzEzy48LxCW6OuP/ThevBfDxPo7FzMyqoNPEHxF/Uc1AzKz/UgVrtvdkn/A6730q0xx/ev39IcCQlr6IuCGvoMzMLD/dJv50mcWJJIn/buAE4CGSG7nMzKyfyXId/2kkV+H8KSLOBCYAO+YalZnlaC7JCqoPAPVp28okS+L/fxGxBdgkaVeSO233zTcsM8vHXGAGSfFcgBVp28m/TLIk/oXpHbjXAIuAx4HH8gzKzPJyMcnS2W1tSPutLLLU6vnbdPM/0sVYdk1vzjKzfmdlhf19rWWa6R2SaaZZwNQqjW0tsiy2vqBlOyKWR8SStn1m1p+MrbC/L3maaXvRaeJPl04cBoyQtHu6ZOIwSfXAnlWL0Mz60CxgaLu+oWl/3jzNtL3oaqrniyRLLO5JMq/f4g3g33OMycxy0zKtcjHJ9M5YqjfdUvQ0U7G2p3IRXd25exVwlaQvR8S/VTEmM8vVVIqZVx9LMr3TUb9VU5aren4o6SuSbksffydpcO6RmVmNKXKaydrKUrLhByQLpf8gbX8euBo4J6+gzKwWtfyWcTbJF7x70xfTTK4TVLmuyjIPiohNwIciYkKbl+6T9ET+oZlZ7ZlKcksQQFOBcZRbV1M9LTdpbZa0X0unpH2BzblGZWZmuelqqqfll6HzgfslPZe264Ez8wzKzGpZU9EBlF5XiX+kpK+n2z8EBgLrSUozHwHcn3NsZmaWg64S/0BgZ9478ydtA+ySW0RmZparrhL/qoj4Zk8PLGkMSc3+PwO2AHMi4qr0buBbSKaMlgNTIuLVno5jZmaV6erL3R5cJLWVTcA3IuJg4GjgXEmHADOBBRExDliQts3MrEq6SvyTenPgiFgVEY+n228CS4G9gJOAxvRtjcDJvRnHzMwq02nij4hX+mqQtLDbEcCjwKiIWJWOsQrYo5N9ZkhaKGnhmjVr+ioUM7PSy7TYem9I2hm4HfhqRLyhjLfMRcQcYA5AQ0NDDd47Z2b9Xd53DUM+dw5nqdXTY2lNn9uBuRFxR9r9sqTR6eujSZZyNDOzKskt8Ss5tf8RsDQivtvmpbuA6en2dGBeXjGYmdm28pzqOYakoNuTkhanff8AzAZulXQ2SSHuz+YYg5lZG01FB7BdyC3xR8RDdH5JaK+uGDIzs57LdY7fzMy2P078ZmYl48RvZlYyTvxmZiXjxG9mVjJO/GZWEnNJigIPSJ/nVnnsR4AHChh7W7mXbDAzK95cYAawIW2vSNvQ28Xes4/9TgFjd8xn/GZWAhfzXtJvsSHtr+WxO+bEb2YlsLLC/loZu2NO/GZWAmMr7K+VsTvmxG9mJTALGNqub2jaX8tjd8yJ38xKYCrJ8h57k5QQ2zttV+PL1SLH7piv6jGzkphKccm2yLG35TN+M7OSceI3MysZJ34zs5Jx4jczKxknfjOzknHiNzMrmdwSv6TrJK2W9Ps2fcMkzZe0LH3ePa/xzcysY3me8V8PfLpd30xgQUSMAxakbTMzq6LcEn9EPAi80q77JKAx3W4ETs5rfDMz61i15/hHRcQqgPR5j87eKGmGpIWSFq5Zs6ZqAZqZ1brt9svdiJgTEQ0R0TBy5MiiwzEzqxnVTvwvSxoNkD6vrvL4ZmalV+3EfxcwPd2eDsyr8vhmZqWX5+WcNwG/AQ6U1CzpbGA2cLykZcDxadvMzKoot7LMEXFGJy9NymtMMzPr3nb75a6ZmeXDid/MrGSc+M3MSsaJ38ysZJz4zcxKxonfzKxknPjNzErGid/MrGSc+M3MSsaJ38ysZJz4zcxKxonfzKxknPjNzErGid/MrGSc+M3MSsaJ38ysZJz4zcxKxonfzKxknPjNzErGid/MrGQKSfySPi3pGUnPSppZRAxmZmVV9cQvaSDw78AJwCHAGZIOqXYcZmZlVcQZ/1HAsxHxXERsBG4GTiogDjOzUhpUwJh7AS+0aTcDH27/JkkzgBlp8y1Jz1QwxghgbY8jrJBUrZG6Hdufu/pjV5U/d6uqffYiP3cH41f6uffuqLOIxN/RH2Ns0xExB5jTowGkhRHR0JN9+zN/7nIp6+eG8n72vvrcRUz1NANj2rTrgJcKiMPMrJSKSPy/BcZJ2kfSDsDpwF0FxGFmVkpVn+qJiE2S/g74FTAQuC4i/tDHw/RoiqgG+HOXS1k/N5T3s/fJ51bENtPrZmZWw3znrplZyTjxm5mVTE0l/rKWgpA0RtL9kpZK+oOk84qOqZokDZT0O0m/KDqWapG0m6TbJD2d/r3/edExVYOkr6X/xn8v6SZJQ4qOKQ+SrpO0WtLv2/QNkzRf0rL0efeeHr9mEn/JS0FsAr4REQcDRwPnluizA5wHLC06iCq7CrgnIg4CJlCCzy9pL+ArQENEHEpyccjpxUaVm+uBT7frmwksiIhxwIK03SM1k/gpcSmIiFgVEY+n22+SJIG9io2qOiTVAZ8Bri06lmqRtCtwLPAjgIjYGBGvFRpU9QwCdpI0CBhKjd4DFBEPAq+06z4JaEy3G4GTe3r8Wkr8HZWCKEXya0tSPXAE8GjBoVTLlcDfA1sKjqOa9gXWAD9Op7iulfS+ooPKW0S8CHwbWAmsAl6PiP8qNqqqGhURqyA52QP26OmBainxZyoFUcsk7QzcDnw1It4oOp68SToRWB0Ri4qOpcoGAR8Ero6II4D19OLX/v4indM+CdgH2BN4n6RpxUbVP9VS4i91KQhJg0mS/tyIuKPoeKrkGGCypOUkU3sfl/TTYkOqimagOSJafqu7jeQHQa37BPB8RKyJiHeBO4CPFBxTNb0saTRA+ry6pweqpcRf2lIQkkQy37s0Ir5bdDzVEhEXRURdRNST/H3fFxE1fwYYEX8CXpB0YNo1CXiqwJCqZSVwtKSh6b/5SZTgS+027gKmp9vTgXk9PVAR1TlzUaVSENurY4DPA09KWpz2/UNE3F1cSJazLwNz05Oc54AzC44ndxHxqKTbgMdJrmT7HTVaukHSTcBEYISkZuAyYDZwq6SzSX4IfrbHx3fJBjOzcqmlqR4zM8vAid/MrGSc+M3MSsaJ38ysZJz4zcxKxonfaoqk+rYVDdu9du32ULyuqxjNqqFmruM3605EnFN0DH1B0qCI2FR0HNZ/+YzfatEgSY2SlqQ164cCSGqS1JBuvyVplqQnJD0iaVT7g0j6p7QuepOk5yR9Je3f6oxd0vmS/qnNGN+T9GBaJ/9Dku5Ia6h/K0OMR0p6QNIiSb9qc4t+k6R/kfQASRlqsx5z4rdadCAwJyLGA28Af9vBe94HPBIRE4AHgb/u5FgHAZ8iKft9WVoTqTsbI+JY4D9Ibqs/FzgU+CtJwzuLMT32vwGnRcSRwHXArDbH3S0ijouI72SIwaxTTvxWi16IiIfT7Z8CH+3gPRuBlhW7FgH1nRzrlxHxTkSsJSmKtc1vBh1oqRH1JPCHdL2Ed0hKK7QUEuwoxgNJfkDMT0tv/CNJscEWt2QY26xbnuO3WtS+DklHdUnejffqlWym8/8L77TZbnnfJrY+aWq//F/LPlva7b+lzTgdxSiSHxSdLaO4vpN+s4r4jN9q0dg2a9CeATzUx8d/GdhD0nBJOwIn9uAYHcX4DDCypV/SYEkf6JOIzdpw4rdatBSYLmkJMAy4ui8PntaC/ybJKme/AJ7uwWG2iTFdMvQ04ApJTwCLKVe9easSV+c0MysZn/GbmZWME7+ZWck48ZuZlYwTv5lZyTjxm5mVjBO/mVnJOPGbmZXM/wf1Cbe/Fm+KnwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "nbins = 10\n", "x = np.linspace(0.5, nbins - 0.5, nbins)\n", "# The background follows a linear shape\n", "b_yields = np.array([ (1 - i/2/nbins) for i in range(0, nbins) ])\n", "b_yields *= b_yields/np.sum(b_yields)\n", "# The signal shape is a peak\n", "s_yields = np.zeros(nbins)\n", "s_yields[3:7] = [ 0.1, 0.4, 0.4, 0.1 ]\n", "# Now generate some data\n", "s = 10\n", "b = 500\n", "s_and_b = s*s_yields + b*b_yields\n", "b_only = b*b_yields\n", "np.random.seed(1) # make sure we always generate the same\n", "data = [ np.random.poisson(s*s_yield + b*b_yield) for s_yield, b_yield in zip(s_yields, b_yields) ]\n", "s30_and_b = 30*s_yields + b*b_yields\n", "plt.bar(x, s30_and_b, color='orange', label='s=30')\n", "plt.bar(x, b_only, color='b', yerr=np.sqrt(b_only), label='b only')\n", "plt.scatter(x, data, zorder=10, color='k')\n", "plt.xlabel('bin number')\n", "plt.ylabel('Total expected yield')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and also the single-bin measurement in bin 4:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa2ElEQVR4nO3deZgU9Z3H8feHQwkKHlwBAQcVRDRilGCiRDGuUVeDkqDxYCUoi0ncFeOaSHYTNWvc4GMOXK+Ix4rRqGyiC4pHECTGW1BBDbq4iIgSOTxQvGD47h9dgwPMUXNUd8/U5/U8/XRVdR0feoZv1/z6V79SRGBmZvnRptQBzMysuFz4zcxyxoXfzCxnXPjNzHLGhd/MLGfalTpAGl27do2KiopSxzAza1Hmz5+/OiK6bbm8RRT+iooK5s2bV+oYZmYtiqTXalruph4zs5xx4TczyxkXfjOznGkRbfxmZjVZv349y5cv5+OPPy51lJLq0KEDvXv3pn379qnWd+E3sxZr+fLldOrUiYqKCiSVOk5JRARr1qxh+fLl9OvXL9U2buoxsxbr448/pkuXLrkt+gCS6NKlS4P+6nHhN7MWLc9Fv0pD3wMXfjOznHEbv5m1Hr9v5rP/U7K5X8n06dP56U9/Sps2bWjXrh2TJ09m2LBhANx///1MmDCByspKxo0bx8SJE5v9+K2/8Df3L4KZlY8B98Gaddntf002IwYcvt9OjJhzA5JY+OJiThw7mpee+AOVlZWc9b1xzPrDlfTu1YMvHTGGESNGMGjQoGY9vpt6zMyaYN26jzjmpHMYfOgp7DPs29xx15/q3Wb77Ttuapdf9+FHm6afeuZF9ujXh90qerPNNu05aeQRTJ8+vdkzt/4zfjOzDN0/53F6fb4rM2+fDMB7az/gB//2ax56dOu/Fk4a+XUmTvgOAHfNfIgfX3wVK1e/w8zbfgPAGytW0adXj03r9+7VgydffKPZM7vwm5k1wRf22p3zLryc8392Bcd+fRhf/coX+c0l59a73chjDmPkMYfx8GPP8NNf/JYH77yamu6BnkWvJRd+M7MmGLDHrsyffTP3znqUH//8Kr4+/EDeeff9es/4qxxy0P7839I3WL3mXXr36s7rb7616bXlb75Fr169mz2zC7+ZWRO8uWIVO+/UmdEn/j3bb9+Rm267h//53S/r3OaVJa+ze7/eSOKZBS/x6afr6bLzDuy4wyAWL1nGq6+9wS49u3P7XbP4/R13NXtmF34zaz2OfLroh3x+0Sv88KL/pE0b0b5dO665rP7ul3+8Zw433zGT9u3b8bkOHbjj+v9AEu3atePKST/iyBPOpnJjJaefMoK999672TOrpjalcjNkyJBo9I1Y3J3TrNVaNOA+9urXtdQxstVlSKrVFi1axF577bXZMknzI2KrHbg7p5lZzmTa1CNpKfA+UAlsiIghknYG7gAqgKXAiRHxTpY5zMzsM8U44z8sIvar9ufGRGB2RPQHZifzZmZWJKVo6jkOmJpMTwWOL0EGM7PcyrrwB/AnSfMljU+W9YiIFQDJc/eMM5iZWTVZd+c8OCLelNQdmCXppbQbJh8U4wH69u2bVT4zs9zJtPBHxJvJ80pJdwFDgbck9YyIFZJ6Aitr2XYKMAUK3TmzzGlmrYO6puv6mFasrrsb+dJlb3LsKT/ghUfuaNbjzn1kPr+86hbuScbwaW6ZNfVI2k5Sp6pp4OvAC8AMYEyy2hig+YeeMzOzWmXZxt8DeETSAuApYGZE3A9MAo6QtBg4Ipk3M2uRNmyoZMxZF7HvISczauz5fPjh1ve+fe75l/nykWPZ95CTGXnaD3nn3bUADB9xJuf/7AqGHjGGAUO/xV8ef3az7TZu3Ej//v1ZtWrVpvk99tiD1atXNylzZoU/IpZExODksXdEXJIsXxMRh0dE/+T57awymJll7eVXXmP8aSNZ+PBtdO60HVff+N9brXPaWRdx6QX/xMKHb+MLe+3Ozy67btNrGzZs4KlZU5l8ybmbLQdo06YNo0eP5tZbbwXgwQcfZPDgwXTt2rSrlX3lrplZE/TZpQcHHzgYgNGjjuaRJxds9vp7az/g3ffe59CDDwBgzEnH8vBjn53Zf/PYrwFwwOCBLH19xVb7P/3007n55psBuPHGGxk7dmyTM7vwm5k1wZbj5Td0/Pxtt2kPQNu2bdmwoXKr1/v06UOPHj2YM2cOTz75JEcffXTjwyZc+M3MmmDZ8r/x+NMLAbjtzgcYlpz9V9mh8/bstGPnTe33v5t2L4cetH+DjjFu3DhGjx7NiSeeSNu2bZuc2cMymzXS8J8Xnuf+pLQ57DP1db/Mwl4D+jH19pmc+S+/oP9uffje2FFbrTP1ygv57nmT+PCjj9lt1134rysuaNAxRowYwdixY5ulmQdc+M3MGq2iby/++ti0etfb7wt78sQD/7XV8rkzrt003bXLjix9dgYAw4cdwPBhB2x6bcGCBQwePJiBAwc2Q2oXfjOzsjZp0iSuueaaTT17moPb+M3MytjEiRN57bXXGDZsWLPt04XfzFqu2EgLuIlg5hp6J0UXfjNrsTp88gpr3t+Q6+IfEaxZs4YOHTqk3sZt/GbWYvV+4yKWcxGrtt0D1ErPY1cuqneVDh060Lt379S7dOE3sxarfeU79Fs2odQxsnVK8/8500o/Is3MrDYu/GZmOePCb2aWMy78ZmY548Jv1gi3PgpPLIY/L4KKCYV5s5bChd+sgW59FMZfD59sKMy/trow7+JvLYULv1kD/ds0+PDTzZd9+GlhuVlL4MJv1kDLarndaW3LzcqNC79ZA/Wt5XantS03Kzcu/GYNdMmJ0HGbzZd13Kaw3KwlcOE3a6BTD4Yp42DbZMCTXbsW5k89uLS5zNLyWD1mjXDqwXDdQ4Vp33rRWhqf8ZuZ5YwLv5lZzrjwm5nljAu/mVnOuPCbmeVMq+/Vo1NzfDNOy9hwAHTq3JKmsNYtTmn+ffqM38wsZzIv/JLaSnpW0j3J/M6SZklanDzvlHUGMzP7TDHO+CcA1W8TPxGYHRH9gdnJvJmZFUmmhV9Sb+AY4Ppqi48DpibTU4Hjs8xgZmabq/XLXUl3A7V+MxoRI1LsfzLwI6BTtWU9ImJFso8VkrrXcvzxwHiAvn37pjiUmZmlUdcZ/y+BXwGvAh8B1yWPD4AX6tuxpGOBlRExvzHBImJKRAyJiCHdunVrzC7MzKwGtZ7xR8SfASRdHBGHVHvpbkkPp9j3wcAISX8PdAA6S7oFeEtSz+Rsvyewsgn5zcysgdK08XeTtFvVjKR+QL2n4BHx44joHREVwEnAnIgYDcwAxiSrjQGmNzi1mZk1WpoLuH4AzJW0JJmvAM5swjEnAdMknQEsA05owr7MzKyB6i38EXG/pP7AwGTRSxHxSUMOEhFzgbnJ9Brg8IbFNDOz5lJXr55v1vLS7pKIiDszymRmZhmq64z/G3W8FoALv5lZC1RXr56xxQxiZmbFUW+vHkk9JN0g6b5kflDyxayZmbVAabpz3gQ8APRK5v8XOCejPGZmlrE0hb9rREwDNgJExAagMtNUZmaWmTSFf52kLiTj9kj6MvBepqnMzCwzaS7gOpfC1ba7S3qUwlW7ozJNZWZmmUlzAdczkg4F9gQEvBwR6zNPZmZmmajrAq6vRcScGi7kGuALuMzMWq66zvgPBeZQ84VcvoDLzKyFqusCrguTyXER4V48ZmatRJpePa9IukzSoMzTmJlZ5tIU/n0pXLR1vaQnJI2X1DnjXGZmlpF6C39EvB8R10XEQRTun3shsELSVEl7ZJ7QrGzNTR5mLUuasXraShoh6S7gcgr34d0NuBu4N+N8ZmbWzNJcwLUYeAi4LCIeq7b8D5IOqWUbMzMrU2kK/74R8UFNL0TE2c2cx8zMMpamjb/Gom9mZi1Tml49ZmbWirjwm5nlTF1j9Zxb14YR8evmj2NmZlmr68vdTsnznsCXKAzNDIWxex7OMpSZmWWnrrF6fgYg6U/A/hHxfjJ/EfDfRUlnZmbNLk0bf1/g02rznwIVmaQxM7PMpenH/zvgqeTK3QBGAjdnmsrMzDKT5g5cl0i6D/hqsmhsRDybbSwzM8tK2u6cHYG1EXE5sFxSvwwzmZlZhtIM0nYhcD7w42RRe+CWLEOZmVl20pzxjwRGAOsAIuJNPuvqaWZmLUyawv9pRASFL3aRtF2aHUvqIOkpSQskvSipqnvozpJmSVqcPO/U+PhmZtZQaQr/NEnXAjtK+kfgQeD6FNt9AnwtIgYD+wFHSfoyMBGYHRH9gdnJvJmZFUmaXj2/lHQEsJbCVbwXRMSsFNsFUDWyZ/vkEcBxwPBk+VQKtzA6v6HBzcyscdJ8uXtpRMyKiB9GxHkRMUvSpWl2nty96zlgJTArIp4EekTECoDkuXst246XNE/SvFWrVqX+B5mZWd3SNPUcUcOyo9PsPCIqI2I/oDcwVNI+aYNFxJSIGBIRQ7p165Z2MzMzq0ddo3N+D/g+sLukhdVe6gQ8VvNWNYuIdyXNBY4C3pLUMyJWSOpJ4a8BMzMrkrrO+H9PYSTO6clz1eOAiDi1vh1L6iZpx2T6c8DfAS9RGOVzTLLamGT/ZmZWJHWNzvke8J6ky4G3q43O2UnSgUl7fV16AlMltaXwATMtIu6R9DiFnkJnAMuAE5rlX2JmZqmkGaTtGmD/avPrali2lYhYCHyxhuVrgMMbkNHMzJpRmi93lXTNBCAiNpLuA8PMzMpQmsK/RNLZktonjwnAkqyDmZlZNtIU/u8CBwFvAMuBA4HxWYYyM7PspLlydyVwUhGymJlZEaS5cneApNmSXkjm95X0k+yjmZlZFtI09VxHYSz+9bCpt47/AjAza6HSFP6OEfHUFss2ZBHGzMyyl6bwr5a0O5+Nxz8KWJFpKjMzy0ya/vhnAVOAgZLeAF4F6h2ywczMylOaXj1LgL9L7rzVpmroBjMza5nS9OrpIuk/gb8AcyVdLqlL9tHMzCwLadr4bwdWAd8CRiXTd2QZyszMspOmjX/niLi42vzPJR2fUR4zM8tYmjP+hySdJKlN8jgRmJl1MDMzy0aawn8mhZuyfJI8bgfOlfS+pLVZhjMzs+aXpldPp2IEMTOz4kjTq+eMLebbSrowu0hmZpalNE09h0u6V1JPSV8AnqBww3UzM2uB0jT1nCLp28DzwIfAyRHxaObJzMwsE2maevoDE4A/AkuBf5DUMeNcZmaWkTRNPXcDF0TEmcChwGLg6UxTmZlZZtJcwDU0ItYCJDdd/5WkGdnGMjOzrKQ54/+cpBsk3Q8gaRBwSLaxzMwsK2kK/03AA0DPZP5/gXMyymNmZhlLU/i7RsQ0YCNARGwAKjNNZWZmmUlT+NclwzBX3YHry8B7maYyM7PMpPly91xgBrC7pEeBbhSGZzYzsxYozQVcz0g6FNgTEPByRKzPPJmZmWUizRl/Vbv+ixlnMTOzIkjTxm9mZq1IZoVfUh9JD0laJOlFSROS5TtLmiVpcfK8U1YZzMxsa7U29Ujav64NI+KZeva9AfiX5DuCTsB8SbOA7wCzI2KSpInAROD8hsU2M7PGqquN/1d1vBbA1+racUSsAFYk0+9LWgTsAhwHDE9WmwrMxYXfzKxoai38EXFYcx1EUgXwReBJoEfyoUBErJDUvZZtxgPjAfr27dtcUczMci9Vrx5J+wCDgA5VyyLi5pTbbk9hSOdzImKtpFTBImIKMAVgyJAhkWojMzOrV72FP7nN4nAKhf9e4GjgEaDewi+pPYWif2tE3JksfktSz+RsvyewspHZzcysEdL06hkFHA78LSLGAoOBbevbSIVT+xuARRHx62ovzQDGJNNjgOkNSmxmZk2Spqnno4jYKGmDpM4UztB3S7HdwcA/AM9Lei5Z9q/AJGBachP3ZcAJDY9tZmaNlabwz5O0I3AdMB/4AHiqvo0i4hEKQzzU5PC0Ac3MrHmlGavn+8nkb5ObsXSOiIXZxjIzs6ykudn67KrpiFgaEQurLzMzs5alrit3OwAdga7JsApVzTadgV5FyGZmZhmoq6nnTAq3WOwFVB+eYS1wVYaZzMwsQ3VduXs5cLmkf46IK4qYyczMMpSmV8+1ks4GDknm5wLX+mYsZmYtU5rCfzXQPnmGQt/8a4BxWYUyM7Ps1PXlbrvkzltfiojB1V6aI2lB9tHMzCwLdXXnrLpIq1LS7lULJe0GVGaayszMMlNXU09V983zgIckLUnmK4CxWYYyM7Ps1FX4u0k6N5m+FmgLrKMwNPMXgYcyzmZmZhmoq/C3BbZn8/F2tk+eO2WWyMzMMlVX4V8REf9etCRmZlYUdX25m+5WWWZm1qLUVfg9dLKZWStUa+GPiLeLGcTMzIojza0XzcysFXHhNzPLGRd+M7OcceE3M8sZF34zs5xx4TczyxkXfjOznHHhNzPLGRd+M7OcceE3M8sZF34zs5xx4TczyxkXfjOznHHhNzPLmcwKv6QbJa2U9EK1ZTtLmiVpcfK8U1bHNzOzmmV5xn8TcNQWyyYCsyOiPzA7mTczsyLKrPBHxMPAljdzOQ6YmkxPBY7P6vhmZlazYrfx94iIFQDJc/faVpQ0XtI8SfNWrVpVtIBmZq1d2X65GxFTImJIRAzp1q1bqeOYmbUaxS78b0nqCZA8ryzy8c3Mcq/YhX8GMCaZHgNML/LxzcxyL8vunLcBjwN7Slou6QxgEnCEpMXAEcm8mZkVUbusdhwRJ9fy0uFZHdPMzOpXtl/umplZNlz4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McqYkhV/SUZJelvSKpImlyGBmlldFL/yS2gJXAUcDg4CTJQ0qdg4zs7wqxRn/UOCViFgSEZ8CtwPHlSCHmVkutSvBMXcBXq82vxw4cMuVJI0HxiezH0h6uZHH6wqsbuS2WXKuhnGuhnGuhinXXEhNyrZrTQtLUfhVw7LYakHEFGBKkw8mzYuIIU3dT3NzroZxroZxroYp11yQTbZSNPUsB/pUm+8NvFmCHGZmuVSKwv800F9SP0nbACcBM0qQw8wsl4re1BMRGyT9E/AA0Ba4MSJezPCQTW4uyohzNYxzNYxzNUy55oIMsiliq+Z1MzNrxXzlrplZzrjwm5nlTKsq/JLaSnpW0j01vDZQ0uOSPpF0XhnlOlXSwuTxmKTBZZLruCTTc5LmSRpWDrmqrfMlSZWSRpVDLknDJb2XvF/PSbqgHHJVy/acpBcl/bkcckn6YbX36oXkZ7lzGeTaQdLdkhYk79fYYmRKkWsnSXcl/yefkrRPU45Vin78WZoALAI61/Da28DZwPHFDJSoK9erwKER8Y6koyl8kbPVBW0lyDUbmBERIWlfYBowsAxyVQ37cSmFDgLFVGcu4C8RcWwR81SpNZekHYGrgaMiYpmk7uWQKyIuAy5LMn4D+EFEvF3qXMBZwF8j4huSugEvS7o1GWWglLn+FXguIkZKGkhh2JvDG3ugVnPGL6k3cAxwfU2vR8TKiHgaWF9muR6LiHeS2ScoXNdQDrk+iM+++d+OGi6yK0WuxD8DfwRWFiMTpM5VdClynQLcGRHLoPD/oExyVXcycFu2iQpS5AqgkyQB21M4YdxQBrkGUTgZIyJeAiok9Wjs8VpN4QcmAz8CNpY4x5Ymkz7XGcB9mab5zGTqySVppKSXgJnA6eWQS9IuwEjgt0XKU2Uy9f8cv5I0Edwnae/ixKo31wBgJ0lzJc2XdFqZ5AJAUkfgKAof5MUwmbpzXQnsReGi0ueBCRFRjJpSX64FwDcBJA2lMBRDo08SW0Xhl3QssDIi5pc6S3UNySXpMAqF//xyyRURd0XEQArNYxeXSa7JwPkRUZl1niopcz0D7BoRg4ErgP8pk1ztgAMonE0eCfxU0oAyyFXlG8CjxWjmSZnrSOA5oBewH3ClpNqa9oqZaxKFD/DnKPzF+yxN+UskIlr8A/gFhaEglgJ/Az4Ebqll3YuA88opF7Av8H/AgHLKtcU2rwJdS50rybE0eXxAobnn+FLnqmGbpWXyfk0ELqo2fwNwQqlzVVv3LuCULPM08P2aCXy12vwcYGipc22xvpJ1Ozf6mMV4w4v5AIYD99TxetEKf5pcQF/gFeCgcnq/gD347AK//YE3qubL4eeYrHMTMKpM3q/PV3u/hgLLyuH9otBsMZvCmX9H4AVgn1LnSl7bgUIb+nbF/BnW835dU/VBCfRIfu8z/QBPmWtHYJtk+h+Bm5tynNbWq2czkr4LEBG/lfR5YB6Fb8w3SjoHGBQRa0uZC7gA6AJcXfg+iQ1RolECt8j1LeA0SeuBj4BvR/JbV+JcZWOLXKOA70naQOH9Oqkc3q+IWCTpfmAhhfbj6yPihVLnShaNBP4UEetKkafKFrkuBm6S9DyFM+vzI6IkwzVvkWsv4GZJlcBfKTQLN37fJfrdNDOzEmkVX+6amVl6LvxmZjnjwm9mljMu/GZmOePCb2aWMy781qpIqpBUY3dFSddLGlTsTDXkqDWjWTG06n78ZtVFxLhSZ2gOktpFROYDh1nr5TN+a43aSZqajF3+h2QgMJKByoYk0x9IuiQZVO2JmkY6lHSRpBuT7ZZIOjtZvtkZu6TzJF1U7Ri/kfSwpEUq3DfgTkmLJf08RcYDJP05GVDtAUk9q+33P1QYT39CRu+b5YQLv7VGewJTImJfYC3w/RrW2Q54IgqDqj1M4TL4mgykMHDXUOBCSe1THP/TiDiEwgii0ymM8b4P8B1JXWrLmOz7CgrDUBwA3AhcUm2/O0bEoRHxqxQZzGrlwm+t0esR8WgyfQtQ093DPgWq7nQ0H6ioZV8zI+KT5LL9lRTGb6nPjOT5eeDFiFgREZ8AS4A+dWTck8IHxKxkFMafsPnQu3ekOLZZvdzGb63RluOQ1DQuyfpqY+lUUvv/hU+qTVett4HNT5o61LLNxi2231jtODVlFIUPiq/UkqWkY9pY6+EzfmuN+kqqKp4nA4808/7fArpL6iJpW6Axt1usKePLQLeq5ZLaF/GGLpYjLvzWGi0CxkhaCOxMYajdZhMR64F/B56k0Fz0UiN2s1XGKNzXdRRwqaQFFG4IclCzhDarxqNzmpnljM/4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxy5v8BuBtRklY7Ry0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sb4 = s30_and_b[4]\n", "b4 = b_only[4]\n", "n4 = data[4]\n", "plt.bar(x[4], sb4, color='orange', label='s=30')\n", "plt.bar(x[4], b4, color='b', yerr=np.sqrt(b_only[4]), label='b only')\n", "plt.scatter(x[4], n4, zorder=10, color='k')\n", "plt.xlabel('bin number')\n", "plt.ylabel('Total expected yield')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We introduce a new ingredient: the background level $b_4$ is now not a fixed number, but only known up to an uncertainty, say $\\sigma_b = \\pm 2$. We assume here that it is Gaussian-distributed, so that the full likelihood is\n", "$$\n", "L(s, b) = P(n_4; s_4 + b_4) G(b^{\\text{aux}}, b_4, \\sigma_b)\n", "$$\n", "\n", "where the first term is the Poisson PDF we had before, and the Gaussian is a new term giving the PDF for $b_4$. The central value of the $b_4$ constraint is $b^{\\text{aux}}$, which is a kind of observable: we can assume that our knowledge of $b_4$ comes from a separate \"auxiliary\" experiment which observed $b^{\\text{aux}}$ with an uncertainty $\\sigma_b$, and we use this to constrain $b_4$ in the experiment here. \n", "\n", "In practice the constraint on $b_4$ can also come from Monte-Carlo simulation, an independent theory prediction, etc., but in this case we can anyway use the same formalism, and set $b^{\\text{aux}}$ to the nominal prediction and $\\sigma_b$ to its uncertainty.\n", "\n", "Here we'll set $b^{\\text{aux}}$ to our original background value, which is what we take as nominal:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "baux = b_only[4]\n", "sigmab = 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define the profile likelihood, with $s$ as our parameter of interest, and $b$ the nuisance parameter. \n", "\n", "First we define $\\lambda(s,b) = -2\\log L(s,b)$, where $L(s,b)$ is the product of two terms:\n", "* The usual Poisson PDF for the one-bin counting process\n", "* The Gaussian representing the constraint on $b$, with an observable $b_{\\text{aux}}$, a mean $b$, and a width $\\sigma_b$.\n", "\n", "Secondly, we define the *profiled* version of the above, where $b$ is set to its best-fit value. To define it, we first must compute the best-fit $\\hat{\\hat {b}}(s)$ for a given $s$, by minimizing $\\lambda(s,b)$ for a fixed value of $s$. Then the profile likelihood is defined as $\\lambda(s) = -2\\log L(s,\\hat{\\hat {b}}(s))$.\n", "\n", "Finally, the $t(s)$ function is defined as usual as $t(s) = \\lambda(s) - \\lambda(\\hat{s})$, after computing the best-fit $\\hat{s}$ which minimizes $\\lambda(s)$. \n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD4CAYAAAAtrdtxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQS0lEQVR4nO3db4xcV33G8e+DTYQo0JDGAWM7dQC3qlu1JdqaVCCKCCDbRDat2ipRKRG8sKI2bWiLwGnUvkFq+SMBDY2IIopwRKoICRBWahRCmoqqUmjWCUkUTMhiATE2xKAKkPIitfLri7mmm2HWO3N2Zmcn+/1IVzv3nHPv/R3Ni2fvvXNnUlVIktTiOdMuQJI0uwwRSVIzQ0SS1MwQkSQ1M0QkSc02TruA1XbhhRfW9u3bp12GJM2Uo0eP/rCqNvW3r7sQ2b59O/Pz89MuQ5JmSpLvDGr3cpYkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqdnUQyTJ7iSPJllIcnBAf5Lc2PU/lOTSvv4NSR5IcsfqVS1JgimHSJINwE3AHmAncFWSnX3D9gA7uuUA8PG+/uuAYxMuVZI0wLTPRHYBC1V1vKqeAm4H9veN2Q/cWj33Aucn2QyQZCvwFuATq1m0JKln2iGyBXh80fqJrm3YMR8F3gM8fa6DJDmQZD7J/OnTp1dUsCTp/007RDKgrYYZk+QK4ImqOrrcQarqlqqaq6q5TZs2tdQpSRpg2iFyAti2aH0rcHLIMa8B9iX5Nr3LYG9I8unJlSpJ6jftELkP2JHkkiTnAVcCh/vGHAbe3n1K6zLgx1V1qqqur6qtVbW92+7fq+ptq1q9JK1zG6d58Ko6k+Ra4E5gA/DJqnokyTVd/83AEWAvsAA8CbxjWvVKkp4pVf23IJ7d5ubman5+ftplSNJMSXK0qub626d9OUuSNMMMEUlSM0NEktTMEJEkNTNEJEnNDBFJUjNDRJLUzBCRJDUzRCRJzQwRSVIzQ0SS1MwQkSQ1M0QkSc0MEUlSM0NEktTMEJEkNTNEJEnNDBFJUjNDRJLUzBCRJDUzRCRJzQwRSVIzQ0SS1MwQkSQ1M0QkSc0MEUlSM0NEktTMEJEkNTNEJEnNDBFJUrOph0iS3UkeTbKQ5OCA/iS5set/KMmlXfu2JPckOZbkkSTXrX71krS+TTVEkmwAbgL2ADuBq5Ls7Bu2B9jRLQeAj3ftZ4C/qapfAy4D/nzAtpKkCZr2mcguYKGqjlfVU8DtwP6+MfuBW6vnXuD8JJur6lRV3Q9QVT8FjgFbVrN4SVrvph0iW4DHF62f4OeDYNkxSbYDrwK+Ov4SJUlLmXaIZEBbjTImyQuAzwLvqqqfDDxIciDJfJL506dPNxcrSXqmaYfICWDbovWtwMlhxyR5Lr0Aua2qPrfUQarqlqqaq6q5TZs2jaVwSdL0Q+Q+YEeSS5KcB1wJHO4bcxh4e/cprcuAH1fVqSQB/gU4VlUfXt2yJUkAG6d58Ko6k+Ra4E5gA/DJqnokyTVd/83AEWAvsAA8Cbyj2/w1wJ8CDyf5Wtf2t1V1ZBWnIEnrWqr6b0E8u83NzdX8/Py0y5CkmZLkaFXN9bdP+3KWJGmGGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmg0dIkl+Iclzute/kmRfkudOrjRJ0lo3ypnIV4DnJdkC3E3vt84/NYmiJEmzYZQQSVU9CfwB8LGq+n1g52TKkiTNgpFCJMnvAn8C/FvXtnH8JUmSZsUoIXIdcD3w+ap6JMnLgXsmU5YkaRYMfSZRVV+hd1/k7Ppx4C/Prif5WFX9xXjLkyStZeP8iO9rxrgvSdIM8DkRSVIzQ0SS1GycIZIx7kuSNANGeWL9j5Zp+6exVCRJmhmjnIlcf662qvrUiquRJM2UZT/im2QPsBfYkuTGRV0vAs5MqjBJ0to3zHMiJ4F5YB9wdFH7T4G/mkRRkqTZsOzlrKp6sKoOAa+sqkOLls9V1f+cHZfksy0FJNmd5NEkC0kODuhPkhu7/oeSXDrstpKkyRr6nkhV/e8yQ14+6sGTbABuAvbQ+zLHq5L0f6njHmBHtxwAPj7CtpKkCRrnR3yrYZtdwEJVHa+qp4Dbgf19Y/YDt1bPvcD5STYPua0kaYKm/bDhFuDxResnurZhxgyzLQBJDiSZTzJ/+vTpFRctSeqZ9sOGg7bpP6NZasww2/Yaq26pqrmqmtu0adOIJUqSlrJsiCS5u/v7gWWGvrfh+CeAbYvWt9L7NNgwY4bZVpI0QcOciWxO8nvAviSvSnLp4uXsoKr6UsPx7wN2JLkkyXnAlcDhvjGHgbd3n9K6DPhxVZ0acltJ0gQN85zI3wMH6f2n/+G+vgLe0HrwqjqT5FrgTmAD8MnuB6+u6fpvBo7Qe9hxAXiS3m+7L7ltay2SpNGlargPVSX5u6p634Trmbi5ubman5+fdhmSNFOSHK2quf72UX7Z8H1J9gGv65r+o6ruGFeBkqTZM8q3+P4jvd9Z/3q3XNe1SZLWqaHPRIC3AL9dVU8DJDkEPMDgb/eVJK0Doz4ncv6i1784xjokSTNolDORfwAeSHIPvQf9XodnIZK0rg0VIkmeAzwNXAb8Dr0QeW9VfX+CtUmS1rihQqSqnk5ybVV9Bh/okyR1RrkncleSdyfZluSCs8vEKpMkrXmj3BN5J70n1P+sr33k3xGRJD07jBIiO+kFyGvphcl/AjdPoihJ0mwYJUQOAT8BbuzWr+ra/njcRUmSZsMoIfKrVfVbi9bvSfLguAuSJM2OUW6sP9B9FTsASV4N/Nf4S5IkzYpRzkReTe93Pb7brV8MHEvyMFBV9Ztjr06StKaNEiK7J1aFJGkmjfJV8N+ZZCGSpNkz6hcwSpL0M4aIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKnZ1EIkyQVJ7kryWPf3xUuM253k0SQLSQ4uav9Qkm8keSjJ55Ocv2rFS5KA6Z6JHATurqodwN3d+jMk2QDcBOwBdgJXJdnZdd8F/Eb3i4rfBK5flaolST8zzRDZDxzqXh8C3jpgzC5goaqOV9VTwO3ddlTVl6rqTDfuXmDrZMuVJPWbZoi8pKpOAXR/LxowZgvw+KL1E11bv3cCXxx7hZKkcxrlN9ZHluTLwEsHdN0w7C4GtFXfMW4AzgC3naOOA8ABgIsvvnjIQ0uSljPREKmqNy7Vl+QHSTZX1akkm4EnBgw7AWxbtL4VOLloH1cDVwCXV1WxhKq6BbgFYG5ubslxkqTRTPNy1mHg6u711cAXBoy5D9iR5JIk5wFXdtuRZDfwXmBfVT25CvVKkvpMM0TeD7wpyWPAm7p1krwsyRGA7sb5tcCdwDHgM1X1SLf9PwMvBO5K8rUkN6/2BCRpvZvo5axzqaofAZcPaD8J7F20fgQ4MmDcKydaoCRpWT6xLklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGZTC5EkFyS5K8lj3d8XLzFud5JHkywkOTig/91JKsmFk69akrTYNM9EDgJ3V9UO4O5u/RmSbABuAvYAO4Grkuxc1L8NeBPw3VWpWJL0DNMMkf3Aoe71IeCtA8bsAhaq6nhVPQXc3m131keA9wA1wTolSUuYZoi8pKpOAXR/LxowZgvw+KL1E10bSfYB36uqB5c7UJIDSeaTzJ8+fXrllUuSANg4yZ0n+TLw0gFdNwy7iwFtleT53T7ePMxOquoW4BaAubk5z1okaUwmGiJV9cal+pL8IMnmqjqVZDPwxIBhJ4Bti9a3AieBVwCXAA8mOdt+f5JdVfX9sU1AknRO07ycdRi4unt9NfCFAWPuA3YkuSTJecCVwOGqeriqLqqq7VW1nV7YXGqASNLqmmaIvB94U5LH6H3C6v0ASV6W5AhAVZ0BrgXuBI4Bn6mqR6ZUrySpz0QvZ51LVf0IuHxA+0lg76L1I8CRZfa1fdz1SZKW5xPrkqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRmhogkqZkhIklqZohIkpoZIpKkZoaIJKmZISJJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmqWqpl3DqkpyGvjOtOtocCHww2kXsYrW23zBOa8XszrnX66qTf2N6y5EZlWS+aqam3Ydq2W9zRec83rxbJuzl7MkSc0MEUlSM0Nkdtwy7QJW2XqbLzjn9eJZNWfviUiSmnkmIklqZohIkpoZImtEkguS3JXkse7vi5cYtzvJo0kWkhwc0P/uJJXkwslXvTIrnXOSDyX5RpKHknw+yfmrVvyIhnjfkuTGrv+hJJcOu+1a1TrnJNuS3JPkWJJHkly3+tWPbiXvcde/IckDSe5YvarHoKpc1sACfBA42L0+CHxgwJgNwLeAlwPnAQ8COxf1bwPupPcw5YXTntOk5wy8GdjYvf7AoO3XwrLc+9aN2Qt8EQhwGfDVYbddi8sK57wZuLR7/ULgm2t9ziuZ76L+vwb+Fbhj2vMZZfFMZO3YDxzqXh8C3jpgzC5goaqOV9VTwO3ddmd9BHgPMCuflljRnKvqS1V1pht3L7B1suU2W+59o1u/tXruBc5PsnnIbdei5jlX1amquh+gqn4KHAO2rGbxDVbyHpNkK/AW4BOrWfQ4GCJrx0uq6hRA9/eiAWO2AI8vWj/RtZFkH/C9qnpw0oWO0Yrm3Oed9P7LW4uGmcNSY4ad/1qzkjn/TJLtwKuAr46/xLFa6Xw/Su8fwKcnVN/EbJx2AetJki8DLx3QdcOwuxjQVkme3+3jza21Tcqk5tx3jBuAM8Bto1W3apadwznGDLPtWrSSOfc6kxcAnwXeVVU/GWNtk9A83yRXAE9U1dEkrx93YZNmiKyiqnrjUn1JfnD2VL47xX1iwLAT9O57nLUVOAm8ArgEeDDJ2fb7k+yqqu+PbQINJjjns/u4GrgCuLy6C8tr0DnnsMyY84bYdi1ayZxJ8lx6AXJbVX1ugnWOy0rm+4fAviR7gecBL0ry6ap62wTrHZ9p35Rx6S3Ah3jmTeYPDhizEThOLzDO3rz79QHjvs1s3Fhf0ZyB3cDXgU3Tnssy81z2faN3PXzxTdf/HuU9X2vLCucc4Fbgo9Oex2rMt2/M65mxG+tTL8CleyPgl4C7gce6vxd07S8Djiwat5fep1W+BdywxL5mJURWNGdggd415q91y83TntM55vpzcwCuAa7pXge4qet/GJgb5T1fi0vrnIHX0rsU9NCi93bvtOczyfd40T5mLkT82hNJUjM/nSVJamaISJKaGSKSpGaGiCSpmSEiSWpmiEiSmhkikqRm/wfhKzJSmZv60wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def logL_sb(s_val,b_val) : \n", " pass # -2 log L here\n", "\n", "def prof_logL_s(s_val) :\n", " # ==> Compute the best-fit bhat using scipy.optimize.minimize_scalar\n", " pass # return the profile likelihood lambda(s) here\n", "def prof_t_s(s_val) :\n", " # ==> find the best-fit shat that minimizes prof_logL_s\n", " pass # return t(s) = prof_logL_s(s) - prof_logL_s(shat) here instead\n", "\n", "s_vals = np.arange(-10,20,1)\n", "plt.plot(s_vals, [ prof_t_s(s_val) for s_val in s_vals ])\n", "\n", "plt.ylabel('prof_t_s');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This should look quite familiar from our previous exercises: we can extract the best-fit value of $s$ from the minimum value, and the uncertainties from the crossings with $t=1$. Just to check, we can compare with what we would get if we took a fixed background as before: " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZOklEQVR4nO3dfZiVdb3v8fd3A0a6bavA8aCo4MlUNEEdEQ+ctF0kCge3Wmdk25Nt9eKcdNc+x9TUSqt9ku1DZacCMjZaPu1MyXzKMMmoIRpElCeTB40JCoRLwfAB5Hv+WItpHOaGGWbWrIF5v65rXTO/+/e77/X9zXDxmXv97nWvyEwkSWrJ31S7AElS12VISJIKGRKSpEKGhCSpkCEhSSrUs9oFdKS+ffvmwIEDq12GJO1W5s6d+1Jm9mupb48KiYEDB1JfX1/tMiRptxIRLxb1+XKTJKmQISFJKmRISJIK7VFrEpJ2L5s3b6ahoYHXX3+92qV0C71792bAgAH06tWr1fsYEpKqpqGhgX333ZeBAwcSEdUuZ4+Wmaxbt46GhgYGDRrU6v18uUlS1bz++uv06dPHgOgEEUGfPn3afNZmSEiqKgOi8+zKz9qQkCQVMiQkqRNNnz6dRYsW7XTcpEmTuP3221t93JkzZzJ27Nj2lNYiF64lqZNs2bKF6dOnM3bsWAYPHrzDsRMmTOikqnbMMwlJ3dYLL7zAscce29i+8cYbufbaawE47bTTuOKKKxg2bBjvec97+NWvfgXAW2+9xWWXXcZ73/tejjvuOL71rW8BMHfuXE499VROPPFETj/9dFavXt14nKuuuopTTz2ViRMn8sADD/C5z32OoUOHsmzZMr73ve9x0kknMWTIEM4991w2bdoEwLXXXsuNN964w1qa27BhA2effTaDBw9mwoQJbN26td0/I88kJHUJ1/10IYtWbejQYw4+6F186b8fs8v7b9myhTlz5vDwww9z3XXXMWPGDKZMmcKKFSuYN28ePXv2ZP369WzevJlLL72Un/zkJ/Tr14977rmHq6++mqlTpwLw8ssv88tf/hKA559/nrFjx/LhD38YgP3224+LLroIgGuuuYbvf//7XHrppa2qpbk5c+awaNEiDjvsMEaPHs19993X+Dy7ypCQpALnnHMOACeeeCIvvPACADNmzGDChAn07Fn67/OAAw5gwYIFLFiwgFGjRgGls43+/fs3Hqe2trbwORYsWMA111zDyy+/zKuvvsrpp5/e6lqaGzZsGIcffjgA48ePZ9asWYaEpD1De/7i31U9e/Z820syzd9D8I53vAOAHj16sGXLFqD0prTml5JmJscccwx1dXUtPs8+++xTWMMnP/lJpk+fzpAhQ5g2bRozZ85scVxLtTTXvK6OuLzYNQlJ3daBBx7ImjVrWLduHW+88QYPPvjgTvf50Ic+xKRJkxr/o16/fj1HHnkka9eubQyJzZs3s3Dhwhb333fffdm4cWNje+PGjfTv35/Nmzdzxx13tGs+c+bMYcWKFWzdupV77rmHkSNHtut4YEhI6sZ69erFF7/4RU4++WTGjh3LUUcdtdN9LrzwQg499FCOO+44hgwZwp133slee+3FvffeyxVXXMGQIUMYOnQov/nNb1rc/7zzzuOGG27g+OOPZ9myZXzlK1/h5JNPZtSoUa16/h055ZRTuPLKKzn22GMZNGgQZ599druOBxCZ2e6DdBU1NTXphw5Ju4/Fixdz9NFHV7uMbqWln3lEzM3MmpbGeyYhSSpkSEiSChkSkqRChoQkqZAhIUkqZEhIkgoZEpK6tVtuuYWjjz6a888/nwceeIDrr7++3ccsum33tGnTuOSSS9p9/M7kbTkkdWvf+c53eOSRRxo/93ncuHFVrqhrqeiZRERMjYg1EbGgoD8i4paIWBoRz0TECc36e0TEvIjY+XvlJamNJkyYwPLlyxk3bhxf//rX3/aX/llnndX4oT+TJ0/m/PPPB+Cxxx7jlFNO4YQTTuAjH/kIr776KgCPPvooRx11FCNHjuS+++4rfM6VK1cyevRojjzySK677roKz7D9Kn0mMQ34f0DRxyudARxRfpwMfLf8dZvPAIuBd1WuREldxr+P2X7bMf8Awy6CNzfBHR/Zvn/oP8Lx58Nf1sF/fPztfRc8tMOnmzRpEo8++ihPPPEEffv2Zdq0aY19U6ZMYcSIEQwaNIibbrqJ2bNn89JLL/HVr36VGTNmsM8++zBx4kRuvvlmLr/8ci666CJ+8Ytf8O53v3uHd32dM2cOCxYsYO+99+akk05izJgx1NS0+GbnLqGiZxKZ+SSwfgdDzgJuz5LZwH4R0R8gIgYAY4BbK1mjJLXkwAMP5Mtf/jLvf//7uemmmzjggAOYPXs2ixYtYsSIEQwdOpTbbruNF198kSVLljBo0CCOOOIIIoKPfvSjhccdNWoUffr04Z3vfCfnnHMOs2bN6sRZtV211yQOBlY2aTeUt60GvgFcDuy7owNExMXAxQCHHnpoRYqU1El29Jf/XnvvuH+fPjs9c2irZ599lj59+rBq1SqgdEvwUaNGcdddd71t3NNPP93q23JX4nbelVTtq5ta+ulkRIwF1mTm3J0dIDOnZGZNZtb069ev4yuU1C3NmTOHRx55hHnz5nHjjTeyYsUKhg8fzq9//WuWLl0KwKZNm/j973/PUUcdxYoVK1i2bBnAdiHS1M9//nPWr1/Pa6+9xvTp0xkxYkSnzGdXVTskGoBDmrQHAKuAEcC4iHgBuBv4+4j4YeeXJ6k7euONN7jooouYOnUqBx10EDfddBOf+tSnGtctxo8fz3HHHcfw4cNZsmQJvXv3ZsqUKYwZM4aRI0dy2GGHFR575MiRfOxjH2Po0KGce+65XXo9AjrhVuERMRB4MDOPbaFvDHAJcCalBetbMnNYszGnAZdl5vYXHTfjrcKl3Yu3Cu98bb1VeEXXJCLiLuA0oG9ENABfAnoBZOYk4GFKAbEU2ARcUMl6JEltU9GQyMzxO+lP4NM7GTMTmNlxVUmSWqvaaxKSurk96dMxu7pd+VkbEpKqpnfv3qxbt86g6ASZybp16+jdu3eb9qv2+yQkdWMDBgygoaGBtWvXVruUbqF3794MGDCgTfsYEpKqplevXo031lPX5MtNkqRChoQkqZAhIUkqZEhIkgoZEpKkQoaEJKmQISFJKmRISJIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqZAhIUkqZEhIkgoZEpKkQoaEJKmQISFJKmRISJIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqVBFQyIipkbEmohYUNAfEXFLRCyNiGci4oTy9kMi4omIWBwRCyPiM5WsU5LUskqfSUwDRu+g/wzgiPLjYuC75e1bgP+TmUcDw4FPR8TgCtYpSWpBRUMiM58E1u9gyFnA7VkyG9gvIvpn5urMfKp8jI3AYuDgStYqSdpetdckDgZWNmk30CwMImIgcDzw284rS5IE1Q+JaGFbNnZG/C3wY+CzmbmhxQNEXBwR9RFRv3bt2gqVKUndU7VDogE4pEl7ALAKICJ6UQqIOzLzvqIDZOaUzKzJzJp+/fpVtFhJ6m6qHRIPAB8vX+U0HHglM1dHRADfBxZn5s3VLVGSuq+elTx4RNwFnAb0jYgG4EtAL4DMnAQ8DJwJLAU2AReUdx0BfAx4NiKeLm+7KjMfrmS9kqS3q2hIZOb4nfQn8OkWts+i5fUKSVInqvbLTZKkLsyQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklSo1SEREftExN+Uv39PRIyLiF6VK02SVG1tOZN4EugdEQcDj1P6POpplShKktQ1tCUkIjM3AecA38rMs4HBlSlLktQVtCkkIuIU4HzgofK2nh1fkiSpq2hLSHwG+Dxwf2YujIjDgScqU5YkqSto9ZlAZj5JaV1iW3s58M/b2hHxrcy8tGPLkyRVU0deAjuiA48lSeoCfJ+EJKmQISFJKtSRIREdeCxJUhfQlndcf2Qn277ZIRVJkrqMtpxJfH5H2zJzWrurkSR1KTu9BDYizgDOBA6OiFuadL0L2FKpwiRJ1dea90msAuqBccDcJts3Av9SiaIkSV3DTkMiM+cD8yPizszcXDQuIn6cmec22zYVGAusycxjW9gnKK1lnAlsAj6ZmU+V+0aX+3oAt2bm9a2fliSpI7R6TWJHAVF2eAvbpgGjd7DPGcAR5cfFwHcBIqIH8O1y/2BgfER4M0FJ6mQdeQlsbrehdCuP9TvY5yzg9iyZDewXEf2BYcDSzFyemW8Cd5fH7tDytX/hR/UrAdj81lZqJ9dx/7wGAF578y1qJ9fx0/mrANjw+mZqJ9fx6ILVAKz/y5vUTq5jxqI/A7Bm4+vUTq5j5nNrAFj18mvUTq5j1vMvAfCHdZuonVzH7OXrAFi29lVqJ9cx98XSdJ/700ZqJ9cxf+XLACxc9Qq1k+tYuOoVAOavfJnayXU896eNAMx9cT21k+tYtvZVAGYvX0ft5Dr+sG4TALOef4nayXWsevk1AGY+t4bayXWs2fg6ADMW/ZnayXWs/8ubADy6YDW1k+vY8Hop2386fxW1k+t47c23ALh/XgO1k+vY/NZWAH5Uv5LayXWNP8u75vyB82+d3dj+Qd0LfGLqnMb21FkruPC23zW2pzy5jAk/+Ourkd+ZuZRL7nyqsX3L48/z2bvnNbZvfuw5LvvR/Mb2xEeX8Pn7nmls/+tDi/jC9AWN7et+upDrfrqwsf2F6Qv414cWNbY/f98zTHx0SWP7sh/N5+bHnmtsf/buedzy+PON7UvufIrvzFza2J7wg7lMeXJZY/vC237H1FkrGtufmDqHH9S90Ng+/9bZ3DXnD43t2sl1/tvz3x6we/7b25Fqv5nuYGBlk3ZDeVvR9u1ExMURUR8R9Zs37+xkR5LUFpG53QnArh0oYl5mHt/C9oHAgwVrEg8BX8vMWeX248DllF66Oj0zLyxv/xgwbGc3EKypqcn6+vp2z0WSupOImJuZNS317fRMovwfNxExcSdDr9iF2hqAQ5q0B1C6mqpouySpE7XmEtj+EXEqMC4i7qbZ7Te2XY2UmY/twvM/AFxSPu7JwCuZuToi1gJHRMQg4I/AecA/7sLxJUnt0JqQ+CJwJaW/5m9u1pfA3xftGBF3AacBfSOiAfgS0AsgMycBD1O6/HUppUtgLyj3bYmIS4CfUboEdmpmLtzuCSRJFdXqNYmI+EJmfqXC9bSLaxKS1HY7WpNoyyfTfSUixgHvK2+amZkPdkSBkqSuqS13gf0apc+5XlR+fKa8TZK0h2r1mQQwBhiamVsBIuI2YB4t3x1WkrQHaOub6fZr8v3fdWAdkqQuqC1nEv8XmBcRT1C6DPZ9eBYhSXu0VoVERPwNsBUYDpxEKSSuyMw/VbA2SVKVtSokMnNrRFySmf9B6Q1wkqRuoC1rEj+PiMsi4pCIOGDbo2KVSZKqri1rEp+i9A7r/9Vse0ufIyFJ2gO0JSQGUwqIkZTC4lfApEoUJUnqGtoSErcBG4Bbyu3x5W3/o6OLkiR1DW0JiSMzc0iT9hMRMb9wtCRpt9eWhet5ETF8WyMiTgZ+3fElSZK6iracSZwMfDwitn246qHA4oh4FsjMPK7Dq5MkVVVbQmJ0xaqQJHVJbblV+IuVLESS1PW09QZ/kqRuxJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBWqeEhExOiIeC4ilkbElS307x8R90fEMxExJyKObdL3LxGxMCIWRMRdEdG70vVKkv6qoiERET2AbwNnAIOB8RExuNmwq4Cny59s93Hgm+V9Dwb+GajJzGOBHsB5laxXkvR2lT6TGAYszczlmfkmcDdwVrMxg4HHATJzCTAwIg4s9/UE3hkRPYG9gVUVrleS1ESlQ+JgYGWTdkN5W1PzgXMAImIYcBgwIDP/CNwI/AFYDbySmY9VuF5JUhOVDoloYVs2a18P7B8RTwOXAvOALRGxP6WzjkHAQcA+EfHR7Z4g4uKIqI+I+rVr13Zo8ZLU3VU6JBqAQ5q0B9DsJaPM3JCZF2TmUEprEv2AFcAHgRWZuTYzNwP3Af+1+RNk5pTMrMnMmn79+lVoGpLUPVU6JH4HHBERgyJiL0oLzw80HRAR+5X7AC4EnszMDZReZhoeEXtHRAAfABZXuF5JUhM9K3nwzNwSEZcAP6N0ddLUzFwYERPK/ZOAo4HbI+ItYBHwT+W+30bEvcBTwBZKL0NNqWS9kqS3i8zmSwS7r5qamqyvr692GZK0W4mIuZlZ01Kf77iWJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSoYqHRESMjojnImJpRFzZQv/+EXF/RDwTEXMi4tgmfftFxL0RsSQiFkfEKZWuV5L0VxUNiYjoAXwbOAMYDIyPiMHNhl0FPJ2ZxwEfB77ZpO+bwKOZeRQwBFhcyXolSW9X6TOJYcDSzFyemW8CdwNnNRszGHgcIDOXAAMj4sCIeBfwPuD75b43M/PlCtcrSWqi0iFxMLCySbuhvK2p+cA5ABExDDgMGAAcDqwF/j0i5kXErRGxT/MniIiLI6I+IurXrl1biTlIUrdV6ZCIFrZls/b1wP4R8TRwKTAP2AL0BE4AvpuZxwN/AbZb08jMKZlZk5k1/fr168jaJanb61nh4zcAhzRpDwBWNR2QmRuACwAiIoAV5cfeQENm/rY89F5aCAlJUuVU+kzid8ARETEoIvYCzgMeaDqgfAXTXuXmhcCTmbkhM/8ErIyII8t9HwAWVbheSVITFT2TyMwtEXEJ8DOgBzA1MxdGxIRy/yTgaOD2iHiLUgj8U5NDXArcUQ6R5ZTPOCRJnSMymy8R7L5qamqyvr6+2mVI0m4lIuZmZk1Lfb7jWpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUqHIzGrX0GEiYi3wYrXr2AV9gZeqXUQnc87dg3PePRyWmf1a6tijQmJ3FRH1mVlT7To6k3PuHpzz7s+XmyRJhQwJSVIhQ6JrmFLtAqrAOXcPznk355qEJKmQZxKSpEKGhCSpkCHRSSLigIj4eUQ8X/66f8G40RHxXEQsjYgrW+i/LCIyIvpWvur2ae+cI+KGiFgSEc9ExP0RsV+nFd8GrfidRUTcUu5/JiJOaO2+XdWuzjkiDomIJyJicUQsjIjPdH71u6Y9v+dyf4+ImBcRD3Ze1R0gM310wgP4N+DK8vdXAhNbGNMDWAYcDuwFzAcGN+k/BPgZpTcM9q32nCo9Z+BDQM/y9xNb2r/aj539zspjzgQeAQIYDvy2tft2xUc759wfOKH8/b7A7/f0OTfp/9/AncCD1Z5PWx6eSXSes4Dbyt/fBvxDC2OGAUszc3lmvgncXd5vm68DlwO7y9UG7ZpzZj6WmVvK42YDAypb7i7Z2e+Mcvv2LJkN7BcR/Vu5b1e0y3POzNWZ+RRAZm4EFgMHd2bxu6g9v2ciYgAwBri1M4vuCIZE5zkwM1cDlL/+pxbGHAysbNJuKG8jIsYBf8zM+ZUutAO1a87NfIrSX2ldTWvqLxrT2rl3Ne2Zc6OIGAgcD/y240vscO2d8zco/YG3tUL1VUzPahewJ4mIGcB/bqHr6tYeooVtGRF7l4/xoV2trVIqNedmz3E1sAW4o23VdYqd1r+DMa3Ztytqz5xLnRF/C/wY+GxmbujA2ipll+ccEWOBNZk5NyJO6+jCKs2Q6ECZ+cGivoj487bT7fIp6JoWhjVQWnfYZgCwCvgvwCBgfkRs2/5URAzLzD912AR2QQXnvO0YnwDGAh/I8gu7XcwO69/JmL1asW9X1J45ExG9KAXEHZl5XwXr7EjtmfOHgXERcSbQG3hXRPwwMz9awXo7TrUXRbrLA7iBty/i/lsLY3oCyykFwrbFsWNaGPcCu8fCdbvmDIwGFgH9qj2XHcxxp78zSq9FN13QnNOW33dXe7RzzgHcDnyj2vPorDk3G3Mau9nCddUL6C4PoA/wOPB8+esB5e0HAQ83GXcmpSs+lgFXFxxrdwmJds0ZWErpNd6ny49J1Z5TwTy3qx+YAEwofx/At8v9zwI1bfl9d8XHrs4ZGEnpZZpnmvxez6z2fCr9e25yjN0uJLwthySpkFc3SZIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqZAhIVVYROwTEQ9FxPyIWBARtdWuSWot790kVd5oYFVmjgGIiL+rcj1Sq3kmIVXes8AHI2JiRPy3zHyl2gVJrWVISBWWmb8HTqQUFl+LiC9WuSSp1Xy5SaqwiDgIWJ+ZP4yIV4FPVrkkqdUMCany3gvcEBFbgc3A/6xyPVKreRdYSVIh1yQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJU6P8DGXcX0cIeNIYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def t_sb(s, b) :\n", " # ==> compute the best-fit s by minimizing logL_sb(s, b) at a fixed b\n", " pass # return t(s,b) = logL_sb(s,b) - logL_sb(shat,b) here instead\n", "\n", "s_vals = np.arange(-10,20,1)\n", "plt.plot(s_vals, [ prof_t_s(s_val) for s_val in s_vals ], label='uncertain b')\n", "plt.plot(s_vals, [ t_sb(s_val, b4) for s_val in s_vals ], linestyle='--', label='fixed b')\n", "plt.xlabel('s')\n", "plt.ylabel('prof_t_s')\n", "plt.axhline(1, linestyle=':')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hopefully everything should work as expected : we get the same result as before, but the size of the confidence interval went up since now the background adds a new source of uncertainty. As before, we could compute the crossings with $t=1$ to get the precise values of the uncertainties in each case.\n", "\n", "To see a bit better what is happening, we can check the values of the best-fit $b$ for different $s$ values:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "b4 = 41.29032258064516\n" ] } ], "source": [ "print('b4 = ', b4)\n", "# ==> compute the best-fit b value for s = 2.7, -10 and 20." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the best-fit point $s=2.7$, the best-fit background should be equal to the nominal one. However at a lower $s$ value $s=-10$, you should get a higher best-fit $b$: this reflects the fact that the background is pushed slightly up to compensate for the fact that the assumed signal is lower than what the data prefers. \n", "\n", "And if we push $s$ up, you should see $b$ going down for the same reason.\n", "\n", "These adjustments makes it more difficult to exclude values of $s$ away from $\\hat{s}$, so that the allowed range of $s$ gets larger -- and therefore the uncertainties get larger as well." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Profiling the multi-bin example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the single-bin example, we could only release $b$ if we also introduced an external constraint in the measurement. This is because the single-bin example only measures one number, so it cannot be used to determine both $s$ and $b$.\n", "\n", "In multi-bin analyses, we count events in multiple bins so it is possible to simultaneously measure multiple parameters. In our example, we can measure $s$ in the central, signal-rich region, while extracting $b$ from the side regions where there is no signal.\n", "\n", "So we can keep the same likelihood definition as in the previous lecture (without the Gaussian term) and simply make $b$ a free (nuisance) parameter." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZOklEQVR4nO3dfZiVdb3v8fd3A0a6bavA8aCo4MlUNEEdEQ+ctF0kCge3Wmdk25Nt9eKcdNc+x9TUSqt9ku1DZacCMjZaPu1MyXzKMMmoIRpElCeTB40JCoRLwfAB5Hv+WItpHOaGGWbWrIF5v65rXTO/+/e77/X9zXDxmXv97nWvyEwkSWrJ31S7AElS12VISJIKGRKSpEKGhCSpkCEhSSrUs9oFdKS+ffvmwIEDq12GJO1W5s6d+1Jm9mupb48KiYEDB1JfX1/tMiRptxIRLxb1+XKTJKmQISFJKmRISJIK7VFrEpJ2L5s3b6ahoYHXX3+92qV0C71792bAgAH06tWr1fsYEpKqpqGhgX333ZeBAwcSEdUuZ4+Wmaxbt46GhgYGDRrU6v18uUlS1bz++uv06dPHgOgEEUGfPn3afNZmSEiqKgOi8+zKz9qQkCQVMiQkqRNNnz6dRYsW7XTcpEmTuP3221t93JkzZzJ27Nj2lNYiF64lqZNs2bKF6dOnM3bsWAYPHrzDsRMmTOikqnbMMwlJ3dYLL7zAscce29i+8cYbufbaawE47bTTuOKKKxg2bBjvec97+NWvfgXAW2+9xWWXXcZ73/tejjvuOL71rW8BMHfuXE499VROPPFETj/9dFavXt14nKuuuopTTz2ViRMn8sADD/C5z32OoUOHsmzZMr73ve9x0kknMWTIEM4991w2bdoEwLXXXsuNN964w1qa27BhA2effTaDBw9mwoQJbN26td0/I88kJHUJ1/10IYtWbejQYw4+6F186b8fs8v7b9myhTlz5vDwww9z3XXXMWPGDKZMmcKKFSuYN28ePXv2ZP369WzevJlLL72Un/zkJ/Tr14977rmHq6++mqlTpwLw8ssv88tf/hKA559/nrFjx/LhD38YgP3224+LLroIgGuuuYbvf//7XHrppa2qpbk5c+awaNEiDjvsMEaPHs19993X+Dy7ypCQpALnnHMOACeeeCIvvPACADNmzGDChAn07Fn67/OAAw5gwYIFLFiwgFGjRgGls43+/fs3Hqe2trbwORYsWMA111zDyy+/zKuvvsrpp5/e6lqaGzZsGIcffjgA48ePZ9asWYaEpD1De/7i31U9e/Z820syzd9D8I53vAOAHj16sGXLFqD0prTml5JmJscccwx1dXUtPs8+++xTWMMnP/lJpk+fzpAhQ5g2bRozZ85scVxLtTTXvK6OuLzYNQlJ3daBBx7ImjVrWLduHW+88QYPPvjgTvf50Ic+xKRJkxr/o16/fj1HHnkka9eubQyJzZs3s3Dhwhb333fffdm4cWNje+PGjfTv35/Nmzdzxx13tGs+c+bMYcWKFWzdupV77rmHkSNHtut4YEhI6sZ69erFF7/4RU4++WTGjh3LUUcdtdN9LrzwQg499FCOO+44hgwZwp133slee+3FvffeyxVXXMGQIUMYOnQov/nNb1rc/7zzzuOGG27g+OOPZ9myZXzlK1/h5JNPZtSoUa16/h055ZRTuPLKKzn22GMZNGgQZ599druOBxCZ2e6DdBU1NTXphw5Ju4/Fixdz9NFHV7uMbqWln3lEzM3MmpbGeyYhSSpkSEiSChkSkqRChoQkqZAhIUkqZEhIkgoZEpK6tVtuuYWjjz6a888/nwceeIDrr7++3ccsum33tGnTuOSSS9p9/M7kbTkkdWvf+c53eOSRRxo/93ncuHFVrqhrqeiZRERMjYg1EbGgoD8i4paIWBoRz0TECc36e0TEvIjY+XvlJamNJkyYwPLlyxk3bhxf//rX3/aX/llnndX4oT+TJ0/m/PPPB+Cxxx7jlFNO4YQTTuAjH/kIr776KgCPPvooRx11FCNHjuS+++4rfM6VK1cyevRojjzySK677roKz7D9Kn0mMQ34f0DRxyudARxRfpwMfLf8dZvPAIuBd1WuREldxr+P2X7bMf8Awy6CNzfBHR/Zvn/oP8Lx58Nf1sF/fPztfRc8tMOnmzRpEo8++ihPPPEEffv2Zdq0aY19U6ZMYcSIEQwaNIibbrqJ2bNn89JLL/HVr36VGTNmsM8++zBx4kRuvvlmLr/8ci666CJ+8Ytf8O53v3uHd32dM2cOCxYsYO+99+akk05izJgx1NS0+GbnLqGiZxKZ+SSwfgdDzgJuz5LZwH4R0R8gIgYAY4BbK1mjJLXkwAMP5Mtf/jLvf//7uemmmzjggAOYPXs2ixYtYsSIEQwdOpTbbruNF198kSVLljBo0CCOOOIIIoKPfvSjhccdNWoUffr04Z3vfCfnnHMOs2bN6sRZtV211yQOBlY2aTeUt60GvgFcDuy7owNExMXAxQCHHnpoRYqU1El29Jf/XnvvuH+fPjs9c2irZ599lj59+rBq1SqgdEvwUaNGcdddd71t3NNPP93q23JX4nbelVTtq5ta+ulkRIwF1mTm3J0dIDOnZGZNZtb069ev4yuU1C3NmTOHRx55hHnz5nHjjTeyYsUKhg8fzq9//WuWLl0KwKZNm/j973/PUUcdxYoVK1i2bBnAdiHS1M9//nPWr1/Pa6+9xvTp0xkxYkSnzGdXVTskGoBDmrQHAKuAEcC4iHgBuBv4+4j4YeeXJ6k7euONN7jooouYOnUqBx10EDfddBOf+tSnGtctxo8fz3HHHcfw4cNZsmQJvXv3ZsqUKYwZM4aRI0dy2GGHFR575MiRfOxjH2Po0KGce+65XXo9AjrhVuERMRB4MDOPbaFvDHAJcCalBetbMnNYszGnAZdl5vYXHTfjrcKl3Yu3Cu98bb1VeEXXJCLiLuA0oG9ENABfAnoBZOYk4GFKAbEU2ARcUMl6JEltU9GQyMzxO+lP4NM7GTMTmNlxVUmSWqvaaxKSurk96dMxu7pd+VkbEpKqpnfv3qxbt86g6ASZybp16+jdu3eb9qv2+yQkdWMDBgygoaGBtWvXVruUbqF3794MGDCgTfsYEpKqplevXo031lPX5MtNkqRChoQkqZAhIUkqZEhIkgoZEpKkQoaEJKmQISFJKmRISJIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqZAhIUkqZEhIkgoZEpKkQoaEJKmQISFJKmRISJIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqVBFQyIipkbEmohYUNAfEXFLRCyNiGci4oTy9kMi4omIWBwRCyPiM5WsU5LUskqfSUwDRu+g/wzgiPLjYuC75e1bgP+TmUcDw4FPR8TgCtYpSWpBRUMiM58E1u9gyFnA7VkyG9gvIvpn5urMfKp8jI3AYuDgStYqSdpetdckDgZWNmk30CwMImIgcDzw284rS5IE1Q+JaGFbNnZG/C3wY+CzmbmhxQNEXBwR9RFRv3bt2gqVKUndU7VDogE4pEl7ALAKICJ6UQqIOzLzvqIDZOaUzKzJzJp+/fpVtFhJ6m6qHRIPAB8vX+U0HHglM1dHRADfBxZn5s3VLVGSuq+elTx4RNwFnAb0jYgG4EtAL4DMnAQ8DJwJLAU2AReUdx0BfAx4NiKeLm+7KjMfrmS9kqS3q2hIZOb4nfQn8OkWts+i5fUKSVInqvbLTZKkLsyQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklSo1SEREftExN+Uv39PRIyLiF6VK02SVG1tOZN4EugdEQcDj1P6POpplShKktQ1tCUkIjM3AecA38rMs4HBlSlLktQVtCkkIuIU4HzgofK2nh1fkiSpq2hLSHwG+Dxwf2YujIjDgScqU5YkqSto9ZlAZj5JaV1iW3s58M/b2hHxrcy8tGPLkyRVU0deAjuiA48lSeoCfJ+EJKmQISFJKtSRIREdeCxJUhfQlndcf2Qn277ZIRVJkrqMtpxJfH5H2zJzWrurkSR1KTu9BDYizgDOBA6OiFuadL0L2FKpwiRJ1dea90msAuqBccDcJts3Av9SiaIkSV3DTkMiM+cD8yPizszcXDQuIn6cmec22zYVGAusycxjW9gnKK1lnAlsAj6ZmU+V+0aX+3oAt2bm9a2fliSpI7R6TWJHAVF2eAvbpgGjd7DPGcAR5cfFwHcBIqIH8O1y/2BgfER4M0FJ6mQdeQlsbrehdCuP9TvY5yzg9iyZDewXEf2BYcDSzFyemW8Cd5fH7tDytX/hR/UrAdj81lZqJ9dx/7wGAF578y1qJ9fx0/mrANjw+mZqJ9fx6ILVAKz/y5vUTq5jxqI/A7Bm4+vUTq5j5nNrAFj18mvUTq5j1vMvAfCHdZuonVzH7OXrAFi29lVqJ9cx98XSdJ/700ZqJ9cxf+XLACxc9Qq1k+tYuOoVAOavfJnayXU896eNAMx9cT21k+tYtvZVAGYvX0ft5Dr+sG4TALOef4nayXWsevk1AGY+t4bayXWs2fg6ADMW/ZnayXWs/8ubADy6YDW1k+vY8Hop2386fxW1k+t47c23ALh/XgO1k+vY/NZWAH5Uv5LayXWNP8u75vyB82+d3dj+Qd0LfGLqnMb21FkruPC23zW2pzy5jAk/+Ourkd+ZuZRL7nyqsX3L48/z2bvnNbZvfuw5LvvR/Mb2xEeX8Pn7nmls/+tDi/jC9AWN7et+upDrfrqwsf2F6Qv414cWNbY/f98zTHx0SWP7sh/N5+bHnmtsf/buedzy+PON7UvufIrvzFza2J7wg7lMeXJZY/vC237H1FkrGtufmDqHH9S90Ng+/9bZ3DXnD43t2sl1/tvz3x6we/7b25Fqv5nuYGBlk3ZDeVvR9u1ExMURUR8R9Zs37+xkR5LUFpG53QnArh0oYl5mHt/C9oHAgwVrEg8BX8vMWeX248DllF66Oj0zLyxv/xgwbGc3EKypqcn6+vp2z0WSupOImJuZNS317fRMovwfNxExcSdDr9iF2hqAQ5q0B1C6mqpouySpE7XmEtj+EXEqMC4i7qbZ7Te2XY2UmY/twvM/AFxSPu7JwCuZuToi1gJHRMQg4I/AecA/7sLxJUnt0JqQ+CJwJaW/5m9u1pfA3xftGBF3AacBfSOiAfgS0AsgMycBD1O6/HUppUtgLyj3bYmIS4CfUboEdmpmLtzuCSRJFdXqNYmI+EJmfqXC9bSLaxKS1HY7WpNoyyfTfSUixgHvK2+amZkPdkSBkqSuqS13gf0apc+5XlR+fKa8TZK0h2r1mQQwBhiamVsBIuI2YB4t3x1WkrQHaOub6fZr8v3fdWAdkqQuqC1nEv8XmBcRT1C6DPZ9eBYhSXu0VoVERPwNsBUYDpxEKSSuyMw/VbA2SVKVtSokMnNrRFySmf9B6Q1wkqRuoC1rEj+PiMsi4pCIOGDbo2KVSZKqri1rEp+i9A7r/9Vse0ufIyFJ2gO0JSQGUwqIkZTC4lfApEoUJUnqGtoSErcBG4Bbyu3x5W3/o6OLkiR1DW0JiSMzc0iT9hMRMb9wtCRpt9eWhet5ETF8WyMiTgZ+3fElSZK6iracSZwMfDwitn246qHA4oh4FsjMPK7Dq5MkVVVbQmJ0xaqQJHVJbblV+IuVLESS1PW09QZ/kqRuxJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBWqeEhExOiIeC4ilkbElS307x8R90fEMxExJyKObdL3LxGxMCIWRMRdEdG70vVKkv6qoiERET2AbwNnAIOB8RExuNmwq4Cny59s93Hgm+V9Dwb+GajJzGOBHsB5laxXkvR2lT6TGAYszczlmfkmcDdwVrMxg4HHATJzCTAwIg4s9/UE3hkRPYG9gVUVrleS1ESlQ+JgYGWTdkN5W1PzgXMAImIYcBgwIDP/CNwI/AFYDbySmY9VuF5JUhOVDoloYVs2a18P7B8RTwOXAvOALRGxP6WzjkHAQcA+EfHR7Z4g4uKIqI+I+rVr13Zo8ZLU3VU6JBqAQ5q0B9DsJaPM3JCZF2TmUEprEv2AFcAHgRWZuTYzNwP3Af+1+RNk5pTMrMnMmn79+lVoGpLUPVU6JH4HHBERgyJiL0oLzw80HRAR+5X7AC4EnszMDZReZhoeEXtHRAAfABZXuF5JUhM9K3nwzNwSEZcAP6N0ddLUzFwYERPK/ZOAo4HbI+ItYBHwT+W+30bEvcBTwBZKL0NNqWS9kqS3i8zmSwS7r5qamqyvr692GZK0W4mIuZlZ01Kf77iWJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSoYqHRESMjojnImJpRFzZQv/+EXF/RDwTEXMi4tgmfftFxL0RsSQiFkfEKZWuV5L0VxUNiYjoAXwbOAMYDIyPiMHNhl0FPJ2ZxwEfB77ZpO+bwKOZeRQwBFhcyXolSW9X6TOJYcDSzFyemW8CdwNnNRszGHgcIDOXAAMj4sCIeBfwPuD75b43M/PlCtcrSWqi0iFxMLCySbuhvK2p+cA5ABExDDgMGAAcDqwF/j0i5kXErRGxT/MniIiLI6I+IurXrl1biTlIUrdV6ZCIFrZls/b1wP4R8TRwKTAP2AL0BE4AvpuZxwN/AbZb08jMKZlZk5k1/fr168jaJanb61nh4zcAhzRpDwBWNR2QmRuACwAiIoAV5cfeQENm/rY89F5aCAlJUuVU+kzid8ARETEoIvYCzgMeaDqgfAXTXuXmhcCTmbkhM/8ErIyII8t9HwAWVbheSVITFT2TyMwtEXEJ8DOgBzA1MxdGxIRy/yTgaOD2iHiLUgj8U5NDXArcUQ6R5ZTPOCRJnSMymy8R7L5qamqyvr6+2mVI0m4lIuZmZk1Lfb7jWpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJUyJCQJBUyJCRJhQwJSVIhQ0KSVMiQkCQVMiQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUqHIzGrX0GEiYi3wYrXr2AV9gZeqXUQnc87dg3PePRyWmf1a6tijQmJ3FRH1mVlT7To6k3PuHpzz7s+XmyRJhQwJSVIhQ6JrmFLtAqrAOXcPznk355qEJKmQZxKSpEKGhCSpkCHRSSLigIj4eUQ8X/66f8G40RHxXEQsjYgrW+i/LCIyIvpWvur2ae+cI+KGiFgSEc9ExP0RsV+nFd8GrfidRUTcUu5/JiJOaO2+XdWuzjkiDomIJyJicUQsjIjPdH71u6Y9v+dyf4+ImBcRD3Ze1R0gM310wgP4N+DK8vdXAhNbGNMDWAYcDuwFzAcGN+k/BPgZpTcM9q32nCo9Z+BDQM/y9xNb2r/aj539zspjzgQeAQIYDvy2tft2xUc759wfOKH8/b7A7/f0OTfp/9/AncCD1Z5PWx6eSXSes4Dbyt/fBvxDC2OGAUszc3lmvgncXd5vm68DlwO7y9UG7ZpzZj6WmVvK42YDAypb7i7Z2e+Mcvv2LJkN7BcR/Vu5b1e0y3POzNWZ+RRAZm4EFgMHd2bxu6g9v2ciYgAwBri1M4vuCIZE5zkwM1cDlL/+pxbGHAysbNJuKG8jIsYBf8zM+ZUutAO1a87NfIrSX2ldTWvqLxrT2rl3Ne2Zc6OIGAgcD/y240vscO2d8zco/YG3tUL1VUzPahewJ4mIGcB/bqHr6tYeooVtGRF7l4/xoV2trVIqNedmz3E1sAW4o23VdYqd1r+DMa3Ztytqz5xLnRF/C/wY+GxmbujA2ipll+ccEWOBNZk5NyJO6+jCKs2Q6ECZ+cGivoj487bT7fIp6JoWhjVQWnfYZgCwCvgvwCBgfkRs2/5URAzLzD912AR2QQXnvO0YnwDGAh/I8gu7XcwO69/JmL1asW9X1J45ExG9KAXEHZl5XwXr7EjtmfOHgXERcSbQG3hXRPwwMz9awXo7TrUXRbrLA7iBty/i/lsLY3oCyykFwrbFsWNaGPcCu8fCdbvmDIwGFgH9qj2XHcxxp78zSq9FN13QnNOW33dXe7RzzgHcDnyj2vPorDk3G3Mau9nCddUL6C4PoA/wOPB8+esB5e0HAQ83GXcmpSs+lgFXFxxrdwmJds0ZWErpNd6ny49J1Z5TwTy3qx+YAEwofx/At8v9zwI1bfl9d8XHrs4ZGEnpZZpnmvxez6z2fCr9e25yjN0uJLwthySpkFc3SZIKGRKSpEKGhCSpkCEhSSpkSEiSChkSkqRChoQkqZAhIVVYROwTEQ9FxPyIWBARtdWuSWot790kVd5oYFVmjgGIiL+rcj1Sq3kmIVXes8AHI2JiRPy3zHyl2gVJrWVISBWWmb8HTqQUFl+LiC9WuSSp1Xy5SaqwiDgIWJ+ZP4yIV4FPVrkkqdUMCany3gvcEBFbgc3A/6xyPVKreRdYSVIh1yQkSYUMCUlSIUNCklTIkJAkFTIkJEmFDAlJUiFDQpJU6P8DGXcX0cIeNIYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "s_vals = np.arange(0,50,2)\n", "\n", "def logL_sb(s_val, b_val) :\n", " return -2*sum( [ np.log(scipy.stats.poisson.pmf(n, s_val*s_yield + b_val*b_yield)) for n, s_yield, b_yield in zip(data, s_yields, b_yields) ] )\n", "\n", "def prof_logL_s(s_val) :\n", " # ==> Define the profile likelihood using the best-fit b as above\n", " pass\n", "\n", "def prof_t_s(s_val) :\n", " # ==> Define t(s) = lambda(s) - lambda(shat) using the best-fit s value shat as before\n", " pass\n", "\n", "def t_sb(s_val, b_val) :\n", " # ==> Define the non-profile likelihood t(s,b) = lambda(s,b) - lambda(shat,b) for a fixed b (i.e. no systematics on b)\n", " pass\n", "\n", "# Plot the results\n", "plt.plot(s_vals, [ prof_t_s(s_val) for s_val in s_vals ], label='uncertain b')\n", "plt.plot(s_vals, [ t_sb(s_val, b_val=428) for s_val in s_vals ], linestyle='--', label='fixed b')\n", "plt.xlabel('s')\n", "plt.ylabel('prof_t_s')\n", "plt.axhline(1, linestyle=':')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So it works as before (just ignore the warning!). Note that the central value is also shifted: the best-fit $b$ is now different from the one we had fixed before ($b=428$ instead of our assumed $b=500$), so the signal also gets a different value. \n", "\n", "This illustrates an important advantage of models with nuisance parameters: since the parameters adapt to the data, they can avoid biases in the results coming from faulty assumptions -- the results are therefore more reliable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Bayesian methods: simple case\n", "\n", "The methods we have seen so far are *frequentist*. This in effect means that there are 2 difference classes of parameters\n", "* Observables, which are randomly given by the experiment.\n", "* Parameters, which always have a unique, fixed value -- which is known, or more often unknown.\n", "\n", "There is another class of methods, *Bayesian* methods, which have a different take on measurement. In this approach, everything is random:\n", "* Observables are random as before\n", "* Parameters are also considered random, with the associated PDF used to quantify the uncertainty on their values.\n", "\n", "So for instance, for a Gaussian measurement:\n", "* The frequentist approach is to say $s = 2 \\pm 3$, meaning that the interval $[-1, 5]$ has a 68% chance of containing the true value of $s$\n", "* The Bayesian approach is to say that $s \\sim G(2, 3)$, where the Gaussian encodes the same information.\n", "\n", "One thing to note is that in the frequentist setting, \"$[-1, 5]$ has a 68% chance of containing the true value of $s$\" is a statement on the interval itself, which is defined for a given experiment: another experiment would find another interval, and 68% of these intervals are guaranteed to contain the true value. For the Bayesian approach, the value of $s$ is itself random, and the experiment is just a way to gather more information on its distribution.\n", "\n", "The fundamental tool of Bayesian inference is Bayes' theorem. It can be written as\n", "$$\n", "P(s) \\propto L(s; n) \\pi(s)\n", "$$\n", "where \n", "* $\\pi(s)$ is the *prior* PDF -- a PDF giving the knowledge of $s$ that we had before the experiment\n", "* $L(s;n)$ is the usual likelihood\n", "* $P(s)$ is the *posterior* PDF on $s$, i.e. the distribution that encodes our knowledge of $s$ including the experimental information.\n", "\n", "A limit on $s$ can then be set by integrating the posterior $P(s)$ and finding $s_{up}$ so that $\\int^{s_{up}} P(s) ds = 95\\%$.\n", "\n", "Consider the example we had before:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa2ElEQVR4nO3deZgU9Z3H8feHQwkKHlwBAQcVRDRilGCiRDGuUVeDkqDxYCUoi0ncFeOaSHYTNWvc4GMOXK+Ix4rRqGyiC4pHECTGW1BBDbq4iIgSOTxQvGD47h9dgwPMUXNUd8/U5/U8/XRVdR0feoZv1/z6V79SRGBmZvnRptQBzMysuFz4zcxyxoXfzCxnXPjNzHLGhd/MLGfalTpAGl27do2KiopSxzAza1Hmz5+/OiK6bbm8RRT+iooK5s2bV+oYZmYtiqTXalruph4zs5xx4TczyxkXfjOznGkRbfxmZjVZv349y5cv5+OPPy51lJLq0KEDvXv3pn379qnWd+E3sxZr+fLldOrUiYqKCiSVOk5JRARr1qxh+fLl9OvXL9U2buoxsxbr448/pkuXLrkt+gCS6NKlS4P+6nHhN7MWLc9Fv0pD3wMXfjOznHEbv5m1Hr9v5rP/U7K5X8n06dP56U9/Sps2bWjXrh2TJ09m2LBhANx///1MmDCByspKxo0bx8SJE5v9+K2/8Df3L4KZlY8B98Gaddntf002IwYcvt9OjJhzA5JY+OJiThw7mpee+AOVlZWc9b1xzPrDlfTu1YMvHTGGESNGMGjQoGY9vpt6zMyaYN26jzjmpHMYfOgp7DPs29xx15/q3Wb77Ttuapdf9+FHm6afeuZF9ujXh90qerPNNu05aeQRTJ8+vdkzt/4zfjOzDN0/53F6fb4rM2+fDMB7az/gB//2ax56dOu/Fk4a+XUmTvgOAHfNfIgfX3wVK1e/w8zbfgPAGytW0adXj03r9+7VgydffKPZM7vwm5k1wRf22p3zLryc8392Bcd+fRhf/coX+c0l59a73chjDmPkMYfx8GPP8NNf/JYH77yamu6BnkWvJRd+M7MmGLDHrsyffTP3znqUH//8Kr4+/EDeeff9es/4qxxy0P7839I3WL3mXXr36s7rb7616bXlb75Fr169mz2zC7+ZWRO8uWIVO+/UmdEn/j3bb9+Rm267h//53S/r3OaVJa+ze7/eSOKZBS/x6afr6bLzDuy4wyAWL1nGq6+9wS49u3P7XbP4/R13NXtmF34zaz2OfLroh3x+0Sv88KL/pE0b0b5dO665rP7ul3+8Zw433zGT9u3b8bkOHbjj+v9AEu3atePKST/iyBPOpnJjJaefMoK999672TOrpjalcjNkyJBo9I1Y3J3TrNVaNOA+9urXtdQxstVlSKrVFi1axF577bXZMknzI2KrHbg7p5lZzmTa1CNpKfA+UAlsiIghknYG7gAqgKXAiRHxTpY5zMzsM8U44z8sIvar9ufGRGB2RPQHZifzZmZWJKVo6jkOmJpMTwWOL0EGM7PcyrrwB/AnSfMljU+W9YiIFQDJc/eMM5iZWTVZd+c8OCLelNQdmCXppbQbJh8U4wH69u2bVT4zs9zJtPBHxJvJ80pJdwFDgbck9YyIFZJ6Aitr2XYKMAUK3TmzzGlmrYO6puv6mFasrrsb+dJlb3LsKT/ghUfuaNbjzn1kPr+86hbuScbwaW6ZNfVI2k5Sp6pp4OvAC8AMYEyy2hig+YeeMzOzWmXZxt8DeETSAuApYGZE3A9MAo6QtBg4Ipk3M2uRNmyoZMxZF7HvISczauz5fPjh1ve+fe75l/nykWPZ95CTGXnaD3nn3bUADB9xJuf/7AqGHjGGAUO/xV8ef3az7TZu3Ej//v1ZtWrVpvk99tiD1atXNylzZoU/IpZExODksXdEXJIsXxMRh0dE/+T57awymJll7eVXXmP8aSNZ+PBtdO60HVff+N9brXPaWRdx6QX/xMKHb+MLe+3Ozy67btNrGzZs4KlZU5l8ybmbLQdo06YNo0eP5tZbbwXgwQcfZPDgwXTt2rSrlX3lrplZE/TZpQcHHzgYgNGjjuaRJxds9vp7az/g3ffe59CDDwBgzEnH8vBjn53Zf/PYrwFwwOCBLH19xVb7P/3007n55psBuPHGGxk7dmyTM7vwm5k1wZbj5Td0/Pxtt2kPQNu2bdmwoXKr1/v06UOPHj2YM2cOTz75JEcffXTjwyZc+M3MmmDZ8r/x+NMLAbjtzgcYlpz9V9mh8/bstGPnTe33v5t2L4cetH+DjjFu3DhGjx7NiSeeSNu2bZuc2cMymzXS8J8Xnuf+pLQ57DP1db/Mwl4D+jH19pmc+S+/oP9uffje2FFbrTP1ygv57nmT+PCjj9lt1134rysuaNAxRowYwdixY5ulmQdc+M3MGq2iby/++ti0etfb7wt78sQD/7XV8rkzrt003bXLjix9dgYAw4cdwPBhB2x6bcGCBQwePJiBAwc2Q2oXfjOzsjZp0iSuueaaTT17moPb+M3MytjEiRN57bXXGDZsWLPt04XfzFqu2EgLuIlg5hp6J0UXfjNrsTp88gpr3t+Q6+IfEaxZs4YOHTqk3sZt/GbWYvV+4yKWcxGrtt0D1ErPY1cuqneVDh060Lt379S7dOE3sxarfeU79Fs2odQxsnVK8/8500o/Is3MrDYu/GZmOePCb2aWMy78ZmY548Jv1gi3PgpPLIY/L4KKCYV5s5bChd+sgW59FMZfD59sKMy/trow7+JvLYULv1kD/ds0+PDTzZd9+GlhuVlL4MJv1kDLarndaW3LzcqNC79ZA/Wt5XantS03Kzcu/GYNdMmJ0HGbzZd13Kaw3KwlcOE3a6BTD4Yp42DbZMCTXbsW5k89uLS5zNLyWD1mjXDqwXDdQ4Vp33rRWhqf8ZuZ5YwLv5lZzrjwm5nljAu/mVnOuPCbmeVMq+/Vo1NzfDNOy9hwAHTq3JKmsNYtTmn+ffqM38wsZzIv/JLaSnpW0j3J/M6SZklanDzvlHUGMzP7TDHO+CcA1W8TPxGYHRH9gdnJvJmZFUmmhV9Sb+AY4Ppqi48DpibTU4Hjs8xgZmabq/XLXUl3A7V+MxoRI1LsfzLwI6BTtWU9ImJFso8VkrrXcvzxwHiAvn37pjiUmZmlUdcZ/y+BXwGvAh8B1yWPD4AX6tuxpGOBlRExvzHBImJKRAyJiCHdunVrzC7MzKwGtZ7xR8SfASRdHBGHVHvpbkkPp9j3wcAISX8PdAA6S7oFeEtSz+Rsvyewsgn5zcysgdK08XeTtFvVjKR+QL2n4BHx44joHREVwEnAnIgYDcwAxiSrjQGmNzi1mZk1WpoLuH4AzJW0JJmvAM5swjEnAdMknQEsA05owr7MzKyB6i38EXG/pP7AwGTRSxHxSUMOEhFzgbnJ9Brg8IbFNDOz5lJXr55v1vLS7pKIiDszymRmZhmq64z/G3W8FoALv5lZC1RXr56xxQxiZmbFUW+vHkk9JN0g6b5kflDyxayZmbVAabpz3gQ8APRK5v8XOCejPGZmlrE0hb9rREwDNgJExAagMtNUZmaWmTSFf52kLiTj9kj6MvBepqnMzCwzaS7gOpfC1ba7S3qUwlW7ozJNZWZmmUlzAdczkg4F9gQEvBwR6zNPZmZmmajrAq6vRcScGi7kGuALuMzMWq66zvgPBeZQ84VcvoDLzKyFqusCrguTyXER4V48ZmatRJpePa9IukzSoMzTmJlZ5tIU/n0pXLR1vaQnJI2X1DnjXGZmlpF6C39EvB8R10XEQRTun3shsELSVEl7ZJ7QrGzNTR5mLUuasXraShoh6S7gcgr34d0NuBu4N+N8ZmbWzNJcwLUYeAi4LCIeq7b8D5IOqWUbMzMrU2kK/74R8UFNL0TE2c2cx8zMMpamjb/Gom9mZi1Tml49ZmbWirjwm5nlTF1j9Zxb14YR8evmj2NmZlmr68vdTsnznsCXKAzNDIWxex7OMpSZmWWnrrF6fgYg6U/A/hHxfjJ/EfDfRUlnZmbNLk0bf1/g02rznwIVmaQxM7PMpenH/zvgqeTK3QBGAjdnmsrMzDKT5g5cl0i6D/hqsmhsRDybbSwzM8tK2u6cHYG1EXE5sFxSvwwzmZlZhtIM0nYhcD7w42RRe+CWLEOZmVl20pzxjwRGAOsAIuJNPuvqaWZmLUyawv9pRASFL3aRtF2aHUvqIOkpSQskvSipqnvozpJmSVqcPO/U+PhmZtZQaQr/NEnXAjtK+kfgQeD6FNt9AnwtIgYD+wFHSfoyMBGYHRH9gdnJvJmZFUmaXj2/lHQEsJbCVbwXRMSsFNsFUDWyZ/vkEcBxwPBk+VQKtzA6v6HBzcyscdJ8uXtpRMyKiB9GxHkRMUvSpWl2nty96zlgJTArIp4EekTECoDkuXst246XNE/SvFWrVqX+B5mZWd3SNPUcUcOyo9PsPCIqI2I/oDcwVNI+aYNFxJSIGBIRQ7p165Z2MzMzq0ddo3N+D/g+sLukhdVe6gQ8VvNWNYuIdyXNBY4C3pLUMyJWSOpJ4a8BMzMrkrrO+H9PYSTO6clz1eOAiDi1vh1L6iZpx2T6c8DfAS9RGOVzTLLamGT/ZmZWJHWNzvke8J6ky4G3q43O2UnSgUl7fV16AlMltaXwATMtIu6R9DiFnkJnAMuAE5rlX2JmZqmkGaTtGmD/avPrali2lYhYCHyxhuVrgMMbkNHMzJpRmi93lXTNBCAiNpLuA8PMzMpQmsK/RNLZktonjwnAkqyDmZlZNtIU/u8CBwFvAMuBA4HxWYYyM7PspLlydyVwUhGymJlZEaS5cneApNmSXkjm95X0k+yjmZlZFtI09VxHYSz+9bCpt47/AjAza6HSFP6OEfHUFss2ZBHGzMyyl6bwr5a0O5+Nxz8KWJFpKjMzy0ya/vhnAVOAgZLeAF4F6h2ywczMylOaXj1LgL9L7rzVpmroBjMza5nS9OrpIuk/gb8AcyVdLqlL9tHMzCwLadr4bwdWAd8CRiXTd2QZyszMspOmjX/niLi42vzPJR2fUR4zM8tYmjP+hySdJKlN8jgRmJl1MDMzy0aawn8mhZuyfJI8bgfOlfS+pLVZhjMzs+aXpldPp2IEMTOz4kjTq+eMLebbSrowu0hmZpalNE09h0u6V1JPSV8AnqBww3UzM2uB0jT1nCLp28DzwIfAyRHxaObJzMwsE2maevoDE4A/AkuBf5DUMeNcZmaWkTRNPXcDF0TEmcChwGLg6UxTmZlZZtJcwDU0ItYCJDdd/5WkGdnGMjOzrKQ54/+cpBsk3Q8gaRBwSLaxzMwsK2kK/03AA0DPZP5/gXMyymNmZhlLU/i7RsQ0YCNARGwAKjNNZWZmmUlT+NclwzBX3YHry8B7maYyM7PMpPly91xgBrC7pEeBbhSGZzYzsxYozQVcz0g6FNgTEPByRKzPPJmZmWUizRl/Vbv+ixlnMTOzIkjTxm9mZq1IZoVfUh9JD0laJOlFSROS5TtLmiVpcfK8U1YZzMxsa7U29Ujav64NI+KZeva9AfiX5DuCTsB8SbOA7wCzI2KSpInAROD8hsU2M7PGqquN/1d1vBbA1+racUSsAFYk0+9LWgTsAhwHDE9WmwrMxYXfzKxoai38EXFYcx1EUgXwReBJoEfyoUBErJDUvZZtxgPjAfr27dtcUczMci9Vrx5J+wCDgA5VyyLi5pTbbk9hSOdzImKtpFTBImIKMAVgyJAhkWojMzOrV72FP7nN4nAKhf9e4GjgEaDewi+pPYWif2tE3JksfktSz+RsvyewspHZzcysEdL06hkFHA78LSLGAoOBbevbSIVT+xuARRHx62ovzQDGJNNjgOkNSmxmZk2Spqnno4jYKGmDpM4UztB3S7HdwcA/AM9Lei5Z9q/AJGBachP3ZcAJDY9tZmaNlabwz5O0I3AdMB/4AHiqvo0i4hEKQzzU5PC0Ac3MrHmlGavn+8nkb5ObsXSOiIXZxjIzs6ykudn67KrpiFgaEQurLzMzs5alrit3OwAdga7JsApVzTadgV5FyGZmZhmoq6nnTAq3WOwFVB+eYS1wVYaZzMwsQ3VduXs5cLmkf46IK4qYyczMMpSmV8+1ks4GDknm5wLX+mYsZmYtU5rCfzXQPnmGQt/8a4BxWYUyM7Ps1PXlbrvkzltfiojB1V6aI2lB9tHMzCwLdXXnrLpIq1LS7lULJe0GVGaayszMMlNXU09V983zgIckLUnmK4CxWYYyM7Ps1FX4u0k6N5m+FmgLrKMwNPMXgYcyzmZmZhmoq/C3BbZn8/F2tk+eO2WWyMzMMlVX4V8REf9etCRmZlYUdX25m+5WWWZm1qLUVfg9dLKZWStUa+GPiLeLGcTMzIojza0XzcysFXHhNzPLGRd+M7OcceE3M8sZF34zs5xx4TczyxkXfjOznHHhNzPLGRd+M7OcceE3M8sZF34zs5xx4TczyxkXfjOznHHhNzPLmcwKv6QbJa2U9EK1ZTtLmiVpcfK8U1bHNzOzmmV5xn8TcNQWyyYCsyOiPzA7mTczsyLKrPBHxMPAljdzOQ6YmkxPBY7P6vhmZlazYrfx94iIFQDJc/faVpQ0XtI8SfNWrVpVtIBmZq1d2X65GxFTImJIRAzp1q1bqeOYmbUaxS78b0nqCZA8ryzy8c3Mcq/YhX8GMCaZHgNML/LxzcxyL8vunLcBjwN7Slou6QxgEnCEpMXAEcm8mZkVUbusdhwRJ9fy0uFZHdPMzOpXtl/umplZNlz4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McqYkhV/SUZJelvSKpImlyGBmlldFL/yS2gJXAUcDg4CTJQ0qdg4zs7wqxRn/UOCViFgSEZ8CtwPHlSCHmVkutSvBMXcBXq82vxw4cMuVJI0HxiezH0h6uZHH6wqsbuS2WXKuhnGuhnGuhinXXEhNyrZrTQtLUfhVw7LYakHEFGBKkw8mzYuIIU3dT3NzroZxroZxroYp11yQTbZSNPUsB/pUm+8NvFmCHGZmuVSKwv800F9SP0nbACcBM0qQw8wsl4re1BMRGyT9E/AA0Ba4MSJezPCQTW4uyohzNYxzNYxzNUy55oIMsiliq+Z1MzNrxXzlrplZzrjwm5nlTKsq/JLaSnpW0j01vDZQ0uOSPpF0XhnlOlXSwuTxmKTBZZLruCTTc5LmSRpWDrmqrfMlSZWSRpVDLknDJb2XvF/PSbqgHHJVy/acpBcl/bkcckn6YbX36oXkZ7lzGeTaQdLdkhYk79fYYmRKkWsnSXcl/yefkrRPU45Vin78WZoALAI61/Da28DZwPHFDJSoK9erwKER8Y6koyl8kbPVBW0lyDUbmBERIWlfYBowsAxyVQ37cSmFDgLFVGcu4C8RcWwR81SpNZekHYGrgaMiYpmk7uWQKyIuAy5LMn4D+EFEvF3qXMBZwF8j4huSugEvS7o1GWWglLn+FXguIkZKGkhh2JvDG3ugVnPGL6k3cAxwfU2vR8TKiHgaWF9muR6LiHeS2ScoXNdQDrk+iM+++d+OGi6yK0WuxD8DfwRWFiMTpM5VdClynQLcGRHLoPD/oExyVXcycFu2iQpS5AqgkyQB21M4YdxQBrkGUTgZIyJeAiok9Wjs8VpN4QcmAz8CNpY4x5Ymkz7XGcB9mab5zGTqySVppKSXgJnA6eWQS9IuwEjgt0XKU2Uy9f8cv5I0Edwnae/ixKo31wBgJ0lzJc2XdFqZ5AJAUkfgKAof5MUwmbpzXQnsReGi0ueBCRFRjJpSX64FwDcBJA2lMBRDo08SW0Xhl3QssDIi5pc6S3UNySXpMAqF//xyyRURd0XEQArNYxeXSa7JwPkRUZl1niopcz0D7BoRg4ErgP8pk1ztgAMonE0eCfxU0oAyyFXlG8CjxWjmSZnrSOA5oBewH3ClpNqa9oqZaxKFD/DnKPzF+yxN+UskIlr8A/gFhaEglgJ/Az4Ebqll3YuA88opF7Av8H/AgHLKtcU2rwJdS50rybE0eXxAobnn+FLnqmGbpWXyfk0ELqo2fwNwQqlzVVv3LuCULPM08P2aCXy12vwcYGipc22xvpJ1Ozf6mMV4w4v5AIYD99TxetEKf5pcQF/gFeCgcnq/gD347AK//YE3qubL4eeYrHMTMKpM3q/PV3u/hgLLyuH9otBsMZvCmX9H4AVgn1LnSl7bgUIb+nbF/BnW835dU/VBCfRIfu8z/QBPmWtHYJtk+h+Bm5tynNbWq2czkr4LEBG/lfR5YB6Fb8w3SjoHGBQRa0uZC7gA6AJcXfg+iQ1RolECt8j1LeA0SeuBj4BvR/JbV+JcZWOLXKOA70naQOH9Oqkc3q+IWCTpfmAhhfbj6yPihVLnShaNBP4UEetKkafKFrkuBm6S9DyFM+vzI6IkwzVvkWsv4GZJlcBfKTQLN37fJfrdNDOzEmkVX+6amVl6LvxmZjnjwm9mljMu/GZmOePCb2aWMy781qpIqpBUY3dFSddLGlTsTDXkqDWjWTG06n78ZtVFxLhSZ2gOktpFROYDh1nr5TN+a43aSZqajF3+h2QgMJKByoYk0x9IuiQZVO2JmkY6lHSRpBuT7ZZIOjtZvtkZu6TzJF1U7Ri/kfSwpEUq3DfgTkmLJf08RcYDJP05GVDtAUk9q+33P1QYT39CRu+b5YQLv7VGewJTImJfYC3w/RrW2Q54IgqDqj1M4TL4mgykMHDXUOBCSe1THP/TiDiEwgii0ymM8b4P8B1JXWrLmOz7CgrDUBwA3AhcUm2/O0bEoRHxqxQZzGrlwm+t0esR8WgyfQtQ093DPgWq7nQ0H6ioZV8zI+KT5LL9lRTGb6nPjOT5eeDFiFgREZ8AS4A+dWTck8IHxKxkFMafsPnQu3ekOLZZvdzGb63RluOQ1DQuyfpqY+lUUvv/hU+qTVett4HNT5o61LLNxi2231jtODVlFIUPiq/UkqWkY9pY6+EzfmuN+kqqKp4nA4808/7fArpL6iJpW6Axt1usKePLQLeq5ZLaF/GGLpYjLvzWGi0CxkhaCOxMYajdZhMR64F/B56k0Fz0UiN2s1XGKNzXdRRwqaQFFG4IclCzhDarxqNzmpnljM/4zcxyxoXfzCxnXPjNzHLGhd/MLGdc+M3McsaF38wsZ1z4zcxy5v8BuBtRklY7Ry0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.bar(x[4], sb4, color='orange', label='s=30')\n", "plt.bar(x[4], b4, color='b', yerr=np.sqrt(b_only[4]), label='b only')\n", "plt.scatter(x[4], n4, zorder=10, color='k')\n", "plt.xlabel('bin number')\n", "plt.ylabel('Total expected yield')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and assume we have a flat prior on $s$, $\\pi(s) = k$. Flat priors are often used in Bayesian applications, since they provide a relatively unbiased starting point.\n", "(Note that technically the flat distribution is not well defined if the range of $s$ is infinite, but this won't crea te problems here).\n", "\n", "Now we can derive the upper limit:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAObElEQVR4nO3dX4ild33H8fenuw3UPzWhGUV3I92W1bgtpugYRfonVlqz8WIRvEi0DQ3CsmDE3pSElv4Bb+pFQcTosoQleOPeGOxaYtPSoimkqZmFGHeVyLjSZFwhGxULEZpu/PbinLbTyWzOszvPmbM73/cLBuZ5zm/P9zfZz3z2mfNnkqpCkrTz/dyiNyBJ2h4WviQ1YeFLUhMWviQ1YeFLUhMWviQ1MbPwkxxP8myS0xe5PUk+nWQ1yZNJ3jb+NqXxmW11M+QK/wHg1pe5/SCwf/pxGPjc1rclbYsHMNtqZGbhV9UjwI9eZskh4PM18RhwbZLXj7VBaV7MtrrZPcJ97AGeWXe8Nj33g40LkxxmcqXEK1/5yrffeOONI4yXXurUqVPPVdXSFu/GbOuKs5Vsj1H42eTcpr+voaqOAccAlpeXa2VlZYTx0ksl+fcx7maTc2ZbC7WVbI/xKp014IZ1x3uBcyPcr7RoZls7yhiFfxK4c/qKhncBP6mql/zIK12FzLZ2lJkP6ST5AnALcH2SNeAvgZ8HqKqjwEPAbcAq8FPgrnltVhqT2VY3Mwu/qu6YcXsBHx1tR9I2MdvqxnfaSlITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITgwo/ya1JnkqymuTeTW5/TZIvJ/lGkjNJ7hp/q9K4zLW6mVn4SXYB9wEHgQPAHUkObFj2UeBbVXUTcAvwN0muGXmv0mjMtToacoV/M7BaVWer6gXgBHBow5oCXp0kwKuAHwEXRt2pNC5zrXaGFP4e4Jl1x2vTc+t9BngLcA74JvDxqvrZxjtKcjjJSpKV8+fPX+aWpVGMlmsw27o6DCn8bHKuNhy/D3gCeAPwG8BnkvziS/5Q1bGqWq6q5aWlpUvcqjSq0XINZltXhyGFvwbcsO54L5MrnvXuAh6siVXge8CN42xRmgtzrXaGFP7jwP4k+6ZPWN0OnNyw5mngvQBJXge8GTg75kalkZlrtbN71oKqupDkbuBhYBdwvKrOJDkyvf0o8AnggSTfZPKj8j1V9dwc9y1tiblWRzMLH6CqHgIe2nDu6LrPzwG/P+7WpPky1+rGd9pKUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1YeFLUhMWviQ1Majwk9ya5Kkkq0nuvciaW5I8keRMkq+Nu01pfOZa3eyetSDJLuA+4PeANeDxJCer6lvr1lwLfBa4taqeTvLaOe1XGoW5VkdDrvBvBlar6mxVvQCcAA5tWPMh4MGqehqgqp4dd5vS6My12hlS+HuAZ9Ydr03Prfcm4LokX01yKsmdm91RksNJVpKsnD9//vJ2LI1jtFyD2dbVYUjhZ5NzteF4N/B24P3A+4A/T/Kml/yhqmNVtVxVy0tLS5e8WWlEo+UazLauDjMfw2dy5XPDuuO9wLlN1jxXVc8Dzyd5BLgJ+M4ou5TGZ67VzpAr/MeB/Un2JbkGuB04uWHN3wK/lWR3klcA7wS+Pe5WpVGZa7Uz8wq/qi4kuRt4GNgFHK+qM0mOTG8/WlXfTvL3wJPAz4D7q+r0PDcubYW5Vkep2viw5fZYXl6ulZWVhczWzpfkVFUtL2K22dY8bSXbvtNWkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpqw8CWpCQtfkpoYVPhJbk3yVJLVJPe+zLp3JHkxyQfH26I0H+Za3cws/CS7gPuAg8AB4I4kBy6y7pPAw2NvUhqbuVZHQ67wbwZWq+psVb0AnAAObbLuY8AXgWdH3J80L+Za7Qwp/D3AM+uO16bn/leSPcAHgKMvd0dJDidZSbJy/vz5S92rNKbRcj1da7Z1xRtS+NnkXG04/hRwT1W9+HJ3VFXHqmq5qpaXlpYGblGai9FyDWZbV4fdA9asATesO94LnNuwZhk4kQTgeuC2JBeq6ktjbFKaA3OtdoYU/uPA/iT7gO8DtwMfWr+gqvb9z+dJHgD+zm8KXeHMtdqZWfhVdSHJ3UxepbALOF5VZ5Icmd4+8/FN6UpjrtXRkCt8quoh4KEN5zb9hqiqP9r6tqT5M9fqxnfaSlITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNWHhS1ITFr4kNTGo8JPcmuSpJKtJ7t3k9g8neXL68WiSm8bfqjQuc61uZhZ+kl3AfcBB4ABwR5IDG5Z9D/idqnor8Ang2NgblcZkrtXRkCv8m4HVqjpbVS8AJ4BD6xdU1aNV9ePp4WPA3nG3KY3OXKudIYW/B3hm3fHa9NzFfAT4ymY3JDmcZCXJyvnz54fvUhrfaLkGs62rw5DCzybnatOFyXuYfGPcs9ntVXWsqparanlpaWn4LqXxjZZrMNu6OuwesGYNuGHd8V7g3MZFSd4K3A8crKofjrM9aW7MtdoZcoX/OLA/yb4k1wC3AyfXL0jyRuBB4A+r6jvjb1ManblWOzOv8KvqQpK7gYeBXcDxqjqT5Mj09qPAXwC/BHw2CcCFqlqe37alrTHX6ihVmz5sOXfLy8u1srKykNna+ZKcWlQ5m23N01ay7TttJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJakJC1+SmrDwJamJQYWf5NYkTyVZTXLvJrcnyaentz+Z5G3jb1Ual7lWNzMLP8ku4D7gIHAAuCPJgQ3LDgL7px+Hgc+NvE9pVOZaHQ25wr8ZWK2qs1X1AnACOLRhzSHg8zXxGHBtktePvFdpTOZa7ewesGYP8My64zXgnQPW7AF+sH5RksNMrpQA/jPJ6Uva7XiuB55rNHeRsxc1980zbh8t13DFZNt89Zg9K9sXNaTws8m5uow1VNUx4BhAkpWqWh4wf3SLmu3XvL1zZy3Z5Nxl5RqujGybrx6zB2T7ooY8pLMG3LDueC9w7jLWSFcSc612hhT+48D+JPuSXAPcDpzcsOYkcOf0VQ3vAn5SVS/5sVe6gphrtTPzIZ2qupDkbuBhYBdwvKrOJDkyvf0o8BBwG7AK/BS4a8DsY5e9661b1Gy/5itk7hxzPXP2HJmvHrMve26qNn1IUpK0w/hOW0lqwsKXpCbmXviLevv6gLkfns57MsmjSW4aY+6Q2evWvSPJi0k+uF1zk9yS5IkkZ5J8bYy5Q2YneU2SLyf5xnT20MfDZ809nuTZi73ufYH5mtuvZVhUtheV66Gz55HtHZfrqprbB5Mnw74L/ApwDfAN4MCGNbcBX2Hymud3Af+2TXPfDVw3/fzgGHOHzl637p+ZPDH4wW36mq8FvgW8cXr82m38e/5T4JPTz5eAHwHXjDD7t4G3Aacvcvui8jX63EVme1G5XmS2d2Ku532Fv6i3r8+cW1WPVtWPp4ePMXmN9RiGfM0AHwO+CDy7jXM/BDxYVU8DVNV2zi7g1UkCvIrJN8aFrQ6uqkem93UxC8nXnOYOmj2nbC8q10NnzyPbOy7X8y78i701/VLXzGPueh9h8q/lGGbOTrIH+ABwdKSZg+YCbwKuS/LVJKeS3LmNsz8DvIXJG5e+CXy8qn420vyt7m0e9zmPuZdzv2Nle1G5HjSb+WR7x+V6yK9W2IpR374+8tzJwuQ9TL4pfnOLMy9l9qeAe6rqxcmFwbbN3Q28HXgv8AvAvyZ5rKq+sw2z3wc8Afwu8KvAPyb5l6r6jy3OHmNv87jPecy9pPsdOduLyvXQ2fPI9o7L9bwLf1FvXx90n0neCtwPHKyqH25x5qXMXgZOTL8prgduS3Khqr4057lrwHNV9TzwfJJHgJuArRb+kNl3AX9dkwcgV5N8D7gR+PoWZ4+xt3nc57x+LcOisr2oXA+dPY9s77xcb/XJhRlPPOwGzgL7+L8nPX5tw5r38/+ffPj6Ns19I5N3UL57u7/mDesfYJwnbYd8zW8B/mm69hXAaeDXt2n254C/mn7+OuD7wPUj/Tf/ZS7+5Nai8jX63EVme1G5XmS2d2KuRwnDjE3fxuRf2e8CfzY9dwQ4Mv08TP5HFN9l8hjY8jbNvR/4MZMfx54AVrbra96wdsxvjJlzgT9h8mqG08Afb+Pf8xuAf5j+HZ8G/mCkuV9g8uuK/4vJVc9HrpB8zWXuIrO9qFwvMts7Ldf+agVJasJ32kpSExa+JDVh4UtSExa+JDVh4UtSExa+JDVh4UtSE/8NbVa/mhVCdpUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def p_s(s_val) :\n", " return scipy.stats.poisson.pmf(n4, s_val + b4)\n", "\n", "import scipy.integrate\n", "\n", "def integral_s(s_val) :\n", " return scipy.integrate.quad(p_s, -20, s_val)[0]\n", "\n", "s_vals = np.linspace(-20,30,50)\n", "plt.subplot(121)\n", "# ==> Plot p(s) as a function of s\n", "plt.subplot(122)\n", "# ==> Plot integral_s(s) as as function of s\n", "# ==> Use scipy.optimize.root_scalar to compute the Bayesian 95% CL upper limit by finding the point where the integral reaches 95%." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Bayesian methods: nuisance parameters\n", "\n", "Models usually also include nuisance parameters. For instance we considered a one-bin model where the background was only defined up to an uncertainty, so that the full likelihood is\n", "$$\n", "L(s, b) = P(n_4; s_4 + b_4) G(b^{\\text{aux}}, b_4, \\sigma_b)\n", "$$\n", "\n", "The way to treat this in a Bayesian way is to integrate over the nuisance parameters: formally\n", "$$\n", "P(s) \\propto \\int L(s, b; n) \\pi(s) G(b) db\n", "$$\n", "So that the constraint becomes the prior PDF of $b$. Using an uncertainty of $5$ on the background, and a central value equal to $b_4$, we have" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def p_sb(s_val, b_val) :\n", " return scipy.stats.poisson.pmf(n4, s_val + b_val)*scipy.stats.norm.pdf(b_val, loc=b4, scale=5)\n", "\n", "def integral_intb(s_val) :\n", " # ==> Marginalize p_sb over b (integrate over b at fixed s)\n", " pass\n", "\n", "def integral_sb(s_val) :\n", " # ==> Return the integral of integral_intb from -infinity (actually -25 should do) to s_val \n", " pass\n", "\n", "s_vals = np.linspace(-20,80,51)\n", "# ==> Plot p_s(s) and integral_intb(s) on the same plot.\n", "# You should see that the latter is a smeared version of the former, reflecting that it is\n", "# a convolution of p_s with an uncertainty on b\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# ==> Plot integral_sb as a function of s (takes a while, since it is a double integral on s and b!)\n", "# ==> Compute the Bayesian limit including b systematics, by finding the point where integral_sb = 0.95" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should find that the limit in this case is slightly larger than in the no-systematics case, as one can expect." ] }, { "cell_type": "code", "execution_count": null, "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }