vasilymat 22 мая в 17:40

Что у меня за распределение

Средний

14 мин

850

Блог компании Т-БанкМашинное обучение*Искусственный интеллектСтатистика в IT

Туториал

Привет! Меня зовут Василий, я ML инженер центра ML-экспертизы в обслуживании. В Т-Банке мы прогнозируем нагрузку на колл-центр: сколько придет обращений, какой длительности и некоторые другие характеристики.

Нужно уметь переводить эту нагрузку в количество людей на линии. Для этого мы реализовали симулятор колл-центра. Для работы симулятора нужно понимать, как распределены те или иные случайные величины, и иметь численные характеристики соответствия наблюдаемых значений теоретическим распределениям. Это называется задачей (критериями) согласия: к какому теоретическому распределению вероятностей принадлежит данная выборка.

«Используй Колмогорова — Смирнова, да и все тут», — скажут мне. В принципе, да, но нет. Мы пойдем чуть дальше и попытаемся разобраться, как все устроено и какие есть ограничения. Расскажу, почему нельзя просто так применять критерии согласия, к каким это приводит ошибкам и как с этим быть.

Основные определения

В прошлой статье мы концентрировались на сравнении двух средних выборок с помощью разных критериев, для этого используется p-value. Для определения, к какому классу распределения отнести нашу выборку, хотелось бы иметь что-то подобное.

Правильно определить класс распределения необходимо, чтобы наше имитационное моделирование колл-центра (или другого процесса) было корректным. И тут нам важно ввести несколько понятий, которые будут важны для нас на протяжении всей статьи: гипотезы, мощность критерия, равная $1-\beta$ , и подходы к критериям согласия.

Гипотезы могут быть простыми и сложными. Например, простая гипотеза звучит так: выборка принадлежит стандартному нормальному распределению, то есть распределению с матожиданием 0 и дисперсией 1.

Сложная гипотеза звучит так: распределение принадлежит нормальному распределению с неизвестными параметрами. Эти параметры еще нужно определить по самой выборке, причем всегда с конечной точностью. Именно с такими задачами мы сталкиваемся чаще всего, и именно они порождают большие сложности.

Мощность критерия, которая равна $1-\beta$ , где $\beta$ показывает вероятность совершить ошибку второго рода — принять нулевую гипотезу, когда она неверна. Чем более критерий мощный, тем лучше он будет замечать неверность нашей нулевой гипотезы.

Например, есть два критерия (функции, которым скармливается выборка). Наша нулевая гипотеза: выборка получена из логнормального распределения. А на самом деле на вход мы подаем гамма-распределение. Более мощный критерий будет реже соглашаться с нулевой гипотезой.

Среди критериев согласия можно выделить два подхода. Первый основывается на функции распределения вероятности. Самый известный из критериев этого типа — критерий Колмогорова — Смирнова.

Второй подход основывается на плотности вероятности и частотах. Самый известный тут — тест $\chi^2,$ он требует группировки экспериментальных данных.

Мы опишем оба эти подхода, а потом сравним с методом численного решения задачи — методом Монте-Карло. «Почему с Монте-Карло?» — спросите вы. Потому что общего теоретического решения задачи не существует. Но обо всем по порядку.

Критерий Колмогорова — Смирнова

Критерий Колмогорова — Смирнова используется для непрерывных и обычно одномерных распределений. Он сравнивает эмпирическую и теоретическую функции распределения. И выглядит так:

— это эмпирическая функция распределения, — теоретическая.

Мы берем две функции распределения и находим максимальное различие между ними. Пусть не смущает, что тут мы ищем только одну точку (а не какое-то среднее).

Функция распределения кумулятивная, это и есть уже нечто усредненное. Чтобы получить различие, нужно. чтобы для каждого значения различия накапливались.

В Scipy тест реализован функцией stats.ks_1samp.

Супремум — это, к сожалению, не максимум. Например, отрезок [0,1) не имеет максимума, но супремум у него 1.

