Bayesian mixture modelling of the high-energy photon counts collected by the Fermi Large Area Telescope

Abstract

Identifying as yet undetected high-energy sources in the $γ$ -ray sky is one of the declared objectives of the Fermi Large Area Telescope (LAT) Collaboration. We develop a Bayesian mixture model which is capable of disentangling the high-energy extra-galactic sources present in a given sky region from the pervasive background radiation. We achieve this by combining two model components. The first component models the emission activity of the single sources and incorporates the instrument response function of the Fermi $γ$ -ray space telescope. The second component reliably reflects the current knowledge of the physical phenomena which underlie the $γ$ -ray background. The model parameters are estimated using a reversible jump MCMC algorithm, which simultaneously returns the number of detected sources, their locations and relative intensities, and the background component. Our proposal is illustrated using a sample of the Fermi LAT data. In the analysed sky region, our model correctly identifies 116 sources out of the 132 present. The detection rate and the estimated directions and intensities of the identified sources are largely unaffected by the number of detected sources.

Keywords

bivariate exponential distribution γ-ray photon instrument response function mixture model reversible jump MCMC

1 Introduction

The Large Area Telescope (LAT) is the primary instrument onboard the Fermi γ -ray space telescope (Atwood et al., 2009). The main objective of the Fermi LAT Collaboration (http://fermi.gsfc.nasa.gov/) is to provide a catalogue of sources present in the $γ$ -ray sky. Since the Fermi LAT began its survey on 11 August 2008, the Collaboration has released several catalogues which are accessible on-line. The latest version of the catalogue, the Third High-Energy Source Catalog (3FHL), lists the $3 034$ $γ$ -ray sources which have been detected until now (Ajello et al., 2017). The data collected by the Fermi LAT are so-called events in the high-energy band. Each event corresponds to a high-energy photon, or $γ$ -ray, which interacts with the space telescope. At every hit, several measurements are taken. We will use the sky directions of the incoming photons, their energy and the quality of the reconstructed events, the so-called event type, to pinpoint high-energy extra-galactic emitting sources.

The main difficulty of identifying high-energy sources in astrophysical count maps is to separate the signal stemming from the sources from the background component which spreads over the whole sky. The method used to detect the sources listed in the 3FHL catalogue consists of two steps. The first step uses an image analysis technique based on wavelet transforms. Broadly speaking, two wavelet transforms (see, e.g., Damiani et al., 1997; Starck and Pierre, 1998; Ciprini et al., 2007) translate the incoming signal into a shape which is suitable for detecting a large variety of sources. These can range from point sources to highly extended sources, and their signal may be smeared by the detector of the LAT. The second step quantifies the statistical significance of a possible source detection using the likelihood ratio statistic (Mattox et al., 1996).

In this article we abandon the commonly used two-step technique described above and consider a Bayesian mixture model instead. We follow previous work by Jones et al. (2015) on the identification of X-ray sources and extend their model into two directions. First, the uncertainty of the measurements is modelled using the instrument response functions provided by the Fermi LAT Collaboration. Second, we translate the simulation-based background model developed by the same Collaboration into a workable parametric formulation. As in practice, we do not have any knowledge about the origin of a single photon, these two model components are combined into a finite mixture model with an unknown number of components. The number of sources, their locations and relative intensities, and the background component are estimated using a tailored reversible jump MCMC algorithm.

Our newly developed model was tested on a collection of $γ$ -ray photons recorded by the Fermi LAT during a period of around 7.2 years of operation with promising results. We were able to identify almost 90% of the sources present in the analysed sky region and to reconstruct their directions with high precision. A careful sensitivity analysis revealed that the parameter estimates of the identified sources remain stable once a sufficiently large number of components is specified. Although the purpose of the analysis is to disentangle high-energy emitting sources from the prominent background radiation, our model can easily be adapted to other types of highly contaminated data in the two-dimensional space where the set of locations at which a spatial process is observed are of interest. As such, our model belongs to the wider class of models known in spatial statistics as Spatial Point Processes (SPPs). The approach we follow differs from the better-known branch of spatial statistics called geostatistics (or point-referenced data) which tries and reconstructs a spatial process that varies continuously in space but is observed only at a few points. A general reference for SPPs is Gaetan and Guyon (2010, Chap. 3), which also well explains the difference between the two approaches.

The article, is organized as follows. Section 2 gives a general overview of the Fermi LAT and the datasets considered in our work. The description of the proposed statistical model and of how we fit it are presented in Sections 3 and 4. Section 5 outlines the main results of our analysis, which are then discussed in Section 6.

2 The Fermi LAT data

The Fermi LAT is a wide field-of-view pair-conversion telescope which covers the energy range from 20 MeV to more than 300 GeV (Atwood et al., 2009). Figure 1 reproduces the whole $γ$ -ray sky in galactic coordinates and at energies larger than 1 GeV according to the data collected by the Fermi LAT over a period of five years. We can clearly identify extra-galactic point and small-area $γ$ -ray sources. The brighter the colour, the brighter are the emitting sources. The brilliant horizontal stripe in the centre of the figure corresponds to the Galactic plane which is mainly formed by the $γ$ -ray photons emitted by our galaxy, the Milky Way, whose centre hosts a supermassive black hole. Indeed, the count map shown in Figure 1 is a mixture of photons which are emitted from both, known and as yet undetected high-energy sources and the diffuse γ -ray background.

Figure 1:

Whole sky map at γ -ray wavelengths and energies larger than 1 GeV based on data accumulated by the LAT over a period of five years of operation (Image Credit: NASA/DOE/Fermi LAT Collaboration). The region framed in white represents the area analysed in this article

Our primary dataset consists of the $26 208$ events recorded by the Fermi LAT during a period of around $7.2$ years of operation and which fall in the subregion framed in white in Figure 1. This region covers the portion of the sky which lies between $[- 10^{\circ}, - 80^{\circ}]$ of galactic longitude and $[10^{\circ}, 80^{\circ}]$ of galactic latitude, respectively. The detected γ -ray photons cover an energy range between $10$ GeV and $1 000$ GeV. Figure 2 plots the nonparametric kernel estimate of the photon density distribution. The various spikes most likely correspond to known and as yet unrevealed high-energy emitting sources. The irregularly shaped $γ$ -ray background which spreads over the entire area is clearly visible.

Figure 2:

Nonparametric kernel estimate of the photon density based on the $γ$ -ray count map accumulated by the Fermi LAT over a period of 7.2 years. The map is expressed in galactic coordinates and refers to the sky region $[- 10^{\circ}, - 80^{\circ}] \times [10^{\circ}, 80^{\circ}]$ , that is, to the area framed in white in Figure 1. The spikes represent potential high-energy extra-galactic emitting sources. The irregularly shaped $γ$ -ray background is clearly visible

Several measurements are taken by the $γ$ -ray detector onboard the Fermi spacecraft for every collected event. Our purpose is to model the directions, expressed in galactic coordinates, of the incoming $γ$ -ray photons and to use the gained insight to identify known and possibly new extra-galactic high-energy emitting sources. We will also consider the energy content of the incoming photons and the associated event type. The latter variable characterizes the quality of the detection, that is, it specifies how accurately the detector identified an incoming photon. There exist four different classes of event types which, in increasing order of reliability, are called PSF0, PSF1, PSF2 and PSF3. We will integrate both pieces of information, energy content and event type, into the model which reconstructs the direction of the high-energy photons. The details are given in Section 3.1.

As mentioned in the introduction, the detection of high-energy emitting sources is hampered by the presence of the photon emitting background which spreads over the whole sky. Jones et al. (2015) model it with a uniform distribution, which is fine for the X-ray photons considered in their work, but inadequate to capture the rather irregular background structure of $γ$ -ray photons. The contaminating effect in particular dominates the galactic part of the sky map and involves non negligible difficulties when we try to discover new sources whose directions lie in this area. This is why in this article we restrict our attention to subregions of the sky whose latitudes lie above $10^{\circ}$ . In Section 3.2, we will present a suitable model to describe the behaviour of the $γ$ -ray background. This model is based on a careful analysis of the intensity maps provided by the Fermi LAT Collaboration. The details are given in the Appendix.

3 Model definition

Denote by $x_{i} \in (180^{\circ}, - 180^{\circ})$ and $y_{i} \in (- 90^{\circ}, 90^{\circ})$ the observed galactic longitude and galactic latitude of photon $i$ with $i = 1, \dots, n$ . Let furthermore $E_{i} \in R^{+}$ be the energy content of an observed photon. Each photon was either emitted by an extra-galactic source or is part of the diffuse background radiation which spreads over the entire observed area. We will hence model the directions of the photons using two different models depending on which situation occurs. These are the single source model $f_{j} (x_{i}, y_{i} ∣ E_{i}, μ_{j})$ , described in Section 3.1, which characterizes how photons scatter around source $j$ , and the background model $g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ)$ of Section 3.2. As in practice, we do not have any knowledge about the origin of a single photon, these two model components are combined into a finite mixture model of the form

