basic circular stats

circular_data_functions.py

"""
=========================================================
circular data analysis functions
=========================================================

"""

# Authors : Anne Kosem and Alexandre Gramfort
# License : Simplified BSD

from math import sqrt, atan2, exp
import numpy as np


def circular_meanangle(phases):  # circular mean of phase distribution "phases"
    ang = np.mean(np.exp(1j * phases))
    return atan2(ang.imag, ang.real)


def circular_PLV(phases):  # Phase Locking Value (Lachaux et al., 1998)
    return abs(np.sum(np.exp(1j * phases)) / len(phases))


def rayleigh_test(phases):  # test if phase distribution is non-uniform; TRUE if pval<0.05 (Fisher, 1995)
    ph = phases.reshape(-1)
    n = len(ph)
    w = np.ones(n)
    r = np.sum(w * np.exp(1j * ph))
    r = abs(r) / np.sum(w)
    R = n * r
    pval = exp(sqrt(1 + 4 * n + 4 * (n ** 2 - R ** 2)) - (1 + 2 * n))
    return pval


def bootstrap_conf(x, n_permutations=10000, seed=0):  # gives 95% confidence interval [ll; ul] of the mean of phase distribution x (Fisher, 1995)
    n_samples = len(x)
    rng = np.random.RandomState(seed)

    C = np.empty(n_permutations, dtype=x.dtype)
    for ind in range(n_permutations):  # surrogate
        C[ind] = circular_meanangle(x[rng.randint(len(x), size=(n_samples))])  # circular mean for each surrogate

    C = np.sort(C)
    indl = 0.025 * n_permutations
    indu = 0.975 * n_permutations
    ll = C[indl]
    ul = C[indu]
    return C, ll, ul

if __name__ == '__main__':

    import pylab as pl

    fig = pl.figure()
    ax = fig.add_subplot(1, 1, 1, polar=True)

    theta = np.pi / 4.
    kappa = 1.  # concentration parameter, decreasing variance with increasing kappa

    # generate sample set following von mises distribution with mean theta and
    # concentration parameter kappa
    sample = np.random.vonmises(theta, kappa, 200)
    pl.hist(sample, bins=20, range=(-np.pi, np.pi), alpha=0.2, normed=True, color='b')
    pl.yticks([])

    # mean angle
    m = circular_meanangle(sample)
    pl.plot([m, m], [0, 1], 'b', linewidth=3)

    # confidence interval
    _, ll, ul = bootstrap_conf(sample)
    pl.plot([ll, ll], [0, 1], '--b', linewidth=2)
    pl.plot([ul, ul], [0, 1], '--b', linewidth=2)

    # rayleigh test against uniformity, PLV
    pval = rayleigh_test(sample)
    plv = circular_PLV(sample)
    pl.xlabel('probability of distribution uniformity = %G, PLV = %G' % (pval, plv))
    pl.show()