Теория вероятностей и статистика для инженера

Предыдущие статьи трека были про мир, в котором всё определено: множество либо содержит элемент, либо нет; программа либо останавливается, либо нет; язык либо регулярен, либо нет. Реальная инженерия так не устроена. Запрос иногда отваливается по таймауту. Кэш иногда промахивается. Новая версия чекаута иногда конвертит лучше, а иногда вам просто повезло с выборкой. Диск, скорее всего, доживёт до конца квартала.

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

Практическая ценность здесь необычно высокая. Почти каждая инженерная катастрофа среднего масштаба, которую я видел, содержала внутри ошибку в вероятностном рассуждении: алерт с «точностью 99%», который на самом деле в пяти случаях из шести ложный; A/B-тест, остановленный на третий день, «потому что уже видно»; SLO по среднему времени ответа, за которым прячется десятипроцентный хвост; ретрай, который превратил единичный сбой в лавину, потому что «вероятность отказа мала» считалась в предположении независимости.

Эта статья закрывает такие дыры. Формулы будут строгие, но каждая — с интуицией и работающим кодом.

Карта статьи

Часть 1. Вероятностное пространство

Интуиция: мера на множестве исходов

Вероятность — это мера, ровно в том же смысле, в каком длина отрезка или площадь фигуры есть мера. У вас есть множество всех возможных исходов, и вы раздаёте кускам этого множества числа от 0 до 1 так, чтобы куски не «создавали массу из ничего»: если разрезать кусок на непересекающиеся части, сумма их мер равна мере целого. Полная масса нормирована на единицу.

Всё остальное — техника вокруг этой одной идеи.

Строгое определение (Колмогоров, 1933)

Вероятностное пространство — тройка (Ω, F, P), где

  Ω  — непустое множество элементарных исходов (sample space);
  F  — сигма-алгебра подмножеств Ω (события):
         (1) Ω ∈ F;
         (2) A ∈ F  ⇒  Ω \ A ∈ F               (замкнутость к дополнению);
         (3) A₁, A₂, ... ∈ F  ⇒  ∪ᵢ Aᵢ ∈ F      (замкнутость к счётному объединению);
  P  — функция F → [0, 1] такая, что
         (A1) P(Ω) = 1;
         (A2) P(A) ≥ 0 для всех A ∈ F;
         (A3) для попарно непересекающихся A₁, A₂, ...:
                P(∪ᵢ Aᵢ) = Σᵢ P(Aᵢ)             (счётная аддитивность).

Всё. Три аксиомы. Из них выводится вообще всё остальное:

P(∅) = 0                              (взять Aᵢ = ∅ в A3)
P(Ω \ A) = 1 − P(A)                   (A и его дополнение разбивают Ω)
A ⊆ B  ⇒  P(A) ≤ P(B)                 (монотонность)
P(A ∪ B) = P(A) + P(B) − P(A ∩ B)     (формула включений-исключений)
P(∪ᵢ Aᵢ) ≤ Σᵢ P(Aᵢ)                   (union bound — рабочая лошадка в CS)

Последнее неравенство, union bound, вы будете использовать чаще всех остальных: «вероятность, что хоть что-нибудь пошло не так, не превосходит суммы вероятностей отдельных поломок». Оно не требует ни независимости, ни каких-либо предположений — только аксиомы. Именно поэтому оно работает в анализе рандомизированных алгоритмов, где зависимости между событиями неописуемы.

Зачем нужна сигма-алгебра, а не просто «все подмножества»

Естественный вопрос: почему бы не разрешить в качестве событий любые подмножества Ω? Для конечного или счётного Ω так и делают — F = 2^Ω. Но на непрерывном Ω = [0, 1] это невозможно: существуют неизмеримые множества (конструкция Витали, использующая аксиому выбора), которым нельзя приписать длину, не сломав аддитивность. Про множества и аксиому выбора — Теория множеств: от наивной к аксиоматической.

Практический вывод для программиста: сигма-алгебра — это не бюрократия, а явное указание, какая информация доступна. В теории случайных процессов фильтрация F₀ ⊆ F₁ ⊆ F₂ ⊆ ... буквально означает «что мы знаем к моменту t». Это ровно то же, что версионирование состояния в event-sourced системе: события накапливаются, знание только растёт.

Событие против исхода

Частая путаница. Исход — точка ω ∈ Ω. Событие — множество исходов A ∈ F. Вероятность определена на событиях, а не на исходах. Для непрерывных распределений вероятность любого конкретного исхода равна нулю: P(X = 3.14159...) = 0, что не значит «невозможно». Отсюда классическая ошибка в коде: сравнивать сгенерированное вещественное число на равенство бессмысленно, работать надо с интервалами.

Часть 2. Условная вероятность и Байес

Определение и интуиция

Для P(B) > 0:   P(A | B) = P(A ∩ B) / P(B)

Интуиция: вы сужаете вселенную до B и заново нормируете массу. P(A | B) — какая доля B при этом занята A.

Независимость: A и B независимы, если P(A ∩ B) = P(A)·P(B). Эквивалентно P(A | B) = P(A) при P(B) > 0 — знание B не меняет ваши шансы на A.

Осторожно: попарная независимость не влечёт независимость в совокупности. Классический контрпример: бросаем две честные монеты, A = «первая орёл», B = «вторая орёл», C = «результаты совпали». Любые две из этих трёх независимы, но P(A ∩ B ∩ C) = 1/4 ≠ 1/8. В распределённых системах этот подвох приносит реальные деньги: три реплики «независимо» отказывают с вероятностью 0.1% каждая, вы считаете 10⁻⁹ — а они стоят в одной стойке с одним питанием, и корреляция отказов делает реальную цифру ближе к 10⁻³.

Формула полной вероятности и Байес

Пусть B₁, ..., Bₙ — разбиение Ω (попарно не пересекаются, объединение = Ω), P(Bᵢ) > 0.

Полная вероятность:  P(A) = Σᵢ P(A | Bᵢ) · P(Bᵢ)

Байес:               P(Bₖ | A) = P(A | Bₖ) · P(Bₖ) / Σᵢ P(A | Bᵢ) · P(Bᵢ)

