Steve: A Hierarchical Bayesian Model for Supernova Cosmology

Traditional ${\chi }^{2}$ analyses, such as those 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), and Betoule et al. (2014), minimize the difference in distance modulus between the observed distance modulus attained from Equation (1) and the cosmologically predicted distance modulus. The ${\chi }^{2}$ function minimized is

$$\begin{eqnarray}&&{\chi }^{2}={\left({\mu }_{\mathrm{obs}}-{\mu }_{C}\right)}^{\dagger }{C}_{\mu }^{-1}({\mu }_{\mathrm{obs}}-{\mu }_{C}),\end{eqnarray} \tag{ 2 }$$

where ${C}_{\mu }^{-1}$ is an uncertainty matrix that combines statistical and systematic uncertainties (see Brout et al. 2019 for a review of these uncertainties for the DES supernova analysis). The predicted ${\mu }_{C}$ is defined as

$$\begin{eqnarray}&&{\mu }_{C}=5\mathrm{log}\left[\displaystyle \frac{{d}_{L}}{10\,\mathrm{pc}}\right],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{d}_{L}=(1+z)\displaystyle \frac{c}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&E(z)=\sqrt{{{\rm{\Omega }}}_{m}{\left(1+z^{\prime} \right)}^{3}+{{\rm{\Omega }}}_{k}{\left(1+z^{\prime} \right)}^{2}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}{\left(1+z^{\prime} \right)}^{3(1+w)}},\end{eqnarray} \tag{ 5 }$$

where d_L is the luminosity distance for redshift z given a specific cosmology; H₀ is the current value of Hubble's constant in $\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$; and ${{\rm{\Omega }}}_{m}$, ${{\rm{\Omega }}}_{k}$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ 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. First, formulating a ${\chi }^{2}$ likelihood results in a loss of model flexibility by assuming Gaussian uncertainty. Second, the method of creating a covariance matrix relies on computing partial derivatives; 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. While 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. Third, the computational efficiency is dependent on both creating the global covariance matrix and 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. While 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 1922), is accounted for by correcting the determined ${\mu }_{\mathrm{obs}}$ from the data or, equivalently, adding a distance bias to the ${\mu }_{C}$ 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 ${\chi }^{2}$ 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 ${\chi }^{2}$ likelihood.

Zoom In Zoom Out Reset image size

To avoid the limitations of the traditional approaches, several recent methods have adopted ABC, where supernova samples are forward-modeled in parameter space and compared to observed distributions. Weyant et al. (2013) provided an introduction to ABC methods for supernova cosmology in the context of the SDSS-II results (Sako et al. 2018) and flat ΛCDM cosmology, while Jennings et al. (2016) demonstrated their superABC method on simulated first-season DES 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 modeling directly from these simulations, the data do 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 ${\chi }^{2}$ 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, modeling 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 (m_B, x₁, and c). Jennings et al. (2016) presented two metrics, which are used to measure the distance between the forward-modeled and observed populations and are minimized in fitting. The first metric compares forward-modeled summary statistic populations (denoted the "Tripp" metric), and the second utilizes 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.

Sitting between the traditional model's simplicity and the complexity of forward modeling lies BHM. Hierarchical models utilize 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 interdependence that 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 while still maintaining most of the computational benefits that come from having a tractable likelihood. Mandel et al. (2009, 2011, 2017) constructed a hierarchical model that they applied to supernova light-curve fitting. March et al. (2011) derived a hierarchical model and simplified it by analytically marginalizing over nuisance parameters to provide increased flexibility with reduced uncertainty over the traditional method, but they did not incorporate bias corrections. March et al. (2014) and Karpenka (2015) improved upon this by incorporating redshift-dependent magnitude corrections from Perrett et al. (2010) to remove bias and validate on 100 realizations of SNLS-like simulations. The recent BAHAMAS model (Shariff et al. 2016) builds on this and reanalyzes the JLA data set (using redshift-dependent bias corrections from Betoule et al. 2014) while including extra freedom in the correction factors α and β, finding evidence for redshift dependence on β. Ma et al. (2016) performed a reanalysis of the JLA data set 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 data set (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 validated their model to be free of significant biases using fits to 30 realizations of supernova data sets that are constructed to mimic the UNITY framework. The well-known Bayesian estimation applied to multiple species (BEAMS) methodology from Kunz et al. (2007) has been extended and applied in several works (Hlozek et al. 2012), most recently to include redshift uncertainty for photometric redshift application as zBEAMS (Roberts et al. 2017) and 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 ${\chi }^{2}$ models.

While 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 supernova simulations with either ΛCDM or flat ΛCDM cosmology. However, quantifying the biases on wCDM cosmology simulations with realistic simulations is becoming critically important as precision supernova cosmology comes into its own, and the focus shifts from determination of ${{\rm{\Omega }}}_{m}$ to both ${{\rm{\Omega }}}_{m}$ and w.

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 while incorporating the flexibility of hierarchical models. We reduce our dependence on an assumed scatter model in simulations by not utilizing bias corrections in an effort to provide a valuable cross-check on analysis methodologies that utilize scatter model–dependent bias corrections.

4.1.1. Observables

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, c is the true color of the supernova, while $\hat{c}$ is the observed color, which, as it has a measurement error, is different from c.

For the moment, let us consider a single supernova and its classic summary statistics, m_B, x₁, and c. For convenience, let us define $\eta \equiv \{{m}_{B},{x}_{1},c\}$. A full treatment of the summary statistics would involve determining $p(\hat{y}| \eta )$, where $\hat{y}$ represents the observed light-curve 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 $\hat{\eta }$ along with a 3 × 3 covariance matrix ${{ \mathcal C }}_{\eta }$ describing the uncertainty on $\hat{\eta }$. Using this simplification, our latent variables are given by

$$\begin{eqnarray}&&p(\hat{\eta }| \eta )\sim { \mathcal N }(\hat{\eta }| \eta ,{{ \mathcal C }}_{\eta }).\end{eqnarray} \tag{ 6 }$$

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

4.2.1. Type Ia Population

Unlike many previous formalisms that utilize M_B as a singular number and model magnitude scatter on the apparent magnitude m_B, we incorporate this scatter into the underlying rest-frame population by having a population in absolute magnitude space. This is mathematically equivalent; however, it allows us to model the underlying population and intrinsic scatter distinctly. To denote this difference, we refer to the mean of our absolute magnitude population with $\langle {M}_{B}\rangle$.

