MetaTrader 5 / 統計と分析
English Русский 中文 Español Deutsch Português 한국어 Français Italiano Türkçe
未知の確率密度関数のカーネル密度推定

未知の確率密度関数のカーネル密度推定

MetaTrader 5 — 統計と分析 | 27 10月 2015, 13:29
2 397 0
Victor
Victor

はじめに

MQL5 パフォーマンスの向上と PC 生産性の確かな上昇により MetaTrader 5 プラットフォームユーザーはかなり洗練された高度な数学的手法をマーケット分析に応用することができるようになりました。こういった手法は経済学、計量経済学、統計学に属しますが、いずれにせよそれらを利用するには確立密度関数を扱う必要があります。

多くの一般的手法は利用されるモデルにおいてデータ分布の正規性またはエラーの正規性を推定することを基にして発展してきました。また、使用されるモデルの多様なコンポーネントの分布を知り分析結果を推定することが必要な場合も多くあります。どちらの場合も未知の確率密度関数を推定することができる『ツール』（理想的には汎用的なもの）を作成するタスクがあります。

本稿では、未知の確率密度関数の推定に用いる汎用的アルゴリズムを多少認識しながらクラスを作成してみようと思います。着想としては推定を行うにあたり、外部手法は使用しない、というものでした。すなわちすべては MQL5 のみによって実現される、ということです。ただ、最終的に元の考えはある程度変更されました。確率密度関数の可視的推定タスクは2つの個別部分で構成されているのは明らかです。

すなわち、推定自体の計算とそれを可視化することです。たとえばグラフやダイアグラムとして表示することです。もちろん、計算は MQL5を用いて行い、一方可視化は HTML ページを作成し、それをウェブブラウザに表示することで行う必要があります。その解は最終的にベクトル形式でグラフィックを取得するために利用されました。

ただ、計算と結果表示部は個別に行われるので、読者のみなさんは別の方法を使って可視化することはもちろん可能です。また、グラフィック的なものを含み、様々なライブラリが MetaTrader 5（知る限り作業進行中です）に登場することを期待しています。MetaTrader 5 がグラフやダイアグラムを構築する高度な方法を提供する場合、提案の解決法で表示部分を変更するのは難しいことではありません。

シーケンス確率密度関数の推定のための真の汎用アルゴリズムを作成することは未達成の目標だとわかったことをあらかじめ知っておく必要があります。提案の解決法は高度に特殊化されたものではありませんが、完全に汎用的とも呼べません。問題は密度推定にあたっては最適化基準はひじょうに難しいということです。たとえば、正規分布や指数分布といった釣鐘型分布に対してです。

よって推定された分布に関して事前情報があれば、各指定ケースに対する適切な解決法を選択することはつねに可能です。ただ、それにもかかわらず、推定された密度の性質については何も知識がないということが前提です。そのような方法は確かに推定の質に影響しますが、まったく異なる密度を推定する可能性を証明することで 効果があると期待しましょう。

マーケットデータ分析非定常シーケンスを扱う必要がよくあることにより、短中シーケンス密度の推定にもっとも着目します。これは使用される推定方法の選択判断をする重大な局面です。

ヒストグラムと Pスプライン曲線 は100万以上の値を持つひじょうに長いシーケンスに対して効果的に使用できます。しかし、10～20の値を持つシーケンスに対してヒストグラムを効果的に構築しようとすると問題がいくつか発生します。よって、これより先、およそ10～10,000の値を持つシーケンスに主に着目していきます。


1. 確率密度関数推定手法

今日、きわめて多く、多かれ少なかれ人気の確率密度関数推定手法が知られています。たとえば、『確率密度推定』、『確率密度』、『密度推定』などのキーワードを使うとインターネットで簡単に見つけることができます。ただ、残念なことにその中でベストなものを選ぶことはできていません。すべてにメリット、デメリットがあります。

従来、密度推定にはヒストグラムが用いられています [1]。ヒストグラム（平滑化されたものを含め）を使用することで質の高い確率密度推定が可能となりますが、それはただ長いシーケンスを扱う場合のみです。先に述べたように、短いシーケンスを大きな数値のグループに分割するのは不可能で、一方2～3本のバーで構成されたヒストグラムではそのようなシーケンスの確率密度分布を説明することができません。よって、ヒストグラムは対象外としました。

その他ひじょうに有名な推定手法はカーネル密度推定です [2]。カーネル平滑化利用の考えは [3] によく示されています。よって、デメリットはあるものの、その手法を選択しました。この手法の実装に関連するいくつかの側面について以下で端的に述べます。

またひじょうに興味深い密度推定手法について述べる必要があります。それは『期待値最大化』アルゴリズムを用いたものです [4]。このアルゴリズムを使用することで、シーケンスをたとえば正規分布を持つコンポーネントに分割することが可能となります。個別のコンポーネントパラメータが決定されたのち、取得した分布曲線を合計することで密度推定値を得ることが可能です。この手法は [5] で述べられています。その他多くのもの同様、本稿執筆中はこの手法はまだ実装されておらず、検証もされていませんでした。様々なソースに述べられている密度推定手法の多くは実際にそのすべてを検証できないようになっています。

実装に選択したカーネル密度推定に話を進めていきます。


2. 確率密度関数のカーネル密度推定

確率密度関数のカーネル密度推定はカーネル平滑化手法を基にしています。この手法の原理はたとえば [6]、[7] に見ることができます。

カーネル平滑化の基本的考えは極めてシンプルです。MetaTrader 5 ユーザーは移動平均（MA）インディケータに詳しいものです。このインディケータは内部で平滑化される加重値によってシーケンスに従いスライドするウィンドウとして簡単に描くことができます。そのウィンドウは長方形、指数またはその他の形態をしています。カーネル平滑化において同様のスライドウィンドウを簡単に確認することができます（たとえば [3]）。ただこの場合は対称です。

カーネル平滑化において頻繁に使用されるウィンドウ例は [8] に見ることができます。カーネル平滑化にゼロ次回帰が使用されている場合、ウィンドウ（カーネル）に影響を与えたシーケンスの加重値はMAにおけるように単に平滑化されています。フィルタリングの問題に取り組む際、同種のウィンドウ関数アプリケーションを確認することができます。ただ、今は同じプロシージャがすこし異なる形で紹介されます。振幅周波数特性 および位相周波数特性が使用されます。またカーネル（ウィンドウ）はフィルタインパルス特性と呼ばれます。

こういった例はしばしば一つのことが異なる方法で表現可能なことを示しています。それはもちろん数学的道具に寄与します。が、その類の問題を議論する場合、混乱も招きます。

カーネル密度推定がすでに述べたカーネル平滑化と同様の原理を使っていても、そのアルゴリズムはやや異なります。

ある時点で密度推定を定義する式を見ます。

ここで

  • x - 長さ n を持つシーケンス
  • K - 対称的なカーネル
  • h - レンジ、平滑化パラメータ

ガウスカーネルのみのちに密度推定に使用されます。

以下のように上記の式からポイント X における密度はポイント X における値とシーケンス値の差により定義された量に対するカーネル値の合計として計算されます。それ以上に密度計算に使用される X ポイントはシーケンス自体の値とは一致しない可能性があります。

