Clone wiki

cosmosis / demos / Demo9

Demo 9: Bayesian evidence with Multinest; scatter and custom plots


Bayesian evidence is a method for model selection - it tells us how likely an entire cosmological model is, rather than just a particular choice of parameters in that model.

In this demo we will use the sampler Multinest to calculate the Bayesian evidence. As a by-product we can also construct all the usual parameter constraint plots, or get a set of posterior samples. We'll use the same likelihood and pipeline as demo 5, the JLA SDSS supernova sample.

We'll also show you how to add additional custom plots as well as the standard ones that are produced by postprocess. That way you can use the cosmosis system that understands how to load and interpret the output files with your own custom plotting.

Run the multinest sampler like this:

cosmosis demos/demo9.ini

This should be a little faster than the emcee example - multinest is clever!

After a preamble about the parameters and supernova data files, you should see something like this:

 MultiNest v3.7
 Copyright Farhan Feroz & Mike Hobson
 Release June 2014

 no. of live points =  500
 dimensionality =    5
 Starting MultiNest
 generating live points
 live points generated, starting sampling
Acceptance Rate:                        0.994575
Replacements:                                550
Total Samples:                               553
Nested Sampling ln(Z):            -130739.431350
Importance Nested Sampling ln(Z):    -376.884711 +/-  0.999095

The last part will repeat until convergence, which should take a few minutes (about 15,000 samples are needed). At the end you'll get a print out of some final statistics, which should look something like this:

Saving 8231 samples
 ln(ev)=  -356.16416856202284      +/-  0.14586312956518516     
 Total Likelihood Evaluations:        15063
 Sampling finished. Exiting MultiNest

The output file output/demo9.txt will contain this:

#cosmological_parameters--omega_m       supernova_params--deltam        supernova_params--alpha supernova_params--beta  supernova_params--m     prior   like    post    weight
# ... [SNIP extra info] ...
0.39559108912944796     -0.7934505939483643     0.1350554633140564      2.812785506248474       -20.956758320331573     2.407945608651872       -11881.874933654408     -11879.466988045757     0.0
0.2770343512296677      -0.7635279893875122     0.1248539662361145      2.390855610370636       -20.933012008666992     2.407945608651872       -11538.145579291855     -11535.737633683204     0.0
0.24225930273532867     -0.8283501863479614     0.14984524011611938     3.0957254767417908      -20.94237756729126      2.407945608651872       -10925.908083130655     -10923.500137522004     0.0
0.36388089656829836     -0.36679399013519287    0.13541900396347045     2.115124434232712       -20.971498131752014     2.407945608651872       -10691.700011258814     -10689.292065650163     0.0
0.16651640236377716     -0.7750347852706909     0.13898276090621947     3.1154564321041107      -20.92144203186035      2.407945608651872       -10221.547816996723     -10219.139871388072     0.0
# ...

