Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://hea-www.harvard.edu/AstroStat/BEHR/park_etal_apj.pdf Äàòà èçìåíåíèÿ: Sun Jun 11 15:44:34 2006 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 01:57:21 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: annular

submitted to the Astrophysical Journal

Bayesian Estimation of Hardness Ratios: Mo deling and Computations
Taeyoung Park1 , Vinay L. Kashyap2 , Aneta Siemiginowska2 , David A. van Dyk3 , Andreas Zezas2 , Craig Heinke4 , & Bradford J. Wargelin2 Department of Statistics Harvard University One Oxford Street, Cambridge, MA 02138 tpark@stat.harvard.edu Smithsonian Astrophysical Observatory, 60 Garden Street, Cambridge, MA 02138 kashyap@head.cfa.harvard.edu aneta@head.cfa.harvard.edu azezas@head.cfa.harvard.edu bwargelin@head.cfa.harvard.edu Department of Statistics University of California, 364 ICS Bldg One, Irvine, CA 92697-1250 dvd@ics.uci.edu Northwestern University, 2131 Tech Drive, Evanston, IL 60208 cheinke@northwestern.edu ABSTRACT 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 metho ds of error propagation fail, and the estimates of spectral hardness become unreliable. Here we develop a rigorous statistical treatment
4 3 2 1


­2­ 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 metho d is Bayesian in nature, and thus can be generalized to carry out a multitude of source-population­based analyses. We verify our metho d with simulation studies, and compare it with the classical metho d. We apply this metho d 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. Subject headings: metho ds: statistical ­ stars: flare ­ X-rays: binaries

1.

Introduction

The shape of the spectrum of an astronomical source is highly informative as to the physical pro cesses 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, Simple Ratio, Color, Fractional Difference, R C HR S H log
10

S H (1 )

H -S H +S

where S and H are the source counts in the two bands, called the soft and hard passbands. The simple formulae above are mo dified for the existence of background counts and instrumental effective areas.
In the context of Chandra data, the passbands often used are called soft (0.3-0.9 keV), medium (0.92.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.
a


­3­ Spectral colors in optical astronomy, defined by the standard optical filters (e.g., UBVRIJK, 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 (Kim et al. 2005). 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, Silverman et al. 2005). 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 metho d 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


­4­ counts are used. This background subtraction is not a go o d solution, especially in the low counts regime (see van Dyk et al. 2001). In general the Gaussian assumption present in the classical metho d (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 mo del-based approach follows a pattern of ever more sophisticated statistical metho ds that are being developed and applied to solve outstanding quantitative challenges in empirical Astronomy and Astrophysics. Bayesian mo del-based metho ds 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 Mo dern 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 metho ds applied to an equal diversity of problems in astronomy and astrophysics. Here, we discuss the classical metho d and present the fully Bayesian metho d for calculating the hardness ratio for counts data.b

1.1.

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 (BS ) and the hard band (BH ) collected in an area
A Fortran and C-based program which calculates the hardness ratios following the methods described in this paper is available for download from http://hea-www.harvard.edu/AstroStat/BEHR/
b


­5­ of r times the source region, the hardness ratio is generalized to S - BS /r H - BH /r S- C = log10 H- (H - BH /r HR = (H - BH /r R=

BS /r BH /r ) - (S - BS /r ) ) + (S - BS /r )

(2 )

The adjusted counts in the background exposure area are directly subtracted from those in the source exposure area. The above equations can be further mo dified 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., R = C = S - BS /r H - BH /r 1 ln(10) 2 HR =
2 2 2 2 S + BS /r 2 H + BH /r 2 + (S - BS /r )2 (H - BH /r )2

2 2 2 2 S + BS /r 2 H + BH /r 2 + (S - BS /r )2 (H - BH /r )2

2 2 2 2 (H - BH /r )2(S + BS /r 2 ) + (S - BS /r )2 (H + BH /r 2)

[(H - BH /r ) + (S - BS /r )]2

(3 )

where XS , XH , BS , and BH are typically approximated with the Gehrels prescription (Gehrels 1986) X X + 0.75 + 1 , (4 ) where X are the observed counts, and the deviation from X is due to identifying the 68% Gaussian (1 ) deviation with the 16th - 84th 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.
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).
c


­6­ Moreover, it has been demonstrated (van Dyk et al. 2001) that using background subtracted source counts (which is the case in the metho d 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 metho d 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 metho d 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 metho d can be very powerful for classifying faint sources, but is not suited for comparisons of spectral properties between observations with different instruments. In this metho d, 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 metho d is heightened in several astrophysically meaningful cases. In §2, we describe the mo del structure we adopt. In §3, we carry out simulations to compare our metho d with the classical metho d. In §4, we outline various applications of our metho d; 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 metho d such as the adopted prior distributions, as well as the advantages and limitations. A summary of the metho d and its advantages is presented in §6. A detailed mathematical derivation of the posterior probability distributions of the hardness ratios and the computational metho ds we use to calculate them are described in Appendix A, and a detailed comparison of the effects of priors is described in Appendix B.

2.

Modeling the Hardness Ratios 2.1. 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 metho d (§1.1) generally fails to pro duce reliable estimates (or upper and lower limits) in such cases, and indeed, it is sometimes not even


­7­ 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 mo del them as an inhomogeneous Poisson pro cess. 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 writed S Poisson(eS · (S + S )) and H Poisson(eH · (H + H )), where S and H are independent data points, S and H are the expected soft counts intensities, S and H are the expected soft and hard background co and eS and eH are correction factors that take into account variations in exposure times, and other instrumental effects.e The observed background mo deled as independent Poisson random variables, BS Poisson(r · eS · S ) and BH Poisson(r · eH · H ) (5 )

and hard source unts intensities, effective areas, counts are also (6 )

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. The intensities and may have other units, determined by how the instrument efficiency is supplied. Thus, if eS is given in [ct s cm2 ph-1 ], then S will have units [ph s-1 cm-2 ]. Alternately, eS could describe the exposure time in [s], in which case S will have units of [ct s-1 ]. In the case where S and H themselves have the units [ct], then eS and eH 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 mo del 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 S Simple Ratio, R = , H
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). 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.
e d


­8­ S H H - S Fractional Difference, HR = . H + S Color, C= log
10

, (7 )

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 mo deling, as in Equations 5 and 6.

2.2.

Bayesian Approach

As a new metho d 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, van Dyk et al. 2001, Kashyap et al. 2002, 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 likeliho o d, e.g., p(S, BS |S , S )) to derive the probability distribution of the mo del parameters (the posterior distribution, e.g., p(S , S |S, BS )) via Bayes' theorem,f i.e., p(S , S |S, BS ) = and similarly for p(H , only be computed up to density, i.e., p(S , S |S, parameter are based on statistics; see §2.2.1 for p(S , S )p(S, BS |S , S ) , p(S , S )p(S, BS |S , S )dS dS (8 )

