ARD-регрессия на MQL5: автоматический отбор признаков и случайные признаки Фурье

Введение

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

Традиционные методы снижения размерности и отбора признаков имеют заметные ограничения. Например, метод главных компонент (PCA) хорошо сжимает пространство признаков, выбирая направления с максимальной дисперсией. Однако он не учитывает целевую переменную и ориентируется только на признаки. В результате PCA часто сохраняет шум с высокой дисперсией и может отбросить полезные признаки, которые имеют относительно низкую дисперсию, но при этом сильно связаны с целевой переменной. Кроме того, главные компоненты представляют собой линейные комбинации исходных переменных, что затрудняет их интерпретируемость.

Более гибким и элегантным решением является ARD-регрессия (Automatic Relevance Determination — Автоматическое определение релевантности). ARD — это полноценный байесовский подход к линейной регрессии, который автоматически оценивает степень важности (релевантности) каждого входного признака в процессе обучения модели.

В этой статье мы подробно разберём математический аппарат ARD-регрессии и реализуем модель ARDRegression на языке MQL5, после чего сравним точность вычислений и скорость работы нашей реализации с эталонным модулем из библиотеки scikit-learn.

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

Что такое ARD-регрессия?

Чтобы понять концепцию ARD-регрессии, стоит начать с классической линейной регрессии. В стандартной постановке мы ищем линейную зависимость между вектором признаков x=[x1, x2, …, xD] и целевой переменной y:

где,

• D – размерность пространства признаков (количество предикторов),
• w = [w1, w2, ..., wD] — вектор весов (коэффициентов) модели,
• w0 — смещение,
• ϵ ∼ N(0, α^−1) — случайный гауссовский шум измерения с точностью α (alpha) (обратной дисперсией α = 1 / σ^2)

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

Байесовский подход кардинально меняет взгляд на природу параметров. В байесовской регрессии коэффициенты wj перестают быть константами и рассматриваются как случайные величины. До получения данных мы задаем наши априорные представления о параметрах регрессии. В частности, в гребневой байесовской регрессии (Bayesian Ridge) на все веса накладывается одно и то же априорное нормальное распределение с общим параметром регуляризации λ(лямбда):

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

Главная идея ARD-регрессии заключается в том,  что на каждый вес wj накладывается индивидуальное априорное нормальное распределение со своим собственным параметром точности  λj:

Каждый гиперпараметр λj выступает в роли индивидуального регулятора силы соответствующего признака:

• Если признак xj информативен, то в процессе обучения модель уменьшает значение λj → 0. Это расширяет априорное распределение, делая его плоским. Высокая дисперсия разрешает весу wj свободно отклоняться от нуля и принимать значения, необходимые для точного объяснения целевой переменной y.
• Если признак xj не оказывает влияния на прогноз целевой переменной, тогда модель увеличивает параметр λj → ∞. Дисперсия 1/λj стремится к нулю, и распределение параметра wj схлопывается в точке 0. Это принудительно зануляет вес wj, исключая признак j из финального уравнения регрессии.

В итоге оптимизация вектора λ превращается в автоматический процесс прореживания пространства признаков.

Обучение модели: EM-алгоритм и оценка параметров

Обучение ARD-регрессии строится на EM-подобном алгоритме (Expectation-Maximization), где вектор весов w играет роль непрерывных скрытых переменных. В отличие от классического EM, здесь на M-шаге производится прямая максимизация маргинального логарифма правдоподобия (Log Evidence).

Алгоритм состоит из двух чередующихся шагов: 

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

На E-шаге мы фиксируем текущие значения гиперпараметров шума α и индивидуальных точностей λ. Наша цель — найти апостериорное распределение весов p(w∣y, X, α, λ). Так как априорное распределение весов и функция правдоподобия являются гауссовскими, апостериорное распределение также будет нормальным N(μw, Σ). Благодаря этому вычисление апостериора сводится к вычислению ковариационной матрицы и вектора средних.

• Ковариационная матрица весов Σ (степень неопределенности коэффициентов):

• Математическое ожидание весов μw (текущие точечные значения коэффициентов):

Практический нюанс реализации: Если число признаков D больше числа наблюдений N, прямое обращение матрицы ковариации Σ размером D×D становится вычислительно дорогим. В коде класса для этого случая используется тождество Вудбери (Woodbury Matrix Identity), позволяющее сводить операцию к обращению матрицы размером N×N.

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

На M-шаге найденные апостериорные характеристики из E-шага используются для обновления гиперпараметров λ и α. В классическом EM-алгоритме M-шаг максимизирует математическое ожидание полного логарифма правдоподобия. В случае ARD-регрессии обновление гиперпараметров производится с помощью прямой максимизации маргинального логарифма правдоподобия (evidence maximization). Этот подход эквивалентен M-шагу. Введение коэффициента эффективности признаков γⱼ ускоряет сходимость по сравнению со стандартным EM-алгоритмом. Чтобы предотвратить численную нестабильность и переобучение самих гиперпараметров, над ними задаются слабоинформативные априорные Гамма-распределения с параметрами (α1,α2) и (λ1,λ2).

• Обновление точностей весов λj:

где, γj  = 1 − λjΣjj — мера эффективности j-го признака.

• Обновление точности шума α:

где, SSE = ||y − Xμw||^2 — сумма квадратов остатков регрессии (Sum of Squared Errors).

• Отбор признаков (Pruning) и критерий остановки обучения

После каждого M-шага выполняется отбор признаков. Если гиперпараметр λj превышает заданный порог threshold_lambda, признак признается нерелевантным:

• Коэффициент wj принудительно обнуляется (wj = 0),
• Признак j исключается из матрицы X на следующих итерациях  для ускорения расчетов.

Итерационный процесс продолжается до тех пор, пока изменение вектора весов || w_old  − w_new||  не станет меньше заданного значения tol, либо пока не будет достигнуто максимальное число итераций max_iter.

Практическая реализация: класс ARDRegression

