2012/05/31

GSLの乱数と確率密度関数を使ってみる。

GSLの乱数と確率密度関数は、条件を付ければラッパーを使わなくてもそのまま使えそうなので試してみました。
※条件とは、GSLには、多くの乱数器がありますが、初期設定の乱数器を使用するという条件です。

ヘッダーの作成

代表的なものだけピックアップし、ヘッダーファイルを作成しました。

//+------------------------------------------------------------------+
//|                                                     gsl_rand.mq4 |
//|                        Copyright 2012, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright 2012, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
#property show_inputs

#import "gsl.dll"
//乱数の基本設定
   //環境設定の取得
   int gsl_rng_env_setup ();
   //型T の乱数発生器のインスタンスを生成して、そのハンドルを返す。
   //T:gsl_rng_env_setupで取得したハンドル。
   int gsl_rng_alloc (int T);
   //乱数発生器を初期化
   //r:gsl_rng_allocで取得したハンドル
   //s:初期化種
   void gsl_rng_set (int r, datetime s);
   //乱数発生器のインスタンス r に割り当てられているメモリを解放する。
   void gsl_rng_free (int r);

//以下 r:gsl_rng_allocで取得したハンドル
//乱数

   //0<=x<1の乱数を返す。
   double gsl_rng_uniform (int r);
   //0<x<1の乱数を返す。
   double gsl_rng_uniform_pos (int r);

//分布に従う乱数と確率密度関数
  
   //正規分布乱数
   double gsl_ran_gaussian (int r, double sigma);
   //正規分布の確率密度関数
   double gsl_ran_gaussian_pdf (double x, double sigma);
   
   //平均値mu の指数分布にしたがう乱数を返す。
   double gsl_ran_exponential (int r, double mu);
   //x における平均値mu の指数分布の確率密度関数
   double gsl_ran_exponential_pdf (double x, double mu);
   
   //幅a のラプラス分布にしたがう乱数を返す。
   double gsl_ran_laplace (int r, double a);
   //x における幅a のラプラス分布の確率密度関数
   double gsl_ran_laplace_pdf (double x, double a);
   
   //底がa、指数がb のべき指数分布にしたがう乱数を返す。
   double gsl_ran_exppow (int r, double a, double b);
   //x における底がa、指数がb のべき指数分布の確率密度関数
   double gsl_ran_exppow_pdf (double x, double a, double b);
   
   //係数a のコーシー分布にしたがう乱数を返す。
   double gsl_ran_cauchy (int r, double a);
   //x における係数a のコーシー分布の確率密度関数
   double gsl_ran_cauchy_pdf (double x, double a);
   
   //ガンマ分布にしたがう乱数を返す。
   double gsl_ran_gamma (int r, double a, double b);
   //引数がa とb のときの、x におけるガンマ分布の確率密度関数
   double gsl_ran_gamma_pdf (double x, double a, double b);
   
   //a からb の間に一様に分布する乱数を返す。
   double gsl_ran_flat (int r, double a, double b);
   //a からb の間の一様分布の確率密度関数
   double gsl_ran_flat_pdf (double x, double a, double b);
   
   //対数正規分布にしたがう乱数を返す。
   double gsl_ran_lognormal (int r, double zeta, double sigma);
   //パラメータ値がzeta とsigma で与えられたときの、対数正規分布の確率密度関数
   double gsl_ran_lognormal_pdf (double x, double zeta, double sigma);
   
   //自由度n のカイ二乗分布にしたがう乱数を返す。
   double gsl_ran_chisq (int r, double nu);
   //自由度n のカイ二乗分布の確率密度関数
   double gsl_ran_chisq_pdf (double x, double nu);
   
   //自由度nu1、nu2 のF 分布にしたがう乱数を返す。
   double gsl_ran_fdist (int r, double nu);
   //自由度nu1、nu2 のF 分布の確率密度関数
   double gsl_ran_fdist_pdf (double x, double nu);
   
   //t 分布にしたがう乱数を返す。
   double gsl_ran_tdist (int r, double nu);
   //自由度nu のt 分布の確率密度関数
   double gsl_ran_tdist_pdf (double x, double nu );
   
   //ベータ分布にしたがう乱数を返す。
   double gsl_ran_beta (int r, double a, double b );
   //パラメータ値がa とb で与えられたとき、ベータ分布の確率密度関数
   double gsl_ran_beta_pdf (double x, double a, double b);
   
   //ロジスティック分布にしたがう乱数を返す。
   double gsl_ran_logistic (int r, double a);
   //係数a が与えられたとき<鴻Wスティック分布の確率密度関数
   double gsl_ran_logistic_pdf (double x, double a);
#import


使い方

先回のコードに移植してみました。※変更ヶ所は、赤字部分です。
gsl.dllとgslcblas.dllをterminal.exeと同じ位置に置いています。

//+------------------------------------------------------------------+
//|                                                SMCFilter_gsl.mq4 |
//|                        Copyright 2012, MetaQuotes Software Corp. |
//|                                        http://www.metaquotes.net |
//+------------------------------------------------------------------+
#property copyright "Copyright 2012, MetaQuotes Software Corp."
#property link      "http://www.metaquotes.net"
#property indicator_separate_window
//#property indicator_chart_window
#include <gsl_rand.mqh>

#property indicator_buffers 2
#property indicator_color1  Yellow
#property indicator_width1  1
#property indicator_color2  Red
#property indicator_width2  1

