Processos gaussianos em machine learning (Parte 1): modelo de classificação em MQL5

Introdução

Damos continuidade ao nosso estudo do modelo de machine learning conhecido como processos gaussianos (PG). No artigo anterior, analisamos em detalhe o problema de regressão, no qual o objetivo principal era prever valores contínuos. Hoje, vamos tratar de um tema muito mais complexo: a classificação. Sua principal dificuldade está no fato de que o processo de inferência (inference) para classificação em processos gaussianos não possui solução analítica, o que exige o uso de métodos aproximados, como a aproximação de Laplace.

Para resolver essa tarefa complexa de forma eficiente, desenvolveremos uma biblioteca modular de processos gaussianos em MQL5. Essa abordagem permitirá estruturar o código, dividindo o modelo de PG em componentes independentes, e fornecerá uma base sólida para melhorias e expansões futuras. Essa biblioteca se tornará uma ferramenta universal tanto para tarefas de regressão quanto de classificação.

Nesta primeira parte do artigo, analisaremos em detalhes a teoria da classificação com PG, incluindo a matemática que serve de base para os métodos aproximados. Também apresentaremos a classe principal da biblioteca, GaussianProcess, que reunirá todos os componentes do modelo, e a classe GPOptimizationObjective, responsável pela integração com a biblioteca de otimização Alglib.

Classificação

A classificação é uma tarefa de machine learning que consiste em atribuir um objeto a uma das categorias predefinidas. Por exemplo, no setor financeiro, a classificação pode ajudar a prever se o preço de uma ação vai subir ou cair com base em dados históricos.

Neste artigo, vamos nos concentrar na classificação binária, na qual o objeto pertence a uma de duas classes, por exemplo, "alta" (+1) ou "queda" (-1). Diferentemente de métodos como máquinas de vetores de suporte (SVM) ou árvores de decisão, que retornam apenas o rótulo da classe, os PG permitem obter uma previsão probabilística. Por exemplo, o modelo pode indicar que a probabilidade de alta de uma ação é de 75%. Esse tipo de informação é especialmente valioso no trading, quando o grau de confiança na previsão ajuda a tomar decisões mais fundamentadas, permitindo filtrar sinais pouco confiáveis. 

Infelizmente, resolver o problema de classificação com o uso de PG é significativamente mais complexo do que no caso da regressão. Isso está relacionado ao tipo de verossimilhança utilizado: 

• Na regressão, em geral, utiliza-se a verossimilhança gaussiana. A combinação dos PG, como distribuição a priori da função, com a verossimilhança gaussiana permite obter analiticamente a distribuição a posteriori, o que simplifica todos os cálculos.
• Na classificação, em que os alvos são rótulos discretos de classe, a verossimilhança gaussiana não é adequada. Em seu lugar, utiliza-se, por exemplo, a verossimilhança logística. Isso leva ao fato de que a distribuição a posteriori também deixa de ser gaussiana e não possui solução analítica.

Como consequência, precisamos recorrer a métodos complexos de inferência aproximada (approximate inference). A ideia central desses métodos é aproximar a verdadeira distribuição a posteriori não gaussiana por uma distribuição gaussiana centrada em sua moda.  No âmbito deste artigo, vamos nos concentrar na aproximação quadrática de Laplace (Laplace Approximation), pois ela é uma das abordagens mais simples e eficazes para obter uma aproximação gaussiana da distribuição a posteriori. 

Para a classificação binária, a ideia fundamental da previsão com PG é bastante simples.  Começamos com a distribuição a priori das funções latentes f(x). Imagine que o PG gera não uma única função, mas um conjunto infinito de funções possíveis, cada uma das quais é uma potencial dependência "latente" nos dados. Em seguida, cada uma dessas possíveis realizações da função latente f(x) é "passada" pela função logística, ou sigmoide. A sigmoide transforma qualquer número real, isto é, o valor de f(x), em uma probabilidade entre 0 e 1, que será nossa probabilidade a priori π(x) de pertencer à classe "+1":

