Two Dimensional Clustering of GammaRay Bursts using durations and hardness
Abstract
GammaRay Bursts (GRBs) have been conventionally bifurcated into two distinct categories: “short” and “long” with durations less than and greater than two seconds respectively. However, there is a lot of literature (although with conflicting results) regarding the existence of a third intermediate class. To investigate this issue, we extend a recent study by Kulkarni & Desai 2017a on classification of GRBs to two dimensions by incorporating the GRB hardness in addition to the observed durations. We carry out this unified analysis on GRB datasets from four detectors, viz. BATSE, RHESSI, Swift (observed and intrinsic frame), and FermiGBM. We consider the duration and hardness features in logscale for each of these datasets and determine the bestfit parameters using Gaussian Mixture Model. This is followed by information theoretic criterion (AIC and BIC) to determine if a threecomponent fit is favored compared to a twocomponent one or viceversa. For BATSE, we find that both AIC and BIC show preference for three components with decisive significance. For Fermi and RHESSI, both AIC and BIC show preference for two components, although the significance is marginal from AIC, but decisive using BIC. For Swift dataset in both the observed and rest frame, we find that three components are favored according to AIC with decisive significance, and two are preferred with BIC with marginal to strong significance.
Department of Computer Science and Engineering, NIT Raipur, Chhatisgarh492010, India \@footnotetextDepartment of Physics, University of Florida, Gainsville, FL 32611, USA \@footnotetextDepartment of Physics, IIT Hyderabad, Kandi, Telangana502285, India \@footnotetextDepartment of Computer Science and Engineering, IIT Hyderabad, Kandi, Telangana502285, India
Keywords GRB Classification; Bayesian Information Criterion; Akaike Information Criterion
1 Introduction
Gammaray bursts (GRBs) are shortduration energetic cosmic explosions with prompt emission between keVGeV energies, which are been continously detected at the rate of about one per day (Zhang et al., 2016a; Schady, 2017). The first convincing case for bifurcating the GRB population into two categories was made from an analysis of the BATSE data (Kouveliotou et al., 1993), and led to establishing the conventional classification of GRBs into short (T90 2 s) and long (T90 2 s) classes, where T90 is the time which encompasses 90% of the burst’s fluence, and is usually used as a proxy for the duration of a GRB.
Despite the conventional wisdom of only two distinct GRB classes, multiple groups have argued over the years for the existence of an intermediate class of GRBs in between the short and long bursts, using T90 as the criterion for classification. The first such claim for an intermediateduration GRB class, with T90 in the range 2â10s in the BATSE dataset was put forward by Horváth (1998) and Mukherjee et al. (1998) and subsequently confirmed by the analysis of the complete BATSE dataset (Horváth, 2002; Chattopadhyay et al., 2007a). However, recently this was disputed by Zitouni et al. (2015), who found that two distributions fit the BATSE T90 data much better compared to three components. Evidence for a third lognormal component was also found in Swift/BAT data (Horváth et al., 2008; Zhang & Choi, 2008; Huja et al., 2009; Horváth et al., 2010; Horváth & Tóth, 2016; Zitouni et al., 2015; Tarnopolski, 2016b). However, these results have been disputed by other authors, who found that the T90 distribution prefers two component (Zhang et al., 2016b). Most recently, Kulkarni & Desai (2017a) carried out a unified classification of the T90 distributions for the GRB datasets from BATSE, Fermi, Swift, and BeppoSax and found that among these, only for Swift GRBs in the observed frame is the evidence for three classes marginally significant at about . However, when the same analysis is done for the Swift GRBs in the intrinsic GRB frame, two components are preferred. For all other datasets, evidence for three components is either very marginal or disfavored.
Extension of studies on GRB classification using both duration and hardness (defined as the ratio of fluences between two different energy bands) have also not reached a common consensus. Horváth et al. (2006) and Chattopadhyay et al. (2007a) argued for three components in the BATSE GRB data using twodimensional clustering in T90hardness and T90fluence planes respectively. Most recently, Chattopadhyay & Maitra (2017) have argued for more than three components in the BATSE data by clustering in six dimensions. For Swift data, Veres et al. (2010) showed using multiple clustering techniques that three components are favored in the two dimensional (T90)  (hardness) plane and the intermediate class has overlap with Xray flashes. However, these results are in conflict with more recent analysis by Yang et al. (2016), who showed by applying two dimensional GMM models on T90 and hardness ratio on Swift GRBs, that the data favor only two components instead of three.
To resolve this imbroglio, we use twodimensional clustering in the hardness vs T90 plane to find out the optimum number of GRB classes, similar to studies done in Yang et al. (2016). We then uniformly apply this method to the latest available data from all the GRB detectors, for which both T90 and hardness have been provided (or can be inferred from the catalogs). These detectors include BATSE, FermiGBM, RHESSI, and Swift. For model comparison, we use two widely used information theoretic criterion, viz. Akaike Information Criterion and Bayesian Information Criterion. We note that in Yang et al. (2016), only Bayesian Information criterion has been used to evaluate the optimum number of components.
Both of these informationcriterion based model comparison techniques have been applied to a variety of problems in astrophysics and particle physics, including in the classification of GRBs (Shi et al., 2012; Desai & Liu, 2016; Desai, 2016; Kulkarni & Desai, 2017a; Ganguly & Desai, 2017; Kulkarni & Desai, 2017b) and references therein.
The outline of this paper is as follows. In Sect. 2 we discuss the methodology used to obtain the bestfit parameters for the mean GRB duration and hardness, along with their covariances after positing two and three classes of GRBs. In Sect. 3, we discuss various techniques used for model comparison. We then present our results for the various GRB datasets in Sect. 4, including a very brief comparison with previous results. We conclude in Sect. 5.
2 Parameter Estimation
2.1 Datasets
Herein, we consider the GRB datasets available from BATSE^{1}^{1}1http://gammaray.msfc.nasa.gov/batse/grb/catalog/current 4B catalog (Paciesas et al., 1999), Swift^{2}^{2}2http://swift.gsfc.nasa.gov/archive/grbtable (Lien et al., 2016), FermiGBM^{3}^{3}3http://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst (Narayana Bhat et al., 2016) and RHESSI^{4}^{4}4https://heasarc.gsfc.nasa.gov/W3Browse/all/rhessigrb.html (Řípa et al., 2009a). For these detectors, spectral information is available. The number of GRBs analyzed for the model comparison are 1939 from BATSE, 991 from Swift, 1376 from Fermi, and 427 from RHESSI. Among these, Swift and Fermi detectors are still online and detecting on the order of about one new GRB per day. We did not consider other catalogs such as those from BeppoSax, INTEGRAL etc, as they either contained a very small sample of GRBs or didn’t have any publicly available data on hardness of the observed bursts.
2.2 Hardness Calculation
Spectral hardness (or hardness ratio) of GRBs is defined as the ratio between the GRB fluence in different energy bands. For most of the datasets analyzed , we use as the ratio between the keV and the keV bands.
Usually, the hardness can be trivially calculated from the ratio of the fluences provided in the catalogs. However, for the most recent Swift catalog (Lien et al., 2016), the fluence was not provided in the public catalogs. So we calculate the hardness from the spectral fits to the Swift GRBs, for which the coefficients were made publicly available. The results from the power law fits are provided in the Swift spectral catalog. The hardness ratio in the observed frame is calculated as,
(1) 
where represents the photon flux at a given energy. Two functional forms for the powerlaw fits have been posited for the Swift catalog: a simple power law as well as the cutoff power law,
(2)  
(3) 
where is the power law photon index and is the peak energy. In Yang et al. (2016), both of the above power laws have been used in the hardness calculation. Here, we shall calculate the hardness using only a simple power law.
2.3 Fitting method
We follow the same procedure as in Zhang et al. (2009) to find the optimum number of components using twodimensional clustering. We use the Gaussian mixture models applied to log (T90) and log (hardness) by varying the number of components and finding the bestfit parameters for each component using the EM algorithm (Dempster et al., 1977). For a given probability density function , where denotes the observed datapoints, represents the parameters used to define the density function, being the total number of GRBs in our study and the probabilities associated with each of the Gaussian distributions (mixing proportions), the loglikelihood can be defined as:
(4) 
and the probability distribution function for a univariate Gaussian as,
(5) 
This can generalized to a bivariate distribution as,
(6) 
where is the correlation, is the covariance of the two variables and is the mean and is the mean .
3 Model Comparison
The comparison of models on the basis of the bestfit likelihood is not the optimum way to do hypothesis testing or to select the preferred model after finding the bestfit parameters for each model. Even though the value of the likelihood increases, addition of extra free parameters leads to overfitting. Therefore, the additional free parameters need to be penalized so as to avoid getting a bad result. To address these issues, a number of both frequentist and Bayesian modelcomparison techniques have been used over the past decade to determine the best model which fits the observational data (Liddle, 2004, 2007; Lyons, 2016). Here, we use information criterion based tests such as Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) for model comparison, since these are straightforward to compute from the likelihoods. AIC/BIC have also been previously used for GRB classification by a number of authors (Mukherjee et al., 1998; Tarnopolski, 2016a, b; Yang et al., 2016; Zhang et al., 2016b). More information about AIC and BIC and its application to astrophysical problems can be found in Liddle (2004, 2007); Shi et al. (2012), etc.
3.0.1 Aic
The Akaike Information Criterion (AIC) is used for model comparison, when we need to penalize for any additional free parameters to avoid overfitting. A preferred model in this test is the one with the smaller value of AIC between the two hypothesis. The AIC is given by,
(7) 
where is the number of free parameters in the model and is the likelihood. The second term favors models with high value of likelihood, while the first term penalizes models which uses large number of parameters. Models with large number of parameters might have a high likelihood but will over fit on the data. The AIC defined in Eq. 7 is good when the ratio is very large i.e. 40 (Burnham & Anderson, 2004). For smaller values, a first order correction is sometimes used (Ganguly & Desai, 2017). However, as all our datasets have a ratio of greater than 40, we don’t have to correct for this. The absolute value of AIC is usually not of interest. The goodness of fit between two hypothesis (A) and (B) is described by the difference of the AIC values and is given by,
(8) 
where  correspond to the AIC values for the hypothesis A and B. Burnham & Anderson (2004) have provided qualitative strength of evidence rules to assess the significance of a model based on the AIC values between the two models. If AIC, then it is considered as strong evidence against the model with higher AIC and AIC is considered as decisive evidence against the model with higher AIC (Liddle, 2007). Values of AIC correspond to weak evidence. ^{5}^{5}5To avoid any ambiguity in our representation of our results, we have consistently kept the 3Gaussian model as the null hypothesis which simplifies the analysis and makes a positive value of , favor the 3Gaussian and a negative value favors the 2Gaussian.
3.0.2 Bic
The Bayesian Inference Criterion (BIC) is also used for penalizing the use of extra parameters. As in the case of AIC, the model with the smaller value of BIC is the preferred model. The penalty in the BIC test is harsher than that in the case of AIC and is given by,
(9) 
The logarithmic term and the number of free parameters act as a very harsh measure needed for the BIC test. The goodness of fit used for hypothesis testing between two models and is given by,
(10) 
Similar to AIC, the model with lower value of BIC is favored. To assess the significance of a model, strength of evidence rules have also been proposed based on BIC (Kass & Raftery, 1995), which are approximately the same as those for AIC.
4 Results
We apply all the techniques discussed in the previous sections to GRB datasets from multiple detectors. For data from each of the GRB detectors, we find the mean value of T90 and its standard deviation by varying the total number of components from one to five, followed by maximizing the likelihood in Eq. 4 for each of the hypothesis. For these bestfit parameters, we then implement the information criterion based modelcomparison techniques outlined in Sect. 3.
The GMM and the corresponding parameter estimation using the EM algorithm are implemented using the sklearn.mixture module of the python library Scikitlearn. Covariance types full and tied are used for generating the model with the number of components ranging from one to five and we choose the covariance type, which yields the maximum value of likelihood. For all the detectors except Fermi, the maximum value is achieved for covariance type equal to full. Note that our main goal is to try to ascertain whether a threecomponent fit is favored compared to a twocomponent one or viceversa. Therefore, we are agnostic to the value for the number of components for which we get the minimum value of AIC or BIC, although we do report its value for all the detectors.
We now present our results for the GRB datasets from BATSE, RHESSI, FermiGBM and Swift.
4.1 Batse
The current BATSE GRB (Paciesas et al., 1999) catalog contains 2036 GRBs detected between 1991 and 2000, of which we have used 1939 for our categorization purposes. We then apply the parameter estimation procedure outlined in Sect. 2.3. While fitting for two components, we find that 681 and 1355 GRBs belong to the short and long category respectively. When we fit for three components we find a total of 574, 886, and 499 GRBs in the short, intermediate, and long categories respectively. A complete summary of the results on applications of GMM, including the bestfit parameters and their covariance matrices are shown in Table 1. The AIC and BIC plots as a function of the number of components can be found in Fig. 1. Here, both AIC and BIC prefer three components. The AIC and BIC values in both cases cross the threshold of 10, needed for decisive evidence. The ellipses for the three components are shown in Fig. 2.
There is more than 20 years of literature on classification of BATSE GRBs. So, we only compare our results to a few selected papers, where both T90 and hardness (or other fluence related parameters) are used for classification. Results of classification of BATSE GRBs using only T90 are summarized in Kulkarni & Desai (2017a). The first cogent case for three GRB classes in BATSE data using spectral information, was made by Chattopadhyay et al. (2007b), who used two multivariate clustering methods using means partitioning and Dirichlet mixture modeling using fluence vs T90 to argue for three components. However, no estimate of the significance was made. Around the same time, an analysis similar to this using GMM in the log(T90) log(hardness) plane was done by Horváth et al. (2006), who find that three components were favored using frequentist model comparison by evaluating chisquare probability with the addition of the third component. They also found an anticorrelation between the duration and the hardness. Most recently, Chattopadhyay & Maitra (2017) used GMM based analysis on six different variables (two durations, peak flux, total fluence, and two spectral hardness ratios), five types of bursts are preferred. Our analysis finds evidence for three components with decisive significance.
AIC  BIC  

