Лет 5 назад я усиленно пытался вникнуть в тервер и статы: книги, статьи, вебсёрфинг. Даже написал несколько статей: раз, два, три. Вообще, в планах было написать довольно большой цикл статей, что бы подсветить какие-то самые сложные вещи, да и самому в них разобраться - совместить полезное с полезным, так сказать. Однако, в какой-то момент я решил, что полученных знаний достаточно для новых проектов и ушел в работу. Работал. Работал. Работал.

Но знаете, все это время меня никак не покидало ощущение, что в процессе своего знакомства с наукой о случайностях я наткнулся на какой-то скелет в шкафу. И знаете что? Мне кажется настало время вытащить этот скелет на всеобщее обозрение.

Есть в тестировании гипотез что-то подозрительное

Почему проверка гипотез, в том числе и статистическая, так важна? Все просто - в конечном счете, качество нашей жизни прямо пропорционально качеству и количеству наших гипотез которые мы можем проверять. Четкие, обоснованные и проверяемые допущения ведут к эффективным решениям, а размытые и ложные к череде ошибок и разочарований. Не формулировать никаких гипотез - это тоже гипотеза, которая состоит в молчаливом и безрассудном допущении, что всё останется по-старому и любые наши действия не будут иметь никаких последствий. Не формулировать никаких гипотез - это самая рискованная ставка из всех возможных.

Гиганты, на плечах которых мы стоим понимали это. Основная проблема связанная с проверкой статистических гипотез, заключается в установке некоторых правил или точнее критериев в соответствии с которыми выдвинутая гипотеза принимается или отвергается. Такие критерии впервые были сформулированы Джоном фон Нейманом и Эгоном Пирсоном. Давайте кратко рассмотрим основные аспекты теории, которая стоит за этим критерием.

Допустим у нас есть монетка в честности которой мы хотим убедиться. Мы подкинули монетку 100 раз и у нас выпало 57 орлов и 43 решки. Основная гипотеза - H_{0} состоит в том, что монетка честная, т.е. вероятность выпадения орла p_{0}=0.5, а конкурирующая гипотеза - H_{1} состоит в том, что монетка нечестная и вероятность выпадения орла равна p_{0}=0.65.

У нас есть 100 наблюдений и очевидно, что если мы подбросим ту же самую монетку еще 100 раз, то результат может оказаться другим. Но разные результаты окажутся более или менее вероятными. Что бы принять гипотезу H_{0} и отвергнуть H_{1} нам необходимо простое правило, какие результаты наблюдений считать допустимыми для H_{0}, а какие критическими. Поскольку результаты наблюдений представляют собой количество "успехов" - выпадений орла в 100 независимых экспериментах (подбрасывании монетки), то распределение результатов подчиняется биномиальному закону. Построим две функции вероятности, которые соответствуют распределению результатов для наших гипотез:

  • f_{0}(x) = P(X = x \; | \; H_{0}) = B(x, n=100, p_{0}=0.5);

  • f_{1}(x) = P(X = x \; | \; H_{1}) = B(x, n=100, p_{1}=0.65).

В результате мы увидим следующее:

Код для отрисовки
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()

На графике видно, что функции перекрываются и очевидно, что ошибки связанные с принятием той или иной гипотезы находятся именно в этой области перекрытия. Итак, что бы принять H_{0} нам необходимо два множества: множество допустимых результатов наблюдений для H_{0} и множество недопустимых результатов для H_{0}. У нас выпало 57 орлов, тогда пусть все результаты, в которых количество орлов будет больше этого значения окажутся в критической области X \geqslant 58, а все остальные результаты в допустимой области. Так мы (довольно произвольно) определили некоторое значение k=58, которое является пограничным для двух множеств. Попадание результатов наблюдений в критическую область означает отклонение H_{0} и принятие H_{1}. Однако, такое отклонение может быть ошибочным, потому что даже при верности H_{0}, т.е. даже при условии, что монетка честная, вероятность попадания в критическую область будет составлять около 0.07. Традиционно, такую ошибку принято обозначать греческой буквой \alpha и называть ошибкой первого рода или уровнем критической области, но чаще всего - уровнем значимости.