In addition to absolute magnitude, the underlying supernova population is also characterized 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 modeled 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; 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 $\langle {x}_{1}^{i}\rangle$ and $\langle {c}^{i}\rangle$. Both the color and stretch means are modeled 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 m_B and c; however, pathological fitting complications required us to simplify our treatment. We parameterize the color skewness ${\alpha }_{c}$ by sampling ${\delta }_{c}={\alpha }_{c}/\sqrt{1+{\alpha }_{c}^{2}}$, which itself is given a uniform prior ${ \mathcal U }(0,0.98)$ that allows ${\alpha }_{c}$ to span positive values in the realm of physical plausibility, as determined from constraints in Scolnic & Kessler (2016). We sample ${\delta }_{c}$ in log space for efficiency in sampling close to zero. The width of the population, represented by the vector $\{{\sigma }_{{M}_{B}},{\sigma }_{{x}_{1}},{\sigma }_{c}\}$, is subject to Cauchy priors with mean zero and width 1, following recommendations from the Stan user guide.

The only constant between survey populations is the absolute magnitude $\langle {M}_{B}\rangle$. We model the color skewness and redshift-dependent means and width of the color and skew distributions individually for each survey. The probability for a supernova to have the true values M_B, x₁, and c given the underlying population is thus given as

$$\begin{eqnarray}\begin{array}{rcl}P({M}_{B},{x}_{1},c,z| \theta ) & = & { \mathcal N }({M}_{B}| \langle {M}_{B}\rangle ,{\sigma }_{{M}_{B}})\,{ \mathcal N }({x}_{1}| \langle {x}_{1}(z)\rangle ,{\sigma }_{{x}_{1}})\\ & & \times \,{{ \mathcal N }}^{\mathrm{skew}}(c| \langle c(z)\rangle ,{\sigma }_{c},{\alpha }_{c}),\end{array}\end{eqnarray} \tag{ 7 }$$

where $\theta =\{\langle {M}_{B}\rangle ,\langle {x}_{1}(z)\rangle ,\langle c(z)\rangle ,{\sigma }_{{m}_{B}},{\sigma }_{{x}_{1}},{\sigma }_{c},{\alpha }_{c}\}$ for legibility.

4.2.2. Outlier Populations

With the fitted summary statistics $\hat{\eta }$ being biased and their uncertainty underreported, we face a significant challenge utilizing 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 utilize the results of computing the bias from both models.

We model the extra dispersion only in color, and we do so by adding independent uncertainty on the color observation $\hat{c}$. We note that extra dispersion in magnitude ${\hat{m}}_{B}$ (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 modeling additional color dispersion. Tests of incorporating extra dispersion on stretch also show that stretch is less biased than color and causes negligible bias in cosmology.

As shown in Kessler et al. (2013), 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 ${\kappa }_{0}+{\kappa }_{1}z$ to our observed color uncertainty (in quadrature). The κ parameters are highly degenerate with the width of the intrinsic color population ${\sigma }_{c}$. We subject them to Cauchy priors centered on zero and with width 0.05, where κ is bounded between 0 and 0.05. 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. (2013). As such, our combined covariance on the observation $\hat{\eta }$ is given by ${{ \mathcal C }}_{\mathrm{tot}}={{ \mathcal C }}_{\eta }+\mathrm{DiagMatrix}\left[0,0,{({\kappa }_{0}+{\kappa }_{1}z)}^{2}\right]$.

Fully covariant extra dispersion on $\{{m}_{B},{x}_{1},c\}$ (rather than just dispersion on c) was also tested by modeling 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 modeled with redshift independence.

4.4.1. Cosmology

4.4.2. Supernova Standardization Parameters

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 SNe Ia shows $\gt 5\sigma$ 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 ${\rm{\Delta }}M=\gamma { \mathcal H }(\mathrm{log}(M)-10))$, where ${ \mathcal H }$ is the Heaviside step function, M 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,

$$\begin{eqnarray}&&{\rm{\Delta }}M=\delta (0)\left[\displaystyle \frac{1.9\left(1-\tfrac{\delta (0)}{\delta (\infty )}\right)}{0.9+{10}^{0.95z}}+\displaystyle \frac{\delta (0)}{\delta (\infty )}\right],\end{eqnarray} \tag{ 8 }$$

where $\delta (0)$ represents the correction at redshift zero, and $\delta (\infty )$ is a parameter allowing the behavior to change with increasing redshift. We take flat priors on $\delta (0)$ and $\delta (0)/\delta (\infty )$. Finally, we assume that the observed mass probability $\hat{m}$ supplied to the model is perfectly determined and thus set $P(\hat{m}| m)=\delta (\hat{m}-m)$.

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 $\tfrac{d{\hat{m}}_{B}}{d{{ \mathcal Z }}_{i}}$, $\tfrac{d{\hat{x}}_{1}}{d{{ \mathcal Z }}_{i}}$, and $\tfrac{d\hat{c}}{d{{ \mathcal Z }}_{i}}$ for each supernova light-curve fit, where ${{ \mathcal Z }}_{i}$ represents the ith 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 $(3\times {N}_{\mathrm{sys}})$ matrix containing the numerical derivatives for the ${N}_{\mathrm{sys}}$ systematics and multiplying it with the row vector containing the offset for each systematic. By scaling the gradient matrix to represent the shift over 1σ 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 the Hubble Space Telescope (HST) Calspec calibration uncertainty, 10 SALT2 model systematic uncertainties, a dust systematic, a global redshift bias systematic, and the systematic peculiar velocity uncertainty. A comprehensive explanation of all systematics is given in Brout et al. (2019); see Table 4 for details. This gives 13 global systematics shared by all surveys that apply globally to all supernova summary statistics plus two systematics per band in each survey. Systematics with known correlations are shifted together to produce covariant deviations, and we thus assume that the numerical derivatives input into our model represent independent systematics. Full details can be found in Brout et al. (2019). With $\eta \equiv \{{m}_{B},{x}_{1},c\}$, our initial conditional likelihood for our observed summary statistics shown in Equation (6) becomes

$$\begin{eqnarray}&&P\left(\hat{\eta },\displaystyle \frac{\partial \hat{\eta }}{\partial {{ \mathcal Z }}_{i}}| \eta ,\delta {{ \mathcal Z }}_{i},{C}_{\eta }\right)={ \mathcal N }\,\left(\hat{\eta }+\delta {{ \mathcal Z }}_{i}\displaystyle \frac{\partial \hat{\eta }}{\partial {{ \mathcal Z }}_{i}}| \eta ,{{ \mathcal C }}_{\eta }\right).\end{eqnarray} \tag{ 9 }$$

4.4.5. Selection Effects

