HFD Algorithm Appears to be Wrong

Hi @thosvarley,

Actually there is some mistakes pointed out by the developers of the pyrem module (http://gilestrolab.github.io/pyrem/pyrem.univariate.html#pyrem.univariate.hfd). Personnally I had a bad time checking some of the pyEEG functions against the original implementation.

For Higuchi, I'm not even sure of the correction in pyrem because there is something wrong with the indexes. I went for this slightly modified version from pyrem, looking at the original paper cited on pyrem's website and using the appendix of a review about the Higuchi's Fractal Dimension application:

def hfd_pyrem(a, k_max):
    r"""[...] (see pyrem sources)
    """
    L = []
    N = a.size

    # [...]
    for k in range(1, k_max + 1): # k = 1, 2, 3, ..., k_max (k_max depends on N) <-- correction of the indexes
        Lk = 0
        for m in range(1, k + 1): # m = 1, 2, 3, ..., k <-- correction of the indexes
            #we pregenerate all idxs
            idxs = np.arange(1, int(np.floor((N - m) / k) + 1), dtype=np.int32) # i = 1, 2, ..., int((N-m)/k) --> idxs = i

            # Sum of L_m(k)
            # Add (...) - 1 to correct the indexes <-- correction
            Lk += (np.sum(np.abs(a[(m + idxs * k) - 1] - a[(m + k * (idxs - 1)) - 1])) * (N - 1) / (int((N - m) / k)* k)) / k

        L.append(np.log(Lk / (m + 1))) # Unchanged
    # Add ``list(map(lambda x: [np.log(1.0/x), 1], np.arange(1, k_max + 1)))``
    return np.linalg.lstsq(list(map(lambda x: [np.log(1.0/x), 1], np.arange(1, k_max + 1))), L, rcond=None)[0][0]

I know it's an old thread but if it helps, it would makes me happy.

Vincent

forrestbao / pyeeg

HFD Algorithm Appears to be Wrong #27