Вся математика байесовского вывода — от итерационного EM-алгоритма до динамического отбора признаков  — объединена в класс ARDRegression.

//+------------------------------------------------------------------+
//| Класс ARDRegression                                             |
//+------------------------------------------------------------------+
class ARDRegression
  {
private:
   int               m_max_iter;            // Максимальное количество итераций EM-алгоритма
   double            m_tol;                 // Порог сходимости по изменению весов
   double            m_alpha_1,m_alpha_2;   // Параметры Гамма-распределения для точности шума alpha(α)
   double            m_lambda_1,m_lambda_2; // Параметры Гамма-распределения для точностей весов lambda(λ)
   double            m_threshold_lambda;    // Порог отбора (если λ_j > threshold -> вес = 0)
   bool              m_fit_intercept;       // Нужно ли обучать параметр смещения (w0)
   bool              m_compute_score;       // Нужно ли вычислять целевую функцию (Log Evidence/Marginal Log-likelihood)
   bool              m_verbose;             // Вывод сообщений о сходимости алгоритма обучения в журнал

   bool              UpdateSigmaAndCoef(const matrix &X_cent, const vector &y_cent, bool &keep_lambda[]);

public:
   vector            coef_;      // Вектор итоговых коэффициентов (w_j)
   double            intercept_; // Итоговое смещение (w_0)
   double            alpha_;     // Итоговая точность шума (1 / sigma^2) (Обратная дисперсия шума σ)
   vector            lambda_;    // Вектор точностей признаков
   matrix            sigma_;     // Ковариационная матрица апостериорного распределения весов
   int               n_iter_;    // Число итераций до сходимости
   double            scores_[];  // История значений Evidence/Marginal Log-likelihood на каждой итерации
   vector            X_offset_;  // Средние значения столбцов матрицы X (для центрирования признаков)
   double            y_offset_;  // Среднее значение целевой переменной y (для центрирования отклика)

                     ARDRegression(int      max_iter = 300,
                 double   tol = 1e-3,
                 double   alpha_1 = 1e-6,
                 double   alpha_2 = 1e-6,
                 double   lambda_1 = 1e-6,
                 double   lambda_2 = 1e-6,
                 double   threshold_lambda = 1e4,
                 bool     fit_intercept = true,
                 bool     compute_score = false,
                 bool     verbose = false);

   bool              fit(const matrix &X, const vector &y);
   vector            predict(const matrix &X, vector &y_std);
  };

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

• max_iter (по умолчанию 300): Максимальное количество итераций EM-алгоритма,
• tol (по умолчанию  0.001): Критерий остановки обучения,

• Гипераприоры Гамма-распределения (alpha_1, alpha_2, lambda_1, lambda_2)

Эти параметры задают априорные значения точности шума α и точностей весов λj через Гамма-распределения Gamma(α1,α2) и Gamma(λ1,λ2). В 99% практических задач оставляйте эти параметры равными 1e-6 как они заданы по умолчанию. Эти значения задают так называемые неинформативные (плоские) априорные распределения, предоставляя возможность данным говорить самим за себя.

• threshold_lambda (по умолчанию 10000): параметр для механизма отбора признаков. Уменьшение (до 100 – 1000) делает отбор признаков более агрессивным. Модель будет сильнее прореживать веса, оставляя только самые сильные сигналы,

• fit_intercept (по умолчанию true): Вычисление смещения. Всегда оставляйте true, если ваши входные данные X не центрированы. Если вы предварительно отцентрировали данные (mean=0), можно поставить false,

• compute_score (по умолчанию false): Расчет целевой функции (Marginal Log-likelihood / Log Evidence) на каждой итерации. Поставьте true, если вы хотите проанализировать график сходимости модели. Результаты вычислений сохраняются в массив scores_,
• verbose (по умолчанию false): Вывод детального лога обучения в журнал.

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

1. Метод обучения: fit(const matrix &X, const vector &y)

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

