# Bayesian Approach for Counting Experiment Statistics applied to a Neutrino Point Source Analysis

###### Abstract

In this paper we present a model independent analysis method following Bayesian statistics to analyse data from a generic counting experiment and apply it to the search for neutrinos from point sources. We discuss a test statistic defined following a Bayesian framework that will be used in the search for a signal. In case no signal is found, we derive an upper limit without the introduction of approximations. The Bayesian approach allows us to obtain the full probability density function for both the background and the signal rate. As such, we have direct access to any signal upper limit. The upper limit derivation directly compares with a frequentist approach and is robust in the case of low-counting observations. Furthermore, it allows also to account for previous upper limits obtained by other analyses via the concept of prior information without the need of the ad hoc application of trial factors. To investigate the validity of the presented Bayesian approach, we have applied this method to the public IceCube 40-string configuration data for 10 nearby blazars and we have obtained a flux upper limit, which is in agreement with the upper limits determined via a frequentist approach. Furthermore, the upper limit obtained compares well with the previously published result of IceCube, using the same data set.

Vrije Universiteit Brussel, Dienst ELEM

Pleinlaan 2, B-1050 Brussels, Belgium

Key words: Neutrino Astronomy, Neutrino Telescopes, Active Galactic Nuclei, Bayesian Statistics.

## 1 Introduction

Bayesian approaches are gaining more and more popularity in scientific analyses [1, 2, 3, 4, 5]. In this paper we discuss a formalism following Bayesian inference [6, 7] to analyse signals in a generic counting experiment and for illustration we apply it to a point source analysis using data from a neutrino telescope.

A frequentist approach is based on the long run relative frequency of occurrences in identical repeats of an experiment. Consequently, this can only provide the probability for a certain outcome under the assumption of a specific hypothesis. On the other hand, the Bayesian approach allows us to directly compute the probability of any particular hypothesis or parameter value based on observations. One of the great strengths of Bayesian inference is the ability to incorporate relevant prior information in the analysis. This provides a mechanism for a statistical learning process that automatically takes previous results into account.

Concerning a signal detection we will analyse the performance of the presented Bayesian approach and we will compare it to methods frequently found in literature. One of the standard test statistics is the frequentist method developed by Li and Ma [8]. While this method is based on the detection of an excess of number of events in a given angular window, the presented Bayesian method is sensitive to the distribution of the angular distance between the arrival direction of the observed events and the source position.

In case no signal is observed, upper limits for the signal flux from the investigated point sources are determined. In the light of a flux upper limit determination we will discuss a method following the Bayesian framework and a comparison will be made to the standard frequentist methods developed by Feldman and Cousins [9] and Rolke et al. [10]. The Feldman-Cousins method introduces the likelihood ratio as an ordering principle, when determining the acceptance interval from which one derives the confidence interval. This construction resolves the issue of empty confidence intervals when the interval is entirely in the non-physical region. The Feldman-Cousins construction also resolves the “flip-flop” problem, which can lead to under coverage. This arises when one wants to report an upper limit if the evidence for a signal is below a certain threshold, but a central confidence interval otherwise. The procedure outlined by Rolke et al. can deal with several nuisance parameters which include uncertainties (for example background expectations or signal efficiencies) by means of the Profile Likelihood method. However, this method has the disadvantage that large samples are assumed and thus applying it to cases with low counting statistics might lead to under coverage.

The Bayesian method presented here will be applied to the study of Active Galactic Nuclei (AGN), which are, together with Gamma Ray Bursts (GRBs), among the leading candidates for the sources of ultra-high energy cosmic rays (UHECRs), as outlined in recent reviews such as [11, 12]. On time scales of order of minutes or even seconds GRBs are transient phenomena, whereas AGN in general can be regarded as steady sources over time periods even exceeding several years. In case hadronic acceleration takes place in these objects, also an accompanying high-energy neutrino flux is expected due to the decay of secondary particles produced in the interactions of accelerated hadrons with the ambient photon field or matter [13].

