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.

//+------------------------------------------------------------------+
//| Класс Гауссовского Процесса                                      |
//+------------------------------------------------------------------+
class GaussianProcess
{
private:
    IKernel*      m_kernel;       // указатель на выбранное ядро
    ILikelihood*  m_likelihood;   // указатель на выбранную функцию правдоподобия
    IInference*   m_inference;    // указатель на выбранный метод инференса
    matrix        m_X_train;      // Тренировочные входные данные Nxd
    vector        m_y_train;      // Тренировочные целевые данные Nx1

    GPInferenceResult m_last_inference_result; // Структура сохраняющая последние результаты инференса
    
    int m_last_termination_type;  // Код завершения операции оптимизации
    int m_last_iterations_count;  // Количество итераций выполненных оптимизатором
    double m_last_nlml_value;     // Итоговое значение NLML после оптимизации
    
private:
    // Вспомогательная функция для численного интегрирования 
    double CalculateNumericalProbability(double mu_f_star, double sigma_f_star_diag, LogitLikelihood *likelihood);

public:
    // Конструктор класса 
    GaussianProcess(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
    const matrix &X_train, const vector &y_train);
    // Статический метод для создания обьекта GaussianProcess с проверкой входных параметров
    static GaussianProcess* Create(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
    const matrix &X_train, const vector &y_train);
    
    // Деструктор
    ~GaussianProcess();
    
    // --- Методы для получения состояния модели ---
    // Возвращает результаты последней операции инференса
    GPInferenceResult GetLastInferenceResult() const;
    // Возвращает тип завершения последней оптимизации гиперпараметров
    int GetLastTerminationType() const;
    // Возвращает количество итераций, выполненных при последней оптимизации гиперпараметров
    int GetLastIterationsCount() const;
    // Возвращает значение отрицательного логарифма маргинального правдоподобия после оптимизации
    double GetLastNLML() const;
    // Возвращает указатель на используемое ядро
    IKernel* GetKernel() const;
    // Возвращает текущие значения всех оптимизируемых гиперпараметров.
    vector GetCurrentHyperparameters();
  
    // --- Методы обучения и настройки ---
    // Запускает полный процесс обучения модели, включая оптимизацию гиперпараметров
    bool Fit();
    // Выполняет одиночный шаг инференса без оптимизации гиперпараметров
    bool PerformInference();
    // Устанавливает тренировочные данные для модели 
    void SetTrainingData(const matrix& X, const vector& y);
    // Устанавливает заданные гиперпараметры для ядра и функции правдоподобия
    void SetHyperparameters(const vector &params);
    // Метод, вызываемый оптимизатором для вычисления целевой функции (NLML)
    double CalculateNLMLObjective(const vector &hyperparameters);
    
    // --- Метод выполняет прогноз для новых тестовых данных
    // Параметр predictmode определяет метод расчета вероятностей для классификации (PROBIT, NUM_INTEGR, MONTE_CARLO)
    bool Predict(const matrix &X_test, GPPredictionResult &result, PredictMode mode = PROBIT);
    
    // --- Вспомогательные методы ---
    // Статический метод для генерации выборок из априорного ГП
    static bool SamplePriorGP(const matrix &x, IKernel* kernel, int num_samples, matrix &f_samples,
                              bool plot_samples = false, int plot_display_seconds = 10);
    //--- Метод для вывода в журнал итоговых значений гиперпараметров
    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                                                     |
//+------------------------------------------------------------------+
GaussianProcess* GaussianProcess::Create(IKernel* kernel, ILikelihood* likelihood, IInference* inference,
const matrix &X_train, const vector &y_train)
{
    // 1. Проверка на NULL-указатели
    if (kernel == NULL || likelihood == NULL || inference == NULL) {
        Print("ERROR: Kernel, Likelihood, or Inference pointer is NULL");
        return NULL;
    } 
    // 2. Проверка корректности входных данных X_train и y_train 
    if (X_train.Rows() == 0 || y_train.Size() == 0 || X_train.Rows() != y_train.Size()) {
    Print("ERROR: Invalid training data dimensions");
    return NULL;
    }
    // 3. Проверка совместимости likelihood и 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. Если все проверки пройдены, создаем объект
    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                               |
//+------------------------------------------------------------------+
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).

//+------------------------------------------------------------------+
//| Метод для обучения модели                                        |
//+------------------------------------------------------------------+
bool GaussianProcess::Fit()
{
    // Создаем обьект GPOptimizationObjective, передавая ему указатель на текущий объект GaussianProcess
    // Этот указатель попадает в приватное поле класса m_gp с помощью которого мы вызываем метод
    // CalculateNLMLObjective,чтобы получить значение NLML для текущего набора гиперпараметров
    GPOptimizationObjective objective_func(GetPointer(this));
    CNDimensional_Rep frep; 
    CObject Obj;
    vector initial_hyperparams = GetCurrentHyperparameters(); // Получаем стартовые значения гиперпараметров
    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[]; // массив указателей на интерфейс IKernel
    
    // Логика получения ядер для установки границ 
    // Этот блок кода определяет, с каким типом ядра мы имеем дело,
    // и заполняет массив kernels_to_process соответствующими указателями:
    if (dynamic_cast<SumKernel*>(m_kernel) != NULL) {          // Проверяем, является ли текущее ядро m_kernel обьектом SumKernel 
        SumKernel* sum_k = dynamic_cast<SumKernel*>(m_kernel); // Если да, то приводим тип m_kernel к типу SumKernel*, 
        sum_k.GetKernels(kernels_to_process); // и вызываем метод GetKernels(), который заполняет массив kernels_to_process всеми ядрами, входящими в сумму
    } else if (dynamic_cast<ProductKernel*>(m_kernel) != NULL) { // Аналогичная логика, если ядро является обьектом ProductKernel
        ProductKernel* prod_k = dynamic_cast<ProductKernel*>(m_kernel);
        prod_k.GetKernels(kernels_to_process);
    } else {
        ArrayResize(kernels_to_process,1); // Если ядро не является ни суммой, ни произведением (т.е. это не композитное ядро), 
        kernels_to_process[0] = m_kernel; // то массив kernels_to_process просто содержит указатель на m_kernel.
    }

   // Этот цикл проходит по каждому базовому ядру, найденному в массиве kernels_to_process,
   // и устанавливает для его гиперпараметров начальный масштаб s, нижнюю bndl и верхнюю 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++;    
                } 
            }           
    }
    
