In the last post, we used the scipy minimize function to find the best-fitting model. However, the best-fitting model is not very useful if we don't know the uncertainties attached to it.
What we really need is a probability for each possible model which could describe our Universe. The probability for model $Y$ with parameters x should reflect the probability that we would observe the BOSS dataset if the Universe were described by $Y(x)$. Such a probability distribution is called a likelihood and a very powerful way to get such a likelihood is a Monte Carlo Markov Chain (MCMC).
What is an MCMC? An MCMC is a chain of parameter sets which has the nice property that the appearance of each parameter in the chain also represents its likelihood distribution. This means one can take the chain, produce a histogram of parameter x, which directly provides the (marginalised) likelihood for parameter x (marginalised just means that you integrate over the likelihood of all other parameters).
There are many algorithms which can produce such a chains, but the most basic and well-known is the Metropolis-Hastings algorithm. Here we will just focus on the Python implementation of such a procedure, but if you have never encountered MCMC it is certainly worth learning about it.
We will use the PocoMC package. First, we have to define a likelihood function, which can make use of the calc_chi2() function we defined earlier
def get_loglike(parameters, data, templates, func): return -0.5*calc_chi2(parameters, data, templates, func)
import pickle def runMCMC(boss_data): labels = ['b', 'A1', 'A2', 'A3', 'A4', 'A5', 'alpha', 'sigmaNL'] if os.path.isfile("samples.pickle"): with open('samples.pickle', 'rb') as handle: samples = pickle.load(handle) return samples, labels k, pk, pk_nw = get_Pk(0.3) pk_nw_interp = interp1d(k, pk_nw) krange = np.arange(0.001, 0.5, 0.001) os_model = get_oscillation(krange, k, pk, pk_nw) prior = pc.Prior([ uniform(loc=0.01, scale=4.0), # [0.01, 4] uniform(loc=-10.0, scale=20.0), # [-10, 10] uniform(loc=-100.0, scale=200.0), # [-100, 100] uniform(loc=-10000.0, scale=20000.0), # [-10000, 10000] uniform(loc=-100000.0, scale=200000.0), # [-100000, 100000] uniform(loc=-100000.0, scale=200000.0), # [-100000, 100000] norm(loc=1.0, scale=0.2), # Normal(1.0, 0.2) uniform(loc=0.1, scale=20.0), # [0.1, 20] ]) # Extra safety: ensure prior.n_dim is a Python int if hasattr(prior, "n_dim"): prior.n_dim = int(prior.n_dim) # Initialise sampler sampler = pc.Sampler( prior=prior, likelihood=get_loglike, likelihood_args=(boss_data, {'noBAO': pk_nw_interp, 'os_model': os_model}, get_shifted_model), random_state=0, n_effective=1024, n_active=512, ) # Start sampling sampler.run() samples, logl, logp = sampler.posterior(resample=True) with open('samples.pickle', 'wb') as handle: pickle.dump(samples, handle, protocol=pickle.HIGHEST_PROTOCOL) return samples, labels
import corner def plot_chain(samples, labels=None, filename="corner.png"): samples = np.asarray(samples) if samples.ndim != 2: raise ValueError("samples must have shape (n_samples, n_dim)") if labels is None: labels = [f"p{i}" for i in range(samples.shape[1])] mean = samples.mean(axis=0) std = samples.std(axis=0, ddof=1) fig = corner.corner( samples, labels=labels, quantiles=[0.16, 0.5, 0.84], show_titles=True, title_fmt=".3f", ) fig.savefig(filename, dpi=150, bbox_inches="tight") plt.show() return fig
 |
Figure 2: 2D likelihood distributions for the 7 parameters in our fit. The 1D marginalised likelihood is given on the right together with the $68\%$ confidence levels. |
Now what does that mean in terms of cosmological parameters? Well as mentioned before, the only parameter of interest here is the scaling parameter $\alpha$. This parameter is defined as
\[
\alpha = \frac{D_V}{D^{\rm fid}_V}
\]where the the distance $D_V$ is defined as
\[
D_V = \left[(1+z)^2D^2_A\frac{cz}{H(z)}\right]^{1/3}
\]The superscribed $^{\rm fid}$ just means that this is the same quantity in the fiducial cosmology, which is the one we defined in the cosmo class in nbodykit.
The distance $D_V(z)$ is an average of the angular diameter distance $D_A(z)$ and the Hubble parameter $H(z)$
\[
\begin{align}
D_A(z) &= \frac{c}{1+z}\int^z_0\frac{dz'}{H(z')}\\
H(z) &= H_0\left[(\Omega_cdm+\Omega_b)(1+z)^3 + \Omega_{\Lambda}\right]^{1/2}
\end{align}
\]With this we can now turn the constraint on $\alpha$ into a constraint on $D_V$
from astropy.cosmology import Planck15 import astropy Hz = Planck15.H(zmean) DC = (1.+zmean)*Planck15.angular_diameter_distance(zmean) c_km_s = astropy.constants.c/1000. DVfid = ( DC**2*(zmean*c_km_s/Hz) )**(1./3.) print("DV = ", mean[6]*DVfid, "+/-", std[6]*DVfid)
>> DV = 2154.4702196107396 +/- 56.825103044911536
The next step would be to combine this constraint with measurements from the Cosmic Microwave Background (CMB) to get constraints on actual cosmological parameters. But for now, I will leave it here.
Thanks a lot for making it to the end of this series. If you are interested in going further, you should consider a few extensions to this analysis:
- Have a look at density field reconstruction
- Consider an anisotropic analysis to constrain $D_A(z)$ and $H(z)$ independently
- Do you think the polynomial model we use is the best choice? Can you come up with a better model?
- Use the simulated BOSS catalogues to test whether our procedure is robust
The final code for this analysis is available on GitHub.
cheers
Florian
