ランダムフォレストで遊ぶ

はじめに

簡単だけど性能がよく、様々な実装が公開されていてマジでパナいと噂の、ランダムフォレストで遊んでみる。

ランダムフォレストとは

  • Breimanによって発展改良された、複数の相関の低い決定木を組み合わせる集団学習の一つ
    • 詳細な紹介や内容は「参考」を参照
  • これ自体は、枠組み(フレームワーク)的な感じが強い
  • 単純な場合、以下のようなパラメータがある
    • 決定木の個数
    • 決定木で使用する学習データの割合
    • 決定木の種類
      • 決定木の深さの制限
      • 決定木の各ノードで使用する判別関数・基準
      • 決定木で使用する素性の割合
      • など
  • 各決定木間の相関が低くなるよう、いろんなところにランダム性を取り入れている
    • 逆に相関が高い場合は、みんな同じような結果を出力しやすいので、みんな間違えてると意味がない
  • また、各決定木は独立しているので、並列処理できる

いろんな実装

本家Breimanによる実装(Fortran)、OpenCV、scikit-learnなどたくさん公開されているし、簡単に使える。

ただし、「Random Forests」(とそれに近い単語)は計算ソフトの商標登録されているらしい。

コード

中身の理解のために、簡単なコードを書いて挙動を追ってみた。
(効率などは考えていない、いつもどおり適当なコード・・・)

パラメータは以下で試してみた。

  • 決定木の個数はB
  • 決定木の学習データの個数はN
  • 決定木の深さ制限はDまで
  • 素性は学習データに出てくる種類の1/3だけ使う
  • 各ノードでは、各素性1つについて、ランダムに閾値を試し、情報利得の大きいものを採用する
    • ランダムな閾値は、th_itr回、th_minからth_maxの間から一様に値を決めて試す
#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <map>
#include <cmath>

//xorshift
// 注意: longではなくint(32bit)にすべき
unsigned long xor128(){
  static unsigned long x=123456789, y=362436069, z=521288629, w=88675123;
  unsigned long t;
  t=(x^(x<<11));
  x=y; y=z; z=w;
  return w=(w^(w>>19))^(t^(t>>8));
}
//[0,1)の一様乱数
// 注意: int_maxぐらいで割った方がよい
double frand(){
  return xor128()%10000000000/static_cast<double>(10000000000); 
}


//決定木
class DecisionTree {
  const static int INF = 999999999;
  
  int M; //素性の種類数
  int m; //この決定木で使用する素性数
  double p; //使用する素性の割合(m = M*p)
  int D; //決定木の深さ制限
  int th_itr; //しきい値の試行回数(th_itr回乱数でしきい値を決め、確認する)
  double th_min, th_max; //ランダムにしきい値を選ぶときの値の範囲
  std::vector< std::pair<int,double> > nodes; //素性,しきい値
  std::vector< std::pair<int,int> > nodes_flag; //フラグ(1:葉,-1:中間ノード),クラスラベル
  
  std::vector< int > feats; //素性のリスト(シャッフルし、前からm個だけを利用する)
  
  //[min_val,max_val)の一様乱数
  double rand_range(double min_val, double max_val){
    return (max_val - min_val) * frand() - min_val;
  }

  //一番数の多いクラスのラベルを返す
  int major_class(const std::vector<int>& t){
    std::map<int,int> cnt;
    for(int i=0; i<t.size(); i++){
      cnt[ t[i] ]++;
    }
    int best = -1, best_cnt = -1;
    for(std::map<int,int>::iterator itr=cnt.begin(); itr!=cnt.end(); ++itr){
      if(best_cnt < itr->second){
        best = itr->first;
        best_cnt = itr->second;
      }
    }
    return best;
  }

  //エントロピー
  double info(const std::vector<int>& t){
    std::map<int,int> cnt;
    for(int i=0; i<t.size(); i++){
      cnt[ t[i] ]++;
    }
    double ret = 0.0;
    for(std::map<int,int>::iterator itr=cnt.begin(); itr!=cnt.end(); ++itr){
      double p = itr->second / (double)(t.size());
      ret += p * log2(p);
    }
    return -ret;
  }

