{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.special import gammaincc\n",
"from scipy.stats import lognorm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Limits with uncertain background an efficiencies
\n",
"The model is $N=\\alpha S + \\beta B$
\n",
"$\\alpha$ and $\\beta$ have lognormal pdfs of mean 1 and unc. $\\delta\\alpha$ and $\\delta\\beta$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Here are counts and backgrounds and uncertainties\n",
"N = [0, 2, 0, 2, 0, 2]\n",
"B = [0.5, 0.5, 2.0, 2.0, 3.5, 3.5]\n",
"dalpha = 0.2\n",
"dbeta = 0.4\n",
"Ntoys = 10000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To calculate the frequentist limit use result\n",
"$\\sum_{n=0}^N \\frac{\\mu^n e^{-\\mu}}{N!} = \\frac{1}{N!}\\Gamma(N+1, \\mu)$
\n",
"where $\\Gamma(a,b) = \\int_b^{\\infty} t^{a-1}e^{-t} dt$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbeUlEQVR4nO3df5AV5Z3v8fc3iEAUGRNmCZnBjFvqqOAqgoqB2ky0jEZNsCJa3jVZUDdU/JlYeyuypnLDpawUm0qtK67RpaIXTMyql929ElDvTZRTBguJg9EERSOxRhkKo6LMMgKJY33vH+eZsefMmZk+Z/qcPqfP51U1Nf3r9DxPT/OZ5umnnzZ3R0REsuVjaRdARESSp3AXEckghbuISAYp3EVEMkjhLiKSQYelXQCApqYmP+6449IuRmref/99jjjiiLSLkapGPwaNXn/QMSin/tu2bXvH3ZuLrauJcJ82bRqdnZ1pFyM1uVyOjo6OtIuRqkY/Bo1ef9AxKKf+Zvb6cOvULCMikkGxwt3MmsxsnZm9bGY7zOxsM/uEmf3CzF4N348O25qZrTKznWb2WzM7vbJVEBGRQnGv3O8AHnf3E4FTgR3AMuAJdz8eeCLMA3wROD58LQXuTrTEIiIyqlHb3M1sCvDXwBIAd/8z8GczWwh0hM3WAjngFmAhcL/nxzV4Jlz1T3f3PYmXXkQaxgcffEB3dzeHDh1KuygVMWXKFHbs2FF03cSJE2ltbWX8+PGx9xfnhuqxwNvA/zKzU4FtwDeBaZHAfhOYFqZbgF2Rz3eHZYPC3cyWkr+yp7m5mVwuF7vQWdPb29vQ9Qcdg0avP4x+DI488kimTZtGS0sLZla9glXJhx9+yLhx44Ysd3d6enp44YUX6O3tjb2/OOF+GHA6cKO7bzWzO/ioCab/h7uZlTQCmbuvBlYDtLe3u+6Sd6RdjFQ1+jFo9PrD6Mdgx44dtLa2ZjLYAfbv38/kyZOLrps8eTK9vb3MnTs39v7itLl3A93uvjXMryMf9n80s+kA4ftbYf1uYEbk861hmYjImGQ12EdTTr1HDXd3fxPYZWbtYdG5wEvAemBxWLYYeCRMrwf+NvSamQf0qL1dRKS64j7EdCPwgJkdDrwGXEX+D8PDZnYN8Dpwedj2UeBCYCdwIGwrIpKo+SufZPe+g4ntr6VpEk8vO6esz7a1tdHZ2cnUqVPHtE2SYoW7uz8PFGvsObfItg5cP7ZiSb0r/Ic3ln84IsXs3neQrpUXJba/tmUbE9tXLdATqlIR/f/w+r+SvMISSdMll1zCnDlzmDlzJqtXrx60rqurixNPPJErr7ySk046iUWLFnHgwIGB9XfeeSenn346p5xyCi+//DIAv/71rzn77LNZsGABn/3sZ3nllVcSKafCXUSkBPfddx/btm2js7OTVatWsXfv3kHrX3nlFa677jp27NjBUUcdxY9+9KOBdVOnTuW5557j2muv5Yc//CEAJ554Ir/61a/YvHkzK1as4NZbb02knAp3EZESrFq1ilNPPZV58+axa9cuXn311UHrZ8yYwfz58wH46le/yubNmwfWfeUrXwFgzpw5dHV1AdDT08Nll13GWWedxc0338yLL76YSDkV7iIiMeVyOX75y1+yZcsWXnjhBWbPnj3kidnCbovR+QkTJgAwbtw4+vr6APjud7/L5z//ebZu3crPf/7zxJ7AVbiLiMTU09PD0Ucfzcc//nFefvllnnnmmSHbvPHGG2zZsgWAn/3sZyxYsGDUfba0tACwZs2axMpaE+O5i4iUqqVpUqI9XFqaJo26zQUXXMA999zDSSedRHt7O/PmzRuyTXt7O3fddRdXX301J598Mtdee+2I+/z2t7/N4sWLWbFiBV/60pfKLn8hhbuI1KU0utZOmDCBxx57bMjy/vbz3t5eDjvsMH76058Ouw3A3LlzB8bROfvss/n9738/MPzAbbfdlkhZFe5SddE+8Or/LlIZCnepuujDJ1l7cEQaW1tbG9u3b0+7GIDCXaok2j4ap21TRMZG4S5VoaYXkepSV0gRkQzSlbvUFN1sFUmGwl1qim62Smy3nwI9byS3vynHwM2/G3GTrq4uLr744tg3TdesWcMXvvAFPv3pTydRwpIo3EWkPvW8Act7ktvf8inJ7StYs2YNs2bNSiXc1eYuqervRbPk8fdpW7ZRPWmk5vX19Q0Z0nfbtm187nOfY86cOZx//vns2bOHdevW0dnZyZVXXslpp53GwYMHWbFiBWeccQazZs1i6dKl5F9/URkKd0nV08vOoWvlRay54Ai6Vl6kNnapeYVD+t51113ceOONrFu3jm3btnH11Vfzne98h0WLFjF37lweeOABnn/+eSZNmsQNN9zAs88+y/bt2zl48CAbNmyoWDnVLCMiUoLCIX2///3vs337ds477zwAPvzwQ6ZPn170s5s2beIHP/gBBw4c4N1332XmzJmJjicTpXAXESlB4ZC+kydPZubMmQMjQQ7n0KFDXHfddXR2djJjxgyWL1+e2PC+xahZRkSkBIVD+s6bN4+33357YNkHH3ww8MKNyZMns3//foCBIJ86dSq9vb2sW7euouXUlbuI1KcpxyTbw2XKMbE2KxzS98Ybb+T888/npptuoqenh76+Pr71rW8xc+ZMlixZwje+8Q0mTZrEli1b+PrXv86sWbP41Kc+xRlnnJFc2YtQuItIfRqlT3oltLW1DbzYOuq0007jqaeeGrL80ksv5dJLLx2Yv+222xIb0nc0apYREckgXbnLmGi4AJHapHCXMdFwAVJN7j6kt0ojKOdhJzXLiEhdmDhxInv37q3oU521yN3Zu3cvEydOLOlzsa7czawL2A98CPS5+1wz+wTwENAGdAGXu/t7lv+zegdwIXAAWOLuz5VUKqlLSb+Qo3B/avJpbK2trXR3d/P222+nXZSKOHTo0LABPnHiRFpbW0vaXynNMp9393ci88uAJ9x9pZktC/O3AF8Ejg9fZwF3h++ScUmHb3R/avKR8ePHc+yxx6ZdjIrJ5XLMnj07sf2NpVlmIbA2TK8FLoksv9/zngGazKz4s7giIlIRca/cHfh/ZubAv7r7amCau+8J698EpoXpFmBX5LPdYdmeyDLMbCmwFKC5uZlcLldWBbKgt7e3ruufRNnjHIN6PkajqfdzIAmNfgySrn/ccF/g7rvN7C+AX5jZoF787u4h+GMLfyBWA7S3t3tHR0cpH8+UXC5H3db/8Y2JlH3UY5DQz6lVdX0OJKTRj0HS9Y8V7u6+O3x/y8z+EzgT+KOZTXf3PaHZ5a2w+W5gRuTjrWGZSNmiN1f753WDVWR4o4a7mR0BfMzd94fpLwArgPXAYmBl+P5I+Mh64AYze5D8jdSeSPONSFkKg1w3WEVGFufKfRrwn+HBgcOAn7n742b2LPCwmV0DvA5cHrZ/lHw3yJ3ku0JelXipRURkRKOGu7u/BpxaZPle4Nwiyx24PpHSiYhIWfSEqohIBincRUQySOEuIpJBCncRkQxSuIuIZJDCXUQkg/SyDilJ9M1LkMzQviKSPIW7lCT65iURqV1qlhERySCFu4hIBincRUQySOEuIpJBCncRkQxSbxmpS9GXd+jFHSJDKdylLkXDXC/uEBlKzTIiIhmkcBcRySCFu4hIBincRUQySOEuIpJBCncRkQxSuIuIZJDCXUQkgxTuIiIZpHAXEckghbuISAbFDnczG2dmvzGzDWH+WDPbamY7zewhMzs8LJ8Q5neG9W0VKruIiAyjlCv3bwI7IvP/CNzu7scB7wHXhOXXAO+F5beH7UQqpn+EyLZlG5m/8sm0iyNSE2KFu5m1AhcBPw7zBpwDrAubrAUuCdMLwzxh/blhe5GKeHrZOXStvIiulRexe9/BtIsjUhPiDvn7z8C3gclh/pPAPnfvC/PdQEuYbgF2Abh7n5n1hO3fie7QzJYCSwGam5vJ5XLl1SADent766r+lShrksegno5lv3o7Byqh0Y9B0vUfNdzN7GLgLXffZmYdSf1gd18NrAZob2/3jo7Edl13crkcdVP/xzdWpKyJHYMKla/S6uocqJBGPwZJ1z/Olft84MtmdiEwETgKuANoMrPDwtV7K7A7bL8bmAF0m9lhwBRgb2IlFhGRUY0a7u7+D8A/AIQr9//u7lea2f8GFgEPAouBR8JH1of5LWH9k+7uiZdcqmb+yicH2rJbmialXBoRiWMsr9m7BXjQzG4DfgPcG5bfC/zEzHYC7wJXjK2Ikrbd+w7StfKitIshIiUoKdzdPQfkwvRrwJlFtjkEXJZA2UREpEx6QlVEJIMU7iIiGaRwFxHJIIW7iEgGKdxFRDJI4S4ikkEKdxGRDFK4i4hkkMJdRCSDFO4iIhmkcBcRySCFu4hIBincRUQySOEuIpJBCncRkQxSuIuIZNBY3sQkGVbWq/VuPwV63shPTzkGbv5dhUonIqNRuEtRZb1ar+cNWN6Tn14+JflCiUhsCnfJlJamSbQt2zho/ull56RYIpF0KNxlbAqbYlJWGOTRoBdpJAp3Gd1IbenRphgRqRkKdxldNMBvP2Vwe3oNXK2LyFAKdylq84SbYPnf5GeiAR63B8yUY4b+EVDvGZGqUbhLUa32ztiaWwqDPKXeM9EbrLq5Ko1E4S7pCu35HQA5Er/Cj4a5bq5KI1G4S7pCe34ul6Ojo2PwFb4eihIpm8JdapceihIp26hjy5jZRDP7tZm9YGYvmtn/DMuPNbOtZrbTzB4ys8PD8glhfmdY31bhOoiISIE4A4f9CTjH3U8FTgMuMLN5wD8Ct7v7ccB7wDVh+2uA98Ly28N2IiJSRaM2y7i7A71hdnz4cuAcIPSVYy2wHLgbWBimAdYB/2JmFvYjjSraNbKc9nN1rRQpSaw2dzMbB2wDjgPuAv4A7HP3vrBJN9ASpluAXQDu3mdmPcAngXcK9rkUWArQ3NxMLpcbU0XqWW9vb83VvwOSLdPsOz/ad27hwL77f07/MYj+3EFliHy+cB+lqLXj3K8Wz4Fqa/RjkHj93T32F9AEbAIWADsjy2cA28P0dqA1su4PwNSR9nvCCSd4I9u0aVPaRRjqe0dVZ99heuAYFFk36j5i+swtG0r+TLXU5DlQZY1+DMqpP9Dpw+RqSS/rcPd9IdzPBprMrP/KvxXYHaZ3h7AnrJ8C7C3nD4+IiJQnTm+ZZjNrCtOTgPOAHeRDflHYbDHwSJheH+YJ658Mf2FERKRK4rS5TwfWhnb3jwEPu/sGM3sJeNDMbgN+A9wbtr8X+ImZ7QTeBa6oQLlFSqahCKSRxOkt81tgdpHlrwFnFll+CLgskdJJ4ynsVZMgDUUgjURPqEpe9FF/oNun0ppGOdS9USQRCnfJK3jpxoJlG+lKrzQiMkYl9ZYREZH6oCt3qb4k2tXH+sSrSMYp3KX6kgji6D40YqTIEGqWERHJIIW7iEgGKdxFRDJIbe5S/3RzVWQIhXsjKXhQaZCEnwatKt1cFRlC4d5ICh5UEpHsUpu7iEgGKdxFRDJI4S4ikkFqc8+66E3Uer5pmrDNE26C5X/z0QL1spGMUbhnXaPdRI12i+yfLxLarfbO4OOiXjaSMQp3yZbCIB8htKMv7OiaWKkCiaRD4S4Nq2vlRR/NLE+tGCIVoXAXAOavfJLd+w4OzLc0TUqxNCIyVgp3AWD3voODr2RFpK4p3KVxqOeQNBCFuzSORus5JA1N4S7ZlsQr/UTqkMJdsk0PJkmDUrhnkdqWRRqewj2L1LYs0vBGHTjMzGaY2SYze8nMXjSzb4blnzCzX5jZq+H70WG5mdkqM9tpZr81s9MrXQkRERkszqiQfcDfu/vJwDzgejM7GVgGPOHuxwNPhHmALwLHh6+lwN2Jl1pEREY0ari7+x53fy5M7wd2AC3AQmBt2GwtcEmYXgjc73nPAE1mNj3pgouIyPBKanM3szZgNrAVmObue8KqN4FpYboF2BX5WHdYtieyDDNbSv7KnubmZnK5XIlFz47e3t5E698BZe0vzd9B0segVB00dv1rQaMfg8Tr7+6xvoAjgW3AV8L8voL174XvG4AFkeVPAHNH2vcJJ5zgjWzTpk3J7vB7R5X8kc/csiHZMpQo8WNQqjKOWZJSr38NaPRjUE79gU4fJldjXbmb2Xjg34EH3P0/wuI/mtl0d98Tml3eCst3AzMiH28Ny6RSol0fQd0fy7CHZqaHh53y0ztTLpHI2Iwa7mZmwL3ADnf/p8iq9cBiYGX4/khk+Q1m9iBwFtDjHzXfSCWU2fUxOhJko48CGQ3z6Xpxh2RAnCv3+cDXgN+Z2fNh2a3kQ/1hM7sGeB24PKx7FLgQ2AkcAK5KssCSHI0EKZJdo4a7u28GbJjV5xbZ3oHrx1guEREZgzj93EVEpM4o3EVEMkhjy9QrDQ4mIiNQuNcrDQ5WUW3LNgL5XkRPLzsn5dKIlE7hLlJEfy+i/pAXqTcKd5FCkbc3bZ4wFVB3Uak/CneRQpG3N7XqgSapUwp3kRF0+9TBAT/lGL26T+qCwl1kBAv+tGrwU7y6kpc6oX7uIiIZpHAXEckgNcvUkzE+uBQdBRI0EmQcLU2TBnWH7Jo4wsaFvx+1zUuKFO71ZIwPLmkUyNINeYBpeWS62Dj6/b8ftc1LyhTuIuXSU8JSw9TmLiKSQbpyFynBoH7vGrBNapjCXaQEQ/q9DycyhIFurkoaFO4ilRANc91clRSozV1EJIMU7iIiGaRwFxHJILW517JiD8mIiMSgcK9lekhGRMqkcBcpQXSsGb1fVWqZwl2kBNEw1/tVpZYp3EUqLfpAU7F1esBJKmDUcDez+4CLgbfcfVZY9gngIaAN6AIud/f3zMyAO4ALgQPAEnd/rjJFF6kTI4W3HnCSConTFXINcEHBsmXAE+5+PPBEmAf4InB8+FoK3J1MMUUawO2n5MO+/+v2U9IukdSxUa/c3f0pM2srWLwQ6AjTa4EccEtYfr+7O/CMmTWZ2XR335NYiUWyqrB3lK7qZQzKbXOfFgnsN4FpYboF2BXZrjssGxLuZraU/NU9zc3N5HK5MotS/3p7e4vWvwMSPy61epyHOwa1bqxl7gj76O3tHbK/Dmr391UJ9XoOJCXp+o/5hqq7u5l5GZ9bDawGaG9v946OjrEWpW7lcjmK1j9H8eUliL5ar6Vp0pj3VynDHoNa9vjGsZc5l/8d9/+jHrS/3Nh///WkLs+BBCVd/3LD/Y/9zS1mNh14KyzfDcyIbNcalklcY3xPaiG9Wk+kMZUb7uuBxcDK8P2RyPIbzOxB4CygR+3tJdJTqXWj8OXZeqhJakmcrpD/Rr75b6qZdQPfIx/qD5vZNcDrwOVh80fJd4PcSb4r5FUVKLNITSgMcj3UJLUkTm+Z/zbMqnOLbOvA9WMtlIiIjI2G/BURySANPyCSpjA0QUf/dJF1A9MapkBKoHAXSVMI7KLd4PQeVhkDhXsaCl7C0QH5Z3xBL+SQ4nQVLyVSuKehoLtjoz+8ITHoKl5KpHDPoMKnUkWk8SjcM0hPpaajam9pKhwfXs00UoTCXSQhVXtLU2GQR4O+cPgKhX7DUriLVEBq71qN3s9R23xDU7iLVIDetSppU7iL1LvCbpJjpaadTFC4i1RYxUePHC58y73xqqadTFC4i1RYaqNHFgZ5/ztaQVfkDUDhLtIoynkQSk/G1i2Fu0iKBj1w9syTtfeyDz0ZW7cU7tWS8OvzpH4VdpPsWnkRuVyOJY+/X71C6Io88xTu1VLB1+dFr/5AQw7Uupq4OtcVeeYp3DNAww1IKtRlsqYp3EUaXbEuk3Goy2RNU7hXktrZpR7oijuTFO6VVKV2drWxi0ghhXudUju7VN1IzTeF60bah/6nUBUKd5EaMdIwBYX/U0ulx81IoRw3sNU2XzUK9yQVvBs16XZ2NcVkW2Fgz1/55JD+8FDno0yOcIXfAYPfJTzcHwz10olF4T5WhSdahdrYQU0xjaYm+sMnbYQgHvQu4eg4ODA4xEfqpaPgH6BwH6sK3jQVaVgjvW1qJOqeOUDhLiL1K24f/WJNphm/qq9IuJvZBcAdwDjgx+6+shI/J+s0rIAUU3jjtXBdJptzhnshyUgBXfiZ6P+wG2D448TD3czGAXcB5wHdwLNmtt7dX0r6Z1VN4V/9qArfNFUbuxQaKbwLb8JGt62JHjflKid84/buKWzfj4oG/2g5UGy7uDeGKeGGckyVuHI/E9jp7q8BmNmDwEKgvsK9ijdK/z53gL2PD+0VIVKqwjAv7FrZf24Vrosj7h+LYv/jLPUPSRL7iG2kEC28wh8uB4bbbrQ/HJH9jXhDuQzm7mPawZAdmi0CLnD3vwvzXwPOcvcbCrZbCiwNs7OA7YkWpL5MBd5JuxApa/Rj0Oj1Bx2Dcur/GXdvLrYitRuq7r4aWA1gZp3uPjetsqSt0esPOgaNXn/QMUi6/h9LakcRu4EZkfnWsExERKqkEuH+LHC8mR1rZocDVwDrK/BzRERkGIk3y7h7n5ndAPxf8l0h73P3F0f52Oqky1FnGr3+oGPQ6PUHHYNE65/4DVUREUlfJZplREQkZQp3EZEMqmq4m9kFZvaKme00s2VF1k8ws4fC+q1m1lbN8lVajPovMbO3zez58PV3aZSzUszsPjN7y8yKPtNgeavC8fmtmZ1e7TJWUoz6d5hZT+T3/z+qXcZKMrMZZrbJzF4ysxfN7JtFtsn6ORDnGCRzHrh7Vb7I31z9A/CXwOHAC8DJBdtcB9wTpq8AHqpW+Wqk/kuAf0m7rBU8Bn8NnA5sH2b9hcBjgAHzgK1pl7nK9e8ANqRdzgrWfzpwepieDPy+yL+BrJ8DcY5BIudBNa/cB4YlcPc/A/3DEkQtBNaG6XXAuWZmVSxjJcWpf6a5+1PAuyNsshC43/OeAZrMbHp1Sld5Meqfae6+x92fC9P7gR1AS8FmWT8H4hyDRFQz3FuAXZH5boZWamAbd+8DeoBPVqV0lRen/gCXhv+OrjOzGUXWZ1ncY5RlZ5vZC2b2mJnNTLswlRKaXGcDWwtWNcw5MMIxgATOA91QrS0/B9rc/a+AX/DR/2KkMTxHfqyQU4E7gf+TbnEqw8yOBP4d+Ja7/1fa5UnDKMcgkfOgmuEeZ1iCgW3M7DBgCrC3KqWrvFHr7+573f1PYfbHwJwqla1WNPTQFe7+X+7eG6YfBcab2dSUi5UoMxtPPtQecPf/KLJJ5s+B0Y5BUudBNcM9zrAE64HFYXoR8KSHOwwZMGr9C9oWv0y+Pa6RrAf+NvSYmAf0uPuetAtVLWb2qf57TGZ2Jvl/n1m5uCHU7V5gh7v/0zCbZfociHMMkjoPqjYqpA8zLIGZrQA63X09+Ur/xMx2kr/xdEW1yldpMet/k5l9GegjX/8lqRW4Aszs38j3BJhqZt3A94DxAO5+D/Ao+d4SO4EDwFXplLQyYtR/EXCtmfUBB4ErMnRxAzAf+BrwOzN7Piy7FTgGGuMcIN4xSOQ80PADIiIZpBuqIiIZpHAXEckghbuISAYp3EVEMkjhLiKSQQp3EZEMUriLiGTQ/wdOfee9GdTgJwAAAABJRU5ErkJggg==\n",
"text/plain": [
"