One large difference between traditional 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 possible events we observe are drawn from the distribution predicted by the underlying theoretical model and that those events, given that 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 S, as

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=P(\mathrm{data}| \theta ,S).\end{eqnarray} \tag{ 10 }$$

As our model so far describes components of a basic likelihood $P(\mathrm{data}| \theta )$, and we wish to formulate a function $P(S| \mathrm{data},\theta )$ that describes the chance of an event being successfully observed, we rearrange the likelihood in terms of those functions and find

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data}| \theta )}{\int P(S| D,\theta )P(D| \theta ){dD}},\end{eqnarray} \tag{ 11 }$$

where the denominator represents an integral over all potential data D, and θ represents top-level parameters. In the case that our selection effects are best characterized by latent variables instead of data, we can add them to our formulation, and our likelihood becomes

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=\displaystyle \frac{\int P(S| L,\theta )P(\mathrm{data}| L)P(L| \theta ){dL}}{\iint P(S| L,\theta )P(L| \theta ){dL}},\end{eqnarray} \tag{ 12 }$$

where L represents our latent parameters. This is derived in Appendix A.1. To evaluate the effect of our selection effects, we need to evaluate both the selection effect terms in the numerator and the integral in the denominator. The numerator represents the probability that we caught the supernova and it was selected into the cosmology sample. The integral represents our global selection efficiency at a location in parameter space, rather than the probability of our data being selected into our sample. As θ represents the vector of all top-level model parameters, and L represents a vector of all latent parameters, 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 a 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 x₁ or c), such that the denominator of Equation (12), now denoted d for simplicity, is given as

$$\begin{eqnarray}&&d=\int \left[\int P(S| {m}_{B})P({m}_{B}| z,\theta ){{dm}}_{B}\right]P(S| z)P(z| \theta ){dz},\end{eqnarray} \tag{ 13 }$$

where $P({m}_{B}| z,\theta )$ can be expressed by translating the underlying M_B, x₁, and c population to m_B given the 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 samples the $P(S| z)P(z| \theta )$ distribution reasonably well, and that the survey color and stretch populations can be treated as Gaussian for the purposes of evaluating $P({m}_{B}| z,\theta )$. We found that discarding the color population skewness entirely resulted in highly biased population recovery (see Figure 12 to see the populations), so we instead characterize the skew-normal color distribution with a Gaussian that follows the mean and variance of a skew normal, with a mean given by $\langle c(z)\rangle +\sqrt{\tfrac{2}{\pi }}{\sigma }_{c}{\delta }_{c}$ and variance ${\sigma }_{c}^{2}(1-2{\delta }_{c}^{2}/\pi )$. This shifted Gaussian approximation for color completely removes the unintended bias when simply discarding skewness. This shift was not required for the stretch population, so it was left out of 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 the resulting population recovery can be found in Appendix A.3.

Zoom In Zoom Out Reset image size

The population $P({m}_{B}| z,\theta )$ becomes ${ \mathcal N }({m}_{B}| {m}_{B}^{* }(z),{\sigma }_{{m}_{B}}^{* })$, where

$$\begin{eqnarray}&&{m}_{B}^{* }(z)=\langle {M}_{B}\rangle +\mu (z)-\alpha \langle {x}_{1}(z)\rangle +\beta \langle c(z)\rangle ,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\sigma }_{{m}_{B}}^{* 2}={\sigma }_{{M}_{B}}^{2}+{\left(\alpha {\sigma }_{{x}_{1}}\right)}^{2}+{\left(\beta {\sigma }_{c}\right)}^{2}.\end{eqnarray} \tag{ 15 }$$

What then remains is determining the functional form of $P(S| {m}_{B})$. 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

$$\begin{eqnarray}&&P(S| {m}_{B})={{\rm{\Phi }}}^{c}({m}_{B}| {\mu }_{\mathrm{CDF}},{\sigma }_{\mathrm{CDF}}),\end{eqnarray} \tag{ 16 }$$

where ${{\rm{\Phi }}}^{c}$ is the complementary cumulative distribution function and ${\mu }_{\mathrm{CDF}}$ and ${\sigma }_{\mathrm{CDF}}$ 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 that suffer from saturation and thus rejection of low-redshift supernovae, or for groups of surveys treated together (as is common with low-redshift surveys), we find that a skew normal is a better analytic form, taking the form

$$\begin{eqnarray}&&P(S| {m}_{B})={{ \mathcal N }}^{\mathrm{Skew}}({m}_{B}| {\mu }_{\mathrm{Skew}},{\sigma }_{\mathrm{Skew}},{\alpha }_{\mathrm{Skew}}).\end{eqnarray} \tag{ 17 }$$

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 that 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 ${\chi }^{2}$ of the analytic fit reached 1, 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 parameterized selection effects in our model by including the covariance matrix of selection effect uncertainty. Formally, we include deviations from the determined mean selection function parameters with parameter vector $\delta S$ and apply a normal prior on this parameter as per the determined uncertainty covariance matrix. While 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 Equation (13) and replace it with a correction for each observed supernova. For the error function (denoted with the subscript "CDF") and skew-normal selection functions (denoted with a subscript "Skew"), respectively, the correction per SN Ia becomes