2  (0.032,2.08)  658  12709.3  12770.5  88.1  54.6  
(3.56,1.11)  1281  
3  (0.58,2.30)  574  12621.2  12715.9  
(2.81,0.92)  886  
(3.72,1.28)  499 
4.2 Rhessi
The RHESSI catalog (Řípa et al., 2009a), which we analyzed contains 1939 GRBs detected between 2002 and 2008. The hardness was obtained by calculating the ratio of the fluence in the (1201500 keV) range and (25120 keV) range. A tabular summary of the results from the GMM algorithm, including the bestfit parameters and their covariance matrices are shown in Table 2. For a twocomponent fit, we find that 107 and 320 GRBs belong to the short and long categories. For a threecomponent fit there are 164, 72, and 191 belonging to short, intermediate, and long categories respectively. The AIC and BIC plots as a function of the number of components can be found in Fig. 3. The ellipses for the three components are shown in Fig. 2. Note that the bestfit maximum likelihood we obtain is negative. Therefore, we obtain negative values of AIC and BIC for and a negative value of AIC for . The scatter plot of log(hardness) vs log (T90) along with the contours can be shown in Fig. 4. Here, both AIC and BIC prefer two components. The BIC value in both the cases crosses the threshold of 10, needed for decisive evidence. However the significance from AIC value is very marginal.
Previously, classification of RHESSI GRBs was done by using both T90 only and from a clustering analysis in the T90 vs hardness plane (Řípa et al., 2009b). The model comparison was done by comparing the difference in likelihoods. Řípa et al. (2009b) found that including both durations and hardness shows a preference for three components, with the probability that the third peak been only a fluctuation to be approximately 0.13% (Řípa et al., 2009b). Our analysis also prefers two components, although the significance is decisive only when evaluated using BIC.
AIC  BIC  

