{
 "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare limits Bayesian (two diff. prior) and frequentist <br>\n",
    "on counting experiments with background <br>\n",
    "No uncertainty on the background.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are counts and backgrouns\n",
    "N = [0,     2,    0,   2,   0, 2]\n",
    "B = [0.5, 0.5,  2.0, 2.0, 3.5, 3.5]\n",
    "\n"
   ]
  },
  {
   "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)$ <br> <br>\n",
    "\n",
    "where $\\Gamma(a,b) = \\int_b^{\\infty} t^{a-1}e^{-t} dt$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 0 and B = 0.5\n",
      "Frequentist                  limit =  2.5\n",
      "Bayesian flat prior          limit =  3.0\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  2.4\n",
      " \n",
      "N = 2 and B = 0.5\n",
      "Frequentist                  limit =  5.8\n",
      "Bayesian flat prior          limit =  5.8\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  5.1\n",
      " \n",
      "N = 0 and B = 2.0\n",
      "Frequentist                  limit =  1.0\n",
      "Bayesian flat prior          limit =  3.0\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  2.7\n",
      " \n",
      "N = 2 and B = 2.0\n",
      "Frequentist                  limit =  4.3\n",
      "Bayesian flat prior          limit =  4.8\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  4.3\n",
      " \n",
      "N = 0 and B = 3.5\n",
      "Frequentist                  limit =  -0.5\n",
      "Bayesian flat prior          limit =  3.0\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  2.7\n",
      " \n",
      "N = 2 and B = 3.5\n",
      "Frequentist                  limit =  2.8\n",
      "Bayesian flat prior          limit =  4.3\n",
      "Bayesian 1/sqrt(S+B) prior   limit =  3.9\n",
      " \n"
     ]
    }
   ],
   "source": [
    "for n,b in zip(N,B):\n",
    "\n",
    "    print(\"N =\",n, \"and B =\",b)\n",
    "    # Frequentist first.\n",
    "    s   = -b\n",
    "    ds  = 0.01\n",
    "    while gammaincc(n+1, s+b) > 0.05:\n",
    "        s = s + ds\n",
    "    print('Frequentist                  limit = ', round(s,1))\n",
    "        \n",
    "    # Now Bayesian. A bit brute force \n",
    "    # (25 is about infty for our purposes)\n",
    "    sarray        = np.linspace(0, 25, 2500)\n",
    "    integrand     = np.power((sarray+b),n) * np.exp(-(sarray+b))\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",
    "    integrand     = np.power((sarray+b),n) * np.exp(-(sarray+b)) / np.sqrt(sarray+b)\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+B) prior   limit = ', round(s,1))\n",
    "    \n",
    "    print(\" \")\n",
    "    \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}