機械学習におけるガウス過程（第1回）：MQL5における分類モデル

はじめに

機械学習モデルであるガウス過程(GP)について、引き続き理解を深めていきます。前回の記事では、連続値を予測することを主な目的とした回帰問題について詳しく検討しました。今回は、さらに複雑なテーマである分類を扱います。分類における最大の難しさは、ガウス過程で分類する場合、推論に閉形式解が存在しないことです。そのため、ラプラス近似などの近似手法を利用する必要があります。

この複雑な問題を効果的に解決するため、今回はMQL5でモジュール型のガウス過程ライブラリを開発します。このアプローチにより、GPモデルを独立したコンポーネントに分割してコードを構造化でき、今後の改良や拡張に向けた強固な基盤を構築できます。このライブラリは、回帰と分類の両方に利用できる汎用的なツールとなることを目指します。

記事の前半では、GP分類の理論について詳しく検討します。そこでは、近似手法の基礎となる数学についても扱います。また、ライブラリの主要クラスであるGaussianProcessについて紹介します。このクラスはモデルを構成するすべてのコンポーネントを統合します。さらに、Alglib最適化ライブラリとの連携を担うGPOptimizationObjectiveクラスについても説明します。

分類

分類とは、機械学習において、ある対象をあらかじめ定義されたカテゴリーのいずれかに割り当てるタスクです。たとえば金融分野では、過去のデータに基づいて株価が上昇するか下降するかを予測するために分類を利用できます。

本記事では、対象が2つのクラスのいずれかに属する二値分類に焦点を当てます。たとえば、「上昇」を+1、「下降」を-1として分類します。サポートベクターマシン(SVM)や決定木などの手法では、基本的にクラスラベルだけが出力されます。一方、GPでは確率的な予測が可能です。たとえば、モデルが「株価が上昇する確率は75%」と判断することができます。このような情報は取引において特に価値があります。予測に対する確信度を利用することで、より適切な判断を下し、信頼性の低いシグナルを除外できるからです。 

残念ながら、GPを利用した分類問題の解決は、回帰よりもはるかに複雑です。その理由は、使用する尤度関数の種類にあります。 

• 回帰では、通常ガウス尤度を使用します。GPを関数の事前分布として考えた場合、GPとガウス尤度を組み合わせることで、事後分布を解析的に求めることができます。そのため、計算が大幅に簡単になります。
• 一方、分類では目的変数が離散的なクラスラベルであるため、ガウス尤度は適していません。その代わりに、たとえばロジット尤度を使用できます。しかし、この場合、事後分布もガウス分布にはならず、閉形式解も存在しません。

その結果、複雑な近似推論手法を利用する必要があります。これらの手法の基本的な考え方は、真の非ガウス事後分布を、その最頻値を中心とするガウス分布で近似することです。この記事では、近似された事後分布を得るための手法として、最も単純かつ効果的なアプローチの1つであるラプラス近似に焦点を当てます。 

二値分類におけるGPベースの予測の基本的な考え方は、非常にシンプルです。まず、潜在関数f(x)の事前分布から始めます。GPが1つの関数だけを生成するのではなく、無限個の可能な関数を生成すると考えてみましょう。それぞれの関数は、データに存在する可能性のある「潜在的な」依存関係を表します。そして、それぞれの潜在関数f(x)の実現値がロジスティック関数（シグモイド関数）に入力されます。シグモイド関数は、f(x)の値である任意の実数を0から1の間の確率へ変換します。これが、+1クラスに属する事前確率π(x)となります。

ここで重要なのは、πがfの決定論的な関数であるという点です。しかし、f自体は確率的なものであり、GPからサンプリングされた値です。そのため、πも確率的な関数になります。この概念は、1次元の入力空間Xについて、図1と図2に明確に示されています。

図1：潜在関数f(x)の実現例

図1は潜在関数の1つの実現（サンプル関数）を示しています。指定されたカーネルのハイパーパラメータに対応した、この関数の典型的な挙動を見ることができます。

図2：シグモイドによって変換された同じ関数

図2では、同じ関数f(x)にロジスティック関数（シグモイド関数）を適用した結果を示しています。

このようにして、π(x)=σ(f(x))というクラス所属確率の事前分布を得ることができます。ただし、この段階ではまだ学習データyは考慮されていません。yの観測がなければ、この事前分布は経験的な裏付けを欠く初期仮定にとどまり、どの仮定が妥当かをモデルは判断できません。 

