Steve

# Steve: A hierarchical Bayesian model for Supernova Cosmology

S. R. Hinton11affiliationmark: , T. M. Davis11affiliationmark: , A. G. Kim22affiliationmark: , D. Brout33affiliationmark: , C. B. D’Andrea33affiliationmark: , R. Kessler44affiliationmark: 55affiliationmark: , J. Lasker44affiliationmark: 55affiliationmark: , C. Lidman66affiliationmark: , E. Macaulay77affiliationmark: , A. Möller88affiliationmark: 66affiliationmark: , M. Sako33affiliationmark: , D. Scolnic55affiliationmark: , M. Smith99affiliationmark: , R. C. Wolf1010affiliationmark: , S. Allam1111affiliationmark: , J. Annis1111affiliationmark: , S. Avila77affiliationmark: , E. Bertin1212affiliationmark: 1313affiliationmark: , D. Brooks1414affiliationmark: , D. L. Burke1515affiliationmark: 1616affiliationmark: , A. Carnero Rosell1717affiliationmark: 1818affiliationmark: , M. Carrasco Kind1919affiliationmark: 2020affiliationmark: , J. Carretero2121affiliationmark: , M. Childress99affiliationmark: , C. E. Cunha1515affiliationmark: , L. N. da Costa1818affiliationmark: 2222affiliationmark: , C. Davis1515affiliationmark: , J. De Vicente1717affiliationmark: , D. L. DePoy2323affiliationmark: , P. Doel1414affiliationmark: , T. F. Eifler2424affiliationmark: 2525affiliationmark: , B. Flaugher1111affiliationmark: , P. Fosalba2626affiliationmark: 2727affiliationmark: , J. Frieman1111affiliationmark: 55affiliationmark: , J. García-Bellido2828affiliationmark: , E. Gaztanaga2626affiliationmark: 2727affiliationmark: , D. W. Gerdes2929affiliationmark: 3030affiliationmark: , R. A. Gruendl1919affiliationmark: 2020affiliationmark: , J. Gschwend1818affiliationmark: 2222affiliationmark: , G. Gutierrez1111affiliationmark: , W. G. Hartley1414affiliationmark: 3131affiliationmark: , D. L. Hollowood3232affiliationmark: , K. Honscheid3333affiliationmark: 3434affiliationmark: , E. Krause2424affiliationmark: , K. Kuehn3535affiliationmark: , N. Kuropatkin1111affiliationmark: , O. Lahav1414affiliationmark: , M. Lima3636affiliationmark: 1818affiliationmark: , M. A. G. Maia1818affiliationmark: 2222affiliationmark: , M. March33affiliationmark: , J. L. Marshall2323affiliationmark: , F. Menanteau1919affiliationmark: 2020affiliationmark: , R. Miquel3737affiliationmark: 2121affiliationmark: , E. Morganson2020affiliationmark: , R. L. C. Ogando1818affiliationmark: 2222affiliationmark: , A. A. Plazas2525affiliationmark: , E. Sanchez1717affiliationmark: , V. Scarpine1111affiliationmark: , R. Schindler1616affiliationmark: , M. Schubnell3030affiliationmark: , S. Serrano2626affiliationmark: 2727affiliationmark: , I. Sevilla-Noarbe1717affiliationmark: , M. Soares-Santos3838affiliationmark: , F. Sobreira3939affiliationmark: 1818affiliationmark: , E. Suchyta4040affiliationmark: , G. Tarle3030affiliationmark: , D. Thomas77affiliationmark: , V. Vikram4141affiliationmark: , and Y. Zhang1111affiliationmark:
(DES Collaboration)
###### Abstract

We present a new Bayesian hierarchical model (BHM) named Steve for performing type Ia supernova (SN Ia) cosmology fits. This advances previous works by including an improved treatment of Malmquist bias, accounting for additional sources of systematic uncertainty, and increasing numerical efficiency. Given light curve fit parameters, redshifts, and host-galaxy masses, we fit Steve simultaneously for parameters describing cosmology, SN Ia populations, and systematic uncertainties. Selection effects are characterised using Monte-Carlo simulations. We demonstrate its implementation by fitting realisations of SN Ia datasets where the SN Ia model closely follows that used in Steve. Next, we validate on more realistic SNANA simulations of SN Ia samples from the Dark Energy Survey and low-redshift surveys (DES Collaboration et al., 2018). These simulated datasets contain more than SNe Ia, which we use to evaluate biases in the recovery of cosmological parameters, specifically the equation-of-state of dark energy, . This is the most rigorous test of a BHM method applied to SN Ia cosmology fitting, and reveals small -biases that depend on the simulated SN Ia properties, in particular the intrinsic SN Ia scatter model. This -bias is less than on average, less than half the statistical uncertainty on . These simulation test results are a concern for BHM cosmology fitting applications on large upcoming surveys, and therefore future development will focus on minimising the sensitivity of Steve to the SN Ia intrinsic scatter model.

cosmology: supernovae

## 1. Introduction

Two decades have passed since the discovery of the accelerating universe (Riess et al., 1998; Perlmutter et al., 1999). Since that time, the number of observed type Ia supernovae (SN Ia) has increased by more than an order of magnitude, with contributions from modern surveys at both low redshift (Bailey et al., 2008; Freedman et al., 2009; Hicken et al., 2009a; Contreras et al., 2010; Conley et al., 2011), and higher redshift (Astier et al., 2006; Wood-Vasey et al., 2007; Frieman et al., 2008; Balland et al., 2009; Amanullah et al., 2010; Chambers et al., 2016; Sako et al., 2018). Cosmological analyses of these supernova samples (Kowalski et al., 2008; Kessler et al., 2009b; Conley et al., 2011; Suzuki et al., 2012; Betoule et al., 2014; Rest et al., 2014; Scolnic et al., 2017) have been combined with complementary probes of large scale structure and the CMB. For a recent review, see Huterer & Shafer (2018). While these efforts have reduced the uncertainty on the equation-of-state of dark energy () by more than a factor of two, it is still consistent with a cosmological constant and the nature of dark energy remains an unsolved mystery.

In attempts to tease out the nature of dark energy, active and planned surveys are continually growing in size and scale. The Dark Energy Survey (DES, Bernstein et al., 2012; Abbott et al., 2016) has discovered thousands of type Ia supernovae, attaining both spectroscopically and photometrically identified samples. The Large Synoptic Survey Telescope (LSST, Ivezic et al., 2008; LSST Science Collaboration et al., 2009) will discover tens of thousands of photometrically classified supernovae. Such increased statistical power demands greater fidelity and flexibility in modelling supernovae for cosmological purposes, as we will require reduced systematic uncertainties to fully utilise these increased statistics (Betoule et al., 2014; Scolnic et al., 2017).