$$\begin{eqnarray}&&{d}_{\mathrm{CDF}}={{\rm{\Phi }}}^{c}\left(\displaystyle \frac{{m}_{B}^{* }-{\mu }_{\mathrm{CDF}}}{\sqrt{{\sigma }_{{m}_{B}}^{* 2}+{\sigma }_{\mathrm{CDF}}^{2}}}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{rcl}{d}_{\mathrm{Skew}} & = & 2{ \mathcal N }\,\left(\displaystyle \frac{{m}_{B}^{* }-{\mu }_{\mathrm{Skew}}}{\sqrt{{\sigma }_{{m}_{B}}^{* 2}+{\sigma }_{\mathrm{Skew}}^{2}}}\right)\\ & & \times \,{\rm{\Phi }}\left(\displaystyle \frac{\mathrm{sign}({\alpha }_{\mathrm{Skew}})({m}_{B}^{* }-{\mu }_{\mathrm{Skew}})}{\tfrac{{\sigma }_{{m}_{B}}^{* 2}+{\sigma }_{\mathrm{Skew}}^{2}}{{\sigma }_{\mathrm{Skew}}^{2}}\sqrt{\tfrac{{\sigma }_{\mathrm{Skew}}^{2}}{{\alpha }_{\mathrm{Skew}}^{2}}+\tfrac{{\sigma }_{{m}_{B}}^{* 2}{\sigma }_{\mathrm{Skew}}^{2}}{{\sigma }_{{m}_{B}}^{* 2}+{\sigma }_{\mathrm{Skew}}^{2}}}}\right),\end{array}\end{eqnarray} \tag{ 19 }$$

and is incorporated into our likelihood. Note that the above efficiencies utilize the common form of the normal distribution rather than the conditional probability notation found previously in this work. This is illustrated in Figure 4. Our corrections for the DES spectroscopic data utilize the CDF functional form, with the combined low-redshift surveys being modeled with the skew-normal efficiency. Further details on this choice are given in Section 5.2.

Zoom In Zoom Out Reset image size

Having laid out each individual aspect of the model, the relationships between variables, and our treatment of uncertainty, here we summarize the relationships in our model mathematically. In this summary of $P(\mathrm{data}| \theta )$, we leave out sample selection for simplicity. Referring to Equation (12), the relationships with $P(\mathrm{data}| \theta )$ are as follows:

$$\begin{eqnarray}&&\begin{array}{l}P(\mathrm{data}| \theta )\\ \propto \,\int d\mu \,{{dm}}_{B}\,{{dx}}_{1}\,{dc}\,{dz},{dm}\ P(\hat{z}| z)P(\hat{m}| m)\\ \times \,P({\hat{m}}_{B},{\hat{x}}_{1},\hat{c},{\hat{{ \mathcal C }}}_{\eta }^{-1}| {m}_{B},{x}_{1},c,\delta {{ \mathcal Z }}_{i},{\kappa }_{0},{\kappa }_{1})\\ \times \,P({x}_{1},c| \langle {x}_{1}^{i}\rangle ,\langle {c}^{i}\rangle ,{\sigma }_{{x}_{1}},{\sigma }_{c},z)\\ \times \,P({m}_{B}| \langle {M}_{B}\rangle ,\mu ,\delta (0),\delta (\infty ),{\sigma }_{{M}_{B}},{\sigma }_{{x}_{1}},{\sigma }_{c},\alpha ,\beta ,z,m)\\ \times \,P(\mu | z,{{\rm{\Omega }}}_{m},w).\end{array}\end{eqnarray} \tag{ 20 }$$

The denominator of Equation (12) is then given by either Equation (18) or (19), depending on the survey; similarly, $P(S| \mathrm{data},\theta )$ is given by Equation (17) or (16), respectively. Combining all of these gives us our full model likelihood with selection effects accounted for.

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 ${ \mathcal U }{(0.0004,0.01)}^{0.5}$, and for the DES-like survey, this is ${ \mathcal U }{(0.008,1.0)}^{0.33}$. For the low-redshift survey, this is equivalent to sampling y = z from 0.02 to $\sqrt{0.1}$, and for the high-redshift survey, this is equivalent to sampling $y={z}^{2}$ from 0.2 to 1.0. These distributions are arbitrary, and this test has been performed with various flat and power-law distributions.
• 2.  
  Draw a random mass probability from ${ \mathcal U }(0,1)$ and calculate the mass-brightness correction using $\delta (0)\,=0.08$, $\delta (0)/\delta (\infty )=0.5$, and Equation (8).
• 3.  
  Draw an absolute magnitude, stretch, and color from the respective distributions ${ \mathcal N }(-19.3,0.1)$, ${ \mathcal N }(0,1)$, and ${ \mathcal N }(0,0.1)$.
• 4.  
  Calculate $\mu (z)$ given the drawn redshift and cosmological parameters ${{\rm{\Omega }}}_{m}=0.3$, $w=-1$ under flat ΛCDM cosmology. Use this to determine the true apparent magnitude of the object m_B using ${m}_{B}=\mu +{M}_{B}-\alpha {x}_{1}+\beta c$.
• 5.  
  Determine if the SN Ia is detected using the detection probability $P(S| {m}_{B})={{ \mathcal N }}^{\mathrm{skew}}(13.72,1.35,5.87)$ for the low-redshift survey (numeric values obtained by fitting to existing low-redshift data). For the DES-like survey, accept with probability $P(S| {m}_{B})={{\rm{\Phi }}}^{C}(23.14,0.5)$. Repeat from step 1 until we have a supernova that passes. We use realistic values for the selection probability to ensure that our model is numerically stable with highly skewed selection functions.
• 6.  
  Add independent, Gaussian observational error onto the true m_B, x₁, and c using Gaussian widths of 0.04, 0.2, and 0.03, respectively (following the mean uncertainty for DES-like SNANA simulations). Add extra color uncertainty in quadrature of ${\kappa }_{0}+{\kappa }_{1}z$, where ${\kappa }_{0}={\kappa }_{1}=0.03$.

The selection function parameters (a skew normal for low-redshift and a complementary error function for high-redshift) are all given an independent uncertainty of 0.01 (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-z 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.

Zoom In Zoom Out Reset image size

We test four models: flat ΛCDM, flat wCDM, ΛCDM, and flat wCDM with a prior ${{\rm{\Omega }}}_{m}\sim { \mathcal N }(0.3,0.01)$, with the latter included to allow sensitive tests on bias for w. To achieve statistical precision, we fit 100 realizations of supernova data sets. Cosmological parameters are recovered without significant bias. Combined posterior surfaces of all 100 realization fits for ΛCDM are shown in Figure 6, and fits for flat wCDM are shown in Figure 7. By utilizing the Stan framework and several efficient parameterizations (discussed further in Appendix B), fits to these simulations of 2000 supernovae take only on order of a single CPU hour to run.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To investigate biases in the model in fine detail, we look for systematic bias in ${{\rm{\Omega }}}_{m}$ in the flat ΛCDM cosmology test and bias in w for the flat wCDM test with strong prior ${{\rm{\Omega }}}_{m}\sim { \mathcal N }(0.3,0.01)$. This allows us to investigate biases without the investigative hindrances of non-Gaussian or truncated posterior surfaces. The strong prior on ${{\rm{\Omega }}}_{m}$ cuts a slice through the traditional "banana" posterior surface in the w–${{\rm{\Omega }}}_{m}$ plane of Figure 7. Without making such a slice, the variation in w is larger due to a shift along the degeneracy direction of the "banana." By focusing the slice at an almost fixed ${{\rm{\Omega }}}_{m}$, we can see the variation in the mean value of w 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 realizations, the precision of our determination of the mean simulation result is approximately a tenth of the quoted scatter (as a degree of non-Gaussianity of our fits will make this relationship inexact). The deviation from truth values is below this threshold, and no significant bias is detected in either the flat ΛCDM or flat wCDM models. With these simple data, we also correctly recover underlying supernova populations, which can be seen in Figure 12.

Table 2.  Cosmological Parameters for Toy Supernova Data

Model	${{\rm{\Omega }}}_{m}\,\langle \mu \rangle ,\langle \sigma \rangle$ (scatter)	$w\,\langle \mu \rangle ,\langle \sigma \rangle$ (scatter)
Flat ΛCDM	$0.301,0.015(0.012)$	⋯
Flat wCDM	⋯	$-1.00,0.042(0.030)$

Note. Cosmological parameters determined from the surfaces of 100 fits to independent realizations of toy supernova data. As described in the main text, each data set comprised 1000 low-redshift and 1000 high-redshift supernovae. For each chain, we record the mean and standard deviation and then show the average mean and average standard deviation in the table. The scatter introduced by simulation variance (the standard deviation of the 100 mean parameter values) is shown in parentheses. Model bias would appear as shifts away from the simulation values of ${{\rm{\Omega }}}_{m}=0.3$, $w=-1$.

Download table as:  ASCIITypeset image

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 data sets that mimic expected observational data in as many ways as possible. To this end, we test using simulations (using the SNANA package) that follow the observational schedule and observing conditions for the DES and low-z surveys, where the low-z sample is based on observations from CfA3 (Hicken et al. 2009a, 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. (2019). The primary differences from the toy data of the previous section are the different underlying color and stretch, inclusion of spectroscopic data and light-curve cuts, and inclusion of intrinsic dispersion models.

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 analyze the effects of two different scatter models via simulations: the G10 and C11 models described in Section 3. The G10 model dispersion has 70% contribution from coherent variation and 30% from chromatic variation, while 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. (2013).

In addition to the improvements in testing multiple scatter models, we also include peculiar velocities for the low-z sample and our full treatment of systematics as detailed in Brout et al. (2019). 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 respective widths of 0.1 and 1; however, the cosmological parameters were not impacted by the choice of the underlying population, and we continue using only the SK16 population for computational efficiency. The selection effects were quantified by comparing all of 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-z 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. Uncertainty on the analytic selection function is incorporated into our fits, mitigating the imperfection of our analytic form by allowing it to vary in our fits.

Zoom In Zoom Out Reset image size

Table 3.  Supernova Population Distributions

Model	$\langle {x}_{1}\rangle$	${\sigma }_{{x}_{1}}$	$\langle c\rangle$	${\sigma }_{c}$
SK16 low-z	0.55 and −1.5	${}_{-1.0}^{+0.45}$ and ${}_{-0.5}^{+0.5}$	−0.055	${}_{-0.023}^{+0.15}$
SK16 DES	0.973	${}_{1.472}^{+0.222}$	−0.054	${}_{0.043}^{+0.101}$

Note. The SK16 low-z stretch distribution is formed as the sum of two bifurcated Gaussians, with the mean and spread of each component given, respectively.

Download table as:  ASCIITypeset image

Each realization of simulated SN Ia light curves contains the SALT2 light-curve fits and redshifts to 128 low-z 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 wCDM model with matter prior.

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

Zoom In Zoom Out Reset image size

Table 4.  Realistic Simulation Standardization Parameters

Model	$\alpha -{\alpha }_{\mathrm{True}}$	$\beta -{\beta }_{\mathrm{True}}$	$\langle {M}_{B}\rangle -\langle {M}_{B}{\rangle }_{\mathrm{True}}$	${\sigma }_{{m}_{B}}^{\mathrm{DES}}$	${\sigma }_{{m}_{B}}^{\mathrm{low}-z}$
G10 Stat + Syst	0.022 [0.009 (0.008)]	0.34 [0.19 (0.18)]	−0.002 [0.028 (0.015)]	0.070 [0.022 (0.018)]	0.073 [0.025 (0.022)]
G10 Stat	0.000 [0.008 (0.008)]	0.33 [0.20 (0.17)]	0.001 [0.016 (0.013)]	0.069 [0.023 (0.019)]	0.072 [0.026 (0.023)]
C11 Stat + Syst	0.002 [0.009 (0.007)]	−0.04 [0.15 (0.13)]	0.014 [0.030 (0.018)]	0.024 [0.016 (0.011)]	0.029 [0.020 (0.014)]
C11 Stat	0.000 [0.008 (0.007)]	−0.05 [0.16 (0.13)]	0.006 [0.016 (0.015)]	0.025 [0.016 (0.012)]	0.027 [0.020 (0.015)]

Note. Standardization parameters and base intrinsic scatter parameter results for the 100 fits to G10 and C11 simulations. We show the average parameter mean and average standard deviation, respectively, with the simulation scatter shown in brackets, such that each cell shows $\langle \mu \rangle [\langle \sigma \rangle (\mathrm{scatter})]$. The width of the intrinsic scatter (${\sigma }_{{m}_{B}}$) does not have an input truth value, as it is determined from the scatter model.

Download table as:  ASCIITypeset image

Table 5 shows that the G10 model is consistent with $w=-1$, while the C11 model shows evidence of bias on w, scattering high. However, their deviation from the truth value represents a shift of approximately $0.5\sigma$ when taking into account the uncertainty on fits to w. The bias is subdominant 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 or 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 5.  Realistic Simulation Determination of w

Model	$w\ \langle \mu \rangle [\langle {\sigma }_{w}\rangle$ (scatter)]	w Bias
G10 Stat + Syst	−0.998 [0.097 (0.073)]	(0.02 ± 0.07)σ
G10 Stat	−1.008 [0.080 (0.068)]	(−0.10 ± 0.08)σ
C11 Stat + Syst	−0.945 [0.098 (0.077)]	(0.55 ± 0.08)σ
C11 Stat	−0.948 [0.079 (0.066)]	(0.65 ± 0.08)σ

Note. Investigating the combined 100 fits to the G10 and C11 simulations, fitting with both statistics only and when including systematics. The quoted value for w represents the average mean of the fits, with the average uncertainty being shown in brackets and the simulation scatter (the standard deviation of the mean of 100 fits) shown in parentheses. The bias significance represents our confidence that the deviation in the mean w away from −1 is not due to statistical fluctuation.

Download table as:  ASCIITypeset image

Table 4 shows a clear difference in both β and ${\sigma }_{{m}_{B}}$ across the G10 and C11 simulations. As expected, the C11 simulations recover a far smaller intrinsic magnitude scatter, giving a result of approximately 0.025 when compared to the result of 0.070 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 $\beta \approx 3.76$ close to the input truth value of 3.8; however, the β recovery for the G10 simulations is biased high, finding $\beta \approx 3.44$ with an input of 3.1. Interestingly, w 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 $\sqrt{100}$ 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 ${\hat{m}}_{B}$, ${\hat{x}}_{1}$, and $\hat{c}$ output from SALT2 light-curve fitting), in addition to an incorrect reported uncertainty on the summary statistics. To confirm this, we run two tests. In the first, we replace the SALT2-fitted ${\hat{m}}_{B}$, ${\hat{x}}_{1}$, and $\hat{c}$ with random numbers drawn from a Gaussian centered on the true SALT2 m_B, x₁, and c values with a covariance as reported by initial light-curve fits. With this test, both the G10 and C11 fits recover $w=-1.00$ exactly. Our second test aims to test our model while allowing biases in the summary statistics not caused by intrinsic scatter through. That is, the first test ascertained that biases in the summary statistics are the cause of cosmological bias. It is thus important to determine the source of those biases, whether they are from the intrinsic scatter model or another aspect of the simulation. To this end, we test a set of 100 simulations generated using an intrinsic dispersion model of only coherent magnitude scatter. We find $w=-1.00$, showing that the source of the biases in the summary statistics is the underlying intrinsic scatter model. From this, the main challenge of improving our methodology is to handle the fact that the 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 supernova data sets 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 standardization parameters and populations.

Difficulty in providing an adequate parameterization for realistic intrinsic dispersion and the simplification of Malmquist bias to only apparent magnitude also lead to biases in the population parameters. To determine if the underlying population mischaracterization was a cause of cosmological bias, we ran fits wherein the underlying population was fixed to the known distributions used for the simulation. These fits did not change the bias in the C11 simulation. We conclude that the biased population recovery is not the cause of cosmological bias. 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 parameterizations of the underlying supernova population.

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

Zoom In Zoom Out Reset image size

Table 6.  Reduced Parameter Correlations with w

Parameter	G10 Stat+Syst	C11 Stat+Syst
${{\rm{\Omega }}}_{m}$	−0.19	−0.21
α	−0.17	−0.20
β	−0.29	−0.23
$\langle {M}_{B}\rangle$	0.68	0.66
$\delta (0)$	0.00	0.00
$\delta (\infty )/\delta (0)$	0.00	0.00
${\sigma }_{{{\rm{m}}}_{B}}^{0}$	0.04	0.07
${\sigma }_{{{\rm{m}}}_{B}}^{1}$	0.23	0.18
${\sigma }_{{x}_{1}}^{0}$	0.04	0.03
${\sigma }_{{x}_{1}}^{1}$	0.05	0.01
${\sigma }_{c}^{0}$	0.01	0.11
${\sigma }_{c}^{1}$	0.08	0.04
${\alpha }_{c}^{0}$	−0.04	0.04
${\alpha }_{c}^{1}$	0.03	0.01
${\kappa }_{c0}^{0}$	−0.10	−0.05
${\kappa }_{c0}^{1}$	−0.20	−0.17
${\kappa }_{c1}^{0}$	−0.05	−0.01
${\kappa }_{c1}^{1}$	−0.01	0.01
$\langle {x}_{1}^{0}\rangle$	−0.01	−0.05
$\langle {x}_{1}^{1}\rangle$	−0.02	0.02
$\langle {x}_{1}^{2}\rangle$	−0.04	−0.04
$\langle {x}_{1}^{3}\rangle$	−0.03	−0.06
$\langle {x}_{1}^{4}\rangle$	−0.06	−0.06
$\langle {x}_{1}^{5}\rangle$	0.04	0.02
$\langle {x}_{1}^{6}\rangle$	0.04	0.04
$\langle {x}_{1}^{7}\rangle$	0.08	0.03
$\langle {c}^{0}\rangle$	−0.05	−0.12
$\langle {c}^{1}\rangle$	0.11	0.03
$\langle {c}^{2}\rangle$	0.11	0.06
$\langle {c}^{3}\rangle$	0.14	0.04
$\langle {c}^{4}\rangle$	−0.11	−0.11
$\langle {c}^{5}\rangle$	−0.15	−0.08
$\langle {c}^{6}\rangle$	−0.12	−0.13
$\langle {c}^{7}\rangle$	−0.12	−0.06
$\delta [{\mathrm{SALT}}_{0}]$	0.05	0.05
$\delta [{\mathrm{SALT}}_{1}]$	−0.01	0.02
$\delta [{\mathrm{SALT}}_{2}]$	−0.10	−0.09
$\delta [{\mathrm{SALT}}_{3}]$	−0.03	−0.03
$\delta [{\mathrm{SALT}}_{4}]$	0.08	0.09
$\delta [{\mathrm{SALT}}_{5}]$	0.01	0.02
$\delta [{\mathrm{SALT}}_{6}]$	0.05	0.07
$\delta [{\mathrm{SALT}}_{7}]$	−0.11	−0.10
$\delta [{\mathrm{SALT}}_{8}]$	0.01	0.02
$\delta [{\mathrm{SALT}}_{9}]$	0.02	0.02
$\delta [{\mathrm{MWE}}_{B-V}]$	0.03	0.02
$\delta [{\text{}}{HST}\ \mathrm{Calib}]$	−0.07	−0.07
$\delta [{v}_{\mathrm{pec}}]$	0.00	−0.01
$\delta [\delta z]$	0.01	0.00
$\delta [{\rm{\Delta }}g]$	0.05	0.11
$\delta [{\rm{\Delta }}r]$	0.16	0.10
$\delta [{\rm{\Delta }}i]$	−0.16	−0.18
$\delta [{\rm{\Delta }}z]$	−0.26	−0.26
$\delta [{\rm{\Delta }}{\lambda }_{g}]$	0.16	0.20
$\delta [{\rm{\Delta }}{\lambda }_{r}]$	0.05	0.06
$\delta [{\rm{\Delta }}{\lambda }_{i}]$	0.00	−0.01
$\delta [{\rm{\Delta }}{\lambda }_{z}]$	0.09	0.07

Note. Correlations determined from the combined 100 simulation fits. Correlations for the low-z band systematics and the latent parameters representing selection function uncertainty are not shown but have negligible correlation. Zero superscripts indicate the DES, and a superscript 1 represents the low-z survey.

Download table as:  ASCIITypeset image

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

With the increased flexibility of BHMs over traditional models, we expect to find an increased uncertainty on parameter inference. This increased uncertainty is one of the strengths of hierarchical models, as it represents a more thorough accounting of model uncertainty. To characterize the influence of the extra degrees of freedom in our model, we analyze the uncertainty on w 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 on w is ±0.08. 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.

Table 7.  w Error Budget

Feature	Parameters	${\sigma }_{w}$	Cumulative
Cosmology	${{\rm{\Omega }}}_{m}$, w	0.051	0.051
Standardization	α, β, $\langle {M}_{B}\rangle$, $\delta (0)$, $\delta (\infty )/\delta (0)$	0.046	0.068
Intrinsic scatter	${\kappa }_{0}$, ${\kappa }_{1}$	0.020	0.071
Redshift-independent populations	${\sigma }_{{M}_{B}}$, ${\sigma }_{c}$, ${\sigma }_{{x}_{1}}$, ${\alpha }_{c}$	0.022	0.074
Redshift-dependent populations	$\langle {c}_{i}\rangle$, $\langle {x}_{1,i}\rangle$	0.030	0.080
Systematics	$\delta {{ \mathcal Z }}_{i}$, $\delta S$	0.054	0.096

Note. Error budget determined from analyzing uncertainty on simulation data while progressively enabling model features. We start from the top of the table, only varying cosmological parameters ${{\rm{\Omega }}}_{m}$ and w, and then progressively unlock parameters and let them fit as we progress down the table. The cumulative uncertainty shows the total uncertainty on w on the fit for all, where the ${\sigma }_{w}$ term is derived by taking the quadrature difference in cumulative uncertainty as we progress.

Download table as:  ASCIITypeset image

We compare the results of our model against those of the BBC+CosmoMC method (Kessler & Scolnic 2017). This method has been used in prior analyses, such as the Pantheon sample analysis of Scolnic et al. (2018), 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 that were also run on the same set of 200 validation simulations and compare the recovered w 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.

Zoom In Zoom Out Reset image size

Table 8.  w Bias Comparison

Model	G10	C11	(G10 + C11)
Steve $\langle w\rangle$	−0.998 ± 0.007	−0.945 ± 0.007	−0.972 ± 0.006
BBC $\langle w\rangle$	−1.044 ± 0.006	−0.978 ± 0.007	−1.010 ± 0.005a
Δ $\langle w\rangle$	+0.044 ± 0.006	+0.033 ± 0.006	+0.038 ± 0.004
Δ ${\sigma }_{w}$	0.057 ± 0.004	0.062 ± 0.004	0.060 ± 0.003

Note. We characterize the bias on w using the 100 simulations for the G10 scatter model and 100 simulations for the C11 scatter model. We also show the results when combining the G10 and C11 models into a combined set of 200 simulations. The mean w value for our method and BBC are presented, along with the mean when averaging the difference between our method and BBC for each individual simulation. Averages are computed giving each simulation sample the same weight. In the model, Δ represents Steve – BBC. The final row shows the scatter between Steve and BBC for the different simulations.

^aThis value is computed with each simulation having the same weight. It disagrees with the value provided in Brout et al. (2019, Table 10, row 3), which uses inverse variance–weighted averages. We do not utilize this weight because the variance is correlated with the value of w due to the ${{\rm{\Omega }}}_{m}$ prior applied in the fitting process. We note that if inverse variance weighting is applied to both data sets, they both shift by ${\rm{\Delta }}w\approx 0.005$, and thus the predicted difference between the BBC method and Steve remains the same.

Download table as:  ASCIITypeset image

As shown in Brout et al. (2019), 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 C11. As such, we expect the BBC method to have a w 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 the BBC method and Steve are sensitive to the intrinsic scatter model, finding differences of ∼0.066 and 0.053, respectively, in w when varying the scatter model. The BBC method finds w biased low for G10 and high for C11 (by about ±0.03), so taking the average result only results in a small bias of −0.01 in w. Our method shows a small improvement in the insensitivity to the intrinsic scatter model (having a decrease in difference in w between the G10 and C11 models), finding no bias for G10 but a w biased high for C11. This decrease in error is not statistically significant, as we have statistical uncertainty of ∼0.01 for 100 simulation realizations. The average bias over the two scatter models is +0.028, representing a larger bias than the BBC method.

When comparing both the G10 and C11 sets of simulations independently, our model differs from BBC in its average prediction of w by +0.044 and +0.033, 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 w when run on the DES-SN3YR sample. When using Planck priors, our uncertainty on w is reduced compared to using our simulation Gaussian prior on ${{\rm{\Omega }}}_{m}$, shrinking the average w difference from 0.06 to 0.04. After factoring this into our uncertainty, we expect our BHM method to, on average, recover ${w}_{\mathrm{BHM}}\,={w}_{\mathrm{BBC}}+0.04\pm 0.04$.

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

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

$$\begin{eqnarray}&&P(\theta | \mathrm{data})\propto P(\mathrm{data}| \theta )P(\theta ),\end{eqnarray} \tag{ 21 }$$

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 $P(\theta )$. Of primary interest to us is the likelihood of observing the data given our parameterized model, ${ \mathcal L }\equiv P(\mathrm{data}| \theta )$. 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 that 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 S, as

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=P(\mathrm{data}| \theta ,S).\end{eqnarray} \tag{ 22 }$$

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 characterize a function $P(S| \mathrm{data},\theta )$, 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 utilizing that $P(x,y,z)=P(x| y,z)P(y,z)=P(y| x,z)P(x,z)$:

$$\begin{eqnarray}&&P(\mathrm{data}| S,\theta )P(S,\theta )=P(S| \mathrm{data},\theta )P(\mathrm{data},\theta )\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&P(\mathrm{data}| S,\theta )=\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data},\theta )}{P(S,\theta )}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data}| \theta )P(\theta )}{P(S| \theta )P(\theta )}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data}| \theta )}{P(S| \theta )},\end{eqnarray} \tag{ 26 }$$

