1. Joe Zuntz Avatar Joe Zuntz
  2. Untitled project
  3. cosmosis

Wiki

Clone wiki

cosmosis / demos / Demo5

View History

Demo 5: Running an MCMC analysis of supernova data

Running

In this example we will run a full MCMC analysis on data from the SDSS-II/SNLS3 Joint Light-curve Analysis (JLA). See http://arxiv.org/abs/1401.4064 for more details about this data.

Run the sampler EITHER like this:

#!bash

cosmosis demos/demo5.ini

OR if you are feeling adventurous and have MPI set up: 

#!bash

mpirun -n 4 cosmosis --mpi demos/demo5.ini

This will take a few minutes to run, depending on your machine. On a Mac you may be asked if python 2.7 should accept incoming network connections - click Allow. You should produce a file output/demo5.txt with 25600 samples:

#!bash

#cosmological_parameters--omega_m       supernova_params--deltam        supernova_params--alpha supernova_params--beta  supernova_params--m     prior   post
# ... [SNIP - all the parameters you ran with are saved here, along with other info like who to cite]
0.29320063566282856     0.0417506764974888      0.13676152572956485     3.1018523606461814      -18.999731746108633     2.407945608651872       -370.0034237651477
0.2934489095354855      0.04092649325195843     0.1367657088117204      3.103591688420887       -18.996000448689696     2.407945608651872       -371.16109978382906
0.29352400617815466     0.041134195900012244    0.1368214157523774      3.1032833459437397      -19.00611268644998      2.407945608651872       -367.2620090287608
0.29383036064910295     0.0402216284783307      0.13676043848518205     3.1033285316438346      -18.997494359334798     2.407945608651872       -370.22872079684635
# ...

(N.B. The actual numbers you get will be different to the above because it is a random process.) Now let's generate some plots to visualize the chain:

#!bash
postprocess  --burn 5000  -o plots -p demo5 demos/demo5.ini
or 
#!bash
./cosmosis/plotting/marginal_densities.r -b 10000 -o plots/ output/demo5.txt -p demo5 --verbose --fill

Note that the R script produces plots similar to those of the Python script, but does not yet produce textual statistical summaries.

You will get output which will include summary information like this:

#!text
    Samples after cutting: 20600

Marginalized mean, std-dev:
    cosmological_parameters--omega_m = 0.298218 ± 0.0341413 
    supernova_params--deltam = -0.0696383 ± 0.0239679 
    supernova_params--alpha = 0.141714 ± 0.00663996 
    supernova_params--beta = 3.1088 ± 0.0764405 
    supernova_params--m = -19.0489 ± 0.0240827 
    prior = 2.40795 ± 4.44089e-16 
    post = -342.583 ± 1.58725 


Marginalized 1D peak, 68% asymmetric error bars:
    cosmological_parameters--omega_m = 0.286266 + 0.0446974 - 0.0227785 
    supernova_params--deltam = -0.0676792 + 0.0227357 - 0.0244075 
    supernova_params--alpha = 0.142259 + 0.00637198 - 0.00692953 
    supernova_params--beta = 3.11276 + 0.0646293 - 0.0821559 
    supernova_params--m = -19.0499 + 0.0245902 - 0.0230999 
    prior = nan + nan - nan 
    post = -341.484 + 0.876216 - 1.82348 

...

Best likelihood:
    cosmological_parameters--omega_m = 0.304464
    supernova_params--deltam = -0.0649367
    supernova_params--alpha = 0.140684
    supernova_params--beta = 3.10063
    supernova_params--m = -19.0464
    prior = 2.40795
    post = -340.11

...

Also the plots directory should also now contain a collection of 1D and 2D histogram plots in pleasant shades of blue, like this one:

2D_supernova_params--deltam_cosmological_parameters--h0.png

You should also find some files with the statistics above saved in a more structured form.

Understanding

The sampling

We are now using a proper MCMC sampler, called emcee. This sampler is described in http://arxiv.org/abs/1202.3665 .

Here's how we tell cosmosis to use it in demos/demo5.ini:

#!ini
[emcee]
walkers = 64
samples = 1000
nsteps = 100
You can find out more about the meaning of these parameters by looking in the ini file itself.

Setting up a cosmosis sampling run involves choosing (or writing) the modules that you want to run, making sure that each produces the quantities that later ones need, and then choosing what likelihoods to take out at the end.

We specify our pipeline like this in the same ini file:

#!ini

[pipeline]
; We use two likelihoods, the JLA (for high redshift) and
; Riess 2011 to anchor H0, which is otherwise degenerate
; with the nuisance parameter M
modules = consistency camb jla riess11
values = demos/values5.ini
likelihoods = jla riess
We take the sum of two likelihoods, jla and riess to get our constraints.

Once again each module in the list above is specified by an ini file section, with options in it. The JLA module, for example, takes quite a lot of parameters pointing to the data files it needs:

#!ini

[jla]
file = cosmosis-standard-library/supernovae/jla_v3/jla.so
data_dir = cosmosis-standard-library/supernovae/jla_v3/data
data_file = jla_lcparams.txt
scriptmcut = 10.0
mag_covmat_file = jla_v0_covmatrix.dat
; ...

The plotting

Our chain analysis program postprocess uses a technique called kernel density estimation (KDE) to smooth the 1D and 2D histograms and form smooth edges.

We use the --burn flag to burn the first 5000 samples. We also specified the output directory and prefix.

You can get more information on postprocess by running it with the -h flag.

Updated