H |H, BH ). Because the posterior distribution in Equation 8 needs a normalizing constant, we often deal with an unnormalized posterior BS ) p(S , S )p(S, BS |S , S ). Bayesian statistical inferences for a the posterior distribution, so that we consider its various summary details. called et al. same for a

Because the underlying likeliho o d functions are Poisson distributions, we assign soconjugate -prior distributions for both source and background intensities (van Dyk 2001); a conjugate prior distribution ensures that a posterior distribution follows the parametric form as the prior distribution, e.g., a -distribution is a conjugate prior
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.


­9­ Poisson likeliho o d.g That is, we assign independent -prior distributions for the source and background counts and , S (S1 , S2 ) and H (H1 , H2 ) , S (S3 , S4 ) and H (H3 , H4 ) , (9 ) (1 0 )

where the values of are calibrated according to our prior knowledge, or lack thereof, about the parameters; see §5.1 for our discussion on cho osing 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 pro duct of the individual marginal posterior distributions,h written as (see Equation A8): p(S , H |S, H, BS , BH ) = p(S |S, BS )p(H |H, BH ) . (1 1 )

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: p(R|S, H, BS , BH ) dR = dR p(C|S, H, BS , BH ) dC = dC d H
H H

p(RH , H |S, H, BS , BH ) ,

(1 2 ) (1 3 )

d H
H

H

10C ln(10) p(10C H , H |S, H, BS , BH ) ,

( 1 - H R ) ( 1 + H R ) p , S, H, BS , BH , 2 2 2 (1 4 ) 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 p(HR|S, H, BS , BH ) dHR = dHR d
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 p(X ) = 1 X ()
-1 - X g

e

,

and has a mean of / , and a variance of / 2 for , > 0. 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, BS ) = p(S , S |S, BS ) d S .
h


­ 10 ­ exclusively to the latter). We use both metho ds of integration because neither has the clear advantage over the other in our case; see §5.4 and Appendix A. 2.2.1. 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 lo cation and metho ds 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 lo cation are the mean, median, and mo de(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 mo de 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 mo de, we repeatedly bisect the Monte Carlo draws and cho ose 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 mo de. 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 pro duct of an abscissa with their probabilities (i.e., the dot pro duct 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 mo de 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 go o d 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


­ 11 ­ density is never lower than that outside the interval. The HPD-interval always contains the mo de, and thus serves as an error bar on it. For a symmetric, unimo dal 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 intervalsi 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 metho d and the HPD-interval for the quadrature.) The equal-tail posterior interval in this case is computed by cho osing 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 cho osing 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.

3.

Verification: Simulation Study

In order to compare the classical metho d with our Bayesian metho d, 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 metho ds 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
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.
i


­ 12 ­ 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., MSE() = E[( - )2 ] = ^ ^ Var() + [E() - ]2 A metho d that constructs shorter intervals with the same coverage rate and pro duces 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: C
ASE

I : hardness ratios for high counts, i.e., S = H = 30, S = H = 0.1; (1 5 a )

C

ASE

I I : hardness ratios for low counts, i.e., S = H = 3, S = H = 0.1. (15b)

In both cases i.e., we take written with S , H , S , H actual extent

we adopt a background-area­to­source-area ratio of r = 100 (see Equation 6), the observed counts in the background region to be 10. Note that these are reference to the counts observed in the "source region", i.e., the units of are all [ct (source area)-1 ], and that we have set eS = eH = 1 here. The 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 metho d, but in CASE II, the Gaussian assumptions made in the classical metho d fail. This is illustrated in Figure 1, where we compare the two metho ds 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 metho d, 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 metho ds.


­ 13 ­

1.5

0.6

High counts
1.0

Low counts

pdf

0.5

pdf
0 1 2 3 4

0.0

0.0 0

0.2

0.4

1

2

3

4

R
1.2

R

3

pdf

pdf

2

1

-1.0

-0.5

0.0

0.5

1.0

0.0

0

0.4

0.8

High counts

Low counts

-1.0

-0.5

0.0

0.5

1.0

C
3.0

C

High counts

0.8

Low counts

pdf

2.0

pdf
0.0 1.0 -1.0 -0.5 0.0 0.5 1.0

0.0 -1.0

0.4

-0.5

0.0

0.5

1.0

HR

HR

Fig. 1.-- Comparison of the shape of posterior distribution for hardness ratios R, C, and HR with the classical Gaussian approximation for high counts (CASE I; left column of the figure) and for low counts (CASE II; right column of the figure). 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).


