DISENTANGLING OVERLAPPING ASTRONOMICAL SOURCES USING SPATIAL AND SPECTRAL INFORMATION

Here we use finite mixture distributions to model several overlapping sources of photons in a high-energy image. Finite mixture distributions are a useful class of statistical models for data that are drawn from a mixture of several subpopulations; these models are finite in that the (possibly unknown) number of subpopulations is a finite positive integer. (See Titterington et al. 1985 and McLachlan & Peel 2004 for comprehensive discussion of finite mixture distributions, and, for example, Mukherjee et al. 1998 for a previous application in astronomy.) We take a Bayesian perspective that allows joint inference for the parameters that describe the photon sources (e.g., their number, intensities, and locations), the basic shape of their spectra, and the probability that any particular photon originated from each source, given its recorded location and energy.

Performing inference jointly on the image and spectra improves the precision of the fitted parameters, and also provides more coherent measures of uncertainty than would be available if the spatial and spectral data were analyzed separately. Furthermore, unlike other methods for overlapping sources, our approach quantifies uncertainty about the number of sources. Whether we are ultimately interested in spatial or spectral aspects of sources, identifying the correct number of sources is clearly fundamental. Consequently, a coherent measure of the uncertainty associated with the fitted number of sources is critical to the appropriate interpretation of the fitted parameters of the individual sources.

In some applications inference for the number of sources may seem unnecessary because the sources are clearly identifiable. For instance, the XMM-Newton observation of FK Aqr and FL Aqr analyzed in Section 6 has relatively weak background noise, and the sources overlap only moderately. In such cases, the main advantage of the proposed method is that it precisely quantifies the uncertainties associated with the positions, intensities and spectral shapes of the sources. As already mentioned, finite mixture analysis also yields, for each observed photon, the probability that it originated from each inferred source (or the background). In this way we do not deterministically assign photons to sources, but rather properly assess the uncertainty of their origins. This is in contrast to other methods, such as those based on source regions, which deterministically assign photons to nearby sources, and therefore do not properly quantify uncertainties in fitted source parameters.

There is a potential for overfitting in finite mixture models if the number of sources is unknown. This is mitigated when substantial prior information regarding the shape of the PSF or the number of sources is available, or both. In practice, we have detailed information about the PSF, and hence know exactly what the distribution of the recorded photon locations should be for each source. (For point sources this is trivial, but even for extended sources one can easily convolve the source model with the PSF.) Even if the PSF varies across the field, the shape of the photon scatter is completely determined by the location of the source. With this complete knowledge of the PSF, there is only a small risk of overfitting, even with limited prior information regarding the number of sources and their spectral shapes. Indeed, our results do not strongly depend on the choice of prior distribution for the number of sources (see Section 5.1).

Our method is designed for analyzing images composed of an unknown number of point sources that are contaminated with background. However, it can be applied to extended sources, with some modifications to account for spatial variations in intensity and spectra. We also mention that the success of our method depends partly on our ability to use spectra to distinguish point sources from the background, which is possible because a typical X-ray point source spectrum is more peaked than the background. Because of this, we are able to use basic models that capture the rough spectral shape in order to exploit spectral information while conserving computational resources. In the X-ray band, this approach offers substantial improvements over analyses using only spatial data without the cost of precisely modeling the spectra. However, the utility of the method in other wavelength bands will depend somewhat on the nature of the spectra typical of those bands.

The remainder of the paper is organized into seven sections. Section 2 develops the statistical model for isolated sources in the context of high-energy data sets, and describes how these models are combined in the case of multiple sources. Section 3 uses a motivating example to illustrate the method and the benefits of incorporating spectral models for the sources. The beginning of Section 4 gives a brief review of Bayesian inference. The remainder of Section 4 describes the details of the proposed Bayesian analysis and computational approach. Section 5 presents two simulation studies. The first illustrates that inference for the number of sources is insensitive to the choice of prior distribution, and the second more thoroughly studies the advantages of using the spectral data. Sections 6 and 7 present the results of our analysis of observations from the XMM-Newton and Chandra X-Ray Observatory. The XMM observation is of the apparent visual binary FK Aqr and FL Aqr, and the Chandra observation is of approximately 14 sources from near the center of the Orion nebula. We summarize in Section 8 and computational details are in the appendices. Our Bayesian Separation of Close Sources (BASCS) software is available on GitHub at https://github.com/astrostat/BASCS.

High-energy detectors record directional coordinates $({x}_{i},{y}_{i})$ and energy ${E}_{i}$ for each detected photon, where $i=1,...,n$ indexes the photons. As mentioned, in practice, the PI channel is used to quantify energy. We denote the full set of spatial and spectral information for $n$ detected photons by $(x,y,E)$. These observed quantities are subject to the effects of the PSF and the spectral Redistribution Matrix Function (RMF). We explicitly account for PSF effects in our model, but model the observed spectra, the convolution of the source spectra and the RMF. This strategy does not allow us to fit source spectral models, but does allow us to leverage spectral data to separate the sources. Even though all the attributes are recorded digitally and are binned quantities, we treat them as continuous variables for simplicity, since this binning is at scales that heavily over-sample the PSF.

Each photon is assumed to originate from one of several point sources or the background, but its exact origin is unknown. Furthermore, the number of point sources contributing photons to the data, their locations, intensities, and spectral distributions are all unknown. We assume background is distributed uniformly across the image, its strength and spectral distribution are often not known.

To introduce notation and our model in the simplest case, we first suppose that the data consist of photons from a single source, with no background contamination. Statistical models specify a distribution for the observed data conditional on a number of typically unknown parameters; we discuss parameter fitting Section 4.3. In the current case, given the unknown position of the source, the detected photons are assumed to be dispersed according to a PSF. That is,

$$\begin{eqnarray}&&({x}_{i},{y}_{i})| \mu \sim {\rm{PSF}}\;{\rm{centered}}\;{\rm{at}}\;{\rm{}}\mu \end{eqnarray} \tag{ 1 }$$

for $i=1,...,n$, where $\mu =({\mu }_{x},{\mu }_{y})$ is the unknown position of the source.⁵ We use the same 2D King profile⁶ in all the simulations and data analyses presented, see Read et al. (2010)⁷ and King (1962). The King profile density, shown in Figure 15 in Appendix C, has heavy tails and is essentially a bivariate Cauchy distribution. Specific parameter values are detailed in Appendix C. More generally, although our method assumes that the PSF is known given $\mu$, it may vary with $\mu$. Furthermore, the PSF may be any function which can be quickly evaluated analytically or numerically. Even in cases where computationally expensive evaluations are required our method is feasible if the PSF is first tabulated.

An important feature of our overall approach is that it also utilizes the spectral data to better assess the likely origin of each photon (when background or more than one source is present). With this end in mind, we propose a simple and computationally practical model for the basic shape of the source spectrum. In particular, we model photon energies using a gamma distribution,⁸

$$\begin{eqnarray}&&{E}_{i}| \alpha ,\gamma \sim \;{gamma}\;(\alpha ,\alpha /\gamma )\end{eqnarray} \tag{ 2 }$$