As such, considerable resources are aimed at developing more sophisticated supernova cosmology analyses. The role of simulations mimicking survey observations has become increasingly important in determining biases in cosmological constraints and validating specific supernova models. First used in SNLS (Astier et al., 2006) and ESSENCE analyses (Wood-Vasey et al., 2007), and then refined and improved for SDSS (Kessler et al., 2009b), simulations are a fundamental component of modern supernova cosmology. Betoule et al. (2014) quantise and correct observational bias using simulations, and more recently Scolnic & Kessler (2016) and Kessler & Scolnic (2017) explore simulations to quantify observational bias in SN Ia distances as a function of multiple factors to improve bias correction. Approximate Bayesian computation methods also make use of simulations, trading traditional likelihoods and analytic approximations for more robust models with the cost of increased computational time (Weyant et al., 2013; Jennings et al., 2016). Bayesian Hierarchical models abound (Mandel et al., 2009; March et al., 2011; March et al., 2014; Rubin et al., 2015; Shariff et al., 2016; Roberts et al., 2017), and either use simulation-determined distance-corrections to correct for biases, or attempt to find analytic approximations for effects such as Malmquist bias to model the biases inside the BHM itself.

In this paper, we lay out a new hierarchical model that builds off the past work of Rubin et al. (2015). We include additional sources of systematic uncertainty, including an analytic formulation of selection efficiency which incorporates parametric uncertainty. We also implement a different model of intrinsic dispersion to both incorporate redshift-dependent scatter and to increase numerical efficiency, allowing our model to perform rapid fits to supernovae datasets.

Section 2 is dedicated to a quick review of supernova cosmology analysis methods, and Section 3 outlines some of the common challenges faced by analysis methods. In Section 4 we outline our methodology. Model verification on simulated datasets is given in Section 5, along with details on potential areas of improvement. We summarise our methodology in Section 6.

## 2. Review

Whilst supernova observations take the form of photometric time-series brightness measurements in many bands and a redshift measurement of the supernova (or its assumed host), most analyses do not work from these measurements directly. Instead, most techniques fit an observed redshift and these photometric observations to a supernova model, with the most widely used being that of the empirical SALT2 model (Guy et al., 2007, 2010). This model is trained separately before fitting the supernova light curves for cosmology (Guy et al., 2010; Betoule et al., 2014). The resulting output from the model is, for each supernova, an amplitude (which can be converted into apparent magnitude, ), a stretch term , and color term , along with a covariance matrix describing the uncertainty on these summary statistics (). As all supernovae are not identical, an ensemble of supernovae form a redshift-dependent, observed population of , and , where the hat denotes an observed variable.

This ensemble of , and represents an observed population, which – due to the presence of various selection effects – may not represent the true, underlying supernova population. Accurately characterising this underlying population, its evolution over redshift, and effects from host-galaxy environment, is one of the challenges of supernova cosmology. Given some modelled underlying supernova population that lives in the redshift-dependent space (absolute magnitude of the supernova, traditionally in the Bessell band), , and , the introduction of cosmology into the model is simple – it translates the underlying population from absolute magnitude space into the observed population in apparent magnitude space. Specifically, for any given supernova our map between absolute magnitude and apparent magnitude is given by the distance modulus:

 μobs=mB+αx1−βc−MB+ΔM⋅m+other % corrections, (1)

where is the mean absolute magnitude for all SN Ia given , is the stretch correction (Phillips, 1993; Phillips et al., 1999), and is the color correction (Tripp, 1998) that respectively encapsulate the empirical relation that broader (longer-lasting) and bluer supernovae are brighter. refers to the host-galaxy mass correlation discussed in Section 4.4.3. The distance modulus is a product of our observations, however a distance modulus can be precisely calculated given cosmological parameters and a redshift. The ‘other corrections’ term often includes bias corrections for traditional analyses. Bias corrections can take multiple forms, such as a redshift-dependent function (Betoule et al., 2014) or a 5D function of , , , and (Kessler & Scolnic, 2017; Scolnic et al., 2017).

Traditional analyses such as that found in Riess et al. (1998); Perlmutter et al. (1999); Wood-Vasey et al. (2007); Kowalski et al. (2008); Kessler et al. (2009b); Conley et al. (2011); Betoule et al. (2014), minimise the difference in distance modulus between the observed distance modulus attained from equation 1, and the cosmologically predicted distance modulus. The function minimised is

 χ2=(μobs−μC)†C−1μ(μobs−μC), (2)

where is an uncertainty matrix that combines statistical and systematic uncertainties (see Brout et al. (2018) for a review of these uncertainties for the DES supernova analysis). The predicted is defined as

 μC = 5log[dL10pc], (3) dL = (1+z)cH0∫z0dz′E(z′), (4) E(z) = √Ωm(1+z′)3+Ωk(1+z′)2+ΩΛ(1+z′)3(1+w) (5)

where is the luminosity distance for redshift given a specific cosmology, is the current value of Hubble’s constant in and , , and represent the energy density terms for matter, curvature and dark energy respectively.

The benefit this analysis methodology provides is speed – for samples of hundreds of supernovae or less, efficient matrix inversion algorithms allow the likelihood to be evaluated quickly. The speed comes with several limitations. Firstly, formulating a likelihood results in a loss of model flexibility by assuming Gaussian uncertainty. Secondly, the method of creating a covariance matrix relies on computing partial derivatives and thus any uncertainty estimated from this method loses information about correlation between sources of uncertainty. For example, the underlying supernova color population’s mean and skewness are highly correlated, however this correlation is ignored when determining population uncertainty using numerical derivatives of population permutations. Whilst correlations can be incorporated into a covariance matrix, it requires human awareness of the correlations and thus methods that can automatically capture correlated uncertainties are preferable. Thirdly, the computational efficiency is dependent on both creating the global covariance matrix, and then inverting a covariance matrix with dimensionality linearly proportional to the number of supernovae. As this number increases, the cost of inversion rises quickly, and is not viable for samples with thousands of supernovae. A recent solution to this computational cost problem is to bin the supernovae in redshift bins, which produces a matrix of manageable size and can allow fitting without matrix inversion at every step. Whilst binning data results in some loss of information, recent works tested against simulations show that this loss does not result in significant cosmological biases (Scolnic & Kessler, 2016; Kessler & Scolnic, 2017).

Selection efficiency, such as the well known Malmquist bias (Malmquist K. G., 1922) is accounted for by correcting the determined from the data, or equivalently, adding a distance bias to the prediction. Specifically, Malmquist bias is the result of losing the fainter tail of the supernova population at high redshift. An example of Malmquist bias is illustrated in Figure 1, which simulates supernovae according to equation (1). Simulations following survey observational strategies and geometry are used to calculate the expected bias in distance modulus, which is then subtracted from the observational data. When using traditional fitting methods such as that found in Betoule et al. (2014), these effects are not built into the likelihood and instead are formed by correcting data. This means that the bias uncertainty is not captured fully in the distribution, and subtle correlations between cosmological or population parameters and the bias correction is lost. Recent developments such as the BBC method (Kessler & Scolnic, 2017) incorporate corrections dependent on and , improving their capture of uncertainty on bias corrections in the likelihood.

### 2.2. Approximate Bayesian Computation

To avoid the limitations of the traditional approaches, several recent methods have adopted Approximate Bayesian Computation, where supernova samples are forward modelled in parameter space and compared to observed distributions. Weyant et al. (2013) provides an introduction into ABC methods for supernova cosmology in the context of the SDSS-II results (Sako et al., 2014) and flat CDM cosmology, whilst Jennings et al. (2016) demonstrates their superABC method on simulated first season Dark Energy Survey samples. In both examples, the supernova simulation package SNANA (Kessler et al., 2009a) is used to forward model the data at each point in parameter space.

