Глава 7.2. Избранные темы: Дерево Гальтона-Уотсона
Процесс Гальтона-Ватсона — это стохастическая модель, описывающая эволюцию популяции, в которой каждый индивидуум размножается независимо в соответствии с одним и тем же распределением вероятностей. Модель была впервые представлена Ирене-Жюлем Бьенеме в 1845 году, а позднее независимо разработана Фрэнсисом Гальтоном и Генри Уильямом Ватсоном в 1874 году для изучения вымирания фамилий.
Первоначальная мотивация исходила из интереса Гальтона к генеалогическим древам и вопросу о том, сохранятся ли аристократические семьи с течением времени. Ватсон и Гальтон вывели вероятность вымирания фамилии, обнаружив, что фамилии в конечном итоге исчезают. С тех пор модель нашла применение не только в генеалогии, но и в биологии (динамика популяций, деление клеток), физике (ветвящиеся процессы в ядерных реакциях), информатике (анализ древовидных структур данных) и экономике (ветвление фирм и организаций).
Размер популяции в следующем поколении зависит только от текущего размера популяции и закона размножения, поэтому последовательность размеров популяции образует цепь Маркова.
Определения и результаты
Определение 1 (Процесс Гальтона-Ватсона) Пусть \(X_0, X_1, X_2, \ldots\) — последовательность случайных величин (представляющих размеры популяции в каждом поколении). Мы говорим, что \((X_t)_{t \ge 0}\) образует процесс Гальтона-Ватсона, если:
- \(X_0 = 1\) (мы начинаем с одного предка).
- Для каждого \(t \ge 1\), \(X_{t+1} = \sum_{n=1}^{X_t} Z_{t,n}\), где \(Z_{t,n}\) представляет собой число потомков \(n\)-го индивидуума в поколении \(t\).
- Случайные величины \((Z_{t,n})_{t \ge 0, n \ge 1}\) являются н.о.р. на \(\mathbb{N}\): \(\mathbb{P}(Z_{t,n} = k)=p_k\) с \(\sum_{k\ge 0} p_k=1\).
Распределение вероятностей \((p_k)_{k \ge 0}\) называется распределением потомства процесса.
Примечание. Процесс Гальтона-Ватсона — это цепь Маркова на пространстве состояний \(S = \mathbb{N}\), где состояние \(0\) является поглощающим (как только популяция вымирает, она исчезает навсегда). Если \(X_t = 0\), то \(X_{t+1} = X_{t+2} = \cdots = 0\). При условии \(X_t=n\), величина \(X_{t+1}\) представляет собой сумму \(n\) н.о.р. случайных величин с распределением \((p_k)\). Следовательно, определяя \((p^{(n)}_k)\) как \(n\)-ю свёртку \((p_k)\), имеем вероятности перехода для \(X_t\): \[ \mathbb{P}(X_{t+1}=k\mid X_t=n)=p^{(n)}_k \] Заметим, что \((p^{(n)})\) определяется соотношениями \[ \begin{aligned} p^{(1)}_k=p_k \\ p^{(n+1)}_k= \sum_{h} p^{(n)}_{h}p_{k-h} \end{aligned} \]
В частности, важно охарактеризовать вероятность1 \(\zeta=\mathbb{P}(\tau_0<\infty)\), вероятность того, что процесс вымрет. Хотя это не особо полезно для вычислений, можно заметить, что для того чтобы процесс вымер в момент времени \(t\), все индивидуумы в момент времени \(t-1\) должны иметь \(0\) потомков: \[ \zeta= \sum_{t=1}^\infty \mathbb{P}(\tau_0=t)=\sum_{t=1}^\infty \sum_{k=0}^\infty \mathbb{P}(\tau_0=t\mid X_{t-1}=k) \mathbb{P}(X_{t-1}=k)= \sum_{t,k=1}^\infty \mathbb{P}(X_{t-1}=k) p_0^k \] Значение \(\zeta\) определяется значениями распределения потомства \((p_k)\), поэтому мы хотим понять, в частности, когда \(\zeta=1\), а когда \(\zeta<1\). Определим функцию \[ \phi(z) := \sum_{k=0}^\infty p_k z^k, \qquad |z|\le 1 \] которая является производящей функцией распределения потомства. В частности, \[ \phi(0)=p_0, \phi(1)=1, \phi'(1)= \mathbb{E}[Z] = \sum_{k=0}^\infty k p_k \]
Лемма 1 Пусть \(\phi_t(z):=\sum_{k=0}^\infty \mathbb{P}(X_t=k) z^k=\mathbb{E}\left[z^{X_t}\right]\), для \(z\in [0,1]\). Тогда выполняется \[ \begin{cases} & \phi_0(z)=z \\ & \phi_1(z)=\phi(z) \\ & \phi_{t+1}(z)= \phi(\phi_t(z)) \end{cases} \] Более того, \(\phi_t\) является неотрицательной, неубывающей, выпуклой и аналитической при \(0\le z<1\).
Доказательство. Тождества для \(\phi_0\) и \(\phi_1\) следуют из определения. При условии значения \(Z_{1}\) имеем \[ \begin{aligned} \phi_{t+1}(z) & =\sum_{k=0}^\infty \mathbb{P}(X_{t+1}=k) z^k= \sum_{k=0}^\infty \sum_{j=0}^\infty \mathbb{P}(X_{t+1}=k\mid X_{1}=j) \mathbb{P}(X_1=j) z^k \\ & =\sum_{j=0}^\infty p_j \left(\sum_{k=0}^\infty \mathbb{P}(X_t=k) z^k\right)^j=\phi(\phi_t(z)) \end{aligned} \] Альтернативно, используя условное математическое ожидание, \(\mathbb{E}[z^{X_{t+1}}]=\mathbb{E}[\mathbb{E}[z^{X_{t+1}} \mid X_1]]=\mathbb{E}[\mathbb{E}[z^{X_{t}}]^{X_1}]=\phi(\phi_{t}(z))\), поскольку при условии \(X_1\), \(X_{t+1}\) является суммой \(X_1\) независимых процессов Гальтона-Ватсона.
Будучи степенным рядом с неотрицательными коэффициентами, \(\phi\) является неотрицательной, неубывающей, выпуклой и аналитической при \(|z|< 1\). Следовательно, её степени \(\phi\circ\phi\circ \ldots \circ \phi\) также неотрицательны, неубывающие, выпуклые и аналитические.
Теорема 1 Пусть \((X_t)_{t \ge 0}\) — процесс Гальтона-Ватсона с распределением потомства, имеющим производящую функцию вероятностей \(\phi\) и среднее \(\mu = \mathbb{E}[Z] = \phi'(1)\in [0,\infty]\). Тогда \(\zeta\) является наименьшей неподвижной точкой \(\phi\), то есть наименьшим решением уравнения \[ z = \phi(z), \qquad z\in [0,1] \] В частности:
- Если \(p_0=0\), тривиально \(\zeta=0\).
- Если \(p_0>0\) и \(\mu \le 1\), то популяция в конечном итоге вымирает п.н., то есть \(\zeta=1\).
- Если \(\mu > 1\), то \(\zeta<1\) и \[ \mathbb{P}(\lim_{t\to \infty} X_t=\infty) =1-\zeta \]
Доказательство. Если \(X_1=k\), то вымирание произойдёт тогда и только тогда, когда вымрут \(k\) (независимых) деревьев, происходящих от каждого из \(k\) индивидуумов в момент времени \(1\), то есть с вероятностью \(\zeta^k\). Таким образом, при условии \(X_1\): \[ \zeta=\sum_{k=0}^\infty \mathbb{P}(\tau_0<\infty\mid X_1=k) \mathbb{P}(X_1=k)= \sum_{k=0}^\infty p_k \zeta^k=\phi(\zeta) \] То есть \(\zeta\) является неподвижной точкой \(\phi\). Нам нужно проверить, что она является минимальной.
Пусть \(\zeta_t = \mathbb{P}(X_t = 0)=\phi_t(0)\) — вероятность вымирания к поколению \(t\). Тогда \(\zeta_t\) неубывающая и \(\zeta = \lim_{t \to \infty} \zeta_t\), в силу непрерывности вероятностных мер на монотонных последовательностях. Пусть \(\xi\) — любая другая неподвижная точка, \(\xi=\phi(\xi)\), тогда, поскольку \(\phi_t(z)\) неубывает по \(z\): \[ \zeta_t=\phi_t(0)\le \phi_t(\xi)=\xi \] Следовательно, переходя к пределу, \(\zeta\le \xi\).
Если \(p_0=0\), то \(\phi(0)=0\) и, следовательно, \(\zeta=0\) является минимальной неподвижной точкой, что доказывает пункт 1.
Заметим, что \(\xi=1\) всегда является неподвижной точкой, так как \(1=\phi(1)\). Однако, поскольку \(\phi\) неубывающая и выпуклая, предположение \(p_0>0\) влечёт, что существует не более одной другой неподвижной точки. Действительно, из элементарного математического анализа:
- если \(\phi'(1)>1\), существует ровно одна другая неподвижная точка \(\zeta<1\).
- если \(\phi'(1)\le 1\), то \(1\) является единственной неподвижной точкой, следовательно \(\zeta=\xi=1\).
Далее проверим, что \(\mathbb{P}(\lim_{t\to \infty} X_t=\infty)=1-\zeta\). Для каждого \(t\) такого, что \(X_t\le k\), существует вероятность не менее \(p_0^k\) вымирания на следующем шаге. Поэтому \[ \mathbb{P}(\tau_0<\infty \mid |\{t\st X_t\le k\}|=\infty)=1 \] Следовательно, на множестве \(\{\tau_0=\infty\}\) имеем, что \(X_t\le k\) лишь конечное число раз п.н. для любого \(k\), и, таким образом, \(X_t\) расходится п.н. на \(\{\tau_0=\infty\}\).
Базовые примеры
Пример 1 (Бинарное деление) Рассмотрим бактерию, имеющую вероятность \(p_0\in (0,1)\) умереть до репликации, вероятность \(p_2\in (0,1)\) реплицироваться в две точные копии (митоз) и вероятность \(p_1=1-p_0-p_2\) претерпеть митоз, при котором одна из двух копий умирает немедленно (так что чистое потомство имеет размер \(1\)). Распределение потомства имеет \(p_k=0\) для \(k\ge 3\) и \(\mu=1-p_0+p_2\). Чтобы найти вероятность вымирания \(\zeta\), нам нужно решить уравнение \[ z = \phi(z) \equiv p_0 + p_1 z + p_2 z^2 \equiv p_0+(1-p_0-p_2) z +p_2 z^2 \] Решениями являются \(z=1\) и \(z=p_0/p_2\). Следовательно \[ \zeta= \begin{cases} p_0/p_2 & \text{если $p_0<p_2$} \\ 1 & \text{если $p_0 \ge p_2$} \end{cases} \] и действительно \(\zeta=1\) тогда и только тогда, когда \(\mu \le 1\).
Пример 2 (Пуассоновское распределение потомства) Предположим, что число потомков подчиняется распределению Пуассона с параметром \(\lambda > 0\): \[ p_k = \frac{e^{-\lambda} \lambda^k}{k!}, \quad k \in \mathbb{N} \] Производящая функция вероятностей имеет вид: \[ \phi(z) = \sum_{k=0}^\infty z^k \frac{e^{-\lambda} \lambda^k}{k!} = e^{\lambda(z-1)} \]
Среднее значение равно \(\mu = \lambda\), поэтому:
- Если \(\lambda \le 1\), вымирание происходит с вероятностью 1.
- Если \(\lambda > 1\), вероятность вымирания \(\zeta\) является наименьшим решением уравнения \(e^{\lambda(z-1)}=z\).
Упражнение 1 Рассмотрим процесс Гальтона-Ватсона, где каждый индивидуум имеет случайное число потомков в соответствии с геометрическим распределением с параметром \(p \in (0,1)\), т.е. \(p_k = p(1-p)^k\) для \(k \ge 0\). Вычислите вероятность вымирания и ожидаемое значение распределения потомства в терминах \(p\).
Имеем \(\phi(z) = \sum_{k=0}^\infty p(1-p)^k z^k = \frac{p}{1 - (1-p)z}\) и \(\mu = \phi'(1) = (1-p)/p\). Решение уравнения \(\phi(z)=z\) приводит к корням \(1\) и \(\frac{p}{1-p}\). Таким образом \[ \zeta = \begin{cases} 1 & \text{если $ p \ge 1/2$} \\ \frac{p}{1-p} & \text{если $p < 1/2$} \end{cases} \] что согласуется со значением \(\mu\).
Упражнение 2 Рассмотрим случайное блуждание \(Y_t\) на \(\mathbb{Z}\) с \(p_{x,x+1}=1-p_{x,x-1}=p \in (0,1/2]\), при условии, что оно начинается в положительной области, то есть предположим \(Y_0=0\), \(Y_1=+1\) для простоты. Пусть \(\tau_0^+\) будет, как обычно, моментом первого возвращения в \(0\), который конечен п.н., так как \(p\le 1/2\), и определим \[ X_n:=| \{0\le t < \tau_0^+ \st Y_{t}=n, Y_{t+1}=n+1 \}| \] как количество раз, когда \(Y_t\) пересекает ребро \((n,n+1)\) до возвращения в \(0\).
Докажите, что \((X_n)\) является процессом Гальтона-Ватсона, найдите распределение потомства и вероятность вымирания \(\zeta\).
Зафиксируем \(n\ge 1\) и пусть \(\sigma_l\) будет моментом \(l\)-го пересечения \((n,n+1)\), таким образом \(Y_{\sigma_l-1}=n\) и \(Y_{\sigma_l}=n+1\). Момент остановки \[ \xi_l:=\inf \{t>\sigma_l : Y_t=n\} \in \mathbb{N} \] является первым возвращением на уровень \(n\) после \(k\)-го пересечения. Число пересечений \((n+1,n+2)\), происходящих в открытом интервале \((\sigma_k,\xi_k)\), равно \[ Z_{n,l}:=|\{\sigma_l <t<\xi_l : Y_t=n+1, Y_{t+1}=n+2\}| \] В силу строгого марковского свойства, величины \(Z_{n,l}\) независимы, и, таким образом, условия Определение 1 выполняются.
Мы интерпретируем пересечение \((n, n+1)\) как «родителя». Начиная с \(n+1\), случайное блуждание должно либо шагнуть в \(n+2\) (создавая потомка), либо шагнуть в \(n\) (умирая). Поскольку \(p_{x,x+1}=p\), вероятность шагнуть в \(n+2\) раньше \(n\) просто равна \(p\) (один шаг вправо), а вероятность шагнуть в \(n\) раньше \(n+2\) равна \(1-p\) (один шаг влево). Если блуждание попадает в \(n+2\), оно возвращается в \(n+1\) п.н., допуская дальнейшие попытки создания «потомства». Таким образом, \(Z_{n,l}\) подчиняется геометрическому распределению с параметром успеха \(p\): \[ \mathbb{P}(Z_{n,l} = k) = p^k (1-p), \quad k \ge 0 \] Таким образом, имеем \(\phi(z)=(1-p)/(1-pz)\), \(\mu=\phi'(1)=p/(1-p)\le 1\), так как \(p\le 1/2\). Следовательно, \(\zeta=1\) для всех \(p\le 1/2\) (если \(p>1/2\), то приведенный выше аргумент неверен, так как возвращение в \(0\) из \(Y_1=1\) не является достоверным).
Моменты и асимптотическое поведение
Далее мы рассмотрим поведение на больших временах, в частности, в надкритическом случае \(\mu>1\).
Утверждение 1 Пусть \((X_t)_{t \ge 0}\) — процесс Гальтона-Ватсона с распределением потомства \((p_k)\).
- Для \(\mu = \mathbb{E}[Z]=\sum_k k p_k\) выполняется \(\mathbb{E}[X_t] = \mu^t\) для \(t \ge 0\).
- Предполагая, что \(\sigma^2 := \operatorname{Var}(Z) = \sum_k (k-\mu)^2 p_k <\infty\), для \(t \ge 0\) выполняется \[ \operatorname{Var}(X_t)= \begin{cases} \sigma^2 \mu^{t-1} \frac{\mu^t - 1}{\mu - 1} & \text{если $\mu \neq 1$} \\ t\sigma^2 & \text{если $\mu = 1$} \end{cases} \]
Доказательство.
\(\mathbb{E}[X_t]=\phi_t'(1)\). По индукции, используя Лемма 1, \(\phi_1'(1)=\phi'(1)=\mu\). \(\phi_{t+1}'(1)= \phi'(\phi_{t}(1)) \phi_{t}'(1)\). Но \(\phi_{t}(1)=1\) и \(\phi_{t}'(1)=\mu^t\) по предположению индукции.
Для дисперсии рассуждаем аналогично, используя \(\operatorname{Var}(X_t)= \mathbb{E}[X_t^2]-\mathbb{E}[X_t]^2 = \phi_t''(1)+\phi_t'(1)-\phi_t'(1)^2\) и применяя индукцию.
Теорема 2 Существует случайная величина \(M\ge 0\) с \(\mathbb{E}[M]<\infty\) и \(\mathbb{P}(M=0)\ge \zeta\), такая что предел \[ \lim_{t\to \infty} X_t \mu^{-t}=M \] существует п.н.
Доказательство. По определению \(X_t\) \[ \mathbb{E}[X_t \mid (X_{t-1},X_{t-2},\ldots,X_0)]= \mu X_{t-1} \] Следовательно, \(M_t:=X_t \mu^{-t}\) является неотрицательным мартингалом. По теореме Дуба о сходимости, предел \(M\) последовательности \(M_t\) существует п.н. и интегрируем. Когда \(\tau_0<\infty\), \(X_t\) в конечном итоге обнуляется, поэтому \(M=0\) на \(\{\tau_0<\infty\}\).
Примечание. Предыдущая теорема устанавливает, что \(X_t \mu^{-t}\) стабилизируется п.н. (возможно, в \(0\)) при \(t\to \infty\): \[ X_t= \mu^{-t}(M+o_t(1)) \] при \(t\to \infty\). Это асимптотическое поведение можно визуализировать с помощью численного моделирования: каждый раз, когда мы сэмплируем процесс Гальтона-Ватсона, мы получаем возможно (если \(\mu>1\)) другой предел \(M\) (т.е. предел \(M\) случаен). Таким образом, для каждой симуляции \(\omega\) мы строим случайную величину \(M(\omega)\), просто беря предел. Если мы хотим оценить распределение \(M\), мы можем просто численно оценить \(\mathbb{P}(M\in A)\), используя метод Монте-Карло: мы подсчитываем долю раз (для достаточно больших \(T\)), когда \(X_T/\mu^T\) попадает в множество \(A\).
Численная справедливость Теорема 2 понятна сразу2.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| components: [viewer, editor]
#| viewerHeight: 800
#| editorHeight: 120
#| layout: vertical
#| id: galtonwatson
from shiny import App, render, ui, reactive, req
import matplotlib.pyplot as plt
import numpy as np
import theme
# ========= 1. UI DEFINITION =========
app_ui = theme.layout(
# 1. Sidebar definition (First argument)
ui.sidebar(
ui.h4("Simulation Controls"),
ui.input_slider(
"n_sims",
"Number of Samples:",
min=10, max=100, value=50, step=10
),
ui.input_slider(
"n_gen",
"Generations (N):",
min=5, max=80, value=40, step=5
),
ui.input_select(
"dist",
"Offspring Distribution:",
choices=["Poisson", "Binomial", "Geometric"]
),
ui.panel_conditional(
"input.dist === 'Poisson'",
ui.input_slider(
"poisson_lambda",
"Expectation λ:",
min=0.1, max=3.0, value=1.5, step=0.1
)
),
ui.panel_conditional(
"input.dist === 'Binomial'",
ui.input_slider(
"binom_n",
"Number of trials n:",
min=1, max=10, value=2, step=1
),
ui.input_slider(
"binom_p",
"Success probability p:",
min=0.0, max=1.0, value=0.5, step=0.01
)
),
ui.panel_conditional(
"input.dist === 'Geometric'",
ui.input_slider(
"geom_p",
"Success probability p:",
min=0.3, max=1.0, value=0.5, step=0.02
)
),
ui.input_action_button("run", "Run Simulation", class_="btn-success"),
title="Galton-Watson Process"
),
# 2. Main Content
ui.output_ui("stats_display"),
ui.output_plot("gw_plot", height="500px"),
title="Galton-Watson Process Convergence"
)
# ========= 2. SERVER LOGIC =========
def server(input, output, session):
theme.setup_plot_style()
sim_results = reactive.Value(None)
def get_dist_params():
"""Extract parameters and compute the expectation μ."""
dist = input.dist()
if dist == "Poisson":
lam = input.poisson_lambda()
return {"type": "Poisson", "lam": lam, "mu": lam}
elif dist == "Binomial":
n = input.binom_n()
p = input.binom_p()
return {"type": "Binomial", "n": n, "p": p, "mu": n * p}
else: # Geometric
p = input.geom_p()
# Mean of Geometric defined on {0, 1, ...} is (1-p)/p
mu = (1.0 - p) / p
return {"type": "Geometric", "p": p, "mu": mu}
def calc_theoretical_extinction(params):
"""
Calculate theoretical extinction probability by solving s = G(s).
Uses Newton-Raphson for Poisson/Binomial and exact formula for Geometric.
"""
mu = params['mu']
if mu <= 1.0:
return 1.0
# Analytic solution for Geometric distribution
if params['type'] == "Geometric":
return 1.0 / mu
# Newton-Raphson for other distributions: f(s) = G(s) - s = 0
s = 0.5 # Initial guess
for _ in range(20):
if params['type'] == "Poisson":
lam = params['lam']
gs = np.exp(lam * (s - 1))
g_prime = lam * gs
elif params['type'] == "Binomial":
n, p = params['n'], params['p']
base = 1 - p + p * s
gs = base ** n
g_prime = n * p * (base ** (n - 1))
else:
return 1.0
f_val = gs - s
f_prime = g_prime - 1
if abs(f_val) < 1e-7:
break
s_new = s - f_val / f_prime
s = s_new
return max(0.0, min(1.0, s))
@reactive.Effect
@reactive.event(input.run, ignore_init=True)
def run_simulation():
n_sims = input.n_sims()
n_gen = input.n_gen()
params = get_dist_params()
# Initialize Population Matrix (using float64)
X = np.zeros((n_sims, n_gen + 1), dtype=np.float64)
X[:, 0] = 1.0
current_pop = np.ones(n_sims, dtype=np.float64)
# Thresholds to switch from discrete sampling to Normal approximation
BINOM_THRESHOLD = 2000000000 # ~2e9
GEOM_THRESHOLD = 1e14
# Safety limit for float precision to stop drawing lines
SAFE_FLOAT_LIMIT = 1e250
try:
for g in range(1, n_gen + 1):
# Mask: only simulate for active populations (positive and finite)
active_mask = (current_pop > 0) & np.isfinite(current_pop)
if not np.any(active_mask):
break
active_pop = current_pop[active_mask]
new_counts = np.zeros_like(active_pop)
if params["type"] == "Poisson":
lam_vec = active_pop * params["lam"]
new_counts = np.random.poisson(lam_vec)
elif params["type"] == "Binomial":
n_total = active_pop * params["n"]
p = params["p"]
large_mask = n_total > BINOM_THRESHOLD
small_mask = ~large_mask
if np.any(small_mask):
n_small = n_total[small_mask].astype(np.int32)
new_counts[small_mask] = np.random.binomial(n_small, p)
if np.any(large_mask):
n_large = n_total[large_mask]
mu_large = n_large * p
std_large = np.sqrt(n_large * p * (1 - p))
sample = np.random.normal(mu_large, std_large)
# Ensure positivity after Normal approximation
new_counts[large_mask] = np.maximum(0.0, np.round(sample))
else: # Geometric
large_mask = active_pop > GEOM_THRESHOLD
small_mask = ~large_mask
if np.any(small_mask):
n_small = active_pop[small_mask]
new_counts[small_mask] = np.random.negative_binomial(n_small, params["p"])
if np.any(large_mask):
k = active_pop[large_mask]
p = params["p"]
mu_geom = k * (1.0 - p) / p
var_geom = k * (1.0 - p) / (p * p)
std_geom = np.sqrt(var_geom)
sample = np.random.normal(mu_geom, std_geom)
# Ensure positivity after Normal approximation
new_counts[large_mask] = np.maximum(0.0, np.round(sample))
# Global positivity enforcement for safety
new_counts = np.maximum(0.0, new_counts)
# Check for numerical danger zone
overflow_mask = new_counts > SAFE_FLOAT_LIMIT
if np.any(overflow_mask):
new_counts[overflow_mask] = np.nan
current_pop[active_mask] = new_counts
X[:, g] = current_pop
# Compute Martingale M_n = X_n / mu^n
mu = params["mu"]
if mu == 0:
M_n = np.zeros_like(X)
M_n[:, 0] = 1
else:
scalings = mu ** np.arange(n_gen + 1)
M_n = X / scalings
# Compute Statistics (ignoring NaNs from overflowed lines)
# Extinction: check last valid value or check if 0 present
# We assume NaN means "exploded", so not extinct.
final_vals = X[:, -1]
extinct_sim = np.sum(final_vals == 0) / n_sims
extinct_theo = calc_theoretical_extinction(params)
sim_results.set({
"M_n": M_n,
"extinction_prob": extinct_sim,
"extinction_theo": extinct_theo,
"mean_final_M": np.nanmean(M_n[:, -1]),
"std_final_M": np.nanstd(M_n[:, -1]),
"params": params,
"n_gen": n_gen,
"n_sims": n_sims
})
except Exception as e:
sim_results.set({"error": str(e)})
@output
@render.ui
def stats_display():
results = sim_results.get()
if not results:
return ui.p("Click 'Run Simulation' to see results.", class_="text-center")
if "error" in results:
return ui.div(ui.h4("Error"), ui.p(str(results['error'])), class_="alert alert-danger")
p = results["params"]
return ui.card(
ui.card_header(ui.h4("Simulation Results", class_="m-0")),
ui.row(
ui.column(4,
ui.h5("Configuration", class_="border-bottom pb-2"),
ui.tags.table(
ui.tags.tr(ui.tags.td(ui.strong("Type:")), ui.tags.td(p['type'])),
ui.tags.tr(ui.tags.td(ui.strong("Expectation (μ):")), ui.tags.td(f"{p['mu']:.3f}")),
class_="table table-sm table-borderless mb-0"
)
),
ui.column(4,
ui.h5("Extinction Prob.", class_="border-bottom pb-2"),
ui.tags.table(
ui.tags.tr(ui.tags.td(ui.strong("Theoretical:")), ui.tags.td(f"{results['extinction_theo']:.4f}")),
ui.tags.tr(ui.tags.td(ui.strong("Simulated:")), ui.tags.td(f"{results['extinction_prob']:.3f}")),
class_="table table-sm table-borderless mb-0"
)
),
ui.column(4,
ui.h5(f"Martingale (M_{results['n_gen']})", class_="border-bottom pb-2"),
ui.tags.table(
ui.tags.tr(ui.tags.td(ui.strong("Empirical mean:")), ui.tags.td(f"{results['mean_final_M']:.3f}")),
ui.tags.tr(ui.tags.td(ui.strong("Empirical std:")), ui.tags.td(f"{results['std_final_M']:.3f}")),
class_="table table-sm table-borderless mb-0"
)
),
class_="align-items-start"
)
)
@output
@render.plot
def gw_plot():
results = sim_results.get()
req(results)
if "error" in results: return
M_n = results["M_n"]
n_gen = results["n_gen"]
n_sims = results["n_sims"]
params = results["params"]
fig, ax = plt.subplots(figsize=(10, 6))
line_width = 1.0 if n_sims > 50 else 1.5
line_alpha = 0.5 if n_sims > 50 else 0.7
x_axis = np.arange(n_gen + 1)
# Matplotlib automatically stops drawing line segments containing NaN
ax.plot(x_axis, M_n.T, alpha=line_alpha, linewidth=line_width)
ax.axhline(y=1.0, color='#e74c3c', linestyle='--', linewidth=1.0, label='E[M_n] = 1')
ax.set_xlabel("Generation (n)", fontsize=11)
ax.set_ylabel(r"$M_n = X_n / \mu^n$", fontsize=11)
crit_type = "Supercritical" if params['mu'] > 1 else "Critical" if params['mu'] == 1 else "Subcritical"
ax.set_title(f"Convergence of {params['type']} Process ({crit_type}, μ={params['mu']:.2f})", fontsize=13)
ax.legend(loc='upper right', frameon=True)
ax.grid(True, linestyle=':', alpha=0.6)
if params['mu'] > 1:
# Robust scaling using nanpercentile to ignore overflowed paths
final_vals = M_n[:, -1]
if np.any(np.isfinite(final_vals)):
y_max = np.nanpercentile(final_vals, 98) * 1.2
if y_max > 0:
ax.set_ylim(0, y_max)
return fig
app = App(app_ui, server)
## file: theme.py
from shiny import ui
import matplotlib.pyplot as plt
def layout(sidebar, *args, **kwargs):
"""
Standard layout wrapper for all Shinylive apps.
Args:
sidebar: A ui.sidebar() object.
*args: Main content elements.
**kwargs: Additional arguments for ui.page_sidebar (e.g., title, fillable).
"""
return ui.page_sidebar(
sidebar,
# Standard CSS across all apps
ui.tags.style("""
:root {
--bs-body-font-family: "Source Sans 3", sans-serif;
--bs-font-monospace: "Source Code Pro", monospace;
}
body {
font-family: var(--bs-body-font-family) !important;
}
.bslib-sidebar-layout .sidebar-content {
font-size: 14px;
}
/* Ensure plots blend in if needed */
.shiny-plot-output {
background-color: transparent;
}
"""),
*args,
**kwargs
)
def setup_plot_style():
"""
Configures Matplotlib to match the web theme (Source Sans 3).
Call this at the top of your server function or global scope.
"""
# Attempt to use Source Sans Pro/3 if available, else fallback to sans-serif
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["Source Sans 3", "Source Sans Pro", "DejaVu Sans", "sans-serif"],
"axes.titlesize": 13,
"axes.labelsize": 11,
"xtick.labelsize": 10,
"ytick.labelsize": 10,
"figure.figsize": (8, 5),
"axes.grid": True,
"grid.alpha": 0.3,
"grid.linestyle": ":",
# Transparent backgrounds to match web theme
"figure.facecolor": "none",
"axes.facecolor": "none",
"savefig.facecolor": "none"
})
Когда \(\mu = 1\) (критический случай) и вымирание происходит с вероятностью \(1\), можно изучать поведение процесса при условии невымирания. Этот результат можно доказать аналитически, используя Лемма 1.
Теорема 3 (Условная предельная теорема) Рассмотрим критический процесс Гальтона-Ватсона (\(\mu= 1\), \(\sigma^2 < \infty\)). Тогда: \[ \mathbb{P}(X_t > 0) \sim \frac{2}{\sigma^2 t} \quad \text{при $t \to \infty$} \]
Кроме того, при условии \(\{\tau_0 >t\}\), нормированный размер популяции \(X_t/t\) сходится по распределению к экспоненциальному распределению: \[ \left(X_t/t \mid \tau_0>t\right) \xrightarrow{d} \mathrm{exp}\left(\frac{2}{\sigma^2}\right) \]
Расширения и современные применения
Обобщения
Процессы Гальтона-Ватсона служат фундаментальными моделями с многочисленными обобщениями, например:
- Многотипные процессы Гальтона-Ватсона: Каждый индивидуум может относиться к одному из нескольких типов, с распределениями потомства, зависящими от типа. Теперь процесс представляет собой векторную цепь Маркова.
- Ветвящиеся процессы с непрерывным временем: Размножение происходит в случайные моменты времени в соответствии с процессом Пуассона, что приводит к цепям Маркова с непрерывным временем.
- Пространственные процессы: Каждому индивидууму приписывается положение в пространстве, и распределение потомства зависит от этого положения.
Хотя исторические корни процесса Гальтона-Ватсона лежат в генеалогии, модель эффективно описывает любую систему, где «объекты» производят свои «копии» независимо.
Эпидемиология и распространение вирусов
В современной эпидемиологии распределение потомства представляет собой число вторичных инфекций, вызванных одним инфицированным индивидуумом. Среднее значение \(\mu\) соответствует базовому репродуктивному числу, \(R_0\). В отличие от стандартных моделей SIR, которые предполагают детерминированное смешивание, ветвящиеся процессы учитывают стохастичность ранних вспышек. Они особенно полезны для моделирования событий «суперраспространения», где распределение потомства \((p_k)\) имеет тяжелые хвосты (например, \(p_0\) велико, но \(p_{k}\) для больших \(k\) существенно отлично от нуля). Это помогает рассчитать вероятность того, что вирус, занесенный в новую популяцию, вымрет естественным образом или вызовет эпидемию.
Полимеразная цепная реакция (ПЦР)
ПЦР — это метод молекулярной биологии, используемый для амплификации одной копии сегмента ДНК в миллионы копий. В идеальном сценарии это бинарный процесс Гальтона-Ватсона с \(p_2=1\) (детерминированное удвоение). В реальности эффективность репликации составляет \(p < 1\). Процесс моделируется как ветвящийся процесс, где каждая цепочка ДНК реплицируется с вероятностью \(p\) или неудачно завершается с вероятностью \(1-p\). Теория ветвящихся процессов используется для оценки дисперсии конечного числа копий, что критически важно для точности количественной ПЦР (кПЦР). Базовая модель Гальтона-Ватсона должна быть уточнена для учета двух физических реалий:
- Насыщение: Ветвление не может расти экспоненциально вечно. По мере приближения популяции \(X_t\) к предельной емкости \(K\) (исчерпание реагентов), среднее число потомства \(\mu(X_t)\) падает ниже \(1\).
- Неопределенность параметров: Эффективность репликации \(p\) неизвестна точно и может зависеть от условий работы оборудования.
Мы моделируем это, помещая априорное распределение на эффективность, например \(p \sim \mathrm{Beta}(\alpha, \beta)\), и предполагая плотностно-зависимое распределение потомства. Для размера популяции \(X_t\), эффективная эффективность становится \(p_{eff}(X_t) = p \cdot (1 - X_t/K)\).
Если мы запускаем \(T\) циклов, для успешного обнаружения необходимо, чтобы конечное количество цепочек \(X_T\) превысило определенное значение \(\bar{x}\). Важный вопрос: Сколько циклов \(T\) мы должны запустить, чтобы обеспечить обнаружение? Другими словами, нас интересует момент первого достижения \(\tau := \inf \{t : X_t \ge \bar{x}\}\). Даже если при больших значениях \(X_t\) имеет место более детерминированное поведение (в силу закона больших чисел), ранняя фаза ветвления и неопределенность в \(p\) вносят вклад в дисперсию случайной величины \(\tau\).
Если \(\pi=\mathrm{Beta}(\alpha, \beta)\) — это априорное распределение для \(p\), вероятность неудачи (ложноотрицательного результата) на цикле \(T\) составляет: \[ \mathbb{P}(\text{False Negative}) = \int_0^1 \mathbb{P}(\tau > T \mid p) \pi(p) dp \] В качестве примера мы можем оценить Интегральную Функцию Распределения (CDF) для \(\tau\). Мы отбираем выборку \(p\) в соответствии с \(\pi\), и для каждого отобранного значения запускаем симуляцию процесса Гальтона-Ватсона с распределением потомства, зависящим от \(X_t\). Для фиксированного допустимого риска ложноотрицательного результата \(\varepsilon>0\), мы затем выбираем время выполнения \(T^*\) такое, что \(\mathbb{P}(\tau \le T^*) \ge 1 - \epsilon\) (например, 99%).
Обратите внимание, что дисперсия в “экспоненциальной фазе” (рассматриваемая как горизонтальное распределение кривых) напрямую переходит в неопределенность количественной оценки. Этот горизонтальный сдвиг соответствует случайной величине \(M\) в мартингейльном пределе \(X_t \approx M \mu^t\). В режиме насыщения определение \(M\) позволяет нам “вычислить в обратном порядке” начальную вирусную нагрузку.
Важно здесь построить логарифмический график:
- Пока \(X_t\) мала по сравнению с предельной емкостью \(X_t \ll K\), мы находимся в классическом режиме Гальтона-Ватсона (распределение потомства постоянно в каждом поколении). Таким образом, имеем \[ \log X_t \simeq \log M + t \log \mu \] Мы видим это ясно на логарифмическом графике ниже, где мы видим по различным симуляциям “ширину” \(\log M\), в то время как (после самых первых итераций) происходит линейный (для \(\log X_t\), экспоненциальный для \(X_t\)) рост с коэффициентом, не зависящим от выборки.
- По мере роста \(X_t\) мы видим, что рост выходит на плато, поскольку мы приближаемся к насыщению и распределение потомства теряет эффективность.
В идеале априорное распределение (значения параметров \(\alpha\) и \(\beta\)) предоставляется производителем оборудования для репликации. Мы можем видеть, как ожидаемое поведение проявляется довольно ясно в симуляции Гальтона-Ватсона. Фактически мы получаем значения, очень похожие на те, что используются в лабораториях для обнаружения низких вирусных нагрузок (поскольку мы начинаем всего с одной копии)3.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| components: [viewer, editor]
#| viewerHeight: 1080
#| editorHeight: 120
#| layout: vertical
#| id: pcr-simulation
from shiny import App, render, ui, reactive, req
import matplotlib.pyplot as plt
import numpy as np
import theme
# ========= 1. UI DEFINITION =========
app_ui = theme.layout(
ui.sidebar(
ui.h4("PCR Parameters"),
ui.input_slider("n_sims", "Number of Runs:", 50, 500, 100, step=50),
ui.h5("Efficiency Prior", class_="mt-3"),
ui.input_slider("alpha", "Alpha (Successes):", 1, 100, 40),
ui.input_slider("beta", "Beta (Failures):", 1, 50, 7),
ui.output_plot("prior_plot", height="150px"),
ui.h5("Reaction Limits", class_="mt-3"),
ui.input_slider("log_capacity", "Carrying Capacity (10^K):", 8, 12, 9, step=0.5),
ui.input_slider("log_threshold", "Detection Threshold (10^x):", 4, 8, 7, step=0.5),
ui.input_action_button("run", "Simulate PCR", class_="btn-primary w-100 mt-3"),
title="PCR Cycle Optimizer"
),
ui.card(
ui.card_header(ui.h4("Simulation Results")),
ui.output_plot("trace_plot", height="400px"),
ui.output_plot("cdf_plot", height="250px"),
ui.output_ui("recommendation_ui")
)
)
# ========= 2. SERVER LOGIC =========
def server(input, output, session):
theme.setup_plot_style()
# Store simulation results
results = reactive.Value(None)
@output
@render.plot
def prior_plot():
"""Small helper plot to visualize the Beta prior in sidebar"""
a, b = input.alpha(), input.beta()
x = np.linspace(0, 1, 200)
# The Beta PDF is: x^(a-1) * (1-x)^(b-1) / B(a,b)
# It is not needed to compute B(a,b) since the box is small
# showing scaling is not important
# B(a,b) = gamma(a)*gamma(b) / gamma(a+b)
# beta_constant = (np.gamma(a) * np.gamma(b)) / np.gamma(a + b)
# Calculate PDF manually
# Add epsilon to x to avoid 0^negative errors, if sliders allowed <1
y = (x**(a - 1) * (1 - x)**(b - 1))
# / beta_constant
fig, ax = plt.subplots(figsize=(4, 1.5))
ax.plot(x, y, color='blue', lw=2)
ax.fill_between(x, y, alpha=0.2, color='blue')
ax.set_title(f"Prior E[p] = {a/(a+b):.2f}", fontsize=10)
ax.set_yticks([])
ax.set_xlim(0, 1)
# Transparent background for sidebar integration
fig.patch.set_alpha(0)
ax.patch.set_alpha(0)
plt.tight_layout()
return fig
@reactive.Effect
@reactive.event(input.run, ignore_init=False) # Run once on startup
def run_simulation():
# 1. Setup Parameters
N = input.n_sims()
a, b = input.alpha(), input.beta()
K = 10 ** input.log_capacity()
threshold = 10 ** input.log_threshold()
T_max = 45
# 2. Vectorized Simulation
# Sample unique efficiency 'p' for each run from Prior
p_base = np.random.beta(a, b, N)
# Initialize population matrix (Sims x Time)
X = np.zeros((N, T_max))
X[:, 0] = 1.0
current_pop = np.ones(N)
# Thresholds for approximations
BINOM_CUTOFF = 1000
for t in range(T_max - 1):
# Calculate density-dependent efficiency
# Clip at 0 to prevent negative efficiency if overshot
p_eff = p_base * (1.0 - current_pop / K)
p_eff = np.maximum(0.0, p_eff)
# Update step
new_copies = np.zeros(N)
# A. Small numbers: Use Binomial (exact branching)
mask_small = current_pop < BINOM_CUTOFF
if np.any(mask_small):
n_small = current_pop[mask_small].astype(int)
p_small = p_eff[mask_small]
new_copies[mask_small] = np.random.binomial(n_small, p_small)
# B. Large numbers: Use Normal Approximation (faster)
mask_large = ~mask_small
if np.any(mask_large):
n_large = current_pop[mask_large]
p_large = p_eff[mask_large]
mu = n_large * p_large
std = np.sqrt(n_large * p_large * (1.0 - p_large))
new_copies[mask_large] = np.random.normal(mu, std)
new_copies = np.maximum(0.0, new_copies)
current_pop += new_copies
X[:, t+1] = current_pop
# 3. Analyze Hitting Times
has_crossed = X >= threshold
# Handle cases that never crossed
did_cross = np.any(has_crossed, axis=1)
# argmax returns first True index
crossing_times = np.argmax(has_crossed, axis=1).astype(float)
crossing_times[~did_cross] = np.nan
results.set({
"X": X,
"times": crossing_times,
"threshold": threshold,
"K": K,
"did_cross": did_cross
})
@output
@render.plot
def trace_plot():
res = results.get()
req(res)
X = res["X"]
threshold = res["threshold"]
K = res["K"]
fig, ax = plt.subplots(figsize=(8, 5))
# Plot a subset of traces (max 100 to avoid clutter)
n_plot = min(100, X.shape[0])
ax.plot(X[:n_plot].T, color='blue', alpha=0.15, lw=1)
# Plot Threshold and Capacity
ax.axhline(threshold, color='red', linestyle='--', linewidth=2, label='Detection Threshold')
ax.axhline(K, color='green', linestyle=':', linewidth=2, label='Saturation (K)')
ax.set_yscale('log')
ax.set_ylim(0.8, K * 2)
ax.set_ylabel("DNA Copies (Log Scale)")
ax.set_xlabel("Cycle Number")
ax.set_title("Stochastic PCR Amplification Traces")
ax.legend(loc='upper left')
ax.grid(True, which="both", ls="-", alpha=0.2)
return fig
@output
@render.plot
def cdf_plot():
res = results.get()
req(res)
times = res["times"]
did_cross = res["did_cross"]
valid_times = times[did_cross]
fig, ax = plt.subplots(figsize=(8, 3))
if len(valid_times) == 0:
ax.text(0.5, 0.5, "No detection occurred", ha='center')
return fig
# Empirical CDF
sorted_times = np.sort(valid_times)
y_vals = np.arange(1, len(sorted_times) + 1) / len(times) # Normalize by TOTAL simulations
ax.step(sorted_times, y_vals, where='post', color='darkorange', lw=2)
# Add 99% line
target_prob = 0.99
idx = np.searchsorted(y_vals, target_prob)
if idx < len(sorted_times):
t_safe = sorted_times[idx]
ax.axvline(t_safe, color='black', linestyle='--')
ax.axhline(target_prob, color='black', linestyle=':')
ax.text(t_safe + 0.5, 0.5, f"T* = {int(t_safe)} cycles\n(99% Confidence)",
verticalalignment='center')
ax.set_ylabel("Detection Prob.")
ax.set_xlabel("Cycle Number")
ax.set_ylim(0, 1.05)
ax.set_xlim(left=0)
ax.grid(True, alpha=0.3)
ax.set_title("CDF of Detection Time (τ)")
plt.tight_layout()
return fig
@output
@render.ui
def recommendation_ui():
res = results.get()
req(res)
times = res["times"]
did_cross = res["did_cross"]
if np.sum(did_cross) == 0:
return ui.div("Signal never reached threshold.", class_="alert alert-warning")
# Calculate T* for 99%
sorted_times = np.sort(times[did_cross])
n_total = len(times)
# We need the 99th percentile of the WHOLE distribution
idx_99 = int(n_total * 0.99)
if idx_99 >= len(sorted_times):
msg = f"Only {len(sorted_times)/n_total:.1%} of runs succeeded. Cannot guarantee 99% detection."
cls = "alert alert-danger"
else:
t_star = int(sorted_times[idx_99])
msg = ui.markdown(f"**Recommendation:** Run protocol for **{t_star} cycles** to ensure < 1% false negative rate.")
cls = "alert alert-success"
return ui.div(msg, class_=cls)
app = App(app_ui, server)
## file: theme.py
from shiny import ui
import matplotlib.pyplot as plt
def layout(sidebar, *args, **kwargs):
"""
Standard layout wrapper for all Shinylive apps.
Args:
sidebar: A ui.sidebar() object.
*args: Main content elements.
**kwargs: Additional arguments for ui.page_sidebar (e.g., title, fillable).
"""
return ui.page_sidebar(
sidebar,
# Standard CSS across all apps
ui.tags.style("""
:root {
--bs-body-font-family: "Source Sans 3", sans-serif;
--bs-font-monospace: "Source Code Pro", monospace;
}
body {
font-family: var(--bs-body-font-family) !important;
}
.bslib-sidebar-layout .sidebar-content {
font-size: 14px;
}
/* Ensure plots blend in if needed */
.shiny-plot-output {
background-color: transparent;
}
"""),
*args,
**kwargs
)
def setup_plot_style():
"""
Configures Matplotlib to match the web theme (Source Sans 3).
Call this at the top of your server function or global scope.
"""
# Attempt to use Source Sans Pro/3 if available, else fallback to sans-serif
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["Source Sans 3", "Source Sans Pro", "DejaVu Sans", "sans-serif"],
"axes.titlesize": 13,
"axes.labelsize": 11,
"xtick.labelsize": 10,
"ytick.labelsize": 10,
"figure.figsize": (8, 5),
"axes.grid": True,
"grid.alpha": 0.3,
"grid.linestyle": ":",
# Transparent backgrounds to match web theme
"figure.facecolor": "none",
"axes.facecolor": "none",
"savefig.facecolor": "none"
})
Анализ алгоритмов
В информатике процессы Гальтона-Ватсона появляются в вероятностном анализе алгоритмов и структур данных. Например, высота случайных бинарных деревьев поиска или глубина рекурсии алгоритма Quicksort могут быть проанализированы путем вложения этих дискретных структур в ветвящиеся процессы с непрерывным временем. «Потомки» здесь соответствуют рекурсивным подзадачам, генерируемым на этапе разбиения.
Теория переноса нейтронов
Хотя это и более старое приложение, оно остается жизненно важным для ядерной безопасности. В ядерном реакторе нейтрон, сталкивающийся с ядром, может вызвать деление, высвобождая больше нейтронов. Это процесс Гальтона-Ватсона, где \(\mu\) — коэффициент критичности. - \(\mu < 1\): Докритический (реакция затухает). - \(\mu = 1\): Критический (стабильная выходная мощность). - \(\mu > 1\): Надкритический (мощность растет экспоненциально, потенциально опасно). Расчеты безопасности опираются на вычисление вероятности вымирания (прекращение цепной реакции деления) в сравнении с экспоненциальным разгоном.