{'fraction_non_alphanumeric': 0.033710199115700765, 'fraction_numerical': 0.011691741594706075, 'mean_word_length': 4.60304264693657, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Estimation of Hardness Ratios: Modeling and Computations

A commonly used measure to summarize the nature of a photon spectrum is the so-called Hardness Ratio, which compares the number of counts observed in different passbands. The hardness ratio is especially useful to distinguish between and categorize weak sources as a proxy for detailed spectral fitting. However, in this regime classical methods of error propagation fail, and the estimates of spectral hardness become unreliable. Here we develop a rigorous statistical treatment of hardness ratios that properly deals with detected photons as independent Poisson random variables and correctly deals with the non-Gaussian nature of the error propagation. The method is Bayesian in nature, and thus can be generalized to carry out a multitude of source-population--based analyses. We verify our method with simulation studies, and compare it with the classical method. We apply this method to real world examples, such as the identification of candidate quiescent Low-mass X-ray binaries in globular clusters, and tracking the time evolution of a flare on a low-mass star.


Introduction
The shape of the spectrum of an astronomical source is highly informative as to the physical processes at the source. But often detailed spectral fitting is not feasible due to various constraints, such as the need to analyze a large and homogeneous sample for population studies, or there being insufficient counts to carry out a useful fit, etc. In such cases, a hardness ratio, which requires the measurement of accumulated counts in two or more broad passbands, becomes a useful measure to quantify and characterize the source spectrum. a A hardness ratio is defined as either the ratio of the counts in two bands called the soft and hard bands, or a monotonic function of this ratio. We consider three types of hardness ratios, where S and H are the source counts in the two bands, called the soft and hard passbands. The simple formulae above are modified for the existence of background counts and instrumental effective areas.
a In the context of Chandra data, the passbands often used are called soft (0.3-0.9 keV), medium (0.9-2.5 keV), and hard (2.5-8.0 keV) bands. Sometimes the soft and medium bands are merged into a single band, also called the soft band (e.g., 0.5-2 keV). In all cases, note that the energies refer to the nominal detector PHA or PI channel boundaries and not to the intrinsic photon energies. The choice of energy or PI channel range is flexible and situational.
Spectral colors in optical astronomy, defined by the standard optical filters (e.g., UB-VRIJK, U-B, R-I, etc), are well known and constitute the main observables in astronomy. They have been widely applied to characterize populations of sources and their evolution. The hardness ratio concept was adopted in X-ray astronomy for early X-ray detectors (mainly proportional counter detectors) which had only a limited energy resolution. The first application of the X-ray hardness ratio was presented in the X-ray variability studies with SAS-3 observations by Bradt et al. (1976). Zamorani et al. (1981) investigated the X-ray hardness ratio for a sample of quasars observed with the Einstein X-ray Observatory. They calculated each source intensity in two energy bands (1.2-3 keV and 0.5-1.2 keV) and discussed a possible evolution of the hardness ratio (the ratio of both intensities) with redshift for their X-ray sample of 27 quasars. Similar ratios have been used in surveys of stars (Vaiana et al. 1981), galaxies (e.g., Kim, Fabbiano, & Trinchieri 1992) and X-ray binaries (Tuohy et al. 1978) studied with the Einstein and earlier observatories. In the case of X-ray binaries in particular, they were used to define two different classes of sources (Z-track and atoll sources; Hasinger & van der Klis, 1989) depending on their time evolution on HR v/s HR diagrams. Since then the concept of X-ray hardness ratio has been developed and broadly applied in a number of cases, most recently in the analysis of large data samples from the Chandra X-ray Observatory (Weisskopf et al. 2000) and XMM-Newton (Jansen et al. 2001).
Advanced X-ray missions such as Chandra and XMM-Newton allow for the detection of many very faint sources in deep X-ray surveys (see review by Brandt & Hasinger 2005). For example, in a typical observation of a galaxy, several tens and in some cases even hundreds of sources are detected, most of which have fewer than 50 counts, and many have only a few counts. Similar types of sources are detected in the ChaMP serendipitous survey . They provide the best quality data to date for studying source populations and their evolution. The most interesting sources are the sources with the smallest number of detected counts, because they have never been observed before. Further, the number of faint sources increases with improved sensitivity limits, i.e., there are more faint sources in deeper observations. Because these faint sources have only a few counts, hardness ratios are commonly used to study properties of their X-ray emission and absorption (e.g., Alexander et al. 2001, Brandt et al. 2001a, Brandt et al. 2001b, Giacconi et al. 2001. With long observations the background counts increase, so the background contribution becomes significant and background subtraction fails to work. Still the background contribution must be taken into account to evaluate the source counts. The sensitivity of X-ray telescopes is usually energy dependent, and thus the classical hardness ratio method can only be applied to observations made with the same instrument. Usually the measurements carry both detection errors and background contamination. In a typical hardness ratio calculation the background is subtracted from the data and only net counts are used. This background subtraction is not a good solution, especially in the low counts regime (see van Dyk et al. 2001). In general the Gaussian assumption present in the classical method (see § § 1.1 and 3) is not appropriate for faint sources in the presence of a significant background. Therefore, the classical approach to calculating the hardness ratio and the errors is inappropriate for low counts. Instead, adopting a Poisson distribution as we do here (see below), hardness ratios can be reliably computed for both low and high counts cases. The Bayesian approach allows us to include the information about the background, difference in collecting area, effective areas and exposure times between the source and background regions.
Our Bayesian model-based approach follows a pattern of ever more sophisticated statistical methods that are being developed and applied to solve outstanding quantitative challenges in empirical Astronomy and Astrophysics. Bayesian model-based methods for high-energy high-resolution spectral analysis can be found for example in Kashyap & Drake (1998), van Dyk et al. (2001), van Dyk and Hans (2002), Protassov et al. (2002), van Dyk and Kang (2004), Gillessen & Harney (2005), and Park et al. (2006, in preparation). More generally, the monographs on Statistical Challenges in Modern Astronomy (Feigelson and Babu, 1992, 2003, Babu and Feigelson, 1997 and the special issue of Statistical Science devoted to Astronomy (May 2004) illustrate the wide breadth of statistical methods applied to an equal diversity of problems in astronomy and astrophysics.
Here, we discuss the classical method and present the fully Bayesian method for calculating the hardness ratio for counts data. b

The Classical Method
A conventional hardness ratio is calculated as set out in Equations 1, where S and H are the "soft" and "hard" counts accumulated in two non-overlapping passbands. In general, an observation will be contaminated by background counts, and this must be accounted for in the estimate of the hardness ratios. The background is usually estimated from an annular region surrounding the source of interest, or from a suitable representative region on the detector that is reliably devoid of sources. The difference in the exposure time and aperture area of source and background observations are summarized by a known constant r for which the expected background counts in the source exposure area are adjusted. With the background counts in the soft band (B S ) and the hard band (B H ) collected in an area of r times the source region, the hardness ratio is generalized to The adjusted counts in the background exposure area are directly subtracted from those in the source exposure area. The above equations can be further modified to account for variations of the detector effective areas by including them in the constant r, in which case the constants for the two bands will be different. c Errors are propagated assuming a Gaussian regime, i.e., where σ X S , σ X H , σ B S , and σ B H are typically approximated with the Gehrels prescription (Gehrels 1986) σ X ≈ √ X + 0.75 + 1 , where X are the observed counts, and the deviation from √ X is due to identifying the 68% Gaussian (1σ) deviation with the 16 th − 84 th percentile range of the Poisson distribution. In addition to the approximation of a Poisson with a faux Gaussian distribution, classical error-propagation also fails on a fundamental level to properly describe the variance in the hardness ratios. Because R is strictly positive, its probability distribution is skewed and its width cannot be described simply by σ R . The fractional difference HR is better behaved, but the range limits of [−1, +1] cause the coverage rates (see §3) to be incorrect. The color C is the best behaved, but the error range will be asymmetrical when the errors are large. c Gross differences in detector sensitivity, e.g., between different observations, or for sources at different offaxis angles, can be accounted for with additional multiplicative factors modifying the background-subtracted counts. This is however not statistically meaningful and can lead to unreliable estimates and error bars in the small counts regime. Instead, we apply such scaling factors directly to the source intensities (see Equation 5).
Moreover, it has been demonstrated ) that using background subtracted source counts (which is the case in the method outlined above) gives a biased estimate of the source intensity, especially for very faint sources. Furthermore, although we can account for gross differences in the detector sensitivity between observations, this is not done in a statistically correct manner, but instead simply by rescaling the background-subtracted counts in each band using additional pre-computed correction factors. Finally, the classical hardness ratio method cannot be applied when a source is not detected in one of the two bands (which is quite common since CCD detectors are more sensitive in the soft band).
An alternative method based on the median (or some quantile) of the distribution of the energies of the detected counts in a source was suggested by Hong et al. (2004). This method can be very powerful for classifying faint sources, but is not suited for comparisons of spectral properties between observations with different instruments. In this method, because the spectral parameter grids are highly non-linear and multi-valued in quantile-quantile space, it is also necessary that the spectral shape of the source being analyzed be known in order to interpret the numerical values.
We have therefore developed fully Bayesian approaches to appropriately compute hardness ratios and their errors, by properly accounting for the Poissonian nature of the observations. The value of our method is heightened in several astrophysically meaningful cases. In §2, we describe the model structure we adopt. In §3, we carry out simulations to compare our method with the classical method. In §4, we outline various applications of our method; we demonstrate its use in identifying candidate quiescent Low-mass X-ray binary sources in a globular cluster and in studying the evolution of a stellar X-ray flare. In §5, we discuss some nuances of our method such as the adopted prior distributions, as well as the advantages and limitations. A summary of the method and its advantages is presented in §6. A detailed mathematical derivation of the posterior probability distributions of the hardness ratios and the computational methods we use to calculate them are described in Appendix A, and a detailed comparison of the effects of priors is described in Appendix B.

The Poisson Case
Hardness ratios are often used to describe faint sources with very few counts; it is not uncommon for either or both of the soft and hard counts to be small enough that the source is undetected in that band. The classical method ( §1.1) generally fails to produce reliable estimates (or upper and lower limits) in such cases, and indeed, it is sometimes not even possible to apply it, such as when the background is large and simple subtraction results in a negative estimate for the source counts. Thus, instead of the Gaussian assumptions on the detected counts, we directly model them as an inhomogeneous Poisson process.
Since the observed counts are the sum of counts due to the source and the background, and because the sum of two independent Poisson random variables is a Poisson random variable with the sum of two Poisson intensities, we may write d S ∼ Poisson(e S · (λ S + ξ S )) and H ∼ Poisson(e H · (λ H + ξ H )), where S and H are independent data points, λ S and λ H are the expected soft and hard source counts intensities, ξ S and ξ H are the expected soft and hard background counts intensities, and e S and e H are correction factors that take into account variations in effective areas, exposure times, and other instrumental effects. e The observed background counts are also modeled as independent Poisson random variables, where ξ is scaled by a known correction factor r, to account for differences in source and background areas and sensitivities. We implicitly assume throughout that the observed data are discrete photon events, or counts. , then e S and e H are dimensionless, and usually taken to be unity. This allows us to compute hardness ratios of fluxes and count rates as well as counts. We place no restriction on the units of the model intensities (other than that λ and ξ have the same units); they are entirely determined by the problem at hand. Thus, our analysis is valid for any instrumental configuration or combination, and the only requirement is that the observed data be described by a Poisson distribution.
Given the expected source counts intensities (λ S and λ H ), it is legitimate for the hardness ratio to be rewritten as d Note that throughout this paper, Greek letters indicate unobserved quantities and Roman letters denote observed data values (except for R, C, and HR, which are used to represent the hardness ratio model parameters in the Poisson case). e Note that in general the effective area changes within an adopted passband, so the chosen values of the correction factor are averages that are calculated after making some prior assumptions about the underlying source spectrum. Their values may change depending on the analysis results.
Color, C = log 10 λ S λ H , These quantities are characteristics of the spectrum in question, rather than of the observed data, like the quantities in Equation 2. As such, the rewritten hardness ratio is of more direct scientific interest. Note also that while these quantities are defined entirely in terms of expected source intensities in the different passbands, the background is implicitly taken into account via direct modeling, as in Equations 5 and 6.

Bayesian Approach
As a new method of evaluating hardness ratios, we adopt a Bayesian approach. Bayesian analysis is explained in detail in numerous articles in the astronomical literature (e.g., Gregory & Loredo 1992, Kashyap & Drake 1998, Protassov et al. 2002, and here we describe it only in rudimentary terms. In this approach, the prior knowledge for a parameter (the prior probability distribution, e.g., p(λ S , ξ S )) is combined with the information in the observed data (the likelihood, e.g., p(S, B S |λ S , ξ S )) to derive the probability distribution of the model parameters (the posterior distribution, e.g., p(λ S , ξ S |S, B S )) via Bayes' theorem, f i.e., and similarly for p(λ H , ξ H |H, B H ). Because the posterior distribution in Equation 8 needs only be computed up to a normalizing constant, we often deal with an unnormalized posterior density, i.e., p(λ S , ξ S |S, B S ) ∝ p(λ S , ξ S )p(S, B S |λ S , ξ S ). Bayesian statistical inferences for a parameter are based on the posterior distribution, so that we consider its various summary statistics; see §2.2.1 for details.
Because the underlying likelihood functions are Poisson distributions, we assign so-called conjugate γ-prior distributions for both source and background intensities ); a conjugate prior distribution ensures that a posterior distribution follows the same parametric form as the prior distribution, e.g., a γ-distribution is a conjugate prior for a f The notation p(A|B) is to be read as "the probability that A is true given that B is true." Thus, except for prior distributions, all the probability distributions that we use are conditional distributions.
Poisson likelihood. g That is, we assign independent γ-prior distributions for the source and background counts λ and ξ, where the values of ψ are calibrated according to our prior knowledge, or lack thereof, about the parameters; see §5.1 for our discussion on choosing a prior distribution.
We describe the application of Bayes' Theorem and the derivation of the posterior probability distribution functions for the interesting parameters in Appendix A. Here, we note that since λ S and λ H are independent of each other, their joint posterior probability distribution is the product of the individual marginal posterior distributions, h written as (see Equation A8): Then the posterior distribution of each type of hardness ratio can be computed by transforming λ S and λ H to the appropriate variables and marginalizing over the resulting nuisance variable, as follows: with ω = λ S + λ H . The integrals can be computed using either Monte Carlo integration or Gaussian quadrature (throughout this paper, we use the term "quadrature" to refer g Formally, these are semi-conjugate prior distributions in that they are conjugate to a particular conditional distribution. For example, the γ-prior distribution on λ S is conjugate to the Poisson distribution of the source counts among the observed soft counts, denoted η S in Appendix A1. The γ(α, β) distribution is a continuous distribution on the positive real line with probability density function and has a mean of α/β, and a variance of α/β 2 for α, β > 0.
h The marginal posterior distribution of interest is computed by integrating out the so-called nuisance parameters out of their joint posterior distributions. For example, p(λ S |S, exclusively to the latter). We use both methods of integration because neither has the clear advantage over the other in our case; see §5.4 and Appendix A.

Finding Point Estimates and Error Bars
Bayesian statistical inferences are made in terms of probability statements, and thus we consider various summary statistics for a posterior distribution. Here we briefly describe the standard summaries of location and methods of computing error bars; a fuller description along with detailed examples can be found in Park et al. (2006, in preparation). Commonly used summaries of location are the mean, median, and mode(s) of the distribution: The posterior mean is the posterior expectation of the parameter, the posterior median is the point that divides the posterior probability evenly such that 50% of the probability resides below and above its value, and the posterior mode is the most likely value of the parameter given the data.
When a posterior simulation of size N is obtained with a Monte Carlo sampler (see §A.1), the posterior mean is the simple average of the sample draws and the posterior median is the [N/2] th draw after sorting the draws in increasing order. To compute the posterior mode, we repeatedly bisect the Monte Carlo draws and choose the half with more draws until the range of the chosen half is sufficiently narrow. The midpoint of the resulting narrow bin approximates the posterior mode.
With quadrature (see §A.2), we can obtain equally spaced abscissas and the corresponding posterior probabilities. Thus, the posterior mean is computed as the sum of the product of an abscissa with their probabilities (i.e., the dot product of the vector of abscissa with the corresponding probability vector). The posterior median is computed by summing the probabilities corresponding to the ordered abscissa one-by-one until a cumulative probability of 50% is reached. The posterior mode is simply the point among the abscissa with the largest probability.
Unlike with point estimates above, there is no unique or preferred way to summarize the variation of a parameter. Any interval that encompasses a suitable fraction α of the area under the probability distribution qualifies as an estimate of the variation. Of these, two provide useful measures of the uncertainty: the equal-tail posterior interval, which is the central interval that corresponds to the range of values above and below which lies a fraction of exactly (α/2) of the posterior probability, is a good measure of the width of the posterior distribution; and the highest posterior density (HPD) interval, which is the range of values that contain a fraction (1 − α) of the posterior probability, and within which the probability density is never lower than that outside the interval. The HPD-interval always contains the mode, and thus serves as an error bar on it.
For a symmetric, unimodal posterior distribution, these two posterior intervals are identical. The equal-tail interval is invariant to one-to-one transformations and is usually easier to compute. However, the HPD-interval always guarantees the interval with the smallest length among the intervals with the same posterior probability. Neither of these intervals is necessarily symmetric around the point estimate, i.e., the upper and lower bounds may be of different lengths. For consistency, we refer to such intervals as posterior intervals; others also refer to them as confidence intervals i or credible intervals.
For a Monte Carlo simulation of size N, we compute either the equal-tail interval or an interval that approximates the HPD interval. (Unless otherwise stated, here we always quote the equal-tail interval for the Monte Carlo method and the HPD-interval for the quadrature.) The equal-tail posterior interval in this case is computed by choosing the [(α/2)N] th and the [(1 − α/2)N] th draws as the boundaries. An approximate HPD-interval is derived by comparing all intervals that consist of the [X] th and the [X + (1 − α)N] th draws and choosing that X which gives the shortest length among them.
When the posterior density is computed by the quadrature, we split parameter space into a number of bins and evaluate the posterior probability at the midpoint of each bin. In this case, a 100(1 − α)% HPD-interval can be computed by beginning with the bin with the largest posterior probability and adding additional bins down to a given value of the probability density until the resulting region contains at least a fraction (1 − α) of the posterior probability.

Verification: Simulation Study
In order to compare the classical method with our Bayesian method, we carried out a simulation study to calculate coverage rates of the classical and posterior intervals. Given pre-specified values of the parameters, source and background counts were generated and then used to construct 95% classical and posterior intervals of each hardness ratio using the methods discussed in this article. From the simulated data we calculated the proportion of the computed intervals that contain the true value of the corresponding hardness ratio. (This is the coverage rate of the classical and probability intervals.) In the ideal case, 95% of the i Technically, confidence intervals are defined in terms of how frequently they contain the "true" value of the estimand, rather than in terms of a posterior probability distribution. simulated intervals would contain the true value of each hardness ratio. Besides the coverage rate, the average length of the intervals and the mean square error of point estimates were also computed and compared. The mean square error of the point estimateθ of θ is defined as the sum of the variance and squared bias for an estimator, i.e., method that constructs shorter intervals with the same coverage rate and produces a point estimate with a lower mean square error is generally preferred. The entire simulation was repeated with different magnitudes of the source intensities, λ S and λ H . Intrinsically, we are interested in the following two prototypical cases: CASE I : hardness ratios for high counts, i.e., CASE II : hardness ratios for low counts, i.e., In both cases we adopt a background-area-to-source-area ratio of r = 100 (see Equation 6), i.e., we take the observed counts in the background region to be 10. Note that these are written with reference to the counts observed in the "source region", i.e., the units of λ S , λ H , ξ S , ξ H are all [ct (source area) −1 ], and that we have set e S = e H = 1 here. The actual extent of the source area is irrelevant to this calculation.
This simulation study illustrates two typical cases, i.e., high counts and low counts sources: CASE I represents high counts sources for which Poisson assumptions tend to agree with Gaussian assumptions; CASE II represents low counts sources where the Gaussian assumptions are inappropriate. In CASE I, the posterior distributions of the hardness ratios agree with the corresponding Gaussian approximation of the classical method, but in CASE II, the Gaussian assumptions made in the classical method fail. This is illustrated in Figure 1, where we compare the two methods for a specific and particular simulation: we assume S = H = 30 in CASE I and S = H = 3 in CASE II, compute the resulting posterior distributions of hardness ratios using the Bayesian method, and compare it to a Gaussian distribution with mean and standard deviation equal to the classical estimates. The right panels in Figure 1 show that there is a clear discrepancy between the two methods. . The solid lines represent the posterior distributions of the hardness ratios and the dashed lines a Gaussian distribution with mean and standard deviation equal to the classical estimates. A prior distribution that is flat on the real line (φ = 1, see §5.1) is adopted for the Bayesian calculations. Note that as the number of counts increase, the two distributions approach each other. At low counts the two distributions differ radically, with the Classical distributions exhibiting generally undesirable properties (broad, and extending into unphysical regimes). In order to compare the statistical properties of estimates based on the Bayesian and classical methods, we simulate test data under both CASE I and CASE II. Table 1 presents the results of 200 simulated sources for each of the two cases. The data are used to compute point estimates and 95% error bars by using the classical method, Monte Carlo sampler, and numerical integration by quadrature. The Bayesian methods use non-informative prior distributions on λ and ξ; in particular, λ ∼ γ(1, 0) and ξ ∼ γ(0.5, 0), see §5.1. The results of CASE I indicate that all three methods are comparable although the quadrature yields slightly shorter intervals and point estimates with a smaller mean square error. Figure 2 illustrates that the resulting intervals from the methods are almost identical. On the other hand, the results of CASE II appear in Figure 3 and confirm that the Bayesian methods outperform the classical method because the mean length of the classical intervals are wider than the posterior intervals, and the classical estimates exhibit much larger mean square error. For example, the hardness ratio HR by definition is limited to the range [−1, +1], so that the maximum length of the intervals should be two. However, the mean length of simulated  Equation 15b). In this, the low counts case, the the Bayesian methods dramatically outperform the classical method.
classical intervals is 3.7, which is clearly larger than needed. Since quadrature tends to yield estimates with smaller mean square errors and shorter posterior intervals while maintaining high coverage rates, it is generally preferred to the Monte Carlo method. However, this comes at a price: because the summation inside the posterior density is over the detected source counts, quadrature is computationally more expensive as the source counts increase; on the other hand, the Monte Carlo method is much faster regardless of the number of source counts. Thus, we recommend using quadrature with relatively low counts data (e.g., S < 20 or H < 20) and the Monte Carlo method with higher counts.
Because the Bayesian inference for each hardness ratio is based on its posterior distribution, we consider these two point estimates, the posterior mode and the posterior mean. Table 1 presents a comparison of the two proposed estimates in terms of their mean square errors. In the case of HR, the posterior mean seems to be more robust to the shape of a posterior distribution than a posterior mode because both left and right limits of HR are bounded. For R and C, on the other hand, the posterior mode appears to perform better. Thus, the posterior mode is preferable for both R and C and the posterior mean for HR, as shown in Table 1.

Quiescent Low-Mass X-ray Binaries
Low-mass X-ray binaries (LMXBs) are binary systems composed of a neutron star or black-hole and a low mass (< 1 − 2 M ⊙ ) donor. These binaries are generally transient sources and spend most of their time in quiescence (qLMXB's). However, even when not active they emit X-rays due to thermal emission from the surface of the neutron star, and sometimes due to an additional hard component associated with residual low-level accretion (e.g., Campana et al. 1998, Brown, Bildsten, & Rutledge 1998. Identification and study of these quiescent systems is very important since they provide a direct census of the overall LMXB population (and therefore the population of neutron stars), and also constrain their duty cycles. Although only the identification of their optical counterparts can unambiguously classify an X-ray source as a qLMXB, soft X-ray spectra or hardness ratios consistent with those expected for thermal emission by neutron stars are very useful for identifying candidate qLMXB's (Heinke et al. 2003).
We apply the methods described in §2 on a long (40 ks) Chandra ACIS-S observation of the globular cluster Terzan-5 (ObsID 3798; PI: R.Wijnands). A total of 49 sources were detected within the half-mass radius of the cluster (for details on the observation, reduction, and analysis, see Heinke et al. 2006, in preparation). We performed the source photometry in the 0.5-2.0 keV (soft) and 2.0-6.0 keV (hard) bands. In Figure 4 we present a color-luminosity diagram of all the detected sources (note that the color here is scaled by a factor 2.5 to be similar to the definition of an optical colors such as B − V ). The posterior modes of the color for all detected sources, computed here using the Jeffrey's prior (φ = 1/2), along with the 95% highest posterior density intervals, are shown. Those sources which are not amenable to analysis via the classical method (weak sources with background-subtracted source counts in at least one band being < 5) are marked with open circles, while the remainder are marked with filled circles. Thermally emitting neutron stars are expected to have X-ray colors in the range [0.5, 2.5] (indicated by dashed lines in the figure). The softer limit (dashed curve to the right; color ≈ 2) corresponds to a purely thermal emission spectrum, and the harder limit (dashed curve to the left; color ≈ 0.5) corresponds to a 20% contribution from a power-law of index 1.5. From this figure it is clear that the Bayesian method, apart from providing more accurate confidence intervals for the source colors, allows us to estimate hardness ratios for the sources which were undetected in one of the two bands. We find 15 sources whose error bars overlap this region and are thus candidate qLMXB's. Note that 8 of these candidates could not have been so classified using the classical method because there are too few source counts. Conversely, many sources that lie at the edge of the region of interest would have been mistakenly classified as candidate qLMXB's with the classical method if the overly large error bars it produces had been used (see discussion of coverage rates in Table 1). Fig. 4.-The X-ray "color-magnitude" diagram for globular cluster Terzan-5. For each of the detected sources, a form of color C = log 10 (λ S /λ H ), scaled by a factor 2.5 for similarity with optical colors, is plotted along the abscissa, and the log 10 (Luminosity [ergs s −1 ]) in the 0.5-6 keV band, along the ordinate. The modes of the posterior probability functions (computed using quadrature) for all sources are shown, along with the 95% confidence intervals (the error bars on the luminosity are not plotted). The sources which would result in unreliable estimates if the color were computed with the classical method (background subtracted source counts < 5 in at least one band) are marked with open circles (red), and the remainder are marked with filled circles. The dashed (blue) lines represent the color-luminosity tracks expected for neutron stars with black-body spectra of different temperatures (the open triangles correspond to temperatures of log 10 (T [K]) = 6.3, 6.2, 6.1 and 6.0 from top to bottom), and varying contribution from a power-law component (slope Γ = 1.5; pure blackbody for the rightmost curve, 10% and 20% contribution from a power-law component in the 0.5-6 keV band luminosity for the middle and left curves respectively). For more details on these spectral models see Heinke et al. (2006).

Evolution of a Flare
In the previous section ( §4.1), we demonstrated the usefulness of using the full Bayesian machinery in classifying weak sources in comparison to the classical method. Here, we shall go further and show that new avenues of analysis and interpretation are opened up when we take advantage of the full probabilistic description available to us. Such analysis methods are difficult to implement correctly and generally cannot be justified under the classical method. As a specific illustrative example, we consider the time evolution of a large X-ray flare on a low-mass dwarf, Ross 154 (Wargelin et al. 2006). Stellar coronae consist of magnetically confined hot plasma (T ∼ 1 − 10 MK) in the outer atmospheres of late-type stars, and emit optically thin thermal emission composed of Bremsstrahlung continuum and atomic line emission components. Detailed analysis of these line-rich spectra yield valuable clues regarding the temperature and emission structure on the star, its composition, and the processes that drive coronal heating (see e.g., Rosner, Golub, & Vaiana, 1985). Of special importance are flares, which result from impulsive energy input into the corona by the reconnection of highly strained magnetic flux tubes, and provide additional information on the dynamics of the corona. Analysis of flare spectra during decay yields constraints on the physical site and size of the flare loops and on the heating process (see, e.g., van den Oord & Mewe 1989, Pallavicini, Tagliaferri, & Stella 1990, Serio et al. 1991, Schmitt & Favata 1999. A useful technique developed explicitly to analyze hydrodynamically evolving loops is to track their evolution in density-temperature space (Reale, Serio, & Peres 1993, Reale et al. 1997, Reale et al. 2004). This makes possible the modeling of flares to derive physically meaningful parameters such as the length of the loop, the heating function, etc. Unfortunately, a comprehensive analysis along these lines requires that complex model spectra be fit to the data at each stage of the flare evolution, and because fitting requires a large number of counts (generally > 500) to produce reliable results, the resolution at which the flare evolution can be studied is strongly limited. This is especially so during the rising phase of flares, where the physically most interesting processes are of short duration and the observed counts are also small. However, hardness ratios can be constructed at relatively finer time binning resolution and can be used as a proxy for the spectral fit parameters. In the following, we apply the Bayesian method described above (see §2) to Chandra data of a stellar X-ray flare, and demonstrate its usefulness.
A large flare was observed on Ross 154, a nearby (2.93 pc) active M dwarf (M4.5Ve), during a 57 ks observation with the Chandra/ACIS-S detector (ObsID 2561; PI: B.Wargelin). During the flare, the counting rate increased by over a factor of 100. Because of strong pileup in the core of the PSF, only those data within an annular region surrounding the center can be used. (Details of the reduction and analysis are given in Wargelin et al. 2006.) The light-curve of the source during the flare is shown in Figure 5, along with light-curves for smaller passbands (soft, 0.25-1 keV; medium, 1-2 keV; hard, 2-8 keV). Note that the hard band light-curve peaks at an earlier time than the softer bands, an impression that is borne out by the evolution of the color, C (lower panel of Figure 5). The same data can also be represented in a color-intensity diagram, as in Figure 6, which shows scatter plots of the color and intensity at each time bin along with their computed errors (connected by a thin line to indicate the time sequence). The advantage of doing so is that color is inversely correlated with plasma temperature, and the intensity approximately tracks the plasma density. j A hysteresis-like evolution is discernible in the color-intensity track, even at a time resolution as small as 40 s. However, the error bars on the data points are large enough to obscure this structure at small time bin sizes, whereas at larger bin sizes, the number of data bins are smaller and it becomes difficult to discern the totality of the evolutionary track, even though the pattern is persistent as bin size and starting times are varied. Mainly, the choices of the bin size and the phase (i.e., the starting point) is completely arbitrary and the choice of no single value can be properly justified. However, Bayesian analysis provides a simple method to work around this difficulty, by simply considering the bin size as a nuisance parameter of the problem and to marginalize over it.
Thus, in order to ameliorate the effects of choosing a single time bin size and a starting time, we have carried out a series of calculations using the Monte Carlo method: we obtain different light-curves by varying the phase of the binning (i.e, the starting time), and for each light-curve thus obtained, we obtain 2000 samples of (λ S , λ H ) for each data point. We hold the time bin size fixed during this process. After cycling through all the phases, the samples of (λ S , λ H ) are all combined to form a single image (the shaded images in Figure 6). This procedure, also called "cycle-spinning" (Esch 2003), allows us to discern the pattern of the evolutionary track clearly. The longer the source spends at any part of the color-intensity diagram, the darker the shade in the images, and the statistical error in the calculation is reflected by the gradient of the shading. Further, all the cycle-spun images thus constructed, for bin sizes ranging from 40 s to 400 s, were then averaged to produce a single coherent image, thus effectively marginalizing over the bin sizes; such an image includes all the errors associated with the analysis (Figure 7).
Note that this type of analysis is not possible using the classical method. For instance, in the Bayesian case, the samples are drawn from the posterior probability distributions of S and H. No such sampling technique is available in the classical case, k and naively using j X-ray luminosity, L X ∝ EM · P (T ), where EM = n 2 e V is the emission measure of coronal plasma with an electron density n e occupying a volume V , and P (T ) is the power emitted by the plasma at temperature T . P (T ) is a weak function of T over the temperature range 6-20 MK expected in this source, and assuming that the volume of emission remains constant, the observed intensity is proportional to the square of the electron density. For a more detailed modeling of the L X − T conversion, see Wargelin et al. (2006). k New counts realizations can be obtained by bootstrapping from the observed counts, but such sampling will be quantized and will render the shaded images unusable at lower counts intensities. the Gaussian approximations as in Equation 4 will overestimate the errors and wash out the signal (see §3). Furthermore, the concept of marginalizing over the nuisance parameter of time bin sizes has no analog in the classical system, and cannot be justified; and cyclespinning leads to correlations between the estimates at different bins since to the statistical independence of the data will be lost.
A sharp initial hardening of the spectrum, coincident with the onset of a small optical flare (Wargelin et al. 2006) is visible at the beginning of the time sequence (points 1-3; see Figures 5,6). A point worth noting is that the Bayesian analysis prevents overinterpretation of the results: at time point 3, the standard analysis indicates that the spectrum is at its hardest, suggesting either a non-thermal origin to the emission or significant transient absorption, possibly due to a coronal mass ejection, but the shaded image demonstrates that the statistical errors are sufficiently large that a more mundane low-density thermal expansion is quite plausible. This stage is followed by a rapid softening (points 3-8), interpreted as the thermalization of the non-thermal hard X-ray flare. Then the star, which lies in the lower central portion of the color-intensity diagram before the flare, moves to the upper left as the flare is set off (the spectrum hardens as intensity increases; points 8-11), turns right at flare peak (softens as the deposited energy cascades to lower temperatures; points 11-12), and then moves down the right flank back to its original state (decays in temperature and intensity; points 12-19). The physical consequences of this analysis are discussed in detail by Wargelin et al. Fig. 7.-Ross 154 color-intensity track. As in Figure 6, but for an image obtained by averaging cycle-spun images at various bin sizes ranging from 40 s to 400 s. The persistence of the evolutionary track of the flare on Ross 154 is quite clear in this visualization. The points corresponding to each bin of the light-curve obtained at a bin size of 150 s are superimposed on the image, with a white solid line connecting them in temporal sequence. A clear pattern of the flaring plasma evolution can be seen, first rising in temperature and intensity, softening at the peak, and then decaying in temperature and intensity.

Informative and Non-informative Prior Distributions
A major component of Bayesian analysis is the specification of a prior distribution, the probability distribution assigned a priori to the parameters defining the model. This is often deemed controversial because of the apparent subjectivity in the assignment: if different people assign different priors and obtain different results, how is it possible to evaluate them? The answer to this conundrum is to realize that all problems always include a bias brought to the analysis by the researcher, either in the choice of statistic, or even in the choice of analysis method, and the specification of a prior distribution codifies and makes explicit this bias. Furthermore, if indeed prior information about the value or the range of the parameter is available, the prior allows a simple and principled method by which this information can be included in the analysis. If no prior information is available, a so-called non-informative prior distribution must be chosen in order to not bias the result.
In our case, if there is a strong belief as to the value or range of the hardness ratio, we can incorporate this via an informative prior distribution, encoded as a γ-distribution (see van Dyk et al. 2001). In addition, the analysis of previously acquired data will produce a posterior probability distribution that can be used as an informative prior on future observations of the same source.
When no prior information is available, we normally use a non-informative prior distribution. Since a Poisson intensity is necessarily positive, three types of non-informative prior distributions immediately present themselves: when X|θ ∼ Poisson(θ), 1. a non-informative prior distribution on the original scale, 2. a Jeffrey's non-informative prior distribution, 3. a non-informative prior distribution under a log transformation, where I θ = E[−∂ 2 log p(X|θ)/∂θ 2 |θ] is the expected Fisher information (Casella & Berger 2002). The first non-informative prior distribution (Equation 16a) is flat from 0 to ∞; the second (Equation 16b) is proportional to the square root of the Fisher information; and the third (Equation 16c) is flat on the whole real line under a log transformation. The functional forms of these prior distributions can be generalized into p(θ) ∝ θ φ−1 , i.e., θ ∼ γ(φ, 0) where φ is an index that varies from 0 to 1: the first corresponds to φ = 1, the second to φ = 1/2; and the third to φ = 0 + (where this notation indicates that φ > 0, but arbitrarily close to zero, for reasons described below). We note that these non-informative prior distributions are improper, i.e., are not integrable. If an improper prior distribution causes a posterior distribution to be improper as well, then no inferences can be made using such non-integrable distributions. In our case, as long as φ is strictly positive, a posterior distribution remains proper, i.e., integrable. Hence, in the third case, we adopt values of φ that are strictly positive but close in value to 0, e.g., φ = 0.01 or φ = 0.1.
Note that while these prior distributions are non-informative, in the sense that in most cases they do not affect the calculated values of the hardness ratios (see Appendix B), they do codify some specific assumptions regarding the range of values that it is believed possible. For instance, if a large dynamic range is expected in the observations, a flat prior distribution in the log scale is more appropriate than a flat prior distribution in the original scale. The choice of the prior distribution is dictated by the analysis and must be reported along with the results.

R versus C versus HR
At low counts, the posterior distribution of the counts ratio, R, tends to be skewed because of the Poissonian nature of data; R only takes positive values. The color, C = log 10 R, is a log transformation of R, which makes the skewed distribution more symmetric. The fractional difference hardness ratio, HR = (1−R)/(1+R), is a monotonically decreasing transformation of R, such that HR → +1 as R → 0 (i.e., a source gets harder) and HR → −1 as R → ∞ (i.e., a source gets softer). The monotonic transformation results in a bounded range of [−1, +1] for HR, thereby reducing the asymmetry of the skewed posterior distribution. R and HR are bounded on one side or two sides, while C is unbounded due to the log transformation.
The posterior distribution of any hardness ratio becomes more symmetric as both soft and hard source intensities increase. Regardless of the size of the intensities, however, the color has the most symmetric posterior distribution among the popular definitions of a hardness ratio. Figure 8 illustrates the effect of the magnitude of the source intensities on the symmetry of the posterior distribution of each hardness ratio; the posterior distribution of C is confirmed to have the most symmetric posterior distribution. In the figure, we fix R = 2 and the soft and hard intensities are determined by beginning with λ S = 2 and λ H = 1 (in units of counts (source area) −1 ) and increasing the intensities by a factor of 5 in each subsequent column. We assume no background contamination in each simulation. Figure 8 also illustrates the effect of the use of different indices for the non-informative prior distribution on the resulting posterior distribution of each hardness ratio. In the case of R, the posterior distribution becomes diffuse and more skewed toward 0 as the value of φ decreases; as expected, this is also the case for C. And in the case of HR, the posterior distribution has more mass near the boundaries (i.e., ±1) as the value of φ decreases, thereby exhibiting a U-shaped density in the case of very faint sources. As the intensities increase, however, the choice of a prior distribution has less effect on the posterior distribution of any of the hardness ratios.
Finally, as pointed out in §3, for R and C, the mode of the posterior probability distribution is a robust estimator, whereas for HR, the mean of the distribution is a better choice. (Notice that the posterior mode of HR in the low count scenario is −1 when φ = 0 + in Figure 8.)

Advantages
A significant improvement that our method provides over the classical way of computing hardness ratios is that we use the correct Poisson distribution throughout and do not make any limiting (high counts) Gaussian assumptions. Thus, while the classical method works well only with high counts data and fails to give reliable results with low counts, our method is valid in all regimes. Indeed, because the observed counts are non-negative integers, it is not appropriate to model low counts with Gaussian distributions which are defined on all the real numbers.
Because our methods are based on Poisson assumptions in a fully model based statistical approach, we need not rely on plug-in estimates of the hardness ratio. Instead we can evaluate the full posterior probability distribution of the hardness ratio, which provides reliable estimates and correct confidence limits even when either or both soft and hard counts are very low. In particular, our method is not limited to "detectable" counts, and can properly calculate upper and lower bounds in cases where the source is detected in only one passband. For high counts, our method gives more precise error bars for the hardness ratios than the classical method (as defined by the coverage; see Table 1 and Appendix A), although both methods yield similar results.
Further, a priori information can be embedded into our model and can be updated by data, thereby producing more accurate results as we collect more data.

Limitations
Unlike the quantile-width method (Hong et al. 2004), the calculation of hardness ratios is limited by necessity to pre-selected non-overlapping passbands. In general, we must use different band definitions in order to achieve the maximum discriminating power for soft and hard sources. A priori, with no knowledge of the character of the source populations, there is no way to ensure that the chosen passbands are the most efficient for splitting the spectrum and obtaining maximum leverage to distinguish spectral model parameters. However, this restriction allows for a more uniform analysis and comparison of results across different instruments and telescopes.
Because the Bayesian method does not allow for a simple analytic solution similar to standard error propagation in the Gaussian case, the computational methods used to determine the probability distributions become important. We have implemented a Monte Carlo integration method (the Gibbs sampler, see Appendix A.1) and a method based on numerical integration (quadrature, see Appendix A.2). The Gibbs sampler is based on Monte Carlo simulation and hence causes our estimates to have simulation errors in addition to their true variability; to reduce the levels of simulation errors, the Markov chain must be run with large enough iterations. On the other hand, the quadrature precisely computes the posterior distribution as long as the number of bins is large enough; however, its computation becomes expensive as the counts become large because of the binomial expansion in Equation A5. In general, the Gibbs sampler is much faster than the method based on quadrature, but care must be taken to ensure that the number of iterations is sufficient to determine the posterior mode and the required posterior interval with good precision. We recommend using the Gibbs sampler for high counts and the quadrature for low counts.

Summary
We have developed a method to compute hardness ratios of counts in non-overlapping passbands that is fully Bayesian and properly takes into account 1. the Poissonian nature of the data, 2. background counts, allowing for differences in collecting areas, effective areas, and exposure times between the source and background regions, and 3. both informative and non-informative prior information.
The algorithm is implemented in two forms: a Monte Carlo based Gibbs sampler that samples from the posterior probability distributions of the counts intensities, and another that numerically integrates the joint probability distribution of the counts intensities using quadrature.
We carry out comparison tests between our method and the classical method, and show that the Bayesian method works even in the low counts regime where the classical method fails, and matches the results of the classical method at high counts where a Gaussian approximation of the errors is appropriate.
We apply the method to real world test cases of 1. determining candidate quiescent Low-Mass X-ray Binaries in the Terzan-5 cluster, where we show that it can be applied even in cases where the classical method is entirely inapplicable, such as those where the source is detected in only one passband, and 2. tracking the time evolution of a flare on the dMe star Ross 154, where we demonstrate the flexibility of the method in being able to work around the conflicting limitations imposed by time bin sizes and time resolution.
Additional applications of these methods include spectral classification of the X-ray sources detected in deep surveys and other galaxies, and spectral variability of AGNs and X-ray binaries. The physical properties of the source influence the source intensities λ S and λ H ; given a source spectral model, the two quantities are related. Thus, variations in the hardness ratios reflect changes in the source properties (e.g., sources with different temperatures have different colors). A relation between source properties, physical model parameters and hardness ratio is a natural extension of this work and will be presented in the future.
We thank Yue Wu for assistance with code development, and Dong-Woo Kim and Peter Edmonds for useful discussions. The authors gratefully acknowledge funding for this project partially provided by NSF grants DMS-01-04129, DMS-04-38240, and DMS-04-06085 and by NASA Contracts NAS8-39073 and NAS8-03060 (CXC), and NASA grant NAG5-13056. This work is a product of joint work with the California-Harvard astrostatistics collaboration (CHASC) whose members are TP, VK, AS, DvD, and AZ, and J. Chiang, A. Connors, D.  (Tanner & Wong 1987) accomplish Monte Carlo simulation by constructing a Markov chain with stationary distribution equal to the target posterior distribution. The Markov chain is then run, and, upon convergence delivers the desired Monte Carlo simulation of the posterior distribution. The Gibbs sampler (Geman & Geman 1984) is a special case of MCMC that constructs the Markov chain through alternating conditional simulation of a subset of the model parameters conditioning on the others. A more detailed description of these methods can be found in van Dyk (2003) and the Appendix of van .
In order to construct a Gibbs sampler, we embed the statistical model described in §2.1 into a larger model. Since the observed counts are the sum of the of source counts (η) and the background counts (β), we write S = η S + β S and H = η H + β H . That is, the source counts are modeled as independent Poisson random variables, η S ∼ Poisson(e S · λ S ) and η H ∼ Poisson(e H · λ H ), and the background counts in the source exposure area as independent Poisson random variables, where λ and ξ denote the expected source and background counts in the source region. This implies Equation 5, where the detected photons in the source region are a convolution of the source and background counts, i.e., S = β S + η S and H = β H + η H .
Using the method of data augmentation, we treat the background counts in the source exposure area (β S and β H ) as missing data; the source counts (η S and η H ) are fully determined when β S and β H are known. Intuitively, it is straightforward to estimate the Poisson intensities if detected photons are split into the source and background counts. Based on this setting, the Gibbs sampler produces Monte Carlo draws of λ S and λ H along with stochastically imputing the missing data (β S and β H ). The conditional independence of λ S and λ H implies that a Monte Carlo simulation of each hardness ratio can be obtained by simply transforming that of λ S and λ H . Due to the conditional independence, both soft and hard bands also have exactly the same sampling steps, and thus we illustrate the Gibbs sampler only for the soft band. First, the joint posterior distribution of λ S , ξ S , and β S is given by That is, conditional on the total soft counts (S), the unobserved background counts in the source exposure area (β S ) follows a binomial distribution: Given the current iterates of the parameters, λ S and ξ (t) S , STEP 1 is given by , for t = 1, . . . , T . Next, STEPS 2 and 3 draw the source and background intensities from the γ-distributions. In particular, STEPS 2 and 3 find the next iterates of the intensities using for t = 1, . . . , T . To accomplish one Gibbs iteration, these steps are implemented for the soft and hard bands. After iterating the Gibbs sampler T times, we collect a posterior sample {λ H , t = t 0 + 1, . . . , T } for a sufficiently long burn-in period l t 0 . The analytical calculation to determine burn-in is far from computationally feasible in most situations. However, visual inspection of plots of the Monte Carlo output is commonly used for determining burn-in. More formal tools for determining t 0 , called convergence diagnostics, have been proposed; for a recent review, see Cowles & Carlin (1996). Under a monotone transformation (see Equations 12,13,14) of the posterior samples, (T − t 0 ) Monte Carlo draws for each hardness ratio are obtained, which enables us to find its point estimates and the corresponding error bar; see §2.2.1 for details.
l The term "burn-in" refers to the number of iterations necessary for the Markov chain to converge and begin to sample from the posterior probability distribution.

A.2. Numerical Integration: Gaussian Quadrature
Instead of introducing the missing data, the Bayes' theorem can analytically compute the joint posterior distribution of λ S and λ H , based on which we obtain the posterior distribution of each hardness ratio through quadrature. Because the models for the hard and soft bands are independent and symmetric, we illustrate the computation only for the soft intensity λ S . First, the joint posterior distribution of λ S and ξ S in Equation 8 is written out as Then, the binomial expansion to (λ S + ξ S ) S , i.e., enables us to analytically integrate ξ S out of the joint distribution in Equation A4; and obtain an analytical solution to the marginal posterior distribution of λ S . Therefore, using the binomial expansion, Equation A4 becomes and then ξ S is integrated out of Equation A6, i.e., Here the prior independence of λ S and λ H implies that the joint posterior distribution of these two intensities is decomposed into a product of their marginal posterior distributions, i.e., This assumption of independence between λ S and λ H can be relaxed, and hierarchical modeling for sources to determine clustering properties and spectral parameters can be devised (in preparation).
Using the joint posterior distribution in Equation A8, we compute the analytical solution to the posterior distribution of each hardness ratio as follows: 1. the posterior distribution of R is obtained after integrating λ H out of 2. the posterior distribution of C is obtained after integrating λ H out of (10) To numerically integrate out a parameter in these analytical solutions, we employ quadrature, which precisely evaluates the integral of a real valued function (and has nothing to do with the Gaussian counts statistic assumptions); refer to Wichura (1989) for details of the algorithm. Due to the complex expression, the probability density for each hardness ratio is approximated for each equally-spaced abscissas equally spaced over the finite range of the hardness ratio. Then we summarize the approximate posterior density in order to calculate the inferences, as detailed in §2.2.1.

B. Effect of Non-informative Priors
In order to compare the effect of these non-informative prior distributions, we have carried out simulations for low counts and high counts cases to test the computed coverage Average length of intervals with φ = 0 + intervals with φ = 0 +

Soft band source intensity
Coverage rate of Average length of λ S [counts (source area) −1 ] intervals with φ = 1/2 intervals with φ = 1/2 Coverage rate of Average length of intervals with φ = 1 intervals with φ = 1 against the simulation: Table 3 presents the coverage rate and mean length of the 95% intervals for the color, C. In order to simplify the comparison, we assume that there is no background contamination, and simulate 1000 realizations of source counts for each case.
The tables are laid out as a grid of soft and hard source counts per source area (i.e., λ S versus λ H , see Equation A1). For each (λ S , λ H ) pair, we report the percentage of the simulations that result in 95% posterior intervals that actually contain (λ S , λ H ) along with the mean length of the 95% posterior intervals from the simulations. This calculation is carried out for three choices of the prior index, φ = 0 + , 1/2, 1, corresponding to the top, middle, and bottom elements in each cell of the table.
In these simulations, we theoretically expect that the computed 95% posterior intervals contain the actual value with a probability of 95%. Due to the Monte Carlo simulation errors, we expect most coverage rates to be between 0.93 and 0.97 which are three standard deviations away from 0.95; a standard deviation of the coverage probability for 95% posterior intervals is given by (0.95)(0.05)/1000 = 0.0069 under a binomial model for the Monte Carlo simulation. Table 3 presents the coverage rate and average length of posterior intervals for small and large magnitudes of source intensities. The key to this table is given in Table 2: For each (λ S , λ H ) pair, the posterior intervals are simulated using different prior distribution indices (φ = 0 + , 1/2, 1) and the summary statistics -the coverage rate and the mean lengths of the intervals -are displayed from top (φ = 0 + ) to bottom (φ = 1) within each cell. The same information is shown in graphical form in Figure 10. The use of φ = 0 + tends to yield very wide posterior intervals and under-cover the true X-ray color when the source intensities are low. On the other hand, the other two non-informative prior distribution indices produce much shorter posterior intervals and maintain high coverage rates. With high counts, however, the choice of a non-informative prior distribution does not have a noticeable effect on the statistical properties of the posterior intervals. The same results are summarized in graphical form in Figure 9, where the 95% coverage rates are shown for various cases. The empirical distributions of the coverage rate computed under either low count (λ S = 0.5, 1, 2, 4 or λ H = 0.5, 1, 2, 4) or high count (all of the remaining) scenarios are shown, along with a comparison with the results from the classical method (last column). The distinction between the Bayesian and classical methods is very clear in this figure.
Note that the posterior intervals over-cover relative to our theoretical expectation, but to a lesser degree than the classical intervals. Over-coverage is conservative and preferred to under-coverage, but should be minimal when possible. Considering the low count and high count scenarios, Figure 9 suggests using the Jeffrey's non-informative prior distribution (i.e., φ = 1/2) in general. At low counts, the shape (and consequently the coverage) of the posterior distribution is dominated by our prior assumptions, as encoded in the parameter φ. The coverage rate will also be improved when informative priors are used. Table 3: Coverage of the X-ray Color (C) Using the Bayesian Method with Different (see Table 2 for key) Non-Informative Prior Distribution Indices (φ = 0 + , 1/2, 1).  Fig. 9.-The empirical distributions of coverage rates with different prior distribution indices. The first three columns of the panels are constructed under different types of noninformative prior distributions, φ = 0 + (first), φ = 1/2 (second), φ = 1 (third). The fourth column depicts results from the classical method, and is presented for reference. (In calculating the coverage and the length of the intervals, We exclude simulations where S = 0 or H = 0 since C is undefined in such cases.) The top panels correspond to the low count scenarios λ S = (0.5, 1, 2, 4) or λ H = (0.5, 1, 2, 4) and the bottom panels correspond to the remaining cases. Coverage rate is the fraction of simulated intervals that contain the parameter under which the data were simulated. In this case, we compute the coverage rate based on the 95% posterior intervals for C with different combinations of values for λ S and λ H . The histograms represent the empirical distribution of the coverage rate. The dotted vertical lines represent the expected coverage rate of 95%.  Table 3. The symbols are located on a grid corresponding to the appropriate λ S (absissa) and λ H (ordinate), for φ = 0 + (top), φ = 1/2 (middle), and φ = 1 (bottom). The shading of each symbol represents the coverage, with 100% being lightest and progressively getting darker as the coverage percentage decreases. The sizes of the symbols are scaled as the √ length of interval. Note that for small values of φ, which correspond to a prior expectation of a large dynamic range in the source colors, the intervals are large when the counts are small (i.e., the choice of prior distribution has a large effect), and decrease to values similar to those at larger φ when there are more counts.