Категориальная скрытая марковская модель на языке MQL5

Введение

Модели пространства состояний (State Space Models, SSM) представляют собой фундаментальный математический аппарат для анализа и прогнозирования последовательных данных. Они широко используются в обработке сигналов, генетике, эконометрике и других областях.
В основе любой модели пространства состояний лежит совместное вероятностное распределение наблюдаемой последовательности y1, …, yT. Ключевое предположение SSM состоит в том, что наблюдаемый процесс порождается последовательностью скрытых (латентных) состояний z1, …, zT, динамика эволюции которых во времени описывается переходной моделью.

В зависимости от характера скрытых состояний (дискретные или непрерывные), типа динамики переходов (линейная или нелинейная) и типа самих наблюдений (категориальные, гауссовские, пуассоновские) модели SSM делятся на разные классы. 

Предметом исследования данной статьи является категориальная скрытая марковская модель (Categorical Hidden Markov Model, HMM) — классический представитель класса SSM с дискретными латентными состояниями и категориальным распределением наблюдаемых данных. Переходы между состояниями в ней описываются с помощью матрицы переходных вероятностей.

Например, модель позволяет по последовательности рыночных событий (белая или чёрная свеча) определять скрытые режимы рынка. В одном режиме (назовём его восходящим трендом) белые свечи встречаются значительно чаще, в другом (нисходящий тренд) — чаще чёрные, в третьем (флэт) — их вероятности равны. Сами режимы мы, конечно, не наблюдаем, они скрыты. Задача модели — угадать текущее состояние рынка, анализируя лишь доступный результат (цвет свечей). 

Для проверки подобных гипотез на языке MQL5 был реализован класс CCategoricalHMM. На простых и наглядных примерах мы разберём основные возможности класса:

• вычисление апостериорных вероятностей скрытых состояний (Inference),
• обучение параметров модели с помощью EM-алгоритма (Fit),
• генерацию синтетических данных (Sample),
• выбор оптимального числа скрытых состояний с помощью критериев AIC и BIC (Model Selection),
• онлайн-фильтрацию скрытых состояний в реальном времени (Filter).

Категориальная скрытая марковская модель

Графически структуру скрытой марковской модели принято представлять в виде направленного ациклического графа.

Рис. 1. HMM представленная в виде графа

На приведённой схеме (рис. 1) элементы графа имеют следующее значение:

• узлы (кружки) — случайные величины,
• белые узлы (zt) — скрытые переменные, которые недоступны для прямого наблюдения,
• темные узлы (yt) — наблюдаемые переменные,
• направленные рёбра (стрелки) — вероятностные зависимости между переменными.

Из графа видно, что каждое скрытое состояние zt зависит только от предыдущего состояния zt−1. Это означает, что последовательность скрытых состояний образует марковскую цепь первого порядка. Кроме того, наблюдаемая переменная yt зависит исключительно от текущего скрытого состояния zt и условно независима от всех предыдущих наблюдений и состояний при известном zt. Это важное предположение модели — условная независимость наблюдений. Формально это означает следующее: если текущий режим системы (zt) известен, то предшествующая история наблюдений и прошлые скрытые состояния не дают дополнительной информации для предсказания текущего наблюдения.

Описанная структура представляет собой классическую модель HMM. Однако существуют и её расширения — например, авторегрессионные скрытые марковские модели (AR-HMM). В них добавляются прямые связи между самими наблюдениями (yt-1 → yt), что позволяет учитывать краткосрочную автокорреляцию в данных.

Класс CCategoricalHMM

Ниже представлено объявление класса и краткое описание его основных компонентов. 

//+------------------------------------------------------------------+
//|                                               CategoricalHMM.mqh |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#include <Math\Stat\Uniform.mqh>
#include <Math\Stat\Gamma.mqh>

//+------------------------------------------------------------------+
//|Categorical Hidden Markov Model                                   |
//|Категориальная скрытая марковская модель                          |
//+------------------------------------------------------------------+
class CCategoricalHMM
  {
public:
   // --- Основные параметры модели ---
   int               m_numStates;      // Количество скрытых состояний
   int               m_numEmissions;   // Количество возможных наблюдаемых символов
   matrix            m_TR;             // Transition Matrix  (m_numStates х m_numStates)
   matrix            m_E;              // Emission Matrix    (m_numStates х m_numEmissions)

   //--- Априоры (регуляризация) ---
   matrix            m_pTR;            // Априорные псевдочастоты для переходов
   matrix            m_pE;             // Априорные псевдочастоты для эмиссий

   vector            m_logliks;        // История лог-правдоподобия
private:
   int               m_seed;           // Seed для воспроизводимости генерации ПСЧ
public:

   // --- Конструкторы ---
                     CCategoricalHMM(int numStates, int numEmissions, int seed = -1);
                     CCategoricalHMM(int numStates, int numEmissions,
                   const matrix &pTR, const matrix &pE, int seed = -1);

   // --- Основные методы ---
   bool              Fit(const vector &y, int maxIter = 100, double tolerance = 1e-6,
                         bool verbose = false); // Baum-Welch training

   bool              Inference(const vector &y, const matrix &tr, const matrix &e,
                               matrix &Sdist, double &LL, matrix &Fdist, matrix &bs, vector &s); // Forward-Backward

   bool              Sample(int T, const matrix &tr, const matrix &e, vector &y, vector &z); // Генерация выборок

   // --- Онлайн-фильтрация ---
   vector            Filter(const vector &PrevState, int nextObs);

   // --- Информационные критерии ---
   double            AIC(double logLike) {return -2.0 * logLike + 2.0 * GetParams();}
   double            BIC(double logLike, int n_samples) {return -2.0 * logLike + GetParams() * MathLog(n_samples);}
   int               GetParams(); // Количество параметров модели

   // --- Инициализация матриц параметров случайными значениями ---
   bool              RandomDirichlet(const double &alpha[], double &result[]);
private:
   matrix            RandomizeMatrix(int rows, int cols, const double &alpha[]);
  };

Основные методы класса:

• Fit() — обучение модели. Реализует классический EM-алгоритм Баума-Велша. По заданной последовательности наблюдений находит оптимальные матрицы переходов (m_TR) и эмиссий (m_E), максимизируя логарифм правдоподобия,
• Inference() — вывод апостериорных распределений. Реализует масштабированный алгоритм Forward-Backward. Вычисляет логарифм маргинального правдоподобия (LL), а также матрицы фильтрации (Fdist), сглаживания (Sdist) и обратных вероятностей, необходимые для обучения и прогнозирования,
• Sample() — генерация синтетических данных. Метод Монте-Карло, который по заданным матрицам TR и E генерирует последовательности скрытых состояний и соответствующих им наблюдений,
• Filter() — онлайн-фильтрация. Рекурсивный байесовский фильтр для работы в реальном времени. На основе предыдущего распределения состояний и нового наблюдения обновляет вероятности текущих скрытых состояний.

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

Для решения этой проблемы в классе реализована байесовская регуляризация с помощью априорных псевдочастот (m_pTR и m_pE). Псевдочастоты выполняют сглаживание оценок и предотвращают вырождение состояний. Даже если определённый паттерн не встречался в истории, он сохраняет минимальную вероятность появления в будущем.

Для случайной инициализации матриц используются вспомогательные методы RandomDirichlet() и RandomizeMatrix(). Распределение Дирихле позволяет генерировать строки матриц, в которых сумма элементов строго равна 1, что гарантирует корректность вероятностных распределений с самого начала обучения.

Кроме того, в классе реализованы методы расчёта информационных критериев AIC и BIC. Они позволяют оценивать оптимальное количество скрытых состояний, находя баланс между качеством подгонки модели и её сложностью.

Эксперимент с игральными костями

Для наглядной демонстрации механизма работы категориальной модели рассмотрим эксперимент с игральными костями.

Предположим, у нас есть две кости с гранями от 1 до 6. Первая — симметричная (Fair), где вероятность выпадения каждой грани равна 1/6. Вторая — асимметричная (Loaded), смещённая в сторону выпадения шестёрки.

