110. Introducing doubt in Bayesian statistics 1

This is a guest post by my old friend Pascal Vaudrevange, who was a student of Lev Kofman when we were together in CITA, and is now a postdoc at Case Western University working with Glenn Starkman. Dmitry.

I have been kindly given the opportunity by Dmitry to talk a little bit about the paper “Introducing doubt in Bayesian statistics” (arXiv:0811.2415 physics:data-an) that came out of a visit of Roberto Trotta (then at Oxford, now at Imperial College) here at Case Western Reserve University in collaboration with Glenn Starkman. In this paper, we give a prescription how to evaluate the goodness-of-fit of a single model in a Bayesian framework. (As opposed to comparing different models using Bayesian model comparison which is well established. Although, truth to be told, in effect we also perform a model comparison, but with a fictional best-fit model.)

This will happen in several parts. This first part lays out the difficulties in determining a “good” model using the example of dark energy and mentions some differences between the Bayesian and the frequentist approach to statistics.

First of all, I should mention that in recent years, cosmological observations have delivered data of such quality that it may be possible (at least in principle) to perform precision checks of cosmological models (as in, the error bars on the parameters of any individual model can be made quite small). But it is not clear how to do determine whether the correct model has been found.

Let me use as example the mystery of dark energy which represents about 75% of the energy density in the universe: . The only thing we know for certain about it is its equation of state (the ratio of its pressure and density) today: . If this equation of state were the same at all times (or equivalently at all redshifts), then dark energy would just be Einstein’s famous cosmological constant. However, if this ratio changes over time, there should be some dynamical mechanism that governs it’s change. In other words, if is not constantly -1, dark energy should be some (yet undiscovered) particle. A very exciting prospect indeed!

Now how can we determine this ratio of pressure to density at different times? Without going into too much detail, there are severa observational data sets providing insight into this, ranging from finding distant supernovae to using a feature in the distribution of matter in the universe going by the name of baryon acoustic oscillation. Unfortunately, it is impossible to solve all the equations and obtain w as a function of time given the measured data.

Instead, one uses clusters of fast computers to approach this problem in a slightly different fashion. (In the following we shall focus on using only supernova data for simplicity.)

One picks a model for as a function of redshift and fixes parameter to some numbers, e.g. if we pick the model , the parameter would be . Also, one needs to fix several other parameters governing the evolution of the universe, such as the amount of dark matter, baryonic matter (meaning matter that consists of quarks), the value of the Hubble parameter, etc..

Then one computes the predicted brightness of a supernova at a given distance, and compares this to the observed brightness of the supernova at the given distance. This whole process is repeated for each supernova, and in the end one ends up with the probability for the observed data given the chosen values for the parameters of the model (in this case ). Then one (almost randomly) chooses a different value for the parameter , say , and repeats this whole process. In this way (which is called a Markov-Chain Monte Carlo (or MCMC) process) one will obtain a list of parameter values distributed according to the posterior, i.e. plotting a histogram with those values will give the probability distribution of the data given the parameters.

For a frequentist, this distribution is the same as the probability of the parameter given the data, whereas for a Bayesian, one must first make use of Bayes theorem which relates the probability of the parameters given the data to the probability of the data given the parameters :

(1)

where is the so called prior on the data, and is the prior on the parameters. These two numbers are responsible for the – mildly put – dislike expressed by some people about Bayesian statistics, as one is free to choose them to be any number one likes. The only condition is that this should happen before one peaks at the data set, i.e. before one runs the MCMC.

Thus after first choosing values for and and then performing the Markov Chain Monte Carlo, we are left with the probability distribution for the parameter “” of our model, indicating the value of constant .

But what if changes with time? Well, in this case we can run the MCMC again, but this time with a different model, e.g. . This would allow the equation of state parameter to change linearly with redshift, and so we would obtain from the MCMC the posterior distribution of both and . If the posterior distribution of was sufficiently different from 0, for a frequentist this would be evidence for a non-constant .

In a true Bayesian fashion however, one would have to compute what is called the evidence for both models and given by

(2)
(3)

where , , are the priors, and is the posterior for obtained from the MCMC and is the posterior for model .

The ratio of both evidences will then indicate which of the models is the preferred one (notice that the constant drops out of the ratio): Using the so-called Jeffrey’s scale, this number can then be interpreted as strong or weak evidence for either the model of or linearly changing with .

There is only one problem left: what if was really a function of the form ? In order to definitely determine whether is a constant, one would have to compare this model with all possible parameterizations of as a function of redshift. Somehow this seems like a daunting task… The next part of my post will explain how we (Glenn, Roberto and me) propose to address this problem.

To be continued.