当然ながら、事前仮定の選択は最終的な事後分布や予測結果に大きな影響を与えます。これはベイズ的アプローチの重要な特徴です。関数の事前分布、ひいては最終的なモデルの特性は、研究者が選択するカーネルの種類に依存するからです。 

推論

さて、根拠のある予測をおこなうためには、実際の学習データyを考慮する必要があります。ここで登場するのが推論です。推論の主な目的は、事前の信念を事後の信念へ変換することです。つまり、観測されたデータを考慮して調整された分布へ変換します。分類の場合、この処理は自然に2つの段階に分けられます。

ステップ1：潜在関数f∗の予測分布

最初のステップでは、新しいテスト点x*について、観測された学習データ(X, y)を条件とした潜在関数f*の事後分布p(f*|X, y, x*)を計算します。これは次の積分によって定義されます。

ここで

• p(f*∣X, x*, f)は、学習点Xにおける潜在関数fが与えられたときの、新しいテスト点x*における潜在関数fの条件付き分布です。GPは定義上、同時正規分布に従うため、この分布は常に正規分布になります。
• p(f|X, y)は、学習データにおける潜在関数fの事後分布です。しかし、非線形な尤度関数（シグモイド関数）のため、これはガウス分布にはなりません。

ここで重要なのは、p(f|X,y)が正規分布ではないため、この積分には閉形式解が存在しないということです。したがって、この計算には近似手法が必要になります。 

ステップ2：最終的な予測確率π*

2番目のステップでは、この予測分布を利用して、最終的な予測確率π*を求めます。これは、テスト点xが正のクラス(y* = +1)に属する確率です。

ここでσ(f*)はロジスティック関数（シグモイド関数）であり、潜在関数f*の値を0から1の間の確率に変換します。この積分は、f*の取り得るすべての値について、それぞれの確率を事後予測分布で重み付けして平均していることを意味します。つまり、この1次元積分は、分布p(f|X,y,x)に関する関数σ(f*)の数学的期待値です。

ここでも、ロジット尤度の場合、この積分には閉形式解が存在しません。そのため、ここでも近似手法が必要になります。先に述べておくと、今回のGPライブラリでは、このような近似を3種類実装します。これにより、精度と計算コストに応じて適切な手法を選択できます。

•  プロビット近似
•  数値積分
•  モンテカルロ法

ここまで説明した2つのステップ、つまり潜在関数の事後分布を計算することと、その後に積分して予測確率を求めることが、GPにおけるベイズ推論の一般的な枠組みです。目的とする予測を得るためには、この2つの積分を計算する必要があり、どちらも近似手法を必要とします。 

ラプラス近似

すでに確認したように、分類におけるベイズ推論には、閉形式で解くことのできない積分が含まれています。ラプラス近似は、非ガウス分布であるp(f|X,y)を、ガウス分布q(f|X,y)で近似することによって、この問題を解決します。条件付き分布p(f*∣X, x*, f)もガウス分布であるため、結果として予測分布p(f*∣X, y, x*)もガウス分布になります。これにより、f*の平均と分散を解析的な式として導出でき、その後の計算が大幅に簡単になります。つまり、ラプラス近似の美しさと計算効率の高さは、事後分布と予測の計算をガウス分布に対する演算へと帰着できることにあります。

ただし、ラプラス近似が一種の妥協であることを理解する必要があります。この手法は、閉形式で扱えない問題を計算可能な問題へ変換しますが、その代償として、真の事後分布の形状を完全には表現できません。この正規近似の品質は、真のp(f|X,y)の分布がどれだけ正規分布に近いかに直接依存します。真の分布が正規分布に近ければ近いほど近似精度は高くなり、逆に離れているほど精度は低下します。

もし近似ではなく、真のp(f*∣X, y, x*)の分布そのものを求めたいのであれば、通常はMCMC（マルコフ連鎖モンテカルロ）法を利用します。MCMC法はより正確な推定値を提供できますが、計算コストが非常に高く、実装も困難です。そのためMCMCは、近似推論手法と比較する際のゴールドスタンダードとして利用できます。

では、ラプラス近似についてもう少し詳しく見てみましょう。この近似は、真の事後分布p(f|X,y)の最頻値（事後分布を最大にする点）を中心として構築されます。具体的には、この最頻値の周辺で事後密度の対数を2次のテイラー展開によって近似します。数学的には、事後密度の対数を次のように近似します。

