如何规范 scikit learn 的 kde？

Posted 2023-02-14

【中文标题】如何规范 scikit learn 的 kde？

【英文标题】：How to normalize kde of scikit learn?

【发布时间】：2019-12-17 06:30:01

【问题描述】：

假设我有一个形状为 (100000,1) 的数组，表示在 0 和 1 之间均匀分布的变量 X 的样本。 我想估计这个变量的概率密度，我使用 Scikit-Learn KernelDensity 来做到这一点。

问题是我只得到一个未标准化的结果。概率密度的积分总和不为 1。我应该如何自动归一化？我是不是做错了什么？

def kde_sklearn(data, grid, **kwargs):
    """
    Kernel Density Estimation with Scikit-learn

    Parameters
    ----------
    data : numpy.array
        Data points used to compute a density estimator. It
        has `n x p` dimensions, representing n points and p
        variables.
    grid : numpy.array
        Data points at which the desity will be estimated. It
        has `m x p` dimensions, representing m points and p
        variables.

    Returns
    -------
    out : numpy.array
        Density estimate. Has `m x 1` dimensions
    """
    kde_skl = KernelDensity(**kwargs)
    kde_skl.fit(data)
    # score_samples() returns the log-likelihood of the samples
    log_pdf = kde_skl.score_samples(grid)
    return np.exp(log_pdf) 

X = np.random.uniform(0,1,1000).reshape(-1,1)
X1 = np.linspace(0,1,100)[:,np.newaxis]

kde_sklearn(X,X1,kernel='tophat')

Out[43]: 
array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
       0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])

我希望向量为 1，因为积分的总和应为 1。

【参考方案1】：

这些是每个点的概率 - 如果

X1 = np.linspace(0,1,10000000)[:,np.newaxis]

?

你得到的数组不是随机变量的分布/样本

【讨论】：

我得到了同样的结果。一个 0.5 的向量。

看这些数字不是某些事件的概率，而是某些统计数据的 p 值，因此它不应该和 1

@quester 它们不是 p 值，而是概率密度，概率密度函数在其域上的积分应该是 1。

@quester 另外，这不是一个实际的答案。应该是评论。

【参考方案2】：

问题不在于规范化，我可以从一个例子中看出。假设我运行以下代码，将 KDE 拟合到来自标准正态分布的样本：

import numpy as np
import sklearn.neighbors as sn

# Sample from a standard normal distribution
XX = np.random.randn(1000).reshape(-1, 1)

# Fit a KDE
kde_sklg = sn.KernelDensity()
kde_sklg.fit(XX)

# Get estimated densities
XX1 = np.linspace(-4.0, 4.0, 100)[:, np.newaxis]
gdens = np.exp(kde_sklg.score_samples(XX1))

然后我可以使用梯形规则估计 PDF 下的面积，如下所示：

my_area = 0.0
for i in range(1,gdens.shape[0]):
    my_area += 0.5*(gdens[i] + gdens[i-1])*(XX1[i,0] - XX1[i-1,0])

我得到的估计面积 (my_area) 约为 0.996，非常接近 1。

问题在于您的 KDE 没有处理统一 PDF 中发生在 0 和 1 处的跳转，因此它会将它们涂抹得太多。 KDE 对您的 PDF 的估计值下大约有一半的区域最终会落在那些模糊区域之下。如果您将X1 的值替换为X2 = np.linspace(-1,2,200)[:,np.newaxis]，您可以看到KDE 在[-1,0] 和[1,2 ].

【讨论】：

不错的答案。谢啦 ：）。我将尝试在我的样本中使用更多示例来训练我的模型，我相信拖影应该会消失。

@RaphaelBenezra 我不确定，但您可能需要 [0,1] 区间之外的样本才能正常工作。您可能还想摆弄不同的内核、带宽等。

【参考方案3】：

我认为发布的答案不清楚，因此我提供了另一个答案。

简而言之，integral 的总和为 1，而不是概率。下面我展示了两种获得确实等于 1 的积分的方法。

import numpy as np
from sklearn.neighbors import KernelDensity

np.random.seed(1)

# some uniform data
X = np.random.uniform(-5,5,100).reshape(-1,1)

# grid to be used later0
grid = np.linspace(-5,5,1000)[:,np.newaxis]

# fit using the data
kde = KernelDensity(kernel = 'tophat', bandwidth= 0.5).fit(X)

# get log probailities of the grid
log_dens = kde.score_samples(grid)

# transform log prob to prob
probs = np.exp(log_dens)

# Integrate
print(np.trapz(probs.ravel(), grid.ravel()))
0.9732232232232225

plt.hist(X, density=True, bins=30)
plt.plot(grid.ravel(),probs.ravel())
plt.show()

请注意，获得积分的另一种方法如下，因为我们在定义的网格中具有相同的步骤：

np.sum(probs*np.diff(grid.ravel())[0])
0.9732232232232225