GSLでヒストグラム■■■

前回の内容に、ヒストグラムを追加しました。ファイルはココにまとめてあります。

DLL追加コード

view plainprint?

1. #include ".\include\gsl\gsl_histogram.h"//ヒストグラム  
2.   
3. //+-------------------------------------------------------------------+  
4. //|ヒストグラム (Histgram)  
5. //|input[inp_N]  入力用バッファ  
6. //|out[N]      出力用バッファ  
7. //|inp_N         入力データ数  
8. //|N             階級数  
9. //|min           最小値(含む）  
10. //|max           最大値(含まない）  
11. //+-------------------------------------------------------------------+  
12. __declspec(dllexport) void __stdcall Histgram(double *input,double *out,int inp_N,int N,double min,double max){  
13.     int i;  
14.     //メモリの確保  
15.     gsl_histogram * h = gsl_histogram_alloc(N);  
16.     if(h==NULL)exit(0);  
17.     //階級のセット  
18.     gsl_histogram_set_ranges_uniform(h,min,max);  
19.     //実行  
20.     for(i=0;i<inp_N;i++)gsl_histogram_increment(h,input[i]);  
21.       
22.     //出力  
23.     for(i=0;i<N;i++)out[i]=gsl_histogram_get(h,i);  
24.       
25.   
26.     //メモリの解放  
27.     gsl_histogram_free(h);  
28. }  

#include ".\include\gsl\gsl_histogram.h"//ヒストグラム

//+-------------------------------------------------------------------+
//|ヒストグラム (Histgram)
//|input[inp_N]  入力用バッファ
//|out[N]      出力用バッファ
//|inp_N         入力データ数
//|N             階級数
//|min           最小値(含む）
//|max           最大値(含まない）
//+-------------------------------------------------------------------+
__declspec(dllexport) void __stdcall Histgram(double *input,double *out,int inp_N,int N,double min,double max){
	int i;
	//メモリの確保
	gsl_histogram * h = gsl_histogram_alloc(N);
	if(h==NULL)exit(0);
	//階級のセット
	gsl_histogram_set_ranges_uniform(h,min,max);
	//実行
	for(i=0;i<inp_N;i++)gsl_histogram_increment(h,input[i]);
	
	//出力
	for(i=0;i<N;i++)out[i]=gsl_histogram_get(h,i);
	

	//メモリの解放
	gsl_histogram_free(h);
}

MQLコード

//+------------------------------------------------------------------+
//|                                               ヒストグラム.mq4 　|
//|                                       Copyright ｩ 2012, bighope  |
//|                      http://expertadviser-bighope.blogspot.com/  |
//+------------------------------------------------------------------+
#property copyright "Copyright ｩ 2012, bighope"
#property link      "http://expertadviser-bighope.blogspot.com/"
#property indicator_separate_window
#property indicator_minimum 0
#include <GSL.mqh>
#property indicator_buffers 3
#property indicator_color1 Snow
#property indicator_color2 Yellow
#property indicator_color3 Red
extern string  SymbA   =   "USDJPY";
extern int     TframeA =   PERIOD_H1;
extern string  SymbB   =   "EURUSD";
extern int     TframeB =   PERIOD_H1;
extern string  SymbC   =   "AUDUSD";
extern int     TframeC =   PERIOD_H1;
extern int     Range   =   2000;//データ数
extern int     N       =    100;//階級数
extern double  min     =    0.0;//階級最小値(pips)
extern double  max     =    2000.0;//階級最大値(pips)

double Ap[];
double Bp[];
double Cp[];
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
       SetIndexStyle(0,DRAW_LINE);
       SetIndexStyle(1,DRAW_LINE);
       SetIndexStyle(2,DRAW_LINE);

       SetIndexBuffer(0,Ap);
       SetIndexBuffer(1,Bp);
       SetIndexBuffer(2,Cp);

//----
   return(0);
  }
//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit()
  {
//----

//----
   return(0);
  }