Как и в случае со сравнением средних, максимальное различие функций распределения тоже случайная величина, у которой есть своя функция плотности вероятности. Ее мы и используем для определения p-value.

T-распределение возникает, если дисперсия генеральной совокупности неизвестна. И ее мы вынуждены определять по самой экспериментальной выборке. Знай мы ее, t-распределение просто превратилось бы в гауссово и не имело таких тяжелых хвостов. В гауссово оно и превращается по мере увеличения выборки и уточнения истинной дисперсии.

Как мы помним, сложная гипотеза — это та, в которой параметры распределения нужно определять по выборке. Но вот беда: аналитического решения (то есть своего t-распределения) для него в общем виде нет.

Если мы хотим ответить на вопрос «Является ли эта выборка семплированной из гамма-распределения с параметрами $\alpha = \alpha_0, \beta = \beta_0,$ значения статистики будут зависеть как от размера выборки, так и от значения ее параметров.

Чтобы немного понять масштаб проблемы, представим, что сравниваем два средних. И вот p-value в каком-нибудь t-тесте зависит от того, сравниваем мы 3 и 4 или 30 и 40. Не говоря уже про размер выборки. В t-тесте он идет как количество степеней свободы, а тут так не выходит.

Все это приводит к тому, что нужно проводить моделирование мощности критерия для каждого значения параметров, для каждого распределения и еще для каждого способа нахождения параметров распределения по выборке (хотя для ряда распределений есть точные решения).

Приведем пример. Сгенерируем выборку, например, из гамма-распределения определенного размера. Выполним тест Колмогорова — Смирнова как тест на простую гипотезу, то есть принадлежит ли эта выборка гамма-распределению с параметрами генерации.

rvs = ss.gamma(1.3, 10, 400).rvs(100,random_state=42)
pvalue = ss.ks_1samp( rvs,ss.gamma.cdf
                  	,(1.3, 10, 400) ).pvalue
# p-value = 0.65, а значит, мы уверенно принимаем данную гипотезу
print('KS p-value простой гипотезы: ',pvalue)

Теперь представим, что мы не знаем параметров распределения и сфитируем
их каким-нибудь методом.

def fit_params(rvs):
	# По умолчанию используется метод максимального правдоподобия
	params_mle = ss.gamma.fit(rvs)
	# Можем использовать метод моментов
	params_mm  = ss.gamma.fit(rvs, method='MM')

	print( 'Параметры с максимизацией правдоподобия (MLE): ',  params_mle)
	print( 'KS p-value MLE: ',  ss.ks_1samp( rvs,ss.gamma.cdf,(params_mle ) ).pvalue )
	print( 'Параметры методом моментов (MM): ',  params_mm)
	print( 'KS p-value MM: ',  ss.ks_1samp( rvs,ss.gamma.cdf,(params_mm ) ).pvalue )

	return params_mle, params_mm

params_mle, params_mm =  fit_params(rvs)

>KS p-value простой гипотезы:  0.649
>Параметры с максимизацией правдоподобия (MLE):  (0.099, 26.503, 3.276)
>KS p-value MLE:  3.124e-170
>Параметры методом моментов (MM):  (1.268, 57.891, 352.505)
>KS p-value MM:  0.927

Мы вызвали метод fit, в Scipy у него есть два параметра: MLE — максимизация функции правдоподобия, а также MM — метод моментов. Обычно MLE дает более точный результат и работает дольше. Но может и совершенно некорректно обучиться на малых выборках.

Код для отрисовки графика:

def plot_pdfs(params_mle,params_mm,N=100):
	x = np.arange(0,1500)

	fig = go.Figure()
	fig.add_trace(go.Scatter(x=x, y = ss.gamma(1.3, 10, 400).pdf(x), name = 'Исходное распределение'))
	fig.add_trace(go.Scatter(x=x, y = ss.gamma(*params_mle).pdf(x), name = 'Фитированное MLE'))
	fig.add_trace(go.Scatter(x=x, y = ss.gamma(*params_mm).pdf(x), name = 'Фитированное MM'))

	fig.
