パターン認識

(7) 決定木 ( Decision Tree )

木構造は、探索・ソート・圧縮などの処理で利用される重要なデータ構造です。パターンを分類するのに木構造が利用されることもあり、このときの木構造を「決定木 ( Decision Tree )」といいます。決定木は分類の流れが視覚化されるので比較的解釈がしやすいという特徴があります。今回は、決定木について紹介します。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 決定木 ( Decision Tree )

「決定木 ( Decision Tree )」は学習パターン X = { x1, x2, ... xn } と教師ベクトル t = ( t1, t2, ... tn )T を使って生成されます。まず、最初の根となるノード ( Root Node ) に全ての学習パターンと教師ベクトルを置きます。次に、学習パターンにある同種の要素 { xji } ( i = 1, 2, ... n ) を使っていくつかの集合に学習パターンと教師ベクトルを分類します。このとき、各分類結果の教師信号ができるだけ等しいものどうし集まるようにします。分類方法としては、xji が性別 ( 男性・女性 ) や国名 ( 日本・中国・アメリカなど ) の「名義尺度 ( Nominal Scale )」なら各カテゴリの組み合わせを利用し、年齢や測定値などの連続値である場合は適当なしきい値を決めて分割するなどが考えられます。分類したデータは根の下に子ノード ( Child Node ) としてぶら下げます。もし、ノード内の教師信号がすべて等しくなったり、これ以上学習パターンによる分類ができなくなったりした場合は、分類処理を打ち切り教師ベクトルの代表値を決定します。このときのノードは「葉ノード ( Leaf Node )」になります。そうでなければ、同様の分類処理を子ノードに対して繰り返します。このようにして決定木が生成されます。

実例を見た方が理解しやすいと思いますので、以下のデータを実際に分類してみます。

表 1-1. 20 人の選んだ果物 ( いちご・りんご・バナナ )
No.性別年齢果物
1男性15りんご
2男性18バナナ
3男性18りんご
4男性21いちご
5男性23りんご
6男性26いちご
7男性29いちご
8男性33バナナ
9男性36バナナ
10男性37バナナ
11女性17いちご
12女性18いちご
13女性20いちご
14女性21いちご
15女性21いちご
16女性24いちご
17女性29いちご
18女性33りんご
19女性36いちご
20女性39いちご

上表は、異なる年齢層の男女それぞれ 10 名に、いちご・りんご・バナナから好きな果物を選択してもらった結果を表した架空のデータです。このとき、性別と年齢は学習パターン、果物は教師信号になります。表を見ると、女性の 9 割はいちごを選択しているので、まずは性別で二つの子ノードに分類します。男性の方はバラバラに選んでいるように見えますが、10 代、20 代、30 代と年齢で分類してみると 10 代はりんご、20 代はいちご、30代はバナナが最も多く選ばれています。女性の方は 30 代にひとりだけりんごを選択した人がいますが、30 代でくくればいちごの方が多くなるので、これ以上は分類しないことにします。したがって、分類木は以下のようになります。

図 1-1. データを決定木にした結果
決定木

この例はデータ数が少ないので簡単に分類することができましたが、実際に扱うデータは通常もっと多く、それを目で見て処理することはまずできません。コンピュータで機械的に処理するとしても、何らかのアルゴリズムを用意する必要があります。決定木を生成するためのアルゴリズムにはたくさんの種類があり、それぞれを組み合わせて利用される場合もあります。今回はその中から代表的なものを紹介します。

2) ID3 アルゴリズム ( ID3 Algorithm )

「ID3 アルゴリズム ( Iterative Dichotomiser 3 (ID3) Algorithm )」は 1979 年に「ジョン・ロス・キンラン ( John Ross Quinlan )」によって提案されたアルゴリズムです。学習パターンに複数の種別があったとき ( 前節の例では性別と年齢 )、どれを使って分類するのが最適なのかを検討するためには何らかの指標が必要になります。この指標のことを「不純度 ( Impurity )」といいます。ID3 アルゴリズムでは不純度として「エントロピー ( Entropy )」を用います。
エントロピーは「ハフマン符号化」の中でも紹介したことがあります。全事象 Ω の中で各事象 x が出現する確率を P(x) としたとき、エントロピーは

H(Ω) = Σx{x∈Ω}( -P(x) log2P(x) )

で求めることができます。但し、P(x) = 0 のときは -P(x) log2P(x) = 0 とします。二つの事象の出現確率がそれぞれ 0, 1 だったとき、H = 0 になります。それに対してどちらも 1 / 2 なら H = 1 です。このように、出現確率に偏りがなくなるほどエントロピーは増加します。前節の例において、教師信号である果物の各出現確率はいちごが 12 / 20 = 3 / 5、りんごとバナナが 4 / 20 = 1 / 5 なのでエントロピーは

H = -( 3 / 5 ) log2( 3 / 5 ) - ( 1 / 5 ) log2( 1 / 5 ) - ( 1 / 5 ) log2( 1 / 5 ) = 1.37

になります。性別で分類したとき、果物の各出現確率は男性の場合いちご 3 / 10、りんご 3 / 10、バナナ 4 / 10 なので

Hmale = -( 3 / 10 ) log2( 3 / 10 ) - ( 3 / 10 ) log2( 3 / 10 ) - ( 4 / 10 ) log2( 4 / 10 ) = 1.57

女性の場合いちご 9 / 10、りんご 1 / 10 なので

Hfemale = -( 9 / 10 ) log2( 9 / 10 ) - ( 1 / 10 ) log2( 1 / 10 ) = 0.47

となって、男性・女性の出現比率はどちらも 1 / 2 なので、平均値は

( 1 / 2 )Hmale + ( 1 / 2 )Hfemale = 1.02

となります。よって、分類前のエントロピー H との差は 0.35 という結果が得られます。この差のことを「情報利得 ( Information Gain )」といいます ( 補足 1 )。エントロピーは出現率の偏り度合いを表すものでした。したがって、情報利得は分類によって偏りが増えた、すなわち同じものどうしが固まって分類された度合いを表します。

年齢を使って 10 代、20 代、30 代に分類した場合、各々のエントロピー H1, H2, H3 は次のようになります。

H1 = -( 2 / 5 ) log2( 2 / 5 ) - ( 2 / 5 ) log2( 2 / 5 ) - ( 1 / 5 ) log2( 1 / 5 ) = 1.52

H2 = -( 8 / 9 ) log2( 8 / 9 ) - ( 1 / 9 ) log2( 1 / 9 ) = 0.50

H3 = -( 2 / 6 ) log2( 2 / 6 ) - ( 1 / 6 ) log2( 1 / 6 ) - ( 3 / 6 ) log2( 3 / 6 ) = 1.46

出現率は 10 代 5 / 20 = 1 / 4、20 代 9 / 20、30 代 6 / 20 = 3 / 10 なので、平均値は

( 1 / 4 )H1 + ( 9 / 20 )H2 + ( 3 / 10 )H3 = 1.04

となって、情報利得は 0.33 となります。よって、性別で分類したほうがわずかながら情報利得は大きく、偏りがより大きくなることを意味します。ID3 アルゴリズムでは、このエントロピーから得られる情報利得を使い、その大小によって分類方法を決定します。

生成された子ノードに対して同じような方法でさらに分類を続けていくと、やがて教師信号が全て等しくなり分類が不要になるか、分類に利用できる学習パターンがなくなります。前者の場合は唯一になった教師信号を登録し、後者の場合は教師信号が複数ある可能性があるので代表値を選んで登録して処理を停止します。このとき、教師信号が登録されたノードは葉ノードとなります。教師信号が複数ある場合の代表値の選択方法としては、最も出現率の高いものを抽出するのが一般的です。

ID3 アルゴリズムのサンプル・プログラムを以下に示します。まずはエントロピーと情報利得を計算する部分です。

/*
  GetSum : 変数ごとのカウント値の合計を返す

  count : 変数ごとのカウント値を保持する配列
*/
template< class T >
unsigned int GetSum( const std::map< T, unsigned int >& count )
{
  typedef std::map< T, unsigned int > Map;

  unsigned int sum = 0;
  for ( typename Map::const_iterator i = count.begin() ;
        i != count.end() ; ++i )
    sum += i->second;

  return( sum );
}

/*
  Impurity : 不純度と総数を返す

  テンプレート引数の T は変数の型を表す
  テンプレート引数の Op は不純度を計算する式を表す

  count : 変数をキーとするカウント
  op : 不純度を計算する式
*/
template< class T, class Op >
std::pair< double, unsigned int >
Impurity( const std::map< T, unsigned int >& count, Op op )
{
  // 総数を計算
  unsigned int total = GetSum( count );
  if ( total == 0 )
    return( std::pair< double, unsigned int >( 0, 0 ) );

  // 不純度の計算
  double imp = 0;
  for ( typename std::map< T, unsigned int >::const_iterator i = count.begin() ;
        i != count.end() ; ++i ) {
    double p = static_cast< double >( i->second ) / total;
    if ( p != 0 )
      imp += op( p );
  }

  return( std::pair< double, unsigned int >( imp, total ) );
}

/*
  Entropy : 交差エントロピーを返す

  p : 発生確率
*/
double Entropy( double p )
{
  return( -p * std::log( p ) / std::log( 2 ) );
}

/*
  GetMode : カウント値が最大の変数(最頻値)を返す

  count : 変数ごとのカウント値を保持する配列
*/
template< class T >
T GetMode( const std::map< T, unsigned int >& count )
{
  typedef std::map< T, unsigned int > Map;
  typedef std::pair< T, unsigned int > Pair;

  typename Map::const_iterator i = count.begin();
  Pair max = *i;
  while ( ++i != count.end() ) {
    if ( i->second > max.second )
      max = *i;
  }

  return( max.first );
}

/*
  ImpByEntropy : 不純度に交差エントロピーを用いる場合の関数オブジェクト

  テンプレート引数の T は変数の型を表す
*/
template< class T >
struct ImpByEntropy
{
  // 交差エントロピーと総数を返す
  std::pair< double, unsigned int >
  operator()( const std::map< T, unsigned int >& count ) const
  {
    return( Impurity( count, Entropy ) );
  }

  // 要約値として最頻値を返す
  T summary( const std::map< T, unsigned int >& count ) const
  {
    return( GetMode( count ) );
  }
};

不純度の計算は関数 Impurity で行われます。引数として渡される count は連想配列 map の形式をとり、教師信号をキーとした要素数としています。したがって、あらかじめキーごとに集計をしておかなければなりません。集計のための関数は後ほど紹介します。もう一つの引数 op は不純度の計算式です。ID3 アルゴリズムでは不純度としてエントロピーが用いられますが、他のアルゴリズムでは異なる不純度が利用されるので、使い回しができるように外部から式を渡すことができるようにしてあります。実際の計算式は関数 Entropy にて実装しています。Impurity の戻り値はエントロピーと要素数の総計のペアとしています。総計は処理の中で必要となるのでこのような形となっています。
関数 GetMode は、最大の要素数を持つ教師信号を得るための関数です。教師信号が複数存在する状態で処理を停止なければならないときは、この関数を使って代表値としての教師信号を得ます。最後に、Entoropy を利用して Impurity を呼び出す形で operator() を実装し、GetMode を呼び出す形で summary をメンバ関数に持つクラス ImpByEntropy を定義します。

ID3 アルゴリズム本体の実装は以下のようになります。

/*
  DecisionTreeNode : 決定木のノード
*/
template< class X, class Y >
struct DecisionTreeNode
{
  typename Matrix< X >::size_type index; // 列番号
  Y target;           // 教師信号
  double impurity;    // 不純度
  unsigned int level; // ノードのレベル(高さ)

  // 枝のノードへのリンク
  std::map< X, typename std::vector< DecisionTreeNode >::size_type > branch;
};

/*
  ID3_Algorithm : ID3 アルゴリズムによる決定木の学習

  テンプレート引数の CATEGORY は学習パターン、TARGET は教師信号の型をそれぞれ表す
  テンプレート引数の Imp で不純度の種類を選択することができる
*/
template< class CATEGORY, class TARGET, template< class > class Imp = ImpByEntropy >
class ID3_Algorithm
{
  // 型の定義
  typedef typename Matrix< CATEGORY >::size_type Index; // データ参照用インデックスの型

  typedef DecisionTreeNode< CATEGORY, TARGET > DecisionTreeNode_; // 決定木のノード
  typedef std::vector< DecisionTreeNode_ > DecisionTree; // 決定木
  typedef typename DecisionTree::size_type size_type;    // 決定木の参照用インデックス

  typedef std::map< TARGET, unsigned int > MapTarget;  // 教師信号に対するカウント値配列
  typedef std::map< CATEGORY, MapTarget > MapCategory; // 学習パターン・教師信号に対するカウント値配列

  typedef std::pair< TARGET, unsigned int > PairTarget;  // 教師信号に対するカウント値要素
  typedef std::pair< CATEGORY, MapTarget > PairCategory; // 学習パターン・教師信号に対するカウント値要素

  // メンバ変数
  DecisionTree tree_; // 決定木
  Imp< TARGET > op_;  // 不純度を計算する関数オブジェクト

  // 情報利得を返す
  double getGain( const MapCategory& count, const MapTarget& targetCount );

  // 最大情報利得を持つ列の番号を求める
  Index getMaxCol( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
                   const std::vector< Index >& rowIndex, const std::vector< Index >& colIndex,
                   const MapTarget& targetCount );

  // 分岐するノードを追加する
  void addNode( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
                size_type currentNode,
                const std::vector< Index >& rowIndex,
                const std::vector< Index >& colIndex );

 public:

  // 学習パターン x と教師信号 y を指定して構築
  ID3_Algorithm( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y );

  // 学習パターンを指定して予測値を返す
  template< class Ran >
  std::pair< TARGET, bool > predict( Ran s ) const;
};

/*
  Classify : 学習パターン・教師信号ごとに分類を行う

  テンプレート引数の X は学習パターン、Y は教師信号の型をそれぞれ表す

  x : 学習パターン
  y : 教師信号
  rowIndex : 読み込み対象の x, y に対する行番号
  colIndex : 読み込み対象の x に対する列番号
  count : 学習パターン・教師信号ごとのカウント値を保持する配列
*/
template< class X, class Y >
void Classify
  ( const Matrix< X >& x, const std::vector< Y >& y,
    const std::vector< typename Matrix< X >::size_type >& rowIndex,
    typename Matrix< X >::size_type colIndex,
    std::map< X, std::map< Y, unsigned int > >* count )
{
  typedef std::map< Y, unsigned int > MapY;
  typedef std::map< X, MapY > MapXY;

  typedef std::pair< Y, unsigned int > PairY;
  typedef std::pair< X, MapY > PairXY;

  typedef typename Matrix< X >::size_type Index;

  count->clear();

  for ( typename std::vector< Index >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i ) {
    X c = x[*i][colIndex]; // 学習パターンの値
    // 学習パターンをキーとする要素への参照を取得
    typename MapXY::iterator cit = count->find( c );
    if ( cit == count->end() )
      cit = ( count->insert( PairXY( c, MapY() ) ) ).first;
    // 取得した参照から教師信号をキーとする要素への参照を取得してカウントをインクリメント
    typename MapY::iterator tit = ( cit->second ).find( y[*i] );
    if ( tit == ( cit->second ).end() )
      tit = ( ( cit->second ).insert( PairY( y[*i], 0 ) ) ).first;
    ++( tit->second );
  }
}

