{'fraction_non_alphanumeric': 0.036729597948406474, 'fraction_numerical': 0.016905400617502658, 'mean_word_length': 4.6331495913324465, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
UNITY: Confronting Supernova Cosmology's Statistical and Systematic Uncertainties in a Unified Bayesian Framework

While recent supernova cosmology research has benefited from improved measurements, current analysis approaches are not statistically optimal and will prove insufficient for future surveys. This paper discusses the limitations of current supernova cosmological analyses in treating outliers, selection effects, shape- and color-standardization relations, unexplained dispersion, and heterogeneous observations. We present a new Bayesian framework, called UNITY (Unified Nonlinear Inference for Type-Ia cosmologY), that incorporates significant improvements in our ability to confront these effects. We apply the framework to real supernova observations and demonstrate smaller statistical and systematic uncertainties. We verify earlier results that SNe Ia require nonlinear shape and color standardizations, but we now include these nonlinear relations in a statistically well-justified way. This analysis was primarily performed blinded, in that the basic framework was first validated on simulated data before transitioning to real data. We also discuss possible extensions of the method.


INTRODUCTION
Recent supernova (SN) cosmological measurements have greatly reduced both the statistical and systematic uncertainty in our knowledge of the accelerated expansion of the universe (the latest such efforts are presented in Suzuki et al. 2012;Betoule et al. 2014;Rest et al. 2014). Despite these improvements, the frameworks currently in use are not statistically optimal. As we build larger supernova samples, these frameworks will become increasingly inadequate.
This paper offers an improved technique for deriving constraints on cosmological parameters from SN measurements (peak magnitude, light-curve shape and color, host-galaxy mass). Although this paper uses SN cosmology as an example, researchers from other fields may find this type of framework useful. Our work is particularly relevant for researchers confronting partially known uncertainties, selection effects, correlated measurements, and outliers. The rest of this introduction summarizes the problems with existing methods and describes in general terms the basic requirements for greater accuracy. Section 2 describes the framework in detail, and Section 3 quantitatively examines the performance of the method on simulated data. In Section 4, we demonstrate the performance on real SN observations, then conclude in Section 5 with future directions.
1.1. The Current Approach Current cosmological constraints are derived from time series of photometric measurements in multiple bands (light curves) and spectroscopy. Before obtaining SN distances, the light curves must be fit with a model such as SALT2 (Guy et al. 2007(Guy et al. , 2010Mosher et al. 2014). SALT2 models each photometric observation in the observer frame with a combination of a mean SN spectral energy distribution (SED) (scaled by a normalization parameter), the first component of the SED variation (scaled by a shape parameter x 1 ), and the mean color variation in magnitudes (which is scaled by a color parameter c). The template is also shifted in time to match the observations with a date-of-maximum parameter. (In this work, we restrict our attention to the SALT2 empirical model, which is the best validated and most widely used such model.) Guy et al. (2007Guy et al. ( , 2010 and Mosher et al. (2014) trained the SALT2 model in an initial (separate) step before the light-curve fitting, using a dataset with well-measured spectra and light curves that partially overlaps with the SN data used for cosmological distances. 7 After measuring a light-curve shape parameter (x obs 1i ), a light-curve color parameter (c obs i ), and a light-curve normalization (m obs Bi , the rest-frame B-band flux), 8 we construct a distance modulus estimate as α and β are the light-curve shape standardization co- 7 We could, in principle, incorporate the training of SALT2 and the fitting of the light curves into our proposed model, but the necessary modeling of selection effects would be difficult. 8 Of course, the real observables in astronomy are electrons per pixel for each image. We assume here that the photometry, calibration, and light-curve fitting yield an approximate multivariate Gaussian likelihood for each SN. efficient and color standardization coefficient, respectively. They quantify the empirical relations that SNe with broader rest-frame optical light curves or bluer restframe colors are intrinsically more luminous (Phillips 1993;Tripp 1998). δ captures residual luminosity correlations with host-galaxy mass, discussed more in Section 2.2. M B is the absolute magnitude in the rest-frame B-band (for a given H 0 ). All of these coefficients are nuisance parameters for cosmological purposes.
We can then construct a χ 2 to use for cosmological fitting (for illustrative purposes only, we assume here that all of the SN observations are uncorrelated): Bi (α, β, δ) − µ(z i , cosmo)) 2 σ 2 µ obs B (α, β) + σ samp 2 (α, β) + σ ext 2 (z) . (2) σ 2 µ obs B (α, β) captures the measurement uncertainties from the light-curve fit. SALT2 handles k-corrections implicitly, as it is a rest-frame spectral model that fits photometry in the observer frame. Limitations of the SALT2 model are important to add, such as its inability to simultaneously model intrinsic color variation and extinction with one c parameter (discussed more below). While SALT2 does include a simple estimate of the dispersion around its model, including this dispersion does not give a χ 2 per degree of freedom of 1 when performing cosmological fits. σ samp 2 (α, β) is a term that captures this sample-dependent unexplained variance in the SN distribution (sometimes called "intrinsic dispersion"). Finally, σ ext 2 (z) captures dispersion due to gravitational lensing and incoherent peculiar velocities. In current cosmological analyses (those since Kowalski et al. 2008), an iteration is performed between estimating the standardization coefficients, estimating σ samp 2 (α, β), and rejecting outlying SNe. 9

Limitations of Current Work
There are several ways in which SN datasets are imperfect. Outliers, selection effects, nonlinear correlations, partially known uncertainties, and heterogeneity each complicate analyses. The proposed framework will address each of these complications, but first we discuss the limitations of current work.
Obtaining either spectra or very high-quality light curves are the only ways to ensure a transient source is a SN Ia. Even at moderate redshifts, both of these techniques are observationally expensive, and non-Ia SNe will inevitably contaminate the sample. A similar issue arises if SNe are of type Ia, but are peculiar, or the redshift is incorrect. The analysis should thus accommodate some amount of non-Ia SNe, which have dissimi-9 There is a still some amount of confusion in the literature concerning Equation 2, and we address it now to contrast against the proposed framework for, e.g., fitting nonlinear {x 1 , c} standardizations. Finding the α and β values that minimize this χ 2 is not the same as minimizing the dispersion of the Hubble diagram. Nor does it eliminate the correlation between {x 1 , c} and Hubble diagram residuals. The residuals at the best-fit α and β are expected to show a correlation with the observed {x 1 , c}, as uncertainties will preferentially scatter the observed values of {x 1 , c} away from the distribution means and it is only the values without this scatter that should be uncorrelated with Hubble diagram residuals. lar colors, decline rates, and absolute magnitudes. The iterative outlier fit described above converges well for the sorts of contaminating distributions we expect when the samples are relatively pure (Kowalski et al. 2008). When the samples are large (several hundred), or impure ( 85%)-conditions the field is beginning to face-the outliers can dominate the other sources of uncertainty in the fit. Kunz et al. (2007) presented a powerful Bayesian technique for simultaneously modeling the distributions of normal SNe and outliers, but does not confront many complexities of the data, including the luminosity standardizations and selection effects.
Selection effects, the tendency for surveys to select against the faintest SNe, sculpt the observed distribution of SNe (this is frequently referred to as Malmquist bias, Malmquist 1922). If not taken into account, this selection will bias both the cosmological parameter estimation and the standardization coefficients. There are SNto-SN variations in the selection efficiency, even within the same survey and for the same redshift (e.g., due to seeing, host-galaxy contamination, or moonlight). These are deterministic but difficult-to-model sources of variation. Noise also plays a role, as a SN right at the detection threshold may stochastically scatter above or below it. Thus, the detection efficiency transitions to zero smoothly as a function of apparent brightness. In most SN cosmological analyses (based on Equation 2), the treatment of selection effects is performed outside the statistical framework, in the form of survey simulations followed by an ad-hoc redshift-dependent adjustment of the m obs Bi to approximately cancel the estimated bias (Kessler et al. 2009;Conley et al. 2011). Bayesian analyses face a related challenge: selection effects also influence the distribution of light-curve width and color with redshift, which can amplify selection effects if not modeled (Wood-Vasey et al. 2007). However, these frequentist and Bayesian analyses both require knowledge of the true population distribution (as a function of redshift) to get accurate results.
The x 1 and c standardizations in Equation 1 are linear. However, nonlinear decline-rate and color standardizations are statistically justified (Amanullah et al. 2010;Scolnic et al. 2014a,b). 10 Our Union2 result (Amanullah et al. 2010) was derived from subdividing the sample by the best-fit latent variables (see Section 2.1). The subdivisions showed very similar cosmological results, but subdivision tests are statistically weak (statistical uncertainties on the difference are ∼ 2 times larger than the uncertainty on the whole sample). Including these nonlinear standardizations in the fit would be preferable for evaluating their impact and would remove any bias created by the assumption of linearity. Unfortunately, neither the Amanullah et al. (2010), nor the Scolnic et al. (2014b) frameworks were able to incorporate these nonlinear standardizations in a fully selfconsistent way.
Type Ia Hubble diagrams show more dispersion in distance than can be explained with measurement uncertainty alone. As noted above, the framework must in-10 Others have seen a trend towards lower color-standardization coefficients for even redder SNe, c 0.5 (Mandel et al. 2011;Burns et al. 2014), but virtually all the SNe used for cosmology are bluer than this. clude a model of the unexplained dispersion, and should not assume that all of this dispersion is in the independent variable (magnitude). In other words, σ samp 2 is really σ samp 2 (α, β). 11 SN variation (beyond lightcurve width and color), is presently only crudely modeled by current light-curve fitters and especially affects the measured color (Wang et al. 2009;Chotard et al. 2011;Foley & Kasen 2011;Foley 2012;Saunders et al. 2015), resulting in color measurements that mix extinction and SN differences in only partially understood ways. When using Equation 2, there is no way to fit for the unexplained dispersion simultaneously with the other parameters. Improvement is needed here, as the unexplained dispersion interacts with the modeling of selection effects, the outlier rejection, and fitting α and β. Currently, we are constrained to use ad-hoc methods, such as performing many fits with a randomized unexplained-dispersion matrix, and computing the distributions of cosmological parameters fit to fit (Marriner et al. 2011). The impact of these variations is currently subdominant to the total uncertainty (Kessler et al. 2013;Mosher et al. 2014), but this may not be true in future samples and a precise way to evaluate the cosmological impact is desirable. March et al. (2011) took a step in the right direction by using a Bayesian hierarchical model to simultaneously fit for cosmological parameters and unexplained dispersion, but their model cannot accommodate calibration uncertainties, outliers, nonlinear standardizations, or selection effects.
Future SN cosmology analyses will confront additional difficulties, such as heterogeneity. Even very homogeneous SN imaging surveys (such as the Dark Energy Survey, Flaugher 2005; Bernstein et al. 2012) will gain heterogeneity as expensive observational resources, such as near-IR measurements or high-quality spectroscopy, will only be available for a subset of objects. A frequentist, "object-by-object" uncertainty analysis such as Equation 2 cannot make efficient use of this information. 12

Desired Properties for a New Approach
In light of these problems, this paper offers an improved technique incorporating a more sophisticated, Bayesian model of the data. Properly making use of heterogeneous information (like measurements that are only available for a subset of objects) requires a model of the SN population in which the parameters of the distribution are treated as unknowns. In addition, the model of the unexplained dispersion should allow for uncertainty in both size and functional form. We discuss the details of our procedure for parameterizing these possibilities in Section 2.7. Only a Bayesian framework can accommodate more exotic possibilities, such as very large numbers of nuisance parameters (each possibly non-zero and hav-11 The SALT2 light-curve fitter already incorporates a significant amount of model uncertainty in both the color and magnitude measurements. 12 The Union2.1 analysis (Suzuki et al. 2012) already uses a simple Bayesian model for the host-galaxy mass standardization, where samples which lack host-mass measurements are given priors derived from other similar samples, but the remainder of that analysis uses a frequentist framework.
ing initially unknown size), as it can make use of a hierarchical prior around zero. 13 One could imagine many possibilities for these parameters (e.g., color-standardization coefficient, β, as a function of host-galaxy-inclination angle); we do not pursue these in this work, but we do note that they could be built into the general framework.
However, not all of the improvements we propose require Bayesian statistics; some are due to the improvements to the data likelihood (and could thus fit into a frequentist framework). First, we introduce a modified approach to fitting for residual correlations with host-galaxy mass, described in Section 2.2. Second, like Kunz et al. (2007), we handle outliers with a mixture model, described in Section 2.3, which offers improved robustness. Third, to account for selection effects, our proposed framework uses a probabilistic model of selection as described in Section 2.4, modifying the classical likelihood to include selection directly. This cleanly estimates and marginalizes the hyperparameters of the true population distributions simultaneously with other parameters, propagating all uncertainties. Fourth, in Section 2.5, we discuss our approach to fitting for the changing SN independent variable distribution with redshift (population drift) simultaneously with other parameters. Finally, our new framework can accommodate nonlinear standardizations (Section 2.6). A frequentist framework can, in principle, include nonlinear standardizations, but it is computationally challenging. 14 Both the old and new techniques can handle correlations due to systematic uncertainties (such as correlated peculiar velocities and photometric calibration uncertainties), assuming the sizes of these uncertainties are known. 15 The modifications to the old framework to include correlations are detailed in Kowalski et al. (2008); Amanullah et al. (2010); Guy et al. (2010); Conley et al. (2011). We implement these correlations with nuisance parameters, as discussed in Section 2.8. The model must also allow for population drift (whether due to selection effects, or changes in the SN population with redshift), or risk a significant bias on cosmological parameters. Population drift can also be handled well with either framework, as Equation 2 handles each object independently (if the measurements are uncorrelated object to object) and our new framework solves for population drift simultaneously with other parameters (Section 2.5).
There are two disadvantages to our approach. First, the more sophisticated model requires more CPU time (now measured in many hours, rather than minutes). Second, it is not necessarily possible to generate a unique 13 In that kind of Bayesian hierarchal model, the width of the prior around zero is incorporated as a hyperparameter in the fit. The prior will smoothly weaken as the statistical evidence builds that some of these nuisance parameters are non-zero. If there is no evidence that any of these parameters are non-zero, the prior will remain clustered around zero.
14 Note that the likelihoods for the latent variables can easily become multi-modal for nonlinear standardizations, as a curved line may approach a datapoint more than once. For this reason, it is more practical to approach the problem with a technique like MCMC than a derivative-based minimizer. The lead author has had success in a frequentist framework only by testing each local χ 2 minimum for each latent variable and recording the best one. 15 There is a benefit to a Bayesian approach if the systematic uncertainties have unknown sizes, and these sizes must be marginalized over.
distance modulus for each SN as in Equation 1. The true-independent-variable priors (Section 2.5) must be treated in bins smaller than the other parameters of interest (e.g., bins of 0.1 in redshift for cosmological parameters). However, no one set of assumptions will exactly cover all use cases (e.g., tests for isotropy). Approximate sets of SN distances are possible, but we leave this for future work.
2. PROPOSED METHOD Our proposed framework surpasses previous analysis efforts by bringing together components that simultaneously address each of the limitations discussed above. We call our framework UNITY (Unified Nonlinear Inference for Type-Ia cosmologY). The required model needs to contain thousands of nonlinear fit parameters (as motivated in the following sections), which poses a problem for many statistical techniques. For example, there is no practical way to analytically marginalize over every parameter except the parameter(s) of interest. We thus must draw random samples from the posterior distribution (Monte Carlo sampling), and use those samples to estimate the posteriors. A natural tool for this sampling is Hamiltonian Monte Carlo (HMC). HMC samples efficiently, with short correlation lengths sample to sample, even for large numbers of fit parameters. To this end, we use Stan (Hoffman & Gelman 2011, 2014Betancourt 2013) through the PyStan interface (Stan Development Team 2014). Stan automatically chooses a mass matrix, and speeds up sampling efficiency even further by using a variant of HMC sampling called the "no-U-turn" sampler. Stan also incorporates automatic, analytic differentiation for computing the gradient of the log-likelihood, making the implementation of the model simple and readable. We show a Probabilistic Graphical Model of our framework in Fig Table 1. Gull (1989) offers an excellent discussion of linear regression with error bars in both dependent and independent variables. We briefly summarize here. Consider the case of fitting a straight line (slope A, intercept B) in two dimensions (y, the dependent variable vs x, the independent variable) with uncertainties in both x and y (σ x and σ y , assumed uncorrelated and Gaussian for simplicity). The generative model for an observed y ("y obs ") is

Parameterization of the True Position on the Standardization Relation
and similarly for x obs where x true and y true are the value for the measurement if no uncertainty is present, and y true = A x true + B. When there is no significant uncertainty in x, Equation 4 is unnecessary, and we have a simple least-squares regression. When uncertainty is present in both variables, we can substitute for y true , but x true remains. For a fit with two independent variables, there are two latent variables representing the true position on the line (now a plane). The same logic holds for more than one observation, in which case there are (Number of Observations)×(Number of Independent Variables) latent variables. For the standard x 1 and c standardizations, this results in 2N SNe additional parameters. 16 In a frequentist analysis, if measurement uncertainties are believed to be Gaussian and the x 1 and c standardizations are linear, the likelihood can be analytically maximized for each of these parameters (Kowalski et al. 2008). This technique enables the proper handling of these nuisance parameters without explicitly including them in the fit. Similarly, in a Bayesian analysis, if the measurement uncertainties are Gaussian (or a Gaussian mixture), the standardizations are linear, and the priors on these parameters are flat or Gaussian (or a Gaussian mixture), these nuisance parameters can be analytically marginalized (Gull 1989), providing a similar computational efficiency boost (as was done in March et al. 2011). In this work, we violate these assumptions; therefore, we must keep these 2N SNe nuisance parameters in the fit (x true 1i and c true i ).

Host-Galaxy Environment Standardization
Relation Kelly et al. (2010) first found evidence that Hubble residuals in current light-curve fitters are correlated with host-galaxy environment. This finding later reached high statistical significance in larger samples Lampeitl et al. 2010;Childress et al. 2013). The current method for including this effect in cosmological analyses fits two separate absolute magnitudes for low-mass-hosted (M * < 10 10 M ⊙ ) and high-masshosted galaxies Suzuki et al. 2012). Rigault et al. (2013) found evidence that most of this effect is due to the age of the progenitor system (see also indirect evidence in Childress et al. 2013). This effect, confirmed independently by Rigault et al. (2015), implies that SNe hosted by high-stellar-mass galaxies may become more like low-stellar-mass-hosted SNe at higher redshift (when progenitors were young in all galaxies). However, newer versions of SALT2 combined with different sample selections may not show as strong an age effect (Rigault et al. 2015;Jones et al. 2015).
We therefore do not assume a constant massstandardization (δ); instead we use a modified version of the model in Rigault et al. (2013) (similar to the model of Childress et al. 2014) Those authors' proposed host-mass-standardization evolution predicts the mass-standardization coefficient will approach zero at high-redshift; we instead assume it smoothly approaches a possibly non-zero quantity, δ(∞). We take a flat prior on δ(∞)/δ(0) from 0 to 1, allowing An edge from one node (e.g., Ωm) to another (e.g., m true Bi ) means that the latter is conditional upon the former (e.g., m true Bi is conditional on Ωm). The enclosed nodes represent variables that are sampled in the MCMC. Global parameters are in orange nodes (single parameters) and red nodes (the set of systematic uncertainty parameters). Green nodes enclose the hyperparameters (parameters of a prior distribution) of the latent-variable priors, and the singly outlined black nodes show those latent variables. Blue nodes show sample-dependent quantities. Finally, the outlined nodes show the observed light-curve fits.
} as the light-curve fit and unexplained dispersion have correlated uncertainties.) i ranges over each SN, j ranges over each SN sample, k ranges over the coefficients in redshift within a sample, and l ranges over each systematic uncertainty (e.g., calibration). Note that the m true Bi are completely determined by other parameters and are not true fit parameters. We fix the selection effect parameters, m cut j , σ cut j , a cut j , b cut j , and the outlier distribution width σ outl , so these are represented filled nodes. Note.
-The final column displays section references for the parameters in our model.
the mass standardization to be constant or declining with redshift, and spanning all of the claims in the literature.

