Ну и вообще, если нас интересует, с какой вероятностью на отрезке [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 неудач" равноправны.

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

В правой части написано совместная плотность для двух независимых случайных величин; одна с плотностью (1/Г(a)) x^{a-1} e^{-x},

другая — с такой же с заменой a на b. Да — такие распределения (логично) называются гамма-распределениями.

Так вот, когда мы этот интеграл считали, мы нарезали всё на отрезки x+y=s. И оказывалось, что интеграл по одному такому отрезку — это
s^{a+b-1}*B(a,b)/(Г(a)Г(b)).
А второе интегрирование тогда даёт
Г(a+b)*B(a,b)/(Г(a)Г(b)),
и из равенства 1 и получается искомый ответ.
Так вот, на вероятностном языке это переговаривается так:
1) свёртка двух гамма-распределений (a.k.a. распределение суммы независимых случайных величин с такими распределениями) — это гамма-распределение с суммой параметров. Потому что мы, собственно, и получили плотность для s=x+y как гамма-распределение с параметром a+b.
2) условное распределение x при условии x+y=s — это растянутое в s раз одно и то же бета-распределение на отрезке [0,1]: распределение на отрезке [0,1] с плотностью
1/B(a,b) x^{a-1} (1-x).

Так вот — возвращаясь к случаю, когда a=m+1 и b=n+1 целые. В этом случае x распределено как m+1-я точка в пуассоновском процессе с единичной интенсивностью, а y — как n+1-я.

Но тогда s=x+y и впрямь распределено как m+n+2-я точка; и чему бы ни равнялось s — условное распределение первых m+n+1-й точки, как мы уже видели, такое же, как при "равномерном бросании" (только нужно упорядочить набор).

И вот ответ на второй вопрос о том, как связаны два взгляда: в классическом подходе тот интеграл, который возникает после замены и перехода на отрезок [0,1] — это (спасибо пуассоновскому взгляду) интеграл от плотности, которую мы получаем для распределения m+1-й по возрастанию точки из m+n+1 равномерно и независимо выбираемых. А эта плотность пропорциональна x^m*(1-x)^n, с коэффициентом пропорциональности, равном числу вариантов "меньшие m точек, m+1-я, большие n",
(m+n+1)!/(m! n!) —
и это и есть наше исходное вероятностное рассуждение.