variations [@Nobili2009; @Yasuda2010], K-correction and template calibration [@Nugent2002; @Hsiao2007], flux calibration [@Faccioli2011], inhomogeneities (peculiar velocities [@Hui2006; @Davis2010], gravitational lensing [@Wang2005b; @Kronborg2010], etc). We would also like to call attention to the particularly clear, informative and up-to-date generic review articles on SNe Ia by @Kirshner2010, @Howell2011, and @Goobar2011. Simultaneously with this physical endeavor from first principles, we should also exercise our best statistical consistent practices to analyze and mine the data and to test the robustness of our inferences. It is to the latter that our paper is devoted.

In contrast to the traditional $\chi^2$ approach present in most of the literature, in this work we discuss a different approach to parameter fitting based on the likelihood. We call attention to the fact that, when we want to estimate both the covariance and the mean of a Gaussian process, the ordinary (uncorrected) $\chi^2$ approach (or any iterative recipe therefrom, for that matter) cannot be straightforwardly applied, lest we might lose a nontrivial term in the objective function to be extremized. This is necessary when the covariance itself does depend on free parameters of the underlying model. Recently, similar criticisms to the traditional $\chi^2$ method were presented by Kim [@Kim2011] (see also 