Выборочные методы MCMC — Алгоритм Метрополиса-Гастингса

Введение

Методы Монте-Карло по схеме марковских цепей (MCMC) — это класс выборочных алгоритмов, позволяющих извлекать выборки из сложного целевого распределения p(x). Цель MCMC — получить набор выборок, которые точно представляют это целевое распределение, чтобы можно было оценить средние значения, дисперсии и другие характеристики этого распределения. Суть MCMC заключается в построении специальной цепи Маркова, стационарное распределение которой совпадает с целевым распределением.  

Эти методы нашли широкое применение в байесовском выводе, машинном обучении и других областях, где требуется аппроксимация апостериорных распределений. В отличие от детерминированных методов, таких как вариационный вывод или аппроксимация Лапласа, MCMC-алгоритмы обладают уникальным свойством: при правильной настройке они в пределе гарантируют получение выборок, точно соответствующих целевому распределению. Среди множества алгоритмов MCMC особое место занимает алгоритм Метрополиса-Гастингса (МГ) — фундаментальный метод, лежащий в основе многих современных подходов.

В данной статье мы рассмотрим алгоритм Метрополиса-Гастингса, начиная с его теоретических основ и ключевых концепций. Затем мы представим реализацию алгоритма на MQL5 в виде класса MHSampler, а также разберем простые примеры его применения для одномерных и многомерных распределений, включая варианты Random Walk и Independent МГ. 

Нужно отметить, что при использовании методов MCMC особое внимание следует уделять диагностике и настройке алгоритма, поскольку его корректная работа требует тщательного анализа сходимости цепи с использованием статистических и визуальных инструментов.

Алгоритм Метрополиса

Базовый алгоритм Метрополиса (Metropolis, 1953) генерирует выборки из целевого распределения, предлагая на каждом шаге переход из текущего состояния x в новое состояние x′ с вероятностью q(x′∣x), где q называется вспомогательным (предлагающим) распределением (proposal distribution). Главное требование к вспомогательному распределению — оно должно гарантировать, что цепь не застрянет в одной области, то есть должна существовать возможность перехода из любой точки пространства, которое занимает целевое распределение, в любую другую точку. Для этого достаточно, чтобы вероятность q(x′∣x) была ненулевой для всех возможных состояний x и x′.

Предложенное состояние x' принимается или отклоняется по критерию, обеспечивающему, чтобы цепь в долгосрочной перспективе посещала состояния пропорционально p(x). Если состояние x' принято, оно становится новым состоянием, в противном случае цепь остается в текущем состоянии (происходит дублирование семпла).

