Глава 7.4. Избранные темы: Методы Монте-Карло
Методы Монте-Карло для марковских цепей (MCMC)
Одним из наиболее важных приложений цепей Маркова является моделирование случайных величин (принимающих значения в \(S\)). Пусть дана целевая мера \(m\in \mathcal{P}(S)\), и мы заинтересованы в моделировании случайной величины \(X\) с распределением \(m\), то есть \(\mathbb{P}(X=x)=m_x\) для всех \(x\in S\). Возможная стратегия состоит в том, чтобы запустить цепь Маркова (в идеале неприводимую, апериодическую) на \(S\) с инвариантной мерой \(m\). При моделировании цепи в течение длительного времени \(T\), распределение \(X_T\) будет сходиться экспоненциально быстро к целевой мере \(m\). Конечно, важно (и крайне нетривиально) точно оценить скорость сходимости распределения \(X_T\). Самое главное, при заданной целевой мере \(m\), разработать такие вероятности перехода, которые гарантируют быструю сходимость и низкую вычислительную стоимость. Последнее свойство можно формализовать математически разными способами (например, избегание вычисления некоторых сложных величин, количество обращений к \(m_x\) за шаг, или просто то, что алгоритм моделирования использует ближайших соседей относительно некоторой естественной структуры на графе, где определена целевая мера \(m\)), ни один из которых не является полностью удовлетворительным.
Пример 1 (Имитация отжига) Предположим, мы хотим найти глобальный минимум функции стоимости \(U: S \to \mathbb{R}\). Во многих приложениях (например, в комбинаторной оптимизации) пространство \(S\) очень велико, а у \(U\) много локальных минимумов, что делает алгоритмы локального или детерминированного поиска неэффективными.
Идея имитации отжига состоит в том, чтобы превратить эту задачу оптимизации в задачу моделирования (сэмплирования). Для параметра \(\beta > 0\) (интерпретируемого как обратная температура) определим меру Гиббса: \[ m^{\beta}_x = \frac{1}{Z^\beta} e^{-\beta U(x)}, \qquad Z^\beta = \sum_{y \in S} e^{-\beta U(y)} \] При \(\beta \to \infty\) (температура стремится к нулю) мера \(m^{\beta}\) сосредотачивается на множестве глобальных минимумов \(U\). Если мы сможем получать выборку из \(m^{\beta}\) для очень больших \(\beta\), мы найдем минимум с большой вероятностью (вероятность упустить его экспоненциально мала по \(\beta\)).
Однако прямое моделирование непрактично, например, из-за константы \(Z^\beta\) (вычисление \(Z^\beta\) подразумевает вычисление всех \(U(x)\)). Вместо этого мы запускаем цепь Маркова, разработанную так, чтобы \(m^{\beta}\) была её инвариантной мерой. Будем считать, что \(S\) имеет структуру графа со степенью, ограниченной \(D\) (то есть ни у одной точки нет более \(D\) соседей). Это мотивировано тем фактом, что, хотя \(|S|\) во многих случаях велико, оно имеет структуру (связного, неориентированного) графа с малой степенью (см. Пример 2 ниже; или рассмотрите случай большого ящика \(S= \{-N,\ldots,N\}^d\) в \(\mathbb{Z}^d\), тогда \(|S|=(2N+1)^d\), но \(D=2d\)). Тогда мы возьмём \[ p^\beta_{x,y}= \begin{cases} e^{-\beta (U(y) - U(x))^+}/D & \text{если $x,y$ --- соседи,} \\ 1- \!\!\!\! \sum\limits_{{z \text{ сосед } x}}\!\!\!\! e^{-\beta (U(z) - U(x))^+}/D & \text{если $y=x$,} \\ 0 & \text{иначе.} \end{cases} \tag{1}\] где \(a^+\) обозначает положительную часть вещественного числа \(a\).
Этот выбор соответствует тому, что на данном шаге, когда \(X_t=x\), мы выбираем с вероятностью \(1/D\) заданного кандидата в соседи \(y\) для \(x\). Если \(U(y)\le U(x)\), то \(X_{t+1}=y\). Если \(U(y)>U(x)\), то с вероятностью \(e^{-\beta (U(y)-U(x))}\) мы полагаем \(X_{t+1}=y\), а в противном случае не двигаемся, \(X_{t+1}=x\). Идея состоит в том, что эта процедура, с одной стороны, благоприятствует переходам в сторону меньших значений \(U\), а с другой стороны, позволяет (пусть и с малой вероятностью) «выбираться из ям» вокруг локальных минимумов \(U\). Инвариантность \(m\) для этой случайной динамики является следствием более общего утверждения Утверждение 1.
На практике в так называемом алгоритме имитации отжига эту цепь моделируют с медленно растущим \(\beta\) (так что \(\beta=\beta_t\)) по мере увеличения времени \(t\), что теоретически позволяет цепи избегать локальных минимумов и оседать в глобальном минимуме. Поскольку увеличение \(\beta\) означает снижение температуры со временем, название этого влиятельного метода происходит от термической обработки металлов.
Пример 2 (Задача коммивояжёра) Рассмотрим \(N\) городов, помеченных \(1, \dots, N\), и матрицу \((c_{i,j})_{1\le i,j \le N}\), представляющую стоимость путешествия между городом \(i\) и городом \(j\). Коммивояжёр должен посетить все \(N\) городов и вернуться домой, и он должен решить, каков оптимальный порядок посещения городов. Это означает, что мы ищем перестановку \(\sigma \in \mathcal{S}_N\), которая минимизирует общую стоимость.
Здесь пространство состояний \(S = \mathcal{S}_N\) — это множество всех перестановок \(N\) городов. Функция стоимости \(U \colon S \to \mathbb{R}\) — это полная стоимость тура (мы понимаем суммы по модулю \(N\), так что \(N+1=1\)): \[ U(\sigma) = \sum_{i=1}^{N} c_{\sigma(i), \sigma(i+1)} \] Чтобы применить метод из Пример 1, мы должны определить структуру графа на \(S\). Попробуем разные стратегии.
- Две перестановки связаны, если одна может быть получена из другой путем обмена двух индексов. Если обозначить через \(\pi^{i,j}\sigma\) транспозицию, меняющую местами элементы на позициях \(i\) и \(j\), то соседями \(\sigma\) являются \(D=\binom{N}{2}\) перестановок \((\pi^{i,j}\sigma)_{i< j}\).
- Две перестановки связаны, если одна может быть получена из другой путем обмена двух соседних индексов. Это означает, что соседями \(\sigma\) являются все перестановки, полученные из \(\sigma\) с помощью смежной транспозиции, а именно \(D=N\) перестановок \((\pi^{i,i+1} \sigma)_{i}\).
- Если \(c_{i,j}=c_{j,i}\) симметрична, две перестановки связаны, если одна может быть получена из другой путем обращения порядка посещения подсегмента индексов. Если обозначить через \(\epsilon^{i,j}\sigma\) перестановку, полученную обращением пути между индексом \(i\) и \(j\), то соседями \(\sigma\) являются \(D=\binom{N}{2}\) перестановок \((\epsilon^{i,j}\sigma)_{i< j}\).
Теперь, как и в Пример 1, мы фиксируем большую константу \(\beta>0\) (которая управляет вероятностью принятия ходов «в гору»), и моделируем цепь Маркова с вероятностями перехода, заданными Уравнение 1. Решающее вычислительное преимущество локальных динамик a., b., c. состоит в том, что нам никогда не нужно вычислять полную стоимость \(U(\sigma)\), а только разность \(U(\cdot)\) между соседними перестановками. Она легко вычисляется (см. Рисунок 1) как \[ \delta U:= \begin{cases} U(\pi^{i,j}\sigma) - U(\sigma) = (c_{\sigma(i-1),\sigma(j)} + c_{\sigma(j),\sigma(i+1)} + c_{\sigma(j-1),\sigma(i)}+ c_{\sigma(i),\sigma(j+1)}) - (c_{\sigma(i-1),\sigma(i)} + c_{\sigma(i),\sigma(i+1)} + c_{\sigma(j-1),\sigma(j)}+c_{\sigma(j),\sigma(j+1)}) & \text{для случая a.} \\ U(\pi^{i,i+1}\sigma) - U(\sigma) = (c_{\sigma(i-1),\sigma(i+1)} + c_{\sigma(i+1),\sigma(i)} + c_{\sigma(i),\sigma(i+2)}) - (c_{\sigma(i-1),\sigma(i)} + c_{\sigma(i),\sigma(i+1)}+c_{\sigma(i+1),\sigma(i+2)}) & \text{для случая b.} \\ U(\epsilon^{i,j}\sigma) - U(\sigma) = (c_{\sigma(i-1),\sigma(j)} + c_{\sigma(i),\sigma(j+1)}) - (c_{\sigma(i-1),\sigma(i)} + c_{\sigma(j),\sigma(j+1)}) & \text{для случая c.} \end{cases} \tag{2}\] Таким образом, чтобы оценить Уравнение 1, нам никогда не нужно вычислять \(Z^\beta\) или даже \(U\), только комбинацию из соответственно \(8\), \(6\), \(4\) слагаемых стоимости (для случаев a., b., c.) (независимо от \(N\)). Алгоритм тогда таков:
- Начните с произвольной перестановки, например, тождественной или случайной.
- Выберите два случайных индекса \(i,j\) (для динамики a.) или один случайный индекс \(i\) (для динамики b., так как \(j=i+1\) в этом случае).
- Вычислите значение \(\delta U\) в Уравнение 2.
- Если \(\delta U \le 0\), то выполните шаг (а именно, поменяйте местами \(i\) и \(j\) в случае a., поменяйте местами \(i\) и \(i+1\) в случае b., обратите путь между \(i\) и \(j\) в случае c.). Если \(\delta U>0\), то выполните шаг с вероятностью \(e^{-\delta U}\), а в противном случае не делайте никаких изменений.
- Повторяйте с пункта 2.
Конечно, критический вопрос заключается в том, сколько раз мы должны повторять шаги 2.-5. Это проблема, представляющая большой практический и теоретический интерес. Ответ критически зависит от \(N\), \(\beta\), \(U\) и допустимости возможных ошибок. Во многих ситуациях можно просто взять «достаточно большое \(T\)» и избежать математических сложностей.
В примерах a., b., c. выше, довольно ясно, чего ожидать:
- Случаи a. и c. относительно похожи, у каждой точки есть \(N(N-1)/2\) соседей, и им обычно требуется меньше шагов, чтобы обойти пространство состояний и «найти» точки минимума \(U\). Обычно ожидается, что c. работает лучше, чем a. Динамика b. напротив допускает только \(N\) соседей и имеет более высокий шанс застрять в «локальных минимумах» \(U\) (а именно, чтобы найти путь меньшей стоимости, сначала нужно поменять местами несколько соседних индексов, при которых стоимость возрастает, но эти переходы маловероятны при больших \(\beta\)). В симуляциях мы видим, что b. действительно является алгоритмом, для которого \(U\) убывает медленнее всего с итерациями.
- Однако каждый отдельный шаг b. намного проще вычислительно, и если мы построим график \(U(X_t)\) не как функцию количества шагов, а как функцию процессорного времени, затраченного на моделирование, разрыв между b. и другим алгоритмом значительно сократится.
- В конечном счёте, мы видим, что увеличение \(\beta\) по ходу моделирования является лучшей стратегией, и, безусловно, для работы b. нам просто нужно настраивать \(\beta \equiv \beta_t\) так, чтобы она росла, но медленно по сравнению с a. и b.
#| '!! 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: 1600
#| editorHeight: 100
#| layout: vertical
#| id: tsp_metropolis
import time
import numpy as np
import matplotlib.pyplot as plt
from shiny import App, render, ui, reactive, req
import theme
# ==============================================================================
# 1. CONFIGURATION & DEFINITIONS
# ==============================================================================
CHECK_INTERVAL = 200 # Check time/limit every N steps
MAX_PLOT_POINTS = 2000 # Max points to render per line (downsampling)
COST_FUNCTIONS = {
"Squared Euclidean (L2^2)": lambda x, y: np.sum((x - y)**2),
"Euclidean (L2)": lambda x, y: np.sqrt(np.sum((x - y)**2)),
"Manhattan (L1)": lambda x, y: np.sum(np.abs(x - y)),
}
DEFAULT_N = 30
DEFAULT_SEED = 123
DEFAULT_ITER = 5000
DEFAULT_TIME = 5.0
# ==============================================================================
# 2. CORE LOGIC
# ==============================================================================
def generate_cities(n, seed):
rng = np.random.default_rng(seed)
return rng.random((n, 2))
def compute_cost_matrix(cities, metric_func):
n = len(cities)
mat = np.zeros((n, n))
for i in range(n):
for j in range(n):
if i != j:
mat[i, j] = metric_func(cities[i], cities[j])
return mat
def compute_total_cost(sigma, cost_mat):
n = len(sigma)
cost = 0.0
for i in range(n):
cost += cost_mat[sigma[i], sigma[(i + 1) % n]]
return cost
def run_metropolis(cities, cost_mat, strategy, beta_range, limit_mode, limit_val):
n = len(cities)
rng = np.random.default_rng()
sigma = np.arange(n)
rng.shuffle(sigma)
current_cost = compute_total_cost(sigma, cost_mat)
best_sigma = sigma.copy()
best_cost = current_cost
# Pre-allocate estimate
est_size = int(limit_val) if limit_mode == 'iter' else 10000
costs_hist = [current_cost]
beta_start, beta_end = beta_range
start_time = time.time()
step = 0
keep_running = True
while keep_running:
step += 1
# Annealing progress
if limit_mode == 'iter':
progress = min(1.0, step / limit_val)
else:
now = time.time()
progress = min(1.0, (now - start_time) / limit_val)
beta = beta_start + (beta_end - beta_start) * progress
# Proposal
i = rng.integers(0, n)
j = 0
delta_u = 0.0
# Strategy A: Random Swap
if strategy == 'a':
j = rng.integers(0, n)
while i == j: j = rng.integers(0, n)
i_prev, i_next = (i - 1) % n, (i + 1) % n
j_prev, j_next = (j - 1) % n, (j + 1) % n
if i_next == j: # Adjacent i, i+1
p_im1, p_i, p_j, p_jp1 = sigma[i_prev], sigma[i], sigma[j], sigma[j_next]
delta_u = (cost_mat[p_im1, p_j] + cost_mat[p_j, p_i] + cost_mat[p_i, p_jp1]) - \
(cost_mat[p_im1, p_i] + cost_mat[p_i, p_j] + cost_mat[p_j, p_jp1])
elif j_next == i: # Adjacent j, j+1
p_jm1, p_j, p_i, p_ip1 = sigma[j_prev], sigma[j], sigma[i], sigma[i_next]
delta_u = (cost_mat[p_jm1, p_i] + cost_mat[p_i, p_j] + cost_mat[p_j, p_ip1]) - \
(cost_mat[p_jm1, p_j] + cost_mat[p_j, p_i] + cost_mat[p_i, p_ip1])
else:
p_im1, p_i, p_ip1 = sigma[i_prev], sigma[i], sigma[i_next]
p_jm1, p_j, p_jp1 = sigma[j_prev], sigma[j], sigma[j_next]
term_add = cost_mat[p_im1, p_j] + cost_mat[p_j, p_ip1] + cost_mat[p_jm1, p_i] + cost_mat[p_i, p_jp1]
term_sub = cost_mat[p_im1, p_i] + cost_mat[p_i, p_ip1] + cost_mat[p_jm1, p_j] + cost_mat[p_j, p_jp1]
delta_u = term_add - term_sub
# Strategy B: Adjacent Swap
elif strategy == 'b':
j = (i + 1) % n
i_prev = (i - 1) % n
j_next = (j + 1) % n
p_im1, p_i, p_j, p_jn = sigma[i_prev], sigma[i], sigma[j], sigma[j_next]
delta_u = (cost_mat[p_im1, p_j] + cost_mat[p_j, p_i] + cost_mat[p_i, p_jn]) - \
(cost_mat[p_im1, p_i] + cost_mat[p_i, p_j] + cost_mat[p_j, p_jn])
# Strategy C: 2-Opt (Reverse Segment)
elif strategy == 'c':
j = rng.integers(0, n)
# Avoid neighbors for 2-opt to be meaningful (though math holds)
while i == j or j == (i - 1) % n or j == (i + 1) % n:
j = rng.integers(0, n)
i_prev = (i - 1) % n
j_next = (j + 1) % n
p_im1, p_i = sigma[i_prev], sigma[i]
p_j, p_jn = sigma[j], sigma[j_next]
delta_u = (cost_mat[p_im1, p_j] + cost_mat[p_i, p_jn]) - \
(cost_mat[p_im1, p_i] + cost_mat[p_j, p_jn])
# Acceptance
if delta_u <= 0 or rng.random() < np.exp(-beta * delta_u):
current_cost += delta_u
if strategy == 'a' or strategy == 'b':
sigma[i], sigma[j] = sigma[j], sigma[i]
elif strategy == 'c':
if j > i:
sigma[i:j+1] = sigma[i:j+1][::-1]
else:
segment = np.concatenate((sigma[i:], sigma[:j+1]))
rev_segment = segment[::-1]
len_first = n - i
sigma[i:] = rev_segment[:len_first]
sigma[:j+1] = rev_segment[len_first:]
if current_cost < best_cost:
best_cost = current_cost
best_sigma = sigma.copy()
costs_hist.append(current_cost)
# Check Limits
if step % CHECK_INTERVAL == 0:
elapsed = time.time() - start_time
if limit_mode == 'iter' and step >= limit_val:
keep_running = False
elif limit_mode == 'time' and elapsed >= limit_val:
keep_running = False
total_time = time.time() - start_time
# Create axes
iter_axis = np.arange(len(costs_hist))
time_axis = np.linspace(0, total_time, len(costs_hist))
return {
"costs": np.array(costs_hist),
"best_sigma": best_sigma,
"best_cost": best_cost,
"time_axis": time_axis,
"iter_axis": iter_axis
}
# ==============================================================================
# 3. UI DEFINITION
# ==============================================================================
app_ui = theme.layout(
ui.sidebar(
ui.h4("Configuration"),
# Compact Row 1: N and Seed
ui.row(
ui.column(6, ui.input_numeric("n_cities", "N (Cities)", value=DEFAULT_N, min=4, max=200)),
ui.column(6, ui.input_numeric("seed", "Seed", value=DEFAULT_SEED))
),
ui.input_select("cost_metric", "Cost Function", choices=list(COST_FUNCTIONS.keys())),
ui.hr(),
ui.h5("Stop Condition"),
# Compact Row 2: Mode and Limit
ui.row(
ui.column(6, ui.input_select("stop_mode", "Stop By", choices={"iter": "Iterations", "time": "Time (s)"})),
ui.column(6, ui.input_numeric("stop_val", "Limit", value=DEFAULT_ITER))
),
ui.hr(),
ui.h5("Annealing Schedules"),
ui.p("Start → End values for Beta", class_="text-muted small mb-2"),
# Compact Rows for Betas
ui.row(
ui.column(6, ui.input_numeric("beta1_s", "Set 1 Start", value=1.0, step=0.5)),
ui.column(6, ui.input_numeric("beta1_e", "End", value=10.0, step=1.0)),
),
ui.row(
ui.column(6, ui.input_numeric("beta2_s", "Set 2 Start", value=0.1, step=0.1)),
ui.column(6, ui.input_numeric("beta2_e", "End", value=50.0, step=5.0)),
),
ui.row(
ui.column(6, ui.input_numeric("beta3_s", "Set 3 Start", value=0.01, step=0.01)),
ui.column(6, ui.input_numeric("beta3_e", "End", value=100.0, step=10.0)),
),
ui.input_action_button("run", "Run Simulation", class_="btn-primary w-100 mt-3"),
width=350
),
ui.tags.style("""
.card-header { font-weight: bold; }
.shiny-plot-output { margin-bottom: 20px; }
"""),
ui.card(
ui.card_header("1. Best Tours Found (Per Method)"),
ui.output_plot("plot_tours", height="300px")
),
ui.card(
ui.card_header("2. Convergence: Strategy A (Random Swap)"),
ui.output_plot("plot_conv_a", height="250px")
),
ui.card(
ui.card_header("3. Convergence: Strategy B (Adjacent Swap)"),
ui.output_plot("plot_conv_b", height="250px")
),
ui.card(
ui.card_header("4. Convergence: Strategy C (Sub-segment Reverse)"),
ui.output_plot("plot_conv_c", height="250px")
),
title="TSP Metropolis Strategies"
)
# ==============================================================================
# 4. SERVER LOGIC
# ==============================================================================
def server(input, output, session):
theme.setup_plot_style()
sim_data = reactive.Value(None)
# Dynamic Defaults for Limit based on Mode
@reactive.Effect
@reactive.event(input.stop_mode)
def _update_defaults():
mode = input.stop_mode()
if mode == 'iter':
ui.update_numeric("stop_val", value=DEFAULT_ITER)
else:
ui.update_numeric("stop_val", value=DEFAULT_TIME)
@reactive.Effect
@reactive.event(input.run)
def run_simulation_handler():
n = input.n_cities()
seed = input.seed()
metric_name = input.cost_metric()
metric_func = COST_FUNCTIONS[metric_name]
stop_mode = input.stop_mode()
stop_val = input.stop_val()
betas = [
(input.beta1_s(), input.beta1_e()),
(input.beta2_s(), input.beta2_e()),
(input.beta3_s(), input.beta3_e())
]
cities = generate_cities(n, seed)
cost_mat = compute_cost_matrix(cities, metric_func)
strategies = ['a', 'b', 'c']
results_storage = {}
with ui.Progress(min=0, max=9) as p:
count = 0
for strat in strategies:
strat_results = []
for b_range in betas:
p.set(count, message=f"Strat {strat.upper()} ({b_range})...")
res = run_metropolis(cities, cost_mat, strat, b_range, stop_mode, stop_val)
res['beta_range'] = b_range
strat_results.append(res)
count += 1
results_storage[strat] = strat_results
sim_data.set({
"cities": cities,
"results": results_storage,
"metric": metric_name,
"stop_mode": stop_mode
})
def plot_convergence(strat_key, ax):
data = sim_data.get()
if not data: return
res_list = data['results'][strat_key]
mode = data['stop_mode']
for res in res_list:
y = res['costs']
x = res['time_axis'] if mode == 'time' else res['iter_axis']
# Downsampling Logic
if len(y) > MAX_PLOT_POINTS:
# Calculate stride to get roughly MAX_PLOT_POINTS
stride = int(np.ceil(len(y) / MAX_PLOT_POINTS))
y = y[::stride]
x = x[::stride]
beta_txt = f"β: {res['beta_range'][0]}→{res['beta_range'][1]}"
ax.plot(x, y, label=beta_txt, linewidth=1.2, alpha=0.8)
ax.set_ylabel("Total Cost U(σ)")
ax.set_xlabel("Time (s)" if mode == 'time' else "Iterations")
ax.grid(True, linestyle=':', alpha=0.6)
ax.legend(fontsize='small', loc='upper right')
@output
@render.plot
def plot_tours():
data = sim_data.get()
if not data: return
cities = data['cities']
results = data['results']
fig, axes = plt.subplots(1, 3, figsize=(10, 3.5), constrained_layout=True)
titles = {'a': "A: Swap", 'b': "B: Adj Swap", 'c': "C: 2-Opt"}
for idx, strat in enumerate(['a', 'b', 'c']):
ax = axes[idx]
res_list = results[strat]
best_run = min(res_list, key=lambda x: x['best_cost'])
sigma = best_run['best_sigma']
tour = cities[np.concatenate((sigma, [sigma[0]]))]
ax.plot(tour[:, 0], tour[:, 1], 'o-', markersize=3, linewidth=1, color='#2c3e50')
ax.plot(tour[0, 0], tour[0, 1], 'o', color='#e74c3c', markersize=5)
ax.set_title(f"{titles[strat]}\nCost: {best_run['best_cost']:.3f}", fontsize=10)
ax.axis('off')
ax.set_aspect('equal')
return fig
@output
@render.plot
def plot_conv_a():
fig, ax = plt.subplots(figsize=(8, 3), constrained_layout=True)
plot_convergence('a', ax)
return fig
@output
@render.plot
def plot_conv_b():
fig, ax = plt.subplots(figsize=(8, 3), constrained_layout=True)
plot_convergence('b', ax)
return fig
@output
@render.plot
def plot_conv_c():
fig, ax = plt.subplots(figsize=(8, 3), constrained_layout=True)
plot_convergence('c', ax)
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"
})
В конечном счёте, метод сводится к следующему подходу, называемому MCMC (Markov Chain Monte Carlo — методы Монте-Карло для марковских цепей). При заданной \(m \in \mathcal{P}(S)\), найти некоторые «вычислительно оптимальные» вероятности перехода \((p^t_{x,y})_{t\ge 0, x,y\in S}\), такие что \(\sum_{y} m_y p^{t}_{y,x} = m_x\), для всех \(t\in \mathbb{N}, x\in S\). Мы сосредоточимся на однородном по времени случае, когда \(p^t_{x,y}\equiv p_{x,y}\) не зависит от \(t\).
Приведем несколько примеров, предполагая (без потери общности), что \(m_x>0\) для всех \(x\in S\).
- \(p_{x,y}=m_y\), тривиальное точное моделирование, дающее цепь Маркова, которая является просто последовательностью н.о.р. величин. В этом случае распределение \(X_t\) сходится к \(m\) всего за один шаг.
- \(p_{x,y}=\varepsilon \sqrt{m_y/m_x} \ind{x\neq y} + (1- \varepsilon \sum_{z\neq x} \sqrt{m_z/m_x}) \ind{x=y}\), известный как глобальная тепловая баня (global heat bath). Здесь \(0< \varepsilon \le \inf_{x,y} \sqrt{m_x/m_y}\), так что \(p_{x,y}\) действительно является вероятностью перехода с обратимой, и, следовательно, инвариантной мерой \(m\).
- \(p_{x,y}=A_{x,y} \sqrt{m_y/m_x}\), обобщая предыдущие примеры. Здесь \(A_{x,y}\) — достаточно малая симметричная матрица, в том смысле, что \[ \begin{aligned} A_{x,y}=A_{y,x} \ge 0, \qquad x,y \in S, \qquad A_{x,x} \sqrt{m_x}=1-\sum_{z\neq x} S_{x,z} \sqrt{m_z}, \qquad x\in S \end{aligned} \tag{3}\] Ясно, что \(m_x p_{x,y}=A_{x,y} \sqrt{m_x m_y}=m_y p_{y,x}\), так что \(m\) является обратимой, и, следовательно, инвариантной.
В типичных приложениях \(S\) велико, а \(m_x\) — какая-то сложная мера. Первые два примера выше крайне непрактичны. В первом случае мы просто моделируем \(m\) детерминированно. Однако мгновенная сходимость подразумевает очень высокую вычислительную стоимость единственного необходимого шага. В глобальной тепловой бане на каждом шаге нам нужно знать отношения \(\sqrt{m_y/m_x}\) для всех \(x,y\in S\), и несмотря на то, что это, как правило, вычислить проще, чем \(m_x\) (см. Пример 2), это всё же требует доступа к \(|S|^2\) отношениям значений инвариантной меры. Третий общий случай более интересен (и будет дополнительно обобщен ниже), но пока мало что говорит нам о хорошем выборе \(A_{x,y}\).
Динамика Метрополиса — Гастингса
Задача становится намного интереснее как математически, так и практически, если предположить, что задана некоторая дополнительная структура. В частности, предположим, что даны некоторые исходные вероятности перехода \(q_{x,y}\). Осторожно: мы не предполагаем здесь, что \(m\) инвариантна для \(q_{x,y}\); скорее, \(q_{x,y}\) следует воспринимать как вытекающие из некоторой особенности, например геометрических свойств, рассматриваемой задачи. Если \(S\) имеет структуру графа, например, \(q_{x,y}\) могут быть вероятностями перехода простого случайного блуждания на \(S\).
Утверждение 1 (Динамика Метрополиса — Гастингса) Пусть \(S\) конечно, и предположим, что даны:
- Целевая мера \(m\in \mathcal{P}(S)\).
- Ядро предложения (proposal kernel) \(Q=(q_{x,y})_{x,y}\), которое является марковским переходным ядром.
- Кососимметричное вихревое ядро \(R=(r_{x,y})_{x,y}\), такое что \(r_{x,y}=-r_{y,x}\) и \(\sum_y r_{x,y}=0\).
Предположим, что \(r_{x,y} \ge - m_y q_{y,x}\) для всех \(x,y\in S\) (в частности, так как \(R\) кососимметрично, \(- m_y q_{y,x}\le r_{x,y}\le m_x q_{x,y}\)). Тогда \[ p_{x,y}:= \begin{cases} \min(q_{x,y},\frac{r_{x,y}+m_y q_{y,x}}{m_x}) & \text{если $y\neq x$} \\ 1- \sum_{z\neq x} \min(q_{x,z},\frac{r_{x,z}+m_z q_{z,x}}{m_x}) & \text{если $y=x$} \end{cases} \tag{4}\] определяет допустимые марковские вероятности перехода, называемые (обобщённым) ядром Метрополиса — Гастингса, которое допускает \(m\) в качестве инвариантной вероятности.
Предположим более того, что \(q_{x,y}\) неприводимо и что \(q_{x,y}>0\), когда \(q_{y,x}>0\). Тогда \(p_{x,y}\) неприводимо, и, следовательно, \(m\) является единственной инвариантной вероятностью.
Доказательство. Свойство нормировки \(\sum_y p_{x,y}=1\) следует непосредственно из определения \(p_{x,x}\). Далее проверим, что \(p_{x,y}\ge 0\). Это как раз условие \(r_{x,y} \ge - m_y q_{y,x}\) для внедиагональных членов. В то же время \(p_{x,x}\ge 0\), так как мы можем ограничить сумму с отрицательным знаком величиной \(\sum_{z\neq x} q_{x,z}=1-q_{x,x}\le 1\).
Далее проверим, что \(m\) инвариантна. Действительно, \[ \begin{aligned} \sum_{y\neq x} m_x p_{x,y} & = \sum_{y\neq x} \min(m_x q_{x,y},r_{x,y}+m_y q_{y,x}) = \sum_{y\neq x} \min(m_x q_{x,y}- r_{x,y},m_y q_{y,x}) + \sum_{y\neq x} r_{x,y} \\ & = \sum_{y\neq x} \min(m_x q_{x,y}+r_{y,x}, m_y q_{y,x}) =\sum_{y\neq x} m_y p_{y,x} \end{aligned} \] где при переходе от первой ко второй строке мы использовали \(r_{x,y}=-r_{y,x}\) и \(\sum_{y\neq x} r_{x,y}=0\) (так как \(r_{x,x}=0\)).
Алгоритмы Метрополиса — Гастингса работают очень хорошо во многих ситуациях, однако они имеют некоторые слабые стороны:
- Если наша задача состоит в том, чтобы сойтись к инвариантной мере как можно быстрее, наличие шанса \(p_{x,x}>0\) не двигаться вообще вряд ли оптимально.
- При увеличении \(\beta\) наличие локальных минимумов \(U\) может резко увеличить время сходимости цепи Маркова.
Упражнение 1 Проверьте, что если \(r=0\), то \(m\) обратима для вероятностей перехода \(p_{x,y}\), определённых в Уравнение 4.
Упражнение 2 Проверьте, что динамика Метрополиса — Гастингса включает в себя вышеупомянутые примеры: тривиальное точное моделирование, глобальную тепловую баню и симметричную динамику Уравнение 3. (Покажите, что для конкретного выбора \(r_{x,y}\), \(q_{x,y}\) получаются эти динамики).
Динамики Глаубера и Кавасаки
Динамика Глаубера — это частный случай MCMC, часто используемый в статистической физике для спиновых систем (например, модели Изинга), соответствующий «одноузловой тепловой бане». В частности, рассмотрим случай, когда \(S\) — это пространство-произведение, скажем \(S = \{-1, +1\}^V\), для некоторого конечного графа \(V\). Мы интерпретируем \(S\) как спиновое пространство: каждая точка в \(V\) имеет признак, скажем спин, который в данный момент времени может быть либо \(+1\), либо \(-1\). Более того, \(V\) сам по себе является графом, а именно, у нас есть понятие «соседей» данной вершины \(v\in V\). Целевая мера — это мера Гиббса \(m(\sigma) \propto e^{-\beta H(\sigma)}\), где Гамильтониан — это функция \(H\colon S\to \mathbb{R}\), например \(H(\sigma) = - \sum_{(u,v)\in E} J_{u,v} \sigma_u \sigma_v\). На каждом шаге динамика происходит следующим образом:
- Выбрать вершину \(v \in V\) равномерно случайно.
- Заменить спин \(\sigma_v\) новым значением \(\sigma'_v \in \{-1, +1\}\), выбранным из условного распределения спина в \(v\) при условии конфигурации его соседей.
Вероятность перехода для переворота спина в узле \(v\) (изменения \(\sigma\) на \(\sigma^v\)) задаётся формулой: \[ p_{\sigma,\sigma^v} = \frac{1}{|V|} \frac{e^{-\beta H(\sigma^v)}}{e^{-\beta H(\sigma)} + e^{-\beta H(\sigma^v)}} = \frac{1}{|V|} \frac{1}{1 + e^{\beta (H(\sigma^v) - H(\sigma))}} \] Эта динамика обратима относительно меры Гиббса. Ключевая особенность состоит в том, что разность энергий \(H(\sigma^v) - H(\sigma)\) зависит только от соседей \(v\), что делает обновление очень дешёвым вычислительно (\(O(1)\), если степень графа ограничена).
В то время как динамика Глаубера позволяет изменяться полной «магниченности» (сумме спинов), динамика Кавасаки предназначена для получения выборки из меры Гиббса при условии фиксированного числа спинов \(+1\) (канонический ансамбль). Например, мы можем думать в этом случае, что \(S=\{0,1\}^V\), где \(0\) представляет пустое место, а \(1\) — место, занятое частицей. Нам нужна динамика, сохраняющая общее число частиц. Например:
- Выбрать пару вершин \((u,v)\), соединённых ребром в графе.
- Предложить поменять местами значения \(\sigma\) в \(u\) и \(v\).
- Принять обмен с вероятностью, определённой правилом Метрополиса: \[ p_{\sigma, \sigma^{u,v}} = \min\left(1, e^{-\beta(H(\sigma^{u,v}) - H(\sigma))}\right) \] где \(\sigma^{u,v}\) — это конфигурация \(\sigma\), в которой значения в \(u\) и \(v\) поменяны местами.
Если \(\sigma_u = \sigma_v\), обмен ничего не меняет, и энергия постоянна. Если они различаются, частица перемещается. Поскольку частицы только перемещаются и никогда не создаются и не уничтожаются, общее количество спинов \(+1\) сохраняется на протяжении всей эволюции цепи.
Точное моделирование
Наконец, мы обсудим несколько иную ситуацию, где мы хотим получить выборку из неизвестной инвариантной меры данной цепи Маркова, так что вероятности перехода в этом случае даны. Мы предполагаем, что вероятности перехода удовлетворяют следующему условию \[ Z:= \sum_y \inf_{x\in S} p_{x,y} >0 \tag{5}\] Другими словами, существует хотя бы одна точка, в которую цепь Маркова может перепрыгнуть с равномерно положительной вероятностью независимо от начальной точки \(x\). Хотя это может показаться очень ограничительным, даже если \(P=(p_{x,y})_{x,y}\) не удовлетворяет Уравнение 5, возможно, некоторая степень \(P^t\) будет удовлетворять. И в этом случае мы можем просто применить ту же технику к вероятности перехода за \(t\) шагов. В частности, это сработает для любой апериодической цепи на конечном или счётном пространстве состояний.
Примечание. Предполагая Уравнение 5, мы можем записать \[ p_{x,y}= (1-\varepsilon) q_{x,y}+ \varepsilon \nu_y \tag{6}\] для некоторых \(\varepsilon>0\), \(\nu\in \mathcal{P}(S)\) и марковских вероятностей перехода \(Q=(q_{x,y})_{x,y}\). Например, мы можем взять \(\varepsilon=Z\) и \[ q_{x,y}= (p_{x,y}-\varepsilon \nu_y)/(1-\varepsilon) \]
Заметим, что если \(p_{x,y}\) неприводима (или, более общо, рассматривая только сообщающийся класс, содержащий носитель \(\nu\)), то \(P\) индуцирует положительно возвратную цепь Маркова (так как на каждом шаге мы имеем вероятность \(\varepsilon \nu_y>0\) вернуться в некоторую точку \(y\) в носителе \(\nu\)). В частности, цепь Маркова допускает единственное инвариантное распределение \(m\in \mathcal{P}(S)\).
Утверждение 2 (Алгоритм точного моделирования) Предположим, неприводимая цепь Маркова на конечном или счётном пространстве состояний \(S\), с вероятностями перехода \(P=(p_{x,y})\), удовлетворяет Уравнение 5, и пусть \(m\) — её инвариантная вероятность. Для любого разложения \(p_{x,y}= (1-\varepsilon) q_{x,y}+ \varepsilon \nu_y\) как в Уравнение 6 выполняется следующее.
Пусть \(\tau\) — геометрическая случайная величина с параметром \(\varepsilon\), то есть \(\mathbb{P}(\tau=k)=\varepsilon (1-\varepsilon)^k\) для \(k\in \mathbb{N}\). Пусть \(\mathbf{Y}:=(Y_t)\) — цепь Маркова с вероятностями перехода \(Q=(q_{x,y})\) и начальным распределением \(\nu\). Тогда \(Y_\tau\) имеет распределение \(m\).
Доказательство. Пусть \(\mu\) — распределение \(Y_\tau\), тогда для \(y\in S\) \[ \begin{aligned} \mu_y = \sum\nolimits_{k=0}^\infty \mathbb{P}(Y_\tau=y| \tau=k) \mathbb{P}(\tau=k)= \sum\nolimits_{k=0}^\infty \varepsilon (1-\varepsilon)^k \mathbb{P}(Y_k=y) \end{aligned} \] Следовательно, \(\mu =\varepsilon \sum_k (1-\varepsilon)^k \nu Q^k\). Нам нужно проверить \(\mu=\mu P\), так как мы уже знаем, что инвариантная мера единственна, и единственным решением этого уравнения является \(\mu=m\).
В операторной записи \(P= \varepsilon R+(1-\varepsilon) Q\), где \(R_{x,y}=\nu_y\) для всех \(y\in S\). Заметим, что \(\lambda R= \nu\) для всех \(\lambda \in \mathcal{P}(S)\). Следовательно, \[ \begin{aligned} \mu P & = \varepsilon \mu R + (1-\varepsilon) \mu Q= \varepsilon \nu + (1-\varepsilon) \left(\varepsilon \sum\nolimits_{k=0}^\infty (1-\varepsilon)^k \nu Q^k \right) Q= \varepsilon \left(\nu +\sum\nolimits_{k=0}^\infty (1-\varepsilon)^{k+1} \nu Q^{k+1} \right) \\ & = \varepsilon \left( \nu (1-\varepsilon)^{0} \nu Q^{0}+ \sum\nolimits_{k=1}^\infty (1-\varepsilon)^{k} \nu Q^{k} \right)=\mu \end{aligned} \]
Это означает, что если \(\nu\) и \(q_{x,y}\) доступны для вычислений, у нас есть очень практичный способ получить выборку из \(m\):
- Смоделируйте геометрическую случайную величину \(\tau\) с параметром успеха \(\varepsilon\) (например, подсчитайте количество неудач до первого успеха при подбрасывании монеты с \(\varepsilon\)).
- Смоделируйте точку \(x_0 \in S\) с распределением \(\nu\).
- Смоделируйте \(\tau\) шагов цепи Маркова, начинающейся в \(x_0\), с вероятностями перехода \((q_{x,y})\).
- Положение цепи Маркова после \(\tau\) шагов является выборкой из \(m\).
Этот подход является точным (perfect), поскольку мы не утверждаем никакой сходимости в долгосрочной асимптотике. Это просто точное сэмплирование, хотя количество итераций, необходимых для моделирования \((Y_t)\), случайно и (с малой вероятностью) может быть очень большим.