Bayesian Inference II: Fitting a Straight Line

Phil Marshall, Daniel Foreman-Mackey and Dustin Lang

Astro Hack Week, New York, September 2015

This is AstroHackWeek/inference/straightline.ipynb

Goals:

• Set up and carry out a simple Bayesian inference, characterizing a simple posterior PDF

• Compare brute force, analytic and MCMC sampled results

• Check models, with posterior predictive distributions of test statistics

• Use the Bayesian Evidence to compare various models, and understand its properties

• Suggest some simple hacks to get some practice with Bayesian inference

The Data

• Let's generate a simple dataset: observations of $y$ with reported uncertainties $\sigma_y$, at given $x$ values.

• Yours will look different to mine, as your random seed will be different. That's good - we can do some (completely unrealistic) comparisons over datasets.

Least-squares Straight Line Fit

• An industry standard: find the slope and intercept of a straight line that minimizes the mean square residual. Since the predicted data depend linearly on the parameters, the least squares solution can be found by a matrix inversion and multiplication, conveneniently packed in numpy.linalg.

• Suppose that the $y$ values are related to the $x$ values by $y = m x + b$, where $x$ and $y$ are "vectors" (numpy arrays). We might say that we want to find $m$ and $b$ that minimize

$W^2 = \sum_k \frac{(y_k - m x_k - b)^2}{\sigma_k^2}$

• Here, $1/\sigma_k^2$ is a sensible weighting factor that reduces the influence of $y_k$ values with large uncertainty $\sigma_k$

• Differentiating $W^2$ with respect to $m$ and $b$ and setting the result to zero leads to a matrix equation that can be solved - this is what np.linalg.lstsq is doing.

Probabilistic Fitting

• The Bayesian answer to the question "What are the values of $b$ and $m$ in a straight line model for the above data?" is the posterior PDF for those parameters, ${\rm Pr}(m,b|y)$

• The posterior PDF doesn't tell us the best-fit parameter combination - it gives us the means to simulate all plausible parameter combinations, given the data.

• Let's derive this posterior PDF, starting with a probabilistic graphical model.

• The PGM illustrates the forward problem of simulating a dataset $\{y_k\}$ by drawing random numbers from conditional probability distributions, which are depicted by the nodes. The conditional dependences are depicted by the edges.

• The PGM is a complete description of the forward problem: it is a representation of the joint PDF for all our variables, ${\rm Pr}(m,b,\{y_k\},\{y^{\rm true}_k\}\,|\,\{x_k\},\{\sigma_k\})$, factorized such that we can simulate data:

${\rm Pr}(m,b,\{y_k\},\{y^{\rm true}_k\}\,|\,\{x_k\},\{\sigma_k\},H_1) = \prod_k {\rm Pr}(y_k\,|\,y^{\rm true}_k,\sigma_k,H_1) \;{\rm Pr}(y^{\rm true}_k\,|\,x_k,m,b,H_1) \;\;\;{\rm Pr}(m,b\,|\,H_1)$

Questions:

1. Where did that product come from? What assumption does it encode?

2. Why are the $\{x_k\}$ and the $\{\sigma_k\}$ always on the right hand side of the bar, but $\{y_k\}$ is sometimes on the left?

3. What is the meaning of "$H_1$"? Hint: notice that it too is always on the right hand side of the bar.

4. What functional form does the conditional PDF ${\rm Pr}(y^{\rm true}_k\,|\,x_k,m,b,H_1)$ have?

5. What functional form do you think we should assume for the sampling distribution, ${\rm Pr}(y_k\,|\,y^{\rm true}_k,x_k,\sigma_k,H_1)$? Hint: imagine you were generating mock data.

6. What functional form should we assume for ${\rm Pr}(m,b\,|\,H_1)$?

Talk to your neighbor(s) for 5 minutes about these, note down some brief answers, and be prepared to provide them when called upon.

The Inverse Problem

• Let's simplify our joint probability expression a bit, by marginalizing out the latent variables $\{y^{\rm true}_k\}$, and absorbing the $\{x_k\}$ and $\{\sigma_k\}$ into $H_1$.