f (x_{i}, y_{i} ∣ E_{i}, θ_{K}, K) = \sum_{j = 1}^{K} ω_{j} f_{j} (x_{i}, y_{i} ∣ E_{i}, μ_{j}) + ω_{0} g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ) .

(3.1)

The vector $ω = (ω_{0}, ω_{1}, \dots, ω_{K})$ contains the mixing proportions $ω_{j} \in (0, 1)$ , with $\sum_{j = 0}^{K} ω_{j} = 1$ . We can interpret these as the relative intensities of the various identified sources $ω_{i}$ , for $i = 1, \dots, K$ , and as the relative intensity $ω_{0}$ of the background. That is, they determine the proportion of photons emitted from an identified source and the proportion of photons which are to be assigned to the background. The vector $μ_{j} = (μ_{jx}, μ_{jy}) \in (180^{\circ}, - 180^{\circ}) \times (- 90^{\circ}, 90^{\circ})$ contains the direction of the $j$ th source, while $(σ_{x}, σ_{y}, δ)$ parametrize the background distribution, with $σ_{x}$ and $σ_{y}$ positive values and $δ \in (0, 1)$ . Further details of the two components of model (3.1) follow in Sections 3.1 and 3.2.

Our model is hence characterized by a set $θ_{K} = {ω, {μ_{j}}_{j = 1}^{K}, σ_{x}, σ_{y}, δ}$ of $3 (K + 1)$ parameters which need to be estimated. Recall, furthermore, that the number $K \in N$ of undetected sources is itself supposed to be unknown and needs to be estimated. We will make inference on $(θ_{K}, K)$ using Bayes’ rule. The a priori available information on these parameters is summarized in Section 3.3.

3.1 The single source model

Suppose that photon $i$ was emitted by the $j$ th high-energy source with sky direction $μ_{j}$ . From the physical point of view, the emission activity of the source is point-like. However, we cannot express it by a point mass in $μ_{j}$ as this neglects all instrumental aspects which need to be accounted for; see Appendix for the details. We will model the direction $(x_{i}, y_{i})$ of photon $i$ which is emitted by source $j$ by a so-called point spread function (PSF). This is a class of functions which describe the instrumental response of optical measurement devices. In case of the Fermi LAT, the PSF describes how the incoming photons spatially scatter around $μ_{j}$ , the direction of the emitting source, as a function of their energy content $E_{i}$ and of the geometry of the detector. In particular, we will adopt the PSF developed by the Fermi LAT Collaboration (Ackermann et al., 2013).

Let $t_{d} (\cdot ∣ b, C)$ be the density of a bivariate Student t distribution with location parameter $b$ , scale matrix $C$ and $d$ degrees of freedom. The density $f_{j} (\cdot ∣ \cdot)$ in Formula (3.1) becomes

f_{j} (x_{i}, y_{i} ∣ E_{i}, μ_{j}) = p_{i} t_{ν_{i, core}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, core}^{2} I_{2}) + (1 - p_{i}) t_{ν_{i, tail}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, tail}^{2} I_{2}),

(3.2)

where $I_{2}$ is the $2 \times 2$ identity matrix. That is, the Fermi LAT PSF is a mixture of two bivariate Student t densities $t_{ν_{i, core}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, core}^{2} I_{2})$ and $t_{ν_{i, tail}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, tail}^{2} I_{2})$ , respectively. Both distributions are centred at $μ_{j}$ , but differ in their scales ${\tilde{σ}}_{i}^{2}$ and their degrees of freedom $ν_{i}$ to characterize the scattering for high and low quality detections. The scale parameters and degrees of freedom depend on the energy content $E_{i}$ of the photon and on its event type. They are calculated as given below using the information provided by the Fermi LAT Collaboration for each single detection.

For each photon $i$ define the scale factor

S_{p} (E_{i}) = \sqrt{{\{c_{0, i} (\frac{E_{i}}{100 MeV})^{- 0.8}\}}^{2} + c_{1, i}^{2}},

where the constants $c_{0, i}$ and $c_{1, i}$ depend on the event type of the photon and are listed in Table 1. Each photon is furthermore characterized by two additional sets of constants, $(σ_{i, core}, γ_{i, core}, N_{i, core})$ and $(σ_{i, tail}, γ_{i, tail}, N_{i, tail})$ , whose values depend again on the event type, but also on the energy content $E_{i}$ of the photon and on how steeply the photon hits the surface of the LAT. These constants can be retrieved from the Science Tools repositories (https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/) of the Fermi LAT project. We define the parameters $({\tilde{σ}}_{i, core}^{2}, ν_{i, core})$ as follows

