Procesos gaussianos en el aprendizaje automático (Parte 1): Modelo de clasificación en MQL5

Introducción

Continuamos nuestro estudio del modelo de aprendizaje automático: los procesos gaussianos (GP). En el artículo anterior, examinamos con detalle un problema de regresión cuyo objetivo principal era predecir valores continuos. Hoy abordaremos un tema mucho más complejo: la clasificación. Su principal dificultad radica en que el proceso de inferencia para la clasificación en procesos gaussianos no tiene una solución analítica, lo cual requiere el uso de métodos aproximados como la aproximación de Laplace.

Para resolver eficazmente este complejo problema, desarrollaremos una biblioteca modular de procesos gaussianos en MQL5. Este enfoque nos permitirá estructurar el código separando el modelo GP en componentes independientes y ofrecerá una base sólida para futuras mejoras y extensiones. Esta biblioteca se convertirá en una herramienta universal tanto para tareas de regresión como de clasificación.

En esta primera parte del artículo, examinaremos con detalle la teoría de la clasificación GP, incluyendo las bases matemáticas de los métodos aproximados. También presentaremos la clase principal de la biblioteca, GaussianProcess, que unificará todos los componentes del modelo, y la clase GPOptimizationObjective, que se encarga de la comunicación con la biblioteca de optimización Alglib.

Clasificación

La clasificación es una tarea de aprendizaje automático que consiste en asignar un objeto a una de las categorías predefinidas. Por ejemplo, en finanzas, la clasificación puede ayudar a predecir si el precio de una acción subirá o bajará usando como base los datos históricos.

En este artículo, nos centraremos en la clasificación binaria, donde un objeto pertenece a una de dos clases, como "crecimiento" (+1) o "caída" (-1). A diferencia de métodos como las máquinas de vectores de soporte (SVM) o los árboles de decisión, que solo producen una etiqueta de clase, los procesos gaussianos permiten realizar una predicción probabilística. Por ejemplo, el modelo puede indicar que la probabilidad de que una acción suba es del 75%. Este tipo de información resulta especialmente valiosa en el trading, donde el grado de confianza en una predicción ayuda a tomar decisiones informadas, permitiendo filtrar las señales poco fiables. 

Desafortunadamente, resolver un problema de clasificación utilizando GP resulta mucho más complejo que la regresión. Esto está relacionado con el tipo de verosimilitud usada: 

• En regresión, se suele utilizar la función de verosimilitud gaussiana. La combinación de la GP (como distribución a priori de la función) y la verosimilitud gaussiana permite obtener la distribución posterior analíticamente, lo que simplifica todos los cálculos.
• Para una clasificación en la que los objetivos son etiquetas de clase discretas, la verosimilitud gaussiana no resulta adecuada. En cambio, se usa, por ejemplo, la verosimilitud logit. Esto da como resultado que la distribución posterior tampoco es gaussiana y no tiene una solución analítica.

En consecuencia, debemos recurrir a métodos complejos de inferencia aproximada. La idea básica de estos métodos consiste en aproximar la verdadera distribución posterior no gaussiana usando una distribución gaussiana centrada en su moda.  En el marco de este artículo nos centraremos en la aproximación cuadrática de Laplace (Laplace Approximation), ya que es uno de los enfoques más simples y eficientes para obtener una aproximación gaussiana de la distribución a posteriori. 

Para la clasificación binaria, la idea subyacente de la predicción GP es bastante simple.  Comenzaremos con una distribución previa de las funciones latentes f(x). Imagine que el proceso gaussiano genera no solo una función, sino un conjunto infinito de funciones posibles, cada una de las cuales es una dependencia "oculta" potencial en los datos. Luego, cada una de estas realizaciones potenciales de la función latente f(x) se "pasa" a través de la función logística (sigmoide). La función sigmoide transforma cualquier número real (valor de f(x)) en una probabilidad entre 0 y 1, que será nuestra probabilidad previa π(x) de pertenecer a la clase “+1”:

Debemos señalar que π es una función determinista de f, pero dado que f en sí misma es estocástica (aleatoria, una muestra del proceso gaussiano), entonces la función π también se vuelve estocástica. Este concepto se ilustra claramente en las fig. 1 y 2 para un espacio de entrada unidimensional X.

Figura 1. Realización de la función latente f(x).

La figura 1 muestra una posible realización de la función latente, mostrando el comportamiento típico de la función correspondiente a los hiperparámetros del kernel dados.

Figura 2. La misma función, pero transformada mediante la función sigmoide.

La fig. 2 muestra el resultado de aplicar la función logística (sigmoide) a la misma función f(x):

De este modo, obtenemos una distribución de probabilidad a priori de pertenencia a la clase π(x)=σ(f(x)), que en esta etapa aún no considera los datos de entrenamiento y. Sin las observaciones de y, esta distribución previa sigue siendo solo nuestra hipótesis inicial, sin respaldo de evidencia empírica; sin ella, el modelo carece de información sobre cuáles de sus supuestos iniciales eran correctos y cuáles requieren revisión. 

Naturalmente, la elección de supuestos a priori influye sustancialmente en los resultados finales a posteriori. Esta es una característica clave del enfoque bayesiano, ya que precisamente de la decisión del investigador sobre el tipo de núcleo dependen las propiedades de la distribución a priori de las funciones y, por lo tanto, el modelo final. 

Proceso de inferencia

Por ello, para hacer predicciones fundamentadas, necesitamos considerar los datos de entrenamiento reales y. Aquí es donde entra en juego el proceso de inferencia. Su objetivo principal es transformar nuestras creencias a priori en creencias a posteriori, es decir, creencias ajustadas para considerar los datos observados. En clasificación, este proceso se divide naturalmente en dos pasos secuenciales.

Paso 1: Distribución predictiva de la función latente f* 

En el primer paso, calculamos p(f*|X, y, x*), la distribución posterior de la función latente f* para un nuevo punto de prueba x*, dados los datos de entrenamiento observados (X, y). Esta se define mediante la siguiente integral:

dónde,

• p(f*∣X, x*, f) es la distribución condicional de la función latente f* en un nuevo punto de prueba x*, dadas las funciones latentes f en los puntos de entrenamiento X. Esta distribución es siempre normal, ya que el GP por definición tiene una distribución normal conjunta,
• p(f|X, y) es la distribución posterior de las funciones latentes f en los datos de entrenamiento. Debido a la función de verosimilitud no lineal (sigmoide), no es gaussiana.

Debemos señalar que, dado que p(f|X, y) no es normal, esta integral no tiene solución analítica. Y esto significa que necesitaremos métodos aproximados para calcularlo. 

Paso 2: Probabilidad predictiva final π*