${\rm Pr}(m,b,\{y_k\}|H_1) = \int\, \prod_k {\rm Pr}(y_k\,|\,y^{\rm true}_k,\sigma_k,H_1) \;{\rm Pr}(y^{\rm true}_k\,|\,x_k,m,b,H_1) \;\;\;{\rm Pr}(m,b\,|\,H_1) \; d y^{\rm true}_k $

$\longrightarrow {\rm Pr}(m,b,\{y_k\}|H_1) = \prod_k {\rm Pr}(y_k\,|\,(m x_k + b),\sigma_k,H_1) \;\;\;{\rm Pr}(m,b\,|\,H_1) $

• Now we just factorize the left hand side and rearrange (like Bayes did), to get:

${\rm Pr}(m,b|\{y_k\},H_1) = \frac{1}{Z} \prod_k {\rm Pr}(y_k\,|\,(m x_k + b),H_1) \;\;\;{\rm Pr}(m,b\,|\,H_1) $

• Here, $Z = {\rm Pr}(\{y_k\}\,|\,H_1)$, and is known as the "evidence" for $H_1$. It also goes by the name "Fully Marginalized Likelihood" because

$Z = {\rm Pr}(\{y_k\}\,|\,H_1) = \int \prod_k {\rm Pr}(y_k\,|\,(m x_k + b),\sigma_k,H_1) \;{\rm Pr}(m,b\,|\,H_1)\;dm\,db$

and $\prod_k {\rm Pr}(y_k\,|\,(m x_k + b),\sigma_k,H_1)$ is the joint likelihood for the parameters. ${\rm Pr}(m,b\,|\,H_1)$ is the prior PDF for the parameters.

Evaluating the posterior PDF on a grid

• For simple, 2-dimensional parameter spaces like this one, evaluating on a grid is not a bad way to go. We'll see that the least squares solution lies at the peak of the posterior PDF - for a certain set of assumptions about the data and the model.

• To compute the evidence $Z$ and normalize the posterior PDF, we need to be able to integrate the likelihood function over the prior - on a grid, we can do this numerically, to normlize the posterior PDF.

• Un-normalized, the likelihood function can take very small values, so we usually work with the log likelihood instead.

Exercise:

With your neighbors, upgrade my log likelihood function below so that it runs faster, and you get a plot of the posterior PDF for $m$ and $b$ without having to wait.

MCMC Sampling

• In problems with higher dimensional parameter spaces, we need a more efficient way of approximating the posterior PDF - both when characterizing it in the first place, and then when doing integrals over that PDF (to get the marginalized PDFs for the parameters, or to compress them in to single numbers with uncertainties that can be easily reported).

• In most applications it's sufficient to approximate a PDF with a (relatively) small number of samples drawn from it; MCMC is a procedure for drawing samples from PDFs.

• The following piece of code implements the Metropolis algorithm, whose derivation can be read in many textbooks. You can read a brief justification aimed at Phil's grad student self in Chapter 2 of Phil's thesis, for example.

• This looks pretty good: no plateauing, or drift.

• A more rigorous test for convergence is due to Gelman & Rubin, and involves comparing the intrachain variance with the inter-chain variance in an ensemble. It's worth reading up on.

• Foreman-Mackey & Hogg recommend looking at the autocorrelation length, and whether it stablizes during the run.

• Note that we kept all the samples, and did not discard any at the beginning of the run as "burn-in". Nor did we thin the chain at all. Both of these are good ideas.

• Now let's look at the samples in parameter space, overlaid on our gridded posterior.

Multi-dimensional Visualization

• For a better view of a set of posterior samples, let's make a corner plot that shows all 1D and 2D marginalized distributions.

Exercise:

Now that you have seen the Metropolis sampler in action, let's try and break it.

With your neighbor(s), see if you can change the control parameters (ie. not the priors) to make the sampler fail in various ways - and be ready to report back in 5 minutes.

Anyway:

With a well-tuned sampler:

• It looks like we made a nice, precise measurement!

• Our next step is to check our model for accuracy.

Model Checking

