zBEAMS: A unified solution for supernova cosmology with redshift uncertainties

TL;DR

zBEAMS is a hierarchical Bayesian method that accurately estimates cosmological parameters from supernova data with uncertain redshifts, photometric redshift errors, and host galaxy misidentification, without spectra.

Contribution

It introduces a unified Bayesian framework that marginalizes over redshift uncertainties and contamination, improving cosmological estimates from photometric supernova surveys.

Findings

Unbiased cosmological parameter estimation with photometric data.

Effective correction for host galaxy misidentification.

Demonstrates robustness through supernova simulations.

Abstract

Supernova cosmology without spectra will be an important component of future surveys such as LSST. This lack of supernova spectra results in uncertainty in the redshifts which, if ignored, leads to significantly biased estimates of cosmological parameters. Here we present a hierarchical Bayesian formalism -- zBEAMS -- that addresses this problem by marginalising over the unknown or uncertain supernova redshifts to produce unbiased cosmological estimates that are competitive with supernova data with spectroscopically confirmed redshifts. zBEAMS provides a unified treatment of both photometric redshifts and host galaxy misidentification (occurring due to chance galaxy alignments or faint hosts), effectively correcting the inevitable contamination in the Hubble diagram. Like its predecessor BEAMS, our formalism also takes care of non-Ia supernova contamination by marginalising over the unknown supernova type. We illustrate this technique with simulations of supernovae with photometric redshifts and host galaxy misidentification. A novel feature of the photometric redshift case is the important role played by the redshift distribution of the supernovae.

Equations

$$\mu (z)=m-M=5{\text{log}}_{10}\left(\frac{d_L}{1Mpc}\right)+25,$$

$$d_L(z)=\frac{c(1+z)}{H_0\sqrt{-{\Omega }_k}}\text{sin}\left(H_0\sqrt{-{\Omega }_k}\int \frac{d{z}^{\prime }}{H(z^{\prime })}\right),$$

$$H(z)=H_{0}\bigg{(}\Omega_{m}(1+z)^{3}+\Omega_{\rm DE}(1+z)^{3(1+w)}+\Omega_{k}(1+z)^{2}\bigg{)}^{1/2},$$

$$P(\mathbf{\theta }|\{D_i\})=\frac{P(\{D_i\}|\mathbf{\theta })P(\mathbf{\theta })}{P(\{D_i\})},$$

$$P(\{D_{i}\}|\boldsymbol{\theta})=\prod_{i}^{N}\,\frac{1}{\sqrt{2\pi{\sigma^{i}_{\mu}}^{2}}}\rm{exp}\bigg{(}-\frac{(\mu_{\rm obs}^{i}-\mu(z^{i},\boldsymbol{\theta}))^{2}}{2{\sigma^{i}_{\mu}}^{2}}\bigg{)}\,.$$

$$P(\mathbf{\theta },\mathbf{\phi }|\mathbf{C})=P(\mathbf{\theta }|\mathbf{\phi },\mathbf{C})P(\mathbf{\phi }|\mathbf{C}).$$

$$P(\mathbf{\theta }|\mathbf{C})=\int P(\mathbf{\theta },\mathbf{\phi }|\mathbf{C})d\mathbf{\phi }.$$

$$P(\mathbf{\theta }|D)$$

$$\propto$$

$$P(\mathbf{\theta }|D)$$

$$\propto$$

$$P(\mathbf{\theta }|D)$$

$$\propto$$

$$P(\mathbf{\theta }|D)$$

$$P(z)$$

$$P(\mathbf{\theta }|D)$$

$$P(\mathbf{\theta }|D)\propto P(\mathbf{\theta })\int \frac{1}{2\pi \sqrt{det|C|}}\exp \left(-\frac{1}{2}{\Delta }^T{C}^{-1}\Delta \right)P(z)dz,$$

$$P(\mathbf{\theta }|D)\propto P(\mathbf{\theta })\int \frac{1}{2\pi {\sigma }_z{\sigma }_{\mu }}\exp \left(-\frac{(z_{\mathrm{obs}}-z{)}^2}{2{\sigma }_z^2}\right)\exp \left(-\frac{({\mu }_{\mathrm{obs}}-\mu (z,\mathbf{\theta }){)}^2}{2{\sigma }_{\mu }^2}\right)P(z)dz.$$

$$P(\mathbf{\theta }|D)\propto P(\mathbf{\theta })\int \exp \left(-\frac{(z_{\mathrm{obs}}-z{)}^2}{2{\sigma }_z^2}-\frac{({\tilde{\mu }}_{\mathrm{obs}}-{\mu }^{\prime }z{)}^2}{2{\sigma }_{\mu }^2}\right)\exp \left(-\frac{(z-\bar{z}{)}^2}{2{\sigma }_p^2}\right)dz,$$

$$P(\mathbf{\theta }|D)\propto \frac{P(\mathbf{\theta })}{C^{1/2}}\exp \left(-\frac{1}{2C}\left[{\sigma }_{\mu }^2(z_{\mathrm{obs}}-\bar{z}{)}^2+{\sigma }_p^2({\mu }_{\mathrm{obs}}-\mu (z_{\mathrm{obs}}){)}^2+{\sigma }_z^2({\tilde{\mu }}_{\mathrm{obs}}-{\mu }^{\prime }\bar{z}{)}^2\right]\right),$$