ここで

• q(f∣X, y)は、事後分布p(f|X,y)に対するガウス近似
• f_hat = argmax(f) p(f|X, y)は、事後分布の最頻値
• A = −∇∇ log p(f|X, y)|f=f_hatは、最頻点における事後分布の負の対数に対するヘッセ行列

まず、ラプラス近似をおこなうためには、潜在関数fの最も確からしい値、つまり最頻値f_hatを求める必要があります。事後分布p(f|X,y)を求めるには、ベイズの定理を利用します。すでに分かっているように、ベイズの定理は事後分布と、尤度p(y|f)、事前分布p(f|X)、周辺尤度p(y|X)を次のように関連付けます。

fについてp(f|X,y)を最大化する場合、正規化定数p(y|X)を知る必要はありません。これはp(y|X)がfに依存せず、したがって最大値の位置に影響を与えないためです。そのため、尤度と事前分布の積p(y|f)p(f|X)に比例する非正規化事後分布を使って計算できます。

さらに、非常に小さな確率値を扱うことによる数値的な問題を避けるため、この非正規化事後分布の対数を取ります。対数の性質により、確率の積はそれぞれの対数の和になります。

ここで得られるΨ(f)が、潜在関数の最頻値を求めるためにニュートン法で最大化する目的関数となります。ニュートン法では、Ψ(f)をfについて1階および2階微分する必要があります。

この式をfについて微分すると、次の式が得られます。

ここで

• W = −∇∇ log p(y|f)は、対数尤度の負のヘッセ行列です。この行列は対角行列になります。

勾配とヘッセ行列を計算したら、ニュートン法を使用して反復的に最頻値を求めます。

各反復において、ニュートン法は勾配とヘッセ行列によって決定される方向へ現在の最頻値の推定値を更新します。そして、収束するまでこの処理を繰り返します。 

最頻値が求まったら、近似されたガウス分布の共分散行列を計算できます。この行列は、最頻点f_hatで評価したΨ(f)の負のヘッセ行列の逆行列に等しくなります。

したがって、ガウス近似における共分散行列Σは次のようになります。

これで、ラプラス近似の第1段階、つまり事後分布の正規近似を求める処理が完了します。

ラプラス近似による予測

q(f|X,y)が得られたら、推論の第2段階である新しいテスト点x*に対する予測へ進むことができます。ここでは、予測分布p(f∗∣X, y, x∗)を求めます。ラプラス近似によってp(f|X,y)がガウス分布、すなわちq(f|X,y)として表現され、さらにp(f∗∣X, x∗, f)もガウス分布であるため、結果としてp(f∗∣X, y, x∗)もガウス分布になります。これにより、その事後平均と分散を解析的に求めることができます。

新しいテスト点x*における潜在関数f*の平均値(mu_f_star)は次のように計算されます。

新しいテスト点xにおける潜在関数fの分散Var(f*)(Sigma_f_star)は、次の式で計算されます。

予測分布の平均と分散が得られたので、最後に目的とするクラス所属確率π*を計算できます。 

この式が、ラプラス近似に基づくGP分類器における確率予測の核心となります。

ここで、クラス確率を単純にσ(E[f∗])として計算していないことに気付いたかもしれません。つまり、f*の事後平均をそのままシグモイド関数に入力する方法です。この方法はMAP予測（最大事後確率予測）と呼ばれ、もちろん有効な手法です。

しかし、MAP予測σ(E[f*])を計算すると、f*に存在する不確実性を無視することになります。単純にf*の中心的な推定値、つまり平均だけを取り出し、それを確率に変換しているからです。一方、E[σ(f*)]を計算する場合、これは積分に対応しますが、f*の分布全体の形状を考慮します。そのため、特にfの不確実性が大きい場合、つまりV[f*]の分散が大きい場合や、f*の分布が非対称である場合には、より正確で意味のある予測確率を得ることができます。この方法は平均予測確率と呼ばれます。

この違いを理解することには、重要な実践的意味があります。

• 目的が単純な二値クラスラベル、たとえば「買い」または「売り」、「+1」または「-1」を得ることだけであれば、より単純なMAP予測で十分な場合があります。これは、計算コストの高い平均予測と同じクラスラベルが得られるためです。
• しかし、確率そのものが重要である場合には、平均予測確率E[σ(f*)]の方が正確です。これはモデルの不確実性を完全に考慮するためです。 