Non-Ia or Other Outlier Contamination
For non-Ia contamination, we use the mixture-model framework of Kunz et al. (2007). This framework models the observed distribution around the modeled mean as a sum of Gaussians, where at least one Gaussian (the normal Ia distribution) is tightly clustered, and an outlying distribution is comparatively dispersed. Although the assumption of a Gaussian contaminating distribu-tion is a strong one, it makes little difference in practice. As the outlying distribution is broader than the inlying distribution, any outlying point will be treated as an outlier. For the relatively pure spectroscopically confirmed SN datasets in use today, modeling the outlying distribution accurately has little impact on the rest of the parameters in the model. Because of this, and as a test of the framework, we perform our fits assuming an outlier distribution that is a unit multinormal (σ outl = 1) in {m B , x 1 , c} (centered on the Ia distribution for that redshift), which is different from the simulated data we generate to test the framework (Section 3).
The relative normalizations of the core and outlying distributions can be chosen object by object (from spectroscopic or other classification evidence), or set to be the same for every object. We fit for the fraction of outliers assuming it is the same for all objects (f outl ), and place a broad log-normal prior on this quantity of −3 ± 0.5 (an outlier fraction of exp(−3) = 0.05 plus or minus 50%). These assumptions work well with Union2.1, as discussed in Section 4.7, but of course they can be adjusted for other datasets.
If the luminosity distribution of SNe Ia turns out to be significantly non-Gaussian (for example, the bimodal model of Rigault et al. 2013), additional Gaussian components can be added (with redshift-dependent normalizations) to give smaller uncertainties and capture possible population drift. We leave this and more complex models of the non-Ia distribution for future work, but these extensions fit easily into this framework.