В тоже самое время, когда мы определили критическую область мы так же создали возможность для появления еще одной ошибки. Если, на самом деле монетка нечестная, т.е. на самом деле верна H_{1}, возникает некоторая вероятность того, что результаты наблюдений за сотней бросков могут оказаться меньше выбранной границы k. Поскольку эти значения находятся в допустимой области для H_{0}, то мы можем принять H_{0}, когда на самом деле верна гипотеза H_{1} и конкретно в нашем случае вероятность такой ошибки составляет около 0.06. Традиционно, такую ошибку принято обозначать греческой буквой \beta и называть ошибкой второго рода. Величину 1 - \beta называют мощностью критической области (или просто мощностью).

Из вышеприведенного графика так же следует, что если для некоторой критической области мы получаем наименьшие значения \alpha и \beta, то такая критическая область является для нас наиболее предпочтительной. Это утверждение может показаться самоочевидным, но представьте, что вы незнакомы с теорией вероятностей, как бы вы принимали решение?

Выбор критической области позволяет сделать сколь угодно малыми либо \alpha либо \beta. Например, если мы уменьшим \alpha, то у нас увеличится \beta. В чем легко убедиться, взглянув на нижеследующий график.

Код для отрисовки
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()

При фиксированном объеме выборки мы не можем одновременно уменьшить и \alpha и \beta. Единственный способ снизить общую долю ложных выводов заключается в увеличении объема выборки.

Представленный пример с монеткой очень прост. В общем случае гипотезы могут быть крайне сложными, а функции распределения вероятностей не всегда возможно задать аналитически. Глядя на графики, можно подумать, что нам достаточно:

  • выбрать приемлемое для нас значение \alpha;

  • дальше, исходя из выражения P(X \geqslant k) \leqslant \alpha, найти значение k - границу критической области;

  • ну а значение \beta появится как-бы само собой.

В общем случае это не так. Очень часто для заданного уровня \alpha существует очень много критических областей и выбор какой-то конкретной критической области должен основываться на следующем принципе, выдвинутом Нейманом и Пирсоном: ограничивая наш выбор областями, для которых вероятность \alpha равна заданной величине, мы должны выбрать ту область из этого множества, для которой вероятность \beta минимальна.

Если следовать этому принципу, то при фиксированном объеме выборки вероятность \beta будет однозначной функцией от \alpha. Следовательно, задавая количество наблюдений, на которых основана проверка гипотезы, мы можем произвольно выбрать одну из величин \alpha или \beta. И наоборот, задавая произвольно одну из величин \alpha или \beta, мы можем найти необходимое количество наблюдений.

Выбор объема выборки и уровня \alpha (или \beta) хоть и могут быть произвольными, но на самом деле очень сильно зависят от контекста. Дело в том, что в разных ситуациях ошибки первого и второго рода могут иметь разную важность, например, если ошибка первого рода обходится нам в один миллион рублей, а ошибка второго рода в сто тысяч рублей, тогда мы постараемся минимизировать \alpha, даже если это приведет к увеличению \beta. Однако, все наши старания будут ограничиваться объемом выборки, поскольку каждое отдельное наблюдение так же имеет некоторую стоимость.

Кажется, что на этом все. Однако, теория Неймана-Пирсона несколько глубже чем можно подумать.

В определенной степени принято считать, что в математической статистике есть два лагеря: лагерь фреквентистов - сторонников частотных методов в статистике и лагерь байесовистов - сторонников байесовских методов вывода. Причем, вторые, как правило, довольно пренебрежительно отзываются о первых (иногда даже слишком). Например, Кэмерон Дэвидсон-Пайлон в своей книге "Вероятностное программирование на Python: байесовский вывод и алгоритмы." пишет о фреквентизме и A/B-тестеровании следующее:

