{ "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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# pick alphas and beta, make sure that lognormal picks make sense\n", "np.random.seed(12345)\n", "alpha = np.random.lognormal(0, np.log(1+dalpha), Ntoys)\n", "beta = np.random.lognormal(0, np.log(1+dbeta), Ntoys)\n", "ax = plt.subplot(111)\n", "bins = np.linspace(max(0,1-4*max(dalpha,dbeta)), 1+4*max(dalpha,dbeta), 101)\n", "ax.hist(alpha, bins, label='alpha', histtype='step')\n", "ax.hist(beta, bins, label='beta', histtype='step')\n", "ax.grid()\n", "ax.set_xlim(bins[0],bins[-1])\n", "ax.legend()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "N = 0 B = 0.5 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 2.6\n", "CLs limit = 3.1\n", "Bayesian flat prior limit = 3.2\n", "Bayesian 1/sqrt(s) prior limit = 1.7\n", " \n", "N = 2 B = 0.5 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 6.1\n", "CLs limit = 6.1\n", "Bayesian flat prior limit = 6.3\n", "Bayesian 1/sqrt(s) prior limit = 5.0\n", " \n", "N = 0 B = 2.0 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 1.1\n", "CLs limit = 3.1\n", "Bayesian flat prior limit = 3.2\n", "Bayesian 1/sqrt(s) prior limit = 1.7\n", " \n", "N = 2 B = 2.0 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 4.5\n", "CLs limit = 5.1\n", "Bayesian flat prior limit = 5.3\n", "Bayesian 1/sqrt(s) prior limit = 3.5\n", " \n", "N = 0 B = 3.5 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 0.0\n", "CLs limit = 3.1\n", "Bayesian flat prior limit = 3.2\n", "Bayesian 1/sqrt(s) prior limit = 1.7\n", " \n", "N = 2 B = 3.5 dalpha = 0.2 dbeta = 0.4\n", "Frequentist limit = 3.1\n", "CLs limit = 4.5\n", "Bayesian flat prior limit = 4.8\n", "Bayesian 1/sqrt(s) prior limit = 2.9\n", " \n" ] } ], "source": [ "for n,b in zip(N,B):\n", "\n", " print(\"N =\",n, \"B =\",b, \"dalpha =\", dalpha, \"dbeta =\", dbeta)\n", " # Start with S=0 and increase it in steps.\n", " # For each S:\n", " # Find average p-value for different alpha,beta\n", " # stop once we reach the 0.05\n", " # Note: we start at S=0 to avoid numerical issues \n", " # since the 2nd parameter of gammaincc must be >-0 \n", " s = 0\n", " ds = 0.01\n", " pav = 999.\n", " while pav > 0.05:\n", " temp = gammaincc(n+1, alpha*s+beta*b)\n", " pav = np.average(temp)\n", " s = s + ds\n", " print('Frequentist limit = ', round(s,1))\n", " \n", " # CL_S now\n", " s = 0\n", " ds = 0.01\n", " pav = 999.\n", " while pav > 0.05:\n", " temp = gammaincc(n+1, alpha*s+beta*b)/gammaincc(n+1, beta*b)\n", " pav = np.average(temp)\n", " s = s + ds\n", " print('CLs limit = ', round(s,1))\n", " \n", " # Bayesian. MC integration. Flat prior.\n", " # (50 is about infty for our purposes)\n", " sarray = np.linspace(0,50,5000)\n", " integrand = np.zeros( len(sarray) )\n", " for d1,d2 in zip(alpha,beta):\n", " sPlusB = d1*sarray + d2*b\n", " integrand = integrand + np.power(sPlusB, n) * np.exp(-sPlusB) \n", " totalIntegral = integrand.sum() # up to factor ds\n", " integral = np.cumsum(integrand) / totalIntegral\n", " index = (integral>0.95).argmax() # index of limit\n", " s = sarray[index]\n", " print('Bayesian flat prior limit = ', round(s,1))\n", " \n", " # Bayesian. MC integration. 1/sqrt(s) prior\n", " # (50 is about infty for our purposes)\n", " sarray = np.linspace(0.0001,50,5000)\n", " integrand = np.zeros( len(sarray) )\n", " for d1,d2 in zip(alpha,beta):\n", " sPlusB = d1*sarray + d2*b\n", " integrand = integrand + np.power(sPlusB, n) * np.exp(-sPlusB)/np.sqrt(sarray) \n", " totalIntegral = integrand.sum() # up to factor ds\n", " integral = np.cumsum(integrand) / totalIntegral\n", " index = (integral>0.95).argmax() # index of limit\n", " s = sarray[index]\n", " print('Bayesian 1/sqrt(s) prior limit = ', round(s,1))\n", "\n", " print(' ')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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.9.5" } }, "nbformat": 4, "nbformat_minor": 4 }