Suppose you have a set of N₁ measurements obtained with one detector, and another set of N₂ measurements obtained with a second detector. And let’s say you wanted something as simple as an estimate of the mean of the quantity (say the intensity) being measured. Let us further stipulate that the measurement errors of each of the points is similar in magnitude and neither instrument displays any odd behavior. How does one combine the two datasets without appealing to subjective biases about the reliability or otherwise of the two instruments?

We’ve mentioned this problem before, but I don’t think there’s been a satisfactory answer.

The simplest thing to do would be to simply pool all the measurements into one dataset with N=N₁+N₂ measurements and compute the mean that way. But if the number of points in each dataset is very different, the simple combined sample mean is actually a statement of bias in favor of the dataset with more measurements.

In a Bayesian context, there seems to be at least a well-defined prescription: define a model, compute the posterior probability density for the model parameters using dataset 1 using some non-informative prior, use this posterior density as the prior density in the next step, where a new posterior density is computed from dataset 2.

What does one do in the frequentist universe?

[Update 9/30] After considerable discussion, it seems clear that there is no way to do this without making some assumption about the reliability of the detectors. In other words, disinterested objectivity is a mirage.

The motivation of writing this posting was originated to Vinay’s recommendation: Practical Statistics for Astronomers (J.V.Wall and C.R.Jenkins), which provides many statistical insights and caveats that astronomers tend to ignore. Without looking at the error distribution and the properties of data, astronomers jump into chi-square and correlation. If someone reads the book, he/she will be careful on adopting statistics of common practice in astronomy, developed many decades ago, and founded on strong assumptions, not compatible with modern data sets. The book addresses many concerns that have been growing in my mind for astronomers and introduces various statistical methods applicable in astronomy.

The view points of astronomers without in-class statistics education but with full readership of this book, would be different from mine. The book mentioned unbiasedness, consistency, closedness, and robustness of statistics, which normally are not discussed nor proved in astronomy papers. Therefore, those readers may miss the insights, caveats, and contents-between-the-lines of the book, which I care about. To reduce such gap, as for quick and easy understanding of classical statistics, I recommend Cartoon Guide to Statistics (Larry Gonick, Woollcott Smith Business & Investing Collins) as a first step. This cartoon book enhances fundamentals in statistics only with fun and a friendly manner, and provides everything that rudimentary textbooks offer.

If someone wants to know beyond classical statistics (so called frequentist statistics) and likes to know popular Bayesian statistics, astronomy professor Phil Gregory’s Bayesian Logical Data Analysis for the Physical Sciences is recommended. If one likes to know little bit more on the modern statistics of frequentists and Bayesians, All of Statistics (Larry Wasserman) is recommended. I realize that textbooks for non-statistics students are too thick to go through in a short time (The book for senior engineering students at Penn State I used for teaching was Probability and Statistics for Engineering and the Sciences by Jay. L Devore, 4th and 5th edition and it was about 600 pages. The current edition is 736 pages). One of well received textbooks for graduate students in electrical engineering is Probability, Random Variables and Stochastic Processes (A. Papoulis & S.U. Pillai). I remember the book offers a rather less abstract definition of measure and practical examples (Personally, Hermite polynomials was useful from the book).

For a casual reading about statistics and its 20th century history, The Lady Tasting Tea: How Statistics Revolutionized Science in the Twentieth Century (D. Salsburg) is quite nice.

Statistics is not just for best fit analysis and error bars. It is a wonderful telescope extracts correct information when it is operated carefully to the right target by the manual. It gets rid of atmospheric and other blurring factors when statistics is understood righteously. It is not a black box nor a magic, as many people think.

The era of treating everything gaussian is over decades ago. Because of the central limit theorem and the delta method (a good example is log-transformation), many statistics asymptotically follows the normal (gaussian) distribution but there are various families of distributions. Because of possible bias in the chi-square method, the error bar cannot guarantee the appointed coverage, like 95%. There are also nonparametric statistics, known for robustness, whereas it may be less efficient than statistics of distribution family assumption. Yet, it does not require model assumption. Also, Bayesian statistics works wonderfully if correct information on priors, suitable likelihood models, and computing powers for hierarchical models and numerical integration are provided.

Before jumping into the chi-square for fitting and testing at the same time, to prevent introducing bias, exploratory data analysis is required for better understanding data and for seeking a suitable statistic and its assumptions. The exploratory data analysis starts from simple scatter plots and box plots. A little statistical care for data and good interests in the truth of statistical methods are all I am asking for. I do wish that these books could assist the realization of my wishes.

—————————————————————————-
[1.] Most of links to books are from amazon.com but there is no personal affiliation to the company.

[2.] In addition to the previous posting on chi-square, what is so special about chi square in astronomy, I’d like to mention possible bias in chi-square fitting and testing. It is well known that utilizing the same data set for fitting, which results in parameter estimates so called in astronomy best fit values and error bars, and testing based on these parameter estimates brings out bias so that the best fit is biased from the true parameter value and the error bar does not match the aimed coverage. See the problem from Aneta’s an example of chi2 bias in fitting x-ray spectra

[3.] More book recommendation is welcome.

If the bootstrap is an automatic processor for frequentist inference, then MCMC is its Bayesian counterpart.

Sometime in my second year of studying statistics, I said that bootstrap and MCMC are equivalent but reflect different streams in statistics. The response to this comment was ‘that’s nonsense.’ Although I forgot details of the circumstance, I was hurt and didn’t try to prove myself. After years, the occasion immediately floats on the surface upon seeing this sentence.

In general, Cosmologists are Bayesians and High Energy Physicists are Frequentists.

I thought it was opposite. By the way, if you crave for more papers, click

• [astro-ph:0712.2544]
  RHESSI Microflare Statistics II. X-ray Imaging, Spectroscopy & Energy Distributions I. G. Hannah et.al.

• [stat.AP;0712.2708]
  The Banff Challenge: Statistical Detection of a Noisy Signal A. C. Davison & N. Sartori

• [astro-ph:0712.2898]
  A study of supervised classification of Hipparcos variable stars using PCA and Support Vector Machines P.G. Willemsen & L. Eyer

• [astro-ph:0712.2961]
  The frequency distribution of the height above the Galactic plane for the novae M. Burlak

• [astro-ph:0712.3028]
  A practical guide to Basic Statistical Techniques for Data Analysis in Cosmology L. Verde

• [astro-ph:0712.3049]
  ZOBOV: a parameter-free void-finding algorithm M. C. Neyrinck

• [stat.CO:0712.3056]
  Gibbs Sampling for a Bayesian Hierarchical Version of the General Linear Mixed Model A. A. Johnson & G L. Jones

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 Y_i(X_i) {i=1..N_bin} and fit a model f(A) with parameters A (usually by minimizing the chisquare (χ²) 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 Y_ij(X_i), {i=1..N_bin, j=1..N_obs}, and deriving Â_j for each of them, and computing sigma_A²=variance( {Â_j} ). By the central limit theorem, as N_obs increases, the distribution of { Â_j } will tend closer to a Gaussian. So, in this case, 3*sigma_A will asymptotically become equal to 3sigma_A.

However, Astronomers are not afforded this luxury of having multiple observations. With dreary regularity, we have to make do with N_obs=1. In that case, one takes recourse in the behavior of the χ², and in particular, the likelihood of obtaining a certain χ² 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 χ² 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 χ² function (computed by stepping through the values of one parameter and varying all the other parameters to minimize the χ² at that point) increases by 1. Similarly, the 3sigma error is computed by checking at what values the χ² increases by 9. As long as the assumptions of local linearity and the Gaussian measurement error hold (i.e., the χ² 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.