Зачастую анализ после эксперимента производится с помощью так называемой проверки статистической гипотезы, например проверки различия средних значений (difference of means test) или проверки различия долей (difference of proportions test). Это требует использования такого запутанного понятия, как «Z-оценка», и еще более запутанных малопонятных «p-значений» (даже не спрашивайте, что это такое). Если вы слушали курс статистики, то вас, возможно, обучали этой методике (хотя вовсе не факт, что вы в ней разобрались). И если мы с вами в этом схожи, то их вывод отнюдь не доставляет вам радости. Если так — отлично. Байесовский подход к решению этой задачи гораздо более естественен. (Дэвидсон-Пайлон Кэмерон. Вероятностное программирование на Python: байесовский вывод и алгоритмы. — СПб.: Питер, 2019. с. 62).

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

  • Байесовские методы - это просто формализация предубеждений. Вы получаете на выходе то, во что изначально верили. Это не наука, а самооправдание. Настоящая наука должна основываться на объективных данных, а не на чьих-то личных верованиях.

  • Байесовство - это не статистика, это вера. Вы, байесовцы, просто верите в свой априор, а потом находите данные, которые подтверждают вашу веру. Это слепая религия, а не наука.

  • Байесовские методы — это интересная академическая игрушка для простых задач с игральными костями, урнами с шарами и монетками. Но для реальной науки они бесполезны и непрактичны.

  • Если вы не получили желаемый результат, вы всегда можете поиграть с априорным распределением, пока не получите его. Это легализованный p-hacking на стероидах.

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

Зачем такое длинное отступление?

изображение.png
изображение.png

Противоборство фреквентистов и байесовистов, это яркий пример того, что даже в научном сообществе могут царить разногласия, инерционность и консерватизм. Можно ли утверждать, что сейчас такого нет? Помните в самом начале я говорил о скелете в шкафу? Так вот, на самом деле я даже не уверен в том - скелет ли это (потому что долгие годы вся окружающая действительность убеждала меня в том, что это не так). Но давайте не будем забегать наперед. Мы еще не закончили с теорией Неймана-Пирсона.

Давайте снова вернемся к нашей монетке, но на этот раз взглянем на результат с точки зрения байесовизма. Запишем теорему Байеса для наших гипотез в пропорциональной форме:

P(H_{0} \; | \; X=x) \propto P(H_{0}) \times P(X=x \; | \; H_{0}),P(H_{1} \; | \; X=x) \propto P(H_{1}) \times P(X=x \; | \; H_{1}).

Пропорциональная форма теоремы Байеса удобна тем, что мы можем ее использовать для сравнения двух гипотез. Для этого достаточно просто разделить один апостериор на другой:

\frac{P(H_{0} \; | \; X=x)}{P(H_{1} \; | \; X=x)} = \frac{P(H_{0})}{P(H_{1})} \times \frac{P(X=x \; | \; H_{0})}{P(X=x \; | \; H_{1})}.

Если это отношение равно 3, то можно сказать, что гипотеза H_{0} объясняет наблюдаемое количество орлов в 3 раза лучше чем H_{1} и наоборот, если отношение равно 1/3, то H_{1} объясняет наблюдения в три раза лучше чем H_{0}. Справедливости ради, надо сказать, что до сих пор есть люди которые считают именно такое соотношение достаточным критерием для принятия и отклонения гипотез, хотя конечно же это неверно.

Отношение апостериорных вероятностей состоит из двух частей:

O(H_{0} : H_{1}) = \frac{P(H_{0})}{P(H_{1})}

и

\Lambda(H_{0} : H_{1} \; | \; X=x) = \frac{P(X=x \; | \; H_{0})}{P(X=x \; | \; H_{1})}.