É importante observar que π é uma função determinística de f, mas, como a própria f é estocástica, aleatória, uma amostra de um PG, então a função π também se torna estocástica. Esse conceito é ilustrado de forma clara nas A Fig. 1 e 2 para um espaço de entrada unidimensional X.

Fig. 1 mostra apenas uma possível realização da função latente

A Fig. 1 mostra apenas uma possível realização da função latente, demonstrando o comportamento típico da função correspondente aos hiperparâmetros definidos do kernel.

Fig. 2 demonstra o resultado da aplicação da função logística à mesma função f(x)

A Fig. 2 demonstra o resultado da aplicação da função logística, ou sigmoide, à mesma função f(x):

Assim, obtemos a distribuição a priori da probabilidade de pertencimento à classe π(x)=σ(f(x)), que, nesta etapa, ainda não leva em conta os dados de treinamento y. Sem as observações y, essa distribuição a priori permanece apenas como nossa hipótese inicial, sem respaldo em fatos empíricos; sem elas, o modelo não dispõe de informações sobre quais de suas suposições iniciais estavam corretas e quais precisam ser revistas. 

Naturalmente, a escolha das suposições a priori influencia de forma significativa os resultados a posteriori finais. Essa é uma característica-chave da abordagem bayesiana, pois são justamente as decisões do pesquisador sobre o tipo de kernel que determinam as propriedades da distribuição a priori das funções e, consequentemente, do modelo final. 

Processo de inferência (Inference)

Portanto, para fazer previsões fundamentadas, precisamos levar em conta os dados reais de treinamento y. É exatamente aqui que entra em cena o processo de inferência (Inference). Seu objetivo principal é transformar nossas crenças a priori em crenças a posteriori, isto é, ajustadas com base nos dados observados. Para a classificação, esse processo se divide naturalmente em duas etapas sequenciais.

Passo 1: Distribuição preditiva da função latente f∗

Na primeira etapa, calculamos p(f|X, y, x), a distribuição a posteriori da função latente f para um novo ponto de teste x, condicionada aos dados de treinamento observados (X, y). Ela é definida pela seguinte integral:

onde:

• p(f∣X, x, f) é a distribuição condicional da função latente f em um novo ponto de teste x, condicionada às funções latentes f nos pontos de treinamento X. Essa distribuição é sempre normal, pois, por definição, o PG possui distribuição normal conjunta,
• p(f|X, y) é a distribuição a posteriori das funções latentes f nos dados de treinamento. Devido à função de verossimilhança não linear, a sigmoide, ela não é gaussiana.

É importante observar: como p(f|X, y) não é normal, essa integral não admite solução analítica. Isso significa que, para calculá-la, precisaremos de métodos aproximados. 

Passo 2: Probabilidade preditiva final π*

Na segunda etapa, usamos essa distribuição preditiva para formar a previsão probabilística final π, a probabilidade de o ponto de teste x pertencer à classe positiva (y* = +1):

Aqui, σ(f) é a função logística, ou sigmoide, que transforma o valor da função latente f em uma probabilidade entre 0 e 1. Já a própria integral significa que fazemos a média dessas probabilidades sobre todos os valores possíveis de f, ponderados por sua distribuição preditiva a posteriori. Em essência, essa integral unidimensional é a esperança matemática da função σ(f) em relação à distribuição p(f|X, y, x).

Novamente, para a verossimilhança logística, essa integral não admite solução analítica. Por isso, também aqui precisaremos de métodos aproximados. Adiantando um pouco, diremos que, em nossa biblioteca GP, foram implementadas três dessas aproximações, o que permite escolher o método mais adequado em função das exigências de precisão e do custo computacional:

•  aproximação probit,
•  integração numérica,
•  método de Monte Carlo.

Essas duas etapas que acabamos de descrever, o cálculo da distribuição a posteriori da função latente e a integração subsequente para obter a probabilidade preditiva, representam o esquema geral da inferência bayesiana em PG. São exatamente essas duas integrais que precisamos calcular para obter a previsão desejada, e ambas exigem o uso de métodos aproximados. 

Aproximação de Laplace