// --- Добавляем границы и масштабы для параметров правдоподобия (если они есть) ---
// LogitLikelihood не имеет гиперпараметров, поэтому для него этот блок будет пропущен
// GaussianLikelihood имеет 1 параметр (sigma)
if (m_likelihood.GetNumHyperparameters() > 0) {
    if (param_idx + m_likelihood.GetNumHyperparameters() <= num_params) {
        s[param_idx] = 1.0;           // Масштаб
        bndl[param_idx] = 1e-10;      // Нижняя граница 
        bndu[param_idx] = 1e3;        // Верхняя граница
        param_idx++;
    } 
}
    CMinBLEICStateShell state;
    CMinBLEICReportShell rep; // объект, который будет содержать отчёт о результатах оптимизации
    //-----------------------  критерии остановки оптимизатора
    double epsg = 0.0001;     //Точность по градиенту (0 означает, что остановка по градиенту отключена)
    double epsf = 0.0000;     //Точность по значению функции    
    double epsw = 0.0000;     //Точность по параметрам 
    //-------------------------   
    double epso = 0.00001;    //Параметры для внешних  условий сходимости в BLEIC
    double epsi = 0.00001;    //Параметры для внутренних условий сходимости в BLEIC
    CAlglib::MinBLEICCreate(theta, state);       // инициализация оптимизатора. Она создаёт начальное состояние state для MinBLEIC, используя начальные значения гиперпараметров из массива theta  
    CAlglib::MinBLEICSetBC(state, bndl, bndu);   // Устанавливает нижние (bndl) и верхние (bndu) границы для каждого параметра
    CAlglib::MinBLEICSetScale(state, s);         //Устанавливает масштабы (s) для каждого параметра. Это может помочь оптимизатору эффективнее работать с параметрами разных порядков величины.
    CAlglib::MinBLEICSetInnerCond(state,epsg,epsf,epsw);    
    CAlglib::MinBLEICSetOuterCond(state, epso, epsi);    
    CAlglib::MinBLEICOptimize(state, objective_func, frep, 0, Obj); // запускаем процесс оптимизации    
    CAlglib::MinBLEICResults(state, theta, rep); // отчет оптимизации
    
    m_last_termination_type = rep.GetTerminationType();
    m_last_iterations_count = rep.GetInnerIterationsCount();
    m_last_nlml_value = objective_func.GetNLML(); // Получаем финальное 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.
//-------------------------------------------------------------------------------------  
// Определяем успешность оптимизации на основе TerminationType
    bool success = true;     
    if (m_last_termination_type < 0)
    {
        Print("Ошибка: Оптимизация GP завершилась неудачно. Тип завершения: ", m_last_termination_type);
        success = false;
    } 
    // Обновляем гиперпараметры модели после оптимизации
    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.