O(H_{0} : H_{1}) называется априорными шансами и выражает нашу степень уверенности в той или иной гипотезе, еще до рассмотрения каких бы то ни было данных. Поскольку, до появления данных о подбрасывании нашей монетки, мы не можем делать никаких предпочтений относительно H_{0} или H_{1}, то получается что P(H_{0}) = P(H_{1}). В этом случае остается только \Lambda(H_{0} : H_{1} \; | \; X=x) - эта часть отношения, иногда называемая коэффициентом Байеса, является отношением правдоподобий и выражает отношение вероятностей получения данных если бы обе гипотезы были верными.

Очевидно, что если гипотеза H_{0} объясняет данные лучше чем H_{1}, то она будет "толкать" \Lambda вверх, а если хуже, то вниз. В этом легко убедиться, если смоделировать несколько выборок подбрасываний честной монетки:

Код для отрисовки
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()

Обратите внимание, что чем лучше H_{0} объясняет данные тем сильнее она толкает отношение правдоподобия вверх. Но немного подозрительно, что этот рост не является монотонным.

Если мы взглянем на график изменения отношения правдоподобия при верной гипотезе H_{1}, то увидим схожую картину, только теперь все будет убывать:

Код для отрисовки
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()

На самом деле Нейман и Пирсон показали, что область состоящая из всех возможных выборок X для которых выполняется неравенство

\Lambda(H_{0} : H_{1} \; | \; X) \geqslant k

является наиболее мощной критической областью для проверки гипотезы H_{0} относительно гипотезы H_{1}. При этом величина k подбирается точно так же как и раньше, из условия необходимого уровня \alpha. Сейчас этот факт известен как лемма Неймана-Пирсона.

Интересно, что писал Сам Джон фон Нейман об этом открытии:

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()

Из-за логарифмирования \Lambda данный график нельзя считать абсолютно достоверным, однако, он хорошо отражает суть леммы Неймана-Пирсона. И вот что странно, эта лемма концентрирует наше внимание только на том где мы оказались, но не на том как мы туда попали. В самом деле, если на первом графике точка оказалась выше 4 то мы считаем монетку честной, а если на втором графике точка оказалась меньше -4, то мы считаем ее нечестной. Но давайте взглянем на траектории этих точек, некоторые из них пересекли указанные границы гораздо раньше. Должны ли мы все равно продолжать эксперимент до 100 бросков, если граница была пересечена, скажем, после 25-го или 80-го броска?

В лемме Неймана-Пирсона ничего не сказано о размере выборки, но формально она выглядит так:

\frac{f_{0}(x_{1})f_{0}(x_{2})...f_{0}(x_{n})}{f_{1}(x_{1})f_{1}(x_{2})...f_{1}(x_{n})} \geqslant k,

где f_{0}(x) = f(x, \theta_{0}) - это распределение нашей случайной величины при \theta_{0} = p_{0}, а f_{1}(x) = f(x, \theta_{1}) - это распределение при \theta_{1} = p_{1} (конкретно в нашем случае - это просто распределение Бернулли). Ну а n - это и есть объем выборки. Так, что да, если n наблюдений обеспечивают заданные \alpha и \beta, то мы обязаны сделать n наблюдений. Давайте взглянем на то как будет меняться \beta теста с ростом числа бросков при \alpha = 0.05.

Код для отрисовки
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()

То есть количество бросков жестко задано уровнем \alpha и соответствующей этому уровню минимальной вероятностью \beta.

Но согласитесь, после увиденного все это кажется немного нелогичным. Что бы это продемонстрировать можно привести простой пример. Допустим, мы тестируем гипотезу которая состоит в том, что монетка либо абсолютно честная H_{0}: p_{0} = 0.5, либо абсолютно нечестная H_{1}: p_{1} = 1. В таком случае, если на каком-то броске такой монетки с номером m выпадает решка, то в дальнейших наблюдениях не будет никакой необходимости. Согласен, пример спорный.

изображение.png