取引では、単純な二値クラスラベル、「買い」または「売り」だけでは十分ではありません。私たちには、確率によって表現される細かな確信度が必要です。確率値を用いることで、取引シグナルをフィルタリングすることができます。成功確率が0.51のシグナルは、ランダムな推測をわずかに上回る程度であるため、成功確率が0.60のシグナルと比べて価値は大きく異なります。この確率を利用することで、取引を開始するための閾値を設定できます。たとえば、成功確率が0.55または0.60を超えた場合のみ取引を開始する、といった条件を設定できます。これによって、偽シグナルの数を減らすことができます。

ラプラス近似における周辺尤度

GPによる分類の推論メカニズムを理解したところで、次に「モデルをどのように調整すれば最適な予測が得られるのか」という問題が生じます。その答えとなるのが周辺尤度(LML)です。これは、モデルのハイパーパラメータθを最適化するための目的関数です。これを計算しなければ、データを最も適切に説明するパラメータを見つけることはできません。 

ここで、B

最適化する目的関数を定義したら、次に重要となるのが、ハイパーパラメータθに対する偏微分の計算です。これは、解析的勾配を使用して最適化するために必要です。解析的勾配を利用することで、数値微分を利用する場合と比較して計算を数倍高速化できます。また、解析的勾配によって、オプティマイザはNLML目的関数の最小値に向かって、より効率的かつ正確に移動できます。 

LMLの勾配は、明示的な部分と暗黙的な部分から構成されます。

明示的な部分は次の式で計算します。

ここでの主な問題は、各ハイパーパラメータについてカーネル行列Kの微分を計算することです。選択したカーネル関数に対する微分処理の実装については、記事の第2部で扱います。

暗黙的な部分は2つの因子から構成されます。暗黙的な部分の第1因子は、次の式を使って求めます。

この式を計算するためには、尤度の対数の3階微分を計算する必要があります。

暗黙的な部分のもう一方の因子は、次のように計算します。 

最後に、NLMLはハイパーパラメータの推定だけでなく、異なるモデルの比較にも利用できることを指摘しておきます。たとえば、異なる種類のカーネルを使用したモデルを比較できます。NLMLの値が低いモデルほど優れていると考えられます。これは、より高い周辺尤度を意味し、モデルが観測されたデータをより適切に説明していることを示すためです。

さらに、GPでは周辺尤度を利用してハイパーパラメータを最適化することにより、データへの適合度とモデルの複雑さのトレードオフを自動的に処理できます。NLMLは複雑すぎるモデルに自然なペナルティを与えるため、過学習を抑制できます。このため、ニューラルネットワークの学習などでおこなうような、過学習を防ぐための明示的な停止条件を必ずしも設定する必要がありません。NLML自体の最適化によって、最適なバランスが探索されるためです。これは、ガウス過程におけるベイズ的アプローチの主要な利点の1つです。

ガウス過程ライブラリ

ここまで必要な理論的概念をすべて確認したので、次に実装へ進みます。今回の主な目標は、予測タスクにおいて信頼性の高いツールとして利用できる汎用的なGPライブラリをMQL5で構築することです。このライブラリはモジュール型アーキテクチャを採用します。GPモデルを独立して差し替え可能なコンポーネントに分割することで、機能を容易に拡張でき、保守性も高めることができます。ライブラリは、以下の主要な機能を考慮して開発します。

• カーネル選択の柔軟性：既存の共分散カーネルを容易に接続できるだけでなく、SumKernelやProductKernelによって複数のカーネルを組み合わせ、より複雑な依存関係をモデル化できること。
• さまざまな尤度関数への対応。
• 分類および回帰問題における、さまざまな事後分布推論手法への対応。
• 汎用性：回帰問題と二値分類問題の両方を解決できる汎用的なライブラリであること。
• ハイパーパラメータ最適化：解析的勾配を利用することで、学習処理の速度と精度を向上させること。Alglibライブラリとの統合により、モデルのハイパーパラメータを効率的に最適化できるようにします。

それでは、ライブラリの構造を詳しく見ていきましょう。ライブラリは、主に6つのコンポーネントから構成され、それぞれが特定の機能を実装します。

• GaussianProcessクラスはライブラリの中心となるハブです。モデルの初期化からハイパーパラメータの最適化、新しいデータに対する予測の実行まで、GPモデルのライフサイクル全体を管理します。

