SAMP is an adapted version of Monte Python, a Markov chain Monte Carlo (MCMC) sampler for Bayesian analysis of cosmological data. Its main purpose is to provide reproducible results for the work presented in arXiv:1709.06022 [astro-ph.CO]. The related necessary likelihoods for polynomial and spline fits to supernovae and galaxy differential ages data are included. After installation, a montepython command is available to run the implemented MCMC chains and analyze them (see montepython -h for usage help).
SAMP also allows to use Python libraries as an alternative to installing the CLASS code (which still remains the recommended interface for difficult computations such as power spectra) to evaluate cosmological quantities. For instance, Astropy can be used to compute cosmological distances. Let's take the compressed likelihood for strong lensing time-delays data given in equation (26) of Suyu et al. (2013). This likelihood is simple enough and already available in the original Monte Python (and from SAMP main directory under montepython/likelihoods/timedelay/), where it is interfaced with CLASS to compute distances, so the two implementations can be compared.
| On this page: |
First, let's write the data points. Open a data/timedelay.txt file and add the following text. (File locations are relative to the SAMP main directory.)
# Data on quasar time delays. # # Columns are [zd, zs, lambda_d, mu_d, sigma_d] # with zd = lens redhsift # zs = source redhsift # # Other three params refer to the compressed lognormal PDF of # the time-delay distance. # # B1608+656 (Suyu et al 2010) 0.6304 1.394 4000.0 7.053 0.2282 # RXJ1131-1231 (Suyu et al 2013) 0.295 0.658 1425.6 6.4993 0.19377
Let's create a new folder containing the likelihood (let's call it pytimedelay to distinguish it from the existing CLASS-based implementation timedelay).
$ mkdir montepython/likelihoods/pytimedelay/
The new likelihood folder must contains two files, pytimedelay.data and __init__.py. Let's start with pytimedelay.data, which in our case lists the path to the default data directory and the data file name. It also lists the parameters varied in the chain, namely the Hubble constant H0 (its value is in km/s/Mpc) and the matter density parameter today Omega_m. The prefix is the same as the likelihood directory.
# Timedelay data.
pytimedelay.data_directory = data.path['data']
pytimedelay.fname = 'timedelay.txt'
pytimedelay.use_nuisance = ['H0', 'Omega_m']
Next, the file __init__.py contains the computation of the likelihood. While adding such an involved computation not related to package initialization to an __init__.py file is unusual for Pythonic code, this is the standard way to add a likelihood in Monte Python. More specifically, the file must contain a class with the same name as the likelihood folder and that inherits from the Likelihood Monte Python class. An initialization method __init__() and a loglkl() method are required. The arguments passed to those methods are Monte Python objects and paths that allow to retrieve information about the MCMC state and configuration. Here is the complete code.
"""Compressed time-delay likelihood, eq. 26 of Suyu et al. (2013)
[arXiv:1208.6010]. This module requires Astropy.
"""
import os
from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
import numpy as np
from montepython.likelihood_class import Likelihood
class pytimedelay(Likelihood):
def __init__(self, path, data, command_line):
# Always required by Monte Python.
Likelihood.__init__(self, path, data, command_line)
# Read data points.
datafile = os.path.join(self.data_directory, self.fname)
td_data = np.loadtxt(datafile).transpose()
self.zd, self.zs, self.lambda_d, self.mu_d, self.sigma_d = td_data
def loglkl(self, cosmo, data):
# Current point in MCMC parameter space.
H0 = (data.mcmc_parameters['H0']['current']
* data.mcmc_parameters['H0']['scale'])
Omega_m = (data.mcmc_parameters['Omega_m']['current']
* data.mcmc_parameters['Omega_m']['scale'])
# Call Astropy. Distances are in Mpc.
cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m, Tcmb0=2.725)
Dd = cosmo.angular_diameter_distance(self.zd)
Ds = cosmo.angular_diameter_distance(self.zs)
Dds = (((1. + self.zs) * Ds - (1 + self.zd) * Dd)
/ ( 1. + self.zs))
# Time delay distance dt = Dt / 1 Mpc.
dt = (1 + self.zd) * Dd * Ds / Dds / u.Mpc
# Compute likelihood.
if np.all(np.greater(dt, self.lambda_d)):
lkls = - ((np.log(dt - self.lambda_d) - self.mu_d) ** 2 /
2. / self.sigma_d ** 2)\
- (np.log(np.sqrt(2. * np.pi) * (dt - self.lambda_d)
* self.sigma_d))
lkl = np.sum(lkls)
else:
lkl = data.boundary_loglike
return lkl
The file starts by importing the necessary modules. Note in particular the Astropy class used to compute cosmological quantities for a flat ΛCDM cosmological model. The __init__() method first initializes Monte Python Likelihood class (this is always required). Then we read the data points from the data file specified above.
The loglkl() method recovers the values of the fit parameters for the current point in the MCMC parameter space (rescaled by a value that the user may have specified in the parameter file, see below). Then it computes the angular diameter distances (to the deflector, to the source, and between the deflector and the source) needed for the time delay distance, rescaled to 1 Mpc. All these relations are given in equations (1) and (2) of Suyu et al. (2013), and equation (26) is used to return the logarithm of the likelihood. If a particular point in parameter space leads to an imaginary logarithm, we reject that point assigning a value at the likelihood boundaries. (Note that the data object used to retrieve the latter information and the value of the varying parameters is an object internal to Monte Python, it does not contain the observed data points.)
Once the likelihood is created, install the updated code:
$ python setup.py build $ python setup.py install --user
This creates a command montepython available from the shell. Let's switch to a new directory to run and analyze the chains:
$ mkdir ~/pytimedelay_chains && cd ~/pytimedelay_chains
To run the chains, two files specifying the MCMC configuration are needed. Since CLASS is not required for the new likelihood, the first one is simply an empty file named default.conf:
$ touch default.conf
Then the newly create likelihood name, the cosmological parameters varied in the chains, and the number of accepted steps before writing the chain are specified into a file called pytimedelay.param:
#------Experiments to test (separated with commas)-----
data.experiments=['pytimedelay']
#------ Parameter list -------
# data.parameters[class name] = [mean, min, max, 1-sigma, scale, role]
data.parameters['H0'] = [70., 0., None, 3., 1, 'nuisance']
data.parameters['Omega_m'] = [0.3, 0., 1., 0.1, 1, 'nuisance']
#------ Mcmc parameters ----
data.write_step=5
When declaring the parameters we also have specified initial guesses for the mean, limiting values (i.e., flat priors), standard deviation, re-scale factor and a role. When using CLASS, the latter can take the values 'cosmo', 'derived' (in both cases the parameter is read by CLASS) or 'nuisance' (if the parameter is only varied in the likelihood code and not in CLASS). However, when not relying on CLASS all parameters must be declared as nuisance ones (even though strictly speaking they describe the cosmological model) as a technical implementation detail. (This may change in the future. So far the dependence on CLASS is bypassed in a trivial way through a simple dummy class replacing the classy Python wrapper to CLASS.)
Finally, run the chains and analyze the result:
$ montepython run -p pytimedelay.param -o ./output -N 10000 $ montepython info ./output
This creates several files containing the chains, information about mariginalized minimum credible intervals and plots. The 68% interval should be about (71.5 +/- 2.6) km/s/Mpc for the Hubble constant, while the matter density parameter won't be constrained by the data (that's expected).
The combination of Numpy and Astropy allows to write the likelihood expression in a compact vectorial form without introducing, e.g., for loops. Apart from the more readable code, the main advantage is that, thanks to Numpy optimizations, the execution time is still competitive despite relying on Python modules only (and it will scale nicely increasing the number of data points). In this simple example, requiring only the computation of angular diameter distances, the Astropy implementation is actually much faster than calling CLASS for every point in parameter space (a chain of 10000 points required 23s on an Intel Core i7-2620M CPU @ 2.70GHz, compared to 15m of the original CLASS code-based implementation; acceptance rate, which may slightly influence execution time, was about 0.3 for both cases). Of course, if involved computations such as power spectra or non-standard cosmologies are required, CLASS is the most efficient interface.