$$C={\sigma }_{\mu }^2({\sigma }_p^2+{\sigma }_z^2)+{\mu }^{\prime 2}{\sigma }_z^2{\sigma }_p^2.$$

$$P(\mathbf{\theta }|D)$$

$$\propto$$

$$P(\mathbf{\theta }|D)$$

$$P(\mathbf{\theta }|D)\propto P(\mathbf{\theta })\left[P_{\mathrm{Ia}}P(D|z_{\ast },{\mu }_{\mathrm{Ia}}(z),\tau =\mathrm{Ia},\mathbf{\theta })+(1-{\mathrm{P}}_{\mathrm{Ia}})\mathrm{P}(D|{\mathrm{z}}_{\ast },{\mu }_{\cancel{\mathrm{Ia}}}(\mathrm{z}),\tau =\cancel{\mathrm{Ia}},\mathbf{\theta })\right].$$

$$P(D|z_{\ast },{\mu }_{\mathrm{Ia}}(z),\tau =\mathrm{Ia},\mathbf{\theta })=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{Ia}}}\exp \left(\frac{-({\mu }_{\mathrm{obs}}-{\mu }_{\mathrm{Ia}}({\mathrm{z}}_{\ast },\mathbf{\theta }){)}^2}{2{\sigma }_{\mathrm{Ia}}^2}\right).$$

$$P(z)=\sum _{\gamma }P(z|\gamma )P(\gamma ).$$

$$P(\mathbf{\theta }|D)$$

$$\propto$$

Full text

1 Introduction

Studies of Type Ia supernovae led to the dark energy breakthrough and modern concordance cosmology, but one can argue that they have been supplanted by Baryon Acoustic Oscillations [2, 3] as the most precise way of constraining cosmology today. To be competitive in the era of LSST [4], EUCLID and the SKA, supernova cosmology faces several big challenges. One is that better control of systematics is required and a number of sophisticated approaches are being developed to improve the control of systematics (e.g. [5, 6, 7, 8, 9, 10, 11, 12]). Another critical problem is that next-generation supernova surveys will be severely spectroscopy limited: LSST will deliver over $10^{5}$ Type Ia Supernova (SNIa) candidates with photometric lightcurves only. The lack of spectroscopy introduces a number of challenges. First, the true identity of any candidate without spectroscopic follow-up is ambiguous - photometric colours only provide a probability for an object to be a SNIa, as opposed to a Type Ibc or II supernova or other transient [13]. Secondly, the precise redshifts of the supernovae are unknown. Photometric redshifts are fairly good if the candidates are known to be SNIa, yielding RMS errors of $\sigma_{z}\sim 0.04(1+z)$, depending on exact assumptions [4, 14, 15]. The problem is that we are exactly in the case where we are not sure whether each candidate is a SNIa or not, and the photometric redshift error is much larger if the object is not a SNIa [16], precisely because they are not standard candles.

A promising approach that dates back to the SDSS II supernova survey [17, 18, 19], is to obtain spectroscopic redshifts for the host galaxies of the supernova candidates and use this as a proxy for the supernova redshift. This will be particularly attractive in the era of big redshift surveys such as 4MOST, SKA and Euclid, where huge numbers of galaxy redshifts will be known. This has the potential to help remove biases [19] and yield improved constraints [20].

However, even this approach has a serious problem: identifying the host galaxy is also not unambiguous. The supernova can appear to lie in between two or more galaxies or may live in a host that is too faint to be detectable (“hostless”) (see figure 1). In general we therefore can assume that instead we have probabilities for the supernova to belong to each of the nearby galaxies on the sky or to be hostless. Current matching algorithms can accurately match the correct host galaxy about $91\%$ of the time when applied to data, potentially increasing to $97\%$ by using machine learning techniques [21]. However, even a $3\%$ contamination may cause significant biases on cosmological parameter inference and must be dealt with.

Our goal with this paper is to develop a single formalism that solves this set of problems, simultaneously handling both the contamination from non-Ia supernovae and contamination from incorrect host identification and redshift assignment. The formalism we desire will be rigorous without resorting to cuts, which are statistically suboptimal.

In laying out the solution we will develop intuition by building systematically to the complete solution. In section 2 we review inference using standard cosmology. In section 3 we consider the case of a single, known Type Ia supernova with a photometric redshift derived from the supernova lightcurve. In section 3.1 we extend the analysis to include host galaxy redshifts. In section 4 we review the BEAMS formalism to handle contamination from non-Ia. In section 5 we present the combined formalism to handle all contamination and finally in section 6 we present an illustrative set of simulations demonstrating our method.

2 Standard Supernova Cosmology Inference

Traditionally, supernova cosmological analysis proceeds with a sample of spectroscopically confirmed type Ia supernovae with well-measured redshifts. The goal is to determine the posterior distribution, $P(\boldsymbol{\theta}|\{D_{i}\})$, over the cosmological parameters, $\boldsymbol{\theta}$, given the dataset which we denote $\{D_{i}\}$.111In much of our discussion we will use $D$ to denote a single supernova so for clarity we use $\{D_{i}\}$ to make explicit the case where we are considering multiple supernovae. Here, the data are the redshift, $z_{i}$, and distance modulus, $\mu_{i}$, for each supernova. The distance modulus for a type Ia supernova is usually estimated from its observed light curve using the SALT2 model [22].

The distance modulus is defined as:

[TABLE]