По определенным правилам (рис. 2) на каждом шаге выбирается одна из костей, производится бросок, результат записывается, после чего кубик может остаться или быть заменён. Задача заключается в следующем: имея только последовательность выпавших чисел (наблюдения), определить, какой именно кубик использовался в каждом броске.

Переведём условия эксперимента на язык моделей пространства состояний:

• Скрытая переменная (z) принимает два дискретных значения z ∈ {1, 2} где 1 — симметричная кость(Fair), а 2 — смещённая(Loaded),
• Наблюдаемая переменная(y) принимает шесть дискретных значений y ∈ {1, 2, 3, 4, 5, 6}, соответствующих выпавшим граням.

Динамику смены кубиков описывает матрица переходов состояний (TR). На рис. 2 стрелками показана вероятность переходов, где вероятность остаться в первом состоянии равна 0.95, а вероятность перейти во второе состояние 0.05. Аналогично, если мы находимся в данный момент времени во втором состоянии, то вероятность остаться в нем равна 0.9, а вернуться в первое состояние — 0.1.

Рис. 2. Матрица переходов TR и матрица эмиссий E

Вторым ключевым компонентом модели является матрица эмиссий (E). Она задаёт категориальное распределение вероятностей выпадения каждой грани для каждого скрытого состояния:

• для честного кубика (Fair): распределение равномерное — вероятность выпадения любой из шести граней равна 1/6,
• для смещённого кубика (Loaded): вероятность выпадения шестёрки равна 0.5, в то время как вероятности остальных пяти граней равны 0.1.

Процесс вывода скрытых состояний для данного примера демонстрируется в скрипте Inference.mql5.

Вывод скрытых состояний (Inference)

//+------------------------------------------------------------------+
//|                                                    Inference.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <CategoricalHMM\CategoricalHMM.mqh>
#include <Graphics\Graphic.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int n_states = 2;    // Количество скрытых состояний
   int n_emissions = 6; // Количество возможных исходов (количество классов)

// --- 1. Создаем модель (seed=99 для повторяемости) ---
   CCategoricalHMM model(n_states, n_emissions, 99);

// --- Задаем параметры модели ---
// --- 2. Матрица переходов TR ---
   matrix tr = {{0.95, 0.05},
        {0.1,  0.9}
     };
// --- 3.Матрицу эмиссий (E) ---
   matrix e = {{1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6},
        {0.1,   0.1,   0.1,   0.1,   0.1,   0.5}
     };

// --- 4. Генерируем последовательность наблюдений и скрытых состояний ---
   vector y;
   vector z;
   int n_samples = 300; // Количество генерируемых выборок
   model.Sample(n_samples, tr, e, y, z);

   PrintFirst(z,"First 5 states",5);       // First 5 states:     [1 1 1 1 1]
   PrintFirst(y,"First 5 emissions",5);    // First 5 emissions:  [6 1 5 2 3]

// --- БЛОК ВИЗУАЛИЗАЦИИ ---
   int plotLimit = 300; // Рисуем первые 300 точек
   double x[], y_[];
   ArrayResize(x, plotLimit);
   ArrayResize(y_, plotLimit);
   for(int i = 0; i < plotLimit; i++)
     {
      x[i] = i;
      y_[i] = z[i];
     }
   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,"HMM_States_Plot",0,0,0,int(width),int(height));
   CCurve *curve = graphic.CurveAdd(x, y_, ColorToARGB(clrDodgerBlue, 255), CURVE_LINES,"States");
   curve.LinesStyle(STYLE_SOLID);
   curve.LinesWidth(2);
   graphic.BackgroundMainColor(ColorToARGB(clrBlack,255));
   graphic.BackgroundMainSize(24);
   graphic.BackgroundMain("States over time");
   graphic.XAxis().Name("Time");
   graphic.XAxis().NameSize(18);
   graphic.YAxis().Name("State (1=Fair, 2=Loaded)");
   graphic.YAxis().NameSize(18);
   graphic.YAxis().Min(1);
   graphic.YAxis().Max(2);
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(5*1000);
   ChartSetInteger(0,CHART_SHOW,true);
   graphic.Destroy();
   ChartRedraw(0);

// --- Эмпирическая вероятность наблюдения=6 для симметричной(Fair) и смещенной (Loaded) кости ---
   PrintFormat("Empirical probability '6' (Fair): %.3f", CalculateSixFreq(y, z, 1));
   PrintFormat("Empirical probability '6' (Loaded): %.3f", CalculateSixFreq(y, z, 2));

// --- Inference ---
   matrix Filtering, Smoothing, bs;
   double marginal_loglik;
   vector s;

   if(model.Inference(y,tr,e,Smoothing,marginal_loglik,Filtering,bs,s))
     {
      Print("---Inference---");
      PrintFormat("Логарифм правдоподобия (marginal_loglik): %.10f", marginal_loglik);
      Print("\nМатрица Smoothing  (апостериорные вероятности):");
      PrintMatrixFormatted(Smoothing);
      Print("\nМатрица Filtering :");
      PrintMatrixFormatted(Filtering);
     }
//MatrixToCSV("FilteringDist.csv",Filtering);
//MatrixToCSV("SmoothingDist.csv",Smoothing);
//VectorToCSV("TrueStates.csv", z);

  }
//+------------------------------------------------------------------+

//+-------------------------------------------------------------------+
//|Возвращает частоту выпадения шестёрки                              |                                                                |
//+-------------------------------------------------------------------+
double CalculateSixFreq(const vector &r, const vector &s, int stateTarget)
  {
   int count = 0, total = 0;
   for(int i=0; i<(int)r.Size(); i++)
     {
      if(s[i] == stateTarget)
        {
         total++;
         if(r[i] == 6)
            count++;
        }
     }
   return (total > 0) ? (double)count/total : 0;
  }
//+------------------------------------------------------------------+

Создаем экземпляр класса, в конструктор которого передаются параметры: numStates = 2 (количество скрытых состояний) и numEmissions = 6 (количество возможных исходов). Для воспроизводимости результатов фиксируется начальное значение генератора псевдослучайных чисел (seed = 99). Далее задаются истинные параметры модели — матрицы переходов и эмиссий. В данном примере мы предполагаем, что параметры модели нам известны заранее.

Осуществим выборку 300 наблюдений и скрытых состояний из этой модели с помощью метода Sample. На вход метод принимает количество желаемых наблюдений и матрицы параметров. На выходе мы получаем два синхронизированных вектора:

• вектор истинных скрытых состояний z (какой кубик бросали),
• вектор наблюдаемых исходов y (какое число выпало).

На рис. 3 представлена полученная последовательность состояний.

Рис.3. Последовательность скрытых состояний

Прежде чем переходить к выводу скрытых состояний, проверим качество сгенерированной выборки. Для этого рассчитаем эмпирические вероятности выпадения шестёрки в каждом из скрытых состояний: 

• Empirical probability '6' (Fair): 0.173,
• Empirical probability '6' (Loaded): 0.468

Как видно, даже на относительно небольшой выборке в 300 наблюдений эмпирические частоты довольно близки к теоретическим значениям (0.167 и 0.5 соответственно). Это подтверждает корректность работы метода Sample().