Simulations provide great flexibility and freedom in how to treat the systematic uncertainties and selection effects associated with supernova surveys. By using forward modelling directly from these simulations, data does not need to be corrected, analytic approximations do not need to be applied, and we are free to incorporate algorithms that simply cannot be expressed in a tractable likelihood such as those found in traditional analyses from Section 2.1. This freedom comes with a cost – computation. The classical method’s most computationally expensive step in a fit is matrix inversion. For ABC methods, we must instead simulate an entire supernova population – drawing from underlying supernova populations, modelling light curves, applying selection effects, fitting light curves and applying data cuts. This is an intensive process.

One final benefit of ABC methods is that they can move past the traditional treatment of supernovae with summary statistics (, , and ). Jennings et al. (2016) presents two metrics, which are used to measure the distance between the forward modelled population and observed population, and are minimised in fitting. The first metric compares forward modelled summary statistic populations (denoted the ‘Tripp’ metric) and the second utilises the observed supernova light curves themselves, moving past summary statistics. However, we note that evaluation of systematic uncertainty was only performed using the Tripp metric.

### 2.3. Bayesian Hierarchical Models

Sitting between the traditional models simplicity and the complexity of forward modelling lies Bayesian hierarchical models (BHM). Hierarchical models utilise multiple layers of connected parameters, with the layers linked via well defined and physically motivated conditional probabilities. For example, an observation of a parameter from a population will be conditioned on the true value of the parameter, which itself will be conditioned on the population distribution of that parameter. We can thus incorporate different population distributions, and parameter inter-dependence which cannot be found in traditional analyses where uncertainty must be encapsulated in a covariance matrix. Unlike ABC methods, which can model arbitrary probability distributions, BHM methods are generally constrained to representing probabilities using analytic forms.

With the introduction of multiple layers in our model, we can add more flexibility than a traditional analysis whilst still maintaining most of the computational benefits that come from having a tractable likelihood. Mandel et al. (2009); Mandel et al. (2011); Mandel et al. (2017) construct a hierarchical model that they apply to supernova light-curve fitting. March et al. (2011) derive a hierarchical model and simplify it by analytically marginalising over nuisance parameters to provide increased flexibility with reduced uncertainty over the traditional method, but do not incorporate bias corrections. March et al. (2014) and Karpenka (2015) improve upon this by incorporating redshift-dependent magnitude corrections from Perrett et al. (2010) to remove bias, and validate on 100 realisations of SNLS-like simulations. The recent BAHAMAS model (Shariff et al., 2016) builds on this and reanalyses the JLA dataset (using redshift dependent bias corrections from Betoule et al., 2014), whilst including extra freedom in the correction factors and , finding evidence for redshift dependence on . Ma et al. (2016) performed a reanalysis of the JLA dataset within a Bayesian formulation, finding significant differences in and values from the original analysis from Betoule et al. (2014). Notably, these methods rely on data that is bias corrected or the methods ignore biases, however the UNITY framework given by Rubin et al. (2015) incorporates selection efficiency analytically in the model, and is applied to the Union 2.1 dataset (Suzuki et al., 2012). The assumption made by the UNITY analysis is that the bias in real data is perfectly described by an analytic function. They validate their model to be free of significant biases using fits to thirty realisations of supernova datasets that are constructed to mimic the UNITY framework. The well known BEAMS (Bayesian estimation applied to multiple species) methodology from Kunz et al. (2007) has been extended and applied in several works (Hlozek et al., 2012), most lately to include redshift uncertainty for photometric redshift application as zBEAMS (Roberts et al., 2017) and to include simulated bias corrections in Kessler & Scolnic (2017). For the latter case, by inferring biases using Bayesian models, sophisticated corrections can be calculated and then applied to more traditional models.

Whilst there are a large number of hierarchical models available, none of them have undergone comprehensive tests using realistic simulations to verify each models’ respective bias. Additionally, testing has generally been performed on supernovae simulations with either CDM cosmology or Flat CDM cosmology. However, quantifying the biases on CDM cosmology simulations with realistic simulations is becoming critically important as precision supernovae cosmology comes into its own, and focus shifts from determination of to both and .

The flexibility afforded by hierarchical models allows for investigations into different treatments of underlying supernova magnitude, color and stretch populations, host-galaxy corrections, and redshift evolution, each of which will be discussed further in the outline of our model below. Our model is designed to increase the numerical efficiency of prior works whilst incorporating the flexibility of hierarchical models. We reduce our dependence on an assumed scatter model in simulations by not utilising bias-corrections in an effort to provide a valuable cross-check on analysis methodologies which utilise scatter-model-dependant bias corrections.

## 3. Challenges in Supernova Cosmology

The diverse approaches and implementations applied to supernova cosmology are a response to the significant challenges and complications faced by analyses. In this Section, we outline several of the most prominent challenges.

Forefront among these challenges is our ignorance of the underlying type Ia intrinsic dispersion. Ideally, analysis of the underlying dispersion would make use of an ensemble of time-series spectroscopy to characterise the diversity of type Ia supernovae. However this data is difficult to obtain, and recent efforts to quantify the dispersion draw inference from photometric measurements. The underlying dispersion model is not a solved problem, and we therefore test against two dispersion models in this work. The first is based on the Guy et al. (2010, hereafter denoted G10) scatter model, the second is sourced from Chotard et al. (2011, hereafter denoted C11). As the SALT2 model does not include full treatment of intrinsic dispersion, each scatter model results in different biases in , , and when fitting the SALT2 model to light curve observations, and results in increased uncertainty on the summary statistics that is not encapsulated in the reported covariance . These two scatter models are currently assumed to span the possible range of scatter in the underlying supernova population. We have insufficient information to prefer one model over the other, and thus we have to account for both possible scatter models.

The underlying supernova population is further complicated by the presence of outliers. Non-type Ia supernovae often trigger transient follow-up in surveys and can easily be mistaken for type Ia supernovae and represent outliers from the standardised SN Ia sample. This contamination is not just a result of non-SN Ia being observed, but can also arise from host galaxy misidentification causing incorrect redshifts being assigned to supernovae. Different optimizations to the host-galaxy algorithm can result in misidentification of the host at the 3% to 9% level (Gupta et al., 2016), resulting in a broad population of outliers. For spectroscopic surveys, where both supernova type and redshift can be confirmed through the supernova spectra, this outlier population is negligible. However, for photometric surveys, which do not have the spectroscopic confirmation, it is one of the largest challenges; how to model, fit, and correct for contaminants.

Finally, one of the other persistent challenges facing supernova cosmology analyses are the high number of systematics. Because of the rarity of SN Ia events, datasets are commonly formed from the SN Ia discoveries of multiple surveys in order to increase the number of supernovae in a dataset. However, each survey introduces additional sources of systematic error, from sources within each survey such as band calibration, to systematics introduced by calibration across surveys. Peculiar velocities, different host environments, and dust extinction represent additional sources of systematic uncertainty which must all be modelled and accounted for.