update_layout( yaxis=dict(range=[0, .005]), title=f'Cравнени pdf, {N} значений в выборке' )

	fig.show()

plot_pdfs(params_mle,params_mm)

С одной стороны, MLE у нас дал просто некорректные результаты, и этого можно не заметить, если что-то массово тестировать. С другой стороны, MM обучился на имеющейся выборке, из-за чего завысил p-value. Он сравнивал выборку не с распределением с исходными параметрами, а с теми, которые лучше подходят к полученной выборке, тем самым максимизируя p-value

Ради интереса увеличим количество значений в выборке и посмотрим, как это все работает.

plot_pdfs( *fit_params( ss.gamma(1.3, 10, 400).rvs(300,random_state=42) ), N=300 )

Параметры с максимизацией правдоподобия (MLE): (1.38, 12.85, 363.73).
KS p-value MLE: 0.994.

Параметры методом моментов (MM): (1.82, −47.59, 309.99).
KS p-value MM: 0.880.

Хорошо видно, что MLE лучше обучился на данных, чем MM. И в обоих случаях значения p-value завышены.

Критерии Крамера — Мизеса — Смирнова и Андерсона — Дарлинга строятся почти как критерий Колмогорова — Смирнова. На вход подаются теоретическая и эмпирическая функции распределения. Только статистика имеет более сложный вид, чем супремум модуля разности.

Тест Андерсона — Дарлинга — это тот же тест Мизеса, только хвостам функции распределения присвоены большие веса. В Scipy для них реализованы функции stats.cramervonmises и stats.anderson. Они также не свободны от распределения для сложных гипотез.

Обычно мощности критерия соотносятся как Андерсон — Дарлинг → Мизес → Колмогоров — Смирнов.

Тесты типа

Если предыдущие критерии работали с функцией распределения (cumulative distribution function), то критерии $\chi2$ работают с частотами.

Критерий $\chi2$ сравнивает вероятность попадания значения в определенный интервал. Грубо говоря, сравнивает теоретическую и эмпирическую гистограммы.

Нужно сгруппировать данные по бинам, чтобы получить эмпирическую гистограмму. Это можно сделать разными способами. В отличие от получения ecdf — эмпирической функции распределения, которая всегда получается только одним способом. Поэтому такие тесты называются тестами на группированных данных.

Для дискретной случайной величины (например, для Пуассона) ничего дополнительно группировать не надо. Для расчета критерия находится расстояние $\chi2$ как:

$\chi^2 = \sum \frac{(O-E)^2}{E}$

Поделив числитель и знаменатель на 1/N^2 , получим выражение для вероятностей:

$\chi^2 = N\sum \frac{(p_{O}-p_{E})^2}{p_{E}}$

Эта величина имеет статистику $\chi2$ c степенями свободы (n — количество бинов, на которые разбиваются данные). Зная эти две величины, мы можем сказать, какова вероятность получить конкретные данные из конкретного распределения, то есть посчитать p-value.

Группировка данных — это отдельная магия. Тут решается классическая задача точности и разброса по типу принципа неопределенности. Есть способы сделать разбиение, минимизируя потери в информации Фишера, но останавливаться на этом мы не будем.

В литературе говорится, что при оптимальной группировке такие тесты обладают большей мощностью, чем тесты на негруппированных данных.

$\chi2$ — самый известный, но далеко не единственный тест по группированным данным. Рассмотрим еще несколько: G-тест, Фримена — Тьюки и другие.

По группированным данным придумана еще масса статистик. Например, в G-тесте статистика считается так:

$G = 2\sum O_i*ln\bigg(\frac{O_i}{E_i}\bigg)$