//+------------------------------------------------------------------+
//| Генерация последовательности состояний и наблюдений              |
//| согласно заданным матрицам переходов и эмиссий.                  |
//|--------------- IN -----------------------------------------------|
//| int T           - длина генерируемой последовательности          |
//| matrix &tr      - матрица переходных вероятностей                |
//| matrix &e       - матрица вероятностей эмиссий                   |
//|--------------- OUT ----------------------------------------------|
//| vector &y       - сгенерированная последовательность наблюдений  |
//| vector &z       - сгенерированная последовательность состояний   |
//+------------------------------------------------------------------+
bool CCategoricalHMM::Sample(int T, const matrix &tr, const matrix &e,
                             vector &y, vector &z)
  {
   if(T <= 0)
      return false;

   int numStates = (int)tr.Rows();
   int numEmissions = (int)e.Cols();
   y.Resize(T);
   z.Resize(T);

//+------------------------------------------------------------------+
//|1. Подготовка кумулятивных распределений (CDF)                    |
//|Это позволяет быстро семплировать из категориального распределения|
//+------------------------------------------------------------------+
   matrix trc = tr.CumSum(1);   // кумулятивная сумма по строкам (переходы)
   matrix ec  = e.CumSum(1);    // кумулятивная сумма по строкам (эмиссии)

//--- Нормализация кумулятивных сумм (приводим к диапазону [0, 1]) ---
   vector tr_last_col = trc.Col(numStates - 1);
   vector e_last_col = ec.Col(numEmissions - 1);

   for(int r = 0; r < numStates; r++)
     {
      trc.Row(trc.Row(r) / tr_last_col[r], r);
      ec.Row(ec.Row(r) / e_last_col[r], r);
     }

//+------------------------------------------------------------------+
//|2. Генерация случайных чисел                                      |
//+------------------------------------------------------------------+
   double statechange[];
   double randvals[];
   MathRandomUniform(0.0, 1.0, T, statechange);
   MathRandomUniform(0.0, 1.0, T, randvals);

//+------------------------------------------------------------------+
//|3. Генерация последовательности                                   |
//+------------------------------------------------------------------+
   int currentState = 0; // Начинаем с состояния 1
   for(int i = 0; i < T; i++)
     {
      // --- Семплирование следующего состояния ---
      double stateVal = statechange[i];
      int nextState = 0;

      for(int s = numStates - 2; s >= 0; s--)
        {
         if(stateVal > trc[currentState,s])
           {
            nextState = s + 1;
            break;
           }
        }

      // --- Семплирование наблюдения (эмиссии) ---
      double emitVal = randvals[i];
      int emit = 0;

      for(int em = numEmissions - 2; em >= 0; em--)
        {
         if(emitVal > ec[nextState][em])
           {
            emit = em + 1;
            break;
           }
        }

      //--- Сохраняем (формат 1...N) ---
      y[i] = emit + 1;
      z[i] = nextState + 1;
      currentState = nextState;
     }
   return true;
  }

Зачем нужен этап генерации синтетических данных? В контролируемом эксперименте мы получаем не только последовательность наблюдений, но и вектор истинных скрытых состояний z — то есть информацию о том, какой кубик (Fair или Loaded) использовался на каждом шаге. Далее мы скрываем вектор z от модели и передаём в Inference() только наблюдения y (300 значений). Задача алгоритма — по наблюдениям и известным матрицам переходов и эмиссий восстановить наиболее вероятную последовательность скрытых состояний. Сравнивая полученные результаты с истинными значениями, мы можем наглядно оценить качество и возможности категориальной модели.

Метод Inference() возвращает распределение вероятностей скрытых состояний. В зависимости от объёма используемых данных это может быть фильтрующее или сглаженное распределение. 

//+-----------------------------------------------------------------------+
//| Функция Inference реализует алгоритм Forward-Backward(scaled)         |
//| и вычисляет апостериорные вероятности скрытых состояний               |
//|--------------- IN ----------------------------------------------------|
//|vector &y       - последовательность наблюдений (эмиссий)              |
//|matrix &tr      - матрица переходных вероятностей                      |
//|matrix &e       - матрица вероятности эмиссий (Local evidence p(yt|zt))|
//|--------------- OUT ---------------------------------------------------|
//|matrix &Sdist   - апостериорные вероятности (Smoothing) p(z_t|y_{1:T}) |
//|double &LL      - логарифм маргинального правдоподобия  p(y_{1:T}|θ)   |
//|matrix &Fdist   - фильтрующие вероятности (Filtering) p(z_t|y_{1:t})   |
//|matrix &bs      - conditional likelihood p(y_{t+1:T}|z_t)              |
//|vector &s       - вектор коэффициентов масштабирования (Evidence)      |
//+-----------------------------------------------------------------------+
bool CCategoricalHMM::Inference(const vector &y, const matrix &tr, const matrix &e,
                                matrix &Sdist, double &LL, matrix &Fdist, matrix &bs, vector &s)
  {
   int originalL = (int)y.Size();
   if(originalL == 0)
      return false;

   int numStates = (int)tr.Rows();
   int numSymbols = (int)e.Cols();

//--- Проверки размеров ---
   if(tr.Rows() != tr.Cols())
      return false;
   if(e.Rows() != numStates)
      return false;

   int L = originalL + 1;
   vector seq_mod(L);

//--- Добавляем фиктивное наблюдение в начало последовательности для упрощения реализации Forward-Backward ---
   seq_mod[0] = numSymbols + 1;
   for(int i = 0; i < originalL; i++)
      seq_mod[i+1] = y[i];

//--- Подготовка матриц ---
   Fdist.Init(numStates, L);
   Fdist.Fill(0);
   bs.Init(numStates, L);
   s.Init(L);
   s.Fill(0);

// --- FORWARD PASS ---
   Fdist[0,0] = 1.0; // стартуем из состояния 1
   s[0] = 1.0;
//--- Fdist - filtering distribution p(z_t | y_{1:t}) ---
   for(int t = 1; t < L; t++)
     {
      int obs = (int)MathRound(seq_mod[t]) - 1;
      if(obs < 0 || obs >= numSymbols)
         return false;

      for(int state = 0; state < numStates; state++)
        {
         double transitionSum = (Fdist.Col(t-1) @ tr.Col(state));
         Fdist[state,t] = e[state,obs] * transitionSum;
        }

      //--- Нормализация (масштабирование) ---
      s[t] = Fdist.Col(t).Sum(); // EVIDENCE
      if(s[t] > 0)
        {
         vector v_norm = Fdist.Col(t) / s[t];
         Fdist.Col(v_norm, t);
        }
     }

// --- BACKWARD PASS ---
   bs.Fill(1.0); // bs - conditional likelihood p(y_{t+1:T} | z_t)
   for(int t = L - 2; t >= 0; t--)
     {
      int obs_next = (int)MathRound(seq_mod[t+1]) - 1;

      for(int state = 0; state < numStates; state++)
        {
         vector next_val = bs.Col(t+1) * e.Col(obs_next);
         double transitionSum = (tr.Row(state) @ next_val);

         if(s[t+1] > 0)
            bs[state,t] = transitionSum / s[t+1];
        }
     }

// --- Расчет маргинального правдоподобия модели для всей последовательности данных p(y_{1:T}|θ) при заданных параметрах θ = {tr,e} ---
   LL = 0;
   for(int i = 0; i < L; i++)
     {
      if(s[i] > 0)
         LL += MathLog(s[i]);
     }

//+-----------------------------------------------------------------------------+
//|                        Smoothing           Filtering        bs(Backward)    |
//|Posterior inference p(z_t | y_{1:T}) ∝ p(z_t | y_{1:t}) × p(y_{t+1:T} | z_t) |
//|Вычисляем апостериорные вероятности состояний при условии всей {1:T}         |
//|последовательности наблюдений. Благодаря масштабированию результат после     |
//|поэлементного умножения уже нормирован (сумма по состояниям = 1)             |
//+-----------------------------------------------------------------------------+
   Sdist = Fdist * bs;

// --- Удаляем фиктивную первую колонку ---
   matrix finalSdist;
   finalSdist.Init(numStates, originalL);
   for(int c = 0; c < originalL; c++)
     {
      vector v_col = Sdist.Col(c+1);
      finalSdist.Col(v_col, c);
     }
   Sdist = finalSdist;

   return true;
  }

Фильтрация (Filtering distribution)

В теории скрытых марковских моделей фильтрация — это вычисление распределения вероятностей скрытых состояний в момент времени t при условии, что нам известны наблюдения только до момента t включительно:

p(z_t | y_{1:t})

Это распределение вычисляется с помощью алгоритма прямого прохода (Forward Algorithm) и сохраняется в матрицу Fdist. Строки матрицы соответствуют скрытым состояниям, а столбцы — моментам времени, то есть для каждого шага рассчитывается своё отдельное распределение.

