0

我想检查概率是否来自经验 CDF 指定的分布。kstest给出我认为错误的 p 值；怎么了？

我编写了一个测试函数来验证 p 值。我正在比较来自两个相同分布的样本数组，并检查从kstest和ks_2samp函数获得的 p 值。由于原假设为真（分布相同），p 值必须均匀分布在 [0,1] 上，换句话说，我必须看到错误发现率等于使用的 p 值阈值。但是，这仅适用于ks_2samp函数给出的 p 值。

from scipy import stats
from statsmodels.distributions.empirical_distribution import ECDF

def test():
    num_runs = 1000
    detected_kstest= 0
    detected_ks_2samp = 0

    for _ in range(num_runs):
        data1 = stats.norm.rvs(size=100, loc=1.0, scale=1.0)
        data2 = stats.norm.rvs(size=100, loc=1.0, scale=1.0)

        ecdf = ECDF(data1)
        p_threshold = 0.05

        _, p_val = stats.kstest(data2, ecdf)
        if p_val < p_threshold:
            detected_kstest += 1

        _, p_val = stats.ks_2samp(data1, data2)
        if p_val < p_threshold:
            detected_ks_2samp += 1

    print(f'FDR for p-value threshold {p_threshold} : kstest: {detected_kstest / num_runs}, ks_2samp: {detected_ks_2samp / num_runs}')

输出是

FDR for p-value threshold 0.05 : kstest: 0.287, ks_2samp: 0.051

我希望两个 fdr 值都接近 0.05，但是 给出的值kstest很奇怪（太高 - 换句话说，kstest经常坚持认为数据来自不同的分布）。

我错过了什么吗？

更新

正如下面所回答的，原因是kstest它不能很好地处理由小样本生成的 ecdf……唉，我必须通过也不是很大的样本来生成经验 CDF。现在，作为一种快速解决方法，我使用一些“混合”方法：

def my_ks_test(data, ecdf, ecdf_n=None):
    n = data.size
    sorted_data = np.sort(data)
    data_cdf = np.searchsorted(sorted_data, sorted_data, side='right')/(1.0 * n)

    data_cdf_by_ecdf = ecdf(sorted_data)

    d = np.max(np.absolute(data_cdf - data_cdf_by_ecdf))

    if ecdf_n is None:
        en = np.sqrt(n)
    else:
        en = np.sqrt(n * ecdf_n/float(n + ecdf_n))

    try:
        p_val = stats.distributions.kstwobign.sf((en + 0.12 + 0.11 / en) * d)
    except:
        p_val = 1.0

    return p_val

因此，它可以将用于生成 ECDF 的样本数量作为参数。也许这并不完全是严格的，到目前为止，这是我能想到的最好的。在大小均为 100 的 data1 和 data2 上进行测试时，它给出

FDR for p-value threshold 0.05 : kstest: 0.268, ks_2samp: 0.049, my_ks_test: 0.037

4

2 回答 2

0

您计算的 ECDF近似于正态分布，但如果您使用该 ECDF 从实际正态分布中测试足够大的样本，kstest将检测到该样本不是来自 ECDF。毕竟，ECDF 不是正态分布。

显然，100 的样本大小（来自实际的正态分布）足够大，kstest通常会检测到这些样本不是来自与基于 的 ECDF 关联的分布data1。

如果data1在保持大小data2固定不变的情况下增加大小，您最终会得到您期望的结果。通过增加 的大小data1，可以增加 ECDF 逼近实际正态分布的程度。

当我将创建更改data1为

        data1 = stats.norm.rvs(size=5000, loc=1.0, scale=1.0)

这是我得到的：

In [121]: test()                                                                                     
FDR for p-value threshold 0.05 : kstest: 0.048, ks_2samp: 0.0465

In [122]: test()                                                                                     
FDR for p-value threshold 0.05 : kstest: 0.0515, ks_2samp: 0.0475

In [123]: test()                                                                                     
FDR for p-value threshold 0.05 : kstest: 0.0515, ks_2samp: 0.05
于 2019-06-19T15:36:28.523 回答
0

所以我认为原因是 ECDF 函数产生了一个阶跃函数并且不做任何插值。kstest 忠实地将分布与这个“看起来很奇怪”的阶跃函数进行比较，如果没有进行更正以考虑到我们实际上正在处理阶跃函数（kstest 的“Smirnov”部分；这就是双面 ks-test 所做的）。

于 2019-06-20T08:29:36.950 回答