LETTER TO THE EDITOR
The Htest probability distribution revisited: improved sensitivity
O. C. de Jager^{1,3}  I. Büsching^{1,2}
1  Unit for Space Physics, School of
Physics, NorthWest University, 2520 Potchefstroom, South Africa
2 
Centre for High Performance Computing (CHPC), CSIR Campus, 15 Lower Hope St. Rosebank, Cape Town, South Africa
3 
South African Department of Science & Technology and National Research Foundation
Reseach Chair: Astrophysics & Space Science, South Africa
Received 5 March 2010 / Accepted 28 May 2010
Abstract
Aims. To provide a significantly improved probability distribution for the Htest for periodicity in Xray and ray arrival times, which is already extensively used by the ray
pulsar community. Also, to obtain an analytical probability
distribution for stacked test statistics in the case of a search for
pulsed emission from an ensemble of pulsars where the significance per
pulsar is relatively low, making individual detections insignificant on
their own. This information is timely given the recent rapid discovery
of new pulsars with the FermiLAT t ray telescope.
Methods. Approximately 10^{14} realisations of the Hstatistic (H) for random (white) noise is calculated from a random number generator for which the repetition cycle is 10^{14}. From these numbers the probability distribution P(>H) is calculated.
Results. The distribution of H is found to be exponential with parameter
so that the cumulative probability distribution
.
If we stack independent values for H, the sum of K such values would follow the ErlangK distribution with parameter
for which the cumulative probability distribution is also a simple analytical expression.
Conclusions. Searches for weak pulsars with unknown pulse
profile shapes in the FermiLAT, Agile or other Xray data bases should
benefit from the Htest since it is known to be powerful
against a broad range of pulse profiles, which introduces only a single
statistical trial if only the Htest is used. The new
probability distribution presented here favours the detection of weaker
pulsars in terms of an improved sensitivity relative to the previously
known distribution.
Key words: methods: statistical  pulsars: general
1 Introduction
When searching for a periodic signal in Xray or ray arrival times dominated by noise, we may either perform a blind search for ray pulsars as demonstrated by the FermiLAT Collaboration (Abdo et al. 2009), or, search for such a signal where the frequency parameters have been prescribed by contemporary radio data (Weltevrede 2010). Following the folding of say N arrival times t_{1}, ...,t_{N} modulo the pulsar spin parameters, we arrive at a set of phases , i=1, ...N. However, in the case of blind searches, Atwood et al. (2006) introduced a time differencing technique in which case the number of trial periods is significantly reduced.de Jager et al. (1989, hereafter DSR) reviewed the general class of Beran statistics (Beran 1969), from which the most general test statistics such as Pearson's , Rayleigh and Z^{2}_{m} statistics are derived, and from within this class they derived the well known Htest for Xray and ray Astronomy.
The probability distribution of the Htest statistic as given by DSR was derived from Monte Carlo simulations employing 10^{8} simulations. The computational power and random number simulators on typical IBM machines during the 1980's had limited ranges of applicability and the Htest suffered accordingly. For values of H<23 we found that the probability distribution was exponential with parameter (or 0.4), whereas a hard tail developed for H>23, which resulted in a significant compromise in sensitivity.
The old version of the Htest probability distribution is already extensively used by e.g. the FermiLAT Collaboration for pulsar searches (e.g. Abdo et al. 2009; and Weltevrede et al. 2010), and from this paper it will become clear that the significances assigned to pulsar detections (or nondetections) may be too conservative, so that some pulsars may be missed, especially when many trial periods are involved, so that large values of the Hstatistic are required for a significant detection. In this Letter we notify the community that all previous published significances from the Htest should be reassessed, based on the new improved distribution presented below. Before we do so, we first briefly review the origin and properties of the Htest.
2 The Beran class of test statistics  towards the Htest
Let
be the pulsar phase measured on the interval ,
so that a full rotation corresponds to .
Assuming noise (e.g. from cosmic rays) are also present such that the pulsed fraction is .
Let
be the observed lineofsight
pulse profile in the absence of noise. The case p=0 then corresponds to no signal (pure noise), whereas p=1 corresponds to no noise (pure pulsed signal). The observed light curve
can then be represented as a mixing of the noise and signal distributions
(see Eq. (2) of DSR 1989)
(1) 
The general form of the Beran statistic is given by Beran (1969) in the form (see also DSR 1989)
It is thus clear that the Beran statistic measures the integrated squared distance between the pulse profile and uniformity, so that if , then as well, or, if (i.e. a flat uniform distribution), then as well. Thus, we would reject uniformity if exceeds a chosen positive critical value.
It was noted by DSR that by replacing
with a density estimator
based on the observed folded phases ,
we retrieve test statistics specific to the kernel of
the density estimator.
Selecting the Fourier series
estimator
with m harmonics (see Eq. (5) of DSR) as representative of the light curve,
results in the wellknown Z^{2}_{m} given the proper normalisation (Eq. (7) of DSR)
(3) 
Since scales with p^{2} (Eq. ( 2)), it is then clear that the statistic Z^{2}_{m} should scale as p^{2}N. The quantity is then also the approximate Gaussian significance of the point source on the skymap if the background level is well known.
The main problem raised by DSR is that we do not know apriori the optimal number of harmonics to select. In the case of the popular Z^{2}_{m} test, the optimal number of harmonics would depend on X and the pulse profile shape.
The rationale behind the Htest was to find a consistent estimator for where the number of harmonics m is not chosen subjectively by the observer, since selecting a number of m values for Z^{2}_{m}, until a pleasing result is obtained, may lead to a false detection, given the difficulty to keep track of the number of trials involved.
Hart (1985) derived a technique to obtain the optimal number of harmonics M such that the
mean integrated squared error (MISE) between the Fourier series estimator
and the
true unknown light curve shape is a minimum (see DSR for a review). Thus,
would give the best Fourier series representation of the true pulse profile, given the constraints
imposed by the statistics and inherent pulsed fraction. Since the minimum of the MISE involves the finding of
a maximum quantity over all harmonics m, DSR redefined this maximal (optimal) quantity in terms of the
socalled Htest statistic (named after Hart):
(4) 
Uniformity would correspond to low values of M and also low values of H, whereas large values of Z^{2}_{M} (relatively strong pulsed signal) corresponds to large values of H. Subtracting 4M and adding 4 (to ensure positive values of H only) effectively limits the variance of H in the absence of a pulsed signal. It is clear that the H test represents a rescaled version of the Z^{2}_{m} test: if M=1, then we retrieve the wellknown Rayleigh test. DSR have shown that this test is powerful against a wide range of alternatives, which makes this test attractive if the pulse profile shape is not known a priori.
To give an example of how H and M depend on the abovementioned significance given real data, we selected archival photon arrival times from the public FermiLAT data base from the Vela pulsar direction, with events selected according to the FermiLAT point spread function, but with energies >500 MeV. Folding these phases with the given contemporary spin parameters of Vela reveals a pulsed fraction of , i.e. nearly 100% and from N=600 events (4 days integration) we already obtain M=20. By adding randomly generated events (i.e. uniformly distributed pulse phases) to the signal events we effectively reduce p and hence X, so that H and M should also decrease accordingly. Figure 1 shows how H and M relates to X: For X>3 (stronger than a signal) we see that H scales with X^{2}, with H=1.9X^{2}, whereas the optimal M is already >10 for X>7.
Figure 1: The dependence of H and M on the approximate skymap significance X for ray pulse phases above 500 MeV from the Vela pulsar as described in the text. The errors represent standard deviations based on 20 independent FermiLAT Vela pulsar data sets, each of length up to 4 days. The horisontal line represent the upper limit of M=20, whereas the dashed line represent a fit of the form H=1.9X^{2} for X>3. 