/*
  Classify : 教師信号ごとに分類を行う

  テンプレート引数の Y は教師信号の型を表す

  y : 教師信号
  rowIndex : 読み込み対象の y に対する行番号
  count : 教師信号ごとのカウント値を保持する配列
*/
template< class Y >
void Classify
  ( const std::vector< Y >& y,
    const std::vector< typename std::vector< Y >::size_type >& rowIndex,
    std::map< Y, unsigned int >* count )
{
  typedef std::map< Y, unsigned int > MapY;
  typedef std::pair< Y, unsigned int > PairY;
  typedef typename std::vector< Y >::size_type Index;

  count->clear();

  for ( typename std::vector< Index >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i ) {
    typename MapY::iterator tit = count->find( y[*i] );
    if ( tit == count->end() )
      tit = ( count->insert( PairY( y[*i], 0 ) ) ).first;
    ++( tit->second );
  }
}

/*
  ID3_Algorithm< CATEGORY, TARGET, Imp >::getGain : 情報利得を返す

  count : 学習パターンを第一キー、教師信号を第二キーとするカウント
  targetCount : 教師信号ごとのカウント
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
double ID3_Algorithm< CATEGORY, TARGET, Imp >::getGain
  ( const MapCategory& count, const MapTarget& targetCount )
{
  typedef std::pair< double, unsigned int > Pair;

  // 教師信号ごとに和をとった不純度
  Pair targetImp = op_( targetCount );

  // 学習パターンごとに不純度を求め、全体との比率を掛ける
  double imp = 0;
  for ( typename MapCategory::const_iterator i = count.begin() ;
        i != count.end() ; ++i ) {
    Pair p = op_( i->second );
    imp += p.first * p.second / targetImp.second;
  }

  return( targetImp.first - imp );
}

/*
  ID3_Algorithm< CATEGORY, TARGET, Imp >::getMaxCol : 最大情報利得を持つ列の番号を求める

  x : 学習パターン
  y : 教師信号
  rowIndex : 読み込み対象の x, y に対する行番号
  colIndex : 読み込み対象の x に対する列番号
  targetCount : 教師信号ごとのカウント値を保持する配列
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
typename ID3_Algorithm< CATEGORY, TARGET, Imp >::Index
ID3_Algorithm< CATEGORY, TARGET, Imp >::getMaxCol
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    const std::vector< Index >& rowIndex, const std::vector< Index >& colIndex,
    const MapTarget& targetCount )
{
  // 学習パターン・教師信号ごとのカウント
  MapCategory count;

  // 最初の列で最大情報利得を初期化
  typename std::vector< Index >::const_iterator c = colIndex.begin(); // 列番号
  Classify( x, y, rowIndex, *c, &count );
  double maxG = getGain( count, targetCount ); // 最大情報利得
  Index maxCol = *c; // 最大情報利得を持つ列の番号

  // 最大情報利得を持つ列の探索
  while ( ++c != colIndex.end() ) {
    Classify( x, y, rowIndex, *c, &count );
    double g = getGain( count, targetCount );
    if ( g > maxG ) {
      maxG = g;
      maxCol = *c;
    }
  }

  return( maxCol );
}

/*
  ID3_Algorithm< CATEGORY, TARGET, Imp >::addNode : 分岐するノードを追加する

  x : 学習パターン
  y : 教師信号
  currentNode : 現在のノード
  rowIndex : 読み込み対象の x, y に対する行番号
  colIndex : 読み込み対象の x に対する列番号
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
void ID3_Algorithm< CATEGORY, TARGET, Imp >::addNode
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    size_type currentNode,
    const std::vector< Index >& rowIndex, const std::vector< Index >& colIndex )
{
  if ( rowIndex.empty() ) return;

  MapTarget targetCount; // 教師信号ごとのカウント
  Classify( y, rowIndex, &targetCount );

  tree_[currentNode].impurity = op_( targetCount ).first;

  // 教師信号が単一なら分岐しないので target に登録して終了
  if ( targetCount.size() < 2 ) {
    tree_[currentNode].target = ( targetCount.begin() )->first;
    return;
  }

  // 対象列がないなら分岐できないので教師信号の代表値を target に登録して終了
  if ( colIndex.empty() ) {
    tree_[currentNode].target = op_.summary( targetCount );
    return;
  }

  // 情報利得が最大となる列番号を取得する
  Index maxCol = getMaxCol( x, y, rowIndex, colIndex, targetCount );

  // 学習パターンが等しい行の番号列に分ける
  std::map< CATEGORY, std::vector< Index > > splitRow;
  for ( typename std::vector< Index >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i ) {
    typename std::map< CATEGORY, std::vector< Index > >::iterator it = splitRow.find( x[*i][maxCol] );
    if ( it == splitRow.end() )
      it = splitRow.insert( std::pair< CATEGORY, std::vector< Index > >( x[*i][maxCol], std::vector< Index >() ) ).first;
    ( it->second ).push_back( *i );
  }

  // 行が分割されないなら分岐できないので教師信号の代表値を target に登録して終了
  if ( splitRow.size() < 2 ) {
    tree_[currentNode].target = op_.summary( targetCount );
    return;
  }

  tree_[currentNode].index = maxCol;

  // 新たな colIndex を生成する
  std::vector< Index > newColIndex( colIndex.begin(), colIndex.end() );
  newColIndex.erase( std::find( newColIndex.begin(), newColIndex.end(), maxCol ) );

  // 分岐したノードを生成する
  for ( typename std::map< CATEGORY, std::vector< Index > >::const_iterator it = splitRow.begin() ;
        it != splitRow.end() ; ++it ) {
    tree_.push_back( DecisionTreeNode_() );
    tree_.back().level = tree_[currentNode].level + 1;
    tree_[currentNode].branch[it->first] = tree_.size() - 1;
    addNode( x, y, tree_.size() - 1, it->second, newColIndex );
  }
}

/*
  GenerateSerialNumber : 連番生成用関数オブジェクト
*/
class GenerateSerialNumber
{
  // インクリメントしながら値を返す
  template< class T >
  class Inc_
  {
    T index;

  public:

    Inc_() : index( T() ) {}
    T operator()() { return( index++ ); }
  };

 public:

  // 連番を生成する
  template< class In >
  void operator()( In s, In e )
  {
    typedef typename std::iterator_traits< In >::value_type value_type;

    Inc_< value_type > inc;
    std::generate( s, e, inc );
  }
};

/*
  ID3_Algorithm< CATEGORY, TARGET, Imp > コンストラクタ

  x : 学習パターン
  y : 教師信号
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
ID3_Algorithm< CATEGORY, TARGET, Imp >::ID3_Algorithm
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y )
{
  ErrLib::CheckLinearModel( x, y.begin(), y.end() );

  tree_.push_back( DecisionTreeNode_() );
  tree_.back().level = 0;

  // 読み込み対象の行列を連番で初期化
  std::vector< Index > rowIndex( y.size() );
  GenerateSerialNumber()( rowIndex.begin(), rowIndex.end() );
  std::vector< Index > colIndex( x.cols() );
  GenerateSerialNumber()( colIndex.begin(), colIndex.end() );

  addNode( x, y, 0, rowIndex, colIndex );
}

/*
  Predict : 学習パターンを指定して予測値を返す

  指定した学習パターンが決定木に存在しなかった場合は戻り値のペアの第二引数は false になる
*/
template< class CATEGORY, class TARGET, class Ran >
std::pair< TARGET, bool >
Predict( const std::vector< DecisionTreeNode< CATEGORY, TARGET > >& tree, Ran s )
{
  typedef std::vector< DecisionTreeNode< CATEGORY, TARGET > > DecisionTree; // 決定木
  typedef typename DecisionTree::size_type size_type; // 決定木の参照用インデックス

  size_type i = 0;
  while ( ! tree[i].branch.empty() ) {
    typename std::map< CATEGORY, size_type >::const_iterator it = tree[i].branch.find( s[tree[i].index] );
    if ( it == tree[i].branch.end() )
      return( std::pair< TARGET, bool >( TARGET(), false ) );
    i = it->second;
  }

  return( std::pair< TARGET, bool >( tree[i].target, true ) );
}

/*
  ID3_Algorithm< CATEGORY, TARGET, Imp >::predict : 学習パターンを指定して予測値を返す

  テンプレート引数の Ran はランダムアクセス反復子であることを想定している
  指定した学習パターンが決定木に存在しなかった場合は戻り値のペアの第二引数は false になる

  s : 学習パターンの開始位置
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
template< class Ran >
std::pair< TARGET, bool >
ID3_Algorithm< CATEGORY, TARGET, Imp >::predict( Ran s ) const
{
  return( Predict( tree_, s ) );
}

DecisionTreeNode は決定木のノードを表す構造体です。重要なメンバ変数は、分割対象となる学習パターンのある列の番号 index と、教師信号を保持する target、学習パターン別に子ノードの位置を保持する branch の 3 つです。ノードが葉ならば、必要となるのは target のみで、branch が空であることで葉であるかどうかを判別するものとします。葉以外のノードは target は不要になり、index に分割対象の学習パターンの列番号を保持し、学習パターンの内容をキーとして子ノードの位置を map 形式の branch へ登録します。決定木は DecisionTreeNode を要素とする可変長配列 vector の形とするので、子ノードの位置としては添字を持たせることにします。なお、テンプレート引数の X は学習パターンの、Y は教師信号のデータ型をそれぞれ表します。

決定木の生成はクラス ID3_Algorithm の構築時に行うようにします。構築時の引数は学習パターンからなる行列の x と教師信号からなるベクトル y の二つです。行列を表す Matrix はサンプル・プログラムの中では実装していませんが、x[行番号][列番号] で任意の要素が取得できることだけ知っていれば理解できると思います。
クラス ID3_Algorithm は 3 つのテンプレートを持っています。CATEGORY は学習パターンの、TARGET は教師信号のデータ型をそれぞれ表しています。最後の Imp は少し見慣れない書き方をしています。

template< class > class Imp

これは、Imp がテンプレート引数を 1 つ持っていることを示しています。メンバ変数としては Imp< TARGET > op_ が定義されているので、テンプレートの型は TARGET であることがわかります。したがって、宣言時には class の後の型の指定は省略することができます。この Imp の実体は、先に紹介した関数オブジェクト ImpByEntropy ですが、この関数は切り替えることが可能な仕様になっています。もし、テンプレート引数として単に class Imp と書いた場合、インスタンス化するときに例えば次のように記述しなければなりません。

ID3_Algorithm< string, int, ImpByEntropy< int > > id3( x, y );

ImpByEntropy の後のテンプレート引数は冗長であり、誤りが発生する原因にもなります。template< class > class Imp としておくことによって、以下のように記述することができます。

ID3_Algorithm< string, int, ImpByEntropy > id3( x, y );

これならテンプレート引数が何であったのか悩む必要もありません。

コンストラクタの中で呼び出される関数 ErrLib::CheckLinearModel は、渡された学習パターンと教師信号の要素数が一致しているかなどをチェックするための関数ですが、サンプル・プログラムの中では実装はされていません。

実際の決定木の生成はメンバ関数 addNode が行います。まず、y にある教師信号をキーとして各要素数を関数 Classify で集計します。集計結果は教師信号をキー、要素数を値とする map 形式 map< TARGET, unsigned int > で保持されます。Classify にはこの他に対象の行番号を保持する配列 rowIndex が渡されます。要素の分割が進むと集計対象の行番号は変化するので、どの行が集計対象になるのかを引数として渡す必要があります。同様に、一度分割に利用した学習パターンはそれ以降は常に単一となり分割する意味がないため除外できます。対象の列番号を保持する変数として colIndex もあることに注意して下さい。
集計の結果、教師信号が単一であった場合はこれ以上ノードを分割しても意味がないため、現在のノードは葉ノードとして扱い、target に教師信号を登録して処理を終了します。また、colIndex が空ならばこれ以上分割できないので教師信号を登録して終了します。このとき、教師信号の選択肢は複数ある可能性があるため、Imp クラスが持つ summary メンバ関数を使って代表値を取得して登録します。先述の通り、ImpByEntropy の場合は要素数が最大の教師信号が選択されます。

情報利得が最大になる学習パターンの列番号は getMaxCol が取得します。getMaxCol ではまず最初に Classify を使って学習パターン・教師信号の組み合わせごとに各列の要素数 count を集計します。このとき、集計結果を保持するコンテナは map< CATEGORY, map< TARGET, unsigned int > > という形式になります。各組み合わせの要素数は count[学習パターン][教師信号] で取得することができるので、連想配列版の二次元配列と考えれば理解しやすいと思います。こうして集計された結果はメンバ関数 getGain に渡されて情報利得が計算されます。getGain に渡される引数は、count の他に、addNode で集計された教師信号ごとの要素数 targetCount もあることに注意して下さい。これは、情報利得を計算するのに分割前のエントロピーの値も必要となるためです。また、Imp クラスの operator() は不純度の他に要素数の総計もペアで戻り値として返すのでした。要素数は比率計算に必要になるため、このような仕様になっています。getGain で情報利得を求めたら、それが現在の最大値を超えているかをチェックして、超えていたらその値と列番号を保持します。こうして情報利得が最大となる列番号を取得することができます。

情報利得が最大になる学習パターンが決まったら、学習パターンが等しい行どうしで配列を生成します。これは新たな rowIndex として利用されます。ここで行が分割されず一つしか存在しなかったら、分割はできないので教師信号の代表値を target に登録して処理を終了します。さらに、今回の分割で利用した学習パターンの列番号を colIndex から除外します。あとは、生成した新しい rowIndex と colIndex を使って addNode を再帰的に呼び出して処理を続けます。

前節の例をサンプル・プログラムで処理したときの決定木を以下に示します。各ノードには、そのときのエントロピーの値と、葉ノードには選択された教師信号が記されています。

図 2-1. ID3 アルゴリズムによる決定木の生成
ID3 アルゴリズムによる決定木

ID3 アルゴリズムは処理の性格上、学習パターンとして名義尺度しか利用することができず、連続値には対応していません。上記の処理においても、年齢はあらかじめ 10 代、20 代、30 代と連番を振ってから実行しています。また、エントロピーの計算は発生頻度から得られることから、教師信号に連続値を利用することも意味を持ちません。連続値に対して処理を行いたい場合は別のアルゴリズムを検討する必要があります。

3) CART ( Classification And Regression Tree )

「CART ( Classification And Regression Tree )」は 1984 年に「レオ・ブレイマン ( Leo Breiman )」などによって発表されたアルゴリズムです。その名の通り、「分類木 ( Classification Tree )」と「回帰木 ( Regression Tree )」の両方に対応しており、教師信号が前者は名義尺度の、後者は連続値の場合にあたります。また、学習パターンが連続値であっても対応できるという特徴もあります。

まず最初に、学習パターンと教師信号がともに名義尺度である場合の処理方法について説明します。学習パターンに複数の種別があったとき、ID3 アルゴリズムでは単純にそれぞれの種別に分解して集計する方法を取っていましたが、CART の場合は 2 つの集合に分割する組み合わせの中で最も情報利得の大きいものを選択します。例えば、A, B, C の 3 つの種別があったとき、これらを 2 つの集合に分ける組み合わせは

A | B C
B | C A
C | A B

の 3 種類です。それぞれの組み合わせで出現率を計算して情報利得を求め、その最大値を持つ組み合わせを採用するわけです。このようなアルゴリズムのため、CART は ID3 アルゴリズムと異なりノードが必ず二分割となります。

情報利得の計算方法も ID3 アルゴリズムとは少し異なります。ID3 アルゴリズムでは出現率からエントロピーを求めたのに対し、CART は「ジニ不純度 ( Gini Impurity )」を用います。全事象 Ω の中で各事象 x が出現する確率を P(x) としたとき、ジニ不純度は以下の式で求められます。

G(Ω) = Σx{x∈Ω}( P(x)[ 1 - P(x) ] )

エントロピーのときと同様に、二つの事象の出現確率がそれぞれ 0, 1 だったとき G = 0 なのに対し、どちらも 1 / 2 なら G = 1 / 2 なので、出現率に偏りがなくなるほどジニ不純度は増大します。前節の例において、果物の各出現率 ( いちご 3 / 5、りんごとバナナ 1 / 5 ) からジニ不純度は

G = ( 3 / 5 )( 1 - 3 / 5 ) + ( 1 / 5 )( 1 - 1 / 5 ) + ( 1 / 5 )( 1 - 1 / 5 ) = 0.56

で、性別で分類したとき

Gmale = ( 3 / 10 )( 1 - 3 / 10 ) + ( 3 / 10 )( 1 - 3 / 10 ) + ( 4 / 10 )( 1 - 4 / 10 ) = 0.66
Gfemale = ( 9 / 10 )( 1 - 9 / 10 ) + ( 1 / 10 )( 1 - 1 / 10 ) = 0.18

なので、情報利得は 0.56 - (1/2) x 0.66 - (1/2) x 0.18 = 0.14 となります ( 補足 1 )。

以下にサンプル・プログラムを示します。

/*
  GiniCoef : ジニ不純度を返す

  p : 発生確率
*/
double GiniCoef( double p )
{
  return( p * ( 1 - p ) );
}