Алгоритм Метрополиса является частным случаем Метрополиса-Гастингса, где вспомогательное распределение q(x'|x) симметрично, то есть вероятность предложения перехода из x в x′ равна вероятности обратного перехода: q(x′∣x)=q(x∣x′).

Благодаря этой симметрии, критерий принятия кандидата x′ значительно упрощается:

A = min (1, p*(x')/p*(x))

где p*(x) — ненормированная плотность целевого распределения, а само распределение задается как p(x)=p*(x)/Zp.

Нормирующая константа Zp может быть неизвестна, что делает алгоритм особенно удобным для задач, где аналитическое вычисление Zp затруднено, как это часто бывает в байесовском выводе.

Если плотность в новой точке x′ выше, чем в текущей точке x (т.е., p*(x′) > p*(x)), то коэффициент принятия A=1, и переход в x′ происходит всегда (с вероятностью 1). Это позволяет цепи быстро взбираться к пикам распределения, обеспечивая поиск областей высокой вероятности.

Если плотность снижается (p*(x′) < p*(x)), переход в x′ происходит с вероятностью A = p*(x′)/p*(x). Этот механизм позволяет цепи спускаться с максимумов распределения и исследовать склоны и долины, распределения обеспечивая эффективный поиск всего пространства, что не позволяет цепи навсегда застрять в одной точке.

Именно этот баланс между безусловным принятием "вверх" и вероятностным принятием "вниз" гарантирует, что цепь сохраняет пропорции целевого распределения, проводя больше времени там, где плотность выше.

Полученная таким образом последовательность выборок образует марковскую цепь, где каждое следующее состояние зависит от предыдущего. Это приводит к автокорреляции между выборками, что может снижать их репрезентативность (из-за высокой автокорреляции соседние выборки содержат почти одинаковую информацию).

Например, если вы собрали 1000 сильно коррелированных выборок, они могут содержать ту же информацию, что и всего 50 независимых выборок. Ваши 1000 выборок не представляют 1000 уникальных наблюдений о распределении, а дают лишь 50 эффективных наблюдений. Как следствие оценки среднего, дисперсии и других характеристик распределения, основанные на этих коррелированных выборках, будут иметь более высокую дисперсию (ошибку), чем если бы вы использовали то же количество независимых выборок. Поэтому для уменьшения корреляции применяют так называемый период прогрева (burn-in), отбрасывающий начальные итерации, и прореживание (thinning), сохраняющее только каждую i-ю выборку.

Рис.1 Алгоритм Метрополиса

Марковские цепи

Методы Монте-Карло по схеме марковских цепей  используют марковские цепи для генерации выборок из сложных распределений. Чтобы понять, как это работает, разберем два ключевых свойства марковских цепей, которые должны выполняться для построения MCMC-алгоритма: их "память" и способность находить нужное распределение.

1. Марковское свойство

Марковская цепь — это последовательность состояний x(1), x(2), … где каждое следующее состояние x(t+1) зависит только от текущего x(t), и не зависит от всей предыдущей истории цепи x(1), …, x(t−1):

P(x(t+1) |x(t), …, x(1)) = P(x(t+1) |x(t))

Это свойство означает, что для прогнозирования будущего состояния достаточно знать только настоящее.

2. Сходимость к целевому распределению

Цель алгоритмов MCMC — построить марковскую цепь, которая со временем (после начального периода прогрева, или burn-in) порождает выборки из целевого распределения p(x). Это достигается за счет специального правила переходов, которое удовлетворяет условию детального равновесия (Detailed Balance):

p(x)T(x′∣x) = p(x′)T(x∣x′)

где:

• p(x) и p(x′) — плотности целевого распределения,
• T(x′∣x) — полная вероятность перехода из x в x'.

В алгоритме Метрополиса-Гастингса это важное условие обеспечивается выбором подходящего коэффициента принятия A, который гарантирует, что цепь со временем "посещает" состояния пропорционально p(x), даже если само распределение задано лишь частично (без нормирующей константы Zp). При этом полная вероятность перехода T(x′∣x) строится как произведение вероятности предложения q(x′∣x) и вероятности принятия A(x′∣x).

Алгоритм Метрополиса-Гастингса

Алгоритм Метрополиса-Гастингса (Hastings, 1970), обобщает алгоритм Метрополиса, позволяя использовать несимметричные предлагающие распределения, то есть q(x'|x) ≠ q(x|x'). Это делает алгоритм более гибким и применимым к широкому классу задач, где симметричные распределения не подходят (например, при использовании логнормальных или экспоненциальных предложений).

Для компенсации асимметрии в критерий принятия A, вводится поправка Гастингса, которая корректирует смещение, создаваемое несимметричным распределением. Вероятность принятия кандидата x' определятся как:

A = min (1, p*(x')q(x|x') / p*(x)q(x'|x))

где:

• p*(x) — ненормированная плотность целевого распределения,
• q(x|x') / q (x'|x) — поправка Гастингса, учитывающая асимметрию предлагающего распределения.

Для вычисления A, как и прежде, не требуется знание нормирующей константы Zp в распределении вероятности р(x) = р*(x)/Zp, поскольку она сокращается.

Влияние вспомогательного распределения на производительность

Выбор  q(x'|x) существенно влияет на эффективность алгоритма. В непрерывных пространствах часто используется нормальное распределение, центрированное в текущем состоянии x(t), что приводит к поведению, известному как Random Walk Metropolis-Hastings (RWMH). При этом ключевую роль играет параметр дисперсии (или шага цепи), который определяет масштаб предлагаемых переходов:

• Малая дисперсия: доля принятых переходов(acceptance rate)  будет высокой, но цепь будет двигаться медленно. Это приведет к длительному времени автокорреляции выборок, поскольку цепь медленно исследует пространство.
• Большая дисперсия: доля принятых переходов будет низкой, поскольку большинство предложенных шагов будут попадать в области низкой вероятности. Это приводит к частым отклонениям и неэффективному использованию вычислений.

На практике оптимальная производительность достигается, когда доля принятия находится в диапазоне 20–40%. Это требует тщательной настройки параметров q, чтобы сбалансировать скорость исследования пространства и вероятность принятия кандидатов.

Рис.2 Алгоритм Метрополиса-Гастингса

Алгоритм Метрополиса-Гастингса реализован в MQL5 с использованием объектно-ориентированного подхода. Основной класс MHSampler требует реализации трех абстрактных классов: LogPDF(логарифм целевой плотности ), PropRND (генератор нового состояния)  и LogPropPDF(логарифм плотности перехода). Благодаря этому можно легко настроить алгоритм для любого целевого и вспомогательного распределения. 

Класс LogPDF

Класс LogPDF определяет интерфейс для вычисления логарифма плотности вероятности целевого распределения p*(x). Метод  LogPdf(const vector &x),  принимает вектор x (точку в пространстве состояний) и возвращает lnp*(x). 

Использование логарифма плотности — стандартная практика в MCMC, так как это предотвращает численные ошибки при работе с малыми вероятностями, особенно в многомерных задачах. Класс может быть реализован для любого целевого распределения, например, нормального или смеси распределений. В байесовском выводе этот класс реализуется для вычисления логарифма ненормированной апостериорной плотности, что выражается суммой логарифмов: lnp*(x) = ln(Правдоподобие)+ln(Априорное распределение).

Класс PropRND

Класс PropRND отвечает за генерацию случайных состояний-кандидатов x′ из вспомогательного распределения q(x′∣x). Метод PropRnd(const vector &x), принимает текущее состояние x и возвращает новый вектор x′, сгенерированный из распределения q(x′∣x).

Например, для Random Walk MH пользователь может реализовать PropRND с гауссовым распределением, центрированным в текущем состоянии x, а для Independent MH — с фиксированным распределением, не зависящим от x.

Класс LogPropPDF

Класс LogPropPDF отвечает за вычисление логарифма плотности предлагающего  распределения q(x'∣x). Этот класс не генерирует новых состояний, а только оценивает плотность вероятности перехода.  Метод LogPropPdf(const vector &x_to, const vector &x_from) возвращает lnq(x'∣x), где x' — предложенное новое состояние, а x — текущее состояние.

Этот класс используется только для несимметричных предлагающих распределений. В таких случаях ln⁡q(x'∣x) необходим для вычисления поправки Гастингса в критерии принятия. Для симметричных случаев этот класс не используется. 

Класс MHSampler

Метод MHSample выполняет алгоритм Метрополиса или Метрополиса-Гастингса, принимая следующие параметры:

• start — начальное состояние x0(вектор),
• samples — матрица для хранения полученных выборок,
• accept_rate — переменная для хранения доли принятых предложений,
• log_pdf — указатель на объект класса LogPDF для вычисления ln p(x),
• prop_rnd — указатель на объект класса PropRND для генерации кандидатов,
• log_prop_pdf — указатель на объект класса LogPropPDF (или NULL для симметричного случая),
• params — структура MH_Params, содержащая параметры:
• nsamples — количество сохраняемых выборок,
• burnin — период прогрева — количество начальных итераций, которые отбрасываются для достижения стационарного распределения,
• thin — интервал прореживания — сохраняется только каждая i-я выборка для уменьшения корреляции,
• symmetric — флаг, указывающий, является ли вспомогательное распределение симметричным, то есть использовать упрощенную формулу (Метрополис) или полную формулу (МГ).

На каждой итерации, начиная с начального состояния x0, выполняются следующие шаги:

1. Генерация кандидата: x' ∼ q(x'∣x0) с помощью объекта prop_rnd.
2. Критерий принятия (вычисляется логарифм отношения принятия r):

• Для симметричного случая (Метрополис):

r = lnp*(x') − lnp*(x0)

• Для несимметричного случая (Метрополис-Гастингс):

r = [ lnp*(x') + lnq(x0∣x') ] − [ lnp*(x0) + lnq(x'∣x0) ]

     3. Принятие/Отклонение: генерируется равномерное случайное число u ∼ U(0,1). Шаг принимается, если выполняется условие ln u ≤ min (r, 0). Если условие выполнено, текущее состояние обновляется x0 = x', иначе x0 остается неизменным.

     4. Сбор выборок: после периода прогрева burnin и с учетом прореживания (thin) текущее состояние x0 записывается в матрицу samples.

В конце цикла всегда вычисляется доля принятия (accept_rate), которая является основным диагностическим инструментом для оценки эффективности выбранного вспомогательного распределения.

#include <Math\Stat\Uniform.mqh>

//+------------------------------------------------------------------+
//| Абстрактный класс для целевой лог-плотности P(x)                 |
//+------------------------------------------------------------------+
class LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) = 0;
  };

//+------------------------------------------------------------------+
//| Абстрактный класс для предлагающей лог-плотности q(x'|x)         |
//| (Для несимметричного случая)                                     |
//+------------------------------------------------------------------+
class LogPropPDF
  {
public:
   // LogPropPdf(x_to, x_from) возвращает q(x_to | x_from)
   virtual double    LogPropPdf(const vector &x_to, const vector &x_from) = 0;
  };

//+------------------------------------------------------------------+
//| Абстрактный класс для генератора кандидатов q(x'|x)              |
//+------------------------------------------------------------------+
class PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) = 0;
  };

//+------------------------------------------------------------------+
//| Структура параметров для MHSampler                               |
//+------------------------------------------------------------------+
struct MH_Params
  {
   int               nsamples;     // кол-во выборок
   int               burnin;       // период прогрева
   int               thin;         // прореживание
   bool              symmetric;   // Симметричное вспомогательное распределение q(x',x) = q(x,x')
  };

//+------------------------------------------------------------------+
//| Класс MHSampler, реализующий алгоритм Метрополиса-Гастингса      |
//+------------------------------------------------------------------+
class MHSampler
  {
public:
   //+------------------------------------------------------------------+
   //| Основная функция MH-семплирования                                |                                      |
   //+------------------------------------------------------------------+
   bool              MHSample(
      const vector &start,
      matrix &samples,
      double &accept_rate,
      LogPDF *log_pdf,
      PropRND *prop_rnd,
      LogPropPDF *log_prop_pdf, // Может быть NULL, если symmetric = true
      const MH_Params &params
   )
     {
      // --- 1. Проверки и инициализация ---
      if(params.nsamples <= 0 || params.thin <= 0)
        {
         Print("Ошибка: nsamples и thin должны быть > 0");
         return false;
        }
      if(!params.symmetric && log_prop_pdf == NULL)
        {
         Print("Ошибка: Для несимметричного вспомогательного распределения требуется log_prop_pdf.");
         return false;
        }

      int dim = (int)start.Size();
      int total_steps = params.nsamples * params.thin + params.burnin;
      vector x0 = start; // Текущее значение
      double accepted_count = 0;
      int sample_idx = 0;

      // Инициализация матрицы сэмплов
      samples.Resize(params.nsamples, dim);

      // --- 2. Основной Цикл Метрополиса-Гастингса ---
      for(int i = 1 - params.burnin; i <= params.nsamples * params.thin; i++)
        {

         // 2.1. Генерация кандидата: x' ~ q(x'| x0)
         vector x_new = prop_rnd.PropRnd(x0);

         // 2.2. Вычисление лог-вероятностей
         double log_pdf_new = log_pdf.LogPdf(x_new);
         double log_pdf_x0 = log_pdf.LogPdf(x0);

         // 2.3. Вычисление лог-коэффициента принятия A
         double r;
         if(params.symmetric)
           {
            // Random Walk MH (Symmetric q(x,x')=q(x',x)): r = log(P(x')/P(x0))
            r = log_pdf_new - log_pdf_x0;
           }
         else
           {
            // General MH (Asymmetric q): r = log( P(x')q(x0|x') / P(x0)q(x'|x0) )
            // q(x_new | x0) - Вероятность прямого перехода
            double log_prop_x0_to_xnew = log_prop_pdf.LogPropPdf(x_new, x0); // q(x' | x0)
            // q(x0 | x_new) - Вероятность обратного перехода
            double log_prop_xnew_to_x0 = log_prop_pdf.LogPropPdf(x0, x_new); // q(x0| x')

            r = (log_prop_xnew_to_x0 + log_pdf_new) - (log_prop_x0_to_xnew + log_pdf_x0);
           }

         // 2.4. Принятие/Отклонение
         int err;
         double U = MathLog(MathRandomUniform(0.0, 1.0,err));
         //Шаг принимается, если log(U)≤min(0,r), что является эквивалентом
         //классического условия U≤A, где коэффициент принятия: A = min(1, exp(r))
         if(U <= MathMin(0.0, r))
           {
            // Принятие: x0 = x'
            x0 = x_new;
            accepted_count++;
           }
         // Отклонение: x0 остается без изменений

         // 2.5. Сэмплирование (учет burnin и thin)
         if(i > 0 && i % params.thin == 0)
           {
            if(sample_idx < params.nsamples)
              {
               samples.Row(x0, sample_idx);
               sample_idx++;
              }
           }
        }

      // 3. Расчет доли принятия (Acceptance Rate)
      accept_rate = accepted_count / total_steps;

      return true;
     }
  };

Одномерный случай с симметричным вспомогательным распределением (Random Walk Metropolis)

Пример №1

Этот пример демонстрирует реализацию алгоритма Метрополиса со случайным блужданием (Random Walk Metropolis, RWM)— где предлагающее распределение симметрично.

Целевым распределением выбрано одномерное стандартное нормальное распределение N(0,1). Кандидаты генерируются с использованием симметричного вспомогательного распределения:

x' = x + u, где u ∼ Uniform(−δ, δ)

Используются следующие параметры:

• начальная точка x0 = 1,
• период прогрева burnin = 500,
• прореживание thin = 10,
• δ = 4, определяет ширину равномерного распределения U(−δ, + δ).

Если шаг цепи δ слишком большой, доля принятия будет низкой, и цепь будет часто "застревать" из-за отклонения кандидатов. Если шаг цепи δ слишком мал, доля принятия будет высокой, но цепь будет двигаться медленно, что увеличивает корреляцию выборок.

Настройка параметра δ позволяет сбалансировать скорость исследования пространства и вероятность принятия кандидатов. Пользователю необходимо подобрать оптимальное значение дельта δ, чтобы получить независимую (или почти независимую) выборку.

#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

//+------------------------------------------------------------------+
//| Реализация LogPDF для N(0, 1)                                    |
//+------------------------------------------------------------------+
class NormalLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      // Log(N(x | 0, 1)) = -0.5*x^2 - 0.5*log(2*pi)
      return -0.5 * x[0] * x[0] - 0.5 * MathLog(2.0 * M_PI);
     }
  };

//+------------------------------------------------------------------+
//| Реализация PropRND (Random Walk)                                 |
//+------------------------------------------------------------------+
class RndWalkPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double delta = 4; // Диапазон Random Walk [-4, 4]
      double x0 = x[0];
      int err;
      double u = MathRandomUniform(-delta, delta,err) ; // Uniform(-d, d)
      vector x_new(1);
      x_new[0] = x0 + u; // x' = x + u
      return x_new;
     }
  };

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
// Параметры
   vector start_point(1);
   start_point[0] = 1.0;

// Инициализация объектов
   NormalLogPDF target_pdf;
   RndWalkPropRND proposal_rnd;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = true; // Random Walk - симметричное предложение

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Запуск Random Walk MH для N(0, 1)...");

// Запуск MH
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd,
                       NULL, // log_prop_pdf не нужен, т.к. symmetric = true
                       params))
     {
      Print("=========================================");
      PrintFormat("Вероятность принятия: %.2f%%", accept_rate * 100.0);

      vector mean = samples.Mean(0);
      PrintFormat("Среднее значение (ожидается 0): %.4f", mean[0]);

      Print("Первые 5 сэмплов:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Сэмпл %d: %.4f", i + 1, samples[i][0]);
        }
     }

