kstest дает странные p-значения - PullRequest
Новая
       35

kstest дает странные p-значения

0 голосов
/ 19 июня 2019

Я хочу проверить, исходит ли вероятность из распределения, определенного эмпирическим CDF.kstest дает p-значения, которые я считаю неправильными;в чем дело?

Я написал тестовую функцию для проверки p-значений.Я сравниваю выборочные массивы из двух идентичных распределений и проверяю p-значения, полученные из функций kstest и ks_2samp.Так как нулевая гипотеза верна (распределения идентичны), p-значения должны быть равномерно распределены на [0,1], другими словами, я должен видеть частоту ложных открытий, равную используемому порогу p-значения.Однако это относится только к значениям p, заданным функцией ks_2samp. 

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.Может быть, это не совсем естественно с точки зрения натуры, на данный момент это лучшее, что я могу придумать.При тестировании на data1 и data2 обоих размеров 100 он дает

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

Ответы [ 2 ]

0 голосов
/ 20 июня 2019

Так что я думаю, что причина в том, что функция ECDF производит пошаговую функцию и не выполняет никакой интерполяции. kstest добросовестно сравнивает распределение с этой «странно выглядящей» степ-функцией и, конечно, обнаруживает разницу, если не вносятся исправления, чтобы учесть, что мы на самом деле имеем дело с степ-функцией (часть «kstest» Смирнова); это то, что делает двусторонний тест ks).

0 голосов
/ 19 июня 2019

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
...