//+------------------------------------------------------------------+
//| Метод обучения модели                                            |
//+------------------------------------------------------------------+
bool ARDRegression::fit(const matrix &X, const vector &y)
  {
   ulong n_samples  = X.Rows();
   ulong n_features = X.Cols();

   matrix X_cent = X;
   vector y_cent = y;

//--- Центрирование данных
   X_offset_.Resize(n_features);
   if(m_fit_intercept)
     {
      y_offset_ = y.Mean();
      y_cent = y - y_offset_;

      for(ulong j = 0; j < n_features; j++)
        {
         X_offset_[j] = X.Col(j).Mean();
         X_cent.Col(X_cent.Col(j) - X_offset_[j], j);
        }
     }
   else
     {
      X_offset_.Fill(0.0);
      y_offset_ = 0.0;
     }

   coef_ = vector::Zeros(n_features);

//--- Инициализация параметров
   double eps = 2.2204460492503131e-16;
   alpha_ = 1.0 / (y_cent.Var() + eps);
   lambda_ = vector::Ones(n_features);

   bool keep_lambda[];
   ArrayResize(keep_lambda, (int)n_features);
   ArrayFill(keep_lambda, 0, (int)n_features, true);

   vector coef_old_ = vector::Zeros(n_features);
   ArrayFree(scores_);

   for(int iter_ = 0; iter_ < m_max_iter; iter_++)
     {
      //--- 1. Подсчет активных признаков
      ulong keep_count = 0;
      for(ulong j = 0; j < n_features; j++)
         if(keep_lambda[j])
            keep_count++;

      if(keep_count == 0)
        {
         coef_.Fill(0.0);
         sigma_.Resize(0, 0);
         n_iter_ = iter_ + 1;
         break;
        }

      //--- 2. Обновление sigma_ и coef_  (E-шаг)
      if(!UpdateSigmaAndCoef(X_cent, y_cent, keep_lambda))
        {
         PrintFormat("ARD Error: Aborted at iteration %d due to numerical solver failure.", iter_);
         return false;
        }

      //--- 3. Обновление гиперпараметров (M-шаг)
      double sse_ = MathPow(y_cent - X_cent @ coef_, 2).Sum();

      vector lambda_keep = SelectByMask(lambda_, keep_lambda);
      vector gamma_ = 1.0 - lambda_keep * sigma_.Diag();

      vector coef_select = SelectByMask(coef_, keep_lambda);
      vector lambda_keep_new = (gamma_ + 2.0 * m_lambda_1) / (coef_select * coef_select + 2.0 * m_lambda_2);
      lambda_keep_new.Clip(1e-10, 1e+10); // Ограничиваем сверху и снизу от переполнения
      ScatterByMask(lambda_, keep_lambda, lambda_keep_new);

      //--- Числитель и знаменатель
      double num = n_samples - gamma_.Sum() + 2.0 * m_alpha_1;
      double den = sse_ + 2.0 * m_alpha_2;

      //--- Защита от деления на 0
      if(den < 1e-12)
         alpha_ = 1e+10;
      else
        {
         alpha_ = num / den;
         // Ограничиваем диапазон [1e-10, 1e+10]
         alpha_ = MathMax(1e-10, MathMin(1e+10, alpha_));
        }
      //  alpha_ = (n_samples - gamma_.Sum() + 2.0 * m_alpha_1) / (sse_ + 2.0 * m_alpha_2);

      //--- 4.Отбор признаков (Pruning)
      keep_count = 0;
      for(ulong j = 0; j < n_features; j++)
        {
         keep_lambda[j] = (lambda_[j] < m_threshold_lambda);
         if(!keep_lambda[j])
            coef_[j] = 0.0;
         else
            keep_count++;
        }

      //--- 4.1. Расчет целевой функции  (Log Evidence / Objective function)
      if(m_compute_score || m_verbose)
        {
         double s = (m_lambda_1 * MathLog(lambda_) - m_lambda_2 * lambda_).Sum();
         s += m_alpha_1 * MathLog(alpha_) - m_alpha_2 * alpha_;

         int sign;
         double logdet_sigma;
         double Signed_Determinant = sigma_.SLogDet(sign,logdet_sigma);
         // Проверяем корректность знака детерминанта (SPD матрица должна иметь sign > 0)
         if(sign <= 0)
           {
            if(m_verbose)
               PrintFormat("ARD Warning: non-positive determinant for sigma_ at iteration %d (sign=%d). Score calculation skipped.", iter_, sign);
           }
         else
           {
            s += 0.5 * (logdet_sigma + (double)n_samples * MathLog(alpha_) + MathLog(lambda_).Sum());
            s -= 0.5 * (alpha_ * sse_ + (lambda_ * coef_ * coef_).Sum());
            int score_size = ArraySize(scores_);
            ArrayResize(scores_, score_size + 1);
            scores_[score_size] = s;
           }
        }
      //--- 5. Проверка сходимости
      if(iter_ > 0)
        {
         if(MathAbs(coef_old_ - coef_).Sum() < m_tol)
           {
            if(m_verbose)
               PrintFormat("Converged after %d iterations", iter_);
            n_iter_ = iter_ + 1;
            break;
           }
        }

      coef_old_ = coef_;

      if(keep_count == 0)
        {
         n_iter_ = iter_ + 1;
         break;
        }

      n_iter_ = iter_ + 1;
     }

//--- 6. Финальный пересчет (с самыми последними alpha_ и lambda_)
   if(!UpdateSigmaAndCoef(X_cent, y_cent, keep_lambda))
     {
      Print("ARD Error: Numerical solver failure on final step");
      return false;
     }

//--- Сохранение итогового смещения
   if(m_fit_intercept)
      intercept_ = y_offset_ - X_offset_.Dot(coef_);
   else
      intercept_ = 0.0;

   return true;
  }

//+------------------------------------------------------------------+
//| Метод для обновления sigma_ и coef_ по текущим alpha_ и lambda_  |
//+------------------------------------------------------------------+
bool ARDRegression::UpdateSigmaAndCoef(const matrix &X_cent, const vector &y_cent, bool &keep_lambda[])
  {
   ulong n_samples  = X_cent.Rows();
   ulong n_features = X_cent.Cols();

//--- 1. Формируем X_keep и lambda_keep через функции-маски
   matrix X_keep      = SelectColsByMask(X_cent, keep_lambda);
   vector lambda_keep = SelectByMask(lambda_, keep_lambda);

   ulong keep_count = X_keep.Cols();

   if(keep_count == 0)
     {
      coef_.Fill(0.0);
      sigma_.Resize(0, 0);
      PrintFormat("ARD Warning: All features pruned (keep_count == 0) at iteration %d. Model cannot be fitted.", n_iter_);
      return true;
     }

//--- 2. Расчет апостериорной ковариационной матрицы весов sigma_ (Direct vs Woodbury)
   if(n_samples >= n_features) // N > D
     {
      matrix gram = X_keep.Transpose().MatMul(X_keep);
      matrix sigma_inv = matrix::Identity(keep_count, keep_count);
      sigma_inv.Diag(lambda_keep);
      sigma_inv = sigma_inv + (gram * alpha_);
      //sigma_ = sigma_inv.PInv();

      matrix residuals;
      vector singular_values;
      long rank;
      //--- Решаем A * X = I через GELSD (SVD Divide & Conquer)
      if(!sigma_inv.LeastSquaresSolutionDC(matrix::Eye(keep_count, keep_count), -1.0, sigma_, residuals, singular_values, rank))
        {return false;}
     }
   else // D > N
     {
      //---  Формируем вектор/матрицу inverse lambda
      vector inv_lambda(keep_count);
      for(ulong j = 0; j < keep_count; j++)
         inv_lambda[j] = 1.0 / lambda_keep[j];

      //---  X_keep_inv = X_keep * inv_lambda (поэлементное умножение каждого столбца j на inv_lambda[j])
      matrix X_keep_inv(n_samples, keep_count);
      for(ulong i = 0; i < n_samples; i++)
        {
         for(ulong j = 0; j < keep_count; j++)
            X_keep_inv[i][j] = X_keep[i][j] * inv_lambda[j];
        }

      //---  Вычисляем инвертируемую матрицу размером (n_samples, n_samples)
      matrix eye_samples = matrix::Identity(n_samples, n_samples);
      matrix inner_mat = (eye_samples / alpha_) + X_keep_inv.MatMul(X_keep.Transpose());
      //matrix inner_inv = inner_mat.PInv();

      matrix residuals;
      vector singular_values;
      long rank;
      matrix inner_inv;
      if(!inner_mat.LeastSquaresSolutionDC(matrix::Eye(n_samples, n_samples), -1.0, inner_inv, residuals, singular_values, rank))
        {return false;}
      //--- inner_inv @ X_keep_inv -> размер (n_samples, keep_count)
      matrix temp = inner_inv.MatMul(X_keep_inv);

      //--- inv_X_T = inv_lambda_col * X_keep.T -> размер (keep_count, n_samples)
      matrix inv_X_T(keep_count, n_samples);
      matrix X_T = X_keep.Transpose();
      for(ulong i = 0; i < keep_count; i++)
        {
         for(ulong j = 0; j < n_samples; j++)
            inv_X_T[i][j] = X_T[i][j] * inv_lambda[i];
        }

      //--- sigma_ = - (inv_X_T @ temp) -> размер (keep_count, keep_count)
      sigma_ = inv_X_T.MatMul(temp) * (-1.0);

      //--- Добавление диагонали 1.0 / lambda
      for(ulong j = 0; j < keep_count; j++)
         sigma_[j][j] += inv_lambda[j];
     }

//--- 3. Расчет апостериорного среднего весов coef_
   vector temp_dot = X_keep.Transpose().MatMul(y_cent);
   vector coef_keep = temp_dot.MatMul(sigma_) * alpha_;

   coef_.Fill(0.0);
   ScatterByMask(coef_, keep_lambda, coef_keep);
   
   return true;
  }