Kilometer-scale neutrino detectors have the sensitivity to measure the predicted neutrino flux. The IceCube Neutrino Telescope [14], located at the geographic South Pole, was completed at the end of 2010 and has been taking data since 2006. In the Northern Hemisphere a kilometer-scale detector called KM3Net is being proposed for construction [15]. There are smaller telescopes, one located in the Mediterranean Sea, named ANTARES which has been collecting data since 2008 [16] and another one in Lake Baikal, Russia, which has been operating its NT configuration since 2005 [17]. These telescopes detect high energy neutrinos ( GeV) by observing Cherenkov radiation in ice or water from secondary particles produced in neutrino interactions.

In the following we introduce the basis of the Bayesian formalism and present the test statistic that we will use to assess consistency of the data with the null hypothesis, being an isotropic distribution of the events. In Section 4, we describe the procedure to determine upper limits in case no significant signal is observed. Subsequently, in Section 5, we apply the method to IceCube public data [18] investigating the neutrino flux from the ten closest blazars (a subclass of AGN). Finally, a summary of our results and conclusions are provided in Section 6.

## 2 Bayesian Formalism

We denote by the probability for hypothesis and to be true under the condition that is true. The product rule [1], , directly yields Bayes’ theorem [19]:

(1) |

Bayes’ theorem is extremely powerful in hypothesis testing. Consider a hypothesis , some observed data and prior information . Bayes’ theorem can then be rewritten as:

(2) |

where

From Eq. (2) it is seen that the Bayesian formalism automatically provides a learning process. The first step is to encode our state of knowledge before analyzing the data into a prior probability . This is then converted into a posterior probability when a new data set is analysed.

We will apply this approach to the analysis of astrophysical point sources. This method will enable us to obtain the full Probability Density Function (PDF) for the source rate (i.e. the number of signal events per unit of time arriving from the source) from which we can derive the corresponding flux or any upper limit.

## 3 Assessment of Significance

As outlined above, Bayesian inference allows us to make statements about the probability of various hypotheses in the light of obtained data. Following the developments described in [20], we quantify the degree to which data support a certain hypothesis and as such make an assessment of the significance. To quantify our degree of belief in a certain hypothesis , one can use the so called evidence [1]

(3) |

where indicates hypothesis to be false. Due to the , this evidence is in a decibel scale. One can calculate [20] that the evidence for any alternative to hypothesis based on the data and prior information is constrained by the observable

(4) |

Thus, provides a reference to quantify our degree of belief in .

Let us now consider an experiment where the probabilities corresponding to the various outcomes on successive trials are independent and stationary. Such experiments belong to the so called Bernoulli class [1]. The probability of observing occurrences of each outcome after trials is therefore given by the multinomial distribution:

(5) |

In terms of the observable we then obtain [20] for each :

(6) |

which is an exact expression.

In case the data are represented in histogram form, the above can be applied if is the total number of entries, represents the number of bins, is the number of entries in bin and is the probability for an entry to fall in bin . Once the various are known, the -value corresponding to a certain observed distribution can easily be obtained using Eq. (3).

In our investigation of neutrino point sources, the total number of events will populate a histogram according to the angular difference of each neutrino arrival direction with respect to the source location in the sky. The values are determined assuming an isotropic background and taking into account the solid angle effect. By generating randomly the same total number of events\@footnotemark\@footnotetextThe number has to be fixed due to the first term in the definition of , Eq. (3). This implies that is sensitive to the shape of the distribution of events within the distance but not to the total number of events, i.e. if there is an excess of signal events that follows an isotropic distribution, the observable will not enable us to distinguish it from background.following an isotropic distribution, we obtain the distribution of for an isotropic background. Comparing the value of for the data with the obtained background -distribution, we determine the P-value or significance of our measurement. A source detection is claimed if the consistency of the data with the null hypothesis has a P-value smaller than , corresponding to a effect in case of a positive single sided Gaussian distribution.

It is important to note that as is sensitive to the distribution of the data in the histogram, it is binning-dependent. To avoid this binning effect in the final P-value, the bin size is chosen such that there is only one or zero event per bin. In this special, quasi-unbinned approach, Eq. (3) would simplify in the following expression

(7) |

where the sum is running over the total number of events .

