Demo 5: Running an MCMC analysis of supernova data
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:
OR if you are feeling adventurous and have MPI set up:
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 demo5.txt with 25600 samples:
#cosmological_parameters--omega_m cosmological_parameters--h0 supernova_params--deltam supernova_params--alpha supernova_params--beta LIKE # ... [SNIP - all the parameters you ran with are saved here, along with other info like who to cite] 0.375349838504 0.543069181465 0.328431291151 0.124554370328 3.28955217497 -2227.10121001 0.277715836563 0.81440479776 6.95767929396 0.152362132558 8.60894281101 -27083.8222028 0.396752581802 0.664902793534 1.88367270414 0.127254401707 4.48364725723 -8182.54328716 0.215070227128 0.762109999036 5.29185049673 0.139453531438 6.46507792112 -26799.1925617 0.214404110344 0.781725633573 2.1142025756 0.141179467329 8.25565740506 -3161.17547089 # ...
(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:
postprocess --burn 5000 -o plots -p demo5 demos/demo5.ini
./cosmosis/plotting/marginal_densities.r -b 10000 -o plots/ 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:
Samples after cutting: 20600 Marginalized mean, std-dev: cosmological_parameters--omega_m = 0.29728 ± 0.0328477 cosmological_parameters--h0 = 0.738122 ± 0.0249374 supernova_params--deltam = 0.040323 ± 0.0768359 supernova_params--alpha = 0.137666 ± 0.00633048 supernova_params--beta = 3.10945 ± 0.0814108 like = -345.562 ± 1.55596 Marginalized median, std-dev: cosmological_parameters--omega_m = 0.296391 ± 0.0328477 cosmological_parameters--h0 = 0.73781 ± 0.0249374 supernova_params--deltam = 0.0412422 ± 0.0768359 supernova_params--alpha = 0.137889 ± 0.00633048 supernova_params--beta = 3.10694 ± 0.0814108 like = -345.255 ± 1.55596 Best likelihood: cosmological_parameters--omega_m = 0.292377 cosmological_parameters--h0 = 0.734679 supernova_params--deltam = 0.0270577 supernova_params--alpha = 0.137495 supernova_params--beta = 3.09374 like = -343.068
Also the plots directory should also now contain a collection of 1D and 2D histogram plots in pleasant shades of blue, like this one:
You should also find some files with the statistics above saved in a more structured form.
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:
[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:
[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:
[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 ; ...
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.