整合2D核密度估计

Question

整合2D核密度估计

pythonintegrationkernel-densityprobability-density

10

我有一个包含点的 x,y 分布，通过 scipy.stats.gaussian_kde 得到 KDE。这是我的代码和输出结果（可以从这里获取 x,y 数据）：

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()

result

红点的坐标为(x1, y1)，与2D图中的每个点一样，它具有由f（核或KDE）给出的相关值，该值介于0和0.42之间。假设f(x1, y1) = 0.08。

我需要将f与在x和y上的积分限进行积分，这些积分限由那些f评估为小于f(x1, y1)的区域给出，即：f(x, y)<0.08。

据我所见，python可以通过数值积分对函数和一维数组进行积分，但我没有看到任何可以让我对二维数组（f核）进行数值积分的东西。此外，我甚至不确定如何识别给定条件（即：f(x, y)小于某个给定值）给出的区域。

这是否可以完成？

- Gabriel

3个回答

3

目前，它是可用的。

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

（注：此处为代码示例，无法翻译。）

- Pavel Prochazka

1

一种直接的方法是将其整合（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 (0,100) 上绘制 x,y 线条，计算时间会很长。

注意：我使用了 sklearn 中的 kde，但我相信你也可以将其转换为其他形式。

使用原始问题中定义的内核：

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)

- ZK Zhao

@Gabriel 我用一个完整的例子来改了它，但是省略了一些 import。在 Python 中，import 对我来说就是一场灾难。 - ZK Zhao

是的，您可以直接使用我问题中定义的“kernel”，而不是使用“sklearn.neighbors.KernelDensity()”。 - Gabriel

1

@Gabriel 另外，我认为你不应该使用 nquad 而不是 integrate.nquad。这会使代码变得不够表达。nquad 用于表示这是一个 n-d 积分而不是一个 1d 积分。 - ZK Zhao

cqcn1991 请查看 jcrudy 回答的评论部分（在上面）。我没有注意到你的答案使用了不同的积分限制。原问题要求在域 f(x,y)<(f(x1,y1) 范围内积分，而不是 (-inf,x1) & (-inf,x1)。 - Gabriel

1

@Gabriel 对不起，我犯了一个错误，完全误解了问题。我在这里给出的答案与你的问题毫无关系。对不起。 - ZK Zhao

显示剩余4条评论

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jcrudy · Accepted Answer

以下是使用蒙特卡罗积分的一种方法。这种方法速度较慢，并且解决方案存在随机性。误差与样本量的平方根成反比，而运行时间与样本量成正比（其中样本量指蒙特卡罗样本（在下面的示例中为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。