Selection Effects
We present the details of our selection model in Appendix B, but outline the important points here. The standard method for incorporating a selection cut is to truncate the data likelihood at the cut and divide by the selection efficiency (e.g., Gelman et al. 2013, see also Kelly 2007 for a discussion of selection effects and nondetections in the context of linear regression). In SN cosmology, the truncation is not sharp, but is instead probabilistic, as discussed in Section 1.2. We assume that the observation likelihood is truncated by an error function. Far from the selection limit, a SN is found or missed with probability one or zero; for SNe near the selection limit, the probability transitions smoothly. An error function reasonably matches the efficiency curves of e.g., Surveys also do not select only on one measured variable. We assume our cut is a plane in three-dimensional space, spanning the dependent variable and both independent variables. In our example, this is magnitude (SALT2 m B ), shape (SALT2 x 1 ), and color (SALT2 c); SNe with m B + a cut x 1 + b cut c > m cut are less than 50% likely to be found (and more than 50% likely to be found if < m cut ). The width of the cut is σ cut ; SNe observed at m cut ± σ cut are considered 16%(+) or 84%(−) likely to be found. m B and c are the primary variables responsible for selection effects in SN searches. For example, SNe found in the rest-frame B-band have a limit in m B (b cut = 0), while SNe selected in the rest-frame V -band have a limit in m B − c ≈ m V or b cut = −1. We note that for selection in just m B , bluer (more-negative c) and slower-declining SNe (larger x 1 ) will be selected, as these correlate with brighter m B . That is, the only effect that a magnitude-based {m B , c} selection ignores is that slower-declining SNe are more likely to be found, irrespective of the maximum brightness, as they stay above the detection threshold longer. A simple simulation shows that this effect is very small compared to 17 An error-function truncation rapidly approaches 100% efficiency on the bright side of the cut. Real surveys are not as ideal, and will asymptote short of 100%. As this asymptote has little brightness dependence, it does not impact the selected population. selection on {m B , c}, even for cadences as large as ten rest-frame days. Another bias related to selection effects is the bias due to larger uncertainties on fainter SNe (e.g., Kowalski et al. 2008). Simple simulations show that our Bayesian framework has much less susceptibility to this bias, and the uncertainty bias is much smaller than the one due to missing faint SNe altogether.