//+-------------------------------------------------------------------+
//| Метод, который будет вызываться оптимизатором для вычисления NLML |
//+-------------------------------------------------------------------+
double GaussianProcess::CalculateNLMLObjective(const vector &hyperparameters)
{
    //  Устанавливаем все гиперпараметры (ядра и правдоподобия)
    SetHyperparameters(hyperparameters);    
    //  Вызываем инференс, который вычислит 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("Ошибка Инференса !");
        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. 

//+------------------------------------------------------------------+    
//| Метод для установки гиперпараметров                              |
//+------------------------------------------------------------------+
void GaussianProcess::SetHyperparameters(const vector &params)
{
//+------------------------------------------------------------------+    
//Это вызов полиморфного метода SetHyperparameters у объекта, на который указывает m_kernel.
//Поскольку m_kernel является указателем базового типа (IKernel*), вызов SetHyperparameters
//будет перенаправлен на конкретную реализацию этого метода в производном классе ядра,
//к которому относится m_kernel.Например,если m_kernel на самом деле  указывает на объект
//RBFKernel,будет вызвана RBFKernel::SetHyperparameters(params).Если это SumKernel, 
//будет вызван метод SumKernel::SetHyperparameters(params), и так далее.    
//+------------------------------------------------------------------+
    int kernel_params_count = m_kernel.GetNumHyperparameters();
    int likelihood_params_count = m_likelihood.GetNumHyperparameters();
    // Устанавливаем параметры ядра
    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);    
    // Устанавливаем параметры правдоподобия
    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. 

//+------------------------------------------------------------------+
//| Метод предсказания для регрессии и классификации                 |
//+------------------------------------------------------------------+    
bool GaussianProcess::Predict(const matrix &X_test, GPPredictionResult &result,PredictMode predict_mode)
{     
    // 1. Проверка, что модель была обучена
    if (!m_last_inference_result.success) {
        Print("Error: Predict - Inference results not available");
        return false;
    }
    // 1.1 Проверка совпадения количества признаков
    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 и K_ss вычисляются независимо от типа инференса/правдоподобия
    matrix K_s = m_kernel.Compute(m_X_train, X_test);            
    matrix K_ss = m_kernel.Compute(X_test, X_test);
    
    // --- 3. Логика для расчета mu_f_star и Sigma_f_star (общая для обоих типов задач) ---
    //------------------------- Алгоритм 2.1 GPML----------------------------------------
    if (m_inference.GetName() == "ExactInference") {
        // Для 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") {
    //------------------------- Алгоритм 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) ---
if (m_likelihood.GetName() == "GaussianLikelihood") {
    // --- 4.1. Регрессия (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; // Для гауссовского правдоподобия mu_y_star = mu_f_star

    } else if (m_likelihood.GetName() == "LogitLikelihood") {
        // --- 4.2. Классификация (LogitLikelihood) ---
          // Убедимся, что m_likelihood является LogitLikelihood для доступа к методу sigmoid
        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];  //среднее апостериорного распределения q(f*|X,y,x*)
            double sigma_f_star_diag_i = result.Sigma_f_star[i, i]; // дисперсия апостериорного распределения 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) Численное интегрирование ---------------------------------------
            else if (predict_mode == NUM_INTEGR) {  
                result.predicted_probabilities[i] = CalculateNumericalProbability(
                    mu_f_star_i,
                    sigma_f_star_diag_i,
                    logit
                );} 
   // ----------------------3) Метод Монте Карло ---------------------------------------              
            else if (predict_mode == MONTE_CARLO) {      
                // Количество сэмплов для Монте-Карло
                int num_samples = 10000;                
                ArrayResize(mc_samples_array, num_samples); 
                double std_dev_f_star_i = MathSqrt(sigma_f_star_diag_i);
                // Генерируем num_samples значений из 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]);
                }
            //Чтобы получить ожидаемую вероятность  p(y*=+1|X,y,x*)
            //мы вычисляем среднее арифметическое всех полученных значений σ(f_sample*). 
            //По закону больших чисел, когда num_samples достаточно велико, 
            //это среднее будет хорошей аппроксимацией истинного значения интеграла     
                result.predicted_probabilities[i] = sum_sigmoid_samples / num_samples;
            }    
            // Предсказанные метки (+1 или -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

//+------------------------------------------------------------------+
//| Класс для целевой функции оптимизатора Alglib                    |
//+------------------------------------------------------------------+
class GPOptimizationObjective : public CNDimensional_Grad
{
private:  
    GaussianProcess* m_gp; // указатель на объект GaussianProcess
    double nlml;           // Отрицательное логарифмическое правдоподобие
public:
    // Конструктор 
    GPOptimizationObjective(GaussianProcess* gp_instance) : m_gp(gp_instance), nlml(0.0) {}
    double GetNLML() { return nlml; }
    ~GPOptimizationObjective() {}

    // Метод Grad, который будет вызван оптимизатором
    virtual void Grad(CRowDouble &w, double &func,CRowDouble &grad, CObject &obj) override {     
        // Преобразуем CRowDouble в vector для передачи в GP
        vector hyperparameters(w.Size());
        for(int i = 0; i < (int)w.Size(); i++) { 
           hyperparameters[i] = w[i];        
        }
        
        // Вызываем метод GP для вычисления 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;
        }            
        // Заполняем grad элементами из 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.