The Type Ia Supernova Color-Magnitude Relation and Host Galaxy Dust: A Simple Hierarchical Bayesian Model

Conventional Type Ia supernova (SN Ia) cosmology analyses currently use a simplistic linear regression of magnitude versus color and light curve shape, which does not model intrinsic SN Ia variations and host galaxy dust as physically distinct effects, resulting in low color-magnitude slopes. We construct a probabilistic generative model for the dusty distribution of extinguished absolute magnitudes and apparent colors as the convolution of a intrinsic SN Ia color-magnitude distribution and a host galaxy dust reddening-extinction distribution. If the intrinsic color-magnitude ($M_B$ vs. $B-V$) slope $\beta_{int}$ differs from the host galaxy dust law $R_B$, this convolution results in a specific curve of mean extinguished absolute magnitude vs. apparent color. The derivative of this curve smoothly transitions from $\beta_{int}$ in the blue tail to $R_B$ in the red tail of the apparent color distribution. The conventional linear fit approximates this effective curve near the average apparent color, resulting in an apparent slope $\beta_{app}$ between $\beta_{int}$ and $R_B$. We incorporate these effects into a hierarchical Bayesian statistical model for SN Ia light curve measurements, and analyze a dataset of SALT2 optical light curve fits of 248 nearby SN Ia at z<0.10. The conventional linear fit obtains $\beta_{app} \approx 3$. Our model finds a $\beta_{int} = 2.3 \pm 0.3$ and a distinct dust law of $R_B = 3.8 \pm 0.3$, consistent with the average for Milky Way dust, while correcting a systematic distance bias of $\sim 0.10$ mag in the tails of the apparent color distribution. Finally, we extend our model to examine the SN Ia luminosity-host mass dependence in terms of intrinsic and dust components.