Open with DEXTER 
This figure is also quite useful to see what typical H and M values we may expect for a typical Velalike pulsar if we know the strength of the signal as derived from a point source on the skymap, and assuming the excess is due to such a pulsed signal.
3 The revised probability distribution of H for uniformity
The calculations were performed on the Institutional Cluster of the NorthWest University, Potchefstroom campus. To parallelize the computations, we used the RngStream package (L'Ecuyer 2002), which guarantees independent, nonoverlapping substreams of random numbers. The repetition cycle for this random number generator is , which is certainly large enough for our purposes. A total of samples were calculated in this way.
In the case of large statistics (N>100) we do not need to simulate individual arrival times directly (see the approximate correction factors in Figs. 1 and 2 of DSR in the case of low statistics  N<100), so that we only need to simulate the Z^{2}_{m} statistics directly: since Z^{2}_{m} is the sum of m statistics, we can simulate a statistic directly from a uniformly distributed random number by taking the transformation and adding these numbers to give Z^{2}_{m}. This speeds up the process considerably. A total of values of H were simulated in this way and the results are shown in Fig. 2. The distribution is everywhere consistent with an exponential distribution with parameter (or 0.4), except for H>70 where a downturn relative to the 0.4 index is possibly seen. DSR arrived at the same precise value of for small values of H (i.e. less than 23) since the random number generator used by DSR did not yet reach the limit of its random cycle for the number of simulations required to reach (with a relatively small error on the corresponding probability) and should therefore reveal the same result as ours for H<23.
It is thus clear that the probability distribution of the
Hstatistic
can be conservatively described by the simple formula
(5) 
Figure 2: The distribution of the Hstatistic derived from Monte Carlo simulations with the best fit model for the cumulative probability (solid line) and the old version of the probability distribution given by DSR shown as a dashed line. 

Open with DEXTER 
4 Incoherent stacking
Suppose we analysed the data from K pulsars, or, K independent observations of the same pulsar, where the effect of e.g. unrecorded glitches excluded the possibility of analysing all the data as one single coherent set of arrival times. In this case we want to see if there is a net signal represented by K values of the Hstatistic. Suppose we arrive at a set of K such values of H, given by H_{i}, i=1, ...,K.
Since we have shown (to the probability level of 10^{14}) that the
Hstatistic follows that of an exponential distribution with parameter
,
we can stack these test statistics through the sum
(6) 
which is known to follow the ErlangK distribution with parameter , so that the significance (or probability for uniformity) of such stacking is given by the simple analytical expression (see e.g., Leemis & McQueston 2008, on univariate distribution relationships, which includes the Erlang distribution as the sum of independent exponential variables with parameter .)
(7) 
5 Conclusions
For M=1 we retrieve the wellknown Rayleigh test, with the exception that the parameter for the exponential distribution has been reduced from (for the Rayleigh test) to for the Htest. This slight loss of sensitivity is the effect of the number of trials taken implicitly into account as a result of the search through m within the Htest.
Finally, it is clear that the corrected distribution of the Hstatistic follows a simple exponential with parameter and evaluation of results for H>23 (i.e. the 10^{4} significance level) would yield more significant results compared to the old distribution presented by DSR. For example, for H=50 a probability of is typically quoted in the literature, whereas the true probability for uniformity is actually  already a factor 20 smaller.
A FermiLAT example of the Vela pulsar (>500 MeV) shows clearly values for M up to 20 for ``skymap'' significances , whereas H scales with as expected for Berantype tests. The scaling H=1.9X^{2} can be used to predict Htest statistics for Velalike pulsars above 500 MeV if we assume the excess on the skymap is all pulsed.
The hard tail of the distribution beyond H>23 presented by DSR probably arose from the repetition cycle of the random number generator used in those days, so that the same fluctuations at large H values were repeated given the finite cycle length of the generator used. In this case we however used a generator with a cycle time much longer than 10^{14}, in which case we did not see the repetition of outliers as a result of a finite cycle length. Confirmation of the possible break (i.e. downturn) in the probability distribution at H>70 requires extensive simulations beyond and is beyond the scope of this paper.
AcknowledgementsThe authors are grateful for partial financial support granted to them by the South African National Research Foundation (NRF) and by the Meraka Institute as part of the funding for the South African Centre for High Performance Computing (CHPC).
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 840 [NASA ADS] [CrossRef] [Google Scholar]
 Atwood, W. B., Ziegler, M., Johnson, R. P., & Baughman, B. M. 2006, ApJ, 652, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Beran, R. J. 1969, Ann. Math. Stat., 40, 1196 [CrossRef] [Google Scholar]
 Buccheri, R., & de Jager, O. C. 1989, in NATO ASI Workshop on Timing Neutron Stars, ed. H. Ögelman, & E. P. J. van den Heuvel (Dordrecht: Kluwer), 95 [Google Scholar]
 de Jager, O. C., Swanepoel, J. W. H., & Raubenheimer, B. C. 1989, A&A, 221, 180 (DSR) [NASA ADS] [Google Scholar]
 Hart, J. D. 1985, J. Statist. Comput. Simul., 21, 95 [CrossRef] [Google Scholar]
 Leemis, L. M., & McQueston, J. T. 2008, American Statistician, 62, 45 [CrossRef] [Google Scholar]
 L'Ecuyer, P., Simard, R., Chen, E. J., & Kelton, W. D. 2002, Operations Res., 50, 1073 [CrossRef] [Google Scholar]
 Weltevrede, P., Abdo, A. A., Ackermann, M., et al. 2010, ApJ, 708, 1426 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Figure 1: The dependence of H and M on the approximate skymap significance X for ray pulse phases above 500 MeV from the Vela pulsar as described in the text. The errors represent standard deviations based on 20 independent FermiLAT Vela pulsar data sets, each of length up to 4 days. The horisontal line represent the upper limit of M=20, whereas the dashed line represent a fit of the form H=1.9X^{2} for X>3. 

Open with DEXTER  
In the text 
Figure 2: The distribution of the Hstatistic derived from Monte Carlo simulations with the best fit model for the cumulative probability (solid line) and the old version of the probability distribution given by DSR shown as a dashed line. 

Open with DEXTER  
In the text 
Copyright ESO 2010