Don't worry about the zero weight! That's just for the early samples. (If you're used to Metropolis-Hasting chains, most of your expectations for the structure of this file will be wrong!) Let's generate some plots and stats, including some extra custom ones:

postprocess -o plots -p demo9 demos/demo9.ini --extra demos/

You'll get some summary stats, this time including the evidence value:

Marginalized mean, std-dev:
    cosmological_parameters--omega_m = 0.295753 ± 0.0336598
    cosmological_parameters--h0 = 0.738363 ± 0.0241566
    supernova_params--deltam = 0.0410873 ± 0.074628
    supernova_params--alpha = 0.136857 ± 0.00635137
    supernova_params--beta = 3.10989 ± 0.081146
    like = -345.526 ± 1.51149
    weight = 0.000348284 ± 0.000128698

Marginalized median, std-dev:
    cosmological_parameters--omega_m = 0.295503 ± 0.0336598
    cosmological_parameters--h0 = 0.738528 ± 0.0241566
    supernova_params--deltam = 0.0444039 ± 0.074628
    supernova_params--alpha = 0.136782 ± 0.00635137
    supernova_params--beta = 3.10814 ± 0.081146
    like = -345.222 ± 1.51149
    weight = 0.000390824 ± 0.000128698

Best likelihood:
    cosmological_parameters--omega_m = 0.299398
    cosmological_parameters--h0 = 0.740779
    supernova_params--deltam = 0.0486396
    supernova_params--alpha = 0.136235
    supernova_params--beta = 3.08996
    like = -343.073
    weight = 0.00589475

Bayesian evidence:
    log(Z) = -356.776 ± 0.145863


The same plots as before will now be generated, in plots/demo9_*, but with two custom ones that we specified in One is a scatter plot using color to show a third variable; you can easily add more such plots for MCMC or Multinest runs. The other is more specific to multinest: an annotated variant of figure 2 from Allison & Dunkley:




The sampling

The multinest code and algorithm are described in Feroz, Hobson & Bridges. There's also a nice summary and comparison to emcee and Metropolis in Allison & Dunkeley.

Here's how we tell cosmosis to use multinest:

sampler = multinest


In this case the sampler did not need the full maximum 50,000 iterations to reach a good level of accuracy - it suceeded long before then. live_points is the size of the multinest set of points; increasing it will get you more samples from the posterior at the end at the cost of more time. A few hundred is typical. If we had set the outfile_root parameter multinest would have given us a large number of different output files of its own - if you're a mulitnest guru you can use those directly.

The rest of our ini file is the same as demo 5, since we have not changed the pipeline.

The plotting

The chain files produced by multinest are a little more complicated than with emcee - they are not equally weighted posterior samples but have an extra "weight" column attached. The postprocess code knows how to handle these weights and uses them in the kernel density estimation process to produce the constraints.

Extra plots

The --extra flag told cosmosis to look in demos/ for a python script containing extra plots to make. You can create any new postprocessing step - it could be a plot, some statistics, or something else completely new.

postprocess searches the file you give it for any new post-processing steps, which must be python classes that inherit from a base post-processing step clas. There are more details in comments.

The first plot we make with this file is the scatter plot. The cosmosis plotting tools already contain the code for this plot; all we have to do is tell it which columns to use for plotting. We define it like this:

class ColorScatter(plots.MultinestColorScatterPlot):
    x_column = "supernova_params--deltam"
    y_column = "cosmological_parameters--h0"
    color_column = "cosmological_parameters--omega_m"
    scatter_filename = "scatter_dm_h0_omm"

The first three lines show the columns to use for the x, y, and color data respectively. They also provide the base filename to save to.

If you want a number of color scatter plots you could create multiple versions of this: ColorScatter1, ColorScatter2, etc.

The example is specific to multinest. It shows how to create a completely new plot by writing the "run" method:

class NestPlot(plots.Plots, plots.MultinestPostProcessorElement):
    def run(self):
        #Get the columns we need, in reduced form 
        #since this is not MCMC
        weight = self.weight_col()
        weight = weight/weight.max()
        #Sort the final 500 samples, since they 
        #are otherwise un-ordered
        weight[-500:] = np.sort(weight[-500:])

        #x-axis is just the iteration number
        x = np.arange(weight.size)

        #Choose a filename and make a figure for your plot.
        filename = self.filename("nest_weight")
        figure = self.figure(filename)

        #Do the plotting, and set the limits and labels
        pylab.plot(x, weight/weight.max())
        #pylab.xlim(-353, -341.5)
        pylab.xlabel("Iteration number $i$")

        #Add some helpful text, on the title and the plot itself
        pylab.title("Normalised posterior weight $p_i$", size='medium')
        pylab.text(500, 0.4, "Increasing likelihood,\n volume still large", size="small")
        pylab.text(5500,0.65, "Volume\ndecreasing", size='small')
        pylab.text(6500, 0.12, "Final $n_\mathrm{live}$\npoints", size='small')

        #Return the list of filenames we made
        return [filename]

In the run method you put ordinary pylab commands to draw whatever you want, except that the figure (matplotlib plotting space) should be generated using the commands shown below - we ask our parent class for a full filename to use (including the prefix and directory) and for the figure to generate.

To get the data for the plot we use method inherited from cosmosis - the reduced_col method gets a particular named column (and removes the burn-in in MCMC or preliminary samples in multinest).

Finally we return a list of any files we have created.