[Feature] fractal_tmf (Multifractal Nonlinearity)

neurokit2/complexity/__init__.py

@@ -52,6 +52,7 @@
 from .fractal_psdslope import fractal_psdslope
 from .fractal_sda import fractal_sda
 from .fractal_sevcik import fractal_sevcik
+from .fractal_tmf import fractal_tmf
 from .information_fisher import fisher_information
 from .information_fishershannon import fishershannon_information
 from .information_gain import information_gain
@@ -200,6 +201,7 @@
     "fractal_sevcik",
     "fractal_mandelbrot",
     "fractal_mfdfa",
+    "fractal_tmf",
     "fractal_nld",
     "fractal_psdslope",
     "fractal_sda",

neurokit2/complexity/complexity_hjorth.py

@@ -11,14 +11,21 @@ def complexity_hjorth(signal):
     NeuroKit returns complexity directly in the output tuple, but the other parameters can be found
     in the dictionary.
 
-
     * The **activity** parameter is simply the variance of the signal, which corresponds to the
       mean power of a signal (if its mean is 0).
 
       .. math::
 
         Activity = \\sigma_{signal}^2
 
+    * The **mobility** parameter represents the mean frequency or the proportion of standard
+      deviation of the power spectrum. This is defined as the square root of variance of the
+      first derivative of the signal divided by the variance of the signal.
+
+      .. math::
+
+        Mobility = \\frac{\\sigma_{dd}/ \\sigma_{d}}{Complexity}
+
     * The **complexity** parameter gives an estimate of the bandwidth of the signal, which
       indicates the similarity of the shape of the signal to a pure sine wave (for which the
       value converges to 1). In other words, it is a measure of the "excessive details" with
@@ -29,14 +36,6 @@ def complexity_hjorth(signal):
 
         Complexity = \\sigma_{d}/ \\sigma_{signal}
 
-    * The **mobility** parameter represents the mean frequency or the proportion of standard
-      deviation of the power spectrum. This is defined as the square root of variance of the
-      first derivative of the signal divided by the variance of the signal.
-
-      .. math::
-
-        Mobility = \\frac{\\sigma_{dd}/ \\sigma_{d}}{Complexity}
-
     :math:`d` and :math:`dd` represent the first and second derivatives of the signal, respectively.
 
     Hjorth (1970) illustrated the parameters as follows:

neurokit2/complexity/entropy_range.py

@@ -3,19 +3,20 @@
 
 
 def entropy_range(signal, dimension=3, delay=1, tolerance="sd", approximate=False, **kwargs):
-    """**Range Entropy (RangeEn)**
+    """**Range Entropy (RangEn)**
 
-    Introduced by Omidvarnia et al. (2018), RangeEn refers to a modified form of SampEn (or ApEn).
+    Introduced by Omidvarnia et al. (2018), Range Entropy (RangEn or RangeEn) refers to a modified
+    form of SampEn (or ApEn).
 
     Both ApEn and SampEn compute the logarithmic likelihood that runs of patterns that are close
     remain close on the next incremental comparisons, of which this closeness is estimated by the
     Chebyshev distance. Range Entropy uses instead a normalized "range distance", resulting in
-    modified forms of ApEn and SampEn, **RangeEn (A)** (*mApEn*) and **RangeEn (B)** (*mSampEn*).
+    modified forms of ApEn and SampEn, **RangEn (A)** (*mApEn*) and **RangEn (B)** (*mSampEn*).
 
-    However, the RangeEn (A), based on ApEn, often yields undefined entropies (i.e., *NaN* or
-    *Inf*). As such, using RangeEn (B) is recommended instead.
+    However, the RangEn (A), based on ApEn, often yields undefined entropies (i.e., *NaN* or
+    *Inf*). As such, using RangEn (B) is recommended instead.
 
-    RangeEn is described as more robust to nonstationary signal changes, and has a more linear
+    RangEn is described as more robust to nonstationary signal changes, and has a more linear
     relationship with the Hurst exponent (compared to ApEn and SampEn), and has no need for signal
     amplitude correction.
 