Priors on True Values of Independent Variables
As this is a Bayesian framework, we must select priors on the true x 1 and c latent parameters (see Gull 1989 for a discussion of Gaussian priors, and Kelly 2007 for a Gaussian-mixture prior). These priors must be chosen very carefully. If the prior mean is wrong, then every distance will also be incorrect in a correlated direction. The variance of the prior has an impact as well. If the prior variance is larger than the population variance, then the true latent parameters will be scattered about the mean more than they should, and the slope of the line will be biased towards zero. The converse will bias the slope of the lines away from zero. The mean and variance of the prior are the most important parameters to estimate accurately, thus Gaussians are normally adequate. In SN cosmology, these priors must also be redshift-dependent as the SN population can drift with redshift.
The optimal way to ensure the proper size and redshiftdependence of the priors is to fit for the prior parameters (the "hyperparameters") simultaneously with every other parameter. We selected skew normal distributions for the c priors (allowing the distribution to be skewed), and Gaussians for the x 1 priors. The prior must be able to vary in redshift more rapidly than the cosmological fit in order to not introduce a bias. What we propose here more than meets this mild requirement, 18 but there is no harm in allowing the hyperparameters to mimic more closely the redshift dependence of the x 1 and c distributions. For each sample, we fit for the mean of the distributions as a function of redshift with a linear spline. We use up to four spline nodes (the x * 1jk and c * jk , for x 1 and c, respectively), equally spaced in redshift over the range of a sample, with linear interpolation between these nodes. 19 We take non-informative flat priors on the means of the distributions and on the log of the standard deviations (the standard deviations for each sample j are R x1 j and R c j ). 20 For the shape parameter α S−N of the skew normals, we take a flat prior on δ S−N = α S−N / 1 + (α S−N ) 2 , which is also allowed to vary in redshift in the same way as the distribution mean. This prior forces the skew normal to approach a Gaussian for samples with few objects. (The superscript "S−N" here is used to distinguish these skew-normal variables from the SN standardization coefficients.) For simplicity, we assume that the true x 1 and c distributions are uncorrelated (the observed distributions of x 1 18 In contrast, the constant-in-redshift priors of March et al. (2011) do not meet this requirement.
19 Specifically, we use two nodes (a single line segment) for datasets with fewer than 30 SNe; we use three for 30-39 and four when the number of SNe is at least 40. This improves the robustness of the determination of the independent variable standard deviations.
20 For a skew normal distribution, the mean and standard deviation are not the same quantities as µ and σ. and c show little correlation). We do note that ignoring correlations will bias the fit if any significant (|ρ| 0.2) correlations are present.

