format: text

text:
# Analysis of Energy Spectra with Low Photon Counts via Bayesian Posterior Simulation

## 1 Introduction

The ever increasing power and sophistication of today’s high energy instruments is opening a new realm of high quality data that is quickly pushing beyond the capabilities of the “classical” data-analysis methods in common use. In this paper, we present an innovative implementation of state-of-the-art statistical methods for fitting high resolution spectra from the Chandra X-ray Observatory. The common “folk-wisdom” of how to bin data, subtract background counts, propagate errors, and, for example, estimate the significance of a spectral line profile are unreliable and can lead to unacceptable results (see Loredo, 1993; Nousek, 1993; Feigelson & Babu, 1997; Siemiginowska et al., 1997; Zimmerman, 1997, for discussion). For example, binning data sacrifices the resolution of the instrument, subtracting background can lead to negative counts with unpredictable results, and statistical black boxes such as the $`\chi ^2`$ and Cash statistics (Lampton et al., 1976; Cash, 1979) although often useful, may not be equipped to answer standard questions (e.g., van Dyk & Protassov, 2000). Some authors have suggested solutions to such problems which involve ad hoc adaptations of commonly used methods (e.g., Gehrels, 1986; Collura et al, 1987; Mighell, 1999). Unfortunately, when such solutions are not rooted in a theoretical framework, they have no justification beyond problems which are more-or-less the same as the simulation studies which justify them and we are often forced into additional ad hoc adaptations. This approach is difficult to justify in light of modern statistical methods which address reasonable model assumptions directly. Thus, in recent years, astrophysicists have increasingly turned to likelihood based (e.g., Lucy, 1974; Cash, 1979; Schmitt, 1985; Sciortino & Micela, 1992) and Bayesian methods (e.g., Bijaoui, 1971; Richardson, 1972; Gregory & Loredo, 1992; Loredo, 1993; Connors, 1997; Siemiginowska, 1997; Kashyap & Drake, 1998; Freeman et al., 1999). The primary purpose of this paper is to illustrate how Bayesian methods can provide practical answers to outstanding real problems which standard methods are not able to handle. The methods described here are equipped with readily available parameter estimates, credible intervals, error bars, model checking techniques, methods for combining information from multiple sources, etc.—all within a flexible theoretical framework and without reliance on asymptotic Gaussian approximations.

We illustrate Bayesian data analysis via two detailed examples. The analysis of Quasar S5 0014+813 offers a straight forward introduction to our methods and the extremely low count Hyprid-Chromosphere Supergiant Star $`\alpha `$ TrA observation shows how we tackle a previously intractable analysis. Together, these examples demonstrate the power of Bayesian methods to handle highly structured models designed to reflect the structure in both the source spectrum and the data collection process. Our methods avoid the binning of counts and thus the sacrificing of high-resolution information required by standard data analysis methods. The analysis of Quasar S5 0014+813 is consistent with the available standard analysis which relies on extra binning and the removal of the high-energy low-count tail.

We emphasize that we model not only the source spectrum, but also other stochastic components of data collection and the instrument such as background contamination and instrument response. Generally we refer to our stochastic representation of the entire process as the (statistical) model. For clarity we refer to the spectral or physical model as the source model and to the model for the observed (PHA111Pulse Height Amplitude, originally in proportional counters, the number of electrons produced by a photon, hence the amplitude of the current pulse registered by the detector electronics. The term now refers to the measure of the energy deposited on the detector (as opposed to the true energy).) counts as the observed-data model. In our detailed example, we develop a model and algorithms for spectral analysis of high-energy (or other) data using a Poisson222Recall, a random variable $`X`$ is said to follow a Poisson distribution with parameter or intensity $`\lambda `$ if $`\mathrm{Pr}(X=x)=e^\lambda \lambda ^x/x!`$. In this case $`\mathrm{E}(X)=\lambda `$ and we often write $`X\stackrel{\mathrm{d}}{}\mathrm{Poisson}(\lambda )`$ (read as $`X`$ is distributed as Poisson with intensity $`\lambda `$). This representation conditions on the intensity parameter, $`\lambda `$, which in turn may vary. process for photon arrivals. We allow for: stochastic instrument response via a photon redistribution matrix; the absorption of photons; the effective area333 The effective area of the telescope is the fraction of the true geometric area that the telescope presents to sky. This varies with energy. of the telescope; and background contamination of the source. In particular, we model information on background emissions as the realization of a second Poisson process (cf. Loredo, 1993), thereby eliminating the need to directly subtract off the background counts and the rather embarrassing resulting problem of negative photon counts. The source energy spectrum is modeled as a mixture of several (Gaussian) line profiles and a Generalized Linear Model (GLM)444In a GLM we assume a transformation (e.g., log) of the model is linear in a set of independent variables. We emphasize that this is not equivalent to transforming the data and proceeding with linear regression. A generalized linear model utilizes the likelihood of the assumed model which may not be Gaussian. (E.g., the assumed model may be Poisson.) See Section 3.3 for details.  (e.g., McCullagh & Nelder, 1989) which accounts for the continuum. GLMs have become the standard statistical method for incorporating information contained in independent variables (as in regression) into many non-Gaussian models and are thus an obvious but innovative choice in this setting.

In addition to several Markov chain Monte Carlo (MCMC) algorithms, we describe and use data augmentation, an important statistical method for Bayesian (and other) analysis. Data augmentation is an elegant computational construct allowing us to take advantage of the fact that if it were possible to collect additional data, statistical analysis would be greatly simplified. This is true regardless of why the so-called “missing data” are not observed. For example, if we were able to record the counts due to background contamination in addition to total counts in each bin, it would, of course, be a trivial task to account for the background. There is a large class of powerful statistical methods designed for “missing data” problems. With the insight that “true” values of quantities recorded with measurement error can be regarded as “missing data,” these methods can usefully be applied to almost any astrophysical problem. In particular, we can treat the true image (before instrument response), the absorbed photon counts, and unbinned energies as “missing data” to account for instrument response, absorption, and binning, respectively.

This introduction fore shadows the tone of the paper: in the process of developing new Bayesian methods, we describe and utilize state-of-the-art statistical reasoning, methods, and algorithms, whereby we explore a larger statistical framework for our problem of interest. Although we introduce many new tools and use terminology that may be unfamiliar, we endeavor to write in a manner accessible to astrophysicists and believe the resulting methods justify the required interdisciplinary work. Table 1 indexes terminology used in the paper that may be unfamiliar to some readers. To aid in the translation from standard statistical notation to standard astrophysical usage, many equations have been written according to both standards when notation is introduced. One notational convention is worthy on mention; we use superscripts to identify model components (e.g., background or absorption).

The paper is organized as follows. After a brief overview of the fundamentals of Bayesian analysis in Section 2, we lay out our hierarchical statistical model555A hierarchical (statistical) model is formulated in terms of unobserved quantities, which are themselves statistically modeled. For example, we may assume photons first arrive at a detector according to a Poisson process and them are randomly redistributed according to a photon redistribution matrix. A hierarchal model separates these two random process into two levels of a structured model. which summarizes the photon collection process and parameterizes many relevant aspects of the energy spectrum in Section 3. Two examples which aim to illustrate a typical data analysis and the advantages of Bayesian methods in this setting are given in Section 4. Section 5 contains brief concluding remarks. Finally, in two appendices, we outline such important general MCMC methods as data augmentation, the Gibbs sampler, the Metropolis-Hastings algorithm, and judging convergence using multiple chains. We also describe in detail how we use these algorithms to fit the hierarchical source model of Section 3.

## 2 Bayesian Analysis

In this section we outline several important methodological and computational issues involved with Bayesian analysis using a simple model that accounts for background in a simplified Poisson process to motivate and illustrate ideas. Our introduction is brief and we encourage interested readers to consult one of the several high-quality recent texts on the subject such as Gelman, Carlin, Stern, & Rubin (1995), Carlin & Louis (1996), Gilks, Richardson, & Spiegelhalter (1996), and Sivia (1996). In Section 3, we show how the ideas developed here can be used for a detailed spectral analysis.

### 2.1 Prior, Sampling, and Posterior Distributions

Bayesian probability analysis is fundamentally based on one simple result known as Bayes’s Theorem which allows us to update a probability distribution based on new data or other information. In particular, knowledge about a (vector) model parameter, $`𝜽`$, is summarized by a probability distribution, $`p(𝜽)`$, such that $`\mathrm{Pr}(𝜽R)=_Rp(𝜽)𝑑𝜽`$ for any region666The notation $`\mathrm{Pr}(𝜽R)`$ represents the probability that $`𝜽`$ is in the region $`R`$ and is computed as the integral of the probability distribution of $`𝜽`$ over $`R`$. $`R`$. Bayes’s Theorem states

$$p(𝜽|𝐘,𝐈)=\frac{p(𝐘|𝜽,𝐈)p(𝜽|𝐈)}{p(𝐘|𝐈)},$$
(1)
where $`𝐘`$ are the observed data or other new information pertaining to $`𝜽`$ and $`𝐈`$ represents any initial information known before $`𝐘`$ is observed. Here, $`p(𝜽|𝐈)`$ represents our knowledge prior to observing $`𝐘`$ and is called the prior distribution. The sampling distribution or likelihood, $`p(𝐘|𝜽,𝐈)`$ represents the likelihood of the data given the model parameters, and $`p(𝜽|𝐘,𝐈)`$ represents our updated knowledge regarding $`𝜽`$ after observing $`𝐘`$ and is called the posterior distribution. Finally, $`p(𝐘|𝐈)`$ represents the unconditional distribution of $`𝐘`$ and acts as the normalizing constant for $`p(𝜽|𝐘,𝐈)`$. The functional form of Bayes’s Theorem describes how our prior knowledge should be updated in light of information contained in the data. The likelihood or sampling distribution is the basis for many standard statistical techniques, while the prior and posterior distributions are specific to Bayesian analysis.

To illustrate Bayes’s Theorem, suppose we have observed counts, $`Y`$, contaminated with background in a (source) exposure and have observed a second exposure of pure background. Throughout this section, we assume the source exposure is $`\tau ^S`$ minutes and the pure background exposure is $`\tau ^B`$ minutes with both exposures using the same area of the detector. (We generally use superscripts to represent photon “sources,” e.g., background or source. Occasionally we use superscripts for powers; for clarity we place powers outside parentheses.) To model the source exposure, we assume $`Y`$ follows a Poisson distribution with intensity $`\lambda ^B+\lambda ^S`$, where $`\lambda ^B`$ and $`\lambda ^S`$ represent the expected counts during the source exposure due to background and source respectively. Thus, the likelihood is

$$p(Y|\lambda ^B,\lambda ^S,𝐈)=e^{(\lambda ^B+\lambda ^S)}(\lambda ^B+\lambda ^S)^Y/Y!\mathrm{for}Y=0,1,2,\mathrm{}$$
(2)
We wish to estimate $`\lambda ^S`$ and treat $`\lambda ^B`$ as a nuisance parameter, a parameter that is of little interest, but must be included in the model. As is detailed below, an important advantage of Bayesian methods is their ability to handle nuisance parameters by computing the marginal posterior distribution of the parameters of interest. The name “marginal” distribution originates with two-way tables of counts where the table margins sum over one of the variables to give the distribution of the other variable alone (e.g., its marginal distribution). Likewise the marginal distribution of the parameter of interest is computed by integrating over (i.e., averaging over) the nuisance parameter.