Говоря применительно к нашему эксперименту на каждом шаге t алгоритм смотрит на историю выпавших цифр от самого первого броска до текущего и рассчитывает, какова вероятность того, что прямо сейчас кубик честный или смещенный. Важная особенность фильтрации — она не использует будущие данные. Именно поэтому данный метод является основным инструментом для принятия решений в режиме реального времени.

Рис. 4. Фильтрация для смещенной кости

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

Сглаживание (Smoothing distribution)

Сглаживание — это вычисление апостериорных вероятностей скрытого состояния в момент времени t с учётом всех имеющихся наблюдений,  как прошлых, так и будущих:

p(z_t | y_{1:T})

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

Это распределение вычисляется с помощью алгоритма прямого и обратного прохода (Forward-Backward) и сохраняется в матрицу Sdist. Как и в случае с фильтрацией, строки этой матрицы соответствуют скрытым состояниям, а столбцы — моментам времени.

На рис. 5 представлен график сглаженных вероятностей. За счет того, что алгоритм при расчете распределения в точке t уже знает, что произошло в будущем (на шагах t + 1, t + 2 и так далее), график сглаживания выглядит более стабильным и точным по сравнению с фильтрацией, практически безошибочно восстанавливая интервалы, когда использовалась смещенная кость.

Рис. 5. Сглаженные вероятности для смещенной кости

Итак, на текущий момент мы научились создавать модель, задавать её параметры, генерировать синтетические данные и выполнять инференс: вычислять фильтрующие и сглаженные распределения скрытых состояний. Однако в реальных задачах истинные параметры матриц переходов и эмиссий нам неизвестны. В распоряжении имеется только последовательность наблюдений. Поэтому следующим шагом мы переходим к задаче обучения параметров модели по историческим данным.

Обучение параметров модели

Рассмотрим совместное распределение вероятностей для скрытых и наблюдаемых величин на отрезке времени от 1 до T:

Где компоненты распределения определяются следующим образом:

• y1:T — {y1, y2, …, yT} — последовательность дискретных наблюдаемых величин,
• z1:T — {z1, z2, …, zT} — последовательность дискретных скрытых состояний,
• p(z1) — распределение вероятностей начального скрытого состояния (стартовый вектор π),
• p(zt | zt-1) — вероятность перехода из состояния j в момент времени t-1 в состояние k в момент времени t (элемент матрицы переходов TR),
• p(yt | zt) — вероятность получить конкретное наблюдение yt = i при условии, что скрытая переменная находится в состоянии zt = k (элемент матрицы эмиссий E).

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

Мы хотим найти параметры θ = (TR, E), которые лучше всего объясняют наблюдаемые данные, то есть максимизируют правдоподобие p(y1:T∣θ). Чтобы найти это правдоподобие,  необходимо просуммировать совместную вероятность p(y1:T, z1:T|θ) по всем возможным последовательностям скрытых состояний z1:T = (z1, z2, …, zT):

Сделать это прямым перебором невозможно: даже для последовательности длиной всего 100 наблюдений и двух состояний существует 2^100 скрытых траекторий (комбинаций смены состояний). Для эффективного решения этой задачи применяется итерационный EM-алгоритм (Expectation-Maximization), который в контексте скрытых марковских моделей известен как алгоритм Баума-Велша. Вместо полного перебора всех траекторий алгоритм использует принцип динамического программирования. В частности, алгоритм прямого хода (Forward) позволяет эффективно суммировать вероятности всех возможных путей на каждом шаге, вычисляя маргинальное правдоподобие (или его логарифм) за линейное время.

//+------------------------------------------------------------------+
//| Fit — Обучение модели (Baum-Welch Training)                      |
//| Параметры:                                                       |
//|   maxIter       - максимальное количество итераций (100)         |
//|   tolerance     - критерий сходимости (по умолчанию 1e-6)        |
//|   verbose       - выводить подробную информацию о ходе обучения  |
//|                                                                  |
//| Обученные параметры сохраняются во внутренних переменных:        |
//|   m_TR      — матрица переходов                                  |
//|   m_E       — матрица эмиссий                                    |
//|   m_logliks — история лог-правдоподобия                          |
//+------------------------------------------------------------------+
bool CCategoricalHMM::Fit(const vector &y,int maxIter, double tolerance, bool verbose)
  {
   int originalL = (int)y.Size();
   int numStates = m_numStates;
   int numEmissions = m_numEmissions;

//--- Инициализация матрицы переходов, если не задана ---
   if(m_TR.Rows() == 0)
     {
      double default_alpha[] = {1.0}; // по умолчанию используем равномерное распределение Дирихле
      m_TR = RandomizeMatrix(numStates, numStates, default_alpha);
      if(verbose)
        {
         Print("TR matrix was empty, initialized with Dirichlet (alpha=1.0)");
         Print(m_TR);
         // Print(m_TR.Sum(1)); //Проверяем, что сумма вероятностей по строкам матрицы равна 1
        }
     }
   else
      if(m_TR.Rows() != numStates || m_TR.Cols() != numStates)
        {
         PrintFormat("Error: Wrong TR size. Expected %dx%d", numStates, numStates);
         return false;
        }

//--- Инициализация матрицы эмиссий, если не задана ---
   if(m_E.Rows() == 0)
     {
      double default_alpha[] = {1.0};
      m_E = RandomizeMatrix(numStates, numEmissions, default_alpha);
      if(verbose)
        {
         Print("E matrix was empty, initialized with Dirichlet (alpha=1.0).");
         Print(m_E);
         // Print(m_E.Sum(1));
        }
     }
   else
      if(m_E.Rows() != numStates || m_E.Cols() != numEmissions)
        {
         PrintFormat("Error: Wrong E size. Expected %dx%d", numStates, numEmissions);
         return false;
        }

   m_logliks.Init(maxIter);
   m_logliks.Fill(0);

//--- Псевдочастоты (априоры) ---
   matrix pseudoTR = (m_pTR.Rows() > 0) ? m_pTR : matrix::Zeros(numStates, numStates);
   matrix pseudoE  = (m_pE.Rows() > 0)  ? m_pE  : matrix::Zeros(numStates, numEmissions);
   matrix ExpNjk = matrix::Zeros(numStates, numStates);   // E[Njk] - матожидание
   matrix ExpNk = matrix::Zeros(numStates, numEmissions); // E[Nk] - матожидание

   double loglik = 1.0;
   bool converged = false;

//--- Подготовка расширенной последовательности ---
   vector seq_ext;
   seq_ext.Init(originalL + 1);
   seq_ext[0] = 0;
   for(int i = 0; i < originalL; i++)
      seq_ext[i+1] = y[i];

   int iter = 0;
   for(iter = 0; iter < maxIter; iter++)
     {
      double oldLL = loglik;
      loglik = 0;

      matrix oldGuessE = m_E;
      matrix oldGuessTR = m_TR;

      // --- Baum-Welch EM ---
      matrix Fdist, bs, Smoothing;
      vector scale;
      double logL;

      if(!Inference(y, m_TR, m_E, Smoothing, logL, Fdist, bs, scale))
         return false;

      loglik += logL;

      matrix logf = MathLog(Fdist);
      matrix logb = MathLog(bs);
      matrix logE = MathLog(m_E);
      matrix logTR = MathLog(m_TR);

      // --- Expectation Step ---
      // --- 1. КСИ (ξ) ---
      for(int k = 0; k < numStates; k++)
        {
         for(int j = 0; j < numStates; j++)
           {
            for(int t = 0; t < originalL; t++)
              {
               int current_obs = (int)MathRound(seq_ext[t+1]) - 1;
               //--- прошлое(альфа)    переход      эмиссия                будущее ---
               double sumLog = logf[k][t] + logTR[k][j] + logE[j][current_obs] + logb[j][t+1];
               //---  MathExp(sumLog) =  p(y1:T,zt,zt+1) ---
               //---  КСИ (ξ) = p(zt,zt+1| y1:T) ---
               ExpNjk[k][j] += MathExp(sumLog) / scale[t+1]; // НАКОПЛЕНИЕ КСИ (ξ)
              }
           }
        }

      // --- 2. ГАММА (γ) ---
      for(int k = 0; k < numStates; k++)
        {
         for(int s_val = 1; s_val <= numEmissions; s_val++)
           {
            double sum_exp = 0;
            for(int t = 0; t < (int)seq_ext.Size(); t++)
              {
               if(MathRound(seq_ext[t]) == s_val) // Если символ совпал с текущим искомым
                 {
                  sum_exp += MathExp(logf[k][t] + logb[k][t]); // НАКОПЛЕНИЕ ГАММЫ (γ)
                 }
              }
            ExpNk[k][s_val-1] += sum_exp;
           }
        }

      // --- Maximization Step ---
      //--- Обновляем параметры модели ---
      for(int i = 0; i < numStates; i++)
        {
         double eSum = ExpNk.Row(i).Sum(); // Обновление матрицы эмиссий E
         if(eSum > 0)
            m_E.Row(ExpNk.Row(i) / eSum, i);

         double tSum = ExpNjk.Row(i).Sum(); // Обновление матрицы переходов TR
         if(tSum > 0)
            m_TR.Row(ExpNjk.Row(i) / tSum, i);
         else
           {
            m_TR.Row(vector::Zeros(numStates), i);
            m_TR[i][i] = 1.0; // Состояние поглощения, если нет выходов
           }
        }

      m_logliks[iter] = loglik;

      // --- Проверка сходимости ---
      if(verbose)
        {
         if(iter == 0)
           {
            Print("Log Likelihood and Relative Changes in Log Likelihood, Transition Matrix and Emission Matrix");
            Print("   Iteration      LogLik      Log Lik Rel     Transition      Emission");
           }
         else
           {
            double relLL = MathAbs(loglik - oldLL) / (1.0 + MathAbs(oldLL));
            double relTR = (m_TR - oldGuessTR).Norm(MATRIX_NORM_INF) / numStates;
            double relE = (m_E - oldGuessE).Norm(MATRIX_NORM_INF) / numEmissions;
            PrintFormat("     %6d   %12g   %12g  %12g  %12g", iter + 1,loglik, relLL, relTR, relE);
           }
        }

      if(iter > 0 && (MathAbs(loglik - oldLL) / (1.0 + MathAbs(oldLL))) < tolerance)
        {
         if((m_TR - oldGuessTR).Norm(MATRIX_NORM_INF) / numStates < tolerance &&
            (m_E - oldGuessE).Norm(MATRIX_NORM_INF) / numEmissions < tolerance)
           {
            converged = true;
            if(verbose)
               PrintFormat("Converged after %d iterations.", iter + 1);
            break;
           }
        }
      //--- Сброс к псевдосчетчикам или нулю если мы их не используем ---
      ExpNjk = pseudoTR;
      ExpNk  = pseudoE;
     }

   if(!converged)
     {
      if(verbose)
        {
         PrintFormat("Algorithm did not converge with tolerance %g in %d iterations.",tolerance, maxIter);
        }
     }

   m_logliks.Resize(iter < maxIter ? iter + 1 : maxIter); // Удаляем неиспользованные хвосты в logliks

   return true;
  }