En el segundo paso, utilizamos esta distribución predictiva para formar una predicción probabilística final π*, la probabilidad de que un punto de prueba x* pertenezca a la clase positiva (y* = +1):

Aquí σ(f*) es la función logística (sigmoide), que transforma el valor de la función latente f* en una probabilidad entre 0 y 1. La integral en sí significa que promediamos estas probabilidades sobre todos los valores posibles de f* ponderados por su distribución predictiva posterior. En esencia, esta integral unidimensional supone la esperanza matemática de la función σ(f*) con respecto a la distribución p(f*|X, y, x*).

Nuevamente, para la verosimilitud logit, esta integral no tiene solución analítica. Por ello, también aquí necesitaremos métodos aproximados. Más adelante veremos que nuestra biblioteca GP implementa tres de estas aproximaciones, lo cual le permite elegir el método apropiado según los requisitos de precisión y los costos computacionales:

•  aproximación probit,
•  integración numérica,
•  método de Monte Carlo.

Estos dos pasos que acabamos de describir —el cálculo de la distribución posterior de la función latente y luego la integración para obtener la probabilidad predictiva— representan el marco general para la inferencia bayesiana en GP. Estas son las dos integrales que debemos calcular para obtener la predicción deseada, y ambas requieren el uso de métodos aproximados. 

La aproximación de Laplace

Como ya hemos comprobado, la inferencia bayesiana para la clasificación produce integrales que no se pueden resolver de forma cerrada. La aproximación de Laplace resuelve este problema aproximando la distribución no gaussiana p(f∣X, y) con la distribución gaussiana q(f∣X, y). Como la distribución condicional p(f*∣X, x*, f) también es gaussiana, la distribución predictiva resultante p(f*∣X, y, x*) también se vuelve gaussiana. Esto nos permite derivar fórmulas analíticas para la media y la varianza de f*, lo cual simplifica significativamente los cálculos posteriores. De este modo, la belleza y la eficiencia computacional de la aproximación de Laplace residen en su capacidad para reducir el cálculo de la distribución posterior y las predicciones a operaciones con distribuciones gaussianas.

Debemos comprender que la aproximación de Laplace es un compromiso. Esta hace que un problema intratable en forma analítica sea computacionalmente tratable, pero a costa de la precisión en la representación de la forma real de la distribución a posteriori. La calidad de esta aproximación normal depende directamente de qué tan cerca esté la distribución real p(f∣X, y) de ser normal. Cuanto más cerca esté, más precisa será la aproximación, y viceversa.

Si nos interesa la distribución real de p(f*∣X, y, x*), y no su aproximación, entonces los métodos MCMC (Markov Chain Monte Carlo) se suelen usar para ello. Si bien el método MCMC puede ofrecer estimaciones más precisas, resulta computacionalmente muy costoso y difícil de implementar. El MCMC puede utilizarse como método de referencia para compararlo con métodos de inferencia aproximada.

Ahora veamos más de cerca qué es la aproximación de Laplace. Dicha aproximación se construye en torno a la moda (máximo) de la verdadera distribución posterior p(f∣X, y). Usa un desarrollo de Taylor de segundo orden del logaritmo de la densidad posterior alrededor de este modo. Matemáticamente, aproximamos el logaritmo de la densidad posterior de la siguiente manera:

dónde,

• q(f∣X, y) — aproximación gaussiana para la distribución posterior p(f∣X, y),
• f_hat = argmax(f) p(f|X, y) — moda de la distribución posterior,
• A = −∇∇ log p(f|X, y)|f=f_hat — hessiano del logaritmo negativo de la distribución posterior en el punto de moda.

En primer lugar, para realizar la aproximación de Laplace, necesitamos encontrar el valor más probable de la función latente f, es decir, la moda f_hat. Para obtener la probabilidad posterior p(f∣X, y), usamos la regla de Bayes. Ya sabemos que esta regla relaciona la distribución posterior con la verosimilitud p(y∣f), la distribución a priori p(f∣X) y la verosimilitud marginal p(y∣X) de la siguiente manera:

Para maximizar p(f∣X, y) con respecto a f, no necesitamos conocer la constante de normalización p(y∣X), ya que no depende de f y, por consiguiente, no afecta la posición del máximo. Por ello, podemos trabajar con una distribución posterior no normalizada, que es proporcional al producto de la verosimilitud y la distribución previa p(y∣f)p(f∣X).

Para simplificar los cálculos y evitar problemas numéricos con valores de verosimilitud muy pequeños, tomamos el logaritmo de esta distribución posterior no normalizada. Debido a la propiedad de los logaritmos, el producto de verosimilitud se convierte en la suma de sus logaritmos:

Ψ(f) es la función objetivo que maximizaremos usando el método de Newton para encontrar la moda de la función latente. El método de Newton requiere calcular las primeras y segundas derivadas de Ψ(f) con respecto a f.

Al derivar esta ecuación con respecto a f, obtenemos:

dónde,

• W = −∇∇ log p(y|f) — hessiano negativo de la log-verosimilitud, que es una matriz diagonal.

Una vez que hemos calculado el gradiente y el hessiano, encontramos la moda de manera iterativa utilizando el método de Newton:

En cada iteración, el método de Newton actualiza nuestra estimación del modo actual en la dirección determinada por el gradiente y el hessiano hasta que se alcanza la convergencia. 

Una vez que tenemos la moda, podemos calcular la matriz de covarianza de la distribución gaussiana aproximada.  Esta matriz es igual al hessiano inverso negativo de Ψ(f) calculado en el punto de moda f_hat.

Por ello, la matriz de covarianza Σ de nuestra aproximación gaussiana es:

Con esto se completa la descripción del primer paso de la aproximación de Laplace: encontrar una aproximación normal de la distribución posterior.

Predicción en la aproximación de Laplace

Una vez hemos obtenido q(f∣X, y), podemos proceder al segundo paso de la inferencia: la predicción para nuevos puntos de prueba x∗. En esta etapa queremos encontrar la distribución predictiva p(f∗∣X, y, x∗). Debido a la aproximación de Laplace que hizo que p(f∣X, y) fuera gaussiana (en la forma q(f∣X, y)), y al hecho de que p(f∗∣X, x∗, f) también es gaussiana, la distribución predictiva resultante p(f∗∣X, y, x∗) también se volverá gaussiana, lo que nos permite obtener analíticamente su media y varianza posteriores.

El valor medio de la función latente f* para un nuevo punto de prueba x* (mu_f_star) se calcula como:

La varianza de la función latente Var(f*) para un nuevo punto de prueba x* (Sigma_f_star) se calcula como:

Ahora que tenemos la media y la varianza de la distribución predictiva, finalmente podemos calcular la probabilidad deseada de pertenecer a la clase π*: 

