//+------------------------------------------------------------------+ //| CategoricalHMM.mqh | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #include #include //+------------------------------------------------------------------+ //|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[]); }; //+------------------------------------------------------------------+ //| Конструктор 1: Без псевдочастот (априоры по умолчанию = 0) | //+------------------------------------------------------------------+ CCategoricalHMM::CCategoricalHMM(int numStates, int numEmissions, int seed = -1): m_numStates(numStates), m_numEmissions(numEmissions), m_seed(seed) { m_pTR.Init(0, 0); m_pE.Init(0, 0); // --- Инициализация генератора ПСЧ --- if(seed != -1) MathSrand(seed); // фиксированный seed для воспроизводимости else MathSrand(GetTickCount()); } //+------------------------------------------------------------------+ //| Конструктор 2: С заданными априорными псевдочастотами | //+------------------------------------------------------------------+ CCategoricalHMM::CCategoricalHMM(int numStates, int numEmissions, const matrix &pTR, const matrix &pE,int seed = -1): m_numStates(numStates), m_numEmissions(numEmissions), m_seed(seed) { // --- Инициализация псевдочастот для переходов --- if(pTR.Rows() == (ulong)m_numStates && pTR.Cols() == (ulong)m_numStates) m_pTR = pTR; else { PrintFormat("HMM Error: Invalid pTR size (%dx%d). Expected (%dx%d). Pseudo-counts reset to zero", pTR.Rows(), pTR.Cols(), m_numStates, m_numStates); m_pTR = matrix::Zeros(m_numStates, m_numStates); } // --- Инициализация псевдочастот для эмиссий --- if(pE.Rows() == (ulong)m_numStates && pE.Cols() == (ulong)m_numEmissions) m_pE = pE; else { PrintFormat("HMM Error: Invalid pE size (%dx%d). Expected (%dx%d). Pseudo-counts reset to ZERO", pE.Rows(), pE.Cols(), m_numStates, m_numEmissions); m_pE = matrix::Zeros(m_numStates, m_numEmissions); } // --- Инициализация генератора --- if(seed != -1) MathSrand(seed); // фиксированный seed для воспроизводимости else MathSrand(GetTickCount()); } //+------------------------------------------------------------------+ //| 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; } //+-----------------------------------------------------------------------+ //| Функция 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; } //+------------------------------------------------------------------+ //| Генерация последовательности состояний и наблюдений | //| согласно заданным матрицам переходов и эмиссий. | //|--------------- 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; } //+------------------------------------------------------------------+ //| Количество параметров модели | //+------------------------------------------------------------------+ int CCategoricalHMM:: GetParams() { int n_transitions = m_numStates * (m_numStates - 1); // Параметры матрицы переходов int n_emissions = m_numStates * (m_numEmissions - 1); // Параметры матрицы эмиссий return n_transitions + n_emissions; } //+------------------------------------------------------------------+ //| Генерация случайной матрицы с помощью Dirichlet Distribution | //+------------------------------------------------------------------+ matrix CCategoricalHMM::RandomizeMatrix(int rows, int cols, const double &alpha[]) { matrix mat(rows, cols); double result[]; //--- Проверяем, один ли параметр alpha для всех столбцов --- bool single_alpha = (ArraySize(alpha) == 1); double local_alpha[]; ArrayResize(local_alpha, cols); if(single_alpha) ArrayFill(local_alpha, 0, cols, alpha[0]); for(int i = 0; i < rows; i++) { if(!single_alpha) ArrayCopy(local_alpha, alpha); if(RandomDirichlet(local_alpha, result)) { vector row_vec; row_vec.Assign(result); mat.Row(row_vec, i); } } return mat; } //+-----------------------------------------------------------------+ //| Генерация случайного вектора из распределения Дирихле | //+-----------------------------------------------------------------+ //| Dirichlet Distribution — это многомерное обобщение Beta- | //| распределения. Используется для генерации случайных вероятностей| //| (строк матрицы переходов tr или матрицы эмиссий e) | //+-----------------------------------------------------------------+ bool CCategoricalHMM::RandomDirichlet(const double &alpha[], double &result[]) { int n = ArraySize(alpha); if(n < 1) { Print("RandomDirichlet: alpha array is empty"); return false; } ArrayResize(result, n); double sum = 0.0; double gamma_val; int error_code = 0; // --- Генерация через гамма-распределение --- for(int i = 0; i < n; i++) { if(alpha[i] <= 0.0) { Print("RandomDirichlet: alpha[", i, "] must be > 0"); return false; } gamma_val = MathRandomGamma(alpha[i], 1.0, error_code); result[i] = gamma_val; sum += gamma_val; } // --- Нормализация [0, 1] --- for(int i = 0; i < n; i++) result[i] /= sum; // --- Корректировка последней компоненты --- double total = 0.0; for(int i = 0; i < n-1; i++) total += result[i]; result[n-1] = 1.0 - total; return true; } //+------------------------------------------------------------------+ //| Filter — Один шаг онлайн-фильтрации (Forward Filtering Step) | //| Вычисляет апостериорное распределение скрытых состояний | //| p(z_t | y_{1:t}) после поступления нового наблюдения. | //| Реализует рекурсивное обновление байесовского фильтра: | //| p(z_t|y_{1:t}) ∝ p(y_t|z_t) * p(z_t | y_{1:t-1}) | //| p(z_t | y_{1:t-1}) = Σ[p(z_t|z_{t-1}) * p(z_{t-1}|y_{1:t-1})] | //|--------------- IN -----------------------------------------------| //| vector &previousFilterProb - вектор вероятностей состояний | //| с предыдущего шага p(z_{t-1} | y_{1:t-1}) | //| int y - текущее наблюдение | //|--------------- OUT ----------------------------------------------| //| vector - обновлённые вероятности состояний | //| p(z_t | y_{1:t}) — filtering distribution | //+------------------------------------------------------------------+ vector CCategoricalHMM::Filter(const vector &previousFilterProb, int y) { // --- 1. Prediction Step: --- // --- p(z_t | y_{1:t-1}) = Σ[p(z_t|z_{t-1}) * p(z_{t-1}|y_{1:t-1})] = previousFilterProb @ TR --- vector statePrediction = previousFilterProb @ m_TR ; // вектор априорных вероятностей состояний // --- 2. Likelihood p(y_t|z_t): Извлекаем столбец из матрицы эмиссий для текущего наблюдения --- int obsIdx = y - 1; vector likelihoods = m_E.Col(obsIdx); // --- 3. Update Step: Поэлементное умножение --- // --- p(z_t | y_{1:t}) ~ statePrediction * likelihoods = p(z_t|y_{1:t-1}) * p(y_t|z_t) --- vector currentFilterProb = statePrediction * likelihoods; // --- 4. Normalization: Делим на сумму (Evidence) = p(y_t|y_{1:t-1}) --- double evidence = currentFilterProb.Sum(); currentFilterProb = currentFilterProb / evidence; return currentFilterProb; } //+------------------------------------------------------------------+ //+------------------------------------------------------------------+ //| Печатает первые n элементов вектора | //+------------------------------------------------------------------+ void PrintFirst(const vector &vec,string str="", int n = 5) { n = (int)MathMin(n, vec.Size()); vector temp; temp.Copy(vec); temp.Resize(n); Print(str,": ",temp); } //+------------------------------------------------------------------+ //| функция для форматированного вывода матрицы с заданной точностью | //+------------------------------------------------------------------+ void PrintMatrixFormatted(const matrix &m, int c = 4) { for(ulong i=0; i