/*
  ImpByGini : 不純度にジニ不純度を用いる場合の関数オブジェクト

  テンプレート引数の T は変数の型を表す
*/
template< class T >
struct ImpByGini
{
  // ジニ不純度と総数を返す
  std::pair< double, unsigned int >
  operator()( const std::map< T, unsigned int >& count ) const
  {
    return( Impurity( count, GiniCoef ) );
  }

  // 要約値として最頻値を返す
  T summary( const std::map< T, unsigned int >& count ) const
  {
    return( GetMode( count ) );
  }
};

/*
  BinTreeNode : 決定木のノード(二分木)
*/
template< class X, class Y >
struct BinTreeNode
{
  typedef typename Matrix< X >::size_type Index; // 学習パターン参照用インデックスの型

  Index index; // 学習パターンの列番号
  Y target;    // 教師信号(代表値)
  double impurity;    // 不純度
  unsigned int total; // 総カウント値
  std::set< X > lCatg; // 左側の学習パターン
  std::set< X > rCatg; // 右側の学習パターン

  // 枝のノードへのリンク
  typename std::vector< BinTreeNode >::size_type left;
  typename std::vector< BinTreeNode >::size_type right;
};

/*
  Nominal : 名義尺度専用 CART アルゴリズム

  テンプレート引数の CATEGORY は学習パターン、TARGET は教師信号の型をそれぞれ表す
  テンプレート引数の Imp で不純度の種類を選択することができる
*/
template< class CATEGORY, class TARGET, template< class > class Imp = ImpByGini >
class Nominal
{
  // 型の定義
  typedef typename Matrix< CATEGORY >::size_type Index; // データ参照用インデックスの型

  typedef std::map< TARGET, unsigned int > MapTarget;  // 教師信号に対するカウント値配列
  typedef std::map< CATEGORY, MapTarget > MapCategory; // 学習パターン・教師信号に対するカウント値配列

  typedef std::pair< TARGET, unsigned int > PairTarget;  // 教師信号に対するカウント値要素
  typedef std::pair< CATEGORY, MapTarget > PairCategory; // 学習パターン・教師信号に対するカウント値要素

  typedef BinTreeNode< CATEGORY, TARGET > BinTreeNode_; // 二分木のノード
  typedef std::vector< BinTreeNode_ > BinTree;   // 二分木
  typedef typename BinTree::size_type size_type; // 二分木のインデックス

  // メンバ変数
  BinTree tree_;     // 二分木
  Imp< TARGET > op_; // 不純度を計算する関数オブジェクト

  // カテゴリの組み合わせに対する情報利得を求める
  double getGain( const MapCategory& count, const MapTarget& targetCount,
                  const std::set< CATEGORY >& category );

  // カテゴリの組み合わせを求める
  static void getPair( const MapCategory& count, unsigned int flag, bool isLeft,
                       std::set< CATEGORY >* category );

  // 最大情報利得を持つカテゴリの組み合わせを求める
  std::pair< double, unsigned int >
  getMaxPair( const MapCategory& count, const MapTarget& targetCount );

  // 最大情報利得を持つ列の番号とカテゴリの組み合わせを求める
  std::pair< Index, unsigned int >
  getMaxVal( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
             const std::vector< Index >& rowIndex,
             const MapTarget& targetCount );

  // 分岐するノードを追加する
  void addNode( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
                size_type currentNode,
                const std::vector< Index >& rowIndex );

 public:

  // 学習パターン x と教師信号 y を指定して決定木を生成する
  Nominal( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y );

  // 学習パターンを指定して予測値を返す
  template< class Ran >
  std::pair< TARGET, bool > predict( Ran s ) const;
};

/*
  Nominal< CATEGORY, TARGET, Imp >::getGain
  : カテゴリの組み合わせに対する情報利得を求める

  count : 学習パターン・教師信号ごとのカウント値
  targetCount : 教師信号ごとのカウント値
  category : true 側のカテゴリ
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
double Nominal< CATEGORY, TARGET, Imp >::getGain
  ( const MapCategory& count, const MapTarget& targetCount,
    const std::set< CATEGORY >& category )
{
  // category に含まれる/含まれない学習パターンに対する教師信号ごとの和を求める
  MapTarget tSum, fSum;
  for ( typename MapCategory::const_iterator c = count.begin() ;
        c != count.end() ; ++c ) {
    MapTarget& sum = ( category.find( c->first ) != category.end() ) ?
      tSum : fSum;
    for ( typename MapTarget::const_iterator t = ( c->second ).begin() ;
          t != ( c->second ).end() ; ++t ) {
      typename MapTarget::iterator it = sum.find( t->first );
      if ( it == sum.end() )
        it = sum.insert( PairTarget( t->first, 0 ) ).first;
      it->second += t->second;
    }
  }

  std::pair< double, unsigned int > tGini = op_( tSum );
  std::pair< double, unsigned int > fGini = op_( fSum );
  std::pair< double, unsigned int > allGini = op_( targetCount );

  return( allGini.first -
          tGini.first * tGini.second / allGini.second -
          fGini.first * fGini.second / allGini.second );
}

/*
  Nominal< CATEGORY, TARGET, Imp >::getPair
  : カテゴリの組み合わせを求める

  count : 学習パターン・教師信号ごとのカウント値
  flag : 対象のカテゴリを 1 としたビット列
  isLeft : 左側のノードが対象なら true
  category : true 側のカテゴリを保持する配列
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
void Nominal< CATEGORY, TARGET, Imp >::getPair
  ( const MapCategory& count, unsigned int flag, bool isLeft,
    std::set< CATEGORY >* category )
{
  unsigned int bit = ( isLeft ) ? 1 : 0;
  category->clear();
  for ( typename MapCategory::const_iterator i = count.begin() ;
        i != count.end() ; ++i ) {
    if ( ( flag & 1 ) == bit )
      category->insert( i->first );
    flag >>= 1;
  }
}

/*
  Nominal< CATEGORY, TARGET, Imp >::getMaxpair
  : 最大情報利得を持つカテゴリの組み合わせを求める

  count : 学習パターン・教師信号ごとのカウント値
  targetCount : 教師信号ごとのカウント値
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
std::pair< double, unsigned int >
Nominal< CATEGORY, TARGET, Imp >::getMaxPair
  ( const MapCategory& count, const MapTarget& targetCount )
{
  typename MapCategory::size_type maxBit = std::numeric_limits< unsigned int >::digits;
  ErrLib::LE( count.size(), maxBit );

  if ( count.empty() ) return( std::pair< double, unsigned int >( 0, 0 ) );

  // 最初の列で最大情報利得を初期化
  unsigned int flag = 1; // 組み合わせを表すビット列
  std::set< CATEGORY > category;
  std::pair< double, unsigned int > max( 0, flag ); // 最大情報利得とその組み合わせ
  getPair( count, flag, true, &category );
  max.first = getGain( count, targetCount, category );

  // 最大情報利得を持つ組み合わせの探索
  unsigned int limit = 1 >> ( count.size() - 1 );
  while ( ++flag <= limit ) {
    getPair( count, flag, true, &category );
    double g = getGain( count, targetCount, category );
    if ( g > max.first ) {
      max.first = g;
      max.second = flag;
    }
  }

  return( max );
}

/*
  Nominal< CATEGORY, TARGET, Imp >::getMaxCol
  : 最大情報利得を持つ列の番号とカテゴリの組み合わせを求める

  x : 学習パターン
  y : 教師信号
  rowIndex : 読み込み対象の x, y に対する行番号
  targetCount : 教師信号ごとのカウント値を保持する配列
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
  std::pair< typename Nominal< CATEGORY, TARGET, Imp >::Index, unsigned int >
Nominal< CATEGORY, TARGET, Imp >::getMaxVal
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    const std::vector< Index >& rowIndex,
    const MapTarget& targetCount )
{
  // 学習パターン・教師信号ごとのカウント
  MapCategory count;

  // 最初の列で最大情報利得を初期化
  Index c = 0;
  Classify( x, y, rowIndex, c, &count );
  std::pair< double, unsigned int > maxVal = getMaxPair( count, targetCount ); // 最大情報利得とカテゴリの組み合わせ
  Index maxCol = c; // 最大情報利得を持つ列の番号

  // 最大情報利得を持つ列の探索
  while ( ++c < x.cols() ) {
    Classify( x, y, rowIndex, c, &count );
    std::pair< double, unsigned int > val = getMaxPair( count, targetCount );
    if ( val.first > maxVal.first ) {
      maxVal = val;
      maxCol = c;
    }
  }

  return( std::pair< Index, unsigned int >( maxCol, maxVal.second ) );
}

/*
  Nominal< CATEGORY, TARGET, Imp >::addNode
  : 分岐するノードを追加する

  x : 学習パターン
  y : 教師信号
  currentNode : 現在のノード
  rowIndex : 読み込み対象の x, y に対する行番号
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
void Nominal< CATEGORY, TARGET, Imp >::addNode
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    size_type currentNode,
    const std::vector< Index >& rowIndex )
{
  if ( rowIndex.empty() ) return;

  MapTarget targetCount; // 教師信号ごとのカウント
  Classify( y, rowIndex, &targetCount );
  std::pair< double, unsigned int > p = op_( targetCount );
  tree_[currentNode].impurity = p.first;
  tree_[currentNode].total = p.second;
  // 最大数を持つ教師信号を各ノードの target に登録する
  // 剪定を行ったときに葉ノードに変わる可能性があるため必要
  tree_[currentNode].target = op_.summary( targetCount );
  // リンク先はない状態で初期化しておく
  tree_[currentNode].left = tree_[currentNode].right = tree_.max_size();

  // 教師信号が単一なら分岐しないので終了
  if ( targetCount.size() < 2 ) return;

  // 情報利得が最大となる列番号とカテゴリの組み合わせを取得する
  std::pair< Index, unsigned int > maxVal = getMaxVal( x, y, rowIndex, targetCount );
  MapCategory count;
  Classify( x, y, rowIndex, maxVal.first, &count );

  // カテゴリが単一なら分岐できないので終了
  if ( count.size() < 2 ) return;

  // 左右の節のカテゴリを抽出する
  getPair( count, maxVal.second, true, &( tree_[currentNode].lCatg ) );
  getPair( count, maxVal.second, false, &( tree_[currentNode].rCatg ) );

  // 左右の節に行番号を分ける
  std::vector< Index > lRow, rRow;
  for ( typename std::vector< Index >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i ) {
    if ( tree_[currentNode].lCatg.find( x[*i][maxVal.first] ) != tree_[currentNode].lCatg.end() )
      lRow.push_back( *i );
    else
      rRow.push_back( *i );
  }

  tree_[currentNode].index = maxVal.first;

  // 分岐したノードを生成する
  tree_.push_back( BinTreeNode_() );
  tree_[currentNode].left = tree_.size() - 1;
  addNode( x, y, tree_.size() - 1, lRow );
  tree_.push_back( BinTreeNode_() );
  tree_[currentNode].right = tree_.size() - 1;
  addNode( x, y, tree_.size() - 1, rRow );
}

/*
  Nominal< CATEGORY, TARGET, Imp > コンストラクタ
  : 学習パターンと教師信号を指定して決定木を生成する

  x : 学習パターン
  y : 教師信号
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
Nominal< CATEGORY, TARGET, Imp >::Nominal
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y )
{
  ErrLib::CheckLinearModel( x, y.begin(), y.end() );

  tree_.push_back( BinTreeNode_() );

  // 読み込み対象の行列を連番で初期化
  std::vector< Index > rowIndex( y.size() );
  GenerateSerialNumber()( rowIndex.begin(), rowIndex.end() );

  addNode( x, y, 0, rowIndex );
}

/*
  Nominal< CATEGORY, TARGET, Imp >::predict : 学習パターンを指定して予測値を返す

  テンプレート引数の Ran はランダムアクセス反復子であることを想定している
  指定した学習パターンが決定木に存在しなかった場合は戻り値のペアの第二引数は false になる
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
template< class Ran >
std::pair< TARGET, bool >
Nominal< CATEGORY, TARGET, Imp >::predict( Ran s ) const
{
  size_type i = 0;
  while ( tree_[i].left != tree_.max_size() || tree_[i].right != tree_.max_size() ) {
    typename std::set< CATEGORY >::const_iterator it = tree_[i].lCatg.find( s[tree_[i].index] );
    if ( it != tree_[i].lCatg.end() ) {
      i = tree_[i].left;
      continue;
    }
    it = tree_[i].rCatg.find( s[tree_[i].index] );
    if ( it == tree_[i].rCatg.end() )
      return( std::pair< TARGET, bool >( TARGET(), false ) );
    i = tree_[i].right;
  }

  return( std::pair< TARGET, bool >( tree_[i].target, true ) );
}

情報利得を計算する仕組みは ID3 アルゴリズムの場合と同様です。ImpByGini クラスによってジニ不純度を計算に利用することができるようになります。したがって、ID3 アルゴリズムでジニ不純度を使ったり、逆に CART でエントロピーを利用するといったことも可能です。BinTreeNode は二分木用のノードで、子ノードへのリンクは左右二つのメンバ変数 left, right で表されています。また、学習パターンも左右の子ノード用に二つ ( lCatg, rCatg ) 用意されていますが、それぞれ複数個の学習パターンを取りうるので set クラスを使って連想配列の形で登録されるようになっています。

Nominal は CART を使って決定木を生成するためのクラスで、ID3 アルゴリズムと同じく決定木は構築時に生成されるようになっています。学習パターンの組み合わせの中で情報利得が最大となるものを求めるのはメンバ関数の getMaxPair で、考えられる組み合わせを getPair で順次求めながら情報利得を getGain で計算し、最大となる組み合わせを探索する流れになっています。
要素を 2 つに分割する組み合わせはビット列で表すことができます。例えば、要素が A, B, C の 3 種類ある場合に分割する組み合わせは、

ABC → A B | C
001

ABC → C A | B
010

ABC → A | B C
011

のように対応させることができます。ビット列をインクリメントするとこの後 100, 101, 110 と増えていきますが、これらはそれぞれ 011, 010, 001 をビット反転させたものになり組み合わせとしては等しくなるため除外できます ( 符号付き整数での正数と負数の関係と同じ意味になります )。したがって、ビット列の数は全てのパターンのちょうど半分を考えればよく、最後に 000 ( 111 ) は考慮しなくてもいいので、組み合わせの数は 23-1 - 1 = 3 と計算できます。一般に、N 個の要素を 2 つに分割する組み合わせは 2N-1 - 1 個となります。getMaxPair ではビット列を 1 から順にインクリメントしながら getPair を使って組み合わせを求め、getGain で組み合わせごとに再集計して不純度を計算し、情報利得を得ます。情報利得がそれまでの最大値を超えたらビット列を保持し、最終的に得られた最大情報利得とそのときのビット列のペアを戻り値として返します。
ビット列は unsigned int 型で保持しているため、要素数はそのビット数 ( 通常 32 ビット ) に制限されます。制限を超えているかのチェックは ErrLib::LE で行っていますが、この関数については実装内容は省略しています。実際に要素数が 32 個あった場合、仮に一つの組み合わせの処理に 1 msec かかるとすれば全ての組み合わせの処理を行うのに 25 日程度かかる計算になります。それだけたくさんの要素数に対して実際に処理することは無理なので、通常は後述する連続値での処理に切り替えることになります。

getMaxCol では情報利得が最大となる学習パターンの列番号を求めます。各列に対して最大の組み合わせは getMaxPair で得られるので、最終的には情報利得が最大となる学習パターンとその組み合わせの二つを求めることができます。得られた結果を元に行を分解し、子ノードを生成しながら再帰的に処理する部分は ID3 アルゴリズムの場合とほぼ同様です。しかし、今回は一度分割に使われた学習パターンが再度新たな組み合わせで最大情報利得をとる可能性があるため、列番号を外していく処理は行いません。

学習パターンが連続値である場合、組み合わせを使った分類とは異なる手法を用います。まずは学習パターンをソートしておき、ある値をしきい値としてそれより小さい要素と大きい要素に分類します。値の異なる要素の数が N 個あった場合、N - 1 通りの分類の仕方があります。教師信号も連続値ならば、決定木は回帰木であり、重回帰分析と同等の手法と考えることができます。このとき、教師信号の出現頻度はほとんどが 1 となり出現率を使った計算は意味がないため、重回帰分析同様に二乗和誤差を利用して分類前後の差を計算して最大となる分類の仕方を求めます ( 補足 1 )。

学習パターンと教師信号がともに連続値である場合のサンプル・プログラムを以下に示します。

/*
  ImpByVar : 不純度に二乗誤差平均(標本分散)を用いる場合の関数オブジェクト

  テンプレート引数の T は変数の型を表す
*/
template< class T >
struct ImpByVar
{
  // 二乗誤差の平均と総数を返す
  std::pair< double, unsigned int >
  operator()( const std::map< T, unsigned int >& count )
  {
    // 総数と平均を計算
    unsigned int total = 0;
    double mean = 0;
    for ( typename std::map< T, unsigned int >::const_iterator i = count.begin() ;
          i != count.end() ; ++i ) {
      mean += i->first * i->second;
      total += i->second;
    }
    if ( total == 0 )
      return( std::pair< double, unsigned int >( 0, 0 ) );
    mean /= total;

    // 二乗誤差平均の計算
    double var = 0;
    for ( typename std::map< T, unsigned int >::const_iterator i = count.begin() ;
          i != count.end() ; ++i ) {
      var += std::pow( i->first - mean, 2 ) * i->second;
    }
    var /= total;

    return( std::pair< double, unsigned int >( var, total ) );
  }