//+------------------------------------------------------------------+
//|GetVal
//|symb        通貨ペア名
//|timeframe   タイムフレーム
//|shift       Barの位置
//+------------------------------------------------------------------+
double GetVal(string symb,int timeframe,int shift){
   double Hi,Lw;

   Hi = iHigh(symb,timeframe,shift);
   Lw = iLow(symb,timeframe,shift);

return(Hi-Lw);
}
//+------------------------------------------------------------------+
//|GetData
//|symb        通貨ペア名
//|timeframe   タイムフレーム
//|out[Ns]     データの取得配列
//|Ns          データの取得数
//+------------------------------------------------------------------+
void GetData(string symb,int timeframe,double &out[],int Ns){
   for(int i=0;i<Ns;i++)out[i]=GetVal(symb,timeframe,i);
}
//+------------------------------------------------------------------+
//|GetHist
//|symb        通貨ペア名
//|timeframe   タイムフレーム
//|out[Ns]     データの取得配列
//|Ns          データの取得数
//+------------------------------------------------------------------+
void GetHist(string symb,int timeframe,int Ns,double &data[]){
   int i;
   double count,minpt,maxpt,pt;
   double input[];
   double output[];
   ArrayResize(input,Ns);
   ArrayResize(output,N);
   ArrayInitialize(data,EMPTY_VALUE);

   GetData(symb,timeframe,input,Ns);
   pt = MarketInfo(symb,MODE_POINT);
   minpt = min*pt;
   maxpt = max*pt;
   Histgram(input,output,Ns,N,minpt,maxpt);
   for(i=0;i<N;i++)count += output[i];
   if(count==0)count=1.0;
   for(i=0;i<N;i++)data[i]=output[i]/count;
}

//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start()
  {
   int    counted_bars=IndicatorCounted();
//----
   if(counted_bars<1){
      GetHist(SymbA,TframeA,Range,Ap);
      GetHist(SymbB,TframeB,Range,Bp);
      GetHist(SymbC,TframeC,Range,Cp);
   }

//----
   return(0);
  }
//+------------------------------------------------------------------+

作動・確認

デフォルトで作動させると以下の様になります。
※縦軸が、発生確率で横軸が階級です。（チャート右が最少数）(白：USDJPY 1H 黄：EURUSD 1H 赤：AUDUSD 1H)

※お化粧が全くなされていないので見苦しいかもしれませんがご了承ください。＾＾；

USDJPYが特異な特徴と示す？

いろいろ試してみるとUSDJPYだけがなんだか特異な形状を示します。なぜ？
例えば、(白：EURJPY 1H 黄：EURUSD 1H 赤：USDCHF 1H)だと

の様に似たような形になります。
他のペアも試してみましたが、試したペアの中でUSDJPYだけが極端に形状が異なっています。

また、(白：USDJPY 4H 黄：EURUSD 1H 赤：GBPJPY 1H)とした場合は、

となり、似たような形が出現します。

しかし、(白：USDJPY 1D 黄：EURUSD 1D 赤：AUDUSD 1D)となると、

特異な形状が解消されたように見えます。（ちなみに、max=5000で表示)

 

まとめ

全てのペア及びタイムフレームを確認したわけでありませんが、ドル円だけがなぜ特異な形となるのか疑問です。他ペアに比べて突発的事情（イベントやニュース）に敏感に反応（突発的なボラ）し、その他のボラは、低い事を示しているのでしょうか？疑問が残るところです。

GSLで多項式回帰を試してみた■■■

今回は、前回の続きとして、多項式回帰を試してみました。

まえがき

今回試してみたのは、Ｙ＝Ｃ０＋Ｃｎ×Ｘ＾Ｃｎの多項式回帰です。
ｎ値を変えることで以下の様に変化します。
ｎ＝２の時

ｎ＝３の時


ｎ＝４の時


ｎ＝１０の時

サンプル

サンプルコードは、ここです。
※上記ｎ値は、変数名Ｐｗです。

ＤＬＬコード

※前回のコードをそのまま残し、修正しました。

view plainprint?