  //情報利得
  double informationGain(int f, double th, const std::vector<int>& t, const std::vector< std::vector< std::pair<int,double> > >& x){
    std::vector<int> lt, rt;
    std::vector< std::vector< std::pair<int,double> > > lx, rx;
    //素性fの値thで2つにデータを分割する
    for(int i=0; i<t.size(); i++){
      double val = 0.0;
      for(int j=0; j<x[i].size(); j++){
        if(x[i][j].first == f) val = x[i][j].second;
      }
      if(val < th){
        lt.push_back(t[i]);
        lx.push_back(x[i]);
      }else{
        rt.push_back(t[i]);
        rx.push_back(x[i]);
      }
    }
    //情報利得の計算
    double G = info(t); //ここはなくても良さそう
    G -= (lt.size()/(double)t.size()) * info(lt);
    G -= (rt.size()/(double)t.size()) * info(rt);
    return G;
  }

  //ノード番号nでの学習
  // 葉ノードでは、p(t|x)を返すべきだが、ここでは一番多いクラスをそのノードのクラスとする
  void train_n(int n, const std::vector<int>& t, const std::vector< std::vector< std::pair<int,double> > >& x){
    int ln = 2*n+1, rn = 2*n+2; //左の子ノードの番号、右の子ノードの番号

    //深さ制限で、子ノードが作れないので、途中で打ち切る
    if(ln >= nodes.size() || rn >= nodes.size()){
      nodes_flag[n] = std::make_pair(1, major_class(t));
      return;
    }

    //現在のデータセットを2つに分割する
    double best_I = -INF;
    int best_n_f = -INF;
    double best_n_th = -INF;
    for(int f=0; f<m; f++){
      for(int i=0; i<1; i++){
        //検証する素性としきい値
        int n_f = feats[f];
        double n_th = rand_range(th_min,th_max);
        //この条件での情報利得
        double I = informationGain(n_f, n_th, t, x);
        if(best_I < I){
          best_I = I;
          best_n_f = n_f;
          best_n_th = n_th;
        }
      }
    }

    //一番よかった条件を使って分岐を作成
    nodes[n] = std::make_pair(best_n_f, best_n_th);

    //次のノードで学習を繰り返す
    std::vector<int> lt, rt;
    std::vector< std::vector< std::pair<int,double> > > lx, rx;
    for(int i=0; i<t.size(); i++){
      double val = 0.0;
      for(int j=0; j<x[i].size(); j++){
        if(x[i][j].first == best_n_f) val = x[i][j].second;
      }
      if(val < best_n_th){
        lt.push_back(t[i]);
        lx.push_back(x[i]);
      }else{
        rt.push_back(t[i]);
        rx.push_back(x[i]);
      }
    }

    if(lt.size() > 0 && rt.size() > 0){
      train_n(ln, lt, lx);
      train_n(rn, rt, rx);
    }
    else if(lt.size() == 0){
      nodes_flag[rn] = std::make_pair(1, major_class(rt));
    }
    else if(rt.size() == 0){
      nodes_flag[ln] = std::make_pair(1, major_class(lt));
    }
  }

public:
  //決定木のデフォルトパラメータ
  DecisionTree(double p_ = 0.3, int D_ = 10, int th_itr_ = 1, double th_min_ = 0.0, double th_max_ = 1.0):
    p(p_),
    D(D_),
    th_itr(th_itr_),
    th_min(th_min_),
    th_max(th_max_),
    nodes((1<<D)-1),
    nodes_flag((1<<D)-1, std::make_pair(-INF,0))
  {}

  //学習
  void train(const std::vector<int>& t, const std::vector< std::vector< std::pair<int,double> > >& x){
    //素性のリストを作成
    std::map<int,int> feats_;
    for(int i=0; i<t.size(); i++){
      for(int j=0; j<x[i].size(); j++){
        feats_[ x[i][j].first ]++;
      }
    }
    M = feats_.size();
    m = (int)(M * p);
    for(std::map<int,int>::iterator itr=feats_.begin(); itr!=feats_.end(); ++itr){
      feats.push_back(itr->first);
    }
    //相関が低くなるように、素性リストはシャッフルして、前からm個だけをこの決定木では利用する
    std::random_shuffle(feats.begin(), feats.end());
    
    //各ノードで再帰的に学習していく
    train_n(0, t, x);
  }

  //予測
  int predict(const std::vector< std::pair<int,double> >& x){
    int n = 0;
    while(nodes_flag[n].first==-INF){
      double val = 0.0;
      for(int j=0; j<x.size(); j++){
        if(x[j].first == nodes[n].first) val = x[j].second;
      }

      if(val < nodes[n].second){
        n = 2*n+1;
      }else{
        n = 2*n+2;
      }
    }
    return nodes_flag[n].second;
  }
};