Nonlinear Independent Variable standardizations
With the true values of the independent variables explicit in the model, it becomes trivial to have nonlinear standardizations. For this work, we suggest a brokenlinear relationship, allowing red/blue and small-/largex 1 SNe to have different size standardization coefficients (β R /β B and α S /α L , respectively). We take a flat prior on the angle of each line segment, but transform to the average slope and the difference in slopes for display purposes: Although the x 1 and c values of the break could be fit parameters, we do not do this for our primary fits. For the moment, we split the sample at x 1 and c of 0.

Unexplained Dispersion
As the unexplained dispersion is parameterized in the model, it can be marginalized. We do not know what functional form to assume, so we use a flexible parameterization. Each SN sample is allowed its own unexplained dispersion, allowing poorer-quality samples to be naturally deweighted. We must also distribute the unexplained dispersion over m B , x 1 , and c, while accounting for possible correlations.
First, we split the variance of the unexplained dispersion into f mB , f x1 , and f c (the fraction of the unexplained variance in m B , x 1 , and c, respectively), which are constrained to sum to one. Then, we scale each of these by 1, 0.13 −2 , and (−3.0) −2 . The values 0.13 and −3.0 approximately scale out α and −β, respectively, where the negative sign for β corresponds to the sign convention in Equation 1. (Note that using α and −β directly would cancel the α or β dependence when computing a marginalized distance uncertainty for each SN, so this would be inappropriate.) We also scale the variance by the unexplained dispersion for each sample, σ 2samp j . Finally, we form a covariance matrix out of this {m B , x 1 , c} unexplained variance, allowing the off-diagonals to be scaled by parameters ρ samp , as follows We take a non-informative "LKJ" prior on the correlation distribution (Lewandowski et al. 2009), with η = 1, as well as a flat prior on log σ samp j .