Esta fórmula es la base de la predicción de verosimilitud en el clasificador GP basado en Laplace.

Posiblemente haya notado que no hemos calculado las probabilidades de clase simplemente como σ(E[f*]), es decir, sustituyendo la media posterior de f* directamente en la función sigmoide. Este enfoque, conocido como predicción MAP (Maximum A Posteriori Prediction), sin duda puede ser útil.

Sin embargo, al calcular la predicción MAP σ(E[f*]), ignoramos la incertidumbre en f*. Simplemente tomamos el punto central f* (la media) y lo convertimos en una probabilidad. Al calcular E[σ(f*)] (que corresponde a la integración), tenemos en cuenta toda la forma de la distribución de f*. Esto nos ofrece una probabilidad predictiva más precisa y significativa, sobre todo cuando existe una incertidumbre significativa en f* (es decir, una gran varianza de V[f*]), o cuando la distribución de f* es asimétrica. Este enfoque se denomina probabilidad predictiva promedio (averaged predictive probability).

Comprender esta diferencia tiene importantes implicaciones prácticas:

• Si nuestro único objetivo es obtener una etiqueta de clase binaria (por ejemplo, "comprar" o "vender", +1 o -1), entonces utilizar la predicción MAP más simple puede ser suficiente, ya que producirá la misma etiqueta que la predicción promedio, que es computacionalmente más costosa.
• Sin embargo, si nos importan las probabilidades en sí, entonces las predicciones promedio (E[σ(f*)]) siguen siendo más precisas, ya que tienen plenamente en cuenta la incertidumbre del modelo. 

En el trading, una simple etiqueta binaria ("comprar" o "vender") no resulta suficiente. Necesitamos la fina gradación de certeza que proporcionan las probabilidades. El valor de probabilidad nos permite filtrar las señales de trading. Una señal con una probabilidad de éxito de 0,51 (que es solo ligeramente mejor que adivinar al azar) tendrá mucho menos valor que una señal con una probabilidad de 0,60. Esto permite al tráder establecer umbrales para realizar una operación. Por ejemplo, puede decidir que las operaciones solo se abrirán cuando la probabilidad de éxito sea superior a 0,55 o 0,60, reduciendo así el número de señales falsas.

Verosimilitud marginal en la aproximación de Laplace

Ahora que comprendemos el mecanismo de inferencia en GP para la clasificación, surge la pregunta: ¿cómo ajustar nuestro modelo para obtener predicciones óptimas? La respuesta reside en la verosimilitud marginal (Marginal Likelihood, LML). Esta es la función objetivo con la que optimizamos los hiperparámetros θ de nuestro modelo. Sin calcularla, es imposible encontrar los mejores parámetros que expliquen nuestros datos: 

donde, B

Una vez definida la función objetivo a optimizar, el siguiente paso importante es calcular sus derivadas parciales con respecto a los hiperparámetros θ. Esto es necesario porque realizaremos una optimización usando gradientes analíticos. Este método acelera los cálculos varias veces en comparación con los métodos numéricos. Los gradientes analíticos permiten que el optimizador se mueva de forma más eficiente y precisa hacia el mínimo de la función objetivo. 

El gradiente LML consta de una parte explícita y otra implícita:

Fórmula para calcular la parte explícita:

Aquí el problema principal consiste en calcular la derivada de la matriz del núcleo K con respecto a cada hiperparámetro. En la segunda parte del artículo, hablaremos de la implementación de derivadas para la función núcleo seleccionada.

La parte implícita consta de dos multiplicadores. El primer multiplicador de la parte implícita se obtiene mediante la fórmula:

Para calcular esta fórmula, debemos calcular la tercera derivada del logaritmo de la verosimilitud.

El segundo multiplicador de la parte implícita se calcula de la siguiente manera: 

En conclusión, observamos que NLML es necesario no solo para estimar hiperparámetros, sino también para comparar diferentes modelos (por ejemplo, con distintos tipos de núcleos). Los modelos con un valor NLML más bajo se consideran mejores porque implican una mayor probabilidad marginal, lo cual significa que el modelo explica mejor los datos observados.

Además, los GP, al usar la verosimilitud marginal para optimizar los hiperparámetros, abordan automáticamente el equilibrio entre el ajuste de los datos y la complejidad del modelo. NLML penaliza de forma natural los modelos excesivamente complejos, evitando así el sobreajuste. Gracias a esto, no será necesario establecer criterios de parada explícitos para evitar el sobreajuste, como se hace, por ejemplo, al entrenar redes neuronales. La optimización de NLML en sí misma se esfuerza por encontrar el equilibrio óptimo. Esta es una de las principales ventajas del enfoque bayesiano para los procesos gaussianos.

Biblioteca de procesos gaussianos

Ahora que hemos abarcado todos los conceptos teóricos necesarios, pasemos a la implementación práctica. Nuestro principal objetivo es crear una biblioteca de GP universal en MQL5 que sirva como una herramienta fiable para las tareas de previsión. Esta biblioteca tendrá una arquitectura modular, donde el modelo de GP se divide en componentes independientes e intercambiables, lo cual permitirá ampliar fácilmente sus capacidades y garantizará un mantenimiento sencillo. Se desarrollará considerando las siguientes características funcionales clave:

• La flexibilidad en la selección del núcleo: la capacidad de conectar fácilmente los núcleos de covarianza existentes, así como crear sus combinaciones (SumKernel, ProductKernel) para modelar dependencias más complejas;
• La compatibilidad con diversas funciones de verosimilitud;
• El soporte para diversos métodos de inferencia de distribución posterior para problemas de clasificación y regresión;
• La multifuncionalidad: la biblioteca debe ser universal, lo que permitirá resolver problemas tanto de regresión como de clasificación binaria;
• La optimización de hiperparámetros: el uso de gradientes analíticos para mejorar la velocidad y la precisión del proceso de entrenamiento. La integración con la biblioteca Alglib debería garantizar una optimización eficiente de los hiperparámetros del modelo.

Vamos a analizar más de cerca la estructura de la biblioteca. Consta de seis componentes principales, cada uno de los cuales implementa una funcionalidad específica:

• La clase GaussianProcess es el nodo central de la biblioteca, ya que gestiona todo el ciclo de vida de un modelo GP, desde la inicialización y la optimización de hiperparámetros hasta la realización de predicciones sobre nuevos datos.

• La clase GPOptimizationObjective es una clase auxiliar que sirve de "puente" entre nuestra biblioteca y la biblioteca de optimización Alglib. Adapta la función objetivo y su gradiente al formato requerido por Alglib (usando herencia de CNDimensional_Grad).

