А мне хотелось правильным показать, что с помощью подхода выше ("ошибаемся — давайте посмотрим, на сколько, и учтём") задача [численного] вычисления тоже делалась.

С большими трудозатратами, быть может, но — делалась, а не стояла непреодолимой стеной.

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

Продолжим, пожалуй?

В прошлый раз мы закончили с формулой Стирлинга на том, что n!~A \sqrt{n} (n/e)^n — но мы ещё не знаем, что A это корень из 2π.

Так вот — оказывается, что корень из 2π, который в формуле Стирлинга в числителе, это "тот же самый" корень из 2π, который у нормального распределения (того, что появляется в центральной предельной теореме) в знаменателе. И есть два способа это увидеть.

Первый — это заодно и способ саму центральную предельную теорему в простейшем случае увидеть.

А именно, давайте зададимся вполне естественным вопросом: если мы подкидываем честную монетку N=2n раз, какие результаты мы будем наблюдать, и с какой вероятностью?

Понятно, что наиболее вероятный результат это n орлов и n решек: вероятность из N подбрасываний выкинуть M орлов равна C_N^M / 2^N, а "центральный" биномиальный коэффициент C_{2n}^n самый большой.
А как такая вероятность начинает уменьшаться?

Скажем, при 10000 подбрасываний — какие количества орлов естественны, а какие неправдоподобны и должны заставить заподозрить, что монетка гнутая (или ещё что-нибудь не так)?
5025? 5320? 6500?

Давайте заметим, что
(C_N^{a+1}) / (C_N^a) = (N-a)/(a+1).
Потому что если в классе из N человек уже выбрано a дежурных, то есть N-a вариантов, как можно назначить (a+1)-го — и каждый набор из (a+1) дежурного получается (a+1) способами.

Сравним теперь C_{2n}^n и C_{2n}^{n+k}: второй отличается от первого умножением на
n/(n+1) * (n-1)/(n+2) * ... * (n-k+1)/(n+k).

Получилось произведение k сомножителей.

А как я уже говорил, "видишь длинное произведение — прологарифмируй!"

Логарифм отношения этих биномиальных коэффициентов равен сумме логарифмов
ln (1- (2j-1)/(n+j) )
по j от 1 до k.

А ln(1+x) примерно равен x. Поэтому будет сумма от
-(2j-1)/(n+j);
да ещё и, если k<<n, то знаменатель можно заменить просто на n. И получается сумма от (2j-1)/n, а сумма нечётных чисел от 1 до (2k-1) это k^2:

То есть вероятность отклонения на k падает со скоростью e^{-k^2/n}.

В частности, характерные отклонения имеют порядок корня из n — и слишком большие отклонения даже такого порядка мы в жизни никогда не увидим (ибо, например, e^{-25}<10^{-10})

Иными словами —