Как реализовать KS-тест в Python - PullRequest
Новая
       102

Как реализовать KS-тест в Python

0 голосов
/ 13 мая 2019

scipy.stats.kstest(rvs, cdf, N) может выполнить KS-Test для набора данных rvs. Он проверяет, соответствует ли набор данных распределению пригодности, значение которого cdf указано в параметрах этого метода.

Рассмотрим теперь набор данных N=4800 выборок. Я выполнил KDE на этих данных и, следовательно, у меня есть приблизительный PDF. Этот PDF очень похож на бимодальный дистрибутив. При построении расчетного PDF и кривой_подгонки к нему бимодального распределения эти два графика в значительной степени идентичны. Параметры подобранного бимодального распределения: (scale1, mean1, stdv1, scale2, mean2, stdv2): [0.6 0.036 0.52, 0.23 1.25 0.4]

Как я могу применить scipy.stats.kstest, чтобы проверить, является ли мой предполагаемый PDF двухмодальным? В качестве моей нулевой гипотезы я утверждаю, что предполагаемый PDF равен следующему PDF: 

hypoDist = 0.6*norm(loc=0, scale=0.2).pdf(x_grid) + 0.3*norm(loc=1, scale=0.2).pdf(x_grid)
hypoCdf = np.cumsum(hypoDist)/len(x_grid)

x_grid - это просто вектор, содержащий значения x, по которым я оцениваю свой предполагаемый PDF. Таким образом, каждая запись pdf имеет соответствующее значение x_grid. Возможно, мои вычисления hypoCdf неверны. Может быть, вместо деления на len(x_grid), я должен делить на np.sum(hypoDist)?

Задача: cdf параметр kstest нельзя указывать как бимодальный. Я также не могу указать, что это hypoDist.

Если бы я хотел проверить, был ли мой набор данных распределен по Гауссу, я бы написал: 

KS_result = kstest(measurementError, norm(loc=mean(pdf), scale=np.std(pdf)).cdf)
print(KS_result)

measurementError - это набор данных, на котором я выполнил KDE. Это возвращает: statistic=0.459, pvalue=0.0 Для меня это немного раздражает, что значение 0,0

1 Ответ

1 голос
/ 14 мая 2019

Аргумент cdf для kstest может быть вызываемым , который реализует накопительную функцию распределения для распределения, с которым вы хотите проверить свои данные.Чтобы использовать его, вы должны внедрить CDF вашего бимодального дистрибутива.Вы хотите, чтобы распределение было смесью двух нормальных распределений.Вы можете реализовать CDF для этого распределения, рассчитав взвешенную сумму CDF двух нормальных распределений, составляющих смесь.

Вот скрипт, который показывает, как вы можете это сделать.Чтобы продемонстрировать, как используется kstest, скрипт запускает kstest дважды.Сначала он использует образец, который не из дистрибутива.Как и ожидалось, kstest вычисляет очень маленькое значение p для этого первого образца.Затем он генерирует образец, который взят из смеси.Для этого примера значение p не мало. 

import numpy as np
from scipy import stats


def bimodal_cdf(x, weight1, mean1, stdv1, mean2, stdv2):
    """
    CDF of a mixture of two normal distributions.
    """
    return (weight1*stats.norm.cdf(x, mean1, stdv1) +
            (1 - weight1)*stats.norm.cdf(x, mean2, stdv2))


# We only need weight1, since weight2 = 1 - weight1.
weight1 = 0.6
mean1 = 0.036
stdv1 = 0.52
mean2 = 1.25
stdv2 = 0.4

n = 200

# Create a sample from a regular normal distribution that has parameters
# similar to the bimodal distribution.
sample1 = stats.norm.rvs(0.5*(mean1 + mean2), 0.5, size=n)

# The result of kstest should show that sample1 is not from the bimodal
# distribution (i.e. the p-value should be very small).
stat1, pvalue1 = stats.kstest(sample1, cdf=bimodal_cdf,
                              args=(weight1, mean1, stdv2, mean2, stdv2))
print("sample1 p-value =", pvalue1)

# Create a sample from the bimodal distribution.  This sample is the
# concatenation of samples from the two normal distributions that make
# up the bimodal distribution.  The number of samples to take from the
# first distributions is determined by a binomial distribution of n
# samples with probability weight1.
n1 = np.random.binomial(n, p=weight1)
sample2 = np.concatenate((stats.norm.rvs(mean1, stdv1, size=n1),
                         (stats.norm.rvs(mean2, stdv2, size=n - n1))))

# Most of time, the p-value returned by kstest with sample2 will not
# be small.  We expect the value to be uniformly distributed in the interval
# [0, 1], so in general it will not be very small.
stat2, pvalue2 = stats.kstest(sample2, cdf=bimodal_cdf,
                              args=(weight1, mean1, stdv1, mean2, stdv2))
print("sample2 p-value =", pvalue2)

Типичный вывод (числа будут отличаться при каждом запуске скрипта): 

sample1 p-value = 2.8395166853884146e-11
sample2 p-value = 0.3289374831186403

Youможет оказаться, что для вашей проблемы этот тест не работает хорошо.У вас есть 4800 сэмплов, но в вашем коде есть параметры, числовые значения которых имеют только одну или две значащие цифры.Если у вас нет веских оснований полагать, что ваш образец взят из распределения с точно этими параметрами, вполне вероятно, что kstest вернет очень небольшое значение р.

...