Como já vimos, na inferência bayesiana para classificação surgem integrais que não admitem solução em forma fechada. A aproximação de Laplace resolve esse problema ao aproximar a distribuição não gaussiana p(f∣X, y) por uma distribuição gaussiana q(f∣X, y). Como a distribuição condicional p(f∣X, x, f) também é gaussiana, a distribuição preditiva resultante p(f∣X, y, x) também se torna gaussiana. Isso nos permite derivar fórmulas analíticas para a média e a variância de f*, o que simplifica significativamente os cálculos posteriores. Assim, toda a elegância e a eficiência computacional da aproximação de Laplace estão em sua capacidade de reduzir o cálculo da distribuição a posteriori e das previsões a operações com distribuições gaussianas.

É importante entender que a aproximação de Laplace é um compromisso. Ela torna computacionalmente tratável um problema que não admite solução em forma fechada, mas à custa da precisão na representação da forma real da distribuição a posteriori. A qualidade dessa aproximação normal depende diretamente de quão próxima da normal está a verdadeira distribuição p(f∣X, y). Quanto mais próxima ela estiver, mais precisa será a aproximação, e vice-versa.

Se o que nos interessa é a verdadeira distribuição p(f∣X, y, x), e não sua aproximação, então normalmente são usados métodos MCMC (Markov Chain Monte Carlo). Embora o MCMC possa fornecer estimativas mais precisas, ele é computacionalmente muito custoso e complexo de implementar. O MCMC pode ser usado como padrão de referência para comparação com métodos aproximados de inferência.

Agora vamos examinar com mais detalhes o que é a aproximação de Laplace. Essa aproximação é construída em torno da moda, ou máximo, da verdadeira distribuição a posteriori p(f∣X, y). Ela utiliza a expansão de Taylor de segunda ordem do logaritmo da densidade a posteriori em torno dessa moda. Matematicamente, aproximamos o logaritmo da densidade a posteriori da seguinte forma:

onde:

• q(f∣X, y) é a aproximação gaussiana da distribuição a posteriori p(f∣X, y),
• f_hat = argmax(f) p(f|X, y) é a moda da distribuição a posteriori,
• A = −∇∇ log p(f|X, y)|f=f_hat é a Hessiana do logaritmo negativo da distribuição a posteriori no ponto da moda.

Para realizar a aproximação de Laplace, antes de tudo precisamos encontrar o valor mais provável da função latente f, isto é, a moda f_hat. Para obter a posteriori p(f∣X, y), usamos a regra de Bayes. Já sabemos que essa regra relaciona a distribuição a posteriori com a verossimilhança p(y∣f), o a priori p(f∣X) e a verossimilhança marginal p(y∣X) da seguinte forma:

Para maximizar p(f∣X, y) em relação a f, não precisamos conhecer a constante de normalização p(y∣X), pois ela não depende de f e, portanto, não afeta a posição do máximo. Por isso, podemos trabalhar com a distribuição a posteriori não normalizada, que é proporcional ao produto da verossimilhança e do a priori p(y∣f)p(f∣X).

Para simplificar os cálculos e evitar problemas numéricos causados por valores de probabilidade muito pequenos, tomamos o logaritmo dessa distribuição a posteriori não normalizada. Graças à propriedade dos logaritmos, o produto das probabilidades se transforma na soma de seus logaritmos:

Ψ(f) é a função objetivo que vamos maximizar com o método de Newton para encontrar a moda da função latente. O método de Newton exige o cálculo das derivadas de primeira e de segunda ordem de Ψ(f) em relação a f.

Ao diferenciar essa equação em relação a f, obtemos:

onde:

• W = −∇∇ log p(y|f) é a Hessiana negativa do logaritmo da verossimilhança, que é uma matriz diagonal.

Depois de calcularmos o gradiente e a Hessiana, encontramos iterativamente a moda com o método de Newton:

Em cada iteração, o método de Newton atualiza nossa estimativa atual da moda na direção determinada pelo gradiente e pela Hessiana, até que a convergência seja alcançada. 

