Анализ гипотез и очень, ооооооочень странные дела
Лет 5 назад я усиленно пытался вникнуть в тервер и статы: книги, статьи, вебсёрфинг. Даже написал несколько статей: раз, два, три. Вообще, в планах было написать довольно большой цикл статей, что бы подсветить какие-то самые сложные вещи, да и самому в них разобраться - совместить полезное с полезным, так сказать. Однако, в какой-то момент я решил, что полученных знаний достаточно для новых проектов и ушел в работу. Работал. Работал. Работал.
Но знаете, все это время меня никак не покидало ощущение, что в процессе своего знакомства с наукой о случайностях я наткнулся на какой-то скелет в шкафу. И знаете что? Мне кажется настало время вытащить этот скелет на всеобщее обозрение.
Есть в тестировании гипотез что-то подозрительное
Почему проверка гипотез, в том числе и статистическая, так важна? Все просто - в конечном счете, качество нашей жизни прямо пропорционально качеству и количеству наших гипотез которые мы можем проверять. Четкие, обоснованные и проверяемые допущения ведут к эффективным решениям, а размытые и ложные к череде ошибок и разочарований. Не формулировать никаких гипотез - это тоже гипотеза, которая состоит в молчаливом и безрассудном допущении, что всё останется по-старому и любые наши действия не будут иметь никаких последствий. Не формулировать никаких гипотез - это самая рискованная ставка из всех возможных.
Гиганты, на плечах которых мы стоим понимали это. Основная проблема связанная с проверкой статистических гипотез, заключается в установке некоторых правил или точнее критериев в соответствии с которыми выдвинутая гипотеза принимается или отвергается. Такие критерии впервые были сформулированы Джоном фон Нейманом и Эгоном Пирсоном. Давайте кратко рассмотрим основные аспекты теории, которая стоит за этим критерием.
Допустим у нас есть монетка в честности которой мы хотим убедиться. Мы подкинули монетку 100 раз и у нас выпало 57 орлов и 43 решки. Основная гипотеза -
У нас есть 100 наблюдений и очевидно, что если мы подбросим ту же самую монетку еще 100 раз, то результат может оказаться другим. Но разные результаты окажутся более или менее вероятными. Что бы принять гипотезу
; .
В результате мы увидим следующее:
Код для отрисовки
import numpy as np
np.random.seed(999)
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
from scipy.stats import bernoulli, binom
N = 100
n_heads, n_tails = 57, 43
p_h0, p_h1 = 0.5, 0.65
plt.figure(figsize=(12, 6))
k_heads = np.arange(35, 67)
p_k_heads_h0 = binom.pmf(k_heads, N, p_h0)
plt.plot(k_heads, p_k_heads_h0, 'o', c='darkred', ms=4, zorder=11, label=r'$H_{0}:\;p=0.5$')
mask = np.cumsum(p_k_heads_h0) < 0.95
plt.vlines(k_heads[mask], 0, p_k_heads_h0[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h0[~mask], color='r', lw=3, zorder=10,
label=r'$\alpha$' + f' = {p_k_heads_h0[~mask].sum():.2f}')
k_heads = np.arange(50, 80)
p_k_heads_h1 = binom.pmf(k_heads, N, p_h1)
plt.plot(k_heads, p_k_heads_h1, 'o', c='darkblue', ms=4, zorder=11, label=r'$H_{1}:\;p=0.65$')
mask = np.cumsum(p_k_heads_h1) > 0.06
plt.vlines(k_heads[mask], 0, p_k_heads_h1[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h1[~mask], color='b', lw=3, zorder=10,
label=r'$\beta$' + f' = {p_k_heads_h1[~mask].sum():.2f}')
plt.axvline(58, c='k', ls='--', zorder=15, lw=1, label='k = 58')
plt.legend()
plt.xlim(30, 85)
plt.xlabel('Число орлов (X)')
plt.ylabel('Вероятность $P(X=x)$')
plt.title('Критерий Неймана–Пирсона: Ошибки I и II рода')
plt.show()На графике видно, что функции перекрываются и очевидно, что ошибки связанные с принятием той или иной гипотезы находятся именно в этой области перекрытия. Итак, что бы принять
В тоже самое время, когда мы определили критическую область мы так же создали возможность для появления еще одной ошибки. Если, на самом деле монетка нечестная, т.е. на самом деле верна
Из вышеприведенного графика так же следует, что если для некоторой критической области мы получаем наименьшие значения
Выбор критической области позволяет сделать сколь угодно малыми либо
Код для отрисовки
plt.figure(figsize=(12, 6))
N = 100
k_heads = np.arange(35, 67)
p_k_heads_h0 = binom.pmf(k_heads, N, p_h0)
plt.plot(k_heads, p_k_heads_h0, 'o', c='darkred', ms=4, zorder=11, label=r'$H_{0}:\;p=0.5$')
mask = np.cumsum(p_k_heads_h0) < 0.99
plt.vlines(k_heads[mask], 0, p_k_heads_h0[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h0[~mask], color='r', lw=3, zorder=10,
label=r'$\alpha$' + f' = {p_k_heads_h0[~mask].sum():.2f}')
k_heads = np.arange(50, 80)
p_k_heads_h1 = binom.pmf(k_heads, N, p_h1)
plt.plot(k_heads, p_k_heads_h1, 'o', c='darkblue', ms=4, zorder=11, label=r'$H_{1}:\;p=0.65$')
mask = np.cumsum(p_k_heads_h1) > 0.23
plt.vlines(k_heads[mask], 0, p_k_heads_h1[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h1[~mask], color='b', lw=3, zorder=10,
label=r'$\beta$' + f' = {p_k_heads_h1[~mask].sum():.2f}')
plt.axvline(62, c='k', ls='--', zorder=15, lw=1, label='k = 62')
plt.legend()
plt.xlim(30, 85)
plt.xlabel('Число орлов (X)')
plt.ylabel('Вероятность $P(X=x)$')
plt.title('Критерий Неймана–Пирсона: Ошибки I и II рода')
plt.show()При фиксированном объеме выборки мы не можем одновременно уменьшить и
Представленный пример с монеткой очень прост. В общем случае гипотезы могут быть крайне сложными, а функции распределения вероятностей не всегда возможно задать аналитически. Глядя на графики, можно подумать, что нам достаточно:
выбрать приемлемое для нас значение
; дальше, исходя из выражения
, найти значение - границу критической области; ну а значение
появится как-бы само собой.
В общем случае это не так. Очень часто для заданного уровня
Если следовать этому принципу, то при фиксированном объеме выборки вероятность
Выбор объема выборки и уровня
Кажется, что на этом все. Однако, теория Неймана-Пирсона несколько глубже чем можно подумать.
В определенной степени принято считать, что в математической статистике есть два лагеря: лагерь фреквентистов - сторонников частотных методов в статистике и лагерь байесовистов - сторонников байесовских методов вывода. Причем, вторые, как правило, довольно пренебрежительно отзываются о первых (иногда даже слишком). Например, Кэмерон Дэвидсон-Пайлон в своей книге "Вероятностное программирование на Python: байесовский вывод и алгоритмы." пишет о фреквентизме и A/B-тестеровании следующее:
Зачастую анализ после эксперимента производится с помощью так называемой проверки статистической гипотезы, например проверки различия средних значений (difference of means test) или проверки различия долей (difference of proportions test). Это требует использования такого запутанного понятия, как «Z-оценка», и еще более запутанных малопонятных «p-значений» (даже не спрашивайте, что это такое). Если вы слушали курс статистики, то вас, возможно, обучали этой методике (хотя вовсе не факт, что вы в ней разобрались). И если мы с вами в этом схожи, то их вывод отнюдь не доставляет вам радости. Если так — отлично. Байесовский подход к решению этой задачи гораздо более естественен. (Дэвидсон-Пайлон Кэмерон. Вероятностное программирование на Python: байесовский вывод и алгоритмы. — СПб.: Питер, 2019. с. 62).
Но такое отношение байесовистов к фреквентистам (с определенной натяжкой) можно понять. Наиболее известным и резким критиком байесовистов был Рональд Фишер. Его авторитет в первой половине XX века был настолько велик, что его скептицизм во многом затормозил развитие байесовских методов. А его риторика, часто выходившая за рамки научной дискуссии и задала тон для целых поколений статистиков-фреквентистов. Прочувствовать этот тон можно с помощью нижеследующего списка "критических" тезисов:
Байесовские методы - это просто формализация предубеждений. Вы получаете на выходе то, во что изначально верили. Это не наука, а самооправдание. Настоящая наука должна основываться на объективных данных, а не на чьих-то личных верованиях.
Байесовство - это не статистика, это вера. Вы, байесовцы, просто верите в свой априор, а потом находите данные, которые подтверждают вашу веру. Это слепая религия, а не наука.
Байесовские методы — это интересная академическая игрушка для простых задач с игральными костями, урнами с шарами и монетками. Но для реальной науки они бесполезны и непрактичны.
Если вы не получили желаемый результат, вы всегда можете поиграть с априорным распределением, пока не получите его. Это легализованный p-hacking на стероидах.
Сейчас накал страстей заметно поутих по сравнению с тем, что было раньше. Теперь разумное мнение состоит в том, что обе парадигмы дополняют друг друга, а не противоречат. Однако, отголоски этих "статистических войн" можно услышать до сих пор.
Зачем такое длинное отступление?
Противоборство фреквентистов и байесовистов, это яркий пример того, что даже в научном сообществе могут царить разногласия, инерционность и консерватизм. Можно ли утверждать, что сейчас такого нет? Помните в самом начале я говорил о скелете в шкафу? Так вот, на самом деле я даже не уверен в том - скелет ли это (потому что долгие годы вся окружающая действительность убеждала меня в том, что это не так). Но давайте не будем забегать наперед. Мы еще не закончили с теорией Неймана-Пирсона.
Давайте снова вернемся к нашей монетке, но на этот раз взглянем на результат с точки зрения байесовизма. Запишем теорему Байеса для наших гипотез в пропорциональной форме:
Пропорциональная форма теоремы Байеса удобна тем, что мы можем ее использовать для сравнения двух гипотез. Для этого достаточно просто разделить один апостериор на другой:
Если это отношение равно 3, то можно сказать, что гипотеза
Отношение апостериорных вероятностей состоит из двух частей:
и
Очевидно, что если гипотеза
Код для отрисовки
def L_ratio(x, n, p_h0, p_h1):
ratio = binom.pmf(x, n, p=p_h0) / binom.pmf(x, n, p=p_h1)
return ratio
f, ax = plt.subplots(1, 2, figsize=(12, 5))
for j in range(5):
samples = bernoulli.rvs(p=p_h0, size=N)
x = np.cumsum(samples)
n = np.r_[1:101]
ax[0].plot(n, x)
ax[0].axhline(58, color='k', ls='--', label='k = 58' if j==0 else None)
L_x = np.log([L_ratio(x[i], n[i], p_h0, p_h1) for i in range(100)])
ax[1].plot(n, L_x)
ax[1].axhline(np.log(58), color='k', ls='--', label='k = ln(58)' if j==0 else None)
ax[0].legend()
ax[1].legend()
ax[0].set_xlabel('Количество бросков (n)')
ax[0].set_ylabel('Число орлов (x)')
ax[0].set_title('Изменение количества орлов в выборке')
ax[1].set_xlabel('Количество бросков (n)')
ax[1].set_ylabel(r'$\ln(\Lambda$)')
ax[1].set_title(r'Изменение $\Lambda(H_{0} : H_{1} \; | \; X=x)$ (верна $H_0$)')
plt.show()Обратите внимание, что чем лучше
Если мы взглянем на график изменения отношения правдоподобия при верной гипотезе
Код для отрисовки
f, ax = plt.subplots(1, 2, figsize=(12, 5))
for j in range(5):
samples = bernoulli.rvs(p=p_h1, size=N)
x = np.cumsum(samples)
n = np.r_[1:101]
ax[0].plot(n, x)
ax[0].axhline(58, color='k', ls='--', label='k = 58' if j==0 else None)
L_x = np.log([L_ratio(x[i], n[i], p_h0, p_h1) for i in range(100)])
ax[1].plot(n, L_x)
ax[1].axhline(np.log(1/58), color='k', ls='--', label='k = ln(1/58)' if j==0 else None)
ax[0].legend()
ax[1].legend()
ax[0].set_xlabel('Количество бросков (n)')
ax[0].set_ylabel('Число орлов (x)')
ax[0].set_title('Изменение количества орлов в выборке')
ax[1].set_xlabel('Количество бросков (n)')
ax[1].set_ylabel(r'$\ln(\Lambda$)')
ax[1].set_title(r'Изменение $\Lambda(H_{0} : H_{1} \; | \; X=x)$ (верна $H_1$)')
plt.show()На самом деле Нейман и Пирсон показали, что область состоящая из всех возможных выборок
является наиболее мощной критической областью для проверки гипотезы
Интересно, что писал Сам Джон фон Нейман об этом открытии:
I can point to the particular moment when I understood how to formulate the undogmatic problem of the most powerful test of a simple statistical hypothesis against a fixed simple alternative. At the present time [probably 1968], the problem appears entirely trivial and within easy reach of a beginning undergraduate. But, with a degree of embarrassment, I must confess that it took something like half a decade of combined effort of E. S. P. [Egon Pearson] and myself to put things straight.» (источник).
То есть то, что сейчас кажется тривиальным стоило величайшим ученым немалых усилий и времени. Отсюда возникает новый вопрос: может ли что-то невероятно крутое и ценное быть прямо под носом на протяжении многих лет, а ты об этом даже не подозреваешь?
Еще любопытно то, что Нейман упоминал, что опубликованная статья с полученными результатами получила "неожиданно" благоприятную рецензию от Рональда Фишера. Слово "неожиданно" использовалось в упоминаниях видимо потому, что фундаментом леммы Неймана-Пирсона является теорема Байееса, которую Фишер так беспощадно критиковал. Но почему? Может потому что он наконец признал право байесовизма на жизнь? Скорее всего нет. Вероятнее всего к тому моменту он уже считал фреквентизм неуязвимым, считал его состоявшейся нормой. Стандартом.
Чтож, мы опять немного отвлеклись. Давайте снова вернемся к лемме Неймана-Пирсона и еще раз взглянем на то как меняется отношение правдоподобий при броске монетки.
Код для отрисовки
f, ax = plt.subplots(1, 2, figsize=(12, 5))
for j in range(10):
samples = bernoulli.rvs(p=p_h0, size=N)
x = np.cumsum(samples)
n = np.r_[1:101]
L_x = np.log([L_ratio(x[i], n[i], p_h0, p_h1) for i in range(100)])
ax[0].plot(n, L_x, c='0.7')
ax[0].plot(n[-1], L_x[-1], 'o')
ax[0].axhline(np.log(58), color='k', ls='--', label='k = ln(58)' if j==0 else None)
samples = bernoulli.rvs(p=p_h1, size=N)
x = np.cumsum(samples)
n = np.r_[1:101]
L_x = np.log([L_ratio(x[i], n[i], p_h0, p_h1) for i in range(100)])
ax[1].plot(n, L_x, c='0.7')
ax[1].plot(n[-1], L_x[-1], 'o')
ax[1].axhline(np.log(1/58), color='k', ls='--', label='k = ln(1/58)' if j==0 else None)
ax[0].legend()
ax[1].legend()
ax[0].set_xlabel('Количество бросков (n)')
ax[0].set_ylabel(r'$\ln(\Lambda$)')
ax[0].set_title(r'Изменение $\Lambda(H_{0} : H_{1} \; | \; X=x)$ (верна $H_0$)')
ax[1].set_xlabel('Количество бросков (n)')
ax[1].set_ylabel(r'$\ln(\Lambda$)')
ax[1].set_title(r'Изменение $\Lambda(H_{0} : H_{1} \; | \; X=x)$ (верна $H_1$)')
plt.show()Из-за логарифмирования
В лемме Неймана-Пирсона ничего не сказано о размере выборки, но формально она выглядит так:
где
Код для отрисовки
n = np.arange(75, 200)
k = binom.ppf(0.95, n, p_h0)
alpha = binom.sf(k, n, p_h0)
beta = binom.cdf(k, n, p_h1)
plt.plot(n, alpha, label=r'$\alpha$')
plt.plot(n, beta, label=r'$\beta$')
plt.legend()
plt.xlabel('Число бросков (n)')
plt.ylabel('Вероятность')
plt.title(r'Изменение $\beta$ с ростом числа бросков')
plt.show()То есть количество бросков жестко задано уровнем
Но согласитесь, после увиденного все это кажется немного нелогичным. Что бы это продемонстрировать можно привести простой пример. Допустим, мы тестируем гипотезу которая состоит в том, что монетка либо абсолютно честная
Но он показывает, что множество всех выборок с разным объемом
Критическая область для гипотезы
, например, все векторы длинны , состоящие из одних единиц. Критическая область для гипотезы
, например, любой вектор длинны , который начинается с единиц и заканчивается одним нулем. Допустимая область для обеих гипотез, например, все векторы длинны
, состоящие из одних единиц.
По крайней мере, хотя бы в этом примере мы не обязаны выполнять фиксированное и заданное наперед количество бросков, а прервать эксперимент на любом шаге
Выше мы видели, что отношение правдоподобий
когда справедлива гипотеза
Так же, выше мы видели, что если
где
Для удобства прологарифмируем
и обозначим i-й член последовательности через
Для нашего примера с монеткой, мы можем записать следующее:
если выпал орел, и
если выпала решка. Тогда правило по которому определяются три интересующих нас области можно записать вот так
где
Единственная серьезная проблема для нашего правила состоит в том, что мы не знаем ка определять необходимые значения
и критическая область для
Тогда для любой выборки
но конкретно для последнего шага
то вероятность получения любой такой выборки будет по крайней мере в
которое можно переписать как
С помощью аналогичных рассуждений мы можем прийти к тому, что
Выходит, что теперь у нас есть критерий, который больше не требует фиксированных объемов выборок? Давайте проверим на нашей монетке.
Для тестирования по правилу (далее будем называть критерием) Неймана-Пирсона, нам необходимо задать количество бросков и необходимый уровень
Код для отрисовки
plt.figure(figsize=(12, 6))
k_heads = np.arange(35, 67)
p_k_heads_h0 = binom.pmf(k_heads, N, p_h0)
plt.plot(k_heads, p_k_heads_h0, 'o', c='darkred', ms=4, zorder=11, label=r'$H_{0}:\;p=0.5$')
mask = np.cumsum(p_k_heads_h0) < 0.955
plt.vlines(k_heads[mask], 0, p_k_heads_h0[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h0[~mask], color='r', lw=3, zorder=10,
label=r'$\alpha$' + f' = {p_k_heads_h0[~mask].sum():.3f}')
k_heads = np.arange(50, 80)
p_k_heads_h1 = binom.pmf(k_heads, N, p_h1)
plt.plot(k_heads, p_k_heads_h1, 'o', c='darkblue', ms=4, zorder=11, label=r'$H_{1}:\;p=0.65$')
mask = np.cumsum(p_k_heads_h1) > 0.087
plt.vlines(k_heads[mask], 0, p_k_heads_h1[mask], color='0.5', lw=1, zorder=9)
plt.vlines(k_heads[~mask], 0, p_k_heads_h1[~mask], color='b', lw=3, zorder=10,
label=r'$\beta$' + f' = {p_k_heads_h1[~mask].sum():.3f}')
plt.axvline(59, c='k', ls='--', zorder=15, lw=1, label='k = 59')
plt.legend()
plt.xlim(30, 85)
plt.xlabel('Число орлов (X)')
plt.ylabel('$P(X=x)$')
plt.title('Критерий Неймана–Пирсона: Ошибки I и II рода')
plt.show()Итак, если в результате сотни подбрасываний монетки количество орлов превосходит 58, то мы считаем, что монетка нечестная. Теперь посмотрим получится ли у нас прийти к таким же выводам при произвольном количестве подбрасываний, с таким же уровнем и мощностью.
Первым делом вычислим значения для границ
Ну а дальше мы просто повторим тест 200 раз и будем следить, когда отношение вероятностей пересечет ту или иную границу.
Код для отрисовки
plt.figure(figsize=(12, 5))
A, B = np.log([(1 - 0.087)/0.044, 0.087/(1 - 0.044)])
for _ in range(200):
samples = bernoulli.rvs(p=p_h0, size=500)
m = np.arange(1, 501)
m_star = np.cumsum(samples)
z = m_star * np.log(p_h1 / p_h0) + (m - m_star) * np.log((1 - p_h1) / (1 - p_h0))
n_step = np.where((B > z) | (z > A))[0][0] + 1
plt.plot(m[:n_step], z[:n_step], c='C0', alpha=0.3)
plt.plot(m[:n_step][-1], z[:n_step][-1], 'C0o', alpha=0.5)
plt.axhline(A, color='darkred', label=f'A = {A:.2f}')
plt.axhline(B, color='darkgreen', label=f'B = {B:.2f}')
plt.axvline(N, color='darkblue', label=r'$n^{\mathrm{(NP)}}$ = 100')
plt.legend(loc='upper right', framealpha=1)
plt.xlabel('Количество бросков (n)')
plt.ylabel(r'$z_{i}$')
plt.title(r'Последовательный критерий (верна $H_{0}$, $\alpha=0.044,\;\beta=0.087$)')
plt.show()Чтож, как видим большинство остановок теста действительно совершаются задолго до 100-го броска монеты. Тоже самое справедливо и для верной гипотезы
Код для отрисовки
plt.figure(figsize=(12, 5))
A, B = np.log([(1 - 0.087)/0.044, 0.087/(1 - 0.044)])
for _ in range(200):
samples = bernoulli.rvs(p=p_h1, size=500)
m = np.arange(1, 501)
m_star = np.cumsum(samples)
z = m_star * np.log(p_h1 / p_h0) + (m - m_star) * np.log((1 - p_h1) / (1 - p_h0))
n_step = np.where((B > z) | (z > A))[0][0] + 1
plt.plot(m[:n_step], z[:n_step], c='C0', alpha=0.3)
plt.plot(m[:n_step][-1], z[:n_step][-1], 'C0o', alpha=0.5)
plt.axhline(A, color='darkred', label=f'A = {A:.2f}')
plt.axhline(B, color='darkgreen', label=f'B = {B:.2f}')
plt.axvline(N, color='darkblue', label=r'$n^{\mathrm{(NP)}}$ = 100')
plt.legend(loc='upper right', framealpha=1)
plt.xlabel('Количество бросков (n)')
plt.ylabel(r'$z_{i}$')
plt.title(r'Последовательный критерий (верна $H_{0}$, $\alpha=0.044,\;\beta=0.087$)')
plt.show()Однако, нельзя не заметить, что иногда требуемое количество бросков может значительно превышать количество бросков требуемое критерием Неймана-Пирсона, которое отмечено вертикальной линией
Код для отрисовки
def wald_test_m(p_h0, p_h1, A, B, H):
m, m_star, z = 0, 0, 0
h0_count = 0
h1_count = 0
if H == 0:
p = p_h0
else:
p = p_h1
while B < z < A:
x = bernoulli.rvs(p)
m += 1
m_star += x
z = m_star * np.log(p_h1 / p_h0) + (m - m_star) * np.log((1 - p_h1) / (1 - p_h0))
if z < B:
h0_count += 1
if z > A:
h1_count += 1
return m, h0_count, h1_count
f, ax = plt.subplots(1, 2, figsize=(12, 4))
N_samples = 10000
n_steps_H0, h0_count_H0, h1_count_H0 = np.array([wald_test_m(p_h0, p_h1, A, B, H=0) for _ in range(N_samples)]).T
n_steps_H1, h0_count_H1, h1_count_H1 = np.array([wald_test_m(p_h0, p_h1, A, B, H=1) for _ in range(N_samples)]).T
sns.histplot(n_steps_H0, ax=ax[0])
ax[0].axvline(100, color='red', lw=3, label=r'$n^{(\mathrm{NP})}$ = 100')
steps_mean = np.mean(n_steps_H0)
ax[0].axvline(steps_mean, color='green', lw=3, label=r'$E(m)$ = ' + f'{steps_mean:.2f}')
alpha = h1_count_H0.sum() / N_samples
ax[0].text(150, 90, r'$\alpha = \frac{{{}}}{{{}}} = {}$'.format(h1_count_H0.sum(), N_samples, alpha))
ax[0].set_xlabel('Количество бросков в экспериментах (m)')
ax[0].legend()
ax[0].set_title(r'Верна $H_{0}$')
sns.histplot(n_steps_H1, ax=ax[1])
ax[1].axvline(100, color='red', lw=3, label=r'$n^{(\mathrm{NP})}$ = 100')
steps_mean = np.mean(n_steps_H1)
ax[1].axvline(steps_mean, color='green', lw=3, label=r'$E(m)$ = ' + f'{steps_mean:.2f}')
beta = h0_count_H1.sum() / N_samples
ax[1].text(150, 90, r'$\beta = \frac{{{}}}{{{}}} = {}$'.format(h0_count_H1.sum(), N_samples, beta))
ax[0].set_xlabel('Количество бросков в экспериментах (m)')
ax[1].legend()
ax[1].set_title(r'Верна $H_{1}$')
plt.show()Удивительно, но наш последовательный критерий работает в 2 раза быстрее (в среднем), чем критерий Неймана-Пирсона. Просто отвлекитесь и подумайте - что все это значит?
Давайте я вам немного помогу:
клинические испытания новых лекарств и вакцин;
A/B тестирование в крупных и новых сервисах;
контроль качества на производстве;
исследования в психологии, экономике, социологии;
мониторинг кибератак и мошенничества;
...;
мониторинг загрязнения окружающей среды.
Так или иначе, прямо или косвенно скорость и качество тестирования статистических гипотез оказывает влияние на качество вашей жизни и жизни ваших близких, детей, друзей, знакомых. Всех с кем вы связаны шестью рукопожатиями.
Но если это так важно, то почему последовательный критерий нигде не используется (на самом деле используется)? В самом деле, вам часто встречались статьи типа "Смотрите, мы заюзали последовательный критерий и порвали конкурентов!" (на самом деле они есть, но не с таким названием). А вам часто попадалось упоминание о последовательном критерии в учебниках? Лично, я смог найти всего один учебник с таким упоминанием (но это только на русском языке). Может сам последовательный тест с каким-то скрытым подвохом (подвохи есть и серьезные, но они перевешивают выгоду). Тоже нет, как только я смог в нем разобраться я сразу добавил его в модель ценообразования и он отлично улавливает аномальные отклонения спроса (не так хорошо, как на самом деле хотелось бы).
Так в чем же дело?
Это и есть тот самый скелет, мысли о котором не давали мне покоя все это время.
Критерий Вальда - рассвет, закат и снова рассвет
Наверняка вам доводилось слышать множество историй успешного успеха: "Я купил курс и получил высокооплачиваемую работу мечты в IT". Замечали, что иногда в комментариях под подобными историями всегда кто-то здравомыслящий упоминал об ошибке выжившего? Так вот, впервые, как известно, данную ошибку обнаружил венгерский математик Абрахам Вальд во время второй мировой войны. Но величайшая его заслуга вовсе не в этом. Во время второй мировой войны именно Абрахам Вальд в значительной степени развил и сформировал идеи последовательного статистического анализа. Естественно его теорию засекретили еще до того как он получил какие-то результаты - военные сами к нему обратились с просьбой найти метод ускорения анализа статистических гипотез. Почему засекретили?.. твердость материалов, прочность соединений, точность измерений, характеристики веществ: представьте что у вас есть возможность отвечать на вопросы быстрее противника или возможность ошибаться дешевле чем противнику.
Идеи последовательного статистического анализа были так же найдены Аланом Тьюрингом и использовались им для взлома Энигмы. Когда я узнал об этом, то по новому, почти физически почувствовал титаническую мощь этого человека. Ну а вы? Мурашки не побежали по коже?
Так может последовательный анализ не получил широкого распространения из-за секретности? Вряд ли. Теория развитая Вальдом была рассекречена практически сразу же после окончания войны. Хотя, результаты полученные Тьюрингом, были рассекречены значительно позже в 80-х годах. Может быть, если бы их рассекретили раньше, то... вряд ли что-то бы изменилось.
Причины, на мой взгляд, гораздо прозаичнее.
Фреквентизм победил не потому что он лучше, он просто оказался удобным для набиравшей обороты научной бюрократии. А мы все знаем что такое бюрократия - она всегда стремится занять все доступное пространство.
К тому же в дело снова включился Рональд фишер, но на этот раз он начал холиварить уже после смерти Вальда, назвав его неопытным ученым, который написал некомпетентную книгу по статистике. Кстати вот эта книга:
И кстати, замечательная книга. На мой взгляд, это до сих пор единственная лучшая книга по последовательному анализу. Вам никогда не казалось, что иногда в книги по статистике добавляют слишком много формул. Не знаю, может для солидности? Но зачем? Такая замечательная и крайне прикладная наука порой излагается на каком-то инопланетном языке. Кстати, рецензентом этой книги был страховщик, а эти люди понимают, что такое риски и как считать деньги. Мда, когда с тебя потихоньку начинает сыпаться песок начинаешь читать введения, благодарности - все, "от корки до корки" как говорится. Но мы что-то снова отвлеклись.
Рональд Фишер... Рональд Фишер был не просто авторитетом, он был самым настоящим идеологом. Осмелились бы вы идти против такого авторитета, который мог позволить себе публично называть ученых еретиками? Ты берешь и делаешь последовательный вероятностный вывод в процессе сбора данных, но почему то вдруг это называют кощунством, поскольку это нарушает "чистоту" эксперимента. Так или иначе фреквентизм захватил ключевые кафедры, журналы и научные советы. Инакомыслие, просто не пропускалось в публичное поле. Молодой ученый рискнувший использовать последовательный анализ в 50-60х годах, мог запросто лишиться карьеры, получив клеймо... нет не "еретика", а как сказала бы научная бюрократия - "ненадежного методолога".
Насколько я понимаю, есть что-то вроде стандарта в научных публикациях, который полностью "фиксированный": p-value, доверительные интервалы, ANOVA и т.д.
Есть ли другие причины забвения? Да. Подавляющее большинство грантов до сих пор требуют предопределенного плана. Вы говорите: "Мы будем тестировать гипотезу последовательно, пока не станет ясно, что...", а бюрократ слышит: "Шарлатанство, шарлатанство, шарлатанство". Фреквентистский дизайн с фиксированным
Фиксированный дизайн приятен, он создает ощущение чистоты и объективности: протокол написан до начала эксперимента, данные собираются "вслепую", анализ проводится один раз.
Давайте добавим сюда сложность расчетов. Помните как легко мы определили и вычислили границы
Серьезный баг системы - вот он, тот самый скелет в шкафу
Я очень хорошо помню, как пытался вникнуть в суть тестирования гипотез. Пазл никак не хотел складываться. И в какой-то из вечеров я решил просто посерфить интернет. И мне действительно попался хороший сайт дизайнера веб приложений. Интересно что именно дизайнер внес какое-то прояснение в моей голове (дорогой мой человек, не знаю как тебя зовут, но спасибо). Насколько я смог тогда разобраться, человек дизайнил приложение одного банка, но он не просто дизайнил, а еще и тестил разные дизайны. Это был первый реально прикладной пример важности тестирования гипотез в бизнесе, который прибавил мне энтузиазма.
Вы когда-нибудь сталкивались с откровенно идиотскими дизайнерскими решениями? Но вы скорее всего не задумывались о том, сколько компания, выкатившая такое приложение, может сэкономить используя даже самое банальное A\B-тестирование? Это можно посчитать. Некоторые крупные компании, особенно IT-гиганты уже давно все посчитали, в том числе и выгоды от последовательного анализа.
Когда крупная технологическая компания запускает фичу
В фармакологии осознание того, что последовательный анализ может сэкономить сотни миллионов долларов и годы времени, перевесило бюрократическое сопротивление. Теперь регуляторы требуют, что бы в крупных испытаниях использования обязательно использовались промежуточные результаты и аналитика.
Упоминания об использовании последовательного анализа теперь можно найти то тут, то там. Конечно, современные вводные учебники по статистике до сих пор являются памятниками "победы" фреквентизма над байесовизмом, но сейчас определенную ценность начинает набирать именно кругозор специалистов по статистике. Бизнес становится более "горизонтален" - можно спокойно постучаться к руководителю и сказать, что "вот так лучше и быстрее".
Казалось бы, что никакой проблемы теперь нет, а скелет в шкафу - это просто артефакт эпохи зарегулированной всем чем угодно. Теперь можно "спокойно дышать". Но все-таки есть кое что подозрительное и во всем этом тоже. Давайте задумаемся о том сколько данных собирается последовательно и о том какая доля этих данных подвергается последовательному анализу. Можно ли сделать абсолютно все тестирование гипотез последовательным?
Последовательный анализ можно сравнить с гоночным болидом Формулы-1, который невероятно быстр на идеальных трассах (простые гипотезы, независимость данных). Изучая последовательный анализ довольно быстро начинаешь понимать границы его применимости. Скажем так, болид Формулы-1 плохо подходит для поездки в магазин за покупками или бездорожья (составные гипотезы, многомерные данные). За последние 70 лет разработано множество обходных путей и приближенных методов, которые расширяют применимость идей последовательного анализа на более широкий класс задач, но всегда ценой потери той самой математической чистоты и оптимальности, которые были в оригинальной работе Вальда.
Например, если последовательный эксперимент вдруг прервется, то придется придумывать какие-то механизмы поправок. Допустим, если с нашей монеткой из вышеприведенного примера эксперимент пришлось оборвать после 40-го броска и у нас выпало 25 орлов и 15 решек, то необходимо задуматься о каком-нибудь механизме хоть какой-то интерпретации получившихся результатов и выработки необходимых поправок.
Что мы можем сделать? Например, мы можем взглянуть на множество перестановок 25 орлов и 15 решек и то какими могут быть максимальные и минимальные значения отношения правдоподобий:
Код для отрисовки
m = np.arange(1, 41)
S = []
for _ in range(10000):
samples = np.random.choice([0]*15+[1]*25, size=40, replace=False)
m_star = np.cumsum(samples)
z = m_star * np.log(p_h1 / p_h0) + (m - m_star) * np.log((1 - p_h1) / (1 - p_h0))
S.append(z)
S = np.array(S)
plt.figure(figsize=(12, 3))
sns.histplot(np.sort(S, axis=1)[:, 0], binwidth=0.25, label=r'$\max(\Lambda)$')
sns.histplot(np.sort(S, axis=1)[:, -1], binwidth=0.25, label=r'$\min(\Lambda)$')
plt.axvline(A, color='red', lw=3, label='A')
plt.axvline(B, color='green', lw=3, label='B')
plt.axvline((A + B) / 2, color='k', lw=1, ls=':')
plt.legend();Есть ли в этом рисунке что-нибудь, что может склонить наше предпочтение в сторону
В последовательном анализе мы можем практически "воочию" лицезреть порядок и беспорядок. Например, если при тех же 40-ка бросках монеты и тех же гипотезах вдруг выпадет 35 решек и всего 5 орлов и если мы взглянем на отсортированные значения отношений правдоподобия в каждом эксперименте или на все возможные точки траекторий и сколько раз в мы в каждой из них побывали, то получим вот такую картину:
Код для отрисовки
m = np.arange(1, 41)
S = []
for _ in range(10000):
samples = np.random.choice([0]*35+[1]*5, size=40, replace=False)
m_star = np.cumsum(samples)
z = m_star * np.log(p_h1 / p_h0) + (m - m_star) * np.log((1 - p_h1) / (1 - p_h0))
S.append(z)
S = np.array(S)
f, ax = plt.subplots(1, 2, figsize=(12, 4))
for i in range(200):
ax[0].plot(np.sort(S[i]), 'b-', alpha=0.1)
ax[0].axhline(A, color='darkred', label=f'A = {A:.2f}')
ax[0].axhline(B, color='darkgreen', label=f'B = {B:.2f}')
for i in range(40):
ax[1].plot(*np.unique(S[:, i], return_counts=True), 'bo', ms=3)На графике все слишком хорошо и красиво что бы быть правдой. Может с нашими гипотезами что-то не так раз мы наблюдаем такую красоту?
Мы использовали случайные перестановки, что бы увидеть все возможные траектории эксперимента, а что если использовать бутстрап? Наверное, можно даже не пытаться и заведомо сказать, что это глупость.
А что если объединить идеи последовательного анализа с порядковыми статистиками? Можно ли заведомо сказать, что это тоже глупость?
А что если объединить идеи последовательного анализа с Байесовской статистикой? Например, в нашем примере эксперимента с монеткой, который оборвался на 40-ом броске, мы могли бы последовательно строить апостериорное бета распределение и наблюдать за тем как меняется наша степень уверенности. Может мы смогли бы придумать какое-то последовательное правило и... помните как в начале статьи мы отметили, что
Есть мнение, что все предопределено, а наша неопределенность связана с незнанием. Например, мы можем резко снизить свою неопределенность относительно нашей монетки если обзаведемся высокоскоростной камерой. А если мы узнаем, вес монетки, ее размеры, температуру воздуха и т.д., то мы можем свести нашу неопределенность практически к нулю. Может быть и так - все есть эпистемическая случайность, которая исчезает с ростом информации. Демон Лапласа - это зеркало в котором отражается наше несовершенство.
Но может быть, например, квантовые флуктуации вакуума - это и есть та самая "чистая", онтологическая случайность, которая не исчезнет даже если бы мы знали все? Может быть это не просто "шум", а фундаментальный механизм, определяющий саму ткань реальности?
Статистика - это язык для работы с обеими случайностями. Мы используем вероятность как инструмент, как для моделирования нашего незнания о скрытых параметрах (байесовский подход), так и для описания поведения принципиально-стохастических систем. Это не недостаток, а универсальная гениальность вероятностного языка.
Что происходит с квантовыми флуктуациями на границах черных дыр? Сама по себе черная дыра - это мертвое равновесие, символ абсолютного детерминизма и порядка. Но природа берет случайность и порождает на ее границе пару частиц из энергии вакуума, нарушая монолитность черной дыры и заставляя ее испаряться.
Если бы случайности не существовало, то все что снаружи черных дыр можно было бы воспринимать как единый, невероятно сложный, но детерминированный часовой механизм, отсчитывающий время до неминуемой тепловой смерти вселенной.
Случайность - это антиэнтропийный механизм, который вбрасывает микровариации и флуктуации в любые системы. Практически все они исчезают, но некоторые усиливаются и дают начало новым структурам: новым химическим элементам в звездах, новым видам в эволюции, новым идеям в сознании.
Можно ли сказать, что природа защищает свои тайны с помощью случайности и неопределенности и наш смысл жизни состоит в том, чтобы преодолеть этот механизм защиты? Я вижу, что иногда к статистике относятся, не просто как к некой взломщице, которая пытается отключить "сигнализацию" Вселенной. А как к оружию способному уничтожить случайность. Но это не так. Смысл не в том, что бы найти "Истину" с большой буквы, потому что вряд ли это вообще возможно. Смысл жизни мыслящего вида - не в том что бы победить случайность, а втом что бы научиться извлекать из нее пользу, точно так же как это возможно с другими природными стихиями.
Настоящий скелет в шкафу - это глубокая связь между статистикой, теорией игр, принятием решений и оптимизацией, которую не то что бы не принято, а вообще сложно замечать. Абрахам Вальд в своей книге показал, что тестирование гипотез можно воспринимать как игру с природой. Природа загадывает некоторое истинное распределение случайной величины
Насколько большим может быть риск? Преодолели бы мы COVID-19 с меньшими потерями если бы мы не знали об адаптивных схемах дизайна экспериментов. А что насчет, возможной полной резистивности бактерий к антибиотикам? Может быть прямо сейчас к нашей планете движется огромный астероид? Может курение и малоподвижный образ жизни действительно убивают?
Даже когда у нас нет никаких данных, мы можем проранжировать наши гипотезы с помощью априорных шансов. Допустим у нас потерялась зарядка от телефона, где ее искать? Возле кровати за тумбочкой? В доме друзей где мы недавно были в гостях? В холодильнике? Может на Марсе, ведь ее запросто могли украсть марсеане для изучения человеческих технологий? А ведь ее могли украсть и кошки - они кажется давно что-то замышляют!
Любой человек, который принимает решения, оперирует гипотезами. Не важно о чем идет речь: о выборе новой стратегии для крупной корпорации или о том, что съесть сегодня на завтрак. Любой человек, принимающий решения, осознанно или неосознанно говорит на языке вероятностей. Последовательный анализ - это возможность вести с природой диалог. И от качества этого диалога будет зависеть очень многое будь то судьбы корпораций, определяемые выбранными стратегиями или наше самочуствие в течение дня, определяемое тем что мы съели на завтрак.