Как реализовать рассчет критерия Колмогорова-Смирнова для двух выборок в SQL (по аналогии с stats.ks_2samp из Python)?

Question

[SpartakGusev]

Всем привет!

Возникла потребность расчёта критерия Колмогорова-Смирнова на сервере SQL для ускорения работы.
Пытаюсь разобрать логику функции stats.ks_2samp из Python.
До момента расчёта D-статистики всё понятно:
data1 = np.sort(ndarray.flatten(data1))
data2 = np.sort(ndarray.flatten(data2))
n1 = data1.shape[0]
n2 = data2.shape[0]
data_all = np.concatenate([data1, data2])
cdf1 = np.searchsorted(data1, data_all, side='right') / n1
cdf2 = np.searchsorted(data2, data_all, side='right') / n2
d = np.max(np.absolute(cdf1 - cdf2))
en = n1 * n2 / (n1 + n2)

После этого необходимо рассчитать значение p-value,
Для больших выборок (от 10 тыс) и двух-стороннего распределения далее в коде функции вызывается
prob = stats.kstwo.sf(d, np.round(en))

И вот здесь я сломался. Может кто-то объяснить, как работает эта функция sf?
Я понимаю, что функция каким-то образом строит распределение статистики исходя из входных параметров, но как это перевести на человеческий язык я не понял:

def sf( x, *args, **kwds):
"""
Survival function (1 - `cdf`) at x of the given RV.

Parameters
----------
x : array_like
quantiles
arg1, arg2, arg3,... : array_like
The shape parameter(s) for the distribution (see docstring of the
instance object for more information)
loc : array_like, optional
location parameter (default=0)
scale : array_like, optional
scale parameter (default=1)

Returns
-------
sf : array_like
Survival function evaluated at x

"""
args, loc, scale = stats.kstwo._parse_args(*args, **kwds)
x, loc, scale = map(np.asarray, (x, loc, scale))
args = tuple(map(np.asarray, args))
_a, _b = stats.kstwo._get_support(*args)
dtyp = np.find_common_type([x.dtype, np.float64], [])
cond0 = stats.kstwo._argcheck(*args) & (scale > 0)
cond1 = stats.kstwo._open_support_mask(x, *args) & (scale > 0)
cond2 = cond0 & (x <= _a)
cond = cond0 & cond1
output = np.zeros(np.shape(cond), dtyp)
place(output, (1-cond0)+np.isnan(x), stats.kstwo.badvalue)
place(output, cond2, 1.0)
if np.any(cond):
goodargs = argsreduce(cond, *((x,)+args))
place(output, cond, stats.kstwo._sf(*goodargs))
if output.ndim == 0:
return output[()]
return output

До строки place(output, (1-cond0)+np.isnan(x), stats.kstwo.badvalue) более-менее понимал что происходит, потом
вообще закипел. Что такое stats.kstwo.badvalue вообще не смог найти в интернете, хотя скорее всего это как раз построение функции распределения статистики ...

В общем - ПАМАГИТИ!!!

Answer 1

[dmshar]

Но мне нужно получить именно p-value для полученного значения статистики, как в функции Python - Дело в том, что критерий Колмогорова-Смирнова построен на факте соответствия статистики, которая получается при его расчете, с распределением Колмогорова (собственно, поэтому и критерий носит такое двойное имя). Причем связана несколько хитрым образом . Поэтому если вам из каких-то соображений надо пройти этот путь самому, "как функция sf", значит вам надо это распределение построить, а потом стандартной процедурой высчитать p-value, т.е. найти интеграл плотности вероятности, который остался правее этого значения.
Из большой любви к искусству конечно можно и этим заняться, или из-за неверия в готовые решения. А других причин для такого садомазохизма я не вижу.

P.S. И да, это наука, а не инженерия, как программирование. И реверс-инжинирингом ее не изучить.

Comments

[SpartakGusev]

Реверс-инжиниринг здесь не принципиален. Мне необходимо перенести расчёт p-value в SQL.
А способ не так важен.
Главное, что бы результат совпадал с результатом функции Python для простоты тестирования.

Вот как раз с построением распределения и расчётом площади правее статистики у меня проблема.
Если сможете хотя бы подсказать, где это можно изучить, буду очень благодарен!
Статистика, к сожалению, не мой профиль и самостоятельная попытка изучения темы пока не даёт результатов.

[dmshar]

Распределение Колмогорова и его свойста, если я правильно помню, даже в продвинутых инженерных ВУЗах не изучают. На МехМате - да, наверное еще на прикладной математике, может - на системном анализе (и то не факт) и пр. А уж всяких околонаучных веб-блог-статеек мне вообще никогда не попадалось, ни на русском, ни на английском. Поэтому изучить "это" можно, взяв соответствующий учебник по теорверу и матстатистики (причем - серьезный учебник, а не адаптацию для ....), найти там соответствующие главы, а потом лЁгко и непринужденно это дело превратить в программку. По дороге правда возникнут вопросы по интерполяции этой функции
и ее приближения для расчёта, так что придется и с численными методами "немного" разобраться. Ну, на худой конец тут может метод Монет-Карло немного облегчить ситуацию. (Кстати, что-то подобное реализовано в Scipy).

И я не знаю другого пути, увы.
Жаль, что не могу вам помочь, но так есть на самом деле.

[SpartakGusev]

В любом случае спасибо, буду разбираться дальше)

[dmshar]

Спасибо на этом форуме - вам, как новичку подскажу - принято выражать несколько по-другому.

[SpartakGusev]

Угадал?)