Работа алгоритма построена на итерационном выполнении двух шагов: 

1. E-шаг (Expectation — Ожидание)

На этом шаге параметры матриц  m_TR и m_E фиксируются (на первой итерации они инициализируются случайными значениями на основе распределения Дирихле или задаются пользователем). С помощью алгоритма прямого-обратного прохода (Forward-Backward) рассчитываются условные апостериорные вероятности, выступающие в роли достаточных статистик:

Кси ξt(k, j)— вероятность перехода из состояния k в состояние j в момент времени t с учетом всей цепочки наблюдений:

где

• p(zt | y1:t) — фильтрация (forward-переменная) logf[k][t],
• p(zt+1 | zt) — элемент матрицы переходов logTR[k][j],
• p(yt+1 | zt+1) — локальное правдоподобие нового наблюдения, элемент матрицы эмиссий logE[j][current_obs],
• p(yt+2:T | zt+1) — backward-переменная logb[j][t+1]

Гамма (γ) — вероятность нахождения системы в конкретном состоянии на каждом шаге с учетом всей истории (сглаженное распределение/Smoothing)

На основе этих вероятностей накапливаются ожидаемые частоты событий:

• В матрице ExpNjk накапливается математическое ожидание количества переходов между состояниями,
• В матрице ExpNk накапливается математическое ожидание частоты появления каждого символа в каждом состоянии.

Таким образом, вместо реальных (но неизвестных) состояний алгоритм работает с их вероятностными весами. 

2. M-шаг (Maximization — Максимизация)

На этом шаге параметры модели обновляются на основе вычисленных достаточных статистик. Матрицы переходов и эмиссий пересчитываются путём нормировки строк: накопленная взвешенная частота делится на общую сумму по строке.

Если для какого-либо состояния сумма исходящих переходов оказывается равной нулю, оно автоматически переводится в поглощающее состояние (на главной диагонали устанавливается значение 1.0).

Цикл (Считаем вероятности по старым матрицам → Обновляем матрицы по новым вероятностям) повторяется до тех пор, пока относительное изменение лог-правдоподобия (relLL) и параметров матриц (relTR, relE) не станет меньше заданной пользователем точности tolerance. 

Для проверки работы алгоритма сгенерируем ту же последовательность, что использовалась в примере с Inference.mql5 и попробуем восстановить исходные матрицы параметров TR и E (скрипт Learning.mql5).

//+------------------------------------------------------------------+
//|                                                     Learning.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <CategoricalHMM\CategoricalHMM.mqh>
#include <Graphics\Graphic.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int n_states = 2;
   int n_emissions = 6;
   int n_samples = 300;

//--- 1. Создаем модель (seed=99 для повторяемости) ---
   CCategoricalHMM model(n_states, n_emissions, 99);

// --- Задаем параметры модели ---
// --- 2. Устанавливаем матрицу переходов (Tr) ---
   matrix tr = {{0.95, 0.05},
        {0.1,  0.9}
     };

// --- 3.Устанавливаем матрицу эмиссий (E) ---
   matrix e = {{1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6},
        {0.1,   0.1,   0.1,   0.1,   0.1,   0.5}
     };

// --- 4. Генерируем последовательность наблюдений и скрытых состояний ---
   vector y,z;
   model.Sample(n_samples, tr, e, y, z);

// --- 5. Создаем новую модель и обучаем на сгенерированых данных ---
   CCategoricalHMM hmm(n_states, n_emissions);
// --- В качестве начальных параметров зададим: ---
   matrix initialTR = {{0.7, 0.3},
        {0.3, 0.7}
     };
   matrix initialE  = {{0.15, 0.05, 0.1, 0.1, 0.4, 0.2},
        {0.05, 0.15, 0.1, 0.1, 0.3, 0.3}
     };

   hmm.m_TR = initialTR;
   hmm.m_E = initialE;

   double tol = 1e-6;
   if(hmm.Fit(y,100, tol, true))
     {
      Print("\n--- Результаты ---");
      PrintFormat("Final LogLikelihood: %.6f", hmm.m_logliks[hmm.m_logliks.Size() - 1]);
      Print("Estimated TR:\n", hmm.m_TR);
      Print("True TR:\n", tr);
      Print("Estimated E:\n", hmm.m_E);
      Print("True E:\n", e);
     }

