Вообще, центральная предельная теорема утверждает вот что. Пусть если мы рассматриваем сумму большого числа "разумных" (не слишком "больших"), более-менее сравнимых и независимых случайных слагаемых (например, считая количество орлов за десять тысяч подбрасываний монетки).

У такой суммы будет какое-то среднее значение; и будут отклонения от него.

Среднее значение на то и среднее, чтобы в среднем отклонения "в плюс" и "в минус" были одинаковы. Но можно посмотреть на среднеквадратичное отклонение (на дисперсию, для тех, кому знакомо это слово).

Так вот, давайте отклонение от среднего перемасштабируем — поделив его на его среднее квадратичное. (Каковое, при одинаково распределённых слагаемых, растёт как корень из n — что мы, собственно, и наблюдали в нашем случае.)

Центральная предельная теорема утверждает, что — что бы мы ни складывали исходно, в пределе получится одно и то же, _универсальное_, распределение: гауссовское, или нормальное, распределение.

Так вот — мы сейчас увидели, что для случая подбрасывания монетки шанс, что отношение k/sqrt{n} попадёт в какой-то интервал [a,b], можно приблизить как интегральную сумму Римана — которая, в свою очередь, стремится к интегралу:

То есть в пределе тут мы видим распределение, у которого есть "линейная плотность"
1/\sqrt{\pi} * e^{-x^2}.

Правда, мы делили отклонение на корень из n, а если посмотреть чуть аккуратнее, то средний квадрат отклонения равен n/2, соответственно, чтобы получить совсем предел "как в ЦПТ", делить надо было на \sqrt{n/2}.

Если распределение выше в корень из двух раз растянуть — получится распределение с плотностью

И именно сходимость к распределению с такой плотностью и утверждает центральная предельная теорема.

===
В рассуждениях выше уже видно, что A, которое оказывается множителем в формуле Стирлинга, появляется в знаменателе у гауссовой плотности.

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

Оказывается, факториал можно представить и интегралом:
n! = \int_0^{\infty} x^n e^{-x} dx.

Собственно, интеграл в правой части называется гамма-функцией — и в него можно подставлять не только целые значения показателя у x:
Г(a) := \int_0^{\infty} x^{a-1} e^{-x} dx;
так что формула выше утверждает, что n!=Г(n+1).

Интегрируя по частям, легко проверить, что Г(a+1)=a Г(a), и по индукции проверить, что действительно
Г(n+1)=n!.

Так вот — а как бы нам оценить, чему равен этот интеграл?

Сначала посмотрим: а где подынтегральная функция f(x)=x^n e^{-x} принимает максимальное значение?

Запишем её как exp(g(x)), где
g(x)= ln f(x) = (n ln x - x);

тогда мы хотим максимизировать g(x). Её производная равна g'(x)=n/x -1, значит, максимальное значение принимается в точке x_0=n.