{
 "metadata": {
  "name": "ps10"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Counting experiments with unknown background rate\n",
      "\n",
      "In class, we worked out that the posterior pdf on the foreground rate after an ON/OFF counting experiment was (assuming uniform priors for both the foreground and background rate)\n",
      "$$\n",
      "f(s|k_{\\rm on},k_{\\rm off},I) \\propto \\sum_{j=0}^{k_{\\rm on}} \\frac{(k_{\\rm on}+k_{\\rm off}-j)!}{(k_{\\rm on}-j)!j!}\n",
      "  \\left(1+\\frac{T_{\\rm off}}{T_{\\rm on}}\\right)^j(sT_{\\rm on})^j e^{-sT_{\\rm on}}\n",
      "$$\n",
      "Let's define a function to plot this.  Note that we'll be interested in largish values of $k_{\\rm off}$ so we won't literally want to take the factorial, but instead the logarithm of the factorial, using the `gammaln` function which takes $\\ln\\Gamma(n)=\\ln([n-1]!)$.  Also, we're free to include additional constant factors, so we keep the numbers manageable by writing\n",
      "$$\n",
      "f(s|k_{\\rm on},k_{\\rm off},I) \\propto \\sum_{j=0}^{k_{\\rm on}} \\frac{(k_{\\rm on}+k_{\\rm off}-j)!}{k_{\\rm off}!(k_{\\rm on}-j)!j!}\n",
      "  \\left(1+\\frac{T_{\\rm off}}{T_{\\rm on}}\\right)^j(sT_{\\rm on})^j e^{-sT_{\\rm on}}\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from __future__ import division\n",
      "rcParams['figure.figsize'] = (10.0, 8.0)\n",
      "rcParams['font.size'] = 14\n",
      "from scipy.special import gammaln"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fgpost(sTon,kon,koff,Toff_Ton,dsTon=1.):\n",
      "    # this takes as input an array of s*Ton values\n",
      "    # the values of kon and koff\n",
      "    # and the ratio Toff_Ton\n",
      "    # and optionally the spacing of the s*Ton points\n",
      "    # It returnes the pdf for s*Ton\n",
      "    j = arange(kon+1)\n",
      "    lncoeff = gammaln(kon+koff-j+1) - gammaln(koff+1) - gammaln(kon-j+1) - gammaln(j+1) + j * log(1+Toff_Ton)\n",
      "    # that is a kon-element array\n",
      "    # we use numpy array broadcasting to turn this into a 1 x kon array and combine it with a ns x 1 array\n",
      "    # to get a ns x kon array and then sum over the last index\n",
      "    pdf = sum( sTon[:,None]**j[None,:] * exp(lncoeff[None,:]-sTon[:,None]), axis=1 )\n",
      "    return pdf/(dsTon*sum(pdf))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As an example, here is the case where $k_{\\rm on}=7$, $k_{\\rm off}=17$, and $T_{\\rm off}/T_{\\rm on}=3.5$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kon = 7\n",
      "koff = 17\n",
      "Toff_Ton = 3.5\n",
      "bhatTon = koff/Toff_Ton\n",
      "deltabTon = sqrt(koff)/Toff_Ton\n",
      "shatTon = kon - bhatTon\n",
      "deltasTon = sqrt(kon+deltabTon**2)\n",
      "sTonmin = max(0,shatTon-5.*deltasTon)\n",
      "sTonmax = shatTon+5.*deltasTon\n",
      "sTon,dsTon = linspace(sTonmin,sTonmax,1000,retstep=True)\n",
      "print \"Roughly, b*Ton = %g +/- %g and s*Ton = %g +/- %g\" % (bhatTon,deltabTon,shatTon,deltasTon)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(sTon,fgpost(sTon,kon,koff,Toff_Ton,dsTon))\n",
      "xlim([sTonmin,sTonmax])\n",
      "xlabel(r'$sT_{\\rm on}$')\n",
      "ylabel(r'$T_{\\rm on}f(s|k_{\\rm on},k_{\\rm off},I)$');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**EXERCISE** Write a function which calculates the posterior pdf $f(s|b,k_{\\rm on},I)$ for the on-source rate, given a known off-source rate."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**EXERCISE** Plot the pdfs for $s T_{\\rm on}$ for the following cases on the same set of axes:\n",
      "\n",
      "1. $k_{\\rm on}=4$, zero background rate\n",
      "\n",
      "2. $k_{\\rm on}=13$, known background rate of $b=9/T_{\\rm on}$\n",
      "\n",
      "3. $T_{\\rm off}=10 T_{\\rm on}$, $k_{\\rm on}=13$, $k_{\\rm off}=90$\n",
      "\n",
      "4. $T_{\\rm off}=T_{\\rm on}$, $k_{\\rm on}=13$, $k_{\\rm off}=9$\n",
      "\n",
      "In each case, also work out the naive maximum likelihood estimate $\\widehat{s}=\\frac{k_{\\rm on}}{T_{\\rm on}}-\\widehat{b}$ and associated error $\\delta s=\\sqrt{\\frac{k_{\\rm on}}{T_{\\rm on}^2}+(\\delta b)^2}$, where $\\widehat{b}$ is either $\\frac{k_{\\rm off}}{T_{\\rm off}}$ or the known (possibly zero) value of $b$, and $\\delta b$ is either $\\frac{\\sqrt{k_{\\rm off}}}{T_{\\rm off}}$ or zero, as appropriate.  Comment on the relationship of these numbers to the features of the pdfs you've plotted.  (It's okay to eyeball the plots; you don't need to calculate widths and peaks numerically.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**EXERCISE** Consider the case where $T_{\\rm off}=10 T_{\\rm on}$, $k_{\\rm off}=90$, and $k_{\\rm on}=7$.  What are $\\widehat{s}$ and $\\delta s$?  Plot $f(s|k_{\\rm on},k_{\\rm off},I)$.  Comment on the situation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}