kstest 给出奇怪的 p 值
kstest gives strange p-values
我想检查概率是否来自经验 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 ...
las,我必须通过同样不是很大的样本生成经验 CDF。
现在,作为一种快速解决方法,我使用了一些 'hybrid' 方法:
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
您计算的 ECDF 近似于 正态分布,但是如果您使用 实际 正态分布测试足够大的样本ECDF,kstest
会检测到样本不是来自ECDF。毕竟ECDF不是正态分布
显然,样本量 100(来自实际正态分布)足够大,以至于 kstest
经常检测到这些样本不是来自与基于 data1
的 ECDF 关联的分布.
如果您增加 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
所以我认为原因是 ECDF 函数产生了一个 step-function 而没有进行任何插值。 kstest 忠实地将分布与此 'strangely-looking' step-function 进行比较,当然会检测到差异,如果没有进行更正以考虑到我们实际上正在处理 step-function ('Smirnov' kstest 的一部分;这就是 two-sided ks-test 所做的)。
我想检查概率是否来自经验 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 ...
las,我必须通过同样不是很大的样本生成经验 CDF。
现在,作为一种快速解决方法,我使用了一些 'hybrid' 方法:
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
您计算的 ECDF 近似于 正态分布,但是如果您使用 实际 正态分布测试足够大的样本ECDF,kstest
会检测到样本不是来自ECDF。毕竟ECDF不是正态分布
显然,样本量 100(来自实际正态分布)足够大,以至于 kstest
经常检测到这些样本不是来自与基于 data1
的 ECDF 关联的分布.
如果您增加 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
所以我认为原因是 ECDF 函数产生了一个 step-function 而没有进行任何插值。 kstest 忠实地将分布与此 'strangely-looking' step-function 进行比较,当然会检测到差异,如果没有进行更正以考虑到我们实际上正在处理 step-function ('Smirnov' kstest 的一部分;这就是 two-sided ks-test 所做的)。