where $m$ is the apparent magnitude of the object, $M$ is the absolute magnitude of the object and $d_{L}$ is the luminosity distance to the object in Mpc. In a $\Lambda$CDM universe, the luminosity distance is related to the cosmological parameters by:

[TABLE]

where

[TABLE]

and $H_{0}$ is the Hubble constant, $\Omega_{m}$ is the energy density of matter, $\Omega_{\rm DE}$ is the energy density of the dark energy, $w$ is the dark energy equation of state where $w=-1$ corresponds to $\Lambda$, the cosmological constant. Finally $\Omega_{k}$ is the curvature parameter. Collectively, we refer to these cosmological parameters as $\boldsymbol{\theta}$.

The posterior probability distribution over the cosmological parameters is then given by Bayes’ theorem:

[TABLE]

where $P(\{D_{i}\}|\boldsymbol{\theta})$ is the likelihood, $P(\boldsymbol{\theta})$ is the prior and $P(\{D_{i}\})$ is called the Bayesian evidence. Since we are not interested in model selection $P(\{D_{i}\})$ is irrelevant for parameter inference and we drop it in all subsequent analysis.

In the case of uncorrelated Gaussian errors, $\sigma^{i}_{\mu}$, on the estimated distance moduli $\mu_{\rm obs}^{i}$, the likelihood over the $N$ supernovae is:

[TABLE]

A standard cosmological analysis would then proceed by applying a numerical sampling method such as MCMC [23, 24] to determine the full posterior of the cosmological parameters. Care needs to be taken around the parameters involved in estimating the distance modulus, such as considered in [25]. While such caveats can be handled easily by introducing new latent parameters, for simplicity we omit these parameters and assume the distance moduli can be measured directly since it is not core to the problem we are addressing.

3 Inference in the presence of redshift uncertainties

To achieve our goal of dealing with the unknown redshifts and contamination we will need two pieces of Bayesian technology in addition to Bayes theorem: marginalisation and the product rule, which we briefly review here. The Product Rule states that, for any sets of parameters $\boldsymbol{\theta},\boldsymbol{\varphi}$ and ${\bf C}$, we have:

[TABLE]

The second piece of technology is marginalisation over nuisance parameters, $\boldsymbol{\varphi}$:

[TABLE]

The combination of these symbol manipulation techniques will be useful to deal with latent parameters; such as the true (unknown) redshift, $z$, and type $\tau$, of the supernova.

As a step to the fully general case, let us consider a single, known Type Ia supernova with an uncertain redshift, such as that derived from the supernova lightcurve. Very small spectroscopic redshift uncertainties, $\delta z$, have usually been dealt with by converting them into an additional error in the distance modulus, $\delta\mu$, by assuming a model $\mu(z)$ for the conversion, and adding the result in quadrature with the usual $\mu$ error. As we will show in detail later this is not statistically correct (for one since we don’t know the correct distance modulus describing our universe) and fails badly for typical photometric redshift errors.

To proceed rigorously we instead start by using Bayes theorem for the posterior for $\boldsymbol{\theta}$ given some data, $D$, from a single supernova. We then expand the arguments of the posterior to include the true redshift ($z)$ and distance modulus ($\mu$) of the supernova as latent (i.e. nuisance) parameters that we marginalise over, since we don’t know their true values:

[TABLE]

We repeatedly apply the product rule to rewrite this multi-dimensional integral as:

[TABLE]

We now make the simplifying assumption that the distribution of the true redshift is independent of the cosmological parameters222A more accurate assumption would be to posit that it depends on $\boldsymbol{\theta}$ only via the volume of spacetime. This is expected to be a weak dependence for the currently allowed range for $\boldsymbol{\theta}$. and note that since $\mu$ is assumed to be a deterministic function of $z$ and $\boldsymbol{\theta}$ the distribution $P(\mu|z,\boldsymbol{\theta})$ is a delta function333This is only exactly true in the background FLRW model and is not true if one allows for effects such as gravitational lensing, but we will ignore such perturbative effects here as is typical in supernova studies., which allows us to eliminate the $\mu$-integral444Using the standard identity $\int f(x)\delta(x-x_{0})dx=f(x_{0})$:

[TABLE]

This is the expression we were seeking. It expresses the posterior as a marginalisation over the unknown supernova redshift. What is the data $D$? In this case let us assume that we have extracted both an estimate of the redshift and distance modulus, ($z_{\rm obs},\mu_{\rm obs}$), for the supernova from its lightcurve. We now have:

[TABLE]

Assuming that the estimates $z_{\rm obs}$ and $\mu_{\rm obs}$ are independent, we can simplify this to:

[TABLE]

An important novel role in this analysis is played by the redshift prior $P(z)$. This has no parallel in the usual supernova analysis where the redshift is known spectroscopically. In the case of the uncertain supernova redshift we must consider Eddington bias: a supernova discovered with a 4m telescope with a redshift estimate of $z_{\rm obs}=0.75$, is much more likely to lie at a true redshift of $z=0.6$ than at $z=0.9$.

This insight, encoded in the prior $P(z)$, must be included in our inference analysis. Gull [26] was perhaps the first to highlight the importance of the data prior, showing that the use of the wrong prior (i.e. not the true distribution from which the data was drawn) leads to a bias.