The probabilities are then the crucial factor that will differentiate between background events following an isotropic distribution and signal events from a source. The signal events will be located at angular distances around the source position following a Gaussian distribution with its standard deviation given by the experimental angular resolution convoluted with the solid angle effect.

On the other hand, the test statistic developed by Li and Ma [8] is based on the absolute value of the number of events in the on-source () and off-source () angular windows. Consequently, both test statistics behave differently as a function of the considered angular window. As we will show in Section 5, the Li-Ma test statistic performs better than for small angular windows while is more sensitive at larger windows. Note that in Section 5 the complete expression of Eq. (3) is used instead of the simplified Eq. (3) since for our background distributions it occurs that in a small fraction of cases, which we control to be at maximum 5%.

In case the observation does not lead to a significant detection, we determine an upper limit on the signal strength. To obtain such an upper limit we will use an exact analytical expression following a Bayesian approach using a uniform prior. The details of this procedure are provided in the following section.

## 4 Bayesian upper limit determination

To obtain an upper limit for a possible source flux, we first have to determine the upper limit for the source rate. Using Bayesian inference we obtain the full posterior PDF for the source rate and from that we can derive the corresponding flux PDF, via the concept of effective area, as explained hereafter.

In any experiment where events are detected at a known rate and independently of the time since the last event, the PDF for the number of observed events is described by a Poisson distribution. Our approach to obtain an upper limit follows the one outlined in [6, 7], except that we apply the exact analytical expression without any approximation.

The data consist of on-source and off-source measurements, where the off-source data consist only of background events and the on-source data are a mix of background and source (i.e. signal) events. We start with the off-source analysis and subsequently use the obtained information as background prior information in the on-source analysis. This implies that the current study is completely data driven and as such is a model independent search for a possible source signal.

### 4.1 Off-source measurements

Consider the case that in an off-source measurement background events have been recorded over a time interval with a constant background rate . Using Eq.(2) we obtain the posterior background PDF by:

(8) |

In the above equation the likelihood function is given by the Poisson distribution corresponding to the measurement of background events over a time span at a constant rate :

(9) |

Since the integrated PDF amounts to 1 (i.e. ), the normalisation factor appearing in Eq. (4.1) is given by:

(10) |

As mentioned before, is the prior PDF for the background rate. For our analysis we use a uniform prior [7], which is given by

(11) |

The uniform prior attributes the same probability to every value of the background rate within the indicated range, reflecting that we do not favour a particular value of the actual background rate. This prior also has the advantage that the derived upper limits are directly comparable to classical frequentist upper limits [5].

To cover the full range of possible background rates, the minimum value of the rate is taken to be zero and Eq. (4.1) can be written as

(12) |

Using this expression for together with Eqs. (4.1) and (4.1) we obtain an analytical expression for the normalisation factor

(13) |

Solving the above equation (for details see Appendix A1) we get:

(14) |

where is the Incomplete Gamma Function. Substitution of Eqs. (4.1), (4.1) and (4.1) in Eq. (4.1) yields

(15) |

which represents the posterior background rate PDF. In [6, 7] the Incomplete Gamma Function is then approximated to . This approximation is valid for and cannot be applied in general when the number of events and time window are small (see [20] for an example). We therefore work with the complete analytical expression.

### 4.2 On-source measurements

Consider the case that in an on-source measurement events, consisting of signal and background, have been recorded over a time interval with a constant signal rate and background rate . Following Eq. (2) the joint probability of source and background is given by:

(16) |

Using the product rule [1] we can write the above equation as,

(17) |

where is the prior probability for the background rate, which is in our case the posterior background PDF obtained from the off-source measurement reflected in Eq. (4.1). The likelihood function is the Poisson distribution for the combined signal and background rate . The normalisation constant is obtained, as outlined in the previous subsection, by integrating the numerator of Eq. (4.2).

It is important to note that since the source rate and the background rate are independent, we can write . Like for the background case discussed before, we use a uniform prior for , i.e.

(18) |

As mentioned before, the uniform prior attributes the same probability to every value of the signal rate within the indicated range, reflecting that we do not favour a particular value of the actual signal rate and also allows us to directly compare the derived upper limits with the classical frequentist results [5].

