Wilkinson Microwave Anisotropy Probe

The data made available through this page has been updated. The most recent version of this data may be accessed through /product/map/current/

Product Name
Gibbs-based low-l likelihood


Code and data from Eriksen et al. to estimate the low-l TT likelihood using a Blackwell-Rao estimator in conjunction with pre-computed Gibbs samples. This code provides a faster exact likelihood implementation than the brute-force N_side=16 code used in the WMAP v2p2 likelihood code.

We provide an F90 module that implements the Blackwell-Rao estimator based on 4000 Gibbs samples obtained from the three-year WMAP V-band data. Usage information is documented in the file README_br_mod.txt, reproduced below. The necessary data files are:

comm_WMAP3_Vband_fgcorr_sigma.fits -- Sky spectrum samples required for the BR estimator (WMAP3; V1 and V2)

comm_WMAP3_Vband_fgcorr_cls.fits -- Power spectrum samples (WMAP3; V1 and V2)

gibbs_WMAP3_Vband_lmax50.tar.gz -- Sky signal samples (WMAP3; V1 and V2)

Here is a listing of the usage information from the file README_br_mod.txt:

USAGE OF BR_MOD_DIST.F90 Hopefully, the Blackwell-Rao module should be straightforward to use: 1) Include the module at the beginning of your program: use br_mod_dist 2) Read the sigma_l sky samples from a FITS file: integer(i4b) :: unit, lmax, numchain, numsamples real(sp), dimension(:,:,:), pointer :: sigmas gibbs_file = unit = call read_gibbs_chain(gibbs_file, unit, lmax, & & numchain, numsamples, sigmas) You now have the samples stored in a three-dimensional array: sigma(0:lmax, 1:numchain, 1:numsamples) Note that samples from a typical Commander run are stored as sigma_l * l*(l+1)/2pi in muK^2. 3) Initialize the BR module (that is, provide the data you want to include). In this case, we want to use only the four first chains, reject the 50 first samples from each chain, and include only 2 <= l <= lmax from each sample: lmin = 2 call initialize_br_mod(lmin, sigma(2:lmax, & & 1:4, 51:numsamples) 4) Compute the estimator for whatever number of C_l spectra you want. Note that there is an offset to be subtracted from each term in the BR estimator in order to avoid overflow errors, and this is based on the first spectrum. Therefore, the first case should not be too different from the best-fit spectrum, although any reasonable spectrum will usually do the job. The call is performed as follows, compute_br_estimator(cls, lnL) The power spectrum must be passed in the same units as the sigma_l's, typically C_l * l*(l+1)/2pi in muK^2 for Commander based results. 5) Please report all bugs to the author. (But not necessarily every problem with making it compile on some obscure platform.. :-)) SPECIFIC NOTES FOR WMAP3 ANALYSIS This module may be used to replace the WMAP pixel-based likelihood code. The benefits are primarily in computational scaling; the WMAP code scales as Npix^3, while Gibbs sampling (for WMAP) scales as Npix^1.5, but produces the same result. Further, the expensive Gibbs sampling step is only done once (and the basic data products for WMAP3 are provided by us), and the computational cost for computing lnL with the Blackwell-Rao estimator is then basically zero. The required changes to the WMAP3 code are the following: 1) Include br_mod_dist.f90 in the Makefile 2) Set lowl_max in WMAP_3yr_options.f90 to, say, 30. 3) Include a 'use br_tools_mod' statement in WMAP_3yr_likelihood.f90 4) Modify the sub-routine PASS2_WMAP_init: a) Comment out 'if(use_lowl_TT)call setup_for_TT_exact(lowl_max)' b) Initialize the BR module as described above, using l's well above the highest l-max you want to include. If you want lowl_max=30, set, say, lmax_br = 50. Remember to reject a few early samples from the beginning of each chain, although this is not very important for the chains provided from the WMAP3 analysis. 5) Modify PASS2_COMPUTE_LIKELIHOOD a) Define the lowl_cl array to be double-precision. b) Comment out 'call compute_like(...)' c) Allocate the lowl_cl array, "allocate(lowl_cl(2:50))", and fill it with the current power spectrum between l=2 and 30. Set l=31 to 50 to a fiducial model that does not change during the run. d) Call the compute_br_estimator routine: call compute_br_estimator(lowl_cl(2:50), lnL) and set the WMAP likelihood value to minus lnL: Like(ttlowlllike) = -lnL 7) Set tt_hi_l_start = 31. 6) Recompile. (If using CosmoMC, remember to recompile both the wmap_likelihood and source directories.) 7) And finally (and most importantly): RUN THE MODIFIED CODE FOR A CASE PUBLISHED IN THE LITERATURE, AND CHECK THAT YOU GET THE SAME RESULT! :-) Note: There may be changes I've forgotten here, so if you are unable to obtain sensible results with the above prescription, please let me know, so I can update it. Hans Kristian Eriksen h.k.k.eriksen@astro.uio.no July 31, 2006

Back to Product Page

A service of the HEASARC and of the Astrophysics Science Division at NASA/GSFC
Goddard Space Flight Center, National Aeronautics and Space Administration
HEASARC Director: Dr. Alan P. Smale
LAMBDA Director: Dr. Eric R. Switzer
NASA Official: Dr. Eric R. Switzer
Web Curator: Mr. Michael R. Greason