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

Wiki

Clone wiki

cosmosis / example_a

View History

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.
The chain seems to have converged nicely. You will also get a collection of nice 1D and 2D plots, and the special scatter plot we added with the --extra flag:

scatter_omm_h0_ns.png

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
Samples is the total number of proposals in the chain, and nstep is how often we save samples and report acceptance fractions. The covariance matrix specified by covmat is used for proposing.

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