  // 要約値としてメディアン値を返す
  T summary( const std::map< T, unsigned int >& count )
  {
    typedef std::map< T, unsigned int > MapTarget;

    std::vector< T > data;

    for ( typename MapTarget::const_iterator i = count.begin() ;
          i != count.end() ; ++i ) {
      for ( unsigned int j = 0 ; j < i->second ; ++j )
        data.push_back( i->first );
    }

    typename std::vector< T >::size_type sz = data.size();
    if ( sz % 2 == 0 ) {
      if ( sz == 0 ) return( 0 );
      return( ( data[sz / 2 - 1] + data[sz / 2] ) / 2 );
    } else {
      return( data[sz / 2] );
    }
  }

};

/*
  Continuous : 連続値専用 CART アルゴリズム

  テンプレート引数の VALUE は学習パターン、TARGET は教師信号の型をそれぞれ表す
  テンプレート引数の Imp で不純度の種類を選択することができる
*/
template< class VALUE, class TARGET, template< class > class Imp = ImpByGini >
class Continuous
{
  // 型の定義
  typedef typename Matrix< VALUE >::size_type Index; // データ参照用インデックスの型

  typedef std::map< TARGET, unsigned int > MapTarget; // 教師信号に対するカウント値配列
  typedef std::map< VALUE, MapTarget > MapValue;      // 学習パターン・教師信号に対するカウント値配列

  typedef std::pair< TARGET, unsigned int > PairTarget; // 教師信号に対するカウント値要素
  typedef std::pair< VALUE, MapTarget > PairValue;      // 学習パターン・教師信号に対するカウント値要素

  typedef BinTreeNode< VALUE, TARGET > BinTreeNode_; // 二分木のノード
  typedef std::vector< BinTreeNode_ > BinTree;   // 二分木
  typedef typename BinTree::size_type size_type; // 二分木のインデックス

  // メンバ変数
  BinTree tree_;     // 決定木
  Imp< TARGET > op_; // 不純度を計算する関数オブジェクト

  // 学習パターンが VALUE であるデータを gtSum から leSum へ移動する
  static void mvValue( MapTarget* leSum, MapTarget* gtSum, const MapTarget& count );

  // VALUE以下の値と大きい値に学習パターンを分けたときの情報利得を求める
  double getGain( const MapTarget& leSum, const MapTarget& gtSum,
                  const MapTarget& targetCount );

  // 最大情報利得を持つ境界を求める
  std::pair< double, VALUE >
  getMaxBound( const MapValue& count, const MapTarget& targetCount );

  // 最大情報利得を持つ列の番号と境界を求める
  std::pair< Index, VALUE >
  getMaxVal( const Matrix< VALUE >& x, const std::vector< TARGET >& y,
             const std::vector< Index >& rowIndex, const MapTarget& targetCount );

  // 分岐するノードを追加する
  void addNode( const Matrix< VALUE >& x, const std::vector< TARGET >& y,
                size_type currentNode,
                const std::vector< Index >& rowIndex );

 public:

  // 学習パターン x と教師信号 y を指定して決定木を生成する
  Continuous( const Matrix< VALUE >& x, const std::vector< TARGET >& y );

  // 学習パターンを指定して予測値を返す
  template< class Ran >
  std::pair< TARGET, bool > predict( Ran s ) const;
};

/*
  Continuous< VALUE, TARGET, Imp >::mvValue
  : 学習パターンが VALUE であるデータを gtSum から leSum へ移動する

  leSum, gtSum : VALUE 以下/より大きい値に分けた結果を保持する配列
  count : 学習パターンが VALUE であるデータの配列
*/
template< class VALUE, class TARGET, template< class > class Imp >
void Continuous< VALUE, TARGET, Imp >::mvValue
  ( MapTarget* leSum, MapTarget* gtSum, const MapTarget& count )
{
  for ( typename MapTarget::const_iterator i = count.begin() ;
        i != count.end() ; ++i ) {
    typename MapTarget::iterator j = leSum->find( i->first );
    if ( j == leSum->end() )
      j = ( leSum->insert( PairTarget( i->first, 0 ) ) ).first;
    j->second += i->second;
    (*gtSum)[i->first] -= i->second;
  }
}

/*
  Continuous< VALUE, TARGET, Imp >::getGain
  : VALUE以下の値と大きい値に学習パターンを分けたときの情報利得を求める

  leSum, gtSum : VALUE 以下/より大きい値に分けた結果を保持する配列
  targetCount : 教師信号ごとのカウント値
*/
template< class VALUE, class TARGET, template< class > class Imp >
double Continuous< VALUE, TARGET, Imp >::getGain
  ( const MapTarget& leSum, const MapTarget& gtSum,
    const MapTarget& targetCount )
{
  std::pair< double, unsigned int > leImp = op_( leSum );
  std::pair< double, unsigned int > gtImp = op_( gtSum );
  std::pair< double, unsigned int > allImp = op_( targetCount );

  return( allImp.first -
          leImp.first * leImp.second / allImp.second - 
          gtImp.first * gtImp.second / allImp.second );
}

/*
  Continuous< VALUE, TARGET, Imp >::getMaxBound
  : 最大情報利得を持つ境界を求める

  count : 学習パターン・教師信号ごとのカウント値
  targetCount : 教師信号ごとのカウント値
*/
template< class VALUE, class TARGET, template< class > class Imp >
std::pair< double, VALUE >
Continuous< VALUE, TARGET, Imp >::getMaxBound
  ( const MapValue& count, const MapTarget& targetCount )
{
  MapTarget leSum;
  MapTarget gtSum( targetCount );
  double maxG = getGain( leSum, gtSum, targetCount ); // 最大情報利得
  VALUE maxValue = std::numeric_limits< VALUE >::min();

  // 最大情報利得を持つ境界の探索
  typename MapValue::const_iterator it = count.begin();
  for ( ; it != count.end() ; ++it ) {
    mvValue( &leSum, &gtSum, it->second );
    double g = getGain( leSum, gtSum, targetCount );
    if ( g > maxG ) {
      maxG = g;
      maxValue = it->first;
    }
  }
  // 境界は次の値との間をとるようにする
  typename MapValue::const_iterator next = count.find( maxValue );
  it = next++;
  if ( it != count.end() && next != count.end() ) {
    unsigned int sum0 = GetSum( it->second );
    unsigned int sum1 = GetSum( next->second );
    maxValue = ( maxValue * sum0 + next->first * sum1 ) / ( sum0 + sum1 );
  }

  return( std::pair< double, VALUE >( maxG, maxValue ) );
}

/*
  Continuous< VALUE, TARGET, Imp >::getMaxVal
  : 最大情報利得を持つ列の番号と境界を求める

  x : 学習パターン
  y : 教師信号
  rowIndex : 読み込み対象の x, y に対する行番号
  targetCount : 教師信号ごとのカウント値を保持する配列
*/
template< class VALUE, class TARGET, template< class > class Imp >
std::pair< typename Continuous< VALUE, TARGET, Imp >::Index, VALUE >
Continuous< VALUE, TARGET, Imp >::getMaxVal
  ( const Matrix< VALUE >& x, const std::vector< TARGET >& y,
    const std::vector< Index >& rowIndex, const MapTarget& targetCount )
{
  // 学習パターン・教師信号ごとのカウント
  MapValue count;

  // 第一列を最大情報利得を持つ列として初期化する
  Index c = 0;
  Classify( x, y, rowIndex, c, &count );
  std::pair< double, VALUE > maxVal = getMaxBound( count, targetCount ); // 最大情報利得と境界
  Index maxCol = c; // 最大情報利得を持つ列の番号

  // 最大情報利得を持つ列の探索
  while ( ++c < x.cols() ) {
    Classify( x, y, rowIndex, c, &count );
    std::pair< double, VALUE > val = getMaxBound( count, targetCount );
    if ( val.first > maxVal.first ) {
      maxVal = val;
      maxCol = c;
    }
  }

  return( std::pair< Index, VALUE >( maxCol, maxVal.second ) );
}

/*
  Continuous< VALUE, TARGET, Imp >::addNode
  : 分岐するノードを追加する

  x : 学習パターン
  y : 教師信号
  currentNode : 現在のノード
  rowIndex : 読み込み対象の x, y に対する行番号
*/
template< class VALUE, class TARGET, template< class > class Imp >
void Continuous< VALUE, TARGET, Imp >::addNode
  ( const Matrix< VALUE >& x, const std::vector< TARGET >& y,
    size_type currentNode, const std::vector< Index >& rowIndex )
{
  if ( rowIndex.empty() ) return;

  MapTarget targetCount; // 教師信号ごとのカウント
  Classify( y, rowIndex, &targetCount );
  std::pair< double, unsigned int > p = op_( targetCount );
  tree_[currentNode].impurity = p.first;
  tree_[currentNode].total = p.second;
  // 教師信号のメディアン値を各ノードの target に登録する
  // 剪定を行ったときに葉ノードに変わる可能性があるため必要
  tree_[currentNode].target = op_.summary( targetCount );
  // リンク先はない状態で初期化しておく
  tree_[currentNode].left = tree_[currentNode].right = tree_.max_size();

  // 教師信号が単一なら分岐しないので終了
  if ( targetCount.size() < 2 ) return;

  // 情報利得が最大となる列番号と境界を取得する
  std::pair< Index, VALUE > maxVal = getMaxVal( x, y, rowIndex, targetCount );

  // 左右の節に行番号を分ける
  std::vector< Index > lRow, rRow;
  for ( typename std::vector< Index >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i ) {
    if ( x[*i][maxVal.first] <= maxVal.second )
      lRow.push_back( *i );
    else
      rRow.push_back( *i );
  }

  // 左右の節のいずれかが空なら分岐できないので終了
  if ( lRow.empty() || rRow.empty() ) return;

  tree_[currentNode].index = maxVal.first;
  tree_[currentNode].lCatg.insert( maxVal.second );

  // 分岐したノードを生成する
  tree_.push_back( BinTreeNode_() );
  tree_[currentNode].left = tree_.size() - 1;
  addNode( x, y, tree_.size() - 1, lRow );
  tree_.push_back( BinTreeNode_() );
  tree_[currentNode].right = tree_.size() - 1;
  addNode( x, y, tree_.size() - 1, rRow );
}