• GPOptimizationObjectiveクラスは、私たちのライブラリとAlglib最適化ライブラリの間をつなぐ補助クラスです。このクラスは、目的関数とその勾配をAlglibが要求する形式へ変換します。これはCNDimensional_Gradを継承することによって実現されます。

• IKernelインターフェース：さまざまな共分散関数（カーネル）のためのメソッド群を定義します。ここには、RBFKernel、LinearKernel、PeriodicKernel、およびそれらの組み合わせ（SumKernel、ProductKernel）などの実装が含まれます。

• ILikelihoodインターフェース：尤度関数のためのメソッド群を定義します。実装としては、回帰用のGaussianLikelihoodと、二値分類用のLogitLikelihoodがあります。

• IInferenceインターフェース：GPの潜在関数について事後分布を推論するためのメソッドを提供します。現在、ExactInferenceとLaplaceInferenceの2つの実装があります。

• 補助構造物とユーティリティ(StructUtils.mqh)：共通の列挙型、データ構造、ユーティリティ関数がまとめられています。ここには、推論結果や予測結果を格納するデータ構造に加え、データ、行列、そして結果を可視化するためのグラフを扱うために必要な関数が含まれています。

このように、明確に定義されたインターフェースとモジュール型の構造を採用することで、新しいカーネル、推論手法、尤度関数を容易に追加できます。その結果、今後のライブラリ開発を柔軟に進めることができ、機能の拡張と保守も容易になります。 

GaussianProcessクラス

GaussianProcessクラスは、ライブラリの中心となるクラスです。GPモデルの構築、学習、予測に必要なすべてのロジックをカプセル化しています。コンポジションの原則に基づいて設計されており、GaussianProcess自身がカーネル、尤度、推論の機能を直接保持することはありません。代わりに、次の3つの主要なインターフェースを介して、これらのコンポーネントを統合します。

• カーネル(IKernel)
• 尤度関数(ILikelihood)
• 推論方法(IInference)

これにより、基盤となるGaussianProcessクラスを変更することなく、GPモデルをさまざまな予測タスクに柔軟に適応させることができます。

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

クラスの主なメソッドを見ていきましょう。

クラスのインスタンスを生成する方法には、主に2つあります。

• Createメソッド：GaussianProcessオブジェクトを安全に生成するために使用します。このメソッドでは、入力データ（X_train、y_train、各種インターフェースへのポインタ）に対して必要なチェックをおこない、正常に生成できた場合はオブジェクトへのポインタを返し、エラーが発生した場合はNULLを返します。

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

• クラスコンストラクタ：データのチェックをおこなわずに、直接初期化する方法を提供します。データが正しいことを確信している場合は、コンストラクタを使用してオブジェクトを生成できます。

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

• Fit()メソッド：モデル全体の学習処理を開始します。このメソッドでは、MinBleicオプティマイザを使用してカーネルのハイパーパラメータと尤度関数を最適化し、負の対数周辺尤度(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;     
} 

Fit()メソッドの内部では、オプティマイザが効率的に動作するために必要なすべての準備をおこないます。

まず、Alglibが理解できる形式でNLMLの目的関数とその解析的勾配を表現する特殊なオブジェクトobjective_func(GPOptimizationObjective)を生成します。そのコンストラクタには、現在のGaussianProcessオブジェクトへのポインタを渡します。これは、CalculateNLMLObjectiveメソッドを呼び出すために必要です。

次に、モデルのすべてのハイパーパラメータの現在値をthetaハイパーパラメータ配列に取得します。これらの値はカーネルと尤度関数から取得され、最適値を探索する際の初期値として使用されます。各ハイパーパラメータについて、スケール(s)、下限(bndl)、上限(bndu)を指定します。境界値を設定することで、問題のある領域や意味を持たない領域、たとえば負のスケール長や分散などを探索することを防ぎます。スケーリングは、オプティマイザがパラメータを正規化するために使用するもので、特にパラメータの桁が大きく異なる場合に、計算の安定性と収束速度を向上させます。デフォルト値はs=1.0です。

次に、IKernelインターフェースへのポインタを格納するkernels_to_process配列を宣言します。この配列は、ハイパーパラメータを最適化する必要があるすべての基本カーネルのリストを格納するために使用します。単純なカーネル（複合カーネルではない）の場合、この配列にはそのカーネルへのポインタが1つだけ格納されます。一方、SumKernelまたはProductKernelの場合は、その合成を構成するすべてのカーネルへのポインタが格納されます。

