scipy gaussian_kde and circular data

Question

I am using scipys gaussian_kde to get probability density of some bimodal data. However, as my data is angular (it's directions in degrees) I have a problem when values occur near the limits. The code below gives two example kde's, when the domain is 0-360 it under estimates as it cannot deal with the circular nature of the data. The pdf needs to be defined on the unit circle but I can't find anything in scipy.stats suitable to this type of data (von mises distribution is there but only works for unimodal data). Has anyone out there ran into this one before? Is there anything (preferable python based) available to estimate bimodal pdf's on the unit circle?

import numpy as np
import scipy as sp
from pylab import plot,figure,subplot,show,hist
from scipy import stats



baz = np.array([-92.29061004, -85.42607874, -85.42607874, -70.01689348,
               -63.43494882, -63.43494882, -70.01689348, -70.01689348,
               -59.93141718, -63.43494882, -59.93141718, -63.43494882,
               -63.43494882, -63.43494882, -57.52880771, -53.61564818,
               -57.52880771, -63.43494882, -63.43494882, -92.29061004,
               -16.92751306, -99.09027692, -99.09027692, -16.92751306,
               -99.09027692, -16.92751306,  -9.86580694,  -8.74616226,
                -9.86580694,  -8.74616226,  -8.74616226,  -2.20259816,
                -2.20259816,  -2.20259816,  -9.86580694,  -2.20259816,
                -2.48955292,  -2.48955292,  -2.48955292,  -2.48955292,
                 4.96974073,   4.96974073,   4.96974073,   4.96974073,
                -2.48955292,  -2.48955292,  -2.48955292,  -2.48955292,
                -2.48955292,  -9.86580694,  -9.86580694,  -9.86580694,
               -16.92751306, -19.29004622, -19.29004622, -26.56505118,
               -19.29004622, -19.29004622, -19.29004622, -19.29004622])


xx = np.linspace(-180, 180, 181)
scipy_kde = stats.gaussian_kde(baz)              
print scipy_kde.integrate_box_1d(-180,180)

figure()
plot(xx, scipy_kde(xx), c='green')             

baz[baz<0] += 360             
xx = np.linspace(0, 360, 181)
scipy_kde = stats.gaussian_kde(baz)              
print scipy_kde.integrate_box_1d(-180,180)
plot(xx, scipy_kde(xx), c='red')             

Answer 1

Dave's answer isn't correct, because scipy's vonmises doesn't wrap around [-pi, pi].

Instead you can use the following code, which is based on the same principle. It is based on the equations described in numpy.

def vonmises_kde(data, kappa, n_bins=100):
    from scipy.special import i0
    bins = np.linspace(-np.pi, np.pi, n_bins)
    x = np.linspace(-np.pi, np.pi, n_bins)
    # integrate vonmises kernels
    kde = np.exp(kappa*np.cos(x[:, None]-data[None, :])).sum(1)/(2*np.pi*i0(kappa))
    kde /= np.trapz(kde, x=bins)
    return bins, kde

Here is an example

import matplotlib.pyplot as plt
import numpy as np
from numpy.random import vonmises

# generate complex circular distribution
data = np.r_[vonmises(-1, 5, 1000), vonmises(2, 10, 500), vonmises(3, 20, 100)]

# plot data histogram
fig, axes = plt.subplots(2, 1)
axes[0].hist(data, 100)

# plot kernel density estimates
x, kde = vonmises_kde(data, 20)
axes[1].plot(x, kde)

Histogram and kernel density plots