// ランダムフォレスト回帰(random forest regression)by 診断人
// データはhttp://c...content-available-to-author-only...e.com/mathematics/brase/understandable_statistics/7e/students/datasets/mlr/frames/frame.htmlのDenver Neighborhoodsを使用しています。
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <algorithm>
#include <assert.h>
#define INCMSE (1) // 重要度INCMSEを使用するか
#define INCNODEPURITY (1) // 重要度INCNODEPURITYを使用するか
using namespace std;
#define SZ(a) ((int)(a).size())
// 乱数は、xorを使ってますが、メルセンヌツイスターの方がよいかも知れません。
class RandXor
{
public:
RandXor()
{
init();
}
void init()
{
x=123456789;
y=362436069;
z=521288629;
w= 88675123;
}
inline unsigned int random()
{
unsigned int t;
t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) );
}
private:
unsigned int x;
unsigned int y;
unsigned int z;
unsigned int w;
};
static RandXor randxor; // マルチスレッド対応にするなら、木ごとに乱数用オブジェクトを用意して、シードを変えましょう。
typedef double FeatureType;
typedef double AnswerType;
struct TreeNode {
bool leaf; // 葉(=子がない)ならtrue
int level; // ノードの深さ。ルートノードは0
int featureID; // 説明変数ID。x0, x1, x2... の0,1,2の部分
FeatureType value; // 分割する値
AnswerType answer; // ノード内(=領域内)の目的変数yの平均値
vector <int> bags; // ノード内(=領域内)に含まれるデータのID
int left; // 左側の子のノードID
int right; // 右側の子のノードID
#if INCNODEPURITY
double nodePurity; // ノードのPurity。平均2乗誤差が入る
#endif // INCNODEPURITY
TreeNode() {
leaf = false;
level = -1;
featureID = -1;
value = 0;
answer = 0;
left = -1;
right = -1;
#if INCNODEPURITY
nodePurity = 0;
#endif // INCNODEPURITY
}
};
class DecisionTree
{
public:
DecisionTree() { }
// 学習。訓練データをいれて、決定木を作成する。
// features 説明変数(特徴量)x0,x1,x2...
// answers 目的変数y
// nodeSize ノード内に入るデータの数。これを下回れば計算は打ち切られ葉になる。Rの実装では、回帰のときは5、分類のときは1が使われている。
// maxLevel ノードの深さの最大値
// numRandomFeatures 領域を分けるときに試す説明変数(グラフでは軸)の数。Rの実装ではmtryという引数で、回帰のときは(説明変数の数)/3、分類のときは(説明変数の数)^0.5が使われている。
// numRandomPositions 領域を分けるときに試すデータ(グラフでは点)の数
// rootBagSizeRatio 全訓練データの中から、どれだけの割合をルートに入れるか(ブートストラップサンプリング)
DecisionTree(const vector <vector <FeatureType> >& features,
const vector <AnswerType>& answers,
int nodeSize,
int maxLevel,
int numRandomFeatures,
int numRandomPositions,
float rootBagSizeRatio)
{
const int numData = SZ(features); // データの数
const int numFeatures = SZ(features[0]); // 説明変数の数
assert(numData == SZ(answers));
assert(numData>1);
#if INCNODEPURITY
m_incNodePurity.clear();
m_incNodePurity.resize(numFeatures);
#endif // INCNODEPURITY
TreeNode root; // ルートのノード
root.level = 0;
const int rootBagSize = static_cast<int>(numData * rootBagSizeRatio); // ルートに入るデータの数
root.bags.resize(rootBagSize);
#if INCMSE
vector <int> freq(numData); // データ番号別の、ルートに選ばれた個数
#endif // INCMSE
// ブートストラップサンプリング
for (int i = 0; i < rootBagSize; i++)
{
// ここで、同じIDが選ばれる可能性があるが、問題なし。
int row = randxor.random() % numData;
root.bags[i] = row;
#if INCMSE
freq[row]++;
#endif // INCMSE
}
#if INCMSE
// 選ばれなかった
m_oob.clear();
for (int i = 0; i < SZ(freq); ++i)
{
if(freq[i]==0)
{
m_oob.emplace_back(i);
}
}
#endif // INCMSE
#if INCNODEPURITY
double sumAll = 0; // ルートのyの総和
for (int p : root.bags)
{
sumAll += answers[p];
}
const double meanAll = sumAll / SZ(root.bags); // ルートのyの平均
// ルートの平均2乗誤差
root.nodePurity = 0;
for (int p : root.bags)
{
root.nodePurity += lossFunction(answers[p], meanAll);
}
#endif // INCNODEPURITY
// 決定木をつくる。m_nodesに子ノードがどんどん追加されていく幅優先探索。
m_nodes.emplace_back(root);
int curNode = 0; // 現在の決定木のノード番号
while (curNode < SZ(m_nodes))
{
TreeNode &node = m_nodes[curNode];
// 現在のノードに入っている目的変数が、すべて同じかどうかを調べる
// (その場合は、ノードを分ける必要がなくなる)
bool equal = true; // すべて同じならtrue
for (int i = 1; i<SZ(node.bags); i++)
{
if (answers[node.bags[i]] != answers[node.bags[i - 1]])
{
equal = false;
break;
}
}
// 葉になる条件のチェック
if (equal || SZ(node.bags) <= nodeSize || node.level >= maxLevel)
{
// 葉にして子ノードは増やさない。
setLeaf(node, curNode, answers);
continue;
}
// どこで分けるのがベストかを調べる
int bestFeatureID = -1;
int bestLeft = 0, bestRight = 0;
FeatureType bestValue = 0;
double bestMSE = 1e99; // 平均2乗誤差(平均との差の2乗の和)
#if INCNODEPURITY
double bestMSELeft = 1e99; // splitValue未満の平均2乗誤差
double bestMSERight = 1e99; // splitValue以上の平均2乗誤差
#endif // INCNODEPURITY
for (int i = 0; i<numRandomFeatures; i++)
{
// x0,x1,x2...の、どの軸で分けるかを決める
int featureID = randxor.random() % numFeatures;
for (int j = 0; j<numRandomPositions; j++) // どの位置で分けるか
{
const int positionID = randxor.random() % SZ(node.bags);
const FeatureType splitValue = features[node.bags[positionID]][featureID];
double sumLeft = 0; // splitValue未満のyの総和
double sumRight = 0; // splitValue以上のyの総和
int totalLeft = 0; // splitValue未満の個数
int totalRight = 0; // splitValue以上の個数
for (int p : node.bags)
{
if (features[p][featureID] < splitValue)
{
sumLeft += answers[p];
totalLeft++;
}
else
{
sumRight += answers[p];
totalRight++;
}
}
// nodeBagのデータが"未満"か"以上"のどちらかに全部偏ってるので
// 分け方として意味がないので、すぐやめる。
if (totalLeft == 0 || totalRight == 0)
continue;
// 平均との差を使って、平均二乗誤差(MSE)を求める
double meanLeft = totalLeft == 0 ? 0 : sumLeft / totalLeft;
double meanRight = totalRight == 0 ? 0 : sumRight / totalRight;
double mseLeft = 0; // 平均二乗誤差(MSE)
double mseRight = 0; // 平均二乗誤差(MSE)
for (int p : node.bags)
{
if (features[p][featureID] < splitValue)
{
mseLeft += lossFunction(answers[p], meanLeft);
}
else
{
mseRight += lossFunction(answers[p], meanRight);
}
}
if (mseLeft + mseRight < bestMSE)
{
bestMSE = mseLeft + mseRight;
#if INCNODEPURITY
bestMSELeft = mseLeft;
bestMSERight = mseRight;
#endif // INCNODEPURITY
bestValue = splitValue;
bestFeatureID = featureID;
bestLeft = totalLeft;
bestRight = totalRight;
}
}
}
// 左か右にどちらかに偏るような分け方しかできなかった場合は、葉にする
// (すべての分け方を試すわけではないので、こういうことは起こりえます)
if (bestLeft == 0 || bestRight == 0)
{
setLeaf(node, curNode, answers);
continue;
}
// うまく分けれたので、新しい子ノードを2つ追加する
TreeNode left;
TreeNode right;
left.level = right.level = node.level + 1;
node.featureID = bestFeatureID;
node.value = bestValue;
node.left = SZ(m_nodes);
node.right = SZ(m_nodes) + 1;
left.bags.resize(bestLeft);
right.bags.resize(bestRight);
for (int p : node.bags)
{
if (features[p][node.featureID] < node.value)
{
left.bags[--bestLeft] = p;
}
else
{
right.bags[--bestRight] = p;
}
}
#if INCNODEPURITY
left.nodePurity = bestMSELeft;
right.nodePurity = bestMSERight;
m_incNodePurity[bestFeatureID] += node.nodePurity - left.nodePurity - right.nodePurity;
#endif // INCNODEPURITY
m_nodes.emplace_back(left);
m_nodes.emplace_back(right);
curNode++;
}
}
// 平均との差の2乗を求める
// y 値
// mean 平均
// 返り値 平均との差の2乗
double lossFunction(double y, double mean) const
{
return (y - mean)*(y - mean);
}
// 予測
// features テスト用の説明変数x0,x1,x2のセット
// 返り値 目的変数yの予測値
AnswerType estimate(const vector <FeatureType>& features) const
{
// ルートからたどるだけ
const TreeNode *pNode = &m_nodes[0];
while (true)
{
if (pNode->leaf)
{
break;
}
const int nextNodeID = features[pNode->featureID] < pNode->value ? pNode->left : pNode->right;
pNode = &m_nodes[nextNodeID];
}
return pNode->answer;
}
// bagsをクリアする。ランダムフォレストで一番メモリを使うのはbagsで、
// 訓練してしまえばいらないので、メモリ不足のときはクリアしたほうがよい。
void clearBags()
{
for (auto& node : m_nodes)
{
node.bags.clear();
node.bags.shrink_to_fit();
}
}
// OOBデータ
const vector <int>& getOOB() const
{
return m_oob;
}
#if INCNODEPURITY
const vector < double >& getIncNodePurity() const
{
return m_incNodePurity;
}
#endif // INCNODEPURITY
private:
// nodeを葉にして、curNodeを次のノードへ進める
void setLeaf(TreeNode& node, int& curNode, const vector<AnswerType>& answers) const
{
node.leaf = true;
// 回帰の場合は、目的変数yの平均を求める
for (int p : node.bags)
{
node.answer += answers[p];
}
assert(SZ(node.bags) > 0);
if (SZ(node.bags))
{
node.answer /= SZ(node.bags);
}
curNode++;
}
vector < TreeNode > m_nodes; // 決定木のノードたち。m_nodes[0]がルート
vector < int > m_oob; // OOBデータ(OOB Examples。訓練データから外れたもの)
#if INCNODEPURITY
vector < double > m_incNodePurity;
#endif // INCNODEPURITY
};
// ランダムフォレスト回帰
class RandomForest {
public:
RandomForest()
{
clear();
}
// 全部の決定木をクリア
void clear()
{
m_trees.clear();
}
// 学習。訓練データをいれて、決定木を作成する。
// 引数は、DecisionTreeにそのまま渡しているだけなので、そちらのコメントを参照してください。
// numAdditionalTrees 追加する木の数
void train(const vector <vector <FeatureType> >& features,
const vector <AnswerType>& answers,
int nodeSize,
int maxLevel,
int numRandomFeatures,
int numRandomPositions,
float rootBagSizeRatio,
int numAdditionalTrees)
{
for (int i = 0; i < numAdditionalTrees; i++)
{
m_trees.emplace_back(DecisionTree(features, answers, nodeSize, maxLevel, numRandomFeatures, numRandomPositions, rootBagSizeRatio));
}
}
// 予測
// features テスト用の説明変数x0,x1,x2のセット
// 返り値 目的変数yの予測値
AnswerType estimateRegression(const vector <FeatureType> &testFeatures) const
{
if (SZ(m_trees) == 0)
{
return 0.0;
}
// すべての木から得られた結果の平均をとるだけ
FeatureType sum = 0;
for (int i = 0; i < SZ(m_trees); i++)
{
sum += m_trees[i].estimate(testFeatures);
}
return sum / SZ(m_trees);
}
// 決定木のbagsをクリア(メモリ節約用)
// treeID 決定木のID
void clearBags(int treeID)
{
m_trees[treeID].clearBags();
}
#if INCMSE
// OOBエラーを求める
// trainingFeatures 訓練データの説明変数(特徴量)x0,x1,x2...
// trainingAnswers 訓練データの目的変数y
// 返り値 OOBエラー
double calculateOOBErr(const vector <vector <FeatureType> >& trainingFeatures, const vector < AnswerType >& trainingAnswers) const
{
const int numFeatures = SZ(trainingFeatures[0]);
double averageOOBErr = 0.0;
for (int i = 0; i < SZ(m_trees); i++)
{
const DecisionTree& tree = m_trees[i];
double ooberr = 0.0;
// OOBデータを使って、予測。2乗誤差の総和をooberrに入れる。
const vector <int>& oob = tree.getOOB();
for (int x = 0; x < SZ(oob); ++x)
{
const double myAnswer = tree.estimate(trainingFeatures[oob[x]]);
ooberr += (trainingAnswers[oob[x]] - myAnswer)*(trainingAnswers[oob[x]] - myAnswer);
}
if(SZ(oob)>0)
{
ooberr /= SZ(oob);
}
averageOOBErr += ooberr;
}
averageOOBErr /= SZ(m_trees);
return averageOOBErr;
}
// IncMSEを求める
// trainingFeatures 訓練データの説明変数(特徴量)x0,x1,x2...
// trainingAnswers 訓練データの目的変数y
// nPerm 特徴量1つあたりのシャッフル予測の回数
// 返り値 特徴量ごとのIncMSE
vector <double> calculateIncMSE(const vector <vector <FeatureType> >& trainingFeatures, const vector < AnswerType >& trainingAnswers, int nPerm) const
{
const int numFeatures = SZ(trainingFeatures[0]);
vector <double> errimp(numFeatures);
for (int i = 0; i < SZ(m_trees); i++)
{
const DecisionTree& tree = m_trees[i];
double ooberr = 0.0;
// OOBデータを使って、予測。2乗誤差の総和をooberrに入れる。
const vector <int>& oob = tree.getOOB();
for (int x = 0; x < SZ(oob); ++x)
{
const double myAnswer = tree.estimate(trainingFeatures[oob[x]]);
ooberr += (trainingAnswers[oob[x]] - myAnswer)*(trainingAnswers[oob[x]] - myAnswer);
}
srand(i);
for (int featureID = 0; featureID < numFeatures; ++featureID) // すべての特徴量featureIDでループ
{
double ooberrperm = 0.0;
// OOBデータを、featureID番目の特徴量だけをシャッフル
vector <FeatureType> swapped(SZ(oob));
for (int x = 0; x < SZ(oob); ++x)
{
swapped[x] = trainingFeatures[oob[x]][featureID];
}
for (int nploop = 0; nploop < nPerm; ++nploop) // 何回シャッフルしたのを試すか。
{
random_shuffle(swapped.begin(), swapped.end());
// すべてのOOBデータでシャッフルしたものを使う
for (int x = 0; x < SZ(oob); ++x)
{
vector <FeatureType> vf(trainingFeatures[oob[x]]);
vf[featureID] = swapped[x]; // featureID番目の要素だけ差し替え
const double myAnswer = tree.estimate(vf);
ooberrperm += (trainingAnswers[oob[x]] - myAnswer)*(trainingAnswers[oob[x]] - myAnswer);
}
}
const double delta = (ooberrperm / nPerm - ooberr) / SZ(oob);
errimp[featureID] += delta;
}
}
for (int featureID = 0; featureID < numFeatures; ++featureID)
{
errimp[featureID] = errimp[featureID] / SZ(m_trees); // 平均をとる
}
return errimp;
}
#endif
#if INCNODEPURITY
// IncNodePurityを求める
// 返り値 特徴量ごとのIncMSE
vector <double> calculateIncNodePurity() const
{
vector < double > retIncNodePurity;
if (SZ(m_trees) > 0)
{
const int numFeatures = SZ(m_trees[0].getIncNodePurity());
retIncNodePurity.resize(numFeatures);
for (int i = 0; i < SZ(m_trees); i++)
{
const DecisionTree& tree = m_trees[i];
const vector <double>& purity = tree.getIncNodePurity();
for (int featureID = 0; featureID < numFeatures; ++featureID)
{
retIncNodePurity[featureID] += purity[featureID];
}
}
for (int featureID = 0; featureID < numFeatures; ++featureID)
{
retIncNodePurity[featureID] /= SZ(m_trees);
}
}
return retIncNodePurity;
}
#endif // INCNODEPURITY
private:
vector < DecisionTree > m_trees; // たくさんの決定木
};
int main()
{
int numAll; // 全データ数
int numTrainings; // 訓練データ数
int numTests; // テストデータ数
int numFeatures; // 説明変数の数
// y = f(x0,x1,x2,...)
// x0,x1,x2は説明変数です。コード上ではfeatureと命名してます。
// yは目的変数です。コード上ではanswerという命名をしてます。
cin >> numAll >> numTrainings >> numTests >> numFeatures;
assert(numTrainings+numTests<=numAll);
// 全データ
vector < vector <FeatureType> > allFeatures(numAll, vector <FeatureType> (numFeatures));
vector < AnswerType > allAnswers(numAll);
for(int i = 0 ; i < numAll; ++i)
{
for (int k = 0; k < numFeatures; ++k)
{
cin >> allFeatures[i][k];
}
cin >> allAnswers[i];
}
// シャッフル用
vector < int > shuffleTable;
for (int i = 0; i < numTrainings+numTests; ++i)
{
shuffleTable.emplace_back(i);
}
random_shuffle(shuffleTable.begin(), shuffleTable.end());
// 訓練データ
vector < vector <FeatureType> > trainingFeatures(numTrainings, vector <FeatureType>(numFeatures));
vector < AnswerType > trainingAnswers(numTrainings);
for (int i = 0; i < numTrainings; ++i)
{
trainingFeatures[i] = allFeatures[shuffleTable[i]];
trainingAnswers[i] = allAnswers[shuffleTable[i]];
}
// テストデータ
vector < vector <FeatureType> > testFeatures(numTests, vector <FeatureType>(numFeatures));
vector < AnswerType > testAnswers(numTests);
for (int i = 0; i < numTests; ++i)
{
testFeatures[i] = allFeatures[shuffleTable[numTrainings+i]];
testAnswers[i] = allAnswers[shuffleTable[numTrainings+i]];
}
// ランダムフォレストを使って予測
RandomForest* rf = new RandomForest();
// 学習
const int numTrees = 100;
rf->train(trainingFeatures, trainingAnswers, 5, 999, 2, 5, 0.66f, numTrees);
// 予測と結果表示
printf("-----\n");
#if INCMSE
{
vector <double> vd = rf->calculateIncMSE(trainingFeatures, trainingAnswers, 5);
for (int i = 0; i < SZ(vd); ++i)
{
printf("IncMSE x%d=%8.2f\n",i,vd[i]);
}
}
#endif // INCMSE
#if INCNODEPURITY
{
vector <double> vd = rf->calculateIncNodePurity();
for (int i = 0; i < SZ(vd); ++i)
{
printf("IncNodePurity x%d=%8.2f\n",i,vd[i]);
}
}
#endif // INCNODEPURITY
double totalError = 0.0;
for (int i = 0; i < numTests; ++i)
{
const double myAnswer = rf->estimateRegression(testFeatures[i]);
const double diff = myAnswer-testAnswers[i];
totalError += abs(diff);
printf("Test%3d myAnswer=%8.2f testAnswer=%8.2f diff=%8.2f\n", i, myAnswer, testAnswers[i], diff);
}
printf("totalError=%8.2f\n",totalError);
delete rf;
return 0;
}