gist.py

import numpy as np
from mixtape.vmhmm import inverse_mbessel_ratio


def fit_vonmises(data):
    """
    Parameters
    ----------
    data : np.array, shape=(n_samples, n_features)
    
    Returns
    -------
    means : np.array, shape=(n_features,)
    kappas : np.array, shape=(n_features,)
    """
    n_samples, n_features = data.shape
    
    means = np.arctan2(np.mean(np.sin(data), axis=0), np.mean(np.cos(data), axis=0))

    kappas = np.zeros(n_features)
    for i in range(n_features):
        kappas[i] = inverse_mbessel_ratio(np.mean(np.cos(data[:, i] - means[i])))
        
    return means, kappas

def wrap(x):
    return (x + np.pi) % (2 * np.pi ) - np.pi
    

if __name__ == '__main__':
    import scipy.stats.distributions
    import matplotlib.pyplot as pp
    
    
    feature_0 = scipy.stats.distributions.vonmises.rvs(10, np.pi, size=10000)
    feature_1 = scipy.stats.distributions.vonmises.rvs(1, 0, size=10000)
    
    data = np.vstack((wrap(feature_0), wrap(feature_1))).T

    means, kappas = fit_vonmises(data)
    
    x = np.linspace(-np.pi, np.pi, 100)
    
    pp.subplot(2,1,1)
    pp.title('feature 1')
    pp.hist(data[:, 0], bins=100)
    pp.twinx().plot(x, scipy.stats.distributions.vonmises.pdf(x, kappa=kappas[0], loc=means[0]), c='r')
    
    pp.subplot(2,1,2)
    pp.title('feature 2')
    pp.hist(data[:, 1], bins=100)
    pp.twinx().plot(x, scipy.stats.distributions.vonmises.pdf(x, kappa=kappas[1], loc=means[1]), c='r')
    
    pp.show()