which is equal to the likelihood ${ \mathcal L }$. At this point, our derivation depends on whether our selection effects are best modeled as a function of observables or latent parameters. In most cases, selection effects will be best modeled directly as a function of observables, and we would introduce an integral over all possible events D, so we can evaluate $P(S| \theta )$,

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data}| \theta )}{\int P(S,D| \theta ){dD}},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=\displaystyle \frac{P(S| \mathrm{data},\theta )P(\mathrm{data}| \theta )}{\int P(S| D,\theta )P(D| \theta ){dD}}.\end{eqnarray} \tag{ 28 }$$

The second option, wherein selection effects can be more computationally efficiently modeled with the inclusion of latent parameters (for example, in the case where we do not have direct access to the observables upon which our data selection is determined), we can introduce our latent parameters in addition to an integral over all possible data,

$$\begin{eqnarray}&&{ \mathcal L }(\theta ;\mathrm{data})=\displaystyle \frac{\int P(S| \mathrm{data},L,\theta )P(\mathrm{data}| L)P(L| \theta ){dL}}{\iint P(S| D,L,\theta )P(D| L,\theta )P(L| \theta ){dD}\,{dL}},\end{eqnarray} \tag{ 29 }$$

where L represents our latent parameters that model the true underlying values of our observables, such that our data are conditioned directly on them. In this formulation, the selection effects can depend on both data and latent variables.