­ 14 ­

Table 1: Statistical Properties of Bayesian and Classical Metho ds. Metho d Classical metho d CASE I (high counts) Monte Carlo metho d Quadrature Classical metho d CASE II (low counts) Monte Carlo metho d Quadrature Hardness Ratio R C HR R C HR R C HR R C HR R C HR R C HR Coverage R a te 95.0% 99.0% 98.5% 96.0% 96.0% 96.0% 94.5% 96.0% 94.5% 97.5% 100.0% 100.0% 98.0% 98.0% 98.0% 97.0% 99.5% 95.0% Mean Length 1.24 0.54 0.61 1.07 0.44 0.49 1.03 0.43 0.49 192.44 6.02 3.70 15.77 1.54 1.26 8.18 1.51 1.23 Mean Square Error of mo de of mean 0.065 0.013 0.016 0.056 0.069 0.012 0.012 0.016 0.015 0.057 0.069 0.012 0.012 0.016 0.015 93.27 0.27 0.21 85.482 0.113 0.083 20.338 0.112 0.083

0.328 0.078 0.181 0.394 0.074 0.187

In order to compare the statistical properties of estimates based on the Bayesian and classical metho ds, 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 metho d, Monte Carlo sampler, and numerical integration by quadrature. The Bayesian metho ds 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 metho ds 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 metho ds are almost identical. On the other hand, the results of CASE II appear in Figure 3 and confirm that the Bayesian metho ds outperform the classical metho d 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


­ 15 ­
Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-4 -2 0 2 4

100

100

50

50

0

0

-4

-2

0

2

4

0

50

100

150

-4

-2

0

2

4

R

R

R

Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-2 -1 0 1 2

100

100

50

50

0

0

-2

-1

0

1

2

0 -2

50

100

150

-1

0

1

2

C

C

C

Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-1.0 -0.5 0.0 0.5 1.0

100

100

50

50

0

0

-1.0

-0.5

0.0

0.5

1.0

0 -1.0

50

100

150

-0.5

0.0

0.5

1.0

HR

HR

HR

Fig. 2.-- Simulated range of hardness ratios calculated using three different metho ds in the high counts case (CASE I; see Equation 15a) ­ the classical metho d (left), Monte Carlo metho d (middle), and quadrature (right) ­ and for the three types of hardness ratios ­ R (top), C (middle), and HR (bottom). The horizontal lines are the 95% intervals computed for each simulated set of counts, and the vertical white line in each panel represents the true value of the hardness ratio. Note that all the different methods exhibit similar performance in this case.


­ 16 ­
Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-4 -2 0 2 4

100

100

50

50

0

0

-4

-2

0

2

4

0

50

100

150

-4

-2

0

2

4

R

R

R

Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-2 -1 0 1 2

100

100

50

50

0

0

-2

-1

0

1

2

0 -2

50

100

150

-1

0

1

2

C

C

C

Classical method
200 200

Monte Carlo method
200

Quadrature

150

150

Simulation

Simulation

Simulation
-1.0 -0.5 0.0 0.5 1.0

100

100

50

50

0

0

-1.0

-0.5

0.0

0.5

1.0

0 -1.0

50

100

150

-0.5

0.0

0.5

1.0

HR

HR

HR

Fig. 3.-- As in Figure 2, for simulated range of hardness ratios in the low counts case (CASE II; Equation 15b). In this, the low counts case, the the Bayesian metho ds dramatically outperform the classical metho d. 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 metho d. However, this comes at a price: because the summation inside the posterior density is over the detected


­ 17 ­ source counts, quadrature is computationally more expensive as the source counts increase; on the other hand, the Monte Carlo metho d 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 metho d 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 mo de 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 mo de because both left and right limits of HR are bounded. For R and C, on the other hand, the posterior mo de appears to perform better. Thus, the posterior mo de is preferable for both R and C and the posterior mean for HR, as shown in Table 1.

4. 4.1.

Applications

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 asso ciated 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 metho ds 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 mo des of the color