// Сохранение семплов в CSV
   MatrixToCSV("MH_RW/MH_samples.csv", samples);

   plotTrace(samples);
  }

//+------------------------------------------------------------------+
//|График выборок                                                    |
//+------------------------------------------------------------------+
void plotTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_RW_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Trace Plot MH_RW");
   graphic.BackgroundMainSize(16);

   vector v = smp.Col(0);
   double x[];
   v.Swap(x);
   graphic.CurveAdd(x, CURVE_LINES, "Sample");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }
//+------------------------------------------------------------------+

После получения выборок проводится их анализ с использованием Python. Основные инструменты диагностики:

• Трасса сэмплов (Trace Plot): график значений x(t) по итерациям, позволяющий оценить сходимость цепи и степень ее перемешивания.
• Гистограмма: сравнение эмпирического распределения выборок с теоретической целевой плотностью N(0,1).
• Автокорреляционная функция (ACF): показывает уровень корреляции между последовательными выборками. В RWM автокорреляция обычно высока из-за природы случайного блуждания, поэтому необходимо настроить прореживание (thin).

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf

# File path
file_path = "C:/Program Files/MetaTrader 5/MQL5/Files/MH_RW/MH_samples.csv"

# Load data
data = pd.read_csv(file_path, header=None)
samples = data.iloc[:, 0].values

