Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : 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 Êîäèðîâêà:

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 E