Type Ia supernovae (SNe Ia) are extremely bright transient events which can be used as standardizable candles for distance measurements in cosmological scales. In the late 1990s they provided the first evidence for the current accelerated expansion of the Universe and consequently the existence of dark energy. Since then they have been central to every large-scale astronomical survey aiming to shed light on the dark-energy mystery.
At maximum brightness the observed magnitude of an SN Ia can be connected to its distance modulus \(\mu\) through the expression
\[\begin{equation}\label{eq:10.21} m_{\mbox{obs}} = \mu + M - \alpha x_1 + \beta c, \end{equation}\]where \(m_{\mbox{obs}}\) is the observed magnitude, \(M\) is the magnitude, and \(x_1\) and \(c\) are the stretch and color corrections derived from the SALT2 standardization (Guy et al., 2007), respectively. To take into account the effect of the host stellar mass \(M_\star\) on \(M\) and \(\beta\), we use the correction proposed by Conley et al. (2011): \[ M = \left\{ \begin{array}{ll} M & \mbox{if $M_\star<10^{10}M_{\odot}$ } \\ M + \Delta M & \mbox{otherwise} \end{array} \right. \]
Considering a flat Universe, \(\Omega_k = 0\), containing dark energy and dark matter, the cosmological connection can be expressed as \[ \begin{array}{rl} \displaystyle \mu & = 5 \log\left( \frac{d_L}{10 pc} \right) \\ \displaystyle d_L(z) & = (1 + z) \frac{c}{H_0}\int_0^z E(z)^{-1} dz \\ \displaystyle E(z) & = \sqrt{\Omega_m(1 + z)^3 + (1 - \Omega_m)(1 + z)^{3(1 + w)}} \end{array} \]
where \(d_L\) is the luminosity distance, \(c\) the speed of light, \(H_0\) the Hubble constant, \(\Omega_m\) the dark-energy density, and \(w\) the dark-energy equation of state parameter.
We used data provided by Betoule et al. (2014), known as the joint light-curve analysis (JLA) sample. This is a compilation of data from different surveys which contains 740 high quality spectroscopically confirmed SNe Ia up to redshift \(z \sim 1.0\)
Our statistical model will thus have one response variable (the observer magnitude, mobs) and four explanatory variables (the redshift \(z\), the stretch \(x_1\), the color \(c\), and the host galaxy mass \(M_{\mbox host}\)). The complete statistical model can be expressed as \[ \begin{array}{rl} m_{obs,i} & \sim N(\eta_i,\varepsilon) \\ \eta_i & = 25 + 5\log_{10}(dl_i(H_0,w,\Omega_m)) + M(\Delta M) - \alpha x_{1,i} + \beta c_i \end{array} \]
We use conservative priors over the model parameters. \[ \begin{array}{rl} \varepsilon & \sim \mbox{Gamma}(10^{-3},10^{-3}) \\ M & \sim N(-20,5) \\ \alpha & \sim N(0,1) \\ \beta & \sim N(0,10) \\ \Delta M & \sim N(0,1) \\ \Omega_m & \sim \mbox{Uniform}(0,1) \\ H_0 & \sim N(70,5) \end{array} \]
These are not completely non-informative but they do allow a large range of values to be searched without putting a significant probability on nonphysical values.
After carefully designing the model we must face the extra challenge posed by the integral in the luminosity-distance definition. Fortunately, Stan has a built-in ordinary differential equation (ODE) solver which provides a user-friendly solution to this problem.
Stan’s ODE solver is explained in detail in the manual Team Stan (2016, Section 44). Here we would merely like to call attention for a couple of important points.
The ODE takes as input a function which returns the derivatives at a given time, or, as in our case, a given redshift. This function must have, as input, time (or redshift), system state (the solution at \(t = 0\)), parameters (input a list if you have more than one parameter), real data, and integer data, in that order. In cases where your problem does not require integer or real data, you must still input an empty list.
You must define a zero point for time (redshift, z0 = 0) and the state of the system at that point (in our case, E0 = 0). These must be given as input data in the dictionary to be passed to Stan.
library(rstan)## Warning: package 'rstan' was built under R version 3.4.3## Loading required package: ggplot2## Loading required package: StanHeaders## Warning: package 'StanHeaders' was built under R version 3.4.3## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())# Preparation
# set initial conditions
z0 = 0 # initial redshift
E0 = 0 # integral(1/E) at z0
# physical constants
c = 3e5 # speed of light
H0 = 70 # Hubble constant
# Data: JLA sample, Betoule et al., 2014
# http://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html
#
# 1 response (obsy - observed magnitude)
# 5 explanatory variable (redshift - redshift,
# ObsMag - apparent magnitude,
# x1 - stretch,
# color - color,
# hmass - host mass)
data <- read.table("./jla_lcparams.txt",header=T)
# remove repeated redshift
data2<-data[!duplicated(data$zcmb),]
# prepare data for Stan
nobs <- nrow(data2) # number of SNe
index <- order(data2$zcmb) # sort according to redshift
ObsMag <- data2$mb[index] # apparent magnitude
redshift <- data2$zcmb[index] # redshift
color <- data2$color[index] # color
x1 <- data2$x1[index] # stretch
hmass <- data2$m3rdvar[index] # host mass
stan_data <- list(nobs = nobs,
E0 = array(E0,dim=1),
z0 = z0,
c = c,
H0 = H0,
obs_mag = ObsMag,
redshift = redshift,
x1 = x1,
color = color,
hmass = hmass)stan_model="
functions {
/**
* ODE for the inverse Hubble parameter.
* System State E is 1 dimensional.
* The system has 2 parameters theta = (om, w)
*
* where
*
* om: dark matter energy density
* w: dark energy equation of state parameter
*
* The system redshift derivative is
*
* d.E[1] / d.z =
* 1.0/sqrt(om * pow(1+z,3) + (1-om) * (1+z)^(3 * (1+w)))
*
* @param z redshift at which derivatives are evaluated.
* @param E system state at which derivatives are evaluated.
* @param params parameters for system.
* @param x_r real constants for system (empty).
* @param x_i integer constants for system (empty).
*/
real[] Ez(real z,
real[] H,
real[] params,
real[] x_r,
int[] x_i) {
real dEdz[1];
dEdz[1] = 1.0/sqrt(params[1]*(1+z)^3 + (1-params[1])*(1+z)^(3*(1+params[2])));
return dEdz;
}
}
data {
int<lower=1> nobs; // number of data points
real E0[1]; // integral(1/H) at z=0
real z0; // initial redshift, 0
real c; // speed of light
real H0; // hubble parameter
vector[nobs] obs_mag; // observed magnitude at B max
real x1[nobs]; // stretch
real color[nobs]; // color
real redshift[nobs]; // redshift
real hmass[nobs]; // host mass
}
transformed data {
real x_r[0]; // required by ODE (empty)
int x_i[0];
}
parameters{
real<lower=0, upper=1> om; // dark matter energy density
real alpha; // stretch coefficient
real beta; // color coefficient
real Mint; // intrinsic magnitude
real deltaM;
real<lower=0> sigint; // magnitude dispersion
real<lower=-2, upper=0> w; // dark matter equation of state parameter
}
transformed parameters{
real DC[nobs,1]; // co-moving distance
real pars[2]; // ODE input = (om, w)
vector[nobs] mag; // apparent magnitude
real dl[nobs]; // luminosity distance
real DH; // Hubble distance = c/H0
DH = (c/H0);
pars[1] = om;
pars[2] = w;
// Integral of 1/E(z)
DC = integrate_ode_rk45(Ez, E0, z0, redshift, pars, x_r, x_i);
for (i in 1:nobs) {
dl[i] = DH * (1 + redshift[i])*DC[i, 1];
if (hmass[i] < 10) mag[i] = 25 + 5*log10(dl[i]) + Mint - alpha*x1[i] + beta*color[i];
else mag[i] = 25 + 5*log10(dl[i]) + Mint + deltaM - alpha*x1[i] + beta*color[i];
}
}
model {
// priors and likelihood
sigint ~ gamma(0.001, 0.001);
Mint ~ normal(-20, 5.);
beta ~ normal(0, 10);
alpha ~ normal(0, 1);
deltaM ~ normal(0, 1);
obs_mag ~ normal(mag, sigint);
}
"# run MCMC
fit <- stan(model_code = stan_model,
data = stan_data,
seed = 42,
chains = 3,
iter =15000,
cores= 3,
warmup=7500
)## In file included from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
## from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file21fc437d1f30.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
## # define BOOST_NO_CXX11_RVALUE_REFERENCES
## ^
## <command-line>:0:0: note: this is the location of the previous definition
## In file included from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/base.hpp:28:0,
## from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array.hpp:21,
## from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint/util/multi_array_adaption.hpp:29,
## from C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/numeric/odeint.hpp:61,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:13,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/arr.hpp:38,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/prim/mat.hpp:298,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:12,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file21fc437d1f30.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<N>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:42:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:43:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp: In static member function 'static void boost::multi_array_concepts::detail::idgen_helper<0ull>::call(Array&, const IdxGen&, Call_Type)':
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:53:43: warning: typedef 'index_range' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index_range index_range;
## ^
## C:/Users/John/Documents/R/win-library/3.4/BH/include/boost/multi_array/concept_checks.hpp:54:37: warning: typedef 'index' locally defined but not used [-Wunused-local-typedefs]
## typedef typename Array::index index;
## ^
## In file included from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:42:0,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
## from C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
## from file21fc437d1f30.cpp:8:
## C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
## C:/Users/John/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:17: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
## static void set_zero_all_adjoints() {
## ^# Output
# print(fit,pars=c("om", "Mint","w","alpha","beta","deltaM","sigint"),intervals=c(0.025, 0.975), digits=3)From Bayesian Models for Astrophysical Data by Hilbe, de Souza & Ishida, 2016, Cambridge Univ. Press
Code 10.26 Bayesian normal model for cosmological parameter inference from type Ia supernova data in R using Stan.