Если отношение близко к единице, логарифм можно разложить в ряд Тейлора и получить тест $\chi2$ . Получается, G-тест — это обобщение $\chi2$ , когда наблюдаемые и ожидаемые частоты не близки. Поэтому он более устойчив к выбросам.

Статистика Фримена — Тьюки (Freeman-Tukey) считается так:

$T^2 = 4*\sum \big( \sqrt{O_i}-\sqrt{E_i} \big)^2$

Информации, в каких конкретно ситуациях она лучше, я не нашел, но, наверное, где-то предпочтительнее.

Все статистики по группированным данным, которые мы описывали, обобщаются статистикой Кресси — Рида (Cressie-Read), где, подставляя различные значения степени, получаем и $\chi2$ , и G, и Фримена — Тьюки, и другие статистики. В Scipy Кресси — Рида есть функция scipy.stats.power_divergence, где параметром lambda_ можно выбрать один из шести стат-тестов. Если залезть внутрь Scipy, мы увидим там вызов power_divergence.

Возьмем нашу выборку и потестируем ее как простую и сложные гипотезы. В идеале значения p-value в обоих случаях должны быть близки.

def get_pvalues_grupp(dist, params: list, lenght: int, f_obs: np.array):

    pr = dist.cdf(bins,*params)
    probabilites = pr[1:] - pr[:-1]
    f_exp = probabilites * lenght / probabilites.sum()

    for lambda_ in ('pearson','log-likelihood','cressie-read'):
        pvalue = ss.power_divergence(f_obs,f_exp,lambda_=lambda_).pvalue
        print(lambda_+' pvalue: {:.4f}'.format(pvalue) )

    return f_exp


rvs = ss.gamma(1.3, 10, 400).rvs(300,random_state=42)

hist, bins = np.histogram(rvs,bins='auto',density=True)
x_pdf = np.linspace(0,max(rvs),500)
f_obs, bins = np.histogram(rvs,bins='auto')

params_mle, params_mm = fit_params( rvs )

print('\n===Тестируем простую гипотезу===')
get_pvalues_grupp(ss.gamma,(1.3, 10, 400), len(rvs), f_obs);
print('\n===Тестируем сложную гипотезу===')
f_exp = get_pvalues_grupp(ss.gamma,params_mle, len(rvs), f_obs)

''' Рисуем гистограмму и расчитанные pdf '''

fig = make_subplots(rows=1, cols=2
,subplot_titles=("ePDF и PDF", "Ожидаемые и наблюдаемые частоты"))
fig.add_trace(go.Scatter(x=x_pdf, y=ss.gamma(1.3, 10, 400).pdf(x_pdf)
                                                , name='Pdf'), row=1, col=1)
fig.add_trace(go.Bar(x=bins, y=hist, name='Histplot'), row=1, col=1)

fig.add_trace(go.Bar(x=bins, y=f_exp, name='Expected freq'), row=1, col=2)
fig.add_trace(go.Bar(x=bins, y=f_obs, name='Observed freq'), row=1, col=2)
fig.show()

>===Тестируем простую гипотезу=== 
>pearson pvalue: 0.5358
>log-likelihood pvalue: 0.2460 
>cressie-read pvalue: 0.4672  
>
>===Тестируем сложную гипотезу=== 
>pearson pvalue: 0.4786 
>log-likelihood pvalue: 0.2600
>cressie-read pvalue: 0.4334

Метод Монте-Карло

Мы можем вычислить выборочные характеристики, например среднее и дисперсию, для нормального распределения по экспериментальной выборке. Потом нагенерировать выборки той же длины и подсчитать на них нужные статистики, не зная, как эти статистики распределены на самом деле (ту же статистику Колмогорова — Смирнова).

Можно провести много численных экспериментов и посмотреть на получившиеся гистограммы интересующих нас статистик. Давайте сделаем это.

