Предисловие

У блондинки спрашивают: "Какова вероятность, выйдя из дома на улицу, встретить динозавра?". Она отвечает: "50%: или встречу, или не встречу."

Уверен, многие знают этот бородатый анекдот. Казалось бы, ответ блондинки безусловно глупый и высмеивающий её интеллектуальные способности. Но не спешите с выводами, ведь на самом деле, если углубиться, не всё так просто и однозначно. И вот, почему.

Поставленный в анекдоте вопрос не содержит никакой информации об опыте прошлого, которым можно пользоваться для оценки такой вероятности. Конечно, можно додумать ситуацию так называемыми данными по умолчанию, например: в городе проживания, сегодня, и тп. Но ведь с тем же успехом можно и додумать, что никаких сведений нет. То есть, например, человек родился в подземном бункере, жил там лет 30, не выходя ни разу на улицу и не имея никакой возможности получать информацию извне, и тут ему предстоит выйти на улицу, но, чтобы не быть съеденным, он пробует оценить вероятность встретить динозавра, не зная буквально ничего об окружающем мире. Когда не известно абсолютно ничего об оцениваемом событии, то его присутствие или отсутствие оказываются равновероятными, а значит логичнее всего действительно ответить 50%, правда, с поправкой на то, что погрешность такой вероятности также 50%.

Интересные формулировки

Самое интересное начинается, если мы начнём углубляться в эту задачу. Допустим, наш герой из бункера действительно оценил вероятность встретить динозавра как 50%. Он как следует приготовился к этой встрече и вышел. Динозавра не оказалось. На следующий день он вышел вновь - и снова не встретил динозавра. И так n раз. Теперь вопрос можно сформулировать обобщённо: Какова вероятность встретить динозавра, выходя на улицу (n+1)-й раз, если в первые n раз динозавр не встречался? Очевидно, что, чем больше n, тем меньше вероятность встретить динозавра в (n+1)-й раз. Но ведь фактическая вероятность действительно может отличаться от нуля, поэтому динозавра нельзя сбрасывать со счетов. А следовательно нам нужен какой-то способ разумно оценить эту вероятность. Правильный ответ:

$$ p = \frac{1}{n+2} $$

Ниже я приведу математическое обоснование, почему нужно вычислять именно так. А сейчас обратим внимание, что при $ n = 0 $, мы получаем те самые $ p = 50\% $.

Давайте двигаться в формулировках дальше: Какова вероятность встретить динозавра, выходя на улицу (n+1)-й раз, если в первые n раз динозавр встретился k раз? Это уже более более печальный сценарий... Но не следует пока беспокоиться за нашего героя: раз он рассуждает о вероятности динозавра перед (n+1)-м выходом, это значит, что k встреч с ним он уже пережил. Абстрагируясь от нашего сюжета, эта задача выглядит как классическая: событие случилось k раз из n возможностей, какова его вероятность?. Большинство уверенно делит k на n, чтобы оценить вероятность, но так делать не совсем верно, особенно, если значение n не является большим. Правильный ответ:

$$ p = \frac{k+1}{n+2} $$

Очевидно, последняя формулировка является наиболее общей и вбирает в себя выше перечисленные как частные случаи: $ n = 0 $ или $ k = 0 $. Но совсем не очевидно, что такая оценка вероятности действительно лучше, чем наивная $ \frac{k}{n} $. Ниже я постараюсь привести математически строгое обоснование, а также численный эксперимент, который показывает, что оценивать вероятность формулой $ \frac{k}{n} $ неверно, особенно при малых n.

Решение

Решения могут быть даны разные. Я решил подойти со стороны вариационного исчисления. Нам ведь нужно уметь оценивать вероятность p по данным числам n и k. То есть в нашей задаче неизвестной является функция $ f(n, k) $, которая минимизирует функционал. А какой именно - мы сейчас сформулируем.

Так как мы оцениваем вероятность p, то ошибку можно считать как квадрат отклонения от реальной вероятности:

$$ (f(n, k) - p)^{2} $$

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

$$ \sum_{k=0}^{n} C_{n}^{k} p^{k} (1-p)^{n-k} (f(n, k) - p)^{2} $$

Реальная вероятность также может быть любой. Вводим предположение, что она безотносительно n и k является равномерно распределённой. А значит, обобщая ошибку для произвольной вероятности, мы получаем:

$$ E[f] = \int_{0}^{1} dp \sum_{k=0}^{n} C_{n}^{k} p^{k} (1-p)^{n-k} (f(n, k) - p)^{2} $$

Это и есть функционал ошибки. Методом вариационного исчисления можно найти такую функцию $ f(n, k) $, при которой функционал принимает минимальное значение. Чтобы это сделать, сначала упростим наше выражение:

$$ E[f] = \int_{0}^{1} dp \sum_{k=0}^{n} C_{n}^{k} p^{k} (1-p)^{n-k} (f(n, k) - p)^{2} = $$

$$ = \sum_{k=0}^{n} C_{n}^{k} \int_{0}^{1} dp p^{k} (1-p)^{n-k} (f^2 - 2fp + p^2) = $$

$$ = \sum_{k=0}^{n} C_{n}^{k} \left( f^2 \int_{0}^{1} dp p^{k} (1-p)^{n-k} - 2f \int_{0}^{1} dp p^{k+1} (1-p)^{n-k} + \int_{0}^{1} dp p^{k+2} (1-p)^{n-k} \right) = $$

$$ = \sum_{k=0}^{n} C_{n}^{k} \left( f^2 \frac{1}{(n+1) C_n^k} - 2f \frac{1}{(n+2) C_{n+1}^{k+1}} + \frac{1}{(n+3) C_{n+2}^{k+2}} \right) = $$

$$ = \sum_{k=0}^{n} \left( f^2 \frac{1}{n+1} - 2f \frac{k+1}{(n+1)(n+2)} + \frac{(k+1)(k+2)}{(n+1)(n+2)(n+3)} \right) = $$

$$ = \frac{1}{n+1} \sum_{k=0}^{n} \left( f^2 - 2f \frac{k+1}{n+2} + \frac{(k+1)(k+2)}{(n+2)(n+3)} \right) $$

С такой формулой работать уже значительно проще.

Согласно вариационному исчислению, если к оптимальному решению $ f(n, k) $ добавить произвольную малую функцию возмущения $ \delta f(n, k) $, то изменение нашего функционала будет $ O(\delta f^2) $:

$$ E[f + \delta f] - E[f] = O(\delta f^2) $$

$$ \frac{1}{n+1} \sum_{k=0}^{n} \left( 2 f \delta f - 2 \delta f \frac{k+1}{n+2} \right) = O(\delta f^2) $$

$$ \frac{2}{n+1} \sum_{k=0}^{n} \left( f(n, k) - \frac{k+1}{n+2} \right) \delta f(n, k) = O(\delta f^2) $$

Данное равенство должно выполняться для любой функции возмущения $ \delta f(n, k) $, поэтому:

$$ f(n, k) = \frac{k+1}{n+2} $$

Что и требовалось получить.

Хочу ещё раз отдельно подчеркнуть, что мы использовали два принципиальных предположения:

  1. Мы решили в качестве ошибки рассматривать квадрат отклонения (L2).
  2. Мы предположили равномерное распределение реальной вероятности p.

Без этих предположений ответ получился бы другим.

Численный эксперимент

Хочется теперь численно удостовериться, что полученный нами способ оценивать вероятность действительно лучше, чем наивный $ \frac{k}{n} $. Для этого напишем скрипт, который для фиксированного n 10000 раз генерирует случайную равномерно распределённую вероятность p, потом, согласно вероятности p, генерирует количество успешных исходов k, далее производит оценки по функциям $ f_{optimal}(n, k) = \frac{k+1}{n+2} $ и $ f_{naive}(n, k) = \frac{k}{n} $ и накапливает квадратичную ошибку для каждой.

from random import random

def f_optimal(n, k):
    return (k + 1) / (n + 2)

def f_naive(n, k):
    return k / n

if __name__ == "__main__":
    n = 10
    E_optimal = 0.0
    E_naive = 0.0
    for _ in range(10000):
        p = random()  # Random uniform p
        k = sum(int(random() < p) for _ in range(n))  # Generating k according p
        p_optimal = f_optimal(n, k)  # Optimal estimating
        p_naive = f_naive(n, k)  # Naive estimating
        E_optimal += (p_optimal - p)**2
        E_naive += (p_naive - p)**2
    print("E optimal: {:.2f}".format(E_optimal))
    print("E naive: {:.2f}".format(E_naive))
    print("Improvement: {:.2f}%".format(100 * (E_naive / E_optimal - 1)))

Таблица для разных n:

n     E_optimal     E_naive     Imporvement
5 241.81 338.74 40.08%
10 140.77 169.47 20.38%
20 76.89 85.36 11.01%
50 32.03 33.27 3.87%
100 16.39 16.69 1.81%
200 8.34 8.41 0.79%
500 3.36 3.38 0.55%
1000 1.67 1.67 0.28%

Как мы видим, разница колоссальна для небольших n, но нивелируется с ростом n, что логично, ведь обе формулы асимптотически неразличимы при больших n и k.