At this point, we specify a prior distribution which allows us to include a priori knowledge (e.g., “allowed parameter ranges”) from other experiments or other scientific information. One of the primary advantages of Bayesian analysis is a well-defined mechanism for the inclusion of information outside the current data set. In the absence of prior information, we use diffuse or so-called non-informative priors, which are ordinarily flat and have minimal influence on the final analysis. The prior distribution itself may be conveniently parameterized using a set of hyperparameters that can be varied to represent the researcher’s knowledge about the value of the model parameters and the degree of certainty of this knowledge. For example, we use the $`\gamma `$ distribution777The $`\gamma `$ $`(\alpha ,\beta )`$ distribution is a continuous distribution on the positive real line with probability density function $`p(Y)=\beta ^\alpha Y^{\alpha 1}e^{\beta Y}/\mathrm{\Gamma }(\alpha )`$, expected value $`\alpha /\beta `$, and variance $`\alpha /\beta ^2`$ for positive $`\alpha `$ and $`\beta `$. to parameterize prior information for $`\lambda ^B`$ and $`\lambda ^S`$;

$$\lambda ^B|𝐈\stackrel{\mathrm{d}}{}\gamma (\alpha ^\mathrm{B},\beta ^\mathrm{B})\text{and}\lambda ^\mathrm{S}|𝐈\stackrel{\mathrm{d}}{}\gamma (\alpha ^\mathrm{S},\beta ^\mathrm{S}),$$
(3)
i.e.,

$$p(\lambda ^B|𝐈)=(\beta ^B)^{\alpha ^B}(\lambda ^B)^{\alpha ^B1}e^{\beta ^B\lambda ^B}/\mathrm{\Gamma }(\alpha ^B)\text{and}p(\lambda ^S|𝐈)=(\beta ^S)^{\alpha ^S}(\lambda ^S)^{\alpha ^S1}e^{\beta ^S\lambda ^S}/\mathrm{\Gamma }(\alpha ^S),$$
(4)
where the notation $`\stackrel{\mathrm{d}}{}`$ is read “follows the distribution” and $`\lambda ^B`$ and $`\lambda ^S`$ are assumed a priori independent. The $`\gamma `$ prior on $`\lambda ^B`$ is mathematically equivalent to a Poisson likelihood resulting from a count equal to $`\alpha ^B1`$ obtained with an exposure of $`\beta ^B`$ times that of the source exposure. (By mathematically equivalent, we mean the prior on $`\lambda ^B`$ is proportional to a Poisson likelihood as a function of $`\lambda ^B`$.) This leads to a natural choice of $`p(\lambda ^B|𝐈)=\gamma `$ $`(\alpha ^B=Y^B+1,\beta ^B=\tau ^B/\tau ^S)`$, where $`Y^B`$ are the counts from the background exposure. Notice that here (and throughout the paper) we explicitly incorporate information from the background exposure into the analysis via the prior distribution on $`\lambda ^B`$. Thus, the counts from the source exposure, $`Y`$, are treated as the observed data $`𝐘`$ in Equation 1. We refer interested readers to Gelman et al. (1995, Section 2.7) for further discussion and examples of the $`\gamma `$ prior distribution with Poisson data.

The equivalence of the $`\gamma `$ prior for $`\lambda ^S`$ and $`\alpha ^S1`$ counts during an exposure of $`\tau ^S\beta ^S`$ minutes leads to a natural interpretation of the hyperparameters—for a relatively non-informative prior we choose $`\beta ^S`$ much less than one. To illustrate this, we consider two priors, one non-informative and improper888An improper distribution is a distribution that is not integrable and thus is not technically a distribution. One should use improper prior distributions only with great care since in some cases they lead to improper posterior distributions which are uninterpretable., $`p(\lambda ^S|𝐈)^{[1]}=\gamma (1,0)1`$, (illustrated as a dotted line in Figure 1); and one informative, where, let us say we know from other means that three counts are to be expected in the same exposure time, hence $`p(\lambda ^S|𝐈)^{[2]}=\gamma (4,1)`$ (solid line in Figure 1). This choice of informative prior is only an example—$`\gamma (4,1)`$ corresponds to Poisson likelihood resulting from 3 counts with an exposure time equal to the source exposure (ten minutes). This is a rather informative prior distribution and is chosen to illustrate the effect of very informative prior. The non-informative prior contains information equivalent to zero counts in an exposure of zero seconds.

Using Bayes’s Theorem, with $`𝜽=(\lambda ^B,\lambda ^S)`$, we can combine the $`\gamma `$ priors and the likelihood given in Equation 2 to compute the posterior distribution,

$$p(\lambda ^B,\lambda ^S|Y,𝐈)e^{[\lambda ^B(\beta ^B+1)+\lambda ^S(\beta ^S+1)]}(\lambda ^B+\lambda ^S)^Y(\lambda ^B)^{\alpha ^B1}(\lambda ^S)^{\alpha ^S1},$$
(5)
for $`\lambda ^B0,\lambda ^S0`$.

Nuisance parameters such as $`\lambda ^B`$ pose a monumental difficulty for classical statistical analysis which often relies on fixing nuisance parameters at estimated values. Unfortunately, this does not account for uncertainty in their estimates and thus tends to be anti-conservative. (Likewise, floating nuisance parameters or “propagating errors” when computing error bars is essentially a Gaussian assumption, which can lead to unpredictable results when such an assumption is not justified.) The Bayesian solution averages over the posterior distribution (i.e., the uncertainty) of the nuisance parameter by computing the marginal (posterior) distribution of the parameters of interest without Gaussian approximations. For example, the marginal posterior distribution of $`\lambda ^S`$ can be computed by (numerical) integration,

$$p(\lambda ^S|Y,𝐈)=_0^{\mathrm{}}p(\lambda ^B,\lambda ^S|Y,𝐈)𝑑\lambda ^B,$$
(6)
and is illustrated for the two priors for $`\lambda ^S`$ in the second plot of Figure 1, where we assign $`Y=1`$ and $`Y^B=48`$. In this example, direct subtraction of background would leave a “negative count” of $`1`$; no such difficulty occurs with the Bayesian analysis. (See Loredo (1993) for another derivation of the marginal distribution of $`\lambda ^S`$ in this setting.)

Since the source count is small relative to the background count, we expect a small $`\lambda ^S`$. Although this is evident in both posterior distributions in Figure 1, the highly informative prior distribution centered at $`\lambda ^S=3`$ pulls the (solid) posterior towards higher values, thus illustrating the effect of an informative prior distribution. Such sensitivity analyses often play an important part in Bayesian (or other) data analyses, since they investigate the sensitivity of the results to the statistical assumptions (e.g., the choice of prior distribution).

The posterior distributions should be interpreted as probability distributions representing the combined information in the prior and data. For example, a region, $`R`$, such that $`_Rp(𝜽|𝐘,𝐈)𝑑𝜽=\zeta `$ is called a $`\zeta `$-level credible interval (or credible region if $`𝜽`$ is multidimensional) and we can say $`\mathrm{Pr}(𝜽R|𝐘,𝐈)=\zeta `$ (e.g., a 67%, 90%, or 95% credible region). The 90% credible region for the posterior distributions illustrated in Figure 1 are (0.77, 4.24) using the informative prior and (0.04, 3.84) using the non-informative prior. Such probability statements are measures of our information regarding the value of the parameter $`𝜽`$, given the data and prior information. This is in contrast to the more traditional frequentist definition of probability, which defines a probability to be the long term frequency of an event generally involving the data given $`𝜽`$. The posterior distribution is a complete summary of our information, but is often summarized by its mean, $`\widehat{𝜽}=\mathrm{E}(𝜽|𝐘,𝐈)`$ and variance, $`\mathrm{Var}(𝜽|𝐘,𝐈)`$, or its modes and the curvatures at these modes. (The curvatures are most useful when the posterior is (locally) approximately Gaussian, as is asymptotically true under certain regularity conditions; e.g., Gelman et al. (1995) chapter 4.) In the following two sections we describe Monte Carlo methods for computing posterior means and variances and credible regions. To compute posterior modes (e.g., maximum likelihood estimates) van Dyk (2000) develops several Expectation Maximization or EM algorithms for use in astrophysical applications. Posterior modes are often used to compute starting values for the more robust but computationally demanding Monte Carl methods; see Appendix A.2. The EM algorithm gets its name because it iteratively maximizes the expected log posterior distribution of $`𝜽`$ given the augmented data.

Although a detailed description is beyond the scope of this paper, Bayesian methodology is well equipped for problems involving model selection. Methods based on Bayes’s factors, computing the relative posterior probabilities of various competing models, and Bayesian ‘p-values’ are all important and remain areas of active statistical research (e.g., Gregory & Loredo, 1992; van Dyk & Protassov, 2000).

### 2.2 Evaluating the Posterior via Monte-Carlo Sampling

For univariate or small dimensional parameter spaces, we can usually compute the posterior mean, variance, credible regions, and other summaries either analytically or via non-stochastic numerical methods (e.g., Gaussian quadrature or Laplace’s method). In higher dimensions, however, these methods can be difficult to implement partially because of the difficulty in finding the region where the integrand is significantly greater than zero. Thus, we often resort to Monte-Carlo integration. In particular, if we can obtain a sample from the posterior, $`\{𝜽_{[t]},t=1,\mathrm{},T\}`$, Monte Carlo integration approximates the mean of any function, $`g`$, of the parameter with

$$\mathrm{E}[g(𝜽)|𝐘,𝐈]\frac{1}{T}\underset{t=1}{\overset{T}{}}g(𝜽_{[t]}),$$
(7)
where we assume $`\mathrm{E}[g(𝜽)|𝐘,𝐈]`$ exists. For example, $`g(𝜽)=𝜽`$ and $`g(𝜽)=(𝜽\mathrm{E}(𝜽|𝐘,𝐈))(𝜽\mathrm{E}(𝜽|𝐘,𝐈))^{}`$ lead to the posterior mean and variance respectively. Probabilities, such as $`\zeta =\mathrm{Pr}(𝜽R)`$ can be computed using $`g(𝜽)=I\{𝜽R\}`$, where the function $`I`$ takes on value 1 if the condition in curly brackets holds and zero otherwise. Likewise, quantiles of the distribution can be approximated by the corresponding quantiles of the posterior sample. In short, a robust data analysis requires only a sample from the posterior distribution. A general strategy is to first sample from the posterior distribution and then approximate various integrals of interest via Monte Carlo.

### 2.3 Obtaining a Sample from the Posterior

The Monte-Carlo approximation methods depend on our ability to obtain a sample from the posterior distribution. Although in some cases the posterior distribution is a well known distribution and trivial to sample, we must often use sophisticated algorithms to obtain a posterior sample. In Appendix A, we discuss three algorithms which have proven widely applicable in practice, the Data Augmentation algorithm (Tanner & Wong, 1987), the Gibbs sampler (Metropolis et al., 1953), and the Metropolis-Hastings algorithm (Hastings, 1970). All of these algorithms construct a Markov chain with stationary distribution equal to the posterior distribution (e.g., Gelfand & Smith, 1990); i.e., once the chain has reached stationarity, it generates samples which are identically (but not independently) distributed according to the posterior distribution. These samples can then be used for Monte Carlo integration as described above; hence these algorithms are known as Markov chain Monte Carlo or MCMC methods (see Tierney, 1996, for regularity conditions for using Equation 7 with MCMC draws). From the onset then, it is clear that three important concerns when using MCMC in practice are (1) selecting starting values for the Markov chain, (2) detecting convergence of the Markov chain to stationarity, and (3) the effect of the lack of independence in the posterior draws. These issues are addressed in Appendix A.2.