def ks_stat(data
        	,cdf
        	,params):

	x = np.sort(data)
	cdfval = cdf(x,*params)
	n = len(x)
	Dpls = max( np.arange(1.,n+1)/n - cdfval )
	Dmin = max( cdfval - np.arange(0.,n)/n )

	return max([Dmin,Dpls])

def get_ecdf(t):
	# Возвращает эмпирическую функцию распределения и значения переменной
	# x, ecdf
	x, cnt = np.unique(t, return_counts=True)
	cumsum = np.cumsum( cnt )
	return x, cumsum/cumsum[-1]

def bad_monte_carlo(dist
                	,data
                	,params
                	,size=999):

	S = np.zeros(size)
	lendata = len(data)
	# Посчитали статистику между исходными данными и сфитированным распределением
	S_obs = ks_stat(data,dist.cdf,params)

	for i in range(size):
    	rvs = dist.rvs(*params,size=lendata)
    	S[i] = ks_stat(rvs,dist.cdf,params)

	print('Monte-Carlo P-value: {:.3f}'.format( (S >= S_obs).sum()/size ) )

def _anderson_darling(dist, data):
	x = np.sort(data, axis=-1)
	n = data.shape[-1]
	i = np.arange(1, n+1)
	Si = (2*i - 1)/n * (dist.logcdf(x) + dist.logsf(x[..., ::-1]))
	S = np.sum(Si, axis=-1)
	return -n - S

def _cramer_von_mises(dist, data):
	x = np.sort(data, axis=-1)
	n = data.shape[-1]
	cdfvals = dist.cdf(x)
	u = (2*np.arange(1, n+1) - 1)/(2*n)
	w = 1 / (12*n) + np.sum((u - cdfvals)**2, axis=-1)
	return w
# Давайте проверим простую гипотезу
data = ss.gamma(1.3, 10, 400).rvs(100,random_state=42)
print('===Простая гипотеза===')
bad_monte_carlo(ss.gamma
            	,data
            	,(1.3, 10, 400))

kstestp = ss.ks_1samp( data
                  	,ss.gamma.cdf
                  	,(1.3, 10, 400) ).pvalue

print('KS test: {:.3f}'.format( kstestp ) ) 
# Ну и сложную, которую раньше тестировали
params_mm  = ss.gamma.fit(data, method='MM')

print('\n===Сложная гипотеза===')
bad_monte_carlo(ss.gamma
            	,data
            	,params_mm )

kstestp = ss.ks_1samp(data
                  	,ss.gamma.cdf
                  	,params_mm ).pvalue

print('KS test: {:.3f}'.format( kstestp ) )
>===Простая гипотеза===
>Monte-Carlo P-value: 0.643
>KS test: 0.649
>
>===Сложная гипотеза===
>Monte-Carlo P-value: 0.926
>KS test: 0.928

Получили те же числа, только считать дольше? Получается, нет никакого резона в Монте-Карло?

Но мы неправильно его реализовали и сравнивали каждый раз с тем же самым распределением. А, как мы видели выше, это завышает p-value.

Попробуем такой алгоритм:

1. Есть экспериментальная выборка — data. По ней находим вектор параметров распределения params.

2. На каждом шаге генерируем rvs — случайную выборку того же размера.

3. Находим _params, обучаясь на этой новой выборке.

4. Находим статистику T(_params,rvs), то есть смотрим, как себя будут вести статистики на сгенерированных выборках.

Посмотрим на наш пример:

def good_monte_carlo(dist
                	,data
                	,params
                	,size=999
                	,method='MLE'
                	,statistic='ks'
                	,print_pvalue=False):
	'''dit — объект scipy stats.name_of_dist, data — сгенерированная выборка,
	params — параметры распределения (истинные — тогда тестируется простая
	гипотеза, посчитанные на выборке — тогда сложная)
	size — количество симуляций Монте-Карло'''
	S = np.zeros(size)
	lendata = len(data)
	# Посчитали статистику между исходными данными и сфитированным распределением
	if statistic == 'ks':
    	S_obs = ks_stat(data,dist.cdf,params)
	elif statistic == 'ad':
    	S_obs = _anderson_darling(dist(*params),data)
	elif statistic == 'cvm':
    	S_obs = _cramer_von_mises(dist(*params),data)

	for i in range(size):
    	# Сгенерировали случайную выборку исходного размера
    	rvs = dist.rvs(*params,size=lendata)
    	# Нашли параметры распределения
    	_params = dist.fit(rvs,method=method)
    	# Посчитали статистику
    	if statistic == 'ks':
        	S[i] = ks_stat(rvs,dist.cdf,_params)
    	elif statistic == 'ad':
        	S[i] = _anderson_darling(dist(*_params),rvs)
    	elif statistic == 'cvm':
        	S[i] = _cramer_von_mises(dist(*_params),rvs)

	pvalue =  (S >= S_obs).sum() / size
	if print_pvalue:
    	print('Monte-Carlo P-value: {:.3f}'.format( pvalue ) )
	return pvalue
# Давайте проверим простую гипотезу
data = ss.gamma(1.3, 10, 400).rvs(100,random_state=42)
print('===Простая гипотеза===')
good_monte_carlo(ss.gamma
            	,data
            	,(1.3, 10, 400)
            	,size=9999
            	,method='MM'
            	,statistic='ks'
            	,print_pvalue=True)
kstestp = ss.ks_1samp( data ,ss.gamma.cdf,(1.3, 10, 400) ).pvalue
print('KS test: {:.3f}'.format( kstestp ) ) 

# Ну и сложную, которую раньше тестировали
params_mm  = ss.gamma.fit(data, method='MM')
print('\n===Сложная гипотеза===')
good_monte_carlo(ss.gamma
            	,data
            	,params_mm
            	,method='MM'
            	,size=9999
            	,statistic='ks'
            	,print_pvalue=True)
kstestp = ss.ks_1samp(data ,ss.gamma.cdf,params_mm ).pvalue
print('KS test: {:.3f}'.format( kstestp ) )

>===Простая гипотеза===
>Monte-Carlo P-value: 0.522
>KS test: 0.649
>
>===Сложная гипотеза===
>Monte-Carlo P-value: 0.879
>KS test: 0.928

P-value уменьшилось, хотя это и не поражает воображение. Метод реализован в Scipy функцией stats.goodness_of_fit.

Какое у нас все-таки распределение

Попробуем сравнить некоторые наши критерии. Распределение выше примерно соответствует распределению среднего времени диалога в одном из интервалов. Вопрос: это гамма-распределение или, например, логнормальное?

Посмотрим, насколько легко наши критерии отличат один от другого. Сгенерируем несколько выборок и посмотрим на распределение p-value. По умолчанию в Scipy используется MLE-метод.

import plotly.express as px

def statsign(basedist
        	,distname1:str
        	,distname2:str
        	,method1='MLE'
        	,method2='MLE'
        	,N=100
        	,l=300
        	,testtype='simple_ks'
        	,mkparams={} ) -> pd.DataFrame:
	''' basedist — базовое распределение, из которого генерируется выборка
	distname1 — первое распределение для сравнения, method1 — метод фитинга 1
	distname2 — второе распределение для сравнения, method2 — метод фитинга 2
	N — количество элементов в выборке
	l — количество выборок
	testtype: simple_ks — простой К-С тест, mk — Монте-Карло
	mkparams — Монте-Карло-словарь
	'''
	results = np.zeros([l,2])
	dist1 = getattr(ss,distname1)
	dist2 = getattr(ss,distname2)

	for i in tqdm( range(l) ):
    	data = basedist.rvs(N)

    	params1 = dist1.fit(data, method=method1)
    	params2 = dist2.fit(data, method=method2)

    	if testtype == 'simple_ks':
        	results[i,0] = ss.ks_1samp( data ,dist1.cdf,params1 ).pvalue
        	results[i,1] = ss.ks_1samp( data ,dist2.cdf,params2 ).pvalue
    	elif testtype == 'mk':
        	results[i,0] = good_monte_carlo(dist1,data,params1,method=method1,**mkparams)
        	results[i,1] = good_monte_carlo(dist2,data,params2,method=method2,**mkparams)

	return pd.DataFrame({distname1:results[:,0],distname2:results[:,1]})