Treatment of Correlated SN Observations
Many effects result in correlated measurements in SN cosmology. The most notable such effect results from common calibration paths: SNe from a given dataset share systematics such as telescope bandpass and photometric calibration uncertainties. Other effects include correlated uncertainties in Milky Way extinction maps, and correlated peculiar velocities. In standard analyses, these effects are propagated into a covariance matrix (Amanullah et al. 2010;Guy et al. 2010;Conley et al. 2011).
In order to speed up the Monte-Carlo sampling, we leave each correlating factor explicit as a parameter (similar to Kowalski et al. 2008), ∆sys l (where l ranges over all systematic uncertainties), while leaving the data uncorrelated SN to SN. These two approaches coincide exactly for a linear model with Gaussian uncertainties (Amanullah et al. 2010). The ∆sys l parameters capture the deviations of a measured quantity, like a zeropoint or a filter bandpass, from the estimated value. For each quantity, we numerically compute the derivative of the light-curve fit m obs Bi , x obs 1i , c obs i with respect to that quantity, giving each ∂m obs Bi /∂∆sys l , ∂x obs 1i /∂∆sys l , and ∂c obs i /∂∆sys l . This lets us marginalize out ∆sys l , with a Gaussian prior around zero set by the estimated size of each systematic uncertainty.
3. SIMULATED DATA TESTING Our analysis framework must pass through careful testing using simulated data before it can be applied to real data. To this end, we generate thirty simulated datasets that incorporate many characteristics of real data. As our analysis takes the SALT2 light-curve fits as inputs, this is the level at which we generate simulated data.
• Each simulated dataset has 250 SNe, except the highest-redshift, which has 50.
• We generate the x 1 population from a unit normal distribution, centered on zero.
• We draw the population c values from the sum of a Gaussian distribution of width 0.1 magnitudes, and an exponential with rate 1/(0.1 magnitudes). We center the distribution on zero.
• We assume that the unexplained dispersion covariance matrix is correct in SALT2, and that only dispersion in m B (gray dispersion) remains. The statistical model does not have access to this information, and fits for the full unknown matrix, overestimating the uncertainties on x 1 and c, and thus slightly biasing α and β away from zero (see Section 2.5). (This is not a unique problem for our framework; the old technique would have the same bias.) • We assume that the uncertainties on m B , x 1 , and c are 0.05, 0.5, and 0.05, and are uncorrelated. In addition, we take 0.1 magnitudes of unexplained dispersion, 0.093z magnitudes of lensing dispersion (Holz & Linder 2005), and 300 km/s peculiar velocity uncertainty. All are approximated as Gaussian and independent SN to SN.
• α and β are assumed to be constant, with values 0.13 and 3.0, respectively. M B is set to −19.1 and Ω m is set to 0.3 (flat ΛCDM model). Note.
-Average over thirty simulated datasets. These results show an expected bias towards larger values of α and β (discussed in Section 3), but Ωm shows no significant bias. Other simulated data fits are also described in Section 3.
• We assume 3% of the SNe are outliers, and draw their observed distribution centered around the normal Ia m B for that redshift, and around zero in x 1 and c. The spread is 1, 2, and 0.5 in m B , x 1 , and c respectively, and we assume these distributions are Gaussian and uncorrelated.
• Each sample has zeropoint uncertainties of size 0.01, 0.01, 0.01, and 0.02 (highest-redshift sample) magnitudes. The uncertainties are taken to be independent sample to sample.
• The datasets have selection effects in m B , with width 0.2 magnitudes. The selection cuts are chosen for 50% completeness at redshifts 0.08, 0.25, 0.6, and 1.45. Note that this selects from the population x 1 and c distributions in a redshiftdependent way. (We randomly draw from the population distributions and pass them through the simulated selection effects until the required number of SNe are generated.) • We assume the redshift distribution of SNe scales linearly with redshift (starting from the minimum redshift of each sample). This is quantitatively incorrect (the real scaling will depend on the cosmological volume, cosmological time dilation, and the SN rate), but does produce SN samples where many SNe are removed by the magnitude cuts, allowing a good test of our selection effect modeling.
The redshift distribution of the four simulated datasets are shown in Figure 2 before selection effects (gray-tinted histograms) and as observed. We show the observed color distributions in Figure 3. This figure shows outliers in gray (knowledge that the statistical framework does not have access to), and the general trend towards bluer SNe at higher redshift due to selection effects. Figure 4 shows the Hubble diagram residuals from the input cosmology with the best-fit ignoring selection effects (Ω m = 0.32) also shown.
After generating thirty complete sets (with all four SN samples), we supply the generated SALT2 result files to the framework. Table 2 summarizes the key results. As expected, α, β show bias away from zero, but any bias on Ω m is small. We also run some simulated datasets fitting for both Ω m and the Dark Energy equation of state parameter w (assumed to be constant in redshift). 21 We see no evidence of bias on w, and mild evidence of bias on the mean Ω m : 0.274 ± 0.015. 2014, did not publish these SNe). Our cosmological fits include Ω m and assume a flat universe. This fit is qualitatively similar to the assumption of a constant equation-of-state parameter w including a CMB or BAO constraint, in that both fits probe the deceleration parameter q 0 . As fitting Ω m only requires SN data, it is a cleaner analysis for our purposes here.
In order to identify the effect of each feature of the analysis on the inferred results, we incrementally transition from the original Union2.1 frequentist framework to the analysis proposed in this work. The results of each step are shown in Figure 5. We conducted this part of the analysis blinded, using real data only after the code was validated on simulated data. The initial version of UNITY required the unexplained dispersion in m B to be fixed, which the improved selection effect model now presented in Appendix B does not require. With the improvements in place, after a second round of blindingunblinding, we found only a small change in Ω m between the two versions, 0.009 (0.2σ).

Frequentist Union2.1 Analysis
First, we show the results of a frequentist calculation, based on the same assumptions as Suzuki et al. (2012) and using its 580 SNe (top line in Figure 5) with SALT2 light-curve fits. 22 All systematic uncertainties from the covariance matrix are included. To reproduce the Suzuki et al. (2012) results, we include only a redshift-independent host-mass standardization. 23 As a cross-check, we also try a hybrid frequentist/Bayesian model, in which the χ 2 from Equation 2 is converted into a likelihood as e −χ 2 /2 and then M B , α, β, and δ are marginalized over (this is similar to the method used in Knop et al. 2003). We obtain essentially the same results as a purely frequentist fit.

Bayesian Model with Same Data
For the next step, we keep all data the same, but transition to a Bayesian model for the data (the credible intervals are shown as the second line in Figure 5). This model includes the SN population terms (described in Section 2.5 and necessary for a Bayesian analysis) and the Union2.1 systematic uncertainties, but does not include our proposed treatment of selection effects, outliers, multi-dimensional unexplained dispersion, or the redshift-dependent host-mass standardization. Other than the type of inference and these differences, this fit is identical to the first fit. The error bars shrink by 5%, but the central value changes very little. As the data are the same for this fit and the last one, the gain in statistical power comes from the ability of a Bayesian hierarchical model to make better use of heterogeneous information. 24 24 A Bayesian hierarchical model and a frequentist line fit will give exactly the same slope and intercept if: all uncertainties are known and Gaussian, every measurement has the same uncertainties (homoscedasticity), the Bayesian hyperpriors are Gaussian, flat priors are taken on all other parameters, and the slope and intercept are evaluated in the Bayesian model at the hyperparameter posterior maximum. These methods will give different results if the measurement uncertainties are heteroscedastic. For example, if there is a dataset with shape and color measurements but no SN has both measurements, Equation 2 will have no constraining power, but a Bayesian model will. The difference between the Bayesian and frequentist Union2.1 fits is thus a measure of how inhomogeneous the uncertainties are within a given redshift range for each SN sample.

Frequentist Redshift-Dependent Host-Mass
Standardization Next, we return to the frequentist framework, but include the redshift-dependent mass standardization as described in Section 2.2 (third line in Figure 5). We remove the Suzuki et al. (2012) covariance matrix term that corresponds to systematic uncertainty on the mass standardization. The best fit shifts to a higher Ω m (brighter standardized magnitudes on average at high redshift) because the high-host-mass half of the high-redshift SNe has less standardization to fainter magnitudes. This step can be conducted through frequentist or Bayesian inference, and we show both analyses.

Bayesian Model with Redshift-Varying Host-Mass
Standardization In this next fit, we continue with the same data and model as the last subsection (redshift-dependent hostmass standardization), but again transition to Bayesian inference, including the population terms (fourth line of Figure 5). The results are similar to the frequentist results, but the Bayesian fit is more agnostic about the value of δ(∞) (essentially unconstrained), so the fit shifts less than the frequentist one to higher Ω m .