続いて、dynamic_cast演算子を使用して、現在のカーネルm_kernel（GaussianProcessクラスのフィールドであり、ユーザーが選択したカーネルを指しています）がSumKernelまたはProductKernelのインスタンスであるかどうかを確認します。該当する場合は、SumKernelまたはProductKernelへの型キャストをおこない、GetKernels()メソッドを呼び出します。このメソッドによって、和または積を構成するすべてのカーネルがkernels_to_process配列に格納されます。カーネルがSumKernelでもProductKernelでもない場合、つまりRBFKernelなどの通常のカーネルの場合は、kernels_to_process配列にm_kernel自身へのポインタを格納します。

その後、kernels_to_processに格納された各基本カーネルを順番に処理し、そのハイパーパラメータにスケールs、下限bndl、上限bnduを設定します。

最後に、すべてのカーネルのハイパーパラメータの処理が完了したら、尤度関数のハイパーパラメータを処理します。ガウス尤度には1つのパラメータがありますが、ロジット尤度にはパラメータがありません。すべてのパラメータの準備が完了すると、最適化処理が開始されます。

• CalculateNLMLObjective()メソッド：メインのGaussianProcessクラスと外部のAlglibオプティマイザをつなぐ役割を果たします。これは、MinBleicオプティマイザが現在のハイパーパラメータの値を評価するために、GPOptimizationObjectiveクラスを介して繰り返し呼び出す目的関数そのものです。その主な役割は、与えられたハイパーパラメータのセットに対するNLMLの値を返すことです。

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

各反復において、MinBleicオプティマイザは新しいハイパーパラメータのセットを提案します。CalculateNLMLObjective()はまず、このハイパーパラメータのセットを受け取り、SetHyperparameters()メソッドを使用して、カーネル(m_kernel)および尤度関数(m_likelihood)の各オブジェクト内部に対応するパラメータを更新します。これは非常に重要です。というのも、その後のすべてのNLMLの計算は、これらの最新のハイパーパラメータの値に基づいておこなわれる必要があるためです。

ハイパーパラメータの更新後、推論オブジェクト(m_inference)のInfer()メソッドを呼び出します。ここが主要な処理段階であり、事後分布の推定に必要となる複雑な数学的計算がすべて実行されます。

推論結果には、NLMLの値とその勾配が含まれており、これらはGrad関数で使用されます。これらの結果は、クラスのprivateフィールドであるm_last_inference_resultに保存されます。

推論が正常に完了すると、メソッドはNLMLを返します。

• GaussianProcess::SetHyperparameters(const vector&params)メソッド：カーネルのハイパーパラメータと尤度関数の最適化された値を、それぞれ対応するコンポーネントに振り分けて設定する役割を担います。 

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

paramsベクトルには、GPモデルのすべてのハイパーパラメータが一定の順序で格納されています。最初にカーネルのすべてのハイパーパラメータが続きます。複合カーネルの場合は、構成されている各カーネルのハイパーパラメータが含まれます。その後に、尤度関数のパラメータが続きます。このメソッドの重要な特徴は、ポリモーフィズムを利用している点です。同じm_kernel.SetHyperparameters()の呼び出しであっても、実行時にm_kernelが指しているオブジェクトの実際の型に応じて、異なる動作をします。 

• Predict()メソッド：これはまさに、モデルが新しいデータに基づいて予測をおこなうために構築されているものです。 

//+------------------------------------------------------------------+
//| 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 labels (+1 or -1)
            result.predicted_labels[i] = (result.predicted_probabilities[i] >= 0.5) ? 1.0 : -1.0;
        }
    }
    return true;    
    }

予測結果（平均、分散、確率、クラスラベル）は、GPPredictionResult構造体に格納されます。

まず、K*行列とK**行列を計算します。これらの行列は、GPにおける予測の基礎となるものです。新しいテスト点における潜在関数f∗の平均と分散を計算するために使用されます。ここでの処理は、学習時にどの推論手法（ExactInferenceまたはLaplaceInference）を使用したかによって異なります。これは、それぞれの手法が予測式に必要となる異なる要素を提供するためです。ExactInferenceについてはRasmussenとWilliamsのGPMLのアルゴリズム2.1、LaplaceInferenceについてはアルゴリズム3.2に対応しています。