Depois de obter a moda, podemos calcular a matriz de covariância da distribuição gaussiana aproximada.  Essa matriz é igual à Hessiana inversa negativa de Ψ(f), calculada no ponto da moda f_hat.

Assim, a matriz de covariância Σ da nossa aproximação gaussiana é igual a:

Isso conclui a descrição da primeira etapa da aproximação de Laplace: encontrar a aproximação normal da distribuição a posteriori.

Predição na aproximação de Laplace

Depois de obter q(f∣X, y), podemos passar à segunda etapa da inferência: a predição para novos pontos de teste x∗. Nesta etapa, queremos encontrar a distribuição preditiva p(f∗∣X, y, x∗). Graças à aproximação de Laplace, que tornou p(f∣X, y) gaussiana, na forma de q(f∣X, y), e ao fato de que p(f∗∣X, x∗, f) também é gaussiana, a distribuição preditiva resultante p(f∗∣X, y, x∗) também se torna gaussiana. Isso nos permite obter analiticamente sua média e sua variância a posteriori.

O valor médio da função latente f para um novo ponto de teste x (mu_f_star) é calculado como:

A variância da função latente Var(f) para um novo ponto de teste x (Sigma_f_star) é calculada como:

Agora que temos a média e a variância da distribuição preditiva, podemos, enfim, calcular a probabilidade desejada de pertencimento à classe π∗: 

Essa fórmula é o núcleo da predição de probabilidades em um classificador GP com base na aproximação de Laplace.

Talvez você tenha notado que não calculamos as probabilidades de classe simplesmente como σ(E[f∗]), isto é, substituindo diretamente a média a posteriori de f* na função sigmoide. Essa abordagem, conhecida como predição MAP (Maximum A Posteriori Prediction), certamente também é válida.

No entanto, ao calcular a predição MAP σ(E[f∗]), ignoramos a incerteza em f. Simplesmente tomamos o ponto central de f, o valor médio, e o transformamos em probabilidade. Quando calculamos E[σ(f)], o que corresponde à integração, levamos em conta toda a forma da distribuição de f. Isso nos fornece uma probabilidade preditiva mais precisa e mais significativa, especialmente quando existe incerteza considerável em f, isto é, uma grande variância V[f], ou quando a distribuição de f* é assimétrica. Essa abordagem é chamada de predições médias (averaged predictive probability).

Entender essa diferença tem importante significado prático.

• Se o seu único objetivo for obter um rótulo binário de classe, por exemplo, "comprar" ou "vender", +1 ou -1, então o uso da predição MAP mais simples pode ser suficiente, pois ela fornecerá o mesmo rótulo que a predição média, computacionalmente mais custosa.
• No entanto, se as próprias probabilidades forem importantes para você, então as predições médias, E[σ(f*)], ainda serão mais precisas, pois levam plenamente em conta a incerteza do modelo. 

No trading, um simples rótulo binário de classe, "comprar" ou "vender", não é suficiente. Precisamos de uma gradação mais refinada de confiança, e é exatamente isso que as probabilidades fornecem. O valor da probabilidade nos permite filtrar sinais de trading. Um sinal com probabilidade de sucesso de 0.51, o que é apenas ligeiramente melhor do que um palpite aleatório, terá muito menos valor do que um sinal com probabilidade de 0.60. Isso permite ao trader definir valores-limite para entrar em uma operação. Por exemplo, pode-se decidir que operações serão abertas apenas quando a probabilidade de sucesso for superior a 0.55 ou 0.60, reduzindo assim a quantidade de sinais falsos.

Verossimilhança marginal na aproximação de Laplace

Agora que entendemos o mecanismo de inferência em PG para classificação, surge a pergunta: como ajustar nosso modelo para obter previsões ideais? A resposta está na verossimilhança marginal (Marginal Likelihood, LML). Essa é a função objetivo pela qual otimizamos os hiperparâmetros θ do nosso modelo. Sem seu cálculo, é impossível encontrar os melhores parâmetros que explicam nossos dados: 

onde, B