{\tilde{σ}}_{i, core}^{2} = σ_{i, core}^{2} S_{p}^{2} (E_{i}) \frac{γ_{i, core}}{γ_{i, core} - 1}, ν_{i, core} = 2 (γ_{i, core} - 1) .

Table 1:

Values for $c_{0}$ and $c_{1}$ for the four types of event classes detected by the LAT

Event type	$c_{0}$ (Radians)	$c_{1}$ (Radians)
PSF0	$1.53 \times 10^{- 1}$	$5.70 \times 10^{- 3}$
PSF1	$9.64 \times 10^{- 2}$	$1.78 \times 10^{- 3}$
PSF2	$7.02 \times 10^{- 2}$	$1.07 \times 10^{- 3}$
PSF3	$4.97 \times 10^{- 2}$	$6.13 \times 10^{- 4}$

The same expression holds for $({\tilde{σ}}_{i, tail}^{2}, ν_{i, tail})$ , only that $(σ_{i, core}, γ_{i, core})$ is replaced by $(σ_{i, tail}, γ_{i, tail})$ . Finally, we define the mixture weights (https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html) as $p_{i} = σ_{i, core}^{2} / (σ_{i, core}^{2} + N_{i, tail} σ_{i, tail}^{2})$ .

To sum up, the single source model (3.2) depends on the unknown source direction $μ_{j}$ and on a number of event-specific constants, provided by the Fermi LAT for each single photon, which are directly linked to the most complex detection mechanism of the space telescope.

3.2 The background model

The emission activity of the $γ$ -ray background is characterized by the density function $g_{b} (\cdot ∣ \cdot)$ which appears in the second part of Model (3.1). Figure 3 plots the marginal distributions of the galactic coordinates of the $γ$ -ray photons emitted by the background as simulated using the Fermi LAT Science Tools; see Appendix for the details. The figure considers the second quadrant spanned by $[0^{\circ}, - 180^{\circ}] \times [0^{\circ}, 90^{\circ}]$ and an observation period of $7.2$ years. The vertical dashed lines delimit the sky region analysed in Section 5. The values of both the longitude and the latitude exhibit an exponential decay when getting away from the galaxy centre. Simply using two independent exponential distributions does not properly fit the background photon emission. We have to account for the dependence between the longitude and the latitude which is clearly visible in the contour plot of Figure 4.

Figure 3:

Histograms for the marginal distributions of the sky directions for the photons emitted by the $γ$ -ray background as simulated by the Fermi LAT Science Tools for the second quadrant spanned by $[0^{\circ}, - 180^{\circ}] \times [0^{\circ}, 90^{\circ}]$ and an observation period of $7.2$ years. Left: galactic longitude. Right: galactic latitude. The vertical dashed lines delimit the sky region analysed in Section 5

Figure 4:

Top: Contour plot of the simulated background photon density for the second sky quadrant spanned by $[0^{\circ}, - 180^{\circ}] \times [0^{\circ}, 90^{\circ}]$ . Bottom: Image plot of the approximate count map reconstructed using the 30 GIE intensity maps provided by the LAT

In this article, we propose to use the bivariate exponential distribution defined by Gumbel (1960) to model the complex structure of the background component in a simple though theoretically coherent way. We define $g_{b} (\cdot ∣ \cdot)$ as

\begin{matrix} g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ) & = & σ_{x} σ_{y} exp {- σ_{x} (- x_{i}) - σ_{y} y_{i} - δ σ_{x} σ_{y} (- x_{i}) y_{i}} \times \\ \times [{1 + δ σ_{x} (- x_{i})} {1 + δ σ_{y} y_{i}} - δ], \end{matrix}

(3.3)

where $σ_{x} > 0$ and $σ_{y} > 0$ are two scale parameters, while $δ \in (0, 1)$ expresses the correlation between the longitude and the latitude. It is straightforward to show that for $δ = 0$ , the model reduces to the product of two independent exponential distributions with rates $σ_{x}$ and $σ_{y}$ . The motivations behind this model choice are given in Appendix.

Differently to what holds for the PSF, where all scale and shape parameters are fixed and provided by the LAT, here the background parameters need to be estimated.

3.3 Prior belief

In Section 4, we will estimate Model (3.1) using Bayes’ rule. The prior distributions for the model parameters $θ_{K} = {ω, {μ_{j}}_{j = 1}^{K}, σ_{x}, σ_{y}, δ}$ and $K$ are set using the available information on the various physical processes and the current knowledge on the already catalogueued sources.

The first element in $(θ_{K}, K)$ are the mixing proportions $ω$ . Conditional on the number $K$ of sources present in the map, we choose a Dirichlet prior

ω ∣ K \sim Dir (1, \dots, 1)

with concentration parameters equal to 1. This corresponds to assuming a flat distribution over the simplex of dimension K+1. Though some prior belief on the intensity of the background radiation $ω_{0}$ is in principle available, we let the model free to specify the single contributions of the sources and of the background.

The only unknown parameters of the PSFs are the galactic coordinates $μ_{j}$ of the putative sources, for $j = 1, \dots, K$ . Since the universe contains a potentially infinite number of high-energy emitting phenomena at any possible direction, we use as prior for the $μ_{j}$ ’s a uniform distribution

μ_{jx} \sim U (180, - 180), μ_{jy} \sim U (- 90, 90)

which extends over the whole map. Two exponential distributions with rate 1 are assigned as priors to the two scale parameters $σ_{x}$ and $σ_{y}$ of the background model

σ_{x} \sim Exp (1), σ_{y} \sim Exp (1), δ \sim U (0, 1),

(3.4)

while a uniform distribution over $(0, 1)$ will be adopted for the correlation parameter $δ$ . For the ease of computation, these three parameters were further transformed to the real domain as $log σ_{x}$ , $log σ_{y}$ and logit(δ).

Finally, our model requires to set a prior distribution for the unknown number of sources $K$ which we may detect in the map. This is the most tricky choice. We assume a truncated Poisson distribution

K \sim tP (κ ∣ κ_{\min}, κ_{\max}),

where the hyper-parameters $κ$ , $κ_{\min}$ and $κ_{\max}$ are set according to the specific sky area which is being analysed using as scientific input the number of catalogue sources in the region of interest (Ajello et al., 2017). As we will discuss later on, the directions and relative intensities of the detected sources remain stable if a sufficiently large number $K$ of sources is considered.

4 Bayesian model fitting

We will follow the Bayesian paradigm and conduct inference on the model parameters $θ_{K} = {ω, {μ_{j}}_{j = 1}^{K}, σ_{x}, σ_{y}, δ}$ and $K$ by exploring their posterior distribution. This distribution is derived in Section 4.1 up to a normalizing constant as the product of the likelihood function of Model (3.1) and of the prior densities and probabilities specified in Section 3.3. Section 4.2 outlines the reversible jump MCMC algorithm used to fit our model.