ExactInferenceを使用した場合は、事前に計算されたL_K_noisyとalphaを取得します。LaplaceInferenceを使用した場合は、W、L_B、sW、およびf_hat（最頻値）を取得します。どちらの場合も、各テスト点における潜在関数の平均(mu_f_star)と共分散行列(Sigma_f_star)が得られます。

記事の理論部分ですでに説明したように、クラス確率を求めるためには積分を計算する必要があります。しかし、この積分には解析解が存在しないため、近似手法を使用します。

• predict_mode == PROBIT（プロビット近似）：

これは、よく使用される高速な近似手法です。シグモイド関数を、形状が類似している正規分布の累積分布関数で置き換えます。これにより、積分を解析的に計算できるようになります。

• predict_mode == NUM_INTEGR（数値積分）：

このモードでは、補助関数CalculateNumericalProbabilityが呼び出されます。この関数は、f*の範囲を離散的な区間に分割し、それぞれの値を加算することで積分を数値的に近似します。より高い精度が期待できる一方で、計算速度は低下します。 

• predict_mode == MONTE_CARLO（モンテカルロ法）：

これは確率的な手法です。まず、事後分布q(f*∣X, y, x*)から多数のランダムなf*のサンプルを生成します。それぞれのサンプルf*について、σ(f*)を計算します。

これらのσ(f*)の算術平均が、求める確率p(y*=+1|X, y, x*)の近似値となります。この方法は、3つの手法の中で最も計算コストが高くなります。正規分布からサンプルを生成するために、標準ライブラリ関数MathRandomNormalを使用します。

計算された確率に基づいて、上記の各近似手法について予測クラスラベルを決定します。クラス+1に属する確率が0.5以上の場合は+1を予測し、それ以外の場合は-1を予測します。

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]);
        }       
    }      
};

このクラスは、私たちのGaussianProcessクラスと外部のAlglib最適化ライブラリ、具体的にはMinBLEICオプティマイザとの橋渡しをします。

Alglibでは、最適化対象となる目的関数が、特定のインターフェースに従っている必要があります。まさにそのために用意されているのがGPOptimizationObjectiveです。このクラスは、このインターフェースを定義するAlglibの基底クラスCNDimensional_Gradを継承しています。この基底クラスには、GPOptimizationObjectiveクラスが実装すべき仮想メソッドが定義されています。これらのメソッドによって、目的関数の値とその勾配の両方を提供できれば、Alglibのオプティマイザは任意の目的関数を扱えるようになります。 

privateメンバGaussianProcess* m_gpには、私たちのGaussianProcessオブジェクトへのポインタが格納されます。これにより、GPOptimizationObjectiveクラスからCalculateNLMLObjectiveメソッドを呼び出し、必要な計算を実行できます。 

Grad()メソッドは、このクラスで最も重要な部分です。このメソッドはCNDimensional_Gradから継承した仮想メソッドをオーバーライドし、Alglibのオプティマイザによって各反復で呼び出されます。Grad()関数は、Alglibから現在のハイパーパラメータベクトルwを受け取り、目的関数funcの値と、その勾配を格納したgradベクトルを返します。 

結論

ここまでの中間的な成果をまとめてみましょう。

記事の前半では、GP分類モデルを理解するための確かな理論的基盤を構築しました。GPによる二値分類の仕組みと、ラプラス近似法について詳しく検討しました。この手法は、分類問題をオンライン取引の用途において実用的かつ計算効率の高いものにするうえで極めて重要です。精度の高いMCMC法は非常に大きな計算コストを必要とするのに対し、ラプラス近似は十分な精度を維持しながら、より現実的な計算量で処理できます。

理論的な構成を確認した後、実装へと進み、GPライブラリの主要な2つのクラスを設計・実装しました。

• GaussianProcess：GPモデルの構築、学習、予測に必要なすべてのロジックをカプセル化する主要クラス。
• GPOptimizationObjective：仲介役として機能し、ハイパーパラメータ最適化のために、目的関数とその勾配をAlglibライブラリが要求する形式に整えます。

第2部では、以下の内容を実装してライブラリを完成させます。

• 各種カーネルに対応するIKernel、推論手法に対応するIInference、尤度関数に対応するILikelihoodという主要インターフェースの詳細な説明と実装コード。
• 合成データを使用したライブラリの動作例。これにより、その機能を具体的に確認します。
• 取引への実践的な応用。私たちのライブラリを基盤として分類および回帰用のインジケータを開発し、GPを取引判断にどのように活用できるかを実際に示します。