Homework 2 (Due Monday, May 2)

The goal of this exercise is to impiment the Metropolis-Hastings algorithm and the Affine-invariant ensemble MCMC for sampling from the basic 6-parameter LCDM model under idealized data.

First download the simulated bandpowers with the following command.

download("https://github.com/EthanAnderes/STA250CMB.jl/homework/homework2/bandpowers.h5", "bandpowers.h5")

Now load the data into Julia.


In [1]:
using HDF5
bandpowers = h5read("bandpowers.h5", "bandpowers")


Out[1]:
5626-element Array{Float64,1}:
    1.37321e-6
 6991.85      
 1734.07      
  680.7       
  311.04      
  283.329     
   60.2673    
  124.896     
   75.3852    
   27.8935    
   53.5414    
   32.1185    
   29.115     
    ⋮         
    1.50434e25
    1.52602e25
    1.57783e25
    1.56306e25
    1.64049e25
    1.67491e25
    1.72643e25
    1.76589e25
    1.83756e25
    1.84843e25
    1.90897e25
    1.9538e25 

These bandpowers were simulated from the Julia lecture STA250CMB/lectures/julia_lecture5_LCDM_MCMC/julia_lecture5.ipynb. You can see that lecture for more details of the simulation procedure but basically the bandpowers are given by $$ \sigma_{\ell} = \frac{1}{2\ell + 1}\sum_{m=-\ell}^\ell |d_{\ell m}|^2. $$ The data spherical harmonic coefficients $d_{\ell m}$ were simulated from the idealized data model $$ d_{\ell m} = T_{\ell m} + \epsilon_{\ell m} $$ where the spectral density of $\epsilon_{\ell m}$ is given by $\sigma^2 \exp\left(\frac{b^2}{8\log 2}\ell(\ell+1)\right)$ with noise and beam given by \begin{align*} \sigma &= \frac{10.0}{3437.75} \\ b &= 0.0035. \end{align*} Note: the beam full width at half max is approximately 12 arcmins wide.

To complete homework 2, impliment Metropolis-Hastings and Affine-invariant ensemble MCMC for sampling from the posterior on LCDM parameters $(\Omega_bh^2, \Omega_ch^2, \tau, \theta_s, 10^9A_s, n_s)$ given these simulated bandpowers. Use a unform prior for the LCDM parameters. You can choose whichever proposal distribution you like (feel free to try Gaussian updates where the covariance matrix of the Gaussian is estimated from the WMAP chain), so long as the resulting chain works mixes well. For both Markov chains, plot the iterations of the samples, the autocorrelation of the chain (after a sufficient burn-in) and make a triangle pairwise-density plot comparing the posterior contours of your chain to the WMAP chain. See STA250CMB/lectures/julia_lecture5_LCDM_MCMC/julia_lecture5.ipynb for details on downloading the WMAP chain and using getDist for plotting.


In [ ]: