Zintegruj oszacowanie gęstości jądra 2D

Question

Mam rozkład punktów, dla których uzyskuję KDE do . To jest mój kod i jak wygląda wyjściowe (dane x,y można otrzymać od here):

import numpy as np 
from scipy import stats 

# Obtain data from file. 
data = np.loadtxt('data.dat', unpack=True) 
m1, m2 = data[0], data[1] 
xmin, xmax = min(m1), max(m1) 
ymin, ymax = min(m2), max(m2) 

# Perform a kernel density estimate (KDE) on the data 
x, y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] 
positions = np.vstack([x.ravel(), y.ravel()]) 
values = np.vstack([m1, m2]) 
kernel = stats.gaussian_kde(values) 
f = np.reshape(kernel(positions).T, x.shape) 

# Define the number that will determine the integration limits 
x1, y1 = 2.5, 1.5 

# Perform integration? 

# Plot the results: 
import matplotlib.pyplot as plt 
# Set limits 
plt.xlim(xmin,xmax) 
plt.ylim(ymin,ymax) 
# KDE density plot 
plt.imshow(np.rot90(f), cmap=plt.cm.gist_earth_r, extent=[xmin, xmax, ymin, ymax]) 
# Draw contour lines 
cset = plt.contour(x,y,f) 
plt.clabel(cset, inline=1, fontsize=10) 
plt.colorbar() 
# Plot point 
plt.scatter(x1, y1, c='r', s=35) 
plt.show() 

czerwony punkt o współrzędnych (x1, y1) ma (jak w każdym punkcie na wykresie 2D) skojarzony wartość podana przez f (jądro lub KDE) między 0 a 0,42. Powiedzmy, że f(x1, y1) = 0.08.

muszę zintegrować f z limitami integracyjnych w x i y podanym przez tych regionach, gdzie f ocenia się mniej niż f(x1, y1), tj: f(x, y)<0.08.

Na co widziałem python może przeprowadzić integrację funkcji i jednego tablic wymiarowych poprzez integrację numerycznej, ale nie mam nic, że pozwolił mi wykonać całkowanie numeryczne na tablicy 2D (jądro f) widziany Co więcej, nie jestem pewien, w jaki sposób rozpoznałbym regiony podane przez ten konkretny warunek (tj .: f(x, y) mniej niż podana wartość)

Czy można to w ogóle zrobić?

Answer 1

Oto sposób wykonania tego przy użyciu integracji monte carlo. Jest trochę powolny, a rozwiązanie ma losowość. Błąd jest odwrotnie proporcjonalny do pierwiastka kwadratowego z wielkości próbki, podczas gdy czas działania jest wprost proporcjonalny do wielkości próbki (gdzie wielkość próbki odnosi się do próbki Monte Carlo (10000 w moim przykładzie poniżej), a nie do wielkości zbioru danych). Oto prosty kod wykorzystujący obiekt 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 

Otrzymuję około 0,2 jako odpowiedź dla twojego zestawu danych.

Answer 2

Bezpośrednim sposobem jest 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]]) 

Problemem jest to, że jest to bardzo powolny, jeśli robisz to w dużej skali. Na przykład, jeśli chcesz wykreślić linię x,y na x (0,100), obliczenie zajmie dużo czasu.

Uwaga: użyłem kde z sklearn, ale uważam, że możesz również zmienić go na inną.

---

Używanie kernel jak określono w pierwotnym pytaniu:

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)