{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identifikace píku\n",
    "-------------------\n",
    "\n",
    "*F. James, kap. 11.5*\n",
    "\n",
    "\n",
    "chceme rozhodnout, zda na pozadí $b(x,\\theta)$ existuje nějaká \"jemná struktura\" (spektrální pík, změna četnosti v čase...) v intervalu AB\n",
    "\n",
    "Pozadí je fitováno nějakou funkcí (s vyloučením int. AB) s odhadovanými parametry $\\hat\\theta$ (a kovar. maticí $V$), potom očekávaná hodnota pro *nulovou hypotézu* je\n",
    "$$\\hat{b}_{AB} = \\int_A^B b(x,\\hat\\theta) dx$$\n",
    "s rozptylem $\\sigma^2 _{AB} = D^T V D$, kde $$D_i=\\frac{\\partial \\hat{b}_{AB}}{\\partial \\theta_i}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testovací statistika \n",
    "\n",
    "$$T=\\frac{(n_{AB}-\\hat{b}_{AB})^2}{V(n_{AB}-\\hat{b}_{AB})}$$\n",
    "\n",
    "kde varianci můžeme uvažovat podle nulové hypotézy a Poissonovy statistiky rovnu\n",
    "\n",
    "$$V(n_{AB})=E(n_{AB})=\\hat{b}_{AB}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pokud $n_{AB}$ a $\\hat{b}_{AB}$ jsou nekorelovány, \n",
    "platí $V(n_{AB}-\\hat{b}_{AB})=\\hat{b}_{AB} + \\sigma^2 _{AB}$.\n",
    "\n",
    "Pro velká $n_{AB}$ má $T$ rozdělení $\\chi^2(1)$ s následujícími `odds` (převrácená hodnota pravděpodobnosti, že struktura dosáhne významnosti alespoň daného počtu `sigma`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "odds: [(1, 3.1514871), (2, 21.977894), (3, 370.39835), (4, 15787.192), (5, 1744277.9)]\n"
     ]
    }
   ],
   "source": [
    "from scipy import stats\n",
    "d=r_[1:6]\n",
    "print(\"odds:\",zip(d,1/stats.chi(1).sf(d).astype(\"single\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pokud nevíme přesně, kde pík (událost) hledat, je pravděpodobnost náhodného nalezení signálu v daném intervalu (obsahujícím *k* možných intervalů) větší než výše uvedená hodnota $p$\n",
    "\n",
    "$$q=1-(1-p)^k \\approx kp $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pravděpodobnost 3-sigma události, pokud máme na výběr z 50 (časových, spektrálních atp.) binů, je najednou 13%. Pokud ovšem takových událostí zaregistrujeme více, můžeme se ptát, jaká je pravděpodobnost, že (stále při platnosti nulové hypotézy) překročí např. $l$ binů z $k$ významnost s rizikem $p$:\n",
    "\n",
    "$${k \\choose l} p^l (1-p)^{k-l}$$\n",
    "\n",
    "a pravděpod., že jich bude alespoň $l$\n",
    "\n",
    "$$\\sum_{i=l}^{k} {k \\choose l} p^l (1-p)^{k-l} = 1-\\sum_{i=0}^{l-1} {k \\choose l} p^l (1-p)^{k-l}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  3.94395797e-01,   3.75615045e-01,   1.69921092e-01,\n",
       "         4.85488834e-02,   9.82536925e-03,   1.49719912e-03,\n",
       "         1.78237991e-04,   1.69750468e-05,   1.31354528e-06,\n",
       "         8.33997006e-08])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#prvnich 10 clenu rady\n",
    "from scipy import misc\n",
    "k=20\n",
    "l=r_[:10]\n",
    "p=1/22.\n",
    "prob20=misc.comb(k,l)*pow(p,l)*pow(1-p,k-l)\n",
    "prob20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, '0.606'),\n",
       " (1, '0.23'),\n",
       " (2, '0.0601'),\n",
       " (3, '0.0115'),\n",
       " (4, '0.00169'),\n",
       " (5, '0.000197'),\n",
       " (6, '1.84e-05'),\n",
       " (7, '1.4e-06'),\n",
       " (8, '8.8e-08'),\n",
       " (9, '4.56e-09')]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#riziko !nahodneho\" vyskytu alespon l binu pres 2 sigma\n",
    "zip(list(l),[\"%.3g\"%g for g in (1-cumsum(prob20))])"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}