2. Метод прогнозирования: predict(const matrix &X, vector &y_std)

В отличие от стандартных моделей машинного обучения, метод predict возвращает не просто точечные значения, а полноценный вероятностный прогноз: 

• Точечный прогноз (Математическое ожидание y*):

• Оценка неопределенности (стандартное отклонение ystd). Апостериорная дисперсия прогноза учитывает как шум в данных α^−1, так и неопределенность оценки весов Σ:

//+------------------------------------------------------------------+
//| Метод предсказания                                               |
//+------------------------------------------------------------------+
vector ARDRegression::predict(const matrix &X, vector &y_std)
  {
   //--- проверка совместимости размерностей признаков
   if(X.Cols() != coef_.Size())
     {
      PrintFormat("ARD Error [predict]: Dimension mismatch! Model expects %d features, but X has %d cols.",
                  coef_.Size(), X.Cols());
      vector empty(0);
      return empty;
     }
// --- прогноз условного среднего
   vector y_mean = X.MatMul(coef_);
   if(m_fit_intercept)
      y_mean = y_mean + intercept_;

//--- расчет стандартного отклонения
   ulong n_samples  = X.Rows();
   ulong n_features = X.Cols();

//--- Подсчет активных признаков
   bool keep_lambda[];
   ArrayResize(keep_lambda, (int)n_features);
   ulong keep_count = 0;

   for(ulong j = 0; j < n_features; j++)
     {
      keep_lambda[j] = (lambda_[j] < m_threshold_lambda);
      if(keep_lambda[j])
         keep_count++;
     }

//--- если активных признаков нет — возвращаем шум 1 / alpha
   if(keep_count == 0)
     {
      y_std = vector::Ones(n_samples) * MathSqrt(1.0 / alpha_);
      return y_mean;
     }

   matrix X_keep = SelectColsByMask(X, keep_lambda);
   matrix X_sigma = X_keep.MatMul(sigma_);
   y_std = MathSqrt((X_sigma * X_keep).Sum(1) + (1.0 / alpha_));

   return y_mean;
  }

Сравнение результатов обучения с моделью Scikit-Learn

Для проверки математической корректности, точности и численной устойчивости разработанного класса ARDRegression было проведено сравнение с эталонной реализацией из библиотеки scikit-learn (Python).

Сравнение выполнялось на двух тестовых скриптах:

• scikit_learn.py — генерирует синтетический датасет, обучает модель sklearn.linear_model.ARDRegression и выводит результаты вычислений в журнал,
• TestARD.mq5 — MQL5 скрипт, импортирует те же данные (CSV-файлы X_ARD.csv и y_ARD.csv) и обучает нашу реализацию.

Обе модели обучались на идентичных наборах данных в двух режимах:

• Датасет 1 (N > D): 200 наблюдений, 50 признаков. Стандартный режим обучения,
• Датасет 2 (N < D): 40 образцов, 50 признаков. Усеченная выборка для проверки корректности вычислений через тождество Вудбери.

Истинные параметры модели:

• смещение (intercept): w0 = 5.0,
• значимые веса: w1 = 3.5, w2 = −2.0, w3 = 1.2,
• оставшиеся 47 весов = 0 (незначимые признаки),
• шум наблюдений: σ = 0.5 (α = 4.0).

//+------------------------------------------------------------------+
//|                                                      TestARD.mq5 |
//|                                                           Eugene |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Eugene"
#property link      "https://www.mql5.com"
#property version   "1.00"

#include <ARDRegression\ARD.mqh>