# Statistics
print("=========================================")
print("MCMC Samples Analysis (MH_RW, N(0, 1)):")
print(f"Number of samples: {len(samples)}")
print(f"Mean: {np.mean(samples):.4f} (expected 0.0000)")
print(f"Std. deviation: {np.std(samples):.4f} (expected 1.0000)")
print("=========================================")

# Plots
fig, axes = plt.subplots(3, 1, figsize=(8, 8))

# Trace Plot
axes[0].plot(samples)
axes[0].set_title('Trace Plot')
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('X')

# Histogram and PDF
x_range = np.linspace(samples.min(), samples.max(), 100)
pdf_theoretical = (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x_range**2)
axes[1].hist(samples, bins=50, density=True, alpha=0.6, label='Samples')
axes[1].plot(x_range, pdf_theoretical, 'r-', label='PDF N(0, 1)')
axes[1].set_title('Histogram and PDF')
axes[1].legend()

# ACF
plot_acf(samples, lags=50, ax=axes[2], title='ACF')
axes[2].set_xlabel('Lags')
axes[2].set_ylabel('Correlation')

plt.tight_layout()
plt.show()

Рис. 3 Трассировочный график и гистограмма выборок с теоретической плотностью

Многомерный случай с симметричным вспомогательным распределением (RWM)