­ 18 ­ 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 metho d (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 metho d, 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 metho d because there are to o 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 metho d if the overly large error bars it pro duces had been used (see discussion of coverage rates in Table 1).


­ 19 ­

Fig. 4.-- The X-ray "color-magnitude" diagram for globular cluster Terzan-5. For each of the detected sources, a form of color C = log10 (S /H ), scaled by a factor 2.5 for similarity with optical colors, is plotted along the abscissa, and the log10 (Luminosity [ergs s-1 ]) in the 0.5-6 keV band, along the ordinate. The mo des 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 metho d (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-bo dy spectra of different temperatures (the open triangles correspond to temperatures of log10 (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 blackbo dy 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 mo dels see Heinke et al. (2006).


­ 20 ­ 4.2. 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 metho d. 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 metho ds are difficult to implement correctly and generally cannot be justified under the classical metho d. 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 pro cesses 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 lo ops and on the heating pro cess (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 hydro dynamically evolving lo ops 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 mo deling of flares to derive physically meaningful parameters such as the length of the lo op, the heating function, etc. Unfortunately, a comprehensive analysis along these lines requires that complex mo del 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 pro duce 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 pro cesses 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 metho d 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


­ 21 ­ be used. (Details of the reduction light-curve of the source during the smaller passbands (soft, 0.25-1 keV; band light-curve peaks at an earlier out by the evolution of the color, C and analysis are given in Wargelin et al. 2006.) The flare is shown in Figure 5, along with light-curves for medium, 1-2 keV; hard, 2-8 keV). Note that the hard time than the softer bands, an impression that is borne (lower panel of Figure 5).

Fig. 5.-- Upper Panel: Light-curves of Ross 154. The observed count rate in the source is shown as a function of time since the beginning of the observation, for counts in the soft (dotted histogram), medium (dashed histogram), hard (dash-dotted histogram), and the total (solid histogram). A bin size of 150 s was chosen for convenience, and the sequence of bins are marked along the light curve for future reference. Lower Panel: Hardness ratio S evolution of Ross 154. The color, C = log10 H is plotted for a bin size of t = 150 s (solid histogram) and t = 40 s (dotted histogram), along with the asso ciated error bars. Larger values of C indicate a softer spectrum, and vice versa.


­ 22 ­ 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 metho d 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 cho osing a single time bin size and a starting time, we have carried out a series of calculations using the Monte Carlo metho d: 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 pro cess. 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 pro cedure, 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 pro duce a single coherent image, thus effectively marginalizing over the bin sizes; such an image includes all the errors asso ciated with the analysis (Figure 7). Note that this type of analysis is not possible using the classical metho d. 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
X-ray luminosity, LX E M · P (T ), where E M = n2 V is the emission measure of coronal plasma with e an electron density ne 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 LX - T conversion, see Wargelin et al. (2006). 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.
k j


­ 23 ­

Fig. 6.-- Color-intensity evolutionary tracks. Spectral hardness increases to the left and source brightness increases towards the top. The paired color and intensity points (shown separately in Figure 5) are plotted for bin sizes of 40 s (left) and 150 s (right). The error bars for each quantity are also plotted as horizontal and vertical bars. The points are connected by a thin solid line to aid in visualization of the temporal sequence (see also Figure 7). Samples of (S , H ) obtained from their posterior probability distributions for each bin and for various starting time phases are displayed as the shaded image underneath the evolutionary tracks.

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.


­ 24 ­ 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.


­ 25 ­

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.


­ 26 ­ 5. 5.1. Discussion

Informative and Non-informative Prior Distributions

A ma jor component of Bayesian analysis is the specification of a prior distribution, the probability distribution assigned a priori to the parameters defining the mo del. This is often deemed controversial because of the apparent sub jectivity 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 metho d, and the specification of a prior distribution co difies 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 metho d 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, enco ded as a -distribution (see van Dyk et al. 2001). In addition, the analysis of previously acquired data will pro duce 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, p ( ) 1 , 2. a Jeffrey's non-informative prior distribution, p() I
1/2

(1 6 a )

, and

(16b)

3. a non-informative prior distribution under a log transformation, p(log ) 1 or equivalently p() 1 , (16c)

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


­ 27 ­ second (Equation 16b) is proportional to the square ro ot 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 co dify 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.

5.2.

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 = log10 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


­ 28 ­ 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 mo de 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 mo de of HR in the low count scenario is -1 when = 0+ in Figure 8.)

5.3.

Advantages

A significant improvement that our metho d 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 metho d works well only with high counts data and fails to give reliable results with low counts, our metho d is valid in all regimes. Indeed, because the observed counts are non-negative integers, it is not appropriate to mo del low counts with Gaussian distributions which are defined on all the real numbers. Because our metho ds are based on Poisson assumptions in a fully mo del 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 metho d 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 metho d gives more precise error bars for the hardness ratios than the classical metho d (as defined by the coverage; see Table 1 and Appendix A), although both metho ds yield similar results. Further, a priori information can be embedded into our mo del and can be updated by


­ 29 ­

15

15

S=2 H=1
pdf

S = 10 H=5
pdf

15

S = 50 H = 25

10

10

pdf

5

5

0

0

0

1

2 R

3

4

5

0

1

2 R

3

4

5

0 0

5

10

1

2 R

3

4

5

40

40

30

30

pdf

pdf

20

20

pdf

10

10

0

0

-1.0

0.0 C

1.0

2.0

-1.0

0.0 C

1.0

2.0

0

10

20

30

S=2 H=1

S = 10 H=5

40

S = 50 H = 25

-1.0

0.0 C

1.0

2.0

S=2 H=1
pdf

S = 10 H=5
pdf

S = 50 H = 25

6

6

pdf

4

4

2

2

0

0

-1.0

-0.5

0.0 HR

0.5

1.0

-1.0

-0.5

0.0 HR

0.5

1.0

0 -1.0

2

4

6

-0.5

0.0 HR

0.5

1.0

Fig. 8.-- The posterior distributions of hardness ratios calculated for R (top row), C (middle row), and HR (bottom row), with different source intensities (low at left and high at right) and prior distributions. The different curves in each figure represent the posterior probability distributions computed with different non-informative prior distribution indices, i.e., = 0+ (solid), = 1/2 (dashed), and = 1 (dotted). At higher counts (large Poisson intensities; right column of figures), the posterior distributions tends to be symmetric and the effect of the prior distributions is minimal, i.e., the posterior distributions with different indices are nearly the same. At small counts (small Poisson intensities; the left column of figures), the non-symmetric shape of the posterior distribution for each hardness ratio becomes clear as do es the effect of the choice of a non-informative prior distribution.


­ 30 ­ data, thereby pro ducing more accurate results as we collect more data.

5.4.

Limitations

Unlike the quantile-width metho d (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 mo del parameters. However, this restriction allows for a more uniform analysis and comparison of results across different instruments and telescopes. Because the Bayesian metho d do es not allow for a simple analytic solution similar to standard error propagation in the Gaussian case, the computational metho ds used to determine the probability distributions become important. We have implemented a Monte Carlo integration metho d (the Gibbs sampler, see Appendix A.1) and a metho d 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 metho d based on quadrature, but care must be taken to ensure that the number of iterations is sufficient to determine the posterior mo de and the required posterior interval with go o d precision. We recommend using the Gibbs sampler for high counts and the quadrature for low counts.

6.

Summary

We have developed a metho d 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


­ 31 ­ 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 metho d and the classical metho d, and show that the Bayesian metho d works even in the low counts regime where the classical metho d fails, and matches the results of the classical metho d at high counts where a Gaussian approximation of the errors is appropriate. We apply the metho d 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 metho d 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 metho d in being able to work around the conflicting limitations imposed by time bin sizes and time resolution. Additional applications of these metho ds 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 mo del, 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 mo del 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 co de development, and Dong-Wo o Kim and Peter Edmonds for useful discussions. The authors gratefully acknowledge funding for this pro ject 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 pro duct 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. Esch, P. Freeman, H. Kang, X.-L. Meng, J. Scargle, E. Sourlas, and Y. Yu.


­ 32 ­ A. A.1. Computing the Hardness Ratios

Monte Carlo Integration: Gibbs Sampler

Monte Carlo metho ds can be used in Bayesian statistical inference if representative values of the mo del parameters can be simulated from the posterior distribution. Given a sufficient simulation size, these simulated values can be used to summarize and describe the posterior distribution and hence the Bayesian statistical inference. For example, the posterior mean can be computed by averaging these simulated values and posterior intervals are computed using quantiles of the the simulated values. Markov chain Monte Carlo (MCMC) metho ds (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 mo del parameters conditioning on the others. A more detailed description of these metho ds can be found in van Dyk (2003) and the Appendix of van Dyk et al. (2001). In order to construct a Gibbs sampler, we embed the statistical mo del described in §2.1 into a larger mo del. 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 mo deled as independent Poisson random variables, S Poisson(eS · S ) and H Poisson(eH · H ), (A1)

and the background counts in the source exposure area as independent Poisson random variables, S Poisson(eS · S ) and H Poisson(eH · H ), (A2)

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 metho d 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 pro duces Monte Carlo draws of S and H along with sto chastically 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


­ 33 ­ 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 p(S , S , S |S, BS ) p(S |S , S )p(BS |S )p(S |S )p(S )p(S ) 1 S -S +S1 -1 B+S +S3 -1 e- (S - S )!S !

(eS +S2 )S -(eS +r ·eS +S4 )S

.

(A3)

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 (t) (t) parameters, S and S , STEP 1 is given by ST 1 : Draw S
(t+1)

EP

from S |(S , S , S, BS ) Binomial S,

(t)

(t)


(t) S

(t) S (t) S

+

,

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 ST ST
EP

2 : Draw 3 : Draw eS + S4 ),

(t+1) S

from S |(S , S from S |(
(t+1) S

(t)

(t+1)

, S, BS ) Gamma(S - S

(t+1)

+ S1 , eS + S2 ) and
(t+1)

EP

(t+1) S

, S

(t+1)

, S, BS ) Gamma(BS + S

+ S3 , eS + r ·

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 (t) (t) sample {S , H , t = t0 + 1, . . . , T } for a sufficiently long burn-in perio dl t0 . 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 to ols for determining t0 , 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 - t0 ) 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.
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.
l


­ 34 ­ A.2. Numerical Integration: Gaussian Quadrature

Instead of intro ducing 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 mo dels 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 p(S , S |S, BS ) = =
00

p(S )p(S )p(S |S , S )p(BS |S ) p(S )p(S )p(S |S , S )p(BS |S )dS d ( S + ( S +

00

S S1 -1 BS +S3 -1 -(e + ) -(e +r ·e + ) e S S2 S e S S S4 S S S )S S S -1 BS +S3 -1 -(e + ) -(e +r ·e + ) e S S2 S e S S S4 S dS S )S S 1 S

. d S (A4)

Then, the binomial expansion to (S + S )S , i.e.,
S

( S + S ) =
j =0

S

(S + 1 ) j (j + 1)(S - j + 1) S

S -j S

,

(A5)

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
S

p(S , S |S, BS ) =

j =0 S

(S + 1 ) (j + 1)(S - j + 1)

j +S1 -1 S -j +BS +S3 -1 -(eS +S )S -(eS +r ·eS +S )S 2 4 e S e S

(A6)
S
3

j =0

(S - j + BS + S3 ) (S + 1 ) · (j + 1)(S - j + 1) (eS + r · eS + S4 )S -j +BS +

(j + S1 ) · (eS + S2 )j +S

1

and then S is integrated out of Equation A6, i.e.,
S

p(S |S, BS ) =

j =0 S

(S - j + BS + S3 ) 1 · (j + 1)(S - j + 1) (eS + r · eS + S4 )S -j +BS +

S


3

j +S1 -1 -(eS +S )S 2 e S

(A7) .
S
3

j =0

1 (S - j + BS + S3 ) · (j + 1)(S - j + 1) (eS + r · eS + S4 )S -j +BS +

(j + S1 ) · (eS + S2 )j +S

1

Here the prior independence of S and H implies that the joint posterior distribution of these two intensities is decomposed into a pro duct of their marginal posterior distributions, i.e., p(S , H |S, H, BS , BH ) = p(S |S, BS )p(H |H, BH ). (A8)


­ 35 ­ This assumption of independence between S and H can be relaxed, and hierarchical mo deling 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 p(R, H |S, H, BS , BH ) dR dH ( S , H ) = p(S , H |S, H, BS , BH ) d S d (R, H ) = p(RH , H |S, H, BS , BH )H dR dH ,

H

2. the posterior distribution of C is obtained after integrating H out of p(C, H |S, H, BS , BH ) dC dH ( S , H ) = p(S , H |S, H, BS , BH ) d S d H (C, H ) = p(10C H , H |S, H, BS , BH )10C ln(10)H dC dH , and 3. the posterior distribution of HR is obtained after integrating = S + p(HR, |S, H, BS , BH ) dHR d ( S , H ) = p(S , H |S, H, BS , BH ) d S d H (H R , ) (1 - HR) (1 + HR) =p S, H, BS , BH , dH R d . 2 2 2 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
H

out of

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


­ 36 ­

Table 2: Legend Key for Table 3. Hard band source intensity, H [counts (source area)-1 ] Coverage rate of Average length of intervals with = 0+ intervals with = 0 Soft band source intensity S [counts (source area)-1 ]

+

Coverage rate of Average length of intervals with = 1/2 intervals with = 1/2 Coverage rate of intervals with = 1 Average length of 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 mo del 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 pro duce much shorter posterior intervals and maintain high coverage rates. With high counts, however, the choice of a non-informative prior distribution do es not have a


­ 37 ­ 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 metho d (last column). The distinction between the Bayesian and classical metho ds 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 enco ded 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 Metho d with Different (see Table 2 for key) Non-Informative Prior Distribution Indices ( = 0+ , 1/2, 1). 100 0.5 100 100 100 1.0 100 100 100 2.0 99.9 99.1 100 4.0 100 98.4 100 8.0 99.9 97.7 100 16.0 99.9 97.5 100 32.0 99.9 98.2 100 64.0 99.9 98.5 0 % % % % % % % % % % % % % % % % % % % % % % % % .5 23 5 2 22 4 2 21 4 2 18 3 2 13 3 1 11 3 1 10 2 1 10 3 1 .9 .1 .8 .9 .9 .7 .4 .5 .5 .6 .9 .2 .6 .2 .9 .3 .1 .9 .8 .9 .8 .9 .0 .8 9 0 7 5 1 3 2 2 1 0 2 4 6 9 9 7 0 1 2 8 0 9 2 2 100 100 100 100 100 100 100 99.9 99.5 100 100 99.2 100 100 98.7 100 99.8 98.0 100 99.6 97.4 100 99.6 97.4 1 % % % % % % % % % % % % % % % % % % % % % % % % .0 22 4 2 21 4 2 19 4 2 16 3 2 12 3 1 9 2 1 9 2 1 9 2 1 .6 .8 .7 .5 .6 .5 .5 .3 .3 .8 .6 .0 .2 .0 .8 .9 .7 .7 .5 .7 .6 .5 .7 .6 7 8 1 9 9 8 7 0 7 0 8 9 8 4 5 6 9 3 6 1 5 2 2 6 100 99.9 99.0 100 99.7 99.2 99.8 99.7 99.7 99.9 99.6 99.0 99.3 99.7 98.7 99.8 99.2 98.2 99.3 99.0 96.0 99.7 99.0 97.1 2 % % % % % % % % % % % % % % % % % % % % % % % % .0 21 4 2 19 4 2 17 3 2 15 3 1 10 2 1 8 2 1 8 2 1 7 2 1 .1 .5 .5 .6 .3 .3 .5 .8 .1 .1 .3 .8 .9 .6 .6 .4 .3 .4 .0 .2 .4 .8 .2 .4 9 1 0 1 1 7 2 9 4 8 0 7 3 4 3 9 4 9 0 1 0 4 3 0 100 99.9 98.4 100 100 98.9 99.7 99.6 98.8 99.2 99.2 99.1 98.3 97.8 97.5 97.7 97.8 97.2 97.7 97.8 97.1 97.7 98.1 97.8 4 % % % % % % % % % % % % % % % % % % % % % % % % .0 19 3 2 17 3 2 15 3 1 13 2 1 8 1 1 6 1 1 5 1 1 5 1 1 .0 .9 .2 .4 .7 .1 .0 .2 .8 .1 .6 .5 .7 .8 .2 .1 .6 .1 .7 .5 .0 .6 .5 .0 2 5 5 0 7 4 4 7 7 1 4 9 4 6 8 1 0 4 4 1 6 9 2 5 100 100 98.4 100 99.9 98.5 99.6 99.5 97.8 97.9 98.5 97.7 99.2 98.8 98.1 98.6 98.0 96.7 97.3 96.5 96.0 97.7 97.7 96.6
H



S

8 % % % % % % % % % % % % % % % % % % % % % % % %

.0 14 3 2 12 3 1 10 2 1 8 1 1 5 1 1 2 0 0 2 0 0 2 0 0 .0 .3 .0 .2 .0 .8 .5 .5 .6 .8 .8 .2 .1 .2 .0 .4 .9 .8 .3 .8 .7 .1 .8 .7 7 5 2 4 4 5 4 6 0 7 9 9 3 5 1 1 5 3 3 8 5 1 3 1 100 99.9 97.2 100 99.8 97.3 99.3 99.2 98.1 97.9 98.0 97.3 98.2 97.6 97.3 97.9 96.6 96.6 95.7 94.4 94.5 96.3 94.8 94.2

16.0 32.0 64.0 % 11.32 100 % 11.14 100 % 11.11 % 3.09 99.7 % 3.02 100 % 3.03 % 1.90 97.0 % 1.82 98.4 % 1.83 % 10.03 100 % 9.84 100 % 9.51 % 2.80 99.7 % 2.74 99.8 % 2.71 % 1.73 97.7 % 1.67 98.0 % 1.65 % 8.45 99.9 % 7.95 99.8 % 7.82 % 2.33 99.0 % 2.21 99.4 % 2.21 % 1.49 97.3 % 1.40 97.7 % 1.39 % 6.06 97.0 % 5.85 96.8 % 5.63 % 1.60 97.9 % 1.53 96.6 % 1.49 % 1.14 97.1 % 1.08 97.4 % 1.04 % 2.50 97.8 % 2.26 97.4 % 2.14 % 0.96 97.2 % 0.86 97.4 % 0.83 % 0.83 96.3 % 0.75 97.1 % 0.71 % 0.73 96.7 % 0.60 95.5 % 0.54 % 0.65 95.6 % 0.55 94.4 % 0.50 % 0.63 95.3 % 0.54 93.7 % 0.49 % 0.59 97.1 % 0.43 93.8 % 0.37 % 0.55 97.2 % 0.43 94.3 % 0.36 % 0.54 97.2 % 0.42 94.5 % 0.36 % 0.53 94.6 % 0.36 96.5 % 0.29 % 0.50 94.8 % 0.36 96.5 % 0.29 % 0.49 94.9 % 0.36 96.5 % 0.28


­ 39 ­

40

40

40

40
= 1.0

= 0.0

+

= 0.5

30

30

30

20

20

20

10

10

10

0

0

0

90 92 94 96 98

90 92 94 96 98

90 92 94 96 98

0 90 92 94 96 98

Coverage (%)
8 8
= 0.0
+

Coverage (%)
8
= 0.5

Coverage (%)
8
= 1.0

10

20

30

Classical method

Coverage (%)
Classical method

6

6

6

4

4

4

2

2

2

0

0

0

90 92 94 96 98

90 92 94 96 98

90 92 94 96 98

0 90 92 94 96 98

Coverage (%)

Coverage (%)

Coverage (%)

2

4

6

Coverage (%)

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 metho d, 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%.


­ 40 ­

Fig. 10.-- Graphical representation of the data in Table 3. The symbols are lo cated 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.


­ 41 ­ REFERENCES Alexander, D.M., Brandt, W.N., Hornschemeier, A.E., Garmire, G.P., Schneider, D.P., Bauer, F.E., & Griffiths, R.E., 2001, AJ, 122, 2156 Babu, G.J., & Feigelson, E.D., 1997, Statistical Challenges in Mo dern Astronomy, SpringerVerlag:New York Bradt, H., Mayer, W., Buff, J., Clark, G.W., Doxsey, R., Hearn, D., Jernigan, G., Joss, P.C., Laufer, B., Lewin, W., Li, F., Matilsky, T., McClinto ck, J., Primini, F., Rappaport, S., & Schnopper, H., 1976, ApJ, 204, L67 Brandt, W.N., Hornschemeier, A.E., Alexander, D.M., Garmire, G.P., Schneider, D.P., Bro os, P.S., Townsley, L.K., Bautz, M.W., Feigelson, E.D., & Griffiths, R.E., 2001a, AJ, 122, 1 Brandt, W.N., Alexander, D.M., Hornschemeier, A.E., Garmire, G.P., Schneider, D.P., Barger, A.J., Bauer, F.E., Bro os, P.S., Cowie, L.L., Townsley, L.K., Burrows, D.N., Chartas, G., Feigelson, E.D., Griffiths, R.E., Nousek, J.A., & Sargent, W.L.W., 2001b, AJ, 122, 2810 Brandt, W.N., & Hasinger, G., 2005, ARAA, 43, 827 Brown, E.F., Bildsten, L., & Rutledge, R.E., 1998, ApJ, 504, 95 Campana, S., Colpi, M., Mereghetti, S., Stella, L., & Tavani, M., 1998, A&ARv, 8, 279 Casella, G. & Berger, R.L., 2002, Statistical Inference, 2nd edition. Duxbury. Cowles, M.K., & Carlin, B.P., 1996, J. Am. Stat. Assn., 91, 883 Esch, D.N., 2003, Ph.D. Thesis, Department of Statistics, Harvard University Feigelson, E.D., & Babu, G.J., 1992, Statistical Challenges in Mo dern Astronomy, SpringerVerlag:Berlin Heidelberg New York Feigelson, E.D., & Babu, G.J., 2003, Statistical challenges in astronomy. Third Statistical Challenges in Mo dern Astronomy (SCMA III) Conference, University Park, PA, USA, July 18 - 21 2001, Springer:New York Gehrels, N., 1986, ApJ, 303, 336 Geman, S., & Geman, D., 1984, IEEE, 6, 721 Giacconi, R., Rosati, P., Tozzi, P., Nonino, M., Hasinger, G., Norman, C., Bergeron, J., Borgani, S., Gilli, R., Gilmozzi, R., & Zheng, W., 2001, ApJ, 551, 624 Gregory, P.C., & Loredo, T.J., 1992, ApJ, 398, 146 Hasinger, G., & van der Klis, M., 1989, A&A, 225, 79


­ 42 ­ Heinke, C.O., Grindlay, J.E., Edmonds, P.D., Lloyd, D.A., Murray, S.S., Cohn, H.N., & Lugger, P.M., 2003, ApJ, 598, 501 Hong, J.S., Schlegel, E.M., & Grindlay, J.E., 2004, ApJ, 614, 508 Gillessen, S., & Harney, H.L., 2005, A&A, 430, 355 Jansen, F., Lumb, D., Altieri, B., Clavel, J., Ehle, M., Erd, C., Gabriel, C., Guainazzi, M., Gondoin, P., Much, R., Munoz, R., Santos, M., Schartel, N., Texier, D., & Vacanti, G., 2001, A&A, 365, L1 Kashyap, V.L., Drake, J.J., Gudel, M., & Audard, M., 2002, ApJ, 580, 1118 ¨ Kashyap, V.L., & Drake, J.J., 1998, ApJ, 503, 450 Kim, D.-W., Barkhouse, W.A., Colmenero, E.R., Green, P.J., Kim, M., Mossman, A., Schlegel, E., Silverman, J.D., Aldcroft, T., Ivezic, Z., Anderson, C., Kashyap. V., Tananbaum, H., & Wilkes, B.J., 2005, submitted to ApJ Kim, D.-W., Fabbiano, G., & Trinchieri, G., 1992, ApJ, 393, 134 Pallavicini, R., Tagliaferri, G., & Stella, L., 1990, A&A, 228, 403 Protassov, R., van Dyk, D.A., Connors, A., Kashyap, V.L., & Siemiginowska, A., 2002, ApJ, 571, 545 Reale, F., Betta, R., Peres, G., Serio, S., & McTiernan, J., 1997, A&A, 325, 782 Reale, F., Gudel, M., Peres, G., & Audard, M., 2004, A&A, 416, 733 ¨ Reale, F., Serio, S., & Peres, G., 1993, A&A, 272, 486 Rosner, R., Golub, L., & Vaiana, G.S., 1985, ARAA, 23, 413 Schmitt, J.H.M.M., & Favata, F., 1999, Nature, 401, 44 Serio, S., Reale, F., Jakimiec, J., Sylwester, B., & Sylwester, J. 1991, A&A, 241, 197 Silverman, J.D., Green, P.J., Barkhouse, W.A., Kim, D.-W., Aldcroft, T.L., Cameron, R.A., Wilkes, B.J., Mossman, A., Ghosh, H., Tananbaum, H., Smith, M.G., Smith, R.C., Smith, P.S., Foltz, C., Wik, D., & Jannuzi, B.T., 2005, ApJ, 618, 123 Tanner, M.A. & Wong, W.H., 1987, J. Am. Stat. Assn., 82, 528 Tuohy, I.R., Garmire, G.P., Lamb, F.K., & Mason, K.O., 1978, ApJ, 226, L17 Vaiana, G.S., Cassinelli, J.P., Fabbiano, G., Giacconi, R., Golub, L., Gorenstein, P., Haisch, B.M., Harnden, F.R.,Jr., Johnson, H.M., Linsky, J.L., Maxson, C.W., Mewe, R., Rosner, R., Seward, F., Topka, K., & Zwaan, C., 1981, ApJ, 245, 163 van den Oord, G.H.J., & Mewe, R., 1989, A&A, 213


­ 43 ­ van Dyk, D.A., 2003, in Statistical chal lenges in astronomy. Third Statistical Chal lenges in Modern Astronomy (SCMA III) Conference, University Park, PA, USA, July 2001, Eds. E.D.Feigelson, G.J.Babu, New York: Springer, p41 van Dyk, D.A., Connors, A., Kashyap, V.L., & Siemiginowska, A., 2001, ApJ 548, 224 van Dyk, D.A., & Hans, C.M., 2002, in Spatial Cluster Mo delling, Eds. D.Denison & A.Lawson, CRC Press: London, p175 van Dyk, D.A., & Kang, H., 2004, Stat. Sci, 19, 275 Wargelin, B.J., Kashyap, V.L., Drake, J.J., Garc´a-Alvarez, D., & Ratzlaff, P.W., 2006, i submitted to ApJ Weisskopf, M.C., Tananbaum, H.D., Van Speybro eck, L.P., & O'Dell, S.L., 2000, Pro c. SPIE, 4012, 2 Wichura, M.J., 1989, Technical Report No. 257, Department of Statistics, The University of Chicago Zamorani, G., Henry, J.P., Maccacaro, T., Tananbaum, H., Soltan, A., Avni, Y., Liebert, J., Sto cke, J., Strittmatter, P.A., Weymann, R.J., Smith, M.G., & Condon, J.J., 1981, ApJ, 245, 357

A This preprint was prepared with the AAS L TEX macros v5.0.