2  (0.34,0.06)  107  54.3  9.7  4.2  28.5  
(3.00,0.06)  320  
3  (0.45,0.06)  164  50.1  18.8  
(2.71,0.6)  72  
(3.11,0.06)  191 
4.3 Fermi
The most recent FermiGBM catalog (Narayana Bhat et al., 2016) contains 1376 GRBs with tabulated values of durations and hardness from 2008 to 2016. For FermiGBM, the hardness is defined as the ratio of fluence in 50300 keV to that in 1050 keV. A complete summary of the results on applications of GMM for FermiGBM data, including the bestfit parameters and their covariance matrices are shown in Table 3. As stated earlier, this is the only dataset for which the maximum value of the likelihood is obtained for covariance type set to tied. For two components, we find that 1144 and 232 GRBs belong to short and long categories respectively. For three components, 164, 72, and 191 belong to short, intermediate, and long categories respectively. The AIC and BIC plots as a function of the number of components can be found in Fig. 5. The minimum value of AIC is obtained for five components. Since our main goal is to find the optimum solution between two and three components, we do not study the properties of the 5component fit. However, it is possible that the number of GRBs in some of the categories is negligible for , making it similar to the solution. Therefore, in Table 3, we show the AIC results only for and . AIC prefers component over , although with marginal significance given by AIC=1.5. The BIC reaches a minimum at two components and BIC between and is greater than 10. Therefore, BIC prefers two components with decisive significance. The ellipses for the two components are shown in Fig. 6.
This is the first paper, which uses both hardness and T90 for the classification of Fermi GRBs. A summary of previous results on the classification of FermiGBM using durations can be found in Kulkarni & Desai (2017a). All previous classification studies with Fermi GRBs show a preference for two GRBs. Our analysis using twodimensional clustering also confirms the evidence for two components, although the significance is decisive only with BIC.
AIC  BIC  