The algorithms used to fit the models described in Section 3 rely on the method of data augmentation. The term ‘data augmentation’ originated with computational methods designed to handle missing data, but the method is really quite general and often useful when there is no missing data per se. In particular, for Monte Carlo integration we aim to obtain a sample from the posterior distribution, $`p(𝜽|𝐘,𝐈)`$. In some cases, we can augment the model to $`p(𝜽,𝐗|𝐘,𝐈)`$, where $`𝐗`$ may be missing data or any other unobserved quantity (e.g., counts due to background). With judicial choice of $`𝐗`$, it may be much easier to obtain a sample from $`p(𝜽,𝐗|𝐘,𝐈)`$ than directly from $`p(𝜽|𝐘,𝐈)`$. Once we have a sample from $`p(𝜽,𝐗|𝐘,𝐈)`$, we simply discard the sample of $`𝐗`$ to obtain a sample from $`p(𝜽|𝐘,𝐈)`$. In Appendix B.1, we describe how this method can be used for fitting the models described in Section 3.

## 3 Fitting High-Resolution Low-Count Spectra

### 3.1 Model Overview

In this section, we describe a new class of (statistical) structured models which simultaneously describes high-resolution source spectra using Gaussian line profiles and a GLM for the continuum and accounts for background contamination of the image, instrument response, and absorption. The model may easily be generalized to account for different line profiles such as the Lorentzian distribution (e.g., Meng & van Dyk, 1999). The (statistical) model is designed to summarize the distribution of photon energies arriving at a detector, which are recorded as counts in a number of energy channels (e.g., as many as 4096 on Chandra/ACIS). Newly developed detectors have much higher resolution than their predecessors, and thus smaller expected counts per bin. Independent Poisson distributions are therefore more appropriate for the counts than the commonly used Gaussian approximation (e.g., $`\chi ^2`$ fitting). We parameterize the intensity in bin $`j𝒥=\{1,\mathrm{},J\}`$, as the sum of a continuum term and $`K`$ Gaussian lines. That is, the expected true counts per bin for a “perfect” instrument with effective area everywhere equal to the maximum possible effective area999We use the maximum value of the effective area over the spectral energy range of interest in this stage of the analysis. This is only a matter of convenience and the full effective area variations are included in Equation 9. is

$$\begin{array}{ccccc}\hfill \text{model intensity}=[& \hfill \text{continuum}+& \hfill \text{lines}]& \text{absorption,}& \text{for each energy bin}\\ \multicolumn{4}{c}{}& \\ \hfill \lambda _j(𝜽)=[& \hfill dE_jf(𝜽^C,E_j)+& \hfill \underset{k=1}{\overset{K}{}}\stackrel{~}{\lambda }^kp_j(\mu ^k,\upsilon ^k)]& u(𝜽^A,E_j),& \text{for}j𝒥,\end{array}$$
(8)
or more formally,

where $`dE_j`$ is the known width of bin $`j`$, $`f(𝜽^C,E_j)`$ is the expected number of counts per keV per maximum effective area from the continuum and is a function of the continuum parameter, $`𝜽^C`$, $`E_j`$ is the known mean energy in bin $`j`$, $`\stackrel{~}{\lambda }^k`$ are the expected counts per maximum effective area from line $`k`$, $`p_j(\mu ^k,\upsilon ^k)`$ is the probability that a Gaussian random variable with mean $`\mu ^k`$ and variance $`\upsilon ^k`$ falls in bin $`j`$, and $`u(𝜽^A,E_j)`$ is the probability that a photon in bin $`j`$ is not absorbed. Specific forms for the continuum and absorption terms are discussed below in Sections 3.2 and 3.4, respectively. The superscripts on the model parameters ($`𝜽`$) are mnemonic and represent absorption (A), background (B), continuum (C), and the lines $`k=1,\mathrm{},K`$. The collection of parameters, $`𝜽^C,𝜽^k=(\stackrel{~}{\lambda }^k,\mu ^k,\upsilon ^k)\mathrm{for}k𝒦=\{1,\mathrm{},K\}`$, and $`𝜽^A`$ (along with $`𝜽^B`$ defined below) are represented by $`𝜽`$. An artificial example, with power law continuum, two spectral lines, and no absorption appears in the first plot of Figure 2.

Since data collection is degraded by effective area, instrument response, and background contamination (see Figure 2), we model the observed counts as independent Poisson variables with intensity

$$\begin{array}{cccc}\genfrac{}{}{0pt}{}{\text{observed}}{\text{intensity}=}\hfill & \hfill \genfrac{}{}{0pt}{}{\text{instrument}}{\text{response}}\left(\genfrac{}{}{0pt}{}{\text{model}}{\text{intensity}}\genfrac{}{}{0pt}{}{\text{effective}}{\text{area}}\right)+& \text{background,}& \genfrac{}{}{0pt}{}{\text{for each}}{\text{energy channel}}\hfill \\ \multicolumn{4}{c}{}\\ \hfill \xi _l(𝜽)=& \hfill \underset{j𝒥}{}M_{lj}\lambda _j(𝜽)d_j+& \lambda _l^B(𝜽^B),& \text{for}l,\end{array}$$
(9)
or more formally,

$`=\{1,\mathrm{},L\}`$ where the $`L\times J`$ matrix $`𝐌=\{M_{lj}\}`$ represents instrument response — a photon arriving in bin $`j`$ has probability $`M_{lj}`$ of being detected in observed bin $`l`$, $`𝐝=(d_1,\mathrm{},d_J)`$ is the effective area of bin $`j`$, normalized so that $`\mathrm{max}_{j𝒥}d_j=1`$, and $`\lambda _l^B(𝜽^B)`$ is the expected counts due to the background which may be known from calibration in space or parameterized in terms of $`𝜽^B`$. As with $`𝒥`$, $``$ may be any subset of detector bins. Generally the counts are also degraded by pile up (e.g., Knoll, 1989); see Figure 2. Here we ignore pile-up which is justifiable for low intensity or spatially diffuse sources (see the discussion in Section 5).

In the next several sections, we describe the stochastic models for each of the sources of photons in turn. This includes both likelihoods which describe the sampling distribution of the data (parameterized by $`𝜽`$) and prior distributions which allow us to incorporate scientific information about the likely parameter values. As described below, the prior distributions are parameterized using the hyper-parameter $`\mathit{\varphi }`$.

### 3.2 The Continuum

The photon counts due to the continuum are modeled via a GLM (McCullagh & Nelder, 1989), specifically a log linear model. That is, the log expected counts per keV per maximum effective area are assumed to be a linear function of a set of independent variables, $`X_j^C`$ which in turn are typically functions of $`E_j`$; hence the notation, $`f(𝜽^C,E_j)`$. In particular we model the counts in bin $`j`$ due to the continuum, denoted $`Y_j^C`$, as

$$Y_j^C\stackrel{\mathrm{d}}{}\mathrm{Poisson}(\mathrm{dE}_\mathrm{j}\mathrm{f}(𝜽^\mathrm{C},\mathrm{E}_\mathrm{j})),$$
(10)
i.e.101010Here and in the remainder of the paper we suppress the conditioning on the initial information, $`𝐈`$. That is, it should be understood that all distributions implicitly condition on $`𝐈`$.,

$$p(Y_j^C|𝜽^C)=e^{\lambda _j^C}(\lambda _j^C)^{Y_j^C}/(Y_j^C)!\mathrm{with}\lambda _j^C=f(𝜽^C,E_j)dE_j.$$
(11)
Here $`\mathrm{log}f(𝜽^C,E_j)=𝐗_j^C𝜽^C,`$ independently for $`j𝒥,`$ with $`𝜽^C`$ a $`(P^C\times 1)`$ vector parameter, $`𝐗_j^C`$ a $`(1\times P^C)`$ vector of independent variables, and $`P^C`$ the number of parameters in the continuum model. Note that we are explicitly using a Poisson process for the photon counts as opposed to an often poor Gaussian approximation.

The flexible framework of the GLM allows us to adjust the expected counts in bin $`j`$ for any set of independent variables. For example, several standard continuum models are easily available. In particular, a power law model is obtained by setting $`𝐗_j=(1,\mathrm{log}(E_j))`$ for $`j𝒥`$ so that

$$f(𝜽^C,E_j)=e^{\theta _1^C}E_j^{\theta _2^C}=\alpha E_j^\beta \mathrm{for}j𝒥,$$
(12)
where the familiar form of the power law model in the last expression is obtained by identifying $`(\alpha ,\beta )`$ with $`(e^{\theta _1^C},\theta _2^C)`$. It is easy to generalize this to handle more complicated models. A break in the power law (i.e., a change point) can be added at $`E_{}`$ by setting $`𝐗_j=(1,\mathrm{log}(E_j),\mathrm{log}(E_j/E_{})I\{E_j>E_{}\})`$, so that