//RandomForest
class RandomForest {
  int B; //決定木の数
  int M; //各決定木で使用する学習データの数
  std::vector<DecisionTree> dts; //決定木のインスタンス
  
public:
  RandomForest(int B_, int M_):
    B(B_),
    M(M_),
    dts(B_)
  {}

  void train(const std::vector<int>& t, const std::vector< std::vector< std::pair<int,double> > >& x){
    for(int i=0; i<B; i++){
      //ブートストラップサンプリング
      std::vector<int> tt(M);
      std::vector< std::vector< std::pair<int,double> > > xx(M);
      for(int j=0; j<M; j++){
        int r = xor128()%(t.size());
        tt[j] = t[r];
        xx[j] = x[r];
      }

      //決定木を作成
      std::cerr << "training tree: " << i << std::endl;
      dts[i].train(tt, xx);
    }
  }

  int predict(const std::vector< std::pair<int,double> >& x){
    std::map<int,int> ret;

    //各決定木の結果を取得
    for(int i=0; i<B; i++){
      int res = dts[i].predict(x);
      ret[res]++;
    }

    //多数決で結果を決定
    int best = -1, best_cnt = -1;
    for(std::map<int,int>::iterator itr=ret.begin(); itr!=ret.end(); ++itr){
      if(best_cnt < itr->second){
        best = itr->first;
        best_cnt = itr->second;
      }
    }

    return best;
  }
};


//libsvm形式のファイルの読み込み
void load_file(std::ifstream& fin,
               std::vector<int>& t,
               std::vector< std::vector< std::pair<int,double> > >& x){
  int c;
  int name;
  double value;
  std::string line, feat;
  
  while(std::getline(fin, line)){
    std::stringstream ss(line);
    ss >> c;
    t.push_back(c);
    std::vector< std::pair<int,double> > tmp;
    while(ss >> feat){
      std::string::size_type idx = feat.find(":");
      name = atoi(feat.substr(0,idx).c_str());
      value = atof(feat.substr(idx+1).c_str());
      tmp.push_back( std::make_pair(name, value) );
    }
    x.push_back(tmp);
  }
}

int main(){
  int B = 30; //決定木の数
  int N = 15000; //決定木で使用する学習データの数

  RandomForest rf(B, N);

  //データの読み込み用(ファイル指定)
  std::ifstream train_file("./a9a");
  std::ifstream test_file("./a9a.t");


  //学習データの読み込み
  std::vector<int> t;
  std::vector< std::vector< std::pair<int,double> > > x;
  load_file(train_file, t, x);

  //ランダムフォレストの構築
  rf.train(t, x);
  
  //学習データで評価
  {
    int cnt = 0;
    for(int i=0; i<t.size(); i++){
      int res = rf.predict(x[i]);
      if(t[i] == res) cnt++;
    }
    std::cout << "close Acc = " << (cnt/(double)t.size()) << " (" << cnt << "/" << t.size() << ")" << std::endl;
  }

  //テストデータの読み込み
  std::vector<int> testt;
  std::vector< std::vector< std::pair<int,double> > > testx;
  load_file(test_file, testt, testx);

  //テストデータで評価
  {
    int cnt = 0;
    for(int i=0; i<testt.size(); i++){
      int res = rf.predict(testx[i]);
      if(testt[i] == res) cnt++;
    }
    std::cout << "Acc = " << (cnt/(double)testt.size()) << " (" << cnt << "/" << testt.size() << ")" << std::endl;
  }

  return 0;
}

結果

いくつか試してみる。パラメータは適当に、Accを見てみる。
http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/から。

a9a
    • 学習データ : a9a
    • テストデータ : a9a.t
(B,N,D,th_itr,th_min,th_max)=(30,15000,10,1,0,1)
close Acc = 0.838856 (27314/32561)
Acc = 0.835637 (13605/16281)
dna
    • 学習データ : dna.scale
    • テストデータ : dna.scale.t
(B,N,D,th_itr,th_min,th_max)=(30,1500,10,1,0,1)
close Acc = 0.9885 (1977/2000)
Acc = 0.930017 (1103/1186)
mnist
    • 学習データ : mnist.scale
    • テストデータ : mnist.scale.t
(B,N,D,th_itr,th_min,th_max)=(20,20000,10,5,0,1)
close Acc = 0.936367 (56182/60000)
Acc = 0.9275 (9275/10000)

結構いい感じみたいだけど、素性数が多いと、各ノードで試す回数が増えて遅い。
(news20は重すぎて試していない。)
文書分類だと10000次元とか余裕で超えそうな気がするので、工夫が必要。