• How do we know if our model is any good? There are two properties that "good" models have: the first is accuracy, and the second is efficiency. We'll return to th equestion of efficiency below.

• Accurate models generate data that is like the observed data. What does this mean? First we have to define what similarity is, in this context. Visual impression is one very important way, that's usually best done first.

• Test statistics that capture relevant features of the data are another. We'll look at the posterior predictive distributions for the datapoints, and for a particularly interesting test statistic, the reduced chi-squared.

Visual Model Checking

• Plot posterior predictions, in data space.

Posterior Predictive Model Checking

• If our model is the true one, then replica data generated by it should look like the one dataset we have. This means that summaries of both the real dataset, $T(d)$, and the replica datasets, $T(d^{\rm rep})$, should be drawn from the same distribution.

• If the real dataset was not generated with our model, then its summary may be an outlier from the distribution of summaries of replica datasets.

• We can account for our uncertainty in the parameters $\theta$ by marginalizing them out, which can be easily done by just making the histogram of $T(d^{\rm rep}(\theta))$ from our posterior samples, after drawing one replica dataset $d^{\rm rep}$ from the model sampling distribution ${\rm Pr}(d^{\rm rep}\,|\,\theta)$ for each one.

• Then, we can ask: what is the posterior probability for the summary $T$ to be greater than the observed summary $T(d)$? If this is very small or very large, we should be suspicious of our model - because it is not predicting the data very accurately.

${\rm Pr}(T(d^{\rm rep})>T(d)\,|\,d) = \int I(T(d^{\rm rep})>T(d))\,{\rm Pr}(d^{\rm rep}\,|\,\theta)\,{\rm Pr}(\theta\,|\,d)\;d\theta\,dd^{\rm rep}$

• Here $I$ is the "indicator function" - 1 or 0 according to the condition.