//+------------------------------------------------------------------+
//| функция для запуска теста и вывода результатов                   |
//+------------------------------------------------------------------+
void RunARDTest(const matrix &X, const vector &y, string test_name)
  {
   ulong n_samples  = X.Rows();
   ulong n_features = X.Cols();

//--- Создаем модель
   ARDRegression model(300, 1e-3, 1e-6, 1e-6, 1e-6, 1e-6, 1e4, true, true, false);

   uint start_time = GetTickCount();
   model.fit(X, y); // обучение
   uint elapsed_time = GetTickCount() - start_time;

// --- Предсказание и проверка ошибки
   vector y_std;
   vector y_pred = model.predict(X, y_std);
   vector errors = y - y_pred;
   double mse = y_pred.RegressionMetric(y,REGRESSION_MSE);

// --- Вывод результатов в журнал
   PrintFormat("=========================================================");
   PrintFormat(" %s (n_samples=%d, n_features=%d)", test_name, n_samples, n_features);
   PrintFormat("=========================================================");
   PrintFormat("Время выполнения:      %u ms", elapsed_time);
   PrintFormat("Число итераций:        %d", model.n_iter_);
   PrintFormat("Intercept (Истинный 5): %.5f", model.intercept_);
   PrintFormat("Alpha (Точность шума):  %.5f", model.alpha_);
   PrintFormat("MSE предсказания:       %.5f", mse);

   int n_scores = ArraySize(model.scores_);
   if(n_scores > 0)
     {
      PrintFormat("Log Evidence (Marginal Likelihood): %.5f", model.scores_[n_scores - 1]);

      // --- Динамика Log Evidence
      PrintFormat("\n--- Динамика Log Evidence по итерациям ---");
      if(n_scores <= 10)
        {
         for(int it = 0; it < n_scores; it++)
            PrintFormat("Итерация %02d | Log Evidence: %12.5f", it + 1, model.scores_[it]);
        }
      else
        {
         //--- Показываем первые 5
         for(int it = 0; it < 5; it++)
            PrintFormat("Итерация %02d | Log Evidence: %12.5f", it + 1, model.scores_[it]);

         PrintFormat("  ...");

         //--- И последние 3
         for(int it = n_scores - 3; it < n_scores; it++)
            PrintFormat("Итерация %02d | Log Evidence: %12.5f", it + 1, model.scores_[it]);
        }
     }

   PrintFormat("\n--- Стандартное отклонение предсказаний (y_std, первые 5 сэмплов) ---");
   for(ulong i = 0; i < 5 && i < n_samples; i++)
     {
      PrintFormat("Сэмпл %02d | y_pred: %8.5f | y_std: %8.5f", i, y_pred[i], y_std[i]);
     }

   PrintFormat("\n--- Веса coef_ (Истинные: X00=3.5, X01=-2.0, X02=1.2, остальные=0) ---");
   for(ulong i = 0; i < 5 && i < n_features; i++)
     {
      PrintFormat("X%02d | Оценка: %8.5f | Lambda: %10.2f",
                  i, model.coef_[i], model.lambda_[i]);
     }

//---  Pruning
   int active_features = 0;
   for(ulong j = 0; j < n_features; j++)
     {
      if(MathAbs(model.coef_[j]) > 1e-4)
         active_features++;
     }
   PrintFormat("\nОсталось активных признаков с весом > 0: %d из %d\n", active_features, n_features);
  }

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   matrix X_csv;
   CSVtoMatrix("X_ARD.csv", X_csv);
   matrix y_csv;
   CSVtoMatrix("y_ARD.csv", y_csv);
   vector y_full = y_csv.Col(0);

// ТЕСТ 1: Полный датасет (200 x 50) -> Direct Gram Matrix
   RunARDTest(X_csv, y_full, "ТЕСТ 1: FULL DATASET (Direct Gram)");

//--- ТЕСТ 2: Первые 40 элементов  -> Woodbury Identity
   matrix X_sub = X_csv;
   vector y_sub = y_full;

//--- урезаем данные до 40 элементов
   X_sub.Resize(40, 50);
   y_sub.Resize(40);

   RunARDTest(X_sub, y_sub, "ТЕСТ 2: SUBSET DATASET (Woodbury Identity)");
  }

import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import ARDRegression

# Функция для вывода результатов работы модели Scikit-Learn
def RunARDTest(title, model, X, y, elapsed_ms, show_score_history=True):
    y_pred, y_std = model.predict(X, return_std=True)
    mse = np.mean((y - y_pred) ** 2)
    active_count = np.sum(model.coef_ != 0.0)
    n_samples, n_features = X.shape

    print("=========================================================")
    print(f" {title} (n_samples={n_samples}, n_features={n_features})")
    print("=========================================================")
    print(f"Время выполнения:       {elapsed_ms} ms")
    print(f"Число итераций:         {model.n_iter_}")
    print(f"Intercept (Истинный 5): {model.intercept_:.5f}")
    print(f"Alpha (Точность шума):  {model.alpha_:.5f}")
    print(f"MSE предсказания:       {mse:.5f}")
    print(f"Log Evidence (Marginal Likelihood): {model.scores_[-1]:.5f}")

    # Вывод динамики целевой функции по итерациям
    if show_score_history:
        scores = model.scores_
        n_scores = len(scores)
        print("\n--- Динамика Log Evidence по итерациям ---")
        if n_scores <= 10:
            for it, sc in enumerate(scores):
                print(f"Итерация {it+1:02d} | Log Evidence: {sc:12.5f}")
        else:
            # Показываем первые 5 и последние 3 итерации
            for it in range(5):
                print(f"Итерация {it+1:02d} | Log Evidence: {scores[it]:12.5f}")
            print("  ...")
            for it in range(n_scores - 3, n_scores):
                print(f"Итерация {it+1:02d} | Log Evidence: {scores[it]:12.5f}")

    print("\n--- Стандартное отклонение предсказаний (y_std, первые 5 сэмплов) ---")
    for i in range(min(5, n_samples)):
        print(f"Сэмпл {i:02d} | y_pred: {y_pred[i]:8.5f} | y_std: {y_std[i]:8.5f}")

    print("\n--- Веса coef_ (Истинные: X00=3.5, X01=-2.0, X02=1.2, остальные=0) ---")
    for j in range(min(5, n_features)):
        print(f"X{j:02d} | Оценка: {model.coef_[j]:8.5f} | Lambda: {model.lambda_[j]:11.2f}")

    print(f"\nОсталось активных признаков с весом > 0: {active_count} из {n_features}\n")