ここにカーネル密度推定アルゴリズムの基本実装手順があります。

  1. 平均値と入力シーケンス標準偏差の判定
  2. 入力シーケンスの正規化各値から以前に取得した平均値の減算および標準偏差値での分割そのような正規化後、元のシーケンスはゼロ平均を持ち標準偏差は1となります。この正常化は密度計算には必要ありませんが、結果のグラフを統一することができます。 X 尺度のあらゆるシーケンスについては標準偏差ユニット内で表現される値があります。
  3. 正常化されたシーケンス内における最大値、最小値の確認
  4. 配列を2つ作成します。望む数値に対応するそのサイズは結果グラフに表示されます。たとえば、チャートが200ポイントを使用して作成されるならば、配列サイズは各々およそ200の値を持ちます。
  5. 結果を格納するため作成した配列の一つを保有二番目の配列はポイント値を形成するために使用されます。そのポイントに対して密度推定が行われます。このために、前に準備した最大値と最小値の間に等間隔値を200作成し、準備済みの配列にそれらを保存する必要があります。
  6. 前に示されている式を使用し、ステップ4で準備した配列に結果を保存する200のテストポイント（今回の場合）で密度推定を実行します。

以下はこのアルゴリズムのソフトウェアへの実装です。

//+------------------------------------------------------------------+
//|                                                        CDens.mqh |
//|                                                    2012, victorg |
//|                                              https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "2012, victorg"
#property link      "https://www.mql5.com"

#include <Object.mqh>
//+------------------------------------------------------------------+
//| Class Kernel Density Estimation                                  |
//+------------------------------------------------------------------+
class CDens:public CObject
  {
public:
   double            X[];              // Data
   int               N;                // Input data length (N >= 8)
   double            T[];              // Test points for pdf estimating
   double            Y[];              // Estimated density (pdf)
   int               Np;               // Number of test points (Npoint>=10, default 200)
   double            Mean;             // Mean (average)
   double            Var;              // Variance
   double            StDev;            // Standard deviation
   double            H;                // Bandwidth
public:
   void              CDens(void);
   int               Density(double &x[],double hh);
   void              NTpoints(int n);
private:
   void              kdens(double h);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
void CDens::CDens(void)
  {
   NTpoints(200);            // Default number of test points
  }
//+------------------------------------------------------------------+
//| Setting number of test points                                    |
//+------------------------------------------------------------------+
void CDens::NTpoints(int n)
  {
   if(n<10)n=10;
   Np=n;                    // Number of test points
   ArrayResize(T,Np);        // Array for test points
   ArrayResize(Y,Np);        // Array for result (pdf)
  }
//+------------------------------------------------------------------+
//| Density                                                          |
//+------------------------------------------------------------------+
int CDens::Density(double &x[],double hh)
  {
   int i;
   double a,b,min,max,h;

   N=ArraySize(x);                           // Input data length
   if(N<8)                                  // If N is too small
     {
      Print(__FUNCTION__+": Error! Not enough data length!");
      return(-1);
     }
   ArrayResize(X,N);                         // Array for input data
   ArrayCopy(X,x);                           // Copy input data
   ArraySort(X);
   Mean=0;
   for(i=0;i<N;i++)Mean=Mean+(X[i]-Mean)/(i+1.0); // Mean (average)
   Var=0;
   for(i=0;i<N;i++)
     {
      a=X[i]-Mean;
      X[i]=a;
      Var+=a*a;
     }
   Var/=N;                                  // Variance
   if(Var<1.e-250)                           // Variance is too small
     {
      Print(__FUNCTION__+": Error! The variance is too small or zero!");
      return(-1);
     }
   StDev=MathSqrt(Var);                      // Standard deviation
   for(i=0;i<N;i++)X[i]=X[i]/StDev;          // Data normalization (mean=0,stdev=1)
   min=X[ArrayMinimum(X)];
   max=X[ArrayMaximum(X)];
   b=(max-min)/(Np-1.0);
   for(i=0;i<Np;i++)T[i]=min+b*(double)i;    // Create test points
//-------------------------------- Bandwidth selection
   h=hh;
   if(h<0.001)h=0.001;
   H=h;
//-------------------------------- Density estimation
   kdens(h);

   return(0);
  }
//+------------------------------------------------------------------+
//| Gaussian kernel density estimation                               |
//+------------------------------------------------------------------+
void CDens::kdens(double h)
  {
   int i,j;
   double a,b,c;

   c=MathSqrt(M_PI+M_PI)*N*h;
   for(i=0;i<Np;i++)
     {
      a=0;
      for(j=0;j<N;j++)
        {
         b=(T[i]-X[j])/h;
         a+=MathExp(-b*b*0.5);
        }
      Y[i]=a/c;                 // pdf
     }
  }
//--------------------------------------------------------------------

NTpoints() メソッドにより等間隔のテストポイントの必要な数値を設定することができます。それに対して密度推定が実行されます。このメソッドは Density() メソッドを呼ぶ前に呼ぶ必要があります。インプットデータおよび配列値（平滑化パラメータ）を持つ配列へのリンクはメソッドが呼ばれる際にそのパラメータとして Density() メソッドにわたされます。

Density() メソッドは無事に完了した場合ゼロを返します。一方テストポイント値と推定結果はそのクラスの配列 T[] および Y[] にセットされます。

配列サイズは NTpoints()にアクセスするとき設定され、デフォルトでは200の値となっています。

以下のサンプルスクリプトは CDens クラスの使用を示しています。

#include "CDens.mqh"
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int i;
   int ndata=1000;          // Input data length
   int npoint=300;          // Number of test points
   double X[];              // Array for data 
   double T[];              // Test points for pdf estimating
   double Y[];              // Array for result

   ArrayResize(X,ndata);
   ArrayResize(T,npoint);
   ArrayResize(Y,npoint);
   for(i=0;i<ndata;i++)X[i]=MathRand();// Create input data
   CDens *kd=new CDens;
   kd.NTpoints(npoint);    // Setting number of test points
   kd.Density(X,0.22);     // Density estimation with h=0.22
   ArrayCopy(T,kd.T);       // Copy test points
   ArrayCopy(Y,kd.Y);       // Copy result (pdf)
   delete(kd);

// Result: T[]-test points, Y[]-density estimation
  }
//--------------------------------------------------------------------

テストポイントの入力シーケンスと推定結果を格納する配列はこの例に準備されています。そして X[] 配列には例ではランダム関数が書かれます。CDens クラスのコピーはすべてのデータが準備されたあと作成されます。

その後、NTpoints()メソッドの起動によって必要なテストポイント数が設定（ここの場合 npoint=300）されます。ポイントのデフォルト数を変更する必要がない場合は NTpoints() メソッドの起動は除外される可能性があります。

計算された値はDensity() メソッドの起動後、定義済み配列にコピーされます。そして以前に作成された CDens クラスのコピーは削除されます。この例は CDens クラスの連携のみ示しています。取得する結果（配列 T[] および Y[] ）へのそれ以上の処理は行われません。

これら結果がグラフ作成に使用されるということならばT[] 配列からの値を X グラフ目盛上に入れることによって簡単に行えます。一方、Y[] 配列からの適切な値は Y 目盛上に入れる必要があります。


3. 最適レンジ選択

図1は正規分布法則と様々な h レンジ値を持つシーケンスに対する密度推定グラフです。

推定は前述の CDens クラスを用いて行われます。グラフはHTML ページ形式で構築されています。そのようなグラフを作成するメソッドは本稿の末尾にて紹介します。HTML 形式でのグラフおよびダイアグラム作成は [9]に見ることができます。

 図1　さまざまな h レンジ値に対する密度推定

図1　さまざまな h レンジ値に対する密度推定

図1はまた3つの密度推定についての真の正規分布密度曲線（ガウス分布）を表示しています。この場合、もっとも適切な推定結果は h=0.22 であったことが簡単に確認できます。その他2例では明確な『過平滑化』および『過小平滑化』が観測されます。

図1は明らかにカーネル密度推定を利用する際、正しい h レンジ選択が重要であることを示しています。In case h 値が誤って選択された場合、推定値は真の密度に対して大きく外れるか大幅に分散します。