for $i=1,...,n$. Here, $\alpha$ and $\gamma$ are the unknown shape and mean parameters used to describe the basic spectral distribution.⁹ The gamma distribution allows flexible modeling of positive quantities with right skewed distributions.¹⁰ We emphasize that we aim to summarize the essential shape of the spectral distribution, rather than to model the details of emission lines and other spectral features. This is practical because for high-energy missions, the effective areas are typically small at low and high energies, with a broad peak in the middle; the resulting counts spectrum is reasonably modeled by a single- or double-component gamma distribution (particularly since we ignore the RMF). Our goal is to identify sources and divide photons among them, not to carry out detailed spectral analysis. However, our algorithm allows for complex spectral models to be built in if necessary. In addition, and computationally more feasible, once the gamma model has fulfilled its role in separating sources, a more sophisticated spectral model may then be used to draw scientific conclusions about the spectral distributions of the disentangled sources. This final stage will be discussed in Section 7.2.

In practice there are multiple sources and background contamination, hence we introduce a finite mixture model. Let $K$ be a parameter denoting the number of sources and $\mu =({\mu }_{1},...,{\mu }_{K})$ be a $2\times K$ matrix giving the source positions i.e., ${\mu }_{j}=({\mu }_{{jx}},{\mu }_{{jy}})$, for $j=1,...,K$. If we knew the origin of every photon then, we could model the spatial and spectral data associated with each point source as we did in Section 2.2. We thus introduce a new variable ${s}_{i}$ which indicates the source number associated with photon i. Each ${s}_{i}$ takes on a value between 1 and $K$, and we let $s$ denote the vector $({s}_{1},...,{s}_{n})$. Note that ${s}_{i}$ is never actually observed and thus is a latent variable. A latent variable is essentially an unknown parameter which is useful for modeling, but may not be of direct interest in itself. Here, we have introduced ${s}_{i}$ to simplify the model and to facilitate the algorithms used for inference, which are described in Section 4.3.

As a parameter, ${s}_{i}$ is also conditioned on in our spatial model, which now becomes

$$\begin{eqnarray}&&({x}_{i},{y}_{i})| (\mu ,{s}_{i}=j)\sim {\rm{PSF}}\;{\rm{centered}}\;{\rm{at}}\;{\rm{}}{\mu }_{j}\end{eqnarray} \tag{ 3 }$$

for $i=1,...,n$. As an unknown parameter, ${s}_{i}$, plays a role similar to μ; it is "given" in Equation (3). The spectral model can also be straightforwardly generalized to the multiple source case. We have

$$\begin{eqnarray}&&{E}_{i}| ({\alpha }_{j},{\gamma }_{j},{s}_{i}=j)\sim {gamma}({\alpha }_{j},{\alpha }_{j}/{\gamma }_{j})\end{eqnarray} \tag{ 4 }$$

for $i=1,...,n$, where the parameters ${\alpha }_{j}$ and ${\gamma }_{j}$ usually differ among the sources.

In addition to point sources, we must model the background. To this end we extend the set of possible values of ${s}_{i}$ to include 0. Throughout, symbols indexed by 0 refer to the background. We assume that photons originating from the background are uniformly distributed across the image,

$$\begin{eqnarray}&&({x}_{i},{y}_{i})| (\mu ,{s}_{i}=0)\sim {\rm{Uniform}}\end{eqnarray} \tag{ 5 }$$

for i such that ${s}_{i}=0$. Instrument effects may cause the background to be non-uniform, and a refinement would be to model such effects.

The background spectrum is also assumed to be flat over the energy range of the source spectra. That is, it is assumed to have a uniform distribution on $({E}_{\mathrm{min}},{E}_{\mathrm{max}})$, where ${E}_{\mathrm{min}}$ and ${E}_{\mathrm{max}}$ are the minimum and maximum photon energy observed. This is a good approximation because the background spectrum is expected to be less peaked than that of a point source.

So far we have not considered the intensities of the different sources and the background. Naturally there should be a parameter for each source, and one for the background, to specify the intensities. Let ${n}_{j}$ denote the number of photons originating from source j, for $j=0,...,K$ (with zero denoting the background), mathematically,¹¹ ${n}_{j}={\displaystyle \sum }_{i=1}^{n}{I}_{\{{s}_{i}=j\}}$. We can realistically model ${n}_{j}$ as a Poisson variable with some mean m_j, for $j=0,...,K$. Because these Poisson means vary with exposure time, however, the relative intensities, ${w}_{j}={m}_{j}/\sum {m}_{j}$, are of more direct interest. Writing $w=({w}_{0},...,{w}_{K})$, and given $n$, the Poisson model for $({n}_{0},...,{n}_{K})$ yields a Multinomial model,¹²

$$\begin{eqnarray}&&({n}_{0},...,{n}_{K})| w,n,K\sim {\rm{Multinomal}}(n;w),\end{eqnarray} \tag{ 6 }$$

where ${\displaystyle \sum }_{j=0}^{K}{w}_{j}=1$. Under this parameterization, the relative strengths of the sources and background can be succinctly expressed by the vector $w=({w}_{0},...,{w}_{K})$ without further reference to $n$. Accordingly, all inference is performed given n, because its value tells us nothing about the number of sources or their parameters.

To complete our introduction of the model we derive the likelihood function, which is the probability of the data expressed as a function of the parameters. The likelihood tells us what values of the parameters are supported by the data and is a key component for principled statistical inference. Let ${{\mathcal{J}}}_{j}$ be the set of photons originating from source j (including j = 0) and let ${\mathcal{J}}$ be the entire collection of observed photons.¹³ Also, denote the value of the PSF centered at $\mu$ and evaluated at $(x,y)$ by ${f}_{\mu }(x,y)$. Lastly, here and throughout, we let ${\theta }_{j}=\{{w}_{j},{\mu }_{j},{\alpha }_{j},{\gamma }_{j}\}$ denote the parameters associated with source j, for $j=1,...,K$. Similarly, for the background, we let ${\theta }_{0}={w}_{0}$. We let ${{\rm{\Theta }}}_{K}$ denote all the source (and background) specific parameters i.e., ${{\rm{\Theta }}}_{K}=\{{\theta }_{0},...,{\theta }_{K}\}$. The remaining parameters are $K$ and $s$. As already discussed, we treat $n$ as fixed, and impose the constraint that the likelihood is zero unless ${\displaystyle \sum }_{j=0}^{K}{n}_{j}=n$. Combining the different parts of the model yields the full model likelihood

$$\begin{eqnarray} &&{L}_{n}^{{\rm{full}}}({{\rm{\Theta }}}_{K},K)\equiv p(x,y,E| {{\rm{\Theta }}}_{K},K,s,n)\\ &&\propto \displaystyle \prod _{i\in {\mathcal{J}}/{{\mathcal{J}}}_{0}}{f}_{{\mu }_{{s}_{i}}}({x}_{i},{y}_{i}){g}_{{\alpha }_{{s}_{i}},{\gamma }_{{s}_{i}}}({E}_{i}), \end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{g}_{{\alpha }_{{s}_{i}},{\gamma }_{{s}_{i}}}({E}_{i})=\displaystyle \frac{{\alpha }_{{s}_{i}}^{{\alpha }_{{s}_{i}}}}{{\gamma }_{{s}_{i}}^{{\alpha }_{{s}_{i}}}{\rm{\Gamma }}({\alpha }_{{s}_{i}})}{E}_{i}^{{\alpha }_{{s}_{i}}-1}{e}^{-{\alpha }_{{s}_{i}}{E}_{i}/{\gamma }_{{s}_{i}}}.\end{eqnarray} \tag{ 8 }$$