By substituting in Eq. (4.2) the expressions of Eq. (4.2) for , (4.1) for and (4.1) for a total rate we obtain the joint PDF for the source and background rates:

(19) |

However, we are interested in the posterior PDF for the source rate alone, independent of the background. The Bayesian formalism allows us to obtain this posterior PDF by marginalisation [1] of the joint PDF, Eq. (4.2), with respect to the background i.e.

(20) |

Solving the above integral we obtain an exact expression for the source rate posterior PDF (for details see Appendix A2):

(21) |

where and . As explained before we use this exact analytical expression without the approximation used in [6, 7].

In case no significant source signal is observed, Eq. (4.2) allows us to derive any upper limit for the source rate. As an example, the 90% source rate upper limit is given by:

(22) |

## 5 Application of the method to a neutrino point source analysis

The approach presented in the current paper is not exploiting a potential source variability when it comes to providing upper limits to a possible source rate, as outlined in Section 6. Consequently the analysis presented here is tailored for sources with a steady rate within the detection time. Our primary goal is the analysis of AGN which, at the time scale considered here, may be regarded to be steady sources of high-energy neutrinos if hadronic acceleration takes place in these sites. By selecting a small region around every well known source location, several sky patches are defined from which data were collected to search for a possible deviation from the background “noise”.

To validate the analysis procedure described in this report, we use the public data [18] of the muon neutrino candidate events recorded by the IceCube Neutrino Observatory [14, 21] in its 40 string configuration (IC40), that collected data during the season 2008-2009. Our analysis is performed on ten nearby blazars (a special class of AGN with one of the jets pointing in the direction of the Earth) following the approach described in the previous sections. The blazars were selected from the online “Roma BZCAT Multi-frequency Catalogue of Blazars” [22] and are listed in Table 1. These blazars are chosen to be nearby, i.e. with a small redshift, and in such a way that their respective angular windows are not overlapping. We limit ourselves to sources in the Northern hemisphere to reduce the atmospheric muon background for the IceCube measurements.

Blazar Name | Right Ascension | Declination | Redshift | Distance [Mpc] |
---|---|---|---|---|

BZUJ1148+5924 | 177.20983 | 59.41567 | 0.011 | 46.2 |

BZUJ0048+3157 | 12.19642 | 31.95697 | 0.015 | 63 |

BZUJ0319+4130 | 49.95067 | 41.51169 | 0.018 | 75.6 |

BZUJ0709+5010 | 107.39246 | 50.18225 | 0.02 | 84 |

BZUJ0153+7115 | 28.35771 | 71.25181 | 0.022 | 92.4 |

BZUJ1719+4858 | 259.81025 | 48.98042 | 0.024 | 100.8 |

BZUJ1632+8232 | 248.13321 | 82.53789 | 0.025 | 105 |

BZUJ1755+6236 | 268.95183 | 62.61225 | 0.027 | 113.4 |

BZBJ1104+3812 | 166.11379 | 38.20883 | 0.030 | 126 |

BZBJ1653+3945 | 253.46758 | 39.76017 | 0.033 | 139.3 |

### 5.1 Assessment of significance

As outlined in [20], we stack the recorded events within a given angular window centered on each of these ten blazars according to their angular distance from the actual blazar position. As mentioned in the previous sections, we need an expression for the various probabilities of Eq. (3) to derive the -value of the signal. In our case, the background is isotropic and consequently the probabilities have to be consistent with the solid angle effect within the selected cone, i.e. :

(23) |

where is the width of each bin in our stacked histogram and is the size of the angular window.

To determine the angular window size for which our test-statistic is most sensitive, we have generated 266 (being the number of events of the actual observation as outlined hereafter) isotropically distributed events plus a signal of 20 events. These 20 events were generated such that their angular distance from the source position follows a Gaussian distribution (with standard deviation of 1 degree) convoluted with the solid angle effect. The chosen standard deviation of the Gaussian distribution is the angular uncertainty of the IceCube track reconstruction [21, 23]. As we mentioned before, our test statistic is sensitive to the distribution of events as a function of the distance to the source. To avoid the dependence on the position generation in our determination of the most optimal angular window size, we have repeated the procedure 10 times for each value of for the signal events and present the mean P-values in Fig. 1.