最適な h レンジ選択に多くの多様な作業が費やされることになります。h 選択にはシンプルでよく証明されているシルバーマンの経験則がよく利用されます。（[10]参照）。

今回は A がシーケンスの標準偏差値の最低値で 1.34で除算され四分位範囲です。る

入力シーケンスが上記 CDens クラスで標準偏差1であったことから、次のコード部分を使用してこの報告を取り入れるのは極めて簡単です。

  ArraySort(X);
  i=(int)((N-1.0)/4.0+0.5);  
  a=(X[N-1-i]-X[i])/1.34;      // IQR/1.34
  a=MathMin(a,1.0);
  h=0.9*a/MathPow(N,0.2);       // Silverman's rule of thumb

推定値はその形式によって正規に近い確率密度関数を持つシーケンスにあてはまります。

多くの場合、四分位範囲を考慮することにより、判定された密度偏差が正規密度偏差から外れる場合、h 値を低い側へと調整することができます。それによりこの判定手法は十分多用途なものになります。よってこの経験則は初期の基本 h 値判定として利用価値があります。またそれには長々として計算は必要ありません。

漸近的また実験的 h レンジ値推定以外にも入力シーケンス自体の分析を基にした多様な手法があります。この場合は、最適な h 値は未知の密度の初期推定を考慮することで決定します。これらいくつかの手法の有効性比較は [11]、[12]にあります。

さまざまなソースに発表されている検証結果によると、 Sheather-Jones プラグイン法（SJPI）がもっとも効果的レンジ推定手法のひとつです。ではこの手法を選びます。長ったらしい数式で本稿をオーバーロードしないために、この先、手法の特徴の一部だけお話していきます。必要であればこの手法の数学的基礎は [13]を参照することができます。

ガウスカーネルを使用します。これはを持たないカーネルです（たとえば、パラメータが負の無限大から正の無限大に代わるとき特定される、など）。それは上記の式から導かれるため、直接計算法と無制限サポートを使用する N 長のシーケンスのポイント M における密度を推定するとき、 O(M*N) オペレーションが必要となります。各シーケンスポイントで推定を行う必要があれば、この値は O(N*N) オペレーションまで造塊し、計算に費やされる時間はシーケンス長の二乗に比例して増えます。

計算リソースに対するそのような高い要求はカーネル平滑化手法の大きなデメリットの一つです。SJPI メソッドに渡せば、計算に必要な量に関して劣ることが判ります。手短に話すと、まず SJPI メソッドを実装する際、入力シーケンス全体長について2つの密度導関数合計を二度計算する必要があります。

その後、取得した結果を使用して推定計算を繰り返し行います。推定値はゼロとなります。引数を見つけるにはニュートン・ラプソン法 [14] が利用される可能性があります。その引数ではこの関数はゼロです。直接計算する場合は、最適レンジ値を判断するため密度推定計算の10倍多く SJPI 手法が適用されます。

カーネル平滑化および密度カーネル密度推定においては多様な計算のスピードアップ方法があります。われわれの場合ではガウスカーネル法が利用され、その値は4より大きい引数値に対して無視してよいほどに小さいと推測される可能性があります。よって、引数の値が4より大きければ、カーネル値を計算する必要はありません。このため必要な計算数を部分的に減らすことが可能です。そのような推定値を基にしたグラフとすべて計算した推定値を基にしたグラフはほとんど区別できません。

計算をスピードアップする別のシンプルな方法は推定が行われるポイント数を減らすことです。前述のよう M が Nより小さければ、必要とされる処理数は O(N*N) から O(M*N)に減ります。たとえば N=100 000のような長いシーケンスがあったなら M=200 ポイントでのみ推定を行うことができます。. このように計算時間を大幅に削減することが可能なのです。

また、必要な計算数を減らすためにより複合的な手法も利用可能です。たとえば、高速フーリエ変換アルゴリズムや ウェーブレット変換を使用した推定です。入力シーケンスディメンションを減らすことを基にした手法（たとえば『データビニング』 [15]）はひじょうに長いシーケンスに対して問題なく適用することができます。

高速ガウス変換 [16] もまた前述の手法に加え、が売るケーネルが使用される場合、計算のスピードアップに適用可能です。ガウス変換 [17] を基にしたアルゴリズムを使用し、レンジ値推定のSJPI 手法を取り入れます。上記のリンクは両方の手法記述ならびに提案されているアルゴリズムを実装するプログラムコードを含む資料につながっています。


4. Sheather-Jonesプラグイン法(SJPI)

SJPI 法の基礎を決定する数式の場合、本稿の [17] に見られるアルゴリズムの実装における数学的基礎をコピーすることはしません。必要であれば [17]での出版物を検証することができます。

CSJPlugin クラス [17]にある資料を基に作成されました。このクラスは最適な h レンジ値を計算し、ただ一つのパブリックメソッドをインクルードしています。それはダブル CSJPlugin::SelectH（double &px[]、double h、double stdev=1）です。

次の引数が起動時このメソッドに渡されます。

  • double &px[] – 元のシーケンスを含む配列へのリンク
  • ダブル h は最初の h レンジ値で、これに関連しニュートン・ラプソン アルゴリズムを用いて式を解くとき検索が行われます。この値が検索している値に可能な限り近いことが望まれます。これは式を解く時間を大幅に削減し、ニュートン・ラプソン法がその信頼性を失う場合の確率を減らします。シルバーマンの経験則を用いて見つけられる値は初期 h 値として使用されます。
  • double stdev=1 – 元のシーケンス標準偏差値デフォルト値は1です。この場合、このメソッドは前に述べられている CDens クラスと併用するのでデフォルト値を変更する必要はありません。そのクラスでは元のシーケンスがすでに正規化されておりその標準偏差は1です。

正常終了の場合、SelectH() メソッドは取得された最適 h レンジ値を返します。失敗の場合のパラメータとして、初期 h 値が返されます。CSJPlugin クラスのソースコードは CSJPlugin.mqh ファイルにあります。

この実装の失敗に関して明確にする必要があります。

ソースシーケンスは一度に [0,1] 区間に変換され、一方で h レンジの初期値は元のシーケンス変換スケールに比例して正規化されます。再正規化は計算中に取得される最適 h 値に適用されます。

eps=1e-4 – 密度推定とその導関数の計算精度、およびP_UL=500 – アルゴリズムの最大内部反復許容数が CSJPlugin クラスコンストラクタ内に設定されます。詳細は [17]を参照ください。

IMAX=20 – ニュートン・ラプソン法に対する最大反復許容数、および PREC=1e-4 – この手法を用いて式を解く精度が SelectH() メソッドに設定されます。

ニュートン・ラプソン法の従来的使用には各値反復において同一ポイントで対象となる関数とその導関数の計算が必要です。この場合、導関数の計算はその引数に小さいインクリメントを加えることで計算される推定値に置き換わっています。

図2には異なる2つの最適 h レンジ値推定手法を用いた例が示されています。

図2　最適 h レンジ値推定比較

図2　最適 h レンジ値推定比較

図2にはランダムシーケンスに対する2つの推定が示されています。その真の密度は赤いグラフ（Pattern）として表示されています。

推定は長さに 10 000 エレメントあるシーケンスに対して行われました。シルバーマンの経験則を用いる場合このシーケンスに対して0.14 の h レンジ値が取得されました。一方 SJPI 手法を用いた場合は0.07でした。

図2に表示されているこれら2つの h 値に対するカーネル密度推定の結果を検証すると、SJPI 手法適用はシルバーマンの経験則に比べるとずっと興味深い h 推定値を得るのに役立っていることが判ります。見てのとおり、鋭い頂点はよりよく描かれ、また h=0.07の場合、傾斜の底部分でほとんどばらつきがありません。