Словами: апостериорное ∝ правдоподобие × априорное. Формула Байеса — это машина, которая переворачивает условие: у вас есть P(данные | гипотеза) (это то, что даёт модель), а нужно P(гипотеза | данные) (это то, что нужно для решения).

Ручной разбор: почему «точный» детектор врёт

Детектор аномальных деплоев: чувствительность (доля пойманных настоящих поломок) 99%, специфичность (доля правильно пропущенных нормальных) 95%. Реально ломают прод 1% деплоев. Пришёл алерт. Какова вероятность, что деплой правда ломает прод?

Считать в вероятностях неудобно — считайте в натуральных частотах на воображаемых 10 000 деплоях:

P(поломка | алерт) = 99 / (99 + 495) = 99/594 = 1/6 ≈ 16.7 %

Пять из шести срабатываний — ложные. Не потому что детектор плохой, а потому что база мала: здоровых деплоев в 99 раз больше, и даже 5% ложных срабатываний по этой огромной базе дают 495 алертов против 99 настоящих. Это base rate fallacy — самая дорогая ошибка интуиции в инженерии.

def posterior(prior, sensitivity, specificity):
    """P(гипотеза | положительный сигнал) по формуле Байеса."""
    tp = sensitivity * prior                 # истинно-положительные
    fp = (1 - specificity) * (1 - prior)     # ложно-положительные
    return tp / (tp + fp)

print(round(posterior(0.01, 0.99, 0.95), 4))   # 0.1667
print(round(posterior(0.20, 0.99, 0.95), 4))   # 0.8319  — база 20%, всё меняется
print(round(posterior(0.01, 0.99, 0.999), 4))  # 0.9091  — та же база, но FPR в 50 раз ниже

Отсюда два инженерных вывода, каждый из которых стоит вписать в чеклист по мониторингу:

1. Precision алерта определяется не его чувствительностью, а произведением частоты ложных срабатываний на объём трафика. Чтобы алерт на редкое событие был полезен, специфичность должна быть экстремальной: 1 − specificity порядка базовой ставки.
2. Комбинируйте слабые сигналы. Два независимых детектора с FPR 5% дают совместный FPR 0.25%, и апостериорная вероятность прыгает с 17% до 80%. Это ровно то, что делает наивный байесовский спам-фильтр: перемножает правдоподобия по словам. Работа Paul Graham «A Plan for Spam» (https://paulgraham.com/spam.html) — исторически первое массовое инженерное применение этой формулы.

Та же арифметика объясняет, почему ML-классификатор с accuracy 99% на несбалансированных данных может быть бесполезен: если позитивов 1%, тривиальная модель «всегда отрицательный» уже даёт 99%. Поэтому в задачах с редкими классами смотрят на precision/recall и PR-AUC, а не на accuracy.

Часть 3. Случайные величины и распределения

Определение

Случайная величина — измеримая функция X: Ω → R,
то есть такая, что для любого x множество {ω : X(ω) ≤ x} принадлежит F.

Функция распределения (CDF):  F_X(x) = P(X ≤ x) — не убывает, непрерывна справа,
                              F(−∞) = 0, F(+∞) = 1.

Дискретный случай: PMF  p(x) = P(X = x),      Σ p(x) = 1
Непрерывный случай: PDF f(x) = F'(x),          ∫ f(x) dx = 1,  P(a<X≤b) = ∫ₐᵇ f

Ключевая мысль: случайная величина — это не число и не «случайность», а функция, переводящая исход в число. «Случайно» само ω; X совершенно детерминирована. Программистская аналогия: Ω — сид генератора, X — чистая функция от сида.

Плотность f(x) — не вероятность. Она может быть больше единицы (равномерное на [0, 0.1] имеет f = 10). Вероятность — это площадь под плотностью.

Дискретные распределения, которые вы встретите

Биномиальные коэффициенты и подсчёт исходов — из Дискретная математика и комбинаторика.

Пуассон заслуживает отдельного слова. Это предел биномиального при n → ∞, p → 0, np → λ: «очень много возможностей, каждая срабатывает очень редко». Ровно модель входящего трафика от большого числа независимых пользователей. Отсюда важное свойство: у Пуассона Var = E. Если ваши счётчики запросов имеют дисперсию сильно больше среднего — трафик не пуассоновский, он пачками (burst), и вся ёмкостная арифметика, построенная на пуассоновском предположении, занижает пики. Это, кстати, основной результат Leland et al. про самоподобие сетевого трафика (https://dl.acm.org/doi/10.1145/190314.190338).

Непрерывные распределения

Отсутствие памяти у экспоненциального распределения — не безобидная деталь. Оно означает, что «пожившая» сущность не стареет: если время до отказа экспоненциально, то компонент, проработавший год, ровно так же надёжен, как новый. Для реального железа это ложь (есть приработка и износ — «ванна»), но это единственное непрерывное распределение с таким свойством, и потому его берут по умолчанию в моделях очередей и в схемах ретраев. Проверяйте, прежде чем брать.

Тяжёлые хвосты: главная практическая тема раздела

Латентность запросов почти никогда не нормальна. Она right-skewed, с длинным хвостом: типичный запрос — 20 мс, но 1 из 1000 — 3 секунды, потому что попал в GC-паузу, промах кэша и ретрай одновременно.

import numpy as np
rng = np.random.default_rng(0)

# логнормальная латентность: медиана ~50 мс, но хвост длинный
lat = rng.lognormal(mean=np.log(50), sigma=0.9, size=1_000_000)

print(f"среднее : {lat.mean():8.1f} мс")
print(f"медиана : {np.median(lat):8.1f} мс")
for q in (50, 90, 99, 99.9):
    print(f"p{q:<5}  : {np.percentile(lat, q):8.1f} мс")

Среднее оказывается заметно больше медианы, а p99.9 — в разы больше p50. Отсюда правило, которое стоит принять как аксиому эксплуатации: среднее время ответа — почти бесполезная метрика. Пользователь страницы, которая делает 20 бэкенд-вызовов, с вероятностью около 1 − 0.99²⁰ ≈ 18% увидит ваш p99 хотя бы раз. Это и есть «tail at scale» из статьи Dean & Barroso (https://research.google/pubs/pub40801/) — обязательное чтение.

Насколько тяжёл хвост, определяется тем, какие моменты вообще существуют. У Парето с α ≤ 2 дисперсия бесконечна, при α ≤ 1 бесконечно и среднее. Это не абстракция: размеры файлов, длины очередей, объёмы трафика на пользователя часто имеют α в районе 1–2.

Часть 4. Ожидание, дисперсия и их ловушки

E[X] = Σ x·p(x)   или   ∫ x·f(x) dx        (если ряд/интеграл абсолютно сходится)

LOTUS:      E[g(X)] = Σ g(x)·p(x)          (не нужно искать распределение g(X)!)
Линейность: E[aX + bY] = aE[X] + bE[Y]     ВСЕГДА, даже для зависимых X, Y
Дисперсия:  Var[X] = E[(X − E X)²] = E[X²] − (E X)²
            Var[aX + b] = a² Var[X]
            Var[X + Y] = Var[X] + Var[Y] + 2·Cov[X, Y]
Ковариация: Cov[X,Y] = E[XY] − E[X]E[Y]
Корреляция: ρ = Cov[X,Y] / (σ_X σ_Y) ∈ [−1, 1]

Линейность ожидания — самый мощный инструмент во всей вероятности, потому что не требует независимости. Классический пример: сколько в среднем корзин останется пустыми при хешировании n ключей в n корзин? Наивно надо считать сложное совместное распределение; с линейностью — вводим индикаторы Iⱼ = «корзина j пуста», E[Iⱼ] = (1 − 1/n)ⁿ ≈ 1/e, и ответ ≈ n/e ≈ 0.37n немедленно. Именно так анализируют хеш-таблицы, Bloom-фильтры и балансировщики.

Три ловушки

Ловушка 1: E[1/X] ≠ 1/E[X]. Неравенство Йенсена: для выпуклой g верно E[g(X)] ≥ g(E[X]). Практически: средняя пропускная способность (RPS) — не обратная величина к среднему времени ответа. Если половина запросов идёт 10 мс, а половина 1000 мс, среднее время 505 мс, но среднее «запросов в секунду на воркер» вовсе не 1/0.505. Метрики, усредняющие отношения, почти всегда неверны — усредняйте числитель и знаменатель отдельно, потом делите.

Ловушка 2: некоррелированность ≠ независимость. Возьмите X ~ Uniform(−1, 1) и Y = X². Тогда Cov[X,Y] = E[X³] − E[X]E[X²] = 0, корреляция нулевая, но Y полностью определяется X. Корреляция Пирсона видит только линейную связь. В мониторинге это значит: нулевая корреляция двух метрик не доказывает их независимость; смотрите scatter plot или взаимную информацию.

Ловушка 3: дисперсия суммы складывается только при нулевой ковариации. Строите SLO на n шардах, считаете Var суммарной нагрузки как сумму дисперсий — а шарды коррелированы через общий трафик, и реальная дисперсия в разы больше.

import numpy as np
rng = np.random.default_rng(1)

x = rng.uniform(-1, 1, 200_000)
y = x**2
print(np.corrcoef(x, y)[0, 1])   # ≈ 0.00x — «независимы»?
print(np.corrcoef(x**2, y)[0, 1])  # 1.0    — а связь идеальная

Ковариационная матрица — это симметричная положительно полуопределённая матрица, и всё, что вы знаете про такие матрицы из Собственные значения, SVD и матричные разложения, применимо здесь: PCA — это спектральное разложение ковариационной матрицы, а «объяснённая дисперсия» — её собственные значения.

Часть 5. Предельные теоремы

Закон больших чисел

Пусть X₁, X₂, ... — независимые одинаково распределённые, E[X] = μ конечно,
X̄ₙ = (X₁ + ... + Xₙ)/n.

Слабый ЗБЧ:  для любого ε > 0   P(|X̄ₙ − μ| > ε) → 0   при n → ∞   (сходимость по вероятности)
Сильный ЗБЧ: P( lim X̄ₙ = μ ) = 1                                  (сходимость почти наверное)

ЗБЧ — это лицензия на измерение: выборочное среднее сходится к истинному. Но он ничего не говорит о скорости. Именно поэтому вопрос «сколько нужно данных» решается не ЗБЧ, а ЦПТ и неравенствами концентрации.

Центральная предельная теорема

Если дополнительно Var[X] = σ² < ∞, то

    (X̄ₙ − μ) · √n / σ  →  N(0, 1)   по распределению.

Эквивалентно: X̄ₙ ≈ N(μ, σ²/n) при больших n.

Смысл: ошибка среднего убывает как 1/√n, и её форма — нормальная, независимо от исходного распределения. Отсюда фундаментальная экономика измерений: чтобы уменьшить погрешность вдвое, нужно вчетверо больше данных. Это и есть причина, почему A/B-тесты на мелкие эффекты требуют миллионов пользователей.

Дополнительная деталь для инженера: ЦПТ означает, что «шум с многих независимых источников выглядит нормальным». Поэтому время ответа сервиса, складывающееся из десятков стадий, могло бы быть нормальным — но не является, потому что стадии не независимы и имеют тяжёлые хвосты.

Когда ЦПТ врёт

ЦПТ требует конечной дисперсии. Если её нет — теорема неприменима, и всё построенное на ней (доверительные интервалы, t-тесты) разваливается.

import numpy as np
from scipy import stats
rng = np.random.default_rng(42)

def means(sampler, n, reps=20_000):
    return sampler((reps, n)).mean(axis=1)

# лёгкий хвост: экспоненциальное, n=30 — уже почти нормально
m_exp = means(lambda s: rng.exponential(1.0, s), 30)
print("exp,   n=30:   skew =", round(float(stats.skew(m_exp)), 3))    # ≈ 0.36

# тяжёлый хвост: Парето α=1.5, дисперсия БЕСКОНЕЧНА
m_par_30   = means(lambda s: rng.pareto(1.5, s) + 1, 30)
m_par_3000 = means(lambda s: rng.pareto(1.5, s) + 1, 3000)
print("pareto n=30:   max =", round(float(m_par_30.max()), 1))        # ≈ 400
print("pareto n=3000: max =", round(float(m_par_3000.max()), 1),
      " skew =", round(float(stats.skew(m_par_3000)), 1))             # skew ≈ 69 (!)

При n = 3000 асимметрия распределения среднего не уменьшилась, а выросла: увеличив выборку в 100 раз, вы не приблизились к нормальности ни на шаг. Инженерный вывод: прежде чем строить доверительный интервал по среднему, посмотрите на хвост. Если данные — выручка на пользователя, длительность сессии, размер файла, то среднее почти наверняка доминируется несколькими наблюдениями, и надо либо винзорировать, либо переходить к квантилям, либо использовать бутстрэп с осторожностью.

Часть 6. Неравенства концентрации: математика рандомизированных алгоритмов

Это тот раздел вероятности, который чаще всего нужен в CS, и который реже всего попадает в вводные курсы.

Марков (X ≥ 0):        P(X ≥ a) ≤ E[X] / a
Чебышёв:               P(|X − μ| ≥ kσ) ≤ 1/k²
Хёфдинг (Xᵢ ∈ [0,1], независимые):
                       P(|X̄ₙ − μ| ≥ ε) ≤ 2·exp(−2nε²)
Чернов (мультипликативная форма, сумма индикаторов, μ = E[S]):
                       P(S ≥ (1+δ)μ) ≤ exp(−δ²μ / (2+δ))

Разница в силе колоссальна. Марков не требует ничего, кроме неотрицательности, и потому очень слаб. Чебышёв требует дисперсии и даёт полиномиальный хвост. Хёфдинг и Чернов требуют независимости и ограниченности — и дают экспоненциальный хвост.

import numpy as np
rng = np.random.default_rng(3)

n, eps = 1000, 0.05
print("Хёфдинг:", round(2*np.exp(-2*n*eps**2), 4))          # 0.0135

x = rng.binomial(1, 0.5, (200_000, n)).mean(axis=1)
print("реально:", float(np.mean(np.abs(x - 0.5) >= eps)))   # ≈ 0.0015

Граница завышена примерно в 9 раз — и это нормально: неравенства концентрации дают гарантию худшего случая, а не оценку. Их ценность в том, что они справедливы для любого распределения из класса, без предположений о нормальности.

Где это работает в реальном коде:

• Оценка размера выборки при отсутствии нормальности. Хёфдинг даёт: чтобы |X̄ − μ| < ε с вероятностью 1−δ, достаточно n ≥ ln(2/δ) / (2ε²). Для ε = 1%, δ = 5% это n ≈ 18 445. Никаких предположений — годится для приёмочного контроля, для оценки доли ошибок в логах, для семплирования аудита.
• Bloom-фильтры и count-min sketch. Оценки ложных срабатываний — прямое применение Маркова и union bound (Cormode & Muthukrishnan, https://dl.acm.org/doi/10.1016/j.jalgor.2003.12.001).
• Балансировка нагрузки. «Power of two choices»: при случайном выборе одной из двух корзин максимальная загрузка падает с Θ(log n / log log n) до Θ(log log n) — доказательство целиком на границах Чернова (Mitzenmacher, https://ieeexplore.ieee.org/document/963420).
• Многорукие бандиты. Границы UCB — это Хёфдинг, повёрнутый в доверительный интервал.
• Амплификация вероятности в рандомизированных алгоритмах. Класс BPP и повторные прогоны — см. Теория сложности: классы P, NP, PSPACE, редукции и полнота.

Часть 7. Статистика: обратная задача

Вероятность: модель известна, что будет с данными? Статистика: данные известны, какова модель? Вторая задача принципиально сложнее — ответ никогда не единственный, и всё, что можно, — это указать оценку вместе с мерой неопределённости.

Оценки и их качество

Оценка (estimator) θ̂ = θ̂(X₁,...,Xₙ) — функция от выборки; сама случайная величина.

Смещение:        bias(θ̂) = E[θ̂] − θ
Состоятельность: θ̂ₙ → θ по вероятности при n → ∞
Среднеквадратичная ошибка:
        MSE(θ̂) = E[(θ̂ − θ)²] = bias(θ̂)² + Var(θ̂)

Последнее разложение — тот самый bias-variance tradeoff, который в ML подают как отдельную идею. Это не отдельная идея, это тождество из двух строк алгебры. И оно объясняет, почему смещённая оценка иногда лучше несмещённой: регуляризация (ridge, L2) вносит смещение, но так сильно снижает дисперсию, что суммарная ошибка падает.

Метод максимума правдоподобия

Правдоподобие:  L(θ) = Π f(xᵢ | θ)        (для наблюдённых данных как функция θ)
Логарифм:       ℓ(θ) = Σ log f(xᵢ | θ)
MLE:            θ̂ = argmax ℓ(θ)

Разберём руками для Бернулли (конверсия в A/B-тесте):

ℓ(p) = Σ [ xᵢ log p + (1 − xᵢ) log(1 − p) ] = k log p + (n − k) log(1 − p),   k = Σxᵢ

dℓ/dp = k/p − (n − k)/(1 − p) = 0
      ⇒ k(1 − p) = (n − k)p
      ⇒ k = np
      ⇒ p̂ = k/n

Оценка — просто доля успехов. Приятно, когда математика подтверждает здравый смысл. Теперь нормальное распределение:

ℓ(μ, σ²) = −n/2·log(2πσ²) − (1/(2σ²))·Σ(xᵢ − μ)²

∂ℓ/∂μ  = 0  ⇒  μ̂ = x̄
∂ℓ/∂σ² = 0  ⇒  σ̂² = (1/n)·Σ(xᵢ − x̄)²      ← делится на n, не на n−1

И вот здесь ловушка: MLE дисперсии смещена вниз. Причина в том, что вы оцениваете разброс вокруг выборочного среднего, а оно само подогнано под данные и лежит ближе к точкам, чем истинное μ. Поправка Бесселя n−1 (одна степень свободы потрачена на μ̂) даёт несмещённую оценку.

import numpy as np
rng = np.random.default_rng(4)

s = rng.normal(0, 1, (100_000, 5))            # 100k выборок по 5 точек, истинная σ² = 1
print("MLE   (÷n)  :", round(float(np.mean(s.var(axis=1))), 4))          # ≈ 0.80
print("Бессель (÷n−1):", round(float(np.mean(s.var(axis=1, ddof=1))), 4)) # ≈ 1.00

0.80 ≈ (n−1)/n = 4/5 — ровно предсказанное смещение. Отсюда ddof=1 в np.std/np.var, если вы оцениваете разброс генеральной совокупности, и ddof=0, если описываете саму выборку. По умолчанию numpy ставит ddof=0 — это регулярный источник тихих багов в аналитическом коде.

Свойства MLE, ради которых его любят: состоятельность, асимптотическая нормальность и асимптотическая эффективность (достигает нижней границы Крамера–Рао). Свойства, из-за которых надо быть осторожным: смещённость при малых n и неинвариантность дисперсии оценки к перепараметризации.

Доверительные интервалы

Интервал [L(X), U(X)] — доверительный интервал уровня 1−α, если

    P( L(X) ≤ θ ≤ U(X) ) ≥ 1 − α   для всех θ.

Случайны ГРАНИЦЫ, а не θ.

Это определение почти всегда толкуют неправильно. «95% CI = [1.2, 3.4]» не значит «с вероятностью 95% истинное значение лежит в [1.2, 3.4]» — истинное значение либо там, либо нет, вероятность здесь 0 или 1. Значит вот что: процедура, применённая к новым данным много раз, накроет истину в 95% случаев.

Рабочая формула для среднего:

CI = x̄ ± t_{1−α/2, n−1} · s / √n           (s — выборочное СКО с ddof=1)

Для доли (Wald, работает плохо на краях):  p̂ ± z · sqrt(p̂(1−p̂)/n)
Лучше: интервал Уилсона или Клоппера–Пирсона.

Насколько это надёжно на реальных данных? Проверим покрытие честным симулятором:

import numpy as np
from scipy import stats
rng = np.random.default_rng(5)

def coverage(dist, n, reps=20_000):
    hits = 0
    for _ in range(reps):
        if dist == "norm":
            x, true = rng.normal(0, 1, n), 0.0
        else:                                  # логнормальное с тяжёлым хвостом
            x, true = rng.lognormal(0, 1.5, n), np.exp(1.5**2 / 2)
        m, s = x.mean(), x.std(ddof=1)
        h = stats.t.ppf(0.975, n - 1) * s / np.sqrt(n)
        hits += (m - h) <= true <= (m + h)
    return hits / reps

print("нормальное,   n=20 :", coverage("norm", 20))   # ≈ 0.948  — как обещано
print("логнормальное, n=20 :", coverage("ln", 20))    # ≈ 0.766  — катастрофа
print("логнормальное, n=200:", coverage("ln", 200))   # ≈ 0.875  — всё ещё не 0.95

Заявленные 95% превращаются в 76% на скошенных данных при n = 20 и не дотягивают даже при n = 200. Если вы строите CI по выручке на пользователя или по латентности — вы систематически переоцениваете свою уверенность. Лечение: бутстрэп (лучше BCa), логарифмирование, квантильные метрики или устойчивые статистики.

Бутстрэп: когда формулы нет

Идея Эфрона (1979): относитесь к выборке как к генеральной совокупности и пересэмплируйте её с возвращением. Распределение статистики по псевдовыборкам приближает её выборочное распределение. Работает для любой статистики — медианы, p95, отношения метрик, AUC — где аналитической формулы просто не существует.

import numpy as np
rng = np.random.default_rng(6)

data = rng.lognormal(0, 1.0, 500)
B = 20_000

idx  = rng.integers(0, len(data), (B, len(data)))    # B псевдовыборок
meds = np.median(data[idx], axis=1)                  # статистика на каждой

print("медиана:", round(float(np.median(data)), 3))
print("95% CI :", np.round(np.percentile(meds, [2.5, 97.5]), 3))

Сложность: O(B·n) по времени и O(B·n) по памяти в векторизованной версии (для больших n считайте партиями). B = 2000 хватает для CI, B ≥ 10 000 — если нужны хвостовые квантили.

Границы применимости: бутстрэп плохо работает для максимума/минимума, для параметров на границе области, для сильно зависимых данных (там нужен block bootstrap) и при очень малых n. Канон — Efron & Tibshirani, «An Introduction to the Bootstrap» (Chapman & Hall, 1993).

Часть 8. Проверка гипотез и A/B-тесты

Строгие определения

H₀ — нулевая гипотеза («эффекта нет»), H₁ — альтернатива.
Статистика T = T(X) — функция от данных с известным распределением при H₀.

p-value = P( T(X_new) не менее экстремальна, чем T(x_obs) | H₀ верна )

Ошибка I рода  (α): отвергли H₀, когда она верна          — ложная тревога
Ошибка II рода (β): не отвергли H₀, когда верна H₁        — пропустили эффект
Мощность = 1 − β

p-value — это не вероятность того, что H₀ верна. Это P(данные | H₀), а не P(H₀ | данные). Переворот требует Байеса и априорной вероятности, которой у частотного подхода нет. Отсюда: p = 0.04 не значит «96% шансов, что эффект реален» — при низкой априорной вероятности гипотезы доля ложных открытий среди значимых результатов легко превышает 50% (Ioannidis, «Why Most Published Research Findings Are False», https://doi.org/10.1371/journal.pmed.0020124).

Выбор теста

Практическая рекомендация по умолчанию: тест Уэлча, а не Стьюдента. Классический t-тест предполагает равные дисперсии в группах — предположение, которое почти никогда не выполняется, когда лечение меняет и среднее, и разброс. Уэлч не требует равенства и почти ничего не теряет, когда они всё-таки равны (Delacre et al., https://doi.org/10.5334/irsp.82).

Мощность и размер выборки: считать ДО эксперимента

import numpy as np
from scipy.stats import norm

def n_per_arm(p, mde, alpha=0.05, power=0.8):
    """Сколько наблюдений на группу нужно для детекции относительного эффекта mde."""
    z_a, z_b = norm.ppf(1 - alpha/2), norm.ppf(power)
    p2   = p * (1 + mde)
    pbar = (p + p2) / 2
    num  = (z_a*np.sqrt(2*pbar*(1-pbar)) + z_b*np.sqrt(p*(1-p) + p2*(1-p2)))**2
    return num / (p2 - p)**2

print(round(n_per_arm(0.10, 0.05)))   # 57 763  — поймать +5% к конверсии 10%
print(round(n_per_arm(0.10, 0.10)))   # 14 751  — эффект вдвое больше → выборка вчетверо меньше
print(round(n_per_arm(0.02, 0.10)))   # 80 681  — редкая конверсия дороже стоит

Обратите внимание на масштабирование n ~ 1/эффект² — прямое следствие ЦПТ. Если у вас 20 000 пользователей в неделю, а нужно 57 763 на группу, эксперимент на +5% физически невозможен за разумный срок: не запускайте его, а меняйте метрику (более чувствительную) или гипотезу (более крупный эффект).

Самая распространённая ошибка в индустрии — не считать мощность вообще. Тест на 500 пользователей с мощностью 12% даёт p > 0.05 практически всегда, и команда делает вывод «фича не работает». Отсутствие значимости — не доказательство отсутствия эффекта; чтобы утверждать «эффекта нет», нужны тесты эквивалентности (TOST) или узкий доверительный интервал вокруг нуля.

Жизненный цикл эксперимента

Подглядывание: измеряем цену

Запустили тест на 2000 пользователей на группу, но проверяете p-value каждые 200 наблюдений и останавливаетесь, как только увидели p < 0.05. Насколько это плохо? Симуляция при полном отсутствии эффекта:

import numpy as np
from scipy import stats
rng = np.random.default_rng(7)

def false_positive_rate(n=2000, step=200, reps=5000):
    fp_final = fp_peek = 0
    checks = range(step, n + 1, step)          # 10 проверок
    for _ in range(reps):
        a, b = rng.normal(0, 1, n), rng.normal(0, 1, n)   # эффекта НЕТ
        fp_peek  += any(stats.ttest_ind(a[:k], b[:k]).pvalue < 0.05 for k in checks)
        fp_final += stats.ttest_ind(a, b).pvalue < 0.05
    return fp_final/reps, fp_peek/reps

print(false_positive_rate())   # (≈0.052, ≈0.191)

Одна честная проверка в конце: 5.2% ложных срабатываний — примерно обещанные α. Десять подглядываний: 19.1%. Почти каждый пятый «победивший» эксперимент — чистый шум. При непрерывном мониторинге доля стремится к 100%: при достаточном терпении вы дождётесь p < 0.05 всегда.

Правильные решения существуют и внедрены в индустрии: групповое последовательное тестирование с alpha spending (границы O’Brien–Fleming), mixture SPRT и always-valid p-values (Johari et al., Optimizely, https://arxiv.org/abs/1512.04922), а также байесовский подход, для которого подглядывание не является нарушением.

Множественные сравнения

Проверяете 20 метрик — при α = 0.05 вероятность хотя бы одного ложного срабатывания 1 − 0.95²⁰ ≈ 64%. Это union bound, работающий против вас.

import numpy as np
from scipy.stats import false_discovery_control
rng = np.random.default_rng(8)

# 95 метрик без эффекта + 5 с настоящим эффектом
p = np.concatenate([rng.uniform(0, 1, 95), rng.uniform(0, 0.001, 5)])

print("наивно  (p<0.05)   :", int((p < 0.05).sum()))                          # 6  (1 ложное)
print("Бонферрони (p<α/m) :", int((p < 0.05/len(p)).sum()))                   # 3  (потеряли реальные)
print("Бенджамини–Хохберг :", int((false_discovery_control(p) < 0.05).sum())) # 5  (ровно то, что надо)

• Бонферрони контролирует FWER (вероятность хотя бы одной ошибки) — консервативен, годится когда любая ложная тревога дорога (регуляторика, безопасность).
• Бенджамини–Хохберг контролирует FDR (ожидаемую долю ложных среди отвергнутых) — правильный выбор для скрининга сотен метрик и для guardrail-дашбордов (https://doi.org/10.1111/j.2517-6161.1995.tb02031.x).

Часть 9. Байесовский подход

Частотный: θ — фиксированное неизвестное число, случайны данные. Байесовский: θ — случайная величина с распределением, отражающим ваше знание; данные фиксированы после наблюдения.

p(θ | data) ∝ p(data | θ) · p(θ)
апостериорное ∝ правдоподобие × априорное

Сопряжённость — когда апостериорное принадлежит тому же семейству, что и априорное, и обновление сводится к арифметике. Для конверсий это пара Beta–Binomial:

Априорное:      θ ~ Beta(a, b)
Данные:         k успехов из n
Апостериорное:  θ ~ Beta(a + k, b + n − k)

Beta(1,1) — равномерное «не знаю ничего». Beta(30, 970) — «в среднем 3%, и я довольно уверен». Обновление — сложение двух чисел, то есть это можно делать онлайн, на каждом событии, без пересчёта по всей истории.

import numpy as np
rng = np.random.default_rng(9)

a_A, b_A = 1 + 120, 1 + (4000 - 120)     # контроль:  120/4000 = 3.0%
a_B, b_B = 1 + 152, 1 + (4000 - 152)     # вариант B: 152/4000 = 3.8%

sA = rng.beta(a_A, b_A, 200_000)
sB = rng.beta(a_B, b_B, 200_000)

print("P(B лучше A)      :", round(float((sB > sA).mean()), 4))
lift = (sB - sA) / sA
print("медианный подъём  :", f"{np.median(lift):.1%}")
print("95% кредибл-интервал:", [f"{v:.1%}" for v in np.percentile(lift, [2.5, 97.5])])
print("ожидаемая потеря при выборе B:", f"{np.maximum(sA - sB, 0).mean():.6f}")

Что здесь важно: P(B лучше A) — это та самая величина, которую все пытаются вычитать из p-value, но не могут. Кредибл-интервал допускает прямое толкование «с вероятностью 95% истинный подъём лежит здесь» — в отличие от доверительного. Цена: нужно задать априорное, и на маленьких данных оно влияет на ответ (что скорее честность, чем недостаток: слабоинформативное априорное защищает от абсурдных оценок вроде «конверсия 100%» по двум наблюдениям).

Метрика expected loss («сколько конверсии я потеряю в среднем, если приму неверное решение») — практически лучший критерий остановки, чем любой p-value: он выражен в единицах бизнеса и напрямую сравнивается с порогом.

Thompson sampling: эксплуатация вместо тестирования

Классический A/B-тест половину трафика гарантированно тратит на худший вариант. Многорукий бандит перераспределяет трафик по ходу дела.

import numpy as np
rng = np.random.default_rng(10)

true_p = np.array([0.10, 0.12])           # B лучше на 2 п.п., мы этого не знаем
a = np.ones(2); b = np.ones(2)            # Beta(1,1) — без априорных предпочтений
pulls = np.zeros(2)

for _ in range(20_000):
    theta = rng.beta(a, b)                # сэмплируем «веру» в каждую руку
    k = int(np.argmax(theta))             # играем ту, что выглядит лучшей ЭТОТ раз
    r = rng.random() < true_p[k]
    a[k] += r; b[k] += 1 - r
    pulls[k] += 1

print("показов:", pulls)                            # ≈ [580, 19420]
print("потеряно конверсий:", round(pulls[0]*0.02))  # ≈ 12 вместо 200 у равного сплита

Равный сплит на 20 000 показов потерял бы 10 000 · 0.02 = 200 конверсий. Thompson sampling — около 12. Плата: выборка неслучайна во времени, что ломает наивный частотный анализ и усложняет учёт сезонности и дрейфа. Практическое правило: бандиты — для оптимизации (какой креатив показать), классический A/B — для принятия решений о продукте, где нужна чистая несмещённая оценка эффекта. Каноническая ссылка — Russo et al., «A Tutorial on Thompson Sampling» (https://arxiv.org/abs/1707.02038). Стратегические аспекты выбора между вариантами при наличии конкурирующих агентов — в Теория игр: равновесия, механизмы, приложения в распределённых системах.

Часть 10. Инженерная практика: перцентили, которые нельзя усреднять

Одна ошибка в наблюдаемости стоит отдельного раздела, потому что она есть примерно у всех.

Перцентили не аддитивны. Нельзя усреднять p99 по инстансам, по минутам, по шардам. Среднее из двух p99 — это не p99 объединённой выборки, и ошибка не мала.

import numpy as np
rng = np.random.default_rng(11)

a = rng.lognormal(0, 1.0, 100_000)      # быстрый инстанс
b = rng.lognormal(1.5, 1.0, 100_000)    # медленный (деградировавший) инстанс

p99_a = np.percentile(a, 99)
p99_b = np.percentile(b, 99)
print("p99 A:", round(float(p99_a), 2), " p99 B:", round(float(p99_b), 2))    # 10.09 / 44.89
print("среднее p99   :", round(float((p99_a + p99_b)/2), 2))                  # 27.49  ← ЛОЖЬ
print("истинный p99  :", round(float(np.percentile(np.concatenate([a,b]), 99)), 2))  # 34.84

Ошибка 21% — и это всего на двух инстансах. На сотне разнородных подов расхождение легко достигает разов, и всегда в сторону оптимизма: усреднение перцентилей систематически прячет деградировавшие узлы.

Что делать:

• Хранить гистограммы или скетчи, а не готовые перцентили: HDR Histogram (http://hdrhistogram.org/), t-digest (Dunning, https://arxiv.org/abs/1902.04023), DDSketch (https://arxiv.org/abs/1908.10693). Все три мержатся корректно и дают контролируемую относительную ошибку по квантилям.
• В Prometheus: histogram_quantile() над бакетами _bucket — правильно (агрегируются бакеты); summary-квантили от разных инстансов — принципиально не агрегируются, что честно написано в документации (https://prometheus.io/docs/practices/histograms/).
• Помнить про координированное упущение (coordinated omission): если нагрузочный генератор ждёт ответа перед следующим запросом, медленные ответы систематически недопредставлены, и измеренный p99 может занижать реальность на порядок. Это открытие Gil Tene, встроенное в поправки HdrHistogram.
• При сэмплировании трейсов помнить, что равномерный сэмплинг 1% почти не увидит хвост: нужен tail-based sampling, который решает о сохранении трейса после того, как узнал его длительность.

Bonus-арифметика: страница, делающая k независимых бэкенд-вызовов, показывает пользователю время максимума, а не среднего. Вероятность «попасть в p99 хотя бы раз» равна 1 − 0.99^k: для k = 20 это 18%, для k = 100 — 63%. Оптимизация хвоста важнее оптимизации медианы ровно в этой пропорции. Способы борьбы: hedged requests (продублировать запрос после p95-таймаута), tied requests, деградация с частичным ответом.

Типичные заблуждения

«p-value — вероятность того, что нулевая гипотеза верна». Нет, это P(данные не менее экстремальны | H₀). Чтобы получить P(H₀ | данные), нужен Байес и априорная вероятность.

«p > 0.05, значит эффекта нет». Нет, значит данных не хватило, чтобы его отличить от нуля. Смотрите на доверительный интервал: [−0.1%, +0.2%] и [−15%, +18%] оба «незначимы», но говорят совершенно разное.

«95% CI означает, что параметр лежит внутри с вероятностью 0.95». Нет — параметр не случаен. Случайны границы. Хотите прямого толкования — берите байесовский кредибл-интервал.

«Наблюдения независимы». Почти никогда. Пользователь делает несколько сессий, запросы приходят пачками, реплики в одной стойке. Игнорирование зависимости занижает дисперсию, сужает интервалы и раздувает ложные срабатывания. Лечится кластеризованными стандартными ошибками, рандомизацией на уровне пользователя, block bootstrap.

«Больше данных всегда решает». При тяжёлых хвостах среднее не концентрируется — см. эксперимент с Парето. При систематическом смещении (сломанный сплиттер, ботовый трафик) больше данных лишь уменьшает CI вокруг неверного числа.

«Корреляция ноль — значит связи нет». Пирсон видит только линейную связь. Y = X² даёт нулевую корреляцию при полной зависимости.

«Средняя латентность — хорошая SLO-метрика». Среднее скрывает хвост, а хвост — это то, что чувствуют пользователи. Ставьте SLO на квантили и на долю запросов быстрее порога.

«Регрессия к среднему — это эффект нашей фичи». Выберите худшие 10% пользователей, что-то сделайте, измерьте — станет лучше даже без вмешательства, просто потому что экстремальные значения содержат шум. Отсюда обязательная контрольная группа.

«Хи-квадрат подойдёт для чего угодно». У него есть условия: ожидаемые частоты в ячейках не меньше ~5, наблюдения независимы. На маленьких счётчиках — точный тест Фишера.

«Симуляция Монте-Карло даёт точный ответ». Её собственная погрешность ~ 1/√N. Если вы оцениваете вероятность порядка 10⁻⁶, вам нужны миллиарды прогонов или importance sampling.

Мини-итог

• Вероятность — мера на сигма-алгебре событий; три аксиомы Колмогорова порождают всё остальное, включая union bound — главный рабочий инструмент в анализе алгоритмов.
• Байес переворачивает условие. При малой базовой ставке даже очень точный детектор даёт в основном ложные срабатывания: считайте в натуральных частотах.
• Линейность ожидания работает без независимости и решает половину комбинаторных задач в анализе структур данных. Дисперсия складывается только при нулевой ковариации.
• ЦПТ даёт ошибка ~ σ/√n: вчетверо больше данных ради вдвое меньшей погрешности. При бесконечной дисперсии ЦПТ не работает вовсе — проверяйте хвост перед тем, как доверять среднему.
• Хёфдинг и Чернов дают экспоненциальные гарантии без предположений о нормальности; на них стоят Bloom-фильтры, балансировка, бандиты и оценки размера выборки.
• MSE = bias² + variance — это и есть bias-variance tradeoff. MLE дисперсии смещён; ddof=1 не декоративен.
• Доверительный интервал — свойство процедуры, а не конкретного интервала. На скошенных данных его реальное покрытие может быть 76% вместо 95%; бутстрэп и квантильные метрики надёжнее.
• p-value — это P(данные | H₀). Подглядывание превращает 5% ложных срабатываний в 18%; множественные метрики требуют Бонферрони или Бенджамини–Хохберга.
• Байесовский подход даёт прямые ответы («P(B лучше A)», ожидаемая потеря) и не ломается от подглядывания; Thompson sampling экономит трафик там, где нужна оптимизация, а не вывод.
• Перцентили нельзя усреднять. Храните гистограммы/скетчи, помните про coordinated omission и про то, что хвост умножается на число вызовов на страницу.

Источники

• А. Н. Колмогоров, «Основные понятия теории вероятностей», 1933 — первоисточник аксиоматики.
• Sheldon Ross, «A First Course in Probability» — стандартный вводный курс с большим числом разобранных задач.
• Rick Durrett, «Probability: Theory and Examples», 5-е изд. — https://services.math.duke.edu/~rtd/PTE/pte.html (свободный PDF; строгий, измеримо-теоретический уровень).
• Larry Wasserman, «All of Statistics», Springer 2004 — самый плотный мост от вероятности к статистике и ML.
• Mitzenmacher & Upfal, «Probability and Computing», 2-е изд., Cambridge 2017 — неравенства концентрации и рандомизированные алгоритмы именно для CS.
• Efron & Tibshirani, «An Introduction to the Bootstrap», Chapman & Hall 1993.
• Gelman et al., «Bayesian Data Analysis», 3-е изд. — http://www.stat.columbia.edu/~gelman/book/ (свободный PDF).
• Kohavi, Tang, Xu, «Trustworthy Online Controlled Experiments», Cambridge 2020 — практическая библия A/B-тестирования; см. также https://exp-platform.com/
• Johari et al., «Peeking at A/B Tests: Why It Matters and What to Do About It» — https://arxiv.org/abs/1512.04922
• Deng, Xu, Kohavi, Walker, «Improving the Sensitivity of Online Controlled Experiments by Utilizing Pre-Experiment Data» (CUPED) — https://dl.acm.org/doi/10.1145/2433396.2433413
• Benjamini & Hochberg, «Controlling the False Discovery Rate» — https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
• Dean & Barroso, «The Tail at Scale», CACM 2013 — https://research.google/pubs/pub40801/
• Dunning, «The t-digest: Efficient Estimates of Rank Statistics» — https://arxiv.org/abs/1902.04023
• Russo et al., «A Tutorial on Thompson Sampling» — https://arxiv.org/abs/1707.02038
• Документация: scipy.stats — https://docs.scipy.org/doc/scipy/reference/stats.html ; numpy.random Generator — https://numpy.org/doc/stable/reference/random/generator.html ; Prometheus histograms — https://prometheus.io/docs/practices/histograms/
• ASA Statement on p-Values — https://doi.org/10.1080/00031305.2016.1154108 (официальная позиция Американской статистической ассоциации о том, чего p-value не значит).

Что дальше

Мы всюду считали, что арифметика работает точно: 2*np.exp(-2*n*eps**2) даёт настоящее число, сумма ста миллионов латентностей не теряет разрядов, а решение системы нормальных уравнений — то самое решение. Ничего из этого не верно в машинной арифметике. Следующая статья трека — про то, как числа с плавающей точкой на самом деле устроены, почему наивное вычисление дисперсии в один проход может выдать отрицательное значение, что такое обусловленность и устойчивость, и как выбирать численные схемы так, чтобы результат сохранял смысл.

Численные методы: точность, устойчивость, арифметика с плавающей точкой