• La interfaz IKernel define un conjunto de métodos para diversas funciones de covarianza (núcleos). Incluye implementaciones como RBFKernel, LinearKernel, PeriodicKernel y sus combinaciones (SumKernel, ProductKernel).

• La interfaz ILikelihood define un conjunto de métodos para funciones de verosimilitud. Las implementaciones incluyen GaussianLikelihood para regresión y LogitLikelihood para clasificación binaria.

• La interfaz de inferencia proporciona los métodos para inferir la distribución posterior de la función GP latente. Actualmente, están implementados ExactInference y LaplaceInference.

• Las estructuras y utilidades auxiliares (StructUtils.mqh) suponen un conjunto de enumeraciones comunes, estructuras de datos (para resultados de inferencia y predicción) y funciones necesarias para trabajar con datos, matrices y gráficos para visualizar los resultados.

Gracias a esta estructura modular e interfaces bien definidas, podemos añadir fácilmente nuevos núcleos, métodos de inferencia y funciones de verosimilitud, lo cual facilita el desarrollo futuro de la biblioteca. 

La clase GaussianProcess

La clase GaussianProcess es la clase central de la biblioteca. Esta engloba toda la lógica necesaria para construir, entrenar y predecir un modelo de procesos gaussianos. Se ha diseñado según el principio de composición; GaussianProcess no contiene directamente funcionalidades de kernel, verosimilitud o inferencia. En cambio, integra estos componentes usando tres interfaces principales:

• el núcleo (IKernel),
• la función de verosimilitud (ILikelihood),
• el método de inferencia (IInference).

Esto permite adaptar de forma flexible el modelo GP a diversas tareas de predicción sin modificar la clase GaussianProcess subyacente.

//+------------------------------------------------------------------+
//| Gaussian process class                                           |
//+------------------------------------------------------------------+
class GaussianProcess
{
private:
    IKernel*      m_kernel;       // pointer to the selected kernel
    ILikelihood*  m_likelihood;   // pointer to the selected likelihood function
    IInference*   m_inference;    // pointer to the selected inference method
    matrix        m_X_train;      // Training input data Nxd
    vector        m_y_train;      // Training target data Nx1

    GPInferenceResult m_last_inference_result; // Structure storing the latest inference results
    
    int m_last_termination_type;  // Optimization operation completion code
    int m_last_iterations_count;  // Number of iterations performed by the optimizer
    double m_last_nlml_value;     // Final NLML value after optimization
    
private:
    // Auxiliary function for numerical integration
    double CalculateNumericalProbability(double mu_f_star, double sigma_f_star_diag, LogitLikelihood *likelihood);

public:
    // Class constructor
    GaussianProcess(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
    const matrix &X_train, const vector &y_train);
    // Static method for creating a GaussianProcess object with input parameters validation
    static GaussianProcess* Create(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
    const matrix &X_train, const vector &y_train);
    
    // Destructor
    ~GaussianProcess();
    
    // --- Methods for getting the model state ---
    // Return the results of the last inference operation
    GPInferenceResult GetLastInferenceResult() const;
    // Return the completion type of the last hyperparameter optimization
    int GetLastTerminationType() const;
    // Return the number of iterations performed during the last hyperparameter optimization
    int GetLastIterationsCount() const;
    // Return the negative logarithm of the marginal likelihood after optimization
    double GetLastNLML() const;
    // Return the pointer to the kernel in use
    IKernel* GetKernel() const;
    // Return the current values of all hyperparameters being optimized.
    vector GetCurrentHyperparameters();
  
    // --- Training and configuration methods ---
    // Run the full model training process, including hyperparameter optimization
    bool Fit();
    // Perform a single inference step without hyperparameter optimization
    bool PerformInference();
    // Set the training data for the model 
    void SetTrainingData(const matrix& X, const vector& y);
    // Set the given hyperparameters for the kernel and likelihood function
    void SetHyperparameters(const vector &params);
    // Method called by the optimizer to calculate the objective function (NLML)
    double CalculateNLMLObjective(const vector &hyperparameters);
    
    // --- The method performs a prediction for new test data
    // The predictmode parameter determines the method for calculating probabilities for classification (PROBIT, NUM_INTEGR, MONTE_CARLO)
    bool Predict(const matrix &X_test, GPPredictionResult &result, PredictMode mode = PROBIT);
    
    // --- Auxiliary methods ---
    // Static method for generating samples from prior GP
    static bool SamplePriorGP(const matrix &x, IKernel* kernel, int num_samples, matrix &f_samples,
                              bool plot_samples = false, int plot_display_seconds = 10);
    //--- Method for logging the final values of hyperparameters
    void PrintOptimizedKernelParameters();  
}; 

Vamos a analizar los métodos principales de la clase:

Existen dos formas principales de crear una instancia de una clase:

• El método Create: este método se usa para crear de forma segura un objeto GaussianProcess. Este método realiza las comprobaciones necesarias sobre los datos de entrada (X_train, y_train, punteros a las interfaces) y retorna un puntero al objeto o NULL en caso de error.

//+------------------------------------------------------------------+
//| Create method                                                    |
//+------------------------------------------------------------------+
GaussianProcess* GaussianProcess::Create(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
const matrix &X_train, const vector &y_train)
{
    // 1. Check for NULL pointers
    if (kernel == NULL || likelihood == NULL || inference == NULL) {
        Print("ERROR: Kernel, Likelihood, or Inference pointer is NULL");
        return NULL;
    } 
    // 2. Check the validity of X_train and y_train inputs
    if (X_train.Rows() == 0 || y_train.Size() == 0 || X_train.Rows() != y_train.Size()) {
    Print("ERROR: Invalid training data dimensions");
    return NULL;
    }
    // 3. Check the compatibility of 'likelihood' and 'inference'
    string likelihood_name = likelihood.GetName();
    string inference_name  = inference.GetName();   
    if (inference_name == "ExactInference" && likelihood_name != "GaussianLikelihood") {
        Print("ERROR: ExactInference supports only GaussianLikelihood!");
        delete kernel; delete likelihood; delete inference;
        return NULL;
    }    
    // 4. If all checks are passed, create the object
    GaussianProcess* gp_model = new GaussianProcess(kernel, likelihood, inference,X_train, y_train);
    if (gp_model == NULL) {
        Print("ERROR: Failed to create GaussianProcess object");
        delete kernel; delete likelihood; delete inference;
        return NULL;
    }    
    return gp_model;
} 

• Un constructor de clase: este posibilita una forma directa de inicialización, sin ninguna comprobación de datos. Si nuestros datos son fiables, podemos crear un objeto utilizando un constructor.