// --- БЛОК ВИЗУАЛИЗАЦИИ ---
   int plotLimit = (int)hmm.m_logliks.Size();
   double x[], y_[];
   ArrayResize(x, plotLimit);
   ArrayResize(y_, plotLimit);
   for(int i = 0; i < plotLimit; i++)
     {
      x[i] = i;
      y_[i] = hmm.m_logliks[i];
     }
   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,"HMM_States_Plot",0,0,0,int(width),int(height));
   CCurve *curve = graphic.CurveAdd(x, y_, ColorToARGB(clrDodgerBlue, 255), CURVE_LINES,"LL");
   curve.LinesStyle(STYLE_SOLID);
   curve.LinesWidth(2);
   graphic.BackgroundMainColor(ColorToARGB(clrBlack,255));
   graphic.BackgroundMainSize(24);
   graphic.BackgroundMain("EM Training Progress");
   graphic.XAxis().Name("Iteration");
   graphic.XAxis().NameSize(18);
   graphic.YAxis().Name("Log Likelihood");
   graphic.YAxis().NameSize(18);
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(5*1000);
   ChartSetInteger(0,CHART_SHOW,true);
   graphic.Destroy();
   ChartRedraw(0);
  }
//+------------------------------------------------------------------+

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

В реальных задачах, если начальные параметры не заданы, метод Fit() автоматически инициализирует их с помощью распределения Дирихле. Если вы хотите использовать автоматическую случайную инициализацию, но сохранить воспроизводимость результатов от запуска к запуску, обязательно передавайте фиксированное значение seed > 0 в конструктор класса CCategoricalHMM.

Сравним параметры, которые восстановил алгоритм, с истинными значениями:

        Estimated TR:
        [[0.960567347171196,0.03943265282880401]
         [0.1389165412495843,0.8610834587504157]]
        True TR:
        [[0.95,0.05]
         [0.1,0.9]]
        Estimated E:
        [[0.1796760822470783,0.1183101172299737,0.2125226716621954,0.1291808115499933,0.1924550232457491,0.1678552940650102]
         [2.431736724262778e-19,0.1712486186920382,0.09543183463784231,0.02722383599510184,0.03038479513837311,0.6757109155366445]]
        True E: [[0.1666666666666667,0.1666666666666667,0.1666666666666667,0.1666666666666667,0.1666666666666667,0.1666666666666667]
         [0.1,0.1,0.1,0.1,0.1,0.5]]

Как видим, в общем и целом полученные результаты очень близки к истинным значениям параметров. Модель определила, что первое скрытое состояние генерирует символы относительно равномерно, а второе имеет смещение в сторону шестерки (вероятность 0.675 против истинной 0.5). Конечно, при анализе точности необходимо делать поправку на небольшой объем выборки — всего 300 наблюдений. Тем не менее даже на таком малом массиве данных алгоритм продемонстрировал хорошую сходимость.

На рис. 6 представлена кривая обучения — график изменения лог-правдоподобия. Заданный уровень точности tolerance = 1e-6 был успешно достигнут за 88 итераций.

Рис. 6. Кривая обучения маргинального правдоподобия

В заключение ещё раз отметим, что EM-алгоритм гарантирует сходимость лишь к локальному максимуму целевой функции (правдоподобия) и крайне чувствителен к начальной инициализации. Для поиска глобального оптимума и получения наиболее устойчивых оценок рекомендуется использовать стратегию multi-start. Она заключается в том, чтобы запускать процесс обучения несколько раз, с различными случайными стартовыми параметрами матриц, с последующим выбором той модели, которая показала наилучшее значение Log-Likelihood.

Выбор оптимального количества скрытых состояний (Model Selection)

На практике, работая с реальными рыночными котировками, мы никогда не знаем истинного порождающего процесса. Возникает закономерный вопрос: сколько скрытых состояний необходимо задать модели? Если состояний слишком мало, модель не сможет описать сложную структуру данных (недостаточное обучение). Если их слишком много, модель начнет подстраиваться под случайный шум (переобучение) полностью теряя предсказательную способность.

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

• AIC (Информационный критерий Акаике),
• BIC (Байесовский информационный критерий).

Оба критерия балансируют между качеством подгонки модели и её сложностью, штрафуя за увеличение числа параметров. Оптимальной считается модель с минимальным значением AIC или BIC.

Рассмотрим процесс выбора модели (скрипт ModelSelection.mq5).

//+------------------------------------------------------------------+
//|                                               ModelSelection.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <CategoricalHMM\CategoricalHMM.mqh>
#include <Graphics\Graphic.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int n_samples = 5000;
   int n_states_true = 3;
   int n_emissions = 6;

// --- 1. ИНИЦИАЛИЗАЦИЯ ИСТИННОЙ МОДЕЛИ ДЛЯ ГЕНЕРАЦИИ ---
   matrix tr_true =
     {
        {0.8, 0.1, 0.1},
        {0.05, 0.9, 0.05},
        {0.1, 0.1, 0.8}
     };

   matrix e_true =
     {
        {1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6, 1.0/6},
        {0.1, 0.1, 0.1, 0.1, 0.1, 0.5},
        {0.4, 0.1, 0.1, 0.1, 0.1, 0.2}
     };

   vector y, states_true;
   CCategoricalHMM gen_model(n_states_true, n_emissions, 42); // Seed 42
   gen_model.Sample(n_samples, tr_true, e_true, y, states_true);
//VectorToCSV("Sample5000.csv",y);

// --- РАСЧЕТ ИСТИННОГО LOG-LIKELIHOOD ---
   double true_ll = 0;
   matrix p_true, fs_true, bs_true;
   vector s_true;

// --- Вызываем Inference на истинных параметрах ---
   gen_model.Inference(y, tr_true, e_true, p_true, true_ll, fs_true, bs_true, s_true);

   PrintFormat("==================================================");
   PrintFormat("TRUE MODEL LOG-LIKELIHOOD: %.2f", true_ll);
   PrintFormat("==================================================");

// --- 2. ПОДБОР МОДЕЛИ (MODEL SELECTION) ---
   int ns[] = {1,2, 3, 4, 5, 6, 7};
   int n_fits = 10;
   int n_points = ArraySize(ns);
   double x_vals[], aic_vals[], bic_vals[]; // Массивы для графика
   ArrayResize(x_vals, n_points);
   ArrayResize(aic_vals, n_points);
   ArrayResize(bic_vals, n_points);

   Print("==================================================");
   Print("HMM MODEL SELECTION ");
   PrintFormat("True States: %d | Samples: %d", n_states_true, n_samples);
   Print("==================================================");

   for(int i = 0; i < ArraySize(ns); i++)
     {
      int n_curr = ns[i];
      double best_ll = -DBL_MAX;
      int n_params;

      // --- Пробуем несколько случайных стартов для каждого N ---
      for(int fit = 0; fit < n_fits; fit++)
        {
         CCategoricalHMM model(n_curr, n_emissions, fit);
         n_params = model.GetParams();

         // --- Обучаем (20 итераций для быстрой оценки) ---
         if(!model.Fit(y, 20, 1e-4,false))
            continue;

         // --- Получаем Log-Likelihood ---
         double current_ll = 0;
         matrix p, fs, bs;
         vector s;
         if(model.Inference(y, model.m_TR, model.m_E, p, current_ll, fs, bs, s))
           {
            if(current_ll > best_ll)
               best_ll = current_ll;
           }
        }

      // --- Оценка информационных критериев для лучшей попытки в текущем N ---
      CCategoricalHMM eval_model(n_curr, n_emissions);
      double val_aic = eval_model.AIC(best_ll);
      double val_bic = eval_model.BIC(best_ll, n_samples);

      // --- Сохраняем данные для графика ---
      x_vals[i] = n_curr;
      aic_vals[i] = val_aic;
      bic_vals[i] = val_bic;

      PrintFormat("N=%d | LL: %.2f | AIC: %.2f | BIC: %.2f | n_params:%d",
                  n_curr, best_ll, val_aic, val_bic,n_params);
     }
   Print("==================================================");