Depois de definir a função objetivo para a otimização, o próximo passo importante é calcular suas derivadas parciais em relação aos hiperparâmetros θ. Isso é necessário porque faremos a otimização usando gradientes analíticos. Essa abordagem acelera os cálculos várias vezes em comparação com métodos numéricos. Os gradientes analíticos permitem que o otimizador avance de forma mais eficiente e precisa em direção ao mínimo da função objetivo. 

O gradiente de LML consiste em uma parte explícita e uma parte implícita:

A fórmula para calcular a parte explícita:

Aqui, o principal problema é calcular a derivada da matriz de kernel K em relação a cada hiperparâmetro. Trataremos da implementação das derivadas para a função de kernel escolhida na segunda parte do artigo.

A parte implícita consiste em dois fatores. O primeiro fator da parte implícita é obtido pela fórmula:

Para calcular essa fórmula, teremos de calcular a terceira derivada do logaritmo da verossimilhança.

O segundo fator da parte implícita é calculado da seguinte forma: 

Em conclusão, vale destacar que o NLML é necessário não apenas para estimar os hiperparâmetros, mas também para comparar diferentes modelos, por exemplo, com diferentes tipos de kernels. Modelos com menor valor de NLML são considerados melhores, pois isso significa maior verossimilhança marginal, ou seja, que o modelo explica melhor os dados observados.

Além disso, os PG, ao usarem a verossimilhança marginal para otimizar os hiperparâmetros, resolvem automaticamente o problema do equilíbrio entre ajuste aos dados e complexidade do modelo. O NLML penaliza naturalmente modelos excessivamente complexos, evitando o sobreajuste. Graças a isso, não há necessidade de critérios explícitos de parada para prevenir overfitting, como é feito, por exemplo, no treinamento de redes neurais. A otimização do NLML, por si só, busca encontrar o equilíbrio ideal. Essa é uma das principais vantagens da abordagem bayesiana nos processos gaussianos.

Biblioteca de Processos Gaussianos

Agora que analisamos todos os conceitos teóricos necessários, passamos à implementação prática. Nosso principal objetivo é criar uma biblioteca universal de PG em MQL5, que servirá como uma ferramenta confiável para tarefas de previsão. Essa biblioteca terá uma arquitetura modular, em que o modelo de PG é dividido em componentes independentes e intercambiáveis, o que permitirá expandir facilmente suas capacidades e garantir simplicidade de manutenção. Ela será desenvolvida levando em conta as seguintes características funcionais principais:

• Flexibilidade na escolha do kernel: possibilidade de conectar facilmente kernels de covariância já existentes, bem como criar suas combinações (SumKernel, ProductKernel) para modelar dependências mais complexas;
• Suporte a diferentes funções de verossimilhança;
• Suporte a diferentes métodos de inferência da distribuição a posteriori para tarefas de classificação e regressão;
• Versatilidade: a biblioteca deve ser universal, permitindo resolver tanto tarefas de regressão quanto de classificação binária;
• Otimização de hiperparâmetros: uso de gradientes analíticos para aumentar a velocidade e a precisão do processo de treinamento. A integração com a biblioteca Alglib deve garantir a otimização eficiente dos hiperparâmetros do modelo.

Vamos analisar com mais detalhes a estrutura da biblioteca. Ela é composta por seis componentes principais, cada um dos quais implementa uma funcionalidade específica:

• Classe GaussianProcess: é o núcleo central da biblioteca, que gerencia todo o ciclo de vida do modelo de PG, desde a inicialização e a otimização dos hiperparâmetros até a execução de predições em novos dados.

• Classe GPOptimizationObjective: essa classe auxiliar funciona como uma "ponte" entre nossa biblioteca e a biblioteca de otimização Alglib. Ela adapta a função objetivo e seu gradiente ao formato exigido pelo Alglib, por meio de herança de CNDimensional_Grad.

• Interface IKernel: define um conjunto de métodos para diferentes funções de covariância, ou kernels. Inclui implementações como RBFKernel, LinearKernel, PeriodicKernel, bem como suas combinações (SumKernel, ProductKernel).

• Interface ILikelihood: define um conjunto de métodos para funções de verossimilhança. As implementações incluem GaussianLikelihood para regressão e LogitLikelihood para classificação binária.

