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

ScienceFriction 发表于 Dev

科幻小说

我有两个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功能和对象的“封闭形式”，可以提供正确的解决方案。

谢谢！

萨沙

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

数据生成，包括随机抽样
拟合
交叉发现

方法

基本思想是使用一些通用的寻根算法。为此，我们正在使用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()