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.

//+------------------------------------------------------------------+
//| Класс Гауссовского Процесса                                      |
//+------------------------------------------------------------------+
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 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                                                     |
//+------------------------------------------------------------------+
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;
}

• 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                               |
//+------------------------------------------------------------------+
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).

//+------------------------------------------------------------------+
//| Метод для обучения модели                                        |
//+------------------------------------------------------------------+
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 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.

//+-------------------------------------------------------------------+
//| Метод, который будет вызываться оптимизатором для вычисления 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;
}

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. 

//+------------------------------------------------------------------+    
//| Метод для установки гиперпараметров                              |
//+------------------------------------------------------------------+
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);
}

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. 

//+------------------------------------------------------------------+
//| Метод предсказания для регрессии и классификации                 |
//+------------------------------------------------------------------+    
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;    
    }

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

//+------------------------------------------------------------------+
//| Класс для целевой функции оптимизатора 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 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.