4.1 The posterior distribution

A common strategy to write down the likelihood function of a finite mixture model is to highlight the latent group variable which pinpoints the component of the mixture which generated a specific observation (Richardson and Green, 1997). With respect to Model (3.1), write $Z_{i} = 0$ if the generic photon $i$ comes from the background with probability $ω_{0}$ and $Z_{i} = j$ if the photon was emitted by the $j$ th source with probability $ω_{j}$ , for $j \in {1, \dots, K}$ . By construction, we have that $\sum_{i = 1}^{n} 1 (Z_{i} = j) = n_{j}$ , where $1 (\cdot)$ is the indicator function and $n_{j}$ the number of photons emitted by source $j$ . It follows that conditional on the number $K$ of sources,

(n_{0}, \dots, n_{K}) ∣ K \sim Mult (n, ω)

is a multinomial distribution indexed by the intensity vector $ω$ .

Let now $x = (x_{1}, \dots, x_{n})$ and $y = (y_{1}, \dots, y_{n})$ be the vectors of observed longitudes and latitudes and $z = (z_{1}, \dots, z_{n})$ represent the unobserved latent class values. The full data likelihood function

\begin{matrix} L (θ_{K}, K ∣ x, y, z) & = p (x, y ∣ z, θ_{K}, K) p (z ∣ θ_{K}, K) \\ = [\prod_{i : z_{i} = 0} g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ) \prod_{j = 1}^{K} \{\prod_{i : z_{i} = j} f_{j} (x_{i}, y_{i} ∣ μ_{j})\}] \prod_{j = 0}^{K} ω_{j}^{n_{j}} \end{matrix}

(4.1)

is given by the joint density function of $(x, y, z)$ . By Bayes’ theorem, the posterior distribution of the model parameters $(θ_{K}, K)$ is

π (θ_{K}, K ∣ x, y, z) \propto L (θ_{K}, K ∣ x, y, z) π (θ_{K}, K) .

(4.2)

Here,

π (θ_{K}, K) = π (θ_{K} ∣ K) π (K),

represents the joint prior distribution obtained by assuming

π (θ_{K} ∣ K) = π (ω ∣ K) π ({μ_{j}}_{j = 1}^{K} ∣ K) π (σ_{x}, σ_{y}, δ),

(4.3)

that is, that the parameters $ω$ , ${μ_{j}}_{j = 1}^{K}$ and ${σ_{x}, σ_{y}, δ}$ of $θ_{K}$ are conditionally independent given $K$ . See Section 3.3 for the elicitation of the individual densities in (4.3) and for the prior $π (K)$ .

4.2 The reversible jump MCMC algorithm

We implemented a two-step procedure for estimating the parameters $(θ_{K}, K)$ which combines Gibbs sampling with a Metropolis step. In turn, values are generated from the full conditional distributions with density $π (θ_{K} ∣ K, x, y, z)$ and probability $π (K ∣ θ_{K}, x, y, z)$ . At Step 1, the number of sources $K$ is kept fixed, while at Step 2 the value of $θ_{K}$ is given. That is, we alternately estimate $θ_{K}$ for the given number of sources and then update the number of sources for a given set of parameters.

The following two paragraphs describe a single iteration of the algorithm; these need to be repeated for a suitably large number $T$ .

4.2.1 Estimate model parameters for given $K$ (Step 1)

Assign the $n$ observed $γ$ -ray photons to the $K + 1$ mixture components with allocation probabilities proportional to $ω_{j} f_{j} (x_{i}, y_{i} ∣ μ_{j})$ for source $j = 1, \dots, K$ , and $ω_{0} g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ)$ for the background, for $i = 1, \dots, n$ . The model parameters $θ_{K}$ are updated following the order in which they are listed in the joint prior density (4.3).

Update the vector of mixture weights $ω$ by drawing from the Dirichlet distribution Dir(n_0+1,\dots,n_K+1) whose hyper-parameters depend on the numbers ( $n_{0}$ , $n_{1}$ , …, $n_{K}$ ) of photons assigned to the $K$ sources and to the background. This step is straightforward given that the Dirichlet distribution is conjugate to the Multinomial.

Update the location parameters $μ_{j}$ by applying a double data augmentation strategy to the PSF of the single source model. This strategy allows us both to avoid the mixture representation of Model (3.2) and to rewrite the Student t distribution as a product of a bivariate Gaussian and a gamma distribution (Section 18.3 Gelman et al., 2013).

Let $S_{i} \sim Be (p_{i})$ be the random indicator variable which determines whether photon $i$ comes from the $core$ density $t_{ν_{i, core}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, core}^{2} I_{2})$ or from the $tail$ density $t_{ν_{i, tail}} (x_{i}, y_{i} ∣ μ_{j}, {\tilde{σ}}_{i, tail}^{2} I_{2})$ of the two-component PSF mixture. The corresponding contribution to the full data likelihood (4.1) of photon $i$ given that $Z_{i} = j$ is

{\{ϕ_{2} (x_{i}, y_{i} ∣ μ_{j}, \frac{{\tilde{σ}}_{i, core}^{2}}{w_{i}} I_{2}) g (w_{i} ∣ \frac{ν_{i, core}}{2}, \frac{ν_{i, core}}{2})\}}^{1 (S_{i} = 1)} \times

{\{ϕ_{2} (x_{i}, y_{i} ∣ μ_{j}, \frac{{\tilde{σ}}_{i, tail}^{2}}{w_{i}} I_{2}) g (w_{i} ∣ \frac{ν_{i, tail}}{2}, \frac{ν_{i, tail}}{2})\}}^{1 - 1 (S_{i} = 1)},

where $ϕ_{2} (\cdot ∣ \cdot)$ and $g (\cdot ∣ \cdot)$ denote, respectively, the densities of a bivariate Gaussian and of a gamma random variable, and $I_{2}$ is again the identity matrix of order 2. Thus,

x_{i}, y_{i} ∣ S_{i} = 1, Z_{i} = j, w_{i} \sim N (μ_{j}, \frac{{\tilde{σ}}_{i, core}^{2}}{w_{i}} I_{2}),

w_{i} ∣ S_{i} = 1 \sim Ga (ν_{i, core} / 2, ν_{i, core} / 2)

and

x_{i}, y_{i} ∣ S_{i} = 0, Z_{i} = j, w_{i} \sim N (μ_{j}, \frac{{\tilde{σ}}_{i, tail}^{2}}{w_{i}} I_{2}),

w_{i} ∣ S_{i} = 0 \sim Ga (ν_{i, tail} / 2, ν_{i, tail} / 2) .

A direct Gibbs sampler for $μ_{j}$ is achieved in three steps:

