gistfile1.txt

{
 "metadata": {
  "name": "statistics"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Statistics"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## What is this?\n",
      "\n",
      "Here we compare significance and excess sensitivity and limits for an on-off measurement computed with the [Rolke](http://adsabs.harvard.edu/abs/2005NIMPA.551..493R) and [Li & Ma](http://adsabs.harvard.edu/abs/1983ApJ...272..317L) methods.\n",
      "\n",
      "* Given `n_on`, `n_off` and `alpha`, compute `significance(n_on, n_off, alpha)`.\n",
      "* Given `n_off` and `alpha` and `significance`, compute `sensitivity(n_off, alpha, significance)`.\n",
      "* Given `n_on`, `n_off` and `alpha`, compute `excess_limits(n_on, n_off, alpha, significance)`.\n",
      "\n",
      "Note that significance and p-value = confidence level are equivalent via the normal distribution:\n",
      "\n",
      "    p_value = scipy.stats.norm.sf(significance)\n",
      "    significance = scipy.stats.norm.isf(p_value)\n",
      "    confidence_level = p_value\n",
      "\n",
      "## TODOs\n",
      "\n",
      "* In the Rolke TS profile below, why is the minimum at TS = 0, but in the Rolke paper at approx. 9?\n",
      "* Check and document where p_values are one- or two-sided.\n",
      "* Finish implementing the functions\n",
      "* Check against HESS code.\n",
      "* Compute the results with HESSE and MINOS instead of Rolke."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from numpy import log, sqrt, abs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Rolke Implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Implements the formulae on page 6 of the Rolke 2005 paper:\n",
      "http://adsabs.harvard.edu/abs/2005NIMPA.551..493R"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_mu_best(x, y, tau):\n",
      "    return x - y / tau\n",
      "\n",
      "def compute_b_best(y, tau):\n",
      "    return y / tau\n",
      "\n",
      "def compute_b_hat(x, y, tau, mu):\n",
      "    term1 = x + y - (1 + tau) * mu\n",
      "    term2 = term1 ** 2 + 4 * (1 + tau) * y * mu\n",
      "    term3 = 2 * (1 + tau)\n",
      "    return (term1 + sqrt(term2)) / term3\n",
      "\n",
      "def compute_L(x, y, tau, mu, b):\n",
      "    from scipy.stats import poisson\n",
      "    P1 = poisson.pmf(x, mu + b)\n",
      "    P2 = poisson.pmf(y, tau * b)\n",
      "    return P1 * P2\n",
      "\n",
      "def compute_lambda(x, y, tau, mu):\n",
      "    mu_best = compute_mu_best(x, y, tau)\n",
      "    b_best = compute_b_best(y, tau)\n",
      "    b_hat = compute_b_hat(x, y, tau, mu)\n",
      "    a = compute_L(x, y, tau, mu, b_hat)\n",
      "    b = compute_L(x, y, tau, mu_best, b_best)\n",
      "    return a / b\n",
      "\n",
      "def compute_TS(x, y, tau, mu):\n",
      "    lambda_ = compute_lambda(x, y, tau, mu)\n",
      "    return - 2 * log(lambda_)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reproduce Figure 1 in the Rolke paper\n",
      "x, y, tau, cl = 8, 15, 5.0, 0.95\n",
      "mu = np.arange(0, 14, 0.1)\n",
      "TS = compute_TS(x, y, tau, mu)\n",
      "plot(mu, TS);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAD9CAYAAAB6DlaSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X98zXX/x/HHMAkp5FdGRLIx5le6RA1RhMgqvyK/UnRd\nKv1y1VX6wbDLF1Fd5UKUUnGFhCKdSr6ShhTiYvs20SKMGWbb+f7xjiS2s7Nzzvt8znnebze3Rs7n\n88T22vu8Pu8fEW63242IiAS1YrYDiIhIwVSsRUQcQMVaRMQBVKxFRBxAxVpExAFUrEVEHKDAYn34\n8GESEhKIjo4mJiaGdevWBSKXiIicpURBv2HkyJF07tyZBQsWkJOTw7FjxwKRS0REzhKR36KYjIwM\nmjRpwu7duwOZSUREzpHvyDolJYVKlSoxcOBANm/eTLNmzZg6dSqlS5cGICIiIiAhRURCTWEXj+fb\ns87JySE5OZnhw4eTnJxMmTJlGD9+/J9u6NQfzzzzjPUMym8/h/I774eTs7vd3u3wkW+xjoqKIioq\nihYtWgCQkJBAcnKyVzcSERHv5Vusq1atSo0aNdixYwcAq1atokGDBgEJJiIivytwNsi0adPo27cv\n2dnZ1KlTh9mzZwciV0DEx8fbjlAkym+X8tvj5Ozeync2SIEvjojwuv8iIhKuvKmdWsEoIuIAKtYi\nIg6gYi0i4gAq1iIiDqBiLSLiACrWIiIBcvQoeDuBrsB51iIiUnRuN/TuDXfc4d3rNbIWEQmA//wH\ndu82BdsbWhQjIuJnR45ATAy89RbccIN3tVPFWkTEz0aOhMxMmDnT/Nyb2qmetYiIH23YAO+8A99/\nX7TrqGctIuInOTkwdCgkJUHFikW7loq1iIifTJliinS/fkW/lnrWIiJ+kJICLVrAunVQt+4f/592\n3RMRCQJuNwwfDqNG/blQe0vFWkTEx955B/bsgUce8d011QYREfGhQ4egQQOzCOa6687/ezTPWkTE\nsqFDoWRJeOmlC/8ezbMWEbHo889h+fKiz6k+H/WsRUR84ORJuPdeePFFuPRS319fxVpExAfGj4dr\nroEePfxzffWsRUSKaPt2aN0aNm6EGjUK/v2aZy0iEmB5eTBsGDz9tGeF2lsq1iIiRTB7Nhw/DiNG\n+Pc+aoOIiHgpPR1iY2HlSmjc2PPXaZ61iEgA9e5tWh8TJxbudZpnLSISIB98AF9//fuBAv6mYi0i\nUkgZGWajprlzoXTpwNyzwDZIrVq1KFeuHMWLFycyMpL169f//mK1QUQkDN1/vzlYYMYM717vlzZI\nREQELpeLChUqeJdKRCSEfP65aYF8911g7+vR1D2NnkVEzBS9IUNg+nS47LLA3tujkfVNN91E8eLF\nGTZsGEOHDv3D/x8zZsyZj+Pj44mPj/d1RhGRoPDcc2aKXvfuhXudy+XC5XIV6d4F9qz37dtHtWrV\n2L9/Px06dGDatGm0adPGvFg9axEJExs3wi23wObNULVq0a7ll+Xm1apVA6BSpUr06NHjDw8YRUTC\nQU4ODB4MEyYUvVB7K99inZWVxdGjRwE4duwYH3/8MbGxsQEJJiISLCZNgssvhwED7GXIt2ednp5O\nj9/2+8vJyaFv37507NgxIMFERILBzp2QlGQWwERE2Muh5eYiIheQlwft2pkHig8+6LvraotUEREf\nmjEDTpyAv/7VdhKNrEVEzuunnyAuDlwuc1q5L2lkLSLiA263WVI+YoTvC7W3tJGTiMg55s+H3bvh\nvfdsJ/md2iAiImdJT4dGjWDpUmjRwj/30OEDIiJFlJAAV18NiYn+u4cOHxARKYL33oOtW+HNN20n\n+TONrEVEgP37zXmKixdDy5b+vZfaICIiXrrrLrjyysKfp+gNtUFERLywcKHZTe/1120nuTCNrEUk\nrB04YNofCxdCq1aBuafaICIihdSnD1SrZnbWCxS1QURECmHRItiwATZtsp2kYBpZi0hYOngQGjaE\nd9+F1q0De2+1QUREPHT33VCxIkyZEvh7W2mDZGZC2bJFvYqISOB88AH87/+aGSBOUeRd92rUgCNH\nfBFFRMT/Dh0yO+rNnAllythO47kiF+vrroOPPvJFFBER/3vwQejRA2680XaSwilyse7SxbylEBEJ\ndosXw5df+neTJn8p8gPG1FQ3zZqZbQWLF/dlNBER39m/Hxo3tjP741xWToq58kqoXt0060VEgtHp\nk1/69bNfqL3lk2O9unZVK0REgtfbb8P27fDcc7aTeE/FWkRC2t695qHinDlQqpTtNN7zSbFu0cKs\nBtq1yxdXExHxDbcbBg82B982a2Y7TdH4pFgXKwa33qrRtYgElxkzzIPFv//ddpKi80mxBuje3UyL\nEREJBrt3myI9Zw5ERtpOU3Q+2xvk+HGoWtX8BVWs6LN8IiKFlpcHbdua52mPPGI7zZ9Zmbp32sUX\nQ/v28OGHvrqiiIh3pk41Bfuhh2wn8R2PinVubi5NmjSha9eu+f6+224z+8OKiNiybRuMHWuO6Aql\nhXoeFeupU6cSExNDREREvr+vSxf45BPTEhERCbScHOjfH55/HurUsZ3Gtwos1nv27GHZsmUMGTKk\nwB5LxYrQtCmsWuWzfCIiHktMhPLl4b77bCfxvQL3s37ooYdISkriyAX2QR0zZsyZj+Pj47nttngW\nLTKNfRGRQElOhmnTzH8LaAIEnMvlwuVyFeka+c4GWbp0KcuXL+ell17C5XIxadIkPjhrMvX5nmim\npsK115pVQyV0wqOIBEBWlln08uSTZv+PYOfz2SBr165lyZIl1K5dm969e7N69Wr69++f7wVr1TKb\nO33+eaFyiIh47fHHIS4O+va1ncR/PJ5n/dlnn/HPf/6zwJE1wPjx8OOP8PLLvgsqInI+K1bAvfea\nI7rKl7edxjN+n2dd0GyQ03r2hPffh9zcQmURESmUAwfM3h9z5jinUHvLb6ebx8WZZn+bNl5nExG5\nILcbbr8d6taFpCTbaQrH6grGc/XsCQsW+OvqIhLuZs8221u88ILtJIHht5H1tm3QsSP83/+ZXflE\nRHxl1y5o2RI+/RRiY22nKbygGllHR0O5crB+vb/uICLhKCfHTM976ilnFmpv+XXMq1aIiPhaYiKU\nLQt/+5vtJIHltzYIwLffms2ddu8OvhVFIuI869eb1dHJyeagbqcKqjYImLcokZHmL1ZEpCgyM037\n46WXnF2oveXXYh0RoVaIiPjGqFHQqhUkJNhOYoff52kkJJhi7X2zRUTC3ZIlsHIlvPii7ST2+L1Y\nN21qnt5u2eLvO4lIKNq3D4YNg7lzzQyzcOX3Yn26FfLee/6+k4iEmrw8uPtuU6xbt7adxq6ALFfp\n1Qvmz1crREQKJykJTp40c6rDXUCKdbNm5r/ffBOIu4lIKFi/Hv7nf2DePO2NDwEq1hER0KcPvPVW\nIO4mIk535Aj07m22Wa5Z03aa4ODXRTFn274d2rWDtLTQOnFYRHzL7TbzqcuWhVdftZ3GP4JuUczZ\n6teHatXgs88CdUcRcaI33oCNG2HyZNtJgktA98Pr3VutEBG5sB07zOKX+fOhdGnbaYJLwNogAHv2\nQOPG5jDdiy7y9q4iEoqys80KxYEDYcQI22n8K6jbIABRUWa/kBUrAnlXEXGCJ580e34MH247SXAK\n+LEAaoWIyLk++si0PmbN0g6dFxLQNgjAr7/CVVeZlsgll3h7ZxEJFenp0KSJGcTFx9tOExhB3wYB\nqFgRbrgBFi0K9J1FJNjk5cGAATBoUPgUam9ZOR2xTx94+20bdxaRYDJpklkA88wztpMEv4C3QQCO\nHTMPEnbuhEqVvL27iDjZ2rXQo4dZVn7llbbTBJYj2iAAZcrArbeaBwoiEn4OHDAbvM2cGX6F2ltW\nijWYPtWcObbuLiK25OVB//6mWHfpYjuNc1gr1u3bm6fAOpRAJLxMnAgZGTB2rO0kzmKtWBcvbr67\nvv66rQQiEmhffAFTppgWaGSk7TTOYuUB42k7d5rTH/bs0T+cSKjbv98c8/faa9Cpk+00dvn8AeOJ\nEydo2bIlcXFxxMTEMHr06CIFPNfVV5sfy5f79LIiEmROH8/Vr58KtbfyLdalSpXi008/ZdOmTXz7\n7bd8+umnrFmzxqcB7rlHrRCRUJeYaKbsPv+87STOVWDPuvRv+xRmZ2eTm5tLhQoVfBrgzjth9Wrz\nFklEQs9nn8H06aZPreO5vFfgX11eXh5NmzZl165d3H///cTExPzh/48ZM+bMx/Hx8cQXcs1ouXLQ\ntavZF2DkyEK9VESCXHo69O1r3j1Xr247jT0ulwuXy1Wka3j8gDEjI4Obb76Z8ePHnynIRX3AeNrq\n1fDww7BpU5EvJSJBIjcXbrkFWraEF16wnSa4+HUF46WXXsqtt97Khg0bCh2sIPHxcPiwirVIKBk7\nFk6dgrPefEsR5FusDxw4wOHDhwE4fvw4K1eupEmTJr4PUcysaNSDRpHQsHIl/Otfpr2pPrVv5NsG\n2bJlCwMGDCAvL4+8vDzuvvtuHn300d9f7KM2CMDu3XDddWbOdcmSPrmkiFjw449w7bXmgaK2PT0/\nb2qn1UUx52rXDu67z8wQERHnOXkS2rSBO+6As8Z1cg7H7Lp3IffdZ946iYgzjRwJNWrAI4/YThJ6\ngqpYd+8OW7fCDz/YTiIihTVnDrhcMHu2zlH0h6Aq1iVLmmPoX3vNdhIRKYxNm8xoeuFCs3ZCfC+o\netZgHjS2bGkeUlx8sU8vLSJ+cOgQNG9upur16mU7jTM4vmcN5uTz5s1hwQLbSUSkIKc3aOraVYXa\n34KuWIMeNIo4xbhxZkFbUpLtJKEv6NogADk5UKsWLFsGjRr5/PIi4gMrVsCgQbBhA1xxhe00zhIS\nbRAwK56GDIFXX7WdRETOZ+dOc9LTu++qUAdKUI6swaxkbNTIPGgsW9YvtxARLxw9alYb//WvpmUp\nhRcyI2uAqCi44QazZFVEgsPpk8mvvx6GDbOdJrwEbbEG8137lVfAT4N3ESmksWPNHtXTpmnhS6AF\ndbHu2BGOHIG1a20nEZElS8xzpIUL4aKLbKcJP0FdrIsVM3sNTJliO4lIeNu+3Tz0X7AAqlWznSY8\nBe0DxtMyM+HKK+Gbb8x0PhEJrIwMs6r4scfMVD0pOsdvkXohp3fw+uc//X4rETlLXh7cdpsZME2f\nbjtN6AjZYp2aCs2amf9econfbyciv/n732HNGvjkE4iMtJ0mdITU1L2z1aplDiaYM8d2EpHw8eab\nZurswoUq1MHAESNrgC+/hHvuMXtdF3PEtxgR51q3Drp1g9WroWFD22lCT8iOrAFatYLLLoMPP7Sd\nRCS0/fgj9OwJs2apUAcTxxTriAh48EFN4xPxp2PHzAPFhx6CLl1sp5GzOaYNApCdDbVrw/Ll2o1P\nxNfy8iAhwbyDnTlTKxT9KaTbIGCO/RoxAiZPtp1EJPQ8/TTs32+2eFChDj6OGlkDHDwIdevC5s3m\nFGURKbp58+Cpp+Crr6ByZdtpQl/IzrM+1yOPmAMK1L8WKbp168yxXKtXQ2ys7TThIWyK9d695in1\nDz9ApUoBv71IyEhJMTOtXnvNFGwJjJDvWZ92xRVwxx3w4ou2k4g416FD0LmzWaWoQh38HDmyBti1\ny2wus3s3lCtnJYKIY2Vnwy23QOPGemBvQ9iMrAHq1DH7XesUdJHCcbth6FAzyNHmaM6Rb7FOS0uj\nbdu2NGjQgIYNG/JikPUdnnjCjApOnLCdRMQ5nnsOtm41M0CKF7edRjyVbxvk559/5ueffyYuLo7M\nzEyaNWvGokWLiI6ONi+22AY5rWtX03e7/36rMUQc4Y034B//MDNAqla1nSZ8+bwNUrVqVeLi4gAo\nW7Ys0dHR7N271/uEfjB6NEycaKbyiciFuVwwapTZX0eF2nlKePobU1NT2bhxIy1btvzDr48ZM+bM\nx/Hx8cTHx/sqm0datTIbo7/9Ntx9d0BvLeIY27fDXXeZr5MGDWynCT8ulwuXy1Wka3g0GyQzM5P4\n+Hieeuopunfv/vuLg6ANAmYy/7BhsG0blPD4249IeNi3D66/3iwnv+ce22kE/DQb5NSpU/Ts2ZN+\n/fr9oVAHk3btzNLzuXNtJxEJLhkZ0KkTDB6sQu10+Y6s3W43AwYMoGLFikw+z2TMYBlZgzmcoG9f\ns6rxootspxGx7+RJ8/C9fn1zfqI2ZwoePl9uvmbNGm644QYaNWpExG//0omJidxyyy1e39CfOnc2\ne/AOH247iYhdeXnQp4958P7OO5qiF2zCZm+QC/nmG3MU0X//CxdfbDuNiB1utzk8YONG+OgjKFXK\ndiI5V1itYDyfZs3MEvRXXrGdRMSepCRzGvnixSrUoSSkRtYAW7ZAhw5mdF22rO00IoE1d66Z9fHl\nl1C9uu00ciFhP7IGsx9vu3bakU/Cz4oV8Nhj5tg7FerQE3Ija4AdO8y80p07zXlyIqFu7Vro3h0W\nLTILxSS4aWT9m3r1zAnNiYm2k4j43+bN0KOH2fdDhTp0heTIGsxpMo0awddfmxPRRULRjh0QHw/T\npkHPnrbTiKc0sj7LFVfAyJFmoyeRUPTjj+Zh+tixKtThIGRH1gBZWXDNNfDuu/CXv9hOI+I76enQ\npg2MGGEGJeIsGlmfo3RpM+p4+GGzUEAkFBw6ZE5J6tdPhTqchHSxBvMJnZ1tRtciTpeZabZVuOkm\nc4iAhI+QboOc9tlnZsexbdu0okuc68QJs/dNrVowY4Y2ZnIytUEu4MYboUkTLZQR5zpxwkzPq1QJ\nXn1VhTochcXIGswCmVatzHJ0HWkkTnLypJntUaaMOeRWB2w4X9jvuleQ0aPNdKd582wnEfFMdjYk\nJEDJkuZIrshI24nEF1SsC5CVZc6f+/e/oX1722lE8nfqFNx5p/n43XdVqEOJetYFKF3a9K2HDzdv\nLUWC1alT0Ls35OaawwNUqCWsijVA164QHW32/BUJRjk5ZsrpiRPw3numBSISVm2Q0378EZo2ha++\ngjp1bKcR+V1ODvTvbxa+vP++ppqGKrVBPFSzptn394EHtLJRgsepU+bcxF9/VaGWPwvLYg3mjLq0\nNFi40HYSEdPy6NnTPEtZskSFWv4sLNsgp61ZA716mbnX5cvbTiPhKivLHBxQvjy8+aYeJoYDTd3z\nwt/+ZvqDb7xhO4mEo6NHzRLy2rVh5kwoXtx2IgkE9ay9kJgI69aZ45BEAunQIbMfdUwMzJqlQi35\nC/tiXaYMzJ5t5l4fOGA7jYSL/fvNwc6tWsHLL0OxsP9KlILoUwRo3dosQHjgAdtJJBzs22eO4urS\nBSZN0qZM4hkV69+88AJs2mQWIYj4y44dcP31cPfd8PzzKtTiubB/wHi2devMU/nNm6FKFdtpJNR8\n/TV062YGBoMH204jNmk2iA888QRs3QqLF2vUI77z8cdmCfm//20KtoQ3zQbxgeeeg59/hqlTbSeR\nUPHWW6bt8f77KtTivQKL9aBBg6hSpQqxsbGByGNdyZJml7Nx48zbVpGimDzZvFtbvdr0qkW8VWCx\nHjhwICtWrAhElqBRu7aZTnXXXZCRYTuNOJHbDY8/Dq+9ZlbKNmhgO5E4XYHFuk2bNpQPw7XYCQnQ\nqRMMGaLNnqRwTp40O+d9/rkp1DVr2k4koaDIp7mNGTPmzMfx8fHEx8cX9ZJBY9IkuO46c0DpfffZ\nTiNOsH+/Odi2alX45BNz4IWIy+XC5XIV6RoezQZJTU2la9eubNmy5Y8vDsHZIOc6PS/244/NCeki\nF7Jtm1no0quXmUOtVYlyIZoN4gf16sH06Wa09MsvttNIsFq5Em68EZ5+GsaOVaEW39OnlAfuusvM\nke3Z05w2LXK2f/3LTM1bsAAGDLCdRkJVgcW6d+/etGrVih07dlCjRg1mz54diFxB57nn4PLLYcQI\nPXAUIzfXHGIxZYp5kHjDDbYTSSjTCsZCyMw0u6QNGWL2wZbwdfiwebd1/LgZUYfhhCkpAvWs/axs\nWbMMPTHR9CglPG3dCtdea+bjr1ihQi2BoWJdSLVrw/z5ZlS1bZvtNBJoCxeaB4lPPgnTpukILgkc\ntUG8NHcu/OMf8MUXWvQQDnJzzb/3vHmmYDdvbjuROJk3tbPIi2LCVf/+5mSZm282Bfvyy20nEn/5\n5RfzTio72+wXU7my7UQSjtQGKYKHHzb7X3fubA4+ldDz+efQtKkZSa9apUIt9qgNUkRuN9x7L6Sm\nwtKlcNFFthOJL+TlwYQJZqvc2bPNPjEivqLDByzJyYE77zSnU7/9NpRQc8nRDhwwba6MDPMwuUYN\n24kk1GjqniUlSpgN5o8eNQfvnjplO5F466OPoHFjiI0Fl0uFWoKHRtY+dOKEGWFHRMC776ol4iQn\nTsDo0WaBy5w50K6d7UQSyjSytqxUKfPFHhlpHjweP247kXjiu+/MIpe0NHNYsgq1BCMVax8rWdL0\nOcuXN9tlHjtmO5FcSG4uJCVB27Zmj4/33oMKFWynEjk/tUH8JDfX7CHyww9miXqlSrYTydm2bYOB\nA83hADNnmpWpIoGiNkgQKV7cFIEbb4S//MUUbbEvJwcmTjQ75A0YYOZOq1CLE2iSmR8VK2Y2fapb\n1xSHd96BEDr1zHE2boRhw+CSS8xKxFq1bCcS8ZxG1gEweLCZ2nfnnWamgQTW0aOmJ33zzaZYr1yp\nQi3Oo2IdIO3bw2efwbPPwmOPaS52ILjd8J//QEyM2X/6++/NN04duSVOpAeMAbZ/v1kdd+SIVsf5\n0/btZu+WlBRz7NaNN9pOJPI7PWB0gEqV4MMPoVs3sznQ0qW2E4WWgwdh5Eho08a8m9m8WYVaQoOK\ntQXFisHjj5t9kYcPh0ce0UG8RXXqlDkMoH598/HWrTBqlJn3LhIKVKwtat0akpNhxw4zyv76a9uJ\nnCc3F958E6Kj4YMPYPVqePllzWuX0KOedRBwu81skVGjoE8feP55KFPGdqrg5nbD++/D009DuXIw\ndqxZiSjiBOpZO1REBPTtC1u2mFNJYmN1IO+F5OXBokXQooX5pjZhAnz5pQq1hD6NrIPQ8uWml924\nsVlUEx1tO5F92dnm/MOJE80p8088AT16aBqeOJNG1iGiUyezd0Xr1mbl47BhsG+f7VR2HDoEkybB\nVVeZgx2mT4f166FnTxVqCS/6dA9SpUqZWSI//GB6sg0bwpNPmjZJOPjmG7OA5aqrzMeLF8PHH5vp\neBERttOJBJ6KdZCrUMFs45mcbBbUXHONOfNx+3bbyXwvIwNmzYKWLc3IuW5d8+d86y1o1sx2OhG7\n1LN2mF9+MVPTXnnFPGQbOdJsll+8uO1k3jl50vTo580zI+e2bc3Wsp06OffPJFIQ9awLyeVy2Y5Q\naJUrw5gx5jT1a65x8cQTULMmPPqoWa3nhO+dR46YBUGdOrm44gqYPBk6djR/pkWLzKENTijUTvz8\nOZuT8zs5u7cKLNYrVqygfv36XH311UyYMCEQmQLGyf/gF18Ml1zi4ptvzDS/kiXNEvZGjeCZZ2Dt\nWrN3czDIyzNHZ02dCh06QPXq8NprEBHhYuNGs8HV0KHmdB0ncfLnDzg7v5OzeyvfYp2bm8sDDzzA\nihUr2Lp1K2+//Tbbtm0LVDbxUEyMWRSSkmJaJMePw4gRZhXf7bebjYw2bQrckvajR2HNGjMHuksX\nuPxyM81uyxYzJXHvXnOK+LXXmncFIlKwfA8fWL9+PXXr1qXWb5v/9urVi8WLFxOtib9BqVgxs4FR\nmzbm5+np5iSUlSvNvhkpKWbvjLg486NWLYiKMjv/XX554WZZHD0KP/1kCm9amplq+N135sf+/eYb\nyPXXwz33wIwZUK2aP/7EIuEj3weMCxYs4KOPPmLGjBkAvPnmm3z11VdMmzbNvFhzqEREvFLYB4z5\njqwLKsaaCSIiEhj59qyrV69OWlramZ+npaURFRXl91AiIvJH+Rbr5s2bs3PnTlJTU8nOzuadd96h\nW7dugcomIiK/ybcNUqJECaZPn87NN99Mbm4ugwcP1sNFERELCpxn3alTJ3744Qf++9//Mnr06DO/\n7vT512lpabRt25YGDRrQsGFDXnzxRduRCi03N5cmTZrQtWtX21EK7fDhwyQkJBAdHU1MTAzr1q2z\nHalQEhMTadCgAbGxsfTp04eTJ0/ajpSvQYMGUaVKFWJjY8/82sGDB+nQoQP16tWjY8eOHD582GLC\n/J0v/6OPPkp0dDSNGzfm9ttvJyMjw2LCCztf9tMmTZpEsWLFOHjwYIHX8WoFYyjMv46MjGTy5Ml8\n//33rFu3jpdeeslxf4apU6cSExPjyFk5I0eOpHPnzmzbto1vv/3WUe/YUlNTmTFjBsnJyWzZsoXc\n3Fzmz59vO1a+Bg4cyIoVK/7wa+PHj6dDhw7s2LGD9u3bM378eEvpCna+/B07duT7779n8+bN1KtX\nj8TEREvp8ne+7GAGjCtXruTKK6/06DpeFeuz519HRkaemX/tJFWrViUuLg6AsmXLEh0dzd69ey2n\n8tyePXtYtmwZQ4YMcdysnIyMDL744gsGDRoEmHbbpZdeajmV58qVK0dkZCRZWVnk5OSQlZVF9erV\nbcfKV5s2bSh/zhLRJUuWMGDAAAAGDBjAokWLbETzyPnyd+jQgWK/7ZPbsmVL9uzZYyNagc6XHeDh\nhx9m4sSJHl/Hq2L9008/UaNGjTM/j4qK4qeffvLmUkEhNTWVjRs30rJlS9tRPPbQQw+RlJR05pPV\nSVJSUqhUqRIDBw6kadOmDB06lKysLNuxPFahQgVGjRpFzZo1ueKKK7jsssu46aabbMcqtPT0dKpU\nqQJAlSpVSE9Pt5zIe7NmzaJz5862Y3hs8eLFREVF0ahRI49f49VXuhPfdl9IZmYmCQkJTJ06lbJl\ny9qO45GlS5dSuXJlmjRp4rhRNUBOTg7JyckMHz6c5ORkypQpE9Rvwc+1a9cupkyZQmpqKnv37iUz\nM5N58+bZjlUkERERjv26Hjt2LCVLlqRPnz62o3gkKyuLcePG8eyzz575NU++jr0q1qEy//rUqVP0\n7NmTfv3ECaD3AAACFUlEQVT60b17d9txPLZ27VqWLFlC7dq16d27N6tXr6Z///62Y3ksKiqKqKgo\nWrRoAUBCQgLJycmWU3luw4YNtGrViooVK1KiRAluv/121q5daztWoVWpUoWff/4ZgH379lG5cmXL\niQrv9ddfZ9myZY76Zrlr1y5SU1Np3LgxtWvXZs+ePTRr1oxfCjhZxKtiHQrzr91uN4MHDyYmJoYH\nH3zQdpxCGTduHGlpaaSkpDB//nzatWvH3LlzbcfyWNWqValRowY7duwAYNWqVTRo0MByKs/Vr1+f\ndevWcfz4cdxuN6tWrSImJsZ2rELr1q0bc+bMAWDOnDmOGrCAmZGWlJTE4sWLKVWqlO04HouNjSU9\nPZ2UlBRSUlKIiooiOTm54G+Wbi8tW7bMXa9ePXedOnXc48aN8/Yy1nzxxRfuiIgId+PGjd1xcXHu\nuLg49/Lly23HKjSXy+Xu2rWr7RiFtmnTJnfz5s3djRo1cvfo0cN9+PBh25EKZcKECe6YmBh3w4YN\n3f3793dnZ2fbjpSvXr16uatVq+aOjIx0R0VFuWfNmuX+9ddf3e3bt3dfffXV7g4dOrgPHTpkO+YF\nnZt/5syZ7rp167pr1qx55uv3/vvvtx3zvE5nL1my5Jm/+7PVrl3b/euvvxZ4nSKdFCMiIoHhvKkE\nIiJhSMVaRMQBVKxFRBxAxVpExAFUrEVEHEDFWkTEAf4fwiZgdAJ9vXEAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x103568990>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Li & Ma Implementation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#@np.vectorize\n",
      "def loglike_lima(on, off, alpha, signal, do_switch=True):\n",
      "    \"\"\"Taken from Utilities::Statistics::LiMa_LogLikelihood\n",
      "    from the HESS software.\"\"\"\n",
      "    #for _ in [on, off, alpha, signal]:\n",
      "    #    _ = np.asarray(_)\n",
      "    # @todo Explain what the model is for negative signal\n",
      "    # Is there a source in the off region?\n",
      "    if do_switch:\n",
      "        positive = signal > 0\n",
      "        n1 = np.where(positive, on, off)\n",
      "        n2 = np.where(positive, off, on)\n",
      "        a = np.where(positive, alpha, 1.0 / alpha)\n",
      "        signal = np.where(positive, signal, -a * signal)\n",
      "    else:\n",
      "        n1 = on\n",
      "        n2 = off\n",
      "        a = alpha\n",
      "\n",
      "    # Estimate background for the on region\n",
      "    # @note Some of the other functions compute background\n",
      "    # in the off region as the parameter\n",
      "    def case1(n1, a, signal):\n",
      "        \"\"\"Corresponds to Li & Ma formula (7)\"\"\"\n",
      "        # @todo Explain this formula better!\n",
      "        return n1 / (a + 1) - signal / a\n",
      "    def case2(n1, n2, a, signal):\n",
      "        # @todo Explain this formula!?\n",
      "        # Is it the same as nll()?\n",
      "        t = (1 + a) * signal - a * (n1 + n2)\n",
      "        delta = t * t + 4 * a * (1 + a) * signal * n2\n",
      "        return (sqrt(delta) - t) / (2 + 2 * a)\n",
      "    background = np.where(off == 0,\n",
      "                          case1(n1, a, signal),\n",
      "                          case2(n1, n2, a, signal))\n",
      "    #background = case2(n1, n2, a, signal)\n",
      "    np.where(background < 0, 0, background)\n",
      "    \n",
      "    # Compute Poisson nll for given signal / background:\n",
      "    # n1 ~ Pois(signal + background)\n",
      "    # n2 ~ Pois(background / a)\n",
      "    mu1 = signal + background\n",
      "    mu2 = background / a\n",
      "    l1 = np.where(n1 > 0, n1 * log(mu1), 0) - mu1\n",
      "    l2 = np.where(n2 > 0, n2 * log(mu2), 0) - mu2\n",
      "    l = l1 + l2\n",
      "    if 0:\n",
      "        print 'background', background\n",
      "        print 'signal', signal\n",
      "        print 'n1:', n1\n",
      "        print 'mu1:', mu1\n",
      "        print 'n2:', n2\n",
      "        print 'mu2:', mu2\n",
      "        print 'of interest', mu2[n2 > 0]\n",
      "        print 'l1', l1\n",
      "        print 'l2', l2\n",
      "        print 'l', l\n",
      "    return -1\n",
      "    #return -l, -l1, -l2\n",
      "\n",
      "@np.vectorize\n",
      "def excess_error(on, off, alpha, nsig=1,\n",
      "                 up=True, loglike=loglike_lima):\n",
      "    \"\"\"Statistical uncertainties on signal using the profile\n",
      "    likelihood method on a variant of the likelihood function\n",
      "    given in the Li & Ma paper.\n",
      "    \n",
      "    @param on: number of counts in on region\n",
      "    @param off: number of counts in off region\n",
      "    @param alpha: on / off exposure ratio\n",
      "    @param nsig: number of standard deviations of the error\n",
      "    @return: signal limit\n",
      "    \n",
      "    The profile likelihood method finds the signal for which\n",
      "    L(signal) = L0 + 0.5 * nsig, where\n",
      "    L0 = L(signal0) is the best-fit likelihood, which here is\n",
      "    signal0 = on - alpha * off, i.e. the signal.\"\"\"\n",
      "    from scipy.optimize import brentq\n",
      "    # Optimal signal and likelihood\n",
      "    signal0 = on - off * alpha\n",
      "    L0 = loglike(on, off, alpha, signal0)\n",
      "    # @todo Shouldn't dL scale with nsig^2\n",
      "    dL = nsig * 0.5\n",
      "    \n",
      "    # We want to use the brentq method to find a root of\n",
      "    # f(signal) = L0 - L(signal) - dL\n",
      "    # The brentq method needs a start interval\n",
      "    # (signal_low, signal_up)\n",
      "    # containing the solution.\n",
      "    # We use signal_low = signal0 and find a signal_up \n",
      "    # such that f(signal_up) is positive.\n",
      "    # The details of the following algorithm are not important,\n",
      "    # we start with an interval of 1 and increase by factors of 2.\n",
      "    signal = signal0 + 1 if up else signal0 - 1\n",
      "    for _ in range(100): # Avoid infinite loop (should never occur)\n",
      "        L = loglike(on, off, alpha, signal)\n",
      "        if (L0 - L < dL) & np.isfinite(L):\n",
      "            signal = signal0 + 2 * (signal - signal0)\n",
      "        else:\n",
      "            break    \n",
      "\n",
      "    # Now use brentq to find the signal\n",
      "    def f(signal):\n",
      "        L = loglike(on, off, alpha, signal)\n",
      "        return L0 - L - dL\n",
      "    signal_low = signal0 if up else signal\n",
      "    signal_up = signal if up else signal0\n",
      "    # @note The HESS software uses bisect\n",
      "    # (i.e. results will be different because we use brentq) \n",
      "    # with hardcoded xtol = 1e-4 and maxiter = 100\n",
      "    signal = brentq(f, signal_low, signal_up,\n",
      "                    xtol=1e-4, maxiter=100)\n",
      "    \n",
      "    # Return difference to excess\n",
      "    return signal - signal0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 76
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Significance"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_significance_lima(n_on, n_off, alpha):\n",
      "    \"\"\"The Li & Ma significance can be evaluated as a simple formula\"\"\"\n",
      "    n_sum = n_on + n_off\n",
      "    temp = (alpha + 1) / n_sum\n",
      "    l = n_on * log(n_on * temp / alpha)\n",
      "    m = n_off * log(n_off * temp)\n",
      "    sign = np.where(n_on - alpha * n_off > 0, 1, -1)\n",
      "    return sign * sqrt(abs(2 * (l + m)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_significance_rolke(n_on, n_off, alpha):\n",
      "    # Translate symbols used in gamma-ray astronomy to those in the Rolke paper\n",
      "    x, y, tau = n_on, n_off, 1. / alpha\n",
      "    TS = compute_TS(x, y, tau, 0.)\n",
      "    return sqrt(TS)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compare significance for one example -> Perfect agreement\n",
      "n_on, n_off, alpha = 10, 10, 0.1\n",
      "significance_lima = compute_significance_lima(n_on, n_off, alpha)\n",
      "significance_rolke = compute_significance_rolke(n_on, n_off, alpha)\n",
      "print significance_lima\n",
      "print significance_rolke"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.70512718528\n",
        "4.70512718528\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Sensitivity"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_sensitivity_lima(n_off, alpha, significance, guess=1e-3):\n",
      "    \"\"\"Compute sensitivity using the Li & Ma formula for significance.\n",
      "    There is no formula, we have to find the solution numerically.\n",
      "\n",
      "\t@note: in weird cases (e.g. on=0.1, off=0.1, alpha=0.001)\n",
      "\tfsolve does not find a solution.\n",
      "\t@todo: Is it possible to find a better starting point or\n",
      "\timplementation to make this robust and fast?\n",
      "\tvalues < guess are often not found using this cost function.\"\"\"\n",
      "    # It would be better to find a range that always contains the solution\n",
      "    # and to then call brentq, which will be a bit slower, but always converge.\n",
      "    from scipy.optimize import fsolve\n",
      "    def f(n_on):\n",
      "        if n_on >= 0:\n",
      "            return compute_significance_lima(n_on, n_off, alpha) - significance\n",
      "        else:\n",
      "            return 1e100\n",
      "    n_on = fsolve(f, guess)[0]\n",
      "    return n_on - alpha * n_off"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# TODO: x must be int; results don't match\n",
      "def compute_sensitivity_rolke(n_off, alpha, significance):\n",
      "    import ROOT\n",
      "    rolke = ROOT.TRolke()\n",
      "    rolke.SetCLSigmas(significance)\n",
      "    x, y, tau, e = int(alpha * n_off), int(n_off), 1. / alpha, 1\n",
      "    rolke.SetPoissonBkgKnownEff(x, y, tau, e)\n",
      "    low, high = ROOT.Double(0), ROOT.Double(0)\n",
      "    rolke.GetSensitivity(low, high)\n",
      "    return high"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute sensitivity for one example\n",
      "n_off, alpha, significance = 10, 0.1, 3.\n",
      "sensitivity_lima = compute_sensitivity_lima(n_off, alpha, significance)\n",
      "sensitivity_rolke = compute_sensitivity_rolke(n_off, alpha, significance)\n",
      "print sensitivity_lima\n",
      "print sensitivity_rolke"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.81849708379\n",
        "6.5305711766\n"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Limits = Errors"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_limits_rolke(n_on, n_off, alpha, significance):\n",
      "    import ROOT\n",
      "    rolke = ROOT.TRolke()\n",
      "    rolke.SetCLSigmas(significance)\n",
      "    x, y, tau, e = int(n_on), int(n_off), 1. / alpha, 1\n",
      "    rolke.SetPoissonBkgKnownEff(x, y, tau, e)\n",
      "    low, high = ROOT.Double(0), ROOT.Double(0)\n",
      "    rolke.GetLimits(low, high)\n",
      "    return low, high"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def compute_limits_lima(n_on, n_off, alpha, significance):\n",
      "    # TODO: Doesn't work yet\n",
      "    return None, None\n",
      "    down = excess_error(n_on, n_off, alpha, significance, False, loglike_lima)\n",
      "    up = excess_error(n_on, n_off, alpha, significance, True, loglike_lima)\n",
      "    return down, up"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 78
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Look at three examples: Gaussian limit, small positive excess, small negative excess\n",
      "examples = [(100, 0, 1, 1), # should roughly give excess 100 +- 10\n",
      "            (100, 0, 1, 2), # should roughly give excess 100 +- 20\n",
      "            (10, 10, 0.1, 1), # excess = 9\n",
      "            (10, 90, 0.1, 1), # excess = 1\n",
      "            (10, 110, 0.1, 1) # excess = -1\n",
      "           ]\n",
      "for n_on, n_off, alpha, significance in examples:\n",
      "    limits_lima = compute_limits_lima(n_on, n_off, alpha, significance)\n",
      "    limits_rolke = compute_limits_rolke(n_on, n_off, alpha, significance)\n",
      "    print n_on, n_off, alpha, limits_lima, limits_rolke"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 0 1 (None, None) (90.33017981414673, 110.33660955973012)\n",
        "100 0 1 (None, None) (81.31051153599459, 121.35596019830203)\n",
        "10 10 0.1 (None, None) (6.141612917334355, 12.516769078804394)\n",
        "10 90 0.1 (None, None) (0.0, 4.617504851826792)\n",
        "10 110 0.1 (None, None) (0.0, 2.6423695593692913)\n"
       ]
      }
     ],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}