Пример №2

Этот пример демонстрирует применение алгоритма RWM для многомерного пространства с коррелированными переменными.

#include <Math\Stat\Normal.mqh>
#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

//+------------------------------------------------------------------+
//|Класс для двумерного нормального распределения                    |
//+------------------------------------------------------------------+
class BivariateNormalLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      vector mu(2);
      mu[0] = 0;
      mu[1] = 0; // Среднее
      matrix Sigma(2, 2);   // Ковариационная матрица
      Sigma[0][0] = 1;
      Sigma[0][1] = 0.8;
      Sigma[1][0] = 0.8;
      Sigma[1][1] = 1;
      matrix Sigma_inv(2, 2);
      Sigma_inv = Sigma.Inv();
      double det = Sigma.Det();
      vector diff = x - mu;
      double quad_form = diff @ Sigma_inv @ diff;
      return -0.5 * (quad_form + MathLog(det) + 2 * MathLog(2 * M_PI));
     }
  };

//
//+------------------------------------------------------------------+
//|Симметричное вспомогательное распределение (Random Walk)          |
//+------------------------------------------------------------------+
class MultiRndWalkPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double sigma = 1;
      vector delta(x.Size());
      int err;
      for(ulong i = 0; i < x.Size(); i++)
        {
         delta[i] = MathRandomNormal(0, sigma, err);
        }
      return x + delta;
     }
  };

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   vector start_point(2);
   start_point[0] = 1.0;
   start_point[1] = 1.0;
   BivariateNormalLogPDF target_pdf;
   MultiRndWalkPropRND proposal_rnd;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = true;

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Запуск Random Walk MH для двумерного N([0,0], [[1,0.8],[0.8,1]])...");
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd, NULL, params))
     {
      PrintFormat("Вероятность принятия: %.2f%%", accept_rate * 100.0);
      vector mean = samples.Mean(0);
      PrintFormat("Среднее значение (ожидается [0,0]): [%.4f, %.4f]", mean[0], mean[1]);
      Print("Первые 5 сэмплов:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Сэмпл %d: [%.4f, %.4f]", i + 1, samples[i][0], samples[i][1]);
        }

      MatrixToCSV("MH_RW/MH_bivariate_samples.csv", samples);

      plotBivariateTrace(samples);
     }
  }

//+------------------------------------------------------------------+
//| График двумерных выборок                                         |
//+------------------------------------------------------------------+
void plotBivariateTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_Bivariate_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Bivariate Normal Distribution");
   graphic.BackgroundMainSize(16);

   vector x1 = smp.Col(0);
   vector y1 = smp.Col(1);
   double x[],y[];
   x1.Swap(x);
   y1.Swap(y);
   graphic.CurveAdd(x, y, CURVE_POINTS, "Samples");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }

В качестве целевого распределения для простоты выберем двумерное нормальное распределение N(μ, Σ) с параметрами:

• вектор средних μ=[0,0],
• ковариационная матрица Σ= { {1, 0.8}, {0.8, 1} }

Эта матрица отражает сильную корреляцию между двумя переменными, что усложняет процесс семплирования.

Для вычисления логарифма плотности целевого распределения напишем класс BivariateNormalLogPDF, который использует стандартную формулу многомерного нормального распределения:

ln ⁡p(x)= −0.5(x − μ)^T Σ^−1(x − μ) − 0.5ln⁡∣Σ∣ − ln⁡(2π)

Поскольку мы используем симметричный алгоритм RWM, вспомогательное распределение также должно быть симметричным. Оно реализовано в классе MultiRndWalkPropRND. Кандидат x′ генерируется из изотропного (круглого) нормального распределения следующим образом:

x' = x + δ

где,

• δ ∼ N(0,σ^2*I),
• σ^2=1, 
• I — единичная матрица.

Этот выбор соответствует поведению случайного блуждания, центрированного в текущем состоянии x.

Семплирование выполняется с использованием следующих параметров:

