python：找到两个gaussian_kde函数（对象）的交集

Question

我有两个 python gaussian_kde 对象，我想找到交点。有简单的方法吗？

请注意，这些函数没有很好地参数化，见图。

这是一种简单的方法（假设只有一个交叉点，但可以很容易地为范围内的所有交叉点修改它，因为指定的交叉点不超过一个 init_interval):

def find_intersection(kde1, kde2, init_interval=0.01, scope =[0,1], convergence=0.0001):
x_left = scope[0]
x_right = scope[0]+init_interval
while x_right < scope[1]:
    left = kde1(x_left)[0]-kde2(x_left)[0]
    right = kde1(x_right)[0]-kde2(x_right)[0]
    if left*right < 0: #meaning the functions intersected (an odd number of times) in the interval
        if init_interval <= convergence:
            return x_right
        else: 
            return find_intersection(kde1, kde2, init_interval/10, scope=[x_left, x_right])
    else: #no intersection or an even number of intersections in the interval
        x_left = x_right
        x_right+=init_interval
return scope[0]-1 #out of scope means no intersection

对于地块的 KDE，我们得到：

>>>from scipy.stats import gaussian_kde
>>>data1 = d_sp.values()
>>>density1 = gaussian_kde(data1)
>>>data2 = d_xp.values()
>>>density2 = gaussian_kde(data2)
>>>xs = np.linspace(0,.2,200)
>>>print find_intersection(density1, density2) 
0.0403   
>>>print find_intersection(density1, density2, convergence=0.000001)
0.0403

我想知道是否有 "closed form" 利用 KDE 功能和对象可以提供正确的解决方案。

谢谢！

Answer 1

如果没有代码很难提供帮助，但我实现了一个完整的示例，包括：

数据生成包括随机抽样
kde-拟合
找交点

方法

基本思想是使用一些通用的寻根算法。为此，我们使用 scipy.

中的 brentq

代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import brentq
from sklearn.neighbors.kde import KernelDensity

# Generate normal functions
x_axis = np.linspace(-3, 3, 100)
gaussianA = norm.pdf(x_axis, 2, 0.5)  # mean, sigma
gaussianB = norm.pdf(x_axis, 0.1, 1.5)

# Random-sampling from functions
a_samples = norm.rvs(2, 0.5, size=100)
b_samples = norm.rvs(0.1, 1.5, size=100)

# Fit KDE
def kde_sklearn(x, x_grid, bandwidth=0.2, **kwargs):
    """Kernel Density Estimation with Scikit-learn"""
    kde_skl = KernelDensity(bandwidth=bandwidth, **kwargs)
    kde_skl.fit(x[:, np.newaxis])
    # score_samples() returns the log-likelihood of the samples
    log_pdf = kde_skl.score_samples(x_grid[:, np.newaxis])
    return kde_skl, np.exp(log_pdf)

kdeA, pdfA = kde_sklearn(a_samples, x_axis, bandwidth=0.25)
kdeB, pdfB = kde_sklearn(b_samples, x_axis, bandwidth=0.25)

# Find intersection
def findIntersection(fun1, fun2, lower, upper):
    return brentq(lambda x : fun1(x) - fun2(x), lower, upper)

funcA = lambda x: np.exp(kdeA.score_samples([[x]][0]))
funcB = lambda x: np.exp(kdeB.score_samples([[x]][0]))

result = findIntersection(funcA, funcB, -3, 3)

# Plot
f, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
ax1.plot(x_axis, gaussianA, color='green')
ax1.plot(x_axis, gaussianB, color='blue')
ax1.set_title('Original Gaussians')
ax2.plot(x_axis, pdfA, color='green')
ax2.plot(x_axis, pdfB, color='blue')
ax2.set_title('KDEs of subsampled Gaussians')
ax2.axvline(result, color='red')
plt.show()

输出

备注

brentq 似乎是最常见的求根算法（因为它稳定且快速）但根据您的数据，可能需要进行参数调整
可以切换到一些other优化算法
（出于建模目的进行了一些简化；例如，通常 kde-bandwith 选择应该通过交叉验证来获得一些东西比我例子中的更好）

编辑：从 fsolve 切换到 brentq，应该更快更稳定

Answer 2

好吧，不打算写代码，但我想了一些主意

从 PDF 的高斯 KDe 估计可以得到 CDF 的估计，这将是无精度损失步骤，从高斯加权和到误差函数的加权和。

所以你有 CDF1(x) 和 CDF2(x)。您构建 f(x)=CDF1(x)-CDF2(x)。它等于在 xmin 处为 0，在 xmax 处为 0。如果我们要找到 f(x) 的 min/max，我们应该计算 f(x) 的导数并检查导数在哪里等于 0。

所以，f'(x)=(CDF1(x)-CDF2(x))' = PDF1(x) - PDF2(x)，你猜怎么着 - 在 PDF 交点处它等于 0观点。因此，我们可以将搜索交叉点问题转换为搜索 min/max 问题。

所以我建议构建 f(x) 并将其传递给某些 min/max 查找例程，例如 fmin() 或 scipy.optimize 中的 minimize_scalar。

python：找到两个gaussian_kde函数（对象）的交集

python: finding the intersection of two gaussian_kde functions (objects)

python

probability

gaussian

方法

代码

输出

备注