The maximum energy ${E}_{\mathrm{max}}$ and the image area are assumed to be known quantities, rather than parameters to be inferred. They are therefore omitted from the likelihood, as are all terms not involving the parameters. In later sections, we compare analyses under the full model to analyses under the spatial-only model that does not use the spectral information. The likelihood of the spatial-only model is

$$\begin{eqnarray}&&{L}_{n}^{{\rm{sp}}}({{\rm{\Theta }}}_{K}^{{\rm{sp}}},K)\equiv p(x,y| {{\rm{\Theta }}}_{K}^{{\rm{sp}}},K,s,n)\propto \displaystyle \prod _{i\in {\mathcal{J}}/{{\mathcal{J}}}_{0}}{f}_{{\mu }_{{s}_{i}}}({x}_{i},{y}_{i}).\end{eqnarray} \tag{ 9 }$$

The notation ${{\rm{\Theta }}}_{K}^{{\rm{sp}}}=\{{w}_{0},...,{w}_{K};{\mu }_{1},...,{\mu }_{K}\}$ represents the set of spatial parameters. Note that, although $w$ does not explicitly appear in either likelihood, the data does nevertheless constrain $w$ in both cases. In particular, the likelihoods indicate probable values of $s$ which in turn indicate probable values of $w$. Conceptually, our method is to apply Bayes rule, briefly reviewed in Section 4.1, to the likelihoods displayed in Equations (7) and (9) to yield a distribution summarizing our knowledge of the parameters given the data, i.e., the joint posterior distribution.

In some situations the gamma spectral model given by Equation (4) is not sufficiently flexible to capture the spectral shape of the observed sources. For example, Figure 1 shows the observed spectrum of the brightest source in the Chandra observation analyzed in Section 7. In particular, the histogram shows the spectrum using one likely assignment of photons produced during the iterations of our algorithm (see Section 4.3). The dashed red curve shows the maximum likelihood fit of the gamma distribution to the observed spectrum. The gamma does not fit the distribution closely. This causes a problem because inference based on the (misspecified) gamma spectral model will suggest there are two sources instead of one in order to better capture the spectral distribution of the source.

Zoom In Zoom Out Reset image size

To solve this problem, we use a mixture of two gamma distributions for a more general spectral model. That is,

$$\begin{eqnarray} &&{E}_{i}| ({\alpha }_{j1},{\alpha }_{j2},{\gamma }_{j1},{\gamma }_{j2},{\pi }_{j1},{\pi }_{j2},{s}_{i}=j)\\ &&\quad \sim {\pi }_{j1}{gamma}\left({\alpha }_{j1},\displaystyle \frac{{\alpha }_{j1}}{{\gamma }_{j1}}\right)+{\pi }_{j2}{gamma}\left({\alpha }_{j2},\displaystyle \frac{{\alpha }_{j2}}{{\gamma }_{j2}}\right), \end{eqnarray} \tag{ 10 }$$

for $i=1,...,n$, where the parameters ${\pi }_{j}$ and ${\pi }_{j2}=1-{\pi }_{j1}$ are the weights of the two gamma components. When this two gamma mixture spectral model is substituted for the one gamma spectral model in Equation (7) we obtain the following extended full model likelihood

$$\begin{eqnarray} &&{L}_{n}^{{\rm{ext}}}({{\rm{\Theta }}}_{K}^{{\rm{ext}}},K)\equiv p(x,y,E| {{\rm{\Theta }}}_{K}^{{\rm{ext}}},K,s,n)\\ &&\propto \displaystyle \prod _{i\in {\mathcal{J}}/{{\mathcal{J}}}_{0}}\left({f}_{{\mu }_{{s}_{i}}}({x}_{i},{y}_{i})\displaystyle \sum _{l=1}^{2}{\pi }_{{s}_{i}l}{g}_{{\alpha }_{{s}_{{il}}},{\gamma }_{{s}_{{il}}}}({E}_{i})\right). \end{eqnarray} \tag{ 11 }$$

The notation ${{\rm{\Theta }}}_{K}^{{\rm{ext}}}$ denotes $\{{\theta }_{0}^{{\rm{ext}}},...,{\theta }_{K}^{{\rm{ext}}}\}$, where ${\theta }_{j}^{{\rm{ext}}}=\{{w}_{j},{\mu }_{j},{\alpha }_{j1},{\alpha }_{j2},{\gamma }_{j1},{\gamma }_{j2},{\pi }_{j1},{\pi }_{j2}\}$ gives the parameters associated with source j, for $j=1,...,K$, and ${\theta }_{0}^{{\rm{ext}}}={\theta }_{0}$. The solid green curve in Figure 1 shows the extended full model fit of the gamma mixture spectral model. In this example, the mixture of gammas quite closely fits the observed spectrum and generally there did not appear to be unwarranted splitting of sources into two in our numerical studies using this model.

Even greater flexibility of the spectral model could be gained by considering a mixture of more than two gammas, but this was not necessary in our numerical studies. For the XMM data of Section 6, the one-gamma spectral model is sufficient in that, for both of the sources, the maximum likelihood fit of the one-gamma and the two-gamma models resulted in essentially identical fits when using a feasible allocation of photons. In the interest of simplicity, we only use the extended full model when necessary (i.e., in Section 7), and elsewhere use the full model given in Equation (7).

2.4.1. Detecting Spectral Model Inadequacy

The Bayesian perspective provides a coherent approach for combining all available information to infer the unknown model parameters ${{\rm{\Theta }}}_{K}$, $K$, and $s$. First, our knowledge (or lack of knowledge) as to the likely values of the parameters before seeing the current data is quantified using a prior distribution. Once the data are observed, Bayes' Theorem allows us to combine the likelihood and the prior distribution to yield the posterior distribution of the parameters. Recall, the likelihood is the probability of the data given the parameters. The posterior distribution expresses our updated knowledge of the parameters after seeing the data. Bayes' Theorem states that, for generic data and parameter vector $\psi$, the posterior distribution is

$$\begin{eqnarray}&&p(\psi | {\rm{data}})=\displaystyle \frac{p({\rm{data}}| \psi )p(\psi )}{p({\rm{data}})},\end{eqnarray} \tag{ 12 }$$

where $p({\rm{data}}| \psi )$ is the likelihood function and $p(\psi )$ is the prior distribution. The denominator $p({\rm{data}})$ is simply a normalizing constant which ensures the posterior integrates to one. In our case, the data is $(x,y,E)$ and under the full model $\psi =\{{{\rm{\Theta }}}_{K},K,s\}$ so

$$\begin{eqnarray}&&p({{\rm{\Theta }}}_{K},K,s| x,y,E)=\displaystyle \frac{p(x,y,E| {{\rm{\Theta }}}_{K},K,s)p({{\rm{\Theta }}}_{K},K,s)}{p(x,y,E)}.\end{eqnarray} \tag{ 13 }$$

Here, all probabilities are conditional on n but this is suppressed. The likelihood $p(x,y,E| {{\rm{\Theta }}}_{K},K,s)$ is given in Equation (7), and the prior distribution $p({{\rm{\Theta }}}_{K},K,s)$ is described in Section 4.2. Referring back to the illustrative example in Section 3, the marginal posterior distribution of $K$,