## 4. Our Method

We construct our hierarchical Bayesian model Steve with several goals in mind: creation of a redshift-dependent underlying supernova population, treatment of an increased number of systematics, and analytic correction of selection effects, including systematic uncertainty on those corrections. We also desire Steve more computationally efficient than prior works, such that cosmological results from thousands of supernovae are obtainable in the order of hours, rather than days. As this is closest to the UNITY method from Rubin et al. (2015, hereafter denoted R15), we follow a similar model scaffold, and construct the model in the programming language Stan (Carpenter et al., 2017; Stan Development Team, 2017). The primary challenge of fitting hierarchical models is their large number of fit parameters, and Stan, which uses automatic differentiation and the no-U-turn Sampler (NUTS, a variant of Hamiltonian Monte Carlo), allows us to efficiently sample high dimensional parameter space.

At the most fundamental level, a supernova cosmology analysis is a mapping from an underlying population onto an observed population, where cosmological parameters are encoded directly in the mapping function. The difficulty arises both in adequately describing the biases in the mapping function, and in adding sufficient, physically motivated flexibility in both the observed and underlying populations whilst not adding too much flexibility, such that model fitting becomes pathological due to increasing parameter degeneracies within the model. In the following sections, we will describe these layers, mapping functions, and occurrences of these fatal pathologies. Summaries of observables and model parameters are shown in Table 1 for easy reference.

### 4.1. Observed Populations

#### 4.1.1 Observables

Like most of the BHM methods introduced previously, we work from the summary statistics, where each observed supernova has a brightness measurement (which is analogous to apparent magnitude), stretch and color , with uncertainty on those values encoded in the covariance matrix . Additionally, each supernova has an observed redshift and a host-galaxy stellar mass associated with it, , where the mass measurement is converted into a probability of being above solar masses. We also have a probability of each supernova being a type Ia, . Our set of observables input into the Steve is therefore given as , as shown in the probabilistic graphical model (PGM) in Figure 2.

As we are focused on the spectroscopically confirmed supernovae for this iteration of the method, we assume the observed redshift is the true redshift such that . Potential sources of redshift error (such as peculiar velocities) are taken into account not via uncertainty on redshift (which is technically challenging to implement) but instead uncertainty on distance modulus. Similarly, we take the mass probability estimate as correct, and do not model a latent variable to represent uncertainty on the probability estimate. One of the strengths of Steve (and the R15 analysis) is that for future data sets where supernovae have been classified photometrically, and we expect some misclassification and misidentification of the host galaxies, these misclassifications can naturally be modelled and taken into account by introducing additional populations that supernovae have a non-zero probability of belonging to.

#### 4.1.2 Latent Variables for Observables

The first layer of hierarchy is the set of true (latent) parameters that describe each supernova. In contrast to the observed parameters, the latent parameters are denoted without a hat. For example, is the true color of the supernova, whilst is the observed color, which, as it has measurement error, is different from .

For the moment, let us consider a single supernova and its classic summary statistics , , . For convenience, let us define . A full treatment of the summary statistics would involve determining , where represents the observed light curves fluxes and uncertainties. However, this is computationally prohibitive as it would require incorporating SALT2 light curve fitting inside our model fitting. Due to this computational expense, we rely on initially fitting the light curve observations to produce a best fit along with a covariance matrix describing the uncertainty on . Using this simplification, our latent variables are given by

 p(^η|η)∼N(^η|η,Cη). (6)

As discussed in Section 3, the SALT2 model does not include full treatment of intrinsic dispersion, and thus this approximation does not encapsulate the full uncertainty introduced from this dispersion.

### 4.2. Underlying Population

#### 4.2.1 Type Ia population

Unlike many previous formalisms which utilise as a singular number and model magnitude scatter on the apparent magnitude , we incorporate this scatter into the underlying rest-frame population by having a population in absolute magnitude space. To denote this difference, we refer to the mean of our absolute magnitude population with .

In addition to absolute magnitude, the underlying supernova population is also characterised by distributions in color and stretch. We follow the prior work of R15 and model the color population as an independent redshift-dependent skew normal distribution for each survey. For the stretch population, we adopt a redshift-dependent normal distribution, and magnitude dispersion is modelled as a normal distribution. We also tested a skew-normal approach for these parameters, reverting to the normal distributions as they are computationally easier to evaluate and we found no reduction in cosmological bias with the skew-normal distributions for stretch and magnitude. Following R15 we allow the mean color and stretch to vary over redshift, anchoring four equally spaced redshift nodes spanning the redshift range of each survey, linearly interpolating between the nodes to determine the mean stretch and color for a given redshift. These nodes are represented as and . Both the color and stretch means are modelled with normal priors. Initial versions of our model adopted a fully covariant multivariate skew normal (with skewness set to zero only for the magnitude component) to capture correlations between and , however pathological fitting complications required us to simplify our treatment. We parametrise the color skewness by sampling which itself is given a uniform prior that allows to span positive values in the realm of physical plausibility as determined from constraints in Scolnic & Kessler (2016). We sample in log space for efficiency in sampling close to zero. The width of the population, represented by the vector is subject to Cauchy priors.

The only constant between survey populations is the absolute magnitude . We model the colour skewness and the redshift-dependent means and width of the colour and skew distributions individually for each survey. The probability for a supernova to have true values , , given the underlying population is thus given as

 P(MB,x1,c,z|θ) = N(MB|⟨MB⟩,σMB)N(x1|⟨x1(z)⟩,σx1) (7) Nskew(c|⟨c(z)⟩,σc,αc),

where for legibility.

#### 4.2.2 Outlier populations

For future use with photometrically classified surveys, we include a simplistic outlier population model. We follow R15 (and therefore Kunz et al., 2007) by implementing a Gaussian mixture model. For surveys with low rates of contamination, it is not possible to fit a contamination population, and the mean of the outlier population has been fixed to the SN Ia population in prior works. However, with the increased number of contaminants expected in the DES photometric sample, we seek a more physically motivated outlier population. We find that an acceptable parametrisation is to model the outlier population with a mean absolute magnitude of , where is constrained to be positive, or even to be greater than a small positive number to reduce degeneracy between the two populations. We note that this represents the mean brightness of outliers, and so outliers could both be brighter and dimmer than the mean SN Ia absolute magnitude. We set the population width to in , and . The probability of each supernova falling into either population is determined by the observed type Ia probability discussed in the previous section. For the spectroscopic survey, we set this to unity. For the photometric proof-of-concept we provide an accurate probability estimate. Further investigation on the effect of inaccurate estimates will be left for future improvements during the analysis of the DES photometric sample.

### 4.3. Correcting biased summary statistics

With the fitted summary statistics being biased and their uncertainty under-reported, we face a significant challenge utilising these statistics naively in supernova cosmology. We must either correct the observables to remove the biases introduced by the intrinsic dispersion of the underlying population, or incorporate this dispersion into our model. We should also avoid assuming a specific dispersion model – either the G10 or C11 model, or utilise the results of computing the bias from both models.

