Пусть λ — среднее число распадов в единицу времени; его называют интенсивностью процесса. Как устроено распределение (случайного) числа отметок (распадов атомов) на отрезке времени от 0 до 1?

Поделим этот отрезок времени на очень большое число N частей (пусть 1/N это просто временное разрешение нашего самописца) — так, чтобы вероятностью получить два распада на одном отрезке можно было бы пренебречь. Тогда на каждом отрезке мы подкидываем монетку, на которой "есть распад" выпадает с очень маленькой вероятностью λ/N, и "нет распада" с оставшейся.
Поэтому, например, вероятность того, что не было ни одного распада, равна (1-λ/N)^N, что при N->\infty стремится к e^{-λ}.

А если мы будем смотреть не на отрезке длины 1, а на отрезке времени длины t — то вероятность будет e^{-λt}.

К фиксированному отрезку мы сейчас вернёмся — но вообще отсюда можно увидеть, что мы можем сказать про (случайный, как и всё остальное) момент первого распада T_1: вероятность того, что T_1>t, равна e^{-λt}.
Это — экспоненциальное распределение с параметром λ. Оно обладает тем (естественно возникающим из независимости подбрасываний монетки) свойством, что если мы его ждём и за какое-то время t_0 оно ещё не произошло — то при этом условии распределение оставшегося времени ожидания такое же, каким распределение T_1 было исходно. Очень неприятно ждать автобус, если у него такое распределение времени прихода: вроде как ждали-ждали, а если он ещё не пришёл, то ждать остаётся (в смысле распределения) "столько же"!

(Да, тут коллеги ругаются на то, что на картинке выше точки слишком посередине каждого отрезка — так что вот улучшенная непрерывная версия!)

А с какой вероятностью у нас на этом отрезке будет отмечен ровно один распад? У нас N "делений времени", на которых он может произойти; на каждом вероятность, что он произойдёт, это (λ/N), да ещё нужно домножить на вероятность, что больше нигде распадов не будет, (1-λ/N)^{N-1}. Итого

Ну и вообще, если нас интересует, с какой вероятностью на отрезке [0,1] произошло ровно k распадов — то нужно выбрать места, где они происходят, а это можно сделать C_N^k способами, и умножить это на вероятность каждой отдельной "конфигурации": k распадов (с вероятностью (λ/N) каждый) и N-k не-распадов. Итого получаем:

Потому что последний сомножитель это так и есть примерно e^{-λ}, а произведение первых двух можно переписать:

И все сомножители в правой части, кроме первого, это почти 1.

Итак, мы (чуть-чуть рукомахательно — переходя к пределу при N->\infty) выяснили, что вероятность получить k=0,1,2,... распадов за время 1 равна
p_k = (1/k!) λ^k e^{-λ}.

Определение. Такое распределение называется пуассоновским с параметром λ (и обозначается π(λ)).

И вот отсюда пуассоновское распределение вообще появляется — как распределение "редких событий" — числа успехов в большом-большом числе испытаний, каждое из которых даёт очень маленький вклад, а в среднем получается λ. Стандартный пример из учебника — число запросов за минуту к АТС, если абонентов много, а каждый отдельный абонент решает позвонить именно сейчас с довольно небольшой вероятностью.

Кстати — пуассоновский процесс можно рассматривать и не только на прямой. Скажем, представим себе, что наш радиоактивный материал как-то размещён на столе — а мы по прошествии минуты (или часа) отмечаем положения распавшихся атомов (не будем задумываться, как это можно реализовать технически). Тогда мы получим конечный набор точек — такой, что для маленького кусочка стола вероятность того, что точка там есть, это примерно масса лежащего там вещества*среднее число распадов на единицу массы за выбранный интервал времени.

И дальше можно либо опять разрезать стол на маленькие кусочки — или сказать, что для непересекающихся областей расположения точек в них независимы, а количество точек в области A распределено по Пуассону с параметром "масса вещества в A*среднее число распадов за данное время на единицу массы". И ровно так работает общее определение пуассоновского процесса, построенного по неатомарной мере μ — где мы говорим, что число точек в области A распределено как π(μ(A)), а для непересекающихся областей подмножества в них независимы.

Вернёмся к одномерному процессу. Ещё одна вещь, которую мне про него хочется сказать, это что на него (как раз благодаря его одномерности) можно смотреть двумя разными способами: можно так, как мы описали, спрашивая, "сколько отметок в данном интервале". А можно — смотреть на расстояния между отметками (разности между последовательными точками): первая отметка в (случайный) момент T_1=x_1, второй после первой надо ждать второй случайное время T_2=x_2-x_1, ... — и все T_j независимы (ибо то, что уже распалось/не распалось, на то, что будет дальше, не влияет) и распределены как экспоненциальные с параметром λ: вероятность, что T_j больше t, это e^{-λt}.

А теперь — после всех этих блужданий (но, надеюсь, что-то интересное выяснив по пути), давайте вернёмся к исходным вопросам. И начнём с первого из них — как можно объяснить, что
n! = \int_0^{\infty} x^n e^{-x} dx ?

Давайте на n! поделим. Получим тожество
\int_0^{\infty} (1/n!) x^n e^{-x} dx = 1.
А здесь уже явно прослеживается теория вероятностей! Если интеграл чего-то неотрицательного равен 1 — это должна быть плотность распределения какой-то случайной величины. А какой?
На самом деле — мы уже видели почти то, что нам нужно, только с заменой x на λ и n на k. И ещё это была не плотность, а вероятность. Ну да сейчас допилим!

А именно — пусть у нас пуассоновский процесс с единичной интенсивностью. Давайте посмотрим на плотность распределения (k+1)-го по порядку момента распада. Так вот, это в точности
ρ(x)=(1/n!) x^n e^{-x}.
Потому что — чтобы (k+1)-й распад попал в маленький отрезок от x до x+dx, нужно, чтобы
1) в этом отрезке распад произошёл — вот отсюда вероятность dx
2) до того распадов было ровно k — а вот отсюда всё остальное.

Вот мы ответ на первый вопрос и получили: (1/n!) x^n e^{-x} это плотность распределения n+1-й точки для пуассоновского процесса единичной интенсивности — и поэтому (полная вероятность!) совершенно естественно, что интеграл от неё равен 1.

Кстати — если мы k+1-ю точку x_{k+1} мы знаем, то при этом условии набор предыдущих точек получается таким же, как если бы мы на отрезок [0,x_{k+1}] равномерно кинули k независимых точек, а потом бы их упорядочили по возрастанию (вот, кстати, ещё один способ увидеть факториал в знаменателе). Например, потому что при делении на N интервалов и при подкидывании соответствующих монеток все конфигурации "k успехов, N-k неудач" равноправны.

Теперь начинает проявляться и ответ на второй вопрос: как связаны "классический" и "вероятностный" способы смотреть на вычисление бета-функции?
А именно — давайте в начале классического доказательства поделим оба интеграла на соответствующие факториалы (точнее, пока ещё на гамма-функции):