/*
  Continuous< VALUE, TARGET, Imp > コンストラクタ

  x : 学習パターン
  y : 教師信号
*/
template< class VALUE, class TARGET, template< class > class Imp >
Continuous< VALUE, TARGET, Imp >::Continuous
  ( const Matrix< VALUE >& x, const std::vector< TARGET >& y )
{
  ErrLib::CheckLinearModel( x, y.begin(), y.end() );

  tree_.push_back( BinTreeNode_() );

  // 読み込み対象の行列を連番で初期化
  std::vector< Index > rowIndex( y.size() );
  GenerateSerialNumber()( rowIndex.begin(), rowIndex.end() );

  addNode( x, y, 0, rowIndex );
}

/*
  Continuous< VALUE, TARGET, Imp >::predict : 学習パターンを指定して予測値を返す

  テンプレート引数の Ran はランダムアクセス反復子であることを想定している
  指定した学習パターンが決定木に存在しなかった場合は戻り値のペアの第二引数は false になる
*/
template< class CATEGORY, class TARGET, template< class > class Imp >
template< class Ran >
std::pair< TARGET, bool >
Continuous< CATEGORY, TARGET, Imp >::predict( Ran s ) const
{
  size_type i = 0;
  while ( tree_[i].left != tree_.max_size() || tree_[i].right != tree_.max_size() ) {
    CATEGORY bound = *( tree_[i].lCatg.begin() );
    i = ( s[tree_[i].index] <= bound ) ? tree_[i].left : tree_[i].right;
  }

  return( std::pair< TARGET, bool >( tree_[i].target, true ) );
}

ImpByVar は、教師信号をキーとした要素数を使って二乗誤差平均 ( 標本分散 ) を計算するための関数オブジェクトです。連続値の場合、教師信号ごとに要素数を集計するのはかえってデータ量を増やすことになるのですが、先に紹介した不純度計算用のクラス ImpByEntropy や ImpByGini と引数を合わせるためにこのような形にしています。名義尺度のデータの場合、葉ノードに複数の教師信号が混在した場合の代表値は要素数の最も多いもの ( 最頻値 ) を選択していました。連続値の場合は最頻値を代表値にするのはあまり意味がないので、平均値やメディアン値などに切り替える必要があります。ImpByVar では代表値を返すメンバ関数 summary にてメディアン値を算出するようにしてあります。

Continuous クラスは、学習パターンが連続値である場合の CART アルゴリズム本体になります。map はデータがソートされた状態で格納されるので並べ替えを行う必要がありません。したがって、配列の開始位置から順番にデータを区切って集計し、情報利得が最大となる区切り方を求めればよいことになります。その処理を行うのがメンバ関数 getMaxBound です。getMaxBound ではまず最初に教師信号ごとの要素数を格納するための map 型変数 leSum と gtSum を用意し、全ての ( 学習パターンによらない ) 要素数 targetCount をいったん gtSum 側にコピーします ( leSum は空の状態になります )。あとは、学習パターン・教師信号ごとに集計された要素数 count の最初の要素から順に集計値を gtSum から leSum に移す処理を行い、情報利得を求めていきます。集計値を移す処理はメンバ関数 mvValue が行っています。残りの処理は名義尺度の場合とほとんど同一です。

下図は、ID3 アルゴリズムのサンプル・プログラムで処理したときと同じデータを使い、Nominal クラスで生成した決定木です。各ノードには、そのときのジニ不純度の値と、葉ノードには選択された教師信号が記されています。

図 3-1. CART アルゴリズムによる決定木の生成 ( 名義尺度 )
CARTアルゴリズムによる決定木

上図を見ると、女性に対する教師信号は 10 代、20 代、30 代全てにおいて "いちご" が選択されています。女性については 10 人中 9 人がいちごを選んでいるので、年代別に分けなければもう少しシンプルな木構造になるのですが、機械的に処理するとこのような結果になります。この状態を「過学習 ( Overfitting )」 といいます。過学習は、このように無駄な分類をしてしまう他に、異常値に対して過度に反応してしまうような状態もあります。CART では、いったん過学習状態になるまで決定木を生成してから後で不要な部分木を取り除く「事後枝刈り ( Post-pruning )」と呼ばれる手法を用いて過学習に対応しています。

枝刈りには以下の評価関数を用います。

C(T) = Σk{1→|T|}( NkGk ) + α|T|

上式において T は葉ノードを表し、|T| は葉ノードの数、Nk、Gk はそれぞれ k 番目の葉ノードにあるデータの要素数と不純度です。また、α は木構造の複雑さを制御するためのパラメータで、α が大きいほど枝刈り後の決定木は単純になります。

枝刈りをすることによって決定木の評価関数の値は変化します。その中から最小となる枝刈りの仕方を探索して枝刈りを行うことを繰り返します。枝刈り前の決定木の評価関数の値が最小になったら処理を停止します。

枝刈りを行うサンプル・プログラムを以下に示します。

/*
  GetScore : 剪定のための評価関数の計算

  tree : 対象の決定木
  i : 計算対象のノード番号 ( 計算対象をルートとする部分木の評価関数を計算する )
  leaf : 葉ノードとみなすノードの番号
  alpha : 木の複雑さを制御するパラメータ
*/
template< class X, class Y >
double GetScore( const std::vector< BinTreeNode< X, Y > >& tree,
                 typename std::vector< BinTreeNode< X, Y > >::size_type i,
                 typename std::vector< BinTreeNode< X, Y > >::size_type leaf,
                 double alpha )
{
  if ( i == leaf )
    return( tree[i].impurity * tree[i].total + alpha );
  if ( tree[i].left == tree.max_size() && tree[i].right == tree.max_size() )
    return( tree[i].impurity * tree[i].total + alpha );

  double score = 0;
  if ( tree[i].left != tree.max_size() )
    score += GetScore( tree, tree[i].left, leaf, alpha );
  if ( tree[i].right != tree.max_size() )
    score += GetScore( tree, tree[i].right, leaf, alpha );

  return( score );
}

/*
  PruningFromRoot : ルート側から対象を探して剪定する

  tree : 対象の決定木へのポインタ
  alpha : 木の複雑さを制御するパラメータ
*/
template< class X, class Y >
void PruningFromRoot( std::vector< BinTreeNode< X, Y > >* tree, double alpha )
{
  typedef typename std::vector< BinTreeNode< X, Y > >::size_type size_type;

  // 現状の評価関数の値を取得
  double minScore = GetScore( *tree, 0, tree->size(), alpha );
  // 途中のノードを葉ノードとみなしたときの評価関数を計算して最小値を探索
  for ( ; ; ) {
    size_type minIndex = tree->size();
    for ( size_type i = 0 ; i < tree->size() ; ++i ) {
      if ( (*tree)[i].left == tree->max_size() && (*tree)[i].right == tree->max_size() )
        continue;
      double score = GetScore( *tree, 0, i, alpha );
      if ( score < minScore ) {
        minIndex = i;
        minScore = score;
      }
    }

    if ( minIndex >= tree->size() ) break;
    (*tree)[minIndex].left = (*tree)[minIndex].right = tree->max_size();
  }
}

評価関数の計算に必要な要素数と不純度は、CART アルゴリズムの中であらかじめメンバ変数 impurity, total として保持してあります。したがって、評価関数の計算はこれらの値を使って集計を行うだけで済みます。集計は再帰的に行い、葉ノードならば計算結果を返し、途中のノードならば左右の子ノードに対して自分自身を呼び出すようにしています。あとは、評価関数の値の中で最小となる枝刈りの仕方を探索し、対象となるノードを葉ノードに変換するだけです。葉ノードが持つべきデータである教師信号の代表値も CART アルゴリズムの中ですでに保持させているので、葉ノードの変換は子ノードへのリンクを削除するだけで終わります。

関数 PruningFromRoot は全てのノードに対して評価関数を計算するため、決定木が大きくなると処理は重くなっていくことが考えられます。全てのノードに対して再帰的に計算することをやめて、葉ノードを持つノードだけに限定して枝刈りを行うようにすることで処理時間をある程度短縮させることができます。その代償として、得られた枝刈りの結果が最適であるという保証はなくなります。

左右両方の子が葉ノードであるノードに対し、その子ノードを捨てて自身を葉ノードにしたときに評価関数が減少する分は、左右の子ノードの要素数 x 不純度の和と自ノードの要素数 x 不純度との差を求め、さらに α を加えることで得ることができます。この減少分が最大となるノードを探索して枝刈りすることを、最大の減少分が負数となるまで繰り返します。以下にサンプル・プログラムを示します。

/*
  HasLeaves : 対象が葉ノードを持つ場合 true を返す

  tree : node を持つ決定木
  node : チェック対象のノード
*/
template< class X, class Y >
bool HasLeaves( const std::vector< BinTreeNode< X, Y > >& tree,
                const BinTreeNode< X, Y >& node )
{
  typedef BinTreeNode< X, Y > BinTreeNode;

  if ( node.left == tree.max_size() || node.right == tree.max_size() )
    return( false );

  if ( node.left != tree.max_size() ) {
    const BinTreeNode& left = tree[node.left];
    if ( left.left != tree.max_size() || left.right != tree.max_size() )
      return( false );
  }
  if ( node.right != tree.max_size() ) {
    const BinTreeNode& right = tree[node.left];
    if ( right.left != tree.max_size() || right.right != tree.max_size() )
      return( false );
  }

  return( true );
}

/*
  PruningFromLeaf : 葉ノードを持つノードから剪定を行う

  tree : 剪定を行う決定木へのポインタ
  alpha : 木の複雑さを制御するパラメータ
*/
template< class X, class Y >
void PruningFromLeaf( std::vector< BinTreeNode< X, Y > >* tree, double alpha )
{
  typedef std::vector< BinTreeNode< X, Y > > BinTree;

  for ( ; ; ) {
    double maxDelta = -1;
    typename BinTree::iterator maxIt = tree->end();
    for ( typename BinTree::iterator i = tree->begin() ;
          i != tree->end() ; ++i ) {
      if ( ! HasLeaves( *tree, *i ) ) continue;
      // 葉ノードの不純度 x データ数 - 自ノードの不純度 x データ数が
      // 剪定をしたときの減少分
      double delta = (*tree)[i->left].impurity * (*tree)[i->left].total +
        (*tree)[i->right].impurity * (*tree)[i->right].total -
        i->impurity * i->total + alpha;
      if ( delta < 0 ) continue;
      if ( delta > maxDelta ) {
        maxDelta = delta;
        maxIt = i;
      }
    }
    if ( maxDelta < 0 ) break;
    maxIt->left = maxIt->right = tree->max_size();
  }
}

下図は、PruningFromRoot を使い、α = 1 で枝刈りを行った後の決定木です。女性に対する教師信号が全て "いちご" に集約されていることがわかります。ちなみに α = 1 で PruningFromLeaf を使った場合、性別のみで分類された決定木となるため、少なくともサンプルのデータに対しては PruningFromLeaf の方が α の影響を強く受ける傾向にあるようです。α = 0.5 にすると同じ結果が得られるようになりました。

図 3-2. 枝刈り後の決定木
枝刈り後の決定木

下図は「アヤメの品種データ ( iris flower data set )」を使って回帰分析を行った結果です。グラフの横軸は「花びらの長さ ( petal length )」、縦軸は「花びらの幅 ( petal width )」となります。花びらの長さを学習パターン、花びらの幅を教師信号として Continuous と ImpByVar の組み合わせで決定木を生成し、α = 0.1 で枝刈りを行った結果から回帰曲線を描画しています。

図 3-3. 決定木による回帰分析
決定木による回帰分析

4) CHAID ( CHi-square Automatic Interaction Detection )

「CHAID ( CHi-square Automatic Interaction Detection )」は「ゴードン・カス ( Gordon Kass )」によって 1980 年に発表された手法です。今まで紹介した手法と異なる点としては、学習パターンの組み合わせから複数個 ( CART のように二つに限定しない ) にノードを分割することができることと、分割方法の評価にカイ二乗値を利用していることが挙げられます。

学習パターンが ( A1, A2, ... Ak ) の k 種類、教師信号が ( B1, B2, ... Bl ) の l 種類あったとき、各要素は ( Ai, Bj ) ( i = 1, 2, ... k ; j = 1, 2, ... l ) の組み合わせごとに集計することができるので、その集計値を rij とします。

N = Σi{1→k}( Σj{1→l}( rij ) )
ai = Σj{1→l}( rij )
bj = Σi{1→k}( rij )

としたとき、

χ2 = Σi{1→k}( Σj{1→l}( ( rij - aibj / N )2 / ( aibj / N ) ) )

は、Ai と Bj がすべて互いに独立ならば自由度 ( k - 1 )( l - 1 ) の χ2-分布に従います。各変数の関係は以下の表のような形になります。

表 4-1. 分割表
A1A2...Ak計
B1r11r21...rk1b1
B2r12r22...rk2b2
:::::
Blr1lr2l...rklbl
計a1a2...akN