• nsamples=5000, burnin=500, thin=10.
• Начальная точка: [1.0,1.0].

Для оценки качества выборок используется Python-скрипт, который строит:

• Трассы выборок: графики значений x1(t) и x2(t) по итерациям для оценки сходимости и перемешивания цепи,
• Автокорреляционная функция (ACF): показывает корреляцию между последовательными выборками для каждой переменной (x1,x2),
• Точечный график: облако выборок с наложенными доверительными эллипсами (95% и 99%) теоретического распределения N(μ, Σ).

Рис.4 График двумерных семплов и доверительные эллипсы

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from matplotlib.patches import Ellipse
from scipy.linalg import eigh

# File path
file_path = "C:/Program Files/MetaTrader 5/MQL5/Files/MH_RW/MH_bivariate_samples.csv"

# Load data
data = pd.read_csv(file_path, header=None)
samples_x1 = data.iloc[:, 0].values
samples_x2 = data.iloc[:, 1].values

# Statistics
print("=========================================")
print("MCMC Samples Analysis (MH_RW, Bivariate):")
print(f"Number of samples: {len(samples_x1)}")
print(f"Mean X1: {np.mean(samples_x1):.4f} (expected 0.0000)")
print(f"Mean X2: {np.mean(samples_x2):.4f} (expected 0.0000)")
print(f"Std. deviation X1: {np.std(samples_x1):.4f} (expected 1.0000)")
print(f"Std. deviation X2: {np.std(samples_x2):.4f} (expected 1.0000)")
print(f"Covariance: {np.cov(samples_x1, samples_x2)[0, 1]:.4f} (expected 0.8000)")
print("=========================================")

# Ellipse plotting function
def plot_ellipse(mean, cov, ax, n_std=1.96, **kwargs):
    vals, vecs = eigh(cov)
    order = vals.argsort()[::-1]
    vals, vecs = vals[order], vecs[:, order]
    theta = np.degrees(np.arctan2(*vecs[:, 0][::-1]))
    width, height = 2 * n_std * np.sqrt(vals)
    ellipse = Ellipse(xy=mean, width=width, height=height, angle=theta, **kwargs)
    ax.add_patch(ellipse)

# Trace Plots and ACF
fig1, axes = plt.subplots(2, 2, figsize=(10, 8))
fig1.suptitle('MCMC Samples Analysis (Trace and ACF)', fontsize=16)

# Trace Plot for X1
axes[0, 0].plot(samples_x1)
axes[0, 0].set_title('Trace Plot (X1)')
axes[0, 0].set_xlabel('Iteration')
axes[0, 0].set_ylabel('X1')

# Trace Plot for X2
axes[0, 1].plot(samples_x2)
axes[0, 1].set_title('Trace Plot (X2)')
axes[0, 1].set_xlabel('Iteration')
axes[0, 1].set_ylabel('X2')

# ACF for X1
plot_acf(samples_x1, lags=50, ax=axes[1, 0], title='ACF (X1)')
axes[1, 0].set_xlabel('Lags')
axes[1, 0].set_ylabel('Correlation')

# ACF for X2
plot_acf(samples_x2, lags=50, ax=axes[1, 1], title='ACF (X2)')
axes[1, 1].set_xlabel('Lags')
axes[1, 1].set_ylabel('Correlation')

plt.tight_layout(rect=[0, 0.03, 1, 0.95])

# Scatter Plot with Ellipses
fig2, ax = plt.subplots(figsize=(8, 8))
ax.scatter(samples_x1, samples_x2, s=10, alpha=0.5, label='Samples')
plot_ellipse([0, 0], np.array([[1, 0.8], [0.8, 1]]), ax, n_std=1.96, edgecolor='red', facecolor='none', label='95% Ellipse')
plot_ellipse([0, 0], np.array([[1, 0.8], [0.8, 1]]), ax, n_std=3.034, edgecolor='green', facecolor='none', label='99% Ellipse')
ax.set_title('Bivariate Scatter Plot with Ellipses')
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.legend()
ax.grid(True)
ax.set_aspect('equal')

plt.show()

Несимметричное вспомогательное распределение

Пример №3

Этот пример демонстрирует необходимость использования полного алгоритма Метрополиса-Гастингса. Это вызвано тем, что целевое распределение имеет ограниченную (положительную) область определения, что требует использования несимметричного предлагающего распределения q.

В качестве целевого распределения выберем одномерное гамма-распределение Gamma(α=3,β=2). Это распределение определено только для x>0.

#include <Math\Stat\Math.mqh>
#include <Math\Stat\Gamma.mqh>
#include <Math\Stat\Lognormal.mqh>
#include <Graphics\Graphic.mqh>
#include <MCMC\MH.mqh>

// Класс для целевого гамма-распределения
class GammaLogPDF : public LogPDF
  {
public:
   virtual double    LogPdf(const vector &x) override
     {
      double alpha = 3.0; // Параметр формы
      double beta = 2.0;  // Параметр масштаба
      return alpha * MathLog(beta) - MathGammaLog(alpha) + (alpha - 1) * MathLog(x[0]) - beta * x[0];
     }
  };

//+--------------------------------------------------------------------------+
//| Класс для несимметричного вспомогательного распределения (логнормальное) |
//+--------------------------------------------------------------------------+
class LogNormalPropRND : public PropRND
  {
public:
   virtual vector    PropRnd(const vector &x) override
     {
      double sigma = 2;
      double mu = MathLog(x[0]);
      vector x_new(1);
      int err;
      x_new[0] = MathRandomLognormal(mu, sigma, err);
      return x_new;
     }
  };