$$f(𝜽^C,E_j)=\{\begin{array}{cc}e^{\theta _1^C}E_j^{\theta _2^C}\hfill & \text{for }E_jE_{}\hfill \\ e^{\theta _1^C}E_j^{\theta _2^C+\theta _3^C}E_{}^{\theta _3^C}\hfill & \text{for }E_j>E_{}\hfill \end{array}\mathrm{for}j𝒥.$$
(13)
The factor $`E_{}^{\theta _3^C}`$ ensures $`f(𝜽^C,E_j)`$ is continuous at $`E_j=E_{}`$. As a final example, we obtain an exponential continuum representing Bremsstrahlung emission by setting $`𝐗_j^C=(1,E_j)`$ so

$$f(𝜽^C,E_j)=e^{\theta _1^C}e^{E_j\theta _2^C}=\frac{\alpha }{\sqrt{T}}e^{E_j/kT}\mathrm{for}j𝒥,$$
(14)
where $`(\alpha ,T)=(e^{\theta _1^C}\sqrt{k\theta _2^C},1/k\theta _2^C)`$.

It is convenient to assume the prior distribution on $`𝜽^C`$ to is multivariate Gaussian with a diagonal variance matrix. That is, $`\theta _p^C\stackrel{\mathrm{d}}{}\mathrm{N}(\mu _\mathrm{p}^\mathrm{C},\upsilon _\mathrm{p}^\mathrm{C})`$ with $`\mathit{\varphi }^C=\{(\mu _p^C,\upsilon _p^C)`$, for $`p=1,\mathrm{},P^C\}`$. The hyper-parameter, $`\mathit{\varphi }^C`$, is set by the user where $`\mu _p^C`$ is a “best guess” of $`\theta _p^C`$ and $`\upsilon _p^C`$ is a measure (in squared standard deviations) of the error of this “best guess.” Large values of $`\upsilon _p^C`$ reflect little prior information for $`\theta _p^C`$.

### 3.3 Emission Lines

Lines reflect deviation in the smooth spectrum due to the continuum because of photon emissions from various ions present in the source. In particular we model the energies of photons due to line $`k𝒦`$, denoted $`Y_i^k`$ as

$$Y_i^k\stackrel{\mathrm{d}}{}\mathrm{N}(\mu ^\mathrm{k},\upsilon ^\mathrm{k})$$
(15)
i.e.,

$$p(Y_i^k|\mu ^k,\upsilon ^k)=\frac{1}{\sqrt{2\pi \upsilon ^k}}e^{(Y_i^k\mu ^k)^2/2\upsilon ^k}$$
(16)
independently for $`i=1,\mathrm{},N^k`$. Equation (15) represents a line with intensity normalized to one. The total line counts for a perfect instrument (i.e., with effective area everywhere equal to its maximum possible value) are denoted $`(N^1,\mathrm{},N^K)`$ and assumed to be independent Poisson random variables,

$$N^k\stackrel{\mathrm{d}}{}\mathrm{Poisson}(\stackrel{~}{\lambda }^\mathrm{k})\mathrm{independently}\mathrm{for}\mathrm{k}𝒦.$$
(17)
Proper prior information for the lines and the continuum is important for a reasonable fit when the spectral model includes emission lines. In particular, prior information is especially important for relatively weak lines, since it is difficult to distinguish a weak line from a chance fluctuation in the continuum. Luckily such prior information is often scientifically forthcoming in the form of knowledge (e.g., laboratory measurements and physics theory) of probable sizes and locations of the various lines. We begin with the line location and width (actually the variance), $`(\mu ^k,\upsilon ^k)`$ for which priors111111We choose this prior distribution partially because it is the so-called conjugate prior distribution, i.e., the resulting posterior distribution is from the same family as the prior distribution (e.g., Gaussian with updated parameters). This property significantly simplifies model fitting with no cost in terms of the accuracy of parameter estimation. are assigned independently for each line $`k𝒦`$,

$$\upsilon ^k\stackrel{\mathrm{d}}{}\frac{\nu _0^\mathrm{k}\upsilon _0^\mathrm{k}}{\chi _{\nu _0}^2}\mathrm{and}\mathrm{p}\left(\mu ^\mathrm{k}|\upsilon ^\mathrm{k}\right)=\mathrm{N}(\mu _0^\mathrm{k},\frac{\upsilon ^\mathrm{k}}{\kappa _0^\mathrm{k}}),$$
(18)
where $`\chi _{\nu _0}^2`$ is a variable that follows the $`\chi ^2`$ distribution with $`\nu _0`$ degrees of freedom. We interpret the hyper-parameter, $`\mathit{\varphi }^k=(\mu _0^k,\upsilon _0^k,\kappa _0^k,\nu _0^k)`$ using the mean and variance of the distributions in Equation 18. For example,

$$\mathrm{E}(\upsilon ^k|\mathit{\varphi }^k)=\nu _0^k\upsilon _0^k/(\nu _0^k2)\mathrm{for}\nu _0^2>2$$
(19)
and

$$\mathrm{Var}(\upsilon ^k|\mathit{\varphi }^k)=2(\nu _0^k\upsilon _0^k)^2(\nu _0^k2)^2(\nu _0^k4)\mathrm{for}\nu _0^2>4.$$
(20)
(Recall that the units here are keV for means and keV2 for variances.) Thus, the mean and variance of the prior for $`\upsilon ^k`$ may be tuned using $`\upsilon _0^k`$ and $`\nu _0^k`$; a small value of $`\nu _0^k`$ results in a wide, relatively non-informative prior. Since the data are discrete, a priori we cannot allow the standard deviation of the line to become too small (say below the PHA bin width of the bin that contains the center of $`\mu ^k`$) since there is not information in the data about the width of a line which is narrower than one PHA channel. This is accomplished by truncating the prior distribution of $`\upsilon ^k`$. For the prior on $`\mu ^k`$, the mean and variance are given by $`\mu _0^k`$ and $`\kappa _0^k`$; $`\mu _0^k`$ is the most probable location of the $`k`$th line and $`\kappa _0^k`$ calibrates the uncertainty in the location of the $`k`$th line relative to the width of the line.

An alternative interpretation of the priors is in terms of additional hypothetical photons. Heuristically, the effect of the prior on $`\mu ^k`$ if $`\upsilon ^k`$ were known would be the same as $`\kappa _0^k`$ photons all known to be from line $`k`$ and equal to $`\mu _0^k`$. Likewise, the effect of the prior on $`\upsilon ^k`$ is the same as adding $`\nu _0^k`$ photons with average squared deviation from the mean equal to $`\upsilon _0^k`$.

We now turn to the prior distribution on $`\stackrel{~}{\lambda }^k`$ and set $`\stackrel{~}{\lambda }^k\stackrel{\mathrm{d}}{}\gamma (\varphi _1^\mathrm{k},\varphi _2^{\stackrel{~}{\lambda }})`$,121212We choose a $`\gamma `$ prior partially because it is conjugate to the Poisson distribution. (independently) which has mean $`\varphi _1^k/\varphi _2^{\stackrel{~}{\lambda }}`$ and variance $`\varphi _1^k/(\varphi _2^{\stackrel{~}{\lambda }})^2`$. Roughly speaking, the $`\gamma `$ prior contains the same information as $`\varphi _2^{\stackrel{~}{\lambda }}`$ Poisson observations (with exposure equal to the source exposure) with a total of $`\varphi _1^k1`$ counts. Since the data consist of a single observation (for each bin), $`\varphi _2^{\stackrel{~}{\lambda }}`$ can be interpreted as the weight put on the prior relative to the data; $`\varphi _2^{\stackrel{~}{\lambda }}=1`$ induces a prior as influential as the data in the absence of absorption, blurring, background, and lines. Thus, values of $`\varphi _2^{\stackrel{~}{\lambda }}1`$ are typically recommended for non informative priors. The hyperparameters $`\varphi _1^k`$, can be interpreted as the prior relative sizes of the lines. That is, $`\varphi _1^k/_{\stackrel{~}{k}}\varphi _1^{\stackrel{~}{k}}`$ is the prior proportion of line photons from line $`k`$. We define the hyper-parameter for line $`k`$ as $`\mathit{\varphi }^k=(\mu _0^k,\kappa _0^k,\nu ^k,(\sigma _0^2)^k,\varphi _1^k,\varphi _2^{\stackrel{~}{\lambda }})`$; the last element is not indexed by $`k`$ since it is constant for $`k𝒦`$.

### 3.4 Absorption and Correction for Effective Area

From the viewpoint of our statistical algorithm, both the telescope effective area and astrophysical absorption (e.g., absorption due to the ISM) are handled in the same way. These two processes act independently on individual photons and randomly prevent a (energy dependent) proportion of photons from being observed. The only essential statistical difference being that only the absorption process has unknown parameters. In particular, we suppose that the probability a photon is not absorbed (statistically speaking “censored”) by either of these two processes is

$$d_ju(𝜽^A,E_j),\mathrm{where}\mathrm{log}u(𝜽^A,E_j)=𝐗_j^A𝜽^A\mathrm{for}j𝒥,$$
(21)
where $`𝜽^A`$ is a $`(P^A\times 1)`$ parameter, $`𝐗_j^A`$ is a $`(1\times P^A)`$ vector of independent variables, and $`P^A`$ is the number of parameters in the absorption model, $`u(𝜽^A,E_j)`$. As an example, simple exponential absorption can be written in this linear form with $`\theta ^A`$ a scalar and $`𝐗_j^A=1/\mathrm{E}_j`$, i.e., $`u(\theta ^A,E_j)=e^{\theta ^A/E_j}`$. For more complicated absorption models, $`𝐗_j^A`$ typically consists of a tabulated absorption function.

The prior for $`𝜽^A`$ is multivariate Gaussian, $`\theta _p^A\stackrel{\mathrm{d}}{}\mathrm{N}(\mu _\mathrm{p}^\mathrm{A},\upsilon _\mathrm{p}^\mathrm{A})`$, independently for $`p=1,\mathrm{},P^A`$. The prior is interpreted similarly to that for the continuum parameter $`𝜽^C`$. We, however truncate this prior to ensure $`\mathrm{exp}\{𝐗_j^A𝜽^A\}<1`$ for each $`j`$, to ensure the proportion of photons not absorbed is less than 1. With appropriately chosen $`𝐗_j^A`$, this can be accomplished by assuming each component of $`𝜽^A`$ is negative.

### 3.5 Background

We assume the availability of a separate observation containing background counts which can be used to model the background spectrum and correct the source spectrum. Rather than simply subtracting off (a scalar multiple of) the background counts, however, we account for the variation due to the Poisson character of the counts. In particular, we suppose the background count in PHA channel $`l`$ is

$$Y_l^B\stackrel{\mathrm{d}}{}\mathrm{Poisson}(\theta _\mathrm{l}^\mathrm{B})\mathrm{independently}\mathrm{for}\mathrm{l},$$
(22)
where the unobserved quantity, $`Y_l^B`$ is the counts in PHA channel $`l`$ that are due to background. We parameterize the prior for $`\theta _l^B`$ as $`\gamma (\varphi _{l,1}^B,\varphi _{l,2}^B)`$, which we expect to be informative based on a pure background exposure. In particular, a reasonable prior based on $`Y_l^{\mathrm{obs},B}`$ background counts would be $`\gamma (Y_l^{\mathrm{obs},B}+1,\tau )`$, where $`\tau `$ is the background exposure time and area relative to the source exposure time and area. If $`Y_l^{\mathrm{obs},B}=0`$, we generally replace it by some small number, such as 0.25 in the prior. (This value should not be too small, since this would preclude the possibility of background counts in the corresponding channel in the same exposure.)

In an extreme case, when the background is very well determined, e.g., via a very long exposure, we may fix $`\theta _l^B=Y_l^{\mathrm{obs},B}/\tau `$ and discard the prior distribution; here we are effectively setting the prior variance to zero. Note this is not equivalent to subtracting off the background because we still allow for Poisson variability in the background counts that contaminate the source counts. An alternative strategy is to fit a parameterized model to the background. For example, we might assume $`Y_l^B\stackrel{\mathrm{d}}{}\mathrm{Poisson}(\lambda _\mathrm{l}^\mathrm{B}(𝜽^\mathrm{B},\stackrel{~}{\mathrm{E}}_\mathrm{l}))`$, where $`\mathrm{log}\lambda _l^B(𝜽^B,\stackrel{~}{E}_l)=𝐗_l^B𝜽^B`$ for $`l`$ with $`𝐗_l^B`$ a row vector of independent variables depending on the energy of PHA channel $`l,\stackrel{~}{E}_l`$. This allows the background counts to be modeled as a power law, broken power law, or any other log linear model.

## 4 Applications

In this section, we illustrate our methods and algorithms using two datasets. We first analyze ASCA/SIS data of high redshift (z=3.384) quasar S5 0014+813 (Elvis et al., 1994) to illustrate the various summaries available in a relatively straight-forward MCMC analysis. The second analysis involves an extremely low-count stellar coronal source ($`\alpha `$ TrA) and illustrates the power of Bayesian methods to combine information from various sources and quantify the weak information available in this data.

### 4.1 Quasar S5 0014+813

A typical quasar X-ray spectrum can be described by an absorbed power law. A fluorescent iron line (Fe K-$`\alpha `$) emitted at energy between $`6.46.8`$ keV if detected can be a signature of a reflection component and its ionization state (George et al., 2000). Quasar S5 0014+813 (Kühr et al., 1981) at redshift z=3.384 is among the highest X-ray flux quasars known with z $``$ 3. S5 0014+813 was observed with ASCA on 1993 October 29 with an exposure time of 22.8 ks in the SIS0 detector (Elvis et al., 1994). Here we apply our model to this data to illustrate the method and look for signatures of the iron emission line.

The spectral data were extracted with the standard screening criteria (Elvis et al., 1994) and standard response matrices were used (ftp://legacy.gsfc.nasa.gov/caldb/data/asca/sis/cpf/94nov9). We use all of the original 512 PHA instrument channels except the unreliable channels below $`0.5`$ keV and above $`10`$keV. In addition we do not group any channels. (Channels are usually grouped in order to justify the use the default $`\chi ^2`$ techniques with their Gaussian assumptions.) As is allowed with a Poisson model, we instead use only the original PHA bins. The source model included the exponential shape of Galactic absorption (see Section 3.5), and a power law continuum (i.e., Equation 12) with a narrow emission line at 1.45 keV (observed frame; $`6.7`$keV rest frame). We accounted for background using a background Poisson process with intensity equal to the (rescaled) background counts in each PHA channel. Flat priors were used on all model parameters.

To estimate the four model parameters (i.e., the power law, normalization, and exponential absorption parameters and the equivalent width131313The equivalent width is defined as $`\stackrel{~}{\lambda }^k/f(\theta ^C,\mu ^k)`$. of the line) a sample from their posterior distribution was obtained by running three MCMC chains using dispersed starting values. The chains showed excellent mixing (as measured with $`\sqrt{\widehat{R}}`$, see Appendix A.2) after 2000 draws. In the Monte Carlo evaluation, the second half of each of the chains were used along with an additional run of 2000 draws from each chain, for a total of 9000 draws.

Summaries of the model fit appear in Table 2, Figure 3, and Figure 4. The parameter estimates are posterior means computed using a transformation which makes the marginal posterior distributions more symmetric and hence the posterior mean a more informative summary (i.e., ln(Normalization) and sqrt(Equivalent Width)). In particular, if we represent the draws of the normalization parameter as $`\{\theta _{[t]},t=1,\mathrm{},9000\}`$, the point estimate of this parameter was computed as

$$\mathrm{exp}\left\{\frac{1}{9000}\underset{t=1}{\overset{9000}{}}\mathrm{ln}(\theta _{[t]})\right\}=\sqrt[9000]{\underset{t=1}{\overset{9000}{}}\theta _{[t]}},$$
(23)
the geometric mean. The credible intervals are computed using the 0.025 and 0.975 quantiles of the draws and are invariant to (monotonic) transformations. Pairwise credible regions appear in Figure 3. The scatter plots illustrate the regions of highest posterior probability by plotting the Monte Carlo draws — $`\mathrm{Pr}(𝜽R)`$ is approximately equal to the proportion of points in that region. The grey-scale images give Monte Carlo estimates of the (darker) 50% and (lighter) 90% marginal posterior regions. The grainy character of the images is due to the Monte Carlo approximation. Even with this relatively large data set and with the use of transformations, the non-Gaussian character of the posterior is evident. We expect that higher dimensional marginal posterior distributions are even less Gaussian in character. Figure 4 compares the fitted source model corrected for effective area and absorption with the PHA counts and illustrates the estimated continuum and the stability of this estimate.

### 4.2 Hybrid-Chromosphere Supergiant Star $`\alpha `$ TrA

Unlike the simple power-law spectrum of the quasar in the previous section, stellar coronal spectra are complicated by a Bremsstrahlung continuum and the presence of numerous emission lines. Such complex spectra are much more difficult to model, and in addition, the intensity of the Bremsstrahlung continuum drops exponentially at high-energies, resulting in very few counts. Analyzing such spectra is however crucial to the understanding of coronal structure, mechanisms of coronal heating, etc. A case in point is the corona of the Hybrid supergiant star $`\alpha `$ TrA (HD 150798, K4II, B-V=1.44, V=1m.92), which shows evidence of both strong magnetic activity as indicated by X-ray emission (Brown et al., 1991; Kashyap et al., 1994) and stellar outflow seen in absorption profiles (Hartmann et al., 1981). X-ray observations with the ROSAT/PSPC (Kashyap et al., 1994) indicate that its corona is dominated by transient, unstable plasma that is confined by magnetic loops that are closed on short length scales (Rosner et al., 1994). Constraining the maximum temperatures present in the corona is therefore of primary importance. Here we use data obtained with ASCA, at higher energies than ROSAT, to model the spectrum. The low number of counts detected at high energies make this spectrum difficult to analyze by traditional means, and we must bring to bear the full power of a hierarchical Bayesian analysis in order to constrain the maximum temperatures present in the corona.

$`\alpha `$ TrA was observed with ASCA in March 1995 for $``$ 34 ks. During this observation, the source exhibited no flares. The count rate was steady, and corresponded to the quiescent state identified with ROSAT.

We model the high-energy region of the ASCA spectrum (2.5 – 7.5 keV) as a combination of a Bremsstrahlung continuum

$$\frac{Norm}{\sqrt{T}}e^{E/k_BT},$$
(24)
where $`T`$ is the electron temperature, $`E`$ is the energy in keV, and $`k_B`$ is the Boltzmann constant, and $`Norm`$ is a normalization; and a number ($`10`$) of narrow emission lines located at the positions of known strong lines whose widths and locations141414Line location can be known only to the resolution of the instrument, and hence each of the model lines represents the sum of a large number of lines within the resolution element; we find that we only exclude $`<5\%`$ of the line flux by this approximation in the energy range considered. are fixed, but intensities are allowed to vary. As is our convention, the units of $`Norm`$ is counts per keV per maximum effective area. Because of the low counts we fit a power law to the background.

We apply the above model to SIS0 data ($`28`$ ks) and the combined data from the 2 GIS detectors ($`33`$ ks in each). Note that the very low counts present in these data ($`150`$ counts in SIS0, $`300`$ in GIS) preclude any “traditional” analysis: it is only by using the full Bayesian machinery that we can derive useful results from such data.

We carry out the analysis in two steps:

1. Choose highly non-informative priors on the parameters to analyze SIS0 data: for the prior on normalization, $`p\left(\frac{Norm}{\sqrt{T}}\right)`$, such that it lies in between $`10^9`$ and $`10^1`$ with a 90% probability; on temperature, $`p\left(\frac{1}{kT}\right)`$ such that it is always positive and also is nearly flat in the temperature range of interest; and on line intensity, $`p\left(\stackrel{~}{\lambda }^k\right)`$ such that a priori all the lines have the same intensity, and the maximum total counts due to lines is 100 (the total source+background counts for the SIS0 observation is only 154; atomic emission line models indicate that for the temperature and energy range of interest, 100 corresponds to the maximum possible contribution to the spectrum from lines). We choose Gaussian forms for the first two and a Gamma prior for the last; these priors are illustrated as solid lines in Figure 5. Thus,

$`p\left(\mathrm{ln}\left({\displaystyle \frac{Norm}{\sqrt{T}}}\right)\right)`$ $`=`$ $`N(\mu =9.21,\sigma =5.58)`$ (25)
$`p\left({\displaystyle \frac{1}{kT}}\right)`$ $`=`$ $`N(\mu =0.95,\sigma =2.5)`$ (26)
$`p\left(\stackrel{~}{\lambda }^k\right)`$ $`=`$ $`\gamma (\varphi _1^k=0.11,\varphi _2^{\stackrel{~}{\lambda }}=0.0033)\mathrm{for}k=1,\mathrm{}10.`$ (27)
2. Use the posterior resulting from the above step to define more informative priors to analyze GIS data. These priors also correct for the difference in exposure time and average effective area between the SIS0 and the GIS data. The posterior variances from the initial analysis were increased somewhat when computing the priors for the second analysis. These priors are illustrated as solid lines in Figure 6. Thus,

$`p\left(\mathrm{ln}\left({\displaystyle \frac{Norm}{\sqrt{T}}}\right)\right)`$ $`=`$ $`N(\mu =9.97,\sigma =1.74)`$ (28)
$`p\left({\displaystyle \frac{1}{kT}}\right)`$ $`=`$ $`N(\mu =0.41,\sigma =0.43)`$ (29)
$`p\left(\stackrel{~}{\lambda }^k\right)`$ $`=`$ $`\gamma (\varphi _1^k=0.12,\varphi _2^{\stackrel{~}{\lambda }}=0.025)\mathrm{for}k=1,\mathrm{},K.`$ (30)

We ran three Markov Chains in each analysis to obtain draws from the posterior distribution of $`(\mathrm{ln}(Norm/\sqrt{T}),1/kT,\stackrel{~}{\lambda }^1,\mathrm{},\stackrel{~}{\lambda }^{10})`$. In both analyses there was excellent mixing after 6000 draws and we used the second half of each chain for a total of 9000 Monte Carlo draws.

The results of the analysis are shown in Figures 58. (In the figures the parameter $`\mathrm{\Omega }`$ refers to the proportion of source photons from the lines and $`\lambda =_{k=1}^{10}\stackrel{~}{\lambda }^k`$.) We find that the plasma temperature is $`19_{>11}^{<64}\times 10^6`$ K (Figure 5). Such a large value (cf. $`2\times 10^6`$ K in the quiet Solar corona) clearly lends credence to the idea that the corona on $`\alpha `$ TrA is dominated even in quiescence by flare-like events.

As a byproduct of our analysis, we also obtain the flux in the modeled lines relative to the continuum. In principle, this allows us to constrain the metallicity for the first time in the corona of $`\alpha `$ TrA by comparing the observed ratio of the line and continuum fluxes151515Incompleteness in atomic line databases ((see Brickhouse, 1998)) contribute to an error of $`<5\%`$ on the line-to-continuum ratios calculated here. They are negligible compared to the measurement error. with that derived from thermal emission models computed over the same temperature range (Drake & Kashyap, 1998; Kashyap et al., 1998). The photospheric metallicity (Taylor, 1999) is \[Fe/H\] = 0.3, and we derive for the coronal metallicity \[Fe/H\] = $`0.4_{0.6}^{+1.1}`$ where the quoted range represents posterior deviations at $`1\sigma `$. While the uncertainty on our measurement is quite large (it is essentially unbounded at high metallicity), it is encouraging that the corona does not appear to be metal abundance deficient (see Drake, 1996).

## 5 Discussion

The power of the Bayesian methods illustrated here lies in their ability to combine information and to directly model the highly structured hierarchical features of the data—both in a principled manner. These features are illustrated in the $`\alpha `$ TrA example. First, by combining information from several detectors, we are able to extract information from the data regarding the plasma temperature. More generally, Bayesian methods allow for the incorporation of various forms quantifiable prior information through the prior distribution. Of course, results are then conditional on the prior information—if these priors are not trusted, the conclusions cannot be trusted either. On the other hand, if the prior information is accepted as reasonable, the posterior distribution should be accepted as a conglomeration of prior scientific information and the data. Second, the extremely low counts in the $`\alpha `$ TrA data, along with many free parameters (ten emission line intensities and two continuum parameters) illustrate a situation in which methods based on the Gaussian distribution and the central limit theorem are simply without justification. Methods which account for the Poisson (e.g., highly variable) character of the data have a sound mathematical basis and, in contract to standard methods such as $`\chi ^2`$ fitting, are equipped to handle such data.

The hierarchy in the model described in Section 3 can be extended to account for various more complicated features in the data, e.g., absorption lines, pile-up, and joint spatial, spectral, and temporal structure. Dealing with pile-up is perhaps the most important outstanding data-analytic challenge for Chandra. Conceptually, however, there is no difficulty in addressing pile-up in a Bayesian framework. After accounting for other features in the data such as instrument response, background, and absorption, we simply need to separate the observed counts into multiple counts of lower or equal energy based on the (current draw of the) spectral and spatial model. The difficulty lies in computation. Simply enumerating the set of photons that could result in a particular observed event, let alone their relative probabilities is an enormous task. Thus, we believe there is great promise in Monte Carlo techniques which if carefully designed, can automatically exclude numerous possibilities with minute probability. Although there remains much work to be done, Bayesian methods in conjunction with MCMC algorithms offer a practical and innovative solution to many outstanding data-analytic challenges in astrophysics.

The authors gratefully acknowledge funding for this project partially provided by NSF grant DMS-97-05157, the Chandra X-ray Observatory, and by NASA grant NAS8-39073. They also thank Rostislav Protassov and C.J. Shen for their work on the programming of the various algorithms, Richard Edgar and Paul Gorenstein for helpful comments on an earlier draft, and the referee, Jeff Scargle, whose careful reading and detailed suggestions have greatly improved the paper. Appendices

## Appendix A Markov Chain Monte Carlo Methods

### A.1 The Data Augmentation Algorithm

The Data Augmentation algorithm is designed to obtain a sample from the posterior distribution for use in Monte Carlo integration. The strategy of the algorithm is to embed the posterior distribution, $`p(𝜽|𝐘)`$, into a distribution in a large space, $`p(𝜽,𝐘^{\mathrm{mis}}|𝐘)`$. If we can obtain a sample from this second distribution, we need only discard the sampled values of$`𝐘^{\mathrm{mis}}`$ to obtain the desired sample from the posterior. The quantity $`𝐘^{\mathrm{mis}}`$ can be any unobserved quantity; it is referred to as “missing data” for historical reasons. For clarity we denote the observed data $`𝐘^{\mathrm{obs}}`$ and the augmented data $`𝐘^{\mathrm{aug}}=(𝐘^{\mathrm{obs}},𝐘^{\mathrm{mis}})`$. In order to obtain a sample from $`p(𝜽,𝐘^{\mathrm{mis}}|𝐘^{\mathrm{obs}})`$, the Data Augmentation algorithm uses an iterative sampling scheme that samples first $`𝐘^{\mathrm{mis}}`$ conditional on the model parameters and $`𝐘^{\mathrm{obs}}`$ and second samples the model parameters given $`𝐘^{\mathrm{aug}}`$. Clearly, the algorithm is most useful when both of these conditional distributions are easily sampled from. The iterative character of the resulting chain naturally leads to a Markov chain, which we initialize at some starting value, $`𝜽_{[0]}`$. For $`t=1,\mathrm{},T`$, where $`T`$ is dynamically chosen, we repeat the following two steps:

Step 1: Draw $`𝐘_{[t]}^{\mathrm{aug}}`$ from $`p(𝐘^{\mathrm{aug}}|𝐘^{\mathrm{obs}},𝜽_{[t1]})`$,

Step 2: Draw $`𝜽_{[t]}`$ from $`p(𝜽|𝐘_{[t]}^{\mathrm{aug}})`$.

Under certain regularity conditions (see Meyn & Tweedie, 1993; Roberts, 1996; Tierney, 1994, 1996, for details) the stationary distribution of the resulting Markov chain is the desired posterior distribution, i.e., for large $`t`$, $`𝜽_{[t]}`$ approximately follows $`p(𝜽|𝐘^{\mathrm{obs}})`$.

To illustrate the utility of the Data Augmentation algorithm, we return to the simple background contamination model introduced in Section 2.1. The choice of $`𝐘^{\mathrm{aug}}`$ is clear in the example; we set $`𝐘^{\mathrm{aug}}=\{Y,Y^S,Y^B\}`$, where $`Y`$ is the total counts, $`Y^S`$ is the unobserved source counts from the source exposure, and $`Y^B`$ is the counts from the pure background observation. (I.e., we can consider $`Y^S`$ to be the missing data.) With this choice of $`𝐘^{\mathrm{aug}}`$, both $`p(𝐘^{\mathrm{aug}}|𝐘^{\mathrm{obs}},𝜽)`$ and $`p(𝜽|𝐘^{\mathrm{aug}})`$ are easy to sample and thus the Data Augmentation algorithm is easy to use; here $`𝐘^{\mathrm{obs}}=\{Y,Y^B\}`$ and $`𝜽=(\lambda ^B,\lambda ^S)`$. Given some $`𝜽_{[0]}=(\lambda _{[0]}^B,\lambda _{[0]}^S)`$ the two steps of the Data Augmentation algorithm at iteration $`t`$ become

Step 1: Draw $`Y_{[t]}^S`$ from

$$Y^S|𝜽_{[t1]},𝐘^{\mathrm{obs}}\stackrel{\mathrm{d}}{}\mathrm{binomial}(\mathrm{Y},\frac{\lambda ^\mathrm{S}}{\lambda ^\mathrm{B}+\lambda ^\mathrm{S}}),$$
(A1)
i.e.,

$$p\left(Y^S|𝜽_{[t1]},𝐘^{\mathrm{obs}}\right)=\left(\genfrac{}{}{0pt}{}{Y}{Y^S}\right)(Y^S)^{\lambda ^S/(\lambda ^B+\lambda ^S)}(YY^S)^{\lambda ^B/(\lambda ^B+\lambda ^S)}$$
(A2)
Step 2: Draw

$$\lambda _{[t]}^B|𝐘_{[t]}^{\mathrm{aug}}\stackrel{\mathrm{d}}{}\gamma (\alpha ^\mathrm{B}+\mathrm{Y}\mathrm{Y}^\mathrm{S},\beta ^\mathrm{B}+1),$$
(A3)
i.e.,

$$p\left(\lambda _{[t]}^B|𝐘_{[t]}^{\mathrm{aug}}\right)=\left(\lambda _{[t]}^B(\beta ^B+1)\right)^{(\alpha ^B+YY^S)}e^{\lambda _{[t]}^B(\beta ^B+1)}/(\lambda _{[t]}^B\mathrm{\Gamma }(\alpha ^B+YY^S)),$$
(A4)
where $`\alpha ^B`$ and $`\beta ^B`$ are typically chosen using the pure background observation as described in Section 2.1, and

$$\lambda _{[t]}^S|𝐘_{[t]}^{\mathrm{aug}}\stackrel{\mathrm{d}}{}\gamma (\alpha ^\mathrm{S}+\mathrm{Y}^\mathrm{S},\beta ^\mathrm{S}+1).$$
(A5)
In the first step, we stochastically divide the source count into source counts and background counts based on the current values of $`\lambda ^B`$ and $`\lambda ^S`$. In the second step we use this division to update $`\lambda ^B`$ and $`\lambda ^S`$. Markov chain theory tells us the iteration converges to the desired draws from the posterior distribution.

By selecting a starting value and iteratively sampling according to Equations A1, A3, and A5, we obtain a Markov chain which delivers a dependent sample from the posterior distribution upon convergence. In the next section, we use the Data Augmentation algorithm to illustrate the important practical issues of selecting starting values, detecting convergence, and accounting for the dependency in the sample.

### A.2 Starting Values, Convergence, and Multiple Chains

An important and difficult aspect of MCMC methods in practice is ascertaining convergence to stationarity. Since the stationary distribution of the Markov chain is the posterior distribution on interest, we can consider $`\{𝜽_{[t]},t>T_0\}`$ to be a (dependent) posterior sample, which can be used for Monte Carlo integration. Thus, determining $`𝜽_{[0]}`$ and $`T_0`$ is critical for valid inference. There is a large and growing literature on these related subjects and we refer interested readers to recent texts on the subject by Gelman et al. (1995), Carlin & Louis (1996), and Gilks et al. (1996), as well as the review article on convergence by Cowles & Carlin (1996). Here we briefly outline the approach that we find most fruitful.

As proposed by Gelman & Rubin (1992), we suggest running multiple Markov chains with a variety of starting values spread throughout the parameter space. This is a useful procedure since a single chain can appear to have converged when actually it has only settled temporarily in one region of the parameter space. This is illustrated with the Markov chain in Figure 9. The three chains show the draws of a variance parameter for a random-effects model (for details see van Dyk & Meng, 2000). Note that although chain 3 appears relatively stable it is far from convergence during the first 10,000 draws. This is evident when it is compared with the other chains, but less so when we look only at the beginning of chain 3. It is recommended that the starting values for the several chains be spread broadly in the parameter space (relative to the region of high posterior probability). This can often be accomplished by roughly mapping the posterior, for example, using estimates and errors based on the chi square estimates, posterior modes, or maximum likelihood estimates. (See van Dyk (2000) for details on the computation of posterior modes and maximum likelihood estimates for our spectral model.) Once such “over dispersed” starting values are obtained, we can run the several chains until all converge to the same region of the parameter space. (There may be more than one mode in the posterior, in which case the chains may converge to different modes, i.e., different regions of the parameter space.) Gelman & Rubin’s (1992) $`\sqrt{\widehat{R}}`$ statistic measures the relative size of the total variance in the draws of a univariate function of the parameter and the average within chain variance of the same function, i.e.,

$$\sqrt{\widehat{R}}=\sqrt{\frac{\frac{T1}{T}W+\frac{1}{T}B}{W}},$$
(A6)
where $`B`$ is the between chain variance, $`W`$ is the within chain variance, and $`T`$ is the number of draws. If the variance within each chain is as great as the total variance in all the draws, i.e., $`\sqrt{\widehat{R}}`$ is near one, then we can be confident that all the chains have converged to the same region of the parameter space. Typically we compute $`\sqrt{\widehat{R}}`$ using the last half (or two thirds) of each of the chains. Once an acceptable level of $`\sqrt{\widehat{R}}`$ is obtained (say below 1.2) we omit the first half (or third) of the chain in all further analysis. If we have several starting values which cover a large enough region of the parameter space, we can be confident that the chains sample all areas with high posterior probability and thus the Monte Carlo approximations are unbiased estimators of the quantities they estimate.

The variance of the Monte Carlo approximations is a function of the posterior variance of the quantity being approximated, the posterior sample size (i.e., $`TT_0`$), and the autocorrelation function of the Markov chain. Typically Monte Carlo errors are small relative to the posterior variance with several thousand posterior draws and thus are of little consequence. Monte Carlo error can be quantified by repeating the analysis for the first half and second half of the Markov chain and noting if the results are substantively different. See Roberts (1996) and the references therein for details and extensions.

### A.3 The Gibbs Sampler

In this and the next section we describe two additional MCMC methods which are designed to delivers a sample from the posterior distribution and are often useful when the Data Augmentation algorithm is not practical. The Gibbs sampler can be viewed as an extension of the Data Augmentation algorithm in which we wish to sample from $`p(𝜽|𝐘^{\mathrm{obs}})`$, and the vector, $`𝜽`$, can be viewed as a combination of model parameters and “missing data”. (In many instances, there is no “missing data”.) We partition $`𝜽`$ into $`(𝜽_1,\mathrm{},𝜽_P)`$, where $`𝜽_p`$ may be a scalar or vector quantity for each $`p`$. The Gibbs sampler, again starts with some starting value $`𝜽_{[0]}`$ and at iteration $`t`$ samples according to the following conditional distributions:

Step 1: Draw $`(𝜽_1)_{[t]}\stackrel{\mathrm{d}}{}\mathrm{p}(𝜽_1|(𝜽_1)_{[\mathrm{t}]},𝐘^{\mathrm{obs}})`$,

Step 2: Draw $`(𝜽_2)_{[t]}\stackrel{\mathrm{d}}{}\mathrm{p}(𝜽_2|(𝜽_2)_{[\mathrm{t}]},𝐘^{\mathrm{obs}})`$,

$`\mathrm{}`$

Step $`P`$: Draw $`(𝜽_P)_{[t]}\stackrel{\mathrm{d}}{}\mathrm{p}(𝜽_\mathrm{P}|(𝜽_\mathrm{P})_{[\mathrm{t}]},𝐘^{\mathrm{obs}})`$,

where $`(𝜽_p)_{[t]}=((𝜽_1)_{[t]},\mathrm{},(𝜽_{p1})_{[t]},(𝜽_{p+1})_{[t1]},\mathrm{},(𝜽_P)_{[t1]})`$. That is, we draw each component of $`𝜽`$ in turn conditional on the current values of the rest of $`𝜽`$ and the data.

The advantage of the Gibbs sampler over the Data Augmentation algorithm is that in many settings additional conditioning results in simpler draws. The disadvantage is that the resulting Markov chains tend to have higher autocorrelation and are slower to converge to stationarity as $`P`$, the number of steps per iteration, increases.

### A.4 The Metropolis-Hastings Algorithm

As a final extension, we consider the case when one (or more) of the steps in the Gibbs sample involves a conditional distribution that is not easy to sample. The Metropolis-Hastings algorithm (Metropolis & Ulam, 1949; Metropolis et al., 1953; Hastings, 1970) replaces the conditional distribution by some convenient “jumping rule” which approximates the conditional distribution. A proposal draw is sampled according to the jumping rule and is either accepted or rejected (in which case, the Markov chain is fixed at the previous draw) according to a rule that maintains the desired stationary distribution (see, e.g., Gelman et al., 1995, for details).

## Appendix B Details of the MCMC Algorithm

### B.1 Data Augmentation

The algorithms used to fit the model described in Section 3.1 rely on the method of data augmentation. In this section, we detail the layers of the data-augmentation scheme we use. We aim to construct an idealized data set for which model fitting is a relatively easy task. That is, given the augmented data, we can easily sample the model parameters. Likewise, given the model parameters, we can easily sample the augmented data and thus we can construct a Data Augmentation algorithm as described in Appendix A.1. Suppose, for example, a data set uncontaminated by background or instrument response were available. Clearly, model fitting would be easier. We define an even larger data set that contains the unbinned, true, and blurred energies of all photons that would have arrived at the detector if there had been no absorption and if we were using a perfect instrument with the effective area equal to its maximum value over all energies used. This data set also includes a variable indicating absorption and loss to reduced effective area161616 Absorption and effective area are handled together, so we need only one indicator variable, see Section 3.4. and a variable indicating the source of each photon, i.e., background ($`B`$), continuum ($`C`$), and each of the $`K`$ line profiles; the set of sources is denoted $`𝒮=\{B,C,1,2,\mathrm{},K\}`$. This idealized data set is summarized in Table 3.

The data augmentation scheme is illustrated in Figure 10 in which squares and circles represent observed and unobserved (“augmented”) quantities, respectively. Given the model parameter $`𝜽`$, we obtain a sample set of photon energies, $`𝐘{}_{}{}^{^{\mathrm{}}}{}_{}{}^{s}=(Y{}_{}{}^{^{\mathrm{}}}{}_{1}{}^{s},\mathrm{},Y{}_{}{}^{^{\mathrm{}}}{}_{N^s}{}^{s})^{}`$ for $`s𝒮`$ (see the third column of Figure 10) represented the undegraded “augmented” data; $`N^s`$ is the total count for source $`s`$. (As a mnemonic device, more dots in the accent above $`Y`$ signifies further removal of a quantity from actual observable quantities.) Here, $`𝐘{}_{}{}^{^{\mathrm{}}}^k`$ contains the exact energy of all photons attributed to line $`k`$ before absorption, with maximum effective area, and no background contamination. (The background photon energies, $`𝐘{}_{}{}^{^{\mathrm{}}}^B`$, do not appear in Figure 10 because we model the detected counts (e.g., in PHA channels) rather than true counts, see Section 3.6.) The first two columns of Figure 10 represent the hyperparameters and model parameters detailed in Sections 3.2–3.5.

The array of energies represented by $`𝐘{}_{}{}^{^{\mathrm{}}}^s`$ are binned into instrument-specific energy bins to obtain a sample spectrum, $`\ddot{𝐘}^s=(\ddot{Y}_1^s,\mathrm{},\ddot{Y}_J^s)^{}`$ (see the fourth column in Figure 10). In particular,

$$\ddot{Y}_j^s=\underset{i=1}{\overset{N^s}{}}I\{Y{}_{}{}^{^{\mathrm{}}}{}_{i}{}^{s}B_j\},\mathrm{for}j𝒥\mathrm{and}s𝒮,$$
(B1)
where $`B_j`$ is the $`j`$th energy bin. The first plot in Figure 2 illustrates the undegraded counts from the continuum and lines, $`\ddot{Y}_j^+=_{s𝒮B}\ddot{Y}_j^s`$ for the artificial data set, where the notation $`𝒮B`$ indicates set subtraction, i.e., the set $`𝒮`$ with $`B`$ removed. The solid line represents $`\mathrm{E}(\ddot{Y}_j^+|𝜽)`$ which equals the term in square brackets in Equation 8. Due to absorption and effective area, a portion of these photons are not detected. The sample counts after absorption (and accounting for effective area) are depicted in the fifth column of Figure 10, by $`\dot{𝐘}^s=(\dot{Y}_1^s,\mathrm{},\dot{Y}_J^s)^{}`$ with

$$\dot{Y}_j^s=\underset{i=1}{\overset{N^s}{}}I\{Y{}_{i}{}^{^{\mathrm{}}}B_j\}(1Z_i^A)\mathrm{for}j𝒥\mathrm{and}s𝒮;$$
(B2)
as described in Table 4, $`Z_i^A`$ is one if photon $`i`$ is absorbed and is zero otherwise. The second plot in Figure 2 represents $`\dot{𝐘}^s`$ with $`\mathrm{E}(\dot{Y}_j^+|𝜽)=\lambda _j(𝜽)d_j`$ plotted as the solid line; see Equation 8. The next two circles in Figure 10 represent the adding of sources and the blurring (i.e., instrument response) process. In particular, $`\dot{𝐘}^+=_{s𝒮B}\dot{𝐘}^s`$. The blurred data, $`𝐘^+=(Y_1^+,\mathrm{},Y_L^+)^{}`$, is a stochastic function of $`\dot{𝐘}^+`$, (i.e., a multinomial distribution171717The multinomial $`(n,𝐩)`$ distribution is a distribution for nonnegative integer valued random vectors and generalizes the binomial distribution. In particular, a vector randomly selected from this distribution sums to $`n`$ and its expected value is $`np`$, where $`p`$ is a probability vector which sums to one.)

$$𝐘^+\stackrel{\mathrm{d}}{}\underset{\mathrm{j}𝒥}{}\mathrm{multinomial}(\dot{\mathrm{Y}}_\mathrm{j}^+,𝐌_\mathrm{j}),$$
(B3)
where $`\dot{𝐘}^+=(\dot{Y}_1^+,\mathrm{},\dot{Y}_J^+)^{}`$ and $`𝐌_j`$ is the $`j`$th column of $`𝐌`$; $`𝐘^+`$ appears in the third plot of Figure 2 with $`\mathrm{E}(Y_l^+|𝜽)=_{j𝒥}M_{lj}\lambda _j(𝜽)d_j`$. The counts due to background contamination are denoted $`\dot{𝐘}^B=(\dot{Y}_1^B,\mathrm{},\dot{Y}_L^B)^{}`$ and the observed data are denoted $`𝐘^{\mathrm{obs}}=(Y_1^{\mathrm{obs}},\mathrm{},Y_L^{\mathrm{obs}})^{}`$ with $`Y_l^{\mathrm{obs}}=Y_l^B+Y_l^+`$ for $`l`$; $`𝐘^{\mathrm{obs}}`$ is illustrated in the final plot of Figure 2 and has expectation $`\xi (𝜽)`$; (cf. Equation 9).

### B.2 The Algorithms

In this section we present the details of the MCMC algorithm which we use to sample from the posterior distribution for our spectral model. We use an algorithm which alternately draws the “missing data” given the model parameters and the parameters given the “missing data”. Both draws are conditional on the observed photon counts and the prior hyperparameters, $`\mathit{\varphi }=\{\mathit{\varphi }^A,\mathit{\varphi }^s,s𝒮\}`$. In particular, we define two groups, (1) the augmented data, $`𝐘^{\mathrm{aug}}=\{𝐘^{\mathrm{obs}},𝐘^B,𝐘^+,\dot{𝐘},\ddot{𝐘},𝐘{}_{}{}^{^{\mathrm{}}}\}`$, where $`\dot{𝐘}=\{\dot{𝐘}^s,s𝒮B\}`$ are the binned true energies, after absorption and accounting for effective area, $`\ddot{𝐘}=\{\ddot{𝐘}^k,k𝒦\}`$ are the binned true energies and $`𝐘{}_{}{}^{^{\mathrm{}}}=\{𝐘{}_{}{}^{^{\mathrm{}}}{}_{}{}^{k},k𝒦\}`$ are the (unbinned) true energies and (2) $`𝜽=\{𝜽^A,𝜽^s,s𝒮\}`$ consists of the various model parameters. Using Bayes’s Theorem, we are able to derive the necessary conditional distributions, which are described below.

First, we draw $`𝐘^{\mathrm{aug}}`$ from $`p(𝐘^{\mathrm{aug}}|𝐘^{\mathrm{obs}},𝜽)`$; the draw is broken into the following five steps.

Draw 1, Step 1: Independently separate the background counts,

$$Y_l^B|𝐘^{\mathrm{obs}},𝜽\stackrel{\mathrm{d}}{}\mathrm{binomial}(\mathrm{Y}_\mathrm{l}^{\mathrm{obs}},\frac{\theta _\mathrm{l}^\mathrm{B}}{\xi _\mathrm{l}(𝜽)}),\mathrm{for}\mathrm{l}.$$
(B4)
Draw 1, Step 2: Restore the blurred the photons,

$$\dot{𝐘}^+|𝐘^B,𝐘^{\mathrm{obs}},𝜽\stackrel{\mathrm{d}}{}\underset{\mathrm{l}}{}\mathrm{multinomial}(\mathrm{Y}_\mathrm{l}^+,\frac{\mathrm{d}_1\lambda _1(𝜽)\mathrm{M}_{1\mathrm{l}},\mathrm{},\mathrm{d}_\mathrm{j}\lambda _\mathrm{j}(𝜽)\mathrm{M}_{\mathrm{Jl}})}{_{\mathrm{j}𝒥}\mathrm{d}_\mathrm{j}\lambda _\mathrm{j}(𝜽)\mathrm{M}_{\mathrm{jl}}}),$$
(B5)
where $`Y_l^+=Y_l^{\mathrm{obs}}Y_l^B`$.

Draw 1, Step 3: Independently separate the counts into line and continuum counts,

$$(\dot{𝐘}^C,\dot{𝐘}^1,\mathrm{},\dot{𝐘}^K)|\dot{𝐘}^+,𝐘^B,𝐘^{\mathrm{obs}},𝜽\stackrel{\mathrm{d}}{}$$
(B6)
$$\mathrm{multinomial}(\dot{Y}_j^+,\frac{(d_jf(𝜽^C,E_j),\stackrel{~}{\lambda }^1p_j^1\mathrm{},\stackrel{~}{\lambda }^Kp_j^K)}{d_jf(𝜽^C,E_j)+_{k=1}^K\stackrel{~}{\lambda }^kp_j^k})\mathrm{for}j𝒥,$$
(B7)
where $`p_j^k=P_j(\mu ^k,\upsilon ^k)`$.

Draw 1, Step 4: Independently restore the absorbed counts in the lines,

$$\ddot{Y}_j^k|\dot{𝐘},𝐘^B,𝐘^{\mathrm{obs}},𝜽\stackrel{\mathrm{d}}{}\dot{\mathrm{Y}}_\mathrm{j}^\mathrm{k}+\mathrm{Poisson}\left(\stackrel{~}{\lambda }^\mathrm{k}\mathrm{p}_\mathrm{j}^\mathrm{k}(1\mathrm{d}_\mathrm{j}\mathrm{u}(𝜽^\mathrm{A},\mathrm{E}_\mathrm{j}))\right),$$
(B8)
for $`j𝒥`$ and $`k𝒦`$.

Draw 1, Step 5: Independently de-round the photon energies from the lines,

$$Y{}_{}{}^{^{\mathrm{}}}{}_{i}{}^{k}|\ddot{𝐘},\dot{𝐘},𝐘^B,𝐘^{\mathrm{obs}},𝜽\stackrel{\mathrm{d}}{}\mathrm{N}(\mu ^\mathrm{k},\upsilon ^\mathrm{k}),\text{truncated to}\mathrm{B}_\mathrm{j}$$
(B9)
for $`i=1,\mathrm{},\ddot{Y}_+^k`$ and $`k𝒦`$ with $`\ddot{𝐘}_+^k=_{j𝒥}\ddot{𝐘}_j^k`$. We note that this draw is omitted for line $`k`$ if $`\upsilon ^k`$ is fixed at zero, i.e., if line $`k`$ is a delta function. (In this case $`\mu ^k`$ is not fit by the algorithm).

Second, we draw $`𝜽`$ from $`p(𝜽|𝐘^{\mathrm{aug}},\mathit{\varphi })`$, taking advantage of conditional independence among several vector components of $`𝜽`$.

Draw 2, Step 1: Independently draw the background model parameters,

$$\theta _l^B|𝐘^{\mathrm{aug}},\mathit{\varphi }\text{}\stackrel{\mathrm{d}}{}\gamma (\varphi _{\mathrm{l},1}^\mathrm{B}+𝐘_\mathrm{l}^\mathrm{B},\varphi _{\mathrm{l},2}^\mathrm{B}+1).$$
(B10)
Draw 2, Step 2: Draw the variance and mean independently for each line profile

$$\upsilon ^k|𝐘^{\mathrm{aug}},\mathit{\varphi }\text{}\stackrel{\mathrm{d}}{}$$
(B11)
$$\frac{1}{\chi _{\nu _0^k+\ddot{Y}_+^k}^2}[\nu _0^k(\sigma _0^2)^k+\underset{i=1}{\overset{\ddot{Y}_+^k}{}}(Y{}_{i}{}^{^{\mathrm{}}}{}_{}{}^{k}Y{}_{}{}^{^{\mathrm{}}}{}_{+}{}^{k}/\ddot{Y}_+^k)^2+\frac{\kappa _0^k\ddot{Y}_+^k}{\kappa _0^k+\ddot{Y}_+^k}(\mu _0^kY{}_{}{}^{^{\mathrm{}}}{}_{+}{}^{k}/\ddot{Y}_+^k)^2],$$
(B12)
$$\mu ^k|\upsilon ^k,𝐘^{\mathrm{aug}},\mathit{\varphi }\text{}\stackrel{\mathrm{d}}{}\mathrm{N}(\frac{\kappa _0^\mathrm{k}\mu _0^\mathrm{k}+\mathrm{Y}{}_{}{}^{^{\mathrm{}}}_{+}^{\mathrm{l}}}{\kappa _0^\mathrm{k}+\ddot{\mathrm{Y}}_+^\mathrm{k}},\frac{\upsilon ^\mathrm{k}}{\kappa _0^\mathrm{k}+\ddot{\mathrm{Y}}_+^\mathrm{k}}),$$
(B13)
where $`Y{}_{}{}^{^{\mathrm{}}}{}_{+}{}^{k}=_{i=1}^{\ddot{Y}_+^k}Y{}_{}{}^{^{\mathrm{}}}_{i}^{k}`$. Again, this step is omitted for line $`k`$ if it is assumed to be a delta function.

Draw 2, Step 3: Draw the line intensities $`\stackrel{~}{\lambda }^k`$, independently for each $`k`$

$$\stackrel{~}{\lambda }^k|𝐘^{\mathrm{aug}},\mathit{\varphi }\text{}\stackrel{\mathrm{d}}{}\gamma (\ddot{\mathrm{Y}}_+^\mathrm{k}+\varphi _1^\mathrm{k},1+\varphi _{\stackrel{~}{\lambda }}^\mathrm{k})\mathrm{independently}\mathrm{for}\mathrm{k}𝒦.$$
(B14)
Draw 2, Step 4: Draw the parameters for the GLM for the continuum and absorption models. For this final step, we condition only on $`𝐘^{\mathrm{obs}}`$, $`𝐘^B`$, $`𝐘^+`$, and $`\dot{𝐘}^+`$, rather than $`𝐘^{\mathrm{aug}}`$. We expect this substitution to improve the rate of convergence of the sampler. Because the conditional distribution is not from a standard family, we use a Metropolis Hastings step. In particular, we note that

$$\dot{Y}_j^+|𝜽,\stackrel{~}{𝝀}\stackrel{\mathrm{d}}{}\mathrm{Poisson}\left(\mathrm{d}_\mathrm{j}\mathrm{exp}(𝐗_\mathrm{j}^\mathrm{A}𝜽^\mathrm{A})\underset{\mathrm{k}=1}{\overset{\mathrm{K}}{}}\stackrel{~}{\lambda }^\mathrm{k}\mathrm{p}_\mathrm{j}^\mathrm{k}\right)\mathrm{for}\mathrm{k}𝒦$$
(B15)
$$\dot{Y}_j^C|𝜽,\stackrel{~}{𝝀}\stackrel{\mathrm{d}}{}\mathrm{Poisson}\left(\mathrm{dE}_\mathrm{j}\mathrm{d}_\mathrm{j}\mathrm{exp}(𝐗_\mathrm{j}^\mathrm{C}𝜽^\mathrm{C}+𝐗_\mathrm{j}^\mathrm{A}𝜽^\mathrm{A})\right)\mathrm{for}\mathrm{j}𝒥,$$
(B16)
where $`\stackrel{~}{𝝀}=(\stackrel{~}{\lambda }^1,\mathrm{},\stackrel{~}{\lambda }^K)`$. Since given $`(𝜽^1,\mathrm{},𝜽^K)`$ the log of each Poisson parameter differs from a linear combination of $`𝜽^C`$ and $`𝜽^A`$ by a known constant, the conditional posterior mode can easily be computed using a minor modification of an Iteratively Reweighted Least Squares algorithm (e.g., Thisted, 1988). This algorithm can also account for the prior information described in Section 3.2 and reports the curvature of the log posterior at the mode. A multivariate t-distribution with four degrees of freedom with the appropriate mode and (perhaps inflated) curvature can be used as a jumping distribution to generate a proposal for the next sample from the conditional distribution. The relative mass of the jumping distribution and actual conditional distribution of $`𝜽^C`$ and $`𝜽^A`$ at the previous draw and proposed draw are combined to determine if the proposal should be accepted or rejected (in which case the previous draw is reused). Several (5 – 10) proposals are drawn at each iteration. We note that the same procedure can be used to fit a GLM to the background counts (as was done in Section 4.2).

Although the MCMC methods detailed above may seem inhibiting as a whole, each of the required steps are quite simple. The power of the MCMC methods described here (e.g., the Gibbs sampler) lies in their ability to break complicated model fitting tasks into a succession of relatively simple tasks. Our general strategy is to use Bayes’s theorem to derive a posterior distribution which hierarchically accounts for the complexity both in the posited model and in data collection. We then use modern statistical algorithms which devolve model fitting into a sequence of relatively simple steps. We believe this is a powerful strategy for dealing with the ever increasing power and sophistication of today’s astronomical instruments.

## Appendix C Internet Resources

There are several internet sites where one can find papers describing Bayesian methods and related software. The MCMC preprint service,

http://www.mcs.surrey.ac.uk/Personal/S.Brooks/MCMC,

and STATLIB,

http://lib.stat.cmu.edu

are both large general statistical sites that offer various software, preprints, and links that may be of interest to astrophysicists. Three sites (that we know of) aim specifically at the interface of astrophysics and statistics:

http://www.fas.harvard.edu/$``$vandyk/astrostat.html

http://astrosun.tn.cornell.edu/staff/loredo/bayes

http://www.astro.psu.edu/statcodes

Caption for Figure 10: A graphical representation of the data-augmentation scheme. Here $`\mathit{\varphi }`$ represents hyperparameters, $`𝜽`$ model parameters, $`𝐘^{^{\mathrm{}}}`$ true photon energies, $`\ddot{𝐘}`$ binned energies, $`\dot{𝐘}`$ binned true photon energies after absorption accounting for effective area, $`𝐘`$ source counts in PHA channels, $`𝐘^{\mathrm{obs}}`$ the observed counts, $`𝐌`$ the instrument response matrix, $`𝐝`$ the effective area vector, $`𝐗^A`$ and $`𝐗^C`$ independent variables describing absorption and continuum respectively, circles represent unobserved quantities, and squares observed quantities; details of the subscripts and superscripts are given in the text. The figure illustrates the interplay of the various model parameters, hyperparameters, observed quantities, and data augmentation. As an example, the first arrow in the row labeled ‘background’ corresponds to the relationship between the background hyperparameters, $`\mathit{\varphi }^B`$, and the background intensities, $`𝜽^B`$, e.g., $`\theta _l^B\stackrel{\mathrm{d}}{}\gamma (\varphi _{\mathrm{l},1}^\mathrm{B},\varphi _{\mathrm{l},2}^\mathrm{B})`$. The second arrow corresponds to the Poisson nature of the background counts; see Equation 22. The final background arrow illustrates the background contamination of the observed PHA counts.