Improved constraints on cosmological parameters from SNIa data

We present a new method based on a Bayesian hierarchical model to extract constraints on cosmological parameters from SNIa data obtained with the SALT-II lightcurve fitter. We demonstrate with simulated data sets that our method delivers tighter statistical constraints on the cosmological parameters over 90% of the time, that it reduces statistical bias typically by a factor ~ 2-3 and that it has better coverage properties than the usual chi-squared approach. As a further benefit, a full posterior probability distribution for the dispersion of the intrinsic magnitude of SNe is obtained. We apply this method to recent SNIa data, and by combining them with CMB and BAO data we obtain Omega_m=0.28 +/- 0.02, Omega_Lambda=0.73 +/- 0.01 (assuming w=-1) and Omega_m=0.28 +/- 0.01, w=-0.90 +/- 0.05 (assuming flatness; statistical uncertainties only). We constrain the intrinsic dispersion of the B-band magnitude of the SNIa population, obtaining sigma_mu^int = 0.13 +/- 0.01 [mag]. Applications to systematic uncertainties will be discussed in a forthcoming paper.


INTRODUCTION
Since the late 1990s when the Supernova Cosmology Project and the High-Z Supernova Search Team presented evidence that the expansion of the universe is accelerating (Riess et al. 1998;Perlmutter et al. 1999), observations of Type Ia supernovae (SNIa) has been seen as one of our most important tools for measuring the cosmic expansion as a function of time. Since a precise measurement of the evolution of the scale factor is likely the key to characterizing the dark energy or establishing that General Relativity must be modified on cosmological scales, the limited data that our universe affords us must be used to the greatest possible advantage. One important element of that task is the most careful possible statistical analysis of the data. Here we report how an improved statistical approach to SNIa data, with, in particular, more consistent treatment of uncertainties, leads to significant improvements in both the precision and accuracy of inferred cosmological parameters.
The fundamental assumption underlying the past and proposed use of Type Ia supernovae to measure the expansion history is that they are "standardizable candles". Type Ia supernovae occur when material accreting onto a white dwarf from a companion drives the mass of the white dwarf above the maximum that can be supported by electron degeneracy pressure, the Chandrasekhar limit of about 1.4 solar masses. This triggers the collapse of the star and the explosive onset of carbon fusion, in turn powering the supernova explosion. Because the collapse happens at a particular critical mass, all Type Ia supernovae are similar. Nevertheless, variability in several factors, including composition, rotation rate, and accretion rate, can lead to measurable differences in supernova observables as a function of time. Indeed, the intrinsic magnitude of the nearby Type Ia supernovae, the distances to which are known via independent means, exhibit a fairly large scatter. Fortunately, this scatter can be reduced by applying the so-called "Phillips corrections" -phenomenological correlations between the intrinsic magnitude of SNIa and their colour as well as between their intrinsic magnitudes and the time scale for the decline of their luminosity (Phillips 1993;Phillips et al. 1999). Such corrections are derived from multi-wavelength observations of the SNIa lightcurves (i.e., their apparent brightness as a function of time). Fortuitously, they make SNIa into standardizable candles -in other words, having measured the colour and light curve of a SNIa, one can infer its intrinsic magnitude with relatively low scatter, typically in the range 0.1 − 0.2 mag. Observations of SNIa at a range of redshifts can then be used to measure the evolution of luminosity distance as a function of redshift, and thence infer the evolution of the scale factor, assuming that the intrinsic properties of SNIa do not themselves evolve (an assumption that has to be carefully checked).
The SNIa sample which is used to measure distances in the Universe has grown massively thanks to a world-wide observational effort (Astier et al. 2006;Wood-Vasey et al. 2007;Amanullah et al. 2010;Kowalski et al. 2008;Kessler et al. 2009a;Freedman et al. 2009;Contreras et al. 2010;Balland et al. 2009;Bailey et al. 2008;Hicken et al. 2009). Presently, several hundred SNIa have been observed, a sample which is set to increase by an order of magnitude in the next 5 years or so. As observations have become more accurate and refined, discrepancies in their modeling have come into focus. Two main methods have emerged to perform the lightcurve fit and derive cosmological parameter constraints. The Multi-Colour Lightcurve Shape (MLCS) ) strategy is to simultaneously infer the Phillips corrections and the cosmological parameters of interest, applying a Bayesian prior to the parameter controlling extinction. The SALT and SALT-II (Guy et al., 2007) methodology splits the process in two steps. First, Phillips corrections are derived from the lightcurve data; the cosmological parameters are then constrained in a separate inference step. As the supernova sample has grown and improved, the tension between the results of the two methods has increased.
Despite the past and anticipated improvements in the supernova sample, the crucial inference step of deriving cosmological constraints from the SALT-II lightcurve fits has remained largely unchanged. For details of how the cosmology fitting is currently done, see for example (Astier et al. 2006;Kowalski et al. 2008;Amanullah et al. 2010;. As currently used, it suffers from several shortcomings, such as not allowing for rigorous model checking, and not providing a rigorous framework for the evaluation of systematic uncertainties. The purpose of this paper is to introduce a statistically principled, rigorous, Bayesian hierarchical model for cosmological parameter fitting to SNIa data from SALT-II lightcurve fits. In particular the method addresses identified shortcomings of standard chi-squared approaches -notably, it properly accounts for the dependence of the errors in the distance moduli of the supernovae on fitted parameters. It also treats more carefully the uncertainty on the Phillips colour correction parameter, escaping the pontetial bias caused by the fact that the error is comparable to the width of the distribution of its value. We will show that our new method delivers considerably tighter statistical constraints on the parameters of interest, while giving a framework for the full propagation of systematic uncertainties to the final inferences. (This will be explored in an upcoming work.) We also apply our Bayesian hierarchical model to current SNIa data, and derive new cosmological constraints from them. We derive the intrinsic scatter in the SNIa absolute magnitude and obtain a statistically sound uncertainty on its value. This paper is organized as follows: in section 2 we review the standard method used to perform cosmological fits from SALT-II lightcurve results and we describe its limitations. We then present a new, fully Bayesian method, which we test on simulated data in section 3, where detailed comparisons of the performance of our new method with the standard approach are presented. We apply our new method to current SN data in section 4 and give our conclusions in section 5.

Definition of the inference problem
Several methods are available to fit SNe lightcurves, including the MLCS method, the ∆m15 method, CMAGIC, (Wang et al., 2003;Conley et al., 2006) SALT, SALT-II and others. Recently, a sophisticated Bayesian hierarchical method to fit optical and infrared lightcurve data has been proposed by Mandel et al. (2009Mandel et al. ( , 2010. As mentioned above, MLCS fits the cosmological parameters at the same time as the parameters controlling the lightcurve fits. The SALT and SALT-II methods, on the contrary, first fit to the SNe lightcurves three parameters controlling the SN magnitude, the stretch and colour corrections. From those fits, the cosmological parameters are then fitted in a second, separate step. In this paper, we will consider the SALT-II method (although our discussion is equally applicable to SALT), and focus on the second step in the procedure, namely the extraction of cosmological parameters from the fitted lightcurves. We briefly summarize below the lightcurve fitting step, on which our method builds. The rest-frame flux at wavelength λ and time t is fitted with the expression dFrest dλ where M0, M1, CL are functions determined from a training process, while the fitted parameters are x0 (which controls the overall flux normalization), x1 (the stretch parameter) and c (the colour correction parameter). The B-band apparent magnitude m * B is related to x0 by the expression where T B (λ) is the transmission curve for the observer's B-band, and t = 0 is by convention the time of peak luminosity. After fitting the SNIa lightcurve data with SALT-II algorithm, e.g. Kessler et al. (2009a) report the best-fit values for m * B , x1, c, the best-fit redshift z of each SNIa and a covariance matrixĈi for each SN, describing the covariance between m * B , x1, c from the fit, of the formĈ Let us denote the result of the SALT-II lightcurve fitting procedure for each SN as (where i runs through the n SNe in the sample, and measured quantities are denoted by a hat). We assume (as it is implicitly done in the literature) that the distribution ofm * Bi ,x1i,ĉi is a multi-normal Gaussian with covariance matrixĈi. The distance modulus µi for each SN (i.e., the difference between its apparent B-band magnitude and its absolute magnitude) is modeled as: where Mi is the (unknown) B-band absolute magnitude of the SN, while α, β are nuisance parameters controlling the stretch and colour correction (so-called "Phillips corrections"), respectively, which have to be determined from the data at the same time as the parameters of interest. The purpose of applying the Phillips corrections is to reduce the scatter in the distance modulus of the supernovae, so they can be used as almost standard candles. However, even after applying the corrections, some intrinsic dispersion in magnitude is expected to remain. Such intrinsic dispersion can have physical origin (e.g., host galaxies properties such as mass (Kelly et al. 2010;) and star formation rate (Sullivan et al. 2006), host galaxy reddening , possible SNe Ia evolution (Gonzalez-Gaitan et al. 2011), etc) or be associated with undetected or underestimated systematic errors in the surveys. Below, we show how to include the intrinsic dispersion explicitly in the statistical model. Turning now to the theoretical predictions, the cosmological parameters we are interested in constraining are where Ωm is the matter density (in units of the critical energy density), ΩΛ is the dark energy density, w is the dark energy equation of state (taken to be redshift-independent, although this assumption can easily be generalized) and h is defined as H0 = 100hkm/s/Mpc, where H0 is the value of the Hubble constant today 1 . The curvature parameter Ωκ is related to the matter and dark energy densities by the constraint equation In the following, we shall consider either a Universe with non-zero curvature (with Ωκ = 0 and an appropriate prior) but where the dark energy is assumed to be a cosmological constant, i.e. with w = −1 (the ΛCDM model), or a flat Universe where the effective equation of state parameter is allowed to depart from the cosmological constant value, i.e. Ωκ = 0, w = −1 (the wCDM model).
In a Friedman-Robertson-Walker cosmology defined by the parameters C , the distance modulus to a SN at redshift zi is given by where DL denotes the luminosity distance to the SN. This can be rewritten as µi = η + 5 log dL(zi, Ωm, ΩΛ, w), where η = −5 log 100h c + 25 and c is the speed of light in km/s. We have defined the dimensionless luminosity distance (with DL = c/H0dL, where c is the speed of light) with the dark energy density parameter In the above equation, we have been completely general about the functional form of the dark energy equation of state, w(z).
In the rest of this work, however, we will make the further assumption that w is constant with redshift, i.e., w(z) = w. We have defined the function sinn(x) = x, sin(x), sinh(x) for a flat Universe (Ωκ = 0), a closed Universe (Ωκ < 0) or an open Universe, respectively. The problem is now to infer, given data D in Eq. (4), the values (and uncertainties) of the cosmological parameters C , as well as the nuisance parameters {α, β}, appearing in Eq. (5) and any further parameter describing the SNe population and its intrinsic scatter. Before building a full Bayesian hierarchical model to solve this problem, we briefly describe the usual approach and its shortcomings.

Shortcomings of the usual χ 2 method
The usual analysis (e.g., Astier et al. (2006); Kowalski et al. (2008); Kessler et al. (2009a);) defines a χ 2 statistics as follows: where µi is given by Eq. (9) as a function of the cosmological parameters and the "observed" distance modulus µ obs i is obtained by replacing in Eq. (5) the best-fit values for the colour and stretch correction and B-band magnitude from the SALT-II output (denoted by hats). Furthermore, the intrinsic magnitude for each SN, Mi, is replaced by a global parameter M , which represents the mean intrinsic magnitude of all SNe in the sample: where the mean intrinsic magnitude M is unknown. The variance σ 2 µi comprises several sources of uncertainty, which are added in quadrature: where σ fit µi is the statistical uncertainty from the SALT-II lightcurve fit, where Ψ = (1, α, −β) andĈi is the covariance matrix given in Eq.
(3). σ z µi is the uncertainty on the SN redshift from spectroscopic measurements and peculiar velocities of and within the host galaxy; finally, σ int µ is an unknown parameter describing the SN intrinsic dispersion. Further discussions of the unknown σ int µ estimation problem, see (Blondin,Mandel& Kirshner 2011; Kim 2011; Vishwakarma& Narlikar 2011). As mentioned above, ideally σ int µ is a single quantity that encapsulates the remaining intrinsic dispersion in the SNIa sample, folding in all of the residual scatter due to physical effects not well captured by the Phillips corrections. However, observational uncertainties such as the estimation of photometric errors can lead to a variation of σ int µ sample by sample (for which there is a growing body of evidence). While we do not consider the latter scenario in this paper, it is important to keep in mind that describing the whole SN population with a single scatter parameter σ int µ is likely to be an oversimplification. Further error terms are added in quadrature to the total variance, describing uncertainties arising from dispersion due to lensing, Milky Way dust extinction, photometric zero-point calibration, etc. In this work, we do not deal with such systematic uncertainties, though they can be included in our method and we comment further on this below.
The cosmological parameter fit proceeds by minimizing the χ 2 in Eq. (13), simultaneously fitting the cosmological parameters, α, β and the mean intrinsic SN magnitude M . The value of σ int µ is adjusted to obtain χ 2 µ /dof ∼ 1 (usually on a sample-by-sample basis), often rejecting individual SNe with a residual pull larger than some arbitrarily chosen cut-off. It was recognized early that fitting the numerator and denominator of Eq. (13) iteratively leads to a large "bias" in the recovered value of β (Kowalski et al. 2008;Astier et al. 2006;Wang et al. 2006). This has been traced back to the fact that the error on the colour correction parameter ci is as large as or larger than the width of the distribution of values of ci, especially for high-redshift SNe. This is a crucial observation, which constitutes the cornerstone of our Bayesian hierarchical model, as explained below. We demonstrate that an appropriate modeling of the distribution of values of ci leads to an effective likelihood that replaces the χ 2 of Eq. (13). With this effective likelihood and appropriate Bayesian priors, all parameters can be recovered without bias. If instead one adopts a properly normalized likelihood function, i.e., replacing the χ 2 of Eq. (13) with (with the pre-factor L0 chosen so that the likelihood function integrates to unity in data space), marginalization over α, β leads to catastrophic biases in the recovered values (up to ∼ 6σ in some numerical tests we performed). This is a strong hint that the naive form of the likelihood function above is incorrect for this problem. The effective likelihood we derive below solves this problem. The standard approach to cosmological parameters fitting outlined above adopted in most of the literature to date has several shortcomings, which can be summarized as follows: (i) The expression for the χ 2 , Eq. (13), has no fundamental statistical justification, but is based on a heuristic derivation. The fundamental problem with Eq. (13) is that some of the parameters being fitted (namely, α, β) control both the location and the dispersion of the χ 2 expression, as they appear both in the numerator and the denominator, via the (σ fit µi ) 2 term. Therefore, the statistical problem is one of jointly estimating both the location and the variance. We show below how this can be tackled using a principled statistical approach.
(ii) Adjusting σ int µ to obtain the desired goodness-of-fit is problematic, as it does not allow one to carry out any further goodness-of-fit test on the model itself, for obviously the variance has been adjusted to achieve a good fit by construction. This means that model checking is by construction not possible with this form of the likelihood function.
(iii) It would be interesting to obtain not just a point estimate for σ int µ , but an actual probability distribution for it, as well. This would allow consistency checks e.g. among different surveys, to verify whether the recovered intrinsic dispersions are mutually compatible (within errorbars). This is currently not possible given the standard χ 2 method. A more easily generalizable approach is desirable, that would allow one to test the hypothesis of multiple SNe populations with different values of intrinsic dispersion, for example as a consequence of evolution with redshift, or correlated with host galaxy properties. Current practice is to split the full SN sample in subsamples (e.g., low-and high-redshift, or for different values of the colour correction) and check for the consistency of the recovered values from each of the subsamples. Our method allows for a more systematic approach to this kind of important model checking procedure.
(iv) It is common in the literature to obtain inferences on the parameters of interest by minimizing (i.e., profiling) over nuisance parameters entering in Eq. (13). This is in general much more computationally costly than marginalization from e.g. MCMC samples (Feroz et al. 2011). There are also examples where some nuisance parameters are marginalized over, while others are maximised (Astier et al. 2006), which is statistically inconsistent and should best be avoided. It should also be noted that maximisation and marginalization do not in general yield the same errors on the parameters of interest when the distribution is non-Gaussian. From a computational perspective, it would be advantageous to adopt a fully Bayesian method, which can be used in conjunction with fast and efficient MCMC and nested sampling techniques for the exploration of the parameter space. This would also allow one to adopt Bayesian model selection methods (which cannot currently be used with the standard χ 2 approach as they require the full marginalization of parameters to compute the Bayesian evidence).
(v) The treatment of systematic errors is being given great attention in the recent literature (see e.g. Nordin et al. (2008)), but the impact of various systematics on the final inference for the interesting parameters, C , has often been propagated in an approximate way (e.g., Kessler et al. (2009a), Appendix F). As we are entering an epoch when SN cosmology is likely to be increasingly dominated by systematic errors, it would be desirable to have a consistent way to include sources of systematic uncertainties in the cosmological fit and to propagate the associated error consistently on the cosmological parameters. The inclusion in the analysis pipeline of systematic error parameters has been hampered so far by the fact that this increases the number of parameters being fitted above the limit of what current methods can handle. However, if a fully Bayesian expression for the likelihood function was available, one could then draw on the considerable power of Bayesian methods (such as MCMC, or nested sampling) which can efficiently handle larger parameter spaces.
Motivated by the above problems and limitations of the current standard method, we now proceed to develop in the next section a fully Bayesian formalism from first principles, leading to the formulation of a new effective likelihood which will overcome, as will be shown below, the above problems. A more intuitive understanding of our procedure can be acquired from the simpler toy problem described in Appendix B.

The Bayesian hierarchical model
We now turn to developing a Bayesian hierarchical model (BHM) for the SNe data from SALT-II lightcurve fits. The same general linear regression problem with unknown variance has been addressed by Kelly (2007), and applied in that paper to X-ray spectral slope fitting. The gist of our method is shown in the graphical network of Fig. 1, which displays the probabilistic and deterministic connection between variables. The fundamental idea is that we introduce a new layer of so-called "latent" variables -that is, quantities which describe the "true" value of the corresponding variables, and which are obviously unobserved (represented in blue in Fig. 1). In particular, we associate with each SN a set of latent variables (zi, ci, x1i, Mi), which give the true value of the redshift, colour correction, stretch correction and intrinsic magnitude for that SN. For the stretch and colour correction, we have noisy estimates in the form of dataĉi,x1i (with an associated covariance matrix) resulting from the SALT-II fits (solid vertical lines). For the redshift zi, we have noisy estimates given by the measured valuesẑi. As above, all observed quantities are denoted by a hat. Each of the latent variables (ci, x1i, Mi) is modeled as an uncorrelated random variable drawn from an underlying parent distribution, described by a Gaussian with two parameters, a population mean and a variance 2 . The latent intrinsic magnitudes Mi are thought of as realizations from a SN population with overall mean intrinsic magnitude M0 and intrinsic dispersion σ int µ (both to be estimated from the data): This is represented in Fig. 1 by the solid arrow connecting the parameters M0, (σ int µ ) 2 with the latent variables Mi. This represents statistically the physical statement that an intrinsic dispersion (σ int µ ) 2 is expected to remain in the distribution of SNe magnitudes, even after applying the Phillips corrections to the observed magnitudes. The simple Ansatz of Eq. (19) can be replaced by a more refined modeling, which describes the existence of multiple SN populations, for example with intrinsic magnitude correlated with host galaxy properties, for which there is a growing body of evidence (Sullivan et al. 2006;Mandel et al. 2010;. We shall investigate this possibility in a forthcoming work.
The parent distribution of the true colour and stretch corrections is similarly represented by Gaussians, again parameterised each by a mean (c , x ) and a variance ( As above, the quantities c , R 2 c , x , R 2 x have to be estimated from the data. The choice of a Gaussian distribution for the latent variables c and x 1 is justified by the fact that the observed distribution ofĉ andx 1 , shown in Fig. 2 for the actual SNIa sample described in section 3.1 below, is fairly well described by a Gaussian. As shown in Fig. 2, there might be a hint for a heavier tail for positive values ofĉ, but this does not fundamentally invalidate our Gaussian approximation. It would be easy to expand our method to consider other distributions, for example mixture models of Gaussians to describe a more complex population or a distribution with heavier tails, if non-Gaussianities in the observed distribution should make such modeling necessary. In this paper we consider the simple uni-modal Gaussians given by Eq. (20).
While we are interested in the value (and distribution) of the SNe population parameters M0, (σ int µ ) 2 , the parameters 2 In this paper we use the notation: to denote a random variable x being drawn from an underlying Gaussian distribution with mean m and variance σ 2 . In vector notation, m is replaced by a vector m, while σ 2 is replaced by the covariance matrix Σ and the distribution is understood to be the appropriate multidimensional Normal, see Appendix A. entering Eqs. (20) are not physically interesting. Therefore they will be marginalized out at the end. It is however crucial to introduce them explicitly in this step: it is precisely the lack of modeling in the distribution of ci that usually leads to the observed "biases" in the reconstruction of β, see Appendix B for an explicit example with a toy model. For ease of notation, we re-write Eq. (19) and (20) in matrix notation as: Having introduced 3n latent (unobserved) variables (c, x 1 , M ), where n is the number of SNe in the sample, the fundamental strategy of our method is to link them to underlying population parameters via Eqs. (19) and (20), then to use the observed noisy estimates to infer constraints on the population parameters of interest (alongside the cosmological parameters), while marginalizing out the unobserved latent variables.
The intrinsic magnitude Mi is related to the observed B-band magnitudem * B and the distance modulus µ by Eq. (5), which can be rewritten in vector notation as: Note that the above relation is exact, i.e. M , x 1 , c are here the latent variables (not the observed quantities), while m * B is the true value of the B-band magnitude (also unobserved). This is represented by the dotted (deterministic) arrows connecting the variables in Fig. 1.
We seek to determine the posterior pdf for the parameters of interest Θ = {C , α, β, σ int µ }, while marginalizing over the uknown population mean intrinsic magnitude, M0. From Bayes theorem, the marginal posterior for Θ is given by (see e.g. Trotta (2008); Hobson et al. (2010) for applications of Bayesian methods to cosmology): where p(D) is the Bayesian evidence (a normalizing constant) and the prior p(Θ, M0) can be written as We take a uniform prior on the variables C , α, β (on a sufficiently large range so as to encompass the support of the likelihood), as well as a Gaussian prior for p(M0|σ int µ ), since M0 is a location parameter of a Gaussian (conditional on σ int µ ). Thus we take where the mean of the prior (Mm = −19.3 mag) is taken to be a plausible value in the correct ballpark, and the variance (σM 0 = 2.0 mag) is sufficiently large so that the prior is very diffuse and non-informative (the precise choice of mean and variance for this prior does not impact on our numerical results). Finally, the appropriate prior for σ int µ is a Jeffreys' prior, i.e., uniform in log σ int µ , as σ int µ is a scale parameter, see e.g. Box & Tiao (1992). Although the intrinsic magnitude M0 and the Hubble constant H0 are perfectly degenerate as far as SNIa data are concerned, we do not bundle them together in a single parameter but treat them separately with distinct priors, as we are interested in separating out the variability due to the distribution of the SNIa intrinsic magnitude.
We now proceed to manipulate further the likelihood, p(D|Θ, M0) = p(ĉ,x 1 ,m * B |Θ, M0): In the first line, we have introduced a set of 3n latent variables, {c, x 1 , M }, which describe the true value of the colour, stretch and intrinsic magnitude for each SNIa. Clearly, as those variables are unobserved, we need to marginalize over them. In (we have also dropped M0 from the likelihood, as conditioning on M0 is irrelevant if the latent M are given). If we further marginalize over M0 (as in Eq. (28), including the prior on M0), the expression for the effective likelihood, Eq. (32), then becomes: where µi ≡ µi(zi, Θ) and we have defined as well as the 3n × 3n block covariance matrix 3 Finally we explicitly include redshift uncertainties in our formalism. The observed apparent magnitude,m * B , on the left-hand-side of Eq. (35), is the value at the observed redshift,ẑ. However, µ in Eq. (35) should be evaluated at the true (unknown) redshift, z. As above, the redshift uncertainty is included by introducing the latent variables z and integrating over them: where we model the redshift errors p(z|ẑ) as Gaussians: with a n × n covariance matrix: It is now necessary to integrate out all latent variables and nuisance parameters from the expression for the likelihood, Eq. (34). This can be done analytically, as all necessary integral are Gaussian. The detailed steps are described in Appendix C. In summary, the procedure consists of: This leads to the final expression for the effective likelihood of Eq. (C32): where the various vectors and matrices appearing in the above expression are defined in Appendix C. This equation is the main result of our paper. The two remaining nuisance parameters Rc, Rx cannot be integrated out analytically, so they need to be marginalized numerically. Hence, they are added to our parameters of interest and are sampled over numerically, and then marginalized out from the joint posterior.

Generation of simulated SNe data
We now proceed to test the performance of our method on simulated data sets, and to compare it with the usual χ 2 approach. We take as starting point the recent compilation of 288 SNIa from Kessler et al. (2009a), the data set on which we apply our Bayesian hierarchical model in the next section. Kessler et al. (2009a) re-fitted the lightcurves from five different surveys:  (Garnavich et al. 1998;Knop et al. 2003;Riess & Strolger 2004, using both the SALT-II method and the MLCS fitter. In the following, we are exclusively employing the results of their SALT-II fits and use those as the observed data set for the purposes of our current work, as described in the previous section. More refined procedures could be adopted, for example by simulating lightcurves from scratch, using e.g. the publicly available package SNANA (Kessler et al. 2009b). In this paper we chose a simpler approach, consisting of simulating SALT-II fit results in such a way to broadly match the distributions and characteristics of the real data set used in Kessler et al. (2009a). The numerical values of the parameters used in the simulation are shown in Table 1. We adopt a vanilla, flat ΛCDM cosmological model as fiducial cosmology. The Phillips correction parameters are chosen to match the best-fit values reported in Kessler et al. (2009a), while the distributional properties of the colour and stretch correction match the observed distribution of their total SN sample. For each survey, we generate a number of SNe matching the observed sample, and we model their redshift distribution as a Gaussian, with mean and variance estimated from the observed distribution within each survey. The observational error of m * B , c, x 1 is again drawn from a Gaussian distribution whose mean and variance have been matched to the observed ones for each survey. Finally, the pseudo-data (i.e., the simulated SALT-II fits results) are generated by drawing from the appropriate distributions centered around the latent variables. For simplicity, we have set to 0 the off-diagonal elements in the correlation matrix (39) in our simulated data, and neglected redshift errors. None of these assumptions have a significant impact on our results. In summary, our procedure for each survey is as follows: (i) Draw a value for the latent redshift zi from a normal distribution with mean and variance matching the observed ones. As we neglect redshift errors in the pseudo-data for simplicity (since the uncertainty in z is subdominant in the overall error budget), we setẑi = zi.
(ii) Compute µi using the fiducial values for the cosmological parameters C and the above zi from Eq. (8).
(iii) Draw the latent parameters x1i, ci, Mi from their respective distributions (in particular, including an intrinsic scatter σ int µ = 0.1 mag in the generation of Mi).
(iv) Compute m * Bi using x1i, ci, Mi and the Phillips relation Eq. (5). (v) Draw the value of the standard deviations σx 1 i, σc i , σm i , from the appropriate normal distributions for each survey type. A small, zi-dependent stochastic linear addition is also made to σx 1 i, σc i , σm i , to mimic the observed correlation between redshift and error.
As shown in Fig. 3, the simulated data from our procedure have broadly similar distributions to the to real ones. The two notable exceptions are the overall vertical shift observed in the distance modulus plot, and the fact that our synthetic data cannot reproduce the few outliers with large values of the variances (bottom panels). The former is a consequence of the different intrinsic magnitude used in our simulated data (as the true one is unknown). However, the intrinsic magnitude is marginalized over at the end, so this difference has no impact on our inferences. The absence of outliers is a consequence of the fact that our simulation is a pure phenomenological description of the data, hence it cannot encapsulate such fine details. While in principle we could perform outlier detection with dedicated Bayesian procedures, we do not pursue this issue further  in this paper. We stress once more that the purpose of our simulations is not to obtain realistic SNIa data. Instead, they should only provide us with useful mock data sets coming from a known model so that we can test our procedure. More sophisticated tests based on more realistically generated data (e.g., from SNANA) are left for future work.
As mentioned above, in keeping with the literature we only consider either flat Universes with a possible w = −1 (the ΛCDM model), or curved Universes with a cosmological constant (w = −1, the wCDM model). Of course it is possible to relax those assumptions and consider more complicated cosmologies with a larger number of free parameters if one so wishes (notably including evolution in the dark energy equation of state).
Of the parameters listed in Eq. (43), the quantities Rc, Rx are of no interest and will be marginalized over. As for the remaining parameters, we are interested in building their marginal 1 and 2-dimensional posterior distributions. This is done by plugging the likelihood (43) into the posterior of Eq. (28), with priors on the parameters chosen according to Table 2. We use a Gaussian prior on the Hubble parameter H0 = 72 ± 8 km/s/Mpc from local determinations of the Hubble constant (Freedman et al. 2001). However, as H0 is degenerate with the intrinsic population absolute magnitude M0 (which is marginalized over at the end), replacing this Gaussian prior with a less informative prior H0[km/s/Mpc] ∼ U(20, 100) has no influence on our results.
Numerical sampling of the posterior is carried out via a nested sampling algorithm (Skilling 2004(Skilling , 2006Feroz & Hobson 2008;Feroz et al. 2009). Although the original motivation for nested sampling was to compute the Bayesian evidence, the recent development of the MultiNest algorithm (Feroz & Hobson 2008;Feroz et al. 2009) has delivered an extremely powerful and versatile algorithm that has been demonstrated to be able to deal with extremely complex likelihood surfaces in hundreds of dimensions exhibiting multiple peaks. As samples from the posterior are generated as a by-product of the evidence computation, nested sampling can also be used to obtain parameter constraints in the same run as computing the Bayesian evidence. In this paper we adopt the publicly available MultiNest algorithm (Feroz & Hobson 2008) to obtain samples from the posterior distribution of Eq. (28). We use 4000 live points and a tolerance parameter 0.1, resulting in about 8 × 10 5 likelihood evaluations. 4 We also wish to compare the performance of our BHM with the usually adopted χ 2 minimization procedure. To this end, we fit the pseudo-data using the χ 2 expression of Eq. (13). In order to mimic what is done in the literature as closely as possible, we first fix a value of σ int µ . Then, we simultaneously minimize the χ 2 w.r.t. the fit parameters ϑ = {Ωm, Ωκ or w, H0, M0, α, β}, as described below. We then evaluate the χ 2 /dof from the resulting best fit point, and we adjust σ int µ to obtain χ 2 /dof = 1. We then repeat the above minimization over ϑ for this new value of σ int µ , and iterate the procedure. Once we have obtained the global best fit point, we derive 1-and 2-dimensional confidence intervals on the parameters by profiling (i.e., maximising over the other parameters) over the likelihood with χ 2 given by Eq. (13). According to Wilks' theorem, approximate confidence intervals are obtained from the profile  Figure 4. Reconstruction of cosmological parameters from a simulated data set encompassing 288 SNIa, with characteristics matching presently available surveys (including realization noise). Blue regions contain 95% and 68% of the posterior probability (other parameters marginalized over) from our BHM method, the red contours delimit 95% and 68% confidence intervals from the standard χ 2 method (other parameters maximised). The yellow star indicates the true value of the parameters. The left panel assumes w = −1 while the right panel assumes Ωκ = 0. Notice how out method produced considerably less biassed constraints on the parameters.
likelihood as the regions where the χ 2 increases by ∆χ 2 from its minimum value, where ∆χ 2 can be computed from the chi-square distribution with the number of degree of freedoms corresponding to the number of parameters of interest and is given in standard look-up tables.
Obtaining reliable estimates of the profile likelihood using Bayesian algorithms (such as MultiNest) is a considerably harder numerical task than mapping out the Bayesian posterior. However, it has been shown that MultiNest can be successfully used for this task even in highly challenging situations (Feroz et al. 2011), provided the number of live points and tolerance value used are adjusted appropriately. For our χ 2 scan, we adopt 10 4 live points and a tolerance of 0.1. We have found that those values give accurate estimates of the profile likelihood more than 2σ into the tails of the distribution for an 8 dimensional Gaussian toy model (whose dimensionality matches the case of interest here). With these MultiNest settings, we gather 1.5 × 10 5 samples, from which the profile likelihood is derived.
Our implementation of the χ 2 method is designed to match the main features of the fitting procedure usually adopted in the literature (namely, maximisation of the likelihood rather than marginalization of the posterior, and iterative determination of the intrinsic dispersion), although we do not expect that it exactly reproduces the results obtained by any specific implementation. Its main purpose is to offer a useful benchmark against which to compare the performance of our new Bayesian methodology.

Parameter reconstruction
We compare the cosmological parameters reconstructed from the standard χ 2 method and our Bayesian approach in Fig. 4 for a typical data realization. The left-hand-side panel shows constraints in the Ωm − ΩΛ plane for the ΛCDM model, both from our Bayesian method (filled regions, marginalized 68% and 95% posterior) and from the standard χ 2 method (red contours, 68% and 95% confidence regions from the profile likelihood). In the right-hand-side panel, constraints are shown in the w − Ωm plane instead (this fit assumes a flat Universe). In a typical reconstruction, our Bayesian method produced considerably tighter constraints on the cosmological parameters of interest than the usual χ 2 approach. Our constraints are also less biassed w.r.t. the true value of the parameters, an important advantage that we further characterize below.
Our BHM further produces marginalized posterior distributions for all the other parameters of the fit, including the stretch and colour corrections and the intrinsic dispersion of the SNe. The 1D marginal posteriors for those quantities are shown in Fig. 5. The recovered posterior means lie within 1σ of the true values. Notice that we do not expect the posterior mean to match exactly the true value, because of realization noise in the pseudo-data. However, as shown below, our method delivers less biassed estimates of the parameters, and a reduced mean squared error compared with the standard χ 2 approach. The stretch correction α is determined with 8% accuracy, while the colour correction parameter β is constrained with an accuracy better than 3%. A new feature of our method is that it produces a posterior distribution for the SN population intrinsic dispersion, σ int µ (right-hand-side panel of Fig 5). This allows to determined the intrinsic dispersion of the SNIa population to typically about 10% accuracy.

Comparison of long-term performance of the two methods
Encouraged by the good performance of our BHM on single instances of simulated data, we now wish to compare the long-term performance of the two methods for a series of simulated data realizations. We are interested in comparing the average ability of both methods to recover parameter values that are as much as possible unbiased with respect to their true values, as well as to establish the coverage properties of the credible and confidence intervals. Coverage is defined as the probability that an interval contains (covers) the true value of a parameter, in a long series of repeated measurements. The defining property of a e.g. 95% frequentist confidence interval is that it should cover the true value 95% of the time; thus, it is reasonable to check if the intervals have the properties they claim. Coverage is a frequentist concept: intervals based on Bayesian techniques are meant to contain a given amount of posterior probability for a single measurement (with no reference to repeated measurements) and are referred to as credible intervals to emphasize the difference in concept. While Bayesian techniques are not designed with coverage as a goal, it is still meaningful to investigate their coverage properties. To our knowledge, the coverage properties of even the standard χ 2 method (which, being a frequentist method would ideally be expected to exhibit exact coverage) have never been investigated in the SN literature.
We generate 100 realizations of the pseudo-data from the fiducial model of Table 1 as described in section 3.1, and we analyze them using our BHM method and the standard χ 2 approach, using the same priors as above, given in Table 2. For each parameter of interest θ, we begin by considering the relative size of the posterior 68.3% range from out method, σ BHM θ , compared with the 68.3% confidence interval from the χ 2 method, σ χ 2 θ , which is summarized by the quantity S θ which shows the percentage change in errorbar size with respect to the errorbar derived using the χ 2 method A value S θ < 1 means that our method delivers tighter errorbars on the parameter θ. A histogram of this quantity for the variables of interest is shown in Fig. 6, from which we conclude that our method gives smaller errobars on Ωm, ΩΛ and w in almost all cases. The flip side of this result is that the uncertainty on α, β is larger from our method than from the χ 2 approach in essentially all data realization. This is a consequence of the large number of latent variables our method marginalizes over. Tight errorbars are good, but not if they come at the expense of a biased reconstruction. To investigate this aspect, we build the following test statistics from each reconstruction: where θBHM is the posterior mean recovered using our BHM method, θ bf χ 2 is the best-fit value for the parameter recovered using the standard χ 2 approach and θtrue is the true value for that parameter. The meaning of T θ is the following: for a given data realization, if the reconstructed posterior mean from our BHM is closer to the true parameter value than the best-fit χ 2 , then T θ < 0, which means that our method is less biassed than χ 2 . A histogram of the distribution of T θ across the 100 realizations, shown in Fig. 7, can be used to compare the two methods: a negative average in the histogram means that the BHM outperforms the usual χ 2 . For all of the parameters considered, our method is clearly much less biassed than the χ 2 method, outperforming χ 2 about 2/3 of the time. Furthermore, the reconstruction of the intrinsic dispersion is better with our Bayesian method almost 3 times out of 4. We emphasize once more that our methodology also provides an estimate of the uncertainty in the intrinsic dispersion, not just a best-fit value as the χ 2 approach.
We can further quantify the improvement in the statistical reconstruction by looking at the bias and mean squared error (MSE) for each parameter, defined as bias = θ − θtrue (47)  Table 3. Comparison of the bias and mean squared error for our Bayesian method and the usual χ 2 approach. The columns labelled "Improvement" give the factor by which our Bayesian method reduces the bias and the MSE w.r.t. the χ 2 approach.  respectively, where the expectation is taken by averaging over the observed values in our 100 simulated trials,θ = θBHM (θ = θ bf χ 2 ) for the BHM (for the χ 2 approach) and Var is the observed parameter variance. The bias is the expectation value of the difference between estimator and true value, while the MSE measures the average of the squares of the errors, i.e., the amount by which the estimator differs from the true value for each parameter. Obviously, a smaller bias and a smaller MSE imply a better performance of the method. The results for the two methods are summarized in Table 3, which shows how our method reduces the bias by a factor ∼ 2 − 3 for most parameters, while reducing the MSE by a factor of ∼ 2. The only notable exception is the bias of the EOS parameter w, which is larger in our method than in the χ 2 approach.
Finally, in Fig. 8 we plot the coverage of each method for 68% and 95% intervals. Errobars give an estimate of the uncertainty of the coverage result, by giving the binomial sampling error from the finite number of realizations considered, evaluated from the binomial variance as N p(1 − p), where N = 100 is the number of trials and p is the observed fractional coverage. Both method slightly undercover, i.e. the credible region and confidence intervals are too short, although the lack of coverage is not dramatic: e.g., the typical coverage of the 1σ (2σ) intervals from our method is ∼ 60% (90%). Our method shows slightly better coverage properties than the χ 2 method, while producing considerably tighter and less biassed constraints (as demonstrated above). This further proves that the tighter intervals recovered by our method do not suffer from bias w.r.t the true values.

COSMOLOGICAL CONSTRAINTS FROM CURRENT SNIA DATA
We now apply our BHM to fitting real SN data. We use the SALT-II fits result for 288 SNIa from Kessler et al. (2009a), which have been derived from 5 different surveys. Our method only includes statistical errors according to the procedure described in  Figure 7. Histograms of the test statistics defined in Eq. (46), comparing the long-term performance of the two methods for the parameters of interest in the ΛCDM model (left) and the wCDM model (right). A predominantly negative value of the test statistics means that our method gives a parameter reconstruction that is closer to the true value than the usual χ 2 , i.e., less biassed. For the cosmological parameters (top row), our method outperforms χ 2 about 2 times out of 3.  Figure 8. Coverage of our method (blue) and standard χ 2 (red) for 68% (solid) and 95% (dashed) intervals, from 100 realizations of pseudo-data for the ΛCDM model (left) and the wCDM model (right). While both methods show significant undercoverage for all parameters, our method has a comparable coverage to the standard χ 2 , except for w. Coverage values for the intrinsic dispersion σ int µ are not available from the χ 2 method, as it does not produce an error estimate for this quantity. section 2.3, coming from redshift uncertainties (arising from spectroscopic errors and peculiar velocities), intrinsic dispersion (which is determined from the data) and full error propagation of the SALT-II fit results. Systematic uncertainties play an important role in SNIa cosmology fitting, and (although not included in this study) can also be treated in our formalism in a fully consistent way. We comment on this aspect further below, though we leave a complete exploration of systematics with our BHM to a future, dedicated work (March et al. 2011). We show in Fig. 9 the constraints on the cosmological parameters Ωm − ΩΛ (left panel, assuming w = −1) and w − Ωm (right panel, assuming Ωκ = 0) obtained with our method. All other parameters have been marginalized over. In order to be consistent with the literature, we have taken a non-informative prior on H0, uniform in the range [20, 100] km/s/Mpc. The figure also compares our results with the statistical contours from Kessler et al. (2009a), obtained using the χ 2 method. (Notice that we compare with the contours including only statistical uncertainties for consistency.) In Fig. 10 we combine our SNIa constraints with Cosmic Microwave Background (CMB) data from WMAP 5-yrs measurements (Komatsu et al. 2009) and Baryonic Acoustic Oscillations (BAO) constraints from the Sloan Digital Sky Survey LRG sample (Eisenstein et al. 2005), using the same method as Kessler et al. (2009a). The combined SNIa, CMB and BAO statistical constraints result in Ωm = 0.28 ± 0.02, ΩΛ = 0.73 ± 0.01 (for the ΛCDM model) and Ωm = 0.28 ± 0.01, w = −0.90 ± 0.05 (68.3% credible intervals) for the wCDM model. Although the statistical uncertainties are comparable to the results by Kessler et al. (2009a) from the same sample, our posterior mean values present shifts of up to ∼ 2σ compared to the results obtained using the standard . Constraints on the cosmological parameters Ωm, Ω Λ (left panel, assuming w = −1) and w, Ωm (right panel, assuming Ωκ = 0) from our Bayesian method (light/dark blue regions, 68% and 95% marginalized posterior), compared with the statistical errors from the usual χ 2 approach (yellow/red regions, same significance level; from Kessler at al. (2009a)). The yellow star gives the posterior mean from our analysis.
χ 2 approach. This is a fairly significant shift, which can be attributed to our improved statistical method, which exhibits a reduced bias w.r.t. the χ 2 approach. Fig. 11 shows the 1d marginalized posterior distributions for the Phillips correction parameters and for the intrinsic dispersion. All parameters are well constrained by the posterior, and we find α = 0.12 ± 0.02, β = 2.7 ± 0.1 and a value of the intrinsic dispersion (for the whole sample) σ int µ = 0.13 ± 0.01 mag. Kessler et al. (2009a) find values for the intrinsic dispersion ranging from 0.08 (for SDSS-II) to 0.23 (for the HST sample), but their χ 2 method does not allow them to derive an error on those determinations. With our method, it would be easy to derive constraints on the intrinsic dispersion of each survey -all one needs to do is to replace Eq. (19) with a corresponding expression for each survey. This introduces one pair of population parameters (M0, σ int µ ) for each survey. In the same way, one could study whether the intrinsic dispersion evolves with redshift. We leave a detailed study of these issues to a future work.
The value of α found in Kessler et al. (2009a) is in the range 0.10 − 0.12, depending of the details of the assumptions made, with a typical statistical uncertainty of order ∼ 0.015. These results are comparable with our own. As for the colour correction parameter β, constraints from Kessler et al. (2009a) vary in the range 2.46 − 2.66, with a statistical uncertainty of order 0.1 − 0.2. This stronger dependence on the details of the analysis seems to point to a larger impact of systematic uncertainties for β, which is confirmed by evidence of evolution with redshift of the value of β (Kessler et al. (2009a) , Fig. 39). Our method can be employed to carry out a rigorous assessment of the evolution with redshift of colour corrections. A possible strategy would be to replace β with a vector of parameters β1, β2, . . . , with each element describing the colour correction in a different redshift bin. The analysis proceeds then as above, and it produces posterior distributions for the components of β, which allows to check the hypothesis of evolution. Finally, in such an analysis the marginalized constraints on all other parameters (including the cosmological parameters of interest) would automatically include the full uncertainty propagation from the colour correction evolution, without the need for further ad hoc inflation of the errorbars. These kind of tests will be pursued in a forthcoming publication (March et al. 2011).

CONCLUSIONS
We have presented a statistically principled approach for the rigorous analysis of SALT-II SNIa lightcurve fits, based on a Bayesian hierarchical model. The main novelty of our method is that it produces an effective likelihood that propagates uncertainties in a fully consistent way. We have introduced an explicit statistical modeling of the intrinsic magnitude distribution of the SNIa population, which for the first time allows one to derive a full posterior distribution of the SNIa intrinsic dispersion.
We have tested our method using simulated data sets and found that it compares favourably with the standard χ 2 approach, both on individual data realizations and in the long term performance. Statistical constraints on cosmological parameters are significantly improved, while in a series of 100 simulated data sets our method outperforms the χ 2 approach at least 2 times out of 3 for the parameters of interest. We have also demonstrated that our method is less biassed and has better coverage properties than the usual approach.
We applied our methodology to a sample of 288 SNIa from multiple surveys. We find that the flat ΛCDM model is still in good agreement with the data, even under our improved analysis. However, the posterior mean for the cosmological  parameters exhibit up to 2σ shifts w.r.t. results obtained with the conventional χ 2 approach. This is a consequence of our improved statistical analysis, which benefits from a reduced bias in estimating the parameters. While in this paper we have only discussed statistical constraints, our method offers a new, fully consistent way of including systematic uncertainties in the fit. As our method is fully Bayesian, it can be used in conjunction with fast and efficient Bayesian sampling algorithms, such as MCMC and nested sampling. This will allow to enlarge the number of parameters controlling systematic effects that can be included in the analysis, thus taking SNIa cosmological parameter fitting to a new level of statistical sophistication. The power of our method as applied to systematic errors analysis will be presented in a forthcoming, dedicated work.
At a time when SNIa constraints are entering a new level of precision, and with a tenfold increase in the sample size expected over the next few years, we believe it is timely to upgrade the cosmological data analysis pipeline in order to extract the most information from present and upcoming SNIa data. This paper represents a first step in this direction. Wang L., et al., 2003, Astrophys. J., 590, 944 Conley A., 2006, Astrophys. J., 644,1 Blondin S., Mandel K. S., Kirshner, R. P.,2011, Astron. Astrophys.,526,81 Kim A. G.,2011, Pub. Am. Soc. Pacific., 123 900 230 Gopal Vishwakarma R., Narlikar J. V., 2011, R. Astron. Astrop.,10,1195 APPENDIX A: NOTATION For reference, we collect here a few useful formulas relating to Gaussian integrals. Use the notation x ∼ Nx(µ, Σ) to denote that the random variable x is drawn from a Normal distribution of mean µ and inverse covariance matrix Σ, given by In performing Gaussian integrals, it is also convenient to recall that: Finally, the integral of the multidimensional Gaussian of Eq. (A1) over all space is unity.

APPENDIX B: TOY MODEL FOR GENERAL LINEAR FITTING
We can get an intuitive understanding of the central ingredient in our Bayesian hierarchical method by considering a simpler toy model, which highlights the salient features of the problem. As mentioned in section 2.2, several authors have observed in the past that the spread of fit values for the colour correction, c, is of the same order as the size of the statistical uncertainty on the values themselves. This leads to a bias in the best-fit value of β obtained from minimizing the χ 2 in Eq. (13). The quantity β gives the slope of the linear relationship between c and µ, see Eq. (27). In the usual χ 2 approach, the latent (true) c are replaced by the observed value,ĉ, as in Eq. (14). If the statistical uncertainty in the independent variable,ĉ, is as large as the spread of its values, the best-fit slope obtained from minimizing the χ 2 may be biased, as the large uncertainty in the actual location of the independent variable leads to confusion as to the value of the slope. The (Bayesian) solution to this problem is to determine the spread of the independent variable directly from the data, and to marginalize over it with an appropriate prior. This gives the most general solution to the problem of linear fitting with errors in both the independent and dependent variable, as shown by Gull (1989). Recent literature in cosmology and astronomy (D'Agostini 2005; Hogg et al. 2010) addresses linear fitting, but not this general case, which has been treated before in Kelly (2007) (see also Andreon (2006); Andreon & Hurn (2010) for related examples). In this short appendix, we will analyse the toy model of fitting a linear function, and compare the performance of a BHM and the χ 2 approach.

B1 Bayesian linear fitting in the presence of large x and y uncertainties
In this subsection, we give a short review of the results of Gull (1989). The simplest toy model which illustrates the methodology we adopt in our paper is shown by the graphical network of Fig. B1. A linear relationship is assumed between the latent variables xi and yi, described by a slope a and intercept b (which are collectively denoted as a parameter vector θ): The observed values for the dependent and independent variables are denoted by hats (xi,ŷi, i = 1, . . . , N ), and they are obtained from the latent values under the assumption of Gaussian noise (with known variances σ 2 x , σ 2 y ): xi ∼ N (xi, σ 2 x ) andŷi ∼ N (yi, σ 2 y ). This probabilistic relationship is depicted in Fig. B1 by the solid arrows connecting the latent variables to the observed quantities. Assuming that errors are uncorrelated, the joint likelihood is given by 5 p(x,ŷ|x, y, σx, σy, θ) = (4π 2 σ 2 The problem can be made more symmetric by defining rescaled versions of the data: where the variables x0, y0 describe the mean value ofx,ŷ, while Rx, Ry describe their spread. These new variables are related to the old ones (a, b) by b = y0 − ax0 and a = Ry/Rx.
Neglecting the normalization constant, the joint posterior for x, y, x0, y0, Rx, Ry can be written as: p(x, y, x0, y0, Rx, Ry|x,ŷ, σx, σy) ∝ p(x,ŷ|x, x0, y0, Rx, Ry)p(x|x0, y0, Rx, Ry)p(x0, y0, Rx, Ry), where the first term on the r.h.s. is the likelihood of Eq. (B3). The key step is to recognize that the appropriate conditional distribution for the latent x is This describe a prior over x centered around the hyperparameter x0 and with standard deviation Rx. Crucially, both x0 and Rx are unknown, and are explicitly determined from the data in the joint posterior, before being marginalized out at the end. Finally, for the prior p(x0, y0, Rx, Ry) appearing in Eq. (B6) we adopt a uniform prior on x0, y0 (as those are location variables) and a prior uniform in log Rx, log Ry (those being scale variables, as apparent from (B7)). From here, one can eliminate y from Eq. (B6) using Eq. (B1), then trade b for Ry using Eq. (B5), to obtain: (B8) From this expression, the latent x can be marginalized out analytically, as well as the nuisance parameters x0, y0, by using appropriate completions of the square in the Gaussian. After some algebra, one finally obtains p(a, log R|x,ŷ, σx, σy) ∝ (a 2 σ 2 where: In the above,x = ix i/N , and similarly forȳ. The marginal posterior for the slope a is obtained by numerical marginalization of log R from the above expression.

B2 Comparison with the χ 2 approach
If instead of the statistically principled solution found above, one writes down a simple χ 2 expression for the likelihood, including error propagation from the linear relationship (B1) one would obtain (D'Agostini 2005): from where the intercept b is eliminated by maximising over it (profiling). We now compare the reconstruction of the slope parameter a from Eq. (B11) (with b eliminated via profiling) with the result obtained using the Bayesian expression, (B9), marginalized numerically over log R (with a uniform prior on log R) with the help of simulated data. Figs. B2 and B3 show in the upper left panel the simulated data points (N = 300), with the error bars giving the size of the standard deviation in the x and y directions for each datum (i.e., the value of σx, σy). The contour plots depict joint posterior regions for (log a, log R) from the Bayesian expression of Eq. (B9), and confidence regions from the χ 2 likelihood (B11) in the a, b plane. The lower right panel shows the 1d marginalized posterior distribution for log a (blue) and the profile likelihood from the χ 2 approach (black). Fig. B2 shows that the two results are essentially identical when the size of the statistical error is smaller than the spread of values of the data points (in this particular example, by a factor of 2). However, when the statistical uncertainty in the x direction is as large as or larger than the spread of the points (Fig. B3), the χ 2 approach gives a biassed result for the slope parameter a. The Bayesian posterior, on the contrary, is closer to the true value, although it (correctly) shows a larger uncertainty.
With this approximation we can now carry out the multi-dimensional integral of Eq. (40), obtaining  Figure B3. As in Fig. B3, but now with a larger statistical uncertainty in the data compared to their spread. The Bayesian marginal posterior for the slope (blue, lower right panel) peaks much closer to the true value than the χ 2 expression (black). In the lower left panel, the 3σ contour from the Bayesian method lies outside the range of the plot.
where from now on, µ = µ(ẑ) and Σm = ΣC + f Σzf T Strictly speaking, one should integrate over redshift in the range 0 zi < ∞, not −∞ < zi < ∞, which would result in the appearance of Gamma functions in the final result. However, as long as σz i z i 1 (as is the case here), this approximation is expected to be excellent.

C3 Integration over population variables {c , x , M0}
The priors on the population variables b = {c , x , M0} in Eq. (C24) can be written as: Thus Eq. (C24) can be written as: p(ĉ,x 1 ,m * B |Θ) = dRc dRx db |2πΣC | − 1 2 |2πΣP | − 1 2 |2πΣA| 1 2 |2πΣ0| − 1 2 p(Rc)p(Rx) where We can now carry out the Gaussian integral over b in Eq. (C24), obtaining our final expression for the effective likelihood, where we have chosen an improper Jeffreys' prior for the scale variables Rc, Rx: and analogously for Rx. These two remaining nuisance parameters cannot be integrated out analytically, so they need to be marginalized numerically. Hence, Rc, Rx are added to our parameters of interest and are sampled over numerically, and then marginalized out from the joint posterior.