To turn the generalized equations from the previous sections into selection effects relevant for this model, we need to highlight that θ represents only our top-level parameters (${{\rm{\Omega }}}_{m}$, w, α, β, etc.), and that our parameterization of the true underlying values (for example, m_B) takes the form of latent parameters. We thus continue from Equation (29) and write out the denominator d:

$$\begin{eqnarray}&&d=\iint P(S| D,L,\theta )P(D| L,\theta )P(L| \theta ){dD}\,{dL}.\end{eqnarray} \tag{ 30 }$$

In our formulation, we assume that our selection effects can be sufficiently encapsulated by latent parameters. That is, we simplify $P(S| D,L,\theta )\to P(S| L)$. This allows us to isolate our integral $\int P(D| L,\theta )$ in the above equation, integrate it to unity, and come to

$$\begin{eqnarray}&&d=\int P(S| L)P(L| \theta ){dL}.\end{eqnarray} \tag{ 31 }$$

Here we can now move from the generic label L to something specific to our model. Specifically, we assume that our selection effects can be quantified using the true apparent magnitude m_B and redshift z, so $L\to \{{m}_{B},z\}$, or, formally, $P(S| D,L,\theta )=P(S| z)P(S| {m}_{B})$. We do not write out other latent parameters found in the model, as they do not impact selection effects—we would simply find them integrated to unity. Note that this assumption does not capture biases engendered by Poisson noise fluctuations. We advocate that future analyses with higher statistical precision use a precisely determined $P(S| D)$. Writing out the latent variables, our denominator becomes