// --- ВИЗУАЛИЗАЦИЯ ---
   ChartSetInteger(0, CHART_SHOW, false);

   CGraphic graphic;
   ulong width = ChartGetInteger(0, CHART_WIDTH_IN_PIXELS);
   ulong height = ChartGetInteger(0, CHART_HEIGHT_IN_PIXELS);

   if(ObjectFind(0, "HMM_Selection_Plot") < 0)
      graphic.Create(0, "HMM_Selection_Plot", 0, 0, 0, int(width), int(height));
   else
      graphic.Attach(0, "HMM_Selection_Plot");

   CCurve *curveAIC = graphic.CurveAdd(x_vals, aic_vals, ColorToARGB(clrDodgerBlue, 255), CURVE_POINTS_AND_LINES, "AIC");
   curveAIC.LinesStyle(STYLE_SOLID);

   CCurve *curveBIC = graphic.CurveAdd(x_vals, bic_vals, ColorToARGB(clrLimeGreen, 255), CURVE_POINTS_AND_LINES, "BIC");
   curveBIC.LinesStyle(STYLE_SOLID);

   graphic.BackgroundMainColor(ColorToARGB(clrBlack, 255));
   graphic.BackgroundMainSize(24);
   graphic.BackgroundMain("HMM Model Selection (AIC vs BIC)");
   graphic.XAxis().Name("Number of States");
   graphic.XAxis().NameSize(18);
   graphic.YAxis().Name("AIC/BIC Value (Lower is Better)");
   graphic.YAxis().NameSize(18);
   graphic.CurvePlotAll();
   graphic.Update();
   Sleep(10000);
   ChartSetInteger(0, CHART_SHOW, true);
   graphic.Destroy();
   ChartRedraw(0);
  }
//+------------------------------------------------------------------+

Сгенерируем выборку объёмом 5000 наблюдений из заведомо известной модели HMM с тремя состояниями (N = 3):

• Состояние 1 (Честный кубик): Равномерное распределение граней (1/6),
• Состояние 2 (Асимметрия): Смещён в сторону выпадения шестёрок (0.5),
• Состояние 3 (Асимметрия): Смещён в сторону выпадения единицы (0.4).

Затем мы последовательно обучим 7 различных архитектур HMM (с количеством состояний N от 1 до 7). Чтобы снизить влияние локальных максимумов, каждая архитектура обучается по 10 раз (n_fits = 10) из случайных стартовых точек. В качестве финального результата отбирается попытка с максимальным лог-правдоподобием (LL). Количество итераций обучения было ограничено до 20 для ускорения эксперимента.

При увеличении числа состояний log-likelihood закономерно растёт:

• Модель с одним состоянием (N=1): превращает матрицу переходов в скаляр 1×1 со значением 1.0, полностью исключая динамику смены режимов. Матрица эмиссий становится обычным вектором частот. Модель просто подсчитывает базовые вероятности выпадения граней во всём массиве данных, затрачивая на обучение одну итерацию. (LL = −8425.95).
• Переход к двум состояниям (N=2): заметного улучшения практически нет (LL = −8422.48).
• Переход к трём состояниям (N=3): происходит значительный рост LL с −8422 до −8397. Это сигнал, что данная модель лучше объясняет данные.
• Дальнейший рост (N≥4): лог-правдоподобие в основном продолжает расти (как показано на рис. 7), так как избыточные параметры позволяют точнее подгоняться под случайные флуктуации выборки.

Рис. 7. Информационные критерии и лог-правдоподобие для различных скрытых состояний

Однако полагаться только на максимум правдоподобия нельзя — это прямой путь к переобучению. Иначе мы бы просто выбрали модель с 6 скрытыми состояниями. Теперь рассмотрим полученные результаты по информационным критериям (рис.8).

Рис. 8. Лучшие результаты AIC/BIC и LL для 7 моделей

AIC уверенно показывает минимум на истинном значении N=3 (AIC = 16837.14), подтверждая свою способность выявлять реальную сложность процесса. BIC, напротив, отдаёт предпочтение самой простой модели N=1 (BIC = 16894.48), помещая истинную модель (N=3) лишь на третье место. Такое расхождение объясняется двумя причинами:

Во-первых, более жёсткий штраф BIC: этот критерий масштабируется как ln(T), где T — объём выборки. При 5000 наблюдений штраф за каждый дополнительный параметр составляет примерно 8.51. Переход от модели с 1 состоянием к модели с 3 состояниями добавляет 16 параметров. Прироста лог-правдоподобия оказалось недостаточно, чтобы компенсировать столь существенное увеличение сложности модели.

Во-вторых, незавершённость обучения  — ограничение в 20 итераций не позволило моделям с N=3 полностью сойтись к локальному оптимуму. Истинная модель имеет LL ≈ −8376.88, тогда как лучшая обученная — только −8397.57. Разница почти в 21 пункт существенно повлияла на итоговые значения критериев.

Тем не менее главный вывод остаётся ясным: при выборе архитектуры HMM не следует гнаться за максимальным правдоподобием. Правильный подход — использовать информационные критерии AIC и BIC, которые помогают найти разумный компромисс между точностью подгонки модели к данным и её сложностью.

Онлайн-фильтрация скрытых состояний и байесовский фильтр

Онлайн-фильтрация в контексте скрытых марковских моделей — это процесс рекурсивного обновления вероятностей скрытых состояний в режиме реального времени по мере поступления новых наблюдений. В отличие от офлайн-методов (например, сглаживания), выполняющих пакетный анализ всей выборки сразу, онлайн-фильтрация работает последовательно во времени (t = 1, 2, 3, …, T).

Главная задача онлайн-фильтрации — вычисление текущего апостериорного распределения P(zt | y1:t). Этот процесс реализуется рекурсивно и на каждом временном шаге состоит из двух фаз:

• Шаг прогноза (Prediction) — расчёт априорного распределения,
• Шаг фильтрации (Update) — коррекция прогноза с помощью нового наблюдения по формуле Байеса (расчет апостериорного распределения).

В рамках HMM классическая теорема Байеса для шага онлайн-фильтрации имеет следующий вид:

где

• Апостериорная вероятность p(z_t = j | y_{1:t}) — уточнённая (отфильтрованная) вероятность того, что система в момент времени t находится в состоянии j, с учетом всей доступной истории, включая последнее наблюдение y_t,
• Локальное правдоподобие p(y_t | z_t = j) — элемент матрицы эмиссий (E), отражающий вероятность генерации символа y_t при условии, что система находится в скрытом состоянии j,
• Априорная вероятность (Шаг прогноза) p(z_t = j | y_{1:t-1}) — предсказание распределения состояний на шаг t, сформированное на основе информации из прошлого шага (t−1). Вычисляется через матричное умножение предыдущего фильтра на матрицу переходов.

• Маргинальное правдоподобие (Evidence) p(y_t |y_{1:t-1}) — полная вероятность увидеть наблюдение y_t при текущих параметрах модели. Это нормировочная константа, которая гарантирует, что сумма вероятностей всех состояний на текущем шаге будет равна 1.

Для демонстрации процесса онлайн-фильтрации рассмотрим скрипт Filter.mql5.

//+------------------------------------------------------------------+
//|                                                       Filter.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"
#include <CategoricalHMM\CategoricalHMM.mqh>
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//--- 1. параметры модели / Model parameters ---
   matrix tr = {{0.95, 0.05},
        {0.10, 0.90}
     };
   matrix e = {{0.16, 0.16, 0.16, 0.16, 0.16, 0.20},
        {0.10, 0.10, 0.10, 0.10, 0.10, 0.50}
     };
   vector y = {1, 6, 6, 6, 2}; // Последовательность наблюдений / Sequence of observations

   CCategoricalHMM hmm(2, 6);
   hmm.m_TR = tr;
   hmm.m_E = e;