• Interface IInference: fornece métodos para inferência da distribuição a posteriori da função latente do PG. No momento, estão implementados ExactInference e LaplaceInference.

• Estruturas auxiliares e utilitários (StructUtils.mqh): conjunto de enumerações gerais, estruturas de dados, para resultados de inferência e predições, e funções necessárias para o trabalho com dados, matrizes e gráficos para visualização dos resultados.

Graças a essa estrutura modular e a interfaces rigorosamente definidas, poderemos adicionar facilmente novos kernels, métodos de inferência e funções de verossimilhança, o que permitirá evoluir a biblioteca com facilidade no futuro. 

Classe GaussianProcess

A classe GaussianProcess é a classe central da biblioteca. Ela encapsula toda a lógica necessária para construir, treinar e gerar previsões com o modelo de PG. Desenvolvida segundo o princípio de composição, GaussianProcess não contém diretamente a funcionalidade de kernels, verossimilhança ou inferência. Em vez disso, ela reúne esses componentes usando três interfaces principais:

• kernel (IKernel),
• função de verossimilhança (ILikelihood),
• método de inferência (IInference).

Graças a isso, é possível adaptar com flexibilidade o modelo de PG a diferentes tarefas de previsão sem alterar a classe principal GaussianProcess.

//+------------------------------------------------------------------+
//| 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 analisar os principais métodos da classe:

Para criar uma instância da classe, há duas formas principais:

• Método Create: use este método para criar com segurança um objeto GaussianProcess. Esse método realiza as verificações necessárias dos dados de entrada (X_train, y_train, ponteiros para as interfaces) e retorna um ponteiro para o objeto ou NULL em caso de erro.

//+------------------------------------------------------------------+
//| 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;
} 

• Construtor da classe: fornece uma forma direta de inicialização, sem qualquer verificação dos dados. Se você estiver seguro quanto aos seus dados, poderá criar o objeto por meio do construtor.

//+------------------------------------------------------------------+
//| 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){ } 

• O método Fit() inicia o processo completo de treinamento do modelo. No âmbito desse método, ocorre a otimização dos hiperparâmetros do kernel e da função de verossimilhança com o uso do otimizador MinBleic, que minimiza o logaritmo negativo da verossimilhança marginal (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 do método Fit(), preparamos tudo o que é necessário para o funcionamento eficiente do otimizador.

É criado um objeto especial objective_func (GPOptimizationObjective), que representa a função objetivo NLML e seu gradiente analítico em um formato compreensível para o Alglib. Um ponteiro para o objeto GaussianProcess atual é passado ao seu construtor, o que é necessário para chamar o método CalculateNLMLObjective.

Em seguida, obtemos os valores atuais de todos os hiperparâmetros do modelo em um array de hiperparâmetros theta. Esses valores, obtidos do kernel e da função de verossimilhança, servirão como ponto inicial para a busca do ótimo. Para cada hiperparâmetro, são definidos escalas (s), limites inferiores (bndl) e superiores (bndu). Os limites impedem a busca de soluções em regiões incorretas ou sem sentido, por exemplo, comprimentos de escala ou variâncias negativas. O escalonamento é usado pelo otimizador para normalizar os parâmetros, o que aumenta a estabilidade e a velocidade de convergência, especialmente quando as ordens de grandeza dos parâmetros diferem muito. Por padrão, s = 1.0

Em seguida, declaramos um array de ponteiros kernels_to_process para a interface IKernel. Ele será usado para armazenar a lista de todos os kernels básicos cujos hiperparâmetros precisam ser otimizados. Se tivermos um kernel simples, não composto, esse array conterá apenas um elemento, um ponteiro para esse kernel. Se, por outro lado, for um SumKernel ou ProductKernel, então ele armazenará ponteiros para todos os kernels que fazem parte dessa composição.

Em seguida, com o operador dynamic_cast, verificamos se o kernel atual m_kernel, que é um campo da classe GaussianProcess e aponta para o kernel selecionado pelo usuário, é uma instância de SumKernel ou ProductKernel. Se for, ocorre a conversão de tipo para SumKernel ou ProductKernel, e o método GetKernels() é chamado para preencher o array kernels_to_process com todos os kernels que compõem a soma ou o produto de kernels. Se o kernel não for nem uma soma nem um produto, isto é, se for um kernel comum, por exemplo, RBFKernel, então o array kernels_to_process simplesmente conterá um ponteiro para o próprio m_kernel.

Depois disso, percorremos em um loop cada kernel básico encontrado em kernels_to_process e definimos, para seus hiperparâmetros, a escala s e os limites bndl e bndu.

Por fim, depois de todos os hiperparâmetros do kernel, são processados os hiperparâmetros da função de verossimilhança. A verossimilhança gaussiana tem um parâmetro, enquanto a verossimilhança logística não tem parâmetros. Depois que todos os parâmetros são preparados, o processo de otimização é iniciado.

• O método CalculateNLMLObjective() atua como elo de ligação entre a classe principal GaussianProcess e o otimizador externo Alglib. Esta é exatamente a função objetivo que o otimizador MinBleic chama continuamente, por meio da classe GPOptimizationObjective, para avaliar os valores atuais dos hiperparâmetros. Sua principal tarefa é retornar o valor de NLML para um determinado conjunto 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;
}

