AdaBoost

ぐーぐる様にでも聞けばいくらでもあるんじゃないかと思ったけど意外とないんですよね。あれば楽だったんですけど・・・。というわけで主に自分用メモとして（＾＾；

まずそもそも何をやるものかというと、ある物体の特徴量から分類を行なう方法の一つ。例えば、ＣＴ値の分布やエッジ、体表からの距離を特徴量として、骨か血管かの区別を行なう。その為の前準備として正解データ(特徴量の値とその結果)を大量に用意してそれを元に学習を行ない、区別するための関数(入力は特徴量の値で出力がbool)を作り出す。
問題は正解データからどうやって区別するための関数(判別機というらしい)を作るか。その方法を説明する。

まずは正解データの各特徴量における弱判別機を作る。弱判別機というのは正解率が50％超える判別機のこと。なんでもよい。適当に線を引いて「こっから上がAで下がB」とかでもよい。それで50％超えるならそれでOK。そういうのを各特徴量においてそれぞれ作る。
次に、それぞれの正解データについて作った弱判別機に通してみて、一番精度のよい弱判別機を持ってくる。そして、その判別機で判別が失敗している正解データに関して重みを増やす。そしてその重み付きの状態でまた弱判別機を作る。今度は前回間違ったデータに重みがついているので弱判別機は精度がよくなるはず。「前回これが間違ったんだからね。今度間違ったらでかいよ」というわけである。
これを全部の特徴量に関して行なったら終了である。最後に全部の弱判別機を足し合わせたのが学習結果の判別機である。

こういった、弱判別機を何度も修正してやる学習方法をBoostingと言うらしい。AdaBoostはその一種らしい。詳しいことはよくわからん（ぉ

定量的なのを抜粋(謎)

まずは正解データを用意。ｘを特徴量の値、ｙが分類(１と-1)、特徴量の数をn、正解データの数をmとすると

となる。目指す判別機ｆは

が成り立つようなｆを作ることである。
それぞれの正解データに関する重みを作る。最初は正解率とか全然わからんのですべて同じにする。

まずは正解率50％の弱判別機lを作る。正確にはｌはスコアを出す関数で弱判別機はこのスコアの符号である。弱判別機は特徴量の数だけ必要である。弱判別機の製造方法はいろいろあるので割愛（ぉ
それぞれの弱判別機のエラー率(正解率の逆)を導出する。

このエラー率が一番低いｎ(判別機番号・・・ていうか特徴量)について重みを更新する。

以下繰り返し。

とまあこれじゃ全然わからんだろうからサンプルプログラム書いてみました。内容は適当な２つの領域をつくって、ある点がどちらの領域に入るかという内容。特徴量はある方向成分に関する射影成分。つまりは射影成分の値からどちらの領域に含まれるかを判別するというもの。

#include 
#include 

using namespace std;


//XとYの領域
const double g_fDistributionArea = 36.0f;
//正解データの個数
const int g_nCorrectSize = 1000;
//特徴量（射影方向）の数。０〜180度までで等分割
const int g_nFeatureSize = 36;
//量子化の数
const int g_nQuantumSize = 20;



//正解領域
bool IsInCorrectArea( double x, double y ){
	//適当（ぉ
	const double r0 = g_fDistributionArea / 2;
	const double ox0 = g_fDistributionArea / 2;
	const double oy0 = g_fDistributionArea / 2;
	const double r1 = g_fDistributionArea / 8;
	const double ox1 = g_fDistributionArea / 3;
	const double oy1 = g_fDistributionArea / 3;
	const double r2 = g_fDistributionArea / 8;
	const double ox2 = g_fDistributionArea / 3;
	const double oy2 = g_fDistributionArea / 3 * 2;
	const double b0 = ( x - ox0 ) * ( x - ox0 ) + ( y - oy0 ) * ( y - oy0 ) - ( r0 * r0 );
	const double b1 = ( x - ox1 ) * ( x - ox1 ) + ( y - oy1 ) * ( y - oy1 ) - ( r1 * r1 );
	const double b2 = ( x - ox2 ) * ( x - ox2 ) + ( y - oy2 ) * ( y - oy2 ) - ( r2 * r2 );
	return b0 * b1 * b2 < 0;
}

//特徴量生成関数
//０〜180°での射影成分
void CreateFeature( double x, double y, int anFeature[g_nFeatureSize] ){
	for( int nFeature = 0; nFeature < g_nFeatureSize; ++nFeature ){
		const double theta = 3.14159 * nFeature / g_nFeatureSize;
		const double reflect = cos( theta ) * x + sin( theta ) * y;
		//↑結局のところ単位方向ベクトルとの内積
		if( nFeature < g_nFeatureSize / 2 ){
			const double reflectMax = sin( theta ) * g_fDistributionArea + cos( theta ) * g_fDistributionArea;
			anFeature[nFeature] = static_cast< int >( ( reflect / reflectMax ) * g_nQuantumSize );
			//量子化して格納
		}else{
			const double reflectMax = sin( theta ) * g_fDistributionArea;
			const double reflectMin = cos( theta ) * g_fDistributionArea;
			anFeature[nFeature] = static_cast< int >( ( ( reflect - reflectMin ) / ( reflectMax - reflectMin ) ) * g_nQuantumSize );
		}
	}
}