Including this is crucial in any practical application since we are not sure about the true redshift prior, $P(z)$, which will depend both on the cosmological volume (which depends on the cosmological parameters) and the rates for each type of supernova. We can handle this uncertainty by introducing hyperparameters, $\varphi$, into the prior, $P(z,\varphi)$ which can then be fit for or marginalised over. Indeed, one can turn this around and instead marginalise over the cosmological parameters and supernova redshifts to produce a posterior for the rates of the different supernova types as functions of redshift555Numerical marginalisation is achieved simply by selecting and histogramming the parameters of interest from the MCMC chain.. Hence this “annoyance” which cannot be handled by the standard formalism becomes a feature in our analysis.

While this will be important for realistic simulations and analysis of real data it is a conceptually straightforward Bayesian extension and we do not consider it further here.

We now consider two special limiting cases to gain some intuition: (i) Spectroscopic galaxy redshifts but unknown host identity, and (ii) photometric supernova redshifts alone.

3.1 Spectroscopic galaxy redshifts but unknown host galaxy identity

In this subsection we answer the question: what happens if we have spectroscopic redshifts for all galaxies but are not sure which galaxy is the true host?

Since we take our data to be coming from the supernova lightcurve only, we can treat spectroscopic galaxy redshifts as a prior $P(z)$ on the SN redshift. Since we don’t know the true host in general we write this prior as a marginalisation over the potential host galaxies, which we index by $\gamma$:

[TABLE]

In the case of figure 1 this sum would likely just consist of just three terms: the two nearby galaxies and a third invisible galaxy which might be too faint to be detected in the image.

In eq. (3.11) $P(\gamma)$ is the probability that the supernova is hosted by the $\gamma^{\rm th}$ galaxy666Estimated for example using the projected distance from the centre of each galaxy of some other measure., and $P(z|\gamma)$ is the corresponding redshift distribution for the $\gamma^{\rm th}$ galaxy. Then we have:

[TABLE]

In the limit where the galaxy redshifts are known spectroscopically to high precision we can approximate the host redshift distribution as a delta function: $P(z|\gamma)=\delta(z-z_{\gamma})$ where $z_{\gamma}$ is the spectroscopically determined redshift of the $\gamma^{\rm th}$ galaxy. Then we can perform the marginalisation analytically, yielding:

[TABLE]

In other words, the final posterior is a mixture model that simply sums all the posteriors arising from assuming the supernova belongs to each of the potential host galaxies, weighted by the probability $P(\gamma)$ that the supernova belongs to each host. This is an intuitively pleasing and simple result.

3.2 Photometric redshifts

A second limiting subcase that will give some useful intuition is the case where we have no spectroscopic host information but instead only have a photometric redshift estimate from the supernova itself or from the host galaxy in the case where the host is unambiguous, i.e. we have an estimate for $P\left(z_{\rm obs}|z\right)$.

For simplicity assume that the resulting photometric redshift distribution is Gaussian777The generalisation to an arbitrary distribution is in principle simple since we perform the marginalisation numerically. The formalism remains unchanged with the new photometric redshift distribution.. If we assume $z_{\rm obs}$ and $\mu_{\rm obs}$ are correlated, for example if they both come from the lightcurve, then eq. (3.9) becomes:

[TABLE]

where $\Delta=\begin{pmatrix}\mu_{obs}-\mu\\ z_{obs}-z\end{pmatrix}$ and $C=\begin{pmatrix}\sigma_{\mu}^{2}&\sigma_{\mu z}\\ \sigma_{\mu z}&\sigma_{z}^{2}\end{pmatrix}$ is the covariance matrix.

Assuming independent $z_{\rm obs}$ and $\mu_{\rm obs}$ estimates for simplicity, eq. (3.9), reduces to:

[TABLE]

In our simulations we numerically perform the marginalisation in eq. (3.15) without considering correlations between $z_{\rm obs}$ and $\mu_{\rm obs}$, which is not generally true. However such correlations can be included in a straight forward way via modeling the covariance function, potentially with hyperparameters which can also be marginalised over.

Examining eq. (3.15) one might be tempted to Taylor expand $\mu(z,\boldsymbol{\theta})$ around $z_{\rm obs}$. Taking only the linear term in $\Delta z\equiv(z-z_{\rm obs})$ one can, assuming a Gaussian prior $P(z)$ with mean $\bar{z}$ and standard deviation $\sigma_{p}$, perform the redshift marginalisation analytically 888Note that one cannot simply assume an improper uniform prior since this leads to biases with real data while assuming a proper top-hat prior is no longer analytically integrable., giving:

[TABLE]

where $\tilde{\mu}_{\rm obs}=\mu_{\rm obs}-\mu(z_{\rm obs})+\mu^{\prime}z_{\rm obs}$. Completing the square and performing the integral gives:

[TABLE]

where $C$ is equal to :

[TABLE]

This has a simple interpretation: the $z$-error is converted into a $\mu$ error by using $\mu^{\prime}$, the instantaneous slope of the $\mu(z,\boldsymbol{\theta})$ curve, which is then added in quadrature with the existing $\mu$ error, $\sigma_{\mu}^{2}$; all modulated by the redshift prior width $\sigma_{p}$. Since we are assuming a Gaussian likelihood this is the same result as for the corresponding Fisher matrix analysis, see [27]. Note that to do this self-consistently one must change the $\mu$ error for each supernova at every point in the MCMC chain (since changing $\boldsymbol{\theta}$ changes $\mu^{\prime}$).