$$\begin{eqnarray}&&d=\int P(S| z)P(S| {m}_{B})P(z,{m}_{B}| \theta )\,{dz}\,{{dm}}_{B}.\end{eqnarray} \tag{ 32 }$$

We can express $P(z,{m}_{B}| \theta )$ as $P({m}_{B}| z,\theta )P(z| \theta )$, where the first term requires us to calculate the magnitude distribution of our underlying population at a given redshift, and the second term is dependent on survey geometry and supernova rates. We can thus state

$$\begin{eqnarray}&&d=\int \left[\int P(S| {m}_{B})P({m}_{B}| z,\theta ){{dm}}_{B}\right]P(S| z)P(z| \theta ){dz}.\end{eqnarray} \tag{ 33 }$$

By assuming that the distribution $P(S| z)P(z| \theta )$ is well sampled by the observed supernova redshifts, we can approximate the integral over redshift by evaluating

$$\begin{eqnarray}&&\int P(S| {m}_{B})P({m}_{B}| z,\theta ){{dm}}_{B}\end{eqnarray} \tag{ 34 }$$

for each supernova in the data set, i.e., Monte Carlo integration with assumed perfect importance sampling.

As stated in Section 4.4.5, the underlying population in apparent magnitude, when we discard skewness, can be represented as ${ \mathcal N }({m}_{B}| {m}_{B}^{* }(z),{\sigma }_{{m}_{B}}^{* })$, where

