Gaußsche Prozesse für maschinellen Lernen (Teil 2): Implementierung und Testen eines Klassifizierungsmodells in MQL5
17
4
Einführung
Im vorherigen Artikel haben wir die theoretischen Grundlagen von Gauß-Prozessen als bayesschem Machine-Learning-Modell kennengelernt und mit der Erstellung einer GP-Bibliothek in MQL5 begonnen, wobei wir zwei Schlüsselklassen beschrieben haben: GaussianProcess und GPOptimizationObjective.
Hier werden wir die Bibliothek vervollständigen, indem wir uns die Implementierung der wichtigsten Schnittstellen im Detail ansehen: IKernel, ILikelihood und IInference. Danach werden wir die Bibliothek mit synthetischen Daten testen und Indikatoren für Klassifizierung und Regression schreiben und ihren Einsatz im Online-Modus – mit einem erneuten Training des Modells bei jeder neuen Bar – testen.
IKernel-Schnittstelle
Die IKernel-Schnittstelle, die Sie in der Datei Kernels.mqh finden, ist die Grundlage für die Implementierung von Kovarianz-Kernels in unserer Bibliothek. Sie macht das System flexibel und leicht erweiterbar: Sie können neue Arten von Kernels oder deren Kombinationen hinzufügen, ohne die Hauptstruktur des Codes zu ändern.
interface IKernel { // Calculate the covariance matrix between two data sets virtual matrix Compute(const matrix &X1, const matrix &X2) = 0; // Calculate the derivative of the covariance matrix with respect to the given hyperparameter // param_index: index of the hyperparameter by which the derivative is taken (starting from 0) virtual matrix ComputeDerivative(int param_index) = 0; // Return the current values of all kernel hyperparameters virtual vector GetHyperparameters() const = 0; // Set new values for kernel hyperparameters virtual void SetHyperparameters(const vector ¶ms) = 0; // Return the number of kernel hyperparameters virtual int GetNumHyperparameters() const = 0; // Return the kernel string name virtual string GetName() const = 0; };
Die Schnittstelle definiert zwei der wichtigsten Methoden, die der Funktionsweise jedes Kernels zugrunde liegen:
- Compute(const matrix &X1, const matrix &X2) – Berechnung der K-Kovarianzmatrix zwischen zwei Datensätzen.
- ComputeDerivative(int param_index) – Berechnung der Ableitung der Kovarianzmatrix in Bezug auf einen der Kernel-Hyperparameter. param_index gibt den Index des Hyperparameters an, nach dem die Ableitung berechnet wird.
RBFKernel-Klasse (Radialer Basiskernel)
RBFKernel, auch bekannt als gaußscher oder quadratischer exponentieller Kernel, ist einer der am häufigsten verwendeten Kovarianz-Kernels. Er zeichnet sich durch zwei Hyperparameter aus: die Signalvarianz σ_f (Amplitude) und die Längenskala l, die die Glattheit der Funktion bestimmt.
//+------------------------------------------------------------------+ //| Kernel RBF class | //+------------------------------------------------------------------+ class RBFKernel : public IKernel { private: double length; double sigma_f; int n; // Number of rows in X1; int m; // Number of rows in X2; matrix m_K; matrix m_D_sq; // matrix of squared distances ||x_i - x_j||^2 public: //Constructor RBFKernel(double ls, double sf) : length(ls), sigma_f(sf) {} //Destructor ~RBFKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); matrix XX(n, m); if (!XX.GeMM(X1, X2.Transpose(), 1.0, 0.0)) { Print("Error: Failed to calculate XX = X * X^T"); return matrix::Zeros(n, m); } matrix diag1(n, 1); matrix X1_sq = X1 @ X1.Transpose(); diag1.Col(X1_sq.Diag(), 0); matrix diag2(m, 1); matrix X2_sq = X2 @ X2.Transpose(); diag2.Col(X2_sq.Diag(), 0); m_D_sq = matrix::Ones(n, 1) @ diag2.Transpose() + diag1 @ matrix::Ones(1, m) - 2 * XX; m_K = sigma_f * sigma_f * MathExp((-1 * m_D_sq) / (2 * length * length)); return m_K; } matrix ComputeDerivative(int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { dK = (2.0 / sigma_f) * m_K; break; } case 1: { dK = m_K * (m_D_sq / (length * length * length)); break; } } return dK; } vector GetHyperparameters() const override { vector params(2); params[0] = sigma_f; params[1] = length; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 2) { sigma_f = params[0]; length = params[1]; } } int GetNumHyperparameters() const override { return 2; } string GetName() const override { return "RBFKernel"; } };
Glücklicherweise sind die Ableitungen für den RBF-Kernel leicht zu finden:
-
Ableitung nach sigma_f:
- Ableitung nach l:
Für den Parameter length führen wir eine elementweise Multiplikation der K-Matrix mit (D_sq / length^3) durch, wobei D_sq die Matrix der quadrierten euklidischen Abstände ist.
LinearKernel-Klasse (Linearer Kernel)
LinearKernel ist ein einfacher Kovarianz-Kernel, der eine lineare Beziehung zwischen Datenpunkten annimmt. Er wird häufig verwendet, wenn erwartet wird, dass die Funktion durch ein lineares Modell angenähert werden kann, oder als Komponente in komplexeren zusammengesetzten Kernels.
Der lineare Kernel hat einen Hyperparameter: die Signalvarianz σ_l.
//+------------------------------------------------------------------+ //| Linear kernel class | //+------------------------------------------------------------------+ class LinearKernel : public IKernel { private: double sigma_l; int n; int m; matrix m_K; public: //Constructor LinearKernel(double sl): sigma_l(sl){} //Destructor ~LinearKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); m_K.GeMM(X1, X2.Transpose(), sigma_l * sigma_l, 0.0); return m_K; } matrix ComputeDerivative(int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { // Derivative with respect to sigma_l dK = (2.0 / sigma_l) * m_K; break; } } return dK; } vector GetHyperparameters() const override { vector params(1); params[0] = sigma_l; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 1) { sigma_l = params[0]; } } int GetNumHyperparameters() const override { return 1; } string GetName() const override { return "LinearKernel"; } };
ComputeDerivative(int param_index)-Methode:
PeriodicKernel-Klasse (Periodischer Kernel)
PeriodicKernel wird verwendet, um Daten mit wiederkehrenden Mustern oder Saisonalität zu modellieren, wodurch zyklische Abhängigkeiten erfasst werden können.
//+------------------------------------------------------------------+ //| Periodic kernel class | //+------------------------------------------------------------------+ class PeriodicKernel : public IKernel { private: double sigma_f; // Signal variance (oscillation amplitude) double length; // Scale length (smoothness of oscillations) double period; // Oscillation period int n; // Number of rows in X1 int m; // Number of rows in X2 matrix m_K; // Covariance matrix matrix m_D; // D matrix (sum of derivatives with respect to d [2 * sin^2(M_PI*distance/period] ) matrix m_distance_dim[]; // Array of difference matrices for each dimension |x_i_k - x_j_k| int m_d_cols; // Number of dimensions (features) (X.Cols()) public: //Constructor PeriodicKernel(double ls, double sf, double p) : sigma_f(sf), length(ls), period(p){} //Destructor ~PeriodicKernel() {} matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); int d = (int)X1.Cols(); m_d_cols = d; m_D = matrix::Zeros(n, m); ArrayResize(m_distance_dim, d); for (int k = 0; k < d; k++) { vector x1_k = X1.Col(k); vector x2_k = X2.Col(k); // Create a matrix where each column is a x1_k vector // (n x 1) * (1 x m) = (n x m) matrix x1_k_copy = x1_k.Outer(vector::Ones(m)) ; // Create a matrix where each row is the transposed vector of x2_k // (n x 1) * (1 x m) = (n x m) matrix x2_k_copy = vector::Ones(n).Outer(x2_k); // Calculate the matrix of all pairwise absolute differences between the elements of vectors x1_k and x2_k matrix distance = MathAbs(x1_k_copy - x2_k_copy); m_distance_dim[k] = distance; // Cache 'distance' for each dimension (feature) matrix sin_term = MathSin(M_PI * distance / period); sin_term = 2.0 * sin_term * sin_term; // 2 * sin^2(u) m_D += sin_term; // Sum up element by element } // Calculate the K matrix m_K = sigma_f * sigma_f * MathExp(-1 * m_D / (length * length)); return m_K; } matrix ComputeDerivative( int param_index) override { matrix dK(n, m); switch (param_index) { case 0: { // Derivative with respect to sigma_f // dK/d(sigma_f) = (2 / sigma_f) * K dK = (2.0 / sigma_f) * m_K; break; } case 1: { // Derivative with respect to 'length' // dK/d(length) = K * (2 * D / length^3) dK = m_K * ( 2.0 * m_D / (length * length * length) ); break; } case 2: { // Derivative with respect to 'period' // dK/d(period) = K * (-1/l^2) * dD/d(period) // dD/d(period) = Sum{k=1:d} (2 * pi * (x_ik - x_jk) / period^2) * sin(2*pi*(x_ik - x_jk)/period) matrix dD_dp = matrix::Zeros(n, m); for (int k = 0; k < m_d_cols; k++) { // Loop through each dimension (column) of the data matrix distance_k = m_distance_dim[k]; // matrix of absolute differences for the k-th dimension matrix u = M_PI * distance_k / period; // Calculate the argument u_k = pi * |x_ik - x_jk| / period matrix sin_2u = MathSin(2.0 * u); // Calculate sin(2 * u_k) // Use the already calculated u to calculate // -pi * |x_ik - x_jk| / period^2 matrix second_term = (-1*u) / period; // Accumulate dD/d(period) for the current measurement dD_dp += 2.0 * sin_2u * second_term; } // unite all parts of the derivative // dK = K * (-1/length^2) * dD_dp dK = m_K * (-1.0 / (length * length)) * dD_dp; break; } } return dK; } vector GetHyperparameters() const override { vector params(3); params[0] = sigma_f; params[1] = length; params[2] = period; return params; } void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 3) { sigma_f = params[0]; length = params[1]; period = params[2]; } } int GetNumHyperparameters() const override { return 3; } string GetName() const override { return "PeriodicKernel"; } };
Ableitung eines periodischen Kernels nach dem Parameter period:
Die Berechnung der Ableitung nach dem Parameter period ist am schwierigsten, da sich period innerhalb der trigonometrischen Funktion befindet, die wiederum Teil der Exponentialfunktion ist.
Die Ableitungsformel basiert auf der Kettenregel, wobei K von D abhängt und D über die Sinusterme von period abhängt. Die ComputeDerivative-Implementierung berechnet iterativ die dD_dp-Komponente für jede Dimension und berechnet dann die endgültige Ableitung:
Dabei wird dD_dp als Summe der Ableitungen 2 * sin^2(M_PI*Abstand/Periode) für alle Dimensionen (Merkmale) berechnet.
Zusammengesetzte Kernel
Gauß-Prozesse ermöglichen es, einfache Kovarianz-Kernel zu kombinieren, um komplexere Modelle zu erstellen, die verschiedene Aspekte der Daten (z. B. Trend, Periodizität, Rauschen) beschreiben können. Diese Funktionalität wird durch zwei zusammengesetzte Kernel implementiert:
- SumKernel – kombiniert mehrere Kernel durch Summierung ihrer Kovarianzmatrizen.
- ProductKernel – kombiniert Kernel durch elementweise Multiplikation ihrer Kovarianzmatrizen.
Beide Klassen verwalten die Hyperparameter ihrer untergeordneten Kernel, sammeln sie für die Optimierung in einem einzigen gemeinsamen Vektor und verteilen sie anschließend wieder.
Die Klasse SumKernel
//+------------------------------------------------------------------+ //| Class for the sum of kernels | //+------------------------------------------------------------------+ class SumKernel : public IKernel { private: IKernel* kernels[]; public: // Constructor SumKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } //Destructor ~SumKernel() { for (int i = 0; i < ArraySize(kernels); i++) { if (kernels[i] != NULL) { delete kernels[i]; } } ArrayFree(kernels); } matrix Compute(const matrix &X1, const matrix &X2) override { int n = (int)X1.Rows(); int m = (int)X2.Rows(); matrix sum = matrix::Zeros(n,m); for (int i = 0; i < ArraySize(kernels); i++) { sum += kernels[i].Compute(X1, X2); // Sum the matrices from each child kernel } return sum; } matrix ComputeDerivative(int param_index) override { int current_param_offset = 0; // Start at offset 0 for the first kernel for (int i = 0; i < ArraySize(kernels); i++) { // Get the number of hyperparameters of the current child kernel int num_params_current_kernel = kernels[i].GetNumHyperparameters(); // Check if the global param_index belongs to the current child kernel // param_index should be: // - greater than or equal to the current offset // - strictly less than the current offset + the number of parameters of the current kernel if (param_index >= current_param_offset && param_index < current_param_offset + num_params_current_kernel) { // If param_index belongs to this kernel, calculate its local index int local_param_index = param_index - current_param_offset; // Call ComputeDerivative on the found child kernel and // pass it the local index. // Since the derivative of the sum of kernels with respect to the parameter of one kernel is equal to the derivative // of this particular kernel by this parameter, we can immediately return the result return kernels[i].ComputeDerivative(local_param_index); } // If param_index does not belong to the current kernel, // increase the offset to check the next kernel current_param_offset += num_params_current_kernel; } Print("Error: SumKernel::ComputeDerivative - Parameter index ", param_index, " out of bounds"); return matrix::Zeros(1, 1); } //+---------------------------------------------------------------------------------+ //| The function returns the values of all hyperparameters of the composite kernel | //+---------------------------------------------------------------------------------+ vector GetHyperparameters() const override { vector all_params(GetNumHyperparameters()); int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Get the values of the parameters of the current child kernel vector kernel_params = kernels[i].GetHyperparameters(); for (int j = 0; j < (int)kernel_params.Size(); j++) { // Copy the parameters into the common vector all_params[current_idx + j] = kernel_params[j]; } // Update the offset for the next kernel current_idx += (int)kernel_params.Size(); } return all_params; } //+----------------------------------------------------------------------+ //|The function takes a vector of values of all hyperparameters, breaks | //|it into the appropriate parts and assigns them to each child kernel | //+----------------------------------------------------------------------+ void SetHyperparameters(const vector ¶ms) override { int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { // Get the number of hyperparameters of the current child kernel int num_params = kernels[i].GetNumHyperparameters(); vector sub_params(num_params); // Copy the corresponding parameters from the general vector for (int j = 0; j < num_params; j++) { sub_params[j] = params[current_idx + j]; } // Call the SetHyperparameters() method on the current child kernel, // passing it only its own parameters, collected in the sub_params vector kernels[i].SetHyperparameters(sub_params); current_idx += num_params; } } //+--------------------------------------------------------------------------------------+ //| Return the total number of hyperparameters for the entire SumKernel composite kernel | //+--------------------------------------------------------------------------------------+ int GetNumHyperparameters() const override { int total_params = 0; for (int i = 0; i < ArraySize(kernels); i++) { total_params += kernels[i].GetNumHyperparameters(); } return total_params; } string GetName() const override { string name = "Sum("; for (int i = 0; i < ArraySize(kernels); i++) { name += kernels[i].GetName(); if (i < ArraySize(kernels) - 1) name += ","; } name += ")"; return name; } //+------------------------------------------------------------------+ //| Provide external access to the list of child kernels | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { // Copy pointers from the internal kernels array to the passed // output_kernels external array. This provides access to child // kernels in the GaussianProcess::Fit() method, where it is necessary to iterate through // all elementary kernels to set optimization boundaries // of their hyperparameters. ArrayCopy(output_kernels, kernels); } };
- Die Compute-Methode berechnet die endgültige K-Kovarianzmatrix als einfache Summe der von jedem untergeordneten Kernel zurückgegebenen Kovarianzmatrizen. Hier kommt das Prinzip des Polymorphismus zum Einsatz: Obwohl das Kernels-Array Zeiger des Basistyps IKernel* speichert, ruft der Aufruf von kernels[i].Compute(X1, X2) ruft tatsächlich die spezifische Implementierung von Compute für jeden einzelnen untergeordneten Kernel auf. Dies ermöglicht es SumKernel, mit beliebigen von IKernel abgeleiteten Kerneln zu arbeiten.
- Die Grundidee hinter der ComputeDerivative-Methode ist, dass die Ableitung der Kovarianzfunktion der Summe von Kerneln in Bezug auf einen spezifischen Hyperparameter gleich der Ableitung der Kovarianzfunktion nur desjenigen untergeordneten Kernels ist, zu dem der gegebene Hyperparameter gehört. Dies bedeutet, dass ComputeDerivative den entsprechenden untergeordneten Kernel anhand von param_index findet und die Ableitung ausschließlich von diesem Kernel zurückgibt. Ableitungen in Bezug auf die Parameter anderer untergeordneter Kernel werden als null betrachtet.
Die Klasse ProductKernel
//+------------------------------------------------------------------+ //| Class for the product of kernels | //+------------------------------------------------------------------+ class ProductKernel : public IKernel { private: IKernel* kernels[]; matrix m_X1; matrix m_X2; int n; int m; public: // Constructor: copies pointers to child kernels ProductKernel(IKernel* &input_kernels[]) { ArrayResize(kernels, ArraySize(input_kernels)); ArrayCopy(kernels, input_kernels); } // Destructor ~ProductKernel() { for (int i = 0; i < ArraySize(kernels); i++) { if (kernels[i] != NULL){ delete kernels[i]; } } ArrayFree(kernels); } matrix Compute(const matrix &X1, const matrix &X2) override { n = (int)X1.Rows(); m = (int)X2.Rows(); m_X1 = X1; m_X2 = X2; matrix product = matrix:: Ones(n,m); for (int i = 0; i < ArraySize(kernels); i++) { // Polymorphism: call the Compute() method of each child kernel, // regardless of its specific type, and multiply the result element by element product *= kernels[i].Compute(X1, X2); } return product; } matrix ComputeDerivative(int param_index) override { matrix dK_prod(n, m); // Final derivative matrix int current_param_offset = 0; int target_kernel_idx = -1; // Kernel index the hyperparameter belongs to int local_param_index = -1; // Local index of the hyperparameter in this kernel // Step 1: Find the kernel and local index of the hyperparameter for (int i = 0; i < ArraySize(kernels); i++) { int num_params_current_kernel = kernels[i].GetNumHyperparameters(); if (param_index >= current_param_offset && param_index < current_param_offset + num_params_current_kernel) { target_kernel_idx = i; local_param_index = param_index - current_param_offset; break; } current_param_offset += num_params_current_kernel; } if (target_kernel_idx == -1) { Print("Error: ProductKernel::ComputeDerivative - Parameter index ", param_index, " out of bounds"); return matrix::Zeros(1, 1); } // Step 2: Calculate dK_k / d_theta_j for the target kernel matrix dK_target_kernel = kernels[target_kernel_idx].ComputeDerivative(local_param_index); // Step 3: Calculate the K_m product for all other kernels (m != k) matrix other_kernels_product = matrix::Ones(n, m); for (int i = 0; i < ArraySize(kernels); i++) { if (i != target_kernel_idx) { // Skip the target kernel other_kernels_product = other_kernels_product * kernels[i].Compute(m_X1, m_X2); } } // Step 4: Multiply the results element by element dK_prod = dK_target_kernel * other_kernels_product; return dK_prod; } // Collect the hyperparameter values of all child kernels into a single vector vector GetHyperparameters() const override { vector all_params(GetNumHyperparameters()); int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { vector kernel_params = kernels[i].GetHyperparameters(); for (int j = 0; j < (int)kernel_params.Size(); j++) { all_params[current_idx + j] = kernel_params[j]; } current_idx += (int)kernel_params.Size(); } return all_params; } void SetHyperparameters(const vector ¶ms) override { int current_idx = 0; for (int i = 0; i < ArraySize(kernels); i++) { int num_params_kernel = kernels[i].GetNumHyperparameters(); vector sub_params(num_params_kernel); for (int j = 0; j < num_params_kernel; j++) { sub_params[j] = params[current_idx + j]; } kernels[i].SetHyperparameters(sub_params); current_idx += num_params_kernel; } } int GetNumHyperparameters() const override { int total_params = 0; for (int i = 0; i < ArraySize(kernels); i++) { total_params += kernels[i].GetNumHyperparameters(); } return total_params; } string GetName() const override { string name = "Prod("; for (int i = 0; i < ArraySize(kernels); i++) { name += kernels[i].GetName(); if (i < ArraySize(kernels) - 1) name += "*"; } name += ")"; return name; } //+------------------------------------------------------------------+ //| Provide external access to the list of child kernels | //+------------------------------------------------------------------+ void GetKernels(IKernel* &output_kernels[]) const { ArrayCopy(output_kernels, kernels); } };
- Die Compute-Methode berechnet die endgültige K-Kovarianzmatrix als elementweise Produkt der von jedem untergeordneten Kernel zurückgegebenen Kovarianzmatrizen. Wie bei SumKernel wird hier Polymorphismus verwendet, um die Compute-Methode des entsprechenden untergeordneten Kernels aufzurufen.
- Die Methode ComputeDerivative: Um die Ableitung von K in Bezug auf einen Hyperparameter zu berechnen, der zu einem bestimmten K_p-untergeordneten Kernel gehört, wendet ProductKernel die Produktregel (Leibniz-Regel) an. Die Ableitung von K in Bezug auf einen gegebenen Parameter entspricht dem Produkt aus der Ableitung der K_p-Matrix in Bezug auf diesen Parameter und dem elementweisen Produkt der Kovarianzmatrizen aller anderen untergeordneten Kernel, genommen ohne die Ableitung:
ILikelihood-Schnittstelle
Die ILikelihood-Schnittstelle definiert, wie jede Likelihood-Funktion in unserer Bibliothek funktionieren sollte. Sie ermöglicht uns, mit verschiedenen Datentypen zu arbeiten, sei es mit kontinuierlichen Werten (wie bei der Regression) oder diskreten Kategorien (wie bei der Klassifikation).
//+------------------------------------------------------------------+ //|Interface for likelihood functions | //+------------------------------------------------------------------+ interface ILikelihood { // Compute the p(y|f) log-likelihood for latent values f and observed values y virtual double LogLikelihood(const vector &f, const vector &y) = 0; //----------------------------- Derivatives of the logarithm of the likelihood with respect to f ----------------- // Calculate the vector of the first derivative of the log-likelihood with respect to f (dlp/df) virtual vector LogLikelihoodGradient(const vector &f, const vector &y) = 0; // Compute the second derivative (Hessian) matrix of the log-likelihood with respect to f (d^2lp/df df^T) virtual matrix LogLikelihoodHessian(const vector &f, const vector &y) = 0; // Calculate the vector of the third derivative of the log-likelihood with respect to f (d^3lp/df^3) virtual vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) = 0; //----------------------------- Derivatives of the logarithm of the likelihood with respect to the parameters ----------------- // The first derivative of the log-likelihood with respect to the j-th likelihood hyperparameter. virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) = 0; // Hessian derivative of the log-likelihood with respect to the j-th likelihood hyperparameter. virtual matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) = 0; // Return the name of the likelihood function virtual string GetName() const = 0; // Return the current values of the likelihood hyperparameters virtual vector GetHyperparameters() const = 0; // Set new values for the likelihood hyperparameters virtual void SetHyperparameters(const vector ¶ms) = 0; // Return the number of likelihood hyperparameters virtual int GetNumHyperparameters() const = 0; };
GaussianLikelihood-Klasse
GaussianLikelihood implementiert eine Likelihood-Funktion für Regressionsprobleme unter der Annahme von normalverteiltem Rauschen in den beobachteten Daten.
//+------------------------------------------------------------------+ //| Gaussian likelihood class (for regression problems) | //+------------------------------------------------------------------+ class GaussianLikelihood : public ILikelihood { private: double m_noise_sigma; // Noise parameter (standard deviation) public: // Constructor GaussianLikelihood(double initial_noise_sigma) : m_noise_sigma(initial_noise_sigma) {} // Return the current value of the hyperparameter vector GetHyperparameters() const override { vector params(1); params[0] = m_noise_sigma; return params; } // Set a new hyperparameter value void SetHyperparameters(const vector ¶ms) override { if (params.Size() == 1) { m_noise_sigma = params[0]; } } // Return the number of hyperparameters int GetNumHyperparameters() const override { return 1; } // Calculate the log p(y|f) log-likelihood for a Gaussian distribution double LogLikelihood(const vector &f, const vector &y) override { int n = (int)y.Size(); double noise_variance = m_noise_sigma * m_noise_sigma; vector residual = y - f; // Formula for the logarithm of the density of a multivariate normal distribution return -0.5 / noise_variance * (residual @ residual) - 0.5 * n * MathLog(2 * M_PI * noise_variance); } //------------------ Derivatives with respect to f ----------------------- // Calculate the vector of the first derivative of the log-likelihood with respect to f (dlp/df) vector LogLikelihoodGradient(const vector &f, const vector &y) override { // d(log p(y|f))/df = (y - f) / sigma^2 return (y - f) / (m_noise_sigma * m_noise_sigma); } // Compute the second derivative (Hessian) matrix of the log-likelihood with respect to f (d^2lp/dfdf^T) matrix LogLikelihoodHessian(const vector &f, const vector &y) override { // d^2(log p(y|f))/dfdf^T = -1/sigma^2 * I int n = (int)y.Size(); return matrix::Identity(n, n) * (-1.0 / (m_noise_sigma * m_noise_sigma)); } // Calculate the vector of the third derivative of the log-likelihood with respect to f vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) override { // For Gaussian likelihood d^3(log p(y|f))/df^3 = 0 return vector::Zeros((int)y.Size()); } // ------------------------------ Derivatives with respect to the parameter ----------------------------------- // // Calculate the first derivative of the log-likelihood with respect to the j-th likelihood hyperparameter virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { // Gaussian likelihood has only 1 hyperparameter: m_noise_sigma (index 0) if (param_index == 0) { // Formula for d(log p(y|f))/d(sigma_n): // = (y-f)^T(y-f) / sigma_n^3 - N / sigma_n vector y_minus_f = y - f; double sn3 = m_noise_sigma * m_noise_sigma * m_noise_sigma; double term1 = (y_minus_f @ y_minus_f) / sn3; double term2 = y.Size() / m_noise_sigma; // N / sigma_n return term1 - term2; } else { return 0.0; } } // Calculate the Hessian derivative of the log-likelihood with respect to the j-th likelihood hyperparameter matrix LogLikelihoodHessianDerivative(const vector &f, const vector &y, int param_index) override { int n = (int)y.Size(); if (param_index == 0) { // Derivative of Hessian with respect to m_noise_sigma: d/d(sigma_n) [-1/sigma_n^2 * I] = (2/sigma_n^3) * I double derivative_coeff = 2.0 / (m_noise_sigma * m_noise_sigma * m_noise_sigma); return matrix::Identity(n, n) * derivative_coeff; } else { return matrix::Zeros(n, n); } } string GetName() const override { return "GaussianLikelihood"; } };
Schlüsselmethoden:
- Der GaussianLikelihood-Konstruktor initialisiert einen einzelnen Hyperparameter, die Rausch-Standardabweichung.
- LogLikelihood() – berechnet die p(y∣f) Log-Likelihood für die Wahrscheinlichkeitsdichtefunktion der multivariaten Normalverteilung:
- LogLikelihoodGradient() – berechnet die erste Ableitung der Log-Likelihood in Bezug auf die latenten Werte von f. Dies ist ein Vektor, bei dem jedes Element gleich ist zu:
- LogLikelihoodHessian() – berechnet die zweite Ableitung (Hesse-Matrix) der Log-Likelihood in Bezug auf f. Für die Gauß-Likelihood ist die Hesse-Matrix eine Diagonalmatrix:
- LogLikelihoodThirdDerivative() – berechnet die dritte Ableitung. Für die Gauß-Likelihood sind alle Ableitungen über der zweiten gleich null.
- LogLikelihoodGradientParam() – berechnet die erste Ableitung der Log-Likelihood in Bezug auf den Hyperparameter m_noise_sigma:
- LogLikelihoodHessianDerivative() – berechnet die Hesse-Ableitung der Log-Likelihood in Bezug auf den Hyperparameter m_noise_sigma:
Die Klasse LogitLikelihood
LogitLikelihood wird für die binäre Klassifikation verwendet, bei der die beobachteten Zielwerte y die Werte {−1,+1} annehmen. Sie setzt die latente Funktion f über eine Sigmoid-Funktion in Beziehung zur Wahrscheinlichkeit der Klassenzugehörigkeit.
Die LogitLikelihood-Implementierung stellt die Hilfsfunktionen sigmoid() und Softplus() mit Prüfungen auf große/kleine Eingabewerte bereit.
//+------------------------------------------------------------------+ //| Class for Logit likelihood (y {-1, +1}) | //+------------------------------------------------------------------+ class LogitLikelihood : public ILikelihood { public: // logit_sigmoid(x) = 1 / (1 + exp(-x)) double sigmoid(double x) const { if (x > 100.0) return 1.0; // Avoid NaN occurrences if (x < -100.0) return 0.0; return 1.0 / (1.0 + MathExp(-x)); } // Softplus log(1 + exp(x)) function double Softplus(double x) const { if (x > 100.0) return x; // For very large x, log(1 + exp(x)) ~ x if (x < -100.0) return MathExp(x); // For very small x, log(1 + exp(x)) ~ exp(x) return MathLog(1.0 + MathExp(x)); } public: // Constructor (logit likelihood has no hyperparameters of its own) LogitLikelihood() {} vector GetHyperparameters() const override { return vector::Zeros(0); } void SetHyperparameters(const vector ¶ms) override {} // Return the number of likelihood hyperparameters int GetNumHyperparameters() const override { return 0; } // --- Calculate the log-likelihood: log p(y|f) --- // Formula: sum_i (-log(1 + exp(-y_i * f_i))) double LogLikelihood(const vector &f, const vector &y) override { int n = (int)y.Size(); double total_log_likelihood = 0.0; for (int i = 0; i < n; i++) { total_log_likelihood += -Softplus(-y[i] * f[i]); } return total_log_likelihood; } // Compute the gradient of the log-likelihood with respect to f // Formula: d(log p(y_i|f_i))/df_i = y_i * sigma(-y_i * f_i) vector LogLikelihoodGradient(const vector &f, const vector &y) override { int n = (int)y.Size(); vector grad(n); for (int i = 0; i < n; i++) { grad[i] = y[i] * sigmoid(-y[i] * f[i]); // grad[i] = y[i] * (1 - sigmoid(y[i] * f[i])); equivalent to } return grad; } // Compute the Hessian of the log-likelihood of f (the H diagonal matrix) // Formula: d^2(log p(y_i|f_i))/df_i^2 = -sigma(y_i*f_i)(1 - sigma(y_i*f_i)) matrix LogLikelihoodHessian(const vector &f, const vector &y) override { int n = (int)y.Size(); matrix H = matrix::Identity(n, n); for (int i = 0; i < n; i++) { double sig = sigmoid(y[i] * f[i]); H[i][i] = -1*sig * (1.0 - sig); } return H; } // Calculate the third derivative of the log-likelihood with respect to f // Formula: d^3(log p(y_i|f_i))/df_i^3 = y_i * sigma(-y_i f_i) * (1 - sigma(-y_i f_i)) * (1 - 2*sigma(-y_i f_i)) vector LogLikelihoodThirdDerivative(const vector &f, const vector &y) override { int n = (int)y.Size(); vector third_deriv(n); for (int i = 0; i < n; i++) { double sig_neg_yf = sigmoid(-y[i] * f[i]); third_deriv[i] = y[i] * sig_neg_yf * (1.0 - sig_neg_yf) * (1.0 - 2.0 * sig_neg_yf); // equivalent to // double sig_yf = sigmoid(y[i] * f[i]); // third_deriv[i] = -y[i] * sig_yf * (1.0 - sig_yf) * (1.0 - 2.0 * sig_yf); } return third_deriv; } // LogitLikelihood has no hyperparameters virtual double LogLikelihoodGradientParam(const vector &f, const vector &y, int param_index) override { return 0.0; } // LogitLikelihood has no hyperparameters matrix LogLikelihoodHessianDerivative(const vector &f_latent, const vector &y, int param_index) override { int n = (int)y.Size(); return matrix::Zeros(n, n); } string GetName() const override { return "LogitLikelihood"; } };
Schlüsselmethoden:
- LogLikelihood() – berechnet die Log-Likelihood logp(y∣f). Für eine binäre Klassifikation mit y {−1,+1} ist dies die Summe der Logarithmen der binären Cross-Entropy-Likelihood:
- LogLikelihoodGradient() – berechnet die erste Ableitung der Log-Likelihood in Bezug auf f. Jedes Element des Gradientenvektors ist gleich:
- LogLikelihoodHessian() – berechnet die zweite Ableitung (Hesse-Matrix) der Log-Likelihood in Bezug auf f. Für LogitLikelihood ist die Hesse-Matrix eine Diagonalmatrix, deren Diagonalelemente gleich sind:
- LogLikelihoodThirdDerivative() – berechnet die dritte Ableitung der Log-Likelihood in Bezug auf f. Jedes Element des Vektors ist gleich:
Hilfsfunktionen und Strukturen
Die Datei StructUtils.mqh enthält eine Reihe von Aufzählungen, Datenstrukturen und Funktionen, die das Arbeiten mit GP erleichtern.
enum PredictMode { PROBIT = 0, // Probit Approximation NUM_INTEGR = 1, // Numerical integration MONTE_CARLO = 2 // Monte Carlo }; //--- Structure for prediction results struct GPPredictionResult { vector mu_f_star; // Posterior mean for the latent f* function on new data matrix Sigma_f_star; // Posterior covariance for the latent f* function on new data // Regression-only fields (GaussianLikelihood) vector mu_y_star; // Posterior mean of observations y* matrix Sigma_y_star; // Posterior covariance of observations y* // Fields for classification only (LogitLikelihood, ProbitLikelihood, etc.) vector predicted_probabilities; // Probabilities p(y*=+1 | X*, D) vector predicted_labels; // Predicted y* labels (+1 or -1) }; // Structure for the inference result struct GPInferenceResult { double nlml_value; // Negative logarithm of the marginal likelihood vector nlml_gradient; // NLML gradient by kernel and likelihood hyperparameters // Inference results on training data: vector mu_f_train; // Posterior mean of the latent function f on the training data matrix Sigma_f_train; // Posterior covariance of the latent function f on the training data (K^−1+W)^−1 matrix L_K_noisy; // Cholesky decomposition K_noisy = cholesky(K + s2n*I) matrix L_B; // Cholesky decomposition B = I + W^0.5 @ K @ W^0.5 matrix sW; // sW = W^0.5; vector sW_diag; // sW diagonal matrix sW_K; // sW @ K vector alpha; // (K_noisy)^-1 * y matrix H; // Hessian of the log likelihood d^2 log p(y|f)/df^2 (for LaplaceInference) bool success; // Flag indicating success of the inference function };
- Die Enumeration PredictMode definiert verschiedene Modi für die Durchführung von Vorhersagen in einem Klassifikationsmodell.
- Die Struktur GPInferenceResult wird verwendet, um alle während der Inferenz erhaltenen Schlüsselergebnisse zu speichern.
- Die Struktur GPPredictionResult ist darauf ausgelegt, alle Ergebnisse zu speichern, die nach der Durchführung einer Vorhersage für neue (Test-)Daten (X_test) erhalten wurden.
- DiagonalTrace
//+-------------------------------------------------------------------+ //| Compute the trace of the product of two matrices using only | //| diagonal elements of the final matrix | //+-------------------------------------------------------------------+ double DiagonalTrace(const matrix &A, const matrix &B) { ulong m = A.Rows(); ulong n = A.Cols(); ulong n_b = B.Rows(); ulong p = B.Cols(); // Find the minimum size for the diagonal ulong k = MathMin(m, p); double tr = 0.0; for(ulong i = 0; i < k; i++) { // Dot product of the i-th row of A and the i-th column of B tr += A.Row(i)@B.Col(i); } return tr; }
Berechnet die Spur des Produkts zweier Matrizen (Tr(A @ B)), ohne die gesamte Produktmatrix bilden zu müssen. Dies geschieht durch Summieren der Skalarprodukte der Zeilen von Matrix A und der entsprechenden Spalten von Matrix B. Dies ermöglicht eine erhebliche Reduzierung des Rechenaufwands im Vergleich zur vollständigen Matrixmultiplikation.
- cho_solve (OpenBLAS LinearEquationsSolutionTriangular-Wrapper)
//+-----------------------------------------------------------------------+ //| Solve the system A*X = B using the Cholesky decomposition A = LL^T | //| Parameters: | //| c: Lower triangular matrix L from the Cholesky decomposition | //| b: right side of the system (matrix) | //| Return: matrix X - system solution | //+---------------------------------------------------------------------+ matrix cho_solve(const matrix &c, const matrix &b) { // Check if 'c' matrix is lower triangular one if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular."); return matrix::Zeros(c.Rows(), b.Cols()); } // Step 1: Solve L * Y = B // Y - intermediate matrix, the result of solving L^-1 * B matrix Y; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_N, b, Y)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L * Y = B failed. Error Code: %d", GetLastError()); return matrix::Zeros(c.Rows(), b.Cols()); } // Step 2: Solve L^T * X = Y // X - L^T * X = Y system solution, which is equivalent to (L^T)^-1 * Y matrix X; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_T, Y, X)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L^T * X = Y failed. Error Code: %d", GetLastError()); return matrix::Zeros(c.Rows(), b.Cols()); } return X; } //+------------------------------------------------------------------+ //| Solve the system Ax = b using the Cholesky decomposition A = LL^T| //| Parameters: | //| c: Lower triangular matrix L from the Cholesky decomposition | //| b: right side of the system (vector) | //| Return: vector x - system solution | //+------------------------------------------------------------------+ vector cho_solve(const matrix &c, const vector &b) { // Check if 'c' matrix is lower triangular one if (!c.IsLowerTriangular()) { Print("Error: cho_solve - Input matrix 'c' is not lower triangular"); return vector::Zeros(c.Rows()); } // Step 1: Solve L * y = b // y - intermediate vector, the result of solving L^-1 * b vector y; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_N, b, y)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L * y = b failed. Error Code: %d", GetLastError()); return vector::Zeros(c.Rows()); } // Step 2: Solve L^T * x = y // x - final solution of the L^T * x = y system, which is equivalent to (L^T)^-1 * y vector x; if (!c.LinearEquationsSolutionTriangular(EQUATIONSFORM_T, y, x)) { PrintFormat("Error: cho_solve- LinearEquationsSolutionTriangular L^T * x = y failed. Error Code: %d", GetLastError()); return vector::Zeros(c.Rows()); } return x; }
Die Funktion cho_solve ist darauf ausgelegt, Systeme linearer Gleichungen Ax=b (oder AX=B für Matrix B) effizient zu lösen, wobei die Matrix A aus ihrer Cholesky-Zerlegung A=LL^T gewonnen wird. Hierbei ist L eine untere Dreiecksmatrix.
Die direkte Berechnung von A^−1 (A.Inv()) ist ressourcenintensiv und kann numerisch instabil sein. Die Cholesky-Zerlegung hingegen bietet einen wesentlich effizienteren und robusteren Ansatz.
Innerhalb von cho_solve werden zwei Hauptoperationen durchgeführt:
- Vorwärtssubstitution: Das System LY=B (oder Ly=b) wird gelöst, wobei L die Eingabematrix c und Y (oder y) das Zwischenergebnis ist. Da L eine Dreiecksmatrix ist, ist diese Lösung relativ schnell.
- Rückwärtssubstitution: Anschließend wird das System L^TX=Y (oder L^Tx=y) gelöst, wobei L^T die Transponierte der Matrix L ist. Dies ist ebenfalls effizient, da L^T eine obere Dreiecksmatrix ist.
Diese Operationen werden mithilfe der LinearEquationsSolutionTriangular-Funktionen aus der OpenBLAS-Hochleistungsbibliothek für lineare Algebra implementiert. Dies sorgt für erhebliche Rechenbeschleunigungen und macht sie für ressourcenintensive Modelle des maschinellen Lernens wie Gauß-Prozesse unverzichtbar.
IInference-Schnittstelle
interface IInference { // This method performs inference (i.e., computes the posterior distribution of f) // and returns the components needed for NLML and prediction virtual void Infer(const matrix &X, const vector &y, IKernel *kernel, ILikelihood *likelihood,GPInferenceResult &result) = 0; virtual string GetName() const = 0; };
Die Schnittstelle ermöglicht das Umschalten zwischen verschiedenen Inferenzmethoden, wie z. B. exakter Inferenz für Gauß-Likelihood und approximativen Methoden für nicht-Gaußsche Fälle, ohne die zentrale Logik von GaussianProcess zu ändern.
Die Infer-Methode ist zentral für diese Schnittstelle. Sie berechnet die Posterior-Verteilung der latenten Funktion f durch und gibt die für weitere Berechnungen benötigten Komponenten (NLML und Vorhersagen) zurück. Sie akzeptiert die folgenden Argumente:
- X – Matrix der Trainingsmerkmale,
- y – Vektor der Trainingslabels,
- kernel – Zeiger auf das Objekt IKernel,
- likelihood – Zeiger auf das ILikelihood-Objekt,
- result – eine Referenz auf die GPInferenceResult-Struktur, die zum Speichern aller Inferenzergebnisse verwendet wird, einschließlich:
- negativer Logarithmus der marginalen Likelihood (NLML), der für die Hyperparameteroptimierung benötigt wird,
- Vektor der NLML-Gradienten über alle Kernel- und Likelihood-Hyperparameter,
- Hilfsmatrizen und -vektoren: (L_K_noisy, L_B, sW, mu_f_train, Sigma_f_train, alpha), die zur Berechnung der NLML-Gradienten und der anschließenden Vorhersage neuer Daten verwendet werden.
Klasse ExactInference
//+------------------------------------------------------------------+ //| ExactInference For Gaussian likelihood only | //+------------------------------------------------------------------+ class ExactInference : public IInference { public: ExactInference() {} void Infer(const matrix &X, const vector &y, IKernel *kernel, ILikelihood *likelihood,GPInferenceResult &result) override { if(likelihood.GetName() != "GaussianLikelihood") { Print("ExactInference supports only GaussianLikelihood"); result.nlml_value = DBL_MAX; return; } result.success = false; int n = (int)y.Size(); // Calculate K with current kernel parameters matrix K = kernel.Compute(X, X); // Get the noise variance vector likelihood_params = likelihood.GetHyperparameters(); double sigma = likelihood_params[0]; double variance = sigma * sigma; double jitter = 1e-6; matrix K_noisy = K + matrix::Identity(n, n) * (variance + jitter); if(!K_noisy.Cholesky(result.L_K_noisy)) { PrintFormat("Error: Cholesky decomposition failed. Error Code: %d",GetLastError()); result.nlml_value = DBL_MAX; return; } //------------- Algorithm 2.1 GPML --------------------------- //Calculate alpha (K_noisy^-1 * y) - O(N^2) result.alpha = cho_solve(result.L_K_noisy, y); //+------------------------------------------------------------------+ //| NLML = 0.5y^T(K+σ^2*I)^-1*y + 0.5*log|K+σ^2*I| +n/2*Log(2π) | //+------------------------------------------------------------------+ //+------------------------------------------------------------------+ //| NLML = 0.5y^T*alpha + 0.5*log|K+σ^2*I| +n/2*Log(2π) | //+------------------------------------------------------------------+ // --- Calculate the first term in the NLML equation: double data_term = 0.5 * (y @ result.alpha); // --- Calculate the second term: 1/2 * log|K + sigma^2*I| = sum(log(L_ii)) double log_det = MathLog(result.L_K_noisy.Diag()).Sum(); // --- Calculate the third term: n/2 * log(2π) double const_term = 0.5 * n * MathLog(2 * M_PI); //--- NLML result.nlml_value = data_term + log_det + const_term; // -------------------- NLML gradients calculations -------------------------------------- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // Calculate (K_noisy^-1) matrix K_noisy_inv = cho_solve(result.L_K_noisy, matrix::Identity(n, n)); // 1. Gradients on kernel hyperparameters vector kernel_hyperparams = kernel.GetHyperparameters(); // Calculate aatK_noisy_inv once before the loop - O(N^2) matrix aatK_noisy_inv = result.alpha.Outer(result.alpha) - K_noisy_inv; matrix dK_dtheta; for(int i = 0; i < (int)kernel_hyperparams.Size(); i++) { dK_dtheta = kernel.ComputeDerivative(i); // O(N^2) result.nlml_gradient[current_grad_idx] = -0.5 * DiagonalTrace(aatK_noisy_inv, dK_dtheta); current_grad_idx++; } // 2. Gradient over noise hyperparameter // dNLML/d(sigma^2) = 0.5 * (trace(K_noisy^-1) - alpha^T * alpha ) double dNLML_d_sigma2n = 0.5 * (K_noisy_inv.Trace() - result.alpha @ result.alpha); result.nlml_gradient[current_grad_idx] = dNLML_d_sigma2n * (2.0 * sigma); result.success = true; } string GetName() const override { return "ExactInference"; } };
Die Klasse ExactInference ist für die Durchführung exakter Inferenz in GP konzipiert und nur bei Verwendung der Gauß-Likelihood anwendbar. Der Hauptfokus ihrer Implementierung liegt auf der effizienten Berechnung der NLML und ihrer analytischen Gradienten. Analytische Gradienten, die direkt mithilfe der Gleichungen berechnet werden, bieten eine höhere Genauigkeit und beschleunigen den Optimierungsprozess im Vergleich zu numerischen Methoden erheblich.
Der Inferenzalgorithmus besteht aus mehreren wichtigen Schritten:
- Bildung der Kovarianzmatrix Knoisy: Zuerst wird die Kernel-Kovarianzmatrix K für die Trainingsdaten berechnet. Dann werden die Rauschvarianz (Varianz = sigma * sigma) und ein kleiner positiver Jitter-Wert (1e-6) zu den Diagonalelementen von K addiert. Dies bildet die Matrix Knoisy;
- Berechnung des Alpha-Vektors;
- Berechnung der negativen logarithmischen marginalen Likelihood (NLML);
- Berechnung der NLML-Gradienten für die Optimierung.
Zur Berechnung der erforderlichen Gradienten und der NLML selbst werden die inverse Kovarianzmatrix (K+σ2*I)^−1 und der Vektor α=(K+σ2*I)^−1 * y benötigt.
Die Funktion cho_solve hilft dabei, diese Variablen so schnell wie möglich zu berechnen.
Klasse LaplaceInference
LaplaceInference implementiert die Laplace-Approximation, eine approximative Inferenzmethode, die sowohl auf Regressions- als auch auf Klassifikationsprobleme in GP anwendbar ist.
//+------------------------------------------------------------------+ //| LaplaceInference | //+------------------------------------------------------------------+ class LaplaceInference : public IInference { private: int m_max_iterations; double m_tolerance; public: LaplaceInference(int max_iter = 100, double tolerance = 1e-10) : m_max_iterations(max_iter), m_tolerance(tolerance) {} void Infer(const matrix &X_train, const vector &y, IKernel *kernel, ILikelihood *likelihood, GPInferenceResult &result) override { result.success = false; int n = (int)y.Size(); // Calculate K matrix K = kernel.Compute(X_train, X_train); double jitter = 1e-6; K = K + matrix::Identity(n, n) * jitter; vector f = (result.mu_f_train.Size() == n) ? result.mu_f_train : vector::Zeros(n); bool converged = false; matrix W(n, n); // W := -Hessian matrix L_B(n, n); // L := cholesky(I + W^0.5*K*W^0.5), B = I + W^0.5*K*W^0.5 vector b(n); // b := Wf + grad loglike p(y|f) vector a(n); // a := b - W^0.5*L^T\(L\(W^0.5*K*b)) matrix sW = matrix::Zeros(n, n); // W^0.5 vector sW_diag(n); // sqrt(W) vector for storing diagonal elements double prev_lml_value = -DBL_MAX; // For the convergence criterion according to LML double current_lml_value = -DBL_MAX; int iter; // --- Step 1: Finding the f_hat mode using Newton's method (According to 3.1 GPML Algorithm) --- for(iter = 0; iter < m_max_iterations; iter++) { // 1. W := - Hessian (diagonal matrix) W = -1 * likelihood.LogLikelihoodHessian(f, y); // Calculate sW_diag_vector (root of W) sW_diag = MathSqrt(W.Diag()); // 2. L_B := cholesky(I + W^0.5*K*W^0.5) // Calculate sW_K = sW * K result.sW_K.Resize(n, n); for(int i = 0; i < n; i++) { result.sW_K.Row(sW_diag[i] * K.Row(i),i); // K[i,:] * sW_diag_vector[i] } // Calculate B = I + sW_K * sW matrix temp_sWKsW(n, n); for(int j = 0; j < n; j++) { temp_sWKsW.Col(result.sW_K.Col(j)*sW_diag[j],j) ; // sW_K[:,j] * sW_diag_vector[j] } matrix B = matrix::Identity(n, n) + temp_sWKsW; if(!B.Cholesky(L_B)) { PrintFormat("Error:Cholesky decomposition of B failed. Error Code: %d", GetLastError()); result.nlml_value = DBL_MAX; return; } result.L_B = L_B; // Save for use in forecasts // 3. b := W*f + grad loglike p(y|f) b = W.Diag() * f + likelihood.LogLikelihoodGradient(f, y); // 4. a := b - W^0.5*L_B^T\(L_B\(W^0.5*K*b)) a = b - sW_diag*cho_solve(L_B,result.sW_K @ b); // 5. f := K * a (new Newton step) f = K @ a; // --- Calculate the current NLML for the convergence criterion --- double log_likelihood_at_f = likelihood.LogLikelihood(f, y); double sum_log_diag_L_B = MathLog(L_B.Diag()).Sum(); current_lml_value = -0.5 * a @ f + log_likelihood_at_f - sum_log_diag_L_B; result.nlml_value = -current_lml_value; // 6. Convergence check double lml_change = current_lml_value - prev_lml_value; // PrintFormat("Iter %d: LML = %g, delta LML = %g",iter, current_lml_value, lml_change); // Check convergence by LML if((lml_change < m_tolerance)) { // PrintFormat("Converged at iter %d: LML = %g, delta LML = %g",iter, current_lml_value, lml_change); converged = true; break; } prev_lml_value = current_lml_value; } if(!converged) { PrintFormat("Warning: Newton algorithm didn't converge at iteration %d. Final LML: %g (delta: %g, tolerance: %g)", iter, current_lml_value, (current_lml_value - prev_lml_value), m_tolerance); } //---Save to the GPInferenceResult structure result.mu_f_train = f; // Posterior mean of the latent function at training points result.H = likelihood.LogLikelihoodHessian(f, y); // Hessian at point f_hat // --- Calculating analytical gradients NLML Algorithm 5.1 GPML --- result.nlml_gradient.Resize(kernel.GetNumHyperparameters() + likelihood.GetNumHyperparameters()); int current_grad_idx = 0; // ============================================================================== // 1. Gradients on KERNEL hyperparameters // ============================================================================== // Calculate R := W^0.5*L_B^-T*(L_B^-1*W^0.5) result.sW.Diag(sW_diag); matrix R(n,n); matrix cs = cho_solve(L_B,result.sW); for(int i = 0; i < n; i++) { R.Row(sW_diag[i] * cs.Row(i),i); } // Calculate C := L_B^-1 * W^0.5 * K // To do this, solve the linear equation: L_B * C = sW * K relative to C matrix C; if(!L_B.LinearEquationsSolutionTriangular(EQUATIONSFORM_N,result.sW_K, C)) { PrintFormat("Error: LinearEquationsSolutionTriangular L_B * C = W^0.5 * K.Error Code: %d", GetLastError()); result.nlml_value = DBL_MAX; return; } // Calculate A = Sigma_f_train = (K^−1+W)^−1 = K - C^T @ C // Sigma_f_train - posterior covariance matrix of the latent f function on the training data matrix CTC = C.Transpose() @ C; result.Sigma_f_train = K - CTC; //------------------- grad NLML = -(s1 + s2^T*s3) // Calculate s2 (first implicit part) // s2 := -0.5 * diag( diag(K) - diag(C^T*C) ) * ThirdDerivative loglike by f vector third_deriv = likelihood.LogLikelihoodThirdDerivative(f, y); //matrix diag; //diag.Diag(K.Diag() - CTC.Diag()); //vector s2 = -0.5 * diag @ third_deriv; vector s2 = -0.5 * (result.Sigma_f_train.Diag() * third_deriv); // First, we calculate the gradient of the log-likelihood vector log_lik_gradient = likelihood.LogLikelihoodGradient(f, y); int num_kernel_hyperparameters = kernel.GetNumHyperparameters(); for(int j = 0; j < num_kernel_hyperparameters; j++) { // C2 := dK_dtheta_j (derivative of the kernel matrix with respect to the current hyperparameter) matrix dK_dtheta_j = kernel.ComputeDerivative(j); ///------------------------------------------------------------------------------------- // s1 := 0.5*a^T*C2*a - 0.5 *trace (R*C2) // explicit part of the derivative double s1_term1 = 0.5 * (a @(dK_dtheta_j @ a)); // 0.5 * a^T * dK_dtheta_j * a // matrix R_dK_dtheta_j = R @ dK_dtheta_j; // double s1_term2 = -0.5 * R_dK_dtheta_j.Trace(); // -0.5 * trace(R * dK_dtheta_j) double s1_term2 = -0.5 * DiagonalTrace(R,dK_dtheta_j); // more efficient option double s1_explicit_part = s1_term1 + s1_term2; ///-------------------------------------------------------------------------------------- // b := C2 * grad loglike vector b = dK_dtheta_j @ log_lik_gradient; // s3 := b - K * R * b // second implicit part vector s3 = b - K @(R @ b); double implicit_part = s2 @ s3; // implicit part of the derivative // NLML result.nlml_gradient[current_grad_idx] = -(s1_explicit_part + implicit_part); current_grad_idx++; } // ============================================================================== // 2. Gradient for LIKELIHOOD hyperparameters dLogLikelihood/d(param_j) // ============================================================================== if(likelihood.GetNumHyperparameters() > 0) { vector likelihood_hyperparams = likelihood.GetHyperparameters(); for(int j = 0; j < (int)likelihood_hyperparams.Size(); j++) { // Term 1: Derivative of log p(y|f*) with respect to sigma double dlogpy_d_param_j = likelihood.LogLikelihoodGradientParam(f, y, j); // (the method returns dHessian/d(param_j) = d^3(log p)/d(param_j)d(f)^2) matrix dH_param_j = likelihood.LogLikelihoodHessianDerivative(f, y, j); // Convert dH_d_param_j to dW_d_param_j (where W = -H) matrix dW_d_param_j = -1.0 * dH_param_j; double term2_lik = -0.5*DiagonalTrace(result.Sigma_f_train,dW_d_param_j); result.nlml_gradient[current_grad_idx] = -(dlogpy_d_param_j + term2_lik); current_grad_idx++; } } result.success = true; } string GetName() const override { return "LaplaceInference"; }
Konstruktor LaplaceInference(int max_iter = 100, double tolerance = 1e-10): Der Klassenkonstruktor ermöglicht die Initialisierung der Parameter des iterativen Newton-Verfahrens, das zur Bestimmung des Modus der A-posteriori-Verteilung verwendet wird.
- max_iter – maximale Anzahl an Iterationen, nach denen der Algorithmus stoppt, auch wenn keine Konvergenz erreicht wurde.
- tolerance – Konvergenzniveau. Der Algorithmus beendet die Iteration, wenn die Änderung des NLML-Werts zwischen aufeinanderfolgenden Schritten unter diesen Schwellenwert fällt.
Die Methode Infer implementiert zwei Schlüsselalgorithmen aus dem Buch „Gaussian Processes for Machine Learning“ (GPML):
- Algorithmus 3.1 Bestimmung des Modus der A-posteriori-Verteilung und Berechnung der NLML,
- Algorithmus 5.1 Berechnung der NLML-Gradienten aus Hyperparametern.
Abb. 1. Algorithmus 3.1 zur Suche nach f_hat- und LML-Modi
- Der Prozess beginnt mit der Berechnung der K-Kern-Kovarianzmatrix. Für eine schnellere und stabilere Konvergenz wird das Ergebnis von mu_f_train aus der vorherigen Hyperparameter-Optimierungsiteration als Anfangsschätzwert für f verwendet, anstatt jedes Mal beim Nullvektor zu beginnen.
- Unter Verwendung des Newton-Verfahrens wird f_hat iterativ aktualisiert.
- Berechnung der NLML: Nachdem der Modus f konvergiert ist (wenn die Änderung der LML kleiner als m_tolerance wird), übergeben wir den NLML-Wert an den Optimierer.
Abb. 2. Algorithmus 5.1 Berechnung von LML-Gradienten
- Gradienten über Kernel-Hyperparameter beinhalten die Berechnung der Hilfsmatrizen R und C. Die R-Matrix wird unter Verwendung von cho_solve berechnet und die C-Matrix wird unter Verwendung von LinearEquationsSolutionTriangular berechnet.
- Die posteriore Kovarianzmatrix der latenten Funktion wird auf den Trainingsdaten Σf_train = K−C^TC berechnet.
- Der Gradient für jeden Kernel-Hyperparameter (dK_dtheta_j) besteht aus einem expliziten (s1) und impliziten Teilen (s2, s3). Die Berechnung von s1 wird durch die Verwendung von DiagonalTrace optimiert, und s3 beinhaltet ebenfalls cho_solve zum Lösen von Systemen.
- Gradienten in Bezug auf die Likelihood-Hyperparameter werden unter Verwendung der Ableitungen der Log-Likelihood und ihrer Hesse-Matrix in Bezug auf die entsprechenden Hyperparameter berechnet. Zu diesem Zweck werden die Methoden LogLikelihoodGradientParam und LogLikelihoodHessianDerivative des 'likelihood'-Objekts verwendet.
Nachdem wir nun alle grundlegenden Bausteine (Kernels, Likelihoods, Inferenzmethoden) haben, fahren wir damit fort, zu demonstrieren, wie die Bibliothek funktioniert.
Testen der Bibliothek an synthetischen Daten
Das Skript Gpsynthetic.mq5 hilft uns dabei, unsere Bibliothek an einfachen synthetischen Daten zu testen. Dies ermöglicht es uns, sicherzustellen, dass die implementierten Methoden funktionieren und korrekt sind.
// --- Include the GP library --- #include <GP/GP.mqh> enum IntervalType { INTERVAL_F = 0, // Confidence interval of the f* latent function INTERVAL_Y = 1 // Confidence interval of y* observations }; enum Type_inference { Exact = 0, // Exact Inference Laplace = 1 // Laplace Inference }; enum Type_Data { Regression = 0, // Regression Classification = 1 // Classification }; //--- Input parameters input IntervalType interval_type = INTERVAL_F; // Interval type to display input Type_Data DataType = Classification; input Type_inference inf = Laplace ; // Inference type //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { string InpFileName1; string InpFileName2; string InpFileName3; string InpFileName4; if(DataType == Regression) { // Regression data InpFileName1 = "Data_Regression/X_train.csv"; InpFileName2 = "Data_Regression/Y_train.csv"; InpFileName3 = "Data_Regression/X_test.csv"; InpFileName4 = "Data_Regression/Y_test.csv"; } else { InpFileName1 = "Data_Classification/X.csv"; // 200 InpFileName2 = "Data_Classification/y.csv"; InpFileName3 = "Data_Classification/X_star.csv"; InpFileName4 = "Data_Classification/y_star.csv"; } //----------- Dataset matrix x_train, y_train, x_test, y_test; CSVtoMatrix(InpFileName1,x_train); CSVtoMatrix(InpFileName2,y_train); CSVtoMatrix(InpFileName3,x_test); CSVtoMatrix(InpFileName4,y_test); // Create label vectors from matrices vector y_train_ = y_train.Col(0); vector y_test_ = y_test.Col(0); //--- 1. Create kernel objects IKernel* rbf = new RBFKernel(1,1); // IKernel* linear = new LinearKernel(1.0); // IKernel* periodic = new PeriodicKernel(1.0, 1.0, 5.8); //--- 2. Combination of kernels (creation of a compound kernel) // IKernel* functional_kernels_array[] = {rbf, linear, periodic}; // SumKernel* combined_functional_kernel = new SumKernel(functional_kernels_array); //--- 3. Creating a likelihood object ILikelihood* likelihood = NULL; // Declare a pointer of the base type if(DataType == Regression) { likelihood = new GaussianLikelihood(1); // Assign the object } if(DataType == Classification) { likelihood = new LogitLikelihood(); } //--- 4. Create the inference object IInference* inference; // Declare a pointer to the interface // Select the inference type: if(inf == Exact) { inference = new ExactInference(); } else { inference = new LaplaceInference(100, 1e-10); } //--- 5. Create an instance of the Gaussian Process model CREATE_GP_MODEL(gp_model, rbf, likelihood, inference,x_train, y_train_); // CREATE_GP_MODEL(gp_model, combined_functional_kernel, likelihood, inference,x_train, y_train_); //-------------------------------------------------- ulong start_time_fit = GetMicrosecondCount(); //--- 6. Hyperparameter optimization gp_model.Fit(); ulong end_time_fit = GetMicrosecondCount(); double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; Print(inference.GetName()); Print("Fit execution time: ", StringFormat("%.3f", elapsed_time_ms), " ms"); gp_model.PrintOptimizedKernelParameters(); //--- 7. Forecast for test points GPPredictionResult gp_predictions; // create a structure where the prediction results are to be set ulong start_time_predict = GetMicrosecondCount(); gp_model.Predict(x_test, gp_predictions,PROBIT); // MONTE_CARLO,NUM_INTEGR,PROBIT ulong end_time_predict = GetMicrosecondCount(); elapsed_time_ms = (end_time_predict - start_time_predict) / 1000.0; Print("Predict execution time: ", StringFormat("%.3f", elapsed_time_ms), " ms"); // --- 8. Classification results if(DataType == Classification) { DisplayClassificationResults(gp_predictions, y_test_); } //--- 9. Regression Results if(DataType == Regression) { VisualizeGP(x_train, y_train, x_test, y_test_, gp_predictions, gp_model, 15); } //--- 10. Releasing memory delete gp_model; }
Abhängig vom gewählten DataType (Regression oder Klassifikation) bestimmt das Skript die Pfade zu den CSV-Dateien mit Trainings- und Testdaten. Es wird davon ausgegangen, dass sich diese Dateien in den Ordnern Data_Regression oder Data_Classification unter Files befinden.
Als Nächstes werden Kernel-Objekte erstellt. Um eine komplexere Kovarianzfunktion zu erstellen, können kombinierte Kernels wie SumKernel oder ProductKernel verwendet werden.
Anschließend wird abhängig vom gewählten DataType ein Likelihood-Funktionsobjekt (ILikelihood) sowie ein Inferenzobjekt (IInference) erstellt.
Mithilfe des Makros CREATE_GP_MODEL (oder eines ähnlichen Konstruktors) wird ein GP-Modell erstellt, das den gewählten Kernel, die Likelihood-Funktion und die Inferenzmethode mit den Trainingsdaten verknüpft.
Die Methode gp_model.Fit(int maxiter = 20) startet den Modelltrainingsprozess. Die Methode verfügt nun über einen neuen Parameter, der für die Einstellung der maximalen Anzahl von maxiter-Iterationen zuständig ist.
Nach Abschluss des Trainings wird die Methode gp_model.Predict() verwendet, um Vorhersagen für neue (Test-)Daten x_test durchzuführen. Die Vorhersageergebnisse (Mittelwerte, Kovarianzen, Wahrscheinlichkeiten/Labels) werden in der Struktur GPPredictionResult gespeichert.
Für die Klassifikation können Sie den Vorhersagemodus (PROBIT, NUM_INTEGR, MONTE_CARLO) auswählen.
Das Skript verarbeitet und zeigt anschließend die Vorhersageergebnisse an:
- für die Klassifikation wird die Funktion DisplayClassificationResults() aufgerufen, die die Genauigkeitsmetrik und die vorhergesagten Wahrscheinlichkeiten berechnet;
- für Regression ist die Funktion VisualizeGP() für die Darstellung der Ergebnisse zuständig.
Für die Regression werden synthetische Daten verwendet, die durch die Funktion y = sin(x) + 0.5 * x + noise(sigma = 0.1) erzeugt werden. Wir werden nicht näher auf dieses Problem eingehen, da wir es bereits in dem Artikel über Regression besprochen haben. Beachten Sie nur, wie stark sich die Trainingszeit verbessert hat, seit wir auf die Verwendung analytischer Gradienten und Methoden aus der OpenBLAS-Bibliothek umgestiegen sind.
Für die binäre Klassifikation werden Daten verwendet, die das XOR-Problem (exklusives ODER) darstellen. Dies ist ein klassisches Problem des maschinellen Lernens, das die Grenzen linearer Modelle und die Notwendigkeit der Verwendung nichtlinearer Ansätze, wie z. B. GP, aufzeigt.
Erinnern Sie sich daran, dass XOR eine logische Operation ist, die zwei binäre Eingaben (0 oder 1) entgegennimmt und 1 zurückgibt, wenn genau eine der Eingaben 1 ist, und andernfalls 0. XOR-Wahrheitstabelle:
Die Aufgabe der binären XOR-Klassifikation besteht darin, ein Modell zu erstellen, das bei zwei gegebenen Eingaben (A, B) deren XOR-Ausgabe (0 und 1) vorhersagt. Da unser Modell jedoch nur mit den Labels +1 und -1 arbeitet, ersetzen wir das Label 0 durch -1.
Die Daten für das XOR-Problem werden wie folgt generiert. Zweidimensionale Eingabepunkte X (X=[x1,x2]) werden generiert, die zufällig über einen Bereich verteilt sind (z. B. [−4,4] auf jeder Achse). Die y-Labels für diese Punkte werden basierend auf der XOR-Logik bestimmt:
- Wenn x1 und x2 die gleichen Vorzeichen haben (beide positiv oder beide negativ), dann ist ihr Produkt x1⋅x2 positiv. In diesem Fall wird das Label +1 zugewiesen.
- Wenn x1 und x2 unterschiedliche Vorzeichen haben (eines ist positiv, das andere ist negativ), dann ist ihr Produkt x1⋅x2 negativ. In diesem Fall wird das Label −1 zugewiesen.
Somit gehören die Punkte (+,+) und (−,−) zur Klasse +1, und die Punkte (+,−) und (−,+) gehören zur Klasse −1. Für diese Aufgabe verwenden wir nur den RBF-Kernel, da er gut mit nichtlinearen Abhängigkeiten zurechtkommt, wodurch Gauß-Prozesse solche Klassen effektiv trennen können.
Während der Tests habe ich 200 Beobachtungen für das Training und 100 neue Punkte zur Überprüfung der Genauigkeit verwendet. Die Vorhersagegenauigkeit betrug 98 %, was die hohe Effizienz des GP bei der Lösung nichtlinearer Klassifizierungsprobleme demonstriert.
GPRegressor-Indikator: Zeitreihenregression mit GP
Der Indikator ist ein Beispiel für die dynamische Prognose von Finanzzeitreihen unter Berücksichtigung von Unsicherheit. Durch die Verwendung eines gleitenden Fensters zur Bildung des Trainingsdatensatzes und das erneute Training des Modells wird das Prinzip der Anpassung an sich ändernde Marktbedingungen umgesetzt.
//+------------------------------------------------------------------+ //| GPRegressor.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Include the GP library --- #include <GP/GP.mqh> #property indicator_chart_window #property indicator_buffers 3 #property indicator_plots 2 // plot Predicted mean #property indicator_label1 "Predicted mean" #property indicator_type1 DRAW_LINE #property indicator_color1 clrBlue #property indicator_style1 STYLE_SOLID #property indicator_width1 2 //--- plot Confidence interval #property indicator_label2 "Confidence interval" #property indicator_type2 DRAW_FILLING #property indicator_color2 clrLightGray #property indicator_width2 1 //--- Input parameters input int WindowLength = 10; // Window length for training the model input int calculate_bars = 100; // History depth for calculation input double zscore = 1.96; // z value for the interval (1.96 = 95%) input bool ShowRMSE = false; // Show RMSE metric //--- Buffers variables for rendering double ExtPredictionBuffer[]; // Buffer for predicted values double ExtUpperBandBuffer[]; // Buffer for the upper border of the confidence interval double ExtLowerBandBuffer[]; // Buffer for the lower border //--- Global objects GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Global matrix for X_train temporary indices matrix g_X_train; //--- Global variables for RMSE data accumulation vector gp_pred; vector naive_pred; vector y_true; // Flag indicating whether the RMSE has already been calculated (for a one-time calculation) bool rmse_calculated = false; string comment; double gp_rmse,naive_rmse; //+------------------------------------------------------------------+ //| Custom indicator initialization function | //+------------------------------------------------------------------+ int OnInit() { if(Bars(Symbol(), Period()) < WindowLength + calculate_bars) { PrintFormat("Error: Not enough bars to calculate. Minimum %d required. Available: %d", WindowLength + calculate_bars, Bars(Symbol(), Period())); return(INIT_FAILED); } rmse_calculated = false; // Reset the flag on initialization/reinitialization //--- Initialize global vectors for RMSE gp_pred.Resize(calculate_bars-1); naive_pred.Resize(calculate_bars-1); y_true.Resize(calculate_bars-1); //--- indicator buffers mapping SetIndexBuffer(0, ExtPredictionBuffer, INDICATOR_DATA); SetIndexBuffer(1, ExtUpperBandBuffer, INDICATOR_DATA); SetIndexBuffer(2, ExtLowerBandBuffer, INDICATOR_DATA); //--- Create kernel, likelihood, and inference method objects g_kernel = new RBFKernel(1.0, 1.0); if(g_kernel == NULL) { Print("Failed to create RBFKernel object"); return(INIT_FAILED); } g_likelihood = new GaussianLikelihood(0.1); if(g_likelihood == NULL) { Print("Failed to create GaussianLikelihood object"); delete g_kernel; return(INIT_FAILED); } g_inference = new ExactInference(); if(g_inference == NULL) { Print("Failed to create ExactInference object"); delete g_kernel; delete g_likelihood; return(INIT_FAILED); } //--- Initialize the feature matrix X_train (one feature - time) g_X_train = matrix::Zeros(WindowLength, 1); for(int i = 0; i < WindowLength; i++) { g_X_train[i][0] = (double)(i + 1); // X = [1, 2, ..., WindowLength] } //--- Create the GaussianProcess object // Pass g_X_train and empty y_train vector initial_y_train = vector::Zeros(WindowLength); g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference, g_X_train, initial_y_train); if(g_gp_model == NULL) { Print("Failed to create GaussianProcess object"); delete g_kernel; delete g_likelihood; delete g_inference; return(INIT_FAILED); } return(INIT_SUCCEEDED); } //+------------------------------------------------------------------+ //| | //+------------------------------------------------------------------+ int OnCalculate(const int32_t rates_total, const int32_t prev_calculated, const datetime &time[], const double &open[], const double &high[], const double &low[], const double &close[], const long &tick_volume[], const long &volume[], const int32_t &spread[]) { // Check for sufficient number of bars if(rates_total < WindowLength) { Print("Error: Not enough bars to calculate. At least ", WindowLength, " is required, available: ", rates_total); return(0); } int start; if(prev_calculated == 0) { // Determine where to start the calculation. // calculate_bars - history depth to calculate. // MathMax(WindowLength, ...) ensures that we start no sooner than there is a full window of data. start = MathMax(WindowLength, rates_total - calculate_bars); // Initialize all buffers with EMPTY_VALUE ArrayInitialize(ExtPredictionBuffer, EMPTY_VALUE); ArrayInitialize(ExtUpperBandBuffer, EMPTY_VALUE); ArrayInitialize(ExtLowerBandBuffer, EMPTY_VALUE); // Set the start of rendering buffers PlotIndexSetInteger(0, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(1, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(2, PLOT_DRAW_BEGIN, start); } else { if(time[rates_total - 1] == time[prev_calculated - 1]) { // This means that a new bar has NOT appeared. // We are on the same bar, just new ticks have arrived. return(rates_total); } // If the time of the last bar has changed, then a new closed bar has appeared. // Start calculation from prev_calculated to handle all new bars. start = prev_calculated; } for(int i = start; i < rates_total && !IsStopped(); i++) { // Check if there are enough bars to form a window if(i < WindowLength) { Print("Error: Not enough bars to form a window on bar ", i, ". Required: ", WindowLength, ", available: ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Skip this bar } // Form y vector from close[i-1], close[i-2], ..., close[i-WindowLength] // close[i-1] - the newest bar in the window // close[i-WindowLength] - the oldest bar in the window vector y(WindowLength); for(int j = 0; j < WindowLength; j++) { y[WindowLength - 1 - j] = close[i - 1 - j]; } //Print(" y = ", y[WindowLength-1]); // The newest bar in the window //--- Data standardization double mu_raw = y.Mean(); // Print("mu_raw = ", mu_raw); double sigma_raw = y.Std(); // Print("sigma_raw = ", sigma_raw); vector y_standardized = (y - mu_raw) / sigma_raw; //--- Set up training data g_gp_model.SetTrainingData(g_X_train, y_standardized); //--- Reset initial hyperparameters before Fit() vector initial_params(3); initial_params[0] = 1.0; // sigma_f (RBFKernel) initial_params[1] = 1.0; // length_scale (RBFKernel) initial_params[2] = 0.1; // sigma_n (GaussianLikelihood) g_gp_model.SetHyperparameters(initial_params); // vector params = g_gp_model.GetCurrentHyperparameters(); // Print(params); // ulong start_time_fit = GetMicrosecondCount(); //--- Train the model if(!g_gp_model.Fit()) { Print("Error: Failed to optimize GP hyperparameters for bar ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; // Skip the current bar } //ulong end_time_fit = GetMicrosecondCount(); // double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; //Print("Fit execution time: ", StringFormat("%.3f", elapsed_time_ms), " мс"); vector params = g_gp_model.GetCurrentHyperparameters(); double sigma_f = params[0]; // sigma_f (RBFKernel) double length_scale = params[1]; // length_scale (RBFKernel) double sigma_n = params[2]; // sigma_n (GaussianLikelihood) // --- Create a point for forecasting (X_star) // Predict the value for the next step after WindowLength matrix X_star = matrix::Zeros(1, 1); X_star[0][0] = (double)(WindowLength + 1); // --- Perform the prediction GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Error: Failed to get forecast from GP model for bar ", i); ExtPredictionBuffer[i] = EMPTY_VALUE; ExtUpperBandBuffer[i] = EMPTY_VALUE; ExtLowerBandBuffer[i] = EMPTY_VALUE; continue; } // --- double predicted_mean_st = gp_predictions.mu_f_star[0]; // --- Perform the inverse transformation of the predicted mean value // back to the original price scale double predicted_mean = predicted_mean_st * sigma_raw + mu_raw; double sigma_f_star = gp_predictions.Sigma_f_star[0,0]; // dispersion (Σ) of f* signal double predicted_variance_st = gp_predictions.Sigma_y_star[0,0]; // dispersion(Σ) of y* observation // --- Perform the inverse dispersion transformation // The variance is scaled by the square of the standard deviation (sigma_raw) of the data, // used for standardization. double predicted_variance = predicted_variance_st * sigma_raw * sigma_raw; // Calculate the predicted confidence interval double std_predict = MathSqrt(predicted_variance); double upper_bound = predicted_mean + zscore * std_predict; double lower_bound = predicted_mean - zscore * std_predict; // --- Save the results to the indicator buffers for rendering ExtPredictionBuffer[i] = predicted_mean; ExtUpperBandBuffer[i] = upper_bound; ExtLowerBandBuffer[i] = lower_bound; comment = StringFormat("GP Regressor | Mean=%.5f | Upper=%.5f | Lower=%.5f |\n", predicted_mean, upper_bound, lower_bound); comment+= StringFormat("Sigma_f = %.5f | Length_Scale = %.5f | Sigma_n = %.5f\n",sigma_f,length_scale,sigma_n); comment+= StringFormat("mu f* = %.5f | Σ y* = %.5f | Σ f*= %.5f",predicted_mean_st,predicted_variance_st,sigma_f_star); Comment(comment); // Debug output Print("Bar forecast ", i, ": Mean=", DoubleToString(predicted_mean, Digits()), " Upper=", DoubleToString(upper_bound, Digits()), " Lower=", DoubleToString(lower_bound, Digits())); } // --- Calculate RMSE to evaluate the efficiency of the GP forecast model compared to a simple "naive" forecast if(ShowRMSE && !rmse_calculated) { vector returns(calculate_bars-1); for(int j=0; j <calculate_bars-1;j++) { // true value y_true[j] = close[rates_total - calculate_bars + j]; //Print("y_true[j]", y_true[j] ); // GP forecast: gp_pred[j] = ExtPredictionBuffer[rates_total-calculate_bars + j]; //Print("gp_pred[j]", gp_pred[j] ); // "Naive" forecast (price tomorrow = price today): naive_pred[j] = close[rates_total - calculate_bars + j-1]; // Print("naive_pred[j]", naive_pred[j] ); returns[j] = close[rates_total - calculate_bars + j] - close[rates_total - calculate_bars + j-1]; } gp_rmse=gp_pred.RegressionMetric(y_true,REGRESSION_RMSE); naive_rmse=naive_pred.RegressionMetric(y_true,REGRESSION_RMSE); double std_returns = returns.Std(); // standard deviation of Close increments comment = comment + StringFormat("\nRMSE GP: %.5f | RMSE Naive: %.5f | StdDev Returns Close: %.5f ", gp_rmse, naive_rmse, std_returns); rmse_calculated = true; } Comment(comment); //--- Return rates_total for the next call return(rates_total); } //+------------------------------------------------------------------+ //| Custom indicator deinitialization function | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { //--- Clear the comment on the chart Comment(""); //--- Free up the allocated memory for the GP object if(g_gp_model != NULL) { delete g_gp_model; } } //+------------------------------------------------------------------+
Eingaben:
- WindowLength – Länge des Datenfensters (in Bars) für das Training des GP-Modells.
- calculate_bars – Historientiefe für die Berechnung und Zeichnung des Indikators.
- zscore – Anzahl der Standardabweichungen zur Konstruktion des Konfidenzintervalls (z. B. 1,96 für 95 %).
- ShowRMSE – Flag zur Anzeige der Quadratwurzel des mittleren quadratischen Fehlers (RMSE).
Initialisierung (OnInit):
- Globale Objekte. Globale Kernel-Objekte (RBFKernel), Likelihood-Funktionen (GaussianLikelihood) und Inferenzmethoden (ExactInference) werden erstellt und initialisiert.
- X_train-Merkmale. GPRegressor implementiert die Regression über die Zeit unter Verwendung von Zeitindizes (von 1 bis WindowLength) innerhalb des Trainingsfensters als X-Eingaben.
- GaussianProcess-Modell. Das Objekt g_gp_model wird mit den angegebenen Komponenten initialisiert, während y_train anfangs leer ist und dynamisch gefüllt werden soll.
Hauptberechnung (OnCalculate):
- Datengenerierung: An jeder Bar wird die WindowLength der neuesten Schlusskurse extrahiert, um den y-Vektor zu generieren.
- Datenstandardisierung: y wird standardisiert (auf Nullmittelwert und Einheitsvarianz reduziert), um die numerische Stabilität der GP-Operation zu gewährleisten.
- Aktualisierung und Training des Modells:
- Hyperparameter werden vor jedem Training auf die Anfangswerte zurückgesetzt (g_gp_model.SetHyperparameters)
- Das GP-Modell wird auf dem aktuellen Fenster der standardisierten Preise g_gp_model.Fit() trainiert.
- Vorhersage: Um die nächste Bar vorherzusagen, wird ein X_star-Punkt mit einem Zeitindex von WindowLength + 1 erstellt. Die Methode g_gp_model.Predict() führt die Vorhersage durch und gibt den Mittelwert (mu_f_star) und die Varianz (Sigma_y_star) für die standardisierten Daten zurück.
- Destandardisierung: Der vorhergesagte Mittelwert und die Varianz werden unter Verwendung des Mittelwerts (mu_raw) und der Standardabweichung (sigma_raw) des aktuellen Fensters auf die ursprüngliche Preisskala zurücktransformiert.
- Konfidenzintervall: Die oberen und unteren Grenzen des Intervalls werden basierend auf dem transformierten Mittelwert und der Standardabweichung unter Verwendung eines vom Benutzer vorgegebenen zscore-Werts berechnet.
Abb. 3. Regressionsmodellvorhersage und 95%-Konfidenzintervall
Berechnung der RMSE-Metrik:
Der Indikator zeigt RMSE-Statistiken für die GP-Modellprognosen und die naive Prognose (wobei der Preis von morgen = Preis von heute ist) für die letzten calculate_bars einmal pro Benutzeranfrage an. Dies hilft dabei, die aktuelle Modelleffizienz objektiv zu bewerten.
Beispielsweise könnten typische Indikatoren für EURUSD auf einem Minuten-Zeitrahmen sein: RMSE GP (0,00032) | RMSE Naive (0,00029), und die Standardabweichung der Schlusskursänderungen (StdDev Returns Close) beträgt 0,00029.
Die Tatsache, dass der GP-RMSE höher ist als der naive RMSE, deutet darauf hin, dass das GP-Modell im gegebenen analysierten Zeitraum schlechter abschneidet als die einfache naive Prognose. Darüber hinaus ist die Gleichheit von RMSE Naive und StdDev Returns Close ein charakteristisches Merkmal eines Random Walk, bei dem zukünftige Preisänderungen unabhängig von vergangenen sind und der aktuelle Wert die beste Prognose darstellt.
Daraus lässt sich schließen, dass es in diesem Abschnitt der Zeitreihe entweder keine vorhersagbaren Muster gibt oder die aktuelle Konfiguration des GP-Modells (das nur die Zeit als Merkmal verwendet) noch nicht in der Lage ist, nützliche Informationen zu extrahieren, um die naive Prognose zu übertreffen.
GPClassifier-Indikator: Klassifizierung basierend auf Gauß-Prozessen
Der GPClassifier-Indikator demonstriert den Prozess des Erstellens, Trainierens und Vorhersagens eines GP-Modells bei einem binären Klassifizierungsproblem. Die verwendeten Merkmale sind Änderungen der Schlusskurse, die in einem gleitenden Fenster historischer Daten berechnet werden. Die Funktion NormalizeData wird verwendet, um die Merkmalsmatrix zu normalisieren und das Modelltraining zu stabilisieren.
//+------------------------------------------------------------------+ //| GPClassifier.mq5 | //| Eugene | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Eugene" #property link "https://www.mql5.com" #property version "1.00" // --- Include the GP library --- #include <GP/GP.mqh> #property indicator_chart_window #property indicator_buffers 4 #property indicator_plots 2 // plot for predicted class "+1" (up arrows) #property indicator_label1 "Predicted Class +1" #property indicator_type1 DRAW_ARROW #property indicator_color1 clrGreen #property indicator_width1 1 // plot for predicted class "-1" (down arrows) #property indicator_label2 "Predicted Class -1" #property indicator_type2 DRAW_ARROW #property indicator_color2 clrBlack #property indicator_width2 1 //--- Input parameters input int WindowLength = 10; // Window length for training the GP model input int NumLags = 1; // Number of features input int calculate_bars = 100; // History depth for calculation input double ProbabilityThreshold = 0.5; // Probability threshold input int ArrowShiftPx = 10; // Arrow offset in pixels input bool ShowACCURACY = false; // Show accuracy on the chart //--- Buffers variables for rendering double ExtClassUpBuffer[]; // Buffer for displaying UP arrows double ExtClassDownBuffer[]; // Buffer for displaying DOWN arrows double ExtPredictedClassBuffer[]; // Buffer for storing predicted class labels double ExtPredictedProbabilityBuffer[]; // Buffer for storing predicted probabilities //--- Global Gaussian process objects GaussianProcess* g_gp_model = NULL; IKernel* g_kernel = NULL; ILikelihood* g_likelihood = NULL; IInference* g_inference = NULL; //--- Global variables for classification metrics bool accuracy_calculated = false; vector gp_pred,naive_pred,y_true,gp_accuracy,naive_accuracy; string comment; int currentsize; //+------------------------------------------------------------------+ //| Custom indicator initialization function | //+------------------------------------------------------------------+ int OnInit() { accuracy_calculated = false; // Reset the flag on initialization/reinitialization //--- Initialize global vectors for Accuracy gp_pred.Resize(0); naive_pred.Resize(0); y_true.Resize(0); gp_accuracy.Resize(1); naive_accuracy.Resize(1); currentsize = 0; //--- indicator buffers mapping SetIndexBuffer(0, ExtClassUpBuffer, INDICATOR_DATA); // For up arrows SetIndexBuffer(1, ExtClassDownBuffer, INDICATOR_DATA); // For down arrows SetIndexBuffer(2, ExtPredictedClassBuffer, INDICATOR_CALCULATIONS); // to store predicted labels SetIndexBuffer(3, ExtPredictedProbabilityBuffer, INDICATOR_CALCULATIONS); // to store probabilities // --- Set arrow symbols PlotIndexSetInteger(0, PLOT_ARROW, 225); // Up arrow for buffer 0 PlotIndexSetInteger(1, PLOT_ARROW, 226); // Down arrow for buffer 1 // --- Set the arrow offset in pixels --- // For the up arrows (buffer 0), use a negative offset so they are above High PlotIndexSetInteger(0, PLOT_ARROW_SHIFT, -ArrowShiftPx); // negative = up // For the down arrows (buffer 1), use a positive offset so they are below Low PlotIndexSetInteger(1, PLOT_ARROW_SHIFT, ArrowShiftPx); // positive = down // Add a check for the ratio of NumLags and WindowLength if(NumLags <= 0) { Print("Error: NumLags must be positive number"); return(INIT_FAILED); } if(NumLags >= WindowLength) { Print("Error: NumLags (", NumLags, ") cannot exceed or be equal to WindowLength (", WindowLength, ")."); return(INIT_FAILED); } //--- Create kernel, likelihood, and inference method objects g_kernel = new RBFKernel(1.0, 1.0); if(g_kernel == NULL) { Print("Failed to create RBFKernel object"); return(INIT_FAILED); } g_likelihood = new LogitLikelihood(); if(g_likelihood == NULL) { Print("Failed to create LogitLikelihood object"); delete g_kernel; return(INIT_FAILED); } g_inference = new LaplaceInference(); if(g_inference == NULL) { Print("Failed to create ExactInference object"); delete g_kernel; delete g_likelihood; return(INIT_FAILED); } // --- Create the GaussianProcess model // Pass empty matrices for X_train and y_train. // They will be overridden in OnCalculate. g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference,matrix::Zeros(1,1), vector::Zeros(1)); if(g_gp_model == NULL) { Print("Failed to create GaussianProcess object"); delete g_kernel; delete g_likelihood; delete g_inference; return(INIT_FAILED); } return(INIT_SUCCEEDED); } //+------------------------------------------------------------------+ //| Custom indicator iteration function | //+------------------------------------------------------------------+ int OnCalculate(const int32_t rates_total, const int32_t prev_calculated, const datetime &time[], const double &open[], const double &high[], const double &low[], const double &close[], const long &tick_volume[], const long &volume[], const int32_t &spread[]) { // Check for sufficient number of bars if(rates_total < WindowLength + NumLags + 1) { Print("Error: Not enough bars to calculate. Required at least ", WindowLength + NumLags + 1, ", available: ", rates_total); return(0); } int start; if(prev_calculated == 0) { start = MathMax(WindowLength + NumLags + 1, rates_total - calculate_bars); // Initialize all EMPTY_VALUE buffers ArrayInitialize(ExtClassUpBuffer, EMPTY_VALUE); ArrayInitialize(ExtClassDownBuffer, EMPTY_VALUE); // Number of initial bars without rendering and values in DataWindow PlotIndexSetInteger(0, PLOT_DRAW_BEGIN, start); PlotIndexSetInteger(1, PLOT_DRAW_BEGIN, start); } else { if(time[rates_total - 1] == time[prev_calculated - 1]) { return(rates_total); } start = prev_calculated; } // --- Main forecast calculation loop for(int i = start; i < rates_total && !IsStopped(); i++) { // Check if there are enough bars to form X and y window if(i < WindowLength + NumLags + 1) { ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } // --- Step 1: Generate training data (X_train and y_train) matrix X_train = matrix::Zeros(WindowLength, NumLags); vector y_train = vector::Zeros(WindowLength); for(int j = 0; j < WindowLength; j++) { // Calculate the bar index in the 'close' array for the current observation in the window. int idx = i - WindowLength + j; // Formation of features (lags) for the current observation (X_train[j] strings) for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { // Example: // lag_idx = 0: (close[idx - 1] - close[idx - 2]) - increment of the previous bar // lag_idx = 1: (close[idx - 2] - close[idx - 3]) - increment of the bar before last X_train[j,lag_idx] = close[idx - 1 - lag_idx] - close[idx - 2 - lag_idx]; } // Y_train[j]: Binary label for price increment on idx bar. if(close[idx] - close[idx - 1] > 0) { y_train[j] = 1; // Price increased } else { y_train[j] = -1; // Price decreased or remained unchanged } } // Print(y_train); // --- Normalize X_train matrix X_train_norm; vector out_mean; vector out_std; if(!NormalizeData(X_train,X_train_norm,out_mean,out_std)) { Print("Error: Failed to normalize feature matrix ", i); continue; } // --- Step 2: Feed data into the GP model g_gp_model.SetTrainingData(X_train_norm, y_train); // Reset initial hyperparameters before each Fit() training vector initial_params(2); initial_params[0] = 1.0; // sigma_f (RBFKernel) initial_params[1] = 1.0; // length_scale (RBFKernel) g_gp_model.SetHyperparameters(initial_params); // --- Step 3: Train the model // ulong start_time_fit = GetMicrosecondCount(); if(!g_gp_model.Fit()) { Print("Error: Failed to optimize GP hyperparameters for bar ", i); ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } // ulong end_time_fit = GetMicrosecondCount(); // double elapsed_time_ms = (end_time_fit - start_time_fit) / 1000.0; // Print("Fit execution time: ", StringFormat("%.3f", elapsed_time_ms), " ms"); // --- Step 4: Create a point for the X_star forecast // X_star - increment of prices of previous NumLags closed bars, // used to predict the movement of the current bar (i bar). matrix X_star = matrix::Zeros(1, NumLags); for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { X_star[0][lag_idx] = close[i - 1 - lag_idx] - close[i - 2 - lag_idx]; } // X_star should also be normalized using the same means and standard deviations, // as for the corresponding X_train columns. for(int lag_idx = 0; lag_idx < NumLags; lag_idx++) { X_star[0][lag_idx] = (X_star[0][lag_idx] - out_mean[lag_idx]) / out_std[lag_idx]; } // --- Step 5: Perform the prediction GPPredictionResult gp_predictions; if(!g_gp_model.Predict(X_star, gp_predictions)) { Print("Error: Failed to get forecast from GP model for bar ", i); ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; continue; } double predicted_probability = gp_predictions.predicted_probabilities[0]; // Probability for a single test point int predicted_class = (int)gp_predictions.predicted_labels[0]; // Label (+1 or -1) // --- Step 6: Save the results to the indicator buffers for rendering ExtClassUpBuffer[i] = EMPTY_VALUE; ExtClassDownBuffer[i] = EMPTY_VALUE; ExtPredictedClassBuffer[i] = predicted_class; ExtPredictedProbabilityBuffer[i] = predicted_probability; // Add the probability output and the number of features (lags) as a chart comment comment = StringFormat("GP Classifier (Lags: %d) | Probability(UP) = %.10f", NumLags, predicted_probability); Comment(comment); // If the predicted class is "+1" and the probability is above the given threshold if(predicted_class == 1 && predicted_probability > ProbabilityThreshold) { ExtClassUpBuffer[i] = high[i]; } // If the predicted class is "-1" and the probability is above the given threshold else if(predicted_class == -1 && (1.0 - predicted_probability) > ProbabilityThreshold) { ExtClassDownBuffer[i] = low[i]; } Print("Bar forecast ", i, ": Probability(UP)=" + DoubleToString(predicted_probability, 10) + ", Predicted Class=" + IntegerToString(predicted_class), " time: ", time[i]); } // --- Calculate classification accuracy for assessing the efficiency of the GP predictive model in comparison //with a "naive" forecast. if(ShowACCURACY && !accuracy_calculated) { for(int j=0; j <calculate_bars-1;j++) { // index to the current bar int bar_idx = rates_total - calculate_bars + j; if(ExtClassUpBuffer[bar_idx]!= EMPTY_VALUE || ExtClassDownBuffer[bar_idx] != EMPTY_VALUE) { currentsize = (int)gp_pred.Size(); currentsize++; gp_pred.Resize(currentsize,100); naive_pred.Resize(currentsize,100); y_true.Resize(currentsize,100); // True label: count to the current bar if(close[bar_idx] - close[bar_idx - 1] > 0) { y_true[currentsize-1] = 1; // Price increased } else { y_true[currentsize-1] = -1; // Price decreased or remained unchanged } // Print("y_true", y_true[currentsize-1] ); // GP forecast: count to the current bar gp_pred[currentsize-1] = ExtPredictedClassBuffer[bar_idx]; // Print(" gp_pred", gp_pred[currentsize-1]); // "Naive" forecast: the "for tomorrow" bar label is equal to the "today" label if(close[bar_idx - 1] - close[bar_idx - 2] > 0) { naive_pred[currentsize-1] = 1; } else { naive_pred[currentsize-1] = -1; } } } if(!(int)gp_pred.Size()==0) { gp_accuracy=gp_pred.ClassificationMetric(y_true,CLASSIFICATION_ACCURACY); naive_accuracy=naive_pred.ClassificationMetric(y_true,CLASSIFICATION_ACCURACY); comment += StringFormat("\nGP Accuracy: %.2f %% | Naive Accuracy: %.2f %%", gp_accuracy[0], naive_accuracy[0]); comment += StringFormat("\nNumber of Filtered Signals: %d", (int)gp_pred.Size()); } else { comment += StringFormat("\n No Signals for ProbabilityThreshold: %.4f ", ProbabilityThreshold); } accuracy_calculated = true; } Comment(comment); //--- Return rates_total for the next call return(rates_total); } //+------------------------------------------------------------------+ //| Custom indicator deinitialization function | //+------------------------------------------------------------------+ void OnDeinit(const int reason) { // Clear the comment on the chart when deleting the indicator Comment(""); if(g_gp_model != NULL) { delete g_gp_model; } }
Eingaben:
- WindowLength – Länge des Datenfensters für das Training des Modells. Anzahl der Beobachtungen in der Trainingsstichprobe.
- NumLags – Anzahl der Lags (Preisinkremente). Wenn beispielsweise NumLags = 1 ist, ist das einzige Merkmal das Preisinkrement der vorherigen Bar. Wenn NumLags = 3 ist, sind die Merkmale die Preisinkremente der drei vorherigen Bars usw.
- calculate_bars – Anzahl der letzten Bars im Chart, für die der Indikator berechnet und gezeichnet wird.
- ProbabilityThreshold – Wahrscheinlichkeitsschwellenwert. Der Indikator zeigt Signale mit Aufwärts-/Abwärtspfeilen an, wenn die vorhergesagte Wahrscheinlichkeit einer Preisbewegung in die entsprechende Richtung diesen Schwellenwert überschreitet.
- ArrowShiftPx – Versatz der Pfeile in Pixeln von den Bar-Extremen.
- ShowACCURACY – Flag zur Anzeige von Klassifizierungsgenauigkeitsmetriken in einem Chart.
OnInit() – Initialisierung des GP-Modells:
- Erstellen Sie einen RBF-Kernel, eine Likelihood-Funktion und eine Inferenzmethode:
- g_kernel = new RBFKernel(1.0, 1.0)
- g_likelihood = new LogitLikelihood()
- g_inference = new LaplaceInference()
- Erstellen Sie das GaussianProcess-Modell:
- g_gp_model = new GaussianProcess(g_kernel, g_likelihood, g_inference, matrix::Zeros(1,1), vector::Zeros(1));
- Das Objekt g_gp_model wird mit leeren Anfangsmatrizen für die Trainingsdaten X_train und y_train erstellt. Diese Daten werden dynamisch gefüllt und bei jeder Bar in der Funktion OnCalculate aktualisiert.
OnCalculate() – Hauptberechnungsschleife:
1. Überprüfung auf eine ausreichende Anzahl von Bars: Vor Beginn der Berechnungen stellen wir sicher, dass das Chart genügend Bars aufweist, um ein Trainingsfenster zu bilden, wobei WindowLength, NumLags und eine Bar für das vorhergesagte Label berücksichtigt werden.
2. Bildung der Trainingsdaten (X_train und y_train):
- X_train-Merkmale – bei jedem Schritt wird eine Matrix der X_train-Merkmale in einem gleitenden Fenster gebildet. Jede Zeile ist eine Beobachtung, die Spalten sind Preisänderungen (Lags) für die vorherigen NumLags-Bars.
- y_train-Labels – binäre Labels, die die Richtung der Preisbewegung der nächsten Bar angeben: „1“ (Preis gestiegen) oder „-1“ (Preis gefallen oder unverändert geblieben).
3. X_train-Normalisierung – die Merkmalsmatrix wird zur numerischen Stabilität standardisiert. Die in diesem Schritt ermittelten Mittelwerte (out_mean) und Standardabweichungen (out_std) werden gespeichert, um den Vorhersagepunkt zu normalisieren.
4. Aktualisierung und Training des GP-Modells:
- g_gp_model.SetTrainingData(g_X_train, y_standardized); – die Trainingsdaten werden an jeder Bar aktualisiert.
- g_gp_model.SetHyperparameters(initial_params) – die Kernel-Hyperparameter werden für jede neue Bar zurückgesetzt und neu optimiert.
- g_gp_model.Fit() – Start des Trainings auf dem aktuellen Fenster der standardisierten Daten.
5. Erstellung und Normalisierung eines neuen Punktes für die Vorhersage (X_star): Normalisierung unter Verwendung der Vektoren out_mean und out_std.
6. Durchführung der Vorhersage: Die Methode g_gp_model.Predict(X_star, gp_predictions) gibt die vorhergesagte Wahrscheinlichkeit einer Aufwärtsbewegung (predicted_probabilities[0]) und das binäre Label (predicted_labels[0]) zurück.
7. Zeichnen und Ausgabe: Abhängig von der predicted_class und dem ProbabilityThreshold wird ein Aufwärtspfeil oder ein Abwärtspfeil im Chart gezeichnet. Informationen über die Vorhersagewahrscheinlichkeit und die vorhergesagte Klasse werden in den Kommentaren im Chart und im Protokoll angezeigt.
Abb. 4. Vorhersageergebnis des GP-Klassifikator-Indikators (RBF-Kernel)
Berechnung der Genauigkeitsmetriken (ShowACCURACY)
Wenn der Parameter ShowACCURACY aktiviert ist, führt der Indikator eine einmalige Berechnung der Klassifizierungsgenauigkeit für die letzten calculate_bars durch und vergleicht die GP-Modellprognosen mit der naiven Prognose.
Wahre Labels (y_true) sind als die tatsächliche Richtung der Preisbewegung definiert. GP-Vorhersage (gp_pred) – die vom GP-Modell vorhergesagten Labels. Die naive Prognose (naive_pred) geht davon aus, dass die Richtung der nächsten Preisbewegung dieselbe sein wird wie die vorherige Bewegung (wenn der Preis beispielsweise an der vorherigen Bar gestiegen ist, ist die naive Prognose Wachstum).
Die CLASSIFICATION_ACCURACY-Metriken für beide Modelle (gp_accuracy und naive_accuracy) werden verglichen und die Ergebnisse in einem Chartkommentar angezeigt. Dies ermöglicht es uns, die Leistung des GP-Modells im Vergleich zu einem einfachen Benchmark objektiv zu bewerten.
Schlussfolgerung
Heute haben wir unsere Betrachtung von Gauß-Prozessen für Klassifikation und Regression abgeschlossen. Nach einer detaillierten Überprüfung der theoretischen Konstrukte konnten wir eine Basisversion der Bibliothek erstellen, die es uns ermöglicht, sowohl Regressions- als auch Klassifizierungsprobleme zu lösen. Tests mit synthetischen Daten bestätigten, dass die implementierten Komponenten korrekt funktionieren, einschließlich der effektiven Nutzung analytischer Gradienten und der Verwendung der leistungsstarken OpenBLAS-Bibliothek für numerische Stabilität und Geschwindigkeit.
Die Indikatoren GPRegressor und GPClassifier zeigten, dass sich GP praktisch für die dynamische Vorhersage von Finanzzeitreihen in Echtzeit einsetzen lässt, während Metriken wie RMSE (für Regression) und Genauigkeit (für Klassifizierung) eine objektive Bewertung der Modellleistung ermöglichten.
Die Bibliotheksimplementierung ist jedoch alles andere als ideal. Trotz aller unternommener Anstrengungen ist ihre Geschwindigkeit etablierten Lösungen wie scikit-learn unterlegen. Dies liegt teilweise daran, dass wir einen relativ langsamen MinBleic-Optimierer verwendet haben.
Daher könnte eine weitere Verbesserung der Bibliothek mit Folgendem zusammenhängen:
- der Verwendung eines effizienteren Gradientenoptimierers. Insbesondere benötigen wir einen schnellen L-BFGS-Algorithmus, der den aktuellen bei Problemen mit einer großen Anzahl von Hyperparametern in der Effizienz deutlich übertrifft;
- einer robusteren Implementierung des Newton-Verfahrens für die Laplace-Approximation. Das Hinzufügen einer linearen Suche ermöglicht ein effizienteres und zuverlässigeres Finden des Modus der A-posteriori-Verteilung, was die Konvergenz der endgültigen Optimierung verbessert;
- Implementierung von Sparse-Methoden. Im Gegensatz zu aktuellen Arbeiten mit dichten Matrizen ermöglicht die Verwendung dünnbesetzter Repräsentationen der Bibliothek das Training mit deutlich größeren Datenmengen, was neue Möglichkeiten für die Skalierung und Anwendung von GPs eröffnet.
Im Artikel verwendete Programme
| # | Name | Typ | Beschreibung |
|---|---|---|---|
| 1 | GPClassifier.mq5 | Indikator | Klassifiziert die Preisrichtung einen Schritt voraus, basierend auf vorhergesagten Wahrscheinlichkeiten (Online-Training) |
| 2 | GPRegressor.mq5 | Indikator | Prognostiziert Preis und Konfidenzintervall für den nächsten Schritt (Online-Training) |
| 3 | GPsynthetic.mq5 | Skript | Testen des GP-Modells an synthetischen Daten |
| 4 | GP.mqh | Klassenbibliothek | Die Klasse GaussianProcess, die Kernel, Likelihood-Funktion und Inferenzmethode kombiniert, sowie die Klasse GPOptimizationObjective, die für die Optimierung zuständig ist |
| 5 | Inference.mqh | Klassenbibliothek | Das Interface IInference und dessen Implementierungen, die Inferenzmethoden definieren |
| 6 | Kernels.mqh | Klassenbibliothek | Das Interface IKernel und Implementierungen von Kovarianz-Kerneln |
| 7 | Likelihoods.mqh | Klassenbibliothek | Interface ILikelihood und Implementierungen von Likelihood-Funktionen |
| 8 | StructUtils.mqh | Klassenbibliothek | Hilfsfunktionen und Datenstrukturen |
| 9 | DataRegression | CSV | GPsynthetic.mq5 Skriptdaten |
| 10 | DataClassification | CSV | GPsynthetic.mq5 Skriptdaten |
Übersetzt aus dem Russischen von MetaQuotes Ltd.
Originalartikel: https://www.mql5.com/ru/articles/19013
Beigefügte Dateien |
MQL5.zip
(50.8 KB)
Warnung: Alle Rechte sind von MetaQuotes Ltd. vorbehalten. Kopieren oder Vervielfältigen untersagt.
Dieser Artikel wurde von einem Nutzer der Website verfasst und gibt dessen persönliche Meinung wieder. MetaQuotes Ltd übernimmt keine Verantwortung für die Richtigkeit der dargestellten Informationen oder für Folgen, die sich aus der Anwendung der beschriebenen Lösungen, Strategien oder Empfehlungen ergeben.
Community of Scientists Optimization (CoSO): Praxis
Wir kehren zum Thema der Optimierung nach dem Modell einer wissenschaftlichen Gemeinschaft zurück. CoSO sollte nicht als fertige Lösung betrachtet werden, sondern als vielversprechende Forschungsplattform. Bei entsprechender Entwicklung kann CoSO seine Nische in Aufgaben finden, bei denen Anpassungsfähigkeit und Widerstandsfähigkeit gegenüber Veränderungen wichtig sind und die Rechenzeit nicht kritisch ist.
Quantenneuronales Netzwerk in MQL5 (Teil III): Ein virtueller Quantenprozessor auf Basis von Qubits
Der Artikel konzentriert sich auf die Erstellung eines Handelssystems mit einem echten Quantensimulator anstelle von mathematischen Analogien. Das System verwendet 3 virtuelle Qubits, Quantengatter und Superpositionsprinzipien zur Marktanalyse. Es ist als Handels-EA für MetaTrader 5 in MQL5 implementiert. Die wichtigste Errungenschaft ist der Übergang von der Simulation zu echten Quantenprinzipien der Finanzinformationsverarbeitung.
Eine alternative Log-datei mit der Verwendung der HTML und CSS
In diesem Artikel werden wir eine sehr einfache, aber leistungsfähige Bibliothek zur Erstellung der HTML-Dateien schreiben, dabei lernen wir auch, wie man eine ihre Darstellung einstellen kann (nach seinem Geschmack) und sehen wir, wie man es leicht in seinem Expert Advisor oder Skript hinzufügen oder verwenden kann.
Community of Scientists Optimization (CoSO): Theorie
Geheimnisse der effektiven Optimierung von Handelsstrategien in metaheuristischen Ansätzen. Community of Scientists Optimization ist ein neuer populationsbasierter Algorithmus, der von den Mechanismen der wissenschaftlichen Gemeinschaft inspiriert ist. Im Gegensatz zu traditionellen, von der Natur inspirierten Metaphern modelliert CoSO einzigartige Aspekte menschlicher wissenschaftlicher Tätigkeit: das Veröffentlichen von Ergebnissen in Fachzeitschriften, das Konkurrieren um Fördermittel und das Bilden von Forschungsteams.
Sie verpassen Handelsmöglichkeiten:
- Freie Handelsapplikationen
- Über 8.000 Signale zum Kopieren
- Wirtschaftsnachrichten für die Lage an den Finanzmärkte
Registrierung
Einloggen
Sie stimmen der Website-Richtlinie und den Nutzungsbedingungen zu.
Wenn Sie kein Benutzerkonto haben, registrieren Sie sich
Guten Tag,
Guten Tag!
Ist es das, was Sie möchten?
Nun werden nicht mehr die Prognose für den nächsten Schritt, sondern der Mittelwert und die Streuung der Trainingsstichprobe angezeigt.
Vielen Dank, das ist sehr nett von Ihnen.
Können die Modelle der Gaußschen Prozessregression im Online-Modus trainiert werden, sodass sie sich in Echtzeit an den Mittelwert und die Varianz der Daten anpassen können?
Vielen Dank, das ist sehr nett von Ihnen.
Ist es möglich, Regressionsmodelle auf Basis von Gauß-Prozessen im Online-Modus zu trainieren, damit sie sich in Echtzeit an den Mittelwert und die Varianz der Daten anpassen können?