//+------------------------------------------------------------------+
//| GaussianProcess class constructor                                |
//+------------------------------------------------------------------+
GaussianProcess::GaussianProcess(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
const matrix &X_train, const vector &y_train) :
    m_kernel(kernel),
    m_likelihood(likelihood),
    m_inference(inference),
    m_X_train(X_train), 
    m_y_train(y_train), 
    m_last_termination_type(0),
    m_last_iterations_count(0),
    m_last_nlml_value(0.0){ } 

• El método Fit(): inicia el proceso completo de entrenamiento del modelo. Este método optimiza los hiperparámetros del kernel y la función de verosimilitud usando el optimizador MinBleic, que minimiza la verosimilitud marginal logarítmica negativa (NLML).

//+------------------------------------------------------------------+
//| Method for training the model                                    |
//+------------------------------------------------------------------+
bool GaussianProcess::Fit()
{
    // Create the GPOptimizationObjective object passing it the pointer to the current GaussianProcess object
    // This pointer goes into the private field of the m_gp class, with which we call the method
    // CalculateNLMLObjective to get the NLML value for the current set of hyperparameters
    GPOptimizationObjective objective_func(GetPointer(this));
    CNDimensional_Rep frep; 
    CObject Obj;
    vector initial_hyperparams = GetCurrentHyperparameters(); // Get the initial values of the hyperparameters
    double theta[];
    ArrayResize(theta, (int)initial_hyperparams.Size()); 
    VectorToArray(initial_hyperparams,theta);
    int num_params = (int)initial_hyperparams.Size();
    double s[];
    double bndl[];
    double bndu[];
    ArrayResize(s, num_params);
    ArrayResize(bndl, num_params);
    ArrayResize(bndu, num_params);

    int param_idx = 0; 
    IKernel* kernels_to_process[]; // array of pointers to the IKernel interface
    
    // Logic for obtaining kernels to set boundaries 
    // This block of code determines what type of kernel we are dealing with
    // and fills the kernels_to_process array with the corresponding pointers:
    if (dynamic_cast<SumKernel*>(m_kernel) != NULL) {          // Check if the current kernel m_kernel is a SumKernel object  
        SumKernel* sum_k = dynamic_cast<SumKernel*>(m_kernel); // If yes, then we cast the m_kernel type to the SumKernel* type
        sum_k.GetKernels(kernels_to_process); // and call the GetKernels() method, which fills the kernels_to_process array with all the kernels included in the sum
    } else if (dynamic_cast<ProductKernel*>(m_kernel) != NULL) { // Similar logic if the kernel is a ProductKernel object
        ProductKernel* prod_k = dynamic_cast<ProductKernel*>(m_kernel);
        prod_k.GetKernels(kernels_to_process);
    } else {
        ArrayResize(kernels_to_process,1); // If the kernel is neither a sum nor a product (i.e. it is not a composite kernel), 
        kernels_to_process[0] = m_kernel; // then the kernels_to_process array simply contains a pointer to m_kernel.
    }

   // This loop iterates over each base kernel found in the kernels_to_process array
   // and sets its hyperparameters to an initial scale s, a lower bound bndl, and an upper bound bndl
    for(int i = 0; i < ArraySize(kernels_to_process); i++) {
        IKernel* current_k = kernels_to_process[i];
            string kernel_name = current_k.GetName();   
            if (kernel_name == "RBFKernel") {
                if (param_idx + 2 <= num_params) {    
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                } 
            } else if (kernel_name == "LinearKernel") {
                if (param_idx + 1 <= num_params) {
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                } 
            } else if (kernel_name == "PeriodicKernel") {
                if (param_idx + 3 <= num_params) {
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                    s[param_idx] = 1.0; bndl[param_idx] = 1e-3; bndu[param_idx] = 1e3; param_idx++;    
                } 
            }           
    }
    
// --- Add bounds and scales for likelihood parameters (if any) ---
// LogitLikelihood has no hyperparameters, so this block will be skipped for it
// GaussianLikelihood has 1 parameter (sigma)
if (m_likelihood.GetNumHyperparameters() > 0) {
    if (param_idx + m_likelihood.GetNumHyperparameters() <= num_params) {
        s[param_idx] = 1.0;           // Scale
        bndl[param_idx] = 1e-10;      // Lower bound 
        bndu[param_idx] = 1e3;        // Upper bound
        param_idx++;
    } 
}
    CMinBLEICStateShell state;
    CMinBLEICReportShell rep; // object that will contain a report on the optimization results
    //-----------------------  optimizer stopping criteria
    double epsg = 0.0001;     //Gradient precision (0 means gradient stopping is disabled)
    double epsf = 0.0000;     //Precision by function value    
    double epsw = 0.0000;     //Accuracy by parameters 
    //-------------------------   
    double epso = 0.00001;    //Parameters for external convergence conditions in BLEIC
    double epsi = 0.00001;    //Parameters for internal convergence conditions in BLEIC
    CAlglib::MinBLEICCreate(theta, state);       // initialize the optimizer. It creates the initial state for MinBLEIC using the initial hyperparameter values from the theta array.
    CAlglib::MinBLEICSetBC(state, bndl, bndu);   // Set the lower (bndl) and upper (bndu) bounds for each parameter
    CAlglib::MinBLEICSetScale(state, s);         //Sets the scales (s) for each parameter. This can help the optimizer work more efficiently with parameters of different orders of magnitude.
    CAlglib::MinBLEICSetInnerCond(state,epsg,epsf,epsw);    
    CAlglib::MinBLEICSetOuterCond(state, epso, epsi);    
    CAlglib::MinBLEICOptimize(state, objective_func, frep, 0, Obj); // start the optimization    
    CAlglib::MinBLEICResults(state, theta, rep); // optimization report
    
    m_last_termination_type = rep.GetTerminationType();
    m_last_iterations_count = rep.GetInnerIterationsCount();
    m_last_nlml_value = objective_func.GetNLML(); // Get the final NLML
//------------------------------------------------------------------------------------    
//    TerminationType field contains completion code, which can be:
//-8     internal integrity control detected    infinite    or    NAN    values    in
//     function/gradient. Abnormal termination signalled.
//-3     inconsistent constraints. Feasible point is
//     either nonexistent or too hard to find. Try to
//     restart optimizer with better initial approximation
// 1     relative function improvement is no more than EpsF.
// 2     relative step is no more than EpsX.
// 4     gradient norm is no more than EpsG
// 5     MaxIts steps was taken
// 7     stopping conditions are too stringent,
//     further improvement is impossible,
//     X contains best point found so far.
// 8     terminated by user who called minbleicrequesttermination(). X contains
//     point which was "current accepted" when    termination    request    was
//     submitted.
//-------------------------------------------------------------------------------------  
// Determine the success of optimization based on TerminationType
    bool success = true;     
    if (m_last_termination_type < 0)
    {
        Print("Error: GP optimization failed. Completion type: ", m_last_termination_type);
        success = false;
    } 
    // Update the model hyperparameters after optimization
    vector optimized_hyperparams;
    optimized_hyperparams.Assign(theta);
    SetHyperparameters(optimized_hyperparams);   
    return success;     
} 