Bayesian Model with Unexplained Dispersion
Remaining in the Bayesian framework, our next addition is the multi-dimensional unexplained dispersion. We first remove the existing Union2.1 unexplained dispersion (which is only in m B ). By effectively increasing the uncertainties on the color, we increase the colorstandardization coefficient β. Thus the color standardization now moves bluer (bluer due to selection effects) high-redshift SNe fainter, decreasing the fitted Ω m (the fourth line from the bottom in Figure 5).

Bayesian Model with Nonlinear Standardizations
We now take advantage of our explicit x true 1 and c true values to include nonlinear color standardizations parameterized by dα and dβ. We remove the Union2.1 covariance terms that describe color-standardization systematics (between the multivariate unexplained dispersion and the nonlinear standardizations, these systematic uncertainties are likely much lower). The color standardization is now strongly nonlinear (discussed further in Section 4.9), with redder SNe requiring a larger coefficient than bluer SNe. This moves the fitted Ω m (third line from the bottom in Figure 5) back in the opposite direction from the previous step.

Bayesian Model with Outliers
Next, we include in the fit all twelve outlier SNe removed by the Union sigma clipping in Union2.1 (a new total of 592). Instead of excluding these, we add our mixture model for handling outlier SNe. We also remove the Union2.1 systematic uncertainties on outlier rejection. The results are quite similar to the previous step, indicating that the Union sigma-clipping worked well with the 2% contamination that was present (second line from the bottom in Figure 5).

Bayesian Model with Selection Effects
Finally, we model selection effects. For simplicity, we approximate all selection as occurring in rest-frame Bband (a cut j = 0, b cut j = 0), and leave a more detailed analysis to future work. 25 For many of the most important samples in Union2.1, this is not a bad approximation. For example, the Sloan Digital Sky Survey SNe (Holtzman et al. 2008) were selected in g, r, and i-band (Dilday et al. 2008). At redshift ∼ 0.4 (the distant end of the survey), r-band corresponds to rest-frame B-band. The Supernova Legacy Survey SNe were selected in iband (Perrett et al. 2012), which matches rest-frame Bband for the highest-redshift SNe with small distance uncertainties (z ∼ 0.7).
For some surveys, rest-frame B-band selection is a poor approximation. Many of the nearby SNe were selected from unfiltered surveys (approximately rest-frame R-band), but these mostly galaxy-targeted surveys have generally weak selection effects in magnitude (discussed below). Some distant SN surveys (Tonry et al. 2003;Riess et al. 2007;Amanullah et al. 2010;Suzuki et al. 2012;Rubin et al. 2013) had selections that were different (bluer) than rest-frame B-band. At least for the HST-discovered SNe, selection effects are small for most of the redshift range (also discussed more below).
We estimate the selection effects for each sample as follows (summarized in Table 3). Nearby SNe are generally limited by spectroscopic followup, for example the ∼ 18.5 magnitude limit for CfA discussed in Hicken et al. (2009) (although this is unfiltered, we approximate it as R-band). We take this limit as typical. The Calàn/Tololo survey (Hamuy et al. 1996a,b) and the SCP nearby search (Kowalski et al. 2008) extend out to higher redshift; the limiting magnitude in this case is R ∼ 19 (Hamuy & Pinto 1999). We take the magnitude limit for SDSS from Dilday et al. (2008) and the limit for SNLS from Perrett et al. (2012). Together, these samples make up most of the mid-redshift weight. For the other mid-redshift samples, we take a limit of R = 24, judged from the approximate rolloff of the SN population in redshift. For the high-redshift ground-discovered samples, we assume the surveys were 50% complete at z = 1; this gives a limit of 23.8 in I-band. Finally, for the HST-discovered SNe, we take a limit from Barbary et al. (2012) of z-band ∼ 25 (Vega). In all cases, we take the width of the selection to be ±0.5 magnitudes (i.e., a SN 0.5 magnitudes brighter than the mean cut has an 84% chance of being selected; a SN 0.5 magnitudes fainter has a 16% chance). As a cross-check, we coherently shift each estimated magnitude limit fainter by 0.5 magnitudes, representing an extreme limit of how inaccurate our estimations are likely to be. The Ω m credible interval shifts by only 0.006.
We remove the Union2.1 covariance matrix terms for Malmquist bias before computing the fit shown in the bottom line of Figure 5. The Ω m credible interval shifts to lower Ω m , as the distant SNe with significant selection effects are standardized fainter. The central value of this final fit closely matches the original Union2.1 result, but we see that our new, smaller (by 7% compared to the Union2.1 analysis) credible interval shows an increase in statistical power. It is worth reiterating that the list Note.
-Approximate magnitude limits for each sample. We use a SALT2 c = 0 SN template to convert each magnitude limit to the limit in rest-frame B-band magnitude at the specified redshift.
of improvements to make was established before the results were known; we did not set out to simply achieve a similar result to Union2.1. The scatter of the intermediate results generally validates the size of the Union2.1 systematics estimates for these effects that we can now properly include in the model.