この性質を利用した検定法を「(ピアソンの)χ2-検定 ( (Pearson's) Chi-square Test )」といいます。χ2 を使って p 値を求めることができるので、この値が小さければ Ai と Bj は互いに独立ではなく何らかの関連性があると判断することができます。CHAID では、この p 値が大きいもの、すなわち rij が学習パターンまたは教師信号ごとに偏りのないものどうしを一つにまとめ、分類の組み合わせを決定します。

N 個の要素を分割する方法の総数は「ベル数 ( Bell Number )」として知られています ( 補足 2 )。この数は、N が増大するに従い爆発的に増えるので、CHAID では全ての可能性を調べる代わりに次のような方法を用います。

ステップ 1 : 二つの学習パターンの組み合わせに対して χ2 を求め、その中から最小値をとるもの ( すなわち p 値が最大となるもの ) を探索する。p 値がしきい値以上なら二つの学習パターンをマージして単一の複合パターンと見なし、この手順を繰り返す。

ステップ 2 : 三つ以上の学習パターンからなる各複合パターンについて、二分割したときに χ2 が最大となる ( すなわち p 値が最小となる ) 方法を探索し、p 値がしきい値より小さいならば分割を実行する。全てを走査したらステップ 1 に戻る。

最初の節で示した例について、年代別に作成した分割表は以下のようになります。

表 4-2. 年代別分割表
10代20代30代計
いちご28212
りんご2114
バナナ1034
計59620

年代は 10 代、20 代、30 代の 3 つあるので、{ 10 代, 20 代 }、{ 20 代, 30 代 }、{ 30 代, 10 代 } のそれぞれのペアで χ2 を求めます。その結果は以下のようになります。

表 4-3. χ2 ( 対角より右上 ) とその p 値 ( 対角より左下 )
10代20代30代
10代4.131.25
20代0.1276.25
30代0.5350.0439

表の中で、対角より右上に示したのが χ2、左下側がその p 値です。教師信号はいちご・りんご・バナナの 3 つなので、χ2 は自由度 ( 2 - 1 ) x ( 3 - 1 ) = 2 の χ2-分布に従うことになり p 値を求めることができます。3 つのペアの中で { 30 代, 10 代 } のペアの χ2 が最も小さく、したがって p 値が最も大きいので、p 値がしきい値以上であればこの二つをマージすることになります。マージ後は 30 代 + 10 代と 20 代の二つをペアとして同様の処理を行います。χ2 = 6.25 ( p = 0.0452 ) となるので、もししきい値が 0.05 ならば p 値はそれより小さくなりマージは行いませんが、しきい値が 0.01 ならばマージされることになります。いずれの場合もこれ以上マージはできないので処理は完了となります。ステップ 2 の分割処理は、もし 10 代、20 代、30 代の全ての年代が一つにマージされていれば、二分割できる全てのパターン ( { 10 代 + 20 代, 30 代 }、{ 20 代 + 30 代, 10 代 }、{ 30 代 + 10 代, 20 代 } ) に対して χ2 と p 値を求めます。このとき、{ 30 代 + 10 代, 20 代 } についてはマージのときに計算済なので省略することができます。残りの 2 パターンについて計算すると、{ 10 代 + 20 代, 30 代 } で χ2 = 4.92 ( p = 0.0854 )、{ 20 代 + 30 代, 10 代 } で χ2 = 1.78 ( p = 0.411 ) となるので、p 値が最小となるのは { 30 代 + 10 代, 20 代 } の場合でしきい値より p 値が小さいなら分割します。もし分割したならば 3 つ以上にマージされたものは存在しなくなり、分割されなければこれ以上分割されることはないので、いずれの場合もステップ 2 は終了します。

分類の組み合わせは各学習パターンに対して得られるので、得られた組み合わせで再度 χ2 とその p 値を求め、その中から p 値の最も小さいものを使って決定木を生成すればよいことになりますが、CHAID では p 値にさらに補正を掛けるようになっています。
検定を行うときの危険率が α だったとき、検定結果が正しい確率は 1 - α となります。したがってこの検定が互いに独立で n 回繰り返されたとき、正しい確率は ( 1 - α )n です。このとき、検定結果が誤りである確率は

1 - ( 1 - α )n ≅ 1 - ( 1 - nα ) = nα

となるので、全体での危険率を α にするためには一回の検定での危険率を α / n にすればよいことになります。これを「ボンフェローニ補正 ( Bonferroni Correction )」といいます。p 値に対しては、逆に検定回数を掛けることになります。この検定回数としては、学習パターンに含まれる種別を分割したときの分割の仕方を使います。これは「第二種スターリング数 ( Stirling Numbers of the Second Kind )」と呼ばれる数になります ( 補足 2 )。例えば、学習パターンに 10 代、20 代、30 代の 3 つのカテゴリがあって、それを 2 つに分割する場合、第二種スターリング数は 3 になるので、求めた p 値に 3 を掛けることになります。

CHAID では全ての組み合わせにおいて p 値がしきい値以上ならば分割はされないので、全学習パターンについてそのような状態になれば処理を停止することになります。また、この他にも CHAID は以下の停止条件を持っています。

  1. 最大レベル数 ( 木の高さ )
  2. 最小ノードサイズ
  3. 最小分割サイズ

最大レベル数は木の高さの最大値で、設定された最大レベル数に達したノードはそれ以上分割されません。また、最小ノードサイズはノードに含まれる要素数の最小値を表し、指定した最小ノード数より少ない要素数を持ったノードはそれ以上分割しません。最小分割サイズは、分割時に各ノードに割り当てられる要素数の最小値を表し、分割後に指定した最小分割サイズに達しないノードが生じる場合は分割するのをやめます。

CHAID のサンプル・プログラムを以下に示します。まずは χ2 の計算処理部分です。

/*
  MapToMatrix : map 形式のデータ count から行列 res への変換

  yCount : Y のみをキーとしたカウント値
*/
template< class X, class Y >
void MapToMatrix( const std::map< X, std::map< Y, unsigned int > >& count,
                  const std::map< Y, unsigned int >& yCount,
                  Matrix< unsigned int >* res )
{
  typedef std::map< Y, unsigned int > MapY;
  typedef std::map< X, MapY > MapXY;

  res->assign( yCount.size(), count.size() );
  Matrix< unsigned int >::size_type x = 0;
  for ( typename MapXY::const_iterator ix = count.begin() ;
        ix != count.end() ; ++ix ) {
    for ( typename MapY::const_iterator iy = ( ix->second ).begin() ;
          iy != ( ix->second ).end() ; ++iy ) {
      typename MapY::const_iterator i = yCount.find( iy->first );
      if ( i == yCount.end() ) continue; // ありえないが念のため
      Matrix< unsigned int >::size_type y = std::distance( yCount.begin(), i );
      (*res)[y][x] = iy->second;
    }
    ++x;
  }
}

/*
  GetChiSquareP : 分割表 count と対象列 flag を元にカイ二乗値の p 値を計算する

  flag はビット列とし、立っているビットの組み合わせでカウントを集計した上で計算する
*/
double GetChiSquareP( const Matrix< unsigned int >& count, const std::vector< unsigned int >& flag )
{
  Matrix< double > obs( count.rows(), flag.size(), 0 ); // 再集計した結果
  std::vector< double > xSum( count.rows(), 0 ); // X 方向の合計
  std::vector< double > ySum( flag.size(), 0 );  // Y 方向の合計
  double sum = 0; // 総計

  // flag にあるビット列を元にカウントを再集計
  for ( Matrix< unsigned int >::size_type ox = 0 ;
        ox < flag.size() ; ++ox ) {
    Matrix< unsigned int >::size_type x = 0;
    for ( unsigned int f = flag[ox] ; f != 0 ; f >>= 1 ) {
      if ( ( f & 1 ) != 0 ) {
        for ( Matrix< unsigned int >::size_type y = 0 ;
              y < count.rows() ; ++y ) {
          obs[y][ox] += count[y][x];
          xSum[y] += count[y][x];
          ySum[ox] += count[y][x];
          sum += count[y][x];
        }
      }
      ++x;
    }
  }

  // カイ二乗値の計算
  double chiSq = 0;
  for ( Matrix< unsigned int >::size_type y = 0 ; y < obs.rows() ; ++y ) {
    for ( Matrix< unsigned int >::size_type x = 0 ; x < obs.cols() ; ++x ) {
      double exp = xSum[y] * ySum[x] / sum;
      if ( exp == 0 ) continue;
      chiSq += std::pow( obs[y][x] - exp, 2 ) / exp;
    }
  }

  ChiSquareDistribution chiDist( ( obs.rows() - 1 ) * ( obs.cols() - 1 ) );

  return( 1 - chiDist.lower_p( chiSq ) );
}

χ2 の計算時は、各要素数を保持するのに学習パターンを列、教師信号を行とする行列 Matrix を利用します。今までのアルゴリズムと同様に、要素数は最初 map で持たせる仕様なので、map 形式から Matrix 形式に型変換するための関数として MapToMatrix を用意してあります。行列への変換には学習パターン・教師信号ごとの要素数だけあれば可能なので、引数としては count のみでもいいのですが、教師信号ごとの要素数が全学習パターンに対してそろっているとは限らないので ( 要素数がゼロの場合、要素ごと存在しないことになるので )、全教師信号を把握するために、教師信号ごとに集計した配列 yCount も引数として渡してあります。そのため yCount はキーとなる教師信号のみを利用し、要素数そのものは使っていません。
なお、Matrix クラスのメンバ関数 assign は、指定した行・列数で行列を初期化するためのものです。

MapToMatrix で生成した行列を使い、GetChiSquareP で χ2 を計算します。計算方法は少し複雑で、引数 flag で渡されたビット列が対象列を表しています。例えば、flag = { 0001, 1010, 0100 } の場合、行列の 1 列目、2 列目と 4 列目の合計、3 列目を使って計算を行うことになります。χ2 を計算したら、この値が χ2-分布に従うとみなして p 値を求め、その値を返します。p 値の計算には ChiSquareDistribution クラスを利用しています。このクラスの詳細については「確率・統計 (6) 標本分布」の「χ2-分布 ( Chi-square Distribution )」にあるサンプル・プログラムをご覧ください。

次はステップ 1 のマージ処理部分です。

/*
  GetCache : キャッシュからのデータ取得
*/
bool GetCache( const std::map< unsigned int, std::map< unsigned int, double > >& cache,
               unsigned int i1, unsigned int i2, double* p )
{
  typedef std::map< unsigned int, double > Map;
  typedef std::map< unsigned int, Map > TwoDimMap;

  unsigned int key1 = std::min( i1, i2 );
  unsigned int key2 = std::max( i1, i2 );

  TwoDimMap::const_iterator i = cache.find( key1 );
  if ( i == cache.end() ) return( false );
  Map::const_iterator j = ( i->second ).find( key2 );
  if ( j == ( i->second ).end() ) return( false );

  *p = j->second;

  return( true );
}

/*
  SetCache : キャッシュへのデータセット
*/
void SetCache( std::map< unsigned int, std::map< unsigned int, double > >* cache,
               unsigned int i1, unsigned int i2, double p )
{
  unsigned int key1 = std::min( i1, i2 );
  unsigned int key2 = std::max( i1, i2 );

  (*cache)[key1][key2] = p;
}

/*
  GetMaxP : 分割表 count から対象列 flag に対するカイ二乗値の p 値が最大のものを求める

  f1, f2 : flag の中で p 値が最大となる組み合わせを取得する変数へのポインタ
*/
double GetMaxP( const Matrix< unsigned int >& count,
                const std::vector< unsigned int >& flag,
                unsigned int* f1, unsigned int* f2,
                std::map< unsigned int, std::map< unsigned int, double > >* cache )
{
  double maxP = -1;
  double p;
  std::vector< unsigned int > pair( 2 );
  for ( std::vector< unsigned int >::const_iterator i1 = flag.begin() ;
        i1 != flag.end() ; ++i1 ) {
    pair[0] = *i1;
    for ( std::vector< unsigned int >::const_iterator i2 = i1 + 1 ;
          i2 != flag.end() ; ++i2 ) {
      if ( ! GetCache( *cache, *i1, *i2, &p ) ) {
        pair[1] = *i2;
        p = GetChiSquareP( count, pair );
        SetCache( cache, *i1, *i2, p );
      }
      if ( maxP < p ) {
        *f1 = *i1;
        *f2 = *i2;
        maxP = p;
      }
    }
  }

  return( maxP );
}

/*
  MergeCategory : 分割表 count から対象列 flag に対するカイ二乗値の p 値が最大の組み合わせをマージしていく

  threshold : マージ処理を停止する p 値のしきい値(最大値がthresholdより小さくなったら停止)
*/
bool MergeCategory( const Matrix< unsigned int >& count,
                    std::vector< unsigned int >* flag, double threshold,
                    std::map< unsigned int, std::map< unsigned int, double > >* cache )
{
  unsigned int f1, f2;
  std::vector< unsigned int >::size_type preSz = flag->size();

  while ( flag->size() > 1 ) {
    double maxP = GetMaxP( count, *flag, &f1, &f2, cache );
    if ( maxP < threshold )
      break;

    flag->erase( std::find( flag->begin(), flag->end(), f1 ) );
    flag->erase( std::find( flag->begin(), flag->end(), f2 ) );
    flag->push_back( f1 + f2 );
  }

  return( flag->size() < preSz );
}

メインとなる関数は GetMaxP です。引数の count は前述の学習パターン・教師信号ごとの要素数で、flag はマージされた列の番号をビット列で表したものです。例えば flag = { 001, 110 } ならば、1 行目は単独で、2, 3 行目がマージされていることを表しています。処理内容は、flag の中から二つのビット列を選択しながら GetChiSquareP を呼び出して p 値を計算し、その中で最大のものを抽出するというものです。p 値が最大となるビット列は f1 と f2 に保持され、p 値の最大値は戻り値として返されます。
計算前に GetCache でキャッシュに p 値が保持されていないかを確認し、あれば計算をスキップするようにしています。もしキャッシュになかった場合は計算後に SetCache を使ってキャッシュに登録します。

関数 MergeCategory は GetMaxP を呼び出して p 値が最大となるビット列を抽出し、p 値がしきい値以上なら二つのビット列をマージする処理を繰り返すだけの内容となっています。もし、p 値がしきい値より小さければ処理は終了となります。また、これ以上マージできない ( ビット列が一つのみになった ) ときも処理停止となっています。

ステップ 2 の分割処理部分を以下に示します。

/*
  SplitFlag : ビット列 flag を mask のビットに従って f1 と f2 に分ける
*/
void SplitFlag( unsigned int flag, unsigned int mask, unsigned int* f1, unsigned int* f2 )
{
  *f1 = *f2 = 0;
  unsigned int bit = 1;
  while ( flag != 0 ) {
    if ( ( flag & 1 ) != 0 ) {
      if ( ( mask & 1 ) != 0 )
        *f1 |= bit;
      else
        *f2 |= bit;
      mask >>= 1;
    }
    flag >>= 1;
    bit <<= 1;
  }
}

/*
  BitCount : flag の立っているビット数を数える
*/
unsigned int BitCount( unsigned int flag )
{
  unsigned int count = 0;
  for ( ; flag != 0 ; flag >>= 1 ) {
    if ( ( flag & 1 ) != 0 )
      ++count;
  }

  return( count );
}

/*
  GetMinP : 分割表 count からビット列 flag を分割したときにカイ二乗値の p 値が最小になるものを求める

  f1, f2 : flag を分割したときに p 値が最小となるものを取得する変数へのポインタ
*/
double GetMinP( const Matrix< unsigned int >& count,
                unsigned int flag, unsigned int* f1, unsigned int* f2,
                std::map< unsigned int, std::map< unsigned int, double > >* cache )
{
  unsigned int bitCount = BitCount( flag );
  if ( bitCount < 3 ) return( 1 );

  double minP = 2;
  double p;
  std::vector< unsigned int > pair( 2 );
  unsigned int limit = std::pow( 2, bitCount - 1 );
  for ( unsigned int mask = 1 ; mask < limit ; ++mask ) {
    SplitFlag( flag, mask, &pair[0], &pair[1] );
    if ( ! GetCache( *cache, pair[0], pair[1], &p ) ) {
      p = GetChiSquareP( count, pair );
      SetCache( cache, pair[0], pair[1], p );
    }
    if ( p < minP ) {
      *f1 = pair[0];
      *f2 = pair[1];
      minP = p;
    }
  }

  return( minP );
}

/*
  SplitCategory : 分割表 count からカイ二乗値の p 値が最小になるように flag を分割する

  threshold : 分割を判断する p 値のしきい値(最大値がthresholdより小さい場合は分割)
*/
bool SplitCategory( const Matrix< unsigned int >& count,
                    std::vector< unsigned int >* flag, double threshold,
                    std::map< unsigned int, std::map< unsigned int, double > >* cache )
{
  unsigned int f1, f2;
  std::vector< unsigned int >::size_type preSz = flag->size();

  for ( std::vector< unsigned int >::size_type i = 0 ; i < flag->size() ; ) {
    double maxP = GetMinP( count, (*flag)[i], &f1, &f2, cache );
    if ( maxP >= threshold ) {
      ++i;
      continue;
    }

    flag->erase( flag->begin() + i );
    flag->push_back( f1 );
    flag->push_back( f2 );
  }

  return( flag->size() > preSz );
}

分割処理のメインは GetMinP で、前述の GetMaxP と処理の流れはかなりよく似ています。渡されたビット列 flag を二つのビット列に分割する処理は SplitFlag が行い、何番目のビットを f1 側に割り当てるかを mask が保持しています。mask の使い方が少し特殊で、flag の中で立っているビットだけに対して何番目を f1 に登録するか判断するような形になります。例えば flag が 1001 の場合、実際に立っているビットは 2 つです。mask が 10 ならば、flag の 1 番目のビットは mask の 1 番目のビット 0 に対応するので f2 側に割り当てられ、flag の 4 番目のビットは mask の 2 番目のビット 1 に対応するので f1 側に割り当てられます。結果として f1 = 1000, f2 = 0001 となります。

関数 SplitCategory は GetMinP を使って p 値が最小になる分割の仕方を求め、p 値がしきい値より小さければ実際にビット列の分割を行います。分割されたビット列は flag の末尾に追加されるので、さらに分割処理の対象となることに注意して下さい。

以下は、残りの部分を示したものになります。

/*
  StirlingNumber_2nd : 第二種スターリング数を求める

  第二種スターリング数は、c 個のカテゴリを g 個のグループに分類する場合の数に相当する
*/
double StirlingNumber_2nd( double c, double g )
{
  double ans = 0;
  for ( unsigned int i = 0 ; i < g ; ++i )
    ans += std::pow( -1, i ) * std::pow( g - i, c ) /
      ( std::tgamma( i + 1 ) * std::tgamma( g - i + 1 ) );

  return( ans );
}

/*
  OptimizeCategory : 最適なカテゴリ分割を探索する

  threshold : マージ・分割処理を停止する p 値のしきい値
*/
template< class CATEGORY, class TARGET >
double OptimizeCategory( const std::map< CATEGORY, std::map< TARGET, unsigned int > >& count,
                         const std::map< TARGET, unsigned int >& targetCount,
                         std::vector< unsigned int >* flag, double threshold )
{
  Matrix< unsigned int > mat;
  MapToMatrix( count, targetCount, &mat );

  std::map< unsigned int, std::map< unsigned int, double > > cache;

  flag->clear();
  for ( unsigned int i = 0 ; i < mat.cols() ; ++i )
    flag->push_back( 1 << i );

  for ( ; ; ) {
    if ( ! MergeCategory( mat, flag, threshold, &cache ) ) break;
    if ( ! SplitCategory( mat, flag, threshold, &cache ) ) break;
  }

  return( ( flag->size() > 1 ) ? GetChiSquareP( mat, *flag ) * StirlingNumber_2nd( mat.cols(), flag->size() ) : 1 );
}

/*
  GetMinPCategory : 最小 p 値を持つ列の番号とカテゴリ分割を求める

  x : 学習パターン
  y : 教師信号
  rowIndex : 読み込み対象の x, y に対する行番号
  targetCount : 教師信号ごとのカウント値を保持する配列
*/
template< class CATEGORY, class TARGET >
typename Matrix< CATEGORY >::size_type
GetMinPCategory( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
                 const std::vector< typename Matrix< CATEGORY >::size_type >& rowIndex,
                 const std::map< TARGET, unsigned int >& targetCount,
                 std::vector< unsigned int >* minFlag, double threshold )
{
  // 学習パターン・教師信号ごとのカウント
  std::map< CATEGORY, std::map< TARGET, unsigned int > > count;
  // 最適なカテゴリ分割
  std::vector< unsigned int > flag;

  double minP = 2; // 最小 p 値
  typename Matrix< CATEGORY >::size_type minCol = 0; // 最小 p 値を持つ列の番号

  // 最小 p 値を持つ列の探索
  for ( typename Matrix< CATEGORY >::size_type c = 0 ; c < x.cols() ; ++c ) {
    Classify( x, y, rowIndex, c, &count );
    double p = OptimizeCategory( count, targetCount, &flag, threshold );
    if ( p < minP ) {
      minP = p;
      minCol = c;
      minFlag->assign( flag.begin(), flag.end() );
    }
  }

  return( minCol );
}

/*
  FlagToMap : ビット列 flag からカテゴリごとのインデックス index を生成する

  count : 学習パターン・教師信号ごとのカウント ( カテゴリ名だけを抽出して利用する )
*/
template< class CATEGORY, class TARGET >
void FlagToMap( const std::vector< unsigned int >& flag,
                const std::map< CATEGORY, std::map< TARGET, unsigned int > >& count,
                std::map< CATEGORY, typename Matrix< CATEGORY >::size_type >* index )
{
  typedef std::map< CATEGORY, std::map< TARGET, unsigned int > > Map;

  std::vector< CATEGORY > category;
  for ( typename Map::const_iterator i = count.begin() ;
        i != count.end() ; ++i )
    category.push_back( i->first );

  for ( typename Matrix< CATEGORY >::size_type i = 0 ;
        i < flag.size() ; ++i ) {
    typename std::vector< CATEGORY >::size_type b = 0;
    for ( unsigned int f = flag[i] ; f != 0 ; f >>= 1 ) {
      if ( ( f & 1 ) != 0 )
        (*index)[category[b]] = i;
      ++b;
    }
  }
}

/*
  Nominal : 名義尺度専用 CHAID アルゴリズム

  テンプレート引数の CATEGORY は学習パターン、TARGET は教師信号の型をそれぞれ表す
*/
template< class CATEGORY, class TARGET >
class Nominal
{
  // 決定木の型
  typedef std::vector< DecisionTreeNode< CATEGORY, TARGET > > DecisionTree;

  DecisionTree tree_;      // 決定木
  double threshold_;       // p 値に対するしきい値
  unsigned int maxLevel_;  // 最大レベル ( 木の高さ )
  unsigned int minParent_; // 親ノードの最小要素数
  unsigned int minChild_;  // 子ノードの最小要素数

  // 分岐するノードを追加する
  void addNode( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
                typename Matrix< CATEGORY >::size_type currentNode,
                const std::vector< typename Matrix< CATEGORY >::size_type >& rowIndex );

 public:

  // 学習パターン x と教師信号 y を指定して決定木を生成する
  Nominal( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
           double threshold, unsigned int maxLevel, unsigned int minParent, unsigned int minChild );

  // 学習パターンを指定して予測値を返す
  template< class Ran >
  std::pair< TARGET, bool > predict( Ran s ) const;
};

/*
  Nominal< CATEGORY, TARGET >::addNode
  : 分岐するノードを追加する

  x : 学習パターン
  y : 教師信号
  currentNode : 現在のノード
  rowIndex : 読み込み対象の x, y に対する行番号
*/
template< class CATEGORY, class TARGET >
void Nominal< CATEGORY, TARGET >::addNode
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    typename Matrix< CATEGORY >::size_type currentNode,
    const std::vector< typename Matrix< CATEGORY >::size_type >& rowIndex )
{
  if ( rowIndex.empty() ) return;

  std::map< TARGET, unsigned int > targetCount; // 教師信号ごとのカウント
  Classify( y, rowIndex, &targetCount );

  tree_[currentNode].impurity = static_cast< double >( targetCount[GetMode( targetCount )] ) / GetSum( targetCount );

  // 教師信号が単一なら分岐しないので終了
  if ( targetCount.size() < 2 ) {
    tree_[currentNode].target = ( targetCount.begin() )->first;
    return;
  }

  // ノードのレベルが maxLevel_ に達した / データ数の合計が minParent_ より小さい場合は
  // 教師信号の代表値を登録して終了
  if ( ( maxLevel_ != 0 && tree_[currentNode].level >= maxLevel_ ) ||
       ( minParent_ != 0 && GetSum( targetCount ) < minParent_ ) ) {
    tree_[currentNode].target = GetMode( targetCount );
    return;
  }

  // p 値が最小となる列番号とカテゴリ分割を取得する
  std::vector< unsigned int > flag;
  typename Matrix< CATEGORY >::size_type minCol = GetMinPCategory( x, y, rowIndex, targetCount, &flag, threshold_ );

  // カテゴリが分割されない場合は分岐できないので教師信号の代表値を登録して終了
  if ( flag.size() < 2 ) {
    tree_[currentNode].target = GetMode( targetCount );
    return;
  }

  tree_[currentNode].index = minCol;
  std::map< CATEGORY, std::map< TARGET, unsigned int > > count;
  Classify( x, y, rowIndex, minCol, &count );

  // ビット列 flag を添字の形に展開
  FlagToMap( flag, count, &( tree_[currentNode].branch ) );

  // 行番号を分ける
  std::vector< std::vector< typename Matrix< CATEGORY >::size_type > > row( flag.size() );
  for ( typename std::vector< typename Matrix< CATEGORY >::size_type >::const_iterator i = rowIndex.begin() ;
        i != rowIndex.end() ; ++i )
    row[tree_[currentNode].branch[x[*i][minCol]]].push_back( *i );

  // データ数が minChild_ より小さいノードがあった場合は分岐しない
  if ( minChild_ != 0 ) {
    for ( typename std::vector< std::vector< typename Matrix< CATEGORY >::size_type > >::const_iterator r = row.begin() ;
          r != row.end() ; ++r ) {
      if ( r->size() < minChild_ ) {
        tree_[currentNode].target = GetMode( targetCount );
        return;
      }
    }
  }

  // 分岐したノードを生成する
  std::set< CATEGORY > changed;
  for ( typename std::vector< std::vector< typename Matrix< CATEGORY >::size_type > >::size_type r = 0 ;
        r < row.size() ; ++r ) {
    tree_.push_back( DecisionTreeNode< CATEGORY, TARGET >() );
    tree_.back().level = tree_[currentNode].level + 1;
    // branch のリンク先を変更
    for ( typename std::map< CATEGORY, typename DecisionTree::size_type >::iterator b = tree_[currentNode].branch.begin() ;
          b != tree_[currentNode].branch.end() ; ++b ) {
      if ( changed.find( b->first ) != changed.end() ) continue;
      if ( b->second == r ) {
        b->second = tree_.size() - 1;
        changed.insert( b->first );
      }
    }
    addNode( x, y, tree_.size() - 1, row[r] );
  }
}

/*
  Nominal< CATEGORY, TARGET > コンストラクタ

  x : 学習パターン
  y : 教師信号
*/
template< class CATEGORY, class TARGET >
Nominal< CATEGORY, TARGET >::Nominal
  ( const Matrix< CATEGORY >& x, const std::vector< TARGET >& y,
    double threshold, unsigned int maxLevel, unsigned int minParent, unsigned int minChild )
  : threshold_( threshold ), maxLevel_( maxLevel ), minParent_( minParent ), minChild_( minChild )
{
  ErrLib::CheckLinearModel( x, y.begin(), y.end() );

  tree_.push_back( DecisionTreeNode< CATEGORY, TARGET >() );
  tree_.back().level = 0;

  // 読み込み対象の行列を連番で初期化
  std::vector< typename Matrix< CATEGORY >::size_type > rowIndex( y.size() );
  GenerateSerialNumber()( rowIndex.begin(), rowIndex.end() );

  addNode( x, y, 0, rowIndex );
}

/*
  Nominal< CATEGORY, TARGET >::predict : 学習パターンを指定して予測値を返す

  テンプレート引数の Ran はランダムアクセス反復子であることを想定している
  指定した学習パターンが決定木に存在しなかった場合は戻り値のペアの第二引数は false になる
*/
template< class CATEGORY, class TARGET >
template< class Ran >
std::pair< TARGET, bool >
Nominal< CATEGORY, TARGET >::predict( Ran s ) const
{
  return( Predict( tree_, s ) );
}

関数 StirlingNumber_2nd は第二種スターリング数を求めるための関数です。第二種スターリング数 S( c, g ) は以下の式で求めることができます。ここで c は学習パターンの数、g はマージ後にできるグループの数をそれぞれ表します ( 補足 2 )。

S( c, g ) = Σi{0→g-1}( ( -1 )i( g - i )c / i!( g - i )! )

OptimizeCategory は前述の MergeCategory と SplitCategory を使って最適な分割の仕方を探索する関数です。戻り値として、最終的に得られた分割の仕方で計算した p 値に第二種スターリング数を掛けた値を返します。この関数を使って各学習パターンごとに最適な分割の仕方を調べ、その中から p 値が最小となる学習パターンを求めるのが GetMinPCategory です。

CHAID アルゴリズムで決定木を生成するためのクラスは Nominal です。決定木を生成するためのメイン・ルーチンはメンバ関数 addNode で、処理内容は今まで紹介したアルゴリズムと大きな差はなく、分割する学習パターンを決めて子ノードを生成し、生成したノードに対して再帰的に呼び出しを行う流れになっています。木の高さが指定した最大値に達していた場合、ノードにある要素数が指定した最小値より小さい場合、分割後の要素数が指定した最小値より小さい場合はいずれも処理を中断するようにしています。また、指定値はゼロならいずれも判定をしないような仕様にしてあります。

サンプル・プログラムを使って最初の節の例で決定木を生成した結果を以下に示します。p 値に対するしきい値は 0.1 とし、その他のしきい値は全て無効にして処理しています。各ノードに示した数値は p 値ではなく、そのノードに含まれる教師信号の中の最大の要素数を総数で割った値です。

図 4-1. CHAID アルゴリズムによる決定木の生成 ( 名義尺度 )
CHAIDアルゴリズムによる決定木の生成

p 値を小さくするほど学習パターンはマージされやすくなるので、決定木は単純化されていきます。実際、p 値を 0.01 にまで小さくした場合、木は一つのノードだけにまで単純化されてしまいます。

サンプル・プログラムでは学習パターンが名義尺度である場合だけに対応していましたが、参考文献の中では順序尺度での処理方法についても記載されています。順序尺度の場合は、ソートした上で隣接する学習パターンのみマージを許可するような扱いになります。これは、名義尺度の場合のプログラムをほんの少し変更するだけで対応ができます。また、ボンフェローニ補正を行う場合の乗数も第二スターリング数とは異なります。マージは隣接するものどうしでしかできないため、ソート済みの要素に対して各グループの始めの要素を選択する仕方を数えればよく、しかも最初のグループは必ず最初の要素が始めになるので、ひとつ減らした残りの要素・開始点について考えればよいことになります。したがって、c 個の学習パターンを g 個のグループにマージする方法は c-1Cg-1 個になります。実際、4 つの学習パターン a, b, c, d を 3 つのグループにマージする方法は 4-1C3-1 = 3 通りで、{ a, b } { c } { d }、{ a } { b, c } { d }、{ a } { b } { c, d } となります。

参考文献ではもう一つの処理方式についても紹介しています。これは、順序尺度に一つの名義尺度 ( 通常は欠損値や無効値など ) を加えた場合で「浮動予測子 ( Floating Predictor )」と呼ばれています。このときは、一つの名義尺度のみ他の学習パターンとの全ての組み合わせで p 値を求めることになります。また、ボンフェローニ補正乗数は、学習パターンが ( 名義尺度のものも含めて ) c 個あり、これを g 個のグループに分けるとき、まず名義尺度だけが一つのグループである場合のみを考えると、残り c - 1 個の順序尺度で g - 1 個のグループを生成するので c-2Cg-2 通りあり、順序尺度のグループの中に名義尺度がある場合、まず名義尺度を除いて c - 1 個の順序尺度で g 個のグループを生成するのに c-2Cg-1 通りあって、一つの組み合わせに対して名義尺度を入れる方法が g 通りあるので、全部で gc-2Cg-1 となります。したがって、合計で c-2Cg-2 + gc-2Cg-1 となります。3 つの順序尺度 a, b, c と名義尺度 f を 3 つのグループにマージする場合、f を単独のグループにしたとき 4-2C3-2 = 2 通りで { a, b } { c } { f } と { a } { b, c } { f }、f を混在させる場合は { a, f } { b } { c }、{ a } { b, f } { c }、{ a } { b } { c, f } で 3 x 4-2C3-1 = 3 通りとなるので合わせて 5 通りあります。

年代が不明である場合を含めた分割表が次のような表になるとして、χ2 と p 値を計算してみます。但し、年代は順序尺度として扱います。

表 4-2. 年代別分割表 ( 浮動予測子 )
10代20代30代不明計
いちご282315
りんご21148
バナナ10337
計5961030

年代は隣接するペアのみで計算するので { 10 代, 20 代 }、{ 20 代, 30 代 } で χ2 を求めます。名義尺度の "不明" は全てと組み合わせるので { 10 代, 不明 }、{ 20 代, 不明 }、{ 30 代, 不明 } で計算しなければなりません。その結果は以下のようになります。

表 4-3. χ2 ( 対角より右上 ) とその p 値 ( 対角より左下 )
10代20代30代不明
10代4.130.225
20代0.1276.257.04
30代0.04391.07
不明0.8940.02960.587

この場合、p 値が最大になるのは { 10 代, 不明 } となるので、まずはこの二つをマージすることになります。

今回は、サンプル・プログラムの量が比較的多くなりましたが、それでも全ての条件に対応しているわけではありません。CHAID は名義尺度だけの対応であり、CART も学習パターンに名義尺度と連続値が混在している場合には未対応です。いずれの場合もサンプル・プログラムを少し変更するだけで対応可能なので、ぜひチャレンジしてみて下さい。

補足 1) 情報利得が非負数であることの証明

