01-Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Methods in Astronomy: Hands-on Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Big Picture: What Are We After?\n",
    "\n",
    "There are two fundamental types of statistical questions we'll want to answer:\n",
    "\n",
    "#### 1. Model Fitting\n",
    "*Given this Model, what parameters best fit my data?*\n",
    "\n",
    "Examples:\n",
    "\n",
    "- What are the slope and intercept of a line of best-fit?\n",
    "- What are the parameters of the best quadratic fit?\n",
    "- What is the frequency, amplitude, and phase of a sinusoidal fit?\n",
    "- What are the orbital parameters of a planet in this radial velocity data?\n",
    "\n",
    "#### 2. Model Selection\n",
    "\n",
    "*Given two potential Models, which better describes my data?*\n",
    "\n",
    "Examples:\n",
    "\n",
    "- Is there a linear trend in this data?\n",
    "- Does a linear or quadratic fit describe our data better?\n",
    "- Is there a periodic signal in this timeseries?\n",
    "- Does this star have a planet around it? Does this star have two planets around it?\n",
    "\n",
    "Often one of the two models is a *null hypothesis*, or a baseline model in which the effect you're interested in is not observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequentist vs Bayesian Approaches\n",
    "\n",
    "Both the model fitting and model selection problems can be approached from either a *frequentist* or a *Bayesian* standpoint.\n",
    "Fundamentally, the difference between these lies in the **definition of probability** that they use:\n",
    "\n",
    "- **A frequentist probability is a measure *long-run frequency* of (real or imagined) repeated trials.** Among other things, this generally means that observed data can be described probabilistically (you can repeat the experiment and get a different realization) while model parameters are fixed, and cannot be described probabilistically (the universe remains the same no matter how many times you observe it). \n",
    "\n",
    "- **A Bayesian probability is a *quantification of belief*.** Among other things, this generally means that observed data are treated as fixed (you know exactly what values you measured) while model parameters – including the \"true\" values of the data reflected by noisy measurements – are treated probabilistically (your degree of knowledge about the universe changes as you gather more data).\n",
    "\n",
    "Perhaps surprisingly to the uninitiated, scientists and statisticians have been vehemently fighting in favor of one approach or the other for decades.\n",
    "For more a more fleshed-out discussion of these different definitions and their consequences, you can see my [series of blog posts](http://jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/) on the topic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Bayesian Problem Setting\n",
    "\n",
    "Thus the end-goal of a Bayesian analysis is a probabilistic statement about the universe.\n",
    "Roughly we want to measure\n",
    "\n",
    "$$\n",
    "P(science)\n",
    "$$\n",
    "\n",
    "Where \"science\" might be encapsulated in the cosmological model, the mass of a planet around a star, or whatever else we're interested in learning about."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't of course measure this without reference to data, so more specifically we want to measure\n",
    "\n",
    "$$\n",
    "P(science~|~data)\n",
    "$$\n",
    "\n",
    "which should be read \"the probability of the science *given* the data.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we should be explicit that this measurement is not done in a vaccum: generally before observing any data we have *some* degree of background information that informs the science, so we should actually write\n",
    "\n",
    "$$\n",
    "P(science~|~data, background\\ info)\n",
    "$$\n",
    "\n",
    "This should be read \"the probability of the science given the data *and* the background information\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, there are often things in the scientific model that we don't particularly care about: these are known as \"nuisance parameters\". As an example of a nuisance parameter, if you are finding a planet in radial velocity data, the secular motion of the star is *extremely* important to model correctly, but in the end you don't really care about its best-fit value.\n",
    "\n",
    "With that in mind, we can write:\n",
    "\n",
    "$$\n",
    "P(science,nuisance\\ parameters~|~data, background\\ info)\n",
    "$$\n",
    "\n",
    "Where as before the comma should be read as an \"and\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is starting to get a bit cumbersome, so let's create some symbols that will let us express this more easily:\n",
    "\n",
    "$$\n",
    "P(\\theta_S, \\theta_N~|~D, I)\n",
    "$$\n",
    "\n",
    "- $\\theta_S$ represents the \"science\": the set of parameters that we are interested in constraining\n",
    "- $\\theta_N$ represents the \"nuisance parameters\": the set of parameters that are important in the model, but are not particularly interesting for the scientific result.\n",
    "- $D$ represents the \"observed data\"\n",
    "- $I$ represents the information or knowledge you had before observing the data, including whatever made you choose the model you're fitting.\n",
    "\n",
    "Finally, we'll often just write $\\theta = (\\theta_S, \\theta_N)$ as a shorthand for all the model parameters.\n",
    "\n",
    "This quantity, $P(\\theta~|~D,I)$ is called the \"posterior probability\" and determining this quantity is the ultimate goal of a Bayesian analysis.\n",
    "\n",
    "Now all we need to do is compute it!\n",
    "\n",
    "The core problem is this: **We do not have a way to directly calculate** $P(\\theta~|~D,I)$. We often do have an expression for $P(D~|~\\theta,I)$, but these two expressions are **not** equal.\n",
    "\n",
    "$$\n",
    "P(\\theta~|~D,I) \\ne P(D~|~\\theta,I)\n",
    "$$\n",
    "\n",
    "\n",
    "We need to figure out how these two expressions are related."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Nitty Gritty: Thinking about Probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization of Probability\n",
    "\n",
    "We will generally be dealing with probability *densities*, that is, $P(A) dA$ is the probability of a value falling between $A$ and $A + dA$.\n",
    "\n",
    "Probability densities are normalized such that the union of all possible events has a probability of unity; mathematically that criterion looks like this:\n",
    "\n",
    "$$\n",
    "\\int P(A) dA = 1\n",
    "$$\n",
    "\n",
    "Among other things, consider the units implied by this expression: because probability is dimensionless, the units of $P(A)$ must be the inverse of the units of $A$.\n",
    "This can be very useful to keep in mind as you manipulate probabilistic expressions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint Probability\n",
    "\n",
    "When we're talking about two different variables, we can write a joint probability:\n",
    "\n",
    "$$\n",
    "P(A,B)\n",
    "$$\n",
    "\n",
    "This should be read \"probability of $A$ and $B$\".\n",
    "\n",
    "In the case that $A$ and $B$ are independent, $P(A,B) = P(A)P(B)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conditional probability\n",
    "\n",
    "Sometimes we want to express the probability of one variable conditioned on the value of another variable. We can write it like this:\n",
    "\n",
    "$$\n",
    "P(A\\mid B) = \\frac{P(A, B)}{P(B)}\n",
    "$$\n",
    "\n",
    "This should be read \"the probability of $A$ given $B$\".\n",
    "\n",
    "Rearranging this we get the more common expression of this:\n",
    "\n",
    "$$\n",
    "P(A, B) = P(A\\mid B)P(B)\n",
    "$$\n",
    "\n",
    "Note that conditional probabilities are still probabilities, so that the normalization is the same as above:\n",
    "\n",
    "$$\n",
    "\\int P(A\\mid B)dA = 1\n",
    "$$\n",
    "\n",
    "This makes clear that $P(A \\mid B)$ has units equivalent to $1/A$ – and no units related to $B$. Among other things, this means that expressions like $\\int P(A\\mid B)dB$ should raise a very big red flag!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marginalization\n",
    "\n",
    "Given the above two properties, we can quite easily show that\n",
    "$$\n",
    "\\int P(A, B) dA = P(B)\n",
    "$$\n",
    "This is known as *marginalization* – we have integrated *A* out of the joint probability, and we are left with the raw probability of *B*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayes' Rule\n",
    "\n",
    "The definition of conditional probability is entirely symmetric, so we can write\n",
    "\n",
    "$$\n",
    "P(A, B) = P(B, A)\n",
    "$$\n",
    "\n",
    "$$\n",
    "P(A\\mid B)P(B) = P(B\\mid A)P(A)\n",
    "$$\n",
    "\n",
    "which is more commonly rearranged in this form:\n",
    "\n",
    "$$\n",
    "P(A\\mid B) = \\frac{P(B\\mid A)P(A)}{P(B)}\n",
    "$$\n",
    "\n",
    "This is known as *Bayes' Theorem* or *Bayes' Rule*, and is important because it gives a formula for \"flipping\" conditional probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From Bayes' Rule to Bayesian Inference\n",
    "\n",
    "If we replace these labels, we find the usual expression of Bayes' theorem as it relates to model fitting:\n",
    "\n",
    "$$\n",
    "P(\\theta \\mid D) = \\frac{P(D\\mid\\theta)P(\\theta)}{P(D)}\n",
    "$$\n",
    "\n",
    "Technically all the probabilities should all be conditioned on the information $I$:\n",
    "\n",
    "$$\n",
    "P(\\theta \\mid D,I) = \\frac{P(D \\mid \\theta,I)P(\\theta \\mid I)}{P(D \\mid I)}\n",
    "$$\n",
    "\n",
    "Recall $\\theta$ is the model we're interested in, $D$ is the observed data, and $I$ encodes all the prior information, including what led us to choose the particular model we're using."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wherefore the Controversy?\n",
    "\n",
    "When we write Bayes' rule this way, we're all of a sudden doing something controversial: can you see where this controversy lies?\n",
    "\n",
    "Two controversial points:\n",
    "\n",
    "- We have a probability distribution over model parameters. A frequentist would say this is meaningless!\n",
    "\n",
    "- The answer depends on the prior $P(\\theta\\mid I)$. This is the probability of the model parameters without any data: how are we supposed to know that?\n",
    "\n",
    "Nevertheless, applying Bayes' rule in this manner gives us a means of quantifying our knowledge of the parameters $\\theta$ given observed data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploring the Terms in Bayesian Inference\n",
    "\n",
    "We have four terms in the above expression, and we need to make sure we understand them:\n",
    "\n",
    "#### $P(\\theta\\mid D, I)$ is the *posterior*.\n",
    "This is the quantity we want to compute: our knowledge of the model given the data & background knowledge (including the choice of model).\n",
    "\n",
    "#### $P(D\\mid\\theta,I)$ is the *likelihood*.\n",
    "This measures the probability of seeing our data given the model. This is identical to the quantity maximized in frequentist *maximum-likelihood* approaches.\n",
    "\n",
    "#### $P(\\theta\\mid I)$ is the *prior*.\n",
    "This encodes any knowledge we had about the answer before measuring the current data.\n",
    "\n",
    "#### $P(D\\mid I)$ is the *Fully Marginalized Likelihood*\n",
    "You might prefer the acronym *FML* (it's also called the *Evidence* among other things). In the context of model fitting, it acts as a normalization constant and in most cases can be ignored."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside on the FML\n",
    "\n",
    "In general the Fully Marginalized Likelihood (FML) is **very costly** to compute, which makes the acronym doubly appropriate in any situation where you actually need it.\n",
    "To see why it's so costly, consider that the FML can be expressed as an integral using the identities we covered above:\n",
    "$$\n",
    "P(D\\mid I) = \\int P(D\\mid\\theta, I) P(\\theta) d\\theta\n",
    "$$\n",
    "In other words, it is the integral over the likelihood for *all possible values of theta*.\n",
    "When your likelihood is a complicated function of many parameters, computing this integral can become extremely costly (a manifestation of the *curse of dimensionality*).\n",
    "\n",
    "In general, for **model fitting**, you can ignore the FML as a simple normalization term. In **model selection**, the FML can become important."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the Point?\n",
    "\n",
    "At first blush, this might all seem needlessly complicated. Why not simply maximize the likelihood and be done with it? Why multiply by a prior at all?\n",
    "\n",
    "There are a couple good reasons to go through all of this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \"Purity\"\n",
    "\n",
    "Many advocates of the Bayesian approach argue for it's philosophical purity: you quantify knowledge in terms of a probability, then follow the math to compute the answer.\n",
    "The fact that you need to specify a prior might be inconvenient, but we can't simply pretend it away.\n",
    "There are good reasons to think that the Bayesian posterior is just the quantity we wish to compute; in that case we should compute it, however inconvenient!\n",
    "\n",
    "Perhaps the most vocal 20th century proponent of this view was Jaynes; I'd highly suggest looking at his book, *Probability Theory: The Logic of Science* ([PDF here](http://bayes.wustl.edu/etj/prob/book.pdf))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter Uncertainties\n",
    "\n",
    "Whether frequentist or Bayesian, the maximum likelihood \"point estimate\" is only a small part of the picture. What we're really interested in scientifically is the *uncertainty* of the estimates. So simply reporting a point estimate is not appropriate.\n",
    "\n",
    "In frequentist approaches, \"error bars\" are generally computed from *Confidence Intervals*, which effectively measure $P(\\hat\\theta\\mid\\theta)$, rather than $P(\\theta\\mid D)$.\n",
    "It takes some mental gymnastics to relate the confidence interval to the quantity we as scientists have in mind when we say \"uncertainty\".\n",
    "\n",
    "In the Bayesian approach, we are actually measuring $P(\\theta\\mid D)$ from the beginning.\n",
    "\n",
    "For some approachable reading on frequentist vs. Bayesian uncertainties, I'd suggest [The Fallacy of Placing Confidence in Confidence Intervals](http://learnbayes.org/papers/confidenceIntervalsFallacy/), as well as my (rather opinionated) blog post on the topic, [Confidence, Credibility, and why Frequentism and Science do not Mix](http://jakevdp.github.io/blog/2014/06/12/frequentism-and-bayesianism-3-confidence-credibility/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Marginalization of nuisance parameters\n",
    "\n",
    "Additionally, Bayesian approaches offer a very natural way to systematically account for nuisance parameters.\n",
    "\n",
    "Consider that we now have a recipe for computing the posterior, $P(\\theta\\mid D,I)$. Recall that in general our model parameters $\\theta$ contain some parameters of interest (what we called $\\theta_S$, or the \"science\") and some nuisance parameters (what we called $\\theta_N$). That is, we're really measuring\n",
    "\n",
    "$$\n",
    "P(\\theta_S,\\theta_N\\mid D, I)\n",
    "$$\n",
    "\n",
    "What we're really interested in, though, is $P(\\theta_S\\mid D, I)$ which we can compute via a straightforward *marginalization* integral:\n",
    "\n",
    "$$\n",
    "P(\\theta_S\\mid D, I) = \\int P(\\theta_S, \\theta_N\\mid D, I)d\\theta_N\n",
    "$$\n",
    "\n",
    "where the integral is over the entire space of $\\theta_N$.\n",
    "This quantity, the marginalized posterior, is our final result of interest.\n",
    "\n",
    "At first glance, you might think this problem would be nearly as difficult as the computation of the FML above.\n",
    "We'll see later, however, that a feature of the MCMC approaches typical for large Bayesian problems is that such marginalized likelihoods come essentially for free."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Looking forward\n",
    "\n",
    "We now have all the probabilistic machinery we need to do Bayesian inference... next we will get our fingers on our keyboards and fit some Bayesian models!"
   ]
  }
 ],
 "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
}