//+------------------------------------------------------------------+
//|Класс для логарифма плотности предлагаемого распределения         |
//+------------------------------------------------------------------+
class LogNormalLogPropPDF : public LogPropPDF
  {
public:
   virtual double    LogPropPdf(const vector &x, const vector &y) override
     {
      double sigma = 2;
      double mu = MathLog(y[0]);
      int err;
      double log_pdf = MathProbabilityDensityLognormal(x[0], mu, sigma, true, err);
      return log_pdf;
     }
  };

//+------------------------------------------------------------------+
//| График выборок                                                   |
//+------------------------------------------------------------------+
void plotTrace(matrix &smp)
  {
   ChartSetInteger(0, CHART_SHOW, false);
   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);
   graphic.Create(0, "MH_Gamma_MCMC", 0, 0, 0, int(width), int(height));
   graphic.BackgroundMain("Gamma Distribution Samples");
   graphic.BackgroundMainSize(16);

   vector x = smp.Col(0);
   double x_arr[], y_arr[];
   x.Swap(x_arr);
   graphic.CurveAdd(x_arr, CURVE_POINTS, "Samples");
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10 * 1000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }

//+------------------------------------------------------------------+
//| Script program start function                                     |
//+------------------------------------------------------------------+
void OnStart()
  {
   vector start_point(1);
   start_point[0] = 1.5;
   GammaLogPDF target_pdf;
   LogNormalPropRND proposal_rnd;
   LogNormalLogPropPDF proposal_pdf;

   MH_Params params;
   params.nsamples = 5000;
   params.burnin = 500;
   params.thin = 10;
   params.symmetric = false;

   MHSampler sampler;
   matrix samples;
   double accept_rate;

   Print("Запуск MH с несимметричным предложением для Gamma(3, 2)...");
   if(sampler.MHSample(start_point, samples, accept_rate, &target_pdf, &proposal_rnd, &proposal_pdf, params))
     {
      PrintFormat("Вероятность принятия: %.2f%%", accept_rate * 100.0);
      PrintFormat("Среднее значение (ожидается %.4f): %.4f", 3.0/2.0, samples.Mean(0)[0]);
      Print("Первые 5 сэмплов:");
      for(int i = 0; i < 5; i++)
        {
         PrintFormat("Сэмпл %d: %.4f", i + 1, samples[i][0]);
        }
      MatrixToCSV("MH_RW/MH_gamma_samples.csv", samples);
      plotTrace(samples);
     }
  }

Логарифм плотности, реализованный в классе GammaLogPDF (унаследованном от абстрактного класса LogPDF), вычисляется как:

ln⁡p(x) = αln⁡β − ln⁡Γ(α) + (α − 1)ln⁡x − βx

В качестве вспомогательного распределения q(x′∣x) возьмем логнормальное распределение LogNormal(lnx,σ2). Это распределение естественно подходит для генерации положительных кандидатов x′ на основе текущего положительного состояния x, так как оно центрирует логарифм нового предложения на логарифме текущего состояния. Однако оно несимметрично, поэтому, нам необходимо реализовать два класса:

• LogNormalPropRND для генерации кандидатов x′,
• LogNormalLogPropPDF для вычисления логарифма плотности предложения.

Семплирование выполняется со следующими параметрами:

• nsamples=5000, burnin=500, thin=10, symmetric=false,
• Начальная точка: x0=1.5

Для оценки работы алгоритма в данном учебном примере мы сравниваем эмпирические результаты с известными ожидаемыми статистическими характеристиками гамма-распределения:

• ожидаемое среднее: E[X] = α/β = 1.5,
• ожидаемая дисперсия: Var[X] = α/β^2 =0.75

В реальных практических задачах (когда точные характеристики p(x) неизвестны) для диагностики качества цепи используются как обычно трассировочные графики, доля принятия и автокорреляционная функция (ACF). Результаты МQL5-скрипта демонстрируют успешную сходимость, доля принятия находится в оптимальном диапазоне (≈30%−40%), эмпирическое среднее близко к ожидаемому значению.

Рис. 5 График трассы и гистограммы выборок с теоретической плотностью Gamma

Independent Metropolis-Hastings

Пример №4

Independent Metropolis-Hastings (IMH) использует независимое вспомогательное распределение q(x'∣x)=q(x'), которое не зависит от текущего состояния x. Это ключевая особенность IMH, в отличие от Random Walk MH.