1. #pragma once  
2. #define WIN32_LEAN_AND_MEAN   
3. #include <Windows.h>  
4. #include <stdlib.h>  
5. #include <stdio.h>  
6. #include <math.h>//数学用  
7. #include ".\include\gsl\gsl_sort.h"//配列並び替え用  
8. #include ".\include\gsl\gsl_wavelet.h"//ウェーブレット用  
9. #include ".\include\gsl\gsl_multifit.h"//重回帰  
10.   
11. BOOL APIENTRY DllMain(HANDLE hModule,DWORD ul_reason_for_call,LPVOID lpReserved)  
12. {  
13. //----  
14. switch(ul_reason_for_call)  
15. {  
16.  case DLL_PROCESS_ATTACH:  
17.  case DLL_THREAD_ATTACH:  
18.  case DLL_THREAD_DETACH:  
19.  case DLL_PROCESS_DETACH:  
20.  break;  
21. }  
22. //----  
23. return(TRUE);  
24. }  
25.   
26. //+-------------------------------------------------------------------+  
27. //|重回帰（Multifit_liner)Y=c0+c1X+c2X^2+c3X^3...cpwX^pw  
28. //|input[]   入力用バッファ  
29. //|pw        変数量  
30. //|out[pw+1]    出力用バッファ[c0,c1,c2,c3,....cpw.chisq]  
31. //|Ns       対象データ数  
32. //+-------------------------------------------------------------------+  
33. __declspec(dllexport) void __stdcall Multifit_liner(double *input,int pw,double *out,int Ns){  
34.  int i,ipw;  
35.  double    ix,chisq;  
36.  gsl_matrix *X, *cov;  
37.  gsl_vector *y, *w, *c;  
38.    
39.  //配列（ベクトル配列）の確保  
40.  X = gsl_matrix_alloc(Ns, pw);  
41.  y = gsl_vector_alloc(Ns);  
42.  c = gsl_vector_alloc(pw);  
43.  cov = gsl_matrix_alloc(pw, pw);  
44.   
45.  //配列（ベクトル配列）入力  
46.  for(i=0;i<Ns;i++){  
47.   ix = i*0.001;  
48.   gsl_matrix_set (X, i, 0, 1.0);  
49.   for(ipw=1;ipw<pw;ipw++)gsl_matrix_set (X, i, ipw, pow(ix,ipw));   
50.   gsl_vector_set (y, i, input[i]);  
51.  }  
52.     
53.  //処理  
54.  gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (Ns, pw);  
55.  gsl_multifit_linear(X,y,c,cov,&chisq,work);  
56.  gsl_multifit_linear_free (work);  
57.    
58.  //抽出  
59.  for(i=0;i<pw;i++)out[i]=gsl_vector_get(c,(i));  
60.  out[pw]=chisq;  
61.   
62.  //メモリ解放  
63.  gsl_matrix_free(X);  
64.  gsl_matrix_free(cov);  
65.  gsl_vector_free(y);  
66.  gsl_vector_free(c);  
67.   
68. }  
69.   
70.   
71. //+-------------------------------------------------------------------+  
72. //|ウェーブレット変換（DWT）  
73. //|input[]  入力用バッファ  
74. //|out[]    出力用バッファ  
75. //|Ns       対象データ数（２のNs乗)  
76. //|D        ウェーブレット係数（2:Harr 4~daubechies  
77. //|fs       フィルタ  
78. //+-------------------------------------------------------------------+  
79. __declspec(dllexport) void __stdcall DWT(double *input,double *out,int Ns,int D,int fs){  
80.  int i;  
81.  int siz = (int)pow(2.0,Ns);  
82.  if(D%2==1)D++;  
83.  if(D>20)D=20;  
84.    
85.  //配列のコピー  
86.  memcpy(out,input,sizeof(double)*siz);   
87.    
88.  //ウェーブレットの初期設定(ウェーブレット名,D値）  
89.  gsl_wavelet *w;  
90.  if(D==2){  
91.   w = gsl_wavelet_alloc(gsl_wavelet_haar,D);  
92.  }else{  
93.   w = gsl_wavelet_alloc(gsl_wavelet_daubechies, D);  
94.    }  
95.     
96.  //作業領域の確保(数量）  
97.  gsl_wavelet_workspace *work;  
98.  work = gsl_wavelet_workspace_alloc(siz);  
99.  size_t *p = (size_t *)malloc(siz * sizeof(size_t));  
100.  double *abs = (double *)malloc(siz * sizeof (double));  
101.  if (abs == NULL || p== NULL){printf("メモリを確保できません"); exit(0);}  
102.    
103.  //変換(初期設定,データ,進み幅,作業エリア）  
104.  gsl_wavelet_transform_forward(w, out, 1, siz,work);  
105.       
106.  //フィルタ作業  
107.  for(i = 0; i < siz; i++) abs[i] = fabs(out[i]);  
108.  gsl_sort_index(p, abs, 1, siz);  
109.  for(i = 0; (i + fs) < siz; i++) out[p[i]] = 0;  
110.       
111.  //逆変換  
112.  gsl_wavelet_transform_inverse(w, out, 1, siz, work);  
113.   
114.     //メモリ解放  
115.  free(abs);  
116.  free(p);  
117.  gsl_wavelet_workspace_free(work);  
118.  gsl_wavelet_free(w);  
119. }  
120. //----------------------------------------------------------------------------  