2  (0.56,0.30)  1144  7820.6  7862.5  1.5  17.1  
(3.45,0.45)  232  
3  (0.08,0.58)  564  7822.1  7879.6  
(2.86,0.46)  590  
(3.66,0.43)  222 
4.4 Swift
The current Swift GRB (Lien et al., 2016) catalog (as of Jan. 2017) contains 1376 GRBs detected after 2004. Since the integrated fluence is not provided for the complete list of 2017 Swift GRBs (unlike the dataset used by Yang et al. (2016)), we calculate the fluence and hardness using the method in Sect. 2.2. We find that for a twocomponent fit, 273 and 718 belong to short and long categories respectively, whereas for a threecomponent fit 527, 342, and 122 belong to short, intermediate, and long categories respectively. A complete summary of the results on applications of GMM, including the bestfit parameters and their covariance matrices are shown in Table 4. The AIC and BIC plots as a function of the number of components can be found in Fig. 8. Here, AIC prefers two components with decisive evidence (AIC ) and BIC prefers two components with strong evidence (BIC). The ellipses for the three components are shown in Fig. 8.
Similar to BATSE, there is a lot of literature on GRB classification from the Swift catalogs, starting from 2008 or so (Horváth et al., 2008). These classification results from Swift based only on T90 are reviewed in Kulkarni & Desai (2017a) and no consensus has emerged among the different groups on the optimum number of GRB categories in the Swift dataset. The first comprehensive twodimensional classification of GRBs in the log (T90)  log (hardness) plane was done by Veres et al. (2010) using 408 Swift GRBs detected until 2009. Using BIC, they found support for three components with decisive evidence. A similar analysis done around the same time by Horváth et al. (2010), using 325 bursts confirmed these earlier results. However, based on a similar analysis of an updated Swift catalog upto Dec. 2012, containing 300 bursts with measured redshifts, Yang et al. (2016) found that the data prefer two components and BIC between two and three components is about 6.5, corresponding to strong evidence. From our analysis, we reach opposite conclusions when evaluating the significance with AIC and BIC. For BIC, our results agree with Yang et al. (2016), and we also find evidence for two components with strong evidence (BIC). However, AIC prefers three components with decisive significance.
AIC  BIC  