int main()
{
	//とりあえず正解領域の表示
	const int nDistributionArea = static_cast< int >( g_fDistributionArea );
	for( int x = 0; x < nDistributionArea; ++x ){
		for( int y = 0; y < nDistributionArea; ++y ){
			if( IsInCorrectArea( x, y ) ){
				cout << "○";
			}else{
				cout << "×";
			}
		}
		cout << endl;
	}
	cout << endl;

	//正解データとその特徴量を生成
	int aanFeature[g_nCorrectSize][g_nFeatureSize];
	int aCorrectTrue[g_nCorrectSize];
	for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
		//まずは正解データ作成
		const double x = rand() * g_fDistributionArea / RAND_MAX;
		const double y = rand() * g_fDistributionArea / RAND_MAX;
		aCorrectTrue[nCorrectNum] = IsInCorrectArea( x, y ) ? 1 : -1;
		//その正解データの特徴量を格納
		CreateFeature( x, y, aanFeature[nCorrectNum] );
	}


	//重み宣言
	double afWeight[g_nCorrectSize];
	//初期値は全部同じ（^^;
	for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
		afWeight[nCorrectNum] = 1.0 / static_cast< double >( g_nCorrectSize );
	}


	//全体ループ開始
	double aafFinalScore[g_nFeatureSize][g_nQuantumSize];
	::memset( aafFinalScore, 0, g_nFeatureSize * g_nQuantumSize * sizeof( double ) );

	for( unsigned int nLoop = 0; nLoop < g_nFeatureSize; ++nLoop ){

		//現在の重みに対する最適な答えを作る。
		//便宜上配列にしているが関数というのが本当は正しい。
		//f(特徴量,その特徴量の量子化された値)
		double aafScore[g_nFeatureSize][g_nQuantumSize];
		::memset( aafScore, 0, g_nFeatureSize * g_nQuantumSize * sizeof( double ) );

		{
			//まずは重みつき分布を求める
			double aaDistributionWithWeightT[g_nFeatureSize][g_nQuantumSize];
			double aaDistributionWithWeightF[g_nFeatureSize][g_nQuantumSize];
			::memset( aaDistributionWithWeightT, 0, g_nFeatureSize * g_nQuantumSize * sizeof( double ) );
			::memset( aaDistributionWithWeightF, 0, g_nFeatureSize * g_nQuantumSize * sizeof( double ) );

			for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
				for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
					if( aCorrectTrue[nCorrectNum] > 0 ){
						aaDistributionWithWeightT[nFeatureNum][aanFeature[nCorrectNum][nFeatureNum]] += afWeight[nCorrectNum];
					}else{
						aaDistributionWithWeightF[nFeatureNum][aanFeature[nCorrectNum][nFeatureNum]] += afWeight[nCorrectNum];
					}
				}
			}

			//弱判別機生成
			//GentleAdaBoost
			for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
				for( int nQuantumNum = 0; nQuantumNum < g_nQuantumSize; ++nQuantumNum ){
					const double T = aaDistributionWithWeightT[nFeatureNum][nQuantumNum];
					const double F = aaDistributionWithWeightF[nFeatureNum][nQuantumNum];
					if( T + F != 0 ){
						aafScore[nFeatureNum][nQuantumNum] = ( T - F ) / ( T + F );
					}else{
						aafScore[nFeatureNum][nQuantumNum] = 0.0;
					}
				}
			}
		}

		{
			//新しくできた関数のエラー率を計算する。
			double aError[g_nFeatureSize];
			::memset( aError, 0, sizeof( double ) * g_nFeatureSize );

			for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
				for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
					const int nFeature = aanFeature[nCorrectNum][nFeatureNum];
					const double fScore = aafScore[nFeatureNum][nFeature];
					aError[nFeatureNum] += afWeight[nCorrectNum] * exp( -fScore * aCorrectTrue[nCorrectNum] );
				}
			}

			//重み更新
			{
				//エラーが最小の特徴量を検索
				int nMinErrorIndex = 0;
				{
					double fMinErrorNow = aError[0];
					for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
						if( fMinErrorNow > aError[nFeatureNum] ){
							fMinErrorNow = aError[nFeatureNum];
							nMinErrorIndex = nFeatureNum;
						}
					}
				}
				//更新
				for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
					const int nFeature = aanFeature[nCorrectNum][nMinErrorIndex];
					afWeight[nCorrectNum] *= exp( -aafScore[nMinErrorIndex][nFeature] * aCorrectTrue[nCorrectNum] );
				}
				//正規化
				double sum = 0;
				for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
					sum += afWeight[nCorrectNum];
				}
				for( int nCorrectNum = 0; nCorrectNum < g_nCorrectSize; ++nCorrectNum ){
					afWeight[nCorrectNum] /= sum;
				}
			}
		}

		//最終スコア分布更新
		for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
			for( int nQuantumNum = 0; nQuantumNum < g_nQuantumSize; ++nQuantumNum ){
				aafFinalScore[nFeatureNum][nQuantumNum] += aafScore[nFeatureNum][nQuantumNum];
			}
		}

		//表示
		{
			const int nDistributionArea = static_cast< int >( g_fDistributionArea );
			for( int x = 0; x < nDistributionArea; ++x ){
				for( int y = 0; y < nDistributionArea; ++y ){
					//それぞれの特徴量を計算
					int anFeature[g_nFeatureSize];
					CreateFeature( x, y, anFeature );

					//その特徴量のスコアの合計を計算
					double fScore = 0.0;
					for( int nFeatureNum = 0; nFeatureNum < g_nFeatureSize; ++nFeatureNum ){
						fScore += aafFinalScore[nFeatureNum][anFeature[nFeatureNum]];
					}

					if( fScore < 0 ){
						cout << "×";
					}else{
						cout << "○";
					}
				}
				cout << endl;
			}
			cout << endl;
		}

	}


	return 0;
}

まあこれでもまだ全然分からんと思うけど自分用メモなのでOKということで（ぉぃ