---
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<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);   
}
"
```

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