From Fig. 1 it is seen that our test statistic does not perform well on small angular windows. This is due to the fact that the difference of the individual probabilities per bin () for small angular windows is not large enough to distinguish a source-like distribution from an isotropic background (i.e. a larger angular window is needed to see the “shape” of the excess). For an angular window of 4 degrees, we see that the sensitivity obtained with becomes optimal.

We also include in Fig. 1 a comparison to the P-value obtained with the standard Li-Ma method. As mentioned in Section 3, Li-Ma is a test statistic based on the total number of events in the on-source () and off-source () angular windows. Li-Ma performs better than for small angular windows (a large ratio). When comparing the smallest P-values of each test statistic, we see that both are similar.

Fixing to 4 degrees, we use the IceCube public data and we obtain the stacked distribution of events presented in Fig. 2. These stacked data comprise 266 events recorded over a time period of 375.5 days [24]. Using Eqs. (3) and (5.1) for a number of entries , the data represented in Fig. 2 yield dB.

As explained in Section 3, the distribution in the case of an isotropic background is obtained by randomly generating times the same number of events as in the on-source region. The distribution is presented in Fig. 3. Comparison of the actual observation with the background distribution gives a -value of . Consequently, we will proceed to give an upper limit on the signal strength.

### 5.2 Upper limit determination using Uniform priors for Source and Background

#### 5.2.1 Determination of the background rate

To determine the number of events in the off-source region () we perform measurements in 4 regions of the sky, shifted from the various blazars positions only in right ascension, keeping the declination constant due to the declination dependence of the IceCube detection efficiency. The specific IC40 configuration of IceCube is also right ascension dependent, so we make shifts of 180 in right ascension to eliminate the right ascension dependence. The IC40 sample has been taken over a detector live time period of 375.5 days, so that both the exposures for on-source, , and off-source, , amount to 375.5 days. The stacked off-source measurements yield a total of 265 events. The posterior background rate PDF is obtained by substitution of the previously mentioned values of and in Eq. (4.1) and by using a sufficiently large value 0.8 mHz. The resulting background rate PDF is shown in Fig. 4.

#### 5.2.2 Determination of the source rate

The posterior source rate PDF is obtained by inserting the previously mentioned values of , , and in Eq. (4.2) and by using a sufficiently large value Hz. The resulting source rate PDF is shown in Fig. 5.

Using the PDF shown in Fig. 5 and applying Eq. (4.2), we obtain the 90% upper limit for the source rate:

To compare this Bayesian method with the frequentist approach, we have also determined the 90% upper limit for the source rate using the Feldman-Cousins [9] and Rolke et al. [10] methods. The values obtained are the following:

We see that the Bayesian approach is equal to the Rolke et al. method and is more conservative than the Feldman-Cousins method. To further test the upper limit calculation, we have generated source signals (or under-fluctuations) by increasing (or decreasing) the number of events in the on-source region, while keeping the same number of background events as in the data. In Fig. 6 we plot the rate upper limits obtained with the Bayesian and frequentist methods as a function of the difference of the number of events between the on-source and off-source regions (). We also show the actual rate for the case of a positive difference of . The Bayesian upper limits are similar to the Rolke et al. limits for a small difference of but the former is more restrictive when this difference increases. When comparing to the Feldman-Cousins results, the Bayesian limits are more conservative for low and tend to the Feldman-Cousins limits as this difference grows, as expected from the fact that we use an uniform prior.

The decrease of the slope when the difference of the number of events is negative shows that the Bayesian method is better protected against under-fluctuations. This effect is shown in Fig. 7, where we consider the background fluctuations by generating isotropic distributions of events in the sky and compute each time and (with , the generated source events). Fig. 7 shows the computed event rate upper limit for the Bayesian and Feldman-Cousins methods as a function of the generated events. We see that the decrease of the slopes of the upper limit determinations (Fig. 6) result in an upper limit that can fall below the actual generated rate. This problem occurs less often for the Bayesian method because the decrease of the slope is less steep compared to Feldman-Cousins. Moreover, as we assume an uniform prior, the Bayesian limit is equal to the Feldman-Cousins for large over-fluctuations and this translates to Fig. 7 by having the same values for the largest upper limits of the rate.