Но он показывает, что множество всех выборок с разным объемом m можно разделить не на две, а на три не пересекающиеся области:

  1. Критическая область для гипотезы H_{0}, например, все векторы длинны m \geqslant 10, состоящие из одних единиц.

  2. Критическая область для гипотезы H_{1}, например, любой вектор длинны m, который начинается с m-1 единиц и заканчивается одним нулем.

  3. Допустимая область для обеих гипотез, например, все векторы длинны m \leqslant 9, состоящие из одних единиц.

По крайней мере, хотя бы в этом примере мы не обязаны выполнять фиксированное и заданное наперед количество бросков, а прервать эксперимент на любом шаге m \leqslant 9. Но возможно ли такое для каких-то более адекватных примеров? ДА!

Выше мы видели, что отношение правдоподобий \Lambda (коэффициент Байеса) по мере появления новых наблюдений может изменяться с разной скоростью. Мы знаем, что для любого количества наблюдений m, вероятность получения выборки (x_{1}, x_{2},..., x_{m}) задается следующим образом

p_{0m} = f(x_{1}, \theta_{0})...f(x_{m}, \theta_{0}),

когда справедлива гипотеза H_{0}, и выражением

p_{1m} = f(x_{1}, \theta_{1})...f(x_{m}, \theta_{1}).

Так же, выше мы видели, что если \Lambda оказывалась ниже заданной границы, которая равнялась \ln(1/58), то мы отвергали гипотезу H_{0}, а если выше границы равной \ln(58), то мы ее принимали (напомню, график отображает лишь принцип, а не истинное положение вещей). Раз так, то очевидно что для любых p_{0m} и p_{1m}, мы можем определить три интересующих нас области с помощью следующего выражения

B \leqslant \frac{p_{1m}}{p_{0m}} \leqslant A,

где A и B некоторые значения, определяющие границы критических и допустимых областей для обеих гипотез.

Для удобства прологарифмируем

\ln \frac{p_{1m}}{p_{0m}} = \ln \frac{f(x_{1}, \theta_{1})}{f(x_{1}, \theta_{0})} + ... + \ln \frac{f(x_{m}, \theta_{1})}{f(x_{m}, \theta_{0})}

и обозначим i-й член последовательности через z_{i}

z_{i} = \ln \frac{f(x_{i}, \theta_{1})}{f(x_{i}, \theta_{0})}.

Для нашего примера с монеткой, мы можем записать следующее:

z_{i} = \ln \frac{p_{1m}}{p_{0m}},

если выпал орел, и

z_{i} = \ln \frac{1- p_{1m}}{1 - p_{0m}},

если выпала решка. Тогда правило по которому определяются три интересующих нас области можно записать вот так

\ln(B) \leqslant m^{*} \ln \frac{p_{1m}}{p_{0m}} + (m - m^{*}) \ln \frac{1- p_{1m}}{1 - p_{0m}} \leqslant ln(A),

где m^{*} - это число орлов в первых m наблюдениях.

Единственная серьезная проблема для нашего правила состоит в том, что мы не знаем ка определять необходимые значения A и B. Но их можно определить, исходя из следующих соображений. У нас есть критическая область для H_{0}, которая отвергается если

\frac{p_{1m}}{p_{0m}} \geqslant A

и критическая область для H_{1}, которая отвергается если

\frac{p_{1m}}{p_{0m}} \leqslant B.

Тогда для любой выборки (x_{1}, x_{2},..., x_{n}) у которой для всех промежуточных значений m выполняется неравенство

B \leqslant \frac{p_{1m}}{p_{0m}} \leqslant A,

но конкретно для последнего шага n выполняется условие

\frac{p_{1m}}{p_{0m}} \geqslant A,

