И равно оно f(n)=(n/e)^n — то есть мы опять поймали экспоненциальную часть приближения.

Вынесем его за интеграл; останется интеграл от
exp(g(x)-g(x_0)).

Вот как выглядит график подынтегральной функции при n=30:

Правда, напоминает гауссовский колокол выше?

На самом деле — он и есть (чем больше n, тем точнее).

Потому что — посмотрим, как ведёт себя функция g(x)-g(n) рядом с точкой x=n, где её значение максимально.

Значение, собственно, у неё там обращается в 0 — мы его вычли.

Производной тоже нет: это точка максимума.

Значит, первый нетривиальный член ряда Тейлора квадратичный;
g''(n) = -1/n,
поэтому

А тогда

На самом деле, приближение выше хорошо работает не только при маленьких y, но при y порядка корня из n — а потом подынтегральная функция уже становится слишком маленькой.

Поэтому гамма-функция Г(n) примерно равна произведению (n/e)^n на интеграл от e^{-y^2/2n}.

Но e^{-y^2/2n} это почти плотность нормального распределения с дисперсией n — только её ещё нужно поделить на \sqrt{2\pi n}.

Поделить — чтобы интеграл был единицей (ибо это полная вероятность).

А значит, до деления он и равнялся \sqrt{2\pi n} — вот мы и видим, что у плотности нормального распределения он в знаменателе, а в формуле Стирлинга в качестве множителя.

Но всё-таки, а откуда в гамма-функции ЦПТ, почему у нас тот самый гауссов интеграл получился?

Ну, почему в числителе минус y^2, понятно: значение мы вычли, линейного слагаемого у ряда Тейлора в точке максимума не бывает, а вторая производная будет отрицательна. Пополам там тоже тейлоровский. А вот почему модуль второй производной оказался равен именно n, а не, скажем, n^2?

Давайте возьмём распределение с плотностью e^{-x} на [0;+\infty) — его (логично) называют экспоненциальным распределением.