--- title: "Workshop on statistical challenges in astronomy -- Hierarchical models in Stan" author: "John Ormerod" output: html_document: null pdf_document: default fontsize: 12pt geometry: margin=1in toc: yes --- # Example 3: Bayesian normal model for cosmological parameter inference from type Ia supernova data in R using Stan. 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. ## Data 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$ ## Statistical model 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. # Running the Model in R using Stan 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. ```{r,echo=TRUE} library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) ``` ```{r,echo=TRUE} # 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) ``` Specify the model ```{r,echo=TRUE} 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 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 om; // dark matter energy density real alpha; // stretch coefficient real beta; // color coefficient real Mint; // intrinsic magnitude real deltaM; real sigint; // magnitude dispersion real 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); } " ``` Fit the model in stan usin 3 cores. ```{r,echo=TRUE} # run MCMC fit <- stan(model_code = stan_model, data = stan_data, seed = 42, chains = 3, iter =5000, cores= 3, warmup=2500 ) ``` Dsiplay summary statistics of the output. ```{r,echo=TRUE} # Output print(fit,pars=c("om", "Mint","w","alpha","beta","deltaM","sigint"),intervals=c(0.025, 0.975), digits=3) ``` # References 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. * Guy et al., (2007). SALT2: using distant supernovae to improve the use of type Ia supernovae as distance indicators, Astronomy and Astrophysics, 466, 11 * Conley et al. (2011). Supernova Constraints and Systematic Uncertainties from the First 3 Years of the Supernova Legacy Survey. * Betoule et al. (2014). Improved cosmological constraints from a joint analysis of the SDSS-II and SNLS supernova samples. Astronomy and Astrophysics, 568, doi:10.1051/0004-6361/201423413