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
    np.random.seed(6314569)
    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
    C, 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))

    from scipy import stats
    import matplotlib.pyplot as plt
    plt.figure()
    mb,sb = C.mean(), C.std()
    p, e, _ = plt.hist(C, bins=50, normed=True, alpha=0.5)
    x = np.linspace(e[0], e[-1], 101)
    plt.plot(x, stats.norm.pdf(x, loc=mb, scale=sb), linewidth=2)

    plt.plot([ll, ll], [0, 1], '--b', linewidth=2)
    plt.plot([ul, ul], [0, 1], '--b', linewidth=2)
    ll2, ul2 = stats.norm.ppf([0.025, 0.975], loc=m, scale=s)
    plt.plot([ll2, ll2], [0, 1], '--g', linewidth=2)
    plt.plot([ul2, ul2], [0, 1], '--g', linewidth=2)

    import circularstats.stats as cs
    mean, ul3, ll3 = cs.circ_mean(sample, return_ci=True, xi=0.05)
    plt.plot([ll3, ll3], [0, 1], '-r', linewidth=2)
    plt.plot([ul3, ul3], [0, 1], '-r', linewidth=2)
    plt.title('Distribution of Circular Mean Estimate')

    print 'bootstrap   ', ll, ul
    print 'boots-normal', ll2, ul2
    print 'asymptotic  ', ll3, ul3
    '''
    bootstrap    0.650457352545 1.03118480889
    boots-normal 0.650198642054 1.03277388054
    asymptotic   [ 0.6549999] [ 1.03193275]
    '''
    import circularstats.api as css
    rp, rs = css.circ_rtest(sample)
    print 'rayleigh test', rp, pval, rp - pval

    pl.show()