# 1. ГЕНЕРАЦИЯ ДАТАСЕТА (200 x 50)
np.random.seed(42)

n_samples = 200
n_features = 50

X = np.random.randn(n_samples, n_features)

true_w = np.zeros(n_features)
true_w[0] = 3.5
true_w[1] = -2.0
true_w[2] = 1.2

# Вычисляем истинный сигнал 
y_true = np.dot(X, true_w) + 5.0

noise = np.random.normal(0, 0.5, size=n_samples)
y = y_true + noise

# ТЕСТ 1: FULL DATASET (200 x 50) -> N > D
t0 = time.perf_counter()
model = ARDRegression(compute_score=True).fit(X, y)
elapsed = int((time.perf_counter() - t0) * 1000)
RunARDTest("ТЕСТ 1 [SKLEARN]: FULL DATASET (200x50)", model, X, y, elapsed)

# ТЕСТ 2: SUBSET DATASET (40 x 50) -> N < K (Срез первых 40 строк)
X_sub = X[:40, :].copy()
y_sub = y[:40].copy()
t0 = time.perf_counter()
model_sub = ARDRegression(compute_score=True).fit(X_sub, y_sub)
elapsed_sub = int((time.perf_counter() - t0) * 1000)
RunARDTest("ТЕСТ 2 [SKLEARN]: SUBSET DATASET (40x50)", model_sub, X_sub, y_sub, elapsed_sub)

# FIGURE 1: LOG EVIDENCE LEARNING CURVES

fig1, axes1 = plt.subplots(1, 2, figsize=(10, 5))

# --- Plot 1: Full Dataset (200x50) ---
iters1 = np.arange(1, len(model.scores_) + 1)
axes1[0].plot(iters1, model.scores_, 'b-o', linewidth=2, markersize=4, label='Log Evidence')
axes1[0].set_title('Learning Curve (Test 1: 200x50, N > D)', fontsize=11, fontweight='bold')
axes1[0].set_xlabel('EM Iteration')
axes1[0].set_ylabel('Log Evidence')
axes1[0].grid(True, linestyle=':', alpha=0.6)
axes1[0].legend()

# --- Plot 2: Subset Dataset (40x50) ---
iters2 = np.arange(1, len(model_sub.scores_) + 1)
axes1[1].plot(iters2, model_sub.scores_, 'r-s', linewidth=2, markersize=4, label='Log Evidence')
axes1[1].set_title('Learning Curve (Test 2: 40x50, N < D)', fontsize=11, fontweight='bold')
axes1[1].set_xlabel('EM Iteration')
axes1[1].set_ylabel('Log Evidence')
axes1[1].grid(True, linestyle=':', alpha=0.6)
axes1[1].legend()
plt.tight_layout()

# FIGURE 2: WEIGHT ESTIMATES vs TRUE WEIGHTS

fig2, axes2 = plt.subplots(1, 2, figsize=(10, 5))
x_features = np.arange(n_features)

# --- Plot 1: Weight Estimates (200x50) ---
axes2[0].stem(x_features, true_w, linefmt='b-', markerfmt='bo', basefmt='k-', label='True Weights')
axes2[0].stem(x_features, model.coef_, linefmt='r--', markerfmt='rx', basefmt='k-', label='ARD Estimates')
axes2[0].set_title('Weight Estimates (Test 1: 200x50, N > D)', fontsize=11, fontweight='bold')
axes2[0].set_xlabel('Feature Index')
axes2[0].set_ylabel('Weight Value')
axes2[0].grid(True, linestyle=':', alpha=0.6)
axes2[0].legend()

# --- Plot 2: Weight Estimates (40x50) ---
axes2[1].stem(x_features, true_w, linefmt='b-', markerfmt='bo', basefmt='k-', label='True Weights')
axes2[1].stem(x_features, model_sub.coef_, linefmt='r--', markerfmt='rx', basefmt='k-', label='ARD Estimates')
axes2[1].set_title('Weight Estimates (Test 2: 40x50, N < D)', fontsize=11, fontweight='bold')
axes2[1].set_xlabel('Feature Index')
axes2[1].set_ylabel('Weight Value')
axes2[1].grid(True, linestyle=':', alpha=0.6)
axes2[1].legend()
plt.tight_layout()
plt.show()

Сопоставление показало полное совпадение результатов MQL5-реализации и scikit-learn с точностью до 5-го знака после запятой по всем ключевым метрикам: оценкам весов, апостериорной точности alpha, маргинальному правдоподобию (Log Evidence) и неопределенности предсказаний.

На рис. 1 представлена динамика логарифма маргинального правдоподобия (Log Evidence) по итерациям EM-алгоритма.

рис. 1 Кривая обучения ARD-регрессии 

График наглядно иллюстрирует быструю сходимость модели: всего 10–15 итераций. 

На рис. 2 представлено сравнение восстановленных весов w (красный пунктир) с истинными значениями весов wtrue (синяя линия), заданными при генерации данных.

рис. 2 Оценки весов ARD-регрессии

Модель ARD верно определила 3 информативных признака, достаточно точно оценив коэффициенты как при полноразмерной выборке (200×50), так и в условиях дефицита данных (40×50). 

Генерация признаков (Feature Engineering) и преодоление проблемы нелинейности

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

Существует три основных подхода к решению этой задачи.

• Явное преобразование признаков

Самый простой способ — сформировать новые признаки вручную с помощью алгебраических и тригонометрических преобразований:

• полиномиальные степени (x^2, x^3),
• нелинейные преобразования (∣x∣, ln(x), exp(x)),
• гармоники (sin(x), cos(x)),
• парные произведения (xi * xj) и т.д.

