核密度估計（Kernel Density Estimation, KDE）演算法透過樣本估計這些樣本所屬的機率密度函式，是non-parametric方法，也就是在進行估計時無需假設分佈的具體形式。本文只討論單變數（univariate）。
數學表達
給定\(n\)個樣本\(x_i \in \mathcal{R}\)，透過KDE演算法估計\(x\)處的機率密度為
其中，核函式\(K\)和頻寬（bandwidth）\(h\)正是KDE演算法的兩個核心引數。
透過公式(1)易得，如果有\(m\)個\(x\)值要估計的機率密度，那麼計算時間複雜度是\(O(m\times n)\)。
Binning
Binning技術是透過將相鄰樣本分箱到一個box中。分箱之前，我們需要對所有的\(n\)個樣本計算核密度函式值；分箱之後，我們只需要對每個box計算一次即可，然後用它來代表所有被分到這個box內的樣本。
根據[1]，透過binning計算可分為以下三步：
-
透過將原始資料分配給相鄰網格點來對資料進行分箱以獲得網格計數。網格計數可以被認為表示其對應網格點附近的資料量。
Bin the data by assigning the raw data to neighboring grid points to obtain grid counts. A grid count can be thought of as representing the amount of data in the neighborhood of its corresponding grid point. -
計算所需的核心權重。網格點等距意味著不同核權重的數量相對較少。
Compute the required kernel weights. The fact that the grid points are equally spaced means that the number of distinct kernel weights is comparatively small. -
結合網格計數和核權重以獲得核估計的近似值。這本質上涉及一系列離散卷積。
Combine the grid counts and the kernel weights to obtain the approximation to the kernel estimate. This essentially involves a series of discrete convolutions.
步驟1涉及劃分多少個box以及如何將樣本對應到box中，對計算速度和精度都非常關鍵。兩種常見的劃分方式是
-
Simple Binning: 將樣本對應的權重劃分到距離到離其最近的grid中。
If a data point at y has surrounding grid points at x and z, then simple binning involves assigning a unit mass to the grid point closest to y. -
Linear Binning: 將樣本對應的權重按比例劃分到其左右相鄰的兩個grid中。注意給左端點x和右端點z的權重是和它到兩個端點的距離“反過來”的，也就是給x的權重是\((z-y)/(z-x)\)，而給右端點z的權重是\((y-x)/(z-x)\)。
Linear binning assigns a mass of (z - y)/(z - x) to the grid point at x, and (y - x)/(z - x) to the grid point at z.
這裡所說的樣本權重就是每個樣本是否重要性相同，可以聯絡python程式碼中sklearn.neighbors.KernelDensity.fit(sample_weight)引數，預設情況下，每個樣本地位相同，其權重都為1（或者說都是1/n）。
根據[2]，我們可以給Simple Binning一種形式化定義：定義box寬度為\(\delta\)的等距網格（其實也就是對應要劃分多少個box），令第\(j\)個bin的中點表示為\(t_j\)，其包含\(n_j\)個樣本。因此有
不難發現，在上述定義中，核心權重對應為\(n_j\)（或者說\(n_j/nh\)），也就是每個box從所有樣本中分到的權重是多少。
卷積近似
在binning的第3步中，我們提到結合網格計數和核權重本質上是在做卷積運算，在許多文獻中我們也經常可以發現為了加速KDE演算法（或解決高維KDE問題）會提到進行卷積近似操作。
首先明確卷積的概念，參考維基百科的定義，卷積是透過兩個函式\(f\)和\(g\)生成第三個函式的一種數學運算元，設\(f(t)\)和\(g(t)\)是實數\(\mathcal{R}\)上的兩個可積函式，定理兩者的卷積為如下特定形式的積分變換：
將公式(3)和公式(4)起來，就不難建立起卷積近似在KDE中的作用了。
Linear Binning + Deriche Approximation
2021年的一篇short paper[3]，表明the combination of linear binning and a recursive filter approximation by Deriche可以實現both fast and highly accurate的效果。
Luckily, 作者在Github上公佈了開原始碼；Unluckily, 程式碼語言是Javascript；於是我將其翻譯為Python版本並公佈在Github上。
然而，對於其效果我持懷疑態度，透過執行density1d.py-main，100個高斯分佈取樣點擬合得到的機率密度函式影像為
而且，更為關鍵的問題是無法該程式碼無法基於樣本生成點x處的機率密度？即，在沒有人為指定的情況下，程式碼會預設將[data.min()-pad×bandwidth, data.max()+pad×bandwidth]作為考慮範圍，然後將這個區間劃分為bins個區間，即\((x_0,x_1),(x_1,x_2),...,(x_{\text{bins}-1}, x_\text{bins})\)，那麼該方法生成的就是\(x_0,x_1,...\)這些點的機率密度，但如果我希望生成任意一個位置\(x, x\ne x_i\)處的機率密度，目前並沒有實現。
sklearn.neighbors.KernelDensity
在sklearn庫中有KDE演算法的實現，而且其考慮採用KD-tree（預設）或Ball-tree來進行加速。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
# 生成一些示例資料
np.random.seed(0)
data = np.random.normal(loc=0, scale=1, size=1000)
# 建立一個KernelDensity物件
kde = KernelDensity(bandwidth=0.5, kernel='gaussian')
# 用示例資料擬合KDE模型
kde.fit(data[:, None])
# 生成一些測試點來評估KDE模型
x = np.linspace(-5, 5, 1000)
x = x[:, None]
# 使用KDE模型評估密度
log_density = kde.score_samples(x)
# 繪製原始資料和KDE估計的密度曲線
plt.figure(figsize=(10, 6))
plt.hist(data, bins=30, density=True, alpha=0.5, color='blue')
plt.plot(x, np.exp(log_density), color='red', lw=2)
plt.title('Kernel Density Estimation (KDE)')
plt.xlabel('Value')
plt.ylabel('Density')
plt.show()
這是非常簡單的一個例子，簡單到甚至我之前沒有思考為什麼要先kde.fit()才能評估密度kde.score_sample()。透過公式(1)，我們可以發現只要清楚指定的核函式\(K\)和頻寬（bandwidth）\(h\)，那麼就可以計算出任意一點的機率密度，而這兩個引數都是在建立KDE物件時就指定好的。
起初，我以為kde.fit()是在訓練引數\(h\)，甚至有那麼一瞬間讓我覺得困擾我們如何指定超引數\(h\)的值根本不是問題，因為它只是個初始值，而在fit之後會訓練到最好的引數。這種思路明顯是錯的。因為缺乏ground_truth，因此模型本身就缺乏訓練\(h\)的能力，這也是為什麼\(h\)被稱作超引數更合適而非引數。
既然核函式和頻寬都是指定好的，那麼如果我們不fit直接去評估密度呢？
import numpy as np
from sklearn.neighbors import KernelDensity
# 生成一些示例資料
np.random.seed(0)
data = np.random.normal(loc=0, scale=1, size=1000)
# 建立一個KernelDensity物件
kde = KernelDensity(bandwidth=0.5, kernel='gaussian')
# 生成一些測試點來評估KDE模型
x = np.linspace(-5, 5, 1000)
x = x[:, None]
# 使用KDE模型評估密度
log_density = kde.score_samples(x)
Traceback (most recent call last):
File "D:\python\fast_kde\main.py", line 17, in <module>
log_density = kde.score_samples(x)
^^^^^^^^^^^^^^^^^^^^
File "F:\anaconda3\Lib\site-packages\sklearn\neighbors\_kde.py", line 261, in score_samples
check_is_fitted(self)
File "F:\anaconda3\Lib\site-packages\sklearn\utils\validation.py", line 1462, in check_is_fitted
raise NotFittedError(msg % {"name": type(estimator).__name__})
sklearn.exceptions.NotFittedError: This KernelDensity instance is not fitted yet. Call 'fit' with appropriate arguments before using this estimator.
看起來，sklearn要求在呼叫kde.score_samples()必須要kde.fit()。
透過閱讀原始碼，kde.fit()沒有改變任何KDE的引數，但構建了內建物件self.tree_，這個就是我們前面提到的KD-tree或Ball-tree，而在進行評估時也是基於這個tree進行的。我的理解是，KD-tree維護樣本之間的距離資訊，在評估一個新的資料\(x\)處的機率密度時，和它距離太遠的樣本就不計算了？
References
[1] @article{wand1994fast,
title = {Fast Computation of Multivariate Kernel Estimators},
author = {Wand, M. P.},
year = {1994},
journal = {Journal of Computational and Graphical Statistics},
volume = {3},
number = {4},
pages = {433},
doi = {10.2307/1390904}
}
[2] @article{wand1994fast,
title = {Fast Computation of Multivariate Kernel Estimators},
author = {Wand, M. P.},
year = {1994},
journal = {Journal of Computational and Graphical Statistics},
volume = {3},
number = {4},
pages = {433},
doi = {10.2307/1390904}
}
[3] @inproceedings{heer2021fast,
title = {Fast & Accurate Gaussian Kernel Density Estimation},
booktitle = {2021 IEEE Visualization Conference (VIS)},
author = {Heer, Jeffrey},
year = {2021},
pages = {11--15},
publisher = {IEEE},
address = {New Orleans, LA, USA},
doi = {10.1109/VIS49827.2021.9623323},
isbn = {978-1-66543-335-8},
langid = {english}
}