2  (1.38,0.19)  273  4219.5  4273.4  21.7  7.7  
(4.29,0.28)  718  
3  (0.45,0.06)  527  4197.8  4281.1  
(0.40,0.48)  342  
(3.2,0.08)  122 
4.5 Swift Intrinsic
We now redo the same analysis by looking at T90 and hardness in the rest frame of Swift GRBs with measured redshifts. We first explain how the observed T90 and hardness are converted from the observer frame to the rest frame.
4.5.1 Hardness in rest frame
We have outlined the procedure for the calculation of hardness of GRBs in the observer frame in Sect. 2.2. For GRBs with measured redshifts, we can calculate all parameters in the intrinsic or the GRB rest frame as well. The correction needed to account for the redshift is different for energy and T90 duration and is outlined below. The hardness in the rest frame is defined as,
(11) 
where the only difference is that the energies are multiplied by the correction factor of , being the redshift.
4.5.2 T90 in rest frame
In similar fashion to the above correction due to redshift, the correction to T90 is as follows:
(12) 
4.6 Results
The Swift detector contains an ultraviolet/optical telescope, which is used to slew towards the position of GRB in order to detect the afterglow. This has enabled the measurement of redshifts for a large number of bursts from the Swift catalog The current Swift GRB catalog (Lien et al., 2016) (as of Jan. 2017) contains 373 GRBs with measured redshifts, which is the largest among all detectors. A complete summary of the results on applications of GMM, including the bestfit parameters and their covariance matrices are shown in Table 5. For a twocomponent fit, 123 and 250 GRBs belong to short and long categories respectively. For a threecomponent fit, 36, 191, and 146 belong to short, intermediate, and long categories respectively. The AIC and BIC plots as a function of the number of components can be found in Fig. 9. Here, AIC prefers three components with decisive evidence (AIC ) and BIC prefers two components(BIC) with strong evidence. The ellipses for the three components are shown in Fig. 10. Therefore, the results are the same as those using the observed frame.
When a similar analysis was done by Yang et al. (2016), their results agree with those using the data in the observed, which is that two bursts are preferred when both T90 and hardness are used. The BIC between two and three components is about 22 when hardness is calculated using a simple power law and about 3.4 when a mixed spectrum model is used. Our results are consistent with those obtained by considering the GRB variables in the observed frame, viz. three components are favored using AIC with decisive significance, and two are favored using BIC with strong significance.
AIC  BIC  