@@ -39,8 +40,8 @@ def entropy_range(signal, dimension=3, delay=1, tolerance="sd", approximate=Fals
         :func:`complexity_tolerance` to estimate the optimal value for this parameter.
     approximate : bool
         The entropy algorithm to use. If ``False`` (default), will use sample entropy and return
-        *mSampEn* (**RangeEn B**). If ``True``, will use approximate entropy and return *mApEn*
-        (**RangeEn A**).
+        *mSampEn* (**RangEn B**). If ``True``, will use approximate entropy and return *mApEn*
+        (**RangEn A**).
     **kwargs
         Other arguments.
 
@@ -50,7 +51,7 @@ def entropy_range(signal, dimension=3, delay=1, tolerance="sd", approximate=Fals
 
     Returns
     -------
-    RangeEn : float
+    RangEn : float
         Range Entropy. If undefined conditional probabilities are detected (logarithm
         of sum of conditional probabilities is ``ln(0)``), ``np.inf`` will
         be returned, meaning it fails to retrieve 'accurate' regularity information.
@@ -68,16 +69,16 @@ def entropy_range(signal, dimension=3, delay=1, tolerance="sd", approximate=Fals
       signal = nk.signal_simulate(duration=2, sampling_rate=100, frequency=[5, 6])
 
       # Range Entropy B (mSampEn)
-      RangeEnB, info = nk.entropy_range(signal, approximate=False)
-      RangeEnB
+      RangEnB, info = nk.entropy_range(signal, approximate=False)
+      RangEnB
 
       # Range Entropy A (mApEn)
-      RangeEnA, info = nk.entropy_range(signal, approximate=True)
-      RangeEnA
+      RangEnA, info = nk.entropy_range(signal, approximate=True)
+      RangEnA
 
       # Range Entropy A (corrected)
-      RangeEnAc, info = nk.entropy_range(signal, approximate=True, corrected=True)
-      RangeEnAc
+      RangEnAc, info = nk.entropy_range(signal, approximate=True, corrected=True)
+      RangEnAc
 
     References
     ----------

neurokit2/complexity/fractal_dfa.py

@@ -128,7 +128,7 @@ def fractal_dfa(
 
     See Also
     --------
-    fractal_hurst
+    fractal_hurst, fractal_tmf
 
     Examples
     ----------

neurokit2/complexity/fractal_tmf.py

@@ -0,0 +1,125 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import scipy.stats
+
+from ..signal import signal_surrogate
+from .fractal_dfa import fractal_dfa
+
+
+def fractal_tmf(signal, n=40, show=False, **kwargs):
+    """**Multifractal Nonlinearity (tMF)**
+
+    The Multifractal Nonlinearity index (*t*\\MF) is the *t*\\-value resulting from the comparison
+    between the multifractality of the signal (measured by the spectrum width, see
+    :func:`.fractal_dfa`) and the multifractality of linearized
+    :func:`surrogates <.signal_surrogate>` obtained by the IAAFT method (i.e., reshuffled series
+    with comparable linear structure).
+
+    This statistics grows larger the more the original series departs from the multifractality
+    attributable to the linear structure of IAAFT surrogates. When p-value reaches significance, we
+    can conclude that the signal's multifractality encodes processes that a linear contingency
+    cannot.
+
+    This index provides an extension of the assessment of multifractality, of which the
+    multifractal spectrum is by itself a measure of heterogeneity, rather than interactivity.
+    As such, it cannot alone be used to assess the specific presence of cascade-like interactivity
+    in the time series, but must be compared to the spectrum of a sample of its surrogates.
+
+    .. figure:: ../img/bell2019.png
+       :alt: Figure from Bell et al. (2019).
+       :target: https://doi.org/10.3389/fphys.2019.00998
+
+    Both significantly negative and positive values can indicate interactivity, as any difference
+    from the linear structure represented by the surrogates is an indication of nonlinear
+    contingence. Indeed, if the degree of heterogeneity for the original series is significantly
+    less than for the sample of linear surrogates, that is no less evidence of a failure of
+    linearity than if the degree of heterogeneity is significantly greater.
+
+
+    .. note::
+
+        Help us review the implementation of this index by checking-it out and letting us know
+        wether it is correct or not.
+
+    Parameters
+    ----------
+    signal : Union[list, np.array, pd.Series]
+        The signal (i.e., a time series) in the form of a vector of values.
+    n : int
+        Number of surrogates. The literature uses values between 30 and 40.
+    **kwargs : optional
+        Other arguments to be passed to :func:`.fractal_dfa`.
+
+    Returns
+    -------
+    float
+        tMF index.
+    info : dict
+        A dictionary containing additional information, such as the p-value.
+
+    See Also
+    --------
+    fractal_dfa, .signal_surrogate
+
+    Examples
+    ----------
+    .. ipython:: python
+
+      import neurokit2 as nk
+
+      # Simulate a Signal
+      signal = nk.signal_simulate(duration=1, sampling_rate=200, frequency=[5, 6, 12], noise=0.2)
+
+      # Compute tMF
+      @savefig p_fractal_tmf.png scale=100%
+      tMF, info = nk.fractal_tmf(signal, n=100, show=True)
+      @suppress
+      plt.close()
+
+    .. ipython:: python
+
+      tMF  # t-value
+      info["p"]  # p-value
+
+
+    References
+    ----------
+    * Ihlen, E. A., & Vereijken, B. (2013). Multifractal formalisms of human behavior. Human
+      movement science, 32(4), 633-651.
+    * Kelty-Stephen, D. G., Palatinus, K., Saltzman, E., & Dixon, J. A. (2013). A tutorial on
+      multifractality, cascades, and interactivity for empirical time series in ecological science.
+      Ecological Psychology, 25(1), 1-62.
+    * Bell, C. A., Carver, N. S., Zbaracki, J. A., & Kelty-Stephen, D. G. (2019). Non-linear
+      amplification of variability through interaction across scales supports greater accuracy in
+      manual aiming: evidence from a multifractal analysis with comparisons to linear surrogates in
+      the fitts task. Frontiers in physiology, 10, 998.
+
+    """
+    # Sanity checks
+    if isinstance(signal, (np.ndarray, pd.DataFrame)) and signal.ndim > 1:
+        raise ValueError(
+            "Multidimensional inputs (e.g., matrices or multichannel data) are not supported yet."
+        )
+
+    info = {}
+
+    w0 = fractal_dfa(signal, multifractal=True, show=False)[0]["Width"]
+
+    w = np.zeros(20)
+    for i in range(20):
+        surro = signal_surrogate(signal, method="IAAFT")
+        w[i] = fractal_dfa(surro, multifractal=True, show=False)[0]["Width"]
+
+    # Run t-test
+    t, p = scipy.stats.ttest_1samp(w, w0)
+    t = t[0]
+    info["p"] = p[0]
+
+    if show is True:
+        pd.Series(w).plot(kind="density", label="Width of surrogates")
+        plt.axvline(x=w0.values, c="red", label="Width of original signal")
+        plt.title(f"tMF = {t:.2f}, p = {info['p']:.2f}")
+        plt.legend()
+
+    return t, info