On the surface this offers a big computational simplification. For a large dataset of $N$ supernovae and $m$ cosmological parameters, eq. (3.17) only requires one to fit for the $m$ parameters rather than all $N+m$ parameters. Unfortunately in most cases eq. (3.17) is not sufficient for several reasons:

• •

It is only applicable when $P(z_{\rm obs}|z)$ is sufficiently narrow, i.e. $\sigma_{\mu}^{2}\gg\sigma_{z}^{2}\mu^{\prime 2}$. But what does this mean in practise?

• •

As alluded to previously we can only perform the integral in eq. (3.15) analytically for special priors such as an improper uniform ($-\infty<z<\infty$) or Gaussian prior.999This is integrable for a prior $P(z)\propto ze^{-\beta z}$ but this prior doesn’t have enough freedom in general for representing supernova rates in astronomy. These are not appropriate for astronomical surveys and hence assuming one of them will give biased results in general [26].

• •

In the general case we do not know the type of the supernova for sure. This implies that the photometric redshift for the supernova may be very non-Gaussian in general. This in turn implies that we cannot analytically integrate eq. (3.15) even if we assume a Gaussian redshift prior.

In the simulations below we find that this approach is biased. For all these reasons we strongly recommend numerically marginalising over all redshifts via MCMC or suitable sampling technique as we present in subsection 3.2.

4 Contamination from non-Ia supernovae

In the previous section we derived the posterior in the presence of redshift contamination but assumed that we knew the object was a SNIa. Unfortunately contamination of supernova types is also inevitable for photometric surveys. This has been addressed by the BEAMS (Bayesian Estimation Applied to Multiple Species) formalism for the case of spectroscopic redshifts in various papers [28, 29, 6, 17]. Here we highlight the key results using the hierarchical Bayesian approach we have used so far in this paper for redshift contamination as a warm-up to the case with both types of contamination.

In this case we can assume we know the true redshift of the supernova but are not sure of its type, labeled by a discrete variable $\tau$, which we here allow to take two values: $\tau=$ Ia and $\tau=\cancel{\rm Ia}$, the latter corresponding to non-Ia objects.101010If there were more than one class of object that gave useful information about $\boldsymbol{\theta}$ we would need to subdivide the non-Ias into more classes; e.g. if the $\boldsymbol{\theta}$ represented information about star formation rates rather than cosmology.

Previously we marginalised over the unknown latent variables $z$ and $\mu$. Now instead we marginalise over $\tau$, $\mu$ and $z$. Following the same approach taken in eq. (3.6) and thereafter, i.e., using Bayes theorem and the Product Rule to marginalise over the latent variables $\tau$, $\mu$ and $z$, one has

[TABLE]

Since we assume we know the redshift of the supernova, $z_{*}$, perfectly, we can write $P(z|\tau,\boldsymbol{\theta})$ as a delta function $\delta(z-z_{*})$ which allows us to do the $z$-integral, yielding:

[TABLE]

As before, $\mu$ is a deterministic function of $z$, $\tau$ and $\boldsymbol{\theta}$ and hence $P(\mu|z_{*},\tau,\boldsymbol{\theta})$ is also a delta-function, which collapses the $\mu$-integral. We write $\mu_{\rm{Ia}}$ for the expression $\mu(z,\boldsymbol{\theta},\tau=\rm{Ia})$ and $\mu_{\cancel{\rm Ia}}$ for $\mu(z,\boldsymbol{\theta},\tau=\cancel{\rm Ia})$. The $\tau$ ‘integral’ is actually just a sum over the two supernova classes. Here, $P(\tau)$ is the supernova type probability.111111This probability is estimated directly from the lightcurve using templates, inference or machine learning [13, 30]. For notational simplicity we write the probability of being a type Ia supernova as $P_{\rm{Ia}}$ (as for $P(\tau=\rm{Ia})$) and the probability of not being a type Ia supernova as $P_{\cancel{\rm Ia}}=1-P_{\rm{Ia}}$.

Then eq. (4.4) reduces to the usual BEAMS result:

[TABLE]

For a Gaussian likelihood with observed distance modulus $\mu_{\rm obs}$, we can write:

[TABLE]

What about the non-Ia term $\mu_{\cancel{\rm Ia}}$? Following [17] and [28] one can take $\mu_{\cancel{\rm Ia}}=\mu_{\rm{Ia}}+b(z,\varphi)$ where $b(z,\varphi)$ allows both for the redshift evolution of the non-Ia population and the difference in mean intrinsic luminosity compared with SNIa, parameterised by some unknown hyperparameters $\varphi$ which are fit simultaneously with all the other parameters. This allows one to learn about the non-Ia populations. For cosmology, the non-Ia population is essentially uninformative and hence one can simply take the corresponding $\sigma_{\cancel{\rm Ia}}$ to be large in actual analysis.

5 The general case: both type and redshift contamination

We are now ready to consider the general case with both type and redshift contamination. Hence we assume that the type $\tau$ of the supernova is unknown, and that we either have photometric redshift estimate for the supernova or redshift information (either photometric or spectroscopic) of potential hosts galaxies.

As before we will treat the host galaxy information as a prior, $P(z)$, on the SN redshift (see eq. (3.11))121212The supernova rates are now a much less important factor in the prior since the redshift priors arising from the galaxies should be very peaked in comparison.:

[TABLE]

