circular_linear_correlation

circular_linear_correlation.py

def circular_linear_correlation(alpha, x):
    # Authors:  Jean-Remi King <jeanremi.king@gmail.com>
    #           Niccolo Pescetelli <niccolo.pescetelli@gmail.com>
    #
    # Licence : BSD-simplified
    """

    Parameters
    ----------
        alpha : numpy.array, shape (n_angles, n_dims)
            The angular data (if n_dims == 1, repeated across all x dimensions)
        x : numpy.array, shape (n_angles, n_dims)
            The linear data
    Returns
    -------
        R : numpy.array, shape (n_dims)
            R values
        R2 : numpy.array, shape (n_dims)
            R square values
        p_val : numpy.array, shape (n_dims)
            P values

    Adapted from:
        Circular Statistics Toolbox for Matlab
        By Philipp Berens, 2009
        berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
        Equantion 27.47
    """

    from scipy.stats import chi2
    import numpy as np

    def corr(X, Y):
        if X.ndim == 1:
            X = X[:, None]
        if Y.ndim == 1:
            Y = Y[:, None]
        if Y.shape != X.shape:
            raise ValueError('X and Y must have identical shapes.')
        coef = np.nan * np.zeros(X.shape[1])
        for idx, (x, y) in enumerate(zip(X.T, Y.T)):
            coef[idx] = np.corrcoef(x, y)[0, 1]
        return coef

    # computes correlation for sin and cos separately
    alpha = alpha % (2 * np.pi) - np.pi
    rxs = corr(x, np.sin(alpha))
    rxc = corr(x, np.cos(alpha))
    rcs = corr(np.sin(alpha), np.cos(alpha))

    # Adapted from equation 27.47
    R = (rxc ** 2 + rxs ** 2 - 2 * rxc * rxs * rcs) / (1 - rcs ** 2)
    R2 = np.sqrt(R ** 2)

    # Get degrees of freedom
    n = len(alpha)
    pval = 1 - chi2.cdf(n * R2, 2)

    return R, R2, pval