情報利得は必ず非負数となることが証明できます。親ノードにある教師信号 yi が子ノードに分割されたとき、子ノード k での yi の出現率は k にある教師信号の総数に対する比率となることから条件付き確率 P( yi | k ) で表されます。親ノードでの yi の出現率を P( yi )、全体に対する子ノード k の比率を P( k ) とすると、エントロピーを利用したときの情報利得 G は

G = Σi( -P( yi ) log2 P( yi ) ) - Σk( P( k )Σi( -P( yi | k ) log2 P( yi | k ) ) )

となります。P( yi ) = Σk( P( yi ∩ k ) )、P( k )P( yi | k ) = P( yi ∩ k ) であることから上式は

G=Σk( Σi( -P( yi ∩ k ) log2 P( yi ) ) ) - Σk( Σi( -P( yi ∩ k ) log2 P( yi | k ) ) )
=-Σk( Σi( P( yi ∩ k ) log2( P( yi ) / P( yi | k ) ) ) )

と変形できるので、もう一度 P( yi ∩ k ) = P( k )P( yi | k ) に戻して符号を左辺に移し、

-G = Σk( P( k )Σi( P( yi | k ) log2( P( yi ) / P( yi | k ) ) ) )

とします。ここで、log2 x は上に凸な関数であることから「イェンセンの不等式 ( Jensen's Inequality )」より

Σi( P( yi | k ) log2( P( yi ) / P( yi | k ) ) )≤log2 Σi( P( yi | k )P( yi ) / P( yi | k ) )
=log2 Σi( P( yi ) ) = 0

なので -G ≤ 0 すなわち G ≥ 0 となります。

ジニ不純度の場合、

G=Σi( P( yi )[ 1 - P( yi ) ] ) - Σk( P( k )Σi( P( yi | k )[ 1 - P( yi | k ) ] ) )
=Σi( P( yi ) - P( yi )2 ) - Σk( Σi( P( yi ∩ k ) - P( k )P( yi | k )2 ) )
=1 - Σi( P( yi )2 ) - [ 1 - Σk( Σi( P( k )P( yi | k )2 ) ) ]
=Σi( Σk( P( k )P( yi | k )2 ) - P( yi )2 )
=Σi( Σk( P( k )P( yi | k )2 ) - Σk( P( k )P( yi | k ) )2 )

と変形することができるので、P( k ) = pk, P( yi | k ) = qk と表して

Σk( pkqk2 ) - Σk( pkqk )2

の最小値がゼロ以上になることを証明すればよいことになります。上式は

Σk( pkqk2 ) - Σk( pkqk )2
=Σk( pkqk2 ) - Σk( pkqkΣi( piqi ) )
=Σk( pkqk2 - pkqkΣi( piqi ) )
=Σk( pk( 1 - pk )qk2 - pkqkΣi{i≠k}( piqi ) )

と変形できるので、0 < pk < 1 ( pk = 0 の場合その学習パターンは最初から除外できるので pk ≠ 0 とします ) ならば qk2 の係数は正数となって、与式は qk を変数とする楕円体になります。したがって、極値を求めればそれが最小値であり、qi で偏微分すると

2piqi - 2piΣk( pkqk )

なので、これがゼロになるとき

qi = Σk( pkqk )

という結果が得られます。これは、qi が i によらず全て等しいことを意味するので、その値を Q とすると

Σk( pkqk2 ) - Σk( pkqk )2
=Σk( pkQ2 ) - Σk( pkQ )2
=Q2Σk( pk ) - Q2Σk( pk )2 = 0

なので、最小値はゼロとなります。

最後に、標本分散の場合を証明します。標本分散を使った情報利得は

G = ( SS / N - m2 ) - Σk( ( Nk / N )( SSk / Nk - mk2 ) )

と表すことができます。但し、SS は全要素の二乗和 Σi( xi2 )、SSk はクラス k に属する要素の二乗和 Σi{xi∈k}( xi2 )、N は全要素数、Nk はクラス k に属する要素の数、m は全要素の平均 S / N = Σi( xi ) / N、mk はクラス k に属する要素の平均 Sk / Nk = Σi{xi∈k}( xi ) / Nk をそれぞれ表します。上式は

( SS / N - m2 ) - Σk( ( Nk / N )( SSk / Nk - mk2 ) )
=[ SS / N - Σk( SSk / N ) ] + Σk( ( Nk / N )mk2 ) - m2
=Σk( ( Nk / N )( Sk / Nk )2 ) - [ Σk( Sk ) / N ]2
=Σk( Sk2 / NkN ) - [ Σk( Sk ) / N ]2

となり、さらに変形すると

=Σk( Sk2 / NkN ) - Σk( Σi( SiSk ) ) / N2
=Σk( Sk2 / NkN - SkΣi( Si ) / N2 )
=Σk( ( 1 / Nk - 1 / N )Sk2 / N - SkΣi{i≠k}( Si ) / N2 )

なので、Nk < N より Sk2 の係数は正数となって、与式は Sk を変数とする楕円体であり極値が最小値となります。Si による偏微分は

2Si / NiN - 2Σk( Sk ) / N2 = 2Si / NiN - 2S / N2

なので、これがゼロになるとき

Si = SNi / N

であり、このときの情報利得は

Σk( ( SNk / N )2 / NkN ) - ( S / N )2
=( S2 / N3 )Σk( Nk ) - ( S / N )2
=( S / N )2 - ( S / N )2 = 0

となるので最小値はゼロとなります。

補足 2) ベル数 ( Bell Number ) と第二種スターリング数 ( Stirling Numbers of the Second Kind )