これは予想されていたため、SJPI 手法は多くの場合で思い通りの h レンジ推定値となります。ひじょうに速いアルゴリズム [17] CSJPlugin クラス作成に使用されているにもかかわらず、この手法を用いた h 値推定は長いシーケンスに対してはまだ時間がかかりすぎます。

この手法のその他の欠点は 10～30 の値で構成される短いシーケンスに対しては h 値を過剰に推定する傾向にあることです。SJPI 手法を用いた推定値はシルバーマンの経験則による h 推定値を上まわる可能性があります。

以下のルールはこういった欠点をいくらか補うためにさらなる実装において使用されます。

  • シルバーマンの推定は h レンジに対するデフォルトの推定値であり、一方別のコマンドによって SJPI 手法を起動することは可能であります
  • SJPI メソッドを使用する場合、前述の2つの手法から取得されたうちの最低値はつねに h の最終値として使用されます。


5. 境界効果

密度推定において最適レンジ値を可能な限り使用したいという思いが前述のかなり面倒な CSJPlugin クラスの作成につながりました。ただレンジサイズとカーネル平滑化手法の高い資源強度を指定すること以外にも問題がひとつあります。それがいわゆる境界効果です。

問題はシンプルです。それは区間 [-1,1]において決められるカーネルを使用して表示されます。そのようなカーネルは有界サポートを持つカーネルと呼ばれます。指定範囲外ではそれはゼロです（存在しない）。

図3　レンジ境界におけるカーネル減少

図3　レンジ境界におけるカーネル減少

図3に示されるように、最初のケースではその中心に対してコアは区間 [-1,1] に位置する元データを完全に重なっています。カーネルが ずれて（右に）いるのでデータが選択されたカーネル関数を完全に使用するに十分ではない状況が明らかになります。

カーネルはすでに最初のケースよりも少ないデータと重なっています。最悪のケースはカーネルの中心がデータ系列の境界にくるときです。その場合カーネルが重なるデータ量は 50% まで減ります。密度推定に使われるデータ数のそのような減少は推定値の大幅な移動につながり、またレンジ境界付近のポイントにおけるばらつきを増やします。

図3はレンジ境界において有界サポートを持つカーネル（Epanechnikovカーネル）に対する減少例を示しています。無限範囲に定義されるガウスカーネル（非有界サポート）はカーネル密度推定法の実装中に使用されることに注意が必要です。理論的にはそのようなカーネルの切り捨てがつねに行われるべきです。ただ、このカーネル値が大きな引数に対してほとんとゼロであることを考えると、境界効果は有界サポートを持つカーネルに対するものと同じように出現します。

この特徴は図1、図2に表されるケースにおける密度推定結果に影響はない可能性があります。というのは、どちらの場合も推定は分布、境界においてほとんどゼロに下がる確率密度関数に対して行われているからです。

レンジ境界でカーネルを切り捨てる影響を示す正の整 X=1,2,3,4,5,6,…n で構成されるシーケンスを作成します。そのようなシーケンスは確率密度分布の均等法則があります。それはこのシーケンスの密度推定は水平な直線でノンゼロレベルに位置していることを意味します。

図4　分布の均等法則を持つシーケンスに対する密度推定

図4　分布の均等法則を持つシーケンスに対する密度推定

予想どおり、図4は明らかにレンジ境界において密度推定の大幅な移動があることを示しています。程度の差はあってもそのような移動を減らす手法は複数あります。おおまかに次のようなグループに分けることができます。

  • データ反映法
  • データ変換法
  • 疑似データ法
  • 境界カーネル法

反映されたデータを使用する考えはデータを追加することで入力シーケンスが故意に増やされているというもので、それはシーケンスレンジ境界に関するシーケンスの一種の正反射です。そのような増加後、同じポイントに対して密度推定が行われます。ただ元のシーケンスについては故意に追加されたデータは推定にも使用されます。

データ変換を持つ手法はそのレンジ境界に近いシーケンス変換に重点を置きます。たとえば密度推定中レンジ境界に近づくといくばくかデータスケールを変形することのできる対数またはその他の変換を利用することは可能です。

いわゆる疑似データは元シーケンスを故意に拡張する（広げる）ために使用されます。それは元のシーケンス値を基に計算されるものです。そのデータは境界においてその変動を捉え、もっとも適した方法でそれを補完する役割をします。

境界カーネル法に特化した公表文献は数多くあります。これら手法においては のカーネルは境界に近づくにつれその形を変えます。カーネル形態は境界において推定値の移動を補正するために変化します。

レンジ境界に出現するゆがみを補正することに特化した手法、それらの比較、効果判定は [18] に見ることができます。

データ反映法はいくつかの短い実験後使用するために選択されました。その選択はこの手法が密度推定が負の値をもつ場合を示さないためです。また、この手法は複合的数学の計算を必要としません。オペレーションのトータル数はシーケンスの各推定をする必要があるためまだ増加し、意図的にその長さを増やします。

その手法を実装する方法は2とおりあります。まず必要なデータによって元のシーケンスを 補足 し、その過程でサイズを3倍に増やします。そうすると前にCDens クラスが紹介されたときに示されたと同じ方法で推定を行うことができます。その次にインプットデータ配列を拡張しないことが可能です。特定の方法で再びそこからデータを選択します。二番目の方法は実装のために選択されました。

前述の CDens クラスにおいて密度推定は void kdens（ダブル h） 関数に実装されています。境界ばらつき修正を加えてこの関数を変更します。

以下が拡大した関数表記です。

//+------------------------------------------------------------------+
//| Kernel density estimation with reflection of data                |
//+------------------------------------------------------------------+
 void kdens(double h)
  {
  int i,j;
  double a,b,c,d,e,g,s,hh;
  
  hh=h/MathSqrt(0.5);
  s=sqrt(M_PI+M_PI)*N*h;
  c=(X[0]+X[0])/hh;
  d=(X[N-1]+X[N-1])/hh;
  for(i=0;i<Np;i++)
    {
    e=T[i]/hh; a=0;
    g=e-X[0]/hh;   if(g>-3&&g<3)a+=MathExp(-g*g);
    g=e-X[N-1]/hh; if(g>-3&&g<3)a+=MathExp(-g*g);
    for(j=1;j<N-1;j++)
      {
      b=X[j]/hh;
      g=e-b;   if(g>-3&&g<3)a+=MathExp(-g*g);
      g=d-e-b; if(g>-3&&g<3)a+=MathExp(-g*g);
      g=c-e-b; if(g>-3&&g<3)a+=MathExp(-g*g);
      }
    Y[i]=a/s;                                        // pdf
    }
  }

この関数は X[] 配列のソースデータがすでに関数が呼ばれるときソートされていると仮定し実装されています。これにより極端なレンジ値に対応する2つのシーケンスエレメントを他の処理から簡単に除外することができます。この関数を実装する際、可能なときはいつでも数学的処理をループの外に出すようにされました。結果、あまり目立ちませんがこの関数は使用されるアルゴリズムの考えを反映します。

引数の大きい値に対してカーネル値を計算しない場合、計算数を減らすことは可能であることはすでに述べられています。それがまさに上記の関数に実装されたものです。この場合、評価済み密度グラフが作成されたのちのいかなる変更もわかります。

kdens() 関数の変更バージョンを使用するときは、図4に示される密度推定は直線に変換され、そうすると境界のくぼみは完全に消えます。ただ、そのような理想的な修正は境界付近の分布がゼロ勾配の場合にのみ可能です。たとえば水平な直線で表される場合です。

境界付近で推定された分布密度が急上昇したり急降下すると、選択されたデータ反映法は境界効果を完全に調整することができなくなります。これは以下の図で示されます。