$$\begin{eqnarray}&&{m}_{B}^{* }(z)=\langle {M}_{B}\rangle +\mu (z)-\alpha \langle {x}_{1}(z)\rangle +\beta \left(\langle c(z)\rangle +\sqrt{\displaystyle \frac{2}{\pi }}{\sigma }_{c}{\delta }_{c}\right),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\sigma }_{{m}_{B}}^{* }={\sigma }_{{M}_{B}}^{2}+{\left(\alpha {\sigma }_{{x}_{1}}\right)}^{2}+{\left(\beta {\sigma }_{c}\sqrt{1-\displaystyle \frac{2{\delta }_{c}^{2}}{\pi }}\right)}^{2}.\end{eqnarray} \tag{ 36 }$$

Then, modeling $P(S| {m}_{B})$ as either a normal or a skew normal, we can analytically perform the integral in Equation (34) and reach Equations (18) and (19).

In this section, we investigate the effect of approximating the skew normal underlying color distribution as a normal. Specifically, Equations (35) and (36) make the assumption that, for our color distribution, ${{ \mathcal N }}^{\mathrm{Skew}}(\mu ,\sigma ,\alpha )$ is well approximated by ${ \mathcal N }(\mu ,\sigma )$. We sought to improve on this approximation by adjusting the mean and standard deviation of the approximated normal to more accurately describe the actual mean and standard deviation of a skew normal. With $\delta \equiv \alpha /\sqrt{1+{\alpha }^{2}}$, the correct mean and standard deviation are

$$\begin{eqnarray}&&{\mu }_{1}={\mu }_{0}+\sqrt{\displaystyle \frac{2}{\pi }}\delta {\sigma }_{0},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{\sigma }_{1}={\sigma }_{0}\sqrt{1-\displaystyle \frac{2{\delta }^{2}}{\pi }},\end{eqnarray} \tag{ 38 }$$

where we highlight that μ here represents the mean of the distribution, not distance modulus. We can then test the approximation ${{ \mathcal N }}^{\mathrm{Skew}}({\mu }_{0},{\sigma }_{0},\alpha )\to { \mathcal N }({\mu }_{1},{\sigma }_{1})$. Unfortunately, this shift to the mean and standard deviation of the normal approximation, where we treat m_B, x₁, and c as a multivariate skew normal, did not produce stable posterior surfaces. Due to this, we treat the underlying m_B, x₁, and c populations as independent.

We tested a fixed ${\sigma }_{c}$ in the shift correction such that ${\mu }_{1}={\mu }_{0}+\sqrt{2/\pi }\delta k$, where we set k = 0.1 to mirror the width of the input simulation population. This resulted in stable posterior surfaces; however, this introduced recovery bias in several population parameters, so we do not fix ${\sigma }_{c}$. Comparing whether we shift our normal in the approximation or simply discard skewness, Figure 3 shows that the calculated efficiency is significantly discrepant to the actual efficiency if the normal approximation is not shifted. The biases when using shifted or unshifted normal approximations when we fit our model on Gaussian and skewed underlying populations are shown in Figure 12, and only the shifted normal approximation correctly recovers the underlying population parameters.

Zoom In Zoom Out Reset image size