то вероятность получения любой такой выборки будет по крайней мере в A раз больше при верности гипотезы H_{1}, чем при верности H_{0} (это утверждение не кажется таким очевидным, как хотелось бы, но в нем можно убедиться численно, используя лемму Неймана-Пирсона для разных n). Вероятность получить любую такую выборку так же означает, что процесс тестирования гипотез окончится принятием H_{1} и отклонением H_{0}. Как мы помним эта вероятность равна \alpha если H_{0} верна и 1 - \beta когда верна H_{1}. Отсюда вытекает простое неравенство

A \cdot \alpha \leqslant 1 - \beta,

которое можно переписать как

A \leqslant \frac{1 - \beta}{\alpha}.

С помощью аналогичных рассуждений мы можем прийти к тому, что

B \geqslant \frac{\beta}{1 - \alpha}.

Выходит, что теперь у нас есть критерий, который больше не требует фиксированных объемов выборок? Давайте проверим на нашей монетке.

Для тестирования по правилу (далее будем называть критерием) Неймана-Пирсона, нам необходимо задать количество бросков и необходимый уровень \alpha. Исходя из заданных параметров, мы вычисляем соответствующую вероятность \beta и границу критической области k. Пусть количество бросков будет равно, как и раньше 100, а уровень 0.05. Но поскольку речь идет о дискретном распределении, то задать уровень который будет в точности равен 0.05 не получится. Выберем в качестве \alpha ближайшее к 0.05 снизу значение. Тогда получим:

Код для отрисовки
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, то мы считаем, что монетка нечестная. Теперь посмотрим получится ли у нас прийти к таким же выводам при произвольном количестве подбрасываний, с таким же уровнем и мощностью.

Первым делом вычислим значения для границ A и B:

\ln(A) = \ln \frac{1 - \beta}{\alpha} = \ln \frac{1 - 0.087}{0.044} = \ln(22.53) = 3.03 \; ,\ln(B) = \ln \frac{\beta}{1 - \alpha} = \ln \frac{0.087}{1 - 0.044} = \ln(0.091) = -2.4 \; .

Ну а дальше мы просто повторим тест 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-го броска монеты. Тоже самое справедливо и для верной гипотезы H_{1}.

Код для отрисовки
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()

Однако, нельзя не заметить, что иногда требуемое количество бросков может значительно превышать количество бросков требуемое критерием Неймана-Пирсона, которое отмечено вертикальной линией n^{\mathrm{(NP)}}. Вообще было бы любопытно взглянуть на распределение возможного количества бросков монеты совершаемых при последовательном критерии:

Код для отрисовки
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 и т.д.

Есть ли другие причины забвения? Да. Подавляющее большинство грантов до сих пор требуют предопределенного плана. Вы говорите: "Мы будем тестировать гипотезу последовательно, пока не станет ясно, что...", а бюрократ слышит: "Шарлатанство, шарлатанство, шарлатанство". Фреквентистский дизайн с фиксированным N идеален в административной среде: четкий бюджет, четкие сроки, понятный план работ. Бюрократия это еще и конвейер.

Фиксированный дизайн приятен, он создает ощущение чистоты и объективности: протокол написан до начала эксперимента, данные собираются "вслепую", анализ проводится один раз.

Давайте добавим сюда сложность расчетов. Помните как легко мы определили и вычислили границы A и B для нашего эксперимента с монеткой. На самом деле вычисление границ - далеко нетривиальная задача. Аппроксимация предложенная Вальдом, которую мы собственно и использовали, является несколько грубой. В большинстве случаев, определение этих границ возможно только численно. А зачем тратить время на исследования для которых нет нужных вычислительных мощностей. Пожалуй это единственная объективная причина.

Серьезный баг системы - вот он, тот самый скелет в шкафу

Я очень хорошо помню, как пытался вникнуть в суть тестирования гипотез. Пазл никак не хотел складываться. И в какой-то из вечеров я решил просто посерфить интернет. И мне действительно попался хороший сайт дизайнера веб приложений. Интересно что именно дизайнер внес какое-то прояснение в моей голове (дорогой мой человек, не знаю как тебя зовут, но спасибо). Насколько я смог тогда разобраться, человек дизайнил приложение одного банка, но он не просто дизайнил, а еще и тестил разные дизайны. Это был первый реально прикладной пример важности тестирования гипотез в бизнесе, который прибавил мне энтузиазма.