図5　段階的に変化する確率密度

図5　段階的に変化する確率密度関数

図6　直線的に増加する確率密度

図6　直線的に増加する確率密度関数

図5と図6は kdens() 関数の元バージョンを用いて取得される密度推定値（赤の線）およびデータ反映法を実装することによる変化を考慮して取得された推定値（ブルーの線）を表示しています。図5において境界効果は完全に修正されています。一方図6では、境界付近の移動は完全には除去されていません。境界付近で推定された分布密度が急上昇したり急降下すると、その境界付近では反映はややなめらかに現れます。

境界効果の調整のために選択されるデータ反映法は既知の手法の中でベストでもなければ最悪でもありません。この手法によって境界効果がすべての場合において除去されるわけではなくてもそれは十分安定しており実装は簡単です。この手法によりロジカルで予測可能な結果を得ることが可能となります。


6. 最終実装CKDensity クラス

h レンジ値を自動で選択し、前に作成したCDensクラスに対して境界効果を修正する機能を追加します。

下記がそのように変更したクラスのソースコードです。

//+------------------------------------------------------------------+
//|                                                    CKDensity.mqh |
//|                                                    2012, victorg |
//|                                              https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "2012, victorg"
#property link      "https://www.mql5.com"

#include <Object.mqh>
#include "CSJPlugin.mqh"
//+------------------------------------------------------------------+
//| Class Kernel Density Estimation                                  |
//+------------------------------------------------------------------+
class CKDensity:public CObject
  {
public:
   double            X[];            // Data
   int               N;              // Input data length (N >= 8)
   double            T[];            // Test points for pdf estimating
   double            Y[];            // Estimated density (pdf)
   int               Np;             // Number of test points (Npoint>=10, default 200)
   double            Mean;           // Mean (average)
   double            Var;            // Variance
   double            StDev;          // Standard deviation
   double            H;              // Bandwidth
   int               Pflag;          // SJ plug-in bandwidth selection flag
public:
   void              CKDensity(void);
   int               Density(double &x[],double hh=-1);
   void              NTpoints(int n);
   void    PluginMode(int m) {if(m==1)Pflag=1; else Pflag=0;}
private:
   void              kdens(double h);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
void CKDensity::CKDensity(void)
  {
   NTpoints(200);                            // Default number of test points
   Pflag=0;                                  // Not SJ plug-in
  }
//+------------------------------------------------------------------+
//| Setting number of test points                                    |
//+------------------------------------------------------------------+
void CKDensity::NTpoints(int n)
  {
   if(n<10)n=10;
   Np=n;                                    // Number of test points
   ArrayResize(T,Np);                        // Array for test points
   ArrayResize(Y,Np);                        // Array for result (pdf)
  }
//+------------------------------------------------------------------+
//| Bandwidth selection and kernel density estimation                |
//+------------------------------------------------------------------+
int CKDensity::Density(double &x[],double hh=-1)
  {
   int i;
   double a,b,h;

   N=ArraySize(x);                           // Input data length
   if(N<8)                                   // If N is too small
     {
      Print(__FUNCTION__+": Error! Not enough data length!");
      return(-1);
     }
   ArrayResize(X,N);                         // Array for input data
   ArrayCopy(X,x);                           // Copy input data
   ArraySort(X);                             // Sort input data
   Mean=0;
   for(i=0;i<N;i++)Mean=Mean+(X[i]-Mean)/(i+1.0); // Mean (average)
   Var=0;
   for(i=0;i<N;i++)
     {
      a=X[i]-Mean;
      X[i]=a;
      Var+=a*a;
     }
   Var/=N;                                  // Variance
   if(Var<1.e-250)                           // Variance is too small
     {
      Print(__FUNCTION__+": Error! The variance is too small or zero!");
      return(-1);
     }
   StDev=MathSqrt(Var);                      // Standard deviation
   for(i=0;i<N;i++)X[i]=X[i]/StDev;          // Data normalization (mean=0,stdev=1)
   a=X[0];
   b=(X[N-1]-a)/(Np-1.0);
   for(i=0;i<Np;i++)T[i]=a+b*(double)i;      // Create test points

//-------------------------------- Bandwidth selection
   if(hh<0)
     {
      i=(int)((N-1.0)/4.0+0.5);
      a=(X[N-1-i]-X[i])/1.34;               // IQR/1.34
      a=MathMin(a,1.0);
      h=0.9*a/MathPow(N,0.2);                // Silverman's rule of thumb
      if(Pflag==1)
        {
         CSJPlugin *plug=new CSJPlugin();
         a=plug.SelectH(X,h);              // SJ Plug-in
         delete plug;
         h=MathMin(a,h);
        }
     }
   else {h=hh; if(h<0.005)h=0.005;}          // Manual select
   H=h;

//-------------------------------- Density estimation
   kdens(h);

   return(0);
  }
//+------------------------------------------------------------------+
//| Kernel density estimation with reflection of data                |
//+------------------------------------------------------------------+
void CKDensity::kdens(double h)
  {
   int i,j;
   double a,b,c,d,e,g,s,hh;

   hh=h/MathSqrt(0.5);
   s=sqrt(M_PI+M_PI)*N*h;
   c=(X[0]+X[0])/hh;
   d=(X[N-1]+X[N-1])/hh;
   for(i=0;i<Np;i++)
     {
      e=T[i]/hh; a=0;
      g=e-X[0]/hh;   if(g>-3&&g<3)a+=MathExp(-g*g);
      g=e-X[N-1]/hh; if(g>-3&&g<3)a+=MathExp(-g*g);
      for(j=1;j<N-1;j++)
        {
         b=X[j]/hh;
         g=e-b;   if(g>-3&&g<3)a+=MathExp(-g*g);
         g=d-e-b; if(g>-3&&g<3)a+=MathExp(-g*g);
         g=c-e-b; if(g>-3&&g<3)a+=MathExp(-g*g);
        }
      Y[i]=a/s;                               // pdf
     }
  }

密度（double &x[],double hh=-1）法はそのクラスに対しては基本的なものです。最初の引数は x[] 配列へのリンクで、そこでは分析された入力シーケンスが含まれる必要があります。入力シーケンスの長さは8エレメント未満です。二番目（オプション）の引数は h レンジ値を強制的に設定するのに使用されます。

任意の正の数が指定されます。設定された h 値はつねに下方より 0.005 分制限されます。このパラメータが負の値である場合、レンジ値は自動的に選択されます。Density() メソッドは正常終了の場合ゼロを返し、非常終了の場合は負の値を返します。

Density() メソッドの起動と正常終了後、T[] 配列がテストポイント値を持ちます。それに対して密度推定が行われ、また Y[] 配列はこれら推定値を持ちます。X[] 配列は正規化され、ソートされた入力シーケンスを持ちます。このシーケンスの平均値はゼロであり、標準偏差値は1です。

以下の値にはクラスメンバーである変数に含まれます。

  • N – 入力シーケンス長
  • Np – テストポイント数。ここで密度推定が行われます。
  • Mean – 入力シーケンス平均値
  • Var – 入力シーケンスバリアンス
  • StDev – 入力シーケンス標準偏差
  • H – レンジ値（平滑化パラメータ）
  • Pflag – h レンジの最適値を決定するための SJPI メソッド応用フラグ

NTpoints(int n) メソッドはテストポイントを設定するのに使用されます。そこで密度推定が実行されます。テストポイント数の設定は Density() 基本メソッドが呼ばれる前に行う必要があります。NTpoints(int n) メソッドを使用しない場合は、デフォルトポイント数200が設定されます。

PluginMode(int m) メソッドにより最適なレンジ値を判定する SJPI メソッドの使用が可能になり、また可能でなくなります。このメソッド起動時に1に等しいパラメータが使用されると、SJPI メソッドが h レンジを指定するのに適用されます。このパラメータ値が1でない場合、SJPI メソッドは使用されません。PluginMode(int m) メソッドが起動されると、デフォルトで SJPI メソッドは無効となります。


7. グラフ作成

本稿執筆にあたり [9] にあるメソッドを利用いたしました。その文献では highcharts [19] JavaScript ライブラリおよびその 呼び出しすべてのフルスペックを含む既製のHTMLページのアプリケーションが取り上げられています。のちに MQL プログラムが表示されるデータを持つテキストファイルを作成し、そのファイルは上述の HTML に連携します。

その場合、表示されるグラフストラクチャにおいてあらゆる修正をおこなうために実装されたJavaScript コードを変更することでHTMLページを編集することが必要です。これを避けるため、シンプルなインターフェースが開発されました。それにより highcharts ライブラリのJavaScriptメソッドを MQL プログラムから直接呼びだすことができます。

インターフェースを作成する際、highcharts ライブラリのすべての機能に対するアクセスを実装するタスクはありません。本稿執筆中、以前に作成された HTML ファイルの変更を許可しない機能のみ実装されていました。

以下に CLinDraw クラスのソースコードを示します。このクラスは連携に highcharts ライブラリメソッドを提供しています。しかしこのコードは完全なソフトウェア実装とはみなさないことです。グラフィック JavaScript ライブラリを持つインターフェースを作成する方法の一例を示しているにすぎません。

//+------------------------------------------------------------------+
//|                                                     CLinDraw.mqh |
//|                                                    2012, victorg |
//|                                              https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "2012, victorg"
#property link      "https://www.mql5.com"

#include <Object.mqh>
#import "shell32.dll"
int ShellExecuteW(int hwnd,string lpOperation,string lpFile,string lpParameters,
                  string lpDirectory,int nShowCmd);
#import
#import "kernel32.dll"
int DeleteFileW(string lpFileName);
int MoveFileW(string lpExistingFileName,string lpNewFileName);
#import
//+------------------------------------------------------------------+
//| type = "line","spline","scatter"                                 |
//| col  = "r,g,b,y"                                                 |
//| Leg  = "true","false"                                            |
//| Reference: http://www.highcharts.com/                            |
//+------------------------------------------------------------------+
class CLinDraw:public CObject
  {
protected:
   int               Fhandle;           // File handle
   int               Num;               // Internal number of chart line
   string            Tit;               // Title chart
   string            SubTit;            // Subtitle chart
   string            Leg;               // Legend enable/disable
   string            Ytit;              // Title Y scale
   string            Xtit;              // Title X scale
   string            Fnam;              // File name

public:
   void              CLinDraw(void);
   void    Title(string s)      { Tit=s; }
   void    SubTitle(string s)   { SubTit=s; }
   void    Legend(string s)     { Leg=s; }
   void    YTitle(string s)     { Ytit=s; }
   void    XTitle(string s)     { Xtit=s; }
   int               AddGraph(double &y[],string type,string name,int w=0,string col="");
   int               AddGraph(double &x[],double &y[],string type,string name,int w=0,string col="");
   int               AddGraph(double &x[],double y,string type,string name,int w=0,string col="");
   int               LDraw(int ashow=1);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
void CLinDraw::CLinDraw(void)
  {
   Num=0;
   Tit="";
   SubTit="";
   Leg="true";
   Ytit="";
   Xtit="";
   Fnam="CLinDraw.txt";
   Fhandle=FileOpen(Fnam,FILE_WRITE|FILE_TXT|FILE_ANSI);
   if(Fhandle<0)
     {
      Print(__FUNCTION__,": Error! FileOpen() error.");
      return;
     }
   FileSeek(Fhandle,0,SEEK_SET);               // if file exists
  }
//+------------------------------------------------------------------+
//| AddGraph                                                         |
//+------------------------------------------------------------------+
int CLinDraw::AddGraph(double &y[],string type,string name,int w=0,string col="")
  {
   int i,k,n;
   string str;

   if(Fhandle<0)return(-1);
   if(Num==0)
     {
      str="$(document).ready(function(){\n"
          "var lp=new Highcharts.Chart({\n"
          "chart:{renderTo:'lplot'},\n"
          "title:{text:'"+Tit+"'},\n"
          "subtitle:{text:'"+SubTit+"'},\n"
          "legend:{enabled:"+Leg+"},\n"
          "yAxis:{title:{text:'"+Ytit+"'}},\n"
          "xAxis:{title:{text:'"+Xtit+"'},showLastLabel:true},\n"
          "series:[\n";
      FileWriteString(Fhandle,str);
     }
   n=ArraySize(y);
   str="{type:'"+type+"',name:'"+name+"',";
   if(col!="")str+="color:'rgba("+col+")',";
   if(w!=0)str+="lineWidth:"+(string)w+",";
   str+="data:[";
   k=0;
   for(i=0;i<n-1;i++)
     {
      str+=StringFormat("%.5g,",y[i]);
      if(20<k++){k=0; str+="\n";}
     }
   str+=StringFormat("%.5g]},\n",y[n-1]);
   FileWriteString(Fhandle,str);
   Num++;
   return(0);
  }
//+------------------------------------------------------------------+
//| AddGraph                                                         |
//+------------------------------------------------------------------+
int CLinDraw::AddGraph(double &x[],double &y[],string type,string name,int w=0,string col="")
  {
   int i,k,n;
   string str;

   if(Fhandle<0)return(-1);
   if(Num==0)
     {
      str="$(document).ready(function(){\n"
          "var lp=new Highcharts.Chart({\n"
          "chart:{renderTo:'lplot'},\n"
          "title:{text:'"+Tit+"'},\n"
          "subtitle:{text:'"+SubTit+"'},\n"
          "legend:{enabled:"+Leg+"},\n"
          "yAxis:{title:{text:'"+Ytit+"'}},\n"
          "xAxis:{title:{text:'"+Xtit+"'},showLastLabel:true},\n"
          "series:[\n";
      FileWriteString(Fhandle,str);
     }
   n=ArraySize(x);
   str="{type:'"+type+"',name:'"+name+"',";
   if(col!="")str+="color:'rgba("+col+")',";
   if(w!=0)str+="lineWidth:"+(string)w+",";
   str+="data:[";
   k=0;
   for(i=0;i<n-1;i++)
     {
      str+=StringFormat("[%.5g,%.5g],",x[i],y[i]);
      if(20<k++){k=0; str+="\n";}
     }
   str+=StringFormat("[%.5g,%.5g]]},\n",x[n-1],y[n-1]);
   FileWriteString(Fhandle,str);
   Num++;
   return(0);
  }
//+------------------------------------------------------------------+
//| AddGraph                                                         |
//+------------------------------------------------------------------+
int CLinDraw::AddGraph(double &x[],double y,string type,string name,int w=0,string col="")
  {
   int i,k,n;
   string str;

   if(Fhandle<0)return(-1);
   if(Num==0)
     {
      str="$(document).ready(function(){\n"
          "var lp=new Highcharts.Chart({\n"
          "chart:{renderTo:'lplot'},\n"
          "title:{text:'"+Tit+"'},\n"
          "subtitle:{text:'"+SubTit+"'},\n"
          "legend:{enabled:"+Leg+"},\n"
          "yAxis:{title:{text:'"+Ytit+"'}},\n"
          "xAxis:{title:{text:'"+Xtit+"'},showLastLabel:true},\n"
          "series:[\n";
      FileWriteString(Fhandle,str);
     }
   n=ArraySize(x);
   str="{type:'"+type+"',name:'"+name+"',";
   if(col!="")str+="color:'rgba("+col+")',";
   if(w!=0)str+="lineWidth:"+(string)w+",";
   str+="data:[";
   k=0;
   for(i=0;i<n-1;i++)
     {
      str+=StringFormat("[%.5g,%.5g],",x[i],y);
      if(20<k++){k=0; str+="\n";}
     }
   str+=StringFormat("[%.5g,%.5g]]},\n",x[n-1],y);
   FileWriteString(Fhandle,str);
   Num++;
   return(0);
  }
//+------------------------------------------------------------------+
//| LDraw                                                            |
//+------------------------------------------------------------------+
int CLinDraw::LDraw(int ashow=1)
  {
   int i,k;
   string pfnam,to,p[];

   FileWriteString(Fhandle,"]});\n});");
   if(Fhandle<0)return(-1);
   FileClose(Fhandle);

   pfnam=TerminalInfoString(TERMINAL_DATA_PATH)+"\\MQL5\\Files\\"+Fnam;
   k=StringSplit(MQL5InfoString(MQL5_PROGRAM_PATH),StringGetCharacter("\\",0),p);
   to="";
   for(i=0;i<k-1;i++)to+=p[i]+"\\";
   to+="ChartTools\\";                          // Folder name
   DeleteFileW(to+Fnam);
   MoveFileW(pfnam,to+Fnam);
   if(ashow==1)ShellExecuteW(NULL,"open",to+"LinDraw.htm",NULL,NULL,1);
   return(0);
  }

このクラスが適切に処理を行うようにしたければ、JavaScript ライブラリを実装している既製のHTMLページが必要です。グラフを作成するためのフィールドサイズおよび位置を指定することも必要です。そのようなページ実装例は本稿添付のファイルに見ることができます。

AddGraph() メソッドがテキストファイル内で呼ばれると、適切な JavaScript コードが作成されます。そしてテキストファイルが以前に作成された HTML ページにインクルードされます。前述のコードはグラフィックライブラリを参照し、引数として AddGraph() メソッドに渡されたデータを利用してグラフ作成を行います。複数のグラフはこのメソッドが再び呼びだされるときにアウトプットフィールド乗に追加することが可能です。

オーバーロードされたAddGraph() 関数の3バージョンは CLinDraw クラスにインクルードされます。そのバージョンの一つは表示データを持つ配列ただ一つが引数としてそれに変換されることを要求します。それはシーケンスエレメント指数が X 軸に表示されるシーケンスグラフを作成するよう設計されています。

二番目のバージョンは引数として配列を2つ取得します。その配列にはそれぞれ X 軸、Y 軸に表示される値が含まれます。三番目のバージョンは Y 軸に対する定数を設定します。それは水平線を引くのに使用することも可能です。これら関数の残りの引数は似たようなものです。

  • ストリング type – グラフタイプを表示します。値『線』、『スプリーン』、『分散』を持ちます。
  • ストリング name – グラフに割り当てられる名前
  • int w=0 – 線幅オプション引数値 0であれば、デフォルトの幅値2が使用されます。
  • ストリング col="" – 線色およびその不透明値を設定します。オプション引数

LDraw() メソッドは最後に呼びだします。このメソッドは JavaScript コードのアウトプットとデフォルトで \MQL5\Files\ に作成されているテキストファイルにデータを完成します。それからこのファイルはグラフィックライブラリファイルと HTML ファイルを持つディレクトリに移動されます。

LDraw() メソッドは単一のオプション引数を持ちます。引数値が設定されていなかったり1に設定されていると、デフォルトのウェブブラウザが呼びだされ、グラフはデータファイルが適切なディレクトリに移動した後表示されます。引数が他の任意に値であれば、データファイルはまた完全に作成されますがウェブブラウザは自動では呼びだされません。

その他使用可能な CLinDraw クラスメソッドによりグラフのヘッダと軸ラベルを設定することができます。highcharts ライブラリおよびCLinDrawクラスを適用することで生じる問題の詳細については本稿では考察しません。highcharts グラフィックライブラリに関する判りやすい情報は [19]にあります。またグラフやダイアグラムを作成するにあたっての応用例は [9]、[20]を参照ください。

CLinDraw クラスの通常処理に対しては端末での外部ライブラリの使用が許可されています。


8. 密度推定およびファイル整列例

ついに3つのクラスを取得しました。 - CKDensity、CSJPlugin、 CLinDrawです。

以下はクラスの利用法を示したソースコードの例です。

//+------------------------------------------------------------------+
//|                                                  KDE_Example.mq5 |
//|                                                    2012, victorg |
//|                                              https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "2012, victorg"
#property link      "https://www.mql5.com"

#include "CKDensity.mqh"
#include "ChartTools\CLinDraw.mqh"
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   int i,n;
   double a;
   int ndata=1000;                       // Input data length
   double X[];                           // Array for data
   double Z[];                           // Array for graph of a Laplace distribution
   double Sc[];                          // Array for scatter plot

//-------------------------- Preparation of the input sequence
   ArrayResize(X,ndata+1);
   i=CopyOpen(_Symbol,PERIOD_CURRENT,0,ndata+1,X);
   if(i!=ndata+1)
     {
      Print("Not enough historical data.");
      return;
     }
   for(i=0;i<ndata;i++)X[i]=MathLog(X[i+1]/X[i]); // Logarithmic returns
   ArrayResize(X,ndata);

//-------------------------- Kernel density estimation
   CKDensity *kd=new CKDensity;
   kd.PluginMode(1);                     // Enable Plug-in mode
   kd.Density(X);                        // Density estimation 

//-------------------------- Graph of a Laplace distribution
   n=kd.Np;                              // Number of test point
   ArrayResize(Z,n);
   for(i=0;i<n;i++)
     {
      a=MathAbs(kd.T[i]*M_SQRT2);
      Z[i]=0.5*MathExp(-a)*M_SQRT2;
     }
//-------------------------- Scatter plot
   n=kd.N;                               // Data length
   if(n<400)
     {
      ArrayResize(Sc,n);
      for(i=0;i<n;i++)Sc[i]=kd.X[i];
     }
   else                                  // Decimation of data
     {
      ArrayResize(Sc,400);
      for(i=0;i<400;i++)
        {
         a=i*(n-1.0)/399.0+0.5;
         Sc[i]=kd.X[(int)a];
        }
     }

//-------------------------- Visualization
   CLinDraw *ld=new CLinDraw;
   ld.Title("Kernel Density Estimation");
   ld.SubTitle(StringFormat("Data lenght=%i, h=%.3f",kd.N,kd.H));
   ld.YTitle("density");
   ld.XTitle("normalized X (mean=0, variance=1)");

   ld.AddGraph(kd.T,Z,"line","Laplace",1,"200,120,70,1");
   ld.AddGraph(kd.T,kd.Y,"line","Estimation");
   ld.AddGraph(Sc,-0.02,"scatter","Data",0,"0,4,16,0.4");

   ld.LDraw();                      // With WEB-browser autostart 
//   ld.LDraw(0);                        // Without autostart

   delete(ld);
   delete(kd);
  }

確率密度関数の判定を行うデータは KDE_Example.mq5 スクリプトの冒頭に準備されます。これを行うには、現シンボルおよび期間の『オープン』値を X[] 配列にコピーします。そうするとそれらを基に対数利益率が計算されます。

次のステップは CKDensity クラスのコピー作成です。その PluginMode() メソッドを呼び出すこと h レンジ判定にSJPIメソッドを使用することが可能となります。それから X[] 配列で作成されたシーケンスに対し密度推定がが行われます。判定プロセスはこのステップで完了です。以降のステップはすべて取得した密度推定値の可視化への取り組みです。

取得した推定値は任意の有名な分布法則と比較する必要があります。このためにラプラス分布法則に対応する値を Z[] 配列に作成します。正規化されソートされたインプットデータ Sc[] 配列に格納されます。入力シーケンス長が400エレメントを超える場合は、そのすべての値が Sc[]にインクルードされるとは限りません。中にはスキップされるものもあります。よって Sc[] 配列サイズは400エレメントを越えないということです。

表示に必要なデータがすべて出そろったらCLinDraw クラスのコピーが作成されます。次に、ヘッダが指定され、必要なグラフが AddGraph() メソッドを使って追加されます。

最初のグラフはテンプレートとしてのラプラス分布法則グラフです。続いてはインプットデータを使用して取得される密度推定のグラフです。一番下のソースデータグループを表示しているグラフは最後に追加されます。表示に必要なエレメントがすべて決定されたら LDraw() メソッドが呼び出されます。

結果、図7に示されるグラフができます。

図7　対数利益率に対する密度推定値（USDJPY、日次）

図7　対数利益率に対する密度推定値（USDJPY、日次）

図7に示される推定は USDJPYの日次に対する対数利益率の値1,000について行われたものです。見てのとおり、推定値はラプラス分布に相当する曲線にひじょうによく似ています。

密度推定および可視化に必要なファイルはすべて下記のアーカイブ KDEFiles.zip にあります。アーカイブには以下のファイルも含まれています。

  • DensityEstimation\ChartTools\CLinDraw.mqh – highcharts.jsと連携するクラス
  • DensityEstimation\ChartTools\highcharts.js - Highcharts JS v2.2.0 ライブラリ(2012-02-16)
  • DensityEstimation\ChartTools\jquery.js – jQuery v1.7.1 ライブラリ
  • DensityEstimation\ChartTools\LinDraw.htm – 表示用HTML ページ
  • DensityEstimation\CKDensity.mqh – 密度推定を実装するクラス
  • DensityEstimation\CSJPlugin.mqh – 最適レンジ値推定の SJPI メソッドを実装するクラス
  • DensityEstimation\KDE_Example.mq5 – 対数利益率に対する密度推定の例

アーカイブ解凍後、すべてのコンテンツを伴う DensityEstimation\ ディレクトリが端末の Scripts（または Indicators）ディレクトリにセットされる必要があります。それから KDE_Example.mq5をすぐにコンパイルし起動することは可能です。必要であればDensityEstimation は名前を付けなおすこともできます。それがプログラム処理に影響を与えることはありません。

CLinDraw クラスの通常処理に対しては端末での外部ライブラリの使用が許可されていることを繰り返します。

DensityEstimation\ ディレクトリは推定結果を可視化するのに必要なすべてのコンポーネントを持つ ChartTools\ サブディレクトリがインクルードされています。可視化ファイルは個別のサブディレクトリにセットされ、 DensityEstimation\ ディレクトリのコンテンツは簡単に確認されます。サブディレクトリChartTools\ の名前はソースコードで直接指定します。よってそれの名前を付けなおすことはありません。あるいはソースコードに変更を加える必要があります。

テキストファイルは可視化中作成されることはすでに述べました。このファイルにはグラフ作成に必要なデータが含まれています。このファイルはもともと端末の特殊な Files ディレクトリで作成されます。ただそのときそれは前述の ChartTools ディレクトリに移動されます。そのため Files\ あるいはその他いかなる端末上のディレクトリにテスト例処理に関わる痕跡はまったく残されていないことになります。よって本稿からインストールしたファイルを完全に消去するにはすべてのコンテンツと共にディレクトリ DensityEstimation を削除すればよいだけです。


おわりに

いくつかの特殊なメソッドを説明する数式は本稿にインクルードされています。本稿を読みやすくするため意図的にそういたしました。必要であればすべての計算は数多くの文献で見つけることが可能です。その何点かに対するリンクは下記に提供しています。

本稿に記述のあるソースコードは本格的検証を経たものではありません。よってそれらは例にしかすぎず、完全に機能するアプリケーションではないことをお伝えしておきます。


参照資料

  1. Histogram.
  2. Kernel density estimation.
  3. Kernel smoother.
  4. Expectation–maximization algorithm.
  5. А.V. Batov, V.Y. Korolyov, А.Y. Korchagin. EM-Algorithm Having a Large Number of Components as a Means of Constructing Non-Parametric Density Estimations (in Russian).
  6. Hardle W. Applied Nonparametric Regression.
  7. Kernel (statistics).
  8. Charts and diagrams in HTML.
  9. Simon J. Sheather. (2004). Density Estimation.
  10. Max Kohler, Anja Schindler, Stefan Sperlich. (2011). A Review and Comparison of Bandwidth Selection Methods for Kernel Regression.
  11. Clive R. Loader. (1999). Bandwidth Selection: Classical or Plug-In?.
  12. S. J. Sheather, M. C. Jones. (1991). A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation.
  13. Newton's method.
  14. Data binning.
  15. Changjiang Yang, Ramani Duraiswami and Larry Davis. (2004). Efficient Kernel Machines Using the Improved Fast Gauss Transform.
  16. Fast optimal bandwidth estimation for univariate KDE.
  17. R.J. Karunamuni and T. Alberts. (2005). A Locally Adaptive Transformation Method of Boundary Correction in Kernel Density Estimation.
  18. Highcharts JS.
  19. Analysis of the Main Characteristics of Time Series.


MetaQuotes Ltdによってロシア語から翻訳されました。
元の記事: https://www.mql5.com/ru/articles/396

添付されたファイル |
ZIPをダウンロード
kdefiles.zip (82.31 KB)
警告: これらの資料についてのすべての権利はMetaQuotes Ltd.が保有しています。これらの資料の全部または一部の複製や再プリントは禁じられています。

この著者による他の記事

ディスカッションに移動
統計の基礎 統計の基礎
たとえファンダメンタル分析支持者であったとしても、すべてのトレーダーは、特定の統計的な計算を使用し作業を行います。この記事は、統計の基礎、基礎的な要素を紹介し、意思決定における統計の重要性を示します。
MQL5.community での新規記事発表システム MQL5.community での新規記事発表システム
MQL5.community での新規記事発表システムをご紹介します。新規システムでは、いくつかのステップに分けることで記事執筆の全手順を明確で安心して行えるものとしました。各ステップでは役に立つアドバイスとお薦め事項を提供します。それは記事執筆経験の真髄のようなものです。本稿がみなさんの数ある疑問の答えを見つける一助となることを願っています。MQL5.community の訪問者の中でみなさんを人気者にするような興味をひかれる新しい材料をお寄せください。
2013 年第三四半期 MetaTrader AppStore 実績 2013 年第三四半期 MetaTrader AppStore 実績
また四半期が経過したところで、 MetaTrader AppStore の実績を集計することにしました。MetaTrader AppStore は MetaTrader の売買ロボットおよびテクニカルインディケータの最大ストアです。報告対象四半期の終わりまでに「マーケット」には 500 人以上の開発者が 1,200 以上のプロダクツを出しました。
非加法的統計分布構造解析への固有座標法の応用 非加法的統計分布構造解析への固有座標法の応用
応用統計学の主な問題は受け入れる統計仮説の問題ですそれはながらく解決不可能と考えられていました。しかし、固有座標法の登場により状況は一変しました。それはシグナルの構造学にとってすぐれた力強いツールであり、近代的な応用統計学手法を用いることで、可能なこと以上のものを見極めることができるようになります。本稿はこの手法の実践的使用に着目し、MQL5によるプログラムの説明をします。また Hilhorst と Schehrによって紹介される分布例を交えて関数同定の問題も取り上げます。