$$\begin{eqnarray}&&p(K| x,y,E)=\displaystyle \sum _{s}\int p(K,s,{{\rm{\Theta }}}_{K}| x,y,E)d{{\rm{\Theta }}}_{K},\end{eqnarray} \tag{ 14 }$$

is displayed in Figure 3. Given the number of unknown parameters, it is not possible to plot their joint posterior distribution, but we can derive and plot the marginal posterior distribution of any one parameter, as in Equation (14) and Figure 3.

Following the Bayesian approach, we specify prior distributions for each of the unknown parameters. First, the positions of the point sources are a priori assumed to be independently and uniformly distributed across the image. That is,

$$\begin{eqnarray}&&{\mu }_{j}\sim {\rm{Uniform}}\end{eqnarray} \tag{ 15 }$$

for $j=1,...,K$. In principle, informative priors can be used if prior information on source locations is available. For example, we might set ${\mu }_{j}\sim N({\mu }_{j0},{\sigma }_{j0}^{2})$, where $({\mu }_{j0},{\sigma }_{j0})$, for $j=1,...,K$, specifies knowledge of the source locations.¹⁶

Next, the vector $w$, that specifies relative intensities, is given a Dirichlet¹⁷ prior distribution with parameter $(\lambda ,...,\lambda )$. A Dirichlet random variable is a probability vector, i.e., it is a vector with non-negative entries that sum to one. We set $\lambda =1$ throughout. This choice is uniform on the probability vector, but very slightly favors sources of equal size. Indeed, setting $\lambda =1$ means the Dirichlet prior has as much information as a single photon count added to each source (including a single count added to the background).¹⁸ Regarding the realized vector of source and background counts $({n}_{0},...,{n}_{K})$, recall that Equation (6) specifies a Multinomial distribution for $({n}_{0},...,{n}_{K})$, given $w$, $n$, and $K$. Since $({n}_{0},...,{n}_{K})$ is a function of the parameter (or latent variable) $s$, Equation (6) is effectively a prior distribution for $s$.¹⁹

External information about the number of sources is amalgamated into a prior for $K$, which we assume to be Poisson with mean parameter $\kappa$.²⁰ Under the Poisson prior, the fitted value of $K$ is relatively robust to the choice of $\kappa$ because the PSF is completely specified.²¹ Indeed, we show in Section 5.1 that the posterior mode for the number of sources may correctly identify the true value of $K$, even when $\kappa$ is quite different from $K$. Therefore, in practice it is adequate to use the Poisson prior for $K$ with $\kappa$ set to any reasonable guess of the number of sources.

To complete the model specification, we must assign prior distributions for the source spectral distribution parameters ${\alpha }_{j}$ and ${\gamma }_{j}$, for $j=1,...,K$. Typically there is sufficient data to overwhelm these prior distributions. Thus, we are not overly concerned with the exact form of these priors. For concreteness, however, we mention that one set of priors we use is ${\alpha }_{j}\sim {gamma}(2,0.5)$ and ${\gamma }_{j}\sim {\rm{Uniform}}({E}_{\mathrm{min}},{E}_{\mathrm{max}})$, for $j=1,...,K$, where ${E}_{\mathrm{min}}$ is the minimum observed energy.²²

To summarize, our prior distribution for the full model parameters ${{\rm{\Theta }}}_{K}$, $K$ and $s$ is

$$\begin{eqnarray} p({{\rm{\Theta }}}_{K},K,s) & = & p(\mu ,\alpha ,\gamma ,s| K,w)p(w| K)p(K)\\ & \propto & \left(\displaystyle \prod _{j=0}^{K}{\alpha }_{j}{e}^{-0.5{\alpha }_{j}}\right)\displaystyle \prod _{j=0}^{K}{w}_{j}^{{n}_{j}}{\left(\displaystyle \prod _{j=0}^{K}{w}_{j}\right)}^{\lambda -1}\displaystyle \frac{{\kappa }^{K}}{K!}, \end{eqnarray} \tag{ 16 }$$

where $\mu$, $\alpha$ and $\gamma$ denote $({\mu }_{1},...,{\mu }_{K})$, $({\alpha }_{1},...,{\alpha }_{K})$, and $({\gamma }_{1},...,{\gamma }_{K})$, respectively. The second term on the second line of Equation (16) comes from the Multinomial prior distribution for $s$. In the case of the extended full model given in Equation (11), the priors for ${\alpha }_{{jl}},{\gamma }_{{jl}}$, $l=1,2$, are the same as those for ${\alpha }_{j},{\gamma }_{j}$, and the prior for ${\pi }_{j1}$ is a $\mathrm{Beta}(2,2)$ distribution,²³ for $j=1,...,K$. (No prior for ${\pi }_{j2}$ is needed because this parameter is determined by ${\pi }_{j1}$, for $j=1,...,K$.) The prior for the spatial model parameters is Equation (16) without the first term.

Given the likelihood in Equation (7) and the prior distribution in Equation (16), we can apply Bayes' Theorem to obtain the posterior distribution of ${{\rm{\Theta }}}_{K}$, $K$ and $s$ (see Equation (13)). The resulting posterior distribution is a complicated function, which we summarize by the low-dimensional marginal distributions as described in Section 4.1 and their means and standard deviations. These summaries are used to estimate the model parameters and their error bars.

We accomplish the necessary numerical integration, e.g., as in Equation (14), using Monte Carlo methods, a cornerstone of statistical computing (Shao & Ibrahim 2000; Liu 2008; Christian & Casella 2011, p. 49–66). The idea of Monte Carlo algorithms is to simulate values of the generic parameter $\psi$ from the posterior distribution in Equation (12) to obtain a Monte Carlo sample $\{{\psi }^{(1)},...,{\psi }^{(T)}\}$. For example, in Figure 3, the height of the bin centered at $k$ is the proportion of the Monte Carlo draws with ${K}^{(t)}$ equal to $k$, i.e.,

$$\begin{eqnarray}&&P(K=k| x,y,E)\approx \displaystyle \frac{1}{T}\displaystyle \sum _{t=1}^{T}{I}_{\{{K}^{(t)}=k\}},\end{eqnarray} \tag{ 17 }$$

for $k=1,...,K$.

A somewhat unusual feature of our model is that the number of parameters is determined by the value of $K$, the unknown number of sources. This necessarily conditional structure means that it only makes sense to consider the posterior distributions of the other parameters for a given inferred value of $K$ (Park et al. 2008 discuss a somewhat similar conditional inference in the context of locating emission lines). For an illustration of why this is so, consider the intensity ${w}_{3}$ of the "third" source in an image. The parameter ${w}_{3}$ does not have the same interpretation when there are three sources versus four, because what is the "third" source in the first scenario may combine two sources from the latter scenario. In fact, for $K=2$ the parameter ${w}_{3}$ does not even exist. In general, there is no clear relationship between the parameters under scenarios with different values of $K$. This prevents us from considering the unconditional posterior distribution of, say, ${w}_{3}$. Instead, we are interested in posterior summaries given a particular value of $K$, such as $p({w}_{3}| K=k,x,y,E)$. For example, the second row of Table 2 provides an estimate of the posterior mean of ${w}_{2}$ conditional on $K=3$, under the full model,²⁴