extern int    Particles_vol = 1000;//粒子の数
extern int    Range         =  100;//試行範囲
extern int    sigmaperiod   =  100;//偏差取得期間
extern double Set           = 3000.0;//初期化範囲(pips)


double Linebuf[];
double Particles[][2];
double ans[];
int r;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   SetIndexStyle(0, DRAW_LINE);
   SetIndexBuffer(0,Linebuf);
   SetIndexEmptyValue( 0,0);
   
   SetIndexStyle(1, DRAW_LINE);
   SetIndexBuffer(1,ans);
   SetIndexEmptyValue( 1,0);
   
   ArrayResize(Particles,Particles_vol);
   
   //乱数の初期化 
   int T=gsl_rng_env_setup();
   r = gsl_rng_alloc(T);
   gsl_rng_set (r,TimeLocal());
//----
   return(0);
  }
//+------------------------------------------------------------------+
//| Custom indicator deinitialization function                       |
//+------------------------------------------------------------------+
int deinit()
  {
//----
   gsl_rng_free (r);
//----
   return(0);
  }
//+------------------------------------------------------------------+
//|粒子の初期化                                                      |
//+------------------------------------------------------------------+  
void setparticles(double & pobj[][],int range,double set){
   for(int i=0;i<Particles_vol;i++){
      Particles[i][0] = data(range)+(((set*2)/Particles_vol)*i-set)*Point;
      Particles[i][1] = 0; 
    }
 }
//+------------------------------------------------------------------+
//|乱数の成形  0<=random<1                                           |
//+------------------------------------------------------------------+
double random(){
   return(gsl_rng_uniform (r));
}
//+------------------------------------------------------------------+
//|正規分布を計算 mu:平均 sigma:偏差                                |
//+------------------------------------------------------------------+
double normal(double x,double sigma){
   return(gsl_ran_gaussian_pdf(x,sigma));
}
//+------------------------------------------------------------------+
//|正規分布となる乱数の成形 mu:平均 sigma:偏差                       |
//+------------------------------------------------------------------+
double normalrand(double sigma){     
     return(gsl_ran_gaussian (r,sigma));
}
//+------------------------------------------------------------------+
//|リサンプリング q:粒子数                                           |
//+------------------------------------------------------------------+
void resample(int q,double &pobj[][] ){
   
   double wt[];
   double rd;
   int i,j;
   
   //粒子のコピー
   double copyp[][2];
   ArrayResize(copyp,q);
   ArrayCopy(copyp,pobj,0,0,q*2);
   
   //累計重み
   ArrayResize(wt,q);
   wt[0] = pobj[0][2];
   for(i=1;i<q;i++) wt[i]=wt[i-1]+pobj[i][2];
   
   //リサンプリング(ルーレット法)
   for(i=0;i<q;i++){
      rd = random()/Particles_vol;
      for(j=0;j<q;j++){
         if(rd > wt[j]){
            continue;
         }else{
            pobj[i][0] = copyp[j][0];
            pobj[i][1] = 0; 
            break;
         }
      }   
   }
}
//+------------------------------------------------------------------+
//|尤度計算および重み付け                                            |
//+------------------------------------------------------------------+
void wlikelihood(double price,int p,double sigma,double &pobj[][]){
   double alfa[];
   double wsum;
   double g = 1;    //Gの偏微分の値
   int i;
   ArrayResize(alfa,p);
   for(i=0;i<p;i++){
      alfa[i] = normal((price-pobj[i][0]),sigma)*g;//(尤度)
      wsum += alfa[i];
   }
   if(wsum<1)wsum=1;
   //正規化
   for(i=0;i<p;i++)pobj[i][1] = alfa[i]/wsum;
}
//+------------------------------------------------------------------+
//|対象データ                                                        |
//+------------------------------------------------------------------+
double data(int index){
   return(Close[index]);
}
//+------------------------------------------------------------------+
//|偏差の取得                                                        |
//+------------------------------------------------------------------+
double setsigma(int range,int index){
   int      i;
   double sum,av;
   for(i=0;i<range;i++)av += data(index+i);
   av /= range; 
   for(i=0;i<range;i++)sum += MathPow((data(index+i)-av),2);
   return(MathSqrt(sum/range)); 
}
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int start(){
   double v;//乱数
   double sum,sigma,limit;
   int i,j;
   static datetime now = 0;
   int    counted_bars=IndicatorCounted();
   if(counted_bars<1)i=Range; else i=0;
   
   if(now==0)setparticles(Particles,Range,Set);
   while(i>=0){
      //各bar一回のみ実行
      if(now!=Time[i]){
         now=Time[i];
         sigma=setsigma(sigmaperiod,i);
         for(j=0;j<Particles_vol;j++){
            Particles[j][0] = Particles[j][0] + normalrand(sigma);      
         }

         //尤度推定          
         wlikelihood(data(i+1),Particles_vol,sigma,Particles);
         //リサンプリング
         resample(Particles_vol,Particles);   
            
         //加重付平均
         sum=0.0;              
         for(j=0;j<Particles_vol;j++)sum += (Particles[j][0])*Particles[j][1];
      
         Linebuf[i] = sum;
      }
      ans[i] = data(i);//答え
   i--;
   }
   return(0);
  }
//+------------------------------------------------------------------+


まとめ

とりあえず乱数と密度関数で時間を取られることはなくなりそうです。