Благодаря независимости q(x'∣x)=q(x') и q(x∣x')=q(x), логарифм коэффициента принятия задается как:

r = ln⁡(p(x')q(x) / p(x)q(x')) = [ lnp(x′) + lnq(x) ] − [ lnp(x) + lnq(x′) ]

Единственное требование к алгоритму IMH — чтобы q(x′) покрывал поддержку (носитель) целевого распределения p(x).

В качестве целевого распределения выбрано бета-распределение Beta(α=2, β=5), которое определено на единичном интервале x∈[0,1]. Логарифм плотности, реализованный в классе BetaLogPDF, вычисляется с использованием логарифма бета-функции, для численной стабильности:

ln⁡p*(x) = (α−1)ln⁡x + (β−1)ln⁡(1−x) − ln⁡B(α, β)

В качестве вспомогательного распределения q(x′) выберем равномерное распределение Uniform(0,1). Этот выбор идеально соответствует области определения целевого бета-распределения (x∈[0,1]). Он реализован с помощью двух классов:

• UniformPropRND: для генерации кандидатов x′ ∼ Uniform(0,1).
• UniformLogPropPDF: для вычисления lnq(x′).

Важное упрощение: поскольку плотность равномерного распределения q(x)=1 и q(x′)=1 для x∈[0,1], отношение q(x)/q(x′) равно 1, а его логарифм равен 0. Таким образом, в данном частном случае IMH коэффициент принятия упрощается до формы Метрополиса:

r = ln p*(x′) − ln p*(x)

Семплирование выполнялось со следующими параметрами:

• начальная точка: x0 = 0.3
• nsamples=5000, burnin=500, thin=10, symmetric=false (хотя в данном случае, это не имеет значения, так как lnq(x)−lnq(x′)=0).

Ключевое преимущество IMH заключается в том, что он может быстро переходить в любую точку пространства состояний, что, при условии хорошего выбора q(x′), приводит к низкой автокорреляции по сравнению с RWM.

Для оценки успешности семплирования в рамках данного учебного примера мы сравниваем эмпирическое среднее, полученное из выборки, с известным ожидаемым средним значением: E[X]= α + βα ≈ 0.2857.

Эмпирические результаты подтверждают успешное семплирование из целевого Beta(2,5).

Рис.6 График трассы и гистограмма выборок с теоретической плотностью Beta

Выбор алгоритма MCMC MH: RWM vs. IMH

Выбор между Random Walk Metropolis (RWM) и Independent Metropolis-Hastings (IMH) зависит от уровня предварительных знаний о целевом распределении p(x).

RWM является наиболее универсальным выбором и имеет большее практическое значение. 

Кандидат x′ генерируется как "случайное блуждание"  вокруг текущего состояния x. Применяется, когда форма целевого распределения плохо известна, или когда пространство является высокоразмерным и сложным (например, многомодальным). Требует тщательной настройки масштаба шага (δ или σ2) для достижения оптимальной доли принятия, чтобы избежать медленного перемешивания (высокой автокорреляции).

IMH — это более эффективный для генерации независимых выборок, но требующий предварительных знаний алгоритм.

Кандидат x′ генерируется из фиксированного (независимого от x) вспомогательного распределения q(x′). Применяется, когда удается выбрать q(x′), которое хорошо аппроксимирует целевое распределение p(x). При удачном выборе q(x′) достигается очень быстрое перемешивание цепи (низкая автокорреляция). Однако, если q(x′) не покрывает хвосты p(x), цепь может застревать, приводя к неверным результатам.

Заключение

Алгоритм Метрополиса-Гастингса — метод Монте-Карло по схеме марковских цепей (MCMC), позволяет генерировать выборки из сложных распределений, для которых прямое семплирование невозможно. В статье рассмотрены как теоретические основы алгоритма так и его реализация на языке MQL5 в виде класса MHSampler. Этот класс позволяет пользователю, реализуя логику для целевого распределения (LogPDF) и генерации предложений (PropRND, LogPropPDF), легко адаптировать алгоритм к различным вариантам MCMC, включая Random Walk Metropolis и Independent Metropolis-Hastings.

Для наглядной демонстрации работы алгоритма в качестве целевых использовались одномерные и двумерные распределения. Качество генерируемых выборок оценивалось с использованием стандартных инструментов диагностики: графиков трасс, гистограмм и автокорреляционных функций построенных в Python.

Доля принятия выборок (acceptance rate) служила ключевым показателем эффективности цепи и использовалась для оптимизации масштаба предлагающего распределения q(x′∣x) определяемого, например, стандартным отклонением для нормального распределения или шириной интервала для равномерного. Правильный выбор масштаба имеет решающее значение: слишком узкое распределение замедляет исследование пространства состояний, увеличивая автокорреляцию выборок, тогда как слишком широкое снижает долю принятия, что уменьшает вычислительную эффективность.

Хотя алгоритм отличается универсальностью, его эффективность снижается в многомерных задачах и при сильной корреляции между переменными. В таких случаях предпочтительны более специализированные методы, такие как гибридный Монте-Карло (HMC) или выборка по уровням (slice sampling), которые эффективнее используют геометрию целевого распределения и требуют меньше ручной настройки. Тем не менее, благодаря своей простоте, алгоритм Метрополиса-Гастингса остается базовым инструментом широко применяемым в байесовском выводе.

Программы, используемые в статье

#	Имя	Тип	Описание
1	MH.mqh	Включаемый файл	Алгоритм Метрополиса-Гастингса
2	Exmp1_MH_RW.mq5	Скрипт	Пример семплирования одномерного распределения с помощью  RWM
3	Exmp1_Plot_MH_RW.py	Скрипт	Диагностика в Python
4	Exmp2_MH_RW.mq5	Скрипт	Пример семплирования двумерного распределения с помощью алгоритма RWM
5	Exmp2_Plot_MH_RW.py	Скрипт	Диагностика в Python
6	Exmp3_MH_RW.mq5	Скрипт	Пример семплирования одномерного распределения с несимметричным предлагающим распределением
7	Exmp3_Plot_MH_RW.py	Скрипт	Диагностика в Python
8	Exmp4_IMH.mq5	Скрипт	Пример семплирования одномерного распределения с помощью алгоритма IMH
9	Exmp4_Plot_IMH.py	Скрипт	Диагностика в Python