Sample power analysis for firing rate task

spiking_power.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spiking power analysis\n",
    "\n",
    "This notebook illustrates a simple simulation-based spiking analysis\n",
    "\n",
    "We will assume:\n",
    "\n",
    "- a baseline firing period prior to each trial $\\lambda_0$\n",
    "- a baseline period of $T_0$ seconds\n",
    "- a stimulus presentation period on each trial $\\alpha \\lambda_0$\n",
    "- a stimulus period of $T$ seconds\n",
    "- a firing rate response that is more variable than Poisson; this is modeled as a trial-to-trial Gaussian noise on the log firing rate: $\\log \\lambda_i = \\log (\\alpha\\lambda) + \\epsilon$\n",
    "- a paired t-test on firing rates to determine significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "np.random.seed(12345)  # for reproducibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's go through it carefully and plot along the way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define some numbers\n",
    "baseline = 20.  # firing rate in Hz\n",
    "Tbaseline = 1.\n",
    "Tstim = 1.\n",
    "effect = 2  # fractional change from baseline\n",
    "trial_noise_std = 0.0  # trial-to-trial noise as a percentage\n",
    "Ntrials = 500  # number of trials for each stimulus\n",
    "Nexpt = 1000  # number of synthetic experiments to conduct\n",
    "alpha = 0.05  # nominal false positive rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make some fake firing rates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "baseline_rates = np.exp(np.log(baseline) * (1 + trial_noise_std * np.random.randn(Ntrials)))\n",
    "stim_rates = np.exp(np.log(effect * baseline) * (1 + trial_noise_std * np.random.randn(Ntrials)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGvFJREFUeJzt3XuUVeWd5vHvA0qGiiJQhCoFLRAcg87YOHaMy2in7I5E\nzHhZnRlBM+2FjDExJmSyTEQnkYodDThpY9ZqdWIMsWTFKHHZoqttQUcrGk1AbVxeULkYLpZSEotr\nyojAb/7Yu8pjWZdzijoc6vX5rFWLfd6zL+97Xs5z9n733ucoIjAzs3QNqnQFzMysvBz0ZmaJc9Cb\nmSXOQW9mljgHvZlZ4hz0ZmaJc9BbUSSdJ+mhflrXY5Jm9Me6ull/0XWVNFvS/B6e/5qkDZK2Shop\naZukcT3Mf4uk/116rc3Kx0FvHSSdJOlJSZsl/UnSE5KOA4iIOyPitArXr07Sbkk9/r/tQ127vJlE\n0n7APwGfi4hhEdEaEQdGxJoetv21iLi2hG2XjaRfSrqm0vWwytuv0hWwfYOkA4EHgEuA3wBDgJOB\ndytZr05EFsrqdgZpcETs6qft1QIfA14uZmZJgyJidz9tu7dt9Wc7LXHeo7d2/xGIiFgQmXcj4pGI\neBFA0gWSnmifOd+zvkTSCkmtkv654LlBkv5J0kZJqyV9vac9cUkzJC2X9Lakf5N0WDd1/G3+7+Z8\nKOXTeb1+J+kGSX8CZndR1xslrZO0RdLTkk7q7cWQdATwSv5wk6RHCtp9eD79S0k3S/pXSduA+sK9\naEmflbRe0rcltUhqlnRhwTZGSnogr9cSSf9YWO9O9Wk/mpkhaS3w//LyBZLelLRJUpOkSXn5xcCX\ngO/mr9XCvPxgSfdIeivvm28UbONT+euzJV/nj3t7nWxgcNBbuxXALkm3SzpN0vAu5uk8xPEF4Djg\nr4BzJE3Jy78CfB44BvgvwNldLAuApLOAWfk8nwCeAH7dTR3/Jv93WD6UsiR//GlgFTAaaB82Kdze\n0rwuI4A7gd9IGtLNNrKFI1YCR+cPD4qIz3WxXoBzgX+MiAOBJ7tYVS1wIHAI8D+BmyQdlD93M7At\nr/eFwAVdrL+zvwE+Sfb6AjwITMjX8e95+4iInwO/Aq7PX6uzJInsqG0ZcDDwd8BMSafm6/opcGNE\nHJSvc0EvdbEBwkFvAETENuAkYDdwK/CWpIWSPtHDYj+KiG0RsR54DJicl/934KcR8WZEbAHm9LCO\nS/L1rMiHPeYAkyUd2sMynYdumiPi5ojYHREfGmrKx+w358//hGw45sge1t/T9jpve2FE/CHfTlfD\nXDvIPgh2RcS/AduBI/Ojm78Hrs6Pnl4GGnupRwCzI+Kd9m1FxO0R0RYR7wHXAH+VD8N15VPAqIi4\nNq/PGuA2YHr+/HvAREnV+TqX9lIfGyAc9NYhIl6NiBkRcRjwn8j2Qm/sYZGWguk24IB8+hBgfcFz\nhdOd1QE/zYd/WoG3yQJtTAlV72n9SLo8HxraJGkTMAwYVcL6+7xt4O1O4/btr9MngMHA6yWsi8L5\n8yGyOZJWSdoM/JHsteuubXXAmPbXOn8triQ7GgCYQfYB+Eo+lPSFIupjA4BPxlqXImKFpNvJhmFK\n9SYwtuBxd2PukIXbDyOiu+GaD1SrxHIknQx8BzglIpbnZa30cEK3RH39+teNwE6y12lVXtbTUUxX\n2zsPOAP424hYlw8JbeL9tnWu23rgtYjo8mgmIlbn60TSF4F7JI2MiHeKqJftw7xHbwBIOjI/aTgm\nf3wo2fjz7/uwugVkY7+H5GP93+1h3v8LXCXpqHy7B0n6b93Mu5FsaGlCCXU5gGxI4m1JQyRdTTZm\nXqz++kD4gHwv/16gQdJQSZ8Ezi+xLgeSXRW1SdLHgR/xwXBvAQ4veLwU2Cbpu5L+g6TBko6W9NcA\nkr4kqf1oYEu+rr1yFZGVl4Pe2m0jO6m5JL+C5CngeeDybubvvLdY+PjnwOJ8+WeBfwV2FgxhdMwb\nEfeRjcvflQ8/PA90eQ18vmd5LfBkPvRwfBHtWpT/rSAb2mijuCGSrtrV1eNSFS7/DWA42RFQI9mJ\n1J4uZ+287TuAdUAz8CJZnxX6BXB0/lrdm7/+/5XsXMofgbfI+mpYPv9pwEuStgI/AaZ1c97BBhgV\n88MjktaQfcLvBt6LiOMljQDuJhv3WwOck594Q9KVZON9O4GZEbG4LLW3AUHSacAtETG+0nXZl0ma\nA9RExEWVroulpdg9+t1AfUQcGxHte1GzgEfy8b5HyU7qkB+CnwNMAqYCN+eXddlHRD4sMDUfGhgD\nzCYbprAC+XDZf86njwe+jF8nK4Nig15dzHsW718O1kh2HTTAmcBdEbEzv3xrJVDMIbalQ8APgFay\noZuXyMLePuhA4F5J28nuHfg/EfFAhetkCSr2qpsAHpa0C/hZRNxGdojZAhARGyS1X6I1hg+ewGum\ntEvlbIDLx9L94d6LiHgGOKLS9bD0FRv0n4mIN/ObZxZLepX+P0llZmZlUFTQR8Sb+b8bJd1HtrfW\nIqkmIlok1ZKdwYdsD77weuCxedkHSPIHg5lZH0RESec9ex2jl1Ql6YB8+uPAFOAF4H6y7+eA7Ds6\nFubT9wPT82uWxwMTya7f7aqyyf6l3sbZs2dXvA5un9vX1V/q772+KGaPvgb4l3wPfD/gVxGxWNIz\nwAJlPyCxluxKGyJiuaQFwHKyG1Uujb7WzszM9livQR8Rf+T9L6sqLG8FPvfhJSAifkR2l56ZmVWY\n74y1Pqmvr690FcrK7bOUFHVnbFk2LCU9otN+j1jKbTTbF6X+3pNElHgy1t9eaWY9GjduHGvXrq10\nNUo20G/Ir6urY82aNf2yLu/Rl0nqexX20ZHvQVa6Gh853b3ufdmj9xi9mVniHPRmZolz0JuZJc5B\nb2YD2vjx43n00Uf32vYGDRrEa6+9BsDXvvY1rr322r227b7yVTdmVrLa2nG0tJTvSpyamjo2bFhT\ntvXvicKreW655ZYK1qR4DnozK1kW8uW7EqelZd+9NHIgXoHkoRszG/CWLl3K0UcfTXV1NV/+8pfZ\nsWMHmzdv5owzzmD06NFUV1dzxhln0Nz8/hfp3n777UyYMIFhw4YxYcIEfv3rX3c8N2/ePI466iiq\nq6uZOnUq69at63K7F110EVdffTUAv/3tbzn00EO54YYbqKmpYcyYMdx+++0d8+7YsYPLL7+curo6\nDj74YC699FLefXfv/CSvg97MBrw777yThx9+mNWrV/Pqq6/ywx/+kIhgxowZrF+/nnXr1lFVVcVl\nl10GQFtbGzNnzmTRokVs3bqVp556ismTs6/0WrhwIXPmzOG+++5j48aNnHzyyZx77rlF1WPDhg1s\n27aNN954g9tuu42vf/3rbNmyBYArrriCVatW8fzzz7Nq1Sqam5u55ppryvOCdFbBr9qMlJEd11a6\nGmZ7rKv/x9n/7yjjX/HvnXHjxsWtt97a8fjBBx+MiRMnfmi+ZcuWxciRIyMi4s9//nOMGDEi7r33\n3njnnXc+MN/UqVNj3rx5HY937doVVVVVsW7duoiIkBSrV6+OiIgLL7wwvv/970dERFNTU1RVVcWu\nXbs6lh09enQsWbIkIiI+/vGPx2uvvdbx3FNPPRXjx4/vtl3dvQZ5eUl56z16Mxvwxo4d2zFdV1fH\nG2+8wV/+8hcuueQSxo0bx/Dhw/nsZz/L5s2biQiqqqq4++67ueWWWzj44IM544wzWLFiBQBr165l\n5syZjBw5kpEjR1JdXY2kDwz7dKe6uppBg96P1aqqKrZv387GjRtpa2vjuOOO61jv1KlTefvtt/v/\nxeiCg97MBrz169d3TK9du5ZDDjmEH//4x6xcuZKnn36azZs38/jjjwPvn0w99dRTWbx4MRs2bODI\nI4/k4osvBuDQQw/lZz/7Ga2trbS2trJp0ya2b9/OCSec0Of6jRo1iqqqKl566aWO9W7evLljWKfc\nHPRmNuDddNNNNDc309raynXXXce0adPYvn07Q4cOZdiwYbS2ttLQ0NAx/1tvvcX9999PW1sb+++/\nPwcccEDHnvhXv/pVrrvuOpYvXw7Ali1buOeee/aofpK4+OKL+da3vsXGjRsBaG5uZvHixXu03mI5\n6M2sZDU1dYDK9petvziSOO+885gyZQoTJ07kiCOO4Hvf+x4zZ86kra2NUaNGceKJJ3L66ad3LLN7\n925uuOEGxowZw6hRo3j88cc7rok/++yzmTVrFtOnT2f48OEcc8wxPPTQQx/YXil1azdnzhwmTpzI\nCSecwPDhw5kyZUrHcFG5+dsry8TfXmmp8LdXVoa/vdLMzIrmoDczS5yD3swscQ56M7PEOejNzBLn\noDczS5yD3swscQ56M7PEOejN7CPh9NNPZ/78+ZWuRkX4ztgy8Z2xloqu7tCsHVtLS3NL2bZZM6aG\nDa9v6PPyP/jBD1i9ejV33HFHP9Zq7+rPO2P9U4JmVrKW5hZoKOP6G8r3IfJR5KEbMxvQ5s6dy9ix\nYxk2bBiTJk3iwQcf5LrrruPuu+/mwAMP5NhjjwXglFNOYd68eQA0NjZy0kkn8e1vf5sRI0YwceJE\nfv/739PY2Mhhhx1GbW3tgD4a6MxBb2YD1ooVK7jpppt49tln2bp1K4sWLWLSpElcddVVTJs2jW3b\ntrFs2bIul126dCmTJ0+mtbWVc889l+nTp/PMM8+wevVq5s+fz2WXXUZbW9teblF5OOjNbMAaPHgw\nO3bs4MUXX2Tnzp0cdthhjB8/vqhlx48fz/nnn48kpk2bxuuvv87s2bPZf//9OfXUUxkyZAirVq0q\ncwv2Dge9mQ1YEyZM4MYbb6ShoYHRo0dz3nnn8eabbxa1bE1NTcf00KFDgeyXoArLtm/f3r8VrhAH\nvZkNaNOnT+eJJ55g3bp1AFxxxRUl/TjIR4GD3swGrBUrVvDYY4+xY8cOhgwZwtChQxk8eDC1tbWs\nWbOmpMubU74U2pdXmlnJasbUlPUSyJoxNb3PBLz77rvMmjWLV155hf33358TTzyRW2+9lSFDhjB/\n/nyqq6s5/PDDeeaZZ3rdy+/8fEpHBUXfMCVpEPAM8HpEnClpBHA3UAesAc6JiC35vFcCM4CdwMyI\n+NAv4PqGKbOBwT8lWBmV+inBmcDygsezgEci4kjgUeDKvBJHAecAk4CpwM1K6aPRzGyAKSroJY0F\nTgduKyg+C2jMpxuBs/PpM4G7ImJnRKwBVgLH90ttzcysZMXu0f8E+A5QeBxRExEtABGxARidl48B\n1hfM15yXmZlZBfR6MlbSF4CWiHhOUn0Ps5Y8iNfQ0NAxXV9fT319T6s3M/voaWpqoqmpaY/W0evJ\nWEnXAf+D7MTqUOBA4F+AvwbqI6JFUi3wWERMkjQLiIiYmy//EDA7IpZ0Wq9PxpoNAD4ZWxl79WRs\nRFwVEYdFxOHAdODRiPgH4AHgwny2C4CF+fT9wHRJQySNByYCS0uplJmZ9Z89uY5+DrBA0gxgLdmV\nNkTEckkLyK7QeQ+4NOldd7PE1dXVJXVN+UBRV1fXb+vyD4+UiYduzCoj9fdeua+jNzOzAchBb2aW\nOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZ\nJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRm\nZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9\nmVnieg16SR+TtETSMkkvSJqdl4+QtFjSq5IWSTqoYJkrJa2U9LKkKeVsgJmZ9UwR0ftMUlVEtEka\nDDwJfBP4IvB2RFwv6QpgRETMknQU8CvgU8BY4BHgiOi0IUmdi5IiCYCU22i2L0r9vSeJiFApyxQ1\ndBMRbfnkx4D9gADOAhrz8kbg7Hz6TOCuiNgZEWuAlcDxpVTKzMz6T1FBL2mQpGXABuDhiHgaqImI\nFoCI2ACMzmcfA6wvWLw5LzMzswoodo9+d0QcSzYUc7yko8n26j8wW39XzszM9tx+pcwcEVslNQGn\nAS2SaiKiRVIt8FY+WzNwaMFiY/OyD2loaOiYrq+vp76+vpTqmJklr6mpiaampj1aR68nYyWNAt6L\niC2ShgKLgDnAZ4HWiJjbzcnYT5MN2TyMT8aa2V6S+nuvLydji9mjPxholDSIbKjn7oh4UNIfgAWS\nZgBrgXMAImK5pAXAcuA94NKkE93MbB9X1OWVZdmw9+jNrAxSf++V7fJKMzMbuBz0ZmaJc9CbmSXO\nQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJ\nc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ\n4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9m\nlrheg17SWEmPSnpJ0guSvpmXj5C0WNKrkhZJOqhgmSslrZT0sqQp5WyAmZn1TBHR8wxSLVAbEc9J\nOgB4FjgLuAh4OyKul3QFMCIiZkk6CvgV8ClgLPAIcER02pCkzkVJkQRAym002xel/t6TRESolGV6\n3aOPiA0R8Vw+vR14mSzAzwIa89kagbPz6TOBuyJiZ0SsAVYCx5dSKTMz6z8ljdFLGgdMBv4A1ERE\nC2QfBsDofLYxwPqCxZrzMjMzq4Cigz4ftrkHmJnv2Xc+LkrzOMnMbIDbr5iZJO1HFvLzI2JhXtwi\nqSYiWvJx/Lfy8mbg0ILFx+ZlH9LQ0NAxXV9fT319fUmVNzNLXVNTE01NTXu0jl5PxgJIugP4U0R8\nu6BsLtAaEXO7ORn7abIhm4fxyVgz20tSf+/15WRsMVfdfAZ4HHiBbHgmgKuApcACsr33tcA5EbE5\nX+ZK4MvAe2RDPYu7WK+D3sz6XervvbIEfbk46M2sHFJ/75Xl8kozMxvYHPRmZolz0JuZJc5Bb2aW\nOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZ\nJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRm\nZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWuF6D\nXtIvJLVIer6gbISkxZJelbRI0kEFz10paaWklyVNKVfFzcysOMXs0f8S+HynslnAIxFxJPAocCWA\npKOAc4BJwFTgZknqv+qamVmpeg36iPgdsKlT8VlAYz7dCJydT58J3BUROyNiDbASOL5/qmpmZn3R\n1zH60RHRAhARG4DRefkYYH3BfM15mZmZVch+/bSe6MtCDQ0NHdP19fXU19f3U3XMzNLQ1NREU1PT\nHq1DEb1ntKQ64IGIOCZ//DJQHxEtkmqBxyJikqRZQETE3Hy+h4DZEbGki3VGMdseqNpPTaTcRrN9\nUervPUlEREnnPosdulH+1+5+4MJ8+gJgYUH5dElDJI0HJgJLS6mQmZn1r16HbiTdCdQD1ZLWAbOB\nOcBvJM0A1pJdaUNELJe0AFgOvAdcmvRuu5nZAFDU0E1ZNuyhGzMrg9Tfe+UcujEzswHKQW9mljgH\nvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXO\nQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJ\nc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJc9CbmSXOQW9mljgHvZlZ\n4soW9JJOk/SKpBWSrijXdszMrGdlCXpJg4B/Bj4PHA2cK+mT5diWVUZTU1Olq1BWbp+lpFx79McD\nKyNibUS8B9wFnFWmbVkFpB4Ubp+lpFxBPwZYX/D49bzMzMz2Mp+M3UP33nsvc+fOrXQ1zD5ynn/+\neb7yla9UuhoDgiKi/1cqnQA0RMRp+eNZQETE3IJ5+n/DZmYfARGhUuYvV9APBl4F/g54E1gKnBsR\nL/f7xszMrEf7lWOlEbFL0mXAYrLhoV845M3MKqMse/RmZrbv2CsnYyX9QlKLpOcLymZLel3Sv+d/\np+2NuvQ3SWMlPSrpJUkvSPpmXj5C0mJJr0paJOmgSte1L7po3zfy8lT672OSlkhalrdvdl6eSv91\n174k+g+y+3byNtyfP06i79rl7VtW0L6S+26v7NFLOgnYDtwREcfkZbOBbRFxQ9krUEaSaoHaiHhO\n0gHAs2T3DFwEvB0R1+d3Bo+IiFmVrGtf9NC+aSTQfwCSqiKiLT+39CTwTeCLJNB/0G37ppJO//0v\n4DhgWEScKWkuifQddNm+krNzr+zRR8TvgE1dPFXSmeN9UURsiIjn8untwMvAWLIwbMxnawTOrkwN\n90w37Wu/J2LA9x9ARLTlkx8jO28VJNJ/0G37IIH+kzQWOB24raA4mb7rpn1QYt9V+jr6yyQ9J+m2\ngX54BSBpHDAZ+ANQExEtkIUlMLpyNesfBe1bkhcl0X/th8bABuDhiHiahPqvm/ZBGv33E+A7vP/h\nBQn1HV23D0rsu0oG/c3A4RExmew/4IA+hMyHNe4BZuZ7vp07ZkCf9e6ifcn0X0TsjohjyY7Ejpd0\nNAn1XxftO4oE+k/SF4CW/Iizpz3cAdl3PbSv5L6rWNBHxMZ4/wTBz4FPVaoue0rSfmQhOD8iFubF\nLZJq8udrgbcqVb891VX7Uuq/dhGxFWgCTiOh/mtX2L5E+u8zwJmSXgN+DfytpPnAhkT6rqv23dGX\nvtubQS8KPpXyDmj398CLe7Eu/W0esDwiflpQdj9wYT59AbCw80IDyIfal0r/SRrVfugraShwKtl5\niCT6r5v2vZJC/0XEVRFxWEQcDkwHHo2IfwAeIIG+66Z95/el78pyw1Rnku4E6oFqSeuA2cApkiYD\nu4E1wCV7oy79TdJngC8BL+TjoAFcBcwFFkiaAawFzqlcLfuuh/adl0L/AQcDjcq+WnsQcHdEPCjp\nDyTQf3TfvjsS6b+uzCGNvuvO9aX2nW+YMjNLXKWvujEzszJz0JuZJc5Bb2aWOAe9mVniHPRmZolz\n0JuZJc5Bb2aWOAe9mVni/j/q4JWBq8Hb9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e9b8c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(baseline_rates, label='baseline');\n",
    "plt.hist(stim_rates, label='stim');\n",
    "plt.legend();\n",
    "plt.title('Single trial firing rates');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and use these to generate fake spike counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "baseline_counts = stats.poisson.rvs(baseline_rates * Tbaseline)\n",
    "stim_counts = stats.poisson.rvs(stim_rates * Tstim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHi5JREFUeJzt3XuUVeWd5vHvA4JNRVEomioFKRCiQXuMl4mxjcYyRiPa\nXnpWj6LpjkraMVE7ZIyJ4CShSCeImcTE1TFObEUJK95zAWeIoAvLayKakHhBQUAuFnIgIgoSQeE3\nf+xNeSyoy7lxqjbPZ62zOPv27vc9p3jOe96zL4oIzMwsu3pVuwJmZlZZDnozs4xz0JuZZZyD3sws\n4xz0ZmYZ56A3M8s4B73tkqQLJT1YprIekTSuHGW1U36X6yppkqQZZd7/CZJeypt+VdJnyrkPs1I4\n6PdgaUA9KWmDpL9IelzSMQARcWdEnF7l+jVI2i6pw7/TIupa1pNHIuKJiBhdzjKrRdLtkr5T7XpY\nee1V7QpYdUjaF3gAuAy4D+gLnAhsqWa92hBJKKvdFaTeEbFt91XJrOdxj37PdQgQEXFvJLZExMMR\n8QKApIskPb5j5bRnfZmkxZLWS/pJ3rJekn4oaZ2kpZKu6KgnLmmcpIWS3pD0W0nD2qnjo+m/GyS9\nLemTab2ekHSDpL8Ak3ZR1x9LWinpLUnPSDqhKy+IpFpJD0h6M63bo3nLXpU0QdKL6bLbJPVNl50k\naVU7ZY6WtEzS+en0AZLul7Q2fa3+rYP6/E36ui5P6/SYpL3TZWdLeiF9L+ZJ+ljedtslHZw33dpL\n31FXSVdJyklqkXRxuuxS4PPAN9LXe2Y6/xpJr6XzXpJ0cldeT+s+HPR7rsXANkl3SDpd0v67WKft\nEMeZwDHAx4HzJJ2Wzv8fwOeAI4CjgXN3sS0Aks4BJqTr/C3wOHBXO3X8dPpv/4joHxFPp9OfBJYA\ng4Hv7aKu89O6DADuBO7bEcqd+BqwCqhNy762zfILgVOBkcChwDfzlu3UXklHAw8CV0TEPZJE8i1q\nAXAAcAowXtKp7dTnh8BRwHHAQOAbwHZJh6Tt+grJa/hb4AFJO76hdzY0VQ/sCxwI/Ctwk6T9IuI/\ngV8A309f73PSfV0BHBMR/Une5+WdlG/djIN+DxURG4ETgO3ALcBaSTMl/W0Hm10XERsjYhXwCHBk\nOv+/AzdGxOsR8RYwtYMyLkvLWRwR29N1j5R0UAfbtB26aYmIn0bE9ojYaagpHbPfkC7/EbA3STB3\n5j2SAB4REdsi4sk2y/8jIlZHxAaSD5gLOijr08BM4J8j4rfpvE8AgyLie2n5y4FbgbFtN04/FC4B\nvhIRa9JvXb+PiPeA84D/GxHz0mGrHwD9gON3bN5JO7cC/57W4bfAJtp/fbaRDOv9naS9ImJlRLza\nSfnWzTjo92ARsSgixkXEMODvSHp4P+5gk1ze883APunzA0l6wjvschgj1QDcmA45rAfeIOmBDimg\n6h2Vj6Sr06GhNyW9CfQHBnWh3O8DS4G5kpZIuqbN8tfynq8gaXd7LgOejIjH8+Y1AEN2tD2t20SS\nbw9tDSL5gFq2i2UHpvsHkvE3ktekq6/hG+mH7A757+WHRMRS4KtAE5CTdKekA7q4H+smHPQGQEQs\nBu4gCfxCvQ4MzZtub8wdkkC6LCIGpo8BEbFPRPx+V9Vqr7rtFS7pRODrwD+lZQ8A3qbzXi4R8U5E\nXB0RI4GzgavajEfnf+toAFZ3UNyXgGGSbsibtwpY1qbt+0XEWbvY/i/AuyTDRG2tTvef7yA++CDa\nDNTkLavvoJ5t7fTaRsTdEXFi3j47+sZm3ZCDfg8l6dD0B7kh6fRBJEMRvyuiuHtJxpoPTMf6v9HB\nuv8HuFbSYel+95P0T+2su45kaGlXYdeefUiGYN6Q1FfSt0nGozsl6UxJO/a1EXifZOhihyskDZE0\nkGT8/u4OitsInA58WtJ16bz5wEZJ30h/aO0t6XBJ/7XtxmkvfRpwQ/oDbi9Jx0nqQ/J6nynpZEl7\nSbqa5ENhx3u3ALgw3eZ04KSutD+VA/J/yD0k3U9fkiGfv5K8J9aDOOj3XBtJftR8WtJG4CngOeDq\ndtZv29PLn/5PYG66/R+A/we8nzc80LpuRPyGpEd4t6QN6Ta7PAY+Iv5KMhb+ZDrUcWwX2jUnfSwG\nXiXp3XY41JPno8DD6evxJHBTRDyWt/zOtJ1LgFf44Ifgnaqe1v9tkh9vT5c0OX09/oHkt41XgbUk\nr13/dsq5GngeeIZkiGsq0Cv99vXPwE9IPgzPBM6KiPfT7b5K8o3kTZIP71930u789/I24PD09f4V\nyfj81HQ/q0l+/J3YSXnWzaizG49Iuo3kjzMXEUe0WfY14H+T/MC0Pp03ERhH0hsaHxFzK1Fx677S\nXuTNETGi2nUpF0mvAl+MiHnVrotZobrSo7+d5JCqD5E0lKS3siJv3miSIwJGA2OAn6ZHD1iGpcMQ\nY9KhiCHAJOBX1a6XmSU6DfqIeILkK2BbPyL50SvfOcDdEfF+eujYK0BXvm5bzyZgMrCeZOjmRZKw\nzxLfc9N6rKIugSDpbGBVRDzfpsM+hA//mNdCYYfNWQ+UjqVn+gM9Ig7ufC2z7qngoJfUj+SIg/bO\n5jMzs26kmB79SGA48Od0/H0o8Mf0iIgWPnwM9dB03k4k+auwmVkRIqKg3z67enil0gcR8UJE1EfE\nwelRFa8BR0XEWmAWcH56/PIIYBTJscPtVTazj0mTJlW9Dm6f27cnti/LbYsorn/cadBLupPkGOtD\nlFwR8JK2eZ33IbCQ5GSOhcBs4PIotmZmZlYWnQ7dRMSFnSw/uM30dcB17axuZma7mc+MrZDGxsZq\nV6Gi3L6eLcvty3LbitXpmbEV27HkUR0zswJJIgr8Mda3EjSzDg0fPpwVK1Z0vqKVVUNDA8uXLy9L\nWe7Rm1mH0h5ktauxx2nvdS+mR+8xejOzjHPQm5llnIPezCzjHPRWlPr64Ugq6FFfP7za1bYMGjFi\nBPPm7b7bBPTq1Ytly5Jb+X75y1/me99r7/4z3YePurGi5HIrKPTKvbmcb02QFfX1w9O/gcqoq2tg\nzZrlFSu/FPlX7L355purWJOuc9CbWcGK+aAvrPzu2ynoiUcgeejGzHq8+fPnc/jhh1NbW8sXv/hF\ntm7dyoYNGzjrrLMYPHgwtbW1nHXWWbS0fHAx3TvuuIORI0fSv39/Ro4cyV133dW6bNq0aRx22GHU\n1tYyZswYVq5cucv9XnLJJXz7298G4NFHH+Wggw7ihhtuoK6ujiFDhnDHHXe0rrt161auvvpqGhoa\nOOCAA7j88svZsmVLZV6QNhz0Ztbj3XnnnTz00EMsXbqURYsW8d3vfpeIYNy4caxatYqVK1dSU1PD\nlVdeCcDmzZsZP348c+bM4e233+app57iyCOPBGDmzJlMnTqV3/zmN6xbt44TTzyRCy64oEv1WLNm\nDRs3bmT16tXceuutXHHFFbz11lsAXHPNNSxZsoTnnnuOJUuW0NLSwne+853KvCBtVfFSm2E9FxAQ\nBT78nvdEu3rfinv/K/O3Mnz48Ljllltap2fPnh2jRo3aab0FCxbEwIEDIyLinXfeiQEDBsSvfvWr\n+Otf//qh9caMGRPTpk1rnd62bVvU1NTEypUrIyJCUixdujQiIi6++OL41re+FRERzc3NUVNTE9u2\nbWvddvDgwfH0009HRMRHPvKRWLZsWeuyp556KkaMGNFuu9p7DdL5BeWte/Rm1uMNHTq09XlDQwOr\nV6/m3Xff5bLLLmP48OHsv//+nHTSSWzYsIGIoKamhnvuuYebb76ZAw44gLPOOovFixcDsGLFCsaP\nH8/AgQMZOHAgtbW1SPrQsE97amtr6dXrg1itqalh06ZNrFu3js2bN3PMMce0ljtmzBjeeOON8r8Y\nu+CgN7Meb9WqVa3PV6xYwYEHHsgPfvADXnnlFZ555hk2bNjAY489BnzwY+qpp57K3LlzWbNmDYce\neiiXXnopAAcddBA/+9nPWL9+PevXr+fNN99k06ZNHHfccUXXb9CgQdTU1PDiiy+2lrthw4bWYZ1K\nc9CbWY9300030dLSwvr165kyZQrnn38+mzZtol+/fvTv35/169fT1NTUuv7atWuZNWsWmzdvpk+f\nPuyzzz6tPfEvfelLTJkyhYULFwLw1ltvcf/995dUP0lceumlfPWrX2XdunUAtLS0MHfu3JLK7SoH\nvZkVrK6ugQ/uMFr+R1J+10jiwgsv5LTTTmPUqFF89KMf5Zvf/Cbjx49n8+bNDBo0iOOPP54zzjij\ndZvt27dzww03MGTIEAYNGsRjjz3Wekz8ueeey4QJExg7diz7778/RxxxBA8++OCH9ldI3XaYOnUq\no0aN4rjjjmP//ffntNNOax0uqjRfvdKKkvwBF/r++SqIPZGvXlkdvnqlmZl1mYPezCzjHPRmZhnn\noDczyzgHvZlZxnUa9JJuk5ST9FzevO9LeknSnyT9UlL/vGUTJb2SLj+tUhU3M7Ou6UqP/nbgc23m\nzQUOj4gjgVeAiQCSDgPOA0YDY4CfqpCDTs3MrOw6DfqIeAJ4s828hyNiezr5e2DHhSbOBu6OiPcj\nYjnJh8Cx5auumZkVqhxj9OOA2enzIcCqvGUt6Twzs6o644wzmDFjRrWrURUl3WFK0v8C3ouIuzpd\neRfyrz3R2NhIY2NjKdUxs92kfmg9uZZcxcqvG1LHmtfWFL395MmTWbp0KT//+c9b582ePbuDLbqv\n5uZmmpubSyqj6KCXdDFwBvCZvNktwEF500PTebuUH/Rm1nPkWnLQVMHymyr3IdLTtO0ET548ueAy\nujp0s+NqQ8mEdDrwdeDsiMi/F9YsYKykvpJGAKOA+QXXysysi66//nqGDh1K//79GT16NLNnz2bK\nlCncc8897Lvvvhx11FEAnHzyyUybNg2A6dOnc8IJJ3DVVVcxYMAARo0axe9+9zumT5/OsGHDqK+v\n/9C3gZ6uK4dX3gk8BRwiaaWkS4D/APYBHpL0R0k/BYiIhcC9wEKScfvLfeUyM6uUxYsXc9NNN/GH\nP/yBt99+mzlz5jB69GiuvfZazj//fDZu3MiCBQt2ue38+fM58sgjWb9+PRdccAFjx47l2WefZenS\npcyYMYMrr7ySzZs37+YWVUanQzcRceEuZt/ewfrXAdeVUikzs67o3bs3W7du5YUXXqC2tpZhw4Z1\nedsRI0bwhS98AYDzzz+fKVOmMGnSJPr06cOpp55K3759WbJkCUcccUSlqr/b+MxYM+uxRo4cyY9/\n/GOampoYPHgwF154Ia+//nqXtq2rq2t93q9fPyC5E1T+vE2bNpW3wlXioDezHm3s2LE8/vjjrFy5\nEoBrrrmmoJuD7Akc9GbWYy1evJhHHnmErVu30rdvX/r160fv3r2pr69n+fLlBd0wJcs/J5Z0HL2Z\n7ZnqhtRV9BDIuiF1na8EbNmyhQkTJvDyyy/Tp08fjj/+eG655Rb69u3LjBkzqK2t5eCDD+bZZ5/t\ntJffdnmWvhX4VoJWFN9KcM/hWwlWh28laGZmXeagNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjPNx\n9GbWoYaGhkwdU95TNDQ0lK0sH0dvRfFx9GbV4ePozcxsJw56M7OMc9CbGZDcB1ZSWR/1Q+ur3SzD\nY/RWJI/RZ4+k8t8HtinbV4WsBo/Rm5nZThz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcZ0GvaTb\nJOUkPZc3b4CkuZIWSZojab+8ZRMlvSLpJUmnVariZmbWNV3p0d8OfK7NvAnAwxFxKDAPmAgg6TDg\nPGA0MAb4qXw1JDOzquo06CPiCeDNNrPPAaanz6cD56bPzwbujoj3I2I58ApwbHmqamZmxSh2jH5w\nROQAImINMDidPwRYlbdeSzrPzMyqpFzXoy/qHOempqbW542NjTQ2NpapOmZm2dDc3Exzc3NJZXTp\nWjeSGoAHIuKIdPoloDEicpLqgUciYrSkCUBExPXpeg8CkyLi6V2U6Wvd9GC+1k32+Fo3PUMlr3Wj\n9LHDLODi9PlFwMy8+WMl9ZU0AhgFzC+kQmZmVl6dDt1IuhNoBGolrQQmAVOB+ySNA1aQHGlDRCyU\ndC+wEHgPuNzddjOz6vJliq0oHrrJHg/d9Ay+TLGZme3EQW9mlnEOerMK8y36rNrKdRy9mbUj15Ir\n+9h3rilX3gIt09yjNzPLOAe9mVnGOejNzDLOQW9mlnEOeqO+fnjBR32YWc/ho26MXG4FxZzlamY9\ng3v0ZmYZ5x69WU/UGw+hWZc56M16om1U5AJklk0eujEzyzgHvZlZxjnozaxy0t8SfFG36vIYvZlV\nTiV+S8AXdSuUe/RmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDnozs4wrKegl/U9JL0h6TtIvJPWV\nNEDSXEmLJM2RtF+5KmtmZoUrOuglHQj8G3B0RBxBckz+BcAE4OGIOBSYB0wsR0XNzKw4pQ7d9AY+\nImkvoB/QApwDTE+XTwfOLXEfZmZWgqKDPiJWAz8EVpIE/FsR8TBQFxG5dJ01wOByVNTMzIpT9CUQ\nJO1P0ntvAN4C7pP0eXa+VVG7ty5qampqfd7Y2EhjY2Ox1TEzy6Tm5maam5tLKqOUa918FlgWEesB\nJP0aOB7ISaqLiJykemBtewXkB72Zme2sbSd48uTJBZdRyhj9SuA4SX+j5FY3pwALgVnAxek6FwEz\nS9iHmZmVqOgefUTMl3Q/sAB4L/33FmBf4F5J44AVwHnlqKiZmRWnpMsUR8RkoO33iPUkwzpmZtYN\n+MxYM7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPe\nzCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDvoMqa8fjqSCH2aWbSXdYcq6l1xuBRBFbOmw\nN8sy9+jNzDLOQW9mlnEOejOzjHPQm5llnIPezCzjSgp6SftJuk/SS5JelPRJSQMkzZW0SNIcSfuV\nq7JmZla4Unv0NwKzI2I08HHgZWAC8HBEHArMAyaWuA8zMytB0UEvqT9wYkTcDhAR70fEW8A5wPR0\ntenAuSXX0szMilZKj34E8BdJt0v6o6RbJNUAdRGRA4iINcDgclTUzMyKU8qZsXsBRwNXRMSzkn5E\nMmzT9tTMdk/VbGpqan3e2NhIY2NjCdUxM8ue5uZmmpubSyqjlKB/DVgVEc+m078kCfqcpLqIyEmq\nB9a2V0B+0JuZ2c7adoInT55ccBlFD92kwzOrJB2SzjoFeBGYBVyczrsImFnsPszMrHSlXtTsK8Av\nJPUBlgGXAL2BeyWNA1YA55W4DzMzK0FJQR8RfwY+sYtFny2lXDMzKx+fGWtmlnEOejOzjHPQ2260\nd1F3wKqvH17tipv1aL7DVDdVXz88vWNUlmyhmDtg5XK+A5ZZKRz03VRxtwV0IJrZzjx0Yz1A4UM+\nHu4x+4B79NYDFD7k4+Eesw+4R29mlnEOejOzjHPQm5llnIPezCzjHPRmZhnnoDczyzgHvVmqfmh9\nUZdo6OxhVm0+jt4slWvJQVMFCq5EmWYFcI/ezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZx\nDnozs4xz0JuZZVzJQS+pl6Q/SpqVTg+QNFfSIklzJO1XejXNzKxY5ejRjwcW5k1PAB6OiEOBecDE\nMuzDzMyKVFLQSxoKnAHcmjf7HGB6+nw6cG4p+zAzs9KU2qP/EfB1PnxDz7qIyAFExBpgcIn7MDOz\nEhR9UTNJZwK5iPiTpMYOVm33rs5NTU2tzxsbG2ls7KgYM7M9T3NzM83NzSWVUcrVKz8FnC3pDKAf\nsK+kGcAaSXURkZNUD6xtr4D8oDczs5217QRPnjy54DKKHrqJiGsjYlhEHAyMBeZFxL8ADwAXp6td\nBMwsdh9mZla6ShxHPxU4VdIi4JR02szMqqQsNx6JiEeBR9Pn64HPlqNcMzMrnc+MNTPLOAe9mfU8\nvSn7vX3rh9ZXu1UV43vGmlnPs42y34s315Qrb4HdiHv0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aW\ncQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPe\nzCzjHPRmZhnnoDczyzgHvZlZxjnozcwyruiglzRU0jxJL0p6XtJX0vkDJM2VtEjSHEn7la+6ZmZW\nqFJ69O8DV0XE4cDfA1dI+hgwAXg4Ig4F5gETS6+mmZkVq+igj4g1EfGn9Pkm4CVgKHAOMD1dbTpw\nbqmVNDOz4pVljF7ScOBI4PdAXUTkIPkwAAaXYx9mZlacvUotQNI+wP3A+IjYJCnarNJ2ulVTU1Pr\n88bGRhobG0utju0petfDtlyHq0jaTZUxq5zm5maam5tLKqOkoJe0F0nIz4iImensnKS6iMhJqgfW\ntrd9ftCbFWRbDprKXGa5yzMrg7ad4MmTJxdcRqlDN9OAhRFxY968WcDF6fOLgJltNzIzs92n6B69\npE8Bnweel7SAZIjmWuB64F5J44AVwHnlqKiZmRWn6KCPiCeB3u0s/myx5ZqZWXn5zFgzs4xz0JuZ\nZZyD3sws4xz0ZmYZ56CvsPr64Ugq+GFmVi4lnxlrHcvlVtDBycEdcNibWXm4R29mlnEOejOzjHPQ\nm5llnIPezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZgbQm6LOYu/oUT+0vtqtAnxmbJfdffe9\njB9/bbWrYWaVso2y304y19TxfY13Fwd9Fz3zzB9Yu/YfgcsK2Go5cGplKmRm1kUO+oIMAkYVsH4x\n17gxMysvj9GbmWWcg97MLOMc9GZmGeegNzPLuIoFvaTTJb0sabGkayq1HzOzbqsCx+YXoyJH3Ujq\nBfwEOAVYDTwjaWZEvFyJ/XVPzcCQaleiG3gU+vw3PjgCaWCR5bTZ7r0SqtQVrwIjKryPaspy+7pT\n2ypwbH4x5VXq8MpjgVciYgWApLuBc4A9LOg/X+1KdAMvwkffgc9tSaffLLKcvO3WAHeVWK3OLKf7\nhEUlLCe77VtOdttWpEoF/RBgVd70ayThb3uivoL9yljexjKWZbYH8AlTXdSnTx/23vsO9t77iS6t\n/+67i+jb91E2bapwxbq9PvAKMK0/8DbQv4gy2my3ZRvwTjkqZ7ZHUET5z96UdBzQFBGnp9MTgIiI\n6/PW8WmjZmZFiIiCfpWtVND3BhaR/Bj7OjAfuCAiXir7zszMrEMVGbqJiG2SrgTmkhzCeZtD3sys\nOirSozczs+6jKmfGZu1kKkm3ScpJei5v3gBJcyUtkjRHUjmPO9ltJA2VNE/Si5Kel/SVdH5W2re3\npKclLUjbNymdn4n27SCpl6Q/SpqVTmemfZKWS/pz+h7OT+dlqX37SbpP0kvp/8NPFtq+3R70eSdT\nfQ44HLhA0sd2dz3K7HaS9uSbADwcEYcC84CJu71W5fE+cFVEHA78PXBF+n5lon0RsQU4OSKOAo4E\nxkg6loy0L894YGHedJbatx1ojIijImLHYdxZat+NwOyIGA18nOR8pMLaFxG79QEcB/w2b3oCcM3u\nrkcF2tUAPJc3/TJQlz6vB16udh3L1M7fAJ/NYvuAGuBZ4BNZah8wFHgIaARmpfOy1L5Xgdo28zLR\nPpLjipfuYn5B7avG0M2uTqbK4rUCBkdEDiAi1gCDq1yfkkkaTtLr/T3JH1km2pcOaywgOef2oYh4\nhgy1D/gR8HU+fCecLLUvgIckPSPpX9N5WWnfCOAvkm5Ph95ukVRDge3z1St3nx79q7ekfYD7gfER\nsYmd29Nj2xcR2yMZuhkKHCvpcDLSPklnArmI+BPQ0bHXPbJ9qU9FxNHAGSRDiyeSkfeP5MjIo4Gb\n0ja+QzIKUlD7qhH0LcCwvOmh6bysyUmqA5BUD6ytcn2KJmkvkpCfEREz09mZad8OEfE2yUWKTic7\n7fsUcLakZSRXCPqMpBnAmoy0j4h4Pf13HcnQ4rFk5/17DVgVEc+m078kCf6C2leNoH8GGCWpQVJf\nYCwwqwr1KDfx4R7TLODi9PlFwMy2G/Qg04CFEXFj3rxMtE/SoB1HLEjqR3I395fISPsi4tqIGBYR\nB5P8X5sXEf8CPEAG2iepJv22iaSPAKcBz5Od9y8HrJJ0SDrrFOBFCmxfVY6jl3Q6yS/JO06mmrrb\nK1FGku4k+aGrFsgBk0h6FvcBBwErgPMiYkO16lgsSZ8CHiP5zxPp41qSs53vpee3778A00n+FnsB\n90TE9yQNJAPtyyfpJOBrEXF2VtonaQTwa5K/y72AX0TE1Ky0D0DSx4FbgT7AMuASoDcFtM8nTJmZ\nZZx/jDUzyzgHvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ9/8BC/9HZhVhUd0A\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e9c7390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(baseline_counts, label='baseline');\n",
    "plt.hist(stim_counts, label='stim');\n",
    "plt.legend();\n",
    "plt.title('Single trial spike counts');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate firing rates from data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "obs_baseline_rates = baseline_counts / Tbaseline\n",
    "obs_stim_rates = stim_counts / Tstim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-62.191816367 3.41747897325e-237\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "tstat, pval = stats.ttest_rel(obs_baseline_rates, obs_stim_rates)\n",
    "print(tstat, pval)\n",
    "print(pval < alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now put it in a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_expt(baseline, effect, Tbaseline, Tstim, trial_noise_std, Ntrials):\n",
    "    # run a single synthetic experiment, return p-value for paired t-test\n",
    "    # of baseline vs stim spike rate\n",
    "    \n",
    "    # random rates\n",
    "    baseline_rates = np.exp(np.log(baseline) * (1 + trial_noise_std * np.random.randn(Ntrials)))\n",
    "    stim_rates = np.exp(np.log(effect * baseline) * (1 + trial_noise_std * np.random.randn(Ntrials)))\n",
    "    \n",
    "    # random counts\n",
    "    baseline_counts = stats.poisson.rvs(baseline_rates * Tbaseline)\n",
    "    stim_counts = stats.poisson.rvs(stim_rates * Tstim)\n",
    "    \n",
    "    # observed rates\n",
    "    obs_baseline_rates = baseline_counts / Tbaseline\n",
    "    obs_stim_rates = stim_counts / Tstim\n",
    "    \n",
    "    # test\n",
    "    return stats.ttest_rel(obs_baseline_rates, obs_stim_rates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_power(alpha, Nexpt, exptfun, *args):\n",
    "    # calculate power for an experiment defined by false positive rate alpha and repeated Nexpt times\n",
    "    # exptfun returns an object with a pvalue\n",
    "    # args are passed on to exptfun\n",
    "    \n",
    "    pvals = np.array([exptfun(*args).pvalue for _ in range(Nexpt)])\n",
    "    \n",
    "    return np.mean(pvals < alpha)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now calculate power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.68400000000000005"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = 0.05\n",
    "Nexpt = 1000\n",
    "baseline = 20.  # firing rate in Hz\n",
    "Tbaseline = 1.\n",
    "Tstim = 1.\n",
    "effect = 1.2  # fractional change from baseline\n",
    "trial_noise_std = 0.10  # trial-to-trial noise as a percentage\n",
    "Ntrials = 50  # number of trials for each stimulus\n",
    "\n",
    "calc_power(alpha, Nexpt, run_expt, baseline, effect, Tbaseline, Tstim, trial_noise_std, Ntrials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}