Wilkinson Microwave Anisotropy Probe DR2

The data made available through this page has been updated. The most recent version of this data may be accessed through https://lambda.gsfc.nasa.gov/product/wmap/current/

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