$$\begin{eqnarray}&&{\hat{w}}_{2}^{F}(k)=\displaystyle \frac{\sum _{t=1}^{T}{w}_{2}^{(t)}{I}_{\{{K}^{(t)}=k\}}}{\sum _{t=1}^{T}{I}_{\{{K}^{(t)}=k\}}}=0.080.\end{eqnarray} \tag{ 18 }$$

More generally, for each one-dimensional parameter τ, we calculate the Monte Carlo estimate

$$\begin{eqnarray}&&\hat{\tau }(k)=\displaystyle \frac{\sum _{t=1}^{T}{\tau }^{(t)}{I}_{\{{K}^{(t)}=k\}}}{\sum _{t=1}^{T}{I}_{\{{K}^{(t)}=k\}}}.\end{eqnarray} \tag{ 19 }$$

In practice, we choose a value of $k$ at which $K$ has relatively high posterior probability, such as the posterior mode, because otherwise the parameters estimated are unlikely to correspond to properties of real sources. (Indeed, our algorithm does not accurately estimate parameters under unlikely values of $K$.) We may decide to consider several different values if the posterior of $K$ is not concentrated on one value. This can be useful despite the fact that, as we have mentioned, the number and interpretation of the parameters is not consistent across values of $K$.

The most popular method for obtaining the Monte Carlo samples needed for estimates such as that in Equation (18) is Markov chain Monte Carlo (MCMC). This an iterative algorithm in which we generate a new value of the parameters ${\psi }^{(t)}$ at each iteration by drawing from a distribution ${\mathcal{G}}$ that only depends on ${\psi }^{(t-1)}$ (and the data) and not earlier members of the Monte Carlo sample. Continuing for T iterations we obtain a sample $\{{\psi }^{(1)},...,{\psi }^{(T)}\}$ of correlated parameter values, which is usually called an MCMC chain. Appropriate choice of ${\mathcal{G}}$ ensures that the sample mimics the posterior distribution in the sense that as $T\to \infty$ the sample empirical distribution approaches the posterior distribution. In implementation, a draw from an appropriate ${\mathcal{G}}$ is typically achieved through two steps: first a new value of the parameters ${\psi }^{*}$ is proposed, and then this value is either accepted or rejected with some probability.²⁵ The Metropolis–Hastings algorithm (Metropolis et al. 1953 and Hastings 1970) is an example of such an algorithm. The reader is referred to Gelman et al. (2013) for details, including discussion of efficient choices of ${\mathcal{G}}$ and monitoring of convergence to the posterior distribution (which is usually done by running multiple MCMC chains in parallel and checking that their behavior is sufficiently similar based on some criterion).

In standard MCMC algorithms the parameter space being explored is fixed throughout. In our context this means the number of sources would have to be known. We therefore turn to reversible jump Markov chain Monte Carlo (RJMCMC) algorithms (first introduced by Green 1995), which allow configurations with differing numbers of sources to be explored. There have been a number of uses of RJMCMC in other astronomy contexts, for example, Umstätter et al. (2005), Brewer & Stello (2009), Jasche & Wandelt (2013), and Walmswell et al. (2013). In RJMCMC algorithms, so called "jump" steps update the value of $K$, the name referring to a jump between parameter spaces (or "models"). These steps are performed by drawing ${K}^{(t)}$ from a distribution only depending on ${\psi }^{(t-1)}=({{\rm{\Theta }}}^{(t-1)},{K}^{(t-1)},{s}^{(t-1)})$, in the same spirit as ordinary MCMC iterations. Feasible values of the parameters ${{\rm{\Theta }}}^{(t)}$ and ${s}^{(t)}$ must simultaneously be drawn because their dimension and interpretation change with $K$. It is this high dimensional sampling that makes RJMCMC challenging. In RJMCMC algorithms, ${K}^{(t)}$ is only allowed to differ from ${K}^{(t-1)}$ by at most one. This constraint facilitates the proposal of appropriate parameters ${{\rm{\Theta }}}^{(t)}$ and ${s}^{(t)}$; RJMCMC moves between configurations by splitting, combining, creating or destroying sources in the model. The standard RJMCMC algorithm for Gaussian mixtures was introduced in Richardson & Green (1997), and Wiper et al. (2001) illustrated RJMCMC for gamma mixtures. Our BASCS software essentially combines these two algorithms. Additional details are given in Appendices A and B. For the analyses found in Sections 5 and 7 we specify the number of iterations for which our RJMCMC algorithm was run (which depended on the observed convergence rate and run time). A single iteration of our RJMCMC algorithm consists of one proposal to change $K$ and ten MCMC updates of the other parameters, i.e., the number of MCMC iterations is ten times greater than the stated number of RJMCMC iterations. In Section 6 we fix $K$ and use MCMC, and thus directly specify the number of MCMC iterations. Our standard approach is to run ten RJMCMC (or MCMC) chains to allow monitoring of convergence, but for simplicity the final results are always computed using a single chain.

As discussed in Section 4.2, having detailed information about the PSF means our estimates are insensitive to the prior on $K$ (see also Section 5.1). Knowledge of the PSF also aids computation in that it limits the number of feasible configurations, meaning the RJMCMC algorithm does not have to jump across many values of $K$. This keeps the number of iterations until approximate convergence comparatively low. Thus, knowledge of the PSF means that, despite the difficulties that are commonly thought to surround mixture models fit with RJMCMC algorithms, our proposed approach is relatively stable and robust. Nonetheless, when the number of sources is clear, MCMC algorithms should be used because they are computationally preferable to RJMCMC algorithms (see Section 6 for an analysis using an MCMC algorithm). In particular, MCMC algorithms are faster per iteration and fewer iterations are needed to obtain enough samples for a given $K$ value of interest. One further challenge is moderate sensitivity to the spectral model, which is the reason why in some applications the gamma spectral model must be replaced by the gamma mixture spectral model introduced in Section 2.4.

To illustrate robustness to the prior on $K$, we simulated data for a one-source (${K}_{{\rm{true}}}=1$) and a ten-source (${K}_{{\rm{true}}}=10$) reality and drew inference for the number of sources under three different settings of the prior mean $\kappa$ (1, 3, and 10). Ten data sets were simulated under each reality, each consisting of images of 20 by 20 spatial units and spectral data (simulated under the single gamma spectral model). We randomly placed the sources in the central 18 by 18 region of the image, avoiding the edges so that source photons are largely contained within the image. The mean number of photons m_j from source j was chosen randomly from the interval 100 to 500, for $j=1,...,K$. The mean total number of photons from the background in each data set, m₀, was set to 400, an average of 1 photon per unit square. The number of photons from source j (or the background) was then simulated from a Poisson distribution with mean m_j, for $j=0,...,K$. Spatial coordinates for the photons were chosen by sampling from the PSF (or the Uniform distribution in the case of the background). We used the King profile density for the PSF; the same PSF is used for analysis of the data sets in Sections 6 and 7.

To complete the data sets, we simulated spectral data under the single gamma spectral model (and from a Uniform distribution in the case of the background). We drew the spectral distribution parameters ${\alpha }_{j}$ (shape) and ${\alpha }_{j}/{\gamma }_{j}$ (rate), $j=1,...,K$, from the Gaussian distributions $N(3,0.{2}^{2})$ and $N(0.005,{0.001}^{2})$, respectively, truncating both distributions to be strictly positive. The resulting spectral parameters are similar to those fitted for the XMM data set in Section 6. An example simulated data set is shown in the left panel of Figure 4. The right panel shows the true spectral distributions for the same data set.