• If our model is true (and we're just uncertain about its parameters, given the data), we can compute the probability of getting a $T$ less than that observed, where T is the reduced chisq relative to a straight line with some reference $(m,b)$.

• Note that we did not have to look up the "chi squared distribution" - we can simply compute the posterior predictive distribution given our generative model.

• Still, this test statistic looks a little bit strange: it's computed relative to the best fit straight line - which is a sensible reference but somehow not really in the spirit of fitting the data...

• Instead, lets look at a discrepancy measure $T(d,\theta)$ that is a function of both the data and the parameters, and compute the posterior probability of getting $T(d^{\rm rep},\theta) > T(d,\theta)$:

${\rm Pr}(T(d^{\rm rep},\theta)>T(d,\theta)\,|\,d) = \int I(T(d^{\rm rep},\theta)>T(d,\theta))\,{\rm Pr}(d^{\rm rep}\,|\,\theta)\,{\rm Pr}(\theta\,|\,d)\;d\theta\,dd^{\rm rep}$

• The conclusion drawn from the discrepancy is more interesting, in this case. All our $\theta = (m,b)$ samples are plausible, so replica datasets generated by them should also be plausible. The straight line defined by each $(m,b)$ should go through the real data points as readily (on average) as it does its replica dataset.

• Do our posterior predictive $p-$values suggest we need to improve our model? What about the visual check?

• The reduced chi-squared is actually not such an interesting test statistic - because we already used it in the likelihood! Better discrepancy measures would stress-test different aspects of the model - such as the correlations between datapoints, or the behavior of the model in particular places, etc.

Model Expansion

• Maybe I see some curvature in the data - or maybe I have a new astrophysical idea for how this data was generated. Let's try adding an extra parameter to the model, to make a quadratic function:

$y = m x + b + q x^2$

• The coefficient $q$ is probably pretty small (we were originally expecting to only have to use a straight line model for these data!), so I guess we can assign a fairly narrow prior, centered on zero.

• All parameters are again precisely constrained.

• The gradient and intercept $m$ and $b$ are significantly different from before, though...

Checking the Quadratic Model

• How do the two models compare? Which one matches the data better?

All of the above discussion was about model checking - assessment of the ability of a model to make accurate predictions.

• We used replicated data, drawn from the model assuming it was true, and then looked at the probability of getting (a summary of) the data, given this assumption.

• Note the parallel with frequentist hypothesis testing:

Frequentist: assume the null hypothesis $H_0$, and calculate the probability of getting the estimators for the parameters by chance, by studying the distribution of a test statistic over a large ensemble of replica datasets (which under certain assumptions is a standard distribution).

Bayesian: assume the alternative hypothesis $H_1$, and calculate the posterior predictive probability distribution for the test statistic given the data, relative to the one test statistic we observed, marginalizing over the model parameters. In general this distribution will be non-standard, but we can often draw samples from the posterior PDF anyway.

• The Bayesian posterior predictive check is fairly conservative, by comparison with tests that focus on the best-fit parameters. This can be a good thing, to err on the side of caution.

Hack Suggestion: Bayesian Cross Validation with SciKit-Learn

• An interesting alternative to posterior predictive model checking is to follow the machine learning approach, and use cross-validation to assess model prediction accuracy.

• What should the algorithm be? Here's a starter suggestion:

   split data into training set dt and holdback test set dh
      infer Pr(x|dt), drawing posterior samples in x
      predict dhp
      compute score, using dhp and dh  
   combine scores

• What should the "score" be? Pr(dh|dt), marginalizing over x? Machine learning "Accuracy"? How should we combine the scores?

• scikit-learn already has a nice API, with (I guess) a Model class and some very useful cross-validation machinery. Maybe we could write BayesianModel that inherits from Model, and give it a sampler, priors, etc, as well as fit and predict methods (which could return ensembles of samples instead of point estimates).

Model Comparison with the Bayesian Evidence

• The evidence for model $H$, ${\rm Pr}(d\,|\,H)$, enables a form of Bayesian hypothesis testing via the evidence ratio (or "odds ratio", or "Bayes Factor"):

$R = \frac{{\rm Pr}(d\,|\,H_1)}{{\rm Pr}(d\,|\,H_0)}$

• This quantity is similar to a likelihood ratio, but it's a fully marginalized likelihood ratio - which is to say that it takes into account our uncertainty about values of the parameters of each model by integrating over them all.

• As well predictive accuracy, the other virtue a model can have is efficiency. The evidence summarizes all the information we put into our model inferences, via both the data and our prior beliefs. You can see this by inspection of the integrand of the fully marginalized likelihood (#FML) integral:

${\rm Pr}(d\,|\,H) = \int\;{\rm Pr}(d\,|\,\theta,H)\;{\rm Pr}(\theta\,|\,H)\;d\theta$

• The consequence of this is that the Evidence depends on both model accuracy (goodness of fit, ability to predict data) and efficiency (not over-fitting by having implausible flexibility).

• The following figure might help illustrate how the evidence depends on both goodness of fit (through the likelihood) and the complexity of the model (via the prior). In this 1D case, a Gaussian likelihood (red) is integrated over a uniform prior (blue): the evidence can be shown to be given by $E = f \times L_{\rm max}$, where $L_{\rm max}$ is the maximum possible likelihood, and $f$ is the fraction of the blue dashed area that is shaded red. $f$ is 0.31, 0.98, and 0.07 in each case.

The illustration above shows us a few things:

1) The evidence can be made arbitrarily small by increasing the prior volume: the evidence is more conservative than focusing on the goodness of fit ($L_{\rm max}$) alone. Of course if you assign a prior you don't believe, then you should not expect to get out a meaningful answer for ${\rm Pr}(d\,|\,H)$.

2) The evidence is linearly sensitive to prior volume ($f$), but exponentially sensitive to goodness of fit ($L_{\rm max}$). It's still a likelihood, after all.

• The odds ratio can, in principle, be combined with the ratio of priors for each model to give us the relative probability for each model being true, given the data:

$\frac{{\rm Pr}(H_1|d)}{{\rm Pr}(H_0|d)} = \frac{{\rm Pr}(d|H_1)}{{\rm Pr}(d|H_0)} \; \frac{{\rm Pr}(H_1)}{{\rm Pr}(H_0)}$

