The AstroStat Slog » HPD http://hea-www.harvard.edu/AstroStat/slog Weaving together Astronomy+Statistics+Computer Science+Engineering+Intrumentation, far beyond the growing borders Fri, 09 Sep 2011 17:05:33 +0000 en-US hourly 1 http://wordpress.org/?v=3.4 Is 3sigma the same as 3*1sigma? http://hea-www.harvard.edu/AstroStat/slog/2007/3-times-sigma/ http://hea-www.harvard.edu/AstroStat/slog/2007/3-times-sigma/#comments Thu, 08 Feb 2007 20:29:22 +0000 vlk http://hea-www.harvard.edu/AstroStat/slog/2007/3-times-sigma/ Sometime early this year, Jeremy Drake asked this innocuous sounding question in the context of determining the bounds on the amplitude of an absorption line: Is the 3sigma error bar the same length as 3 times the 1sigma error bar?

In other words, if he measured the 99.7% confidence range for his model parameter, would it always be thrice as large as the nominal 1sigma confidence range? The answer is complicated, and depends on who you ask: Frequentists will say “yes, of course!”, Likelihoodists will say “Maybe, yeah, er, depends”, and Bayesians will say “sigma? what’s that?” So I think it would be useful to expound a bit more on this to astronomers, whose mental processes are generally Bayesian but whose computational tools are mostly Frequentist.

In the Frequentist paradigm, one fits a parameterized model to numerous datasets, and computes the errors on the fitted parameters by looking at the ensemble of the parameter values. That is, if you have one dataset Yi(Xi) {i=1..Nbin} and fit a model f(A) with parameters A (usually by minimizing the chisquare (χ2) deviations or something of that sort) you will get one set of best-fit parameters, Â. It is important to note that in this paradigm, the concept of determining an error on  does not exist!. Errors are only determined by fitting f(A) to a series of datasets Yij(Xi), {i=1..Nbin, j=1..Nobs}, and deriving Âj for each of them, and computing sigmaA2=variance( {Âj} ). By the central limit theorem, as Nobs increases, the distribution of { Âj } will tend closer to a Gaussian. So, in this case, 3*sigmaA will asymptotically become equal to 3sigmaA.

However, Astronomers are not afforded this luxury of having multiple observations. With dreary regularity, we have to make do with Nobs=1. In that case, one takes recourse in the behavior of the χ2, and in particular, the likelihood of obtaining a certain χ2 for the observed data given the model. If we assume that the measurement errors are Normally distributed, and that we are looking in that part of the parameter space close to the best-fit solution, then it can be shown that the χ2 varies as a parabolic function of the parameter values. The 1sigma error on the parameter (not necessarily from a linear model; see Cash 1979, for example) is computed by determining the parameter values at which the χ2 function (computed by stepping through the values of one parameter and varying all the other parameters to minimize the χ2 at that point) increases by 1. Similarly, the 3sigma error is computed by checking at what values the χ2 increases by 9. As long as the assumptions of local linearity and the Gaussian measurement error hold (i.e., the χ2 function is parabolic), then we will have 3*1sigma=3sigma, otherwise not.

Bayesians on the other hand compute the marginalized posterior probability density for each parameter and calculate the bounds on the parameter via credible regions that enclose a given percentage of the area under the curve. There is no set recipe for specifying how the bounds are calculated: it could be one-sided, centered on the median, or include the mode and the highest posterior densities. (The last leads to the smallest possible interval.) The percentage levels for which the intervals are calculated are usually matched to the percentages that Gaussian standard deviations provide, e.g., 68% == 1sigma, 99.7% == 3sigma, etc. Thus, when the 68% credible range is mentioned, it could also be referred to as the “Gaussian-equivalent 1sigma region”. However, because the intervals can be defined in many different ways, they do not necessarily match the 1sigma error bar, and for that matter are often asymmetrical. Given that, it is quite unreasonable to expect that 3*1sigma, even on one side, would match the Gaussian-equivalent 3sigma. IF the pdf is Gaussian, AND we choose to use HPD (highest posterior density) intervals, THEN the 99.7% interval would be thrice as long as the 68% interval.

]]>
http://hea-www.harvard.edu/AstroStat/slog/2007/3-times-sigma/feed/ 3