N 個の要素を任意の数にまとめる仕方を数えてみます。要素が 1 個なら明らかに 1 です。2 個の場合、{ a, b } と { a } { b } の 2 通りあります。3 個になると、{ a, b, c }、{ a, b } { c }、{ c, a } { b }、{ b, c } { a }、{ a } { b } { c } の 5 通りです。

このような、N 個の要素をまとめる方法の総数を「ベル数 ( Bell Number )」といい、BN で表します。上記の結果から B1 = 1、B2 = 2、B3 = 5 になります。

N = 4 のときはベル数は更に大きくなります。全ての可能性を列挙する代わりに、4 つの要素 a, b, c, d の最初の要素 a を固定して考えます。まず、a を単独で一つのグループとしたとき、残りの要素 b, c, d をまとめる仕方は B3 通りあります。次に、a を含むグループの要素数が 2 つである場合を考えます。例えば { a, b } の場合、残りの 2 つの要素 c, d をまとめる方法は B2 通りです。a を含むグループの要素数が 2 つである場合は他に { a, c } { a, d } の合わせて 3 通り ( = 3C1 ) あります。このとき、a を含むグループの要素数が 1 つであるときと 2 つであるときでは重複は必ず発生しないことに注意して下さい。次に、a を含むグループの要素数が 3 つである場合、例えば { a, b, c } を考えると、残りは d のみなので 1 通り ( = B1 ) となります。a を含むグループの要素数が 3 つである場合は他に { a, b, d } { a, c, d } の 3 通り ( = 3C2 ) あります。このときも重複は発生しません。最後に、a を含むグループの要素数が 4 つである場合 { a, b, c, d } があるので、合計は

B4 = B3 + 3C1B2 + 3C2B1 + 1 = 15

という結果が得られます。この考え方を用いれば、一般に BN は

BN = BN-1 + Σk{1→N-2}( N-1CN-k-1Bk ) + 1

という漸化式で得られますが、N-1CN-k-1 = N-1Ck であることと、N-1CN-1 = N-1C0 = 1 であることを利用すれば、B0 = 1 とすることによって

BN = Σk{0→N-1}( N-1CkBk )

ともう少しすっきりした形で表すことができます。

N 個の要素に 0 から k - 1 までの番号を付ける方法は kN 通りあり、各々の番号の付け方は N 桁の k 進数で表すことができます。例えば、3 個の要素に 0 から 2 までの番号を付ける方法は

000 001 002 010 011 012 020 021 022
100 101 102 110 111 112 120 121 122
200 201 202 210 211 212 220 221 222

の 33 通りあります。これは、N 個の要素を k 個の異なるグループに割り振る方法の数に該当しますが、場合によっては要素の存在しないグループが生じることになるのでそれらを除外することを考えます。

まず、0 のない数は 1 から k - 1 までの数で構成されているので ( k - 1 )N 個あります。1 から k - 1 まで同様に考えると、少なくとも 1 つの数を含まない N 桁の k 進数は k( k - 1 )N 個以下になります。ところが、この中には重複が発生しています。例えば 4 進数 4 桁の数 1122 は 0 を含まない数であると同時に 3 も含みません。したがって 2 つの重複が発生してしまいます。また、4 進数 4 桁の数 1111 は 0 と 2 と 3 を含まないので 3 つの重複が発生しています。一般に、r 個の数を含まない k 進数は r 個の重複が発生していることになります。
少なくとも 2 つの数を含まない N 桁の k 進数は、k 個の数の中から 2 つを選んで除外した数から生成できるので kC2( k - 2 )N 個以下です。これらのつじつまを合わせるために、まずは kN から k( k - 1 )N を引いて数を 1 つだけ含まない k 進数を除外し、数を 2 つ含まない重複分は kC2( k - 2 )N を加えて調整します。しかし、加えた分の中には 3 つ以上の数を含まない k 進数が重複しています。例えば、5 進数 5 桁の数 11222 は { 0, 3 } { 0, 4 } { 3, 4 } を含まない組み合わせの中にあるので 3 つの重複が発生しています。一般に、r 個の数を含まない k 進数は r 個の中から 2 つの数を選択する組み合わせ分だけあるので rC2 個の重複が発生しています。

この、kCi( k - i )N ( i = 0, 1, 2, ... k - 1 ) を足したり引いたりする操作を繰り返したとき、重複分がどのようになるかを考えます。r 個の数を含まない k 進数は最初一つだけあって、r 個引かれ、rC2 個加えられ... と繰り返されます。やがて i = r になったときに重複分はなくなり ( = rCr 個 )、その後は発生しなくなるので、重複分の合計は

1 - r + rC2 - ... + (-1)rrCr = Σi{0→r}( (-1)irCi )

です。ところが、二項定理から

Σi{0→r}( (-1)irCi ) = Σi{0→r}( rCi(-1)i・1r-i ) = ( 1 - 1 )r = 0

なので重複分は打ち消されます。最後まで残るのは r = k - 1 のときで、これらはただ一つの数からなる k 進数となります。以上のことから、N 個の要素を k 個の異なるグループに空のグループがないように分割する方法の数は

Σi{0→k-1}( ( -1 )ikCi( k - i )N )

となります。k - i → i と変換すると、項は ( -1 )k-ikCk-iiN となって、kCk-i = kCi なので

Σi{1→k}( ( -1 )k-ikCiiN )

と表すこともできます。

先述の例では、少なくとも一つの数を含まない数が 3・23 = 24 個以下であって、実際には

111 112 121 122 211 212 221 222 ( 0 を除いた場合 )
000 002 020 022 200 202 220 222 ( 1 を除いた場合 )
000 001 010 011 100 101 110 111 ( 2 を除いた場合 )

となります。この中で 2 つ重複しているのが 000, 111, 222 の 3 つなので

33 - 24 + 3 = 6

という結果が得られます。実際の値は

012 021 102 120 201 210

となります。

上の例を見てわかるように、各値は区別のあるグループに分割した場合の数を表しています。グループに区別がなく分割の仕方だけを考える場合はその重複分も除外する必要があります。k 個のグループに分割した上でそれを順番に並べる並べ方は k! 通りなので、k! で割った値

Σi{1→k}( ( -1 )k-ikCiiN ) / k!

は N 個の要素を k 個に分割する場合の数を表すことになります。この数を「第二種スターリング数 ( Stirling Numbers of the Second Kind )」といい S( N, k ) で表します。

ベル数の定義から、明らかに以下の式が成り立ちます。

BN = Σk{1→N}( S( N, k ) )
<参考文献>
  1. 「画像認識」 原田 達也 著 (講談社)
  2. 「決定機による判別と予測」
  3. 「CART Algorithm」
  4. 「CHAID and Earlier Supervised Tree Methods」
  5. 「Carnegie Mellon University」-「Geoffrey J. Gordon」-「Homework 4」 ... エントロピーを使った情報利得が非負値であることの証明が載っています。
  6. 「生物測定学研究室」-「生物統計学 3」
  7. Wikipedia
[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system