def plot_hist(pvalues:pd.DataFrame) -> None:

	fig = px.histogram(
    	pvalues,
    	title='Распределение p-value',
    	histnorm='probability',
    	barmode="overlay",
    	nbins = 20 )
	fig.show()
# Сохраним в словарике разные p-value для разных тестов
pvalues_tests = {}
pvalues = statsign(ss.gamma(1.3, 10, 400)
               	,'gamma'
               	,'lognorm'
               	,N=100
               	,l=300)

pvalues = pvalues[pvalues.gamma >0.01] # Исключили неправильный фит
pvalues_tests['ks'] = pvalues
delta = pvalues.gamma - pvalues.lognorm > 0
print( delta.value_counts()/len(delta)  )

>True 	0.560166
>False	0.439834

Мы сгенерировали 300 выборок из гамма-распределения по 100 элементов. Посмотрели, в каком количестве пар p-value гамма-распределения было больше, чем p-value логнормального (значение True). При этом пары, где p-value гамма-распределения было меньше 0,01, выкинули, так как это говорит о некорректной работе MLE-оптимизатора.

Получили примерно одинаковое количество ложных и истинных срабатываний. Причем в большом проценте экспериментов p-value логнормального распределения было в районе 0,9, см. гистограмму.

Теперь используем Монте-Карло-симуляцию в тех же условиях.

pvalues = statsign(ss.gamma(1.3, 10, 400)
               	,'gamma'
               	,'lognorm'
               	,N=100
               	,testtype='mk'
               	,l=300)

pvalues = pvalues[pvalues.gamma > 0.01]
pvalues_tests['monte-karlo'] = pvalues

delta = pvalues.gamma - pvalues.lognorm > 0
print( delta.value_counts()/len(delta)  )
plot_hist(pvalues)

>True 	0.691057
>False	0.308943

При использовании Монте-Карло алгоритм ошибается куда в меньшем количестве случаев в тех же условиях (70 на 30) при том же критерии. Понятно, что нужно получить p-value на данное p-value и, если мы еще раз запустим процедуру, результат будет несколько отличатся. Но общая картина ясна.

Заметим, что с ростом выборки все становится заметно лучше и уже при N = 600 в 95% случаев мы не ошибемся, если выберем распределение с большим p-value без всякого Монте-Карло, просто используя тест Колмогорова — Смирнова как тест на простую гипотезу.

pvalues = statsign(ss.gamma(1.3, 10, 400)
               	,'gamma'
               	,'lognorm'
               	,N=600
               	,l=10000)

pvalues = pvalues[pvalues.gamma >0.01] # Исключили неправильный фит
delta = pvalues.gamma - pvalues.lognorm > 0
print( delta.value_counts()/len(delta)  )

plot_hist(pvalues)

True 	0.944528
False	0.055472

Другими словами, если у нас нет способа увеличить выборку, нужно использовать правильное Монте-Карло. Посмотреть на разные статистики и выбрать ту, которая делает критерий более мощным.

Выводы

Нет единого аналитического подхода для всех объемов выборок и всех распределений. Самый очевидный вариант «из коробки» может легко обмануть.

Для каждой конкретной задачи лучше провести моделирование, насколько мощны ваши критерии в ней. Какой критерий лучше выбрать и какой объем выборки необходим.

Я на этом не останавливался, но существуют специализированные тесты для сложных гипотез, например тест Лиллиефорса для нормального распределения. Если есть возможность, лучше пользоваться ими.

Хабы: