PulsarSearch

Log survival function of the chi-squared distribution.

Examples

julia> using PulsarSearch;

julia> using Distributions;

julia> chi2 = 31;

julia> d = Chisq(2);

julia> isapprox(chi2_logp(chi2, 2), logccdf(d, chi2), atol=0.1)
true

julia> chi2 = Array([5, 32]);

julia> all(isapprox.(chi2_logp.(chi2, 2), logccdf.(d, chi2), atol=0.1))
true

Number of Gaussian sigmas corresponding to tail probability.

This function computes the value of the characteristic function of a standard Gaussian distribution for the tail probability equivalent to the provided p-value, and turns this value into units of standard deviations away from the Gaussian mean. This allows the user to make a statement about the signal such as “I detected this pulsation at 4.1 sigma

Number of Gaussian sigmas corresponding to tail log-probability.

This function computes the value of the characteristic function of a standard Gaussian distribution for the tail probability equivalent to the provided p-value, and turns this value into units of standard deviations away from the Gaussian mean. This allows the user to make a statement about the signal such as “I detected this pulsation at 4.1 sigma

The example values below are obtained by brute-force integrating the Gaussian probability density function using the mpmath library between Nsigma and +inf.

Examples

julia> using PulsarSearch

julia> pvalues = Array([0.15865525393145707, 0.0013498980316301035, 9.865877e-10, 6.22096e-16, 3.0567e-138]);

julia> log_pvalues = log.(pvalues);

julia> sigmas = Array([1, 3, 6, 8, 25]);

julia> all(isapprox(equivalent_gaussian_Nsigma_from_logp.(log_pvalues), sigmas; atol=0.1))
true

Equivalent gaussian sigma for small log-probability.

Return the equivalent gaussian sigma corresponding to the natural log of the cumulative gaussian probability logp. In other words, return x, such that Q(x) = p, where Q(x) is the cumulative normal distribution. This version uses the rational approximation from Abramowitz and Stegun, eqn 26.2.23, that claims to be precise to ~1e-4. Using the log(P) as input gives a much extended range.

The parameters here are the result of a best-fit, with no physical meaning.

Translated from Scott Ransom's PRESTO

Natural log of the Gamma function in its asymptotic limit.

Return the natural log of the gamma function in its asymptotic limit as z->infty. This is from Abramowitz and Stegun eqn 6.1.41.

Translated from Scott Ransom's PRESTO

Asymptotic natural log of incomplete gamma function.

Return the natural log of the incomplete gamma function in its asymptotic limit as z->infty. This is from Abramowitz and Stegun eqn 6.5.32.

Translated from Scott Ransom's PRESTO

Examples

julia> using PulsarSearch

julia> pvalues = Array([0.15865525393145707, 0.0013498980316301035]);

julia> sigmas = Array([1, 3]);

julia> all(isapprox(equivalent_gaussian_Nsigma.(pvalues), sigmas; atol=0.1))
true

Calculate a multi-trial p-value from the log of a single-trial one.

This allows to work around Numba's limitation on longdoubles, a way to vectorize the computation when we need longdouble precision.

Parameters

logp1 : float The natural logarithm of the significance at which we reject the null hypothesis on each single trial. n : int The number of trials

Returns

logpn : float The log of the significance at which we reject the null hypothesis after multiple trials

Calculate a multi-trial p-value from the log of a single-trial one.

This allows to work around Numba's limitation on longdoubles, a way to vectorize the computation when we need longdouble precision.

Parameters

logpn : float The natural logarithm of the significance at which we want to reject the null hypothesis after multiple trials n : int The number of trials

Returns

logp1 : float The log of the significance at which we reject the null hypothesis on each single trial.

Return the detection level for the $Z^2_n$ statistics.

See Buccheri et al. (1983), Bendat and Piersol (1971).

Parameters

n : int, default 2 The $n$ in $Z^2_n$ (number of harmonics, including the fundamental) epsilon : float, default 0.01 The fractional probability that the signal has been produced by noise

Other Parameters

ntrial : int The number of trials executed to find this profile nsummedspectra : int Number of Z_2^n periodograms that are being averaged

Returns

detlev : float The epoch folding statistics corresponding to a probability epsilon * 100 % that the signal has been produced by noise

Calculate the probability of a certain folded profile, due to noise.

Parameters

z2 : float A $Z^2_n$ statistics value n : int, default 2 The $n$ in $Z^2_n$ (number of harmonics, including the fundamental)

Other Parameters

ntrial : int The number of trials executed to find this profile nsummedspectra : int Number of $Z_2^n$ periodograms that were averaged to obtain z2

Returns

p : float The probability that the $Z^2_n$ value has been produced by noise

Calculate the probability of a certain folded profile, due to noise.

Parameters

z2 : float A $Z^2_n$ statistics value n : int, default 2 The $n$ in $Z^2_n$ (number of harmonics, including the fundamental)

Other Parameters

ntrial : int The number of trials executed to find this profile nsummedspectra : int Number of $Z_2^n$ periodograms that were averaged to obtain z2

Returns

p : float The probability that the $Z^2_n$ value has been produced by noise

z_n(phases, n) --> zsq

$Z^2_n$ statistics, à la Buccheri+83, A&A, 128, 245, eq. 2.

Parameters

phase : array of floats The phases of the events

n : int, default 2 Number of harmonics, including the fundamental

Returns

zsq : float The $Z^2_n$ statistic

Examples

julia> using PulsarSearch

julia> z_n([10., 0., 0., 0., 0.], 2)
20.0
julia> z_n(ones(10), 2)
40.0
julia> z_n(Array([0.5]), 2)
0.0

z_n_binned(profile, n) --> zsq

$Z^2_n$ statistic for pulse profiles from binned events

See Bachetti+2021, arXiv:2012.11397

Parameters

profile : array of floats The folded pulse profile (containing the number of photons falling in each pulse bin)

n : int Number of harmonics, including the fundamental

Returns

zsq : float The value of the statistic

Examples

julia> using PulsarSearch

julia> z_n_binned(zeros(10), 2)
0.0

julia> z_n_binned(ones(10), 2) < 1e-30
true

julia> z_n_binned([10., 0., 0., 0., 0.], 2)
40.0

z_n_search(times, n, fmin, fmax [,oversample]) --> freqs, zsq_stat

Calculate the $Z^2_n$ statistics at trial frequencies in photon data.

Parameters

times : array-like the event arrival times

n : int the number of harmonics in $Z^2_n$

Other Parameters

fmin : float minimum pulse frequency to search

fmax : float maximum pulse frequency to search

oversample : float Oversampling factor with respect to the usual 1/T/n rule

Returns

fgrid : array-like frequency grid of the epoch folding periodogram

zsq_stat : array-like the Z^2_n statistics corresponding to each frequency bin.

z_n_search_hist(times, n, fmin, fmax [,oversample]) --> freqs, zsq_stat

Calculate the $Z^2_n$ statistics at trial frequencies in photon data. Pre-bins the data using a histogram. At the moment, it is not faster. It will be after I add the 2-d hist

• shift-and-add. But it allowed to test that z_n_binned works with

binned data.

Parameters

times : array-like the event arrival times

n : int the number of harmonics in $Z^2_n$

Other Parameters

fmin : float minimum pulse frequency to search

fmax : float maximum pulse frequency to search

oversample : float Oversampling factor with respect to the usual 1/T/n rule

Returns

fgrid : array-like frequency grid of the epoch folding periodogram

zsq_stat : array-like the Z^2_n statistics corresponding to each frequency bin.