Other Parameters
We now present the results for important nuisance parameters and their relation to Ω m in the form of 1D and 2D credible regions. In both cases, our credible regions are derived by using Kernel Density Estimation with a Gaussian kernel on the MCMC samples, 26 then solving for the contour level that encloses 68.3% (inner shaded regions) and 95.4% (outer shaded regions) of the posterior. 27 Figure 6 shows the significant degeneracy between both host-mass standardization parameters and Ω m , illustrating that neither one should be neglected. The mean value of the estimated fraction of outliers is similar to the 12/592 found by the Union sigma clipping. This parameter has no significant degeneracy with any parameter for the spectroscopically confirmed datasets in Union2.1, confirming that outlier rejection is not a significant concern at this high level of purity. As an additional crosscheck, instead of assuming the outlier distribution is centered on the SN Ia distribution, we fit for it. Including six parameters in the model for the mean and dispersion in {m B , x 1 , c} (taken to be uncorrelated, and assuming a constant distribution with redshift) leaves the error bar unchanged but shifts the credible region by 0.009 in Ω m . If future versions of the UNITY framework are run on samples with larger contamination, a more flexible outlier parameterization can be matched with the increased number of outliers in the fit (e.g., Hlozek et al. 2012). Figure 7 shows the degeneracy between the fraction of the unexplained variance in x 1 and α, and similarly c/β. Including these unexplained dispersion parameters increases the uncertainty and the mean value for α and β. (Here, α and β represent the mean standardization 26 We draw ∼3,000 samples from sixteen chains. The Gelman & Rubin (1992) R statistics are 1.01, indicating good convergence.
27 Kernel Density Estimation gives biased results when the posterior samples lie against a parameter boundary, as the density will be smoothed including regions with zero density. For all our variables, we reflect the nearby samples through the minimum and maximum samples, creating a virtual dataset on the other side of such boundaries.
coefficient.) The model prefers most of the unexplained dispersion in m B and c, rather than x 1 . We also see statistical evidence of non-zero dα and dβ (recall that dα is the x 1 standardization coefficient for broad-light-curve SNe minus the coefficient for narrow SNe; dβ is similarly defined for red minus blue SNe): dα is −0.089 +0.035 −0.040 and dβ is 2.26 +0.43 −0.38 . 28 Both have a correlation with Ω m , although the correlation with dβ is larger. We note, at least for this compilation of SNe and SALT2-1, that blue SNe still have a non-zero β: 0.56 +0.31 −0.34 , and red SNe have a β of 2.83 +0.18 −0.17 . This value is significantly less than 4.1 (the value expected if all reddening were due to the mean extinction law of the Milky Way diffuse interstellar medium). 29 5. CONCLUSIONS In this work, we propose UNITY, a unified Bayesian model for handling outliers, selection effects, shape and color standardizations, unexplained dispersion, and heterogeneous observations. We demonstrate the method with the Union2.1 SN compilation, and show that our method has 7% smaller uncertainties, but results that are consistent with the Union2.1 analysis. The advantages of UNITY will likely be even larger in upcoming datasets, as scarce followup resources will introduce heterogeneity, and enlarged samples may reduce the cosmologicalparameter statistical uncertainties below the size of our improvements.
There are several future directions for this research. We could allow the unexplained dispersion covariance matrix to vary with other parameters (x 1 , c, redshift, or rest-frame wavelength range). We could decompose the unexplained dispersion into a sum of Gaussians as in the model of Rigault et al. (2013). We could use Gaussian process regression (Rasmussen & Williams 2006) to handle the redshift-dependent priors on x 1 and c. Our model also lets us include more than one light-curve fit for each SN, with the covariance matrix between different light-curve fitters parameterized. We could even make use of the hierarchical model to constrain "exotic possibilities," as sketched in Section 1.3. All of these are straightforward modifications of what we present here, but our proposed model is already superior to current SN cosmological analysis frameworks. We believe these concepts will also be useful for other applications.
P (detect|{m obs Bi , x obs 1i , c obs i }) and P (detect) are described in Equations B2 and B10, respectively. m true Bi is not a free parameter; instead it is completely determined by other parameters and is given by ∆{m B , x 1 , c} are the contributions from the systematic uncertainty terms, e.g., The broken-linear slopes are defined as and δ(z i ) is given by Equation 5. We take the following priors on all parameters: DETAILS OF THE SELECTION-EFFECT MODEL In this section, we detail the approach we take for modeling selection effects in UNITY. We start with a general likelihood containing missing observations in the case when both the number of observed objects (N obs ) and the number of missed objects (N miss ) are exactly known: Bi , x obs 1i , c obs i }|params)P (detect|{m obs Bi , x obs 1i , c obs i }) (1 − P (detect))   . (B1) As described in Section 2.4, we assume that the efficiency smoothly varies with {m obs Bi , x obs 1i , c obs i }: where Φ is the Gaussian CDF. Of course, we do not know any of the parameters for each of the missing observations. P (detect) is thus marginalized over the entire distribution when it is not referencing a specific observation. We now make the counterintuitive approximation that the redshift of each missed SN is exactly equal to the redshift of a detected SN. This approximation is accurate because the SN samples have, on average, enough SNe that the redshift distribution is resonably sampled. We now refactor Equation B1 to consider each detected SN: The combinatoric factor trivially becomes (1 + N miss i ). To minimize the number of parameters in the cosmological fit, we now seek to marginalize over (1 + N miss i ) from 1 to ∞. It is common (e.g., Gelman et al. 2013;Kelly 2007) to take a flat prior on log (N miss i ); for reasons that will become clear in a moment, we force this prior to zero faster than this by multiplying this (1 + N miss i ) −1 prior by a weak exponential decay P (1 + N miss i ) = exp (−ǫ(N miss i + 1)) 1 + N miss where ǫ is a small positive number. Thus, the marginalization becomes: Bi , x obs 1i , c obs i }|params)P (detect|{m obs Bi , x obs 1i , c obs i }) ǫ + P (detect|z i ) This expression highlights the benefit of this prior. If a flat prior on log (N miss i ) is assumed (ǫ = 0), then the likelihood can be poorly behaved in regions of parameter space where the efficiency is poor (even if care is taken in the log likelihood to correctly handle the asymptotes of the log of the efficiency). Running UNITY on Union2.1 with ǫ either 0 or 0.01 makes virtually no difference, but Stan has difficulty with many of the simulated datasets (which we constructed with severe Malmquist bias as a test, see Section 3) for ǫ = 0, so we make ǫ = 0.01 our default choice.
We now need the efficiency as a function of the parameters. We derived the following expression by analytically computing the variance, then expanding to first order in x * 1 (this approximation is good to better than one percent in the variance for our distributions). For simplicity, we use the same expression for c, although the c distribution is skewnormal instead of normal. (Changing the assumed intrinsic color distribution to normal, rather than skew-normal, changes Ω m by only 0.001, so this approximation is not a significant concern.) V mB = C mB i + σ 2 cut + (δ(z i )/2) 2 (B8) + ((α − a cut j ) R x1 ) 2 + (α − a cut j ) dα x * 1 R x1 2 π + (dα R x1 ) 2 π − 2 4π The first line contains the measurement and unexplained dispersion in m B , the dispersion of the selection cut, and the dispersion due to the host-galaxy relation. The second and third lines give the dispersion due to the x 1 − m B and c − m B relations, respectively. Using the same approximations, we derived an expression for the mean of the (unstandardized) magnitude distribution in this case, our approximations are valid to a few hundredths of a magnitude. We then compute the selection efficiency assuming that the distribution of magnitudes is Gaussian This approximation is not completely valid, as there is a tail to faint magnitudes of red SNe, but for the central part of the distribution (which is the part that is cut into by severe Malmquist bias), it is a good approximation.