Note that the obtained rate does not take into account the reconstruction efficiency. The latter is taken into account by converting the source rate upper limit into a flux upper limit by means of the so called Effective Area, , which is defined as

For the current analysis we use the angle averaged Effective Area determined from a simulated spectrum [24], taking into account the observed energy estimate for each individual observed event [18]. The median value corresponds to Acm over the considered energy range. Our analysis is performed on a circular area of 4 centered on each of the 10 sources, representing in total . From the result of and taking the effective area and the size of the on-source region into account we arrive at a 90% upper limit for an signal flux of .

However, this flux upper limit does not take into account the effect of neutrino oscillations. At the source, astrophysical models predict a flavor ratio of . Assuming maximum oscillation we expect to observe at Earth . A tiny fraction of the will produce a muon which might also be detected in IceCube and as such have entered our event sample. However, we will neglect this effect since it is marginal and would require a special simulation which is beyond the scope of this paper. So our final value for the 90% upper limit for a E signal flux is:

For consistency checking, we can compare our limit to the result published by the IceCube Collaboration [24], which was obtained with a different analysis concerning a search for a diffuse high-energy neutrino flux in the full Northern hemisphere: scmsr, which is comparable to our result. In that analysis the same data set was used and since the ten blazars we studied are randomly located in the sky and have not been selected based on any (astro)physical characteristics, the windows around them represent a fair sample of the sky which can be used to compare to a diffuse search.

## 6 Conclusion and Outlook

In this paper we have discussed a statistical method to analyse point sources using data from a neutrino telescope following Bayesian inference. Using the observable , we have indicated how to assess the significance of a possible signal in the data by comparing it to the distribution expected for an isotropic background. We have shown how to obtain upper limits for the corresponding flux in case the observation does not lead to a significant signal detection. Our calculations are similar to [6, 7] but we have made no approximations in the final results and thus this method can be applied to low counting observations.

Applying this method we have analysed the public IceCube 40-string configuration data for 10 nearby blazars located in the Northern sky. From our analysis it was also seen that the on-source data is consistent with an isotropic background only hypothesis. Therefore we have determined a 90% upper limit, which has been compared to the upper limits obtained using the same data set but applying the Feldman-Cousins and the Rolke et al. methods. Simulating a signal from a source, by artificially changing the number of events in the on-source region, we have shown that the Bayesian limits are similar to the Rolke et al. calculations for small difference in the number of events between the on-source and off-source region and tend to the Feldman-Cousins limits as this difference in the number of events increases. We have shown that in the case of under-fluctuations in the background the Bayesian method is better protected.

It is our intention to extend the current method also for non-steady sources like for instance GRBs. Apart from providing a signal significance for discovery [20], this should also provide a mechanism to accurately determine flux upper limits. For flaring sources we do not know the time window in which the neutrinos are emitted. If we take a time window large enough to cover all possible scenarios for neutrino emission, we would obtain a rate which is not the actual one (because of the existence of time intervals with and without neutrino emission within our time window). The proper extension of the method is currently under study.

## Acknowledgements

The authors would like to thank the IceCube Collaboration for providing the public data used in this report to evaluate our analysis method. This research was performed with financial support from the Odysseus programme of the Flemish Foundation for Scientific Research (FWO) under contract number G.0917.09.

## Appendix A1

The normalisation factor of the posterior background rate PDF given in Eq. (4.1) may be written as:

(24) |

The integral part can be expressed as the so-called Incomplete Gamma function given by:

(25) |

Using this expression, we can rewrite Eq. (24) as follows:

(26) | |||||

(27) |

which is the expression reflected in Eq. (4.1).

## Appendix A2

According to Eq. (4.2) the posterior PDF for the source rate alone is given by:

(28) |

where

Substitution of the various expressions given in Section 4 yields:

(29) |

Combination of Eqs. (28) and (29) yields:

Considering the numerator of Eq. (Appendix A2), we obtain:

A | ||||

Using the Newtonian Binomial, , we find:

Using the Incomplete Gamma function, Eq. (25), and simplifying the above equation, we finally find

(31) |

where and .

