MCMC представляет собой семейство алгоритмов, которые позволяют генерировать выборки из сложных вероятностных распределений, которые не имеют аналитического решения, либо оно сложно реализуемо. В основе этих методов лежит элегантное сочетание теории Марковских цепей и принципов Монте-Карло, что позволяет исследовать высокоразмерные пространства параметров с удивительной эффективностью.
В этой статье мы погрузимся в мир методов MCMC, теорию и практику Марковских цепей. Осторожно: статья очень длинная, но прочитать ее стоит. Независимо от того, занимаетесь ли вы финансовым моделированием, байесовским выводом или стохастическими процессами, понимание MCMC значительно расширит ваши аналитические возможности.
Теория Марковских процессов: от случайных блужданий к сложным распределениям
Марковский процесс — это стохастический процесс, обладающий «свойством Маркова», то есть, будущие состояния зависят только от текущего состояния и не зависят от истории процесса. Этот принцип выражается знаменитой фразой: «будущее зависит от прошлого только через настоящее».
Формально, для последовательности случайных величин X₁, X₂, …, Xₙ, процесс является Марковским, если:
P(Xₙ₊₁ = x | X₁ = x₁, X₂ = x₂, …, Xₙ = xₙ) = P(Xₙ₊₁ = x | Xₙ = xₙ).
Это свойство позволяет нам моделировать сложные системы с помощью относительно простых вероятностных правил перехода между состояниями.
Марковские цепи, которые представляют собой дискретные Марковские процессы, характеризуются матрицей переходных вероятностей P, где элемент P(i,j) представляет вероятность перехода из состояния i в состояние j.
Ключевое свойство, которое делает Марковские цепи особенно полезными для MCMC, — это их эргодичность. Эргодическая цепь имеет следующие характеристики:
- Неприводимость — из любого состояния можно достичь любого другого состояния за конечное число шагов;
- Апериодичность — цепь не циклична, то есть, нет фиксированной периодичности возвращения в одно и то же состояние;
- Положительная возвратность — всегда есть ожидаемое время возвращения в любое состояние.
При выполнении этих условий Марковская цепь сходится к своему стационарному распределению независимо от начального состояния. Это свойство становится центральным для методов MCMC, где мы конструируем Марковские цепи таким образом, чтобы их стационарное распределение соответствовало интересующему нас распределению вероятностей.
На практике я часто сталкиваюсь с задачами, где непосредственное вычисление вероятностей или ожиданий по сложным распределениям невозможно. Но используя свойства Марковских цепей, мы можем сгенерировать выборки, которые асимптотически следуют нужному распределению, и затем использовать эти выборки для численного оценивания интересующих нас величин.
Принципы методов Монте-Карло и их ограничения
Методы Монте-Карло представляют собой класс вычислительных алгоритмов, которые используют случайность для решения детерминированных задач. Название происходит от района Монте-Карло в Монако, известного своими казино, что метафорически отражает случайную природу этих методов. Основная идея методов Монте-Карло заключается в многократном генерировании случайных выборок для получения численных результатов.
Классические методы Монте-Карло основаны на законе больших чисел, который гарантирует, что среднее значение большого числа независимых случайных величин сходится к их ожидаемому значению. Например, для оценки интеграла функции f(x) по области D мы можем сгенерировать N случайных точек в D и вычислить среднее значение f в этих точках, что даст нам приближение к интегралу.
Однако, несмотря на всю мощь, классические методы Монте-Карло имеют существенные ограничения:
- Проклятие размерности — эффективность методов быстро падает с увеличением размерности пространства. В высокоразмерных пространствах большинство случайных выборок концентрируются в областях с низкой вероятностью, что делает прямое семплирование неэффективным.
- Сложность генерации выборок из произвольных распределений — для многих сложных распределений не существует прямых методов генерации случайных чисел.
- Вычислительная сложность — в задачах с большим числом параметров требуется огромное количество выборок для достижения приемлемой точности.
- Необходимость знать нормализующую константу — многие распределения известны только с точностью до множителя, что делает прямое семплирование невозможным.
Именно эти ограничения и привели к разработке более совершенных методов, таких как MCMC, которые позволяют эффективно исследовать сложные вероятностные пространства. Вместо того чтобы генерировать независимые выборки, MCMC создает коррелированную последовательность точек, которая, тем не менее, сходится к нужному распределению.
В моей практике я неоднократно сталкивался с ситуациями, когда классические методы Монте-Карло оказывались бессильны перед сложностью задачи. Например, в финансовом моделировании при оценке производных инструментов со множеством зависимых факторов или в байесовской статистике при работе с иерархическими моделями. В таких случаях MCMC становится не просто удобным, а необходимым инструментом для получения надежных результатов.
Основные принципы MCMC: как построить правильную цепь
Сердцем MCMC является идея конструирования Марковской цепи таким образом, чтобы ее стационарное распределение совпадало с целевым распределением, из которого мы хотим получить выборки. Это достигается путем тщательного проектирования правил перехода между состояниями.
Ключевым свойством, которое должна иметь наша Марковская цепь, является детальный баланс. Это условие гарантирует, что для любых двух состояний i и j, поток вероятности из i в j равен потоку из j в i:
π(i)P(i,j) = π(j)P(j,i)
где:
- π(i) — вероятность состояния i в стационарном распределении
- P(i,j) — вероятность перехода из i в j.
Если детальный баланс выполняется для всех пар состояний, то π является стационарным распределением Марковской цепи. Важно отметить, что условие детального баланса является достаточным, но не необходимым для существования стационарного распределения.
Для построения Марковской цепи с нужным стационарным распределением мы используем различные стратегии, наиболее распространенными из которых являются:
- Алгоритмы, основанные на принятии/отклонении предложенных переходов (например, алгоритм Метрополиса-Гастингса).
- Алгоритмы, использующие полные условные распределения (например, семплирование по Гиббсу).
- Гибридные алгоритмы, сочетающие различные подходы (например, Hamiltonian Monte Carlo).
Каждый из этих подходов имеет свои преимущества и ограничения, которые необходимо учитывать при выборе конкретного метода для решения задачи.
В процессе построения MCMC алгоритма необходимо обеспечить:
- Сходимость к стационарному распределению — цепь должна иметь единственное стационарное распределение, которое совпадает с целевым.
- Эффективное исследование пространства состояний — алгоритм должен быстро исследовать области высокой вероятности и не застревать в локальных модах.
- Баланс между исследованием (exploration) и эксплуатацией (exploitation) — необходимо найти компромисс между широким охватом пространства и концентрацией выборок в областях высокой вероятности.
На практике я стараюсь внимательно подходить к выбору начальных значений, периода «прогрева» (burn-in), оценке автокорреляции и других диагностических мер, чтобы убедиться, что мой MCMC алгоритм работает корректно и эффективно для конкретной задачи.
Алгоритм Метрополиса-Гастингса: классический подход к MCMC
Алгоритм Метрополиса-Гастингса (МГ) является, пожалуй, наиболее фундаментальным и широко используемым методом MCMC. Он представляет собой обобщение алгоритма Метрополиса, разработанного в 1953 году, и получил свое современное формулирование благодаря работе Гастингса в 1970 году.
Элегантность алгоритма МГ заключается в его способности генерировать выборки из распределения, известного лишь с точностью до нормализующей константы, что делает его незаменимым инструментом для байесовской статистики.
Основная идея алгоритма МГ состоит в построении Марковской цепи путем предложения новых состояний из распределения предложений (proposal distribution) и принятия или отклонения этих предложений с определенной вероятностью, которая гарантирует сходимость к целевому распределению.
Пусть π(x) — целевая функция плотности вероятности, из которой мы хотим получить выборки, а q(x’ | x) — распределение предложений, определяющее вероятность предложить переход из текущего состояния x в новое состояние x’. Алгоритм МГ состоит из следующих шагов:
1) Инициализация: выбрать начальное состояние x₀.
2) Для каждого шага t = 1, 2, …:
- Сгенерировать кандидата x’ из распределения предложений q(x’ | xₜ₋₁).
- Вычислить вероятность принятия: α = min(1, (π(x’) * q(xₜ₋₁ | x’)) / (π(xₜ₋₁) * q(x’ | xₜ₋₁)))
- С вероятностью α принять x’ и установить xₜ = x’, в противном случае оставить xₜ = xₜ₋₁.
Ключевой момент здесь — вероятность принятия α, которая обеспечивает выполнение условия детального баланса:
π(x)q(x’ | x)α(x, x’) = π(x’)q(x | x’)α(x’, x)
где α(x, x’) — вероятность принятия перехода из x в x’.
Элегантность этого решения в том, что для вычисления α нам не нужно знать нормализующую константу целевого распределения, так как она сокращается в отношении π(x’)/π(x). Это делает алгоритм МГ применимым к широкому классу задач, где целевое распределение известно только с точностью до множителя, что часто встречается в байесовской статистике с ненормированными апостериорными распределениями.
Математическое обоснование алгоритма МГ базируется на теории эргодических Марковских цепей. Можно показать, что при достаточно общих условиях на π и q, построенная цепь является эргодической, и ее стационарное распределение совпадает с π. Это означает, что после достаточного количества шагов («прогрева»), состояния цепи можно рассматривать как выборки из целевого распределения.
Реализация на Python: простой пример
Давайте реализуем алгоритм Метрополиса-Гастингса на Python для решения простой задачи — получения выборок из смеси нормальных распределений. Этот пример позволит нам увидеть, как работает алгоритм МГ на практике, и проиллюстрирует его способность исследовать многомодальные распределения.
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
# Определяем целевое распределение - смесь двух нормальных распределений
def target_distribution(x):
# Ненормированная плотность вероятности
return 0.3 * stats.norm.pdf(x, -3, 1) + 0.7 * stats.norm.pdf(x, 2, 2)
# Реализация алгоритма Метрополиса-Гастингса
def metropolis_hastings(n_samples, proposal_width=1.0):
# Инициализация
samples = np.zeros(n_samples)
x_current = 0 # Начальное состояние
# Счетчик принятых предложений
accepted = 0
for i in range(n_samples):
# Генерация предложения
x_proposal = x_current + np.random.normal(0, proposal_width)
# Вычисление отношения плотностей
p_current = target_distribution(x_current)
p_proposal = target_distribution(x_proposal)
# Вычисление вероятности принятия
acceptance_ratio = p_proposal / p_current
# Принятие или отклонение предложения
if np.random.uniform() < acceptance_ratio:
x_current = x_proposal
accepted += 1
samples[i] = x_current
acceptance_rate = accepted / n_samples
return samples, acceptance_rate
# Функция автокорреляции
def autocorrelation(x, max_lag):
result = np.zeros(max_lag)
x_mean = np.mean(x)
x_var = np.var(x)
for lag in range(max_lag):
if lag == 0:
# Для нулевого лага автокорреляция равна 1
result[lag] = 1.0
else:
# Для остальных лагов используем стандартную формулу
result[lag] = np.mean((x[:-lag] - x_mean) * (x[lag:] - x_mean)) / x_var
return result
# Запуск алгоритма
n_samples = 50000
burn_in = 5000 # Количество шагов "прогрева"
samples, acceptance_rate = metropolis_hastings(n_samples)
# Отбрасываем период "прогрева"
samples = samples[burn_in:]
# Визуализация результатов
plt.figure(figsize=(12, 8))
# Гистограмма полученных выборок
plt.subplot(2, 1, 1)
plt.hist(samples, bins=50, density=True, alpha=0.7, color='#333333')
# Истинное распределение для сравнения
x = np.linspace(-8, 8, 1000)
plt.plot(x, target_distribution(x), 'r-', lw=2, label='Истинное распределение')
plt.title('Сравнение выборок MCMC с целевым распределением')
plt.xlabel('x')
plt.ylabel('Плотность вероятности')
plt.legend()
# График трассировки для оценки перемешивания
plt.subplot(2, 1, 2)
plt.plot(samples[:1000], color='#333333')
plt.title('Трассировка MCMC (первые 1000 выборок после прогрева)')
plt.xlabel('Итерация')
plt.ylabel('Значение')
plt.tight_layout()
plt.savefig('metropolis_hastings_example.png', dpi=300)
plt.show()
print(f"Коэффициент принятия: {acceptance_rate:.2f}")
# Автокорреляция выборок
max_lag = 50
acf = autocorrelation(samples, max_lag)
plt.figure(figsize=(10, 6))
plt.bar(range(max_lag), acf, color='#333333')
plt.title('Автокорреляционная функция выборок MCMC')
plt.xlabel('Лаг')
plt.ylabel('Автокорреляция')
plt.axhline(y=0, linestyle='--', color='#999999')
plt.grid(True, alpha=0.3)
plt.savefig('metropolis_hastings_autocorrelation.png', dpi=300)
plt.show()
# Дополнительная диагностика
print(f"Среднее значение выборок: {np.mean(samples):.3f}")
print(f"Стандартное отклонение выборок: {np.std(samples):.3f}")
print(f"Размер выборки после прогрева: {len(samples)}")
# Эффективный размер выборки
effective_sample_size = len(samples) / (1 + 2 * np.sum(acf[1:]))
print(f"Эффективный размер выборки: {effective_sample_size:.0f}")
Рис. 1: Сравнение эмпирического распределения выборок, полученных методом Метрополиса-Гастингса, с целевым распределением (смесь двух нормальных распределений) и трассировка цепи Маркова для первых 1000 итераций после периода прогрева
Рис. 2: Автокорреляционная функция выборок MCMC. Быстрое убывание автокорреляции указывает на низкую зависимость между соседними выборками, что свидетельствует об эффективности алгоритма
Коэффициент принятия: 0.85
Среднее значение выборок: 0.625
Стандартное отклонение выборок: 2.831
Размер выборки после прогрева: 45000
Эффективный размер выборки: 1032Этот код иллюстрирует основные компоненты алгоритма Метрополиса-Гастингса и включает важные элементы анализа результатов: визуализацию выборок в сравнении с целевым распределением, график трассировки для оценки «перемешивания» цепи и автокорреляционную функцию для оценки независимости выборок.
В реализации алгоритма мы использовали нормальное распределение в качестве распределения предложений (q(x’ | x) ~ N(x, σ²)), что соответствует случайному блужданию в пространстве состояний. Параметр proposal_width (σ) контролирует размер шага и существенно влияет на эффективность алгоритма. Слишком маленькое значение приведет к высокой автокорреляции выборок и медленному исследованию пространства, а слишком большое — к низкому коэффициенту принятия и также неэффективному семплированию.
Обратите внимание на период «прогрева» (burn-in), который мы отбрасываем перед анализом результатов. Это стандартная практика в MCMC, поскольку начальные состояния цепи могут сильно зависеть от выбранного начального значения и еще не представлять стационарное распределение.
График трассировки позволяет визуально оценить качество «перемешивания» цепи — хорошо перемешанная цепь должна быстро переходить между различными областями пространства состояний и не «застревать» в отдельных регионах. Автокорреляционная функция дает количественную оценку зависимости между последовательными выборками — быстрое убывание автокорреляции с увеличением лага указывает на более эффективное семплирование.
Коэффициент принятия (acceptance rate) — еще один важный диагностический показатель. Оптимальное значение зависит от размерности задачи, но часто рекомендуют поддерживать его в диапазоне 0.2-0.5 для случайного блуждания. Слишком низкий коэффициент означает, что цепь редко перемещается и неэффективно исследует пространство, а слишком высокий может указывать на слишком мелкие шаги и высокую автокорреляцию.
Тонкости и оптимизация работы алгоритма
Успешное применение алгоритма Метрополиса-Гастингса требует понимания ряда тонкостей и возможных оптимизаций. За годы работы с MCMC я выработал несколько практических рекомендаций, которые значительно повышают эффективность этого метода.
Настройка распределения предложений
Выбор подходящего распределения предложений q(x’ | x) важен для эффективности МГ. Идеальное распределение предложений должно быть близко к целевому распределению, но при этом легко семплироваться. Вот несколько стратегий:
- Адаптивный МГ — автоматическая настройка параметров распределения предложений на основе истории цепи. Например, можно периодически обновлять ковариационную матрицу для многомерного нормального распределения предложений, используя накопленные выборки.
- Использование распределения предложений, учитывающего структуру целевого распределения. Например, для распределений с тяжелыми хвостами лучше использовать распределение предложений с тяжелыми хвостами, а не нормальное.
- Настройка масштаба предложений для достижения оптимального коэффициента принятия. В академическом сообществе для многомерных распределений оптимальный коэффициент принятия принято считать за 0.234.
def adaptive_metropolis_hastings(n_samples, target_log_density, initial_state, adaptation_interval=100):
d = len(initial_state) # Размерность
samples = np.zeros((n_samples, d))
current_state = initial_state.copy()
# Начальная ковариационная матрица для предложений
cov = np.eye(d) * 0.1
# Счетчик принятий
accepted = 0
for i in range(n_samples):
# Адаптация ковариационной матрицы
if i > 0 and i % adaptation_interval == 0 and i <= n_samples // 2: # Используем накопленные выборки для оценки ковариации if i > adaptation_interval:
cov = np.cov(samples[:i].T) + 1e-6 * np.eye(d)
# Масштабирование для достижения оптимального коэффициента принятия
current_acceptance_rate = accepted / max(1, i)
scale = 1
if current_acceptance_rate < 0.2: scale = 0.9 # Уменьшаем шаг, если принимается слишком мало предложений elif current_acceptance_rate > 0.3:
scale = 1.1 # Увеличиваем шаг, если принимается слишком много предложений
cov *= scale
# Генерация предложения
proposal = np.random.multivariate_normal(current_state, cov)
# Вычисление логарифмов плотностей
log_p_current = target_log_density(current_state)
log_p_proposal = target_log_density(proposal)
# Вычисление логарифма вероятности принятия
log_acceptance_ratio = log_p_proposal - log_p_current
# Принятие или отклонение предложения
if np.log(np.random.uniform()) < log_acceptance_ratio:
current_state = proposal
accepted += 1
samples[i] = current_state
acceptance_rate = accepted / n_samples
return samples, acceptance_rateРабота с многомодальными распределениями
Стандартный МГ может неэффективно исследовать многомодальные распределения, поскольку цепь может «застрять» в одной из мод на долгое время. Для решения этой проблемы применяются следующие подходы:
- Параллельные цепи с обменом (Parallel Tempering или Replica Exchange MCMC) — запускается несколько цепей с различными «температурами», что позволяет более горячим цепям легче преодолевать барьеры между модами, и периодически обмениваются состояниями между цепями;
- Использование смесей распределений предложений, включающих как локальные шаги (для исследования текущей моды), так и глобальные прыжки (для перехода между модами);
- Многозапускные стратегии — запуск нескольких независимых цепей из различных начальных точек и объединение результатов.
Вычислительные оптимизации
Ниже пример кода демонстрирующий работу с логарифмами вероятностей вместо самих вероятностей для предотвращения численной нестабильности:
def metropolis_hastings_log_scale(n_samples, target_log_density, proposal_width=1.0):
samples = np.zeros(n_samples)
x_current = 0
for i in range(n_samples):
x_proposal = x_current + np.random.normal(0, proposal_width)
log_p_current = target_log_density(x_current)
log_p_proposal = target_log_density(x_proposal)
# Логарифм отношения вероятностей
log_acceptance_ratio = log_p_proposal - log_p_current
# Принятие или отклонение в логарифмической шкале
if np.log(np.random.uniform()) < log_acceptance_ratio:
x_current = x_proposal
samples[i] = x_current
return samplesТакже мы можем применять другие методы вычислительной оптимизации:
- Использование градиентной информации (если доступна) для более эффективного исследования пространства. Это приводит нас к более продвинутым методам, таким как Hamiltonian Monte Carlo, которые мы обсудим далее;
- Векторизация вычислений для одновременного обновления нескольких параметров, что особенно эффективно при использовании библиотек, оптимизированных для векторных операций, таких как NumPy.
Диагностика и мониторинг
Для оценки качества работы MCMC алгоритма необходимо анализировать ряд диагностических метрик:
Эффективный размер выборки (Effective Sample Size, ESS) — оценка количества «независимых» выборок, эквивалентных полученной автокоррелированной последовательности:
def effective_sample_size(samples):
n = len(samples)
max_lag = min(n // 5, 1000) # Ограничиваем максимальный лаг
# Вычисляем автокорреляцию
acf = autocorrelation(samples, max_lag)
# Находим лаг, при котором автокорреляция становится отрицательной или близкой к нулю
cutoff = max_lag
for i in range(1, max_lag):
if acf[i] < 0.05: # Порог близкий к нулю
cutoff = i
break
# Вычисляем сумму автокорреляций до точки отсечения
sum_rho = 1 + 2 * np.sum(acf[1:cutoff])
# Эффективный размер выборки
ess = n / sum_rho
return essМетрика Гельмана-Рубина (Gelman-Rubin statistic или potential scale reduction factor, PSRF) — для оценки сходимости нескольких цепей. Значения близкие к 1 указывают на хорошую сходимость:
def gelman_rubin(chains):
# chains - список массивов, каждый представляет отдельную цепь MCMC
m = len(chains) # Количество цепей
n = chains[0].shape[0] # Длина каждой цепи
# Средние значения по каждой цепи
chain_means = np.array([np.mean(chain) for chain in chains])
# Общее среднее
overall_mean = np.mean(chain_means)
# Межцепная дисперсия
B = n * np.sum((chain_means - overall_mean)**2) / (m - 1)
# Внутрицепные дисперсии
W = np.mean([np.var(chain, ddof=1) for chain in chains])
# Оценка дисперсии
var_hat = (n - 1) / n * W + B / n
# Фактор потенциального сокращения шкалы
R_hat = np.sqrt(var_hat / W)
return R_hatМы также можем использовать графический анализ трассировок, автокорреляций и распределений выборок для визуального контроля качества.
Преимущества и недостатки
Алгоритм Метрополиса-Гастингса имеет ряд неоспоримых преимуществ:
- Простота реализации и интуитивная понятность;
- Гибкость в выборе целевого распределения и распределения предложений;
- Теоретическая обоснованность и асимптотическая корректность.
Однако, у него есть и существенные недостатки:
- Низкая эффективность в высокоразмерных пространствах («проклятие размерности»);
- Трудности с исследованием многомодальных распределений;
- Необходимость тонкой настройки параметров для оптимальной производительности;
- Высокая автокорреляция выборок, что снижает эффективный размер выборки.
Эти ограничения привели к разработке более совершенных методов MCMC, которые мы рассмотрим в следующих разделах. Тем не менее, понимание принципов работы алгоритма МГ и умение его правильно настраивать — важный навык для каждого специалиста по анализу данных и статистическому моделированию.
Семплирование по Гиббсу
Семплирование по Гиббсу (Gibbs sampling) представляет собой специальный случай алгоритма Метрополиса-Гастингса, который особенно эффективен для многомерных распределений, где условные распределения каждой переменной при фиксированных остальных переменных легко семплировать. Этот метод назван в честь физика Джозайи Уилларда Гиббса и является одним из наиболее широко используемых алгоритмов MCMC в байесовской статистике и машинном обучении.
В отличие от алгоритма Метрополиса-Гастингса, который обновляет все параметры одновременно или группами, семплирование по Гиббсу последовательно обновляет каждую переменную по отдельности, используя ее полное условное распределение при фиксированных значениях всех остальных переменных. Ключевое преимущество этого подхода заключается в том, что вероятность принятия каждого предложения всегда равна 1, т.е. каждое предложение принимается, что значительно повышает эффективность алгоритма.
Формально, для многомерного распределения p(x₁, x₂, …, xₙ), алгоритм Гиббса последовательно генерирует выборки из условных распределений p(xᵢ | x₁, …, xᵢ₋₁, xᵢ₊₁, …, xₙ) для каждого i от 1 до n. Одна итерация алгоритма состоит из n шагов, на каждом из которых обновляется одна переменная.
Можно показать, что при определенных условиях, в частности, если все условные распределения положительны, полученная Марковская цепь сходится к совместному распределению p(x₁, x₂, …, xₙ). Это следует из того, что каждое условное обновление сохраняет целевое распределение как инвариантное.
Основные отличия семплирования по Гиббсу от алгоритма Метрополиса-Гастингса:
- В Гиббсе обновляются переменные по одной, а не все сразу;
- Каждое предложение принимается с вероятностью 1, нет необходимости в вычислении вероятности принятия;
- Для Гиббса необходимо знать и уметь семплировать из полных условных распределений, в то время как МГ требует только возможности вычислять отношение плотностей вероятности.
Семплирование по Гиббсу особенно эффективно для определенных классов моделей, таких как:
- Иерархические байесовские модели с сопряженными априорными распределениями;
- Графические модели с определенной структурой зависимостей;
- Модели скрытых переменных, такие как смеси распределений или скрытые марковские модели.
В этих случаях условные распределения часто имеют стандартную форму, из которой легко генерировать выборки, что делает алгоритм Гиббса идеальным выбором.
Реализация на Python: байесовская линейная регрессия
Байесовская линейная регрессия представляет собой отличный пример модели, для которой семплирование по Гиббсу особенно эффективно. В этой модели мы хотим оценить параметры линейной регрессии (коэффициенты β и дисперсию шума σ²) с учетом априорных распределений на эти параметры.
Рассмотрим следующую байесовскую модель:
- Данные: y = Xβ + ε, где ε ~ N(0, σ² · I)
- Априорное распределение для β: β ~ N(μ₀, Λ₀⁻¹)
- Априорное распределение для σ²: σ² ~ Inv-Gamma(a₀, b₀)
Для этой модели полные условные распределения имеют следующий вид:
p(β | σ², X, y) = N(μₙ, Λₙ⁻¹)
где:
- Λₙ = X^T X / σ² + Λ₀
- μₙ = Λₙ⁻¹ (X^T y / σ² + Λ₀ μ₀)
p(σ² | β, X, y) = Inv-Gamma(aₙ, bₙ)
где:
- aₙ = a₀ + n/2
- bₙ = b₀ + (y — Xβ)^T (y — Xβ) / 2
Теперь реализуем семплирование по Гиббсу для этой модели:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
import seaborn as sns
from sklearn.preprocessing import StandardScaler
# Загрузка и подготовка данных
def load_stock_data(ticker='MSFT', start_date='2022-01-01', end_date='2023-01-01'):
import yfinance as yf
data = yf.download(ticker, start=start_date, end=end_date)
return data
# Загружаем данные Microsoft
msft_data = load_stock_data('MSFT')
# Подготовка данных для регрессии
def prepare_features(data, lags=5):
# Используем логарифмические доходности
data['log_return'] = np.log(data['Close'] / data['Close'].shift(1))
# Создаем признаки на основе лагированных доходностей
for lag in range(1, lags + 1):
data[f'lag_{lag}'] = data['log_return'].shift(lag)
# Удаляем строки с NaN
data = data.dropna()
# Стандартизация признаков
scaler = StandardScaler()
features = [f'lag_{lag}' for lag in range(1, lags + 1)]
X = scaler.fit_transform(data[features])
y = np.array(data['log_return'])
return X, y, features
X, y, feature_names = prepare_features(msft_data)
# Реализация Гиббса для байесовской линейной регрессии
def gibbs_bayesian_regression(X, y, n_samples=5000, burn_in=1000,
mu_0=None, Lambda_0=None, a_0=0.01, b_0=0.01):
n, p = X.shape
# Инициализация априорных параметров
if mu_0 is None:
mu_0 = np.zeros(p)
if Lambda_0 is None:
Lambda_0 = np.eye(p) * 0.01
# Инициализация цепи
beta = np.zeros(p)
sigma_sq = 1.0
# Хранение выборок
beta_samples = np.zeros((n_samples, p))
sigma_sq_samples = np.zeros(n_samples)
# Предварительные вычисления
X_T_X = X.T @ X
X_T_y = X.T @ y
for i in range(n_samples + burn_in):
# Обновление beta
Lambda_n = X_T_X / sigma_sq + Lambda_0
Lambda_n_inv = np.linalg.inv(Lambda_n)
mu_n = Lambda_n_inv @ (X_T_y / sigma_sq + Lambda_0 @ mu_0)
beta = np.random.multivariate_normal(mu_n, Lambda_n_inv)
# Обновление sigma^2
a_n = a_0 + n / 2
b_n = b_0 + np.sum((y - X @ beta)**2) / 2
sigma_sq = 1 / np.random.gamma(a_n, 1 / b_n)
# Сохранение выборок после прогрева
if i >= burn_in:
beta_samples[i - burn_in] = beta
sigma_sq_samples[i - burn_in] = sigma_sq
return beta_samples, sigma_sq_samples
# Запуск алгоритма Гиббса
n_samples = 10000
burn_in = 2000
beta_samples, sigma_sq_samples = gibbs_bayesian_regression(X, y, n_samples, burn_in)
# Анализ результатов
def plot_gibbs_results(beta_samples, sigma_sq_samples, feature_names):
# Оценка параметров модели (апостериорные средние)
beta_mean = np.mean(beta_samples, axis=0)
sigma_sq_mean = np.mean(sigma_sq_samples)
# Апостериорные интервалы для beta
beta_5 = np.percentile(beta_samples, 5, axis=0)
beta_95 = np.percentile(beta_samples, 95, axis=0)
# Таблица с результатами
results = pd.DataFrame({
'Признак': feature_names,
'Среднее': beta_mean,
'5% квантиль': beta_5,
'95% квантиль': beta_95
})
print("Апостериорные оценки параметров:")
print(results)
print(f"Апостериорная средняя дисперсия шума: {sigma_sq_mean:.6f}")
# Визуализация апостериорных распределений коэффициентов
plt.figure(figsize=(14, 8))
# Трассировки для beta
plt.subplot(2, 2, 1)
for j in range(len(feature_names)):
plt.plot(beta_samples[:1000, j], alpha=0.7, label=feature_names[j])
plt.title('Трассировки для коэффициентов (первые 1000 выборок)')
plt.xlabel('Итерация')
plt.ylabel('Значение')
plt.legend()
# Трассировка для sigma^2
plt.subplot(2, 2, 2)
plt.plot(sigma_sq_samples[:1000], color='black')
plt.title('Трассировка для дисперсии шума')
plt.xlabel('Итерация')
plt.ylabel('Значение')
# Апостериорные плотности для beta
plt.subplot(2, 2, 3)
for j in range(len(feature_names)):
sns.kdeplot(beta_samples[:, j], label=feature_names[j])
plt.title('Апостериорные распределения коэффициентов')
plt.xlabel('Значение')
plt.ylabel('Плотность')
plt.legend()
# Апостериорная плотность для sigma^2
plt.subplot(2, 2, 4)
sns.kdeplot(sigma_sq_samples, color='black')
plt.title('Апостериорное распределение дисперсии шума')
plt.xlabel('Значение')
plt.ylabel('Плотность')
plt.tight_layout()
plt.savefig('gibbs_bayesian_regression.png', dpi=300)
plt.show()
# Автокорреляция для одного из коэффициентов
plt.figure(figsize=(10, 6))
max_lag = 50
acf = np.array([np.corrcoef(beta_samples[:-lag, 0], beta_samples[lag:, 0])[0, 1]
for lag in range(1, max_lag + 1)])
plt.bar(range(1, max_lag + 1), acf, color='#333333')
plt.title('Автокорреляция первого коэффициента')
plt.xlabel('Лаг')
plt.ylabel('Автокорреляция')
plt.axhline(y=0, linestyle='--', color='#999999')
plt.savefig('gibbs_bayesian_regression_autocorr.png', dpi=300)
plt.show()
plot_gibbs_results(beta_samples, sigma_sq_samples, feature_names)
# Предсказания с использованием апостериорных выборок
def predict_with_uncertainty(X_new, beta_samples, sigma_sq_samples, n_pred_samples=1000):
# Количество тестовых точек и выборок
n_test = X_new.shape[0]
n_samples = min(len(beta_samples), n_pred_samples)
# Случайно выбираем подмножество выборок для предсказания
idx = np.random.choice(len(beta_samples), n_samples, replace=False)
beta_subset = beta_samples[idx]
sigma_subset = np.sqrt(sigma_sq_samples[idx])
# Предсказания для каждой выборки параметров
predictions = np.zeros((n_samples, n_test))
for i in range(n_samples):
# Среднее предсказание
mean_pred = X_new @ beta_subset[i]
# Добавляем шум
predictions[i] = mean_pred + np.random.normal(0, sigma_subset[i], size=n_test)
# Статистики по предсказаниям
mean_predictions = np.mean(predictions, axis=0)
lower_bound = np.percentile(predictions, 5, axis=0)
upper_bound = np.percentile(predictions, 95, axis=0)
return mean_predictions, lower_bound, upper_bound, predictions
# Пример предсказания
X_test = X[-30:] # Используем последние 30 точек для тестирования
y_test = y[-30:]
mean_pred, lower_bound, upper_bound, all_predictions = predict_with_uncertainty(X_test, beta_samples, sigma_sq_samples)
# Визуализация предсказаний
plt.figure(figsize=(12, 6))
plt.plot(y_test, 'o-', color='black', label='Фактические значения')
plt.plot(mean_pred, '--', color='blue', label='Среднее предсказание')
plt.fill_between(range(len(y_test)), lower_bound, upper_bound, color='blue', alpha=0.2, label='90% доверительный интервал')
# Добавим несколько случайных траекторий для иллюстрации
for i in range(10):
idx = np.random.randint(0, all_predictions.shape[0])
plt.plot(all_predictions[idx], alpha=0.3, color='gray')
plt.title('Байесовские предсказания с неопределенностью')
plt.xlabel('Временной индекс')
plt.ylabel('Логарифмическая доходность')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('gibbs_bayesian_prediction.png', dpi=300)
plt.show()Апостериорные оценки параметров:
Признак Среднее 5% квантиль 95% квантиль
0 lag_1 -0.000872 -0.003457 0.001683
1 lag_2 -0.002052 -0.004600 0.000466
2 lag_3 -0.002222 -0.004772 0.000308
3 lag_4 0.000082 -0.002495 0.002687
4 lag_5 0.001154 -0.001400 0.003717
Апостериорная средняя дисперсия шума: 0.000584
Рис 3: Результаты байесовской линейной регрессии методом Гиббса для прогнозирования доходности акций Microsoft (трассировки коэффициентов регрессии, трассировки дисперсии шума, апостериорные распределения коэффициентов и дисперсии шума)
Рис. 4: Автокорреляционная функция для первого коэффициента регрессии в алгоритме Гиббса
Рис. 5: Байесовские предсказания логарифмической доходности акций Microsoft с квантификацией неопределенности. Черные точки — фактические значения, синяя пунктирная линия — среднее предсказание, сиреневая область — 90% доверительный интервал. Серые траектории иллюстрируют вариабельность предсказаний из апостериорного распределения параметров
Данный код демонстрирует полный рабочий процесс байесовской линейной регрессии с использованием семплирования по Гиббсу:
- Загрузка и подготовка данных — используем реальные данные о ценах акций Microsoft для прогнозирования логарифмических доходностей на основе их прошлых значений;
- Реализация алгоритма Гиббса для байесовской линейной регрессии — поочередное обновление коэффициентов регрессии β и дисперсии шума σ² из их полных условных распределений;
- Анализ результатов — визуализация апостериорных распределений параметров, трассировок цепи и автокорреляций для оценки качества семплирования;
- Байесовские предсказания — использование апостериорных выборок параметров для генерации предсказаний с оценкой неопределенности.
Ключевые преимущества байесовской регрессии с использованием MCMC:
- Полная апостериорная информация о параметрах модели, включая их распределения, корреляции и доверительные интервалы;
- Естественная количественная оценка неопределенности предсказаний, которая учитывает неопределенность в параметрах модели;
- Возможность включения априорной информации, что особенно полезно при ограниченном количестве данных или высокой размерности модели;
- Устойчивость к переобучению благодаря регуляризации через априорные распределения.
В данном примере мы использовали нормально-обратное гамма сопряжение для регрессионных коэффициентов и дисперсии шума, что позволило нам эффективно применить семплирование по Гиббсу. Эта комбинация является классическим примером условно-сопряженной модели, для которой алгоритм Гиббса показывает наибольшую эффективность.
Преимущества и ограничения метода
Семплирование по Гиббсу имеет ряд существенных преимуществ, которые делают его популярным выбором для многих задач байесовского вывода:
Преимущества:
- Высокая эффективность для условно-сопряженных моделей — при наличии аналитических выражений для полных условных распределений, алгоритм Гиббса обеспечивает 100% принятие предложений, что значительно повышает эффективность по сравнению с алгоритмом Метрополиса-Гастингса;
- Естественная интерпретируемость — поскольку каждый параметр обновляется по отдельности, легко отслеживать их эволюцию и понимать, как они взаимодействуют;
- Модульность — для сложных моделей можно использовать различные стратегии обновления для разных групп параметров, комбинируя метод Гиббса с другими методами MCMC;
- Простота реализации — для многих стандартных моделей условные распределения имеют хорошо известную форму, что упрощает реализацию;
- Эффективное исследование условно низкоразмерных пространств — даже если совместное распределение сложное, условные распределения часто проще и легче семплировать.
Ограничения:
- Необходимость знания полных условных распределений — для применения алгоритма Гиббса необходимо иметь аналитические выражения для всех условных распределений или эффективные методы их семплирования, что не всегда возможно;
- Высокая автокорреляция для сильно коррелированных параметров — при наличии сильных корреляций между параметрами, метод Гиббса может очень медленно исследовать пространство состояний, поскольку он обновляет параметры по одному;
- Медленная сходимость для многомодальных распределений — как и многие другие методы MCMC, подход Гиббса может «застревать» в одной моде многомодального распределения;
- Сложность обработки ограничений — если параметры имеют сложные ограничения, условные распределения могут не иметь стандартную форму;
- Необходимость случайного обновления параметров в некоторых случаях — для некоторых моделей порядок обновления параметров может влиять на свойства сходимости, что требует случайного порядка обновления.
Стратегии преодоления ограничений:
- Блочное обновление — вместо обновления каждого параметра по отдельности, можно обновлять группы сильно коррелированных параметров вместе, используя, например, алгоритм Метрополиса-Гастингса для блока;
- Параметризация моделей — изменение параметризации модели может уменьшить корреляции между параметрами и улучшить сходимость Гиббса;
- Использование иерархических моделей — хорошо спроектированная иерархическая структура может упростить условные распределения и улучшить сходимость;
- Комбинирование с другими методами — для параметров, для которых не известны аналитические условные распределения, можно использовать шаги Метрополиса-Гастингса (получая так называемый Metropolis-within-Gibbs);
- Адаптивные схемы обновления — динамическое изменение стратегии обновления параметров на основе накопленной информации о их взаимозависимостях.
В своей практике я обнаружил, что семплирование по Гиббсу часто является методом первого выбора для многих байесовских иерархических моделей, особенно в сочетании с библиотеками, такими как PyMC, которые автоматизируют многие аспекты вывода. Однако для сложных моделей с высокими корреляциями между параметрами или нестандартными распределениями, часто требуются более продвинутые методы, такие как Hamiltonian Monte Carlo, который мы рассмотрим в следующем разделе.
Hamiltonian Monte Carlo: продвинутый метод для сложных моделей
Hamiltonian Monte Carlo (HMC), или гамильтоновский метод Монте-Карло, представляет собой мощное расширение традиционных методов MCMC, которое использует физические принципы для более эффективного исследования пространства параметров.
Этот метод был впервые предложен физиком Саймоном Дуэйном в 1987 году и изначально назывался «гибридным методом Монте-Карло». В контексте байесовского вывода HMC был популяризирован работами Радфорда Нила и сейчас является методом выбора для многих сложных вероятностных моделей.
Физические основы и интуиция метода
HMC основан на физической аналогии движения частицы в потенциальном поле. Представьте, что параметры модели — это положение частицы в многомерном пространстве, а отрицательный логарифм целевой плотности вероятности — это потенциальная энергия этой частицы. Тогда наша задача состоит в том, чтобы исследовать области с низкой потенциальной энергией (высокой плотностью вероятности).
Традиционные методы MCMC, такие как Метрополис-Гастингс с нормальным распределением предложений, соответствуют случайному блужданию частицы, что часто приводит к неэффективному исследованию пространства. HMC, вместо этого, вводит дополнительные переменные «импульса» (momentum) для каждого параметра и использует законы гамильтоновой механики для определения траектории частицы.
Ключевая идея состоит в том, что частица может накапливать кинетическую энергию, которая позволяет ей преодолевать области с высокой потенциальной энергией (низкой плотностью вероятности) и эффективно исследовать пространство параметров даже в высокоразмерных задачах.
Математическая формулировка
Формально, HMC расширяет пространство состояний, добавляя к каждому параметру модели θ соответствующую переменную импульса p. Вместе они образуют расширенное состояние (θ, p). Определяется гамильтониан H(θ, p), который представляет собой общую энергию системы:
H(θ, p) = U(θ) + K(p)
где:
- U(θ) = -log(π(θ)) — потенциальная энергия, равная отрицательному логарифму ненормированной целевой плотности π(θ);
- K(p) = p^T M^(-1) p / 2 — кинетическая энергия, где M — матрица масс (обычно диагональная или пропорциональная единичной).
Эволюция системы во времени описывается гамильтоновыми уравнениями:
dθ/dt = ∂H/∂p = M^(-1) p
dp/dt = -∂H/∂θ = -∇U(θ) = ∇log(π(θ))
Эти уравнения определяют траекторию в расширенном пространстве (θ, p), вдоль которой общая энергия H сохраняется. Для генерации выборки в HMC используется следующий алгоритм:
- Случайно сгенерировать новый импульс p из распределения, соответствующего кинетической энергии (обычно нормального распределения);
- Начиная с текущего состояния (θ, p), смоделировать эволюцию системы по гамильтоновым уравнениям на протяжении некоторого времени τ, используя численное интегрирование (обычно метод leap-frog);
- В конце траектории получить предложенное состояние (θ’, p’);
- Принять или отклонить предложение с вероятностью Метрополиса: α = min(1, exp(H(θ, p) — H(θ’, p’))).
Важно отметить, что если бы мы могли интегрировать гамильтоновы уравнения точно, то H(θ’, p’) = H(θ, p), и вероятность принятия была бы равна 1. Однако, из-за численных ошибок интегрирования, необходим шаг Метрополиса для обеспечения корректности выборки.
Преимущества HMC с точки зрения физики
- Направленное движение: В отличие от случайного блуждания, HMC использует градиент логарифма целевой плотности (∇log(π(θ))), что позволяет эффективно двигаться в направлении увеличения плотности вероятности;
- Сохранение энергии: Траектории в HMC сохраняют общую энергию, что позволяет предложениям перемещаться на большие расстояния, сохраняя при этом высокую вероятность принятия;
- Адаптация к геометрии: Движение частицы естественным образом адаптируется к геометрии целевого распределения — она ускоряется в областях с крутым градиентом и замедляется в плоских областях;
- Эффективное исследование высокоразмерных пространств: HMC особенно эффективен в высокоразмерных пространствах, где традиционные методы MCMC страдают от «проклятия размерности».
Эта физическая интуиция помогает понять, почему HMC так эффективен в исследовании сложных распределений вероятностей, особенно в контексте байесовского вывода для нелинейных моделей с большим числом параметров.
Реализация с использованием PyMC: практический пример
Сейчас я реализую HMC с использованием библиотеки PyMC для решения практической задачи — байесовской логистической регрессии. Это классический пример нелинейной модели, для которой HMC демонстрирует значительные преимущества по сравнению с более простыми методами MCMC.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import pymc as pm
import arviz as az
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
# Создадим синтетические данные для логистической регрессии
def generate_logistic_data(n_samples=1000, n_features=5, random_seed=42):
np.random.seed(random_seed)
# Истинные коэффициенты
true_beta = np.random.normal(0, 1, n_features)
true_intercept = 0.5
# Генерация признаков с корреляцией
cov_matrix = np.eye(n_features)
# Добавляем корреляцию между некоторыми признаками
cov_matrix[0, 1] = cov_matrix[1, 0] = 0.7
cov_matrix[2, 3] = cov_matrix[3, 2] = 0.5
X = np.random.multivariate_normal(np.zeros(n_features), cov_matrix, n_samples)
# Вычисление линейного предиктора
linear_pred = true_intercept + X @ true_beta
# Вероятность принадлежности к положительному классу
p = 1 / (1 + np.exp(-linear_pred))
# Генерация меток классов
y = np.random.binomial(1, p)
return X, y, true_beta, true_intercept
# Генерация данных
X, y, true_beta, true_intercept = generate_logistic_data(n_samples=2000, n_features=5)
# Разделение на тренировочную и тестовую выборки
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Стандартизация признаков
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# Посмотрим на корреляционную матрицу признаков
plt.figure(figsize=(10, 8))
sns.heatmap(np.corrcoef(X_train.T), annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('Корреляционная матрица признаков')
plt.savefig('feature_correlation.png', dpi=300)
plt.show()
# Байесовская логистическая регрессия с использованием PyMC и HMC
def bayesian_logistic_regression(X, y, n_samples=2000, tune=1000, chains=4, cores=1):
with pm.Model() as logistic_model:
# Априорные распределения для параметров
beta = pm.Normal('beta', mu=0, sigma=1, shape=X.shape[1])
intercept = pm.Normal('intercept', mu=0, sigma=1)
# Линейный предиктор
linear_pred = intercept + pm.math.dot(X, beta)
# Функция связи (логистическая функция)
p = pm.math.sigmoid(linear_pred)
# Правдоподобие
y_obs = pm.Bernoulli('y_obs', p=p, observed=y)
# Семплирование с использованием NUTS (No U-Turn Sampler - улучшенная версия HMC)
trace = pm.sample(n_samples, tune=tune, chains=chains, cores=cores, return_inferencedata=True)
return trace, logistic_model
# Запуск байесовского вывода
trace, model = bayesian_logistic_regression(X_train_scaled, y_train)
# Анализ результатов
def analyze_results(trace, true_beta, true_intercept, X_test, y_test):
# Получение апостериорных выборок
posterior_beta = trace.posterior['beta'].values
posterior_intercept = trace.posterior['intercept'].values
# Reshaping для удобства
n_chains, n_samples, n_features = posterior_beta.shape
posterior_beta_flat = posterior_beta.reshape(n_chains * n_samples, n_features)
posterior_intercept_flat = posterior_intercept.reshape(n_chains * n_samples)
# Апостериорные средние
mean_beta = np.mean(posterior_beta_flat, axis=0)
mean_intercept = np.mean(posterior_intercept_flat)
# Апостериорные 95% доверительные интервалы
ci_beta = np.percentile(posterior_beta_flat, [2.5, 97.5], axis=0)
ci_intercept = np.percentile(posterior_intercept_flat, [2.5, 97.5])
# Таблица результатов
results = pd.DataFrame({
'Истинное значение': np.concatenate([[true_intercept], true_beta]),
'Апостериорное среднее': np.concatenate([[mean_intercept], mean_beta]),
'2.5% квантиль': np.concatenate([[ci_intercept[0]], ci_beta[0]]),
'97.5% квантиль': np.concatenate([[ci_intercept[1]], ci_beta[1]])
}, index=['intercept'] + [f'beta_{i}' for i in range(n_features)])
print("Результаты байесовского вывода:")
print(results)
# Визуализация апостериорных распределений
plt.figure(figsize=(12, 8))
# Визуализация для коэффициентов
for i in range(n_features):
plt.subplot(2, 3, i + 1)
sns.histplot(posterior_beta_flat[:, i], kde=True, color='#333333')
plt.axvline(true_beta[i], color='red', linestyle='--', label='Истинное значение')
plt.axvline(mean_beta[i], color='blue', linestyle='-', label='Апостериорное среднее')
plt.title(f'beta_{i}')
if i == 0:
plt.legend()
# Визуализация для интерсепта
plt.subplot(2, 3, 6)
sns.histplot(posterior_intercept_flat, kde=True, color='#333333')
plt.axvline(true_intercept, color='red', linestyle='--', label='Истинное значение')
plt.axvline(mean_intercept, color='blue', linestyle='-', label='Апостериорное среднее')
plt.title('intercept')
plt.tight_layout()
plt.savefig('posterior_distributions.png', dpi=300)
plt.show()
# Оценка качества предсказаний на тестовой выборке
n_pred_samples = min(500, n_chains * n_samples) # Ограничим количество выборок для предсказаний
pred_indices = np.random.choice(n_chains * n_samples, n_pred_samples, replace=False)
y_pred_proba = np.zeros((n_pred_samples, X_test.shape[0]))
for i, idx in enumerate(pred_indices):
chain_idx = idx // n_samples
sample_idx = idx % n_samples
beta_sample = posterior_beta[chain_idx, sample_idx]
intercept_sample = posterior_intercept[chain_idx, sample_idx]
linear_pred = intercept_sample + X_test @ beta_sample
y_pred_proba[i] = 1 / (1 + np.exp(-linear_pred))
# Средняя вероятность
mean_proba = np.mean(y_pred_proba, axis=0)
# Предсказанные классы
y_pred = (mean_proba > 0.5).astype(int)
# Оценка точности
accuracy = np.mean(y_pred == y_test)
print(f"Точность на тестовой выборке: {accuracy:.4f}")
# ROC-кривая
from sklearn.metrics import roc_curve, auc
fpr, tpr, thresholds = roc_curve(y_test, mean_proba)
roc_auc = auc(fpr, tpr)
plt.figure(figsize=(10, 8))
plt.plot(fpr, tpr, color='#333333', lw=2, label=f'ROC curve (area = {roc_auc:.3f})')
plt.plot([0, 1], [0, 1], color='#999999', lw=1, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC кривая байесовской логистической регрессии')
plt.legend(loc="lower right")
plt.savefig('roc_curve.png', dpi=300)
plt.show()
# Неопределенность в предсказаниях
plt.figure(figsize=(10, 8))
# Сортировка по средней вероятности для лучшей визуализации
sort_indices = np.argsort(mean_proba)
sorted_mean_proba = mean_proba[sort_indices]
sorted_y_test = y_test[sort_indices]
sorted_y_pred_proba = y_pred_proba[:, sort_indices]
# Отображаем только подмножество тестовых примеров для ясности
n_display = min(100, len(y_test))
step = len(y_test) // n_display
display_indices = range(0, len(y_test), step)
# Отображение неопределенности
plt.errorbar(range(len(display_indices)),
sorted_mean_proba[display_indices],
yerr=np.std(sorted_y_pred_proba[:, display_indices], axis=0),
fmt='o', color='#333333', ecolor='#999999', elinewidth=1, capsize=3)
# Отображение истинных меток
for i, idx in enumerate(display_indices):
plt.scatter(i, sorted_y_test[idx], color='red', marker='x', s=50)
plt.axhline(0.5, color='#999999', linestyle='--')
plt.xlabel('Индекс тестового примера')
plt.ylabel('Вероятность положительного класса')
plt.title('Байесовские предсказания с неопределенностью')
plt.savefig('prediction_uncertainty.png', dpi=300)
plt.show()
return results, accuracy, roc_auc
# Анализ результатов HMC
results_hmc, accuracy_hmc, roc_auc_hmc = analyze_results(trace, true_beta, true_intercept, X_test_scaled, y_test)
# Диагностика MCMC
def mcmc_diagnostics(trace, algorithm_name="MCMC"):
# Диаграмма трассировки
az.plot_trace(trace, var_names=['beta', 'intercept'])
plt.suptitle(f'Трассировки для {algorithm_name}')
plt.savefig(f'trace_plot_{algorithm_name.lower()}.png', dpi=300)
plt.show()
# Эффективный размер выборки
ess = az.ess(trace, var_names=['beta', 'intercept'])
print(f"Эффективный размер выборки для {algorithm_name}:")
print(ess)
# Диагностика Гельмана-Рубина (должна быть близка к 1)
rhat = az.rhat(trace, var_names=['beta', 'intercept'])
print(f"\nRhat статистика для {algorithm_name}:")
print(rhat)
# Автокорреляция
az.plot_autocorr(trace, var_names=['beta', 'intercept'])
plt.suptitle(f'Автокорреляция для {algorithm_name}')
plt.savefig(f'autocorrelation_{algorithm_name.lower()}.png', dpi=300)
plt.show()
# Парные графики для оценки корреляций между параметрами
az.plot_pair(trace, var_names=['beta'], kind='hexbin',
divergences=True, figsize=(12, 10))
plt.suptitle(f'Парные графики для {algorithm_name}')
plt.savefig(f'pair_plot_{algorithm_name.lower()}.png', dpi=300)
plt.show()
# Энергетическая диагностика - только для HMC-подобных алгоритмов
if hasattr(trace.sample_stats, 'energy'):
az.plot_energy(trace)
plt.suptitle(f'Энергетическая диагностика для {algorithm_name}')
plt.savefig(f'energy_plot_{algorithm_name.lower()}.png', dpi=300)
plt.show()
else:
print(f"Энергетическая диагностика недоступна для {algorithm_name}")
# Диагностика HMC
print("=== Диагностика HMC ===")
mcmc_diagnostics(trace, "HMC")
# Сравнение с Metropolis-Hastings (для демонстрации преимуществ HMC)
def metropolis_hastings_logistic(X, y, n_samples=2000, tune=1000):
with pm.Model() as mh_model:
# Те же априорные распределения
beta = pm.Normal('beta', mu=0, sigma=1, shape=X.shape[1])
intercept = pm.Normal('intercept', mu=0, sigma=1)
# Линейный предиктор
linear_pred = intercept + pm.math.dot(X, beta)
# Функция связи
p = pm.math.sigmoid(linear_pred)
# Правдоподобие
y_obs = pm.Bernoulli('y_obs', p=p, observed=y)
# Используем Metropolis вместо NUTS
trace_mh = pm.sample(n_samples, tune=tune, step=pm.Metropolis(),
chains=2, cores=1, return_inferencedata=True)
return trace_mh, mh_model
# Запуск Metropolis-Hastings (это может занять больше времени)
trace_mh, model_mh = metropolis_hastings_logistic(X_train_scaled, y_train, n_samples=2000, tune=1000)
# Анализ результатов Metropolis-Hastings
results_mh, accuracy_mh, roc_auc_mh = analyze_results(trace_mh, true_beta, true_intercept, X_test_scaled, y_test)
# Диагностика для Metropolis-Hastings
print("\n=== Диагностика Metropolis-Hastings ===")
mcmc_diagnostics(trace_mh, "Metropolis-Hastings")
# Дополнительная диагностика принятия для Metropolis-Hastings
def metropolis_acceptance_diagnostics(trace_mh):
"""Специальная диагностика для Metropolis-Hastings"""
# Попробуем получить статистики принятия, если они доступны
if hasattr(trace_mh.sample_stats, 'accepted'):
acceptance_rate = np.mean(trace_mh.sample_stats.accepted)
print(f"Коэффициент принятия Metropolis-Hastings: {acceptance_rate:.3f}")
# График коэффициента принятия по итерациям
plt.figure(figsize=(10, 6))
for chain in range(trace_mh.sample_stats.dims['chain']):
accepted = trace_mh.sample_stats.accepted.sel(chain=chain)
# Скользящее среднее для коэффициента принятия
window_size = 100
rolling_acceptance = np.convolve(accepted, np.ones(window_size)/window_size, mode='valid')
plt.plot(rolling_acceptance, label=f'Цепь {chain}', alpha=0.7)
plt.axhline(y=0.234, color='red', linestyle='--',
label='Оптимальный коэффициент (23.4%)')
plt.xlabel('Итерация')
plt.ylabel('Скользящий коэффициент принятия')
plt.title('Эволюция коэффициента принятия Metropolis-Hastings')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('metropolis_acceptance_rate.png', dpi=300)
plt.show()
else:
print("Статистики принятия недоступны для данного trace")
# Диагностика принятия для Metropolis-Hastings
metropolis_acceptance_diagnostics(trace_mh)
# Сравнение эффективности HMC и Metropolis-Hastings
def compare_efficiency(trace_hmc, trace_mh):
# Эффективный размер выборки для обоих методов
ess_hmc = az.ess(trace_hmc, var_names=['beta', 'intercept'])
ess_mh = az.ess(trace_mh, var_names=['beta', 'intercept'])
# Среднее значение ESS
mean_ess_hmc = np.mean([ess_hmc.sel(variable='beta').values.flatten(),
ess_hmc.sel(variable='intercept').values.flatten()])
mean_ess_mh = np.mean([ess_mh.sel(variable='beta').values.flatten(),
ess_mh.sel(variable='intercept').values.flatten()])
# Относительная эффективность
relative_efficiency = mean_ess_hmc / mean_ess_mh
print(f"Средний эффективный размер выборки для HMC: {mean_ess_hmc:.2f}")
print(f"Средний эффективный размер выборки для Metropolis-Hastings: {mean_ess_mh:.2f}")
print(f"Относительная эффективность HMC: {relative_efficiency:.2f}x")
# Сравнение автокорреляции
_, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
# Автокорреляция для первого параметра
max_lag = 20
# HMC
chain_hmc = trace_hmc.posterior['beta'].values.reshape(-1, X.shape[1])[:, 0]
acf_hmc = np.array([np.corrcoef(chain_hmc[:-lag], chain_hmc[lag:])[0, 1]
for lag in range(1, max_lag + 1)])
# Metropolis-Hastings
chain_mh = trace_mh.posterior['beta'].values.reshape(-1, X.shape[1])[:, 0]
acf_mh = np.array([np.corrcoef(chain_mh[:-lag], chain_mh[lag:])[0, 1]
for lag in range(1, max_lag + 1)])
# Визуализация
ax1.bar(range(1, max_lag + 1), acf_hmc, color='#333333', alpha=0.7, label='HMC')
ax1.set_title('Автокорреляция для HMC')
ax1.set_xlabel('Лаг')
ax1.set_ylabel('Автокорреляция')
ax1.axhline(y=0, linestyle='--', color='#999999')
ax2.bar(range(1, max_lag + 1), acf_mh, color='#999999', alpha=0.7, label='Metropolis-Hastings')
ax2.set_title('Автокорреляция для Metropolis-Hastings')
ax2.set_xlabel('Лаг')
ax2.set_ylabel('Автокорреляция')
ax2.axhline(y=0, linestyle='--', color='#999999')
plt.tight_layout()
plt.savefig('autocorrelation_comparison.png', dpi=300)
plt.show()
return relative_efficiency
# Сравнение эффективности
relative_eff = compare_efficiency(trace, trace_mh)
Рис. 6: Корреляционная матрица признаков синтетического датасета. Наблюдается умеренная корреляция между признаками 0-1 (r=0.72) и 2-3 (r=0.53), что создает реалистичные условия для тестирования байесовских методов
Рис. 7: Апостериорные распределения параметров байесовской логистической регрессии. Красные пунктирные линии показывают истинные значения параметров, синие сплошные — апостериорные средние
Рис. 8: ROC-кривая байесовской логистической регрессии на тестовой выборке
Рис. 9: Байесовские предсказания с квантификацией неопределенности. Черные точки с планками погрешностей — средние вероятности и их стандартные отклонения, красные крестики — истинные метки классов. Пунктирная линия показывает порог классификации (0.5)
Рис. 10: Диагностика трассировок MCMC для алгоритма HMC (NUTS)
Рис. 11: Парные графики апостериорных выборок коэффициентов регрессии
Рис. 12: Энергетическая диагностика HMC. Совпадение распределений энергии переходов и стационарной энергии свидетельствует об отсутствии проблем с геометрией целевого распределения
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.74 7 348.60 draws/s 0:00:08 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.60 3 192.27 draws/s 0:00:15 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.88 7 138.86 draws/s 0:00:21 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.77 7 104.46 draws/s 0:00:28 0:00:00
Результаты байесовского вывода:
Истинное значение Апостериорное среднее 2.5% квантиль \
intercept 0.500000 0.446489 0.323415
beta_0 0.496714 0.618128 0.434850
beta_1 -0.138264 -0.328945 -0.509488
beta_2 0.647689 0.686349 0.533827
beta_3 1.523030 1.311625 1.129142
beta_4 -0.234153 -0.189631 -0.312734
97.5% квантиль
intercept 0.572482
beta_0 0.808265
beta_1 -0.151950
beta_2 0.839836
beta_3 1.499253
beta_4 -0.066887
Точность на тестовой выборке: 0.7700
Эффективный размер выборки для HMC:
Size: 88B
Dimensions: (beta_dim_0: 5)
Coordinates:
* beta_dim_0 (beta_dim_0) int64 40B 0 1 2 3 4
Data variables:
beta (beta_dim_0) float64 40B 7.344e+03 7.681e+03 ... 9.654e+03
intercept float64 8B 9.578e+03
Rhat статистика для HMC:
Size: 88B
Dimensions: (beta_dim_0: 5)
Coordinates:
* beta_dim_0 (beta_dim_0) int64 40B 0 1 2 3 4
Data variables:
beta (beta_dim_0) float64 40B 1.001 1.001 1.0 1.001 1.0
intercept float64 8B 1.0 
Рис. 13: Апостериорные распределения параметров байесовской логистической регрессии (Metropolis-Hastings). Сравнение с HMC показывает схожие результаты восстановления параметров, но с большей вариабельностью
Рис. 14: ROC-кривая байесовской логистической регрессии (Metropolis-Hastings) на тестовой выборке
Рис. 15: Байесовские предсказания с квантификацией неопределенности (Metropolis-Hastings)
Рис. 16: Байесовские предсказания с квантификацией неопределенности (Metropolis-Hastings)
Рис. 17: Автокорреляционные функции параметров для Metropolis-Hastings. Медленное убывание автокорреляции свидетельствует о менее эффективном семплировании по сравнению с HMC
Рис. 18: Парные графики апостериорных выборок коэффициентов регрессии (Metropolis-Hastings). Наблюдаются менее гладкие распределения и хаотичные артефакты семплирования по сравнению с HMC
Этот код демонстрирует полный рабочий процесс байесовской логистической регрессии с использованием HMC (реализованного через NUTS — No U-Turn Sampler, улучшенную версию HMC) и Metropolis-Hastings для сравнения.
Обратите внимание на следующие ключевые аспекты:
- Генерация данных: Мы создаем синтетические данные с коррелированными признаками, что делает задачу более сложной для традиционных методов MCMC;
- Реализация моделей: Мы используем PyMC для создания байесовской логистической регрессии и запуска MCMC. PyMC автоматически выбирает NUTS (улучшенную версию HMC) по умолчанию, но мы также явно указываем Metropolis для сравнения;
- Анализ результатов: Мы оцениваем точность модели на тестовой выборке, визуализируем апостериорные распределения параметров и анализируем неопределенность в предсказаниях;
- Диагностика MCMC: Мы используем различные диагностические метрики, включая трассировки, эффективный размер выборки, Rhat-статистику, автокорреляцию и энергетическую диагностику, для оценки качества семплирования;
- Сравнение эффективности: Мы сравниваем эффективность HMC и Metropolis-Hastings с точки зрения эффективного размера выборки и автокорреляции.
HMC обычно демонстрирует значительные преимущества по сравнению с Metropolis-Hastings, особенно для:
- Коррелированных параметров: HMC эффективно исследует области с сильными корреляциями между параметрами;
- Высокоразмерных моделей: Эффективность HMC растет с увеличением размерности, в то время как Metropolis-Hastings страдает от «проклятия размерности»;
- Эффективного размера выборки: HMC обычно обеспечивает гораздо больший эффективный размер выборки на то же количество итераций;
- Автокорреляции: Выборки, полученные с помощью HMC, имеют значительно меньшую автокорреляцию.
Вышеуказанный код представляет собой комплексный пример использования HMC для байесовского вывода, который можно адаптировать для широкого спектра задач в анализе данных и машинном обучении.
Практические рекомендации по использованию HMC
За годы работы с HMC и его вариантами, я выработал ряд практических рекомендаций, которые помогут вам эффективно использовать этот мощный метод:
1. Настройка параметров HMC
HMC имеет два основных параметра, которые критически влияют на его эффективность:
- Размер шага (step size, ε): Контролирует точность численного интегрирования. Слишком большой шаг приводит к численной нестабильности и низкому коэффициенту принятия, слишком маленький — к неэффективному исследованию пространства;
- Количество шагов (number of leapfrog steps, L): Определяет длину траектории. Слишком короткие траектории приводят к случайному блужданию, слишком длинные — к избыточным вычислениям и потенциальному возвращению в окрестность начальной точки.
Современные реализации HMC, такие как NUTS (No U-Turn Sampler), автоматически адаптируют эти параметры, что значительно упрощает использование метода. Тем не менее, важно понимать их влияние на производительность.
2. Рекомендации по предобработке данных
Правильная предобработка данных крайне важна для эффективности HMC:
- Стандартизация или нормализация признаков: Приведение всех признаков к одному масштабу значительно улучшает сходимость HMC;
- Параметризация модели: Хорошо спроектированная параметризация может уменьшить корреляции между параметрами и улучшить геометрию пространства параметров;
- Центрирование данных: Для многих моделей центрирование данных вокруг нуля улучшает сходимость;
- Избегание сильно коррелированных признаков: При наличии сильных корреляций, рассмотрите возможность использования методов снижения размерности или ортогонализации признаков.
3. Диагностика и мониторинг
Важно регулярно проверять качество полученных выборок:
- Трассировки параметров: Визуальный осмотр трассировок для проверки хорошего «перемешивания»;
- Эффективный размер выборки (ESS): Чем выше, тем лучше. Значения ESS ниже 100 обычно указывают на проблемы с семплированием;
- Метрика Gelman-Rubin (Rhat): Должна быть близка к 1 для всех параметров. Значения выше 1.1 указывают на проблемы со сходимостью;
- Энергетическая диагностика: NUTS предоставляет диагностику на основе сохранения энергии, которая помогает выявить проблемы с геометрией целевого распределения;
- Дивергенции: Высокое количество дивергентных переходов в NUTS указывает на проблемы с геометрией пространства параметров, часто требующие репараметризации модели.
4. Репараметризация моделей
Одним из наиболее эффективных способов улучшения работы HMC является репараметризация модели:
- Некоррелированная параметризация: Трансформация параметров для уменьшения корреляций;
- Центрированная vs. нецентрированная параметризация: Для иерархических моделей нецентрированная параметризация часто более эффективна, особенно при малом количестве данных;
- Логарифмическая трансформация: Для параметров с положительными ограничениями, таких как дисперсии или масштабные параметры;
- Трансформация Фишера для корреляций: Для корреляционных матриц использование z-трансформации вместо прямого моделирования корреляций.
5. Вычислительные оптимизации
HMC требует вычисления градиента логарифма целевой плотности, что может быть вычислительно затратно:
- Автоматическое дифференцирование: Используйте библиотеки с эффективным автоматическим дифференцированием, такие как PyMC, TensorFlow Probability или Pyro;
- Векторизация вычислений: Максимальное использование векторизованных операций для ускорения вычислений;
- Параллельные цепи: Запуск нескольких цепей параллельно для эффективного использования многоядерных процессоров;
- Использование GPU: Для очень больших моделей рассмотрите возможность использования GPU для ускорения вычислений.
6. Расширенные варианты HMC
Существует несколько улучшенных вариантов HMC, которые стоит рассмотреть для сложных моделей:
- NUTS (No U-Turn Sampler): Адаптивный вариант HMC, который автоматически определяет оптимальную длину траектории. Это стандартный выбор в большинстве современных библиотек;
- Римановский HMC: Использует информацию о локальной геометрии пространства параметров (через метрику Фишера или гессиан) для более эффективного исследования. Особенно полезен для моделей с сильными нелинейностями;
- Адаптивная HMC: Автоматически настраивает матрицу масс и размер шага во время прогрева для оптимизации эффективности;
- Предобусловленный HMC: Использует предобуславливание для улучшения геометрических свойств целевого распределения.
7. Работа с ограничениями и граничными условиями
Многие статистические модели включают параметры с ограничениями (например, положительные дисперсии, вероятности в диапазоне [0,1]). Существует два основных подхода к работе с такими ограничениями в HMC:
- Репараметризация: Трансформация ограниченных параметров в неограниченное пространство. Например, использование логарифмической трансформации для положительных параметров или логит-трансформации для вероятностей;
- Отражающие границы: Модификация алгоритма HMC для корректной обработки отражений от границ допустимой области.
Большинство современных библиотек, таких как PyMC, автоматически применяют соответствующие трансформации, что упрощает работу с ограниченными параметрами.
Применение MCMC в практических задачах
Байесовская оптимизация портфеля активов
Одной из наиболее интересных областей применения методов MCMC является финансовое моделирование, и в частности, байесовская оптимизация портфеля активов. В отличие от классических подходов, таких как оптимизация Марковица, байесовский подход позволяет естественным образом учитывать неопределенность в оценках ожидаемой доходности и ковариационной матрице, а также включать априорные представления о рынке.
Рассмотрим задачу оптимизации портфеля с использованием байесовского подхода и методов MCMC. В портфеле будут 5 ETF :
- VTI — Vanguard Total Stock Market ETF;
- BND — Vanguard Total Bond Market ETF;
- VEA — Vanguard Developed Markets ETF;
- GLD — SPDR Gold Shares ETF;
- VNQ — Vanguard Real Estate ETF.
Мы будем анализировать их историческую доходность за последние 2 года.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pymc as pm
import arviz as az
from scipy import stats
import yfinance as yf
import pytensor.tensor as pt
# Загрузка данных о ценах активов
def load_portfolio_data(tickers=['VTI', 'BND', 'VEA', 'GLD', 'VNQ'],
start_date='2023-07-01', end_date='2025-07-01'):
data = yf.download(tickers, start=start_date, end=end_date)['Close']
returns = data.pct_change().dropna()
return returns, data
# Загружаем данные
returns, prices = load_portfolio_data()
# Визуализация исторических доходностей
plt.figure(figsize=(12, 6))
(returns + 1).cumprod().plot()
plt.title('Кумулятивная доходность активов')
plt.ylabel('Кумулятивный рост')
plt.grid(True, alpha=0.3)
plt.savefig('cumulative_returns.png', dpi=300)
plt.show()
# Статистика доходностей
returns_stats = pd.DataFrame({
'Среднее (%)': returns.mean() * 100,
'Ст. отклонение (%)': returns.std() * 100,
'Мин. (%)': returns.min() * 100,
'Макс. (%)': returns.max() * 100,
'Skewness': returns.skew(),
'Kurtosis': returns.kurtosis()
})
print("Статистика доходностей активов:")
print(returns_stats)
# Визуализация корреляций
plt.figure(figsize=(10, 8))
sns.heatmap(returns.corr(), annot=True, cmap='coolwarm', vmin=-1, vmax=1)
plt.title('Корреляционная матрица доходностей')
plt.savefig('returns_correlation.png', dpi=300)
plt.show()
# Оптимизированная байесовская модель с упрощенной структурой
def bayesian_returns_estimation(returns, n_samples=1000, tune=500):
n_assets = returns.shape[1]
with pm.Model() as model:
# Упрощенные априорные распределения
mu = pm.Normal('mu', mu=0, sigma=0.02, shape=n_assets)
# Используем готовую упакованную LKJ ковариацию
packed_L = pm.LKJCholeskyCov('packed_L',
n=n_assets,
eta=3.0,
sd_dist=pm.HalfNormal.dist(sigma=0.05))
# Извлекаем компоненты
chol = packed_L[0] # Холецкий фактор
# corr = packed_L[1] # Корреляционная матрица (если нужна)
# stds = packed_L[2] # Стандартные отклонения (если нужны)
# Правдоподобие с использованием Холецкого разложения
observed = pm.MvNormal('observed', mu=mu, chol=chol, observed=returns.values)
# Настройки семплирования
trace = pm.sample(n_samples, tune=tune,
chains=2,
cores=1,
target_accept=0.9,
return_inferencedata=True)
return trace, model
# Запуск байесовской оценки параметров
print("Запуск MCMC семплирования...")
trace, model = bayesian_returns_estimation(returns)
# Упрощенный анализ апостериорных распределений
def analyze_posterior_fast(trace, returns):
asset_names = returns.columns
# Извлечение апостериорных выборок
posterior_mu = trace.posterior['mu'].values.reshape(-1, len(asset_names))
# Извлекаем стандартные отклонения и корреляции из упакованной структуры
posterior_stds = trace.posterior['packed_L_stds'].values.reshape(-1, len(asset_names))
posterior_corr = trace.posterior['packed_L_corr'].values.reshape(-1, len(asset_names), len(asset_names))
# Апостериорные средние для ожидаемых доходностей
mu_mean = np.mean(posterior_mu, axis=0)
mu_std = np.std(posterior_mu, axis=0)
# Таблица с результатами
mu_results = pd.DataFrame({
'Актив': asset_names,
'Апост. ср. доходность (%)': mu_mean * 100,
'Апост. ст. откл. (%)': mu_std * 100,
'Эмпирическая доходность (%)': returns.mean() * 100
}).set_index('Актив')
print("Апостериорные оценки ожидаемых доходностей:")
print(mu_results)
# Строим графики
plt.figure(figsize=(12, 4))
# Апостериорные распределения доходностей
plt.subplot(1, 2, 1)
for i, name in enumerate(asset_names):
plt.hist(posterior_mu[:, i] * 100, alpha=0.6, label=name, bins=20)
plt.axvline(0, color='black', linestyle='--')
plt.title('Апостериорные распределения доходностей (%)')
plt.xlabel('Доходность (%)')
plt.ylabel('Частота')
plt.legend()
plt.grid(True, alpha=0.3)
# Апостериорная корреляционная матрица
plt.subplot(1, 2, 2)
mean_corr = np.mean(posterior_corr, axis=0)
sns.heatmap(mean_corr, annot=True, cmap='coolwarm', vmin=-1, vmax=1,
xticklabels=asset_names, yticklabels=asset_names)
plt.title('Апостериорная корреляционная матрица')
plt.tight_layout()
plt.savefig('posterior_analysis.png', dpi=300)
plt.show()
return mu_results, mean_corr, posterior_mu, posterior_stds, posterior_corr
# Анализ результатов
mu_results, mean_corr, posterior_mu, posterior_stds, posterior_corr = analyze_posterior_fast(trace, returns)
# Байесовская оптимизация портфеля
def portfolio_optimization_full(posterior_mu, posterior_stds, posterior_corr,
risk_aversion=2.0):
n_samples = posterior_mu.shape[0] # Используем все выборки
n_assets = posterior_mu.shape[1]
# Векторизованное вычисление ковариационных матриц
sigma_diag = np.array([np.diag(posterior_stds[i]) for i in range(n_samples)])
cov_matrices = np.array([sigma_diag[i] @ posterior_corr[i] @ sigma_diag[i] for i in range(n_samples)])
# Векторизованная оптимизация
optimal_weights = np.zeros((n_samples, n_assets))
for i in range(n_samples):
mu = posterior_mu[i]
cov = cov_matrices[i]
try:
# Регуляризация для численной стабильности
cov_reg = cov + np.eye(n_assets) * 1e-6
inv_cov = np.linalg.inv(cov_reg)
weights = inv_cov @ mu / risk_aversion
# Применяем ограничения (только длинные позиции)
weights = np.maximum(weights, 0)
if np.sum(weights) > 0:
weights = weights / np.sum(weights)
else:
weights = np.ones(n_assets) / n_assets
except np.linalg.LinAlgError:
weights = np.ones(n_assets) / n_assets
optimal_weights[i] = weights
# Вычисляем метрики для всех выборок
optimal_returns = np.array([np.sum(optimal_weights[i] * posterior_mu[i]) for i in range(n_samples)])
optimal_risks = np.array([np.sqrt(optimal_weights[i] @ cov_matrices[i] @ optimal_weights[i])
for i in range(n_samples)])
# Средние результаты
mean_weights = np.mean(optimal_weights, axis=0)
std_weights = np.std(optimal_weights, axis=0)
return {
'optimal_weights': optimal_weights,
'mean_weights': mean_weights,
'std_weights': std_weights,
'optimal_returns': optimal_returns,
'optimal_risks': optimal_risks,
'mean_portfolio_return': np.mean(optimal_returns),
'mean_portfolio_risk': np.mean(optimal_risks)
}
# Оптимизация портфеля
print("Оптимизация портфеля...")
portfolio_results = portfolio_optimization_full(posterior_mu, posterior_stds, posterior_corr)
# Анализ результатов
def analyze_portfolio_results_fast(portfolio_results, returns, posterior_mu, posterior_stds):
asset_names = returns.columns
mean_weights = portfolio_results['mean_weights']
std_weights = portfolio_results['std_weights']
# Таблица весов
weights_df = pd.DataFrame({
'Актив': asset_names,
'Оптимальный вес (%)': mean_weights * 100,
'Ст. откл. веса (%)': std_weights * 100
}).set_index('Актив')
print("Оптимальные веса портфеля:")
print(weights_df)
# Основные графики
plt.figure(figsize=(15, 8))
# Оптимальные веса
plt.subplot(2, 2, 1)
plt.bar(asset_names, mean_weights * 100, yerr=std_weights * 100,
color='#333333', alpha=0.7, capsize=5)
plt.title('Оптимальные веса портфеля')
plt.ylabel('Вес (%)')
plt.xticks(rotation=45)
plt.grid(True, alpha=0.3)
# Распределение доходности и риска
plt.subplot(2, 2, 2)
optimal_returns = portfolio_results['optimal_returns'] * 100
optimal_risks = portfolio_results['optimal_risks'] * 100
plt.scatter(optimal_risks, optimal_returns, alpha=0.6, color='#333333')
plt.scatter(np.mean(optimal_risks), np.mean(optimal_returns),
color='red', marker='*', s=200, label='Среднее')
plt.title('Риск-доходность портфеля')
plt.xlabel('Риск (%)')
plt.ylabel('Доходность (%)')
plt.legend()
plt.grid(True, alpha=0.3)
# Распределение весов
plt.subplot(2, 2, 3)
optimal_weights = portfolio_results['optimal_weights'] * 100
for i, name in enumerate(asset_names):
plt.hist(optimal_weights[:, i], alpha=0.6, label=name, bins=15)
plt.title('Распределения весов')
plt.xlabel('Вес (%)')
plt.ylabel('Частота')
plt.legend()
plt.grid(True, alpha=0.3)
# Эффективная граница
plt.subplot(2, 2, 4)
risk_aversions = np.linspace(1, 8, 8)
frontier_returns = []
frontier_risks = []
for risk_aversion in risk_aversions:
results = portfolio_optimization_full(posterior_mu, posterior_stds, posterior_corr,
risk_aversion=risk_aversion)
frontier_returns.append(results['mean_portfolio_return'] * 100)
frontier_risks.append(results['mean_portfolio_risk'] * 100)
plt.plot(frontier_risks, frontier_returns, 'o-', color='#333333', label='Эффективная граница')
plt.scatter(np.mean(optimal_risks), np.mean(optimal_returns),
color='red', marker='*', s=200, label='Выбранный портфель')
# Отдельные активы
for i, name in enumerate(asset_names):
mean_return = np.mean(posterior_mu[:, i]) * 100
mean_risk = np.mean(posterior_stds[:, i]) * 100
plt.scatter(mean_risk, mean_return, label=name, s=50)
plt.title('Эффективная граница')
plt.xlabel('Риск (%)')
plt.ylabel('Доходность (%)')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('portfolio_analysis_fast.png', dpi=300, bbox_inches='tight')
plt.show()
# График 1: Кумулятивная доходность портфеля vs отдельных активов
plt.figure(figsize=(12, 6))
portfolio_returns = returns.dot(mean_weights)
cumulative_returns = (1 + portfolio_returns).cumprod()
cumulative_returns.plot(label='Оптимальный портфель', color='#333333', linewidth=3)
(1 + returns).cumprod().plot(alpha=0.7)
plt.title('Кумулятивная доходность портфеля vs отдельных активов', fontsize=14)
plt.ylabel('Кумулятивный рост', fontsize=12)
plt.xlabel('Дата', fontsize=12)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('portfolio_cumulative_returns.png', dpi=300, bbox_inches='tight')
plt.show()
# График 2: Ежемесячная доходность оптимального портфеля
plt.figure(figsize=(14, 6))
monthly_returns = portfolio_returns.resample('M').apply(lambda x: (1 + x).prod() - 1) * 100
colors = ['green' if x >= 0 else 'red' for x in monthly_returns]
bars = plt.bar(range(len(monthly_returns)), monthly_returns, color=colors, alpha=0.7)
plt.title('Ежемесячная доходность оптимального портфеля', fontsize=14)
plt.ylabel('Доходность (%)', fontsize=12)
plt.xlabel('Месяц', fontsize=12)
plt.axhline(y=0, color='black', linestyle='-', alpha=0.5, linewidth=1)
plt.grid(True, alpha=0.3)
# Настраиваем подписи по оси X
tick_positions = range(0, len(monthly_returns), max(1, len(monthly_returns)//12))
tick_labels = [monthly_returns.index[i].strftime('%Y-%m') for i in tick_positions]
plt.xticks(tick_positions, tick_labels, rotation=45)
plt.tight_layout()
plt.savefig('portfolio_monthly_returns.png', dpi=300, bbox_inches='tight')
plt.show()
# График 3: Просадка портфеля
plt.figure(figsize=(12, 6))
cumulative_max = cumulative_returns.cummax()
drawdown = (cumulative_returns / cumulative_max - 1) * 100
drawdown.plot(color='red', linewidth=2)
plt.fill_between(drawdown.index, drawdown, 0, color='red', alpha=0.3)
plt.title('Просадка портфеля', fontsize=14)
plt.ylabel('Просадка (%)', fontsize=12)
plt.xlabel('Дата', fontsize=12)
plt.axhline(y=0, color='black', linestyle='-', alpha=0.5, linewidth=1)
plt.grid(True, alpha=0.3)
# Добавляем аннотацию максимальной просадки
max_dd_value = drawdown.min()
max_dd_date = drawdown.idxmin()
plt.annotate(f'Макс. просадка: {max_dd_value:.1f}%',
xy=(max_dd_date, max_dd_value),
xytext=(max_dd_date, max_dd_value-2),
arrowprops=dict(arrowstyle='->', color='black'),
fontsize=10, ha='center')
plt.tight_layout()
plt.savefig('portfolio_drawdown.png', dpi=300, bbox_inches='tight')
plt.show()
# Отдельный график метрик производительности
plt.figure(figsize=(10, 6))
# Вычисляем метрики
annual_return = (cumulative_returns.iloc[-1] ** (252 / len(cumulative_returns)) - 1) * 100
annual_volatility = portfolio_returns.std() * np.sqrt(252) * 100
sharpe_ratio = annual_return / annual_volatility if annual_volatility > 0 else 0
max_drawdown = drawdown.min()
# Создаем таблицу метрик
performance_metrics = pd.DataFrame({
'Метрика': ['Годовая доходность (%)', 'Годовая волатильность (%)',
'Коэффициент Шарпа', 'Максимальная просадка (%)'],
'Значение': [annual_return, annual_volatility, sharpe_ratio, max_drawdown]
}).set_index('Метрика')
print("\nМетрики производительности оптимального портфеля:")
print(performance_metrics)
# Визуализация метрик
metrics = performance_metrics.index
values = performance_metrics['Значение']
colors = ['green' if v > 0 else 'red' for v in values]
bars = plt.bar(metrics, values, color=colors, alpha=0.7)
plt.title('Метрики производительности оптимального портфеля', fontsize=14)
plt.ylabel('Значение', fontsize=12)
plt.xticks(rotation=45, ha='right')
plt.grid(True, alpha=0.3)
# Добавляем значения на столбцы
for bar, value in zip(bars, values):
height = bar.get_height()
plt.text(bar.get_x() + bar.get_width()/2., height + (0.5 if height >= 0 else -1),
f'{value:.2f}', ha='center', va='bottom' if height >= 0 else 'top')
plt.tight_layout()
plt.savefig('portfolio_metrics.png', dpi=300, bbox_inches='tight')
plt.show()
return weights_df
# Анализ результатов
weights_df = analyze_portfolio_results_fast(portfolio_results, returns, posterior_mu, posterior_stds)
# Быстрая диагностика MCMC
def mcmc_diagnostics_fast(trace):
print("=== Быстрая диагностика MCMC ===")
# Только основные метрики
ess = az.ess(trace, var_names=['mu'])
rhat = az.rhat(trace, var_names=['mu'])
print(f"Средний эффективный размер выборки: {np.mean(ess.mu):.0f}")
print(f"Максимальный R-hat: {np.max(rhat.mu):.3f}")
if np.max(rhat.mu) < 1.1:
print("✓ Сходимость хорошая (R-hat < 1.1)") else: print("⚠ Возможные проблемы со сходимостью (R-hat >= 1.1)")
# Упрощенная трассировка
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
mu_samples = trace.posterior['mu'].values
for i in range(min(3, mu_samples.shape[-1])): # Показываем только первые 3 актива
for chain in range(mu_samples.shape[0]):
plt.plot(mu_samples[chain, :, i], alpha=0.7, label=f'Asset {i}, Chain {chain}' if i == 0 else "")
plt.title('Трассировки (первые 3 актива)')
plt.xlabel('Итерация')
plt.ylabel('Значение')
if mu_samples.shape[-1] > 0:
plt.legend()
plt.subplot(1, 2, 2)
az.plot_ess(trace, var_names=['mu'], kind='local', ax=plt.gca())
plt.title('Локальный эффективный размер выборки')
plt.tight_layout()
plt.savefig('mcmc_diagnostics_fast.png', dpi=300)
plt.show()
# Диагностика
mcmc_diagnostics_fast(trace)
Рис. 19: Историческая кумулятивная доходность ETF активов портфеля за 2 года
Рис. 20: Эмпирическая корреляционная матрица дневных доходностей активов за период 2023-2025. Высокая корреляция между акциями США и развитых рынков (VTI-VEA), слабая корреляция облигаций и золота с рисковыми активами подтверждает их диверсификационную ценность
Статистика доходностей активов:
Среднее (%) Ст. отклонение (%) Мин. (%) Макс. (%) Skewness Kurtosis
Ticker
BND 0.018201 0.373042 -1.205624 1.277082 -0.047779 0.392389
GLD 0.112249 0.992289 -3.568330 3.699125 -0.032621 1.126554
VEA 0.059395 0.947142 -6.329876 7.469883 0.262242 10.991525
VNQ 0.033581 1.147283 -4.339915 5.965661 0.079166 3.169672
VTI 0.074964 1.049094 -5.868308 10.145648 0.937027 19.892745
Запуск MCMC семплирования...
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 0 0.35 15 60.17 draws/s 0:00:24 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 1500 0 0.28 15 31.56 draws/s 0:00:47 0:00:00
Апостериорные оценки ожидаемых доходностей:
Апост. ср. доходность (%) Апост. ст. откл. (%) Эмпирическая доходность (%)
Актив
BND 0.018025 0.016615 0.018201
GLD 0.112731 0.043598 0.112249
VEA 0.059264 0.041446 0.059395
VNQ 0.034087 0.050416 0.033581
VTI 0.074750 0.046479 0.074964
Рис. 21: Байесовский анализ параметров рынка. а) Апостериорные распределения ожидаемых доходностей с полной квантификацией неопределенности; б) Апостериорная корреляционная матрица, усредненная по всем MCMC выборкам, отражающая байесовские оценки взаимосвязей между активами
Оптимизация портфеля...
Оптимальные веса портфеля:
Оптимальный вес (%) Ст. откл. веса (%)
Актив
BND 27.562952 20.593457
GLD 33.868963 15.395348
VEA 4.006936 10.125138
VNQ 2.001456 6.020346
VTI 32.559693 15.955165
Метрики производительности оптимального портфеля:
Значение
Метрика
Годовая доходность (%) 18.908235
Годовая волатильность (%) 9.409020
Коэффициент Шарпа 2.009586
Максимальная просадка (%) -7.271414
=== Быстрая диагностика MCMC ===
Средний эффективный размер выборки: 2070
Максимальный R-hat: 1.006
✓ Сходимость хорошая (R-hat < 1.1)
Рис. 22: Полная байесовская оптимизация портфеля с использованием всех MCMC выборок. а) Оптимальные веса с планками неопределенности, отражающими всю апостериорную информацию; б) Распределение соотношения риск-доходность портфеля; в) Апостериорные распределения весов активов; г) Байесовская эффективная граница с отдельными активами и выбранным портфелем
Рис. 23: Историческая производительность оптимального байесовского портфеля. Оптимальный портфель (черная линия) демонстрирует сбалансированную динамику роста, сочетая высокую доходность акций с защитными свойствами облигаций и альтернативных активов
Рис. 24: Ежемесячная доходность оптимального байесовского портфеля за весь период наблюдения. Зеленые столбцы отражают месяцы с положительной доходностью, красные — с отрицательной. Портфель демонстрирует преобладание положительных месяцев и контролируемую величину просадок
Рис. 25: Временная динамика просадок оптимального байесовского портфеля. Красная заливка показывает периоды снижения стоимости портфеля относительно исторических максимумов
Рис. 26: Ключевые метрики производительности оптимального байесовского портфеля
Рис. 27: Диагностика качества MCMC семплирования для байесовской оптимизации портфеля. Локальный эффективный размер выборки подтверждает высокое качество и независимость выборок для надежных статистических выводов
Этот код демонстрирует полный процесс байесовской оптимизации портфеля с использованием методов MCMC:
- Загрузка и анализ данных — загружаем исторические цены активов, вычисляем доходности и анализируем их статистические свойства;
- Байесовская оценка параметров — используем MCMC для оценки апостериорных распределений ожидаемых доходностей, волатильностей и корреляций;
- Оптимизация портфеля — находим оптимальные веса активов, учитывая неопределенность в оценках параметров;
- Анализ результатов — визуализируем оптимальные веса, их распределения, эффективную границу и оцениваем историческую производительность оптимального портфеля.
Ключевые преимущества байесовского подхода к оптимизации портфеля:
- Учет неопределенности в оценках параметров — вместо точечных оценок мы работаем с распределениями, что позволяет более реалистично оценивать риск;
- Естественная регуляризация через априорные распределения — предотвращает экстремальные аллокации, характерные для классической оптимизации Марковица;
- Гибкость в моделировании — можно использовать более сложные модели доходностей, включая асимметрию, тяжелые хвосты и зависимости, меняющиеся во времени;
- Возможность включения экспертных мнений — априорные распределения могут отражать экспертные прогнозы и убеждения;
- Прямая оценка неопределенности в оптимальных весах — мы получаем не только точечные оценки оптимальных весов, но и их распределения.
Байесовская оптимизация портфеля с использованием MCMC представляет собой мощный инструмент для управления инвестициями, который естественным образом учитывает неопределенность в оценках параметров и позволяет принимать более обоснованные инвестиционные решения.
Моделирование временных рядов с MCMC
Временные ряды представляют собой еще одну важную область, где методы MCMC могут быть особенно полезны. Байесовский подход к моделированию временных рядов позволяет естественным образом учитывать структурные изменения, нелинейности и сложные зависимости в данных.
Рассмотрим применение MCMC для байесовской структурной модели временного ряда, которая позволяет выделить тренд, сезонность и шум:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pymc as pm
import arviz as az
from scipy import stats
import yfinance as yf
from statsmodels.tsa.seasonal import seasonal_decompose
import pytensor.tensor as pt
from numba import jit
import warnings
warnings.filterwarnings('ignore')
# Настройка JAX backend для PyMC (если доступен)
try:
import jax
os.environ["XLA_FLAGS"] = "--xla_force_host_platform_device_count=4"
print("JAX backend доступен для ускорения")
except ImportError:
print("JAX недоступен, используем стандартный backend")
# Загрузка недельных данных временного ряда
def load_time_series_data(ticker='SPY', start_date='2023-06-19', end_date='2025-06-29'):
data = yf.download(ticker, start=start_date, end=end_date, interval='1wk')
return data
# Загружаем недельные данные S&P 500 ETF
spy_data = load_time_series_data()
# Преобразование данных для моделирования
log_price = np.log(spy_data['Close'])
# Визуализация временного ряда
plt.figure(figsize=(12, 6))
log_price.plot(color='#333333')
plt.title('Логарифм цены закрытия SPY (недельные данные)')
plt.ylabel('Log(Цена)')
plt.grid(True, alpha=0.3)
plt.savefig('log_price_series.png', dpi=300)
plt.show()
# Классическая декомпозиция временного ряда для сравнения
decomposition = seasonal_decompose(log_price, model='additive', period=52) # 52 недели в году
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)
decomposition.observed.plot(ax=axes[0], color='#333333')
axes[0].set_title('Наблюдаемый ряд')
axes[0].grid(True, alpha=0.3)
decomposition.trend.plot(ax=axes[1], color='#333333')
axes[1].set_title('Тренд')
axes[1].grid(True, alpha=0.3)
decomposition.seasonal.plot(ax=axes[2], color='#333333')
axes[2].set_title('Сезонность')
axes[2].grid(True, alpha=0.3)
decomposition.resid.plot(ax=axes[3], color='#333333')
axes[3].set_title('Остатки')
axes[3].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('classical_decomposition.png', dpi=300)
plt.show()
# Байесовская структурная модель временного ряда
def bayesian_structural_time_series(data, n_samples=2000, tune=1000):
y = data.values
T = len(y)
t = np.arange(T)
t_norm = t / T
# Параметры модели
n_harmonics = 2
n_knots = 6
knot_locations = np.linspace(0, 1, n_knots)
# Предварительно вычисляем базисные функции
def compute_bspline_basis(x, knots):
n = len(knots)
basis = np.zeros((len(x), n))
for i, k in enumerate(knots):
basis[:, i] = np.exp(-50 * (x - k)**2)
return basis
# Вычисляем базис
basis_matrix = compute_bspline_basis(t_norm, knot_locations)
# Предварительно вычисляем тригонометрические функции
sin_components = np.zeros((T, n_harmonics))
cos_components = np.zeros((T, n_harmonics))
for h in range(n_harmonics):
freq = 2 * np.pi * (h + 1) / 52 # 52 недели в году
sin_components[:, h] = np.sin(freq * t)
cos_components[:, h] = np.cos(freq * t)
with pm.Model() as model:
# Компонент тренда
spline_coef = pm.Normal('spline_coef', mu=0, sigma=1, shape=n_knots)
trend = pm.Deterministic('trend', pt.dot(basis_matrix, spline_coef))
# Сезонные компоненты
sin_coef = pm.Normal('sin_coef', mu=0, sigma=0.1, shape=n_harmonics)
cos_coef = pm.Normal('cos_coef', mu=0, sigma=0.1, shape=n_harmonics)
# Векторизованное вычисление сезонности
seasonal = pm.Deterministic('seasonal',
pt.dot(sin_components, sin_coef) + pt.dot(cos_components, cos_coef))
# Шумовая компонента
sigma = pm.HalfNormal('sigma', sigma=0.1)
rho = pm.Uniform('rho', 0, 1)
# Векторизованная AR(1) модель
mu = trend + seasonal
# Используем GaussianRandomWalk для AR(1) процесса
innovations = pm.Normal('innovations', mu=0, sigma=sigma, shape=T)
# Векторизованная AR(1) структура
def ar1_logp(value, mu, rho, sigma, innovations):
# Первое наблюдение
logp = pm.Normal.logp(value[0], mu[0], sigma / pt.sqrt(1 - rho**2))
# Остальные наблюдения
for i in range(1, T):
logp += pm.Normal.logp(value[i], mu[i] + rho * (value[i-1] - mu[i-1]), sigma)
return logp
# Векторизованная модель наблюдений
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
# Настройки семплирования с попыткой использования JAX
try:
# Пробуем nutpie sampler (JAX-based)
trace = pm.sample(n_samples, tune=tune,
chains=4,
cores=4,
nuts_sampler="nutpie",
target_accept=0.9,
return_inferencedata=True)
except:
# Fallback на стандартный NUTS
trace = pm.sample(n_samples, tune=tune,
chains=4,
cores=4,
target_accept=0.9,
max_treedepth=12,
return_inferencedata=True)
return trace, model
# Запуск байесовской структурной модели
print("Запуск байесовской структурной модели...")
trace, model = bayesian_structural_time_series(log_price)
# Анализ результатов
def analyze_structural_model(trace, data):
# Извлечение компонент с правильной обработкой размерностей
trend_samples = trace.posterior['trend'].values
seasonal_samples = trace.posterior['seasonal'].values
# Извлекаем chains, draws, time_points
n_chains, n_draws, n_timepoints = trend_samples.shape
# Объединяем все цепи и выборки
trend = trend_samples.reshape(n_chains * n_draws, n_timepoints)
seasonal = seasonal_samples.reshape(n_chains * n_draws, n_timepoints)
# Средние значения компонент
mean_trend = np.mean(trend, axis=0)
mean_seasonal = np.mean(seasonal, axis=0)
# Доверительные интервалы
trend_lower = np.percentile(trend, 5, axis=0)
trend_upper = np.percentile(trend, 95, axis=0)
seasonal_lower = np.percentile(seasonal, 5, axis=0)
seasonal_upper = np.percentile(seasonal, 95, axis=0)
# Вычисление остатков
residuals = data.values - (mean_trend + mean_seasonal)
# График 1: Наблюдаемый ряд и модель
plt.figure(figsize=(12, 6))
data.plot(color='#333333', alpha=0.7, label='Наблюдения')
plt.plot(data.index, mean_trend + mean_seasonal, color='blue',
label='Модель (тренд + сезонность)', linewidth=2)
plt.fill_between(data.index,
mean_trend + mean_seasonal + 2*np.std(residuals),
mean_trend + mean_seasonal - 2*np.std(residuals),
color='blue', alpha=0.2, label='±2σ доверительная полоса')
plt.title('Наблюдаемый ряд и байесовская структурная модель')
plt.ylabel('Log(Цена)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('observed_vs_model.png', dpi=300)
plt.show()
# График 2: Тренд
plt.figure(figsize=(12, 6))
plt.plot(data.index, mean_trend, color='#333333', label='Средний тренд', linewidth=2)
plt.fill_between(data.index, trend_lower, trend_upper, color='#333333', alpha=0.2,
label='90% доверительный интервал')
plt.title('Компонента тренда')
plt.ylabel('Тренд')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('trend_component.png', dpi=300)
plt.show()
# График 3: Сезонность
plt.figure(figsize=(12, 6))
plt.plot(data.index, mean_seasonal, color='#333333', label='Средняя сезонность', linewidth=2)
plt.fill_between(data.index, seasonal_lower, seasonal_upper, color='#333333', alpha=0.2,
label='90% доверительный интервал')
plt.title('Сезонная компонента')
plt.ylabel('Сезонность')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('seasonal_component.png', dpi=300)
plt.show()
# График 4: Остатки
plt.figure(figsize=(12, 6))
plt.plot(data.index, residuals, color='#333333', alpha=0.7)
plt.axhline(y=0, color='red', linestyle='--', linewidth=2)
plt.title('Остатки модели (наблюдения - тренд - сезонность)')
plt.ylabel('Остатки')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('residuals.png', dpi=300)
plt.show()
# График 5: Автокорреляция остатков
plt.figure(figsize=(10, 6))
# Векторизованное вычисление автокорреляции
max_lag = 15
autocorr_values = []
for lag in range(1, max_lag + 1):
if lag < len(residuals): corr = np.corrcoef(residuals[:-lag], residuals[lag:])[0, 1] autocorr_values.append(corr) else: autocorr_values.append(0) plt.bar(range(1, len(autocorr_values) + 1), autocorr_values, color='#333333', alpha=0.7) plt.axhline(y=0, color='red', linestyle='-', alpha=0.5) plt.title('Автокорреляция остатков') plt.xlabel('Лаг (недели)') plt.ylabel('Автокорреляция') plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('residuals_autocorrelation.png', dpi=300) plt.show() return { 'mean_trend': mean_trend, 'mean_seasonal': mean_seasonal, 'residuals': residuals, 'trend_ci': (trend_lower, trend_upper), 'seasonal_ci': (seasonal_lower, seasonal_upper) } # Анализ результатов структурной модели results = analyze_structural_model(trace, log_price) # Анализ распределения параметров def analyze_parameters(trace): rho_samples = trace.posterior['rho'].values sigma_samples = trace.posterior['sigma'].values sin_coef_samples = trace.posterior['sin_coef'].values cos_coef_samples = trace.posterior['cos_coef'].values # Объединяем цепи и выборки rho = rho_samples.flatten() sigma = sigma_samples.flatten() sin_coef = sin_coef_samples.reshape(-1, sin_coef_samples.shape[-1]) cos_coef = cos_coef_samples.reshape(-1, cos_coef_samples.shape[-1]) # График 1: Параметр автокорреляции plt.figure(figsize=(10, 6)) sns.histplot(rho, kde=True, color='#333333', alpha=0.7) plt.axvline(np.mean(rho), color='red', linestyle='--', label=f'Среднее: {np.mean(rho):.3f}') plt.title('Апостериорное распределение параметра автокорреляции (ρ)') plt.xlabel('ρ') plt.ylabel('Плотность') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('rho_distribution.png', dpi=300) plt.show() # График 2: Параметр шума plt.figure(figsize=(10, 6)) sns.histplot(sigma, kde=True, color='#333333', alpha=0.7) plt.axvline(np.mean(sigma), color='red', linestyle='--', label=f'Среднее: {np.mean(sigma):.4f}') plt.title('Апостериорное распределение стандартного отклонения шума (σ)') plt.xlabel('σ') plt.ylabel('Плотность') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('sigma_distribution.png', dpi=300) plt.show() # График 3: Коэффициенты синусов plt.figure(figsize=(10, 6)) for i in range(sin_coef.shape[1]): sns.kdeplot(sin_coef[:, i], label=f'Гармоника {i+1}', alpha=0.7) plt.title('Апостериорные распределения коэффициентов синусов') plt.xlabel('Коэффициент') plt.ylabel('Плотность') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('sin_coefficients.png', dpi=300) plt.show() # График 4: Коэффициенты косинусов plt.figure(figsize=(10, 6)) for i in range(cos_coef.shape[1]): sns.kdeplot(cos_coef[:, i], label=f'Гармоника {i+1}', alpha=0.7) plt.title('Апостериорные распределения коэффициентов косинусов') plt.xlabel('Коэффициент') plt.ylabel('Плотность') plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('cos_coefficients.png', dpi=300) plt.show() # Анализ параметров analyze_parameters(trace) # Векторизованное прогнозирование @jit(nopython=True) def generate_forecast_sample(trend_last, trend_slope, seasonal_values, rho_val, sigma_val, n_forecast): """Генерация одной выборки прогноза""" forecast = np.zeros(n_forecast) last_error = 0.0 for t in range(n_forecast): # Линейное продолжение тренда trend_forecast = trend_last + trend_slope * (t + 1) # Повторение сезонности seasonal_forecast = seasonal_values[t % len(seasonal_values)] # AR(1) процесс для ошибок error_forecast = rho_val * last_error + np.random.normal(0, sigma_val) forecast[t] = trend_forecast + seasonal_forecast + error_forecast last_error = error_forecast return forecast def forecast_time_series(trace, data, n_forecast=10): # Разделение данных train_data = data[:-n_forecast] test_data = data[-n_forecast:] print("Обучение модели на тренировочных данных...") trace_train, model_train = bayesian_structural_time_series(train_data) # Извлечение компонент с обработкой размерностей trend_samples = trace_train.posterior['trend'].values seasonal_samples = trace_train.posterior['seasonal'].values rho_samples = trace_train.posterior['rho'].values sigma_samples = trace_train.posterior['sigma'].values # Объединяем цепи и выборки trend_train = trend_samples.reshape(-1, trend_samples.shape[-1]) seasonal_train = seasonal_samples.reshape(-1, seasonal_samples.shape[-1]) rho_train = rho_samples.flatten() sigma_train = sigma_samples.flatten() n_samples = len(rho_train) # Параллельное прогнозирование print("Генерация прогнозов...") forecasts = np.zeros((n_samples, n_forecast)) # Векторизованное вычисление для всех выборок for i in range(n_samples): trend_last = trend_train[i, -1] trend_slope = (trend_train[i, -1] - trend_train[i, -2]) if len(trend_train[i]) > 1 else 0
seasonal_values = seasonal_train[i, -52:] if len(seasonal_train[i]) >= 52 else seasonal_train[i]
forecast = generate_forecast_sample(
trend_last, trend_slope, seasonal_values,
rho_train[i], sigma_train[i], n_forecast
)
forecasts[i] = forecast
# Статистики прогноза
mean_forecast = np.mean(forecasts, axis=0)
lower_forecast = np.percentile(forecasts, 5, axis=0)
upper_forecast = np.percentile(forecasts, 95, axis=0)
# Визуализация прогноза
plt.figure(figsize=(12, 6))
plt.plot(train_data.index, train_data.values, color='#333333',
label='Обучающие данные', linewidth=1)
plt.plot(test_data.index, test_data.values, color='blue',
label='Фактические данные', linewidth=2)
plt.plot(test_data.index, mean_forecast, color='red', linestyle='--',
label='Прогноз', linewidth=2)
plt.fill_between(test_data.index, lower_forecast, upper_forecast,
color='red', alpha=0.2, label='90% доверительный интервал')
plt.title('Байесовский прогноз временного ряда SPY')
plt.ylabel('Log(Цена)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bayesian_forecast.png', dpi=300)
plt.show()
# Метрики качества
mse = np.mean((test_data.values - mean_forecast)**2)
mae = np.mean(np.abs(test_data.values - mean_forecast))
coverage = np.mean((test_data.values >= lower_forecast) &
(test_data.values <= upper_forecast))
print(f"Среднеквадратическая ошибка (MSE): {mse:.6f}")
print(f"Средняя абсолютная ошибка (MAE): {mae:.6f}")
print(f"Покрытие 90% доверительного интервала: {coverage:.2f}")
return {
'mean_forecast': mean_forecast,
'lower_forecast': lower_forecast,
'upper_forecast': upper_forecast,
'mse': mse,
'mae': mae,
'coverage': coverage
}
# Прогнозирование
forecast_results = forecast_time_series(trace, log_price)
# Диагностика MCMC
def mcmc_diagnostics(trace):
print("=== Диагностика MCMC ===")
# R-hat статистики
rhat = az.rhat(trace)
print("R-hat статистики:")
for var in rhat.data_vars:
max_rhat = rhat[var].values.max()
print(f" {var}: {max_rhat:.4f}")
# Эффективный размер выборки
ess = az.ess(trace)
print("\nЭффективный размер выборки:")
for var in ess.data_vars:
min_ess = ess[var].values.min()
print(f" {var}: {min_ess:.0f}")
# Трассировки ключевых параметров
az.plot_trace(trace, var_names=['rho', 'sigma'])
plt.tight_layout()
plt.savefig('trace_diagnostics.png', dpi=300)
plt.show()
# Автокорреляция параметров
az.plot_autocorr(trace, var_names=['rho', 'sigma'], max_lag=50)
plt.tight_layout()
plt.savefig('parameter_autocorr.png', dpi=300)
plt.show()
# Диагностика
mcmc_diagnostics(trace)
Рис. 28: Временной ряд логарифма цены закрытия S&P 500 ETF (SPY) за период 2022-2025. Наблюдается общий восходящий тренд с периодами экстремальной и спокойной волатильности и возможными сезонными паттернами
Рис. 29: Классическая аддитивная декомпозиция временного ряда SPY с годовой сезонностью (52 недели). а) Исходный ряд; б) Выделенный тренд показывает долгосрочную динамику; в) Сезонная компонента с годовым циклом; г) Остатки после удаления тренда и сезонности
Запуск байесовской структурной модели...
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.18 31 32.99 draws/s 0:01:30 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.24 31 39.54 draws/s 0:01:15 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.20 31 35.61 draws/s 0:01:24 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.22 15 36.92 draws/s 0:01:21 0:00:00 
Рис. 30: Апостериорное распределение параметра автокорреляции (ρ) в модели остатков. Распределение показывает степень оставшейся временной зависимости после учета структурных компонент, с красной линией, обозначающей среднее значение
Рис. 31: Апостериорное распределение стандартного отклонения шума (σ). Параметр характеризует уровень непредсказуемой еженедельной волатильности в логарифмах цен SPY после учета всех детерминистических компонент модели
Рис. 32: Апостериорные распределения коэффициентов синусов для двух сезонных гармоник. Первая гармоника отражает основной годовой цикл (52 недели), вторая — полугодовой цикл (26 недель), захватывающие различные аспекты сезонности фондового рынка
Рис. 33: Апостериорные распределения коэффициентов косинусов для сезонных гармоник. Совместно с синусными коэффициентами формируют полное описание фазы и амплитуды внутригодовых циклов в движении цен SPY
Обучение модели на тренировочных данных...
Progress Draws Divergences Step size Grad evals Sampling Speed Elapsed Remaining
─────────────────────────────────────────────────────────────────────────────────────────────────────────────────
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.13 31 31.38 draws/s 0:01:35 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.10 31 31.04 draws/s 0:01:36 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.12 31 29.97 draws/s 0:01:40 0:00:00
━━━━━━━━━━━━━━━━━━━━━━━━━ 3000 0 0.12 31 30.78 draws/s 0:01:37 0:00:00
Генерация прогнозов...
Среднеквадратическая ошибка (MSE): 0.517752
Средняя абсолютная ошибка (MAE): 0.678277
=== Диагностика MCMC ===
R-hat статистики:
cos_coef: 1.0005
innovations: 1.0033
rho: 1.0015
seasonal: 1.0011
sigma: 1.0005
sin_coef: 1.0009
spline_coef: 1.0006
trend: 1.0009
Эффективный размер выборки:
cos_coef: 4200
innovations: 9683
rho: 10869
seasonal: 3331
sigma: 12661
sin_coef: 3217
spline_coef: 3267
trend: 3339
Рис. 34: Прогноз недельных цен SPY с использованием байесовской структурной модели на 10 недель вперед. Красная пунктирная линия — средний прогноз, красная область — 90% доверительный интервал. Синяя линия показывает фактические значения для оценки качества прогнозирования
Рис. 35: Диагностика сходимости MCMC для ключевых параметров (ρ и σ). Четыре стационарные цепи без трендов подтверждают успешную сходимость алгоритма к стационарному распределению
Этот код демонстрирует применение MCMC для байесовской структурной модели временного ряда, которая разделяет ряд на компоненты тренда, сезонности и шума. Основные этапы анализа:
- Загрузка и визуализация данных — работаем с логарифмом цены закрытия индексного ETF S&P 500 (SPY);
- Классическая декомпозиция — для сравнения проводим стандартную декомпозицию временного ряда;
- Байесовская структурная модель — создаем модель, которая включает нелинейный тренд, моделируемый сплайном, сезонную компоненту, представленную суммой гармоник, и шумовую компоненту с автокорреляцией;
- Анализ результатов — визуализируем выделенные компоненты, их доверительные интервалы и распределения параметров модели;
- Прогнозирование — используем модель для прогнозирования будущих значений ряда с оценкой неопределенности.
Несмотря на то, что Байесовский подход к моделированию временных рядов может ошибаться, он все же имеет ряд важных преимуществ:
- Естественная количественная оценка неопределенности — мы получаем не только точечные оценки компонент, но и их полные распределения;
- Гибкость в спецификации модели — можно легко включать дополнительные компоненты, такие как структурные изменения, нелинейные зависимости, или внешние факторы;
- Устойчивость к выбросам — с правильно выбранными распределениями шума модель становится устойчивой к аномальным наблюдениям;
- Естественная обработка пропущенных данных — пропуски можно рассматривать как латентные переменные и моделировать вместе с остальными параметрами;
- Включение априорной информации — можно использовать экспертные знания через априорные распределения параметров;
- Надежные доверительные интервалы для прогнозов — учитывается неопределенность как в параметрах модели, так и в будущих реализациях шума.
Структурные модели временных рядов с MCMC представляют собой мощный инструмент для анализа и прогнозирования временных рядов, который может быть применен к широкому спектру задач, от финансовых рынков до климатических данных и экономических показателей.
Выводы
На протяжении этой статьи мы подробно рассмотрели методы Марковских цепей Монте-Карло (MCMC), начиная с теоретических основ и заканчивая практическими применениями в анализе данных и статистическом моделировании. Теперь, подводя итоги, можно выделить ключевые преимущества MCMC, которые делают их незаменимым инструментом в современном байесовском анализе:
- Универсальность — MCMC применимы к широкому спектру вероятностных моделей, от простых до крайне сложных, без необходимости аналитического решения интегралов;
- Байесовская интерпретация — методы MCMC естественным образом вписываются в байесовскую парадигму, позволяя работать с полными апостериорными распределениями вместо точечных оценок;
- Оценка неопределенности — MCMC предоставляют надежные оценки неопределенности для всех параметров модели и производных величин;
- Гибкость в моделировании — возможность создавать и исследовать сложные иерархические модели, которые лучше отражают реальные процессы;
- Естественная обработка латентных переменных — MCMC позволяют эффективно моделировать ненаблюдаемые величины вместе с параметрами модели;
- Устойчивость к переобучению — байесовский подход с MCMC обеспечивает автоматическую регуляризацию через априорные распределения;
- Адаптивность — современные алгоритмы MCMC, такие как HMC и NUTS, адаптируются к геометрии целевого распределения, обеспечивая эффективное исследование пространства параметров;
- Теоретическая обоснованность — методы MCMC имеют строгое математическое обоснование, гарантирующее их асимптотическую корректность.
Мы рассмотрели три основных алгоритма MCMC — Метрополис-Гастингс, семплирование по Гиббсу и Hamiltonian Monte Carlo, каждый из которых имеет свои преимущества и области применения. Практические примеры показали, как эти методы могут быть успешно применены к различным задачам, включая байесовскую логистическую регрессию, оптимизацию инвестиционного портфеля и структурное моделирование временных рядов. Во всех этих случаях MCMC обеспечили не только точечные оценки параметров, но и полное представление о их распределениях и взаимосвязях.
Методы MCMC являются одним из основных инструментов в анализе финансовых данных, предлагая уникальные возможности для работы со сложными вероятностными моделями. Их способность обеспечивать полное представление о распределениях параметров и взаимосвязях делает их незаменимыми в современной статистике и машинном обучении.