Processos gaussianos em machine learning (Parte 2): Implementação e teste do modelo de classificação em MQL5
58
0
Introdução
No artigo anterior, apresentamos os fundamentos teóricos do modelo bayesiano de machine learning, os processos gaussianos (PG), e também iniciamos o desenvolvimento da biblioteca de PG em MQL5, descrevendo duas classes-chave: GaussianProcess e GPOptimizationObjective.
Hoje concluiremos o desenvolvimento da biblioteca, analisando em detalhes a implementação das interfaces-chave: IKernel, ILikelihood e IInference. Depois disso, testaremos a biblioteca com dados sintéticos e escreveremos indicadores para classificação e regressão, que mostrem sua operação em modo online, com o modelo sendo retreinado a cada nova barra.
Interface IKernel
A interface IKernel, que você encontrará no arquivo Kernels.mqh, é a base para a criação de kernels de covariância em nossa biblioteca. Ela torna o sistema flexível e fácil de estender, pois permite adicionar novos tipos de kernels ou suas combinações sem alterar a estrutura principal do código.
interface IKernel { // Вычисляет ковариационную матрицу между двумя наборами данных virtual matrix Compute(const matrix &X1, const matrix &X2) = 0; // Вычисляет производную ковариационной матрицы по заданному гиперпараметру // param_index: индекс гиперпараметра, по которому берется производная (начиная с 0) virtual matrix ComputeDerivative(int param_index) = 0; // Возвращает текущие значения всех гиперпараметров ядра virtual vector GetHyperparameters() const = 0; // Устанавливает новые значения гиперпараметров ядра virtual void SetHyperparameters(const vector ¶ms) = 0; // Возвращает количество гиперпараметров ядра virtual int GetNumHyperparameters() const = 0; // Возвращает строковое имя ядра virtual string GetName() const = 0; };
A interface define os dois métodos mais importantes que sustentam o funcionamento de qualquer kernel:
- Compute(const matrix &X1, const matrix &X2): calcula a matriz de covariância K entre dois conjuntos de dados.
- ComputeDerivative(int param_index): calcula a derivada da matriz de covariância em relação a um dos hiperparâmetros do kernel. O param_index indica o índice do hiperparâmetro em relação ao qual a derivada é obtida.
Classe RBFKernel (Kernel de base radial)
RBFKernel, também conhecido como kernel gaussiano ou exponencial quadrático, é um dos kernels de covariância mais usados. Ele tem dois hiperparâmetros: a variância do sinal σ_f (amplitude) e a escala de comprimento l (length), que determina a suavidade da função.
//+------------------------------------------------------------------+ //| Класс RBF ядра | //+------------------------------------------------------------------+ class RBFKernel : public IKernel { private: double length; double sigma_f; int n; // Количество строк в X1; int m; // Количество строк в X2; matrix m_K; matrix m_D_sq; // матрица квадратов расстояний ||x_i - x_j||^2 public: //Конструктор RBFKernel(double ls, double sf) : length(ls), sigma_f(sf) {} //Деструктор ~RBFKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); matrix XX(n, m); if (!XX.GeMM(X1, X2.Transpose(), 1.0, 0.0)) { Print("Ошибка: Не удалось вычислить XX = X * X^T"); return matrix::Zeros(n, m); } matrix diag1(n, 1); matrix X1_sq = X1 @ X1.Transpose(); diag1.Col(X1_sq.Diag(), 0); matrix diag2(m, 1); matrix X2_sq = X2 @ X2.Transpose(); diag2.Col(X2_sq.Diag(), 0); m_D_sq = matrix::Ones(n, 1) @ diag2.Transpose() + diag1 @ matrix::Ones(1, m) - 2 * XX; m_K = sigma_f * sigma_f * MathExp((-1 * m_D_sq) / (2 * length * length)); return m_K; } matrix ComputeDerivative(int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { dK = (2.0 / sigma_f) * m_K; break; } case 1: { dK = m_K * (m_D_sq / (length * length * length)); break; } } return dK; } vector GetHyperparameters() const override { vector params(2); params[0] = sigma_f; params[1] = length; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 2) { sigma_f = params[0]; length = params[1]; } } int GetNumHyperparameters() const override { return 2; } string GetName() const override { return "RBFKernel"; } };
Felizmente, para o kernel RBF, as derivadas são fáceis de obter:
Derivada em relação a sigma_f:
- Derivada em relação a l:
Para o parâmetro length, multiplicamos a matriz K elemento a elemento por (D_sq / length^3), sendo D_sq a matriz dos quadrados das distâncias euclidianas.
Classe LinearKernel (Kernel linear)
LinearKernel é um kernel de covariância simples que pressupõe uma dependência linear entre os pontos de dados. Ele é frequentemente aplicado quando se espera que a função seja aproximável por um modelo linear, ou como componente em kernels compostos mais complexos.
O kernel linear possui um hiperparâmetro: a variância do sinal σ_l.
//+------------------------------------------------------------------+ //| Класс для линейного ядра | //+------------------------------------------------------------------+ class LinearKernel : public IKernel { private: double sigma_l; int n; int m; matrix m_K; public: //Конструктор LinearKernel(double sl): sigma_l(sl){} //Деструктор ~LinearKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); m_K.GeMM(X1, X2.Transpose(), sigma_l * sigma_l, 0.0); return m_K; } matrix ComputeDerivative(int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { // Производная по sigma_l dK = (2.0 / sigma_l) * m_K; break; } } return dK; } vector GetHyperparameters() const override { vector params(1); params[0] = sigma_l; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 1) { sigma_l = params[0]; } } int GetNumHyperparameters() const override { return 1; } string GetName() const override { return "LinearKernel"; } };
Método ComputeDerivative(int param_index):
Classe PeriodicKernel (Kernel periódico)
PeriodicKernel é usado para modelar dados com padrões repetitivos ou sazonalidade, permitindo captar dependências cíclicas.
//+------------------------------------------------------------------+ //| Класс для периодического ядра | //+------------------------------------------------------------------+ class PeriodicKernel : public IKernel { private: double sigma_f; // Дисперсия сигнала (амплитуда колебаний) double length; // Длина масштаба (гладкость колебаний) double period; // Период колебаний int n; // Количество строк в X1 int m; // Количество строк в X2 matrix m_K; // Ковариационная матрица matrix m_D; // Матрица D (сумма производных по d [2 * sin^2(M_PI*distance/period] ) matrix m_distance_dim[]; // Массив матриц разностей по каждому измерению |x_i_k - x_j_k| int m_d_cols; // Количество измерений(признаков) (X.Cols()) public: //Конструктор PeriodicKernel(double ls, double sf, double p) : sigma_f(sf), length(ls), period(p){} //Деструктор ~PeriodicKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); int d = (int)X1.Cols(); m_d_cols = d; m_D = matrix::Zeros(n, m); ArrayResize(m_distance_dim, d); for (int k = 0; k < d; k++) { vector x1_k = X1.Col(k); vector x2_k = X2.Col(k); // Создаем матрицу, где каждая колонка - это вектор x1_k // (n x 1) * (1 x m) = (n x m) matrix x1_k_copy = x1_k.Outer(vector::Ones(m)) ; // Создаем матрицу, где каждая строка - это транспонированный вектор x2_k // (n x 1) * (1 x m) = (n x m) matrix x2_k_copy = vector::Ones(n).Outer(x2_k); // Вычисляем матрицу всех попарных абсолютных разностей между элементами векторов x1_k и x2_k matrix distance = MathAbs(x1_k_copy - x2_k_copy); m_distance_dim[k] = distance; // Кешируем distance для каждого измерения(признака) matrix sin_term = MathSin(M_PI * distance / period); sin_term = 2.0 * sin_term * sin_term; // 2 * sin^2(u) m_D += sin_term; // Суммируем поэлементно } // Вычисляем матрицу K m_K = sigma_f * sigma_f * MathExp(-1 * m_D / (length * length)); return m_K; } matrix ComputeDerivative( int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { // Производная по sigma_f // dK/d(sigma_f) = (2 / sigma_f) * K dK = (2.0 / sigma_f) * m_K; break; } case 1: { // Производная по length // dK/d(length) = K * (2 * D / length^3) dK = m_K * ( 2.0 * m_D / (length * length * length) ); break; } case 2: { // Производная по period // dK/d(period) = K * (-1/l^2) * dD/d(period) // dD/d(period) = Sum{k=1:d} (2 * pi * (x_ik - x_jk) / period^2) * sin(2*pi*(x_ik - x_jk)/period) matrix dD_dp = matrix::Zeros(n, m); for (int k = 0; k < m_d_cols; k++) { // Цикл по каждому измерению (колонке) данных matrix distance_k = m_distance_dim[k]; // матрица абсолютных разностей для k-го измерения matrix u = M_PI * distance_k / period; // Вычисление аргумента u_k = pi * |x_ik - x_jk| / period matrix sin_2u = MathSin(2.0 * u); // Вычисление sin(2 * u_k) // Используем уже вычисленное u для расчета // -pi * |x_ik - x_jk| / period^2 matrix second_term = (-1*u) / period; // Накопление dD/d(period) для текущего измерения dD_dp += 2.0 * sin_2u * second_term; } // объединение всех частей производной // dK = K * (-1/length^2) * dD_dp dK = m_K * (-1.0 / (length * length)) * dD_dp; break; } } return dK; } vector GetHyperparameters() const override { vector params(3); params[0] = sigma_f; params[1] = length; params[2] = period; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 3) { sigma_f = params[0]; length = params[1]; period = params[2]; } } int GetNumHyperparameters() const override { return 3; } string GetName() const override { return "PeriodicKernel"; } };
Derivada do kernel periódico em relação ao parâmetro period:
O cálculo da derivada em relação ao parâmetro period é o mais complexo, pois period está dentro de uma função trigonométrica, que, por sua vez, faz parte da exponencial.
A fórmula da derivada se baseia na regra da cadeia, em que K depende de D, e D depende de period por meio dos termos senoidais. A implementação de ComputeDerivative calcula iterativamente o termo dD_dp para cada dimensão e, em seguida, obtém a derivada resultante:
em que dD_dp é calculado como a soma das derivadas de 2 * sin^2(M_PI*distance/period) em todas as dimensões (atributos).
Kernels compostos
Os processos gaussianos permitem combinar kernels de covariância simples para construir modelos mais complexos, capazes de descrever diferentes aspectos dos dados (por exemplo, tendência, periodicidade, ruído). Essa funcionalidade usa dois kernels compostos:
- SumKernel: combina vários kernels somando suas matrizes de covariância.
- ProductKernel: combina kernels multiplicando suas matrizes de covariância elemento a elemento.
Ambas as classes gerenciam os hiperparâmetros de seus kernels filhos, reunindo-os em um único vetor para otimização e, em seguida, distribuindo-os de volta para os kernels filhos.
Classe SumKernel
//+------------------------------------------------------------------+ //| Класс для суммы ядер | //+------------------------------------------------------------------+ class SumKernel : public IKernel { private: IKernel* kernels[]; public: // Конструктор SumKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } //Деструктор ~SumKernel() { for (int i = 0; i < ArraySize(kernels); i++) { if (kernels[i] != NULL) { delete kernels[i]; } } ArrayFree(kernels); } matrix Compute(const matrix &X1, const matrix &X2) override { int n = (int)X1.Rows(); int m = (int)X2.Rows(); matrix sum = matrix::Zeros(n,m); for (int i = 0; i < ArraySize(kernels); i++) { sum += kernels[i].Compute(X1, X2); // Суммируем матрицы от каждого дочернего ядра } return sum; } matrix ComputeDerivative(int param_index) override { int current_param_offset = 0; // Начинаем со смещения 0 для первого ядра for (int i = 0; i < ArraySize(kernels); i++) { // Получаем количество гиперпараметров текущего дочернего ядра int num_params_current_kernel = kernels[i].GetNumHyperparameters(); // Проверяем, принадлежит ли глобальный param_index текущему дочернему ядру // param_index должен быть: // - больше или равен текущему смещению // - строго меньше, чем текущее смещение + количество параметров текущего ядра if (param_index >= current_param_offset && param_index < current_param_offset + num_params_current_kernel) { // Если param_index принадлежит этому ядру, вычисляем его локальный индекс int local_param_index = param_index - current_param_offset; // Вызываем ComputeDerivative для найденного дочернего ядра и // передаем ему локальный индекс. // Поскольку производная суммы ядер по параметру одного ядра равна производной // этого конкретного ядра по этому параметру, мы можем сразу вернуть результат return kernels[i].ComputeDerivative(local_param_index); } // Если param_index не принадлежит текущему ядру, // увеличиваем смещение для проверки следующего ядра current_param_offset += num_params_current_kernel; } Print("Error: SumKernel::ComputeDerivative - Parameter index ", param_index, " out of bounds"); return matrix::Zeros(1, 1); } //+-------------------------------------------------------------------+ //| Функция возвращает значения всех гиперпараметров составного ядра | //+-------------------------------------------------------------------+ vector GetHyperparameters() const override { vector all_params(GetNumHyperparameters()); int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Получаем значения параметров текущего дочернего ядра vector kernel_params = kernels[i].GetHyperparameters(); for (int j = 0; j < (int)kernel_params.Size(); j++) { // Копируем параметры в общий вектор all_params[current_idx + j] = kernel_params[j]; } // Обновляем смещение для следующего ядра current_idx += (int)kernel_params.Size(); } return all_params; } //+----------------------------------------------------------------------+ //|Функция принимает вектор значений всех гиперпараметров, разбирает | //|его на соответствующие части и присваивает их каждому дочернему ядру | //+----------------------------------------------------------------------+ void SetHyperparameters(const vector ¶ms) override { int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Получаем количество гиперпараметров текущего дочернего ядра int num_params = kernels[i].GetNumHyperparameters(); vector sub_params(num_params); // Копируем соответствующие параметры из общего вектора for (int j = 0; j < num_params; j++) { sub_params[j] = params[current_idx + j]; } // Вызываем метод SetHyperparameters() на текущем дочернем ядре, // передавая ему только его собственные параметры, собранные в векторе sub_params kernels[i].SetHyperparameters(sub_params); current_idx += num_params; } } //+---------------------------------------------------------------------------------+ //| Возвращает общее количество гиперпараметров для всего составного ядра SumKernel | | //+---------------------------------------------------------------------------------+ int GetNumHyperparameters() const override { int total_params = 0; for (int i = 0; i < ArraySize(kernels); i++) { total_params += kernels[i].GetNumHyperparameters(); } return total_params; } string GetName() const override { string name = "Sum("; for (int i = 0; i < ArraySize(kernels); i++) { name += kernels[i].GetName(); if (i < ArraySize(kernels) - 1) name += ","; } name += ")"; return name; } //+------------------------------------------------------------------+ //| Предоставляет внешний доступ к списку дочерних ядер | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { // Копируем указатели из внутреннего массива kernels в переданный // внешний массив output_kernels. Это обеспечивает доступ к дочерним // ядрам в методе GaussianProcess::Fit(), где требуется пройтись // в цикле по всем элементарным ядрам для установки границ оптимизации // их гиперпараметров. ArrayCopy(output_kernels, kernels); } };
- Método Compute: calcula a matriz de covariância final K como uma soma simples das matrizes de covariância retornadas por cada kernel filho. Aqui, aplica-se o princípio do polimorfismo: embora o array kernels armazene ponteiros do tipo base IKernel*, a chamada kernels[i].Compute(X1, X2), na prática, chama a implementação específica de Compute para cada kernel filho concreto. Isso permite que SumKernel trabalhe com qualquer kernel herdado de IKernel.
- Método ComputeDerivative: a derivada da função de covariância da soma de kernels em relação a um hiperparâmetro específico corresponde à derivada da função de covariância somente do kernel filho ao qual esse hiperparâmetro pertence. Isso significa que ComputeDerivative encontra o kernel filho correspondente por param_index e retorna somente sua derivada. As derivadas em relação aos parâmetros dos demais kernels filhos são consideradas iguais a zero.
Classe ProductKernel
//+------------------------------------------------------------------+ //| Класс для произведения ядер | //+------------------------------------------------------------------+ class ProductKernel : public IKernel { private: IKernel* kernels[]; matrix m_X1; matrix m_X2; int n; int m; public: // Конструктор: копирует указатели на дочерние ядра ProductKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } // Деструктор ~ProductKernel() { for (int i = 0; i < ArraySize(kernels); i++) { if (kernels[i] != NULL){ delete kernels[i]; } } ArrayFree(kernels); } matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); m_X1 = X1; m_X2 = X2; matrix product = matrix:: Ones(n,m); for (int i = 0; i < ArraySize(kernels); i++) { // Полиморфизм: вызываем метод Compute() каждого дочернего ядра, // независимо от его конкретного типа, и поэлементно умножаем результат product *= kernels[i].Compute(X1, X2); } return product; } matrix ComputeDerivative(int param_index) override { matrix dK_prod(n, m); // Итоговая матрица производной int current_param_offset = 0; int target_kernel_idx = -1; // Индекс ядра, которому принадлежит гиперпараметр int local_param_index = -1; // Локальный индекс гиперпараметра в этом ядре // Шаг 1: Находим ядро и локальный индекс гиперпараметра for (int i = 0; i < ArraySize(kernels); i++) { int num_params_current_kernel = kernels[i].GetNumHyperparameters(); if (param_index >= current_param_offset && param_index < current_param_offset + num_params_current_kernel) { target_kernel_idx = i; local_param_index = param_index - current_param_offset; break; } current_param_offset += num_params_current_kernel; } if (target_kernel_idx == -1) { Print("Error: ProductKernel::ComputeDerivative - Parameter index ", param_index, " out of bounds"); return matrix::Zeros(1, 1); } // Шаг 2: Вычисляем dK_k / d_theta_j для целевого ядра matrix dK_target_kernel = kernels[target_kernel_idx].ComputeDerivative(local_param_index); // Шаг 3: Вычисляем произведение K_m для всех остальных ядер (m != k) matrix other_kernels_product = matrix::Ones(n, m); for (int i = 0; i < ArraySize(kernels); i++) { if (i != target_kernel_idx) { // Пропускаем целевое ядро other_kernels_product = other_kernels_product * kernels[i].Compute(m_X1, m_X2); } } // Шаг 4: Поэлементно перемножаем результаты dK_prod = dK_target_kernel * other_kernels_product; return dK_prod; } // Собирает значения гиперпараметров всех дочерних ядер в один вектор vector GetHyperparameters() const override { vector all_params(GetNumHyperparameters()); int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { vector kernel_params = kernels[i].GetHyperparameters(); for (int j = 0; j < (int)kernel_params.Size(); j++) { all_params[current_idx + j] = kernel_params[j]; } current_idx += (int)kernel_params.Size(); } return all_params; } void SetHyperparameters(const vector ¶ms) override { int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { int num_params_kernel = kernels[i].GetNumHyperparameters(); vector sub_params(num_params_kernel); for (int j = 0; j < num_params_kernel; j++) { sub_params[j] = params[current_idx + j]; } kernels[i].SetHyperparameters(sub_params); current_idx += num_params_kernel; } } int GetNumHyperparameters() const override { int total_params = 0; for (int i = 0; i < ArraySize(kernels); i++) { total_params += kernels[i].GetNumHyperparameters(); } return total_params; } string GetName() const override { string name = "Prod("; for (int i = 0; i < ArraySize(kernels); i++) { name += kernels[i].GetName(); if (i < ArraySize(kernels) - 1) name += "*"; } name += ")"; return name; } //+------------------------------------------------------------------+ //| Предоставляет внешний доступ к списку дочерних ядер | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { ArrayCopy(output_kernels, kernels); } };
- Método Compute: calcula a matriz de covariância final K como o produto elemento a elemento das matrizes de covariância retornadas por cada kernel filho. Assim como em SumKernel, aqui o polimorfismo é usado para chamar o método Compute do kernel filho correspondente.
- Método ComputeDerivative: para calcular a derivada de K em relação ao hiperparâmetro pertencente a um kernel filho específico, K_p, ProductKernel aplica a regra do produto (regra de Leibniz). A derivada de K em relação a esse parâmetro é igual à derivada da matriz K_p em relação a esse parâmetro, multiplicada pelo produto elemento a elemento das matrizes de covariância de todos os demais kernels filhos, sem derivar:
Interface ILikelihood
A interface ILikelihood define como qualquer função de verossimilhança deve operar em nossa biblioteca. Ela nos permite trabalhar com diferentes tipos de dados, tanto com valores contínuos (como na regressão) quanto com categorias discretas (como na classificação).
//+------------------------------------------------------------------+ //|Интерфейс для функций правдоподобия | //+------------------------------------------------------------------+ interface ILikelihood { // Вычисляет логарифм правдоподобия p(y|f) для скрытых значений f и наблюдаемых y virtual double LogLikelihood(const vector &f, const vector &y) = 0; //----------------------------- Производные логарифма правдоподобия по f ----------------- // Вычисляет вектор первой производной логарифма правдоподобия по f (dlp/df) virtual vector LogLikelihoodGradient(const vector &f, const vector &y) = 0; // Вычисляет матрицу второй производной (Гессиан) логарифма правдоподобия по f (d^2lp/df df^T) virtual matrix LogLikelihoodHessian(const vector &f, const vector &y) = 0; // Вычисляет вектор третьей производной логарифма правдоподобия по f (d^3lp/df^3) virtual vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) = 0; //----------------------------- Производные логарифма правдоподобия по параметрам ----------------- // Первая производная логарифма правдоподобия по j-му гиперпараметру правдоподобия. virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) = 0; // Производная Гессиана логарифма правдоподобия по j-му гиперпараметру правдоподобия. virtual matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) = 0; // Возвращает имя функции правдоподобия virtual string GetName() const = 0; // Возвращает текущие значения гиперпараметров правдоподобия virtual vector GetHyperparameters() const = 0; // Устанавливает новые значения гиперпараметров правдоподобия virtual void SetHyperparameters(const vector ¶ms) = 0; // Возвращает количество гиперпараметров правдоподобия virtual int GetNumHyperparameters() const = 0; };
Classe GaussianLikelihood
GaussianLikelihood implementa a função de verossimilhança para tarefas de regressão, assumindo ruído normal nos dados observados.
//+------------------------------------------------------------------+ //| Класс Гауссовского правдоподобия (для задач регрессии) | //+------------------------------------------------------------------+ class GaussianLikelihood : public ILikelihood { private: double m_noise_sigma; // Параметр шума (стандартное отклонение) public: // Конструктор GaussianLikelihood(double initial_noise_sigma) : m_noise_sigma(initial_noise_sigma) {} // Возвращает текущее значения гиперпараметра vector GetHyperparameters() const override { vector params(1); params[0] = m_noise_sigma; return params; } // Устанавливает новое значения гиперпараметра void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 1) { m_noise_sigma = params[0]; } } // Возвращает количество гиперпараметров int GetNumHyperparameters() const override { return 1; } // Вычисляет логарифм правдоподобия log p(y|f) для гауссовского распределения double LogLikelihood(const vector &f, const vector &y) override { int n = (int)y.Size(); double noise_variance = m_noise_sigma * m_noise_sigma; vector residual = y - f; // Формула логарифма плотности многомерного нормального распределения return -0.5 / noise_variance * (residual @ residual) - 0.5 * n * MathLog(2 * M_PI * noise_variance); } //------------------ Производные по f ----------------------- // Вычисляет вектор первой производной логарифма правдоподобия по f (dlp/df) vector LogLikelihoodGradient(const vector &f, const vector &y) override { // d(log p(y|f))/df = (y - f) / sigma^2 return (y - f) / (m_noise_sigma * m_noise_sigma); } // Вычисляет матрицу второй производной (Гессиан) логарифма правдоподобия по f (d^2lp/dfdf^T) matrix LogLikelihoodHessian(const vector &f, const vector &y) override { // d^2(log p(y|f))/dfdf^T = -1/sigma^2 * I int n = (int)y.Size(); return matrix::Identity(n, n) * (-1.0 / (m_noise_sigma * m_noise_sigma)); } // Вычисляет вектор третьей производной логарифма правдоподобия по f vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) override { // Для Гауссовского правдоподобия d^3(log p(y|f))/df^3 = 0 return vector::Zeros((int)y.Size()); } // ------------------------------ Производные по параметру ----------------------------------- // // Вычисляет первую производную логарифма правдоподобия по j-му гиперпараметру правдоподобия virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { // Гауссовское правдоподобие имеет только 1 гиперпараметр: m_noise_sigma (индекс 0) if (param_index == 0) { // Формула для d(log p(y|f))/d(sigma_n): // = (y-f)^T(y-f) / sigma_n^3 - N / sigma_n vector y_minus_f = y - f; double sn3 = m_noise_sigma * m_noise_sigma * m_noise_sigma; double term1 = (y_minus_f @ y_minus_f) / sn3; double term2 = y.Size() / m_noise_sigma; // N / sigma_n return term1 - term2; } else { return 0.0; } } // Вычисляет производную Гессиана логарифма правдоподобия по j-му гиперпараметру правдоподобия matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) override { int n = (int)y.Size(); if (param_index == 0) { // Производная Hessian по m_noise_sigma: d/d(sigma_n) [-1/sigma_n^2 * I] = (2/sigma_n^3) * I double derivative_coeff = 2.0 / (m_noise_sigma * m_noise_sigma * m_noise_sigma); return matrix::Identity(n, n) * derivative_coeff; } else { return matrix::Zeros(n, n); } } string GetName() const override { return "GaussianLikelihood"; } };
Métodos principais:
- Construtor GaussianLikelihood: inicializa o único hiperparâmetro, o desvio padrão do ruído.
- LogLikelihood(): calcula o logaritmo da verossimilhança p(y∣f) para a densidade de probabilidade da distribuição normal multivariada:
- LogLikelihoodGradient(): calcula a primeira derivada do logaritmo da verossimilhança em relação aos valores latentes f. É um vetor em que cada elemento é igual a:
- LogLikelihoodHessian(): calcula a segunda derivada (Hessiano) do logaritmo da verossimilhança em relação a f. Para a verossimilhança gaussiana, a matriz Hessiana é diagonal:
- LogLikelihoodThirdDerivative(): calcula a terceira derivada. Para a verossimilhança gaussiana, todas as derivadas de ordem superior à segunda são iguais a zero.
- LogLikelihoodGradientParam(): calcula a primeira derivada do logaritmo da verossimilhança em relação ao hiperparâmetro m_noise_sigma:
- LogLikelihoodHessianDerivative(): calcula a derivada da Hessiana do logaritmo da verossimilhança em relação ao hiperparâmetro m_noise_sigma:
Classe LogitLikelihood
LogitLikelihood é usada para classificação binária, em que os valores-alvo observados y assumem os valores {−1,+1}. Essa função relaciona a função latente f à probabilidade de pertencimento a uma classe via função sigmoide.
A implementação de LogitLikelihood inclui as funções auxiliares sigmoid() e Softplus(), com tratamento de valores de entrada muito grandes ou muito pequenos.
//+------------------------------------------------------------------+ //| Класс для Логит-правдоподобия (y {-1, +1}) | //+------------------------------------------------------------------+ class LogitLikelihood : public ILikelihood { public: // logit_sigmoid(x) = 1 / (1 + exp(-x)) double sigmoid(double x) const { if (x > 100.0) return 1.0; // Избегаем появления NaN if (x < -100.0) return 0.0; return 1.0 / (1.0 + MathExp(-x)); } // Функция Softplus log(1 + exp(x)) double Softplus(double x) const { if (x > 100.0) return x; // Для очень больших x, log(1 + exp(x)) ~ x if (x < -100.0) return MathExp(x); // Для очень маленьких x, log(1 + exp(x)) ~ exp(x) return MathLog(1.0 + MathExp(x)); } public: // Конструктор (логит-правдоподобие не имеет собственных гиперпараметров) LogitLikelihood() {} vector GetHyperparameters() const override { return vector::Zeros(0); } void SetHyperparameters(const vector ¶ms) override {} // Возвращает количество гиперпараметров правдоподобия int GetNumHyperparameters() const override { return 0; } // --- Вычисляет логарифм правдоподобия: log p(y|f) --- // Формула: sum_i (-log(1 + exp(-y_i * f_i))) double LogLikelihood(const vector &f, const vector &y) override { int n = (int)y.Size(); double total_log_likelihood = 0.0; for (int i = 0; i < n; i++) { total_log_likelihood += -Softplus(-y[i] * f[i]); } return total_log_likelihood; } // Вычисляет градиент логарифма правдоподобия по f // Формула: d(log p(y_i|f_i))/df_i = y_i * sigma(-y_i * f_i) vector LogLikelihoodGradient(const vector &f, const vector &y) override { int n = (int)y.Size(); vector grad(n); for (int i = 0; i < n; i++) { grad[i] = y[i] * sigmoid(-y[i] * f[i]); // grad[i] = y[i] * (1 - sigmoid(y[i] * f[i])); эквивалентно } return grad; } // Вычисляет Гессиан логарифма правдоподобия по f (диагональная матрица H) // Формула: d^2(log p(y_i|f_i))/df_i^2 = -sigma(y_i*f_i)(1 - sigma(y_i*f_i)) matrix LogLikelihoodHessian(const vector &f, const vector &y) override { int n = (int)y.Size(); matrix H = matrix::Identity(n, n); for (int i = 0; i < n; i++) { double sig = sigmoid(y[i] * f[i]); H[i][i] = -1*sig * (1.0 - sig); } return H; } // Вычисляет третью производную логарифма правдоподобия по f // Формула: d^3(log p(y_i|f_i))/df_i^3 = y_i * sigma(-y_i f_i) * (1 - sigma(-y_i f_i)) * (1 - 2*sigma(-y_i f_i)) vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) override { int n = (int)y.Size(); vector third_deriv(n); for (int i = 0; i < n; i++) { double sig_neg_yf = sigmoid(-y[i] * f[i]); third_deriv[i] = y[i] * sig_neg_yf * (1.0 - sig_neg_yf) * (1.0 - 2.0 * sig_neg_yf); // эквивалентно // double sig_yf = sigmoid(y[i] * f[i]); // third_deriv[i] = -y[i] * sig_yf * (1.0 - sig_yf) * (1.0 - 2.0 * sig_yf); } return third_deriv; } // LogitLikelihood не имеет гиперпараметров virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { return 0.0; } // LogitLikelihood не имеет гиперпараметров matrix LogLikelihoodHessianDerivative(const vector &f_latent, const vector &y, int param_index) override { int n = (int)y.Size(); return matrix::Zeros(n, n); } string GetName() const override { return "LogitLikelihood"; } };
Métodos principais:
- LogLikelihood(): calcula o logaritmo da verossimilhança logp(y∣f). Para classificação binária com y em {−1,+1}, isso corresponde à soma dos logaritmos da verossimilhança de entropia cruzada binária:
- LogLikelihoodGradient(): calcula a primeira derivada do logaritmo da verossimilhança em relação a f. Cada elemento do vetor gradiente é igual a:
- LogLikelihoodHessian(): calcula a segunda derivada (Hessiana) do logaritmo da verossimilhança em relação a f. Para LogitLikelihood, a matriz Hessiana é diagonal cujos elementos diagonais são iguais a:
- LogLikelihoodThirdDerivative(): calcula a terceira derivada do logaritmo da verossimilhança em relação a f. Cada elemento do vetor é igual a:
Funções e estruturas auxiliares
O arquivo StructUtils.mqh contém um conjunto de enumerações, estruturas de dados e funções que facilitam a operação com PG.
enum PredictMode { PROBIT = 0, // Probit Approximation NUM_INTEGR = 1, // Numerical integration MONTE_CARLO = 2 // Monte Carlo }; //--- Структура для результатов предсказания struct GPPredictionResult { vector mu_f_star; // Апостериорное среднее для скрытой функции f* на новых данных matrix Sigma_f_star; // Апостериорная ковариация для скрытой функции f* на новых данных // Поля только для регрессии (GaussianLikelihood) vector mu_y_star; // Апостериорное среднее наблюдений y* matrix Sigma_y_star; // Апостериорная ковариация наблюдений y* // Поля только для классификации (LogitLikelihood, ProbitLikelihood и т.д.) vector predicted_probabilities; // Вероятности p(y*=+1 | X*, D) vector predicted_labels; // Предсказанные метки y* (+1 или -1) }; // Структура для результата инференса struct GPInferenceResult { double nlml_value; // Отрицательный логарифм маргинального правдоподобия vector nlml_gradient; // Градиент NLML по гиперпараметрам ядра и правдоподобия // Результаты инференса на обучающих данных: vector mu_f_train; // Апостериорное среднее скрытой функции f на обучающих данных matrix Sigma_f_train; // Апостериорная ковариация скрытой функции f на обучающих данных (K^−1+W)^−1 matrix L_K_noisy; // Разложение Холецкого K_noisy = cholesky(K + s2n*I) matrix L_B; // Разложение Холецкого B = I + W^0.5 @ K @ W^0.5 matrix sW; // sW = W^0.5; vector sW_diag; // диагональ sW matrix sW_K; // sW @ K vector alpha; // (K_noisy)^-1 * y matrix H; // Гессиан логарифма правдоподобия d^2 log p(y|f)/df^2 (для LaplaceInference) bool success; // Флаг успешности функции инференса };
- A enumeração PredictMode define diferentes modos para executar previsões no modelo de classificação.
- A estrutura GPInferenceResult é usada para armazenar todos os resultados principais obtidos durante a inferência.
- A estrutura GPPredictionResult armazena todos os resultados obtidos após a previsão em novos dados de teste (X_test).
- DiagonalTrace
//+-------------------------------------------------------------------+ //| Вычисляет след произведения двух матриц, используя только | //| диагональные элементы итоговой матрицы | //+-------------------------------------------------------------------+ double DiagonalTrace(const matrix &A, const matrix &B) { ulong m = A.Rows(); ulong n = A.Cols(); ulong n_b = B.Rows(); ulong p = B.Cols(); // Находим минимальный размер для диагонали ulong k = MathMin(m, p); double tr = 0.0; for(ulong i = 0; i < k; i++) { // Скалярное произведение i-й строки A и i-го столбца B tr += A.Row(i)@B.Col(i); } return tr; }
Calcula o traço do produto de duas matrizes A e B (Tr(A @ B)) sem precisar formar toda a matriz produto. Para isso, soma os produtos escalares das linhas da matriz A com as colunas correspondentes da matriz B. Essa abordagem permite reduzir significativamente o custo computacional em comparação com a multiplicação completa de matrizes.
- cho_solve (wrapper da OpenBLAS para LinearEquationsSolutionTriangular)
//+------------------------------------------------------------------+ //| Решает систему A*X = B, используя разложение Холецкого A = LL^T | //| Параметры: | //| c: Нижняя треугольная матрица L из разложения Холецкого | //| b: правая часть системы (матрица) | //| Возвращает: матрица X - решение системы | //+------------------------------------------------------------------+ matrix cho_solve(const matrix &c, const matrix &b) { // Проверка, является ли матрица 'c' нижней треугольной if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular."); return matrix::Zeros(c.Rows(), b.Cols()); } // Шаг 1: Решаем L * Y = B // Y - промежуточная матрица, результат решения L^-1 * B matrix Y; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_N, b, Y)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L * Y = B failed. Error Code: %d", GetLastError()); return matrix::Zeros(c.Rows(), b.Cols()); } // Шаг 2: Решаем L^T * X = Y // X - решение системы L^T * X = Y, что эквивалентно (L^T)^-1 * Y matrix X; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_T, Y, X)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L^T * X = Y failed. Error Code: %d", GetLastError()); return matrix::Zeros(c.Rows(), b.Cols()); } return X; } //+------------------------------------------------------------------+ //| Решает систему Ax = b, используя разложение Холецкого A = LL^T | //| Параметры: | //| c: Нижняя треугольная матрица L из разложения Холецкого | //| b: правая часть системы (вектор) | //| Возвращает: вектор x - решение системы | //+------------------------------------------------------------------+ vector cho_solve(const matrix &c, const vector &b) { // Проверка, является ли матрица c нижней треугольной if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular"); return vector::Zeros(c.Rows()); } // Шаг 1: Решаем L * y = b // y - промежуточный вектор, результат решения L^-1 * b vector y; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_N, b, y)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L * y = b failed. Error Code: %d", GetLastError()); return vector::Zeros(c.Rows()); } // Шаг 2: Решаем L^T * x = y // x - конечное решение системы L^T * x = y, что эквивалентно (L^T)^-1 * y vector x; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_T, y, x)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L^T * x = y failed. Error Code: %d", GetLastError()); return vector::Zeros(c.Rows()); } return x; }
A função cho_solve serve para resolver com eficiência sistemas de equações lineares Ax=b (ou AX=B quando B é uma matriz), em que a matriz A é obtida a partir de sua decomposição de Cholesky A=LL^T. Aqui, L é a matriz triangular inferior.
O cálculo direto de A^−1 (A.Inv()) exige muitos recursos e pode ser numericamente instável. A decomposição de Cholesky, por outro lado, oferece uma abordagem significativamente mais eficiente e confiável.
Em cho_solve, são executadas duas operações principais:
- Substituição progressiva: resolve-se o sistema LY=B (ou Ly=b), em que L é a matriz de entrada c e Y (ou y) é o resultado intermediário. Como L é uma matriz triangular, essa resolução é relativamente rápida.
- Substituição regressiva: em seguida, resolve-se o sistema L^TX=Y (ou L^Tx=y), em que L^T é a matriz transposta de L. Isso também é eficiente, pois L^T é uma matriz triangular superior.
Essas operações são implementadas com as funções LinearEquationsSolutionTriangular da biblioteca de álgebra linear de alto desempenho OpenBLAS. Também proporcionam uma aceleração significativa dos cálculos e se tornam indispensáveis para modelos de machine learning com alto consumo de recursos, como os processos gaussianos.
Interface IInference
interface IInference { // Этот метод выполняет инференс (вывод апостериорного распределения f) // и возвращает компоненты, необходимые для NLML и предсказания virtual void Infer(const matrix &X, const vector &y, IKernel *kernel, ILikelihood *likelihood,GPInferenceResult &result) = 0; virtual string GetName() const = 0; };
A interface permite alternar entre diferentes métodos de inferência, como a inferência exata para a verossimilhança gaussiana e métodos aproximados para casos não gaussianos, sem alterar a lógica central de GaussianProcess.
Infer é o método central desta interface. Ele executa a inferência da distribuição a posteriori da função latente f e retorna os componentes necessários para cálculos posteriores (NLML e previsões). Recebe os seguintes argumentos:
- X, matriz de atributos de treinamento,
- y, vetor de rótulos de treinamento,
- kernel, ponteiro para um objeto IKernel,
- likelihood, ponteiro para um objeto ILikelihood,
- result, referência à estrutura GPInferenceResult, usada para armazenar todos os resultados da inferência, incluindo:
- o logaritmo negativo da verossimilhança marginal (NLML), necessário para a otimização dos hiperparâmetros,
- o vetor de gradientes de NLML em relação a todos os hiperparâmetros do kernel e da verossimilhança,
- matrizes e vetores auxiliares: (L_K_noisy, L_B, sW, mu_f_train, Sigma_f_train, alpha), usados como suporte para calcular os gradientes de NLML e para a posterior previsão para novos dados.
Classe ExactInference
//+------------------------------------------------------------------+ //| ExactInference Только для гауссовского правдоподобия | //+------------------------------------------------------------------+ class ExactInference : public IInference { public: ExactInference() {} void Infer(const matrix &X, const vector &y, IKernel *kernel, ILikelihood *likelihood,GPInferenceResult &result) override { if(likelihood.GetName() != "GaussianLikelihood") { Print("ExactInference supports only GaussianLikelihood"); result.nlml_value = DBL_MAX; return; } result.success = false; int n = (int)y.Size(); // Вычисление K с текущими параметрами ядра matrix K = kernel.Compute(X, X); // Получаем дисперсию шума vector likelihood_params = likelihood.GetHyperparameters(); double sigma = likelihood_params[0]; double variance = sigma * sigma; double jitter = 1e-6; matrix K_noisy = K + matrix::Identity(n, n) * (variance + jitter); if(!K_noisy.Cholesky(result.L_K_noisy)) { PrintFormat("Error: Cholesky decomposition failed. Error Code: %d",GetLastError()); result.nlml_value = DBL_MAX; return; } //------------- Algorithm 2.1 GPML --------------------------- //Вычисляем alpha (K_noisy^-1 * y) - O(N^2) result.alpha = cho_solve(result.L_K_noisy, y); //+------------------------------------------------------------------+ //| NLML = 0.5y^T(K+σ^2*I)^-1*y + 0.5*log|K+σ^2*I| +n/2*Log(2π) | //+------------------------------------------------------------------+ //+------------------------------------------------------------------+ //| NLML = 0.5y^T*alpha + 0.5*log|K+σ^2*I| +n/2*Log(2π) | //+------------------------------------------------------------------+ // --- Вычисляем первое слагаемое в формуле NLML : double data_term = 0.5 * (y @ result.alpha); // --- Вычисляем второе слагаемое: 1/2 * log|K + sigma^2*I| = sum(log(L_ii)) double log_det = MathLog(result.L_K_noisy.Diag()).Sum(); // --- Вычисляем третье слагаемое: n/2 * log(2π) double const_term = 0.5 * n * MathLog(2 * M_PI); //--- NLML result.nlml_value = data_term + log_det + const_term; // -------------------- Вычисления градиентов NLML -------------------------------------- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // Вычисление (K_noisy^-1) matrix K_noisy_inv = cho_solve(result.L_K_noisy, matrix::Identity(n, n)); // 1. Градиенты по гиперпараметрам ядра vector kernel_hyperparams = kernel.GetHyperparameters(); // Вычисляем aatK_noisy_inv один раз до цикла - O(N^2) matrix aatK_noisy_inv = result.alpha.Outer(result.alpha) - K_noisy_inv; matrix dK_dtheta; for(int i = 0; i < (int)kernel_hyperparams.Size(); i++) { dK_dtheta = kernel.ComputeDerivative(i); // O(N^2) result.nlml_gradient[current_grad_idx] = -0.5 * DiagonalTrace(aatK_noisy_inv, dK_dtheta); current_grad_idx++; } // 2. Градиент по гиперпараметру шума // dNLML/d(sigma^2) = 0.5 * (trace(K_noisy^-1) - alpha^T * alpha ) double dNLML_d_sigma2n = 0.5 * (K_noisy_inv.Trace() - result.alpha @ result.alpha); result.nlml_gradient[current_grad_idx] = dNLML_d_sigma2n * (2.0 * sigma); result.success = true; } string GetName() const override { return "ExactInference"; } };
A classe ExactInference executa a inferência exata em PG e é aplicável apenas quando se utiliza a verossimilhança gaussiana. Sua implementação se concentra no cálculo eficiente de NLML e de seus gradientes analíticos. Os gradientes analíticos, calculados diretamente a partir das fórmulas, garantem maior precisão e aceleram significativamente a otimização em comparação com métodos numéricos.
O algoritmo de inferência consiste em várias etapas principais:
- Geração da matriz de covariância Knoisy: primeiro, calcula-se a matriz de covariância do kernel K para os dados de treinamento. Em seguida, adiciona-se aos elementos diagonais de K a variância do ruído (variance = sigma * sigma) e um pequeno valor de jitter positivo (1e-6). Assim é formada a matriz Knoisy;
- Cálculo do vetor αlpha;
- Cálculo do logaritmo negativo da verossimilhança marginal (NLML);
- Cálculo dos gradientes de NLML para otimização.
Para calcular os gradientes necessários e o próprio NLML, são necessários a inversa da matriz de covariância (K+σ2I)^−1 e o vetor α=(K+σ2I)^−1 * y.
A função cho_solve ajuda a calcular essas variáveis com a máxima rapidez.
Classe LaplaceInference
LaplaceInference implementa a aproximação de Laplace, um método de inferência aproximada aplicável tanto a tarefas de regressão quanto de classificação em PG.
//+------------------------------------------------------------------+ //| LaplaceInference | //+------------------------------------------------------------------+ class LaplaceInference : public IInference { private: int m_max_iterations; double m_tolerance; public: LaplaceInference(int max_iter = 100, double tolerance = 1e-10) : m_max_iterations(max_iter), m_tolerance(tolerance) {} void Infer(const matrix &X_train, const vector &y, IKernel *kernel, ILikelihood *likelihood, GPInferenceResult &result) override { result.success = false; int n = (int)y.Size(); // Вычисление K matrix K = kernel.Compute(X_train, X_train); double jitter = 1e-6; K = K + matrix::Identity(n, n) * jitter; vector f = (result.mu_f_train.Size() == n) ? result.mu_f_train : vector::Zeros(n); bool converged = false; matrix W(n, n); // W := -Hessian matrix L_B(n, n); // L := cholesky(I + W^0.5*K*W^0.5), B = I + W^0.5*K*W^0.5 vector b(n); // b := Wf + grad loglike p(y|f) vector a(n); // a := b - W^0.5*L^T\(L\(W^0.5*K*b)) matrix sW = matrix::Zeros(n, n); // W^0.5 vector sW_diag(n); // вектор для хранения диагональных элементов sqrt(W) double prev_lml_value = -DBL_MAX; // Для критерия сходимости по LML double current_lml_value = -DBL_MAX; int iter; // --- Шаг 1: Нахождение моды f_hat методом Ньютона (Согласно Алгоритму 3.1 GPML) --- for(iter = 0; iter < m_max_iterations; iter++) { // 1. W := - Hessian (диагональная матрица) W = -1 * likelihood.LogLikelihoodHessian(f, y); // Вычисление sW_diag_vector (корень из W) sW_diag = MathSqrt(W.Diag()); // 2. L_B := cholesky(I + W^0.5*K*W^0.5) // Вычисляем sW_K = sW * K result.sW_K.Resize(n, n); for(int i = 0; i < n; i++) { result.sW_K.Row(sW_diag[i] * K.Row(i),i); // K[i,:] * sW_diag_vector[i] } // Вычисляем B = I + sW_K * sW matrix temp_sWKsW(n, n); for(int j = 0; j < n; j++) { temp_sWKsW.Col(result.sW_K.Col(j)*sW_diag[j],j) ; // sW_K[:,j] * sW_diag_vector[j] } matrix B = matrix::Identity(n, n) + temp_sWKsW; if(!B.Cholesky(L_B)) { PrintFormat("Error:Cholesky decomposition of B failed. Error Code: %d", GetLastError()); result.nlml_value = DBL_MAX; return; } result.L_B = L_B; // Сохраняем для использования в прогнозах // 3. b := W*f + grad loglike p(y|f) b = W.Diag() * f + likelihood.LogLikelihoodGradient(f, y); // 4. a := b - W^0.5*L_B^T\(L_B\(W^0.5*K*b)) a = b - sW_diag*cho_solve(L_B,result.sW_K @ b); // 5. f := K * a (новый ньютоновский шаг) f = K @ a; // --- Вычисление текущего NLML для критерия сходимости --- double log_likelihood_at_f = likelihood.LogLikelihood(f, y); double sum_log_diag_L_B = MathLog(L_B.Diag()).Sum(); current_lml_value = -0.5 * a @ f + log_likelihood_at_f - sum_log_diag_L_B; result.nlml_value = -current_lml_value; // 6. Проверка сходимости double lml_change = current_lml_value - prev_lml_value; // PrintFormat("Iter %d: LML = %g, delta LML = %g",iter, current_lml_value, lml_change); // Проверка сходимости по LML if((lml_change < m_tolerance)) { // PrintFormat("Converged at iter %d: LML = %g, delta LML = %g",iter, current_lml_value, lml_change); converged = true; break; } prev_lml_value = current_lml_value; } if(!converged) { PrintFormat("Warning: Newton algorithm didn't converge at iteration %d. Final LML: %g (delta: %g, tolerance: %g)", iter, current_lml_value, (current_lml_value - prev_lml_value), m_tolerance); } //---Сохраняем в структуру GPInferenceResult result.mu_f_train = f; // Апостериорное среднее скрытой функции в обучающих точках result.H = likelihood.LogLikelihoodHessian(f, y); // Гессиан в точке f_hat // --- Вычисления аналитических градиентов NLML Алгоритм 5.1 GPML --- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // ============================================================================== // 1. Градиенты по гиперпараметрам ЯДРА // ============================================================================== // Вычисление R := W^0.5*L_B^-T*(L_B^-1*W^0.5) result.sW.Diag(sW_diag); matrix R(n,n); matrix cs = cho_solve(L_B,result.sW); for(int i = 0; i < n; i++) { R.Row(sW_diag[i] * cs.Row(i),i); } // Вычисление C := L_B^-1 * W^0.5 * K // Для этого решаем линейное уравнение: L_B * C = sW * K относительно C matrix C; if(!L_B.LinearEquationsSolutionTriangular(EQUATIONSFORM_N,result.sW_K, C)) { PrintFormat("Error: LinearEquationsSolutionTriangular L_B * C = W^0.5 * K.Error Code: %d", GetLastError()); result.nlml_value = DBL_MAX; return; } // Вычисление A = Sigma_f_train = (K^−1+W)^−1 = K - C^T @ C // Sigma_f_train - апостериорная ковариационная матрица скрытой функции f на обучающих данных matrix CTC = C.Transpose() @ C; result.Sigma_f_train = K - CTC; //------------------- grad NLML = -(s1 + s2^T*s3) // Вычисление s2 (первая неявная часть) // s2 := -0.5 * diag( diag(K) - diag(C^T*C) ) * ThirdDerivative loglike по f vector third_deriv = likelihood.LogLikelihoodThirdDerivative(f, y); //matrix diag; //diag.Diag(K.Diag() - CTC.Diag()); //vector s2 = -0.5 * diag @ third_deriv; vector s2 = -0.5 * (result.Sigma_f_train.Diag() * third_deriv); // Предварительно вычислим градиент логарифма правдоподобия vector log_lik_gradient = likelihood.LogLikelihoodGradient(f, y); int num_kernel_hyperparameters = kernel.GetNumHyperparameters(); for(int j = 0; j < num_kernel_hyperparameters; j++) { // C2 := dK_dtheta_j (производная матрицы ядра по текущему гиперпараметру) matrix dK_dtheta_j = kernel.ComputeDerivative(j); ///------------------------------------------------------------------------------------- // s1 := 0.5*a^T*C2*a - 0.5 *trace (R*C2) // явная часть производной double s1_term1 = 0.5 * (a @(dK_dtheta_j @ a)); // 0.5 * a^T * dK_dtheta_j * a // matrix R_dK_dtheta_j = R @ dK_dtheta_j; // double s1_term2 = -0.5 * R_dK_dtheta_j.Trace(); // -0.5 * trace(R * dK_dtheta_j) double s1_term2 = -0.5 * DiagonalTrace(R,dK_dtheta_j); // более эффективный вариант double s1_explicit_part = s1_term1 + s1_term2; ///-------------------------------------------------------------------------------------- // b := C2 * grad loglike vector b = dK_dtheta_j @ log_lik_gradient; // s3 := b - K * R * b // вторая неявная часть vector s3 = b - K @(R @ b); double implicit_part = s2 @ s3; // неявная часть производной // NLML result.nlml_gradient[current_grad_idx] = -(s1_explicit_part + implicit_part); current_grad_idx++; } // ============================================================================== // 2. Градиент по гиперпараметрам ПРАВДОПОДОБИЯ dLogLikelihood/d(param_j) // ============================================================================== if(likelihood.GetNumHyperparameters() > 0) { vector likelihood_hyperparams = likelihood.GetHyperparameters(); for(int j = 0; j < (int)likelihood_hyperparams.Size(); j++) { // Term 1: Производная log p(y|f*) по sigma double dlogpy_d_param_j = likelihood.LogLikelihoodGradientParam(f, y, j); // (метод возвращает dHessian/d(param_j) = d^3(log p)/d(param_j)d(f)^2) matrix dH_param_j = likelihood.LogLikelihoodHessianDerivative(f, y, j); // Преобразуем dH_d_param_j в dW_d_param_j (где W = -H) matrix dW_d_param_j = -1.0 * dH_param_j; double term2_lik = -0.5*DiagonalTrace(result.Sigma_f_train,dW_d_param_j); result.nlml_gradient[current_grad_idx] = -(dlogpy_d_param_j + term2_lik); current_grad_idx++; } } result.success = true; } string GetName() const override { return "LaplaceInference"; }
Construtor LaplaceInference(int max_iter = 100, double tolerance = 1e-10): o construtor da classe permite inicializar os parâmetros do método iterativo de Newton, usado para encontrar a moda da distribuição a posteriori.
- max_iter: número máximo de iterações; ao atingir esse limite, o algoritmo será interrompido, mesmo que a convergência não tenha sido alcançada.
- tolerance: critério de convergência. O algoritmo para de iterar se a variação de NLML entre etapas consecutivas ficar abaixo desse limiar.
O método Infer implementa dois algoritmos fundamentais do livro "Gaussian Processes for Machine Learning" (GPML):
- Algoritmo 3.1: busca da moda da distribuição a posteriori e cálculo da NLML,
- Algoritmo 5.1: cálculo dos gradientes de NLML em relação aos hiperparâmetros.
Fig.1. Algoritmo 3.1 para busca da moda f_hat e LML
- O processo começa com o cálculo da matriz de covariância do kernel K. Como aproximação inicial de f, usa-se o resultado mu_f_train da iteração anterior de otimização dos hiperparâmetros para obter uma convergência mais rápida e estável da iteração anterior de otimização dos hiperparâmetros, em vez de começar a cada vez com um vetor nulo.
- Usando o método de Newton, f_hat é atualizado iterativamente.
- Cálculo da NLML: após a convergência da moda f (quando a variação de LML fica menor que m_tolerance), passamos o valor de NLML ao otimizador.
Fig.2. Algoritmo 5.1 para cálculo dos gradientes de LML
- Os gradientes em relação aos hiperparâmetros do kernel incluem o cálculo das matrizes auxiliares R e C. A matriz R é calculada usando cho_solve, e a matriz C, usando LinearEquationsSolutionTriangular.
- Calcula-se a matriz de covariância a posteriori da função latente para os dados de treinamento Σf_train = K−C^TC.
- O gradiente em relação a cada hiperparâmetro do kernel (dK_dtheta_j) é a soma da parte explícita (s1) e das partes implícitas (s2, s3). O cálculo de s1é otimizado com DiagonalTrace, e s3 também usa cho_solve para resolver sistemas.
- Os gradientes em relação aos hiperparâmetros da verossimilhança são calculados com o uso das derivadas da log-verossimilhança e de sua Hessiana em relação aos hiperparâmetros correspondentes. Para isso, são usados os métodos LogLikelihoodGradientParam e LogLikelihoodHessianDerivative do objeto likelihood.
Agora que temos todos os componentes básicos (kernels, verossimilhanças, métodos de inferência), podemos demonstrar como a biblioteca funciona.
Teste da biblioteca com dados sintéticos
O script Gpsynthetic.mq5 nos ajudará a testar nossa biblioteca com dados sintéticos simples. Isso permitirá verificar se os métodos implementados funcionam e se estão corretos.
// --- Подключаем библиотеку GP --- #include <GP/GP.mqh> enum IntervalType { INTERVAL_F = 0, // Доверительный интервал скрытой функции f* INTERVAL_Y = 1 // Доверительный интервал наблюдений y* }; enum Type_inference { Exact = 0, // Exact Inference Laplace = 1 // Laplace Inference }; enum Type_Data { Regression = 0, // Regression Classification = 1 // Classification }; //--- Входные параметры input IntervalType interval_type = INTERVAL_F; // Тип интервала для отображения input Type_Data DataType = Classification; input Type_inference inf = Laplace ; // Тип инференса //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { string InpFileName1; string InpFileName2; string InpFileName3; string InpFileName4; if(DataType == Regression) { // Данные для регрессии InpFileName1 = "Data_Regression/X_train.csv"; InpFileName2 = "Data_Regression/Y_train.csv"; InpFileName3 = "Data_Regression/X_test.csv"; InpFileName4 = "Data_Regression/Y_test.csv"; } else { InpFileName1 = "Data_Classification/X.csv"; // 200 InpFileName2 = "Data_Classification/y.csv"; InpFileName3 = "Data_Classification/X_star.csv"; InpFileName4 = "Data_Classification/y_star.csv"; } //----------- Dataset matrix x_train, y_train, x_test, y_test; CSVtoMatrix(InpFileName1,x_train); CSVtoMatrix(InpFileName2,y_train); CSVtoMatrix(InpFileName3,x_test); CSVtoMatrix(InpFileName4,y_test); // Создаем векторы меток из матриц vector y_train_ = y_train.Col(0); vector y_test_ = y_test.Col(0); //--- 1. Создание объектов ядер IKernel* rbf = new RBFKernel(1,1); // IKernel* linear = new LinearKernel(1.0); // IKernel* periodic = new PeriodicKernel(1.0, 1.0, 5.8); //--- 2. Комбинация ядер (создание составного ядра) // IKernel* functional_kernels_array[] = {rbf, linear, periodic}; // SumKernel* combined_functional_kernel = new SumKernel(functional_kernels_array); //--- 3. Создание объекта правдоподобия ILikelihood* likelihood = NULL; // Объявляем указатель базового типа if(DataType == Regression) { likelihood = new GaussianLikelihood(1); // Присваиваем объект } if(DataType == Classification) { likelihood = new LogitLikelihood(); } //--- 4. Создание объекта инференса IInference* inference; // Объявляем указатель на интерфейс // Выбираем тип инференса: if(inf == Exact) { inference = new ExactInference(); } else { inference = new LaplaceInference(100, 1e-10); } //--- 5. Создание экземпляра модели Гауссовского Процесса CREATE_GP_MODEL(gp_model, rbf, likelihood, inference,x_train, y_train_); // CREATE_GP_MODEL(gp_model, combined_functional_kernel, likelihood, inference,x_train, y_train_); //-------------------------------------------------- ulong start_time_fit = GetMicrosecondCount(); //--- 6. Оптимизация гиперпараметров gp_model.Fit(); ulong end_time_fit = GetMicrosecondCount(); double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; Print(inference.GetName()); Print("Время выполнения Fit: ", StringFormat("%.3f", elapsed_time_ms), " мс"); gp_model.PrintOptimizedKernelParameters(); //--- 7. Предсказание для тестовых точек GPPredictionResult gp_predictions; // создаем структуру куда будут записаны результаты предсказаний ulong start_time_predict = GetMicrosecondCount(); gp_model.Predict(x_test, gp_predictions,PROBIT); // MONTE_CARLO,NUM_INTEGR,PROBIT ulong end_time_predict = GetMicrosecondCount(); elapsed_time_ms = (end_time_predict - start_time_predict) / 1000.0; Print("Время выполнения Predict: ", StringFormat("%.3f", elapsed_time_ms), " мс"); // --- 8. Результаты Классификации if(DataType == Classification) { DisplayClassificationResults(gp_predictions, y_test_); } //--- 9. Результаты Регрессии if(DataType == Regression) { VisualizeGP(x_train, y_train, x_test, y_test_, gp_predictions, gp_model, 15); } //--- 10. Освобождение памяти delete gp_model; }
Dependendo do DataType selecionado (regressão ou classificação), o script define os caminhos para os arquivos CSV com os dados de treinamento e teste. Pressupõe-se que esses arquivos estejam nas pastas Data_Regression ou Data_Classification, na seção Files.
Em seguida, são criados os objetos dos kernels. Para construir uma função de covariância mais complexa, podem ser usados kernels combinados, como SumKernel ou ProductKernel.
Depois, é criado o objeto da função de verossimilhança (ILikelihood) de acordo com o DataType selecionado, bem como o objeto de inferência (IInference).
Com a macro CREATE_GP_MODEL (ou um construtor semelhante), é criado o modelo de PG, que associa o kernel selecionado, a função de verossimilhança e o método de inferência aos dados de treinamento.
O método gp_model.Fit(int maxiter = 20) inicia o treinamento do modelo. O método agora inclui um novo parâmetro, maxiter, que define o número máximo de iterações.
Após a conclusão do treinamento, o método gp_model.Predict() é usado para gerar previsões em novos dados de teste x_test. Os resultados das previsões (valores médios, covariâncias, probabilidades/rótulos) são salvos na estrutura GPPredictionResult.
Para classificação, é possível selecionar o modo de previsão (PROBIT, NUM_INTEGR, MONTE_CARLO).
Em seguida, o script exibe os resultados da previsão:
- para classificação: é chamada a função DisplayClassificationResults(), que calcula a métrica de acurácia e as probabilidades previstas;
- para regressão: a função VisualizeGP() plota os resultados.
Para regressão, são usados dados sintéticos gerados pela função y = sin(x) + 0.5 * x + noise(sigma = 0.1). Não vamos nos aprofundar nessa tarefa, pois já a analisamos no artigo dedicado à regressão. Apenas observe como o tempo de treinamento melhorou depois que passamos a usar gradientes analíticos e métodos da biblioteca OpenBLAS.
Para classificação binária, são usados dados que representam a tarefa XOR (OU exclusivo). Essa é uma tarefa clássica em machine learning, que demonstra a limitação dos modelos lineares e a necessidade de recorrer a abordagens não lineares, como os PG.
Vamos lembrar que XOR é uma operação lógica que recebe duas entradas binárias (0 ou 1) e retorna 1 se exatamente uma das entradas for igual a 1 e 0 caso contrário. Tabela de verdade XOR:
A tarefa de classificação binária XOR consiste em construir um modelo que, com base em duas entradas (A, B), prevê a saída XOR correspondente (0 ou 1). Mas, como nosso modelo trabalha apenas com os rótulos "+1" e "-1", substituiremos o rótulo "0" por "-1".
Os dados para a tarefa XOR são gerados da seguinte forma. São criados pontos de entrada bidimensionais X (X=[x1,x2]), distribuídos aleatoriamente em certo intervalo (por exemplo, [−4,4] em cada eixo). Os rótulos y desses pontos são definidos com base na lógica XOR:
- Se x1 e x2 tiverem o mesmo sinal (ambos positivos ou ambos negativos), então o produto x1⋅x2 será positivo. Nesse caso, o rótulo atribuído será +1.
- Se x1 e x2 tiverem sinais diferentes (um positivo e o outro negativo), então o produto x1⋅x2 será negativo. Nesse caso, o rótulo atribuído será −1.
Assim, os pontos (+,+) e (−,−) pertencem à classe +1, enquanto os pontos (+,−) e (−,+) pertencem à classe −1. Para esta tarefa, usaremos apenas o kernel RBF, pois ele lida muito bem com dependências não lineares, permitindo que os processos gaussianos separem essas classes com eficiência.
Durante os testes, usei 200 observações para treinamento e 100 novos pontos para verificar a acurácia. A acurácia da previsão foi de 98%, o que demonstra a alta eficiência dos PG ao resolver problemas de classificação não linear.
Indicador GPRegressor: Regressão em função do tempo com PG
O indicador é um exemplo de previsão dinâmica de séries temporais financeiras que considera a incerteza. Usando uma janela deslizante para montar a amostra de treinamento e com retreinamento contínuo do modelo, ele implementa o princípio de adaptação às condições de mercado em constante mudança.
//+------------------------------------------------------------------+ //| GPRegressor.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Подключаем библиотеку GP --- #include <GP/GP.mqh> #property indicator_chart_window #property indicator_buffers 3 #property indicator_plots 2 // plot Predicted mean #property indicator_label1 "Predicted mean" #property indicator_type1 DRAW_LINE #property indicator_color1 clrBlue #property indicator_style1 STYLE_SOLID #property indicator_width1 2 //--- plot Confidence interval #property indicator_label2 "Confidence interval" #property indicator_type2 DRAW_FILLING #property indicator_color2 clrLightGray #property indicator_width2 1 //--- Входные параметры input int WindowLength = 10; // Длина окна для обучения модели input int calculate_bars = 100; // Глубина истории для расчета input double zscore = 1.96; // z-значение для интервала (1.96 = 95%) input bool ShowRMSE = false; // Показывать метрику RMSE //--- Переменные-буферы для отрисовки double ExtPredictionBuffer[]; // Буфер для предсказанных значений double ExtUpperBandBuffer[]; // Буфер для верхней границы доверительного интервала double ExtLowerBandBuffer[]; // Буфер для нижней границы //--- Глобальные объекты GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Глобальная матрица для временных индексов X_train matrix g_X_train; //--- Глобальные переменные для накопления данных для RMSE vector gp_pred; vector naive_pred; vector y_true; // Флаг, указывающий, был ли RMSE уже рассчитан (для однократного расчета) bool rmse_calculated = false; string comment; double gp_rmse,naive_rmse; //+------------------------------------------------------------------+ //| Custom indicator initialization function | //+------------------------------------------------------------------+ int OnInit() { if(Bars(Symbol(), Period()) < WindowLength + calculate_bars) { PrintFormat("Ошибка: недостаточно баров для вычислений. Требуется минимум %d. Доступно: %d", WindowLength + calculate_bars, Bars(Symbol(), Period())); return(INIT_FAILED); } rmse_calculated = false; // Сбрасываем флаг при инициализации/переинициализации //--- Инициализируем глобальные векторы для RMSE gp_pred.Resize(calculate_bars-1); naive_pred.Resize(calculate_bars-1); y_true.Resize(calculate_bars-1); //--- indicator buffers mapping SetIndexBuffer(0, ExtPredictionBuffer, INDICATOR_DATA); SetIndexBuffer(1, ExtUpperBandBuffer, INDICATOR_DATA); SetIndexBuffer(2, ExtLowerBandBuffer, INDICATOR_DATA); //--- Создаем объекты ядра, правдоподобия и метода инференса g_kernel = new RBFKernel(1.0, 1.0); if(g_kernel == NULL) { Print("Failed to create RBFKernel object"); return(INIT_FAILED); } g_likelihood = new GaussianLikelihood(0.1); if(g_likelihood == NULL) { Print("Failed to create GaussianLikelihood object"); delete g_kernel; return(INIT_FAILED); } g_inference = new ExactInference(); if(g_inference == NULL) { Print("Failed to create ExactInference object"); delete g_kernel; delete g_likelihood; return(INIT_FAILED); } //--- Инициализируем матрицу признаков X_train (один признак - время) g_X_train = matrix::Zeros(WindowLength, 1); for(int i = 0; i < WindowLength; i++) { g_X_train[i][0] = (double)(i + 1); // X = [1, 2, ..., WindowLength] } //--- Создаем объект GaussianProcess // Передаем g_X_train и пустой y_train vector initial_y_train = vector::Zeros(WindowLength); g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference, g_X_train, initial_y_train); if(g_gp_model == NULL) { Print("Failed to create GaussianProcess object"); delete g_kernel; delete g_likelihood; delete g_inference; return(INIT_FAILED); } return(INIT_SUCCEEDED); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int OnCalculate(const int32_t rates_total, const int32_t prev_calculated, const datetime &time[], const double &open[], const double &high[], const double &low[], const double &close[], const long &tick_volume[], const long &volume[], const int32_t &spread[]) { // Проверка на достаточное количество баров if(rates_total < WindowLength) { Print("Ошибка: недостаточно баров для вычислений. Требуется минимум ", WindowLength, ", доступно: ", rates_total); return(0); } int start; if(prev_calculated == 0) { // Определяем, откуда начинать расчет. // calculate_bars - это глубина истории для расчета. // MathMax(WindowLength, ...) гарантирует, что начнем не раньше, чем есть полное окно данных. start = MathMax(WindowLength, rates_total - calculate_bars); // Инициализируем все буферы значением EMPTY_VALUE ArrayInitialize(ExtPredictionBuffer, EMPTY_VALUE); ArrayInitialize(ExtUpperBandBuffer, EMPTY_VALUE); ArrayInitialize(ExtLowerBandBuffer, EMPTY_VALUE); // Устанавливаем начало отрисовки буферов PlotIndexSetInteger(0, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(1, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(2, PLOT_DRAW_BEGIN, start); } else { if(time[rates_total - 1] == time[prev_calculated - 1]) { // Это означает, что новый бар НЕ появился. // Мы находимся на том же баре, просто пришли новые тики. return(rates_total); } // Если время последнего бара изменилось, значит, появился новый закрытый бар. // Начинаем расчет с prev_calculated, чтобы обработать все новые бары. start = prev_calculated; } for(int i = start; i < rates_total && !IsStopped(); i++) { // Проверка достаточности баров для формирования окна if(i < WindowLength) { Print("Ошибка: недостаточно баров для формирования окна на баре ", i, ". Требуется: ", WindowLength, ", доступно: ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Пропускаем этот бар } // Формирование вектора y из close[i-1], close[i-2], ..., close[i-WindowLength] // close[i-1] - самый новый бар в окне // close[i-WindowLength] - самый старый бар в окне vector y(WindowLength); for(int j = 0; j < WindowLength; j++) { y[WindowLength - 1 - j] = close[i - 1 - j]; } //Print(" y = ", y[WindowLength-1]); // Самый новый бар в окне данных //--- Стандартизация данных double mu_raw = y.Mean(); // Print("mu_raw = ", mu_raw); double sigma_raw = y.Std(); // Print("sigma_raw = ", sigma_raw); vector y_standardized = (y - mu_raw) / sigma_raw; //--- Устанавливаем обучающие данные g_gp_model.SetTrainingData(g_X_train, y_standardized); //--- Сброс начальных гиперпараметров перед Fit() vector initial_params(3); initial_params[0] = 1.0; // sigma_f (RBFKernel) initial_params[1] = 1.0; // length_scale (RBFKernel) initial_params[2] = 0.1; // sigma_n (GaussianLikelihood) g_gp_model.SetHyperparameters(initial_params); // vector params = g_gp_model.GetCurrentHyperparameters(); // Print(params); // ulong start_time_fit = GetMicrosecondCount(); //--- Обучаем модель if(!g_gp_model.Fit()) { Print("Ошибка: не удалось оптимизировать гиперпараметры GP для бара ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Пропускаем текущий бар } //ulong end_time_fit = GetMicrosecondCount(); // double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; //Print("Время выполнения Fit: ", StringFormat("%.3f", elapsed_time_ms), " мс"); vector params = g_gp_model.GetCurrentHyperparameters(); double sigma_f = params[0]; // sigma_f (RBFKernel) double length_scale = params[1]; // length_scale (RBFKernel) double sigma_n = params[2]; // sigma_n (GaussianLikelihood) // --- Создаем точку для прогноза (X_star) // Предсказываем значение для следующего шага после WindowLength matrix X_star = matrix::Zeros(1, 1); X_star[0][0] = (double)(WindowLength + 1); // --- Выполняем предсказание GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Ошибка: не удалось получить прогноз от GP-модели для бара ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; } // --- double predicted_mean_st = gp_predictions.mu_f_star[0]; // --- Выполняем обратное преобразование предсказанного среднего значения // обратно в исходную ценовую шкалу double predicted_mean = predicted_mean_st * sigma_raw + mu_raw; double sigma_f_star = gp_predictions.Sigma_f_star[0,0]; // дисперсия (Σ) сигнала f* double predicted_variance_st = gp_predictions.Sigma_y_star[0,0]; // дисперсия(Σ) наблюдения y* // --- Выполняем обратное преобразование дисперсии // Дисперсия масштабируется квадратом стандартного отклонения (sigma_raw) данных, // которые были использованы для стандартизации. double predicted_variance = predicted_variance_st * sigma_raw * sigma_raw; // Расчет спрогнозированного доверительного интервала double std_predict = MathSqrt(predicted_variance); double upper_bound = predicted_mean + zscore * std_predict; double lower_bound = predicted_mean - zscore * std_predict; // --- Сохраняем результаты в буферы индикатора для отрисовки ExtPredictionBuffer[i] = predicted_mean; ExtUpperBandBuffer[i] = upper_bound; ExtLowerBandBuffer[i] = lower_bound; comment = StringFormat("GP Regressor | Mean=%.5f | Upper=%.5f | Lower=%.5f |\n", predicted_mean, upper_bound, lower_bound); comment+= StringFormat("Sigma_f = %.5f | Length_Scale = %.5f | Sigma_n = %.5f\n",sigma_f,length_scale,sigma_n); comment+= StringFormat("mu f* = %.5f | Σ y* = %.5f | Σ f*= %.5f",predicted_mean_st,predicted_variance_st,sigma_f_star); Comment(comment); // Отладочный вывод Print("Прогноз для бара ", i, ": Mean=", DoubleToString(predicted_mean, Digits()), " Upper=", DoubleToString(upper_bound, Digits()), " Lower=", DoubleToString(lower_bound, Digits())); } // --- Расчет RMSE для оценки эффективности прогнозной модели ГП по сравнению с простым "наивным" прогнозом if(ShowRMSE && !rmse_calculated) { vector returns(calculate_bars-1); for(int j=0; j <calculate_bars-1;j++) { // истинное значение y_true[j] = close[rates_total - calculate_bars + j]; //Print("y_true[j]", y_true[j] ); // GP прогноз: gp_pred[j] = ExtPredictionBuffer[rates_total-calculate_bars + j]; //Print("gp_pred[j]", gp_pred[j] ); // "Наивный" прогноз (цена завтра = цена сегодня): naive_pred[j] = close[rates_total - calculate_bars + j-1]; // Print("naive_pred[j]", naive_pred[j] ); returns[j] = close[rates_total - calculate_bars + j] - close[rates_total - calculate_bars + j-1]; } gp_rmse=gp_pred.RegressionMetric(y_true,REGRESSION_RMSE); naive_rmse=naive_pred.RegressionMetric(y_true,REGRESSION_RMSE); double std_returns = returns.Std(); // стандартное отклонение приращений Close comment = comment + StringFormat("\nRMSE GP: %.5f | RMSE Naive: %.5f | StdDev Returns Close: %.5f ", gp_rmse, naive_rmse, std_returns); rmse_calculated = true; } Comment(comment); //--- Возвращаем rates_total для следующего вызова return(rates_total); } //+------------------------------------------------------------------+ //| Custom indicator deinitialization function | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { //--- Очищаем комментарий на графике Comment(""); //--- Освобождаем выделенную память для объекта GP if(g_gp_model != NULL) { delete g_gp_model; } } //+------------------------------------------------------------------+
Parâmetros de entrada:
- WindowLength: comprimento da janela de dados (em barras) para treinar o modelo de PG.
- calculate_bars: profundidade do histórico para cálculo e renderização do indicador.
- zscore: quantidade de desvios padrão para construir o intervalo de confiança (por exemplo, 1.96 para 95%).
- ShowRMSE: flag para exibir a raiz do erro quadrático médio (RMSE) da previsão.
Inicialização (OnInit):
- Objetos globais. São criados e inicializados os objetos globais do kernel (RBFKernel), da função de verossimilhança (GaussianLikelihood) e do método de inferência (ExactInference).
- Atributos X_train. GPRegressor implementa regressão em função do tempo, usando os índices de tempo (de 1 a WindowLength) dentro da janela de treinamento como atributos de entrada X.
- Modelo GaussianProcess. O objeto g_gp_model é inicializado com os componentes especificados, enquanto y_train fica inicialmente vazio e será preenchido dinamicamente.
Laço principal de cálculo (OnCalculate):
- Geração dos dados: em cada barra, são extraídos os últimos WindowLength valores de fechamento para montar o vetor y.
- Padronização dos dados: y é padronizado (transformado para média zero e variância unitária) a fim de aumentar a estabilidade numérica dos cálculos com PG.
- Atualização e treinamento do modelo:
Os hiperparâmetros são redefinidos para os valores iniciais antes de cada treinamento (g_gp_model.SetHyperparameters).
O modelo GP é treinado na janela atual de preços padronizados com g_gp_model.Fit().
- Previsão: para prever a próxima barra, é criado o ponto X_star com o índice de tempo WindowLength + 1. O método g_gp_model.Predict() gera a previsão, retornando a média (mu_f_star) e a variância (Sigma_y_star) para os dados padronizados.
- Despadronização: a média e a variância previstas são transformadas de volta para a escala original de preços usando a média (mu_raw) e o desvio padrão (sigma_raw) da janela atual.
- Intervalo de confiança: os limites superior e inferior do intervalo são calculados com base na média e no desvio padrão transformados, usando o valor zscore definido pelo usuário para o número de desvios padrão.
Fig.3. Previsão do modelo regressor e 95% do intervalo de confiança
Cálculo da métrica RMSE:
O indicador exibe uma única vez, a pedido do usuário, as estatísticas de RMSE para as previsões do modelo de PG e da previsão "ingênua" (em que o preço de amanhã = o preço de hoje) para as últimas calculate_bars barras. Essas estatísticas ajudam a avaliar objetivamente a eficiência atual do modelo.
Por exemplo, valores típicos para EURUSD no timeframe de um minuto podem ser: RMSE GP (0.00032) | RMSE Naive (0.00029), enquanto o desvio padrão das variações dos preços de fechamento (StdDev Returns Close) é 0.00029.
O fato de RMSE GP ser maior que RMSE Naive indica desempenho inferior do modelo de PG em comparação com a previsão naive no período analisado. Além disso, a igualdade entre RMSE Naive e StdDev Returns Close é uma característica típica de passeio aleatório, quando as variações futuras de preço não dependem das variações passadas e o valor atual é a melhor previsão.
A partir disso, podemos concluir que, nesse trecho da série temporal, ou não há padrões previsíveis, ou a configuração atual do modelo de PG, que usa apenas o tempo como atributo, ainda não é capaz de extrair informações úteis a ponto de superar a previsão naive.
Indicador GPClassifier: classificação com base em processos gaussianos
O indicador GPClassifier demonstra como criar, treinar e gerar previsões com um modelo de PG em uma tarefa de classificação binária. Como atributos, são usadas as variações dos preços de fechamento, calculadas em uma janela deslizante sobre os dados históricos. A função NormalizeData normaliza a matriz de atributos para estabilizar o treinamento do modelo.
//+------------------------------------------------------------------+ //| GPClassifier.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Подключаем библиотеку GP --- #include <GP/GP.mqh> #property indicator_chart_window #property indicator_buffers 4 #property indicator_plots 2 // plot для предсказанного класса "+1" (стрелки вверх) #property indicator_label1 "Predicted Class +1" #property indicator_type1 DRAW_ARROW #property indicator_color1 clrGreen #property indicator_width1 1 // plot для предсказанного класса "-1" (стрелки вниз) #property indicator_label2 "Predicted Class -1" #property indicator_type2 DRAW_ARROW #property indicator_color2 clrBlack #property indicator_width2 1 //--- Входные параметры input int WindowLength = 10; // Длина окна для обучения модели GP input int NumLags = 1; // Количество признаков input int calculate_bars = 100; // Глубина истории для расчета input double ProbabilityThreshold = 0.5; // Порог вероятности input int ArrowShiftPx = 10; // Смещение стрелок в пикселях input bool ShowACCURACY = false; // Показывать точность на графике //--- Переменные-буферы для отрисовки double ExtClassUpBuffer[]; // Буфер для отображения стрелок UP double ExtClassDownBuffer[]; // Буфер для отображения стрелок DOWN double ExtPredictedClassBuffer[]; // Буфер для хранения предсказанных меток класса double ExtPredictedProbabilityBuffer[]; // Буфер для хранения предсказанных вероятностей //--- Глобальные объекты гауссовского процесса GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Глобальные переменные для метрик классификации bool accuracy_calculated = false; vector gp_pred,naive_pred,y_true,gp_accuracy,naive_accuracy; string comment; int currentsize; //+------------------------------------------------------------------+ //| Custom indicator initialization function | //+------------------------------------------------------------------+ int OnInit() { accuracy_calculated = false; // Сбрасываем флаг при инициализации/переинициализации //--- Инициализируем глобальные векторы для Accuracy gp_pred.Resize(0); naive_pred.Resize(0); y_true.Resize(0); gp_accuracy.Resize(1); naive_accuracy.Resize(1); currentsize = 0; //--- indicator buffers mapping SetIndexBuffer(0, ExtClassUpBuffer, INDICATOR_DATA); // Для стрелок вверх SetIndexBuffer(1, ExtClassDownBuffer, INDICATOR_DATA); // Для стрелок вниз SetIndexBuffer(2, ExtPredictedClassBuffer, INDICATOR_CALCULATIONS); // для хранения предсказанных меток SetIndexBuffer(3, ExtPredictedProbabilityBuffer, INDICATOR_CALCULATIONS); // для хранения вероятностей // --- Устанавливаем символы стрелок PlotIndexSetInteger(0, PLOT_ARROW, 225); // Стрелка вверх для буфера 0 PlotIndexSetInteger(1, PLOT_ARROW, 226); // Стрелка вниз для буфера 1 // --- Устанавливаем смещение стрелок в пикселях --- // Для стрелок вверх (буфер 0), используем отрицательное смещение, чтобы они были над High PlotIndexSetInteger(0, PLOT_ARROW_SHIFT, -ArrowShiftPx); // отрицательное = вверх // Для стрелок вниз (буфер 1), используем положительное смещение, чтобы они были ниже Low PlotIndexSetInteger(1, PLOT_ARROW_SHIFT, ArrowShiftPx); // положительное = вниз // Добавляем проверку на соотношение NumLags и WindowLength if(NumLags <= 0) { Print("Ошибка: NumLags должен быть положительным числом"); return(INIT_FAILED); } if(NumLags >= WindowLength) { Print("Ошибка: NumLags (", NumLags, ") не может быть больше или равен WindowLength (", WindowLength, ")."); return(INIT_FAILED); } //--- Создаем объекты ядра, правдоподобия и метода инференса g_kernel = new RBFKernel(1.0, 1.0); if(g_kernel == NULL) { Print("Failed to create RBFKernel object"); return(INIT_FAILED); } g_likelihood = new LogitLikelihood(); if(g_likelihood == NULL) { Print("Failed to create LogitLikelihood object"); delete g_kernel; return(INIT_FAILED); } g_inference = new LaplaceInference(); if(g_inference == NULL) { Print("Failed to create ExactInference object"); delete g_kernel; delete g_likelihood; return(INIT_FAILED); } // --- Создаем модель GaussianProcess // Передаем пустые матрицы для X_train и y_train. // Они будут переопределяться в OnCalculate. g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference,matrix::Zeros(1,1), vector::Zeros(1)); if(g_gp_model == NULL) { Print("Failed to create GaussianProcess object"); delete g_kernel; delete g_likelihood; delete g_inference; return(INIT_FAILED); } return(INIT_SUCCEEDED); } //+------------------------------------------------------------------+ //| Custom indicator iteration function | //+------------------------------------------------------------------+ int OnCalculate(const int32_t rates_total, const int32_t prev_calculated, const datetime &time[], const double &open[], const double &high[], const double &low[], const double &close[], const long &tick_volume[], const long &volume[], const int32_t &spread[]) { // Проверка на достаточное количество баров if(rates_total < WindowLength + NumLags + 1) { Print("Ошибка: недостаточно баров для вычислений. Требуется минимум ", WindowLength + NumLags + 1, ", доступно: ", rates_total); return(0); } int start; if(prev_calculated == 0) { start = MathMax(WindowLength + NumLags + 1, rates_total - calculate_bars); // Инициализируем все буферы EMPTY_VALUE ArrayInitialize(ExtClassUpBuffer, EMPTY_VALUE); ArrayInitialize(ExtClassDownBuffer, EMPTY_VALUE); // Количество начальных баров без отрисовки и значений в DataWindow PlotIndexSetInteger(0, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(1, PLOT_DRAW_BEGIN, start); } else { if(time[rates_total - 1] == time[prev_calculated - 1]) { return(rates_total); } start = prev_calculated; } // --- Главный цикл расчета прогнозов for(int i = start; i < rates_total && !IsStopped(); i++) { // Проверка достаточности баров для формирования окна X и y if(i < WindowLength + NumLags + 1) { ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } // --- Шаг 1: Формирование обучающих данных (X_train и y_train) matrix X_train = matrix::Zeros(WindowLength, NumLags); vector y_train = vector::Zeros(WindowLength); for(int j = 0; j < WindowLength; j++) { // Расчет индекса бара в массиве close для текущего наблюдения в окне. int idx = i - WindowLength + j; // Формирование признаков (лагов) для текущего наблюдения (строки X_train[j]) for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { // Пример: // lag_idx = 0: (close[idx - 1] - close[idx - 2]) - приращение предыдущего бара // lag_idx = 1: (close[idx - 2] - close[idx - 3]) - приращение позапрошлого бара X_train[j,lag_idx] = close[idx - 1 - lag_idx] - close[idx - 2 - lag_idx]; } // Y_train[j]: Бинарная метка для приращения цены на баре idx. if(close[idx] - close[idx - 1] > 0) { y_train[j] = 1; // Цена выросла } else { y_train[j] = -1; // Цена упала или осталась неизменной } } // Print(y_train); // --- Нормализация X_train matrix X_train_norm; vector out_mean; vector out_std; if(!NormalizeData(X_train,X_train_norm,out_mean,out_std)) { Print("Ошибка: не удалось нормализовать матрицу признаков ", i); continue; } // --- Шаг 2: Подаем данные в GP-модель g_gp_model.SetTrainingData(X_train_norm, y_train); // Сброс начальных гиперпараметров перед каждым обучением Fit() vector initial_params(2); initial_params[0] = 1.0; // sigma_f (RBFKernel) initial_params[1] = 1.0; // length_scale (RBFKernel) g_gp_model.SetHyperparameters(initial_params); // --- Шаг 3: Обучаем модель // ulong start_time_fit = GetMicrosecondCount(); if(!g_gp_model.Fit()) { Print("Ошибка: не удалось оптимизировать гиперпараметры GP для бара ", i); ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } // ulong end_time_fit = GetMicrosecondCount(); // double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; // Print("Время выполнения Fit: ", StringFormat("%.3f", elapsed_time_ms), " мс"); // --- Шаг 4: Создаем точку для прогноза X_star // X_star - это приращения цен предыдущих NumLags закрытых баров, // используемые для предсказания движения текущего бара (бара i). matrix X_star = matrix::Zeros(1, NumLags); for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { X_star[0][lag_idx] = close[i - 1 - lag_idx] - close[i - 2 - lag_idx]; } // X_star также должен быть нормализован, используя те же средние и стандартные отклонения, // что и для соответствующих столбцов X_train. for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { X_star[0][lag_idx] = (X_star[0][lag_idx] - out_mean[lag_idx]) / out_std[lag_idx]; } // --- Шаг 5: Выполняем предсказание GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Ошибка: не удалось получить прогноз от GP-модели для бара ", i); ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } double predicted_probability = gp_predictions.predicted_probabilities[0]; // Вероятность для единственной тестовой точки int predicted_class = (int)gp_predictions.predicted_labels[0]; // Метка (+1 или -1) // --- Шаг 6: Сохраняем результаты в буферы индикатора для отрисовки ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; ExtPredictedClassBuffer[i] = predicted_class; ExtPredictedProbabilityBuffer[i] = predicted_probability; // Добавляем вывод вероятности и количество признаков (лагов) в комментарий на график comment = StringFormat("GP Classifier (Lags: %d) | Probability(UP) = %.10f", NumLags, predicted_probability); Comment(comment); // Если прогнозируемый класс - "+1" и вероятность выше заданного порога if(predicted_class == 1 && predicted_probability > ProbabilityThreshold) { ExtClassUpBuffer[i] = high[i]; } // Если прогнозируемый класс - "-1" и вероятность падения выше заданного порога else if(predicted_class == -1 && (1.0 - predicted_probability) > ProbabilityThreshold) { ExtClassDownBuffer[i] = low[i]; } Print("Прогноз для бара ", i, ": Probability(UP)=" + DoubleToString(predicted_probability, 10) + ", Predicted Class=" + IntegerToString(predicted_class), " time: ", time[i]); } // --- Расчет точности классификации для оценки эффективности прогнозной модели ГП по сравнению //с "наивным" прогнозом. if(ShowACCURACY && !accuracy_calculated) { for(int j=0; j <calculate_bars-1;j++) { // индекс до текущего бара int bar_idx = rates_total - calculate_bars + j; if(ExtClassUpBuffer[bar_idx]!= EMPTY_VALUE || ExtClassDownBuffer[bar_idx] != EMPTY_VALUE) { currentsize = (int)gp_pred.Size(); currentsize++; gp_pred.Resize(currentsize,100); naive_pred.Resize(currentsize,100); y_true.Resize(currentsize,100); // Истинная метка: считаем до текущего бара if(close[bar_idx] - close[bar_idx - 1] > 0) { y_true[currentsize-1] = 1; // Цена выросла } else { y_true[currentsize-1] = -1; // Цена упала или осталась неизменной } // Print("y_true", y_true[currentsize-1] ); // Прогноз GP: считаем до текущего бара gp_pred[currentsize-1] = ExtPredictedClassBuffer[bar_idx]; // Print(" gp_pred", gp_pred[currentsize-1]); // "Наивный" прогноз: метка бара "на завтра" равна метке "сегодня" if(close[bar_idx - 1] - close[bar_idx - 2] > 0) { naive_pred[currentsize-1] = 1; } else { naive_pred[currentsize-1] = -1; } } } if(!(int)gp_pred.Size()==0) { gp_accuracy=gp_pred.ClassificationMetric(y_true,CLASSIFICATION_ACCURACY); naive_accuracy=naive_pred.ClassificationMetric(y_true,CLASSIFICATION_ACCURACY); comment += StringFormat("\nGP Accuracy: %.2f %% | Naive Accuracy: %.2f %%", gp_accuracy[0], naive_accuracy[0]); comment += StringFormat("\nNumber of Filtered Signals: %d", (int)gp_pred.Size()); } else { comment += StringFormat("\n No Signals for ProbabilityThreshold: %.4f ", ProbabilityThreshold); } accuracy_calculated = true; } Comment(comment); //--- Возвращаем rates_total для следующего вызова return(rates_total); } //+------------------------------------------------------------------+ //| Custom indicator deinitialization function | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { // Очищаем комментарий на графике при удалении индикатора Comment(""); if(g_gp_model != NULL) { delete g_gp_model; } }
Parâmetros de entrada:
- WindowLength: comprimento da janela de dados para treinar o modelo. Número de observações na amostra de treinamento.
- NumLags: quantidade de lags (variações de preço defasadas). Por exemplo, se NumLags = 1, o atributo será apenas a variação do preço na barra anterior. Se NumLags = 3, os atributos serão as variações do preço nas três barras anteriores, e assim por diante.
- calculate_bars: define o número das últimas barras no gráfico para as quais o indicador será calculado e plotado.
- ProbabilityThreshold: limiar de probabilidade. O indicador exibe sinais com setas para cima ou para baixo se a probabilidade prevista de movimento do preço na direção correspondente exceder esse limiar.
- ArrowShiftPx: deslocamento das setas em pixels em relação à máxima e à mínima da barra.
- ShowACCURACY: flag para exibir as métricas de acurácia da classificação no gráfico.
OnInit(): inicialização do modelo de PG:
- Criamos o kernel RBF, a função de verossimilhança e o método de inferência:
- g_kernel = new RBFKernel(1.0, 1.0)
- g_likelihood = new LogitLikelihood()
- g_inference = new LaplaceInference()
- Criamos o modelo GaussianProcess:
- g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference, matrix::Zeros(1,1), vector::Zeros(1));
- O objeto g_gp_model é criado com estruturas iniciais vazias para os dados de treinamento X_train e y_train. Esses dados serão preenchidos e atualizados dinamicamente a cada barra dentro da função OnCalculate.
OnCalculate(): laço principal de cálculo:
1. Verificação de barras suficientes: antes de iniciar os cálculos, verificamos se há barras suficientes no gráfico para formar a janela de treinamento, considerando WindowLength, NumLags e uma barra para o rótulo a ser previsto.
2. Geração dos dados de treinamento (X_train e y_train):
- Atributos X_train: a cada etapa, a matriz de atributos X_train é montada em uma janela deslizante. Cada linha é uma observação, e as colunas são as variações de preço (lags) nas NumLags barras anteriores.
- Rótulos y_train: rótulos binários que indicam a direção do movimento do preço na próxima barra: "1" (o preço subiu) ou "-1" (o preço caiu ou não se alterou).
3. Normalização de X_train: a matriz de atributos é padronizada para aumentar a estabilidade numérica. Os valores médios (out_mean) e os desvios padrão (out_std) obtidos nessa etapa são salvos para normalizar o ponto de previsão.
4. Atualização e treinamento do modelo de PG:
- g_gp_model.SetTrainingData(X_train_stdard, y_train): os dados de treinamento são atualizados a cada barra.
- g_gp_model.SetHyperparameters(initial_params): os hiperparâmetros do kernel são redefinidos e otimizados novamente para cada nova barra.
- g_gp_model.Fit(): executa o treinamento na janela atual de dados padronizados.
5. Criação e normalização do novo ponto de previsão (X_star): o ponto é normalizado usando os vetores out_mean e out_std.
6. Execução da previsão: o método g_gp_model.Predict(X_star, gp_predictions) retorna a probabilidade prevista de alta (predicted_probabilities[0]) e o rótulo binário (predicted_labels[0]).
7. Desenho e exibição: dependendo de predicted_class e ProbabilityThreshold, uma seta para cima ou para baixo é desenhada no gráfico. As informações sobre a probabilidade da previsão e a classe prevista são exibidas no comentário do gráfico e no log.
Fig.4. Resultado da previsão do indicador GP Classifier (kernel RBF)
Cálculo das métricas de acurácia (ShowACCURACY)
Quando o parâmetro ShowACCURACY é ativado, o indicador executa um único cálculo da acurácia da classificação para as últimas calculate_bars barras, comparando as previsões do modelo de PG com a previsão "ingênua".
Os rótulos reais (y_true) são definidos como a direção efetiva do movimento do preço. A previsão do modelo de PG (gp_pred) corresponde aos rótulos previstos pelo modelo de PG. A previsão "ingênua" (naive_pred) assume que a direção do próximo movimento do preço será a mesma do movimento anterior, por exemplo, se o preço subiu na barra anterior, a previsão ingênua será de alta.
São comparadas as métricas CLASSIFICATION_ACCURACY de ambos os modelos (gp_accuracy e naive_accuracy), e os resultados são exibidos no comentário do gráfico. Isso permite avaliar objetivamente o desempenho do modelo de PG em relação a um benchmark simples.
Considerações finais
Hoje concluímos a análise do modelo de classificação com processos gaussianos. Depois de examinar em detalhes os fundamentos teóricos, conseguimos criar uma versão básica da biblioteca que permite resolver tanto problemas de regressão quanto de classificação. O teste com dados sintéticos confirmou que os componentes implementados da biblioteca funcionam corretamente, incluindo a aplicação eficiente de gradientes analíticos e da biblioteca de alto desempenho OpenBLAS para estabilidade numérica e velocidade.
Os indicadores GPRegressor e GPClassifier demonstraram as possibilidades práticas de uso dos PG para a previsão dinâmica de séries temporais financeiras em tempo real, e métricas como RMSE (para regressão) e Accuracy (para classificação) permitiram avaliar objetivamente a eficiência do modelo.
Ainda assim, a implementação da biblioteca está longe do ideal. Apesar de todos os esforços empreendidos, sua velocidade de execução é inferior à de soluções consolidadas, como o scikit-learn. Isso se deve em parte ao fato de termos usado o otimizador MinBleic, que não é o mais rápido.
Portanto, as melhorias futuras da biblioteca podem passar por:
- uso de um otimizador baseado em gradientes mais eficiente. Em particular, precisaremos de um algoritmo L-BFGS rápido, que supera significativamente o atual em eficiência para tarefas com grande número de hiperparâmetros;
- implementação mais robusta do método de Newton para a aproximação de Laplace. A adição de busca linear permitirá encontrar a moda da distribuição a posteriori de forma mais eficiente e confiável, melhorando a convergência da otimização resultante;
- implementação de métodos esparsos. Diferentemente da abordagem atual com matrizes completas, o uso de representações esparsas permitirá treinar modelos com volumes de dados significativamente maiores, abrindo novas possibilidades de escalabilidade e aplicação dos PG.
Programas utilizados no artigo
| # | Nome | Tipo | Descrição |
|---|---|---|---|
| 1 | GPClassifier.mq5 | Indicador | Classifica a direção do preço um passo à frente com base em probabilidades previstas (treinamento online) |
| 2 | GPRegressor.mq5 | Indicador | Prevê o preço e o intervalo de confiança um passo à frente (treinamento online) |
| 3 | GPsynthetic.mq5 | Script | Verificação do modelo de processo gaussiano com dados sintéticos |
| 4 | GP.mqh | Biblioteca de classes | Classe GaussianProcess, que combina o kernel, a função de verossimilhança e o método de inferência; classe GPOptimizationObjective, responsável pela otimização |
| 5 | Inference.mqh | Biblioteca de classes | Interface IInference e suas implementações, que definem os métodos de inferência |
| 6 | Kernels.mqh | Biblioteca de classes | Interface IKernel e implementações dos kernels de covariância |
| 7 | Likelihoods.mqh | Biblioteca de classes | Interface ILikelihood e implementações das funções de verossimilhança |
| 8 | StructUtils.mqh | Biblioteca de classes | Funções auxiliares e estruturas de dados |
| 9 | DataRegression | CSV | Dados para o script GPsynthetic.mq5 |
| 10 | DataClassification | CSV | Dados para o script GPsynthetic.mq5 |
Traduzido do russo pela MetaQuotes Ltd.
Artigo original: https://www.mql5.com/ru/articles/19013
Arquivos anexados |
MQL5.zip
(53.23 KB)
Aviso: Todos os direitos sobre esses materiais pertencem à MetaQuotes Ltd. É proibida a reimpressão total ou parcial.
Esse artigo foi escrito por um usuário do site e reflete seu ponto de vista pessoal. A MetaQuotes Ltd. não se responsabiliza pela precisão das informações apresentadas nem pelas possíveis consequências decorrentes do uso das soluções, estratégias ou recomendações descritas.
Ir para discussão
Caminhe em novos trilhos: Personalize indicadores no MQL5
Vou agora listar todas as possibilidades novas e recursos do novo terminal e linguagem. Elas são várias, e algumas novidades valem a discussão em um artigo separado. Além disso, não há códigos aqui escritos com programação orientada ao objeto, é um tópico muito importante para ser simplesmente mencionado em um contexto como vantagens adicionais para os desenvolvedores. Neste artigo vamos considerar os indicadores, sua estrutura, desenho, tipos e seus detalhes de programação em comparação com o MQL4. Espero que este artigo seja útil tanto para desenvolvedores iniciantes quanto para experientes, talvez alguns deles encontrem algo novo.
Previsão no trading e modelos Grey
Este artigo aborda a aplicação de modelos Grey à previsão de séries temporais financeiras. Vamos analisar os princípios de funcionamento dos modelos Grey e as particularidades de sua aplicação a séries financeiras. Também discutiremos as vantagens e limitações do uso desses modelos em trading.
Está chegando o novo MetaTrader 5 e MQL5
Esta é apenas uma breve resenha do MetaTrader 5. Eu não posso descrever todos os novos recursos do sistema por um período tão curto de tempo - os testes começaram em 09.09.2009. Esta é uma data simbólica, e tenho certeza que será um número de sorte. Alguns dias passaram-se desde que eu obtive a versão beta do terminal MetaTrader 5 e MQL5. Eu ainda não consegui testar todos os seus recursos, mas já estou impressionado.
MetaTrader 5 Global Optimizer: Uma Estrutura Profissional para Otimizar EAs por Grupos, Subgrupos e Critérios de Robustez
Apresentamos uma metodologia para transformar a otimização de EAs no MetaTrader 5 em um fluxo organizado e auditável. A automação em Python cria .set e .ini, orquestra otimizações por grupos e subgrupos, compara cada etapa ao baseline e aplica rewind quando necessário. O leitor poderá escolher os melhores parâmetros considerando lucro, estabilidade, drawdown, trades, concentração de resultado e consistência em vários ativos.
Você está perdendo oportunidades de negociação:
- Aplicativos de negociação gratuitos
- 8 000+ sinais para cópia
- Notícias econômicas para análise dos mercados financeiros
Registro
Login
Você concorda com a política do site e com os termos de uso
Se você não tem uma conta, por favor registre-se