We model the extra dispersion only in color, and do so by adding independent uncertainty on the color observation . We note that extra dispersion in magnitude (from coherent scatter) is absorbed completely by the width of the underlying magnitude population (discussed in Section 4.2.1) without introducing cosmological bias, which is not true of the color term, hence the requirement for modelling additional color dispersion. Tests on incorporating extra dispersion on stretch as well show that stretch is less biased than color, and causes negligible bias in cosmology.

As shown in (Kessler et al., 2013b), the extra color dispersion shows heavy redshift dependence, increasing with redshift. This is an artifact of different filters, however as we may be subject to similar effects in our observational data, we decide to incorporate redshift dependence in our extra uncertainty. We thus add to our observed color uncertainty (in quadrature). The parameters are highly degenerate with the width of the intrinsic color population . We subject them to Cauchy priors centered on zero and with width , where is bounded between and . We pick this maximum value to allow extra dispersion without completely subsuming the intrinsic population widths due to the severe degeneracy, where this maximum value easily encapsulates the determined dispersion according to the results of Kessler et al. (2013b). As such, our combined covariance on the observation is given by .

Fully covariant extra dispersion on (rather than just dispersion on ) was also tested, by modelling the dispersion as a multivariate Gaussian, but it showed negligible improvement in recovering unbiased cosmology over just color dispersion, and was far more computationally inefficient. We note here that we model dispersion in magnitude, but this is done at the level of underlying populations, not observed populations. This magnitude dispersion is modelled with redshift independence.

### 4.4. Mapping function

#### 4.4.1 Cosmology

We formulate our model with three different cosmological parameterisations; Flat CDM, Flat CDM, and standard CDM. is given the prior , was treated with and the equation of state was similarly set to a flat prior . For calculating the distance modulus, we fix . If the Hubble constant has a different value, the absolute magnitude is with the other cosmological parameters unaffected.

#### 4.4.2 Supernova Standardisation Parameters

With increasingly large datasets and more nuanced analyses, the choice of how to handle and becomes an important consideration when constructing a model. R15 employs a broken linear relationship for both color and stretch, where different values of and are adopted depending on whether and are positive or negative (although the cut could be placed at a location other than 0). Shariff et al. (2016) instead model as redshift-dependent, testing two phenomenological models; and a second model which effects a rapid but smooth change in at a turnover redshift .

We tested two models with varying against simulated supernova sets; and . See Section 5.2 for details on simulation generation. We found for both models that non-zero values for are preferred even with constant used in simulation, due to severe degeneracy with selection effects. This degeneracy resulted in a significant bias in recovered cosmology. Due to the recovery of non-zero , we continue to adopt the constant and found in traditional analyses. As such, our calculation of distance modulus mirrors that found in Equation (3).

#### 4.4.3 Host Galaxy Environment

There are numerous results showing statistically significant correlations between host-galaxy environment and supernova properties (Kelly et al., 2010; Lampeitl et al., 2010; Sullivan et al., 2010; D’Andrea et al., 2011; Gupta et al., 2011; Johansson et al., 2013; Rigault et al., 2013). The latest sample of over 1300 spectroscopically confirmed type Ia supernovae show evidence for correlation between host mass and luminosity (Uddin et al., 2017). The traditional correction, as employed in analyses such as Suzuki et al. (2012) and Betoule et al. (2014), invokes a step function such that , where is the Heaviside step function, is the galaxy mass in solar masses and represents the size of the magnitude step. The scale of this step function varies from analysis to analysis, and is treated as a fit parameter. In this work we adopt the model used in R15, which follows the work from Rigault et al. (2013), such that we introduce two parameters to incorporate a redshift-dependent host galaxy mass correction:

 ΔM=δ(0)⎡⎢ ⎢⎣1.9(1−δ(0)δ(∞))0.9+100.95z+δ(0)δ(∞)⎤⎥ ⎥⎦, (8)

where represents the correction at redshift zero, and a parameter allowing the behaviour to change with increasing redshift. We take flat priors on and .

#### 4.4.4 Uncertainty Propagation

The chief difficulty with including systematic uncertainties in supernova analyses is that they have difficult-to-model effects on the output observations. As such, the traditional treatment for systematics is to compute their effect on the supernova summary statistics – computing the numerical derivatives , , , where represents the th systematic.

Assuming that the gradients can be linearly extrapolated – which is a reasonable approximation for modern surveys with high quality control of systematics – we can incorporate into our model a deviation from the observed values by constructing a matrix containing the numerical derivatives for the systematics and multiplying it with the row vector containing the offset for each systematic. By scaling the gradient matrix to represent the shift over of systematic uncertainty, we can simply enforce a unit normal prior on the systematic row vector to increase computational efficiency.

This method of creating a secondary covariance matrix using partial derivatives is used throughout the traditional and BHM analyses. For each survey and band, we have two systematics — the calibration uncertainty and the filter wavelength uncertainty. We include these in our approach, in addition to including HST Calspec calibration uncertainty, ten SALT2 model systematic uncertainties, a dust systematic, a global redshift bias systematic, and also the systematic peculiar velocity uncertainty. A comprehensive explanation of all systematics is given in Brout et al. (2018); see Table 4 for details. This gives thirteen global systematics shared by all surveys, plus two systematics per band in each survey. With , our initial conditional likelihood for our observed summary statistics shown in Equation (6) becomes

 P(^η,∂^η∂Zi|η,δZi,Cη)=N(^η+δZi∂^η∂Zi|η,Cη). (9)

#### 4.4.5 Selection Effects

One large difference between traditional methods and BHM methods is that we treat selection effects by incorporating selection efficiency into our model, rather than relying on simulation-driven data corrections. We describe the probability that the events we observe are drawn from the distribution predicted by the underlying theoretical model and that those events, given they happened, are observed and pass cuts. To make this extra conditional explicit, we can write the likelihood of the data given an underlying model, , and that the data are included in our sample, denoted by , as

 L(θ;data) =P(data|θ,S). (10)

As our model so far describes components of a basic likelihood , and we wish to formulate a function that describes the chance of an event being successfully observed, we rearrange the likelihood in terms of those functions and find

 L(θ;data) =P(S|data,θ)P(data|θ)∫P(S|D,θ)P(D|θ)dD, (11)

where the denominator represents an integral over all potential data . This is derived in Appendix A.1 As represents the vector of all model parameters, and represents a vector of all observed variables, this is not a trivial integral. Techniques to approximate this integral, such as Monte-Carlo integration or high-dimensional Gaussian processes failed to give tractable posterior surfaces that could be sampled efficiently by Hamiltonian Monte-Carlo, and post-fitting importance sampling failed due to high-dimensionality (a brief dismissal of many months of struggle). We therefore simplify the integral and approximate the selection effects from their full expression in all of -space, to apparent magnitude and redshift space independently (not dependent on or ), such that the denominator of equation (11), denoted now for simplicity, is given as

 d=∫[∫P(S|mB)P(mB|z,θ)dmB]P(S|z)P(z|θ)dz, (12)

where can be expressed by translating the underlying , , and population to given cosmological parameters. A full derivation of this can be found in Appendix A.2.