Dentro del método Fit(), preparamos todo lo necesario para que el optimizador funcione de manera eficaz.

Luego se crea un objeto especial, objective_func (GPOptimizationObjective), que representa la función objetivo NLML y su gradiente analítico en un formato comprensible para Alglib. Después se transmite a su constructor un puntero al objeto GaussianProcess actual (esto es necesario para llamar al método CalculateNLMLObjective).

A continuación, obtenemos los valores actuales de todos los hiperparámetros del modelo en el array de hiperparámetros theta. Estos valores (obtenidos a partir del núcleo y la función de verosimilitud) nos servirán como punto de partida para la búsqueda del óptimo. Para cada hiperparámetro, se especifican las escalas (s), los límites inferior (bndl) y superior (bndu). Los límites impiden la búsqueda de soluciones en regiones mal condicionadas o sin sentido (por ejemplo, longitudes de escala o varianzas negativas). El optimizador usa el escalado para normalizar los parámetros, lo que mejora la estabilidad y la velocidad de convergencia, especialmente cuando los órdenes de magnitud de los parámetros difieren considerablemente. Valor predeterminado s = 1.0

A continuación, declaramos un array de punteros kernels_to_process a la interfaz Ikernel. Se usará para almacenar una lista de todos los núcleos base cuyos hiperparámetros necesitan ser optimizados. Si tenemos un núcleo simple (no compuesto), entonces este array tendrá un solo elemento: el puntero a dicho núcleo. Si se trata de un SumKernel o un ProductKernel, almacenará los punteros a todos los núcleos que forman parte de esta composición.

A continuación, usando el operador dynamic_cast, comprobamos si el kernel actual m_kernel (que es un campo de la clase GaussianProcess y apunta al kernel seleccionado por el usuario) es una instancia de SumKernel o ProductKernel. En ese caso, se produce una conversión de tipo a SumKernel o ProductKernel, y se llama al método GetKernels(), que llena el array kernels_to_process con todos los kernels que forman parte de la suma o el producto de kernels. Si el kernel no es ni una suma ni un producto (es decir, es un kernel regular, como un RBFKernel), entonces el array kernels_to_process simplemente contiene un puntero al propio m_kernel.

Después de eso, iteramos cada kernel básico encontrado en kernels_to_process y establecemos sus hiperparámetros para escalar s y limitar bndl y bndu.

Finalmente, tras procesar todos los hiperparámetros del kernel, se procesan los hiperparámetros de la función de verosimilitud. La función de verosimilitud gaussiana tiene un parámetro, mientras que la función de verosimilitud logit no tiene parámetros. Una vez preparados todos los parámetros, comienza el proceso de optimización.

• El método CalculateNLMLObjective() actúa como un enlace entre la clase principal GaussianProcess y el optimizador externo Alglib. Esta es precisamente la función objetivo que el optimizador MinBleic llama constantemente (a través de la clase GPOptimizationObjective) para evaluar los valores actuales de los hiperparámetros. Su tarea principal consiste en devolver el valor NLML para un conjunto dado de hiperparámetros.

//+-------------------------------------------------------------------+
//| Method that will be called by the optimizer to calculate NLML     |
//+-------------------------------------------------------------------+
double GaussianProcess::CalculateNLMLObjective(const vector &hyperparameters)
{
    //  Set all hyperparameters (kernels and likelihoods)
    SetHyperparameters(hyperparameters);    
    //  Call the inference function, which will calculate NLML
    m_inference.Infer(m_X_train, m_y_train, m_kernel, m_likelihood,m_last_inference_result);
    if (!m_last_inference_result.success) {    
        Print("Inference Error !");
        return DBL_MAX;    
    }    
    return m_last_inference_result.nlml_value;
}

En cada iteración, el optimizador MinBleic ofrece un nuevo conjunto de hiperparámetros. Lo primero que hace CalculateNLMLObjective() es tomar este conjunto (hiperparámetros) y usar el método SetHyperparameters() para actualizar los parámetros correspondientes dentro de los objetos kernel (m_kernel) y función de verosimilitud (m_likelihood). Esto resulta muy importante porque todos los cálculos NLML posteriores deben basarse en estos valores de hiperparámetros actuales.

Tras actualizar los hiperparámetros, el método llama a Infer() en el objeto de inferencia (m_inference). Este es el paso principal donde se realizan todos los cálculos matemáticos complejos destinados a estimar la posterior distribución.

Los resultados de la inferencia, incluido el valor NLML y sus gradientes (que usará la función Grad), se almacenan en el campo de clase privado m_last_inference_result.

Si la inferencia es exitosa, el método retorna NLML.

• El método GaussianProcess::SetHyperparameters(const vector &params) se encarga de distribuir y establecer valores optimizados de los hiperparámetros del kernel y la función de verosimilitud. 

//+------------------------------------------------------------------+    
//| Method for setting hyperparameters                               |
//+------------------------------------------------------------------+
void GaussianProcess::SetHyperparameters(const vector &params)
{
//+------------------------------------------------------------------+    
//This is a call of the polymorphic SetHyperparameters method on the object pointed to by m_kernel.
//Since m_kernel is a pointer to a base type (IKernel*), calling SetHyperparameters
//will be redirected to a concrete implementation of this method in the derived kernel class
//m_kernel refers to. For example, if m_kernel actually points to an object
//RBFKernel, RBFKernel::SetHyperparameters(params) is called. If this is SumKernel, 
//the SumKernel::SetHyperparameters(params) method is called, and so on.    
//+------------------------------------------------------------------+
    int kernel_params_count = m_kernel.GetNumHyperparameters();
    int likelihood_params_count = m_likelihood.GetNumHyperparameters();
    // Set kernel parameters
    vector kernel_hps(kernel_params_count);
    for(int i = 0; i < kernel_params_count; i++) {
        kernel_hps[i] = params[i];
    }
    m_kernel.SetHyperparameters(kernel_hps);    
    // Set the likelihood parameters
    vector likelihood_hps(likelihood_params_count);
    for(int i = 0; i < likelihood_params_count; i++) {
        likelihood_hps[i] = params[kernel_params_count + i];
    }
    m_likelihood.SetHyperparameters(likelihood_hps);
}

El vector de parámetros contiene todos los hiperparámetros del modelo GP en un orden fijo: primero vienen todos los hiperparámetros del núcleo (o núcleos, si es un núcleo compuesto) y luego vienen los parámetros de la función de verosimilitud. La característica clave de este método es el uso del polimorfismo. La misma llamada a m_kernel.SetHyperparameters() se comporta de forma distinta dependiendo del tipo real del objeto al que apunta m_kernel en tiempo de ejecución. 

