circstat / pycircstat

Toolbox for circular statistics with Python
MIT License
157 stars 42 forks source link

trouble with correlation #79

Open david-klindt opened 2 years ago

david-klindt commented 2 years ago

Hi guys, thanks for this lovely repo!

I am having trouble with circular correlation between two circular variables. As you can see they are strongly correlated, but the correlation function returns nearly 0:

image

Here are the two signals if you want to try it: timeseries: 1, 2

Is there a problem with the code or am I doing sth wrong?

philippberens commented 2 years ago

It's also not what I would expect - I need to see whether I find time to debug it. Is there a chance you can run the Matlab version if it yields the same?

david-klindt commented 2 years ago

That would be great :) Unfortunately, I don't have a Matlab license

huangziwei commented 2 years ago

the current implementation (Jammalamadaka & SenGupta, 2001) subtracts the mean out of the original data, then compute the correlation, which turns out to be problematic for your data, because ts1 is more or less uniformly distributed, while ts2 is bimodal (more accurately, axial).

here's an alternative implementation (Fisher & Lee, 1983; Zar, 1999) without mean subtraction that can give you a r_{aa}=0.94557:

def aacorr(a: np.ndarray, b: np.ndarray) -> float:

    """Fisher & Lee (1983); Zar (1999) eq (27.43)

    a, b: np.array, (n, )
        circular data in radian

    raa: float
         angular-angular correlation 
    """
    aij = np.triu(a[: ,None] - a).flatten()
    bij = np.triu(b[: ,None] - b).flatten()

    num = np.sum(np.sin(aij) * np.sin(bij))
    den = np.sqrt(np.sum(np.sin(aij) ** 2) * np.sum(np.sin(bij) ** 2))

    raa = num / den

    return raa
Macko94 commented 2 years ago

Hi guys, thanks for this lovely repo!

I am having trouble with circular correlation between two circular variables. As you can see they are strongly correlated, but the correlation function returns nearly 0:

image

Here are the two signals if you want to try it: timeseries: 1, 2

Is there a problem with the code or am I doing sth wrong?

Hi David, Make sure you convert your angular values from degrees to radians. I had the same issue and turned out that the problem was wrong units. Cheers,