Вы когда-нибудь сталкивались с откровенно идиотскими дизайнерскими решениями? Но вы скорее всего не задумывались о том, сколько компания, выкатившая такое приложение, может сэкономить используя даже самое банальное A\B-тестирование? Это можно посчитать. Некоторые крупные компании, особенно IT-гиганты уже давно все посчитали, в том числе и выгоды от последовательного анализа.

Когда крупная технологическая компания запускает фичу A, то теперь с гораздо большей вероятностью она используют последовательный дизайн и в реальном времени следит с какой вероятностью эта фича лучше старой фичи B и насколько.

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

Упоминания об использовании последовательного анализа теперь можно найти то тут, то там. Конечно, современные вводные учебники по статистике до сих пор являются памятниками "победы" фреквентизма над байесовизмом, но сейчас определенную ценность начинает набирать именно кругозор специалистов по статистике. Бизнес становится более "горизонтален" - можно спокойно постучаться к руководителю и сказать, что "вот так лучше и быстрее".

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

Последовательный анализ можно сравнить с гоночным болидом Формулы-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();

Есть ли в этом рисунке что-нибудь, что может склонить наше предпочтение в сторону H_{0} или H_{1}? Тут есть центр масс и это что-то вроде весов? Все эти и подробные вопросы не возникли бы, если бы последовательный анализ не имел ограничений. Да и вообще, использование множества перестановок не кажется чем-то обоснованным. Но с другой стороны, это способ более точной калибровки и возможность обнаружения "загрязнений" выборки прямо налету.

В последовательном анализе мы можем практически "воочию" лицезреть порядок и беспорядок. Например, если при тех же 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-ом броске, мы могли бы последовательно строить апостериорное бета распределение и наблюдать за тем как меняется наша степень уверенности. Может мы смогли бы придумать какое-то последовательное правило и... помните как в начале статьи мы отметили, что \alpha и \beta выбираются исходя из контекста, т.е. из стоимости ошибок и наблюдений, но что если контекст достаточно шире? Неужели желание сделать этот контекст частью эксперимента, это тоже глупая идея?

Есть мнение, что все предопределено, а наша неопределенность связана с незнанием. Например, мы можем резко снизить свою неопределенность относительно нашей монетки если обзаведемся высокоскоростной камерой. А если мы узнаем, вес монетки, ее размеры, температуру воздуха и т.д., то мы можем свести нашу неопределенность практически к нулю. Может быть и так - все есть эпистемическая случайность, которая исчезает с ростом информации. Демон Лапласа - это зеркало в котором отражается наше несовершенство.

Но может быть, например, квантовые флуктуации вакуума - это и есть та самая "чистая", онтологическая случайность, которая не исчезнет даже если бы мы знали все? Может быть это не просто "шум", а фундаментальный механизм, определяющий саму ткань реальности?

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

Что происходит с квантовыми флуктуациями на границах черных дыр? Сама по себе черная дыра - это мертвое равновесие, символ абсолютного детерминизма и порядка. Но природа берет случайность и порождает на ее границе пару частиц из энергии вакуума, нарушая монолитность черной дыры и заставляя ее испаряться.

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

Случайность - это антиэнтропийный механизм, который вбрасывает микровариации и флуктуации в любые системы. Практически все они исчезают, но некоторые усиливаются и дают начало новым структурам: новым химическим элементам в звездах, новым видам в эволюции, новым идеям в сознании.