We now apply two further approximations similar to those made in R15 – that the redshift distribution of the observed supernovae reasonably well samples the distribution, and that the survey color and stretch populations can be treated as Gaussian for the purposes of evaluating . We found that discarding the color population skewness entirely resulted in highly biased population recovery (see Figure 12 to see the populations), and so we instead characterise the skew normal color distribution with a Gaussian that follows the mean and variance of a skew normal; with mean given by and variance . This shifted Gaussian approximation for color completely removes the unintended bias when simply discarding skewness. This shift was not required for the stretch population, and so was left out for the stretch population for numerical reasons. The impact of this approximation on the calculated efficiency is shown in Figure 3, and more detail on this shift and resulting population recovery can be found in Appendix A.3.

The population becomes , where

 m∗B(z) = ⟨MB⟩+μ(z)−α⟨x1(z)⟩+β⟨c(z)⟩ (13) σ∗2mB = σ2MB+(ασx1)2+(βσc)2. (14)

What then remains is determining the functional form of . For the treatment of most surveys, we find that the error function, which smoothly transitions from some constant efficiency down to zero, is sufficient. Formally, this gives

 P(S|mB)=Φc(mB|μCDF,σCDF), (15)

where is the complementary cumulative distribution function and and specify the selection function. The appropriateness of an error function has been found by many past surveys (Dilday et al., 2008; Barbary et al., 2010; Perrett et al., 2012; Graur et al., 2013; Rodney et al., 2014). However, for surveys which suffer from saturation and thus rejection of low-redshift supernovae, or for groups of surveys treated together (as is common to do with low-redshift surveys), we find that a skew normal is a better analytic form, taking the form

 P(S|mB)=NSkew(mB|μSkew,σSkew,αSkew). (16)

The selection functions are fit to apparent magnitude efficiency ratios calculated from SNANA simulations, by taking the ratio of supernovae that are observed and passed cuts over the total number of supernovae generated in that apparent magnitude bin. That is, we calculate the probability we would include a particular supernova in our sample, divided by the number of such supernovae in our simulated fields. To take into account the uncertainty introduced by the imperfection of our analytic fit to the efficiency ratio, uncertainty was uniformly added in quadrature to the efficiency ratio data from our simulations until the reduced of the analytic fit reached one, allowing us to extract an uncertainty covariance matrix for our analytic fits to either the error function or the skew normal. This is mathematically identical to fitting the efficiency ratio with a second ‘intrinsic dispersion’ parameter which adds uncertainty to the efficiency ratio data points.

We thus include into our model parametrised selection effects by including the covariance matrix of selection effect uncertainty. Formally, we include deviations from the determined mean selection function parameters with parameter vector , and apply a normal prior on this parameter as per the determined uncertainty covariance matrix. Whilst this uncertainty encapsulates the potential error from the simulations not matching the analytic approximations, it does not cover potential variations of the selection function at the top level — varying cosmology or spectroscopic efficiency. Tests with changing the intrinsic scatter model used in the selection efficiency simulations show that the uncertainty introduced is negligible.

With the well-sampled redshift approximation we can remove the redshift integral in Eq (12) and replace it with a correction for each observed supernova. For the error function (denoted with the subscript ‘CDF’) and skew normal selection functions respectively (denoted with a subscript ‘Skew’), the correction per SN Ia becomes

 dCDF = Φc⎛⎜ ⎜⎝m∗B−μCDF√σ∗2mB+σ2CDF⎞⎟ ⎟⎠ (17) dSkew = 2N⎛⎜ ⎜⎝m∗B−μSkew√σ∗2mB+σ2Skew⎞⎟ ⎟⎠× (18) Φ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝sign(αSkew)(m∗B−μSkew)σ∗2mB+σ2Skewσ2Skew√σ2Skewα2Skew+σ∗2mBσ2Skewσ∗2mB+σ2Skew⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

and is incorporated into our likelihood. This is illustrated in Figure 4. Our corrections for the DES spectroscopic data utilise the CDF functional form, with the combined low redshift surveys being modelled with the skew normal efficiency. Further details on this choice are given in Section 5.2.

## 5. Model Verification

In order to verify our model we run it through stringent tests. First, we validate on toy models, verifying that we recover accurate cosmology when generating toy supernovae data constructed to satisfy the assumptions of the BHM construction. We then validate our model on SNANA simulations based on a collection of low redshift surveys and the DES three-year spectroscopic sample, termed the DES-SN3YR sample

### 5.1. Applied to Toy Spectroscopic Data

We generate simple toy data to validate the basic premise of the model. The data generation algorithm is described below:

1. Draw a redshift from a power law distribution. For the low redshift survey this is , and for the DES-like survey this is .

2. Draw a random mass probability from and calculate the mass-brightness correction using , , and equation (8).

3. Draw an absolute magnitude, stretch and color from the respective distributions , , . Calculate the true absolute magnitude using Equation (1).

4. Calculate given the drawn redshift and cosmological parameters , under Flat CDM cosmology. Use this to determine the true apparent magnitude of the object .

5. Determine if the SN Ia is detected using detection probability for the low redshift survey (numeric values obtained by fitting to existing low redshift data). For the DES-like survey, accept with probability . Repeat from step one until we have a supernova that passes. We use realistic values for the selection probability to ensure our model is numerically stable with highly skewed selection functions.

6. Add independent, Gaussian observational error onto the true , , using Gaussian widths of , , respectively (following the mean uncertainty for DES-like SNANA simulations). Add extra color uncertainty in quadrature of , where .

The selection functions parameters (a skew normal for low-redshift and a complimentary error function for high-redshift) are all given independent uncertainty of (mean and width for the CDF selection function, and mean, width and skewness for the skew normal selection function). Draw from each survey simulation until we have 1000 low- supernovae and 1000 DES-like supernovae, representing a statistical sample of greater power than the estimated 350 supernovae for the DES-SN3YR sample. Sample data for 1000 high and low redshift supernovae are shown in Figure 5, confirming the presence of strong selection effects in both toy surveys, as designed.

We test four models: Flat CDM, Flat CDM, CDM, and Flat CDM with a prior , with the latter included to allow sensitive tests on bias for . To achieve statistical precision, we fit 100 realisations of supernovae datasets. Cosmological parameters are recovered without significant bias. Combined posterior surfaces of all 100 realisations fits for CDM are shown in Figure 6 and fits for Flat CDM are shown in Figure 7. By utilising the Stan framework and several efficient parametrisations (discussed further in Appendix B), fits to these simulations of 2000 supernovae take only on order of a single CPU-hour to run.

To investigate biases in the model in fine detail, we look for systematic bias in in the Flat CDM cosmology test, and bias in for the Flat CDM test with strong prior . This allows us to investigate biases without the investigative hindrances of non-Gaussian or truncated posterior surfaces. The strong prior on cuts a slice through the traditional ‘banana’ posterior surface in the - plane of Figure 7. Without making such a slice, the variation in is larger due to a shift along the degeneracy direction of the ‘banana’. By focusing the slice at an almost fixed , we can see the variation in the mean value of approximately perpendicular to the lines of degeneracy, instead of along them. The results of the analysis are detailed in Table 2, and demonstrate the performance of our model in recovering the true cosmological parameters. As we are using 100 independent realisations, the precision of our determination of the mean simulation result is a tenth of the quoted scatter. The deviation from truth values is below this threshold, no significant bias is detected in either the Flat CDM model or the Flat CDM model. With this simple data, we also correctly recover underlying supernova populations, which can be seen in Figure 12.