These need not be spectroscopic redshifts and $P(z|\gamma)$ may have significant spread.

As before we begin from Bayes rule, introduce and marginalise over latent variables for the unknown true variables, $z,\mu,\tau$, and then repeatedly apply the product rule:

[TABLE]

where the last step accounts for the fact that $\tau$ is a discrete (categorical) variable. We can now substitute for the prior $P(z)$ expressed as the sum over the potential host galaxies, eq. (3.11), which gives:

[TABLE]

Eq. (5.8) is our main result for the posterior arising from a single supernova. It handles both Type and redshift-host contamination. Note that because the true redshift of the supernova is unknown we can allow for redshift-dependence of the type probabilities, $P(\tau|z)$, whereas in previous BEAMS analysis this was taken as a constant.

To compute the full posterior over $N$ supernovae with collective data $\{D_{i}\}$, we now have $2N$ nuisance parameters corresponding to the type and redshift, $(\tau,z)$ for each supernova, in addition to the cosmological parameters $\boldsymbol{\theta}$. We can then either assume the supernovae are independent, and simply multiply the posteriors to get:

[TABLE]

or allow for correlations between the supernovae as in [29].

In practise the marginalisation over redshift required by eq. (5.9) can be achieved efficiently through MCMC by allowing the redshift of each supernova to be a free nuisance parameter that is varied along with the cosmological parameters $\boldsymbol{\theta}$. If the supernovae are correlated then we must also introduce the type of each supernova to include as a parameter. This does not significantly alter the MCMC analysis [29].

At this point one can again ask what the data, $D_{i}$, is for each supernova? At the most basic – and correct – level this would be the lightcurve measurements in various colour bands as a function of time. However, a convenient simplification is to consider $D_{i}=(z_{\rm obs},\mu_{\rm obs})$, assuming that the object is a SNIa. In general this would be inappropriate if more than one type of supernova contained useful information about the cosmological parameters $\boldsymbol{\theta}$. However, since in the non-Ia case the derived redshift and distance modulus give almost no useful cosmological information, one can simply take $P(z_{\rm obs}|z,\tau=\cancel{\rm Ia})$ to be a wide uniform distribution or very wide Gaussian and use $z_{\rm obs},\mu_{\rm obs}$ derived assuming the object is a SNIa [28, 6]. Then the fact that one is using the “wrong” values for $z_{\rm obs},\mu_{\rm obs}$ has no impact. We note that obviously this is only an issue when the type of the supernova is unknown.

5.1 Dealing with Selection Bias

Our eq. (5.9) is not completely general of course. It neglects the fact that real surveys will have censored data due to selection effects, also known as Malmquist bias: the faint end of the population will not be detected, especially at high redshift. If untreated this will bias the cosmological results. One approach is to apply a Malmquist bias correction directly to the data before undertaking the analysis, as done in the analysis of the SDSS-II data in [17]. This is unsatisfactory from a Bayesian point of view however. Another approach, following e.g. [6], is to consider the full population of supernovae, only a subset of which are actually detected by any given SN survey. Introducing a latent label $\epsilon$ which is unity if the SN is detected by the survey and zero otherwise, we can then treat $\epsilon$ as we have the other latent variables and marginalise over it, which will lead to the appearance of the probability distribution $P(\mbox{detect}|D)$ which will depend explicitly on the telescope and survey strategy. In addition, we will need to marginalise over the unknown number of undetected supernovae.

Within our current approach of treating the data to be the distance modulus and redshift, we could simply truncate the likelihood at a given threshold value of $\mu$. This has been implemented in the context of a hierarchical Bayesian model like the one used here by Sereno [31]. However, a simple censoring of distance moduli will not capture the true complexities of selection effects in real surveys. A better approach will be to go back to the lightcurve flux measurements directly but this in turn will depend on the specific lightcurve fitter one uses; see [6] for an analysis of Malmquist bias in the context of the SALT II lightcurve fitter and e.g. [32] section 3.1 for a pedagogical introduction to dealing with selection effects within Bayesian statistics. In such studies a key role is played by demanding that the censored likelihood is still normalised.

We will ignore selection bias in our simulations since to do so correctly would require working at the level of light curves, which we have avoided for simplicity, and the standard approach still fails dramatically even in this idealised case. We now move to test the zBEAMS formalism in different cases.

6 Simulations

Here we perform an illustrative set of catalogue simulations to show how zBEAMS, as described in sections 3 to 5, recovers the correct cosmological parameters by marginalising over unknown supernova types and redshifts. We consider the two cases corresponding to sections 5 and 3.2: one where we have spectroscopic redshift estimates from potential host galaxies, hereafter referred to as the spectroscopic case, and one where we have photometric measurements only, hereafter referred to as the photometric case. Conceptually, it is trivial to combine these cases when dealing with a dataset with mixed spectroscopic and photometric measurements, however we keep the two cases distinct.

In both cases for simplicity we assume that all objects are detected, no matter how faint. We demonstrate that while the standard likelihood, ignoring redshift uncertainty, results in biased cosmological parameter estimates even in this case, zBEAMS correctly marginalises over type and redshift uncertainties, recovering the fiducial cosmology.