for every photon $i$ assigned to the $j$ th source draw from

S_{i} ∣ w_{i}, μ_{j}, x_{i}, y_{i}, \sim Be ({\tilde{p}}_{i}), with

{\tilde{p}}_{i} \propto ϕ_{2} (x_{i}, y_{i} ∣ μ_{j}, \frac{{\tilde{σ}}_{i, core}^{2}}{w_{i}} I_{2}) g (w_{i} ∣ \frac{ν_{i, core}}{2}, \frac{ν_{i, core}}{2}) p_{i};

for every photon $i$ assigned to the $j$ th source draw from

w_{i} ∣ S_{i} = 1, μ_{j}, x_{i}, y_{i} \sim Ga (\frac{ν_{i, core}}{2} + 1, \frac{1}{2} \{ν_{i, core} + \frac{| | r_{i}^{⊤} - μ_{j}^{⊤} | |^{2}}{{\tilde{σ}}_{i, core}^{2}}\})

w_{i} ∣ S_{i} = 0, μ_{j}, x_{i}, y_{i} \sim Ga (\frac{ν_{i, tail}}{2} + 1, \frac{1}{2} \{ν_{i, tail} + \frac{| | r_{i}^{⊤} - μ_{j}^{⊤} | |^{2}}{{\tilde{σ}}_{i, tail}^{2}}\}),

depending on the value of $S_{i}$ at the previous step. Here, $| | r_{i}^{⊤} - μ_{j}^{⊤} | |$ represents the squared Euclidean distance between the observed direction $r_{i} = (x_{i}, y_{i})$ and the direction $μ_{j}$ of source $j$ ;

write $x_{j} = (x_{1}, \dots, x_{n_{j}})$ , $y_{j} = (y_{1}, \dots, y_{n_{j}})$ , $s_{j} = (s_{1}, \dots, s_{n_{j}})$ and $w_{j} = (w_{1}, \dots, w_{n_{j}})$ . Draw from

μ_{j}^{⊤} ∣ s_{j}, w_{j}, x_{j}, y_{j} \sim N_{2} (ψ, ξ),

where

ψ = {(\sum_{i = 1}^{n_{j}} \frac{w_{i}}{{\tilde{σ}}_{i, s_{i}}^{2}})}^{- 1} (\sum_{i = 1}^{n_{j}} \frac{w_{i}}{{\tilde{σ}}_{i, s_{i}}^{2}} r_{i}^{⊤}), ξ = {(\sum_{i = 1}^{n_{j}} \frac{w_{i}}{{\tilde{σ}}_{i, s_{i}}^{2}})}^{- 1}

with ${\tilde{σ}}_{i, S_{i} = 1}^{2} = {\tilde{σ}}_{i, core}^{2}$ and ${\tilde{σ}}_{i, S_{i} = 0}^{2} = {\tilde{σ}}_{i, tail}^{2}$ , guarantees that the values are bounded into the map limits.

Jointly update the transformed background parameters ${log σ_{x}, log σ_{y}, logit (δ)}$ with a Metropolis–Hastings step. The proposed values are sampled from a three-dimensional Student t distribution with scale matrix $c_{prop} \cdot Σ_{prop}$ and degrees of freedom $ν_{prop}$ . Set $Σ_{prop}$ equal to minus the inverse observed Fisher information matrix obtained by fitting the background model (3.3) to the $γ$ -ray photon counts simulated using the Fermi LAT Science Tools (see Appendix). Set the tuning parameters $c_{prop}$ and $ν_{prop}$ to 0.6 and 3 to guarantee an acceptance rate of about 60% and a moderate autocorrelation among the draws.

4.2.2 Update $K$ for given model parameters (Step 2)

We update the size of the finite mixture model (3.1) by drawing from the posterior probability $π (K ∣ θ_{K}, x, y, z)$ of $K$ . This step changes the dimension of the parameter space and is carried out using the reversible jump MCMC algorithm described in Green (1995). According to Richardson and Green (1997) and Jones et al. (2015), we set four different types of possible moves called split, birth, combine and death. The first two increase the dimension of the parameter space by adding a new source to the model; the latter two reduce the size $K$ of the finite mixture by one. The four moves are chosen at random with equal probability. Note that adding a new source to the model or removing a source does not change the background component $g_{b} (x_{i}, y_{i} ∣ σ_{x}, σ_{y}, δ)$ .

In a split move, the algorithm attempts to augment the parameter space passing from a mixture of size $K$ to size $K^{'} = K + 1$ . Given a source $j$ , chosen at random, with location $μ_{j}$ and relative weight $ω_{j}$ , a split implies the original parameters to be separated into two components $(μ_{j_{1}}, μ_{j_{2}})$ and $(ω_{j_{1}}, ω_{j_{2}})$ . This is achieved by first generating three independent draws $u = (u_{1}, u_{2}, u_{3})$ from a Beta(2,2) distribution, as outlined on page 739 of Richardson and Green (1997), and then transforming $(μ_{j}, ω_{j}, u)$ as given there. The new entire set of parameters is called $θ_{K^{'}}$ . The split move is accepted with probability