#pragma once
#define WIN32_LEAN_AND_MEAN 
#include <Windows.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>//数学用
#include ".\include\gsl\gsl_sort.h"//配列並び替え用
#include ".\include\gsl\gsl_wavelet.h"//ウェーブレット用
#include ".\include\gsl\gsl_multifit.h"//重回帰

BOOL APIENTRY DllMain(HANDLE hModule,DWORD ul_reason_for_call,LPVOID lpReserved)
{
//----
switch(ul_reason_for_call)
{
 case DLL_PROCESS_ATTACH:
 case DLL_THREAD_ATTACH:
 case DLL_THREAD_DETACH:
 case DLL_PROCESS_DETACH:
 break;
}
//----
return(TRUE);
}

//+-------------------------------------------------------------------+
//|重回帰（Multifit_liner)Y=c0+c1X+c2X^2+c3X^3...cpwX^pw
//|input[]   入力用バッファ
//|pw        変数量
//|out[pw+1]    出力用バッファ[c0,c1,c2,c3,....cpw.chisq]
//|Ns       対象データ数
//+-------------------------------------------------------------------+
__declspec(dllexport) void __stdcall Multifit_liner(double *input,int pw,double *out,int Ns){
 int i,ipw;
 double    ix,chisq;
 gsl_matrix *X, *cov;
 gsl_vector *y, *w, *c;
 
 //配列（ベクトル配列）の確保
 X = gsl_matrix_alloc(Ns, pw);
 y = gsl_vector_alloc(Ns);
 c = gsl_vector_alloc(pw);
 cov = gsl_matrix_alloc(pw, pw);

 //配列（ベクトル配列）入力
 for(i=0;i<Ns;i++){
  ix = i*0.001;
  gsl_matrix_set (X, i, 0, 1.0);
  for(ipw=1;ipw<pw;ipw++)gsl_matrix_set (X, i, ipw, pow(ix,ipw)); 
  gsl_vector_set (y, i, input[i]);
 }
  
 //処理
 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (Ns, pw);
 gsl_multifit_linear(X,y,c,cov,&chisq,work);
 gsl_multifit_linear_free (work);
 
 //抽出
 for(i=0;i<pw;i++)out[i]=gsl_vector_get(c,(i));
 out[pw]=chisq;

 //メモリ解放
 gsl_matrix_free(X);
 gsl_matrix_free(cov);
 gsl_vector_free(y);
 gsl_vector_free(c);

}


//+-------------------------------------------------------------------+
//|ウェーブレット変換（DWT）
//|input[]  入力用バッファ
//|out[]    出力用バッファ
//|Ns       対象データ数（２のNs乗)
//|D        ウェーブレット係数（2:Harr 4~daubechies
//|fs       フィルタ
//+-------------------------------------------------------------------+
__declspec(dllexport) void __stdcall DWT(double *input,double *out,int Ns,int D,int fs){
 int i;
 int siz = (int)pow(2.0,Ns);
 if(D%2==1)D++;
 if(D>20)D=20;
 
 //配列のコピー
 memcpy(out,input,sizeof(double)*siz); 
 
 //ウェーブレットの初期設定(ウェーブレット名,D値）
 gsl_wavelet *w;
 if(D==2){
  w = gsl_wavelet_alloc(gsl_wavelet_haar,D);
 }else{
  w = gsl_wavelet_alloc(gsl_wavelet_daubechies, D);
   }
  
 //作業領域の確保(数量）
 gsl_wavelet_workspace *work;
 work = gsl_wavelet_workspace_alloc(siz);
 size_t *p = (size_t *)malloc(siz * sizeof(size_t));
 double *abs = (double *)malloc(siz * sizeof (double));
 if (abs == NULL || p== NULL){printf("メモリを確保できません"); exit(0);}
 
 //変換(初期設定,データ,進み幅,作業エリア）
 gsl_wavelet_transform_forward(w, out, 1, siz,work);
    
 //フィルタ作業
 for(i = 0; i < siz; i++) abs[i] = fabs(out[i]);
 gsl_sort_index(p, abs, 1, siz);
 for(i = 0; (i + fs) < siz; i++) out[p[i]] = 0;
    
 //逆変換
 gsl_wavelet_transform_inverse(w, out, 1, siz, work);

    //メモリ解放
 free(abs);
 free(p);
 gsl_wavelet_workspace_free(work);
 gsl_wavelet_free(w);
}
//----------------------------------------------------------------------------

まとめ

適用範囲を広くとりすぎたり、ｎ値を大きくしすぎるとバグが発生する可能性がありますので注意してください。