From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida 

 

you are kindly asked to include the complete citation if you used this material in a publication

Code 10.26 Bayesian normal model for cosmological parameter inference from type Ia supernova data in R using Stan

================================================================================

library(rstan)

# 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
data <- read.table("
https://raw.githubusercontent.com/astrobayes/BMAD/master/data/Section_10p11/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)

# Fit
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
                 )

# Output 
print(fit,pars=c("om", "Mint","w","alpha","beta","deltaM","sigint"),intervals=c(0.025, 0.975), digits=3)

================================================================================

Output on screen:

Inference for Stan model: 92f6251bcff3a6ee78ad0ecef758d250.

3 chains, each with iter=15000; warmup=7500; thin=1;

post-warmup draws per chain=7500, total post-warmup draws=22500.

 

                 mean      se_mean          sd      2.5%        25%         50%       75%       97.5%       n_eff       Rhat

om            0.232           0.001    0.090     0.036       0.173       0.243      0.300        0.380        8406            1

Mint      -19.059           0.000    0.017  -19.093   -19.070    -19.059   -19.048     -19.027      11305            1

w            -0.844            0.002    0.178    -1.232     -0.960      -0.829     -0.710       -0.555        9081            1

alpha       0.119            0.000     0.006     0.106      0.114        0.119      0.123         0.131      18587            1

beta         2.432            0.001     0.071     2.293      2.384        2.432      2.480        2.573       19238            1

deltaM   -0.031            0.000     0.013    -0.055    -0.039       -0.031     -0.022       -0.006       16075           1

sigint      0.159             0.000    0.004      0.151     0.156        0.159       0.162        0.168       20519           1

 

Samples were drawn using NUTS(diag_e) at Fri May 5 15:14:39 2017.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

© 2017 by Emille E. O. Ishida