For all simulations, we assume a flat $\Lambda$CDM universe with a fiducial cosmology given by the latest results from the Planck collaboration [33], i.e., $H_{0}=67.74~{}\rm{km/s/Mpc}$, $\Omega_{m}=0.31$ and $w=-1$. We perform inference over the parameters using Markov Chain Monte Carlo (MCMC) methods, specifically we use the Metropolis-Hastings [23, 24] algorithm for low-dimensional sampling and block Metropolis-Hastings when numerically marginalising over redshift in the photometric case. Detailed descriptions of both cases follow below.

6.1 Spectroscopic Case

For the spectroscopic case, we simulate 1000 SNe from a uniform distribution across a redshift range of $z=0.015-1$. We create two datasets, one pure Ia dataset without any redshift errors for reference, herein referred to as the unbiased dataset, and one with both host galaxy redshift and non-Ia contamination, herein referred to as the biased dataset.

The unbiased SNe data is generated with a dispersion of $0.2~{}\rm{mag}$. For the biased dataset, we assume the host galaxies have spectroscopically confirmed redshifts, but the supernova is observed using photometry, hence the supernova type is not known and it is not always clear which galaxy the supernova belongs to if multiple galaxies lie within a small angular distance of one another. We also considered that supernovae with $z<0.1$ will be identified spectroscopically. We assume $5\%$ type misidentification where the non-Ia population is offset from the Ia population by $2~{}\rm{mag}$, and has a Gaussian dispersion of $1.5~{}\rm{mag}$ (similar to [17]). A more realistic distribution can be used with the same method, provided the form is known. We assume a $9\%$ host misidentification [21], where the misidentified host redshift is drawn from a normal distribution $z\sim\mathcal{N}(z_{\rm{true}},0.1^{2})$.

In figure 2 we plot the distance modulus of a contaminated SNe Ia dataset. The biased dataset we analyse using both the standard likelihood (which does not take redshift error into account) and with the zBEAMS likelihood. We use the zBEAMS posterior in eq. (5.9) to fully marginalise over both type and redshift uncertainties and thus produce unbiased cosmological estimates. In this analysis, we solve for $\Omega_{m}$, $H_{0}$ and $w$ while we assume the parameters of the populations (such as the magnitude offset and standard deviation of the non-Ia population) are known. However, it would be simple to solve for these simultaneously as done in earlier BEAMS papers [17, 29]. We infer the marginalised posterior distribution for $w$ and $\Omega_{m}$ for each of these three instances, and their respective contours are shown in figure 3.

6.2 Photometric Case

Here we consider the case where the redshift of the host galaxy is obtained photometrically. Now to use zBEAMS, the marginalisation over redshift must be performed numerically, as the integral in eq. (3.15) has no analytic solution. We assume for this work that the redshift uncertainties are Gaussian distributed with a standard deviation of $0.04(1+z)$, though any more realistic distribution can be assumed with little change in complexity. We simulate 998 SNe from a redshift distribution given by $P(z)\sim ze^{-\beta z}$ spanning the redshift range $z=0.015-1.4$, with $\beta=3$. Note that this distribution is also the prior on the fitted redshifts and in reality one would need to fit for $\beta$ simultaneously with all other parameters. While we do not address this here, this distribution could be extended to include modeling of instrumental selection effects, in addition to intrinsic supernova rate information. As before, we assumed Gaussian errors with dispersion $0.2~{}\rm{mag}$ in the distance modulus. In figure 4 one can see the magnitude residuals for the observational redshift (main figure) and for the redshifts recovered by the zBEAMS analysis (inset figure).

We then do a MCMC using the standard likelihood (where we only solve for 3 cosmological parameters assuming incorrectly that all the redshifts are correct). The result is plotted in figure 6 in shades of red. Note that the result is clearly biased with respect to the solid black contours which are obtained using the true SNe redshifts. The blue contours in figure 6 show the result when applying zBEAMS to the biased dataset. We used a block Metropolis-Hastings sampling method – affectionately dubbed “Arabian nights” – to fit for 1001 parameters simultaneously (3 cosmological parameters and 998 redshifts), i.e., numerically computing the posterior given by eq. (3.15). The block Metropolis-Hastings proceeds identically to the usual Metropolis-Hastings sampling algorithm, except that parameters are updated in blocks instead of updating all parameters every step. We took the 3 cosmological parameters and each redshift as a block, but experimented with block sizes of 1-10 redshifts.

The block size has little impact on accuracy, as long as the blocks are small enough not to reduce the acceptance ratio, they do however impact on algorithm speed. We assume that each supernova redshift has a prior coming from the host galaxy (or from the supernova lightcurve itself) which we take to be Gaussian centred on the observed redshift with standard deviation of $0.04(1+z_{\rm obs})$. The prior on the overall SNIa redshift distribution was taken to be $P(z)=ze^{-\beta z}$, where we fixed the value of $\beta$ to 3. In a case with real data, one would need to fit for these hyperparameters as well.

We found that the block Metropolis-Hastings recovers the true redshifts of the low-z supernovae ($\sigma_{z}=0.02$ for $z<0.25$) well, with worsening performance as the redshift increases. This is due to two effects: first we assume the photometric redshift error scales with $(1+z)$ and secondly the Hubble diagram progressively flattens out at $z>0.25$ removing the signal that allows MCMC to constrain the redshift. This can be clearly seen in figure 5 where we show the stacked 1D histograms $z_{i,\rm chain}-z_{i,\rm true}$ for all 998 SNIa. It can be seen that while the error increases with redshift, the redshift estimates show no systematic bias. The marginalised posterior distributions for $w$ and $\Omega_{m}$ for each of these instances are represented in the contour plots shown in figure 6. We can see that zBEAMS recovers the correct cosmology, and contours, as desired.