• El método Predict(). En esencia, para eso se crea un modelo: para hacer predicciones basadas en nuevos datos. 

//+------------------------------------------------------------------+
//| Prediction method for regression and classification              |
//+------------------------------------------------------------------+    
bool GaussianProcess::Predict(const matrix &X_test, GPPredictionResult &result,PredictMode predict_mode)
{     
    // 1. Check that the model has been trained
    if (!m_last_inference_result.success) {
        Print("Error: Predict - Inference results not available");
        return false;
    }
    // 1.1 Check the match of the number of features
    if (X_test.Cols() != m_X_train.Cols()) {
        Print("Error: Predict - Number of features in X_test  must match X_train ");
        return false;
    }

    int N_train = (int)m_X_train.Rows();
    int N_test = (int)X_test.Rows();

    // 2. K_s and K_ss are calculated regardless of the type of inference/likelihood
    matrix K_s = m_kernel.Compute(m_X_train, X_test);            
    matrix K_ss = m_kernel.Compute(X_test, X_test);
    
    // --- 3. Logic for calculating mu_f_star and Sigma_f_star (common for both types of problems) ---
    //------------------------- Algorithm 2.1 GPML----------------------------------------
    if (m_inference.GetName() == "ExactInference") {
        // For ExactInference
        matrix L_K_noisy = m_last_inference_result.L_K_noisy;
        vector alpha = m_last_inference_result.alpha;
        
        result.mu_f_star = K_s.Transpose() @ alpha;            

        matrix V(N_train, N_test);    
        if (!L_K_noisy.LinearEquationsSolution(K_s, V)) {
            Print("Error: Predict (Exact) - LinearEquationsSolution failed");
            return false;        
        }          
        result.Sigma_f_star = K_ss - V.Transpose() @ V;

    } else if (m_inference.GetName() == "LaplaceInference") {
    //------------------------- Algorithm 3.2 GPML ----------------------------------------
        matrix W = -1 * m_last_inference_result.H;    
        matrix L_B = m_last_inference_result.L_B;    
        matrix sW = m_last_inference_result.sW;        
        vector f_hat = m_last_inference_result.mu_f_train;    
        vector grad_f_hat = m_likelihood.LogLikelihoodGradient(f_hat, m_y_train);
        // Eq[f*∣X,y,x*]=k(x*)^T K^−1 f_hat = k(x*)^T ∇log p(y∣f_hat) 
        result.mu_f_star = K_s.Transpose() @ grad_f_hat;
        
        matrix SwKs = sW @ K_s;
        matrix V(N_train, N_test);    
        if (!L_B.LinearEquationsSolution(SwKs, V)) {
            Print("Error: Predict (Laplace) - LinearEquationsSolution failed");
            return false;
        }
        // Vq[f*|X, y,x*] = Kss - Ks^T(K + W^-1)^-1 Ks
        result.Sigma_f_star = K_ss - V.Transpose() @ V;
    }    
    
    // --- 4. Likelihood-specific logic (Likelihood) ---
if (m_likelihood.GetName() == "GaussianLikelihood") {
    // --- 4.1. Regression (GaussianLikelihood) ---
    double noise_variance = 0.0;
    vector likelihood_params = m_likelihood.GetHyperparameters();
    if (likelihood_params.Size() > 0) {
        noise_variance = likelihood_params[0] * likelihood_params[0];
    }    
    result.Sigma_y_star = result.Sigma_f_star + matrix::Identity(N_test, N_test) * noise_variance;
    result.mu_y_star = result.mu_f_star; // For Gaussian likelihood mu_y_star = mu_f_star

    } else if (m_likelihood.GetName() == "LogitLikelihood") {
        // --- 4.2. Classification (LogitLikelihood) ---
          // Make sure m_likelihood is a LogitLikelihood to access the sigmoid method
        LogitLikelihood *logit = dynamic_cast<LogitLikelihood*>(m_likelihood);
        if (logit == NULL) {
            Print("Error: Failed to cast m_likelihood to LogitLikelihood in Predict");
            return false;
        }
        
        result.predicted_probabilities.Resize(N_test);
        result.predicted_labels.Resize(N_test);
        double mc_samples_array[]; 
 
        for (int i = 0; i < N_test; i++) {
            double mu_f_star_i = result.mu_f_star[i];  //mean of the posterior distribution q(f*|X,y,x*)
            double sigma_f_star_diag_i = result.Sigma_f_star[i, i]; // variance of the posterior distribution q(f*|X,y,x*)
    //------------------- 1)Probit Approximation----------------------
            if (predict_mode == PROBIT) {
                double k_i = 1.0 / MathSqrt(1.0 + M_PI / 8.0 * sigma_f_star_diag_i);
                result.predicted_probabilities[i] = logit.sigmoid(mu_f_star_i * k_i);}
    // ----------------- 2) Numerical integration ---------------------------------------
            else if (predict_mode == NUM_INTEGR) {  
                result.predicted_probabilities[i] = CalculateNumericalProbability(
                    mu_f_star_i,
                    sigma_f_star_diag_i,
                    logit
                );} 
   // ----------------------3) Monte Carlo Method ---------------------------------------              
            else if (predict_mode == MONTE_CARLO) {      
                // Number of samples for Monte Carlo
                int num_samples = 10000;                
                ArrayResize(mc_samples_array, num_samples); 
                double std_dev_f_star_i = MathSqrt(sigma_f_star_diag_i);
                // Generate num_samples values from N(mu_f_star_i, std_dev_f_star_i)             
                MathRandomNormal(mu_f_star_i, std_dev_f_star_i, num_samples, mc_samples_array);
                double sum_sigmoid_samples = 0.0;
                for (int s = 0; s < num_samples; s++) {
                    sum_sigmoid_samples += logit.sigmoid(mc_samples_array[s]);
                }
            //To get the expected probability p(y*=+1|X,y,x*)
            //we calculate the arithmetic mean of all obtained values σ(f_sample*). 
            //By the law of large numbers, when num_samples is large enough, 
            //this average will be a good approximation of the true value of the integral     
                result.predicted_probabilities[i] = sum_sigmoid_samples / num_samples;
            }    
            // Predicted marks (+1 or -1)
            result.predicted_labels[i] = (result.predicted_probabilities[i] >= 0.5) ? 1.0 : -1.0;
        }
    }
    return true;    
    }

Los resultados de la predicción (media, varianza, probabilidades, etiquetas de clase) se escriben en la estructura GPPredictionResult.