### 5.2. DES SN data validation

Many BHM methods have previously been validated on data constructed explicitly to validate the assumptions of the model. This is a useful consistency check that the model implementation is correct, efficient and free of obvious pathologies. However, the real test of a model is its application to realistic datasets that mimic expected observational data in as many possible ways. To this end, we test using simulations (using the SNANA package) that follow the observational schedule and observing conditions for the DES and low- surveys, where the low- sample is based on observations from CfA3 (Hicken et al., 2009a; Hicken et al., 2009b), CfA4 (Hicken et al., 2012) and CSP (Contreras et al., 2010; Folatelli et al., 2010; Stritzinger et al., 2011). Simulation specifics can be found in Kessler et al. (2018).

Additionally, prior analyses often treated intrinsic dispersion simply as scatter in the underlying absolute magnitude of the underlying population (Conley et al., 2011; Betoule et al., 2014), but recent analyses require a more sophisticated approach. In our development of this model and tests of intrinsic dispersion, we analyse the effects of two different scatter models via simulations, the G10 and C11 models described in Section 3. The G10 models dispersion with 70% contribution from coherent variation and 30% from chromatic variation whilst the C11 model has 25% coherent scatter and 75% from chromatic variation. These two broadband scatter models are converted to spectral energy distribution models for use in simulations in Kessler et al. (2013a).

In addition to the improvements in testing multiple scatter models, we also include peculiar velocities for the low- sample, and our full treatment of systematics as detailed in Brout et al. (2018). Our simulated populations are sourced from Scolnic & Kessler (2016, hereafter SK16) and shown in Table 3. Initial tests were also done with a second, Gaussian population with color and stretch populations centered on zero and with respective width and , however cosmological parameters were not impacted by choice of the underlying population and we continue using only the SK16 population for computational efficiency. The selection effects were quantified by comparing all the generated supernovae to those that pass our cuts, as shown in Figure 8. It is from this simulation that our analytic determination of the selection functions for the low- and DES survey are based. We run two simulations to determine the efficiency using the G10 and C11 scatter models and find no difference in the functional form of the Malmquist bias between the two models.

Each realisation of simulated SN Ia light curves contains the SALT2 light-curve fits and redshifts to 128 low- supernovae, and 204 DES-like supernovae, such that the uncertainties found when combining chains is representative of the uncertainty in the DES-SN3YR sample. As our primary focus is Dark Energy, we now focus specifically on the Flat CDM model with matter prior.

Points of maximum posterior for 100 data realisations are shown in Figure 9. The parameter bounds and biases for are listed in Table 5, and secondary parameters are shown in Table 4.

Table 5 shows that the G10 model is consistent with , whilst the C11 model show evidence of bias on , scattering high. However, their deviation from the truth value represents a shift of approximately when taking into account the uncertainty on fits to . The bias is sub-dominant to both the size of the uncertainty for each fit, and the scatter induced by statistical variance in the simulations. We also note that the simulations do not vary cosmological parameters nor population. As our model does include uncertainty on those values, the simulation scatter is expected to be less than the model uncertainty, and represents a minimum bound on permissible uncertainty values.

Table 4 shows a clear difference in both and across the G10 and C11 simulations. As expected, the C11 simulations recover a far smaller intrinsic magnitude scatter, giving results of approximately when compared to the result of for the G10 simulations. The extra smearing of the C11 model does not result in a significantly biased value compared to the average uncertainty on , with recovery of close to the input truth value of , however the recovery for the G10 simulations is biased high, finding with an input of . Interestingly, -bias is only found for the C11 simulations. A measure of the significance of the parameter bias can be calculated by comparing the bias to a tenth of the scatter (as our Monte-Carlo estimate uncertainty is of the scatter). From this, we can see that most biases are detected with high statistical significance due to the large number of simulations tested against.

We investigate the cosmological bias and find its source to be a bias in the observed summary statistics (i.e. the , , output from SALT2 light curve fitting), in addition to incorrect reported uncertainty on the summary statistics. To confirm this, we run two tests. The first of which, we replace the SALT2-fitted , and with random numbers drawn from a Gaussian centered on the true SALT2 , and values with covariance as reported by initial light curve fits. With this test, both the G10 and C11 fits recover exactly. Our second test aims to test our model whilst allowing biases in the summary statistics not caused from intrinsic scatter through. To this end, we test a set of 100 simulations generated using an intrinsic dispersion model of only coherent magnitude scatter, and also find , showing that the source of the biases in summary statistics is the underlying intrinsic scatter model. From this, the main challenge of improving our methodology is to handle the fact that observational uncertainty reported from fitting the SALT2 model to light curves is incorrect, non-Gaussian and biased. Our current model and techniques can quantify the effect of different scatter models on biasing the observed summary statistics, but being unable to constrain the ‘correct’ (simulated) scatter model in our model fit means we cannot fully correct for the bias introduced by an unknown scatter model.

Unfortunately, adding extra fit parameters to allow for shifting observables washes out our ability to constrain cosmology, and applying a specific bias correction requires running a fiducial simulation (assuming cosmology, population and scatter model), which presents difficulties when trying to account for correlations with population and scatter model. This is compounded by the fact that bias corrections do not in general improve fits (increase the log posterior), and so are difficult to fit parametrically. Works such as Kessler & Scolnic (2017) show that bias corrections can be applied to supernovae datasets that can robustly handle multiple intrinsic scatter models, and future work will center on uniting these methodologies — incorporating better bias corrections that separate intrinsic scatter bias and non-Gaussian summary statistic bias from Malmquist bias, without having to precompute standardisation parameters and populations.

Difficulty in providing an adequate parametrisation for realistic intrinsic dispersion, and the simplification of Malmquist bias to only apparent magnitude also leads to biases in the population parameters. Tests when fixing the population parameters to true does not resolve the cosmological bias observed in the C11 simulations. As the population parameters recovered using the simplistic toy supernova data in the previous section do not exhibit significant bias, future work will focus on intrinsic dispersion and Malmquist bias rather than alternate parametrisations of the underlying supernova population.

Table 6 lists the fit correlations between our model fit parameters (excluding the low- band systematics, and Malmquist bias uncertainty parameters which had negligible correlation), showing (in order) cosmological parameters, standardisation parameters, population width and skewness parameters, intrinsic dispersion parameters, mass-step parameters, population mean parameters, SALT2 model systematics, dust systematic, global HST calibration systematic, peculiar velocity systematic, global redshift systematic and DES band magnitude and wavelength systematics. Figure 10 show the full correlations between all non-systematic model parameters. Other interesting correlations are shown and discussed in Figure 10. The band systematics for DES filters , and also show significant correlation with , highlighting the importance of minimising instrumental uncertainty.

