Repository: StingraySoftware/stingray

Issue:
Calculate signal upper limits from power spectra
It would be useful if the various PowerSpectrum objects could identify the highest peak in the spectrum (within some frequency range) and calculate detection thresholds, as well as signal upper limits (e.g. following [Vaughan et al. 1994](https://ui.adsabs.harvard.edu/abs/1994ApJ...435..362V))

Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9460c449..371f4488 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,15 +1,8 @@
 Changelog
 =========
 
-Not Yet Released
-----------------
-
-`Full list of changes`__
-
-__ https://github.com/StingraySoftware/stingray/compare/v0.3...main
-
-v0.3 (2021-05-31)
------------------
+Unreleased (2021-01-18)
+-----------------------
 
 - Lots of performance improvements
 - Faster simulations
@@ -28,7 +21,7 @@ v0.3 (2021-05-31)
 
 `Full list of changes`__
 
-__ https://github.com/StingraySoftware/stingray/compare/v0.2...v0.3
+__ https://github.com/StingraySoftware/stingray/compare/v0.2...main
 
 v0.2 (2020-06-17)
 -----------------
@@ -47,7 +40,7 @@ v0.2 (2020-06-17)
 
 `Full list of changes`__
 
-__ https://github.com/StingraySoftware/stingray/compare/v0.1.3...v0.2
+__ https://github.com/StingraySoftware/stingray/compare/v0.1.3...v0.2>
 
 v0.1.3 (2019-06-11)
 -------------------
diff --git a/stingray/io.py b/stingray/io.py
index 3bf2785e..6da03caa 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -584,6 +584,7 @@ def load_events_and_gtis(
     if mission_key not in probe_header:
         mission_key = "TELESCOP"
     mission = probe_header[mission_key].lower()
+
     db = read_mission_info(mission)
     instkey = get_key_from_mission_info(db, "instkey", "INSTRUME")
     instr = mode = None
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 18e66df1..00744250 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -10,7 +10,7 @@ import stingray.utils as utils
 from stingray.crossspectrum import AveragedCrossspectrum, Crossspectrum
 from stingray.gti import bin_intervals_from_gtis, check_gtis
 from stingray.largememory import createChunkedSpectra, saveData
-from stingray.stats import pds_probability
+from stingray.stats import pds_probability, amplitude_upper_limit
 from stingray.utils import genDataPath
 
 from .events import EventList
@@ -284,6 +284,82 @@ class Powerspectrum(Crossspectrum):
 
         return pvals
 
+    def modulation_upper_limit(self, fmin=None, fmax=None, c=0.95):
+        """Upper limit on a sinusoidal modulation.
+
+        To understand the meaning of this amplitude: if the modulation is described by:
+
+        ..math:: p = \overline{p} (1 + a * \sin(x))
+
+        this function returns a.
+
+        If it is a sum of sinusoidal harmonics instead
+        ..math:: p = \overline{p} (1 + \sum_l a_l * \sin(lx))
+        a is equivalent to :math:`\sqrt(\sum_l a_l^2)`.
+
+        See `stingray.stats.power_upper_limit`, `stingray.stats.amplitude_upper_limit`
+        for more information.
+
+        Parameters
+        ----------
+        fmin: float
+            The minimum frequency to search (defaults to the first nonzero bin)
+
+        fmax: float
+            The maximum frequency to search (defaults to the Nyquist frequency)
+
+        Other Parameters
+        ----------------
+        c: float
+            The confidence value for the upper limit (e.g. 0.95 = 95%)
+
+        Results
+        -------
+        a: float
+            The modulation amplitude that could produce P>pmeas with 1 - c probability
+
+        Examples
+        --------
+        >>> pds = Powerspectrum()
+        >>> pds.norm = "leahy"
+        >>> pds.freq = np.arange(0., 5.)
+        >>> # Note: this pds has 40 as maximum value between 2 and 5 Hz
+        >>> pds.power = np.array([100000, 1, 1, 40, 1])
+        >>> pds.m = 1
+        >>> pds.nphots = 30000
+        >>> pds.modulation_upper_limit(fmin=2, fmax=5, c=0.99)
+        0.1016...
+        >>> pds.norm = "frac"
+        >>> pds.modulation_upper_limit(2, 5, 0.99)
+        Traceback (most recent call last):
+         ...
+        ValueError: Modulation upper limit is only available in the Leahy...
+        """
+
+        # TODO: add case with norm different from Leahy
+        if self.norm != 'leahy':
+            raise ValueError("Modulation upper limit is only available in the Leahy normalization")
+
+        freq = self.freq
+        fnyq = np.max(freq)
+        power = self.power
+        freq_mask = freq > 0
+        if fmin is not None or fmax is not None:
+            if fmin is not None:
+                freq_mask[freq < fmin] = 0
+            if fmax is not None:
+                freq_mask[freq > fmax] = 0
+        freq = freq[freq_mask]
+        power = power[freq_mask]
+
+        maximum_val = np.argmax(power)
+        nyq_ratio = freq[maximum_val] / fnyq
+
+        # I multiply by M because the formulas from Vaughan+94 treat summed powers, while here
+        # we have averaged powers.
+        return amplitude_upper_limit(
+            power[maximum_val] * self.m, self.nphots, n=self.m, c=c, nyq_ratio=nyq_ratio, fft_corr=True)
+
 
 class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     """
diff --git a/stingray/stats.py b/stingray/stats.py
index ad349e7c..fb867d31 100644
--- a/stingray/stats.py
+++ b/stingray/stats.py
@@ -20,7 +20,17 @@ __all__ = ['p_multitrial_from_single_trial',
            'classical_pvalue',
            'chi2_logp',
            'equivalent_gaussian_Nsigma',
-           'equivalent_gaussian_Nsigma_from_logp']
+           'equivalent_gaussian_Nsigma_from_logp',
+           'power_confidence_limits',
+           'power_upper_limit',
+           'pf_from_ssig',
+           'pf_from_a',
+           'pf_upper_limit',
+           'a_from_pf',
+           'a_from_ssig',
+           'ssig_from_a',
+           'ssig_from_pf',
+           'amplitude_upper_limit']
 
 
 @vectorize([float64(float32),
@@ -776,3 +786,335 @@ def _pavnosigfun(power, nspec):
         m -= 1
 
     return sum
+
+
+def power_confidence_limits(preal, n=1, c=0.95):
+    """Confidence limits on power, given a (theoretical) signal power.
+
+    This is to be used when we *expect* a given power (e.g. from the pulsed
+    fraction measured in previous observations) and we want to know the
+    range of values the measured power could take to a given confidence level.
+    Adapted from Vaughan et al. 1994, noting that, after appropriate
+    normalization of the spectral stats, the distribution of powers in the PDS
+    and the Z^2_n searches is always described by a noncentral chi squared
+    distribution.
+
+    Parameters
+    ----------
+    preal: float
+        The theoretical signal-generated value of power
+
+    Other Parameters
+    ----------------
+    n: int
+        The number of summed powers to obtain the result. It can be multiple
+        harmonics of the PDS, adjacent bins in a PDS summed to collect all the
+        power in a QPO, or the n in Z^2_n
+    c: float
+        The confidence level (e.g. 0.95=95%)
+
+    Results
+    -------
+    pmeas: [float, float]
+        The upper and lower confidence interval (a, 1-a) on the measured power
+
+    Examples
+    --------
+    >>> cl = power_confidence_limits(150, c=0.84)
+    >>> np.allclose(cl, [127, 176], atol=1)
+    True
+    """
+    rv = stats.ncx2(2 * n, preal)
+    return rv.ppf([1 - c, c])
+
+
+def power_upper_limit(pmeas, n=1, c=0.95):
+    """Upper limit on signal power, given a measured power in the PDS/Z search.
+
+    Adapted from Vaughan et al. 1994, noting that, after appropriate
+    normalization of the spectral stats, the distribution of powers in the PDS
+    and the Z^2_n searches is always described by a noncentral chi squared
+    distribution.
+
+    Note that Vaughan+94 gives p(pmeas | preal), while we are interested in
+    p(real | pmeas), which is not described by the NCX2 stat. Rather than
+    integrating the CDF of this probability distribution, we start from a
+    reasonable approximation and fit to find the preal that gives pmeas as
+    a (e.g.95%) confidence limit.
+
+    As Vaughan+94 shows, this power is always larger than the observed one.
+    This is because we are looking for the maximum signal power that,
+    combined with noise powers, would give the observed power. This involves
+    the possibility that noise powers partially cancel out some signal power.
+
+    Parameters
+    ----------
+    pmeas: float
+        The measured value of power
+
+    Other Parameters
+    ----------------
+    n: int
+        The number of summed powers to obtain pmeas. It can be multiple
+        harmonics of the PDS, adjacent bins in a PDS summed to collect all the
+        power in a QPO, or the n in Z^2_n
+    c: float
+        The confidence value for the probability (e.g. 0.95 = 95%)
+
+    Results
+    -------
+    psig: float
+        The signal power that could produce P>pmeas with 1 - c probability
+
+    Examples
+    --------
+    >>> pup = power_upper_limit(40, 1, 0.99)
+    >>> np.isclose(pup, 75, atol=2)
+    True
+    """
+    def ppf(x):
+        rv = stats.ncx2(2 * n, x)
+        return rv.ppf(1-c)
+
+    def isf(x):
+        rv = stats.ncx2(2 * n, x)
+        return rv.ppf(c)
+
+    def func_to_minimize(x, xmeas):
+        return np.abs(ppf(x) - xmeas)
+
+    from scipy.optimize import minimize
+    initial = isf(pmeas)
+
+    res = minimize(func_to_minimize, [initial], pmeas, bounds=[(0, initial * 2)])
+
+    return res.x[0]
+
+
+def amplitude_upper_limit(pmeas, counts, n=1, c=0.95, fft_corr=False, nyq_ratio=0):
+    """Upper limit on a sinusoidal modulation, given a measured power in the PDS/Z search.
+
+    Eq. 10 in Vaughan+94 and `a_from_ssig`: they are equivalent but Vaughan+94
+    corrects further for the response inside an FFT bin and at frequencies close
+    to Nyquist. These two corrections are added by using fft_corr=True and
+    nyq_ratio to the correct :math:`f / f_{Nyq}` of the FFT peak
+
+    To understand the meaning of this amplitude: if the modulation is described by:
+
+    ..math:: p = \overline{p} (1 + a * \sin(x))
+
+    this function returns a.
+
+    If it is a sum of sinusoidal harmonics instead
+    ..math:: p = \overline{p} (1 + \sum_l a_l * \sin(lx))
+    a is equivalent to :math:`\sqrt(\sum_l a_l^2)`.
+
+    See `power_upper_limit`
+
+    Parameters
+    ----------
+    pmeas: float
+        The measured value of power
+
+    counts: int
+        The number of counts in the light curve used to calculate the spectrum
+
+    Other Parameters
+    ----------------
+    n: int
+        The number of summed powers to obtain pmeas. It can be multiple
+        harmonics of the PDS, adjacent bins in a PDS summed to collect all the
+        power in a QPO, or the n in Z^2_n
+    c: float
+        The confidence value for the probability (e.g. 0.95 = 95%)
+    fft_corr: bool
+        Apply a correction for the expected power concentrated in an FFT bin,
+        which is about 0.773 on average (it's 1 at the center of the bin, 2/pi
+        at the bin edge.
+    nyq_ratio: float
+        Ratio of the frequency of this feature with respect to the Nyquist
+        frequency. Important to know when dealing with FFTs, because the FFT
+        response decays between 0 and f_Nyq similarly to the response inside
+        a frequency bin: from 1 at 0 Hz to ~2/pi at f_Nyq
+
+    Results
+    -------
+    a: float
+        The modulation amplitude that could produce P>pmeas with 1 - c probability
+
+    Examples
+    --------
+    >>> aup = amplitude_upper_limit(40, 30000, 1, 0.99)
+    >>> aup_nyq = amplitude_upper_limit(40, 30000, 1, 0.99, nyq_ratio=1)
+    >>> np.isclose(aup_nyq, aup / (2 / np.pi))
+    True
+    >>> aup_corr = amplitude_upper_limit(40, 30000, 1, 0.99, fft_corr=True)
+    >>> np.isclose(aup_corr, aup / np.sqrt(0.773))
+    True
+    """
+
+    uplim = power_upper_limit(pmeas, n, c)
+    a = a_from_ssig(uplim, counts)
+    if fft_corr:
+        factor = 1 / np.sqrt(0.773)
+        a *= factor
+    if nyq_ratio > 0:
+        factor = np.pi / 2 * nyq_ratio
+        sinc_factor = np.sin(factor) / factor
+        a /= sinc_factor
+    return a
+
+
+def pf_upper_limit(*args, **kwargs):
+    """Upper limit on pulsed fraction, given a measured power in the PDS/Z search.
+
+    See `power_upper_limit` and `pf_from_ssig`.
+    All arguments are the same as `amplitude_upper_limit`
+
+    Parameters
+    ----------
+    pmeas: float
+        The measured value of power
+
+    counts: int
+        The number of counts in the light curve used to calculate the spectrum
+
+    Other Parameters
+    ----------------
+    n: int
+        The number of summed powers to obtain pmeas. It can be multiple
+        harmonics of the PDS, adjacent bins in a PDS summed to collect all the
+        power in a QPO, or the n in Z^2_n
+    c: float
+        The confidence value for the probability (e.g. 0.95 = 95%)
+    fft_corr: bool
+        Apply a correction for the expected power concentrated in an FFT bin,
+        which is about 0.773 on average (it's 1 at the center of the bin, 2/pi
+        at the bin edge.
+    nyq_ratio: float
+        Ratio of the frequency of this feature with respect to the Nyquist
+        frequency. Important to know when dealing with FFTs, because the FFT
+        response decays between 0 and f_Nyq similarly to the response inside
+        a frequency bin: from 1 at 0 Hz to ~2/pi at f_Nyq
+
+    Results
+    -------
+    pf: float
+        The pulsed fraction that could produce P>pmeas with 1 - c probability
+
+    Examples
+    --------
+    >>> pfup = pf_upper_limit(40, 30000, 1, 0.99)
+    >>> np.isclose(pfup, 0.13, atol=0.01)
+    True
+    """
+
+    return pf_from_a(amplitude_upper_limit(*args, **kwargs))
+
+
+def pf_from_a(a):
+    """Pulsed fraction from fractional amplitude of modulation.
+
+    If the pulsed profile is defined as
+    p = mean * (1 + a * sin(phase)),
+
+    we define "pulsed fraction" as 2a/b, where b = mean + a is the maximum and
+    a is the amplitude of the modulation.
+
+    Hence, pulsed fraction = 2a/(1+a)
+
+    Examples
+    --------
+    >>> pf_from_a(1)
+    1.0
+    >>> pf_from_a(0)
+    0.0
+    """
+    return 2 * a / (1 + a)
+
+
+def a_from_pf(p):
+    """Fractional amplitude of modulation from pulsed fraction
+
+    If the pulsed profile is defined as
+    p = mean * (1 + a * sin(phase)),
+
+    we define "pulsed fraction" as 2a/b, where b = mean + a is the maximum and
+    a is the amplitude of the modulation.
+
+    Hence, a = pf / (2 - pf)
+
+    Examples
+    --------
+    >>> a_from_pf(1)
+    1.0
+    >>> a_from_pf(0)
+    0.0
+    """
+    return p / (2 - p)
+
+
+def ssig_from_a(a, ncounts):
+    """Theoretical power in the Z or PDS search for a sinusoid of amplitude a.
+
+    From Leahy et al. 1983, given a pulse profile
+    p = lambda * (1 + a * sin(phase)),
+    The theoretical value of Z^2_n is Ncounts / 2 * a^2
+
+    Note that if there are multiple sinusoidal components, one can use
+    a = sqrt(sum(a_l))
+    (Bachetti+2021b)
+
+    Examples
+    --------
+    >>> round(ssig_from_a(0.1, 30000), 1)
+    150.0
+    """
+    return ncounts / 2 * a ** 2
+
+
+def a_from_ssig(ssig, ncounts):
+    """Amplitude of a sinusoid corresponding to a given Z/PDS value
+
+    From Leahy et al. 1983, given a pulse profile
+    p = lambda * (1 + a * sin(phase)),
+    The theoretical value of Z^2_n is Ncounts / 2 * a^2
+
+    Note that if there are multiple sinusoidal components, one can use
+    a = sqrt(sum(a_l))
+    (Bachetti+2021b)
+
+    Examples
+    --------
+    >>> a_from_ssig(150, 30000)
+    0.1
+    """
+    return np.sqrt(2 * ssig / ncounts)
+
+
+def ssig_from_pf(pf, ncounts):
+    """Theoretical power in the Z or PDS for a sinusoid of pulsed fraction pf.
+
+    See `ssig_from_a` and `a_from_pf` for more details
+
+    Examples
+    --------
+    >>> round(ssig_from_pf(pf_from_a(0.1), 30000), 1)
+    150.0
+    """
+    a = a_from_pf(pf)
+    return ncounts / 2 * a ** 2
+
+
+def pf_from_ssig(ssig, ncounts):
+    """Estimate pulsed fraction for a sinusoid from a given Z or PDS power.
+
+    See `a_from_ssig` and `pf_from_a` for more details
+
+    Examples
+    --------
+    >>> round(a_from_pf(pf_from_ssig(150, 30000)), 1)
+    0.1
+    """
+    a = a_from_ssig(ssig, ncounts)
+    return pf_from_a(a)


Test patch:
diff --git a/stingray/tests/test_stats.py b/stingray/tests/test_stats.py
index a968a744..d6a2de5b 100644
--- a/stingray/tests/test_stats.py
+++ b/stingray/tests/test_stats.py
@@ -212,3 +212,7 @@ class TestClassicalSignificances(object):
         logp = z2_n_logprobability(stat, n=n)
         p = z2_n_probability(stat, n=n)
         assert np.allclose(logp, np.log(p))
+
+    def test_power_upper_limit(self):
+        # Use example from Vaughan+94
+        assert np.isclose(power_upper_limit(40, 1, 0.99), 75, rtol=0.1)