Em cada iteração, o otimizador MinBleic propõe um novo conjunto de hiperparâmetros. Antes de tudo, CalculateNLMLObjective() recebe esse conjunto (hyperparameters) e usa o método SetHyperparameters() para atualizar os parâmetros correspondentes dentro dos objetos do kernel (m_kernel) e da função de verossimilhança (m_likelihood). Isso é muito importante, pois todos os cálculos posteriores de NLML devem se basear nesses valores atualizados dos hiperparâmetros.

Depois que os hiperparâmetros são atualizados, o método chama Infer() no objeto de inferência (m_inference). Esse é o passo principal, no qual ocorrem todos os cálculos matemáticos complexos voltados para a estimação da distribuição a posteriori.

Os resultados da inferência, incluindo o valor de NLML e seus gradientes, que serão usados pela função Grad, são armazenados no campo privado da classe m_last_inference_result.

Se a inferência for bem-sucedida, o método retorna o valor de NLML.

• O método GaussianProcess::SetHyperparameters(const vector &params) é responsável por distribuir e definir os valores otimizados dos hiperparâmetros do kernel e da função de verossimilhança. 

//+------------------------------------------------------------------+    
//| 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);
}

O vetor params contém todos os hiperparâmetros do modelo de PG em uma ordem fixa: primeiro vêm todos os hiperparâmetros do kernel, ou dos kernels, se for um kernel composto, e depois os parâmetros da função de verossimilhança. A principal característica desse método é o uso de polimorfismo. A mesma chamada m_kernel.SetHyperparameters() se comporta de maneira diferente dependendo do tipo real do objeto para o qual m_kernel aponta durante a execução do programa. 

• Método Predict(). É para isso, em essência, que o modelo é construído: obter previsões em novos dados. 

//+------------------------------------------------------------------+
//| 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;    
    }

Os resultados das previsões, média, variância, probabilidades e rótulos de classe, são gravados na estrutura GPPredictionResult.

Em primeiro lugar, calculamos as matrizes K e K. Essas matrizes são a base para as previsões em PG. Elas são necessárias para calcular a média e a variância da função latente f nos novos pontos de teste. A lógica aqui depende de qual método de inferência, ExactInference ou LaplaceInference, foi usado durante o treinamento, pois eles fornecem componentes diferentes para as fórmulas de predição (Algoritmo 2.1 para Exact, Algoritmo 3.2 para Laplace, do livro GPML de Rasmussen e Williams).

Se foi usada a inferência exata (ExactInference), são extraídos os valores previamente calculados de L_K_noisy e alpha. Se foi usada a aproximação de Laplace (LaplaceInference), são extraídos W, L_B, sW e f_hat, a moda. Em ambos os casos, o resultado são o valor médio (mu_f_star) e a matriz de covariância (Sigma_f_star) da função latente para cada ponto de teste.

