bayesian_statistics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian statistics\n",
    "\n",
    "Bayesian statistics allows priors (null hypotheses) and posteriors (effect sizes) to take on *distributions*. This allows us to assign degrees of certainty to events and statements. \n",
    "\n",
    "Typically, a harsh contrast is drawn between frequentists (null hypothesis, p-values) and Bayesians, and this is probably because they address different *questions*.\n",
    "\n",
    "The frequentist asks: \n",
    "* do I have sufficient data to reject the null hypothesis -- that there is no difference between my observed data and a straw-man (null) hypothesis?\n",
    "\n",
    "The Bayesian asks:\n",
    "* what can I learn about my data and the underlying system that was used to generate it?\n",
    "\n",
    "The divide arises from different ways of seein the natural world. A frequentist will assume that *there is a true value of a system* and the data we use to estimate the system are noisy. A Bayesian would say that we would do better to think of the data as *fixed*, and assume that the *system is noisy*. \n",
    "\n",
    "This is best illustrated with an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This function is a t-test-like operation using Bayesian statistics\n",
    "def bayes_ttest(groups=None, N=40, show=False):\n",
    "  \"\"\"\n",
    "  Run a Bayesian t-test on sample or true data.\n",
    "  \"\"\"\n",
    "  if groups is None: # Generate some data\n",
    "    group1, group2 = gen_data(N=40)\n",
    "  elif len(groups) != 2:\n",
    "    print('T-test requires only 2 groups, not %i' %len(groups))\n",
    "    return None\n",
    "  else:\n",
    "    group1, group2 = groups\n",
    "  \n",
    "  pooled = np.concatenate((group1, group2)) # Pooled data\n",
    "  # Establish priors\n",
    "  mu1 = pm.Normal(\"mu_1\", mu=pooled.mean(), tau=1.0/pooled.var()/N)\n",
    "  mu2 = pm.Normal(\"mu_2\", mu=pooled.mean(), tau=1.0/pooled.var()/N)\n",
    "  sig1 = pm.Uniform(\"sigma_1\",lower=pooled.var()/1000.0,upper=pooled.var()*1000)  \n",
    "  sig2 = pm.Uniform(\"sigma_2\",lower=pooled.var()/1000.0,upper=pooled.var()*1000)\n",
    "  v = pm.Exponential(\"nu\", beta=1.0/29)\n",
    "  \n",
    "  # Set up posterior distribution\n",
    "  t1 = pm.NoncentralT(\"t_1\", mu=mu1, lam=1.0/sig1, nu=v, value=group1,\n",
    "                      observed=True)\n",
    "  t2 = pm.NoncentralT(\"t_1\", mu=mu2, lam=1.0/sig2, nu=v, value=group2,\n",
    "                      observed=True)\n",
    "  \n",
    "  # Generate the model\n",
    "  model = pm.Model( [t1, mu1, sig1, t2, mu2, sig2, v] ) # Push priors\n",
    "  mcmc = pm.MCMC(model) # Generate MCMC object\n",
    "  mcmc.sample(40000, 10000, 2) # Run MCMC sampler # \"trace\"\n",
    "  \n",
    "  # Get the numerical results\n",
    "  mus1 = mcmc.trace('mu_1')[:]  \n",
    "  mus2 = mcmc.trace('mu_2')[:]  \n",
    "  sigmas1 = mcmc.trace('sigma_1')[:]  \n",
    "  sigmas2 = mcmc.trace('sigma_2')[:]  \n",
    "  nus = mcmc.trace('nu')[:] \n",
    "  diff_mus = mus1-mus2  # Difference in mus\n",
    "  diff_sigmas = sigmas1-sigmas2  \n",
    "  normality = np.log(nus)  \n",
    "  effect_size = (mus1-mus2)/np.sqrt((sigmas1**2+sigmas2**2)/2.)  \n",
    "  print('\\n   Group 1 mu: %.4f\\n   Group 2 mu: %.4f\\n   Effect size: %.4f'\n",
    "        %(mus1.mean(), mus2.mean(), effect_size.mean()))\n",
    "  \n",
    "  if show: # Plot some basic metrics if desired\n",
    "    from pymc.Matplot import plot as mcplot\n",
    "    # mcplot(mcmc) # This plots 5 graphs, only useful as a benchmark.\n",
    "    \n",
    "    # Finally, what can this tell is about the null hypothesis?\n",
    "    # Split distribution\n",
    "    fig2 = plt.figure() \n",
    "    ax2 = fig2.add_subplot(121)\n",
    "    minx = min(min(mus1),min(mus2))  \n",
    "    maxx = max(max(mus1),max(mus2))  \n",
    "    xs = np.linspace(minx,maxx,1000)\n",
    "    gkde1 = stats.gaussian_kde(mus1)  \n",
    "    gkde2 = stats.gaussian_kde(mus2)\n",
    "    ax2.plot(xs,gkde1(xs),label='$\\mu_1$')  \n",
    "    ax2.plot(xs,gkde2(xs),label='$\\mu_2$')  \n",
    "    ax2.set_title('$\\mu_1$ vs $\\mu_2$')\n",
    "    ax2.legend()\n",
    "    \n",
    "    # Difference of mus\n",
    "    ax3 = fig2.add_subplot(122)\n",
    "    minx = min(diff_mus)  \n",
    "    maxx = max(diff_mus)  \n",
    "    xs = np.linspace(minx,maxx,1000)  \n",
    "    gkde = stats.gaussian_kde(diff_mus)\n",
    "    ax3.plot(xs,gkde(xs),label='$\\mu_1-\\mu_2$')  \n",
    "    ax3.legend()\n",
    "    ax3.axvline(0, color='#000000',alpha=0.3,linestyle='--')\n",
    "    ax3.set_title('$\\mu_1-\\mu_2$')\n",
    "    plt.show()\n",
    "  \n",
    "  return\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Actually run it with our test groups, A and B\n",
    "bayes_ttest(groups=[A,B], show=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can see that there is a difference in the means of A and B, but unlike with a simple t-test we get an estimate of the confidence of this effect size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Behind the scenes\n",
    "\n",
    "Bayesian statistics are most commonly run using Markov Chain Monte Carlo simulation methods. \n",
    "\n",
    "MCMC takes a data set and treats it as the operation to act on a prior distribution. MCMC is interested in re-creating the parameters used to generate the dataset. For a normal distribution (as above), this means a $\\mu$ (sample mean) and $\\sigma$ (standard deviation/variance). \n",
    "\n",
    "MCMC (hopefully) converges on parameters by repeatedly applying a mathematically complex but heuristic algorithm. In order to arrive at the posterior distribution, mu and sigma will take on values that approximate that of the data.\n",
    "\n",
    "Prior distributions of $\\mu$ and $\\sigma$ are fed into the MCMC. If we have no reason to suspect a bias, these can be uninformative priors, or uniform, spanning some range that will include the mean of the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(range(80), [1./80. for i in range(80)], lw=2.)\n",
    "plt.fill_between(range(80), [1./80. for i in range(80)], np.zeros(80), \n",
    "                 facecolor='blue', alpha=0.4)\n",
    "plt.ylim([0,0.25])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whereas for $\\sigma$ it makes less sense to choose a uniform prior (though we certainly could). We might instead be inclined to believe $\\sigma$ is closer to 0 than to 80, and so we make its prior an exponential distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(0., 30., 300)\n",
    "E = [np.exp(-(i/5.))/5. for i in x]\n",
    "plt.plot(x, E, linewidth=2.)\n",
    "plt.fill_between(x, E, np.zeros(len(E)), facecolor='b', alpha=0.4)\n",
    "plt.ylim([0,0.25])\n",
    "plt.xlim([0,80])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These distributions are approximated with randomly chosen numbers from the generator functions (beta, binomial, Bernoulli, Poisson, exponential) in a method known as Monte Carlo sampling.\n",
    "\n",
    "Here is a sample random exponential distribution that "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "E_data = np.random.exponential(5, 200)\n",
    "plt.hist(E_data, bins=20, facecolor='b', edgecolor='none', \n",
    "             alpha=0.4, normed=1)\n",
    "plt.ylim([0,0.25])\n",
    "plt.xlim([0,80])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The job of the MCMC sampler is to use the given data to transform the prior, all the while tracking how the parameters ($\\mu$ and $\\sigma$ here) respond to shifts in the posterior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, MCMC flows:\n",
    "\n",
    "1. Prior distribution is subjected to data operations:\n",
    " * a random data point is chosen and can move one unit direction\n",
    " * the random data point moves probabilistically based on which direction has more observed data points\n",
    " * the random data point can move or stay the same, and has no memory of where it's been (Markovian)\n",
    "2. This is repeated several thousand times, and each time creates a distribution\n",
    "\n",
    "At the end, we have a posterior distribution. This distribution is much  more informative than frequentist statistical tests. We can ask more detailed questions about the data, such as\n",
    "\n",
    "* what is the probability of the data having a mean of $\\mu_x$?\n",
    "* what is the likelihood of the data being greater than $x_0%?\n",
    "* what are the credible regions?\n",
    "* what are the credible regions of the effect sizes?\n",
    "\n",
    "Most valuably, we have probabilities for all of these questions and numbers, not just binary yes/no."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## P-substitute\n",
    "\n",
    "We also have a p-value substitute ready at hand. A simple alternative that is absolutely viable is the percent of samples that fall beyond a benchmark. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Simple p-substitute values\n",
    "locs = [1.5, 5.5, 11.5]\n",
    "scales = [1., 3., 5.]\n",
    "fig = plt.figure()\n",
    "for s in range(3):\n",
    "    ax = fig.add_subplot(1,3,s+1)\n",
    "    x = np.random.normal(loc=locs[s], scale=scales[s], size=400)\n",
    "    cnts, bins, patches = ax.hist(x, bins=20, facecolor='gray', \n",
    "                                  edgecolor='white', alpha=0.4)\n",
    "    bincents = 0.5 * (bins[:-1] + bins[1:]) # Bin centers\n",
    "    cols, cnts_sofar = [], 0. # Colors\n",
    "    for b in bincents:\n",
    "        if b < 0.:\n",
    "            cols.append('red')\n",
    "        else:\n",
    "            cols.append('blue')\n",
    "            cnts_sofar += cnts[list(bincents).index(b)]\n",
    "    # Set the facecolors\n",
    "    for c, p in zip(cols, patches): \n",
    "        plt.setp(p, 'facecolor', c)\n",
    "    # Mark the benchmark\n",
    "    ax.plot([0,0], [0,max(cnts)], '--', color='r', linewidth=2.)\n",
    "    ax.set_title('%.2f %% above 0' %float(cnts_sofar/sum(cnts)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sources\n",
    "\n",
    "** Chan, C** *Data analysis with Python -- Practical computing for biologists.* http://people.duke.edu/~ccc14/pcfb/analysis.html\n",
    "\n",
    "**Manly, B** *Randomization, Bootstrap and Monte Carlo methods in Biology.*\n",
    "\n",
    "**Pilon, CD** *Probabilistic programming and Bayesian methods for Hackers.* http://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/tree/master/\n",
    "\n",
    "**Vanderplas, J** *Frequentism and Bayesianism: A Python-driven Primer.* http://jakevdp.github.io/blog/2014/06/14/frequentism-and-bayesianism-4-bayesian-in-python/\n",
    "(Also an associated paper on arXiv: http://arxiv.org/abs/1411.5018)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}