Zoom In Zoom Out Reset image size

For each of the 20 simulated data sets, ten RJMCMC chains were run to assess convergence, but for simplicity only one chain per data set was used in the final analysis.²⁶ The chains were run for 200,000 RJMCMC iterations, the first 100,000 of which formed the convergence period (or burnin) and were discarded. For each data set, the posterior probability of being in state $K=k$ was calculated, using Equation (17), for all feasible values of $k$. Figure 5 summarizes the inference for $K$ under the ten-source (${K}_{{\rm{true}}}=10$, left panels) and one-source (${K}_{{\rm{true}}}=1$, right panels) realities, for $\kappa =1,3$ and 10 (top, middle and bottom panels, respectively). Recall that $\kappa$ is the prior mean number of sources. The 25% and 75% quantiles of the posterior probabilities across the ten data sets are indicated for each value of $K$.

Zoom In Zoom Out Reset image size

Figure 5 shows that, for the ten-source reality, the posterior probability is concentrated around $K=$ 9–11, regardless of which of the three values of $\kappa$ is used. Indeed, the prior probability of ten sources specified by the prior with $\kappa =10$ is nearly 1.25 million times that of the probability specified by the prior with $\kappa =1$. Despite the difference in the prior probability as a function of $\kappa$, the posterior probabilities of $K=10$ are quite consistent; the average (across simulations) differs by only about 0.1 (comparing $\kappa =1$ with $\kappa =10$, see Figure 5). In other words, there is about a 1.25 multiplicative increase in the posterior probability of ten sources when the prior mean is changed from $\kappa =1$ to $\kappa =10$. This modest difference in posterior probability is acceptable as it is unlikely that prior information would allocate the truth 1.25 million to one odds.

There are appreciable differences among the simulated data sets as indicated by the quantiles in Figure 5. This is to be expected because the source positions and intensities are chosen randomly. Some of the simulated data sets have two sources very close to each other, making it hard to determine that they are distinct. In some cases, it is possible to separate these very close sources based on the spectral data (using the full model), i.e., if the spectral data appear to come from two gamma distributions rather than one. However, in other cases it is difficult to separate such nearby sources, even with the spectral data. Indeed, checks confirmed that data sets with sizeable posterior probability at $K=9$ under the full model include overlapping sources that cannot be separated by eye and have similar spectral distributions. Posterior probability at $K$ values of 11 and above appear because chance clusters of photons are sometimes mistaken for separate sources. The precise location of these "ghost" sources, however, is highly erratic across RJMCMC iterations. There is limited evidence for them in the data and thus wide error bars for their "locations" in the posterior distribution.

Inference is also robust to the choice of $\kappa$ under the one-source reality (${K}_{{\rm{true}}}=1$). The posterior mode is clearly $K=1$ for all three values of $\kappa$. Owing to the skewness of the Poisson density, the difference in prior probability of $K=1$ across the different $\kappa$ values is less dramatic than that for $K=10$. When $\kappa =1$ the a priori probability of $K=1$ is around 800 times that when $\kappa =10$. Consequently, the difference in posterior probabilities is also less noticeable. Indeed, the qualitative difference in the posteriors under $\kappa =1$ and $\kappa =10$ is marginal, see Figure 5.

Our key conclusion is that the posterior probability of the true number of sources $K$ seems insensitive to the prior probability assigned to $K$, at least when using the Poisson prior. Consequently, the value of $\kappa$ only needs to be in the region of the true number of sources in order for the fit for $K$ to be reasonable. These conclusions match our intuition that knowing the precise PSF statistically constrains the mixture model sufficiently for the data to drive the fitted values of the parameters. Our simulations are representative of typical data sets, but establishing similar conclusions for smaller datsets may require more studies. A data set could also be larger than those in our simulations, but as $\kappa$ (and $K$) increases, greater Poisson variance means that the absolute deviation of $\kappa$ from the true number of sources has progressively less influence on posterior inferences. (Intuitively, it is more reasonable to a priori suspect 101 sources when there are 110, than to suspect 1 when there are 10.) In our context, prior information typically consists of previous observations, possibly from a different wavelength band. Therefore, it can be assumed that the information is quite reliable and gross prior "misspecification" is unlikely. Clearly, priors other than the Poisson distribution can be considered if a more diffuse prior distribution is desired.

Here we investigate the performance of our model and methods for a range of background intensities, source separations and relative source intensities. We compare the performance of the spatial-only and full models. For simplicity, we simulated data for a two-source (${K}_{{\rm{true}}}=2$) reality. In each simulation, the number of photons from the background and the number from each source were drawn from Poisson distributions with respective means ${m}_{0},{m}_{1},{m}_{2}$. We set ${m}_{2}=1000$ and ${m}_{1}={m}_{2}/r$, for $r=1,2,5,10,50$. We refer to $r$ as the relative intensity of the two sources. To set m₀ and quantify the strength of the simulated background in an astronomically meaningful way we define a source region in terms of the PSF. Specifically, we again use the King profile PSF and define the source region as the region with PSF greater than 10% of its maximum. (The King profile density has no finite moments). We next define $q$ to be the probability that a photon from a source falls within its source region and set the background per source region to be ${m}_{0}={{bqm}}_{2}$, for $b=0.001,0.01,0.1,1$. That is, the mean number of background photons in the faint source region was varied between $1/1000$ and 1 times the mean number of photons from the faint source falling in the same region. As we shall discuss and unsurprisingly, the faint source was difficult to locate in data sets that were simulated with $b=1$ and less so for those simulated with $b=0.001$. Finally, the separation of the two sources was set to be $0.5,1,1.5$, or 2 distance units. These separations can be interpreted using the fact that our source regions are approximately circles of radius 1.

Spectral data was also simulated for source and background photons. An aim of this simulation study aims is to investigate how much using the spectral data improves the fitted parameters. Since sources can only be distinguished by their spectra if their spectra are different, we used different spectra for the two simulated sources; specifically we set ${\alpha }_{1}=3$, ${\gamma }_{1}=600$, ${\alpha }_{2}=6$, and ${\gamma }_{2}=1500$.

In summary, our simulation study consists of a $5\times 4\times 4$ grid of configuration settings ($r=1,2,3,5,10,50$; $b\;=0.001,0.01,0.1,1$; and source separations of $0.5,1,1.5,2$). One hundred data sets were simulated for each of the resulting 80 configurations, and analyzed using first the spatial-only model and then the full model. In particular, for each data set our algorithm was run for 20,000 RJMCMC iterations, the first 10,000 of which formed the convergence period (or burnin) and were discarded.²⁷ The median posterior probability of two sources is shown in Figure 6 for each of the different simulation settings. The left and the right panels correspond to the spatial-only and full models, respectively. We use the median posterior probability across the 100 simulated data sets because in a few simulations the faint source is unusually bright or unusually faint, which noticeably effects the mean posterior probability of two sources. Nevertheless, summaries based on the mean posterior probability are qualitatively very similar, albeit with slightly more noise. We have organized the results by background intensity because in practical applications background is often well determined.

Zoom In Zoom Out Reset image size