Applying the same procedure, we obtain for the denominator:

B | ||||

Which finally gives, after simplification, for the posterior source rate probability density function:

(32) |

## References

- [1] E.T. Jaynes, Probability Theory, Cambridge University Press, Cambridge, 2003.
- [2] E. D. Feigelson and G. J. Babu, Statistical Challenges in Modern Astronomy III, Springer, Heidelberg, 2002.
- [3] C. Arina, Bayes and present dark matter direct search status, J. Phys. Conf. Ser. 375 (2012) 012009 [arXiv:1110.0313].
- [4] G. D’Agostini, Bayesian reasoning in high-energy physics: Principles and applications, CERN-YELLOW-99-03, Geneva, 1999.
- [5] R. D. Cousins, Why isn’t every physicist a Bayesian?, Am. J. Phys. 63 (1995) 398.
- [6] T. J. Loredo, chapter 5.2, Statistical Challenges in Modern Astronomy, E.D. Feigelson and G.J. Babu (eds.), Springer- Verlag, New York, pp. 275-297 (1992).
- [7] P. Gregory, Bayesian Logical Data Analysis for the Physical Sciences, Cambridge University Press, Cambridge, 2010.
- [8] T. P. Li and Y. Q. Ma, Analysis methods for results in gamma-ray astronomy, Astrophys. J. 272 (1983) 317-324.
- [9] G. J. Feldman and R. D. Cousins, A Unified Approach to the Classical Statistical Analysis of Small Signals, Phys. Rev. D 57 (1998) 3873-3889.
- [10] W. A. Rolke et al., Limits and Confidence Intervals in the Presence of Nuisance Parameters, Nucl. Instrum. Meth. A 551, p. 493-503 (2005).
- [11] A. Letessier-Selvon and T. Stanev, Ultrahigh Energy Cosmic Rays, Rev. Mod. Phys 83 (2011) 907-942 [astro-ph:1103.0031].
- [12] K. Kotera, A. V. Olinto, The Astrophysics of Ultrahigh Energy Cosmic Rays, Ann. Rev. Astron. Astrophys. 49 (2011) [astro-ph/1101.4256].
- [13] E. Waxman, J. N. Bahcall, High energy neutrinos from astrophysical sources: An upper bound, Phys. Rev. D 59 023002 (1999)
- [14] F. Halzen and S. R. Klein, IceCube: An Instrument for Neutrino Astronomy, Rev. Sci. Instrum. 81 (2010) 081101 [arXiv:1007.1247v2].
- [15] www.km3net.org
- [16] J. Aguilar et al., the ANTARES Collaboration, ANTARES: the first undersea neutrino telescope, Nucl. Instrum. Meth. A 656 (2011) 11-38.
- [17] V. Aynutdinov et al., the Baikal Collaboration, The Baikal neutrino experiment: Status, selected physics results and perspectives, Nucl. Instrum. Meth. A 588 (2008) 9910.
- [18] http://icecube.wisc.edu/science/data/ic40.
- [19] T. Bayes and Mr. Price, An Essay towards solving a Problem in the Doctrine of Chances, Phil. Trans. 53 (1763) 370Ð418.
- [20] N. van Eijndhoven, On the Observability of High-Energy Neutrinos from Gamma Ray Bursts, Astropart. Phys. 28 (2008) 540-546 [arXiv:astro-ph/0702029v2].
- [21] J. Ahrens et al., IceCube Collaboration, Sensitivity of the IceCube detector to astrophysical sources of high energy muon neutrinos, Astropart. Phys. 20 (2004) 507-532 [arXiv:astro-ph/0305196].
- [22] http://www.asdc.asi.it/bzcat/.
- [23] R. Abbasi et al., IceCube Collaboration, The Shadow of the Moon in Cosmic Rays measured with IceCube, Proc. of the 32nd International Cosmic Ray Conference V4 (2011) 264-267 [arXiv:1111.2741].
- [24] R. Abbasi et al., IceCube Collaboration, A Search for a Diffuse Flux of Astrophysical Muon Neutrinos with the IceCube 40-String Detector, Phys. Rev. D 84 (2011) 082001 [arXiv:1104.5187].