Можно ли сказать, что природа защищает свои тайны с помощью случайности и неопределенности и наш смысл жизни состоит в том, чтобы преодолеть этот механизм защиты? Я вижу, что иногда к статистике относятся, не просто как к некой взломщице, которая пытается отключить "сигнализацию" Вселенной. А как к оружию способному уничтожить случайность. Но это не так. Смысл не в том, что бы найти "Истину" с большой буквы, потому что вряд ли это вообще возможно. Смысл жизни мыслящего вида - не в том что бы победить случайность, а втом что бы научиться извлекать из нее пользу, точно так же как это возможно с другими природными стихиями.

Настоящий скелет в шкафу - это глубокая связь между статистикой, теорией игр, принятием решений и оптимизацией, которую не то что бы не принято, а вообще сложно замечать. Абрахам Вальд в своей книге показал, что тестирование гипотез можно воспринимать как игру с природой. Природа загадывает некоторое истинное распределение случайной величины F, а статистик должен его определить. Он создает множество гипотез {F_{1}, F_{2}, ..., F_{i}}, но в этот же момент неминуемо начинает чем-то рисковать. Риск состоит в том, что неверные решающие правила могут сделать эксперимент слишком дорогим или слишком долгим, но решающие правила так же определяют стоимость последствий от выбора той или инной гипотезы F_{i}.

Насколько большим может быть риск? Преодолели бы мы COVID-19 с меньшими потерями если бы мы не знали об адаптивных схемах дизайна экспериментов. А что насчет, возможной полной резистивности бактерий к антибиотикам? Может быть прямо сейчас к нашей планете движется огромный астероид? Может курение и малоподвижный образ жизни действительно убивают?

Даже когда у нас нет никаких данных, мы можем проранжировать наши гипотезы с помощью априорных шансов. Допустим у нас потерялась зарядка от телефона, где ее искать? Возле кровати за тумбочкой? В доме друзей где мы недавно были в гостях? В холодильнике? Может на Марсе, ведь ее запросто могли украсть марсеане для изучения человеческих технологий? А ведь ее могли украсть и кошки - они кажется давно что-то замышляют!

Любой человек, который принимает решения, оперирует гипотезами. Не важно о чем идет речь: о выборе новой стратегии для крупной корпорации или о том, что съесть сегодня на завтрак. Любой человек, принимающий решения, осознанно или неосознанно говорит на языке вероятностей. Последовательный анализ - это возможность вести с природой диалог. И от качества этого диалога будет зависеть очень многое будь то судьбы корпораций, определяемые выбранными стратегиями или наше самочуствие в течение дня, определяемое тем что мы съели на завтрак.

изображение.png

Комментарии (4)


  1. wiolowan
    27.09.2025 05:46

    Познакомился с последовательным анализом Ввльда в мединституте в 80х, всё время помнил о нём, но в науке всегда спрашивали минимальный обьем выборки и количество наборов реактивов, так что было не до Вальда :)..


  1. GidraVydra
    27.09.2025 05:46

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

    Сколько пишу и рецензирую заявки на гранты, ни разу не видел, даже у самых многословных аппликантов, чтобы в заявке в явном виде указывали такие частности, как метод статанализа экспериментальных даных. Абсолютный максимум, что встречал, это спойлер типа "результаты будут проанализированы методами мат. статистики", без конкретики.

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

    -ставить эксперименты последовательно, каждый раз после этого производя расчеты, хватит экспериментов или надо ещё

    или

    -поставить одновременно в 1.5-2 раза больше экспериментов, но зато точно и заранее известно сколько именно, и обработать результаты один раз после завершения серии

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


  1. 1Fedor
    27.09.2025 05:46

    Если это отношение равно 3, то можно сказать, что гипотеза H_{0} объясняет наблюдаемое количество орлов в 3 раза лучше чем H_{1} и наоборот, если отношение равно 1/3, то H_{1} объясняет наблюдения в два раза лучше чем H_{0}

    Почему в первом случае, когда соотношение 3:1 в три раза, а во втором случае, когда 1:3 в два раза?


    1. uchitel Автор
      27.09.2025 05:46

      Очепятка. Спасибо, поправил.