In images simulated with relative intensity 50 the posterior probability of two sources tends to be low. This is because $r=50$ corresponds to a faint source intensity of ${m}_{2}=20$, while the brighter source has intensity ${m}_{1}=1000$. Thus, the faint source is typically not bright enough to be distinguished from noise; its photons can be adequately explained as a random cluster formed of photons from the brighter source or the background. In this case the posterior probability peaks sharply at $K=1$. The spatial-only model is more likely to mistake a cluster of background photons for a faint source and therefore, in the case of $r=50$ and small source separation, typically gives slightly higher posterior probabilities of two sources than the full model (but the probabilities are still very small). For less extreme relative intensities, using the full model increases the posterior probability of two sources. The improvement is particularly noticeable for relative intensities 5 and 10, regardless of the background strength. The spectral distribution of source counts reduces the plausibility that the faint source is just a cluster of photons from the background or the bright source. When both sources are bright and reasonably separated both the spatial-only and full models give high posterior probability at $K=2$.²⁸

To fit the source parameters, we fix $K$ at its posterior mode value and use Equation (19). Although the fitted parameters of the bright source are always accurate, those for the faint source may be poor, especially if the posterior mode of $K$ is at 1 or if "ghost" sources have appreciable posterior probability. The accuracy of the faint source's fitted parameters essentially follows the pattern seen in Figure 6. When the real faint source is very weak or located too close to the bright source, then a fitted second source (when the posterior mode of $K$ is greater than 1) is likely to be a "ghost" consisting mainly of a cluster of photons from the background or the bright source. In which case, its fitted parameters bear little resemblance to those of the true faint source. This is illustrated in Figure 7, which shows the mean (conditional on $K=2$) posterior locations of the two sources for all 100 data sets under each configuration of simulation settings. Crosses indicate the true locations of the sources. The mean posterior locations of the bright source (red dots) are not always visible in the plots because they are often in the middle of the red crosses. The location of the bright source becomes slightly harder to fit as the intensity of the faint source increases. (This is at least partly because the background intensity is proportional to the faint source intensity.) The size of the dots indicate the posterior probability of two sources.

Zoom In Zoom Out Reset image size

The full model again yields more accurate fits. The fitted locations of the faint source (blue dots) center around its true location (blue crosses) for $r\leqslant 10$, even when the source separation is small. For the spatial-only model there is more scatter. Under both models, when $r=50$ we can see that many of the fitted faint source locations correspond to spurious clusters of photons surrounding the bright source. As the separation increases some of the fitted faint source locations are halfway between the true locations of the two sources. This occurs when the posterior distribution of the faint source x-coordinate is bimodal, a spurious cluster of photons and the real faint source both being supported as possible second sources. In an actual analysis this bimodal behavior would be apparent from inspection of the posterior draws of the source location. For r = 50 and large separation, the full model sometimes accurately fits the faint source location, but the spatial-only model never does. The behavior of the other fitted parameters follows the same pattern illustrated in Figure 7 because the fitted source locations indicate how well photons are allocated to the correct source. This is confirmed by inspecting tables of the mean (or median) squared error of each parameter (not shown).

The number of Monte Carlo samples used in estimating the mean posterior locations (conditional on $K=2$) is determined by the posterior probability of two sources, and thus is indicated by the size of the dots. Very small dots may have non-negligible Monte Carlo error, i.e., the true posterior mean location (conditional on $K=2$) may be somewhat inaccurately approximated. This is because applying Equation (18) for each parameter does not accurately compute the mean of $p({{\rm{\Theta }}}_{K}| K,x,y,E)$ for values of $K$ that have low posterior probability.²⁹ However, in practice, when the number of sources is unknown, it makes sense to only consider values of $K$ with relatively high posterior probability. Furthermore, one typically checks the level of Monte Carlo error for the values of $K$ of interest, by running multiple chains. Large variation in the parameter estimates across the chains indicates high Monte Carlo error. In which case, one should run the chains longer in order to obtain a larger Monte Carlo sample.

For both models, ten RJMCMC chains were run for 150,000 RJMCMC iterations, the first 100,000 of which formed the convergence period (or burnin) and were discarded. The posterior distribution of the number of sources, under the spatial and extended full models, is displayed in Figure 12. The mode of both posteriors is at 14. However, the spatial-only model shows slightly more uncertainty, and some support for 15 sources.

Zoom In Zoom Out Reset image size

As mentioned in Section 4, $K$ determines the number and meaning of the other model parameters and therefore we must condition on a value of $K$ to draw meaningful inferences for them from the RJMCMC output. Figure 11 shows 90% posterior credible regions (blue) for the locations of the sources under the two models, given $K=14$. Each credible region shows an area which has 0.9 posterior probability (given $K=14$) of containing the location of the relevant source, i.e., an integral of the posterior distribution of the source location (given $K=14$) over this area would evaluate to 0.9. The credible regions look to be similar under the two models. The estimated relative intensities also appear in Figure 11 and are also similar, but are slightly lower under the spatial-only model for most sources. This is due to a higher estimate of the relative background intensity under the spatial-only model (0.0053 versus 0.0006 under the extended full model³³ ). Table 5 gives the the posterior mean fit of the source locations and relative intensities under the extended full model for $K=14$. The detected sources are also matched to the source catalog from the Chandra Orion Ultradeep Project (COUP; Getman et al. 2005).

Table 5.  Extended Full Model Fit for the Chandra Observation in Figure 11

COUP #	${\mu }_{{jx}}$		${\mu }_{{jy}}$		Relative Intensity (%)
732	4054.42	(4054.41, 4054.43)	4149.45	(4149.44, 4149.46)	34.59	(34.16, 35.03)
745	4052.83	(4052.81, 4052.84)	4140.67	(4140.66, 4140.68)	28.11	(27.71, 28.51)
689	4069.93	(4069.91, 4069.94)	4175.93	(4175.91, 4175.94)	14.10	(13.79, 14.40)
724	4058.57	(4058.56, 4058.59)	4176.73	(4176.71, 4176.74)	11.43	(11.16, 11.71)
744	4051.53	(4051.50, 4051.55)	4147.57	(4147.55, 4147.60)	7.41	(7.14, 7.68)
765	4045.40	(4045.35, 4045.46)	4181.20	(4181.15, 4181.25)	1.42	(1.32, 1.53)
649	4088.16	(4088.08, 4088.24)	4165.95	(4165.87, 4166.03)	0.57	(0.50, 0.63)
766	4045.36	(4045.27, 4045.45)	4155.18	(4155.10, 4155.25)	0.77	(0.68, 0.87)
788	4043.48	(4043.36, 4043.61)	4155.74	(4155.64, 4155.84)	0.56	(0.47, 0.64)
682	4072.11	(4072.01, 4072.21)	4181.12	(4181.03, 4181.22)	0.46	(0.39, 0.52)
640	4091.73	(4091.53, 4091.92)	4137.42	(4137.26, 4137.59)	0.13	(0.10, 0.16)
664	4081.43	(4081.22, 4081.63)	4159.41	(4159.21, 4159.61)	0.11	(0.08, 0.14)
665	4082.84	(4082.67, 4083.02)	4137.28	(4137.14, 4137.43)	0.15	(0.12, 0.19)
779	4044.39	(4043.86, 4044.60)	4140.72	(4140.43, 4140.90)	0.14	(0.09, 0.18)
Background	⋯	⋯	⋯	⋯	0.06	(0.01, 0.10)

Note. Posterior mean locations and relative intensities (as percentages), with 68% intervals indicated.

Download table as:  ASCIITypeset image

Other observations of Orion suggest that the source circled (in green) in the right panel of Figure 11 is a genuine source. Its location is more uncertain than other sources because it is more difficult to detect. Indeed, with an estimated intensity between 13 and 25 counts, this source is at the edge of detectability of local detection methods, particularly since the estimate of the local background in such methods would be high due to contamination from nearby bright sources. Thus, we expect that more basic approaches would either have failed to find this source, or would only find it by rendering their detection threshold to a point where spurious detections became problematic. Indeed, the reason the spatial-only model gives non-negligible weight to 15 sources (see Figure 12) is that it tends to split sources into two. The problem is that a single empirical PSF may exhibit chance variations that appear to be evidence for multiple PSFs. The spatial-only model also mistakes clusters of background photons for sources. The locations and spectra of these spurious sources show considerable posterior variability. Although any particular instance has low probability, there are multiple instances that together create erroneous support for an additional source. The main advantage of using the spectral information, in this example, is that it mitigates these issues, leading to a greater certainty that there are really 14 sources. Additionally, under the extended full model, the standard deviations of the parameters are almost invariably slightly smaller.

The extended full model only captures the basic shape of the source spectra and we now illustrate how detailed follow-up spectral analysis can incorporate probabilistic event allocations. We perform this analysis for the two overlapping sources that are enclosed in the red box in the right panel of Figure 11 (COUP sources #732 and #744; Getman et al. 2005). Their estimated relative intensities are 0.3459 and 0.0741 under the extended full model. This is a good example to test the probabilistic event allocations, since the sources are close together (separation $\approx 1\buildrel{\prime\prime}\over{.} 7$), each have sufficient counts for a useful spectral fit ($\approx 4350$ and $\approx 910$ counts between $0.5-7$ keV for the bright and faint sources, respectively), and one source is substantially weaker than the other.

As described in Section 2.3, ${s}_{i}$ indicates the source (or background) number associated with photon i. These are unknown parameters (or latent variables) that are updated at each iteration of the RJMCMC sampler. The variability in ${s}_{i}$ indicates the uncertainty in the source of photon i (due to the PSF and uncertainty in the source parameters). We can account for this uncertainty by conducting many spectral analyses, each according to a sampled photon allocation (i.e., sampled values of ${s}_{i}$), and combining the results. We focus on photons with spatial location in the red box in Figure 11 (right panel) and to values of ${s}_{i}$ sampled conditional on $K=14$. Since we are only interested in COUP sources #732 and #744 we ignore any photons that are attributed to one of the other sources (in a given allocation). (The photons in the red box in Figure 11 are attributed to one of the other sources only rarely.)

Based on the photon allocations, we construct a sample of 1000 simulated spectral data sets for both sources, constructed from photon allocations based on every $10\mathrm{th}$ iteration of the RJMCMC algorithm that sets $K=14$ (up to the $10,000\mathrm{th}$ RJMCMC iteration that sets $K=14$). The variability in the source counts across the 1000 iterations is $\pm 17$ for both the bright and faint sources. The specific photons that are allocated to each source also varies, even when the total source counts do not. Each individual spectrum is fit with an absorbed single temperature thermal model (xsphabs*xsapec in CIAO/Sherpa v4.6) fitting the absorption column (${N}_{{\rm{H}}}$), temperature (kT), metallicity (Z), and normalization. A pile-up correction is needed for all spectra for the bright source since the measured count rate of 0.7 counts frame⁻¹ is higher than the threshold at which pile-up becomes significant ($\approx 0.3$ counts frame⁻¹). We use the jdpileup model in Sherpa, fitting the grade migration parameter α and the pile-up strength parameter f (Davis 2001). We call the entire collection of spectral fits the disentangled analysis.

For comparison, we also carry out a spectral analysis of the sources based on a naïve allocation of photons that collects events from within 1'' of the fitted location of each source and assumes that there is no contamination from the other source. The only difference in the spectral model for the naïve and disentangled analyses is in how the effective areas are defined. In the case of the naïve analysis, a correction is made post-facto to the normalization based on how much of the source is expected to be included within the 1'' source photons extraction radius. In the disentangled analysis, the assumed extraction radius for the spectra with allocated events is set to be $2\buildrel{\prime\prime}\over{.} 5$ and the subsequent correction is negligible.

The results of the spectral fits to the disentangled spectra are shown as histograms of best-fit values for ${N}_{{\rm{H}}}$, kT, Z, and model flux computed for each of the 1000 spectra, see Figure 13. In several cases, a bimodal distribution is apparent. This suggests that a multi-temperature component spectrum would be a better fit. The separation of the modes, however, is generally too small to be picked up by typical multi-temperature model fits. Not shown are the pile-up parameters for the bright source, which are consistent between the naïve and disentangled analyses ($(\alpha ,f)=(0.6,0.93)$ for naïve, and $(0.53,0.89)$ for the disentangled spectra), though the former indicates that the pile-up strength is slightly higher. This is to be expected, since the naïve analysis is carried out for photons in the core of the PSF, where naturally pile-up is most significant. The disentangled spectra include photons from the wings, thus reducing the strength of pile-up effects and decreasing the correction needed to the source flux by about 60%.

Zoom In Zoom Out Reset image size

The spread in the histograms in Figure 13 indicates the uncertainty in the best-fit values due to uncertainty in the allocation of photons. The best-fit values from the naïve calculation are shown as solid red vertical lines. The dashed red vertical lines give 68% intervals indicating the statistical errors, due to randomness in the photons emitted and detected, under the naïve analysis. These statistical errors do not account for uncertainty in the photon allocations. The histograms, on the other hand, represent only errors due to uncertainty in the photon allocations, but do not account for statistical errors (due to randomness in photon emission and detection). Because the two sources of error are independent, and because we expect the statistical errors for the disentangled analyses to be similar to those for the naïve analysis, the total errors could be represented by a perturbation of the histograms with σ equal to the statistical errors from the naïve analysis. For these data, with the exception of flux (panels (g) and (h) of Figure 13), the statistical errors dominate the errors due to uncertainty in the photon allocation. Despite this, the disentangled analysis provides reasonable evidence that the absorption column of the faint source (panels (b) of Figure 13) and the flux of the two sources (panels (g) and (h) of Figure 13) are different from the best-fit values under the naïve analysis.

The variability of the true parameters around each of the best fit values recorded in the histograms is expected to be similar to that indicated for the naïve fit. However, we did not calculate these uncertainties because of the large computational cost. For these data, with the exception of flux (panels (g) and (h) of Figure 13), the variability in the true spectral parameters around the best fit values is likely larger than the uncertainty in the best fit values (due to the uncertainty in the allocation of photons). Despite this, the disentangled analysis provides reasonable evidence that the absorption column of the faint source (panel (b) of Figure 13) and the flux of the two sources (panels (g) and (h) of Figure 13) are different to the naïve analysis best-fit values.

Overall, the naïve analysis best-fit values for the fainter source are in greater disagreement with the disentangled analysis than those for the bright source. This is to be expected, since in the naïve analysis, the contamination of the fainter source by the brighter source is larger. Our algorithm effectively removes this contamination. This causes the spectral fit parameter values to change and the measured source flux of the fainter source to decrease. In summary, the observed changes to the spectral model parameters are as would be expected when contamination is reduced and the data quality is improved.