scipy gaussian_kde и круговые данные

Question

Я использую scipys gaussian_kde, чтобы получить плотность вероятности некоторых бимодальных данных. Однако, поскольку мои данные являются угловыми (это направления в градусах), у меня возникает проблема, когда значения происходят вблизи пределов. В приведенном ниже коде приведены два примера kde, когда домен 0-360 под оценками, поскольку он не может справиться с круговым характером данных. PDF необходимо определить на единичном круге, но я не могу найти что-либо в scipy.stats, подходящем для данных этого типа (распределение von mises существует, но работает только для унимодальных данных). Кто-нибудь там сталкивался с этим раньше? Есть ли что-нибудь (предпочтительное на основе python) для оценки бимодальных PDF-файлов на единичном круге?

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

Итак, у меня есть то, что я считаю разумным решением. В основном я использую распределение фон Мизеса в качестве базовой функции для оценки плотности ядра. Код ниже, если он полезен для кого-либо еще.

def vonmises_KDE(data, kappa, plot=None):  

    """  
    Create a kernal densisity estimate of circular data using the von mises 
    distribution as the basis function. 

    """ 

    # imports 
    from scipy.stats import vonmises 
    from scipy.interpolate import interp1d 

    # convert to radians 
    data = np.radians(data) 

    # set limits for von mises 
    vonmises.a = -np.pi 
    vonmises.b = np.pi 
    x_data = np.linspace(-np.pi, np.pi, 100) 

    kernels = [] 

    for d in data: 

     # Make the basis function as a von mises PDF 
     kernel = vonmises(kappa, loc=d) 
     kernel = kernel.pdf(x_data) 
     kernels.append(kernel) 

     if plot: 
      # For plotting 
      kernel /= kernel.max() 
      kernel *= .2 
      plt.plot(x_data, kernel, "grey", alpha=.5) 


    vonmises_kde = np.sum(kernels, axis=0) 
    vonmises_kde = vonmises_kde/np.trapz(vonmises_kde, x=x_data) 
    f = interp1d(x_data, vonmises_kde) 


    if plot: 
     plt.plot(x_data, vonmises_kde, c='red') 

    return x_data, vonmises_kde, f 

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]) 
kappa = 12 
x_data, vonmises_kde, f = vonmises_KDE(baz, kappa, plot=1) 

Answer 2

ответ Дэйва не является правильным, потому что scipy «s vonmises не обернуть вокруг [-pi, pi].

Вместо этого вы можете использовать следующий код, основанный на том же принципе. Он основан на уравнениях, описанных в 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 

Вот пример

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

Answer 3

Вот быстрый приближение к @ kingjr более точный ответ:

def vonmises_pdf(x, mu, kappa): 
    return np.exp(kappa * np.cos(x - mu))/(2. * np.pi * scipy.special.i0(kappa)) 


def vonmises_fft_kde(data, kappa, n_bins): 
    bins = np.linspace(-np.pi, np.pi, n_bins + 1, endpoint=True) 
    hist_n, bin_edges = np.histogram(data, bins=bins) 
    bin_centers = np.mean([bin_edges[1:], bin_edges[:-1]], axis=0) 
    kernel = vonmises_pdf(
     x=bin_centers, 
     mu=0, 
     kappa=kappa 
    ) 
    kde = np.fft.fftshift(np.fft.irfft(np.fft.rfft(kernel) * np.fft.rfft(hist_n))) 
    kde /= np.trapz(kde, x=bin_centers) 
    return bin_centers, kde 

Test (С использованием tqdm для прогресса бара и времени, а также Matplotlib для проверки результатов):

import numpy as np 
from tqdm import tqdm 
import scipy.stats 
import matplotlib.pyplot as plt 

n_runs = 1000 
n_bins = 100 
kappa = 10 

for _ in tqdm(xrange(n_runs)): 
    bins1, kde1 = vonmises_kde(
     data=np.r_[ 
      np.random.vonmises(-1, 5, 1000), 
      np.random.vonmises(2, 10, 500), 
      np.random.vonmises(3, 20, 100) 
     ], 
     kappa=kappa, 
     n_bins=n_bins 
    ) 


for _ in tqdm(xrange(n_runs)): 
    bins2, kde2 = vonmises_fft_kde(
     data=np.r_[ 
      np.random.vonmises(-1, 5, 1000), 
      np.random.vonmises(2, 10, 500), 
      np.random.vonmises(3, 20, 100) 
     ], 
     kappa=kappa, 
     n_bins=n_bins 
    ) 

plt.figure() 
plt.plot(bins1, kde1, label="kingjr's solution") 
plt.plot(bins2, kde2, label="dolf's FFT solution") 
plt.legend() 
plt.show() 

Результаты:

100%|██████████| 1000/1000 [00:07<00:00, 135.29it/s] 
100%|██████████| 1000/1000 [00:00<00:00, 1945.14it/s] 

(1945/135 = 14 раз быстрее)

Для еще большей скорости используйте целую мощность 2 в качестве количества ящиков. Он также масштабируется лучше (то есть он быстро сохраняется со множеством ящиков и большим количеством данных). На моем ПК это в 118 раз быстрее исходного ответа с n_bins = 1024.

Почему это работает?

Произведение FFTs двух сигналов (без нулевого заполнения) равен circular (or cyclic) convolution двух сигналов. A kernel density estimation - это в основном ядро, свернутое с сигналом, который имеет импульс в позиции каждой точки данных.

Почему это не так?

Поскольку я использую гистограмму для равномерного распределения данных, я теряю точное положение каждого образца и использую только центр бункера, к которому он принадлежит. Количество выборок в каждом бункере используется как величина импульса в этой точке. Например: Игнорируя нормализацию на мгновение, если у вас есть бункер от 0 до 1 и два образца в этом бункере, равный 0,1 и 0,2, KDE будет the kernel centred around 0.1 + the kernel centred around 0.2. Аппроксимация будет равна 2x `ядра с центром около 0,5, который является центром бункера.