Интегрируйте двумерную оценку плотности ядра

Вот способ сделать это, используя интеграцию Монте-Карло. Это немного медленно, и в решении есть случайность. Ошибка обратно пропорциональна квадратному корню из размера выборки, в то время как время выполнения прямо пропорционально размеру выборки (где размер выборки относится к выборке Монте-Карло (10000 в моем примере ниже), а не к размеру вашего набора данных. ). Вот простой код, использующий ваш объект kernel.

#Compute the point below which to integrate
iso = kernel((x1,y1))

#Sample from your KDE distribution
sample = kernel.resample(size=10000)

#Filter the sample
insample = kernel(sample) < iso

#The integral you want is equivalent to the probability of drawing a point 
#that gets through the filter
integral = insample.sum() / float(insample.shape[0])
print integral

Я получаю примерно 0,2 в качестве ответа для вашего набора данных.

В настоящее время он доступен

kernel.integrate_box([-np.inf,-np.inf], [2.5,1.5])

Прямой путь к integrate

import matplotlib.pyplot as plt
import sklearn
from scipy import integrate
import numpy as np

mean = [0, 0]
cov = [[5, 0], [0, 10]]
x, y = np.random.multivariate_normal(mean, cov, 5000).T
plt.plot(x, y, 'o')
plt.show()

sample = np.array(zip(x, y))
kde = sklearn.neighbors.KernelDensity().fit(sample)
def f_kde(x,y):
    return np.exp((kde.score_samples([[x,y]])))

point = x1, y1
integrate.nquad(f_kde, [[-np.inf, x1],[-np.inf, y1]])

Проблема в том, что это очень медленно, если вы делаете это в больших масштабах. Например, если вы хотите построить линию x,y в точке x (0,100), расчет займет много времени.

Примечание: я использовал kde из sklearn, но я считаю, что вы также можете изменить его на другую форму.

Используя kernel, как определено в исходном вопросе:

import numpy as np
from scipy import stats
from scipy import integrate

def integ_func(kde, x1, y1):

    def f_kde(x, y):
        return kde((x, y))

    integ = integrate.nquad(f_kde, [[-np.inf, x1], [-np.inf, y1]])

    return integ

# Obtain data from file.
data = np.loadtxt('data.dat', unpack=True)
# Perform a kernel density estimate (KDE) on the data
kernel = stats.gaussian_kde(data)

# Define the number that will determine the integration limits
x1, y1 = 2.5, 1.5
print integ_func(kernel, x1, y1)