Главный недостаток подхода — проклятие размерности. При росте количества исходных признаков (D) число новых признаков растёт очень быстро. Например, при D=50 полиномиальное расширение второй степени уже даёт более 1300 комбинаций, а третьей степени — десятки тысяч признаков.

• Ядерный трюк

В классическом машинном обучении проблему проклятия размерности решают с помощью ядерного трюка (Kernel Trick). В таких моделях, как Kernel Ridge Regression, метод опорных векторов (SVM) или Гауссовские процессы (GP), скалярное произведение признаков в высокоразмерном пространстве заменяется функцией ядра k(x, x′).

Самым популярным и универсальным является RBF-ядро (Radial Basis Function):

где,

• x, x′ ∈ RD — два вектора признаков,
• ∥ x − x′ ∥^2 — квадрат евклидова расстояния между ними,
• γ > 0 — параметр масштаба (гамма), определяющий ширину ядра и чувствительность к расстоянию между точками.

Однако у классических ядерных методов есть серьезный недостаток — они требуют вычисления и хранения матрицы Грама NxN и её обращения, что требует O(N^3) операций. Это сильно ограничивает применимость метода на больших выборках.

• Метод случайных признаков (Random Fourier Features)

Эффективное решение предложил в 2007 году алгоритм Random Fourier Features (Rahimi & Recht). Вместо использования ядерной модели мы явно генерируем новые признаки, которые хорошо аппроксимируют RBF-ядро, и передаем их в обычную линейную модель.

Аппроксимация RBF-ядра выполняется следующим образом:

• Генерация случайных частот — семплируем M случайных векторов из нормального распределения, соответствующего RBF ядру N(0, 2γ),
• Генерация случайных фаз — семплируем M случайных смещений из равномерного распределения Uniform(0, 2*pi),
• Преобразование данных — для любой входной точки признаков x размерности D строится трансформированный вектор признаков φ(x) размерности M:

где,

• x∈RD — вектор исходных признаков,
• M — количество сгенерированных случайных признаков,
• wj — веса частот,
• bj — случайные фазы,
• Sqrt(2/M) — нормировочный коэффициент.

Преобразование cos(Wx + b) структурно похоже на формулу нейрона многослойного персептрона. Только вместо обучения весов W и смещений b через обратное распространение ошибки, мы однократно генерируем их из соответствующих распределений и фиксируем. Дальше обучение проходит на этих новых признаках с помощью любой линейной модели.

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

Реализация генератора случайных признаков: класс CRBFSampler

//+------------------------------------------------------------------+
//| Класс генерации случайных признаков (Random Fourier Features)    |
//| для RBF ядра                                                     |
//+------------------------------------------------------------------+
class CRBFSampler
  {
private:
   double            m_gamma;           // Параметр ширины ядра (gamma)
   ulong             m_n_components;    // Количество генерируемых компонентов (M)
   int               m_seed;            // Seed (-1: случайная генерация, >=0: фиксированная)

   matrix            m_random_weights;  // Матрица весов W: (n_features x n_components)
   vector            m_random_offsets;  // Вектор смещений b: (n_components)
   bool              m_is_fitted;       // Флаг успешности генерации весов
public:
                     CRBFSampler(const double gamma=1.0,
               const ulong n_components=100,
               const int seed=-1);
                    ~CRBFSampler(void) {}

   bool              Fit(const matrix &X);
   matrix            Transform(const matrix &X);
   matrix            FitTransform(const matrix &X);
   bool              IsFitted(void) const { return m_is_fitted; }
  };

//+------------------------------------------------------------------+
//| Конструктор                                                      |
//+------------------------------------------------------------------+
CRBFSampler::CRBFSampler(const double gamma=1.0,
                         const ulong n_components=100,
                         const int seed=-1)
   : m_gamma(gamma),
     m_n_components(n_components),
     m_seed(seed),
     m_is_fitted(false)
  {
// --- Инициализация генератора ПСЧ
   if(m_seed != -1)
      MathSrand((uint)m_seed);
   else
      MathSrand(GetTickCount());
  }

//+------------------------------------------------------------------+
//| Генерация случайных признаков                                    |
//+------------------------------------------------------------------+
bool CRBFSampler::Fit(const matrix &X)
  {
   ulong n_features = X.Cols();

//--- Генерация весов W ~ N(0, 2*gamma) с помощью класса CMVN
   vector Mu = vector::Zeros(m_n_components); // Вектор нулевых средних
   matrix Sigma = matrix::Identity(m_n_components, m_n_components);
   Sigma = Sigma * 2.0 * m_gamma;

// Генерируем матрицу нормальных случайных величин размера [n_features x n_components]
   if(!CMVN::rnd(Mu, Sigma, (int)n_features, m_random_weights))
     {
      Print("Ошибка CRBFSampler::Fit: Не удалось сгенерировать веса через CMVN::rnd");
      return false;
     }

//--- Генерация смещений b ~ Uniform(0, 2*PI)
   m_random_offsets = vector::Random(m_n_components, 0.0, 2.0 * M_PI);

   m_is_fitted = true;
   return true;
  }

//+------------------------------------------------------------------+
//| Преобразование матрицы признаков                                 |
//+------------------------------------------------------------------+
matrix CRBFSampler::Transform(const matrix &X)
  {
   if(!m_is_fitted)
     {
      Print("Ошибка CRBFSampler::Transform: Модель еще не обучена (вызовите Fit)!");
      matrix empty(0,0);
      return empty;
     }

   ulong rows = X.Rows();

//--- 1. Проекция: Projection = X @ W (Размерность: N x M)
   matrix projection = X @ m_random_weights;

//--- 2. Добавление смещения b (векторное сложение по строкам)
   for(ulong i = 0; i < rows; i++)
     {
      projection.Row(MathCos(projection.Row(i) + m_random_offsets), i);
     }

//--- 3. Масштабирование: sqrt(2.0 / n_components)
   double scale = MathSqrt(2.0 / (double)m_n_components);
   projection = projection * scale;

   return projection;
  }