//--- Стартовое состояние  / Initial state ---
   vector currentFilterState = {1.0, 0.0};
   vector cfs = currentFilterState;
   Print("=== ПОСЛЕДОВАТЕЛЬНЫЙ ОНЛАЙН ПРОГНОЗ ===");

   for(int i = 0; i < (int)y.Size(); i++)
     {
      int obs = (int)y[i];

      // --- ШАГ 1: ПРОГНОЗ (до того, как узнали текущее наблюдение) ---
      vector statePrediction = currentFilterState @ hmm.m_TR;
      PrintFormat("Время t=%d:", i+1);
      PrintFormat("  Прогноз состояния : [1]%.4f, [2]%.4f", statePrediction[0], statePrediction[1]);

      vector obsPrediction = statePrediction @ hmm.m_E; // Умножаем прогноз состояний на матрицу эмиссий

      ulong maxObsIdx = obsPrediction.ArgMax(); // Находим самое вероятное наблюдение
      PrintFormat("  Самое вероятное след. наблюдение: %d (с вероятностью %.4f)", maxObsIdx + 1, obsPrediction[maxObsIdx]);

      // --- ШАГ 2: ОБНОВЛЕНИЕ (ФИЛЬТРАЦИЯ) ---
      // --- пришло наблюдение obs, обновляем текущее состояние ---
      currentFilterState = hmm.Filter(currentFilterState, obs);

      PrintFormat("  РЕАЛЬНОСТЬ: Пришло наблюдение %d. Текущее состояние уточнено (Фильтрация): [1]%.4f, [2]%.4f",
                  obs, currentFilterState[0], currentFilterState[1]);
      Print("---------------------------------------");
     }

// --- ВАРИАНТ 1: ПОСЛЕДОВАТЕЛЬНЫЙ (ONLINE) ---
   Print("=== 1. ONLINE FILTER (STEP-BY-STEP) ===");
   currentFilterState = cfs;
   for(int i = 0; i < (int)y.Size(); i++)
     {
      currentFilterState = hmm.Filter(currentFilterState, (int)y[i]);
      PrintFormat("Время t = %d (after y = %d): [1]%.4f, [2]%.4f", i+1, (int)y[i], currentFilterState[0], currentFilterState[1]);
     }

// --- ВАРИАНТ 2: ПАКЕТНЫЙ (Inference) ---
   Print("\n=== 2. BATCH Inference (Filtering Distribution) ===");
   matrix pStates, fs, bs;
   vector s_vec;
   double logP;

   if(hmm.Inference(y, tr, e, pStates, logP, fs, bs, s_vec))
     {
      PrintMatrixFormatted(fs,4);
     }
  }

//-------------Результаты выполнения скрипта ----------------------

// === ПОСЛЕДОВАТЕЛЬНЫЙ ОНЛАЙН ПРОГНОЗ ===
// Время t=1:
//   Прогноз состояния : [1]0.9500, [2]0.0500
//   Самое вероятное след. наблюдение: 6 (с вероятностью 0.2150)
//   РЕАЛЬНОСТЬ: Пришло наблюдение 1. Текущее состояние уточнено (Фильтрация): [1]0.9682, [2]0.0318
// ---------------------------------------
// Время t=2:
//   Прогноз состояния : [1]0.9229, [2]0.0771
//   Самое вероятное след. наблюдение: 6 (с вероятностью 0.2231)
//   РЕАЛЬНОСТЬ: Пришло наблюдение 6. Текущее состояние уточнено (Фильтрация): [1]0.8273, [2]0.1727
// ---------------------------------------
// Время t=3:
//   Прогноз состояния : [1]0.8032, [2]0.1968
//   Самое вероятное след. наблюдение: 6 (с вероятностью 0.2590)
//   РЕАЛЬНОСТЬ: Пришло наблюдение 6. Текущее состояние уточнено (Фильтрация): [1]0.6201, [2]0.3799
// ---------------------------------------
// Время t=4:
//   Прогноз состояния : [1]0.6271, [2]0.3729
//   Самое вероятное след. наблюдение: 6 (с вероятностью 0.3119)
//   РЕАЛЬНОСТЬ: Пришло наблюдение 6. Текущее состояние уточнено (Фильтрация): [1]0.4022, [2]0.5978
// ---------------------------------------
// Время t=5:
//   Прогноз состояния : [1]0.4418, [2]0.5582
//   Самое вероятное след. наблюдение: 6 (с вероятностью 0.3674)
//   РЕАЛЬНОСТЬ: Пришло наблюдение 2. Текущее состояние уточнено (Фильтрация): [1]0.5588, [2]0.4412
// ---------------------------------------
// === 1. ONLINE FILTER (STEP-BY-STEP) ===
// Шаг t = 1 (after y = 1): [1]0.9682, [2]0.0318
// Шаг t = 2 (after y = 6): [1]0.8273, [2]0.1727
// Шаг t = 3 (after y = 6): [1]0.6201, [2]0.3799
// Шаг t = 4 (after y = 6): [1]0.4022, [2]0.5978
// Шаг t = 5 (after y = 2): [1]0.5588, [2]0.4412
//
// === 2. BATCH Inference (Filtering Distribution) ===
//   1.0000 0.9682 0.8273 0.6201 0.4022 0.5588
//   0.0000 0.0318 0.1727 0.3799 0.5978 0.4412

Мы предполагаем, что модель уже предварительно обучена, а матрицы переходов (TR) и эмиссий (E) нам известны. Скрипт имитирует поступление данных в реальном времени на примере последовательности: y = {1, 6, 6, 6, 2}.

Процесс обработки каждого временного интервала разделён на два шага:

• Шаг прогноза (до закрытия текущего бара t): Имея вектор распределения состояний с предыдущего шага currentFilterState (на шаге t=1 он инициализируется начальными вероятностями, например {1.0, 0.0}), мы выполняем его матричное умножение на матрицу переходов. Прогноз скрытого состояния строится до того, как закроется текущий бар и будет получено новое наблюдение. Дополнительно на этом этапе можно рассчитать, какое именно наблюдение (какой исход) является наиболее вероятным на текущем баре.
• Шаг фильтрации (после закрытия бара t): Как только бар сформирован, мы получаем фактическое значение наблюдения y. Используя формулу Байеса, модель выполняет коррекцию — уточняет вероятности состояний для текущего момента t с учетом этого наблюдения.
• Переход к шагу t+1: Полученный апостериорный вектор становится основой для следующего шага. С его помощью мы делаем прогноз на момент времени t+1, и система переходит в режим ожидания нового наблюдения.

Цикл (прогноз → закрытие бара → фильтрация → прогноз на следующий шаг) повторяется для каждого нового наблюдения {t+1, t+2, …}, обновляя оценку скрытого состояния в режиме онлайн.

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

Заключение

В рамках данной статьи мы детально разобрали теоретические основы и практическую реализацию скрытых марковских моделей с категориальными эмиссиями (Categorical HMM) на языке MQL5. В ходе исследования мы последовательно рассмотрели ключевые задачи, которые решает модель:

• вывод скрытых состояний (метод Inference) — расчет апостериорных распределений на основе масштабированного алгоритма прямого и обратного прохода (Forward-Backward),
• обучение параметров (метод Fit) — итерационный EM-алгоритм, включающий стохастическую инициализацию матриц через распределение Дирихле и регуляризацию с помощью априорных псевдочастот,
• выбор модели (Model Selection) — выбор оптимального числа скрытых состояний с помощью информационных критериев AIC и BIC, защищающих модель от переобучения,
• онлайн-фильтрацию (метод Filter) — рекурсивный байесовский фильтр, позволяющий в режиме реального времени обновлять вероятности скрытых состояний по мере поступления новых данных.

Категориальная модель предоставляет трейдеру строгий математический инструмент для распознавания скрытых состояний по наблюдаемым величинам. Понимание принципов работы данной модели закладывает прочный фундамент для дальнейшего изучения и разработки более сложных моделей пространства состояний.

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

#	Имя	Тип	Описание
1	CategoricalHMM.mqh	Класс	Класс категориальной скрытой марковской модели
2	Inference.mq5	Скрипт	Пример вывода скрытых состояний и генерации синтетических данных из заданной модели
3	Learning.mq5	Скрипт	Пример обучения модели
4	ModelSelection.mq5	Скрипт	Пример выбора оптимального числа скрытых состояний с помощью AIC/BIC критериев
5	Filter.mq5	Скрипт	Пример построения рекурсивного байесовского фильтра
6	MQL5.zip	Архив	Все вышеперечисленные файлы Вы сможете найти в архиве MQL5

Проект поддерживается в Algo Forge