α_{K}^{split} = min \{1, \frac{π (θ_{K^{'}}, K^{'} ∣ x, y, z)}{π (θ_{K}, K ∣ x, y, z)} \frac{q (K^{'} \to K) J (θ_{K}, u)}{q (K \to K^{'}) b (u)}\} .

Here, $q (K \to K^{'}) = 1 / (4 K)$ is the probability of augmenting the parameter space with the selected split move, while $q (K^{'} \to K) = 1 / {4 (K + 1)}$ is the probability of the reverse move. The expression $b (u)$ represents the product of three Beta(2,2) densities evaluated at $u$ and $J (θ_{K}, u) = ω_{1} / {u_{1} (1 - u_{1})}$ is the Jacobian of the transformation entailed by the split move.

A combine move reduces the size of the mixture from $K$ to $K^{'} = K - 1$ . The corresponding acceptance probability is just the inverse of $α_{K - 1}^{split}$ when we augment the size of the mixture from $K - 1$ to $K$ . That is, $α_{K}^{comb} = (α_{K - 1}^{split})^{- 1}$ .

The birth and death moves are simpler. For a birth move we generate a random direction $μ^{*}$ from the prior density $π (μ_{j})$ with corresponding intensity $ω^{*} \sim Beta (1, K)$ . Note that the existing intensities need be rescaled by $ω_{j} = ω_{j} (1 - ω^{*})$ , for $j = 0, \dots, K$ so that all weights of the augmented model sum to 1. In a death move, an existing component is randomly chosen and deleted. Again, the remaining weights are rescaled so as to sum up to 1. Note that all existent/remaining directions are maintained in a birth/death move.

5 A real data example

We fitted the model presented in Section 3 to the high-energy photons recorded by the Fermi LAT in the sky region framed in white in Figure 1. As already mentioned, this region is of size $[- 10^{\circ}, - 80^{\circ}] \times [10^{\circ}, 80^{\circ}]$ and contains $26 208$ events whose energy contents range from 10 GeV to 1 000 GeV. The 3FHL catalogue lists $132$ confirmed sources for this same region. Section 5.1 presents the parameter estimates, while in Section 5.2 we validate and benchmark our model by cross-checking the identified sources with the ones listed in the 3FHL catalogue.

5.1 Parameter estimates

Figure 5 plots the posterior distribution of the number $K$ of identified sources obtained from $20 000$ runs of the algorithm described in Section 4. All components of the chain stabilized very quickly. The posterior mode agrees with the upper bound $κ_{\max} = 200$ of the support of $K$ and was visited $1 590$ times after a burn-in of suitable length was discarded. Here, we will illustrate the results obtained by focusing on the minor mode $K = 188$ which was visited $1 290$ times. Indeed, because of as yet uncleared reasons, our algorithm tends to explore the upper bound of the support of $K$ a larger number of times as expected. A careful sensitivity analysis on the number of sources $K$ (results not reported here) revealed, however, that the directions and relative intensities of the detected sources are largely unaffected by the final value chosen for $K$ . That is, for $K = 200$ we re-detect the sources of $K = 188$ , with almost no change in their parameter estimates, plus additional 12 ones.

Figure 5:

Posterior distribution of the number of mixture components $K$ for the analysed sky region. The modal value $K = 200$ appears $1 590$ times. The minor mode $K = 188$ was visited 1 290 times

We hence rerun the algorithm for $10 000$ iterations with the number of components fixed at $K = 188$ . Figure 6 plots the traces of the $10 000$ simulated values for the background parameters $σ_{x}$ , $σ_{y}$ and $δ$ obtained for the model value K=188. The posterior modes and corresponding 95% HPD intervals are: ${\tilde{σ}}_{x} = 0.0218$ $[0.0215 - 0.0222]$ , ${\tilde{σ}}_{y} = 0.0282$ $[0.0279 - 0.0288]$ and $\tilde{δ} = 0.7372$ $[0.6871 - 0.7895]$ . The marginal distribution of the longitude for $γ$ -ray photons generated by the diffuse background is more dispersed than for the latitude as we may have expected from Figure 3. The Bayesian estimate of the correlation parameter $δ$ is quite high. There is a rather strong linear dependence between the two exponential distributions which marginally model the background contamination along its longitude and its latitude.

Figure 6:

Trace plots of the $10 000$ simulated values of the background model parameters obtained for the modal value $K = 188$ : $σ_{x}$ (top), $σ_{y}$ (middle) and $δ$ (bottom). The solid horizontal line at the centre represents the posterior mode; the two dashed lines delimit the $95 %$ highest posterior density (HPD) interval

Table 2 lists the first 10 sources detected by our model with the corresponding estimated longitudes ( $l_{alg}$ ) and latitudes ( $b_{alg}$ ) in Columns 2 and 3. Column 4 reports the estimated mixing proportions $ω$ . All detected sources contribute a total of 6.27% of the detected $γ$ -ray photons with an average intensity of 0.000293, while the intensity of the diffuse background is $ω_{0} = 0.9373$ (HPD interval $[0.9340 - 0.9399]$ ).

Table 2:

Sample list of the sources detected by our model in the analysed sky region. Columns 2 and 4 report the estimated longitudes ( $l_{alg}$ ) and latitudes ( $b_{alg}$ ); Column 4 records the estimated mixing proportions ( $ω$ ). These sources were associated with confirmed sources present in the 3FHL catalogue within an angular distance of $1^{\circ}$ (Column 8). The longitudes ( $l_{3 FHL}$ ) and latitudes ( $b_{3 FHL}$ ) listed in the catalogue are given in Columns 6 and 7. The complete list can be found in the Supplementary Material

Detected source				Confirmed source
${ID}_{alg}$	$l_{alg}$	$b_{alg}$	$ω$	${ID}_{3 FHL}$	$l_{3 FHL}$	$b_{3 FHL}$	distance( $^{\circ}$ )
1	$-$ 79.213	17.062	1.27e-04	1	$-$ 79.214	17.071	0.010
2	$-$ 79.431	21.486	1.76e-04	2	$-$ 79.424	21.503	0.019
3	$-$ 79.339	28.077	1.65e-04	3	$-$ 79.333	28.090	0.014
4	$-$ 75.963	23.756	2.08e-04	4	$-$ 75.954	23.754	0.009
5	$-$ 74.738	21.444	1.17e-04	5	$-$ 74.759	21.463	0.028
6	$-$ 76.516	28.177	3.01e-04	6	$-$ 76.512	28.177	0.003
7	$-$ 70.647	13.130	4.38e-04	7	$-$ 70.624	13.160	0.037
8	$-$ 70.732	22.939	3.21e-04	8	$-$ 70.744	22.951	0.016
9	$-$ 78.489	46.552	2.43e-04	9	$-$ 78.463	46.580	0.034
10	$-$ 73.005	35.191	2.40e-04	10	$-$ 72.991	35.182	0.014

5.2 Model validation

Our model detected a total of 188 high-energy emitting sources. We were able to match 116 of these with one of the 132 confirmed sources listed in the 3FHL catalogue within and angular distance of $1^{\circ}$ . An excerpt is given in Table 2. Columns 6 and 7 of the table report the longitude ( $l_{3 FHL}$ ) and latitude ( $b_{3 FHL}$ ) of the matched confirmed source. Column 8 reports the angular distance between the two sources, the detected and the confirmed one. The complete list of all possible matches is given in the Supplementary Material. This list consists of 127 rows, 11 of which are struck through. These are for the nine detected sources with two possible matches, and the two catalogue sources associated with two detected sources. For these multiple matches, the finally selected matched pair corresponds to the two sources whose angular distance is smaller. All observed angular distances are below $0 . 5^{\circ}$ and most are within the typical values of $0.2 - 0 . 3^{\circ}$ used in Astrophysics.

In all, we were able to identify 116 confirmed sources out of 132 which corresponds to a success rate of 88%; 16 sources remained unidentified. There are 72 detected sources with no correspondence in the 3FHL catalogue. We used the a posteriori available information on their intensities to further discriminate whether these correspond to real $γ$ -ray emitting sources. A threshold value was defined as the median of the estimated mixing proportions $ω_{j}$ for the 116 identified sources. Out of the 72 detected but unmatched sources, 20 sources exhibited posterior modes for $ω_{j}$ which exceeded this threshold value. We qualified them as possible unidentified sources. They are currently under examination as foreseen by the Fermi LAT Collaboration. The 52 discarded sources may be false positives or they may correspond to as yet unidentified sources which are characterized by too faint signals to be detected with sufficiently high power.

6 Discussion

The model we developed in this paper represents a promising resource for identifying extra-galactic high-energy emitting sources using data from the Fermi LAT Collaboration. Indeed, we were able to identify 116 of the 132 $γ$ -ray emitting sources confirmed for the analysed sky region, which corresponds to a detection rate of almost 90%. Of the 72 sources detected by our algorithm, which we were not able to match with a confirmed source, 20 are currently under investigation. However, the same model can easily be adapted to other types of highly contaminated data in the two-dimensional space where the set of locations at which a spatial process is observed are of interest.

Nonetheless, there are still a number of issues with the model and with the algorithm that need be settled.

A first question is why the posterior distribution of $K$ tends to concentrate on the upper bound of its support $κ_{\max}$ . To explore the implications we carried out a sensitivity analysis by inspecting the posterior distributions of the remaining model parameters for a set of values of $K$ between the minor mode $K = 188$ and $K = 200$ . The reconstructed directions and relative intensities of the detected sources remained stable for all these scenarios. The high number of detected sources which found no confirmed counterpart may be a side effect of the parametric background model which is not flexible enough to capture local irregularities and thus generates a large amount of false positives. Some of them, however, may correspond to sources whose signal is too faint to be detected at a sufficiently high significance. Beware that the commonly used two-step procedure which led to the 3FHL catalogue described in the introduction classifies a new discovery using the 5 sigma standard, that is, for $p$ -values smaller than $2.87 \times 10^{- 7}$ (Dorigo, 2016). The heuristic approach used to qualify the newly detected sources based on the estimated intensities needs to be replaced by a formal procedure which accounts also for the further available information on the photons such as the energy content. The ultimate solution would be to directly model the Fermi LAT intensity maps, especially if we wanted to work on a sky portion of the Galactic plane where it is more complicated to disentangle sources from the background.

A further issue is how to incorporate the possible fish eye effect into the current PSF. Indeed, as we used a mixture of two bivariate Student t distributions with equal spread on both directions (longitude and latitude), the scattering of the $γ$ -ray photons around their emitting source is necessarily symmetric. Asymmetric model formulations are available, but the corresponding densities are no longer available in closed form. Future developments must also consider the fine tuning of our algorithm so as to be able to analyse larger portions of the sky.

This article includes a substantial amount of improvements of the model and of the fitting routine with respect to a very first and highly incomplete attempt which was presented at an Italian national conference (Sottosanti et al., 2019). In particular, we replaced the simplistic approximation of the instrument response function based on King's PSF (King, 1962) with the more PSF function proposed in Ackermann et al. (2013), as discussed in Section 3.1. This allowed us to take into account both the energy content of the photon and its event type, which determine the quality of each event. We furthermore carefully studied the GIE component of the $γ$ -ray emitting background. This allowed us to formulate the simple though theoretically coherent parametric model of Section 3.2 for the diffuse background component which incorporates all information provided on it by the Fermi LAT project. On the computational side, we heavily improved the efficiency of the reversible jump MCMC algorithm of Sottosanti et al. (2019) both in terms of methodological innovations and coding; see Section 4. Last but not least, we carried out a sensitivity analysis aimed at assessing the extent to which the parameter estimates of the single source and background models may depend on the number $K$ of detected sources.

APPENDIX

Despite its simple form, the background model (3.3) wraps up the insight gained from a detailed study of the available information on the diffuse $γ$ -ray background. Indeed, the background emission itself has two origins: the isotropic $γ$ -ray background which uniformly spreads over the entire map and the so-called Galactic Interstellar Emission (GIE) which dominates the galactic part of the sky map.

The Fermi LAT Collaboration provides data on the GIE component in the form of $30$ so-called intensity maps measured in ${cm}^{- 2} s^{- 1} {sr}^{- 1} {MeV}^{- 1}$ . These maps are based on an elaborate theoretical model which incorporates the current knowledge of the physical phenomena which determine GIE (Acero et al., 2016). Informally, we can interpret these intensities as expected counts which were normalized with respect to four features of the space telescope. These are the so-called effective area which broadly speaking corresponds to the size (in ${cm}^{2}$ ) of the surface of the LAT, the duration of observation (in seconds, $s$ ), the volume of the sky (in steradians, $sr$ ) covered by the space telescope and the energy contents of the $γ$ -ray photon (in MeV). The $30$ intensity maps correspond to 30 energy bins which range from $58$ to $513 056$ MeV with a step size of $0.136$ measured on the log $_{10}$ scale. Each map divides the whole sky into a spatial grid of $2 880 \times 1 441$ pixels at which the intensities are provided (Gorski et al., 2005).

In Section 3.2, we simulated a count map of the GIE using the intensities of these 30 GIE maps thanks to the support provided by the Fermi LAT Science Tools (https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/). Here, we will present the different steps of an empirical analysis which allowed us to drop the dependence on the features of the space telescope when modelling the diffuse background.

We first converted the $30$ intensity maps into a unique whole sky count map, that is, into the same form our primary dataset is delivered. We can easily account for the energy layer by multiplying each intensity by the central value of the energy bin it refers to. These intermediate results are then multiplied by $4 π / (2.800 \times 1.441)$ , which is the solid angle covered by the Fermi LAT at every pixel, and by the total observation time of 7.2 years (expressed in seconds). An issue remains with the effective area component. This area depends on the geometrical conformation of the space telescope as well as on the capacity of the detector to convert and correctly identify incoming $γ$ -rays. The behaviour of the effective area is described on the Fermi LAT performance page (http://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm) as a function of the energy content of the photon and of its direction expressed in polar coordinates. The four panels of Figure 7 plot the effective area (in $m^{2}$ ) which characterizes the latest release of the Fermi LAT software (Pass8) for the four event types PSF0, PSF1, PSF2 and PSF3. Recall that these correspond to increasing quality of the measurements. Independently of the event type, the higher the energy content of the photon is and the steeper it hits the surface of the LAT, the larger is the effective area. It seems, however, that with the exception of the less accurate event type PSF0 and for incident angles which are steep enough we may approximate the effective area by a constant value for energy levels larger than 10 GeV (that is, 4 on the log $_{10}$ scale), which corresponds to the energy range we are interested in. This motivated us to discard in our analysis all events of type PSF0.

Figure 7:

The effective area (in $m^{2}$ ) for the Pass8 release of the Fermi LAT software. The quality of the measurements, as expressed by the four event types PSF0, PSF1, PSF2 and PSF3, increases from the top left to the bottom right panel. The x-axis reports the energy content of the photon (in units of MeVs) on the log $_{10}$ scale; the y-axis reports the cosine of the incident angle ( $θ$ )

A rough approximation of the expected number of counts is hence provided by summing up the rescaled intensities over the 30 maps on a pixel by pixel basis. The bottom panel of Figure 4 reports the approximate count map we reconstructed according to the steps outlined in the previous paragraphs. This map is in close agreement with the simulated count map of Section 3.2 and shown in the top panel of Figure 4. On this basis, we decided to neglect the influence of the effective area of the space telescope in the formulation of our background model (3.3).

As far as the prior elicitation on the parameters of Model (3.3) goes, an alternative could have been to fit the bivariate exponential model (3.3) by maximum likelihood to the reconstructed count maps and to centre the prior distributions (3.4) at the corresponding maximum likelihood estimates. We deliberately did not choose this solution because of two reasons. First, Model (3.3) embraces both diffuse background components, that is, the GIE and the isotropic one, while the reconstructed count maps are for the GIE. Second, the reconstructed count map depends heavily on how we rescaled it with respect to the effective area. In any case, the fraction of γ -ray photons which are emitted by the background is so high that their likelihood contribution overweights any prior elicitation.

Supplementary material

Supplementary materials for this article are available from http://www.statmod.org/smij/archive.html.

Footnotes

Acknowledgements

We would like to thank the associate editor and an anonymous referee for their most useful comments on a previous version of the article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

This research was supported by SID 2018 grant ‘Advanced statistical modelling for indexing celestial objects’ (BIRD185983) awarded by the Department of Statistical Sciences of the University of Padova. We furthermore acknowledge the financial support by Professor Junhui Fan (grants n. NSFC11733ss001 and n. NSFCU1531245).

References

Acero

Ackermann

Ajello

Albert

Baldini

Ballet

Barbiellini

Bastieri

Bellazzini

Bissaldi

. (2016) Development of the model of galactic interstellar emission for standard point-source analysis of Fermi Large Area Telescope data. The Astrophysical Journal Supplement Series, 223, 26.

Ackermann

Ajello

Allafort

Asano

Atwood

Baldini

Ballet

Barbiellini

Bastieri

Bechtol

. (2013) Determination of the point-spread function for the Fermi Large Area Telescope from on-orbit data and limits on pair halos of active galactic nuclei. The Astrophysical Journal, 765, 54.

Ajello

Atwood

Baldini

Ballet

Barbiellini

Bastieri

Bellazzini

Bissaldi

Blandford

Bloom

. (2017) 3FHL: The third catalog of hard Fermi-LAT sources. The Astrophysical Journal Supplement Series, 232, 18.

Atwood

Abdo

Ackermann

Althouse

Anderson

Axelsson

Baldini

Ballet

Band

Barbiellini

. (2009) The large area telescope on the Fermi gamma-ray space telescope mission. The Astrophysical Journal, 697, 1071.

Ciprini

Tosti

Marcucci

Cecchi

Discepoli

Bonamente

Germani

Impiombato

Lubrano

Pepe

(2007) The first GLAST symposium. In American Institute of Physics Conference Series, edited Ritz

Michelson

Meegan

C. A.

. Vol. 921, pages 546–547. College Park, MD: American Institute of Physics.

Damiani

Maggio

Micela

Sciortino

(1997) A method based on wavelet transforms for source detection in photon-counting detector images. I. Theory and general properties. The Astrophysical Journal, 483, 350.

Dorigo

(2016) Anomaly! Collider Physics and the Quest for New Phenomena at Fermilab. Singapore: World Scientific.

Gaetan

Guyon

(2010) Spatial Statistics and Modeling. Berlin: Springer-Verlag.

Gelman

Carlin

Stern

Dunson

Vehtari

Rubin

(2013) Bayesian Data Analysis, 3rd Edition. Boca Raton, FL: Chapman & Hall/CRC.

10.

Gorski

Hivon

Banday

Wandelt

Hansen

Reinecke

Bartelmann

(2005) HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere. The Astrophysical Journal, 622, 759.

11.

Green

(1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82, 711–732.

12.

Gumbel

(1960) Bivariate exponential distributions. Journal of the American Statistical Association, 55, 698–707.

13.

Jones

Kashyap

Van Dyk

(2015) Disentangling overlapping astronomical sources using spatial and spectral information. The Astrophysical Journal, 808, 137.

14.

King

(1962) The structure of star clusters. I. An empirical density law. The Astronomical Journal, 67, 471.

15.

Mattox

Bertsch

Chiang

Dingus

Digel

Esposito

Fierro

Hartman

Hunter

Kanbach

. (1996) The likelihood analysis of EGRET data. The Astrophysical Journal, 461, 396.

16.

Richardson

Green

(1997) On Bayesian analysis of mixtures with an unknown number of components (with discussion). Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59, 731–92.

17.

Sottosanti

Costantin

Bastieri

Brazzale

(2019) Discovering and locating high-energy extra-galactic sources by Bayesian mixture modelling. In New Statistical Developments in Data Science, edited by Petrucci

Racioppi

Verde

, pages 135–148. New York, NY: Springer International Publishing.

18.

Starck

J-L

Pierre

(1998) Structure detection in low intensity x-ray images. Astronomy and Astrophysics Supplement Series, 128, 397–407.

Bayesian mixture modelling of the high-energy photon counts collected by the Fermi Large Area Telescope

Abstract

Keywords

1 Introduction

2 The Fermi LAT data

Figure 1:

Whole sky map at γ -ray wavelengths and energies larger than 1 GeV based on data accumulated by the LAT over a period of five years of operation (Image Credit: NASA/DOE/Fermi LAT Collaboration). The region framed in white represents the area analysed in this article

Figure 3:

Top: Contour plot of the simulated background photon density for the second sky quadrant spanned by [ 0 ∘ , − 180 ∘ ] × [ 0 ∘ , 90 ∘ ] . Bottom: Image plot of the approximate count map reconstructed using the 30 GIE intensity maps provided by the LAT

4.1 The posterior distribution

4.2.1 Estimate model parameters for given K (Step 1)

4.2.2 Update K for given model parameters (Step 2)

5 A real data example

5.1 Parameter estimates

Figure 5:

Posterior distribution of the number of mixture components K for the analysed sky region. The modal value K = 200 appears 1 590 times. The minor mode K = 188 was visited 1 290 times

6 Discussion

APPENDIX

Supplementary material

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

References

Top: Contour plot of the simulated background photon density for the second sky quadrant spanned by $[0^{\circ}, - 180^{\circ}] \times [0^{\circ}, 90^{\circ}]$ . Bottom: Image plot of the approximate count map reconstructed using the 30 GIE intensity maps provided by the LAT

4.2.1 Estimate model parameters for given $K$ (Step 1)

4.2.2 Update $K$ for given model parameters (Step 2)

Posterior distribution of the number of mixture components $K$ for the analysed sky region. The modal value $K = 200$ appears $1 590$ times. The minor mode $K = 188$ was visited 1 290 times