• Prior probabilities are very difficult to assign in most practical problems (notice that no theorist ever provides them). So, one way to interpret the evidence ratio is to note that:

  • If you think that having seen the data, the two models are still equally probable,
  • then the evidence ratio in favor of $H_1$ is you the odds that you would have had to have been willing to take against $H_1$, before seeing the data.
  • That is: the evidence ratio updates the prior ratio into a posterior one - as usual.

Calculating the Evidence

• The FML is in general quite difficult to calculate, since it involves averaging the likelihood over the prior. MCMC gives us samples from the posterior - and these cannot, it turns out, be reprocessed so as to estimate the evidence stably.

• A number of sampling algorithms have been developed that do calculate the evidence, during the process of sampling. These include:

  • Nested Sampling (including MultiNest and DNest)
  • Parallel Tempering, Thermodynamic Integration
  • ...

• If we draw samples from the prior, we can then estimate the evidence via the usual sum over samples,

${\rm Pr}(d\,|\,H) \approx \sum_k\;{\rm Pr}(d\,|\,\theta,H)$

• Sampling the prior and computing the likelihood at each sample position is called "Simple Monte Carlo", and while it works in certain low-dimensional situations, in general it is very inefficient (at best).

• Still, let's give it a try on our two models, and attempt to compute the Bayes Factor

$R = \frac{{\rm Pr}(d\,|\,{\rm quadratic})}{{\rm Pr}(d\,|\,{\rm straight line})}$

Now that we have the likelihood for each sample, let's check that we did actually sample the posterior well. Here are the corner plots (note the weights, and the need to not plot the points as well as the contours!):

Question:

• Why does $q$ appear in this plot? Does it affect the evidence calculation?

Exercise:

Play around with the Simple Monte Carlo sampler, and see if you can break it.

• Try running it a few times, with the same N: how precise these evidence estimates from 50000 samples?

• We said we believed a priori that $q$ should be small, and set the prior width accordingly. What would happen if we started to doubt ourselves, and re-assigned a broader prior? Try running with some larger values for the Gaussian prior width $q[1]$, and note what happens to the Bayes Factor.

Conclusions

• The Bayesian evidence is qualitatively different from other model assessments. While they focus primarily on prediction accuracy, the evidence is the way in which information from the prior PDF propagates through into our posterior beliefs about the model as a whole.

• There are no mathematical limitations to its use, in contrast to various other hypothesis tests that are only valid under certain assumptions (such as the models being nested). Any two models can be compared and the odds ratio computed.

• Good things to keep in mind when using, or reading about, the evidence for a model:

  • Garbage in, garbage out: if you don't believe your prior, you won't believe your evidence.

  • The evidence is only linearly sensitive to prior volume, but exponentially sensitive to goodness of fit. If you want a clearer sense of which model makes the data more probable, get more data.

  • If you don't believe your priors, the evidence can be a distraction from things you care about more - like accuracy, or cost.

  • If you do believe your priors, the evidence is very useful: it appears in second-level inferences, when parts of your model that you previously considered to be constant now need to be varied and inferred. The FML is only one parameter thaw away from becoming a likelihood.

.

More Hacks!

• Intrinsic Scatter: Implement the log likelihood for the straight line plus intrinsic scatter, and compare it with the other models via the evidence. This is more of an exercise than a hack, but it's good to practice.

• Better Discrepancy Measures: Design and implement some more interesting discrepancy measures $T(d,\theta)$ that do a better job of highlighting the failures of the straight line model.

• Improve This Notebook: The coding in this notebook leaves a lot to be desired. The samplers could be packaged into def's for much easier re-use, and the plots could be combined to allow much easier comparison. Also, who does linear algebra these days? Surely scikit-learn has a linear model we could use. Pull requests are most welcome!

• Simple Monte Carlo: Write a general purpose Simple Monte Carlo sampler whose API mimics something well-used (like emcee), package it up, and compete for business. Can you automate the calculation of the evidence to a desired precision? If you can, I think that will probably guarantee that you have sampled the posterior well.

• Simple Monte Carlo Sci-Kit Learn Model: So that the machine learning community can experiment with scoring with the Bayesian evidence.

.

Copyright (C) Phil Marshall, Daniel Foreman-Mackey & Dustin Lang

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.