For the sample size of the DES low- supernova samples (332 supernova), the bias from intrinsic scatter models is sub-dominant to the statistical uncertainty, as shown in Figure 9. For our full systematics model, the bias represents a deviation between to depending on scatter model, and given that they remain sub-dominant, we will leave more complicated treatment of them for future work.

### 5.3. Uncertainty Analysis

With the increased flexibility of Bayesian hierarchical models over traditional models, we expect to find an increased uncertainty on parameter inference. To characterise the influence of the extra degrees of freedom in our model, we analyse the uncertainty on averaged across 10 nominal simulations of the DES-SN3YR sample with various model parameters allowed to either vary or stay locked to a fixed value. By taking the difference in uncertainty in quadrature, we can infer the relative contribution for each model feature to the uncertainty error budget.

The error budget detailed in Table 7 shows that our uncertainty is dominated by statistical error, as the total statistical uncertainty is on is . With the low number of supernovae in the DES-SN3YR sample, this is expected. We note that the label ‘Systematics’ in Table 7 represents all numerically computed systematics (as discussed in Section 4.4.4) and systematic uncertainty on the selection function.

### 5.4. Methodology Comparison

We compare the results of our model against those of the BBC+CosmoMC method (Kessler & Scolnic, 2017). BBC+CosmoMC has been used in prior analyses, such as the Pantheon sample analysis of Scolnic et al. (2017) and is being used in the primary analysis of the DES-SN3YR sample (DES Collaboration et al., 2018). The BBC method is a two-part process, BBC computes bias corrections for observables, and then the corrected distances are fit using CosmoMC (Lewis & Bridle, 2002). For shorthand, we refer to this combined process as the BBC method hereafter in this paper, as we are concerned with the results of cosmological parameter inference. As a leading supernova cosmology method, it provides a good consistency check as to the current levels of accuracy in recovering cosmological parameters.

To this end, we take the results of the BBC method which were also run on the same set of 200 validation simulations and compare the recovered values to those of our method. The results are detailed in Table 8, and a scatter plot of the simulation results is presented in Figure 11.

As shown in Brout et al. (2018), the BBC method recovers cosmological parameters without bias so long as the intrinsic scatter model is known. As we do not know the correct intrinsic scatter model, the BBC method averages the results when using bias corrections from G10 and from C11. As such, we expect the BBC method to have a -bias in one direction for G10 simulations and the other direction for C11 simulations. These results are consistent with those displayed in Table 8. Both BBC method and Steve are sensitive to the intrinsic scatter model, finding differences of and respectively in when varying the scatter model. The BBC method finds biased low for G10 and biased high for C11 (by about ), so taking the average result only results in a small bias of in . Our method shows a small improvement in the insensitivity to the intrinsic scatter model (having a decrease in difference in between the G10 and C11 models), finding no bias for G10 but a biased high for C11. This decrease in error is not statistically significant as we have statistical uncertainty of for simulation realisations. The average bias over the two scatter models is , representing a larger bias than the BBC method.

When comparing both the G10 and C11 set of simulations independently, our model differs from BBC in its average prediction of by and respectively. For the G10 model this difference is a result of bias in the BBC results, however for the C11 simulations this is a result of both bias from BBC, and a larger bias from our method. These results also allow us to state the expected values for when run on the DES-SN3YR sample. When using Planck priors our uncertainty on is reduced compared to using our simulation Gaussian prior on , shrinking the average -difference from to . After factoring this into our uncertainty, we expect our BHM method to, on average, recover .

Having established that our method exhibits similar shifts in the recovery of compared to BBC, future work will focus on improving the parametrisation of intrinsic scatter model into our framework, with the goal of minimising the effect of the underlying scatter model on the recovery of cosmology.

## 6. Conclusions

In this paper we have outlined the creation of a hierarchical Bayesian model for supernova cosmology. The model takes into account selection effects and their uncertainty, fits underlying populations and standardisation parameters, incorporates unexplained dispersion from intrinsic scatter color smearing and incorporates uncertainty from peculiar velocities, survey calibration, HST calibration, dust, a potential global redshift offset, and SALT2 model uncertainty. Furthermore, our uncertainties in standardisation, population, mass-step and more, being explicitly parametrised in our model, are captured with covariance intact, an improvement on many previous methods. The model has been optimised to allow for hundreds of supernovae to be modelled fully with latent parameters. It runs in under an hour of CPU time and scales linearly with the number of supernovae, as opposed to polynomial complexity of matrix inversion of other methods.

The importance of validating models using high-precision statistics gained by performing fits to hundreds of data realisations cannot be overstated, however this validation is lacking in many earlier BHM models for supernova cosmology. We have validated this model against many realisations of simplistic simulations with well-known and well-defined statistics, and found no cosmological bias. When validating using SNANA simulations, we find evidence of cosmological bias which is traced back to light curve fits reporting biased observables and incorrect covariance. Allowing fully parametrised corrections on observed supernovae summary statistics introduces too many degrees of freedom and is found to make cosmology fits too weak. Allowing simulation based corrections to vary in strength is found to give minor reductions in bias, however the uncertainty on the intrinsic scatter model itself limits the efficacy of the bias corrections. For the data size represented in the DES three-year spectroscopic survey, the determined biases should be sub-dominant to other sources of uncertainty, however this cannot be expected for future analyses with larger datasets. Stricter bias corrections calculated from simulations are required to reduce bias. Ideally, this would include further work on the calculation of intrinsic dispersion of the type Ia supernova population such that we can better characterise this bias.

With our model being validated against hundreds of simulation realisations, representing a combined dataset over more than simulated supernovae, we have been able to accurately determine biases in our model and trace their origin. With the current biases being sub-dominant to the total uncertainty, we now prepare to analyse the DES three-year dataset.

## Acknowledgements

Plots of posterior surfaces and parameter summaries were created with ChainConsumer (Hinton, 2016).

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

## Appendix A Selection Effect Derivation

### a.1. General Selection Effects

When formulating and fitting a model using a constraining dataset, we wish to resolve the posterior surface defined by

 P(θ|data)∝P(data|θ)P(θ), (A1)

which gives the probability of the model parameter values () given the data. Prior knowledge of the allowed values of the model parameters is encapsulated in the prior probability . Of primary interest to us is the likelihood of observing the data given our parametrised model, . When dealing with experiments that have imperfect selection efficiency, our likelihood must take that efficiency into account. We need to describe the probability that the events we observe are both drawn from the distribution predicted by the underlying theoretical model and that those events, given they happened, are subsequently successfully observed. To make this extra conditional explicit, we write the likelihood of the data given an underlying model, , and that the data are included in our sample, denoted by , as:

 L(θ;data)=P(data|θ,S). (A2)

A variety of selection criteria are possible, and in our method we use our data in combination with the proposed model to determine the probability of particular selection criteria. That is, we characterise a function , which colloquially can be stated as the probability of a potential observation passing selection cuts, given our observations and the underlying model. We can introduce this expression in a few lines due to symmetry of joint probabilities and utilising that :

 P(data|S,θ)P(S,θ) = P(S|data,θ)P(data,θ) (A3) P(data|S,θ) = P(S|data,θ)P(data,θ)P(S,θ) (A4) = P(S|data,θ)P(data|θ)P(θ)P