Como já discutimos na parte teórica do artigo, existe o problema de calcular a integral para obter a probabilidade da classe. Por isso, são usadas aproximações:

• predict_mode == PROBIT (aproximação probit):

Esta é uma aproximação rápida frequentemente usada. Ela substitui a função sigmoide pela função de distribuição acumulada da distribuição normal, que tem forma semelhante. Isso permite calcular a integral analiticamente.

• predict_mode == NUM_INTEGR (integração numérica):

Nesse modo, é chamada a função auxiliar CalculateNumericalProbability. Ela aproxima numericamente a integral, dividindo o intervalo de f* em intervalos discretos e somando os valores. Isso pode ser mais preciso, mas também mais lento. 

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

Este é um método estocástico. É gerado um grande número de samples, amostras, aleatórios f a partir da distribuição a posteriori q(f∣X, y, x). Para cada sample f, calcula-se sigma(f*).

A média aritmética de todos esses valores sigma(f) é uma aproximação da probabilidade desejada p(y=+1|X, y, x*). Este é o método computacionalmente mais custoso. Para gerar amostras da distribuição normal, é usada a função da biblioteca padrão MathRandomNormal.

Com base nas probabilidades calculadas, para cada um dos métodos de aproximação acima é tomada a decisão sobre o rótulo de classe previsto. Se a probabilidade de pertencer à classe +1 for maior ou igual a 0.5, prevê-se +1, caso contrário, -1.

Classe 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 classe é o elo de ligação entre a nossa classe GaussianProcess e a biblioteca externa de otimização Alglib, em particular o otimizador MinBLEIC.

O Alglib exige que a função objetivo que ele otimiza siga uma determinada interface. É exatamente para isso que serve GPOptimizationObjective. Ela herda de CNDimensional_Grad, a classe base do Alglib que define essa interface. Essa classe base fornece métodos virtuais que a classe GPOptimizationObjective deve implementar. Esses métodos permitem que os otimizadores do Alglib trabalhem com qualquer função objetivo, desde que ela forneça tanto o valor da função quanto o seu gradiente. 

O membro privado GaussianProcess* m_gp contém um ponteiro para o nosso objeto GaussianProcess. Isso permite que a classe GPOptimizationObjective chame o método CalculateNLMLObjective para executar os cálculos necessários. 

O método Grad() é a parte mais importante dessa classe. Ele sobrescreve o método virtual de CNDimensional_Grad e será chamado pelo otimizador Alglib em cada iteração. A função Grad() recebe do Alglib o vetor atual de hiperparâmetros w e deve retornar o valor da função objetivo func e o vetor de seus gradientes grad. 

Considerações finais

Vamos resumir os resultados intermediários.

Na primeira parte do artigo, estabelecemos uma base teórica sólida para compreender o modelo de classificação com PG. Analisamos em detalhes os princípios de funcionamento dos PG para classificação binária e o método de aproximação de Laplace. Esse método é criticamente importante, pois torna a tarefa de classificação praticamente viável e computacionalmente eficiente para as necessidades do trading online, em contraste com o método MCMC, exato, porém extremamente custoso.

Depois de entender a construção teórica, passamos à implementação prática, projetando e descrevendo duas classes-chave da nossa biblioteca de PG:

• GaussianProcess: a classe principal, que encapsula toda a lógica de construção, treinamento e predição do modelo de PG,
• GPOptimizationObjective: que atua como intermediária, preparando a função objetivo e seu gradiente no formato exigido pela biblioteca Alglib para a otimização dos hiperparâmetros.

Na segunda parte, concluiremos a implementação da biblioteca, apresentando:

• descrição detalhada e código de implementação das interfaces-chave: IKernel, para diferentes kernels, IInference, para métodos de inferência, e ILikelihood, para funções de verossimilhança;
• exemplos de funcionamento da biblioteca com dados sintéticos, para uma demonstração clara de suas capacidades;
• aplicação prática no trading: desenvolveremos indicadores para classificação e regressão com base em nossa biblioteca, mostrando como os PG podem ser usados na tomada de decisões de trading.