Procesos gaussianos en el aprendizaje automático (Parte 2): Implementación y prueba de un modelo de clasificación en MQL5
26
4
Introducción
En el artículo anterior nos familiarizamos con los fundamentos teóricos del modelo bayesiano de aprendizaje automático, los procesos gaussianos (GP), y también empezamos a crear una biblioteca de GP en MQL5, describiendo dos clases clave: GaussianProcess y GPOptimizationObjective.
Hoy completaremos la construcción de la biblioteca, analizando en detalle la implementación de las interfaces clave: IKernel, ILikelihood e IInference. A continuación, probaremos la biblioteca con datos sintéticos y crearemos indicadores de clasificación y regresión que muestren su funcionamiento en modo online, con reentrenamiento del modelo en cada nueva barra.
Interfaz IKernel
La interfaz IKernel, que encontrará en el archivo Kernels.mqh, constituye la base para crear núcleos de covarianza en nuestra biblioteca. Esto aporta flexibilidad al sistema y facilita su ampliación: se pueden añadir nuevos tipos de núcleos o combinaciones de ellos sin modificar la estructura principal del código.
interface IKernel { // Calcula la matriz de covarianza entre dos conjuntos de datos virtual matrix Compute(const matrix &X1, const matrix &X2) = 0; // Calcula la derivada de la matriz de covarianza con respecto a un hiperparámetro dado // param_index: índice del hiperparámetro respecto al cual se calcula la derivada (empezando por 0) virtual matrix ComputeDerivative(int param_index) = 0; // Devuelve los valores actuales de todos los hiperparámetros del núcleo virtual vector GetHyperparameters() const = 0; // Establece nuevos valores para los hiperparámetros del núcleo virtual void SetHyperparameters(const vector ¶ms) = 0; // Devuelve el número de hiperparámetros del núcleo virtual int GetNumHyperparameters() const = 0; // Devuelve el nombre del núcleo como cadena virtual string GetName() const = 0; };
La interfaz define los dos métodos más importantes que constituyen la base del funcionamiento de cualquier núcleo:
- Compute(const matrix &X1, const matrix &X2) — calcula la matriz de covarianza K entre dos conjuntos de datos.
- ComputeDerivative(int param_index) — calcula la derivada de la matriz de covarianza con respecto a uno de los hiperparámetros del núcleo. param_index indica el índice del hiperparámetro respecto al cual se calcula la derivada.
Clase RBFKernel (núcleo de base radial)
RBFKernel, también conocido como núcleo gaussiano o núcleo exponencial cuadrático, es uno de los núcleos de covarianza más utilizados. Se caracteriza por dos hiperparámetros: la varianza de la señal σ_f (amplitud) y la escala de longitud l (length), que determina la suavidad de la función.
//+------------------------------------------------------------------+ //| Clase del núcleo RBF | //+------------------------------------------------------------------+ class RBFKernel : public IKernel { private: double length; double sigma_f; int n; // Número de filas en X1; int m; // Número de filas en X2; matrix m_K; matrix m_D_sq; // matriz de cuadrados de distancias ||x_i - x_j||² public: //Constructor RBFKernel(double ls, double sf) : length(ls), sigma_f(sf) {} //Destructor ~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("Error: No se ha logrado calcular 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"; } };
Afortunadamente, en el caso del núcleo RBF, las derivadas se calculan fácilmente:
-
Derivada con respecto a sigma_f:
- Derivada con respecto a l:
Para el parámetro length, realizamos una multiplicación elemento a elemento de la matriz K por (D_sq / length^3), donde D_sq es la matriz de los cuadrados de las distancias euclídeas.
Clase LinearKernel (núcleo lineal)
LinearKernel es un núcleo de covarianza sencillo que supone una dependencia lineal entre los puntos de datos. Se suele utilizar cuando se espera que una función pueda aproximarse mediante un modelo lineal, o como componente de núcleos compuestos más complejos.
El núcleo lineal tiene un hiperparámetro: la varianza de la señal σ_l.
//+------------------------------------------------------------------+ //| Clase para el núcleo lineal | //+------------------------------------------------------------------+ class LinearKernel : public IKernel { private: double sigma_l; int n; int m; matrix m_K; public: //Constructor LinearKernel(double sl): sigma_l(sl){} //Destructor ~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: { // Derivada con respecto a 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):
Clase PeriodicKernel (núcleo periódico)
PeriodicKernel se utiliza para modelar datos con patrones repetitivos o estacionalidad, lo que permite captar dependencias cíclicas.
//+------------------------------------------------------------------+ //| Clase para el núcleo periódico | //+------------------------------------------------------------------+ class PeriodicKernel : public IKernel { private: double sigma_f; // Varianza de la señal (amplitud de las oscilaciones) double length; // Longitud de la escala (suavidad de las oscilaciones) double period; // Periodo de oscilación int n; // Número de filas en X1 int m; // Número de filas en X2 matrix m_K; // Matriz de covarianza matrix m_D; // Matriz D (suma de derivadas con respecto a d [2 * sin^2(M_PI*distance/period] ) matrix m_distance_dim[]; // Array de matrices de diferencias para cada dimensión |x_i_k - x_j_k| int m_d_cols; // Número de dimensiones (características) (X.Cols()) public: //Constructor PeriodicKernel(double ls, double sf, double p) : sigma_f(sf), length(ls), period(p){} //Destructor ~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); // Creamos una matriz en la que cada columna es un vector x1_k // (n × 1) * (1 × m) = (n × m) matrix x1_k_copy = x1_k.Outer(vector::Ones(m)) ; // Creamos una matriz en la que cada fila es el vector x2_k transpuesto // (n × 1) * (1 × m) = (n × m) matrix x2_k_copy = vector::Ones(n).Outer(x2_k); // Calculamos la matriz de todas las diferencias absolutas por pares entre los elementos de los vectores x1_k y x2_k matrix distance = MathAbs(x1_k_copy - x2_k_copy); m_distance_dim[k] = distance; // Almacenamos en la caché distance para cada dimensión (característica) matrix sin_term = MathSin(M_PI * distance / period); sin_term = 2.0 * sin_term * sin_term; // 2 * sin^2(u) m_D += sin_term; // Sumamos elemento a elemento } // Calculamos la matriz 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: { // Derivada con respecto a sigma_f // dK/d(sigma_f) = (2 / sigma_f) * K dK = (2.0 / sigma_f) * m_K; break; } case 1: { // Derivada con respecto a length // dK/d(length) = K * (2 * D / length^3) dK = m_K * ( 2.0 * m_D / (length * length * length) ); break; } case 2: { // Derivada con respecto a 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++) { // Bucle sobre cada dimensión (columna) de los datos matrix distance_k = m_distance_dim[k]; // matriz de diferencias absolutas para la k-ésima dimensión matrix u = M_PI * distance_k / period; // Cálculo del argumento u_k = pi * |x_ik - x_jk| / period matrix sin_2u = MathSin(2.0 * u); // Cálculo de sin(2 * u_k) // Utilizamos el valor de u ya calculado para el cálculo de // -pi * |x_ik - x_jk| / period^2 matrix second_term = (-1*u) / period; // Acumulación de dD/d(period) para la dimensión actual dD_dp += 2.0 * sin_2u * second_term; } // combinación de todas las partes de la derivada // 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 del núcleo periódico respecto al parámetro period:
El cálculo de la derivada con respecto al parámetro period es el más complejo, ya que period se encuentra dentro de una función trigonométrica que, a su vez, forma parte de una función exponencial.
La fórmula de la derivada se basa en la regla de la cadena, donde K depende de D y D depende de period a través de términos sinusoidales. La implementación de ComputeDerivative calcula iterativamente la componente dD_dp para cada dimensión y, a continuación, calcula la derivada final:
Donde dD_dp se calcula como la suma de las derivadas de 2 * sin^2(M_PI*distance/period) en todas las dimensiones (características).
Núcleos compuestos
Los procesos gaussianos permiten combinar núcleos de covarianza simples para crear modelos más complejos, capaces de describir diversos aspectos de los datos (por ejemplo, la tendencia, la periodicidad o el ruido). Esta funcionalidad está implementada mediante dos núcleos compuestos:
- SumKernel: combina varios núcleos sumando sus matrices de covarianza.
- ProductKernel: combina los núcleos mediante la multiplicación elemento a elemento de sus matrices de covarianza.
Ambas clases gestionan los hiperparámetros de sus núcleos hijos, agrupándolos en un único vector común para su optimización y redistribuyéndolos después.
Clase SumKernel
//+------------------------------------------------------------------+ //| Clase para la suma de núcleos | //+------------------------------------------------------------------+ class SumKernel : public IKernel { private: IKernel* kernels[]; public: // Constructor SumKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } //Destructor ~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); // Sumamos las matrices de cada núcleo hijo } return sum; } matrix ComputeDerivative(int param_index) override { int current_param_offset = 0; // Empezamos con un desplazamiento de 0 para el primer núcleo for (int i = 0; i < ArraySize(kernels); i++) { // Obtenemos el número de hiperparámetros del núcleo hijo actual int num_params_current_kernel = kernels[i].GetNumHyperparameters(); // Comprobamos si el param_index global pertenece al núcleo hijo actual // param_index debe ser: // - es mayor o igual que el desplazamiento actual // - es estrictamente menor que el desplazamiento actual + el número de parámetros del núcleo actual if (param_index >= current_param_offset && param_index < current_param_offset + num_params_current_kernel) { // Si param_index pertenece a este núcleo, calculamos su índice local int local_param_index = param_index - current_param_offset; // Llamamos a ComputeDerivative para el núcleo hijo encontrado y // le pasamos el índice local. // Dado que la derivada de la suma de núcleos respecto a un parámetro de un núcleo es igual a la derivada // de este núcleo concreto respecto a este parámetro, podemos devolver el resultado inmediatamente return kernels[i].ComputeDerivative(local_param_index); } // Si param_index no pertenece al núcleo actual, // Aumentamos el desplazamiento para comprobar el siguiente núcleo current_param_offset += num_params_current_kernel; } Print("Error: SumKernel::ComputeDerivative - Parameter index ", param_index, " out of bounds"); return matrix::Zeros(1, 1); } //+-------------------------------------------------------------------+ //| Retorna los valores de los hiperparámetros del núcleo compuesto | //+-------------------------------------------------------------------+ vector GetHyperparameters() const override { vector all_params(GetNumHyperparameters()); int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Obtenemos los valores de los parámetros del núcleo hijo actual vector kernel_params = kernels[i].GetHyperparameters(); for (int j = 0; j < (int)kernel_params.Size(); j++) { // Copiamos los parámetros en un vector común all_params[current_idx + j] = kernel_params[j]; } // Actualizamos el desplazamiento para el siguiente núcleo current_idx += (int)kernel_params.Size(); } return all_params; } //+----------------------------------------------------------------------+ //|Recibe un vector con valores de hiperparámetros y lo descompone | //|en las partes correspondientes y las asigna a cada núcleo hijo | //+----------------------------------------------------------------------+ void SetHyperparameters(const vector ¶ms) override { int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Obtenemos el número de hiperparámetros del núcleo hijo actual int num_params = kernels[i].GetNumHyperparameters(); vector sub_params(num_params); // Copiamos los parámetros correspondientes del vector común for (int j = 0; j < num_params; j++) { sub_params[j] = params[current_idx + j]; } // Llamamos al método SetHyperparameters() en el núcleo hijo actual, // pasándole únicamente sus propios parámetros, reunidos en el vector sub_params kernels[i].SetHyperparameters(sub_params); current_idx += num_params; } } //+---------------------------------------------------------------------------------+ //| Retorna el número total de hiperparámetros del núcleo compuesto 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; } //+------------------------------------------------------------------+ //| Proporciona acceso externo a la lista de núcleos hijos | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { // Copiamos los punteros del array interno kernels al array pasado // array externo output_kernels. Esto proporciona acceso a los núcleos // hijos en el método GaussianProcess::Fit(), donde es necesario recorrer // en un bucle todos los núcleos elementales para establecer los límites de optimización de // sus hiperparámetros. ArrayCopy(output_kernels, kernels); } };
- El método Compute calcula la matriz de covarianza final K como la suma simple de las matrices de covarianza devueltas por cada núcleo hijo. Aquí se aplica el principio del polimorfismo: a pesar de que el array kernels almacena punteros del tipo base IKernel*, la llamada a kernels[i].Compute(X1, X2) invoca, en realidad, la implementación específica de Compute para cada núcleo hijo concreto. Esto permite que SumKernel funcione con cualquier núcleo derivado de IKernel.
- El método ComputeDerivative: la derivada de la función de covarianza de la suma de núcleos respecto a un hiperparámetro concreto es igual a la derivada de la función de covarianza únicamente del núcleo hijo al que pertenece dicho hiperparámetro. Esto significa que ComputeDerivative encuentra el núcleo hijo correspondiente mediante param_index y devuelve únicamente su derivada. Las derivadas con respecto a los parámetros de otros núcleos hijos se consideran iguales a cero.
Clase ProductKernel
//+------------------------------------------------------------------+ //| Clase para el producto de núcleos | //+------------------------------------------------------------------+ class ProductKernel : public IKernel { private: IKernel* kernels[]; matrix m_X1; matrix m_X2; int n; int m; public: // Constructor: copia los punteros a los núcleos hijos ProductKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } // Destructor ~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++) { // Polimorfismo: llamamos al método Compute() de cada núcleo hijo, // independientemente de su tipo concreto, y multiplicamos el resultado elemento a elemento product *= kernels[i].Compute(X1, X2); } return product; } matrix ComputeDerivative(int param_index) override { matrix dK_prod(n, m); // Matriz final de la derivada int current_param_offset = 0; int target_kernel_idx = -1; // Índice del núcleo al que pertenece el hiperparámetro int local_param_index = -1; // Índice local del hiperparámetro en este núcleo // Paso 1: Encontramos el núcleo y el índice local del hiperparámetro 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); } // Paso 2: Calculamos dK_k / d_theta_j para el núcleo objetivo matrix dK_target_kernel = kernels[target_kernel_idx].ComputeDerivative(local_param_index); // Paso 3: Calculamos el producto K_m para todos los demás núcleos (m != k) matrix other_kernels_product = matrix::Ones(n, m); for (int i = 0; i < ArraySize(kernels); i++) { if (i != target_kernel_idx) { // Omitimos el núcleo objetivo other_kernels_product = other_kernels_product * kernels[i].Compute(m_X1, m_X2); } } // Paso 4: Multiplicamos los resultados elemento a elemento dK_prod = dK_target_kernel * other_kernels_product; return dK_prod; } // Recopila los valores de los hiperparámetros de todos los núcleos hijos en un único vector 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; } //+------------------------------------------------------------------+ //| Proporciona acceso externo a la lista de núcleos hijos | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { ArrayCopy(output_kernels, kernels); } };
- Método Compute: calcula la matriz de covarianza final K como el producto elemento a elemento de las matrices de covarianza devueltas por cada núcleo hijo. Al igual que en SumKernel, aquí se utiliza el polimorfismo para llamar al método Compute del núcleo hijo correspondiente.
- Método ComputeDerivative: para calcular la derivada de K con respecto a un hiperparámetro perteneciente a un núcleo hijo concreto K_p, ProductKernel aplica la regla del producto (regla de Leibniz). La derivada de K respecto a este parámetro es igual al producto de la derivada de la matriz K_p respecto a este parámetro por el producto elemento a elemento de las matrices de covarianza de todos los demás núcleos hijos, tomadas sin derivar:
Interfaz ILikelihood
La interfaz ILikelihood define cómo debe funcionar cualquier función de verosimilitud de nuestra biblioteca. Nos permite trabajar con diferentes tipos de datos, ya sean valores continuos (como en la regresión) o categorías discretas (como en la clasificación).
//+------------------------------------------------------------------+ //|Interfaz para funciones de verosimilitud | //+------------------------------------------------------------------+ interface ILikelihood { // Calcula el logaritmo de la verosimilitud p(y|f) para los valores latentes f y los valores observados y virtual double LogLikelihood(const vector &f, const vector &y) = 0; //----------------------------- Derivadas de la log-verosimilitud con respecto a f ----------------- // Calcula el vector de la primera derivada del logaritmo de la verosimilitud con respecto a f (dlp/df) virtual vector LogLikelihoodGradient(const vector &f, const vector &y) = 0; // Calcula la matriz de la segunda derivada (matriz hessiana) del logaritmo de la verosimilitud con respecto a f (d^2lp/df df^T) virtual matrix LogLikelihoodHessian(const vector &f, const vector &y) = 0; // Calcula el vector de la tercera derivada del logaritmo de la verosimilitud con respecto a f (d^3lp/df^3) virtual vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) = 0; //----------------------------- Derivadas de la log-verosimilitud con respecto a los parámetros ----------------- // Primera derivada del logaritmo de la verosimilitud con respecto al j-ésimo hiperparámetro de la verosimilitud. virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) = 0; // Derivada de la matriz hessiana del logaritmo de la verosimilitud con respecto al j-ésimo hiperparámetro de la verosimilitud. virtual matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) = 0; // Devuelve el nombre de la función de verosimilitud virtual string GetName() const = 0; // Devuelve los valores actuales de los hiperparámetros de la verosimilitud virtual vector GetHyperparameters() const = 0; // Establece nuevos valores para los hiperparámetros de la verosimilitud virtual void SetHyperparameters(const vector ¶ms) = 0; // Devuelve el número de hiperparámetros de la verosimilitud virtual int GetNumHyperparameters() const = 0; };
Clase GaussianLikelihood
GaussianLikelihood implementa la función de verosimilitud para problemas de regresión, suponiendo que el ruido en los datos observados sigue una distribución normal.
//+------------------------------------------------------------------+ //| Clase de verosimilitud gaussiana (para problemas de regresión) | //+------------------------------------------------------------------+ class GaussianLikelihood : public ILikelihood { private: double m_noise_sigma; // Parámetro de ruido (desviación estándar) public: // Constructor GaussianLikelihood(double initial_noise_sigma) : m_noise_sigma(initial_noise_sigma) {} // Devuelve el valor actual del hiperparámetro vector GetHyperparameters() const override { vector params(1); params[0] = m_noise_sigma; return params; } // Establece un nuevo valor para el hiperparámetro void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 1) { m_noise_sigma = params[0]; } } // Devuelve el número de hiperparámetros int GetNumHyperparameters() const override { return 1; } // Calcula el logaritmo de la verosimilitud log p(y|f) para una distribución gaussiana 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; // Fórmula del logaritmo de la densidad de una distribución normal multivariante return -0.5 / noise_variance * (residual @ residual) - 0.5 * n * MathLog(2 * M_PI * noise_variance); } //------------------ Derivadas con respecto a f ----------------------- // Calcula el vector de la primera derivada del logaritmo de la verosimilitud con respecto a 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); } // Calcula la matriz de segundas derivadas (Hessiano) del logaritmo de la verosimilitud con respecto a 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)); } // Calcula el vector de terceras derivadas del logaritmo de la verosimilitud con respecto a f vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) override { // Para la verosimilitud gaussiana, d^3(log p(y|f))/df^3 = 0 return vector::Zeros((int)y.Size()); } // ------------------------------ Derivadas con respecto al parámetro ----------------------------------- // // Calcula la primera derivada del logaritmo de la verosimilitud con respecto al j-ésimo hiperparámetro de la verosimilitud virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { // La verosimilitud gaussiana solo tiene un hiperparámetro: m_noise_sigma (índice 0) if (param_index == 0) { // Fórmula para 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; } } // Calcula la derivada del Hessiano del logaritmo de la verosimilitud con respecto al j-ésimo hiperparámetro de verosimilitud matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) override { int n = (int)y.Size(); if (param_index == 0) { // Derivada del Hessiano respecto a 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 clave:
- El constructor GaussianLikelihood inicializa un único hiperparámetro: la desviación estándar del ruido.
- LogLikelihood() — calcula el logaritmo de la verosimilitud p(y∣f) para la densidad de probabilidad de una distribución normal multivariante:
- LogLikelihoodGradient() — calcula la primera derivada del logaritmo de la verosimilitud con respecto a los valores latentes f. Es un vector en el que cada elemento es igual a:
- LogLikelihoodHessian() — calcula la segunda derivada (Hessiano) del logaritmo de la verosimilitud con respecto a f. Para la verosimilitud gaussiana, el Hessiano es una matriz diagonal:
- LogLikelihoodThirdDerivative() — calcula la tercera derivada. Para la verosimilitud gaussiana, todas las derivadas superiores a la segunda son iguales a cero.
- LogLikelihoodGradientParam() — calcula la primera derivada del logaritmo de la verosimilitud con respecto al hiperparámetro m_noise_sigma:
- LogLikelihoodHessianDerivative() — calcula la derivada del Hessiano del logaritmo de la verosimilitud con respecto al hiperparámetro m_noise_sigma:
Clase LogitLikelihood
LogitLikelihood se utiliza para la clasificación binaria, donde los valores objetivo observados y toman los valores {−1, +1}. Vincula la función latente f con la probabilidad de pertenencia a la clase mediante una función sigmoide.
En la implementación de LogitLikelihood se incluyen las funciones auxiliares sigmoid() y Softplus(), con comprobaciones para valores de entrada grandes o pequeños.
//+------------------------------------------------------------------+ //| Clase para la verosimilitud logit (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; // Evitamos la aparición de NaN if (x < -100.0) return 0.0; return 1.0 / (1.0 + MathExp(-x)); } // Función Softplus: log(1 + exp(x)) double Softplus(double x) const { if (x > 100.0) return x; // Para valores muy grandes de x, log(1 + exp(x)) ≈ x if (x < -100.0) return MathExp(x); // Para valores muy pequeños de x, log(1 + exp(x)) ≈ exp(x) return MathLog(1.0 + MathExp(x)); } public: // Constructor (la verosimilitud logit no tiene hiperparámetros propios) LogitLikelihood() {} vector GetHyperparameters() const override { return vector::Zeros(0); } void SetHyperparameters(const vector ¶ms) override {} // Devuelve el número de hiperparámetros de la verosimilitud int GetNumHyperparameters() const override { return 0; } // --- Calcula el logaritmo de la verosimilitud: log p(y|f) --- // Fórmula: 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; } // Calcula el gradiente del logaritmo de la verosimilitud con respecto a f // Fórmula: 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])); de forma equivalente } return grad; } // Calcula el Hessiano del logaritmo de la verosimilitud con respecto a f (matriz diagonal H) // Fórmula: 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; } // Calcula la tercera derivada del logaritmo de la verosimilitud con respecto a f // Fórmula: 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); // de forma equivalente // 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 no tiene hiperparámetros virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { return 0.0; } // LogitLikelihood no tiene hiperparámetros 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 clave:
- LogLikelihood() — calcula el logaritmo de la verosimilitud logp(y∣f). Para la clasificación binaria con y ∈ {−1, +1}, es la suma de los logaritmos de la verosimilitud basada en la entropía cruzada binaria:
- LogLikelihoodGradient() — calcula la primera derivada del logaritmo de la verosimilitud con respecto a f. Cada elemento del vector gradiente es igual a:
- LogLikelihoodHessian() — calcula la segunda derivada (Hessiano) del logaritmo de la verosimilitud con respecto a f. Para LogitLikelihood, el Hessiano es una matriz diagonal cuyos elementos diagonales son iguales a:
- LogLikelihoodThirdDerivative() — calcula la tercera derivada de la log-verosimilitud con respecto a f. Cada elemento del vector es igual a:
Funciones y estructuras auxiliares
El archivo StructUtils.mqh contiene un conjunto de enumeraciones, estructuras de datos y funciones que facilitan el trabajo con los GP.
enum PredictMode { PROBIT = 0, // Aproximación probit NUM_INTEGR = 1, // Integración numérica MONTE_CARLO = 2 // Monte Carlo }; //--- Estructura para los resultados de la predicción struct GPPredictionResult { vector mu_f_star; // Media a posteriori de la función latente f* para datos nuevos matrix Sigma_f_star; // Covarianza a posteriori de la función latente f* para datos nuevos // Campos exclusivos para la regresión (GaussianLikelihood) vector mu_y_star; // Media a posteriori de las observaciones y* matrix Sigma_y_star; // Covarianza a posteriori de las observaciones y* // Campos destinados exclusivamente a la clasificación (LogitLikelihood, ProbitLikelihood, etc.) vector predicted_probabilities; // Probabilidades p(y*=+1 | X*, D) vector predicted_labels; // Etiquetas predichas y* (+1 o -1) }; // Estructura para el resultado de la inferencia struct GPInferenceResult { double nlml_value; // Logaritmo negativo de la verosimilitud marginal vector nlml_gradient; // Gradiente de la NLML respecto a los hiperparámetros del núcleo y de la verosimilitud // Resultados de la inferencia en los datos de entrenamiento: vector mu_f_train; // Media a posteriori de la función latente f en los datos de entrenamiento matrix Sigma_f_train; // Covarianza a posteriori de la función latente f en los datos de entrenamiento (K^−1+W)^−1 matrix L_K_noisy; // Descomposición de Cholesky K_noisy = cholesky(K + s2n*I) matrix L_B; // Descomposición de Cholesky B = I + W^0.5 @ K @ W^0.5 matrix sW; // sW = W^0.5; vector sW_diag; // diagonal de sW matrix sW_K; // sW @ K vector alpha; // (K_noisy)^-1 * y matrix H; // Hessiano del logaritmo de la verosimilitud d^2 log p(y|f)/df^2 (para LaplaceInference) bool success; // Indicador de éxito de la función de inferencia };
- La enumeración PredictMode define los distintos modos para realizar predicciones en un modelo de clasificación.
- La estructura GPInferenceResult se utiliza para almacenar todos los resultados clave obtenidos durante el proceso de inferencia.
- La estructura GPPredictionResult está diseñada para almacenar todos los resultados obtenidos tras realizar una predicción sobre datos nuevos (de prueba) (X_test)
- DiagonalTrace
//+-------------------------------------------------------------------+ //| Calcula la traza del producto de dos matrices usando únicamente | //| elementos diagonales de la matriz final | //+-------------------------------------------------------------------+ 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(); // Calculamos el tamaño mínimo para la diagonal ulong k = MathMin(m, p); double tr = 0.0; for(ulong i = 0; i < k; i++) { // El producto escalar de la i-ésima fila de A y la i-ésima columna de B tr += A.Row(i)@B.Col(i); } return tr; }
Calcula la traza del producto de dos matrices A y B (Tr(A @ B)) sin necesidad de formar toda la matriz producto. Lo hace sumando los productos escalares de las filas de la matriz A y las columnas correspondientes de la matriz B. Esto permite reducir considerablemente los costos computacionales en comparación con la multiplicación completa de matrices.
- cho_solve (envoltorio de LinearEquationsSolutionTriangular de OpenBLAS)
//+------------------------------------------------------------------+ //| Resuelve el sistema A*X = B usando Cholesky A = LL^T | //| Parámetros: | //| c: Matriz triangular inferior L de Cholesky | //| b: miembro derecho del sistema (matriz) | //| Devuelve: matriz X, solución del sistema | //+------------------------------------------------------------------+ matrix cho_solve(const matrix &c, const matrix &b) { // Comprobación de si la matriz «c» es triangular inferior if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular."); return matrix::Zeros(c.Rows(), b.Cols()); } // Paso 1: Resolvemos L * Y = B // Y: matriz intermedia, resultado de calcular 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()); } // Paso 2: Resolvemos L^T * X = Y // X: solución del sistema L^T * X = Y, lo que equivale a (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; } //+------------------------------------------------------------------+ //| Resuelve el sistema Ax = b usando Cholesky A = LL^T | //| Parámetros: | //| c: Matriz triangular inferior L de Cholesky | //| b: miembro derecho del sistema (vector) | //| Devuelve: el vector x, solución del sistema | //+------------------------------------------------------------------+ vector cho_solve(const matrix &c, const vector &b) { // Comprobación de si la matriz c es triangular inferior if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular"); return vector::Zeros(c.Rows()); } // Paso 1: Resolvemos L * y = b // y: vector intermedio, resultado de calcular 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()); } // Paso 2: Resolvemos L^T * x = y // x: solución final del sistema L^T * x = y, lo que equivale a (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; }
La función cho_solve está diseñada para resolver de forma eficaz sistemas de ecuaciones lineales Ax=b (o AX=B, si B es una matriz), donde la matriz A se obtiene a partir de su descomposición de Cholesky A=LL^T. Aquí, L es una matriz triangular inferior.
El cálculo directo de A^−1 (A.Inv()) requiere muchos recursos y puede ser numéricamente inestable. La descomposición de Cholesky, por el contrario, ofrece un enfoque mucho más eficaz y fiable.
Dentro de cho_solve se llevan a cabo dos operaciones principales:
- Sustitución directa: se resuelve el sistema LY=B (o Ly=b), donde L es la matriz de entrada c e Y (o y) es el resultado intermedio. Dado que L es una matriz triangular, la resolución es relativamente rápida.
- Sustitución inversa: a continuación, se resuelve el sistema L^TX=Y (o L^Tx=y), donde L^T es la matriz transpuesta de L. Esto también es eficaz, ya que L^T es una matriz triangular superior.
Estas operaciones se llevan a cabo mediante las funciones LinearEquationsSolutionTriangular de la biblioteca de álgebra lineal de alto rendimiento OpenBLAS. Esto permite acelerar considerablemente los cálculos, lo que los convierte en indispensables para modelos de aprendizaje automático intensivos en recursos, como los procesos gaussianos.
Interfaz IInference
interface IInference { // Este método realiza la inferencia (inferencia de la distribución a posteriori de f) // y devuelve los componentes necesarios para el NLML y la predicción virtual void Infer(const matrix &X, const vector &y, IKernel *kernel, ILikelihood *likelihood,GPInferenceResult &result) = 0; virtual string GetName() const = 0; };
La interfaz permite alternar entre distintos métodos de inferencia, como la inferencia exacta para la verosimilitud gaussiana y los métodos aproximados para casos no gaussianos, sin modificar la lógica principal de GaussianProcess.
El método Infer es central en esta interfaz. Realiza la inferencia de la distribución a posteriori de la función latente f y devuelve los componentes necesarios para los cálculos posteriores (NLML y predicción). Recibe los siguientes argumentos:
- X — matriz de características de entrenamiento,
- y — vector de etiquetas de entrenamiento,
- kernel — puntero a un objeto IKernel,
- likelihood — puntero a un objeto ILikelihood,
- result — referencia a la estructura GPInferenceResult, que se utiliza para almacenar todos los resultados de la inferencia, incluidos:
- el logaritmo negativo de la verosimilitud marginal (NLML), necesario para la optimización de los hiperparámetros,
- el vector de gradientes del NLML con respecto a todos los hiperparámetros del núcleo y de la verosimilitud,
- las matrices y los vectores auxiliares: (L_K_noisy, L_B, sW, mu_f_train, Sigma_f_train, alpha), utilizados para calcular los gradientes de NLML y para la predicción posterior de nuevos datos.
Clase ExactInference
//+------------------------------------------------------------------+ //| ExactInference Solo para la verosimilitud gaussiana | //+------------------------------------------------------------------+ 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(); // Cálculo de K con los parámetros actuales del núcleo matrix K = kernel.Compute(X, X); // Obtenemos la varianza del ruido 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; } //------------- Algoritmo 2.1 GPML --------------------------- //Calculamos 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π) | //+------------------------------------------------------------------+ // --- Calculamos el primer término de la fórmula del NLML: double data_term = 0.5 * (y @ result.alpha); // --- Calculamos el segundo término: 1/2 * log|K + sigma^2*I| = sum(log(L_ii)) double log_det = MathLog(result.L_K_noisy.Diag()).Sum(); // --- Calculamos el tercer término: n/2 * log(2π) double const_term = 0.5 * n * MathLog(2 * M_PI); //--- NLML result.nlml_value = data_term + log_det + const_term; // -------------------- Cálculo de los gradientes del NLML -------------------------------------- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // Cálculo de (K_noisy^-1) matrix K_noisy_inv = cho_solve(result.L_K_noisy, matrix::Identity(n, n)); // 1. Gradientes con respecto a los hiperparámetros del NÚCLEO vector kernel_hyperparams = kernel.GetHyperparameters(); // Calculamos aatK_noisy_inv una vez antes del bucle: 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. Gradiente con respecto al hiperparámetro del ruido // dNLML/d(sigma^2) = 0.5 * (traza(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"; } };
La clase ExactInference está diseñada para realizar inferencia exacta en procesos gaussianos (GP) y solo es aplicable cuando se utiliza la verosimilitud gaussiana. En su implementación, se ha prestado especial atención al cálculo eficiente del NLML y de sus gradientes analíticos. Los gradientes analíticos, calculados directamente a partir de fórmulas, garantizan una mayor precisión y aceleran considerablemente el proceso de optimización en comparación con los métodos numéricos.
El algoritmo de inferencia consta de varios pasos clave:
- Construcción de la matriz de covarianza Knoisy: primero se calcula la matriz de covarianza K del núcleo para los datos de entrenamiento. A continuación, se añade a los elementos diagonales de K la varianza del ruido (variance = sigma * sigma) y un pequeño valor positivo de jitter (1e-6). Así se forma la matriz Knoisy;
- Cálculo del vector αlpha;
- Cálculo del logaritmo negativo de la verosimilitud marginal (NLML);
- Cálculo de los gradientes del NLML para la optimización.
Para calcular los gradientes necesarios y el propio NLML se requieren la inversa de la matriz de covarianza (K+σ²*I)⁻¹ y el vector α = (K+σ²*I)⁻¹ * y.
La función cho_solve permite calcular estas variables con la mayor rapidez posible.
Clase LaplaceInference
LaplaceInference implementa la aproximación de Laplace, un método de inferencia aproximada aplicable tanto a tareas de regresión como de clasificación en GP.
//+------------------------------------------------------------------+ //| 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(); // Cálculo de 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 := -Hessiano 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 log-verosimilitud 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); // vector para almacenar los elementos diagonales de sqrt(W) double prev_lml_value = -DBL_MAX; // Para el criterio de convergencia del LML double current_lml_value = -DBL_MAX; int iter; // --- Paso 1: Búsqueda de la moda f_hat mediante el método de Newton (según el algoritmo 3.1 de GPML) --- for(iter = 0; iter < m_max_iterations; iter++) { // 1. W := -Hessiano (matriz diagonal) W = -1 * likelihood.LogLikelihoodHessian(f, y); // Cálculo de sW_diag_vector (raíz cuadrada de W) sW_diag = MathSqrt(W.Diag()); // 2. L_B := cholesky(I + W^0.5*K*W^0.5) // Calculamos 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] } // Calculamos 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; // Lo guardamos para utilizarlo en las predicciones // 3. b := W*f + grad log-verosimilitud 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 (nuevo paso de Newton) f = K @ a; // --- Cálculo del NLML actual para el criterio de convergencia --- 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. Comprobación de la convergencia double lml_change = current_lml_value - prev_lml_value; // PrintFormat("Iteración %d: LML = %g, delta LML = %g", iter, current_lml_value, lml_change); // Comprobación de la convergencia según el LML if((lml_change < m_tolerance)) { // PrintFormat("Convergió en la iteración %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); } //---Guardamos en la estructura GPInferenceResult result.mu_f_train = f; // Media a posteriori de la función latente en los puntos de entrenamiento result.H = likelihood.LogLikelihoodHessian(f, y); // Hessiano en el punto f_hat // --- Cálculo de los gradientes analíticos del NLML. Algoritmo 5.1 de GPML --- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // ============================================================================== // 1. Gradientes con respecto a los hiperparámetros del NÚCLEO // ============================================================================== // Cálculo de 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álculo de C := L_B^-1 * W^0.5 * K // Para ello, resolvemos la ecuación lineal L_B * C = sW * K con respecto a 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; } // Cálculo de A = Sigma_f_train = (K^−1 + W)^−1 = K - C^T @ C // Sigma_f_train: matriz de covarianza a posteriori de la función latente f en los datos de entrenamiento matrix CTC = C.Transpose() @ C; result.Sigma_f_train = K - CTC; //------------------- grad NLML = -(s1 + s2^T*s3) // Cálculo de s2 (primera parte implícita) // s2 := -0.5 * diag( diag(K) - diag(C^T*C) ) * tercera derivada de la log-verosimilitud respecto de 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); // Calculemos previamente el gradiente del logaritmo de la verosimilitud 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 (derivada de la matriz del núcleo con respecto al hiperparámetro actual) matrix dK_dtheta_j = kernel.ComputeDerivative(j); ///------------------------------------------------------------------------------------- // s1 := 0.5*a^T*C2*a - 0.5 * traza(R*C2) // parte explícita de la derivada 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 * traza(R * dK_dtheta_j) double s1_term2 = -0.5 * DiagonalTrace(R,dK_dtheta_j); // una variante más eficiente double s1_explicit_part = s1_term1 + s1_term2; ///-------------------------------------------------------------------------------------- // b := C2 * grad log-verosimilitud vector b = dK_dtheta_j @ log_lik_gradient; // s3 := b - K * R * b // segunda parte implícita vector s3 = b - K @(R @ b); double implicit_part = s2 @ s3; // parte implícita de la derivada // NLML result.nlml_gradient[current_grad_idx] = -(s1_explicit_part + implicit_part); current_grad_idx++; } // ============================================================================== // 2. Gradiente con respecto a los hiperparámetros de la VEROSIMILITUD dLogLikelihood/d(param_j) // ============================================================================== if(likelihood.GetNumHyperparameters() > 0) { vector likelihood_hyperparams = likelihood.GetHyperparameters(); for(int j = 0; j < (int)likelihood_hyperparams.Size(); j++) { // Término 1: Derivada de log p(y|f*) con respecto a sigma double dlogpy_d_param_j = likelihood.LogLikelihoodGradientParam(f, y, j); // (el método devuelve dHessian/d(param_j) = d³(log p)/d(param_j)d(f)²) matrix dH_param_j = likelihood.LogLikelihoodHessianDerivative(f, y, j); // Convertimos dH_d_param_j en dW_d_param_j (donde 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"; }
Constructor LaplaceInference(int max_iter = 100, double tolerance = 1e-10): el constructor de la clase permite inicializar los parámetros del método iterativo de Newton, que se utiliza para encontrar la moda de la distribución a posteriori.
- max_iter: el número máximo de iteraciones tras las cuales el algoritmo se detendrá, incluso si no se ha alcanzado la convergencia.
- tolerance: umbral de convergencia. El algoritmo detendrá las iteraciones si la variación del valor del NLML entre pasos consecutivos es inferior a este umbral.
El método Infer implementa dos algoritmos clave del libro «Gaussian Processes for Machine Learning» (GPML):
- Algoritmo 3.1: búsqueda de la moda de la distribución a posteriori y cálculo del NLML,
- Algoritmo 5.1: cálculo de los gradientes del NLML con respecto a los hiperparámetros.
Fig. 1. Algoritmo 3.1 para buscar la moda f_hat y calcular el LML
- El proceso comienza con el cálculo de la matriz de covarianza del núcleo K. Para lograr una convergencia más rápida y estable, se utiliza el resultado mu_f_train de la iteración anterior de optimización de hiperparámetros como aproximación inicial de f, en lugar de empezar cada vez con un vector nulo.
- Mediante el método de Newton, f_hat se actualiza de forma iterativa.
- Cálculo del NLML: una vez que converge la moda f (cuando la variación del LML es inferior a m_tolerance), pasamos el valor del NLML al optimizador.
Fig. 2. Algoritmo 5.1: cálculo de los gradientes del LML
- Los gradientes respecto a los hiperparámetros del núcleo incluyen el cálculo de las matrices auxiliares R y C. La matriz R se calcula utilizando cho_solve, mientras que la matriz C se calcula mediante LinearEquationsSolutionTriangular.
- Luego se calcula la matriz de covarianza a posteriori de la función latente sobre los datos de entrenamiento: Σf_train = K−C^TC.
- El gradiente para cada hiperparámetro del núcleo (dK_dtheta_j) se compone de una parte explícita (s1) y de partes implícitas (s2, s3). El cálculo de s1 se ha optimizado mediante el uso de DiagonalTrace, mientras que s3 también incluye cho_solve para resolver sistemas.
- Los gradientes respecto a los hiperparámetros de la verosimilitud se calculan utilizando las derivadas de la log-verosimilitud y su matriz hessiana con respecto a los hiperparámetros correspondientes. Para ello se utilizan los métodos LogLikelihoodGradientParam y LogLikelihoodHessianDerivative del objeto likelihood.
Ahora que ya disponemos de todos los componentes básicos (núcleos, verosimilitudes y métodos de inferencia), pasamos a mostrar cómo funciona la biblioteca.
Pruebas de la biblioteca con datos sintéticos
El script Gpsynthetic.mq5 nos ayudará a probar nuestra biblioteca con datos sintéticos sencillos. Esto permitirá verificar el funcionamiento y la corrección de los métodos implementados.
// --- Incluimos la biblioteca de GP --- #include <GP/GP.mqh> enum IntervalType { INTERVAL_F = 0, // Intervalo de confianza de la función latente f* INTERVAL_Y = 1 // Intervalo de confianza de las observaciones y* }; enum Type_inference { Exact = 0, // Inferencia exacta Laplace = 1 // Inferencia de Laplace }; enum Type_Data { Regression = 0, // Regresión Classification = 1 // Clasificación }; //--- Parámetros de entrada input IntervalType interval_type = INTERVAL_F; // Tipo de intervalo que se mostrará input Type_Data DataType = Classification; input Type_inference inf = Laplace ; // Tipo de inferencia //+------------------------------------------------------------------+ //| Función de inicio del script | //+------------------------------------------------------------------+ void OnStart() { string InpFileName1; string InpFileName2; string InpFileName3; string InpFileName4; if(DataType == Regression) { // Datos para la regresión 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"; } //----------- Conjunto de datos 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); // Creamos vectores de etiquetas a partir de matrices vector y_train_ = y_train.Col(0); vector y_test_ = y_test.Col(0); //--- 1. Creación de objetos de núcleos IKernel* rbf = new RBFKernel(1,1); // IKernel* linear = new LinearKernel(1.0); // IKernel* periodic = new PeriodicKernel(1.0, 1.0, 5.8); //--- 2. Combinación de núcleos (creación de un núcleo compuesto) // IKernel* functional_kernels_array[] = {rbf, linear, periodic}; // SumKernel* combined_functional_kernel = new SumKernel(functional_kernels_array); //--- 3. Creación de un objeto de verosimilitud ILikelihood* likelihood = NULL; // Declaramos un puntero de tipo base if(DataType == Regression) { likelihood = new GaussianLikelihood(1); // Asignamos un objeto } if(DataType == Classification) { likelihood = new LogitLikelihood(); } //--- 4. Creación de un objeto de inferencia IInference* inference; // Declaramos un puntero a la interfaz // Elegimos el tipo de inferencia: if(inf == Exact) { inference = new ExactInference(); } else { inference = new LaplaceInference(100, 1e-10); } //--- 5. Creación de una instancia del modelo de GP 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. Optimización de hiperparámetros gp_model.Fit(); ulong end_time_fit = GetMicrosecondCount(); double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; Print(inference.GetName()); Print("Tiempo de ejecución de Fit: ", StringFormat("%.3f", elapsed_time_ms), " ms"); gp_model.PrintOptimizedKernelParameters(); //--- 7. Predicción para los puntos de prueba GPPredictionResult gp_predictions; // Creamos una estructura en la que se guardarán los resultados de las predicciones 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("Tiempo de ejecución de Predict: ", StringFormat("%.3f", elapsed_time_ms), " ms"); // --- 8. Resultados de la clasificación if(DataType == Classification) { DisplayClassificationResults(gp_predictions, y_test_); } //--- 9. Resultados de la regresión if(DataType == Regression) { VisualizeGP(x_train, y_train, x_test, y_test_, gp_predictions, gp_model, 15); } //--- 10. Liberación de memoria delete gp_model; }
En función del DataType seleccionado (regresión o clasificación), el script determina las rutas de acceso a los archivos CSV con los datos de entrenamiento y de prueba. Se supone que estos archivos se encuentran en las carpetas «Data_Regression» o «Data_Classification» de la sección «Files».
A continuación, se crean los objetos de núcleo. Para construir una función de covarianza más compleja, se pueden utilizar núcleos compuestos, como SumKernel o ProductKernel.
Luego se crea un objeto de función de verosimilitud (ILikelihood) según el DataType seleccionado, así como un objeto de inferencia (IInference).
Mediante la macro CREATE_GP_MODEL (o un constructor similar) se crea un modelo de GP, que vincula el núcleo seleccionado, la función de verosimilitud y el método de inferencia con los datos de entrenamiento.
El método gp_model.Fit(int maxiter = 20) inicia el proceso de entrenamiento del modelo. En el método se ha añadido un nuevo parámetro, «maxiter», que permite establecer el número máximo de iteraciones.
Una vez finalizado el entrenamiento, se utiliza el método gp_model.Predict() para realizar predicciones sobre los nuevos datos (de prueba) x_test. Los resultados de las predicciones (valores medios, covarianzas, probabilidades/etiquetas) se guardan en la estructura GPPredictionResult.
Para la clasificación, se puede seleccionar el modo de predicción (PROBIT, NUM_INTEGR, MONTE_CARLO).
A continuación, el script muestra los resultados de la predicción:
- para la clasificación: se invoca la función DisplayClassificationResults(), que calcula la métrica de precisión y las probabilidades predichas;
- En el caso de la regresión, la función VisualizeGP() se encarga de representar gráficamente los resultados.
Para la regresión se utilizan datos sintéticos generados por la función y = sin(x) + 0.5 * x + noise(sigma = 0.1). No nos detendremos en detalle en esta tarea, ya que la tratamos en el artículo dedicado a la regresión. Simplemente observe cuánto ha mejorado el tiempo de entrenamiento gracias a que hemos pasado a utilizar gradientes analíticos y métodos de la biblioteca OpenBLAS.
Para la clasificación binaria se utilizan datos que representan un problema XOR (disyunción exclusiva). Se trata de un problema clásico del aprendizaje automático que pone de manifiesto las limitaciones de los modelos lineales y la necesidad de recurrir a enfoques no lineales, como los GP.
Recordemos que XOR es una operación lógica que toma dos entradas binarias (0 o 1) y devuelve 1 si exactamente una de las entradas es 1, y 0 en caso contrario. Tabla de verdad de XOR:
La tarea de clasificación binaria XOR consiste en construir un modelo que, a partir de dos entradas (A, B), prediga su salida XOR (0 o 1). Pero, dado que nuestro modelo solo funciona con las etiquetas «+1» y «-1», sustituiremos la etiqueta «0» por «-1».
Los datos para el problema XOR se generan de la siguiente manera. Se generan puntos de entrada bidimensionales X (X = [x₁, x₂]), distribuidos aleatoriamente en un intervalo determinado (por ejemplo, [−4, 4] en cada eje). Las etiquetas y de estos puntos se determinan a partir de la lógica XOR:
- Si x₁ y x₂ tienen el mismo signo (ambos positivos o ambos negativos), su producto x₁⋅x₂ será positivo. En este caso, se asigna la etiqueta +1.
- Si x1 y x2 tienen signos diferentes (uno positivo y otro negativo), su producto x1⋅x2 será negativo. En este caso, a la etiqueta se le asigna el valor −1.
Por lo tanto, los puntos (+,+) y (−,−) pertenecen a la clase +1, mientras que los puntos (+,−) y (−,+) pertenecen a la clase −1. Para esta tarea utilizaremos únicamente el núcleo RBF, ya que maneja muy bien las dependencias no lineales, lo que permite que los procesos gaussianos separen eficazmente dichas clases.
Durante las pruebas, utilicé 200 observaciones para el entrenamiento y 100 puntos nuevos para comprobar la exactitud. La precisión de la predicción fue del 98 %, lo que demuestra la alta eficacia de los GP en la resolución de problemas de clasificación no lineales.
Indicador GPRegressor: Regresión temporal mediante GP
El indicador es un ejemplo de previsión dinámica de series temporales financieras que tiene en cuenta la incertidumbre. Mediante el uso de una ventana deslizante para generar el conjunto de entrenamiento y el reentrenamiento constante del modelo, aplica el principio de adaptación a las condiciones cambiantes del mercado.
//+------------------------------------------------------------------+ //| GPRegressor.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Incluimos la biblioteca de 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 //--- Trazado del intervalo de confianza #property indicator_label2 "Confidence interval" #property indicator_type2 DRAW_FILLING #property indicator_color2 clrLightGray #property indicator_width2 1 //--- Parámetros de entrada input int WindowLength = 10; // Longitud de la ventana para entrenar el modelo input int calculate_bars = 100; // Profundidad de la historia para el cálculo input double zscore = 1.96; // Valor z para el intervalo (1,96 = 95 %) input bool ShowRMSE = false; // Mostrar la métrica RMSE //--- Variables de búfer para el dibujado double ExtPredictionBuffer[]; // Búfer para los valores predichos double ExtUpperBandBuffer[]; // Búfer para el límite superior del intervalo de confianza double ExtLowerBandBuffer[]; // Búfer para el límite inferior //--- Objetos globales GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Matriz global para los índices temporales X_train matrix g_X_train; //--- Variables globales para acumular datos para el RMSE vector gp_pred; vector naive_pred; vector y_true; // Bandera que indica si el RMSE ya se ha calculado (para calcularlo una sola vez) bool rmse_calculated = false; string comment; double gp_rmse,naive_rmse; //+------------------------------------------------------------------+ //| Función de inicialización del indicador personalizado | //+------------------------------------------------------------------+ int OnInit() { if(Bars(Symbol(), Period()) < WindowLength + calculate_bars) { PrintFormat("Error: barras insuficientes para el cálculo. Se necesitan un mínimo de %d. Disponibles: %d", WindowLength + calculate_bars, Bars(Symbol(), Period())); return(INIT_FAILED); } rmse_calculated = false; // Restablecemos la bandera al inicializar o reinicializar //--- Inicializamos los vectores globales para el RMSE gp_pred.Resize(calculate_bars-1); naive_pred.Resize(calculate_bars-1); y_true.Resize(calculate_bars-1); //--- asignación de los búferes del indicador SetIndexBuffer(0, ExtPredictionBuffer, INDICATOR_DATA); SetIndexBuffer(1, ExtUpperBandBuffer, INDICATOR_DATA); SetIndexBuffer(2, ExtLowerBandBuffer, INDICATOR_DATA); //--- Creamos los objetos del núcleo, la verosimilitud y el método de inferencia 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); } //--- Inicializamos la matriz de características X_train (una característica: el tiempo) 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] } //--- Creamos un objeto GaussianProcess // Pasamos g_X_train y un y_train vacío 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[]) { // Comprobación de que haya suficientes barras if(rates_total < WindowLength) { Print("Error: barras insuficientes para el cálculo. Se necesita un mínimo de ", WindowLength, ", disponibles: ", rates_total); return(0); } int start; if(prev_calculated == 0) { // Determinamos a partir de dónde empezar el cálculo. // calculate_bars es la profundidad del historial que se va a utilizar para el cálculo. // MathMax(WindowLength, ...) garantiza que no empecemos antes de disponer de una ventana de datos completa. start = MathMax(WindowLength, rates_total - calculate_bars); // Inicializamos todos los búferes con el valor EMPTY_VALUE ArrayInitialize(ExtPredictionBuffer, EMPTY_VALUE); ArrayInitialize(ExtUpperBandBuffer, EMPTY_VALUE); ArrayInitialize(ExtLowerBandBuffer, EMPTY_VALUE); // Establecemos el inicio del dibujo de los búferes 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]) { // Esto significa que NO ha aparecido una barra nueva. // Estamos en la misma barra; simplemente han llegado nuevos ticks. return(rates_total); } // Si la hora de la última barra ha cambiado, significa que ha aparecido una nueva barra cerrada. // Empezamos el cálculo a partir de prev_calculated para procesar todas las barras nuevas. start = prev_calculated; } for(int i = start; i < rates_total && !IsStopped(); i++) { // Comprobación de que haya suficientes barras para formar la ventana if(i < WindowLength) { Print("Error: barras insuficientes para formar la ventana en la barra ", i, ". Se necesitan: ", WindowLength, ", disponibles: ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Omitimos esta barra } // Formación del vector y a partir de close[i-1], close[i-2], ..., close[i-WindowLength] // close[i-1] - la barra más reciente de la ventana // close[i-WindowLength] - la barra más antigua de la ventana vector y(WindowLength); for(int j = 0; j < WindowLength; j++) { y[WindowLength - 1 - j] = close[i - 1 - j]; } //Print(" y = ", y[WindowLength-1]); // Barra más reciente de la ventana de datos //--- Estandarización de datos 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; //--- Establecemos los datos de entrenamiento g_gp_model.SetTrainingData(g_X_train, y_standardized); //--- Restablecimiento de los hiperparámetros iniciales antes de 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(); //--- Entrenamos el modelo if(!g_gp_model.Fit()) { Print("Error: no se han podido optimizar los hiperparámetros GP para la barra ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Saltamos la barra actual } //ulong end_time_fit = GetMicrosecondCount(); // double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; //Print("Tiempo de ejecución de Fit: ", StringFormat("%.3f", elapsed_time_ms), " 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) // --- Creamos un punto para la predicción (X_star) // Predecimos el valor para el siguiente paso después de WindowLength matrix X_star = matrix::Zeros(1, 1); X_star[0][0] = (double)(WindowLength + 1); // --- Realizamos la predicción GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Error: no se ha podido obtener el pronóstico para el modelo GP para la barra ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; } // --- double predicted_mean_st = gp_predictions.mu_f_star[0]; // --- Realizamos la transformación inversa del valor medio predicho // de vuelta a la escala de precios original double predicted_mean = predicted_mean_st * sigma_raw + mu_raw; double sigma_f_star = gp_predictions.Sigma_f_star[0,0]; // varianza (Σ) de la señal f* double predicted_variance_st = gp_predictions.Sigma_y_star[0,0]; // varianza (Σ) de la observación y* // --- Realizamos la transformación inversa de la varianza // La varianza se escala por el cuadrado de la desviación estándar (sigma_raw) de los datos, // que se utilizaron para la estandarización. double predicted_variance = predicted_variance_st * sigma_raw * sigma_raw; // Cálculo del intervalo de confianza predicho double std_predict = MathSqrt(predicted_variance); double upper_bound = predicted_mean + zscore * std_predict; double lower_bound = predicted_mean - zscore * std_predict; // --- Guardamos los resultados en los búferes del indicador para su trazado 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); // Salida de depuración Print("Pronóstico para la barra ", i, ": Mean=", DoubleToString(predicted_mean, Digits()), " Upper=", DoubleToString(upper_bound, Digits()), " Lower=", DoubleToString(lower_bound, Digits())); } // --- Cálculo del RMSE para evaluar la eficacia del modelo predictivo de GP en comparación con una predicción «ingenua» sencilla if(ShowRMSE && !rmse_calculated) { vector returns(calculate_bars-1); for(int j=0; j <calculate_bars-1;j++) { // valor verdadero y_true[j] = close[rates_total - calculate_bars + j]; //Print("y_true[j]", y_true[j] ); // Predicción del GP: gp_pred[j] = ExtPredictionBuffer[rates_total-calculate_bars + j]; //Print("gp_pred[j]", gp_pred[j] ); // Predicción «ingenua» (precio de mañana = precio de hoy): 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(); // desviación estándar de los incrementos de 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); //--- Devolvemos rates_total para la siguiente llamada return(rates_total); } //+------------------------------------------------------------------+ //| Función de desinicialización del indicador personalizado | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { //--- Borramos el comentario del gráfico Comment(""); //--- Liberamos la memoria asignada al objeto GP if(g_gp_model != NULL) { delete g_gp_model; } } //+------------------------------------------------------------------+
Parámetros de entrada:
- WindowLength: longitud de la ventana de datos (en barras) para entrenar el modelo GP.
- calculate_bars: profundidad del historial para el cálculo y el trazado del indicador.
- zscore: número de desviaciones estándar para calcular el intervalo de confianza (por ejemplo, 1,96 para el 95 %).
- ShowRMSE: bandera para mostrar la raíz del error cuadrático medio (RMSE) de la predicción.
Inicialización (OnInit):
- Objetos globales. Se crean e inicializan los objetos globales del núcleo (RBFKernel), de la función de verosimilitud (GaussianLikelihood) y del método de inferencia (ExactInference).
- Características X_train. GPRegressor realiza una regresión respecto al tiempo utilizando los índices de tiempo (de 1 a WindowLength) dentro de la ventana de entrenamiento como características de entrada X.
- El modelo GaussianProcess. Se inicializa el objeto g_gp_model con los componentes indicados; y_train está inicialmente vacío y se irá rellenando de forma dinámica.
Ciclo principal de cálculo (OnCalculate):
- Formación de datos: en cada barra se extraen los WindowLength últimos precios de cierre para formar el vector y.
- Estandarización de los datos: y se estandariza (se lleva a media cero y varianza unitaria) para la estabilidad numérica del GP.
- Actualización y entrenamiento del modelo:
- Los hiperparámetros se restablecen a sus valores iniciales antes de cada entrenamiento (g_gp_model.SetHyperparameters)
- El modelo GP se entrena en la ventana actual de precios estandarizados: g_gp_model.Fit().
- Predicción: para pronosticar la siguiente barra, se crea un punto X_star con el índice temporal WindowLength + 1. El método g_gp_model.Predict() realiza la predicción y devuelve la media (mu_f_star) y la varianza (Sigma_y_star) para los datos estandarizados.
- Desestandarización: la media y la varianza predichas se transforman de nuevo a la escala de precios original utilizando la media (mu_raw) y la desviación estándar (sigma_raw) de la ventana actual.
- Intervalo de confianza: los límites superior e inferior del intervalo se calculan a partir de la media y la desviación estándar transformadas, utilizando el valor de zscore (número de desviaciones estándar) especificado por el usuario.
Fig. 3. Predicción del modelo de regresión e intervalo de confianza del 95 %
Cálculo de la métrica RMSE:
El indicador muestra, una sola vez y a petición del usuario, las estadísticas de RMSE para las predicciones del modelo de GP y la predicción «ingenua» (en la que el precio de mañana = el precio de hoy) correspondientes a las últimas calculate_bars. Ayuda a evaluar de forma objetiva la eficacia actual del modelo.
Por ejemplo, los valores típicos para el EURUSD en un marco temporal de un minuto podrían ser: RMSE GP (0,00032) | RMSE Naive (0,00029), y la desviación estándar de los incrementos de los precios de cierre (StdDev Returns Close) sería de 0,00029.
El hecho de que el RMSE de GP sea superior al RMSE Naive indica que el desempeño del modelo de GP es inferior al de una simple predicción «ingenua» en el periodo analizado. Además, la igualdad entre RMSE Naive y StdDev Returns Close es un rasgo característico de un paseo aleatorio, en el que las variaciones futuras del precio no dependen de las pasadas y el valor actual es la mejor predicción.
De ello podemos concluir que, en este tramo de la serie temporal, o bien no existen patrones predecibles, o bien la configuración actual del modelo de GP (que utiliza únicamente el tiempo como característica) aún no es capaz de extraer información útil que le permita superar a una predicción ingenua.
Indicador GPClassifier: clasificación basada en procesos gaussianos
El indicador GPClassifier muestra el proceso de creación y entrenamiento de un modelo de GP, así como la predicción con este, en una tarea de clasificación binaria. Como características se utilizan los incrementos de los precios de cierre, que se calculan en una ventana deslizante de datos históricos. Mediante la función NormalizeData se lleva a cabo la normalización de la matriz de características para estabilizar el entrenamiento del modelo.
//+------------------------------------------------------------------+ //| GPClassifier.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Incluimos la biblioteca de GP --- #include <GP/GP.mqh> #property indicator_chart_window #property indicator_buffers 4 #property indicator_plots 2 // trazado para la clase predicha «+1» (flechas hacia arriba) #property indicator_label1 "Predicted Class +1" #property indicator_type1 DRAW_ARROW #property indicator_color1 clrGreen #property indicator_width1 1 // trazado para la clase predicha «-1» (flechas hacia abajo) #property indicator_label2 "Predicted Class -1" #property indicator_type2 DRAW_ARROW #property indicator_color2 clrBlack #property indicator_width2 1 //--- Parámetros de entrada input int WindowLength = 10; // Longitud de la ventana para el entrenamiento del modelo de GP input int NumLags = 1; // Número de características input int calculate_bars = 100; // Profundidad de la historia para el cálculo input double ProbabilityThreshold = 0.5; // Umbral de probabilidad input int ArrowShiftPx = 10; // Desplazamiento de las flechas en píxeles input bool ShowACCURACY = false; // Mostrar la precisión en el gráfico //--- Variables de búfer para el dibujado double ExtClassUpBuffer[]; // Búfer para mostrar las flechas UP double ExtClassDownBuffer[]; // Búfer para mostrar las flechas DOWN double ExtPredictedClassBuffer[]; // Búfer para almacenar las etiquetas de clase predichas double ExtPredictedProbabilityBuffer[]; // Búfer para almacenar las probabilidades predichas //--- Objetos globales del proceso gaussiano GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Variables globales para las métricas de clasificación bool accuracy_calculated = false; vector gp_pred,naive_pred,y_true,gp_accuracy,naive_accuracy; string comment; int currentsize; //+------------------------------------------------------------------+ //| Función de inicialización del indicador personalizado | //+------------------------------------------------------------------+ int OnInit() { accuracy_calculated = false; // Restablecemos la bandera al inicializar o reinicializar //--- Inicializamos los vectores globales para Accuracy gp_pred.Resize(0); naive_pred.Resize(0); y_true.Resize(0); gp_accuracy.Resize(1); naive_accuracy.Resize(1); currentsize = 0; //--- asignación de los búferes del indicador SetIndexBuffer(0, ExtClassUpBuffer, INDICATOR_DATA); // Para las flechas hacia arriba SetIndexBuffer(1, ExtClassDownBuffer, INDICATOR_DATA); // Para las flechas hacia abajo SetIndexBuffer(2, ExtPredictedClassBuffer, INDICATOR_CALCULATIONS); // para almacenar las etiquetas predichas SetIndexBuffer(3, ExtPredictedProbabilityBuffer, INDICATOR_CALCULATIONS); // para almacenar probabilidades // --- Configuramos los símbolos de las flechas PlotIndexSetInteger(0, PLOT_ARROW, 225); // Flecha hacia arriba para el búfer 0 PlotIndexSetInteger(1, PLOT_ARROW, 226); // Flecha hacia abajo para el búfer 1 // --- Configuramos el desplazamiento de las flechas en píxeles --- // Para las flechas hacia arriba (búfer 0), utilizamos un desplazamiento negativo para que queden por encima de High PlotIndexSetInteger(0, PLOT_ARROW_SHIFT, -ArrowShiftPx); // negativo = hacia arriba // Para las flechas hacia abajo (búfer 1), utilizamos un desplazamiento positivo para que queden por debajo de Low PlotIndexSetInteger(1, PLOT_ARROW_SHIFT, ArrowShiftPx); // positivo = hacia abajo // Añadimos una comprobación de la relación entre NumLags y WindowLength if(NumLags <= 0) { Print("Error: NumLags debe ser un número positivo"); return(INIT_FAILED); } if(NumLags >= WindowLength) { Print("Error: NumLags (", NumLags, ") no puede ser mayor o igual a WindowLength (", WindowLength, ")."); return(INIT_FAILED); } //--- Creamos los objetos del núcleo, de la verosimilitud y del método de inferencia 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); } // --- Creamos el modelo GaussianProcess // Pasamos matrices vacías para X_train e y_train. // Se redefinirán en 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); } //+------------------------------------------------------------------+ //| Función de iteración del indicador personalizado | //+------------------------------------------------------------------+ 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[]) { // Comprobación de que haya suficientes barras if(rates_total < WindowLength + NumLags + 1) { Print("Error: barras insuficientes para los cálculos. Se necesitan un mínimo de ", WindowLength + NumLags + 1, ", disponibles: ", rates_total); return(0); } int start; if(prev_calculated == 0) { start = MathMax(WindowLength + NumLags + 1, rates_total - calculate_bars); // Inicializamos todos los búferes a EMPTY_VALUE ArrayInitialize(ExtClassUpBuffer, EMPTY_VALUE); ArrayInitialize(ExtClassDownBuffer, EMPTY_VALUE); // Número de barras iniciales sin trazado ni valores en 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; } // --- Ciclo principal de cálculo de predicciones for(int i = start; i < rates_total && !IsStopped(); i++) { // Comprobación de que haya suficientes barras para formar la ventana de X e y if(i < WindowLength + NumLags + 1) { ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } // --- Paso 1: Formación de los datos de entrenamiento (X_train e y_train) matrix X_train = matrix::Zeros(WindowLength, NumLags); vector y_train = vector::Zeros(WindowLength); for(int j = 0; j < WindowLength; j++) { // Cálculo del índice de la barra en el array close para la observación actual en la ventana. int idx = i - WindowLength + j; // Formación de características (rezagos) para la observación actual (fila de X_train[j]) for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { // Ejemplo: // lag_idx = 0: (close[idx - 1] - close[idx - 2]) - incremento de la barra anterior // lag_idx = 1: (close[idx - 2] - close[idx - 3]) - incremento de la barra de hace dos barras X_train[j,lag_idx] = close[idx - 1 - lag_idx] - close[idx - 2 - lag_idx]; } // Y_train[j]: Etiqueta binaria para el incremento del precio en la barra idx. if(close[idx] - close[idx - 1] > 0) { y_train[j] = 1; // El precio ha subido } else { y_train[j] = -1; // El precio ha bajado o se ha mantenido sin cambios } } // Print(y_train); // --- Normalización de X_train matrix X_train_norm; vector out_mean; vector out_std; if(!NormalizeData(X_train,X_train_norm,out_mean,out_std)) { Print("Error: no se ha podido normalizar la matriz de características ", i); continue; } // --- Paso 2: Introducimos los datos en el modelo de GP g_gp_model.SetTrainingData(X_train_norm, y_train); // Restablecimiento de los hiperparámetros iniciales antes de cada entrenamiento con 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); // --- Paso 3: Entrenamos el modelo // ulong start_time_fit = GetMicrosecondCount(); if(!g_gp_model.Fit()) { Print("Error: no se han podido optimizar los parámetros GP para la barra ", 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("Tiempo de ejecución de Fit: ", StringFormat("%.3f", elapsed_time_ms), " ms"); // --- Paso 4: Creamos un punto para la predicción X_star // X_star - incrementos de precio de las NumLags barras cerradas anteriores, // que se utilizan para predecir el movimiento de la barra actual (barra 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 también debe normalizarse utilizando las mismas medias y desviaciones estándar, // que se usaron para las columnas correspondientes de 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]; } // --- Paso 5: Realizamos la predicción GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Error: no se ha podido obtener el pronóstico del modelo GP para la barra ", i); ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } double predicted_probability = gp_predictions.predicted_probabilities[0]; // Probabilidad para un único punto de prueba int predicted_class = (int)gp_predictions.predicted_labels[0]; // Etiqueta (+1 o -1) // --- Paso 6: Guardamos los resultados en los búferes del indicador para el trazado ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; ExtPredictedClassBuffer[i] = predicted_class; ExtPredictedProbabilityBuffer[i] = predicted_probability; // Añadimos al comentario del gráfico la probabilidad y el número de características (rezagos) comment = StringFormat("GP Classifier (Lags: %d) | Probability(UP) = %.10f", NumLags, predicted_probability); Comment(comment); // Si la clase prevista es «+1» y la probabilidad es superior al umbral establecido if(predicted_class == 1 && predicted_probability > ProbabilityThreshold) { ExtClassUpBuffer[i] = high[i]; } // Si la clase prevista es «-1» y la probabilidad de caída es superior al umbral establecido else if(predicted_class == -1 && (1.0 - predicted_probability) > ProbabilityThreshold) { ExtClassDownBuffer[i] = low[i]; } Print("Pronóstico para la barra ", i, ": Probability(UP)=" + DoubleToString(predicted_probability, 10) + ", Predicted Class=" + IntegerToString(predicted_class), " time: ", time[i]); } // --- Cálculo de la precisión de la clasificación para evaluar la eficacia del modelo predictivo de GP en comparación //con una predicción «ingenua». if(ShowACCURACY && !accuracy_calculated) { for(int j=0; j <calculate_bars-1;j++) { // índice hasta la barra actual 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); // Etiqueta verdadera: calculamos hasta la barra actual if(close[bar_idx] - close[bar_idx - 1] > 0) { y_true[currentsize-1] = 1; // El precio ha subido } else { y_true[currentsize-1] = -1; // El precio ha bajado o se ha mantenido sin cambios } // Print("y_true", y_true[currentsize-1] ); // Predicción del GP: calculamos hasta la barra actual gp_pred[currentsize-1] = ExtPredictedClassBuffer[bar_idx]; // Print(" gp_pred", gp_pred[currentsize-1]); // Predicción «ingenua»: la etiqueta de la barra «de mañana» es igual a la de «hoy» 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); //--- Devolvemos rates_total para la siguiente llamada return(rates_total); } //+------------------------------------------------------------------+ //| Función de desinicialización del indicador personalizado | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { // Limpiamos el comentario del gráfico al eliminar el indicador Comment(""); if(g_gp_model != NULL) { delete g_gp_model; } }
Parámetros de entrada:
- WindowLength: longitud de la ventana de datos para el entrenamiento del modelo. Número de observaciones en la muestra de entrenamiento.
- NumLags: número de rezagos (incrementos de precios). Por ejemplo, si NumLags = 1, la característica será solo el incremento del precio en la barra anterior. Si NumLags = 3, las características serán los incrementos del precio en las tres barras anteriores, y así sucesivamente.
- calculate_bars: determina cuántas de las últimas barras del gráfico se usarán para realizar el cálculo y el trazado del indicador.
- ProbabilityThreshold: umbral de probabilidad. El indicador muestra señales mediante flechas hacia arriba o hacia abajo cuando la probabilidad predicha de que el precio se mueva en la dirección correspondiente supera ese umbral.
- ArrowShiftPx: desplazamiento de las flechas en píxeles respecto a los extremos de la barra.
- ShowACCURACY: bandera para mostrar las métricas de precisión de la clasificación en el gráfico
OnInit() — inicialización del modelo GP:
- Creamos el núcleo RBF, la función de verosimilitud y el método de inferencia:
- g_kernel = new RBFKernel(1.0, 1.0)
- g_likelihood = new LogitLikelihood()
- g_inference = new LaplaceInference()
- Creamos un modelo GaussianProcess:
- g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference, matrix::Zeros(1,1), vector::Zeros(1));
- El objeto g_gp_model se crea con matrices iniciales vacías para los datos de entrenamiento X_train e y_train. Estos datos se completarán y actualizarán dinámicamente en cada barra dentro de la función OnCalculate.
OnCalculate(): ciclo principal de cálculo:
1. Comprobación de que haya suficientes barras: antes de iniciar los cálculos, nos aseguramos de que haya suficientes barras en el gráfico para formar la ventana de entrenamiento, teniendo en cuenta WindowLength, NumLags y una barra para la etiqueta prevista.
2. Formación de los datos de entrenamiento (X_train e y_train):
- Características X_train: en cada paso, se forma una matriz de características X_train dentro de una ventana deslizante. Cada fila corresponde a una observación, y las columnas representan los incrementos de los precios (rezagos) durante las NumLags barras anteriores.
- Las etiquetas y_train son etiquetas binarias que indican la dirección del movimiento del precio en la siguiente barra: «1» (el precio subió) o «-1» (el precio bajó o no cambió).
3. Normalización de X_train — la matriz de características se estandariza para garantizar la estabilidad numérica. Los valores medios (out_mean) y las desviaciones estándar (out_std) obtenidos en esta etapa se guardan para normalizar el punto de predicción.
4. Actualización y entrenamiento del modelo GP:
- g_gp_model.SetTrainingData(X_train_norm, y_train): los datos de entrenamiento se actualizan en cada barra.
- g_gp_model.SetHyperparameters(initial_params) — los hiperparámetros del núcleo se restablecen y se vuelven a optimizar para cada nueva barra
- g_gp_model.Fit() — iniciamos el proceso de entrenamiento en la ventana actual de datos estandarizados.
5. Creación y normalización de un nuevo punto para la predicción (X_star): se normaliza utilizando los vectores out_mean y out_std.
6. Realización de la predicción: el método g_gp_model.Predict(X_star, gp_predictions) devuelve la probabilidad predicha de movimiento al alza (predicted_probabilities[0]) y la etiqueta binaria (predicted_labels[0]).
7. Trazado y muestra: en función de predicted_class y ProbabilityThreshold, se dibuja en el gráfico una flecha hacia arriba o hacia abajo. La información sobre la probabilidad de la predicción y la clase predicha se muestra como comentario en el gráfico y se registra en el registro.
Fig. 4. Resultado de la predicción del indicador GP Classifier (núcleo RBF)
Cálculo de las métricas de precisión (ShowACCURACY)
Al activar el parámetro ShowACCURACY, el indicador realiza un único cálculo de la precisión de la clasificación para las últimas calculate_bars barras, comparando las predicciones del modelo de GP con una predicción «ingenua».
Las etiquetas verdaderas (y_true) se definen como la dirección real del movimiento del precio. La predicción del GP (gp_pred) corresponde a las etiquetas predichas por el modelo de GP. La predicción «ingenua» (naive_pred) supone que la dirección del siguiente movimiento del precio será la misma que la del movimiento anterior (por ejemplo, si el precio subió en la barra anterior, la predicción «ingenua» predice una subida).
Luego se compara la métrica CLASSIFICATION_ACCURACY de ambos modelos (gp_accuracy y naive_accuracy), y los resultados se muestran en un comentario del gráfico. Esto permite evaluar de forma objetiva la eficacia del modelo de GP frente a un punto de referencia sencillo.
Conclusión
Hoy hemos concluido el análisis del modelo de clasificación de procesos gaussianos. Tras un análisis detallado de los fundamentos teóricos, hemos logrado crear una versión básica de la biblioteca que permite resolver tanto problemas de regresión como de clasificación. Las pruebas realizadas con datos sintéticos han confirmado el correcto funcionamiento de los componentes implementados de la biblioteca, incluida la aplicación eficaz de gradientes analíticos y el uso de la biblioteca de alto rendimiento OpenBLAS para garantizar la estabilidad numérica y la velocidad.
Los indicadores GPRegressor y GPClassifier demostraron las posibilidades prácticas de aplicar procesos gaussianos para la predicción dinámica de series temporales financieras en tiempo real, y métricas como el RMSE (para la regresión) y la precisión (Accuracy, para la clasificación) permitieron evaluar de forma objetiva la eficacia del modelo.
No obstante, la implementación de la biblioteca dista mucho de ser perfecta. A pesar de todos los esfuerzos realizados, su velocidad de ejecución es inferior a la de soluciones reconocidas, como, por ejemplo, scikit-learn. En parte, esto se debe a que utilizamos el optimizador MinBleic, que no es precisamente el más rápido.
Por lo tanto, la mejora futura de la biblioteca podría estar relacionada con:
- el uso de un optimizador de gradiente de mayor rendimiento. En concreto, necesitaremos un algoritmo L-BFGS rápido que supere con creces al actual en eficiencia para problemas con un gran número de hiperparámetros;
- una implementación más robusta del método de Newton para la aproximación de Laplace. La incorporación de una búsqueda lineal permitirá encontrar de forma más eficaz y fiable la moda de la distribución a posteriori, mejorando la convergencia de la optimización final;
- la implementación de métodos dispersos. A diferencia del trabajo actual con matrices completas, el uso de representaciones dispersas permitirá a la biblioteca entrenarse con volúmenes de datos considerablemente mayores, lo que abrirá nuevas posibilidades para el escalado y la aplicación de los GP.
Programas utilizados en el artículo
| # | Nombre | Tipo | Descripción |
|---|---|---|---|
| 1 | GPClassifier.mq5 | Indicador | Clasifica la dirección del precio un paso adelante basándose en las probabilidades predichas (aprendizaje online) |
| 2 | GPRegressor.mq5 | Indicador | Predice el precio y el intervalo de confianza un paso hacia delante (aprendizaje online) |
| 3 | GPsynthetic.mq5 | Script | Comprobación del modelo de GP con datos sintéticos |
| 4 | GP.mqh | Biblioteca de clases | Clase GaussianProcess, que combina el núcleo, la función de verosimilitud y método de inferencia, y clase GPOptimizationObjective, responsable de la optimización |
| 5 | Inference.mqh | Biblioteca de clases | Interfaz IInference y sus implementaciones, que definen los métodos de inferencia |
| 6 | Kernels.mqh | Biblioteca de clases | Interfaz IKernel e implementaciones de núcleos de covarianza |
| 7 | Likelihoods.mqh | Biblioteca de clases | Interfaz ILikelihood e implementaciones de funciones de verosimilitud |
| 8 | StructUtils.mqh | Biblioteca de clases | Funciones auxiliares y estructuras de datos |
| 9 | DataRegression | CSV | Datos para el script GPsynthetic.mq5 |
| 10 | DataClassification | CSV | Datos para el script GPsynthetic.mq5 |
Traducción del ruso hecha por MetaQuotes Ltd.
Artículo original: https://www.mql5.com/ru/articles/19013
Archivos adjuntos |
MQL5.zip
(53.23 KB)
Advertencia: todos los derechos de estos materiales pertenecen a MetaQuotes Ltd. Queda totalmente prohibido el copiado total o parcial.
Este artículo ha sido escrito por un usuario del sitio web y refleja su punto de vista personal. MetaQuotes Ltd. no se responsabiliza de la exactitud de la información ofrecida, ni de las posibles consecuencias del uso de las soluciones, estrategias o recomendaciones descritas.
Predicción en el trading y modelos de Grey
En este artículo se analiza la aplicación de los modelos Grey para la predicción de series temporales financieras. Analizaremos los principios de funcionamiento de los modelos Grey y las particularidades de su aplicación a las series financieras. Discutiremos las ventajas y las limitaciones del uso de estos modelos en el trading.
Automatización de estrategias de trading en MQL5 (Parte 23): Recuperación por zonas con trailing stop y lógica de cestas
En este artículo, mejoramos nuestro sistema de recuperación por zonas mediante la incorporación de órdenes stop dinámicas y funciones de negociación con cestas múltiples. Analizamos cómo la arquitectura mejorada utiliza stops dinámicos para asegurar las ganancias y un sistema de gestión de cestas para gestionar de forma eficiente múltiples señales de negociación. Mediante la implementación y las pruebas retrospectivas, demostramos que se trata de un sistema de negociación más sólido, diseñado para adaptarse al comportamiento del mercado.
Gestor de riesgos para robots de trading (Parte I): archivo include para el control de riesgos en asesores expertos
El trading exige una estricta disciplina sn la gestión de riesgos. El presente trabajo presenta un análisis de las principales causas de los fracasos de los tráders y propone una solución técnica en forma de la clase CEnhancedRiskManager para la plataforma MQL5. Incluye pruebas prácticas con un asesor experto de cuadrícula agresivo.
Algoritmo del átomo artificial — Artificial Atom Algorithm (A3)
Implementación del algoritmo A3 en MQL5: un método metaheurístico de optimización inspirado en los procesos químicos. Con solo dos parámetros configurables, la compacidad y una población reducida proporcionan una alta velocidad de funcionamiento con una calidad suficiente de las soluciones.
Está perdiendo oportunidades comerciales:
- Aplicaciones de trading gratuitas
- 8 000+ señales para copiar
- Noticias económicas para analizar los mercados financieros
Registro
Entrada
Usted acepta la política del sitio web y las condiciones de uso
Si no tiene cuenta de usuario, regístrese
Buenos días,
¡Buenos días!
¿Es esto lo que quieres?
Ahora ya no se muestra la predicción para el siguiente paso, sino la media y la dispersión de la muestra de entrenamiento.
Muchas gracias, es muy amable por tu parte.
¿Se pueden entrenar los modelos de regresión de procesos gaussianos en modo en línea, de modo que se adapten en tiempo real a la media y la varianza de los datos?
Muchas gracias, es un detalle por su parte.
¿Es posible entrenar modelos de regresión basados en procesos gaussianos en modo en línea, de modo que puedan adaptarse a la media y a la dispersión de los datos en tiempo real?