n维数据估计经验分布的累积概率

Posted 2023-03-12

【中文标题】n维数据估计经验分布的累积概率

【英文标题】：Cumulative probability of estimated empirical distribution for n-dimensional data

【发布时间】：2020-08-05 09:37:57

【问题描述】：

问题

我有一个包含 4 个数字特征和 1000 个数据点的数据集。值的分布是未知的（numpy randint 生成统一的整数，但这只是为了说明）。给定新的数据点（4 个数字），我想找出这个特定数据点的累积概率（单个数字）。

import numpy as np

data = np.random.randint(1, 100, size=(1000, 4))
array([[28, 52, 91, 66],
       [78, 94, 95, 12],
       [60, 63, 43, 37],
       ...,
       [81, 68, 45, 46],
       [14, 38, 91, 46],
       [37, 51, 68, 97]])

new_data = np.random.randint(1, 100, size=(1, 4))
array([[75, 24, 39, 94]])

我试过了：

Scipy

可以估计pdf，不知道如何估计累积概率。可能的方法是 monte-carlo sim 或集成 (scipy.integrate.nquad)，这对我的情况来说太慢了 Integrate 2D kernel density estimate。

import scipy.stats
kde = scipy.stats.gaussian_kde(data.T)
kde.pdf(new_data)

Scikit-learn

同上，不知道如何估计累积概率。

from sklearn.neighbors import KernelDensity
model = KernelDensity()
model.fit(data)
np.exp(model.score_samples(new_data))

统计模型

无法存档任何内容，因为它只接受一维数据。

from statsmodels.distributions.empirical_distribution import ECDF
ecdf = ECDF(data[:, 0])
ecdf(new_data[0][0])

问题是，是否有一种快速有效的方法来估计具有提供的 scipy 或 sklearn（最好）模型的 4 维数据点的累积概率？

我是朝着正确的方向前进还是有完全不同的方法来解决这个问题？也许变分自动编码器是要走的路？有没有简单的方法来解决这个问题？

【问题讨论】：

IIUC 你想做类似 ecdf(new_data) 的事情并得到一个结果元组，对吗？

不，我想要的是一个数据点的这 4 个特征值的单个累积概率。例如输入=数组（[[75,24,39,94]]），输出=0.01。

【参考方案1】：

一个点的多元 ecdf 只会计算值小于该点的观测值的分数。

类似下面的东西

np.random.seed(0)
data = np.random.randint(1, 100, size=(1000, 4))
new_data = np.random.randint(1, 100, size=(2, 4))

def ecdf_mv(new_data, data):
    new_data = np.atleast_2d(new_data)
    ecdf = []
    for row in new_data:
        ecdf.append((data <= row).all(1).mean())

    return np.asarray(ecdf)

ecdf_mv(new_data, data)

array([0.039, 0.002])

一些检查：

ecdf_mv(np.ones(4) * 100 / 2, data), 0.5**4
(array([0.067]), 0.0625)

marginal = 100 * np.ones((4, 4)) - 50 * np.eye(4)
ecdf_mv(marginal, data)
array([0.521, 0.515, 0.502, 0.54 ])

在单变量情况下，我们可以对数据进行排序以获得快速算法来计算原始点的 ecdf。 我不知道是否有一种数据结构或算法在计算上比蛮力比较更有效，如果 ecdf 必须在很多点上进行评估。

【讨论】：

谢谢！我明白你的意思。那真的会奏效！但是，如果我的 1000 个数据点的数据集不是同质的，我真的想估计经验概率分布并首先平滑我的模型，然后才使用它来估计新数据点的概率。我该怎么做？

此外，我可以在给定 ecdf 的情况下绘制 X 个数据点，并计算小于该点的值的观察分数。我在想是否有更好的解决方案。

ecdf 不平滑。例如，可以使用类似于 kde 的内核 cdf 来计算平滑 cdf。我对此不是很熟悉，也从未在多变量案例中尝试过。 statsmodels 有 statsmodels.org/dev/generated/… 但我猜它没有经过测试或使用太多。

这正是我想要的。一直在寻找这个几天没有运气。可能是因为正如你所说，它没有经过测试或使用太多（我想知道为什么？）。如果您愿意，请对此做出回答，我会将其标记为最佳答案。非常感谢！

【参考方案2】：

做了一些尝试和错误，我发现了以下内容：

纯 numpy 解决方案（基于 Josef 的）：

import numpy as np

def ecdf_mv(new_data, data):
    rows = np.expand_dims(new_data, axis=1)
    ecdf = (data < rows).all(axis=2).mean(axis=1)
    return np.asarray(ecdf)

这将返回一组点的经验累积分布函数。

如果需要多变量 KDE 上的 CDF，则可以使用以下代码（虽然这要慢得多）：

from statsmodels.nonparametric.kernel_density import KDEMultivariate

def cdf_kde_mv(new_data, data, data_type):
    data_kde = KDEMultivariate(data, var_type=data_type)
    return data_kde.cdf(new_data)

在性能方面，纯 numpy 可能更快，但更占用内存，并且在某些情况下，Josef 使用 for 循环的方法更快。

对于“平滑”的感觉，它不依赖于“ECDF vs CDF over Multivariate KDE”，而是更多地取决于您的样本大小和箱数。在下面的示例中，即使不使用 KDE，可视化看起来也很“流畅”，这是因为样本量对于该数量的 bin 来说足够大。被解释为 smothing 的插值是 matplotlib plot_surface。如果您还需要分析上的平滑结果，而不仅仅是图形，请考虑使用 KDE 方法。

通过一些可视化应用它：

import numpy as np
import matplotlib.pyplot as plt

sample_size = 10_000
bins_count= 20

sample = np.random.multivariate_normal([0, 0], [[1, 0.0], [0.0, 1]] , sample_size)
bins_edges = np.linspace(-4, 4, bins_count + 1)

bins_centers = (bins_edges[:-1] + bins_edges[1:]) / 2
X, Y = np.meshgrid(bins_centers, bins_centers)
x, y = X.ravel(), Y.ravel()

evaluation_points = np.stack([x, y], axis=1)

cummulative_probability = cdf_kde_mv(evaluation_points, sample, "cc")
# Alternatively:
cummulative_probability = ecdf_mv(evaluation_points, sample)
cummulative_probability = cummulative_probability.reshape(bins_count, bins_count)

# Plotting
fig = plt.figure(figsize=(10, 10), constrained_layout=True)
ax = fig.add_subplot(projection='3d')
ax.view_init(elev=20., azim=-135)
ax.plot_surface(X, Y, cummulative_probability, rcount=bins_count, ccount=bins_count, 
                antialiased=False, vmin=0, vmax=cummulative_probability.max(), 
                alpha=0.5, cmap='viridis')
ax.margins(0)
ax.set_zlim(0, 1)
plt.show()