//+------------------------------------------------------------------+
//| Объединенный метод Fit + Transform                               |
//+------------------------------------------------------------------+
matrix CRBFSampler::FitTransform(const matrix &X)
  {
   if(!Fit(X))
     {
      matrix empty(0,0);
      return empty;
     }
   return Transform(X);
  }

Для работы с алгоритмом Random Fourier Features был написан класс CRBFSampler. Он отвечает за генерацию случайных частот и фаз и преобразование исходной матрицы признаков из D-мерного пространства в M-мерное.

В конструкторе задаются три основных параметра:

• n_components (M) количество случайных признаков Фурье. Чем больше M, тем точнее аппроксимация RBF-ядра,
• gamma (γ) — параметр ширины ядра,
• seed — зерно генератора случайных чисел (для воспроизводимости экспериментов).

При малом gamma (< 3) модель строит очень плавную, гладкую кривую. Она видит только общую тенденцию и игнорирует мелкие колебания. С увеличением параметра gamma модель строит извилистую, сильно колеблющуюся кривую, пытаясь подстроиться под каждую точку данных.

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

• fit(const matrix &X) — генерирует матрицу случайных частот и вектор фазовых смещений,
• transform(const matrix &X) — преобразует входную матрицу X размерности N×D в матрицу новых признаков размера N×M,
• FitTransform(const matrix &X) — объединяет метод генерации весов и смещений и трансформации обучающей выборки X.

Если перемножить полученную матрицу сгенерированных признаков X_rff на ее транспонированную (X_rff @ X_rff^T), то результатом вычисления скалярных произведений между всеми парами векторов будет приближенная матрица RBF-ядра (матрица сходства).

Превращаем линейную ARD-регрессию в нелинейную модель

Проверим на практике связку CRBFSampler и ARDRegression. Для  этого создадим учебный пример - скрипт ARD_RBF.mq5

Для эксперимента сгенерируем выборку N=300 точек из сложной нелинейной функции: 

В качестве входного признака использовались равномерно распределённые значения в диапазоне [-4, 4]. Исходный одномерный признак X_raw был преобразован в 400-мерное пространство с помощью CRBFSampler(gamma=2.5, n_components=400). Значение gamma=2.5 выбрано в качестве примера — в реальной работе его рекомендуется подбирать через кросс-валидацию.

После преобразования полученная матрица передается в метод обучения модели ARD-регрессии ard.fit(). Из 400 случайных гармоник ARD отобрал всего 12 действительно значимых признаков, отбросив остальные 388 как бесполезные.

Помимо точечного прогноза модель рассчитывает дисперсию, что позволяет построить 95%-й доверительный интервал (±2σ).

На рис.3 показана аппроксимация нелинейной функции на обучающей выборке [-4,4], а также экстраполяция (прогноз) за пределы выборки [- 4.5, -4] & [4, 4.5]. Красной линией отмечено предсказание ARD-регрессии, голубой — истинная функция, серым пунктиром — границы 95%-го доверительного интервала.

рис. 3 Аппроксимация нелинейной функции с помощью ARD-регрессии и RBF ядра

Время обучения модели на MQL5 составило около 2.3 секунды, тогда как аналогичный код на Python выполняется за 2.8 секунды. Преимущество в скорости достигается за счёт использования функции LeastSquaresSolutionDC. В отличие от метода PInv() эта функция использует алгоритм Divide and Conquer (OpenBLAS/LAPACK). При этом решение системы A⋅X=I через LeastSquaresSolutionDC математически эквивалентно вычислению псевдообратной матрицы Мура — Пенроуза (A+) через SVD-разложение, что гарантирует точность и устойчивость вычислений на плохо обусловленных матрицах.

После обучения CRBFSampler сохраняет внутри себя сгенерированные случайные веса W и смещения b. Поэтому для корректной работы необходимо:

• один раз вызвать FitTransform() (или Fit()) на обучающей выборке,

• для всех новых данных использовать только метод Transform().

Повторный вызов FitTransform() приведёт к генерации новых случайных признаков, и веса обученной ARD-регрессии станут некорректными.

Таким образом связка RBFSampler + ARDRegression превращает обычную линейную  регрессию в нелинейный аппроксиматор функций любой сложности. 

Чтобы запустить скрипт scikit_learn.py, вам понадобится интерпретатор Python и стандартный стек статистических пакетов. О том, как их установить, рассказывается в статье Python + MetaTrader 5: быстрый исследовательский контур для данных, признаков и прототипов.

После установки пакетов поместите файл scikit_learn.py в папку (например, Scripts/ARD).

Запускайте скрипт через MetaEditor:

• в MetaEditor нажмите Файл → Открыть и выберите файл scikit_learn.py из папки ARD,

• после открытия файла нажмите клавишу F7 (или кнопку "Компилировать" на панели инструментов),

• результаты вычислений должны появиться во вкладке Журнал.

Заключение

В этой статье мы подробно разобрали математический фундамент и практическую реализацию ARD-регрессии (Automatic Relevance Determination) на MQL5, а также расширили её возможности до нелинейного моделирования с помощью случайных признаков Фурье (CRBFSampler).

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

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

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

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

#	Имя	Тип	Описание
1	ARDRegression.mqh	Включаемый файл	Содержит класс модели ARD регрессии и класс генерации случайных признаков Фурье (CRBFSampler)
2	scikit_learn.py	Скрипт	Пример обучения эталонной модели ARDregression из библиотеки scikit-learn
3	TestARD.mq5	Скрипт	Пример обучения MQL5 модели ARDregression
4	X_ARD.csv	CSV	Признаки
5	y_ARD.csv	CSV	Целевые значения
6	ARD_RBF.mq5	Скрипт	Пример обучения нелинейной модели ARD регрессии с помощью случайного преобразования признаков Фурье
7	MQL5.zip	Архив	Все вышеперечисленные файлы можно найти в архиве