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 

我想知道是否有一个利用 KDE 函数和对象的“封闭形式”可以提供正确的解决方案。

谢谢!

Answer 1

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

• 数据生成，包括随机采样
• kde 拟合
• 寻找交叉点

方法

基本思想是使用一些通用的求根算法。为此，我们使用 brentq来自 scipy。

代码

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，这应该更快、更稳定