Examining figure 4 we can see the origin of the bias of the standard likelihood. A large number of the data points are more than $3\sigma$ away from the model and some are over $20\sigma$ away. This is an artefact of using the wrong redshifts. The inset shows the same residual Hubble diagram when the data is instead plotted using the mean redshifts recovered from the MCMC chain for each redshift. Very few datapoints are now more than $2\sigma$ from the fiducial model, even at low redshifts where the excursions were the strongest.

This allows zBEAMS to produce unbiased cosmology contours that almost match the contour sizes of the perfect, spectroscopic case. A more standard approach to this same problem might be to significantly increase all the $\mu$-error bars of the points to account for the redshift uncertainties. We find that doing so yields biased results for reasons discussed in section 3.2. Increasing the error bars further might unbias the contours but only at the expense of significantly inflating the associated contours.

It should be noted that accurate sampling in realistic scenarios will not be trivial since we are fitting for a posterior that typically has more unknown parameters than data points. We have explored both block Metropolis-Hastings and other algorithms such as Diffusive Nested Sampling [34] as viable solutions. Hamiltonian Monte Carlo [35] may also be well-suited to the high-dimensionality of this problem.

The code used in this section is available on Github at:

https://github.com/MichelleLochner/zBEAMS.

7 Conclusions and Future Work

Future large surveys such as LSST will likely deliver large numbers ($>10^{5}$) of good supernova candidates without the spectroscopic confirmation that has historically been required to use them for cosmology. This means we will have to do a cosmological analysis with a sample for which both the true type and redshift of the supernovae are unknown. Instead, only the probability distributions of both supernova type and redshift will be available. In particular, the redshifts will be uncertain either because the redshift is only known photometrically or because the true identity of the galaxy hosting the supernova is unsure, even if the redshifts of potential host galaxies are perfectly known.

In this paper we have shown how to achieve unbiased cosmology with such a sample, simultaneously handling both the non-Ia contamination and the problem of imprecise supernova redshifts in a unified framework. Our formalism - zBEAMS - generalises the original BEAMS formalism [28] to handle the redshift uncertainties of the supernovae by employing a hierarchical Bayesian approach. We introduce nuisance parameters for the type and redshift of each supernova in the dataset and then marginalise over these nuisance parameters using numerical sampling. In the special case where the supernova belongs to one of several galaxies each with spectroscopically known redshift, the final posterior is a simple weighted mixture model over the posteriors assuming the supernova is in each of the different potential host galaxies (see eq. (3.13)).

We show in figure 3 that a model with a $9\%$ host misidentification error leads to large bias ($\sim 3\sigma$) using the standard inference approach while zBEAMS removes the bias at essentially no extra computational cost. We also consider the case of photometric uncertainties, where we numerically marginalise over 998 redshift parameters to produce the unbiased contours in figure 6.

There are a number of ways in which this work can be extended in a straight-forward way:

• •

In this analysis we have assumed that the probability of belonging to a given host galaxy $\gamma$, encoded in the terms $P(z|\gamma)$, and the probability of being a given type of supernova, $\tau$, encoded in $P(\tau|z)$, are known apriori. It would be interesting to extend our formalism to allow these to be partially known nuisance parameters that are estimated by available data.

• •

The zBEAMS formalism could be extended to include correlations with host galaxy information, such as host influence on Hubble residuals via stellar mass etc. [36, 37].

• •

We found that the redshift distribution plays an important role for the photometric redshift case. A much more complex model could be used than the one we assumed which could include some systematic effects and allow one to learn something about supernova rates (See Malz & Peters, et al. (in prep.) for upcoming work on this problem in the LSST context).

• •

While we have presented the zBEAMS formalism emphasising its generic nature for any data $D$, in our examples we took $D$ to be the measured distance moduli. It would be useful to develop zBEAMS specifically for the case in which $D$ is the set of lightcurve flux measurements in different bands; i.e. one step further back in the analysis chain.

• •

As discussed in detail in (5.1), realistic supernova surveys censor the true SN population because of the magnitude limits of the telescope and cuts performed during the analysis. Selection effects within a Bayesian framework have already been extensively covered in e.g. [6] and could be incorporated into the zBEAMS likelihood.

These extensions are left to future work.

Acknowledgments

The authors thank Alan Heavens, Nosiphiwo Zwane, Laura Richter and other participants of the 2016 Bayes School in Stellenbosch for early contributions and Rahul Biswas, Rick Kessler, Bob Nichol and Hiranya Peiris for comments on the draft. We also thank the anonymous referee for comments. BB thanks Renée Hlozek and Martin Kunz for long-term collaborations related to BEAMS and Roberto Trotta for discussions. ML thanks Brendon Brewer for help with DNest. JF thanks Mário G. Santos for useful discussions and Marta Spinelli for help with MCMC techniques. This paper originated at the 2016 Bayes School and Workshop in November 2016 funded by NITHeP, and grew out of Alan Heaven’s lectures at the same School. The authors acknowledge support from the NRF, the South African Square Kilometre Array Project, the National Astrophysics and Space Science Programme and AIMS. We acknowledge use of the Hubble ESA Archive for images. This paper and the research it describes were created using scrum methodology.