2012/05/27

モンテカルロフィルタのサンプルを作ってみた。

TwitterのTLに粒子フィルタと言う言葉が出てきたので調べてみました。呼び方がいろいろある様で、モンテカルロフィルタとか、パーティクルフィルタなどと呼ばれているようです。今回は、粒子をベクトル化させない(変化量を与えていない)のでモンテカルロフィルタとしました。

どんなもの?

解りやすくまとめてあるサイトが数多くありますので、参照してください。
手順は、以下の様になります。
  1. モデルの決定
  2. 粒子の初期化(散布)
  3. 予測
  4. 尤度の計算、重み付
  5. リサンプリング
  6. 重み付平均
  7. 3に戻る

作成

今回作成したプログラムの基本設計は、以下のとおりです。また、直前の値からほとんど動かない予測値を想定しています。

【モデル】
システムモデル x[t]=x[-t1]+v (※vは標準正規分布のノイズ)
観測モデル y[t]=x[t]+w (※w=0と設定)
粒子設計 Particles[粒子数][構造] (構造[0:位置情報,1:重み情報])
尤度の計算 標準正規分布を使用
正規分布の偏差 直近データより取得
リターン値 直近の終値を基準にした予測値(黄線)と答えとなる終値(赤線)

【パラメータ】
int    Particles_vol =1000; 粒子数量を設定するパラメータ
int    Range = 100; 試行範囲(粒子数量×試行回数で計算量が増えるので注意)
int    sigmaperiod =100; 標準正規分布のパラメター(偏差)を取得するための計算範囲
double Set = 3000.0; 粒子の初期設定の一様分布の範囲(pipsで記載) -Set<= Particles[][0] <=Set

【コード】

//+------------------------------------------------------------------+
//|                                                    SMCFilter.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
#define PI  3.1415926
#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 ans[];
double Particles[][2]; 
//+------------------------------------------------------------------+
//| 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);
   //乱数の初期化 
   MathSrand(TimeLocal());
   ArrayResize(Particles,Particles_vol);
//----
   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(){
   double rd = MathRand();
   rd/= 32768;
   return(rd);
}
//+------------------------------------------------------------------+
//|正規分布を計算 mu:平均 sigma:偏差                                |
//+------------------------------------------------------------------+
double normal(double x,double sigma){
   double mu = 0;
   return(MathExp(-MathPow((x-mu),2)/(2*sigma*sigma))/(MathSqrt(2*PI)*sigma));
}
//+------------------------------------------------------------------+
//|正規分布となる乱数の成形 mu:平均 sigma:偏差                       |
//+------------------------------------------------------------------+
double normalrand(double sigma){
     double mu = 0;
     double t,u,r1,r2;
     sigma = MathSqrt(sigma);
     t=MathSqrt(-2.0*MathLog(1-random()));
     u=2*PI*random();
     r1=t*MathCos(u);
     return(r1*sigma+mu);
}
//+------------------------------------------------------------------+
//|リサンプリング 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);
   //limit =Set*Point;
   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);
  }
//+------------------------------------------------------------------+

テスト

稼働させると以下の様になります。
5598
今回のテストでは、直前の終値から次の終値は、ほとんど変化しないという条件で作成したため、このような結果となりました。

まとめ

今回は、線形・ガウス(正規分布)のモデルを使用しましたが、この手法は、非線形・非ガウスへ対応ができる手法の様です。(参考)(・・? 参考までに。。