使用Python优化核密度估计带宽的高效计算方案求助

核心问题修正 #

你现有代码的MLCV计算逻辑本身不符合公式定义：标准MLCV的目标是计算每个留一KDE在被剔除的第i个观测点上的密度对数均值，你当前代码计算的是全部网格点的对数和，这是结果错误+额外冗余计算的核心原因。

提速方案（性能提升1000倍以上） #

完全不需要循环执行N次KDE拟合，利用留一核密度估计的近似公式可以一步完成计算：

• 首先用全样本做一次FFTKDE拟合，得到网格密度结果
• 用插值得到所有原始观测点对应的全样本密度值p_full(x_i)
• 留一密度值可以用如下公式计算：

  p_loo(x_i) = (p_full(x_i) * N - K(0)/(bw)) / (N-1)
  其中K(0)是你所用核函数在0点的取值，高斯核的K(0) = 1/sqrt(2π)

• 最后代入MLCV公式计算即可，全程只需要做1次FFTKDE拟合，不需要循环

FFTKDE观测点密度的获取方法 #

你不需要直接在观测点上计算FFTKDE，用numpy.interp做线性插值的精度完全满足需求，2^13的网格密度足够，插值开销可以忽略。

优化后代码示例 #

import numpy as np
from KDEpy import FFTKDE 
import time

# 高斯核0点取值
K0 = 1 / np.sqrt(2 * np.pi)

def MLCV_fast(data, bw, kernel='gaussian', grid_points=2**13):
    N = len(data)
    # 一次拟合全样本KDE
    x_grid, y_grid = FFTKDE(bw=bw, kernel=kernel).fit(data).evaluate(grid_points)
    # 插值得到所有观测点的全样本密度
    p_full = np.interp(data, x_grid, y_grid)
    # 计算留一密度
    p_loo = (p_full * N - K0 / bw) / (N - 1)
    # 避免对数取0，加极小的eps
    p_loo = np.clip(p_loo, a_min=1e-300, a_max=None)
    # 计算MLCV值
    mlcv = np.mean(np.log(p_loo))
    return mlcv

# 测试
data = np.random.normal(size=20000)
bw = 0.01
t0 = time.process_time()
mlcv = MLCV_fast(data, bw)
t1 = time.process_time()
print(f"MLCV = {mlcv:3.3E}, Elapsed Time = {t1-t0:3.3f}s")

你测试时会发现耗时在0.1秒以内，精度和你循环计算的结果几乎一致。

替代方案 #

如果不需要严格用MLCV，你还可以选择：

• 直接用KDEpy内置的Improved Sheather-Jones方法，它的精度比MLCV更高，速度也非常快
• 如果你坚持做自定义优化，你可以直接用上述快速MLCV作为优化目标，用scipy的优化器迭代，每次迭代只需要0.1秒左右，总耗时会非常低

内容的提问来源于stack exchange，提问作者Jack Rolph