En primer lugar, calculamos las matrices K* y K**. Estas matrices suponen la base para las predicciones en GP, y son necesarias para calcular la media y la varianza de la función latente f* en nuevos puntos de prueba. La lógica en este caso depende del método de inferencia (ExactInference o LaplaceInference) utilizado durante el entrenamiento, ya que proporcionan componentes diferentes para las fórmulas de predicción (Algoritmo 2.1 para Exact, Algoritmo 3.2 para Laplace del libro GPML de Rasmussen y Williams).

Si hemos utilizado ExactInference, se recuperan los valores precalculados de L_K_noisy y alpha. Si hemos utilizado LaplaceInference, se extraen W, L_B, sW y f_hat (moda). En ambos casos, el resultado es la media (mu_f_star) y la matriz de covarianza (Sigma_f_star) de la función latente para cada punto de prueba.

Como ya hemos comentado en la parte teórica del artículo, existe un problema al calcular la integral para obtener la probabilidad de clase. Por lo tanto, se usan aproximaciones:

• predict_mode == PROBIT (aproximación Probit):

Esta es una aproximación rápida que se usa con frecuencia. Sustituye la función sigmoide por la función de distribución acumulativa de la distribución normal, que tiene una forma similar. Esto nos permite calcular la integral analíticamente.

• predict_mode == NUM_INTEGR (Integración numérica):

En este modo, se llama a la función auxiliar CalculateNumericalProbability. Esta se aproxima numéricamente a la integral dividiendo el rango de f* en intervalos discretos y sumando los valores. Puede que sea más preciso, pero más lento. 

• predict_mode == MONTE_CARLO (Método de Monte Carlo):

Este es un método estocástico. Se genera un gran número de muestras aleatorias f* partiendo de la distribución posterior q(f*∣X, y, x*). Para cada muestra f*, se calcula sigma(f*).

La media aritmética de todos estos valores sigma(f*) es una aproximación de la probabilidad deseada p(y*=+1|X, y, x*). Este es el método que requiere mayor capacidad de cálculo. Para generar muestras a partir de una distribución normal, se usa la función de la biblioteca estándar MathRandomNormal.

En función de las probabilidades calculadas, se toma una decisión sobre la etiqueta de clase prevista para cada uno de los métodos de aproximación anteriormente mencionados. Si la probabilidad de pertenecer a la clase +1 es mayor o igual a 0,5, entonces se predice +1; de lo contrario, -1.

La clase GPOptimizationObjective

//+------------------------------------------------------------------+
//| Class for the Alglib optimizer objective function                |
//+------------------------------------------------------------------+
class GPOptimizationObjective : public CNDimensional_Grad
{
private:  
    GaussianProcess* m_gp; // pointer to GaussianProcess object
    double nlml;           // Negative log-likelihood
public:
    // Constructor 
    GPOptimizationObjective(GaussianProcess* gp_instance) : m_gp(gp_instance), nlml(0.0) {}
    double GetNLML() { return nlml; }
    ~GPOptimizationObjective() {}

    // Grad method that will be called by the optimizer
    virtual void Grad(CRowDouble &w, double &func,CRowDouble &grad, CObject &obj) override {     
        // Convert CRowDouble to a vector for passing to GP
        vector hyperparameters(w.Size());
        for(int i = 0; i < (int)w.Size(); i++) { 
           hyperparameters[i] = w[i];        
        }
        
        // Call the GP method to calculate NLML
        func = m_gp.CalculateNLMLObjective(hyperparameters);
        nlml = func; 
        
         GPInferenceResult current_result = m_gp.GetLastInferenceResult(); 
        if (!current_result.success ) {
            Print("Warning: GPOptimizationObjective::Grad - Gradient calculation failed");
            for(int i = 0; i < (int)w.Size(); i++) {
                grad.Set(i, DBL_MAX); 
            }
            return;
        }            
        // Fill grad with elements from current_result.nlml_gradient
        for(int i = 0; i < (int)w.Size(); i++) {
        grad.Set(i, current_result.nlml_gradient[i]);
        }       
    }      
};

Esta clase actúa como puente entre nuestra clase GaussianProcess y la biblioteca de optimización externa Alglib (en concreto, el optimizador MinBLEIC).

Alglib requiere que la función objetivo que optimiza se ajuste a una interfaz determinada. Y para eso precisamente sirve GPOptimizationObjective. Esta hereda de CNDimensional_Grad, la clase básica de Alglib que define esta interfaz. Esta clase base ofrece métodos virtuales que la clase GPOptimizationObjective debe implementar. Estos métodos permiten que los optimizadores de Alglib trabajen con cualquier función objetivo, siempre que ofrezca tanto el valor de la función como su gradiente. 

El miembro privado GaussianProcess* m_gp contiene el puntero a nuestro objeto GaussianProcess. Esto permite que la clase GPOptimizationObjective llame al método CalculateNLMLObjective para realizar los cálculos necesarios. 

El método Grad() es la parte más importante de esta clase: sobrescribe el método virtual de CNDimensional_Grad y será llamado por el optimizador Alglib en cada iteración. La función Grad() recibe el vector de hiperparámetros actual w de Alglib y debe retornar el valor de la función objetivo func y el vector de sus gradientes grad. 

Conclusión

Vamos a resumir los resultados provisionales.

En la primera parte del artículo, hemos sentado una sólida base teórica para comprender el modelo de clasificación GP. Asimismo, hemos examinado con detalle los principios de funcionamiento de los procesos gaussianos para la clasificación binaria y el método de aproximación de Laplace. Este método resulta de vital importancia porque hace que el problema de clasificación sea práctico y computacionalmente eficiente para las necesidades del trading en línea, a diferencia del método MCMC, que es preciso pero increíblemente costoso.

Una vez abordados los constructos teóricos, hemos pasado a la implementación práctica, diseñando y describiendo dos clases clave de nuestra biblioteca de GP:

• GaussianProcess: la clase principal que encapsula toda la lógica para construir, entrenar y predecir un modelo GP,
• GPOptimizationObjective: actúa como un proxy, preparando la función objetivo y su gradiente en el formato requerido por la biblioteca Alglib para la optimización de hiperparámetros.

En la segunda parte, completaremos la implementación de la biblioteca proporcionando al lector:

• Una descripción detallada y código de implementación de las interfaces clave: IKernel (para diferentes núcleos), IInference (para métodos de inferencia) e ILikelihood (para funciones de verosimilitud);
• Ejemplos del funcionamiento de la biblioteca con datos sintéticos para demostrar claramente sus capacidades;
• Una aplicación práctica en el trading: desarrollaremos indicadores de clasificación y regresión basados en nuestra biblioteca, demostrando cómo se pueden utilizar los GP para tomar decisiones de trading.