import numpy as npimport matplotlib.pyplot as pltimport networkx as nxfrom scipy import optimizerng = np.random.default_rng(2026)def simulate_gw_poisson(lam, generations=8, population_cap=5000, rng=None):# Simulate Galton-Watson generation sizes for Poisson offspring.if rng isNone: rng = np.random.default_rng() x = [1]for _ inrange(generations):if x[-1] ==0: x.append(0)continue nxt = rng.poisson(lam, size=x[-1]).sum() nxt =min(int(nxt), population_cap) x.append(nxt)return np.array(x, dtype=int)def simulate_gw_tree_poisson(lam, generations=6, max_nodes=250, rng=None):# Simulate a finite genealogical tree.if rng isNone: rng = np.random.default_rng() G = nx.DiGraph() levels = {0: [0]} G.add_node(0, generation=0) next_id =1for gen inrange(1, generations +1): current = []for parent in levels.get(gen -1, []):iflen(G) >= max_nodes:break n_children =int(rng.poisson(lam))for _ inrange(n_children):iflen(G) >= max_nodes:break child = next_id next_id +=1 G.add_node(child, generation=gen) G.add_edge(parent, child) current.append(child) levels[gen] = currentifnot current:breakreturn G, levelsdef layered_positions(levels):# Position nodes by generation. pos = {}for gen, nodes in levels.items(): m =len(nodes)if m ==0:continue xs = np.linspace(-1, 1, m) if m >1else np.array([0.0])for x, node inzip(xs, nodes): pos[node] = (x, -gen)return posfig, axes = plt.subplots(1, 3, figsize=(15, 4.8))for ax, lam inzip(axes, [0.8, 1.1, 2]): local_rng = np.random.default_rng(int(1000* lam) +7) G, levels = simulate_gw_tree_poisson( lam, generations=6, max_nodes=180, rng=local_rng ) pos = layered_positions(levels) nx.draw_networkx_edges( G, pos, ax=ax, arrows=False, width=0.8, alpha=0.6 ) nx.draw_networkx_nodes( G, pos, ax=ax, node_size=55 ) generation_sizes = [len(levels.get(g, [])) for g inrange(7)] ax.set_title(rf"$\lambda={lam}$"+"\n"+rf"$X_n={generation_sizes}$" ) ax.axis("off")plt.suptitle("Три реализации процесса Гальтона–Ватсона", y=1.03, fontsize=15)plt.tight_layout()plt.show()
TipТри режима ветвящегося процесса
Главной характеристикой закона размножения является
\[
\mu=\mathbb E\xi.
\]
Уже на уровне симуляций видно три качественно различных режима:
\(\mu<1\) — докритический;
\(\mu=1\) — критический;
\(\mu>1\) — надкритический.
Особенно интересен критический случай: процесс может совершать большие случайные выбросы, хотя в дальнейшем всё равно вырождается с вероятностью \(1\).
TipНасколько большим может оказаться вымерший процесс?
Даже если
\[
q=1,
\]
общее число частиц
\[
Y=\sum_{n=0}^{\infty}X_n
\]
может оказаться очень большим.
Особенно интересно приближение к критическому режиму
\[
\mu\uparrow1.
\]
Сравним
\[
\mu=0.6,\qquad
\mu=0.9,\qquad
\mu=0.99.
\]
Show Python code
import numpy as npimport matplotlib.pyplot as pltdef simulate_total_progeny_poisson( lam, n_sim=5000, max_generations=10000, max_total=2_000_000, seed=0):""" Simulate the total progeny Y = X_0 + X_1 + X_2 + ... for a Galton-Watson process with Poisson(lam) offspring. Since X_{n+1} | X_n = x ~ Pois(lam * x), we can simulate each generation directly. """ rng = np.random.default_rng(seed) totals = np.empty(n_sim, dtype=int)for i inrange(n_sim): x =1 total =1for _ inrange(max_generations):if x ==0:break# Exact conditional distribution:# X_{n+1} | X_n=x ~ Pois(lam * x) x = rng.poisson(lam * x) total += x# Safety cap for extremely large realizationsif total >= max_total: total = max_totalbreak totals[i] = totalreturn totalsdef empirical_survival(values):""" Empirical survival function P(Y >= y). """ values = np.asarray(values) xs, counts = np.unique(values, return_counts=True)# Number of observations >= x survival_counts = np.cumsum(counts[::-1])[::-1] surv = survival_counts /len(values)return xs, surv# ------------------------------------------------------------# Simulation# ------------------------------------------------------------mus = [0.6, 0.9, 0.99]plt.figure(figsize=(9, 5.5))for i, mu inenumerate(mus): totals = simulate_total_progeny_poisson( mu, n_sim=7000, max_generations=10000, max_total=2_000_000, seed=500+ i ) x, surv = empirical_survival(totals) plt.step( x, surv, where="post", linewidth=2, label=rf"$\mu={mu}$" )# ------------------------------------------------------------# Figure formatting# ------------------------------------------------------------plt.xscale("log")plt.yscale("log")plt.xlabel(r"Общее число частиц $Y$")plt.ylabel(r"$P(Y\geq y)$")plt.title("Хвост распределения общего числа частиц")plt.grid( alpha=0.25, which="both")plt.legend()plt.show()
8 Время до вырождения
Определим
\[
\boxed{
T=\inf\{n\ge0:X_n=0\}.
}
\]
Для докритического и критического процессов
\[
P(T<\infty)=1,
\]
но распределение \(T\) зависит от близости к критической точке.
TipЧто происходит около критической точки?
Рассмотрим
\[
\xi\sim\operatorname{Pois}(\lambda)
\]
при
\[
\lambda=0.5,\quad0.8,\quad0.95,\quad1.
\]
Оценим
\[
P(T>n).
\]
Несмотря на почти верное вырождение во всех четырёх случаях, около \(\lambda=1\) процесс может существовать очень долго.
Show Python code
def simulate_extinction_times_poisson( lam, n_sim=10000, max_generations=500, seed=0): local_rng = np.random.default_rng(seed) times = np.full(n_sim, max_generations +1, dtype=int)for i inrange(n_sim): x =1for n inrange(1, max_generations +1): x =int(local_rng.poisson(lam, size=x).sum()) if x >0else0if x ==0: times[i] = nbreakreturn timesmax_n =120n_grid = np.arange(max_n +1)plt.figure(figsize=(9, 5.5))for i, lam inenumerate([0.5, 0.8, 0.95, 1.0]): T = simulate_extinction_times_poisson( lam, n_sim=12000, max_generations=max_n, seed=900+ i ) survival = np.array([(T > n).mean() for n in n_grid]) plt.plot(n_grid, survival, linewidth=2, label=rf"$\lambda={lam}$")plt.yscale("log")plt.xlabel("Поколение $n$")plt.ylabel(r"$P(T>n)$")plt.title("Время до вырождения около критической точки")plt.grid(alpha=0.25, which="both")plt.legend()plt.show()
9 Интерпретации и приложения
NoteГде возникают ветвящиеся процессы?
Несмотря на простоту модели, процессы Гальтона–Ватсона возникают как естественное приближение во многих задачах.
9.1 Эпидемии
На ранней стадии эпидемии один заражённый человек заражает случайное число новых людей.
Если обозначить это число через \(\xi\), то
\[
\mu=\mathbb E\xi
\]
играет роль эффективного коэффициента воспроизводства.
При \(\mu<1\) цепочка заражений быстро прекращается, а при \(\mu>1\) существует положительная вероятность большой вспышки.
9.2 Распространение информации в сетях
Каждый пользователь может передать сообщение случайному числу новых пользователей.
На ранней стадии, пока повторные контакты и циклы редки, такая динамика естественно приближается ветвящимся процессом.
9.3 Финансовые каскады
В моделях финансового заражения дефолт одного участника может вызвать дефолты нескольких других участников.
Если сеть разрежена и каскад находится на ранней стадии, число новых дефолтов, вызванных одним дефолтом, можно приближённо описывать законом размножения \(\xi\).
Условие
\[
\mathbb E\xi>1
\]
указывает на возможность самоподдерживающегося каскада.
9.4 Machine Learning и Data Science
Ветвящиеся идеи возникают при изучении:
каскадов распространения информации;
случайных деревьев;
локального исследования разреженных графов;
случайных сетей;
алгоритмов поиска, порождающих случайное число новых состояний.
Во многих таких моделях процесс Гальтона–Ватсона служит первым приближением локального поведения.