INTRODUCTION
Type Ia supernova (SN Ia) rest-frame optical light curves have been used as cosmological distance indicators to trace the history of cosmic expansion, detect cosmic acceleration (Riess et al. 1998;Perlmutter et al. 1999), and to constrain the equationof-state parameter w of dark energy (Garnavich et al. 1998;Wood-Vasey et al. 2007;Astier et al. 2006;Kowalski et al. 2008;Hicken et al. 2009b;Kessler et al. 2009a;Freedman et al. 2009;Amanullah et al. 2010;Conley et al. 2011;Sullivan et al. 2011;Rest et al. 2014;Scolnic et al. 2014a;Betoule et al. 2014). Determining supernova distances with high precision and small systematic error is essential to accurate constraints on the cosmic expansion history and the properties of dark energy.
Inferring peak optical absolute magnitudes of SN Ia from distance-independent measures such as their light curve shapes and colors underpins the evidence for cosmic acceleration. Empirical studies show that SN Ia with broader, slower declining optical light curves are more luminous ("broader-brighter") and that SN Ia with redder colors are dimmer. But the "redder-dimmer" colorluminosity relation widely used in cosmological SN Ia analyses masks the fact that it has two separate physical origins. An intrinsic correlation arises from the physics of exploding white dwarfs while interstellar dust in the host galaxy also makes supernovae appear dimmer and redder. The conventional approach of fitting a single (usually linear) function for luminosity vs. color, for a given light curve shape, is too simple. This leads to considerable uncertainty regarding the physical interpretation of the color-luminosity distribution of SN Ia, the confounding of extrinsic host galaxy dust reddening with the intrinsic color variations of SN Ia, and the proper way to use SN Ia color measurements to estimate accurate photometric distances. In this paper, we present a new probabilistic model describing the apparent SN Ia colormagnitude distribution as arising from the combination of intrinsic color-luminosity variations and host galaxy dust reddening and extinction, and apply this to SN Ia data to determine the characteristics of these physical components.
Cosmological analyses of high-z SN Ia data depend on empirical correlations originally observed in samples of nearby low-z SN Ia (Hamuy et al. 1996a;Riess et al. 1999;Jha et al. 2006;Hicken et al. 2009a;Contreras et al. 2010;Stritzinger et al. 2011;Hicken et al. 2012). Light curve fitting methods, including MLCS (Riess et al. 1996a(Riess et al. , 1998Jha et al. 2007), SALT2 (Guy et al. , 2010, SNooPy (Burns et al. 2011), and BayeSN , all make use of the optical luminosity-light curve width correlation (Phillips 1993;Hamuy et al. 1996b;Phillips et al. 1999). However, current approaches conceptually differ on how measured apparent colors are used to infer the SN Ia luminosities and thus estimate the photometric distance. Methods such as MLCS, SNooPy, and BayeSN explicitly model the intrinsic SN Ia light curves and the effects of host galaxy dust extinction as separate components. However, the most popular tool for fitting cosmological SN Ia light curves is currently SALT2, a spectral template model that does not attempt to separate intrinsic SN Ia variations from host galaxy dust effects.
A longstanding puzzle in the analysis of SN Ia light curves is the nature of their apparent color and brightness variations. In principle, they comprise color and luminosity variations intrinsic to the supernovae, as well as reddening and extinction by interstellar dust along the line of sight in their host galaxies. However, the fact that astronomers only observe the combination of these effects poses a challenging inference problem. The function of dust absorption over wavelength (e.g. CCM, Cardelli et al. 1989) is typically parameterized by the ratio of total to selective extinction, This ratio normally has an average value of 3.1 for interstellar dust in the Milky Way (MW) Galaxy, although it can vary between 2.1 and 5.8 (Draine 2003). Schlafly et al. (2016) find that the MW extinction curve is fairly uniform: with a narrow spread σ(R V ) ≈ 0.18. Similar extinction curves have been found in external galaxies; for example, Finkelman et al. (2008Finkelman et al. ( , 2010 found average values of R V ≈ 2.8. Early analyses of SN Ia data estimated unphysically low values of R V 1 (c.f. Branch & Tammann 1992, for a review), although these analyses did not take into account empirical correlations between the luminosity, color and light curve shape of the events. Riess, Press, & Kirshner (1996b) used the first MLCS (Riess, Press, & Kirshner 1996a) method to fit SN Ia optical BV RI light curves and, by minimizing the Hubble diagram scatter, they found a dust R V = 2.6 ± 0.3, consistent with the Milky Way average. They noted that the failure to properly account for intrinsic color-luminosity correlations would cause estimates of the dust extinctionreddening ratio R V to be biased low.
Extending the wavelength range of the observations from the optical to the rest-frame NIR improves the constraints on the dust law and thus helps disentangle intrinsic color variation from dust reddening. For several nearby, highly reddened (peak apparent B − V ≫ 0.6) SN Ia with optical and NIR light curves, R V can be fit precisely and unusually low values of 1.5-1.8 have been reported (Krisciunas et al. 2007;Elias-Rosa et al. 2006Wang et al. 2008). While the origin of these apparently low R V values is poorly understood (Wang 2005;Goobar 2008;Phillips et al. 2013;Johansson et al. 2014;Amanullah et al. 2015), these extremely red and dim objects are not present in the cosmological sample due to selection effects and cuts (B − V < 0.3). Freedman et al. (2009) constructed a SN Ia Hubble diagram using rest-frame i-band magnitudes and B − V colors. By minimizing the Hubble residuals χ 2 , they estimated R V ≈ 1.74 ± 0.27, and suggested that either dust in SN Ia host galaxies has a substantially different extinction law than Milky Way dust, or that there is significant SN Ia intrinsic color-luminosity dispersion, independent of the light curve shape. The latter hypothesis was supported by Folatelli et al. (2010) to explain a discrepancy found when analyzing of the nearby Carnegie Supernova Project (CSP, Contreras et al. 2010) SN Ia sample. When examining the optical-NIR colors they inferred R V = 3.2 ± 0.4 when the extremely red objects (E(B − V ) 1) were excluded. However, when minimizing the Hubble diagram dispersion, they find low values of R V ≈ 1 − 2. Recent analyses of larger optical-NIR nearby samples found that the majority of SN Ia Ia with low reddening appear extinguished by dust with R V closer to 3 Phillips 2012;Burns et al. 2014).
Another strategy for separating intrinsic color variation from dust effects is to measure spectral features that correlate with the intrinsic photometric properties of SN Ia. For example, Foley & Kasen (2011) have correlated the velocity of the Si II λ6355 absorption feature with the peak intrinsic B − V color, and found R V ≈ 2.5. As they controlled for optical decline rate, this is evidence for independent intrinsic color variations. Chotard et al. (2011) modeled the components of the apparent SN Ia spectroscopic and photometric variations depending upon Si II and Ca II H&K equivalent widths, finding that the remainder is well described by a CCM dust reddening law with R V = 2.8 ± 0.3, consistent with the MW value. However, the vast majority of the high-z SN Ia observations currently used for cosmological analysis consist of rest-frame optical photometry and lack the high-quality spectra or rest-frame NIR data used by the above analyses.
In the following subsections, we describe the conventional method, the Tripp formula, for modeling correlations between SN Ia magnitude, color and light curve shapes from optical light curve data, currently used to estimate cosmological distances. SALT2mu is a generalization of this method to account for scatter around the Tripp formula (Marriner et al. 2011). We comment on the drawbacks of these approaches, primarily because they are inadequate for properly accounting for the physically distinct factors of intrinsic SN Ia variation and host galaxy dust underlying the data. We introduce a new statistical model, Simple-BayeSN, that we have developed to address these shortcomings. Simple-BayeSN analyzes the peak apparent magnitude, apparent color, and light curve shape obtained from light curve fits to the SN Ia photometric time series. It models the SN Ia data as arising from a probabilistic generative process combining intrinsic SN Ia variations, host galaxy dust effects, and measurement error. Simple-BayeSN uses a hierarchical Bayesian framework to fit the SN Ia data on the Hubble Diagram, while coherently estimating the parameters driving the underlying effects.

The Tripp Formula
Conventional cosmological SN Ia analysis (e.g. Betoule et al. 2014;Rest et al. 2014) currently proceeds by fitting primarily rest-frame optical light curve data to obtain estimates of peak apparent magnitude m s , apparent color c s and light curve shape x s for each SN s. A simple linear regression model for the absolute magnitude M s as a function of the distance-independent light curve observables is constructed using the Tripp formula 7 : (Tripp 1998), where µ s is the supernova distance modulus. The global coefficients (M 0 , α, β) are found by fitting this relation with measurements (m s ,ĉ s ,x s ) on the Hubble diagram, µ s = µ(z s ), for a sample of SN Ia {s}. The optical "broader-brighter" width-luminosity relation (Phillips 1993) is captured by α, the "redderdimmer" color-luminosity relation by β. The expected absolute magnitude at x s = c s = 0 is M 0 . This equation would be correct if the color-luminosity relation were entirely due to (small amounts) of dust (and the light curve shape dependence truly linear). In that case m s − βc s is the reddening-free Wesenheit magnitude (Madore 1982). However, in the presence of intrinsic color-luminosity variations, this formula is not fundamental: it is just the simplest linear model for absolute magnitude as a function of the observables, light curve shape and apparent color.
The residual scatter 8 around this model that is unaccounted for by measurement error or peculiar velocities is ǫ s res with variance σ 2 res . This contributes to the uncertainty in absolute magnitude M s , and thus the photometric distance modulus µ s , of an individual SN.
In typical usage, M s and m s are peak absolute and apparent magnitudes effectively in rest-frame B-band, and the color corresponds to peak apparent B−V . Hence β is the slope of the change in B-magnitude for a unit change in B − V . This can be compared to the expected slope for normal Milky Way dust extinction A B vs. E(B − V ) reddening, R B ≡ R V + 1 = 4.1. Tripp (1998) and Tripp & Branch (1999) originally found β ≈ 2 using peak apparent B − V colors and ∆m 15 (B) for light curve shape (Phillips 1993). Using the first SALT model to determine optical light curve stretch and color, Guy et al. (2005) and Astier et al. (2006) also found low values β ≈ 1.5 − 2. Conley et al. (2007) fit nearby SN Ia on the Hubble diagram and found that the empirical relation between SN Ia optical luminosity and apparent color, controlling for light curve shape, still required a low value of β ≈ 2. They speculated that this is much less than the normal dust R B ≈ 4 either because the dust in SN Ia hosts is nonstandard, or because the estimated β may actually be measuring some combination of intrinsic color variations (not accounted for by light curve shape x s ) and normal interstellar dust.
The SALT2 spectral template (Guy et al. , 2010 is the most popular (and well-tested, Mosher et al. 2014) model currently applied to fit cosmological SN Ia light 7 This formula is often written with an arbitrary negative sign preceding α (e.g. Guy et al. 2005;Astier et al. 2006). When it is regarded as a linear regression model for absolute magnitude versus the covariates x and c, it is most natural to precede all the regression coefficients by a positive sign. 8 This variance term has been variously called the intrinsic scatter or dispersion σ int (e.g. Astier et al. 2006;Conley et al. 2011;Marriner et al. 2011). However, in their conventional usage, there is no implication that this scatter is solely attributed to physical properties intrinsic to the SN Ia with host galaxy dust subtracted. Scolnic et al. (2014b) instead refers to it as residual scatter: it is the additional variance needed to account for the scatter in the Hubble residuals. We adopt the latter usage to avoid confusion, and reserve intrinsic within our model to conceptually refer to the latent properties of the SN Ia in the absence of host galaxy dust. curve data. Using SALT2 and the Tripp formula to fit SNLS1 and SDSS-II SN Ia data, Guy et al. (2007) and Kessler et al. (2009a) found β ≈ 1.8 − 2.6. For SNLS3, Guy et al. (2010) and Conley et al. (2011) found β ≈ 3.1. Similarly, the combined SDSS+SNLS3 [JLA] analysis of Betoule et al. (2014) obtained β = 3.10 ± 0.08, significantly less than R B = 4.1. Marriner et al. (2011) introduced a more general formalism (SALT2mu) for accounting for the residual scatter around the linear model, Eq.1, and fitting for its regression coefficients. Sensible methods for fitting Eq. 1 take into account the fact that the fitted values (m s ,ĉ s ,x s ) are different from the true, latent (unobserved) values (m s , c s , x s ) that obey Eq. 1. The difference amounts to random measurement error with the SALT2 fit covariance matrix. Marriner et al. (2011) further supposes an additional source of residual scatter between the measured values and the latent values. For example, the measured color of a SN s is decomposed aŝ

Luminosity vs. Color Residual Scatter
where c s mod is the "model" color component that enters into the Tripp model (Eq. 1), and is linearly correlated with the luminosity, c s n is the color measurement error "noise" and c s r is a random "residual" color scatter term (Scolnic et al. 2014b). Similar equations can be written for the scatter in the magnitude and light curve shape components of the data. While the measurement errors are quantified by a covariance matrix estimated from the SALT2 light curve fit, the residual scatter covariance matrix Σ r is a priori unknown and poorly constrained by the data, and some choices must be made regarding its entries. To date, it has been used generally with non-zero entries only for residual magnitude σ 2 mr and color σ 2 cr variances. The special case in which only the residual magnitude variance entry is nonzero corresponds to the conventional assumption that the residual scatter is attributed to unexplained variance in luminosity, σ 2 res = σ 2 mr . With an assumed residual matrix, SALT2mu estimates the M 0 , α, β coefficients by minimizing the Hubble residual χ 2 , modified to incorporate measurement errors and the residual scatter matrix. Applying SALT2mu to SDSS-II data, Marriner et al. (2011) find that attributing the residual scatter only to color led to a larger estimated β ≈ 3.2, still significantly less than R B = 4.1. Scolnic et al. (2014b) analyzed a combined dataset, consisting of SDSS-II, SNLS3 and nearby samples, to examine the dependence of the estimated β on the relative attribution of residual scatter to luminosity or color when analyzing the data within the SALT2mu framework. When "luminosity variation" is assumed (σ mr ≫ σ cr ), the estimated β ≈ 3.2 but when "color variation" is assumed (σ mr ≪ σ cr ), it increases to β ≈ 3.7, closer to the normal MW average. Simulations from a color variation model, in which c mod had a "dust-like" distribution (with a "one-sided" tail to the red, but a sharp edge to the blue), with a MW dust-like true β = 4.1, could better match the pattern of Hubble residuals vs. color seen in the data, compared to a luminosity variation simulation with a conventional β = 3.1. They estimated that a misattribution of the residual scatter to luminosity rather than color variation, could lead to a bias ∆β ≈ −1 in the recovered slope, and a 4% shift in the inferred w.
Recently, Scolnic & Kessler (2016) presented a method to determine the underlying distribution of the latent c mod colors by matching the observed data distribution to realistic forward simulations that incorporate measurement noise, and selection effects, and two SALT2 spectral variation models: luminosity-variation dominated (G10, Guy et al. 2010) and color-variation dominated (C11, based on Chotard et al. 2011). Applying this to the cosmological SN Ia compilation of Scolnic et al. (2015), they uncover a "dust-like" underlying color distribution, with a red tail and blue edge, when simulating with a colordominated C11 model. By matching their simulations to data, they find β = 3.85 using the color-variation model, and β = 3.1 with the luminosity-variation model.
These analyses suggest that the proper modeling of color-luminosity subcomponents is important for understanding the observed SN Ia color-magnitude distribution, the accurate estimation of distances, and inferences for cosmology.

Shortcomings of these approaches
The simplicity of the Tripp formula has enabled its widespread application in cosmological SN Ia analyses. However, in its conventional usage, Eq. 1 is too simplistic. Because neither the magnitude m s nor the color c s that enter into it are corrected for host galaxy dust extinction or reddening, the absolute magnitude M s = m s − µ s is actually the dust-extinguished absolute magnitude M ext s , and c s is the dust-reddened apparent color c app s . However, in reality, the extinguished absolute magnitude results from the dimming of the supernova's intrinsic luminosity by dust extinction M ext The apparent color results from the dust reddening the supernova's intrinsic color: c app The dust extinction is surely correlated with its reddening, and can only be positive. The supernova's intrinsic luminosity may be correlated with its intrinsic color, independently of light curve shape, as speculated by e.g. Conley et al. (2007) and Freedman et al. (2009); at the very least we do not know that this intrinsic correlation is zero, and would like to estimate it. By regressing only the sum M ext s against the sum c app s , the Tripp formula tries to capture all the color-magnitude correlations in a single trend with slope β. As it is a priori highly unlikely that the intrinsic color-magnitude slope would be exactly equal to the dust reddening-extinction law, this parameterization is clearly limited and inadequate. It fails to distinguish between the different physical characteristics of color-luminosity variation intrinsic to the SN Ia vs. reddening-extinction by extrinsic host galaxy dust.
The residual matrix framework of SALT2mu adds some additional degrees of freedom to the conventional Tripp formula. By decomposing the apparent color c s into a model color c mod and residual color c r , and allowing σ cr > 0, one can in effect capture additional colormagnitude variations. However, this residual color scatter does not correlate with luminosity; this additional component essentially has its own β r = 0. Hence, the residual color scatter by itself would not be able to capture a separate non-zero color-luminosity correlation different from β. This suggests that the color-magnitude covariance entry in the residual matrix Σ r could be made non-zero. This raises two challenges: either a numerical value for this residual covariance would have to be set a priori, or it would have to be inferred jointly with the other parameters. Unfortunately, in the former case, it is unclear what value to set it to a priori; in the latter case, the inference of this covariance would likely be highly degenerate with β, unless strong priors were set.
Systematic uncertainties in the treatment of dust and color of SN Ia Ia have important implications for cosmological inference. If the single slope β is actually measuring a combination of intrinsic color-luminosity variation and host galaxy dust reddening-extinction, then different SN Ia subsamples may have different proportions of each. The proportions may even be redshift-dependent, owing to the physical environment of the progenitor systems, host galaxies, or selection effects. Applying one β slope across the entire sample would incur complex systematic biases that propagate into cosmological inferences. Resolving the confusion between the intrinsic variation and extrinsic host galaxy dust effects is imperative for the proper analysis of SN Ia observables.

Simple-BayeSN
To remedy the aforementioned shortcomings of the conventional methodology, we propose a new statistical model, Simple-BayeSN, describing the observed colormagnitude distribution as arising from the probabilistic combination of an intrinsic SN Ia color-magnitude variations and host galaxy dust reddening and extinction. The host galaxy dust is given a physically-motivated distribution, allowing for only positive extinction. The intrinsic color-magnitude slope β int can be different from the reddening-extinction slope R B . The observed data arises from the combination of these effects with measurement error. By fitting this statistical model to the light curve data, the separate physical characteristics of the intrinsic and dust distributions are coherently inferred. As our model uses the same SN Ia light curve measurements that are conventionally used by the Tripp formula, it can be readily applied to current cosmological SN Ia datasets.
We adopt a hierarchical Bayesian, or multi-level modeling, framework to build a structured probability model conceptually describing the multiple random effects that underlie the observed SN Ia. This principled strategy enables us to coherently model and make probabilistic inferences at both the level of an ensemble or population of objects as well as at the level of the constituent individuals (Gelman et al. 2003;Loredo & Hendry 2010;Loredo 2012). Inference with the hierarchical model may be regarded as a probabilistic deconvolution of the observed SN data into the multiple, unobserved, latent random effects generating it (Mandel 2012). Hierarchical Bayesian statistical modeling was first applied to SN Ia analysis by Mandel et al. (2009, who constructed the BayeSN model for optical and NIR SN Ia light curves . Fundamentally, BayeSN models the photometric time series observations of SN Ia as arising from intrinsic light curves and a host galaxy dust extinction across optical and NIR wave-lengths. The training process of BayeSN learns the covariance structure of the intrinsic SN Ia light curve distribution across phase and wavelength, as well as the characteristics of the dust distribution and dust law R V . The latent variables for each SN, and the hyperparameters of the population distributions are inferred coherently from the joint posterior density conditional on the set of light curve data. Mandel et al. (2011) demonstrated with BayeSN that the combination of optical and NIR observations could significantly improve constraints on SN Ia dust and the precision of photometric distances.
The Simple-BayeSN approach distills the core concept of BayeSN: hierarchically modeling the SN Ia data as a combination of intrinsic SN Ia variations, dust effects, and measurement error. However, BayeSN is a complex framework that directly and non-parametrically models the multi-wavelength photometric time series observations. As a simplification, Simple-BayeSN instead employs the outputs from an external light curve model that fits the photometric time series data of individual SN Ia and estimates the three parameters used in conventional analyses: peak apparent magnitude, apparent color, and light curve shape. In this paper, we employ the widely-used SALT2 light curve model, but Simple-BayeSN can generally work with the fit parameters from any external model for apparent SN Ia light curve data. March et al. (2011) constructed the first hierarchical Bayesian model for fitting the cosmological SN Ia Hubble diagram with the SALT2 optical light curve fit parameters. They encapsulated the Tripp formula in a hierarchical linear regression with population distributions for light curve shape and apparent color. The regression coefficients and the population hyperparameters are inferred simultaneously with cosmological parameters in the posterior distribution. This concept was further extended recently by Rubin et al. (UNITY, 2015) and Shariff et al. (BAHAMAS, 2016). However, these Bayesian models inherit the same fundamental, conceptual limitations of the Tripp formula, because they are essentially still regressing extinguished absolute magnitudes directly against apparent light curve parameters, rather than modeling the constituent, and physically distinct, intrinsic and dust components underlying the data. This paper is structured as follows. In §2, we construct a probabilistic generative model for the distribution of extinguished absolute magnitudes and apparent colors as a convolution of the intrinsic SN Ia color-magnitude distribution and the host galaxy dust distribution. We describe the features and generic implications of this model. In §3, we encapsulate this generative model in a hierarchical Bayesian framework, Simple-BayeSN, for analyzing SN Ia magnitudes, colors and light curve shape measurements. We demonstrate inference of the parameters of this model from the SN Ia data via maximum likelihood and Gibbs sampling. In §4 we apply this model to analyze a data set of SALT-II parameters for a sample of 248 nearby SN Ia (z < 0.10). In §4.5, we demonstrate how the host galaxy stellar mass dependence (Kelly et al. 2010) can be included in this new framework. We discuss our results in §5 and conclude in §6. Mathematical and computational details about our methods are described in Appendices §A, B, C, & D.

MOTIVATION: A PROBABILISTIC GENERATIVE MODEL
In this section, we illustrate the essential concepts underlying our statistical approach by constructing a probabilistic generative model for the SN Ia color-magnitude relation. We simulate a SN Ia sample of N SN = 250 SN Ia uniformly distributed between redshifts z = 0.01 and 0.10, and assume a fiducial ΛCDM cosmology of h = 0.72, Ω M = 0.27, Ω Λ = 0.73 and w = −1. We assume that the optical light curve data of each SN Ia is fit with a light curve model, which returns three useful measurements in the rest-frame of the SN: the peak Bband apparent magnitudem s , the peak apparent B − V colorĉ s =ĉ app s , and a light curve shape parameterx s , plus estimates of their fitting uncertainties. (The "hatted" quantities indicate estimated or measured values, which differ from the true values by some random measurement error).
To focus on understanding the effects of intrinsic colormagnitude variations and host galaxy dust, we will assume that the the observables, the apparent magnitude, apparent color and light curve shape, for each SN s are estimated without error, (m s ,ĉ app s ,x s ) = (m s , c app s , x s ), and that peculiar velocities are negligible: σ pec = 0 km s −1 . Hence, conditional on the redshift and cosmological parameters, the distance moduli are known, as are the extinguished absolute magnitudes: M ext s = m s − µ(z s ). In §3, we will construct a statistical model for analyzing observed data accounting for measurement error and peculiar velocities, inference of parameters and prediction of photometric distances.
To simulate the data, we generate values of the latent variables of intrinsic color c int s , intrinsic absolute magnitude M int s , light curve shape x s , and host galaxy dust reddening (or color excess) E s = E(B − V ) for each supernova s. In our usage, latent variables or parameters are those physical quantities for each supernova that are unobserved, but underlie the observed measurements. Intrinsic parameters refer to the properties of the supernova in the absence of host galaxy dust effects and measurement error. We assume simple forms for the population distributions underlying these latent variables. These distributions are governed by a set of hyperparameters, describing the intrinsic SN Ia population distribution and the host galaxy dust distribution. Data or parameters pertaining to individual supernovae are denoted by a subscript s, whereas hyperparameters common to the population are not. We adopt values that are similar to those estimated later in §4 by applying our model to the observed data.

Intrinsic Absolute Magnitudes and Colors
We generate light curve shape parameters {x s } for each SN, drawn from a Gaussian population distribution with mean x 0 = −0.40 and variance σ 2 s } is drawn from an assumed Gaussian population distribution with mean c int 0 = −0.06 mag and variance σ 2 c,int = (0.06 mag) 2 . We assume that the intrinsic absolute magnitude (in rest-frame B-band), M s , is simply related to these quantities by a linear relation:  where M int 0 is the expected intrinsic absolute magnitude for an (x s = 0, c int s = 0) supernova, α is the slope of the Phillips light curve width-luminosity relation, and β int is the slope against intrinsic color. We adopt true values of M int 0 = −19.40 mag, α = −0.15 and β int = 2.2. The random scatter around the mean relation is ǫ int s ∼ N (0, σ 2 int ) with variance σ 2 int = (0.10 mag) 2 . In Figure 1, we show the result of simulating these intrinsic supernova quantities from this generative model in this manner. Since we are focusing on effects in the color-magnitude plane, we control for light curve shape by plotting the light curve shape-corrected intrinsic absolute magnitude M int s − α × x s versus intrinsic color c int s .

Host Galaxy Dust Extinction and Reddening
As the light of the supernova leaves its host galaxy, it passes along the line-of-sight through a random column density of interstellar dust, which absorbs and scatters the light in wavelength-dependent process, resulting in extinction and reddening. The resulting change in rest-frame B − V supernova color is the color excess or dust reddening and is denoted The resulting dimming is quantified by the change in the peak absolute magnitude in B-band by the extinction A B . The dust extinction and color excess are related by the parameter R B , a property of the dust: A s B = R B E s . Since dust only dims and reddens, E s is a positive quantity. We assume that the dust reddening for each SN Ia is randomly drawn from an exponential population distribution with population mean τ = 0.07 mag: E s ∼ Expon(τ ). The exponential distribution for the dust population has been previously used by, e.g. Jha et al. (2007); Mandel et al. (2009Mandel et al. ( , 2014. The effect of dust is to dim the absolute magnitude and redden the color  The slope of the red extinction-reddening vectors is R B = 4.1. The length of each red vector is given by a random value of the reddening Es, drawn from the exponential distribution. We illustrate this for three random SN Ia (large points).
of SN s. We assume for simplicity that the same R B value characterizes the dust in all SN Ia host galaxies.
In Figure 2, we show the distribution of host galaxy dust reddening E s drawn from assumed exponential distribution. We adopt a value of R B = 4.1 and demonstrate the effect of reddening and extinction on the intrinsic magnitudes and color. Each intrinsic (blue) point maps to a red point through the effects of dust. In Fig

Implications of Inference with the Tripp Model
Next, we analyze the data by fitting the linear Tripp formula, as is conventionally done: where the residual scatter about the relation is ǫ t res,s ∼ N (0, (σ t res ) 2 ). To improve clarity, we introduce additional notation compared to Eq. 1. The absolute magnitude and color are actually the dust-extinguished absolute magnitude M ext s and dust-reddened apparent color c app s . The regression coefficients (M ext,t 0 , α t , β t app ) and the residual variance σ t res are labelled by superscripts t (for Tripp) to distinguish them from the hyperparameters of the generative model.
We estimate (M t 0 , α t , β t , σ t res ) from the simulated data via maximum likelihood. The resulting fit in the color-magnitude plane is shown as the black line, with slope β t = 3.23 ± 0.08. The residual scatter around the black line is σ t res = 0.127 mag. We find that when the true intrinsic magnitude-color slope is β int = 2.2, and the true dust law slope is R B = 4.1, the linear Tripp estimator obtains neither. Rather, it fits a value somewhere inbetween the true intrinsic slope and the dust law. For different values of the true hyperparameters, the value that β t app will obtain depends on the values of R B and β int , as well as the intrinsic color dispersion σ c,int , the intrinsic scatter σ int and the average amount of dust extinction R B τ in the host galaxies.
In the very red (positive) tail of the apparent color distribution, there are more red points below the fitted Tripp relation than above. This can also be seen in the blue (negative tail) of the apparent color distribution. This indicates that the linear Tripp model gives a biased trend of the extinguished absolute magnitude with apparent color for supernovae in the tails of the apparent color distribution. Furthermore, if the true mean trend of extinguished absolute magnitude vs. apparent color were actually linear, then the naive linear model would have captured this, resulting in no bias in the tails. This suggests that the true mean trend of M ext vs. c app , under the true generative model, must be non-linear as a function of c app (for a fixed light curve shape x).
We can calculate this trend given the mathematical model for the generative process we have just described. The extinguished absolute magnitudes and apparent colors result from adding dust extinction and reddening to the intrinsic absolute magnitudes and intrinsic colors. Therefore, the "dusty" distribution of extinguished absolute magnitudes and apparent colors is a convolution . Simulation: The intrinsic SN Ia magnitude-color distribution (blue points) is the same as that in Figs. 1 and 2, and has an intrinsic slope β int = 2.2. The effect of host galaxy dust drawn from an exponential distribution Es ∼ Expon(τ ) is to map each blue point (intrinsic magnitude and intrinsic color) to a red point (extinguished magnitude and apparent color) along the dust vector (red arrow) with slope R B = 4.1. When the conventional linear Tripp model is used to fit the relation between extinguished magnitudes and apparent colors, the resulting slope (black line) is β t = 3.23 ± 0.08 with residual scatter σ t res = 0.127 mag. The result of the conventional linear fit returns a value that is neither the true intrinsic slope, nor the dust slope, but a weighted average between the two. of the intrinsic SN Ia distribution and the host galaxy dust distribution 9 . From this convolved distribution, the mean trend of extinguished absolute magnitude with apparent color for a given light curve shape can be computed (see Appendix §A). The trend has the properties that in the blue limit (c app → −∞) its derivative smoothly approaches the intrinsic slope β int and in the red limit (c app → +∞) its slope approaches the dust law R B . The linear Tripp formula approximates an intermediate slope.
In the left panel of Figure 4, we plot the trend of M ext vs. c app , computed under our generative model (Simple-BayeSN or SBayeSN, red line), and the colormagnitude relation modeled by the conventional linear Tripp formula (black line) for an SN with an average light curve shape x = x 0 , using the values fit from the simulation. We see that the apparent color-magnitude slope β t app obtained by the linear Tripp estimator approximates the derivative of the true curve near the middle of the apparent color distribution. In the right panel, we show the difference between the trend under our model and the linear fit. We also overplot the true apparent color population distribution, Eq. A5 of Mandel et al. (2014), using the true values of (c int 0 , σ c,int , τ ). With respect to the nonlinear prediction of the true generative model, the naive linear model systematically overestimates the extinguished absolute luminosities, and thus the distances, of the SN Ia with very blue (negative) or very red (positive) apparent colors. It also slightly systematically underestimates the extinguished absolute luminosities, and thus the distances, of the SN Ia with average apparent color. Although these biases are small relative to the variance σ t res for one supernova, they will not decrease with increasing sample size.

Comparison to the Tripp Formula
In the absence of measurement error, we can substitute Eq. 5 and Eq. 4 into Eq. 3 to obtain By comparing this to Eq. 6, we see that if the dust law R B is equal to the intrinsic color-magnitude slope β int , then our model reduces to the Tripp formula, up to relabelling of the fit parameters, since (R B − β int )E s = 0 for every value of E s , regardless of the distribution of E s . In this case, intrinsic color effects are effectively indistinguishable from dust effects, at least within this model. The regression coefficients should match (M int 0 , α, β int ) = (M ext,t 0 , α t , β t app ) when fit to the data. If the intrinsic color-magnitude slope is differs from the dust law R B = β int , then the regression equation Eq. 7 is (up to relabelling of fit parameters) the Tripp formula plus an extra random variable (R B − β int )E s . If R B > β int , as we find in this paper, then this quantity is always positive. (At other wavelengths, it may be that R < β int , in which case it would always be negative). Since dust only reddens, E s > 0, this term could be regarded as an extra "noise" term that has a nonzero mean and an asymmetric distribution (unlike ǫ int s ,  which is assumed to be Gaussian and symmetric). Neglecting this additional random noise term will result in biased estimates of the other parameters, so that generally (M int 0 , α, β int ) = (M ext,t 0 , α t , β t app ). The simulation above demonstrates the inherently probabilistic nature underlying the apparent colormagnitude distribution. Consider, for example, a blue SN Ia with a well-measured light curve shape and apparent color c app s = −0.10 . There is some probability that it is unaffected by host galaxy dust (E s = 0) and its intrinsic color is also c int s = −0.10. But there is also some probability that its intrinsic color is actually c int s = −0.20, but it suffers from E(B − V ) s = +0.10 mag of host galaxy reddening. Since the intrinsic slope and the dust law are different, β int = R B , these two scenarios should result in different calculations for M ext s , which should in turn yield different distances µ s , given a well-measured apparent magnitude m s . The Tripp formula, and indeed, any statistical model directly modeling a functional form of extinguished magnitude vs. apparent color and light curve shape, would predict the same distance in either of these two scenarios. What is needed is a statistical model that properly weighs these different possibilities by their probabilities and marginalizes over them to produce a probability distribution for the photometric distance.
3. SIMPLE-BAYESN: A SIMPLE HIERARCHICAL BAYESIAN MODEL FOR SN IA In the previous section, we made some simple assumptions about the underlying probabilistic processes generating the SN Ia data, including the possibility that the intrinsic SN Ia magnitude-color slope β int is different from the dust law, R B . To test whether the data is consistent with two different slopes, we need to estimate the parameters governing the underlying processes. To do this, we first construct a hierarchical Bayesian statistical model to describe the distribution of the observed data as a probabilistic convolution of the the intrinsic variations and dust effects. We also include random effects such as measurement error and peculiar velocity uncertainties. We derive the likelihood function and joint posterior probability of the latent variables and hyperparameters. The hierarchical model is fit to the observed data with this model to estimate the hyperparameters using maximum likelihood and Gibbs sampling.
We regard as the observed data the peak apparent magnitudem s , the peak apparent colorĉ app s , and the light curve shapex s for each SN s, obtained from fitting the SN Ia light curve (time series) data, as well as the measured redshift z s . The latent parameters for each SN s are the observable light curve parameters φ s = (m s , c app s , x s ), the dust reddening E s and the distance modulus µ s . The hyperparameters of the dust distribution are the population mean dust reddening and the dust law parameter Θ dust = (τ, R B ). The hyperparameters Θ SN governing the intrinsic SN Ia correlations and distributions are described in §3.4.

Light Curve Fitting Error Likelihood
The observed data for conventional SN Ia analysis consists of broadband optical SN Ia light curves or brightness time series. A statistical model for the SN Ia time series ("light curve fitter") is fit to these data and returns estimates d s of parameters useful for summarizing the light curve, as well as an estimate of their joint estimation uncertainty, W s , for each SN s. The light curve fitter accounts for the redshifting of the SN Ia spectral energy distribution into the observed photometric passbands, and dust extinction from our Milky Way Galaxy. We refer to these light curve parameter estimates as the light curve "data," rather than the original multi-filter time series observations, and consider the error in these parameter estimates from the light curve fit as the "measurement error." Our likelihood function for light curve fitting and estimation models the probability of the data d s = (m s ,ĉ app s ,x s ) T given the latent observable parameters φ s as Gaussian with covariance matrix W s .
For the analysis in this paper, we use SN Ia light curve parameter estimates and error covariance matrices obtained from the SALT2 (Guy et al. 2010) light curve fitter. However, our hierarchical Bayesian statistical model can be used with parameter outputs from any suitable light curve fitter with the above characteristics.
3.2. Redshift-Distance Likelihood The expected theoretical distance modulus at redshift z in a smooth cosmology with parameters The approximation was made in linearizing the distance modulus at the cosmological redshift f (z c s ; Ω) about the observed redshift z s . The uncertainty in µ given the redshift is significant for low-z objects, for which where σ z is the redshift measurement error and σ pec = 200 km s −1 is the peculiar velocity dispersion (Carrick et al. 2015). With P (µ) ∝ 1, this leads to For the analysis of nearby SN Ia in this study, we fixed Ω =Ω = (0.72, 0.27, 0.73, −1) to its concordance values. It is important to remember that, because SN Ia by themselves are only relative distance indicators, absolute magnitudes and distance moduli are only determined up to an overall additive constant. Thus, a change in h will trivially shift all the absolute magnitudes and distances accordingly. However, a correct marginalization (either numerical or analytical) over the absolute magnitude constant M 0 (or the quantity M 0 = M 0 − 5 log 10 h) removes the dependence on h from all other inferences from the SN Ia data.

Latent Variable Equations
Defining the observable parameter vector as φ s = (m s , c app s , x s ) T , we can write this concisely as: where e E ≡ e 2 + R B e 1 and e i is a unit vector along the ith coordinate axis. It is not exactly true that the effect of applying a dust reddening law with a given extinction A V to SN Ia spectra will have linear effect on the magnitudes and colors measured through broadband filters. However, linearity is a good approximation up to moderate extinction values (Jha, Riess, & Kirshner 2007).

Intrinsic SN Ia Population Distribution
We construct a simple model for the joint population distribution of the intrinsic supernova parameters, P (ψ s | Θ SN ), depending on some hyperparameters Θ SN . It is convenient to express this joint distribution in terms of a product of conditional probability densities, each of which we model separately.
We can model each conditional factor with a mean relation and (Gaussian) scatter about the relation: and the marginal population distribution of the light curve shapes x as a Gaussian with mean x 0 and variance σ 2 x . We adopt the most general linear model for the mean intrinsic absolute magnitude trend f M (c int s , x s ; Θ SN ), and the mean intrinsic color trend f c (x s ; Θ SN ). Hence, the model for the intrinsic SN Ia parameters can be written as a hierarchy of linear equations.
where the intrinsic scatter terms are ǫ is similar to the linear Tripp formula but relates the intrinsic absolute magnitude to the light curve shape and intrinsic color.
To summarize, the nine hyperparameters governing the structure of the population distribution for the intrinsic SN Ia parameters are • M int 0 : the intrinsic absolute magnitude constant is the expected intrinsic absolute magnitude for a SN with light curve shape x s = 0 and intrinsic color c int s = 0, • α : the slope of the trend of intrinsic absolute magnitude vs. light curve shape, • β int : the slope of the trend of intrinsic absolute magnitude vs. intrinsic color, • σ 2 int : the intrinsic variance around the mean trend of intrinsic absolute magnitude vs. light curve shape and color, • c int 0 : the expected intrinsic color for a SN Ia with light curve shape x = 0. If α int c = 0, then c int 0 is the population mean intrinsic color, • α int c : the slope of the trend of intrinsic color vs. light curve shape, • σ 2 c,int : the intrinsic variance around the mean trend of intrinsic color vs. light curve shape, • x 0 : the mean of the x light curve shape population distribution, • σ 2 x : the variance of the x light curve shape population distribution.
These model choices can be expanded upon by including additional variables such as spectroscopic indicators or host galaxy properties, or non-linear trends in modeling the population distribution Eq. 15. For example, in §4.5, we modify Eq. 19 so that the intrinsic absolute magnitude constant depends on host galaxy stellar mass.

Host Galaxy Dust Population Distribution
The host galaxy dust reddening E s ≡ E(B − V ) is assumed to be drawn from an exponential population distribution with average τ : E s ∼ Expon(τ ). This has a probability density only on positive redding E s > 0 because dust only causes dimming and reddening:  found that this model describes well the distribution of peak apparent B − V colors of nearby SN Ia up to A V < 1. Futhermore, we assume that the properties of the dust in host galaxies are described by the ratio of extinction to reddening The dust hyperparameters are Θ dust = (τ, R B ): • τ : The population average of the exponential distribution of dust reddening: τ = E s .
• R B : The ratio of A B dust extinction to E(B − V ) reddening.
These simple choices can be expanded upon in the future to allow for the dust distribution or R B to vary within subpopulations. For example, in §4.5, we allow the average reddening τ to depend on host galaxy stellar mass.
3.6. Hyperpriors We need to specify the prior in the hyperparameters, or hyperprior P (Θ). For θ i = σ 2 int , σ 2 c,int , σ 2 x , or τ , we use flat priors on positive values only, θ i ∼ U (0, ∞). For all other hyperparameters we use flat priors on positive and negative values, θ i ∼ U (−∞, ∞). In our application, our posterior inferences are insensitive to using proper uniform hyperpriors over a wide range.

Bayesian Inference and Parameter Estimation
In Appendix §B, we use the assumptions of the model to derive the marginal likelihoods and posterior distribution of all the parameters and hyperparameters given the light curve data and redshifts. The conceptual relationships between the hyperparameters, latent parameters and data are depicted as a probabilistic graphical model in Fig. 16.
We use two methods to estimate the hyperparameters ("train" the model) from the light curve data D = {d s } and redshifts Z = {z s } for a fixed cosmological model Ω =Ω. The first method finds the peak of the log marginal likelihood P (d s | z s ; Θ,Ω) in the 11-dimensional parameter space Θ given by Eq. B4 using a constrained nonlinear optimization algorithm. The constraints are that σ x , σ c , σ int , τ > 0. The uncertainties in the maximum likelihood values Θ MLE are estimated from the inverse Hessian of the negative log likelihood (the Fisher information matrix), although these are only asymptotically lower bounds on the parameter uncertainties.
To explore the parameter space and obtain a more complete measure of the joint parameter uncertainty, we run a Gibbs sampler to generate a Markov chain that converges to the global posterior probability Eq. B2. This generates an MCMC in the 5N SN + 11 dimensional parameter space of {φ s , E s , µ s } and Θ. We have constructed a Gibbs sampling algorithm that exploits the conditional independence properties of the probabilistic graphical model to generate random moves that efficiently explore the parameter space. It is similar to the BayeSN algorithm of Mandel et al. (2009, except here the "data" are the 3 light curve parameter estimates rather than the full multi-wavelength time series data. The Gibbs sampler draws from the full set of conditional posterior densities, and does not require tuning of step sizes (which would be required for Metropolis algorithms). The resulting chains are used to compute numerical summaries of the posterior density. See Appendix §B and Appendix §C for more details.
In addition to using two separate methods for parameter estimation, maximum likelihood and MCMC sampling, we also implemented these algorithms in independent MatLab and Python codes. We tested the internal consistency of these codes by simulating data from the model under reasonable hyperparameter values, and then fitting the mock data with each code to verify that the truth was recovered.  (Contreras et al. 2010;Stritzinger et al. 2011). The light curve cuts applied to this sample are those used by Betoule et al. (2014). We demonstrate the Simple-BayeSN model by analyzing the nearby sample at low-z, where the distances µ(z) are insensitive to cosmological parameters. We will present an analysis of the full-z sample in a forthcoming paper.
In total, there are 248 SN Ia with 0.01 < z < 0.10 in this sample, with high quality light curves mainly from the CfA and CSP surveys. The sample includes 29 SN Ia observed both by CfA3/4 and CSP, and we arbitrarily selected the CSP light curves in these cases. The SALT2 light curve model was fit to the optical light curves. This light curve fitter returns estimates of the optical peak magnitude 10m B , the peak optical colorĉ (corresponding to B − V color), andx 1 ("stretch"), a measure of the optical light curve shape. We assign d s = (m B ,ĉ,x 1 ) T as our light curve data. SALT2 also returns an error covariance matrix for the light curve fit parameters, which we denote W s . The sample is cut to a range of light curve shapes of −3 <x 1 < 3 and a range of colors to −0.30 <ĉ < 0.30. The light curve shape cut corresponds approximately to the normal range in optical decline rate 1.6 ∆m 15 (B) 0.75 (Hicken et al. 2009a;Kessler et al. 2009a). These cuts exclude the very fast declining and red peculiar SN 1991bg-like objects, whose light curves SALT2 was not designed to fit ). The median error for the apparent magnitudesm B is 0.066 mag, the median error for colorsĉ is 0.030 mag, and the median error in the light curve shapex 1 is 0.121. Fig. 5 summarizes these data.

Fitting the Linear Tripp Model
First, we fit the conventional linear Tripp formula Eq. 6 to the data and estimated the parameters. We used a Bayesian linear regression that properly accounts for measurement error in the apparent magnitudesm B and covariates (ĉ app ,x), peculiar velocity uncertainties in the absolute magnitudes, as well as the residual variance around the regression σ 2 res (c.f. Appendix §D). We obtained M ext,t 0 = −19.299 ± 0.010 mag, 10 Kessler et al. (2013) describe m B as "an effective B-band magnitude is defined to be m B = −2.5 log 10 (x 0 ) + 10.635; this is the observed magnitude through an idealized filter that corresponds to the B-band in the rest frame of the SN." The overall flux amplitude of the SALT2 fit is x 0 . They also show with simulations that the SALT2 c parameter approximates the peak B − V color with scatter ≈ 0.02 mag.
The values of α t ≈ −0.15 and β t app ≈ 3 are typical for cosmological SN I analyses with the linear Tripp model (e.g. c.f. Rest et al. (2014)). The value σ t res = 0.124±0.009 mag is also typical, but we caution that this number should not be interpreted as the precision of distance estimates in the Hubble diagram. It corresponds to the portion of the scatter in the Hubble diagram that is unexplained by measurement errors or uncertainties in the SALT2 light curve fits. Hence, the realized scatter of estimated distance moduli in the Hubble diagram will be larger than this, owing the the light curve fit and peculiar velocity errors. In §4.4, we find that the rms scatter in the Hubble diagram is 0.16 mag for this sample.

Fitting the Simple-BayeSN model
Next, we estimated the hyperparameters of the hierarchical model using two numerical methods: optimizing the marginal likelihood and running an MCMC Gibbs sampler. The results (Table 1) from the two estimation methods agree quite well, but the posterior estimates typically provide more conservative uncertainty estimates than those obtained from the Fisher matrix evaluated at the maximum likelihood values.
The slope of the width-luminosity relation α remains the same. The model estimates the mean and standard deviation of the intrinsic color distribution to be c int 0 = −0.06 ± 0.01 and σ c,int = 0.067 ± 0.009. The population average host galaxy dust E(B − V ) reddening τ is non-zero and estimated τ = 0.07 ± 0.01 mag. From our posterior estimates, we estimate the population average host galaxy extinction to be A V = τ (R B − 1) = 0.184 ± 0.025 mag.
In Figure 6, we show a scatter plot of the light curve shape-corrected extinguished absolute magnitudes versus  Figure 7 shows the joint and marginal posterior probability densities for the hyperparameters R B and β int . They indicate that it is highly unlikely β int or R B are equal or that either takes on the value of β app ≈ 3 found by the linear Tripp model. In particular the inference of the host galaxy dust law slope R B = 3.8 ± 0.3 (R V = 2.7 ± 0.3) is consistent with the canonical Milky Way average R B = 4.1 (R V = 3.1). The intrinsic colormagnitude slope β int = 2.33 ± 0.26 is significantly (9σ) different from zero. The posterior mean of the difference R B − β int is positive and more than 3σ away from zero, while the tail probability that the difference is less than zero is 0.6%. Figure 8 shows the marginal posterior probability densities for the hyperparameters σ c,int and τ . Both hyperparameters are well constrained; in particular the population mean of the host galaxy dust distribution is clearly non-zero. Furthermore, the total apparent color variation comprises approximately equal contributions from intrinsic variation (σ c,int ≈ 0.07 mag) and dust reddening (τ ≈ 0.07 mag).

Latent Distributions
We explore the implications of our probabilistic model for the latent intrinsic and dust distributions. We examine the probabilistic properties of the distribution of the dusty latent variables that arise from the combination of intrinsic SN Ia variations and host galaxy dust. These aspects are computed using the fitted values of the hyperparametersΘ (Table 1), and the expressions   In Figure 9, we compare the data against the joint and marginal model probability distributions of extinguished absolute magnitudes and apparent colors described by the fitted model hyperpameters. These model distributions are derived from P (M ext , c app , x| Θ) as given in Eq. A2. The joint distribution of extinguished absolute magnitudes (corrected for light curve shape) and apparent colors accounts for a long redder-dimmer tail attributed to host galaxy dust with a different slope R B = 3.8 ± 0.3 than the intrinsic color magnitude slope β int = 2.3 ± 0.3. The model's marginal apparent color distribution captures the fat, red tail evident in a smoothed kernel density estimate of the apparent color data distribution. The model also captures an asymmetry in the absolute magnitude distribution. In these plots, the data distributions are slightly wider than the underlying fitted model distributions due to the effects of estimation errors in the light curve fit parameters and peculiar velocities.
In Figure 10, we illustrate the model intrinsic and dusty population distributions of the absolute magnitude and colors implied by the fitted hyperparametersΘ. The model intrinsic distributions are shown in blue. The population distribution of host galaxy dust reddening E s is well-described by an exponential distribution with mean τ = 0.07 ± 0.01 (red curve). The convolution of the intrinsic distribution (blue) with the dust distribution (red) yields the distribution of extinguished magnitudes and apparent colors (magenta).
In Figure 11, we further examine the model joint distributions of (light curve shape-corrected) absolute magnitude and color. The isodensity contours containing ap- proximately 95% and 68% of the model population probability are shown as dashed contours. The blue contours are the model's implied joint distribution intrinsic absolute magnitude and intrinsic color. The magenta contours depicts the model distribution of extinguished absolute magnitudes and apparent colors. Also shown are the intrinsic relation between absolute magnitude and intrinsic color, with slope β int (blue solid line), and a line with the slope of the dust law, R B (red solid line). The mean trend of extinguished absolute magnitude (corrected for light curve shape), as a function of apparent color is shown as the magenta curve. This trend is nonlinear with the following limiting properties. In the blue (negative) limit of the apparent color distribution, its derivative asymptotes to that of the intrinsic relation β int . In the red (positive) limit, it acquires the slope of the dust law R B . The curve smoothly transitions between the two limits at an apparent color value in the red (positive) tail of the intrinsic color distribution. SN Ia with apparent colors redder (more positive) than this transition are more likely to have been reddened by host galaxy dust rather than to be intrinsically red deviations from the population mean intrinsic color.

Distance Estimates
With the posterior estimates of the hyperparameters, we can calculate photometric distances from the Simple-BayeSN model and compare them to those obtained from the linear Tripp formula. These are obtained from the probability density P (µ s | d s ;Θ), for the photomet-  ric distance modulus based on the light curve data for each SN s, from which we can compute an expected valueμ s and varianceσ 2 µ,s ( §B.3). This distance uncertainty marginalizes over the tradeoffs between the latent factors of dust reddening-extinction and intrinsic colormagnitude variations. In Figure 12, we show a Hubble diagram of photometric distance moduli estimates obtained under the fit hyperparameters.
We compute the precision-weighted root-mean-squared Hubble residuals of the photometric distance moduli relative to the ΛCDM distance-redshift relation.
where the precision weights are w −1 s =σ 2 µ,s + σ 2 µ|z,s . For the conventional linear Tripp model, scatter of distance modulus estimates about the theoretical distanceredshift relation is about wRMS = 0.16 mag. Simple-BayeSN applied to the same SALT2 parameter estimates yields a minor improvement to wRMS = 0.15 mag. Figure 13 shows the Hubble residuals versus the estimated apparent colorĉ app for each SN for both the Tripp model and Simple-BayeSN. The conventional linear Tripp model on average overestimates the distance moduli for the SN Ia with very blue and very red apparent colors |c app | > 0.2. We fit a simple quadratic curve to the Hubble residuals versus apparent color, 2 p=0 b p (ĉ app ) p , and found b 0 = −0.015 ± 0.012, b 1 = −0.12 ± 0.11, and b 2 = 1.60 ± 0.65 (blue dashed curve). The quadratic coefficient b 2 is 2.5σ from zero. Simple-BayeSN reduces the distance bias (∼ 0.1 mag) in the tails of the appar-ent color distribution by fitting for the effective nonlinear trend generated by the convolution of intrinsic and dust effects. A quadratic fit of the Simple-BayeSN Hubble residuals vs. apparent color yields coefficients consistent with zero.
As shown in §2, if the true intrinsic color-magnitude relation is linear and has a slope different from the dust reddening-extinction vector, then the resulting "dusty" color-magnitude relation will be a smooth curve. Attempting to fit this curve with the linear Tripp model results in distance biases that increase as the SN Ia apparent color deviates from the mean, in either direction. In Figure 14, we plot the Hubble residuals as a function of the absolute deviation of the apparent color from the mean, |∆ĉ app | = |ĉ app −c app |. We also compute the mean Hubble residual within bins of width 0.05 mag. Under the linear Tripp formula, we see that the mean Hubble residual is significantly positive for SN Ia with very red or blue apparent colors |∆ĉ app | > 0.2 relative to the mean. However, the Simple-BayeSN model accounts for the β int = R B effect, and the Hubble residuals show no trend with |∆ĉ app |.
We tested the sensitivity of the model fit to the SN Ia with the most extreme colors, seen on the edges of Fig.  13. We narrowed our apparent color cut to |ĉ app | < 0.25 mag, and refit our model to the remaining 240 SN Ia. We obtained consistent hyperparameter posterior estimates: β int = 2.4±0.2, R B = 3.9±0.4, and R B −β int = 1.4±0.5. This shows that the fit for the intrinsic slope and dust law is not primarily driven by a few very blue or red SN Ia. The mean host galaxy reddening hyperparameter decreased to τ = 0.05±0.01, as expected since the narrower cut removed events with the most dust reddening.

Host Galaxy Mass Dependence
In this section, we investigate the dependence of the Hubble residuals on the host galaxy stellar masses, M stellar . Kelly et al. (2010) found that the Hubble residuals of SN Ia depended on host mass, after accounting for optical light curve shape and color correlations with  Figure 14. Hubble residuals versus the absolute deviation of apparent color from the mean, |∆ĉ app | = |ĉ app −c app | for each SN Ia for both the Tripp model (left) and Simple-BayeSN (right). We also plot the average Hubble residual in each bin of width 0.05 mag. The Tripp model systematically overestimates the distance moduli for the SN Ia with very blue and very red apparent colors |∆ĉ app | > 0.2. The SBayeSN model reduces the distance biases in the tails of the apparent color distribution by accounting for the specific nonlinear trend caused by the convolution of intrinsic and dust effects.
luminosity. This was confirmed in subsequent analyses: the Hubble residuals of SN Ia in more massive galaxies are brighter by about 0.05 − 0.10 mag (Sullivan et al. 2010;Lampeitl et al. 2010;Childress et al. 2013). The SN Ia host galaxies are divided between those with low stellar mass (LM) and high stellar mass (HM). We adopt the same mass split at M stellar = 10 10 M ⊙ used by Betoule et al. (2014).
The 0.01 < z < 0.10 sample contains host mass estimates for 215 SN Ia. We computed the Hubble residuals +3.740 ± 0.370 R B +3.633 ± 0.417 · · · · · · τ LM 0.097 ± 0.027 τ 0.068 ± 0.014 τ HM 0.066 ± 0.013 Note. -Hyperparameter estimates from fitting the Simple-BayeSN model to the 215 SN Ia in the 0.01 < z < 0.10 sample with host stellar mass estimates. The parameters and estimates on the right (left) are obtained by fitting Simple-BayeSN model with (without) host mass step dependence of the intrinsic absolute magnitude offset M int 0 and the host galaxy dust mean reddening τ . Estimates are the posterior mean and standard deviation of the MCMC samples.
obtained via the fitting the linear Tripp model to these SN Ia. The difference between the mean Hubble residuals of the HM and LM subsamples is −0.059 ± 0.028. Next, we computed the Hubble residuals obtained with the Simple-BayeSN model fitted to these SN Ia. The hyperparameter estimates are listed in Table 2. The difference between the mean Hubble residuals of the HM and LM subsamples is slightly smaller, −0.053 ± 0.028, with the Simple-BayeSN model.
Cosmological SN Ia analyses (Sullivan et al. 2011;Betoule et al. 2014) have incorporated this mass dependence into the conventional Tripp formula (Eq. 6) by modifying the absolute magnitude constant M ext,t 0 to be a ("mass step") function of M stellar : In our model, the hyperparameters express the intrinsic properties of the SN Ia and the host galaxy dust distribution separately. To account for a difference in the extinguished absolute magnitude offset between the two host mass classes, one can either make the intrinsic absolute magnitude offset M int 0 a function of host mass, or make the host dust distribution different between the two host mass classes. In fact, the net effect can be a combination of intrinsic SN Ia and dust population differences between the host mass classes. The simplest, sensible way to modify Simple-BayeSN to account for an overall offset in extinguished absolute magnitude is to We have modified our Gibbs sampler to sample the posterior distribution in the expanded parameter set and ran it on our sample to train the model and estimate the hyperparameters. The posterior mean estimates of the hyperparameters are shown in Table 2. We recomputed distances and the Hubble residuals obtained from the trained Simple-BayeSN model with host mass dependence. After including intrinsic M int 0 and dust τ host mass dependences in the model to compute distances, the difference between the mean Hubble residuals of the HM and LM subsamples is −0.005 ± 0.028 mag.
In Figure 15, we show the joint posterior density of the difference in the average host galaxy dust reddening δτ ≡ τ HM − τ LM and the intrinsic absolute magnitude offset δM int 0 = M int 0,HM − M int 0,LM . The peak of this distribution favors a solution that explains the Hubble-residual mass step partially by intrinsic SN Ia δM int 0 = −0.04 mag, and partially by dust distribution δτ ≈ −0.02 mag differences. However, as the size of this low-z sample is small, the uncertainties are wide; (δτ = 0, δM int 0 = 0) falls just inside the 95% highest posterior density contour, and the effects may be sensitive to the choice of the mass split. Future application of this model to a large high-z cosmological sample will yield more robust constraints on these host galaxy dependencies.

DISCUSSION
Our Simple-BayeSN model describes the observed "dusty" distribution of the extinguished absolute magnitudes and apparent colors of SN Ia as arising from the combination of intrinsic SN Ia color-magnitude variations (independent from light curve shape) and the distribution of host galaxy dust reddening-extinction. The probabilistic convolution of these two distributions generically results in a non-linear dusty apparent colormagnitude trend (Figs. 4,11). Previous empirical studies have noted evidence for trends of Hubble residuals with respect to apparent color in the conventional linear Tripp analysis, indicating that SN Ia with bluer or redder apparent colors may follow different effective colormagnitude slopes. Sullivan et al. (2011) analyzed SNLS3 SN Ia data using the linear Tripp formula and found β t app = 3.097 ± 0.094, but noted a possible trend of Hubble residuals blueward of c app < 0.15, suggesting bluer SN preferred a shallower effective color-magnitude slope. Puzzlingly, when they separately analyzed subsamples split at c app = 0, they found each β t app ≈ 3.8 to be higher than the global slope. Ganeshalingam et al. (2013) analyzed the combination of the Lick Observatory Supernova Search (LOSS) SN Ia sample (Ganeshalingam et al. 2010) with other low-z and high-z samples, using SALT2 and the linear Tripp method. They estimated a typical β t app ≈ 3.17 ± 0.08, but found significant positive Hubble residuals at apparent colors c app < −0.10, indicating that a shallower apparent color-magnitude trend (smaller β app ) would be a better fit to the bluer SN Ia. They also noted a slight positive deviation in the average Hubble residuals for c app > 0.3, but it was not statistically significant. Scolnic et al. (2014b) analyzed a compilation of SNLS3, SDSS-II and nearby samples with the linear Tripp relation and found that the Hubble residuals for SN Ia with c app bluer (more negative) and redder (more positive) than zero favored two different slopes, β app = 1.68 ± 0.38 and β app = 3.22 ± 0.29, respectively.
Analyzing the Union 2.1 data compilation, Suzuki et al. (2012) find a shallower β app = 1.3 ± 0.3 for bluer SN (c app < 0.05), and a steeper slope β app = 2.77 ± 0.09 for redder SN (c app > 0.05). Rubin et al. (2015) included a broken-linear apparent color-magnitude relation in their Bayesian UNITY model and found β app = 0.6 ± 0.3 for blue c app < 0, and β app = 2.8 ± 0.2 for red c app > 0, when they applied it to analyze Union 2.1. Notably, the two different apparent blue and red slopes are both much less than the R B = 4.1 of normal MW interstellar dust.
Fitting for different apparent color-magnitude slopes β app for apparently blue and red SN Ia is a useful empirical test for non-linearity of effective color-magnitude trend. However, as we pointed out in §2.4, the incorporation of any (linear or nonlinear) functional form s , x s ), directly between extinguished absolute magnitudes and apparent light curve parameters, within a statistical model for SN Ia data, would inherit the same conceptual flaw as the conventional linear Tripp formula. It would not properly weigh the probabilities of the different random combinations of intrinsic (M int , c int ) and dust (A B , E) that could have produced a given (M ext , c app ), and it would fail to capture the essential probabilistic nature of the latent astrophysical processes underlying the data.
By probabilistically deconvolving the SN Ia data into intrinsic variations and host galaxy dust components, we estimate a dust law with R B = 3.8 ± 0.3, consistent with the normal MW value of R B = 4.1. This is consistent with recent SN Ia analyses that do not use the Tripp formula. Burns et al. (2014), analyzing the optical-NIR colors of nearby SN Ia, found that those with low reddening E(B − V ) < 0.3 are consistent with dust reddening with R B = 4.1 (c.f. Phillips 2012). , accounting for intrinsic SN Ia covariance over phase and wavelength with the BayeSN hierarchical optical-NIR light curve model, found an R B ≈ 3.8 in the limit of low reddening. Correlating spectral equivalent widths with intrinsic photometric variability, Chotard et al. (2011) estimate R B = 3.8 ± 0.3. Sasdelli et al. (2016) developed a new technique to use principal components of SN Ia spectra to predict the intrinsic B light curves and B − V color curves of nearby SN Ia. By comparing these predictions against the observed light curves, they were able to deduce dust reddening and extinction and estimated R B = 3.8 ± 0.3. These studies have leveraged additional observational data (NIR photometry or spectra) to reach their conclusions. In this work, we have obtained a consistent and sensible estimate of R B using only optical photometric data fit with the same SALT2 light curve model used for cosmological SN Ia studies.
Our results for the dust law R B = 3.8 ± 0.3 are in quantitative agreement with the results of Scolnic et al. (2014b), who found a β ≈ 3.7 when attributing the residual scatter to solely to color σ cr within the SALT2mu framework (Marriner et al. 2011). One may be tempted to identify their residual scatter in color c r (Eq. 2) to Simple-BayeSN's intrinsic color c int . However, our approaches are conceptually different. Marriner et al. (2011) introduced the residual scatter matrix for additional dispersion in the light curve parameters that, in effect, acts in the same way as measurement errors around the primary magnitude, light curve shape and color relations of the usual Tripp model. In Simple-BayeSN, the intrinsic SN Ia distribution is fundamental: supernovae have physical intrinsic colors and luminosities before any dust reddening-extinction or measurement occurs ( §2).
There are also important practical differences. Before applying SALT2mu to fit for the Tripp coefficients, it is necessary to first assume the proportions of residual scatter attributed to luminosity and color, since it does not itself estimate them from the data. Further, Scolnic et al. (2014b) assumed the residual scatter in color is uncorrelated with luminosity. SALT2mu minimizes a "χ 2 ", a simplistic method with known frequentist biases in linear regression when dealing with scatter in the covariates (c.f. §D).
In contrast, Simple-BayeSN is a hierarchical Bayesian model that estimates the hyperparameters of the intrinsic SN Ia and dust component populations from the observed data. The host galaxy dust distribution is given a physically motivated form that only allows positive extinction. We find β int = 2.3 ± 0.3, indicating that the intrinsic colors of SN Ia are correlated with their intrinsic luminosity, for a given light curve shape. Hence, the "model" c mod and "residual" c r color components in Eq. 2 in the SALT2mu analysis of Scolnic et al. (2014b) do not exactly correspond to the dust reddening E(B − V ) and intrinsic color c int , respectively, in the Simple-BayeSN framework presented here. However, is it likely that combinations of the former map into combinations of the latter.
This variation may result from the effects of observational viewing angle into asymmetric SN Ia explosions (Kasen et al. 2009;Maeda et al. 2011). Theoretical simulations of multi-dimensional asymmetric delayed-detonation SN Ia explosions by Kasen, Röpke, & Woosley (2009) follow an intrinsic color-luminosity slope of β int = 4.45, controlling for light curve shape (in the absence of dust). Reaching quantitative agreement with the shallower slope found by our statistical model will be a challenge for theoretical models of SN Ia explosion physics.
We have used the SALT2 SN Ia spectral template model to fit the optical photometric time series to obtain estimates of the peak apparent m B , c (B − V color), and light curve shape x 1 for each SN Ia. Using Simple-BayeSN to analyze these derived parameter estimates and find two latent trends, β int and R B . However, internal to SALT2, there is a single, empirically-derived average color law CL(λ), describing the color variation in the SN Ia spectral energy distribution independent from light curve shape. Our results suggest that an improved SN Ia SED model should internally account for two physical sources of chromatic variation: dust reddening and intrinsic color variations.
Previous studies have found a correlation between the peak intrinsic B − V color and the optical decline rate, with a slope of ≈ +0.10 against ∆m15(B) (Phillips et al. 1999;Altavilla et al. 2004;Folatelli et al. 2010). This roughly translates to a slope of ≈ −0.014 against SALT2 stretch x 1 . Using SALT2 fit parameters, our estimate of the intrinsic color-light curve shape slope is α int c = −0.008 ± 0.005 (Eq. 20, Table 1) is not inconsistent with those estimates. However, it may be a shallower slope because internally the SALT2 model only has one color law and the model's peak B − V color has a negligible dependence on light curve shape x 1 (Kessler et al. 2013). Applying Simple-BayeSN to the light curve estimates from an improved fitting method may tighten the constraints on this and other key supernova parameters.
Our estimated R B = 3.8 ± 0.3 may in fact be somewhat biased low due to observational selection effects, which we have not modeled here. For any given apparent color c app , supernovae dimmer than the flux limits of nearby supernova searches will not be found and followed up. Hence, SN Ia with extinguished absolute magnitudes near the bottom portion of Fig. 6 would not be found, and their absence may cause the estimated R B vector to be slightly shallower. Although this mainly affects dust-reddened SN Ia, this bias is not the same as a cut in apparent color. Proper modeling the selection effects of nearby surveys is challenging, but such a correction would to tend to push R B closer to the MW mean 4.1.
Another potentially important interplay between our model and observational biases would occur at high-z, where selection effects are strong. The combination of a nonlinear apparent color-magnitude curve with highz selection biases may contribute to the apparent decrease in the estimated β app with redshift ("β evolution"). Kessler et al. (2009a)'s analysis of the SDSS-II data found that when the SALT2 parameters were fit with the Tripp formula to subsamples binned by redshift, the resulting estimates of β app tended to decrease significantly at z > 0.6. Using an updated SALT2, Guy et al. (2010) also found a decreasing estimated β app at higher z, but did not conclude that it was real, due to the possibility of systematic errors in the estimation of SN Ia color and its uncertainties confounding the determination of β app at high-z. While Betoule et al. (2014) did not report any β evolution, recent analyses of their JLA compilation by Li et al. (2016) and Shariff et al. (2016) have claimed decreasing β app at high-z. For example, Shariff et al. (2016) report an unexplained drop from β app = 3.1 ± 0.1 by ∆β app = −1.1 ± 0.2 at z ≈ 0.66.
We demonstrated in §2 that, if the data is generated from intrinsic variation and dust distributions with different β int and R B , then their probabilistic convolution results in a nonlinear mean apparent color-magnitude trend. When the observed data is fit with the linear Tripp formula, it yields an estimated β t app that approximates the derivative of this curve near the mean of the apparent color distribution (Fig. 4). At high-redshifts, the dimmer, and thus redder and dustier SN Ia, are less likely to be found and included in the sample. This will cause the mean apparent color of the observed sample to shift to bluer (more negative) values from low to high redshift. The observed high-z sample will be, in effect, less dusty causing the apparent color-magnitude curve to be shallower overall. The blue shift in the mean apparent color at high redshift will cause the β t app to approximate the derivative of the underlying curve at a bluer (more negative) mean apparent color, where the derivative is closer to the intrinsic slope β int (Fig. 11). The net result will naturally cause the Tripp estimator of β t app to obtain a lower value at high-z relative to low-z.
In the context of our results, our estimate β t app = 3.01±0.12 (Eq. 26) is consistent with the β t app ≈ 3 found at z 0.7 by Li et al. (2016) and Shariff et al. (2016), and our β int = 2.3 ± 0.3 is consistent with β t app ≈ 2 they find at higher redshifts. This is consistent with our understanding above that the SN Ia at high-z will be on average bluer and less dusty due to selection bias, and thus should adhere closer to the shallow intrinsic colorluminosity trend β int . This explanation rests on the relative proportions of intrinsic color to dust reddening, in the observed sample, changing with redshift due to selection bias, and does not require either the intrinsic slope β int or dust law R B themselves to "evolve" with redshift. Disentangling this generic prediction of our model from the possible systematic errors particular to the fitting of noisy high-z SN Ia light curve data will be a challenge.
Several empirical studies have looked for astrophysical systematic effects by examining the distribution of SN Ia parameters, or the Hubble residuals, as a function of host galaxy properties or local environment of the SN (e.g. Kelly et al. 2010;Lampeitl et al. 2010;Sullivan et al. 2010;Kelly et al. 2015). For example, Sullivan et al. (2011) find a 4.4σ difference in β t app between SN Ia host galaxies with low and high stellar mass. They note that the estimated β t app is likely conflating dust effects with intrinsic variations, and hypothesize that this difference is tied to the relative dustiness of host galaxies as a func-tion of stellar mass and star formation rate. Because of the confounding of intrinsic variation and dust, the Tripp formula is a blunt tool to examine these hypotheses.
In contrast, our framework is well-suited to these kinds of analyses because we have separately modeled the latent components of the data attributed to intrinsic variation and host galaxy dust, and parameterized their properties (e.g. β int , R B , σ c,int , τ ). This provides a richer vocabulary with which to investigate the astrophysical connections between the SN observables and the properties of the progenitor, local environments or host galaxies. For example, in §4.5 we incorporated the host galaxy "mass step" into our model in a new way: we allow a combination of the intrinsic SN Ia brightness and host galaxy dustiness to vary with host stellar mass. This approach may lead to insights into the physical origin of the observed effect, e.g. whether it is caused by the the host galaxy dust properties and/or the physics of SN Ia progenitors. This is just one simple example of how our framework can be used to investigate astrophysical systematics. Future studies could investigate the dependence of our intrinsic and dust parameters on star formation rate, host galaxy stellar mass, metallicity and morphology.
6. CONCLUSION The Tripp formula (Eq. 1) is widely used in conventional analyses of cosmological SN Ia light curve data. However, it is also simplistic: by directly regressing extinguished absolute magnitudes against apparent color, it fails to take into account that both factors comprise the physically distinct effects of intrinsic SN Ia variation and extrinsic host galaxy dust. This shortcoming has led to estimates of the apparent color-magnitude slope β app that are puzzlingly smaller than the normal MW interstellar dust reddening-extinction law R B = 4.1.
To address this, we have constructed a hierarchical Bayesian model (Simple-BayeSN) describing the "dusty" distribution of extinguished absolute magnitudes and apparent colors as arising from two population distributions: an intrinsic SN Ia distribution, and a host galaxy dust distribution. The intrinsic distribution includes the dependence of intrinsic luminosity on light curve shape and intrinsic color (Eq. 19). It allows for a non-zero trend of intrinsic absolute magnitude vs. intrinsic color with slope β int , controlling for light curve shape. It also models intrinsic color variations that can be uncorrelated with light curve shape (Eq. 20). The host galaxy dust distribution includes a law parameter, R B , characterizing the direction of the dust extinctionreddening vector. This model provides a more physical decomposition of the sources of variation underlying the SN Ia data. Inference with this model in effect performs a probabilistic deconvolution of the observed SN Ia measurements to estimate the characteristics of the underlying intrinsic and dust component distributions.
Analyzing the optical light curve fit data from compilation of 248 nearby SN Ia (z < 0.1), we find that fitting the Tripp formula obtains β t app = 3.0 ± 0.1, significantly less than the R B = 4.1 of normal MW dust, but consistent with the findings of recent cosmological analysis. In contrast, Simple-BayeSN, by modeling the data as probabilistic convolution of intrinsic and dust components, finds a non-zero intrinsic color-magnitude slope β int = 2.3 ± 0.3 and a dust law slope of R B = 3.8 ± 0.3. The slope of the dust law is consistent with the average value for normal Milky Way interstellar dust R B = 4.1. The width of the intrinsic B − V color distribution is found to be σ c,int ≈ 0.07 mag, while the average dust E(B − V ) reddening of the sample is τ ≈ 0.07 mag.
Since β int = R B (at 3σ), the convolution of the intrinsic SN Ia color-magnitude distribution with the dust reddening-extinction distribution results in a non-linear curve of extinguished absolute magnitude vs. apparent color (Fig. 4 & 11). The conventional linear Tripp formula approximates this curve near the bulk of the empirical apparent color distribution of the samples. It obtains a linear slope that approximates the slope of a tangent to this curve near the mean apparent color. This results in ∼ 0.1 mag overestimates of photometric distance moduli in the tails of the apparent color distribution. This systematic bias vanishes when we use our model to account for the distinct effects of intrinsic variation and host galaxy dust.
As photometric calibration uncertainties become better understood and controlled (Scolnic et al. 2015), astrophysical systematics caused by incorrect modeling of SN Ia color-luminosity effects will become a major limiting factor for dark energy constraints from SN Ia. Future research will evaluate the relative impacts of these systematics on the constraints on w from the cosmological SN Ia sample. Since the Tripp formula is also used for calibrating nearby SN Ia on the absolute distance scale, similar color-dependent systematic errors may impact the inference of H 0 . However, the effect on the H 0 estimate is likely to be small, as it is mainly influenced by the measurement of the average luminosity of the low-z SN Ia, whereas these biases are most pronounced in the color tails.
There are many potential directions for applying and extending the Simple-BayeSN framework. We will apply this model for the cosmological analysis of SN Ia over the full range of redshifts, accounting for systematics and selection effects, for determination of the cosmological parameters. While we have applied Simple-BayeSN to model intrinsic and dust effects using parameters derived from optical photometric light curve fits, we can incorporate other useful information that is likely to improve inferences. For example, additional measurements from NIR light curves (e.g. Friedman et al. 2015) is likely to improve constraints on dust and distances. We can also extend the model to incorporate spectroscopic information. For example, the spectroscopic expansion velocitycolor relation (VCR, Foley & Kasen 2011) can be modeled by adding velocity as an additional variable that may correlate with intrinsic color and luminosity. Non-Gaussian distributions of the intrinsic parameters can be tested (Mandel et al. 2014). The Simple-BayeSN hyperparameters could be used to explore astrophysical systematics related to the dependence of intrinsic SN Ia properties (e.g. β int and σ c,int ) and dust (τ, R B ) with properties of the host galaxy and local environment. We will incorporate this model into SNANA (Kessler et al. 2009b) to use in realistic simulations of SN Ia surveys.
Precise and accurate SN Ia distance estimates are essential to the success of current or future cosmological surveys such as the Dark Energy Survey (DES), the Large Synoptic Survey Telescope (LSST) survey, and WFIRST. Simple-BayeSN is both a conceptual advance and practical improvement in the proper statistical modeling, inference and understanding of the intrinsic SN Ia variations and host galaxy dust effects underlying these measurements.
We thank Saurabh Jha, Rick Kessler, Pat Kelly, Xiao-Li Meng, David Spergel, and Roberto Trotta for useful discussions. Supernova cosmology at the Harvard College Observatory is supported in part by National Science Foundation grants AST-156854, AST-1211196, and NASA grant NNX15AJ55G. This manuscript is based upon work supported by the National Aeronautics and Space Administration under Contract No. NNG16PJ34C issued through the WFIRST Science Investigation Teams Program. R.J.F. and D.S. were supported in part by NASA grant 14-WPS14-0048. R.J.F.'s UCSC group is supported in part by NSF grant AST-1518052 and the Alfred P. Sloan Foundation. D.S. acknowledges support from KICP and from NASA through Hubble Fellowship grant HST-HF2-51383.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. H.S. was supported by a Marie-Skodowska-Curie RISE (H2020-MSCA-RISE-2015-691164) Grant provided by the European Commission. This work was supported by Grant ST/N000838/1 from the Science and Technology Facilities Council (UK).

A. THE DUSTY SN IA DISTRIBUTION
The joint probability distribution of the dusty SN Ia parameters P (M ext , c app , x| Θ) can be derived from the model assumptions for the intrinsic SN Ia population distribution ( §3.4) and the dust population model ( §3.5). This depends on the their respective hyperparameters: Θ SN , and Θ dust = (R B , τ ). Letψ = (M ext , c app , x s ) T be a vector of the dusty parameters of a SN Ia, modified by extinction and reddening. This is related to the intrinsic SN Ia parameters ψ viaψ = ψ + e E E, where e E = e 2 + R B e 1 . The joint distribution ofψ and dust reddening E is P (ψ, E| Θ) = P (ψ =ψ − e E E| Θ SN , R B )P (E| τ ).
The marginal probability distribution ofψ is a convolution of the intrinsic SN Ia population distribution and the dust distribution for a given set of model hyperparameters Θ = {Θ SN , R B , τ }, obtained by this integral: From this, one can compute the conditional probability P (M ext | c app , x, Θ) and then the mean trend of extinguished absolute magnitude versus apparent color and light curve shape, E(M ext | c app , x, Θ). The joint distribution of the light curve shape-corrected extinguished absolute magnitude M ext − αx and apparent color is mathematically equivalent to evaluating the above expressions with x = 0.

B. BAYESIAN INFERENCE B.1. Global Posterior Probability Density
To estimate all the parameters and hyperparameters, we compute the joint posterior density. For one SN s, the joint probability of the data d s and latent parameters φ s , E s , and µ s given redshift z s , the hyperparameters Θ = (Θ SN , Θ dust ) and cosmological parameters Ω is: where P (Ω) includes constraints from other data (e.g. CMB or BAO). However, we do not do this in this paper, since we are restricting our analysis to the low-z sample with Ω fixed.
B.2. Probabilistic Graphical Model In Figure 16, we display a directed acyclic graph, a probabilistic graphical representation of our hierarchical Bayesian model. Probabilistic graphical models were first used to express hierarchical Bayesian inference with SN Ia by Mandel et al. (2009. The probabilistic graphical model describes how the unknown parameters of individual SN Ia (labelled by index s) and the hyperparameters of the dust and intrinsic SN Ia populations are related to the measured supernova data {d s , z s }, and cosmological parameters Ω. From the intrinsic SN Ia population, described by hyperparameters Θ SN , a vector of intrinsic light curve parameters ψ s is drawn for each SN s. Random values of dust reddening E s and extinction A s B for each SN are drawn from the host galaxy dust population, described by hyperparameters R B , τ . The intrinsic and dust latent variables combined with distance modulus and light curve fitting error yield the estimated light curve parameters. The distance modulus is related to the observed redshift through the cosmological parameters and peculiar velocity error. In this paper, we have fixed the cosmological parametersΩ for the low-z analysis. The open boxes represent unknown parameters and hyperparameters, the shaded boxes represent the observed data, and the arrows indicate relations of conditional probability between the parameters, hyperparameters and data. The arrow colors heuristically represent the multiple sources of randomness and uncertainty underlying the data: intrinsic variation (blue), dust (red), distances and peculiar velocities (green), measurement errors (purple) and cosmological parameters (black). One can read the graph as a describing the generative process for the data.
The probability density of the photometric distance modulus is P (µ s |d s ;Θ) ∝ P (d s | µ s ;Θ)P (µ s ), where the prior P (µ s ) is taken to be flat. Ideally, one would also marginalize over the posterior uncertainty in the hyperparameters Θ, but if they are well-determined, this may not be worth the extra effort.
C. SIMPLE-BAYESN GIBBS SAMPLING ALGORITHM We sketch an MCMC Gibbs sampling algorithm to sample from the global posterior Eq. B2 of the Simple-BayeSN hierarchical model. The purpose of an MCMC algorithm is to generate a sequence of random parameter vectors with a long-run stationary distribution equal to the global posterior. Our Gibbs sampler proceeds by sequentially drawing new parameter values from a full set of conditional posterior densities derived from the global posterior distribution. Gibbs samplers for hierarchical Bayesian models for SN Ia were previously developed by Mandel et al. (2009Mandel et al. ( , 2014 and recently by Shariff et al. (2016).
We begin a chain with randomized initial values for the apparent colors and distances {φ s , µ s }, as well as the population hyperparameters Θ SN , R B , and τ , overdispersed around the maximum likelihood values Θ MLE . We alternate between updating the individual SN parameters {φ s , E s , µ s } conditional on the hyperparameters, and updating the population hyperparameters (Θ SN , R B , τ ) conditional on the current values of the set of individual SN parameters. Steps 1-4 utilize draws from tractable conditional distributions.