2  (1.99,0.92)  123  1782.5  1825.6  19.27  4.1  
(4.3,1.29)  250  
3  (0.23,0.68)  36  1763.2  1829.7  
(2.93,1.08)  191  
(4.63,1.33)  146 
Dataset  

Model preferred  Magnitude  Model preferred  Magnitude  
BATSE  3  88  3  55 
RHESSI  2  4  2  29 
Fermi  2  1.5  2  17 
Swift  3  22  2  8 
Intrinsic Swift  3  19  2  4.1 
5 Conclusions
The main goal of this paper was to find the optimum number of GRB components between two and three categories for a variety of datasets, by carrying out a twodimensional clustering in the T90 vs hardness plane, similar to a recent analysis carried out for the Swift data (Zhang et al., 2016b). We did a comprehensive unified analysis of the T90 versus hardness distributions of GRBs from five different GRB datasets from four spacebased detectors: BATSE, RHESSI, FermiGBM, Swift (observed frame), Swift (intrinsic frame) by fitting the data to two as well as three bivariate normal distributions. We then used two information criterion based statistical tests to ascertain the best model among these two hypotheses. These tests include AIC and BIC model comparison tests. The statistical significance from the information criterion based tests was obtained qualitatively using empirical strength of evidence rules (Shi et al., 2012). Our results are stated below. A tabular summary of all these results are summarized in Table 6.

For BATSE, we find that both AIC and BIC prefer three components with AIC and BIC 10 in both the cases, pointing to decisive evidence for three components.

For RHESSI, we find that both AIC and BIC favor two components. However, the significance from AIC is marginal, but BIC points to decisive significance for two components.

For FermiGBM also both AIC and BIC prefer two components, but with marginal significance for AIC and decisive significance for BIC. We also note that the minimum value of AIC occurs for five components.

When we looked at Swift GRBs in the observer frame, we find that AIC and BIC reach opposite conclusions. AIC prefers three components with decisive significance. BIC shows a prefer for two components, albeit with strong significance. The minimum value of AIC is obtained for five components.

