Wiki
Clone wikicosmosis / example_a
Example A: Simple CosmoMC-like Metropolis-Hastings analysis for WMAP9
Running
In this example we will run a Metropolis-Hastings MCMC on WMAP9 data with camb. If you have a good covariance matrix MH can be extremely efficient. Since it's more efficient to sample on Omega_m h^2 than Omega_m in this space we will also demonstrate a "consistency" module which converts between parameterizations.
We will do this with just a single chain, but compile with camb's open-mp facility to make use of multiple CPU cores (if your computer has them).
First, let's recompile with openmp:
#!bash make clean COSMOSIS_OMP=1 make
Whether you want open-mp or not depends on what kind of jobs you are doing. If using CAMB with a serial sampler like metropolis it will provide a good speed-up on multi-core machines.
This is a longer example that will take about an hour to run on a four-core machine. You can press ctrl-c to finish early if you get bored:
#!bash
cosmosis examples/example_a.ini
#!bash #cosmological_parameters--ommh2 cosmological_parameters--h0 cosmological_parameters--ombh2 cosmological_parameters--tau cosmological_parameters--n_s cosmological_parameters--a_s cosmological_parameters--omega_m cosmological_parameters--omega_b like # [ - SNIP EXTRA INFO - ] 0.128580559015 0.725501943148 0.0224232500497 0.0845032166586 0.975169764135 2.12083201207e-09 0.244285758862 0.0426011575661 -3780.64706738 0.128580559015 0.725501943148 0.0224232500497 0.0845032166586 0.975169764135 2.12083201207e-09 0.244285758862 0.0426011575661 -3780.64706738 0.128580559015 0.725501943148 0.0224232500497 0.104527759608 0.977514119639 2.21135177454e-09 0.244285758862 0.0426011575661 -3780.87309981 0.128580559015 0.725501943148 0.0224232500497 0.104527759608 0.977514119639 2.21135177454e-09 0.244285758862 0.0426011575661 -3780.87309981 0.128580559015 0.725501943148 0.0224232500497 0.104527759608 0.977514119639 2.21135177454e-09 0.244285758862 0.0426011575661 -3780.87309981 0.128580559015 0.725501943148 0.0224232500497 0.104527759608 0.977514119639 2.21135177454e-09 0.244285758862 0.0426011575661 -3780.87309981 0.128580559015 0.725501943148 0.0224232500497 0.104527759608 0.977514119639 2.21135177454e-09 0.244285758862 0.0426011575661 -3780.87309981 0.137305306271 0.696705158166 0.0223658339511 0.0855245750445 0.973073683495 2.20340858708e-09 0.282871549476 0.0460773023046 -3779.04391909 0.137305306271 0.696705158166 0.0223658339511 0.0855245750445 0.973073683495 2.20340858708e-09 0.282871549476 0.0460773023046 -3779.04391909 0.137243536116 0.69503395726 0.0224614818516 0.0792628034488 0.97306092368 2.17243037374e-09 0.284105638081 0.0464971525384 -3778.91171105 0.137243536116 0.69503395726 0.0224614818516 0.0792628034488 0.97306092368 2.17243037374e-09 0.284105638081 0.0464971525384 -3778.91171105 # ...
Now let's generate some plots to visualize the chain:
#!bash
postprocess -o plots -p example_a examples/example_a.ini --extra examples/extra_a.py
You will get some summary statistics, as with the other samplers in demos three, five and nine, for example:
#!text Marginalized mean, std-dev: cosmological_parameters--ommh2 = 0.137262 ± 0.00391691 cosmological_parameters--h0 = 0.693811 ± 0.0195768 cosmological_parameters--ombh2 = 0.0224055 ± 0.00045473 cosmological_parameters--tau = 0.0856337 ± 0.0128334 cosmological_parameters--n_s = 0.969108 ± 0.011548 cosmological_parameters--a_s = 2.1983e-09 ± 6.39297e-11 cosmological_parameters--omega_m = 0.286191 ± 0.0230251 cosmological_parameters--omega_b = 0.0466207 ± 0.00212351 like = -3781.27 ± 1.52141 # etc. ...
You will also get a set of convergence statistics based on the criterion from Dunkley et al (2005):
#!text Dunkely et al (2005) power spectrum test. For converged chains j* > 20.0: cosmological_parameters--ommh2 j* = 25.1 cosmological_parameters--h0 j* = 28.3 cosmological_parameters--ombh2 j* = 29.2 cosmological_parameters--tau j* = 32.8 cosmological_parameters--n_s j* = 32.3 cosmological_parameters--a_s j* = 28.4 The power spectra for this chain suggest good convergence.
Understanding
The sampling
The sampler in this case uses the classic Metropolis-Hastings algorithm, just like CosmoMC. We use a multivariate Gaussian proposal based on a covariance matrix. When your covmat is a good representation of the likelihood space then MH is very effective. Here's how we configure the sampler in examples/example_a.ini:
#!ini [runtime] sampler=metropolis [metropolis] samples=3200 nstep=10 covmat=examples/covmat_a.txt random_start=F
We could have started off our chain in some randomized position in the prior space by setting random_start to T. This is useful when using several chains to compare convergence and burn-in.
Our pipeline has three modules. See the earlier demos for an overview of what modules do and how we specify them here.
#!ini [pipeline] modules = consistency camb wmap values = examples/values_a.ini likelihoods = wmap9 extra_output = cosmological_parameters/omega_m cosmological_parameters/omega_b
we want to sample over Omega_m h^2 and Omega_v h^2 , but our CAMB interface wants the parameters Omega_c, Omega_b, and Omega_lambda. We use the consistency module to fix this. Consistency knows a large number of relations between cosmological parameters, for example that Omega_c + Omega_b = Omega_m. It takes the parameters you provide and applies these rules to deduce all the others. If it can't figure them all out then it tries again, this time adding more assumptions about Omega_nu and Omega_k. It checks for both under- and over-constrained models.
Since we will want to make plots of them later we ask for the Omega_m and Omega_b derived parameters to be included in the output file.
The plotting
We use Kernel Density Estimation to make the 1D and 2D constraint plots. This is a smoothing/interpolation method that adds together Gaussians (in a normalized space) laid on top of each of the MCMC samples.
You can experiment with this smoothing by varying the --factor-kde parameter (the default is 2.0).
By using "--extra" we get extra post-processing steps from the file examples/extra_a.py. This is how we ask for the scatter plot - see demo nine for more details.
Updated