When we repeat the analysis for Swift GRBs in the intrinsic frame, we reach the same conclusions as in the observer frame. AIC prefers three components with decisive significance, where BIC prefers two with strong evidence.
6 Acknowledgements
We are grateful to P. Narayana Bhat for providing us the hardness data for FermiGBM GRBs.
References
 Burnham & Anderson (2004) Burnham, K. P., & Anderson, D. R. 2004, Sociological methods & research, 33, 261
 Chattopadhyay & Maitra (2017) Chattopadhyay, S., & Maitra, R. 2017, Mon. Not. R. Astron. Soc., 469, 3374
 Chattopadhyay et al. (2007a) Chattopadhyay, T., Misra, R., Chattopadhyay, A. K., & Naskar, M. 2007a, Astrophys. J., 667, 1017
 Chattopadhyay et al. (2007b) —. 2007b, Astrophys. J., 667, 1017
 Dempster et al. (1977) Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, Journal of the royal statistical society. Series B (methodological), 1
 Desai (2016) Desai, S. 2016, EPL (Europhysics Letters), 115, 20006
 Desai & Liu (2016) Desai, S., & Liu, D. W. 2016, Astroparticle Physics, 82, 86
 Ganguly & Desai (2017) Ganguly, S., & Desai, S. 2017, Astropart. Phys., C94, 17
 Horváth (1998) Horváth, I. 1998, Astrophys. J., 508, 757
 Horváth (2002) —. 2002, Astron. Astrophys., 392, 791
 Horváth et al. (2010) Horváth, I., Bagoly, Z., Balázs, L. G., et al. 2010, Astrophys. J., 713, 552
 Horváth et al. (2006) Horváth, I., Balázs, L. G., Bagoly, Z., Ryde, F., & Mészáros, A. 2006, Astron. Astrophys., 447, 23
 Horváth et al. (2008) Horváth, I., Balázs, L. G., Bagoly, Z., & Veres, P. 2008, Astron. Astrophys., 489, L1
 Horváth & Tóth (2016) Horváth, I., & Tóth, B. G. 2016, Astrophys. Space Sci., 361, 155
 Huja et al. (2009) Huja, D., Mészáros, A., & Řípa, J. 2009, Astron. Astrophys., 504, 67
 Kass & Raftery (1995) Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 773
 Kouveliotou et al. (1993) Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, Astrophys. J. Lett., 413, L101
 Kulkarni & Desai (2017a) Kulkarni, S., & Desai, S. 2017a, Astrophys. Space Sci., 362, 70
 Kulkarni & Desai (2017b) —. 2017b, ArXiv eprints, arXiv:1708.00605
 Liddle (2004) Liddle, A. R. 2004, Mon. Not. R. Astron. Soc., 351, L49
 Liddle (2007) —. 2007, Mon. Not. R. Astron. Soc., 377, L74
 Lien et al. (2016) Lien, A., Sakamoto, T., Barthelmy, S. D., et al. 2016, Astrophys. J., 829, 7
 Lyons (2016) Lyons, L. 2016, ArXiv eprints, arXiv:1607.03549
 Mukherjee et al. (1998) Mukherjee, S., Feigelson, E. D., Jogesh Babu, G., et al. 1998, Astrophys. J., 508, 314
 Narayana Bhat et al. (2016) Narayana Bhat, P., Meegan, C. A., von Kienlin, A., et al. 2016, Astrophys. J. Suppl. Ser., 223, 28
 Paciesas et al. (1999) Paciesas, W. S., Meegan, C. A., Pendleton, G. N., et al. 1999, Astrophys. J. Suppl. Ser., 122, 465
 Schady (2017) Schady, P. 2017, ArXiv eprints, arXiv:1707.05214
 Shi et al. (2012) Shi, K., Huang, Y. F., & Lu, T. 2012, Mon. Not. R. Astron. Soc., 426, 2452
 Tarnopolski (2016a) Tarnopolski, M. 2016a, Astrophys. Space Sci., 361, 125
 Tarnopolski (2016b) —. 2016b, New Astron., 46, 54
 Řípa et al. (2009a) Řípa, J., Mészáros, A., Wigger, C., et al. 2009a, Astron. Astrophys., 498, 399
 Řípa et al. (2009b) —. 2009b, Astron. Astrophys., 498, 399

Veres et al. (2010)
Veres, P., Bagoly, Z., Horváth, I., Mész
’aros, A., & Balázs, L. G. 2010, Astrophys. J., 725, 1955  Yang et al. (2016) Yang, E. B., Zhang, Z. B., & Jiang, X. X. 2016, Astrophys. Space Sci., 361, 257
 Zhang et al. (2016a) Zhang, B., Lü, H.J., & Liang, E.W. 2016a, ArXiv eprints, arXiv:1611.01948
 Zhang et al. (2009) Zhang, B., Zhang, B.B., Virgili, F. J., et al. 2009, Astrophys. J., 703, 1696
 Zhang & Choi (2008) Zhang, Z.B., & Choi, C.S. 2008, Astron. Astrophys., 484, 293
 Zhang et al. (2016b) Zhang, Z.B., Yang, E.B., Choi, C.S., & Chang, H.Y. 2016b, Mon. Not. R. Astron. Soc., 462, 3243
 Zitouni et al. (2015) Zitouni, H., Guessoum, N., Azzam, W. J., & Mochkovitch, R. 2015, Astrophys. Space Sci., 357, 7