読者です 読者をやめる 読者になる 読者になる

Kaggle体験記(TopCoder機械学習マラソンとの違いなど)

tanzakuさん、今年もありがとうございます! Competitive Programming (その2) Advent Calendar 2016 - Adventarの25日目!

はじめに

Kaggle(https://www.kaggle.com/)も、TopCoder機械学習マラソン(TopCoder)も、広義の競技プログラミング! ということで、今回Kaggleに軽めに参加してみました。これらの違いについて取りあげようと思います。機械学習そのものについては取り上げてません。あしからず。

Kaggleの特徴(TopCoderと異なる点)

予測データのみを提出すれば良い。

  • TopCoderは、マッチにより、予測する実行コードを提出する形式か、予測したデータを提出する形式か違ってきますが、Kaggleは、ほぼ予測データのみ提出です(例外はある)。そういうわけで、ライブラリや計算資源を好きなだけ使えます

参加しやすい

  • ページのレイアウトが分かりやすく、いたれりつくせりで、TopCoderのように管理がされなくなったりということもなさげで(過去問が動かないとか)、全体的に、ちゃんとしてます。
  • とりあえず提出できるサンプルデータも用意されてるので、最初の1歩は踏み出しやすいです。
  • 運営の方から、「新たな検索機能を追加したので、感想をきかせてね」といった連絡がリアルタイムで送られてきて、驚きました。積極的に改善しようという姿勢がうかがえます

マッチの開催中でも、フォーラム内でなら、解法を議論してよい。

  • 解法そのものについて直接言及しているのもあります(こうすると、何点までは行けるとか)
  • データのビジュアライズなど、面倒そうなことをやってくれてる親切な人もいます。
  • ただ、上位の人は、下位の人向けへのヒントをとどめて、本当に上位にいけそうなアイディアを公開するのは控えているようです。

マッチの開催中でも、コードを公開することができる。

  • 後から始めたほうが情報が入る分有利で、公平でないとか、そういう議論もあるようです。
    • そんな人はTopCoder。TopCoderは議論禁止・コード公開禁止、また、コード提出の場合は計算資源も全員いっしょなので、そういう意味で公平な競争がしたい人にはおすすめ
      • ただ、TopCoderのようにコード提出形式のもはじまったようです。
      • 上位の情報は書いてないので、上位の人たちは、別にこのままのルールでも、公平な競争ができてると感じているかも。でも、中位以下の人のモチベーションがどうなるか気になります。

書いたコードを、Kaggleのブラウザ上で実行できる。

  • TopCoderでもできなくはないけど、ファイルも取り扱えたり。folkもできたり、高機能。多くのライブラリも使えます。Dockerのおかげ?
  • 計算資源がしょぼいのを除けば、ブラウザ上だけでなんでもできそう。すごい。

最終提出データを2種類選択できる。

  • TopCoderだと1番最後に提出した1つになるので、万が一間違ったものを提出したことを考えると、時間ギリギリで提出するのは無謀。
  • 1つは過学習したやつ、1つはクロスバリデーションしっかりやった無難なやつといった、作戦もありそうです。

参加者のレベル

  • 語れる立場にありません(泣)。マシな結果を出してからにします。ただ、最近はTopCoder→Kaggleへ参戦して結果を出しているトップクラスの方や(colunさん・Komakiさん・marek.cyganさん・iwiさんなど)、またKaggleからTopCoderへ参戦してる方もいます(fugusukiさんなど)ので、強い人はどこでも強いということだけは言えるでしょう。

Kaggle体験記じゃ

今回、参加したのはこれ。
www.kaggle.com
作業過程20時間をすべて録画していたので、時刻つきで、見直してみました。流れは分かるかと。

  • (00:00:00) なんでもいいので、最低限の解答を提出しよう。自分は大風呂敷を広げて暴走する悪癖があるので、まず簡単な解をめざす!
  • (00:31:54) Kernelsってところで、他の人のコードが見れて、ショック!folkって書いてあるってことは、これを使っていいのか?同じぐらいの得点の人もたくさんいるし、ぱくってそう…。

f:id:shindannin:20161226002752p:plain

  • (00:43:49) そのまま実行できるの?どれぐらい時間かかるんだろう…。

f:id:shindannin:20161226003012p:plain

  • (00:52:37) 提出できるし、Dockerすごい。まぁ、パクってるけど、最低限の解答ということにする。

f:id:shindannin:20161226003205p:plain

  • (02:36:56) 罪悪感をかんじてきたので、せめて自分で実行する。Linuxのほうがインストールが簡単そうだけど、家にLinuxがないので、Amazon Web Service(AWS) EC2で、無料のUbuntuインスタンスを借りてみる。

f:id:shindannin:20161226003920p:plain

  • スクリプトで使われてるパッケージ(NumPy, scipy, pandas, scikit-learn)を、全部インストールしたつもりなのに、動かず、ハマる。python 2で動くのか疑心暗鬼に。原因は、間違ってpython等を古いバージョンに戻したため、バージョンの整合性のせいで動かず、一からインストールしなおしたら、治った。 python 2で普通に全部動く。(2016年12月時点の話ですが、xgboostのインストールで、赤字の部分は古い情報も見かけるので注意。)

#xgboost
git clone --recursive https://github.com/dmlc/xgboost.git
cd xgboost; make; cd python-package; sudo python setup.py install; cd ~

  • でも、やっぱり別のエラーがでる…。スクリプトの文字コードにスペイン語が入ってるのにも関わらず、秀丸はSHIFT-JISと誤認識して、文字化けしてた。おーい…。絶対UTFで。
  • (08:30:59) 今度はMemoryエラー。単に無料インスタンスではメモリ不足のよう。有料で1番安いインスタンスを借りなおし、インストール用のbashスクリプトをちゃんと整備。

f:id:shindannin:20161226054019p:plain

  • (10:37:38) 動くようになる。ただし、単に自分でスクリプトを実行しただけでなので、スコアは同じ。この時点で、527位/1700ぐらい。
  • (11:30:22) スクリプトの中身とxgboostのパラメータなどを理解しようとする。
  • (12:30:57) データセット・スコアの理解につとめる
  • (13:27:10) フォーラムをみる。 Lagと呼ばれる特徴量を追加すると0.03点いけるらしいのでやってみる。*1

f:id:shindannin:20161226063348p:plain

  • (16:04:55) 簡易版のLagを入れてみると、実際0.03にかなり近い得点にいけた!てきと~なことこの上ないけど、上位10%まであと少し。ただ、終了まで3時間しかない…。

f:id:shindannin:20161226065853p:plain

  • ちゃんとLagを入れてみるが、点が下がる。実装ミスかもしれないけど、AWS上だとIDEはさすがに使えないので、デバッグしづらい…
  • 5分程度しか計算してないので、短すぎかも?
  • (18:09:16)もう高いEC2インスタンスC4.8XLarge借りちゃえ。1時間$3ぐらいだし、もう1時間30分しかないので、問題ないでしょう。

f:id:shindannin:20161226071750p:plain

  • うわー、それでも結果伸びない。他人のを借りてきただけど、いろいろちゃんとセットアップしないので、結局損。xgboostしか使ってないのは、そこまで大問題ではないと思うけど、特徴量選択の自動化・ハイパーパラメータ調整・クロスバリデーションも何もやってない。眠い。天啓による調整に頼る(ダメ)。
  • 最後、結構時間ギリギリ(打ち切れるように改良してないのもだめだし、そもそもギリギリまでやってるのがダメ)
  • なんとか間に合って
  • (19:38:59) あぁぁぁぁ、あせって間違ってデータではなく、pythonスクリプトを提出してることに気づく。つい、TopCoderの癖が…。バカすぎる…。

f:id:shindannin:20161226071207p:plain

  • マッチ終了
  • (19:42:19) Post Deadlineと表示される。提出遅れには厳粛に対応、しかし、遅れた結果はしっかり見せる、優しいのか鬼なのか分からないけど、ちゃんとしてる仕様。

f:id:shindannin:20161226071050p:plain

  • 最終結果は、215位/1784(上位13%)でした。TopCoderのように、レート下がるってことはなさげなので、下位に沈んだら放棄する人も多いのかも。
  • 徹夜したので、次の日は即寝落ち。




















  • 2日後、高額なインスタンスにつなげっぱなしだったのに気づく!
  • 結果

f:id:shindannin:20161226072340p:plain









  • …わしはバカすぎじゃ。

Kaggle体験記まとめ

  • 普段からLinuxを使っている人や、AWS EC2やS3を使っている人であれば、もっとすんなり行くと思います。
  • 今回は全部AWS EC2上のUbuntu上で行いましたが、自分のPCにも開発環境を用意したほうがいいと思います。IDEが使えないのが痛い…。
  • Kaggleの流れを知りたいのであれば、最初の1回は今回のように他の人のスクリプトを足がかりにスタートするのは悪くないと思います。
  • ただ、成功例をたどるだけでは、勉強にはならなそう(短時間のわりに良いスコアは得られるかもしれませんが)。試行錯誤したり、アルゴリズムを導入・改良したり、自動化など環境をより便利にしたりってところは、自分でやったほうが、長期的にみると、スコア的にも勉強の観点からも得だと思います。
  • Forumでの解法の事後公開やディスカッションもとても盛んなので、読みましょう。そのあとの復習が大事。
  • 日本人のKaggle参加者の参戦記もたくさん落ちているので、ぜひ読みましょう。

おわりに

来年こそはKaggleもがんばりたい…といいたいところですが、本職が大変になりそうなので、ぜひ、これをみた皆さんは、わしの分まで、Kaggleをやるのじゃぞ!

*1:(同じお客様IDの過去5年間のデータを探し、それも特徴量に加えるという方法です。過去の記事 ランダムフォレストのつかいかた - じじいのプログラミング -> 同種の説明変数を追加する。 -> (1)同種のデータとのかかわりを表す変数を追加する ぽいやつ)

競技プログラミングの2015年の結果、2016年の目標

2015年の結果

SRMに56.6%以上参加 20点 達成
SRMに56.6%以上参加したうえでRating1800以上 20点 失敗
TCO15 SRM Round 3 20点 失敗
マラソンマッチ レッドコーダー 80点 失敗
TCO15 マラソン予選10位以内 1回ごとに 20点 失敗
TCO15 マラソンワールドファイナル出場 300点 失敗
Google Code Jam Round 3 10点 失敗
プロコンTシャツゲット1枚ごとに 5点 2枚
プロコン賞金10万円ごとに 10点 失敗
プロコンでオンサイト決勝出場 50点 失敗
ニコ生放送200枠分(=100時間) 10点 失敗
海外でプロコンオフ会をやる 10点 達成
Codeforcesで栗原市ランキングいり 10点 失敗

2015年の結果 40点でした…。(なお、2014年は60点でした)

2015年の反省

マラソンマッチ レッドコーダー 80点 失敗
TCO15 マラソン予選10位以内 1回ごとに 20点 失敗
TCO15 マラソンワールドファイナル出場 300点 失敗

まずは、ここらへん。TCO15 Marathon Match Round 1の最終提出が間に合わず(もちろん自分が悪い)、レーティング激減して、気持ちが切れてしまった感じがあります…。あと、そもそも時間が作れなかったというのもありました(特に2015年上半期)。
また、自分の場合は、安定した結果を出すために、マッチに参加するよりも、突っ込んだ復習が大事な気がしているのですが、復習のモチベ維持もできていなかった気がします。

2016年はもっと時間がなさそうで、未来は暗い気がしますが、参加するマラソンマッチを絞って、出たやつについては丁寧にやっていこうと思います。

プロコン賞金10万円ごとに 10点 失敗

主に機械学習マッチでチャンスがあると思っていたのですが、そもそもあまり参加ができませんでした。少なくとも解答がでた場合は、結果が悪くても提出するようにはしているのですが、解答まで至らなかったマッチが2つありました…(地震と前立腺がん)。こっちは短時間で準備するノウハウはつかんできた気がするので、今後、通常マラソンよりは結果が出せそうな気はします。

SRMに56.6%以上参加したうえでRating1800以上 20点 失敗

2015年は初めて一時的にRating 1800を越え、最終戦まで1800を越えるチャンスはあったのですが、逃してしまいました。Mediumが解けないのはしょうがないにしても、Easyで落とすケースが目立ちました。というか、Easyが昔ほど簡単ではありません。
まずEasyの出てない回の問題もたまってきたので、Easy問題を練習で解いてEasy安定に努めたいと思います。これなら、時間はそんなにかからないはず。

Google Code Jam Round 3 10点 失敗

大会では、もうちょっとうまくやれた気がします。Code Jam Round 3は決して不可能ではない目標なのですが、ピンチになって難しい問題を解かないとダメな状況になるとほぼ失敗するので、ピンチになってもあせらずに簡単な問題を確実に解いていくようにしたいです。

ニコ生放送200枠分(=100時間) 10点 失敗

結果、放送したのは120枠分(=60時間)程度でした。現状だと100時間はかなり難しいですが、できる範囲で継続してやっていきます。

海外でプロコンオフ会をやる 10点 達成

シンガポールでオフ会をするのは予想以上に大変でした…。「人を集める方法」「人に連絡する方法」「場所探し」「参加費」「何をすればいいのか」、すべてが難しかったです。かなりの時間を使った気がしますが、達成できてよかったです。ただ、これをやらなければ、他に時間をつかえた気もします。

あと、数少ない良かった点としては、「目標を常に覚えていた」ことでしょう。目標がなければ、もっと無為に過ごしていたと思われるので、それよりはマシ。2016年も目標を忘れない!

2016年の得点表と目標

正直なところ、2016年は2015年より忙しくなる要素しかないので、目標は2014年の水準、60点以上にします。あと、Codeforces, Facebook Hacker Cup, Kaggleとバリエーションを増やしました。(これらを増やしても2015年は40点です。)

SRMに50%以上参加 20点
SRMに50%以上参加したうえでRating1800以上 20点
TCO16 SRM Round 3 20点
マラソンマッチ レッドコーダー 80点
マラソンマッチ Rating 2000以上 40点
マラソンマッチ Rating 1800以上 20点
TCO16 マラソン予選10位以内 1回ごとに 30点
TCO16 マラソンワールドファイナル出場 300点
Codeforces Highest Rating 2200以上(薄い橙) 50点
Codeforces Highest Rating 1900以上(紫) 10点
Google Code Jam Round 3 10点
Facebook Hacker Cup Round 2 5点
Facebook Hacker Cup Round 3 20点
プロコンTシャツゲット1枚ごとに 5点
プロコン賞金10万円ごとに 10点
プロコンでオンサイト決勝出場 50点
Kaggle 上位25% 5点
Kaggle 上位10% 10点
ニコ生放送200枠分(=100時間) 10点
海外でプロコンオフ会をやる 10点
Codeforcesで栗原市ランキングいり 10点

なお、得点表は、難度だけではなく、自分がやりたいかどうかも含めて決めています。

ところで…

ほっほー。

2015年の目標を振り返ってみたところ、40点のところをみたら、

30~59 怠惰・堕落(FUJIYAMAに乗る)

と、書いてあったのじゃ。絶対にこんな低い点は取らないと思って、てきとーに書いたのに…。しかも、不幸は重なるもので、甲府に行く用事ができてしまったので、ついでに行ってきたのじゃ…。
f:id:shindannin:20160105160851j:plain

頭がおかしいアトラクションだらけですね。わしは死にたくないのじゃ。
f:id:shindannin:20160105145828j:plain

でも勇気を振り絞って、以下の乗り物に乗ってきたのじゃ。天にも届くようなスロープに、背後にそびえる富士山が、一層の恐怖をかきたてるのじゃ。変な生き物もいたし、地獄だったのじゃ。さすがFUJIYAMA!
f:id:shindannin:20160105152400j:plain

…では、さらばじゃ。ほっほー()

累積和を使う動的計画法

この記事はCompetitve Programming Advent Calendar 2015の23日目の記事です。tanzakuさんに感謝
www.adventar.org


今回は、累積和を使う動的計画法についてです。TopCoderのDiv2上位ぐらいの人向けの難易度です。

問題

AtCoder Typical DP Conestの問題です。
F: 準急 - Typical DP Contest | AtCoder

Problem Statement
ある路線には駅 1 から駅 N までの N 個の駅がある。すぬけ君は、この路線に準急を走らせることにした。
準急は、駅 1 に止まり、{駅 2, ..., 駅 N-1} の部分集合に止まり、駅 N に止まる。
連続する K個以上の駅に止まると、客が飽きてしまうので、そのようなことはしない。
準急の停車駅の組み合わせとして何通り考えられるか、mod 1,000,000,007 で求めよ。


Constraints
2≤K≤N≤1000000

最初の方針(ただし、計算時間が間に合わない)

まず答えがでる実例N=9, K=5で手計算すると、以下のように112通りになります。図は矢印にそって足していくだけです。駅1と駅9は絶対に停まるのには注意してください。
f:id:shindannin:20151224005704p:plain:w700

  • 状態は、
    • 何駅目まで来たか
    • 現在何駅連続で停まっているか
  • 遷移は
    • 次の駅に停まる→連続で停まった回数が+1
    • 次の駅に停まらない→連続で停まった回数が0に戻る

となるので、動的計画法で
dp[何駅目まで来たか][現在何駅連続で停まっているか]
とすれば、コードは以下のようになります。

#include <iostream>
using namespace std;

long long dp[123][123];

int main()
{
	int N,K;
	cin >> N >> K;

	long long  MOD = 1000000007;

	dp[1][1]=1;
	for (int n = 1; n <= N-1; ++n)
	{
		for(int k = 0; k <= K-1; ++k)
		{
			dp[n+1][k+1] += dp[n][k];  // 駅にとまる
			dp[n+1][k+1] %= MOD;
			dp[n+1][0]   += dp[n][k];  // 駅にとまらない
			dp[n+1][0]   %= MOD;
		}
	}

	long long ans = 0;
	for(int k = 1; k <= K-1; ++k )
	{
		ans += dp[N][k];
		ans %= MOD;
	}

	cout << ans << endl;

	return 0;
}

しかし、この問題ではNもKが1000000まで取りえます。計算時間もメモリも間に合いません。

累積和で状態数を減らす

動的計画法の計算量は、ざっくりいうと、以下のように状態数×(1状態あたりの)遷移数になります*1
f:id:shindannin:20151224005712p:plain:w700
つまり、計算量減らして高速化するには、

  1. 状態数を減らす
  2. 遷移数を減らす

のが基本です。今回は、状態数が多いので、状態数を減らしましょう。1連続停車~K-1連続停車までを1つの状態にします*2
f:id:shindannin:20151224005718p:plain:w700
駅に停まらないケース(青線)は、矢印元と矢印先が同じ状態なので、そのまままとめられます。簡単。
f:id:shindannin:20151224005722p:plain:w700
ただし、駅に停まるケース(赤線)はどうでしょう?1番上のK-1回連続(先の例だと4連続)のケースは、駅にこれ以上停まれないので失われています。この部分は差し引かないとダメですが、状態をまとめてしまってはK-1連続停車の回数がわかりません。この部分をどうやって求めるかが、累積和を使った動的計画法のポイントです。
f:id:shindannin:20151224005727p:plain:w700
最初の図をよく見ると、回数が斜め方向に同じになることがわかります。改めて考えてみると、n駅でK-1連続停車した回数は、n-(K-1)駅で0連続停車した回数といっしょです。
f:id:shindannin:20151224005733p:plain:w700
0連続停車した回数は分かるので、これで計算できます。例えば駅5の1連続~K-1回連続は、4+4-1=7と求められます。
f:id:shindannin:20151224005738p:plain:w700
これで計算量はO(N)まで落とせます。コードは以下のようになります。

#include <iostream>
using namespace std;

long long dp[1234567];  // dp[n] = n駅まできて、現在0連続停車
long long sum[1234567]; // sum[n] = n駅まできて、現在1~K-1連続停車

int main()
{
	int N,K;
	cin >> N >> K;

	long long  MOD = 1000000007;

	dp[0]=1;
	sum[1]=1;
	for (int n = 1; n <= N-1; ++n)
	{
		dp[n+1]=sum[n]+dp[n];
		dp[n+1]%=MOD;
		sum[n+1]=sum[n]+dp[n];
		sum[n+1] %= MOD;
		if(n+1>=K)
		{
			sum[n+1]=sum[n+1]-dp[n+1-K]+MOD;
			sum[n+1] %= MOD;
		}
	}

	cout << sum[N] << endl;

	return 0;
}

まとめ

動的計画法の高速化には、状態数を減らすのが方法の1つです。累積和を使って減らすことができます。累積和の部分と累積和ではない部分の関係をどう立式するかで、アイディアが必要になってきます。
ただ、今回の方針は一例にすぎず、同じ問題でも多数のやり方があります。ぜひ「Typical DP 準急」で検索して、他の方針との違いも見てみてください。

応用問題として、ちゃっかりこの問題をあげておきます。累積和だけではありませんが、チャレンジしてみてください。
Contest Page | CodeChef


あす24日クリスマスイブは、競技プログラミングアドベントカレンダー創始者のtanzakuさんと、実際に役立った例を紹介してくださるuvarimeronさんです。楽しみですね。それでは、楽しいクリスマスを!

*1:本当は1遷移あたりの計算量も掛けないといけませんが、O(1)の問題のほうが多いでざっくり

*2:着想は、自分の場合は、1連続停車~K-1連続停車まで、だいたい同じような遷移をしているところから思いつきました。

MSXの思い出 - 別冊テープログイン MOCKY

MSXアドベントカレンダー8日目の記事です。主催者のmsiroさん、ありがとうございます!
www.adventar.org

1985年の別冊テープログインに載っていた、MOCKYというゲームを紹介しようと思います。

Tagoo : MSXソフトウエア検索 : 別冊テープログインMSX GAME BOOK
f:id:shindannin:20151209001724p:plain

このゲームは風船をつくり
f:id:shindannin:20151208235356p:plain:w400
敵をおびき寄せ
f:id:shindannin:20151209000000p:plain:w400
敵がからまったら
f:id:shindannin:20151209000017p:plain:w400
敵をくくりつけて
f:id:shindannin:20151209000030p:plain:w400
敵を空に飛ばして倒すゲームでした
f:id:shindannin:20151209000036p:plain:w400
2人プレイもできます
f:id:shindannin:20151209000045p:plain:w400
片方のプレイヤーが敵にやられて全滅しても
f:id:shindannin:20151209000053p:plain:w400
やられたプレイヤーはお化けになり、ゲームが続けられます
f:id:shindannin:20151209000107p:plain:w400
お化けになったプレイヤーは風船を壊すことができます
f:id:shindannin:20151209000117p:plain:w400
風船を割れるので、邪魔することができます
f:id:shindannin:20151209000128p:plain:w400
しかし、実はお邪魔プレイだけではありません

このゲーム、実は風船が同時に2個しか作れず、変なところに風船を作るとピンチになります
f:id:shindannin:20151209000136p:plain:w400
そこでお化けが、いらない風船を割ってあげると
f:id:shindannin:20151209000143p:plain:w400
風船を新たにつくることができ
f:id:shindannin:20151209000155p:plain:w400
追い詰められてしまっていても、協力して、無事敵を倒せます!
f:id:shindannin:20151209000201p:plain:w400


まとめると、

  • 死んだあともプレイできる。
  • プレイヤー1とプレイヤー2が違う役割となる
  • 死んだ後に、邪魔プレイと協力プレイの選択ができる
    • ゲームをやり直したければ、邪魔プレイ。
    • 性能の差をいかして、協力プレイ

という工夫がありました。その後、ボンバーマンのみそボンなど同じコンセプトのものは現れますが、画期的なゲームデザインでした。

明日9日目はsunflatさんのVPOKEについての記事です。楽しみですね!ではでは!

ボツネタ集(だいたいランダムフォレスト)

どうでもよいまえがき
今回の記事を最後に、機械学習関係の記事は、しばらくお休みにします!

  • 我流すぎて、機械学習の知識が不足していて、記事を書くのに妙に時間がかかる。
  • 最近、機械学習マッチで結果を出してないので、その中で記事を書いても、自分の中で説得力がいまいち。

というわけでしばらくの間は、機械学習について、以下のことに専念します。

  • 本とか論文とか読んでインプット
  • 機械学習マッチで結果を出す

ただ、記事を書こうと思ってたボツネタがあるので、勢いで書きました。雑です。メモ書き程度だと思ってもらえば。

それではスタートじゃ!!!!

ボツネタ1 : ID分解

世の中のIDというものは、複数の情報を持っていることがあります。これを別の特徴量として追加します。数学も機械も関係ない、セコいテクニックですが、意外と実戦的かもしれません。

まず、データを眺めてください。
f:id:shindannin:20150425163540p:plain

自分だったら、x0のidの赤で囲われた桁を、別な特徴量として追加します。紫で囲った部分も念のため追加するかもしれません。なぜでしょう?
f:id:shindannin:20150425163552p:plain

なぜかというと、いろいろと法則に気づくからです。

  • 紫の最初の5桁は、startdayが同じとき、常に同じ値になります。日時にかかわる値だと推測できます。ただ、この5桁とstartdayの大小関係違うので、一応特徴量として追加します。
  • 赤の桁は、完全に独立した値です。0~2の値しかとりません。とてもあやしいので、特徴量として追加します。
  • 下2桁と3桁は常に31なので不要です。
  • 下1桁はcyclesと常に同じ値なので不要です

f:id:shindannin:20150425163601p:plain

まぁ、ただ注意点は、

  • 分解前の元のIDを削るのは、特に慎重になったほうがいいです。気づいてない重要な情報があるかもしれないので。とりあえずは分解したのを追加するだけにとどめておきましょう。
  • つねに同じ値というのも、必ず全部のデータで常に同じ値になってるか確認しましょう。まれに違うことがあって、それが重要度の高い特徴量になってたら。大惨事です。

データを数値でも眺めるのは本当に重要です。機械学習する前に、必ず見るようにしましょう。まぁ、こういうのすら自動化するのが最高ですけどね。

ボツネタ2 : ランダムフォレスト実装例の簡単なバリエーションについて

過去の記事で、ランダムフォレストといってもバリエーションがいろいろあるといって逃げっぱなしなので、「ちょっとは説明しないと…」という罪悪感をかんじていました。

  • 効果 : ★が大きいほど、効果が高そう…といいたいですが、完全に自分の経験上なので、信ぴょう性はあまりないです。
  • 謎度 : ★が大きいほど、めったにみかけないもので、謎度を高くしてます。

C++の実装例を使って、ちょっと説明します。改良はみんなの宿題!

ルートバッグ1.5倍

効果 ★
謎度 ★★★★★
TopCoderのPsyhoさんのコードをみて知ったテクの中でも、もっとも謎なものです。

ブートストラッピングでは、ルートバッグのサイズは元のデータは0.632倍が使われています。
bootstrap - What is the .632+ rule in bootstrapping? - Cross Validated
実装例のrootBagSizeRatioの部分です。
115行目

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
}

rootBagSizeRatioの部分は、調整するにしても、普通は1以下の値を入れようと考えますが、Psyhoさんはここで1.5倍のように、1を超える値を入れてくることもあります。いったいどういう意味があるのか分かりませんが…。
なお、

  • 1.5倍を入れるようにした場合でも、OOBデータもわずかに残ります。重複も許してるので。
  • このままだと計算時間も長くなるので、防ぐために特徴量毎にweightという変数を用意し、何回でてきたかをとっておけば、計算時間は長くなりません。ソースコードの改良は読者の宿題!!

領域を分けるとき、全部の境界を試す

効果 ★★★
謎度

Rのランダムフォレスト実装であるので、そちらが標準だと思います。

ある特徴量の、全部の境界を試してます。自分の実装例では以下のように、numRandomPositions個しか試していません。
224行目

for (int j = 0; j<numRandomPositions; j++)	 // どの位置で分けるか
{
    const int positionID = randxor.random() % SZ(node.bags);
    const FeatureType splitValue = features[node.bags[positionID]][featureID];

ただ、多く分けるほうが必ず良いというわけでもありません。特に、少ない特徴量で似たような木ができる現象の要因になることも。計算時間も長いです。

上記の実装を全部の境界で切るようにしたいなら、

for (int j = 0; j<SZ(node.bags); j++)	 // どの位置で分けるか
{
    const int positionID = j;

とすれば、機能としてはRの実装と一緒になりますが、速度としてはお話にならないです、ダメ!前回の日記で述べた式

  {Σ[i∈P](Yi-E[P])^2}                          - {Σ[i∈L](Yi-E[L])^2}                          - {Σ[i∈R](Yi-E[R])^2}
= {Σ[i∈P]Yi^2} - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2 - {Σ[i∈L]Yi^2} + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2 - {Σ[i∈R]Yi^2} + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... 2乗を展開した
=                - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2                  + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2                  + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... {Σ[i∈P]Yi^2} = {Σ[i∈L]Yi^2} + {Σ[i∈R]Yi^2} で消した。
=                - 2E[P]S[P]         + |P|E[P]^2                  + 2E[L]S[L]         - |L|E[L]^2                  + 2E[R]S[R]         - |R|E[R]^2 ... 総和Sであらわした
=                - 2S[P]^2/|P|       + S[P]^2/|P|                 + 2S[L]^2/|L|       - S[L]^2/|L|                 + 2S[R]^2/|R|       - S[R]^2/|R| ... E[A] = S[A]/|A|で平均Eを消した
=                -  S[P]^2/|P|                                    +  S[L]^2/|L|                                    +  S[R]^2/|R|                    ... 足し算

これで評価するようにすれば、xが小さいほうから高速に全ての分け方を試せます(予めqsortを使います)。まぁ、説明不足なので、本家のソースコードregTree.cの関数findBestSplitを熟読して、実装しましょう。読者の宿題!
なお、後述する平均2乗誤差、3乗誤差や1.5乗誤差を変えるというのをやると、以上の数式展開が使えず、速くできません。

領域を分けるとき、任意の場所で分けられるようにする。

効果 ???
謎度 ★★★
一度も試したことはありません。上記の例ではデータのある場所で分けていますが、現在のノードに入ってるデータの、特徴量xの最大と最小を求めておき、(max-min)*rand + min の位置で切るといった感じの実装で、任意の場所で分けることができます。これをやるとxの値の差で選ばれる確率が変わるようになるので、いろいろと問題が起きそうな気もしますが…。

領域を分けるとき、マージン最大化っぽいことをする

効果 ★★
謎度
過去の日記で、yowaさんからコメントをいただきました。Rの実装でもやっているので、謎ではありません。
ランダムフォレストのつくりかた(C++の実装例つき) - じじいのプログラミング

おそらく日記をみて、「せっかく赤色と青色の領域を分けるのに、なぜ、わざわざギリギリの点上で分けるのか。中間で分けた方がいいのではないか」と思った人もいるのではないかと思います。それです。

改良前
278行目

node.value = bestValue;

改良後はこんなかんじになるでしょうか。

FeatureType nextValue = -BIG;

for (int i = 0; i < SZ(node.bags); ++i)
{
    if( features[node.bags[i]][bestFeatureID] < bestValue)
    {
        nextValue = max(nextValue, features[node.bags[i]][bestFeatureID]);
    }
}
node.value = (bestValue + nextValue) / 2.0;
if (!(node.value > nextValue)) node.value = bestValue;

あまり効果があった経験はないのですが、計算のオーダーが増えるわけでもないので、いれていいと思います。

木の深さに応じ、特徴量を選ぶ数numRandomFeaturesを変える

効果 ★★★
謎度 ★★★

以下のnumRandomFeaturesは、通常は固定ですが、これを木の深さに応じて変えます。

198行目

for (int i = 0; i<numRandomFeatures; i++)
{
    // x0,x1,x2...の、どの軸で分けるかを決める
    int featureID = randxor.random() % numFeatures;

実装は読者の宿題!まぁ、これは簡単でしょう。木の深さを表す変数をノードにもたせておき、numRandomFeatures[depth]みたいにすればよろしい。

これは割と効果があります。特に、後述する似た木ができる現象を避けるのに使えます。つまり、木の浅いところではnumRandomFeaturesの値を少なく、深いところではnumRandomFeaturesを大きめにすると良いです。

MSE以外の誤差を使用する

効果 ★★
謎度 ★★★

311行目

double lossFunction(double y, double mean) const
{
    return (y - mean)*(y - mean);
}

ここを、絶対値の3乗誤差、絶対値の1.5乗誤差、絶対値などに置き換える方法です。

ボツネタ3 : ランダムフォレストで似た木ばかりができると予測性能が落ちる

ランダムフォレストで、同じような木がたくさんできてしまうと、極端に予測性能が落ちることがあります。以下のような条件のとき要注意です。

  • 特徴量が少ない(10個以下)
  • 特徴量に比べて、ノード分割の際に選ぶ数numRandomFeatures(RのmTry)が多い。(目安は、回帰 特徴量の数/3、分類 特徴量の数^0.5)
  • 切る場所を試す数numRandomPositionsが大きい(Rだと全部で切るので、多い)

木をたくさん増やすことである程度解消できますが、限度があります。ちゃんとパラメータを調整して、避けましょう。

以下の図は、特徴量が7個しかないケースで、予測性能が極端に落ちたときの、決定木10本をビジュアライズしたものです。

  • 上側がルートです。
  • 帯の長さが、ノードの属するデータの個数に当たります。
  • 帯の色が、ノードをどの特徴量xで分割したかにあたります。計7色。

これを見ると、一見カラフルで、いろんな木ができてますが、上部の色はだいたい、茶色か水色か黄緑しか来ていないことが分かります。
f:id:shindannin:20150425190854p:plain

念のため、上記の画像を出力するコードも貼っておきます。htmlで出力する手抜きなやつですが…。

    void visualize(int id)
    {
        char path[1001];
        sprintf(path,"../visual/VisualTree.html");
        FILE *fp = fopen(path, id==0 ? "w" : "a");

        vector <string> colors = 
        {
            "Red",
            "Maroon",
            "Yellow",
            "Olive",
            "Lime",
            "Aqua",
            "Teal",
            "Blue",
            "Green",
        };

        if (fp)
        {
            fprintf(fp,"<canvas id=\"sample%d\" width=\"15000\" height=\"300\">\n",id);
            fprintf(fp,"<p></p>\n");
            fprintf(fp,"</canvas>\n");
            fprintf(fp,"<script>\n");
            fprintf(fp,"var canvas = document.getElementById('sample%d');\n",id);
            fprintf(fp,"var context = canvas.getContext('2d');\n");

            float xRatio = 0.5f;

            queue < pair <int, int> > q; // ノード番号, 左座標
            q.push(make_pair(0,0));
            while(!q.empty())
            {
                pair <int,int> now = q.front();
                q.pop();

                const TreeNode& nd = m_nodes[now.first];


                fprintf(fp,"context.fillStyle = \"%s\";\n", colors[(nd.featureID+SZ(colors))%SZ(colors)].c_str());
                fprintf(fp,"context.fillRect(%d,%d,%d,%d);\n", now.second, 10*nd.level, (int)(SZ(nd.bags)*xRatio), 10);

                if(!nd.leaf)
                {
                    q.push(make_pair(nd.left,now.second));
                    q.push(make_pair(nd.right,now.second+(int)(SZ(m_nodes[nd.left].bags)*xRatio)));
                }
            }
            fprintf(fp,"</script>\n");
            fclose(fp);
        }
    }

ボツネタにすらまとめられなかったボツネタ

カテゴリカル変数ごとに、他の特徴量xや目的変数の統計データを求めて、新たな特徴量として追加する。

機械学習は、基本不等式で分けるので、カテゴリカル変数(列挙型)はやっかい。そのままは使えない。ただ集計して特徴量を足すのに使いやすい。カテゴリカル変数x1をだとしたら、x1の値ごとに、他の特徴量x2,x3..やyを集計して(足す・割る・割合を求める・最大・最小)、新たな特徴量を追加する話。なぜボツにしたかというと、
ランダムフォレストのつかいかた - じじいのプログラミング
の「同種のデータとのかかわりを表す変数を追加する」とネタがかぶってるから。ただ、あの記事ではカテゴリカル変数とは言わなかったので、改めて書くか迷った。ただ非常に効果があるテクニックなので、みなさん使ってほしいです。

ブースティングのC++の実装例

C++の実装例で決定木の部分は実装してるので、てきとーなブースティング(木0の予測外れ分を木1のyとツッコむ)というのを実装するのは割と簡単。ただ、まじめなブースティングは、木の重みづけの部分とかも含めて理論になってるっぽく、そこは理解してないので、やめた。また、実戦で使ったことがないのも不安なので、やめた。







以上です。では、さようなら~。ほっほー。

ランダムフォレスト 特徴量の重要度(C++の実装例つき)

はじめに

今回の記事は、alfredplplさんの以下の記事とだいぶかぶっています…。図つきで、とても分かりやすい記事なので、お勧めです。こちらをはじめに読んだほうが良いと思います。alfredplpl.hatenablog.com
alfredplplさんの記事は分類を例にあげてましたので、こちらは回帰を例にあげて説明します。基本的な中身は一緒ですので、比較しながら読んでもらえればと思います。

Rのランダムフォレストの重要度

Rのランダムフォレストには2種類の重要度(importance)があります。この記事で説明するのは、IncMSEとIncNodePurityです。

アプローチ 回帰(regression) 分類(classification)
シャッフルしたOOBデータで予測して求める IncMSE MeanDecreaseAccuracy
決定木のノードの純度(Purity)の増分から求める IncNodePurity MeanDecreaseGini

Rでの求め方

(以下の例はTopCoder機械学習マッチ"ChildStuntedness5"でのデータを使用しました。要ログイン→TopCoder)
(1)randomForestの引数にimportance=TRUEをセットする。(注:以下はただの例ですので、他の引数とかは適宜変更してください。)

forest <- randomForest(GENIQ~., data=data.train, ntree=200, importance=TRUE)

(2) プロットしたいならvarImpPlot、しないならimportanceを呼ぶ。

ret1 = varImpPlot(forest)
ret2 = importance(forest)

結果は以下のようになります。

f:id:shindannin:20150424021604p:plain

  • IncMSEとIncNodePurityは別なので、重要度の値はもちろんのこと、上記のように順位が異なってくる場合もあります
  • 上記の方法ではなく、importance(forest)で重要度を出力すると、IncNodePurityは標準誤差で割られた値になります。*1

IncMSE

C++での実装例

C++の実装例と実行結果をideoneにアップしました。#define INCMSE (1)で使用するしないを切り替えられます。#if INCMSEで囲まれている部分を見れば、どう計算してるか分かると思います。もし、ランダムフォレストそのものの部分が分からない場合は、ランダムフォレストのつくりかた(C++の実装例つき) - じじいのプログラミングを参照してください。

アルゴリズム

入力変数X, 目標変数をyとします。データサンプルは(X1,y1),(X2,y2),...(Xi,yi),...(Xn,yn)のn個あるとします。
またXはベクトルでXi = {xi1, xi2, xi3, ..., xim}とm個の特徴量を持つとします。

以下の処理は、ランダムフォレストの決定木1本1本に対して行います。

  1. ランダムフォレストで、決定木を作るときに訓練に使われなかったデータサンプル(ブートストラップサンプリングで外れたデータ)を使います。これをOOBデータと呼びます(Out-of-bag(OOB) data、OOB examplesとも)。
  2. まずOOBエラーを求めます。OOBデータのXを決定木に入力すると、予測した結果y'を求められます。y'と目標変数yがどれだけズレてるかがOOBエラーです。回帰のときは、MSE(mean square error、平均二乗誤差)をOOBエラーとします。つまりΣ(yi-yi')^2 / kとなります(kはOOBデータのデータ個数)。
  3. m個の特徴量それぞれに対して重要度を求めます。単純にm回ループして、ループ変数を特徴量fとして、以下のことをします。
    1. OOBデータの特徴量fだけをシャッフルします。他の特徴量とyはシャッフルしないので、気を付けてください。
      • X1=( 1.0, 2.1, 3.2) y1=1.2
      • X2=( 0.5, 1.6,-2.1) y2=2.1
      • X3=(-1.2, 1.5, 2.0) y3=3.5
      • ↓シャッフル
      • X1=( 1.0, 1.6, 3.2) y1=1.2
      • X2=( 0.5, 1.5,-2.1) y2=2.1
      • X3=(-1.2, 2.1, 2.0) y3=3.5
    2. シャッフルしたデータを使って、(2)と同じようにOOBエラーを求めます。特徴量fだけシャッフルしたので、特徴量fが重要なのであれば、重要なデータがシャッフルでおかしくなったわけなので予測が悪くなるはずです。予測が悪化した値、つまり、(シャッフルしたデータのOOBエラー) - (元のOOBエラー) を重要度IncMSEとします。

最後に、全ての決定木のIncMSEの総和を求めて、最後に木の本数で割って平均を求めます。

ポイント

  • 計算時間は、かかります。訓練の計算は1回だけで済みますが、予測は、(OOBデータ個数)*(木の本数)*(特徴量の数)*(特徴量1つあたりのシャッフル予測の回数 nPerm)回だけ行うことになります。予測は訓練に比べたら計算時間はずっと短いですが、回数が増えればもちろん時間はかかります。時間がかかりすぎの場合は、上記の4つのどれかを減らす必要があります。
  • 重要度として、それなりに信頼できます。シャッフルする部分以外は、通常のOOBデータを使った予測と変わらないので。
  • マイナスの値もとりえます。予測に影響をほぼ及ぼさない特徴量だと、シャッフルして、たまたま逆に予測精度があがることもあり得ます。もちろん、重要度は最悪という意味です。
  • Rでの変数名は、%IncMSEと%がついていますが、全部足して100%になるわけではありません。

R Projectによる、C言語での実装例

本体はC言語で書かれています。CRAN - Package randomForestからrandomForest_4.6-10.tar.gzをダウンロードしましょう。 src/regrf.c の262-295行と327-338行のあたり、if(varimp){} のスコープ内のコードを読むと分かりやすいと思います。

IncNodePurity

C++での実装例

先程のコードの#if INCNODEPURITYで囲まれている部分で、求めています。

アルゴリズム

ランダムフォレストのノードの分け方を知ってるなら、何も新しい知識はいりません。
決定木を作る際に、各ノードを子2つに分けるベストな分け方は、(親のMSE)-(子2つのMSE)が最大をとるときです((親のMSE)は固定なので、(子2つのMSE)の最小化でもいいです)。数式で書けば以下のようになります。

Aを集合として、i番目の目的変数をYiとして、親の集合Pを、子の集合Lと集合Rに2つに分割する。
要素数 |A|
総和 S[A] = {Σ[i∈A]Yi}
平均 E[A] = S[A]/|A|
と表す。MSE(平均2乗誤差)の変化量は以下のようになる。
  {Σ[i∈P](Yi-E[P])^2} - {Σ[i∈L](Yi-E[L])^2} - {Σ[i∈R](Yi-E[R])^2}

これを選ばれた特徴量別に足していくだけです。最後に木の数で割ります。

ポイント

  • 計算時間は速いです。決定木を作る際に求まるので、訓練をすればいいだけで、予測の時間はかかりません。これはIncNodePurityの最大のメリットです。
  • 重要度としては、あまり信頼できません。確かに分割の際に使われる特徴量は重要とは言えますが、あまり直接的ではないので。ただ、例えば、特徴量が1000個あるときに下位50個を削りたいと言った、ざっくりとした特徴量抽出には使いやすいです。 (実際、1位を取ったマッチToxCast Challenge | EPA TopCoderでは、これを使いました。)
  • 決定木の数が増えても、値は増えません(木の数で割って平均をとってるので)。
  • サンプルデータの数に比例して、値は大きくなります。
  • 経験上、決定木の数が少ない場合でも、かなり値は安定することが多いです。最初は多くの決定木を使って、次に少ない決定木を使って、結果があまり変わらなかったら、それ以降は少ない本数で試すというアプローチが、時間を節約にもなって良いと思います。

R Projectによる、C言語での実装例

先ほどダウンロードした、randomForest_4.6-10.tar.gzを参照してください。

  • src/regrf.c の変数tginiがIncNodePurityになります。
  • src/regTree.cで、tginiを計算していますが、変数crit→critvar→critmax→decsplit→tginiというふうにデータが渡されています。247行目 crit = (suml * suml / npopl) + (sumr * sumr / npopr) - critParent; が重要です。平均2乗誤差を使う場合は式展開すると、以下のように簡単になり、これで合ってます。

  {Σ[i∈P](Yi-E[P])^2}                          - {Σ[i∈L](Yi-E[L])^2}                          - {Σ[i∈R](Yi-E[R])^2}
= {Σ[i∈P]Yi^2} - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2 - {Σ[i∈L]Yi^2} + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2 - {Σ[i∈R]Yi^2} + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... 2乗を展開した
=                - 2E[P]{Σ[i∈P]Yi} + |P|E[P]^2                  + 2E[L]{Σ[i∈L]Yi} - |L|E[L]^2                  + 2E[R]{Σ[i∈R]Yi} - |R|E[R]^2  ... {Σ[i∈P]Yi^2} = {Σ[i∈L]Yi^2} + {Σ[i∈R]Yi^2} で消した。
=                - 2E[P]S[P]         + |P|E[P]^2                  + 2E[L]S[L]         - |L|E[L]^2                  + 2E[R]S[R]         - |R|E[R]^2 ... 総和Sであらわした
=                - 2S[P]^2/|P|       + S[P]^2/|P|                 + 2S[L]^2/|L|       - S[L]^2/|L|                 + 2S[R]^2/|R|       - S[R]^2/|R| ... E[A] = S[A]/|A|で平均Eを消した
=                -  S[P]^2/|P|                                    +  S[L]^2/|L|                                    +  S[R]^2/|R|                    ... 足し算

重要度をどう使うか?

特徴量選択

無意味な特徴量を削ることで、予測精度を上げることができます。また、計算時間も速くなり、多くの木を作れるので、それでも精度を上げることができます。

以下に、自分が行う特徴量選択の方法をまとめてみました。今回の2つは(3)(4)にあたります。

手法 選択性能 計算時間
(1) 特徴量を削りクロスバリデーションで予測するのを、全ての特徴量で試す 良い*2 とても遅い
(2) 特徴量を削りOOBデータを使って予測するのを、全ての特徴量で試す 良い 遅い
(3) IncMSEで重要度の低い特徴量を削ってみる 普通 普通
(4) IncNodePurityで重要度の低い特徴量を削ってみる 悪い 短い
  • かけられる時間やデータ量にもよるのですが、(4)→(3)→(2)→(1)の順に試すのが良いと、自分は思います。基本は(2)になると思います
  • (1)(2)(3)については、MSEではなく、精度の評価関数(scoring function)に置き換えることもできます。例えば、yの大きい順を予測したい順位当ての問題のときは、順位差の絶対値の総和をMSEの変わりに使うといったこともできるでしょう。
  • ただ「重要度を低い特徴量を削る」と簡単に書きましたが、貪欲に1番重要度が低い特徴量を削っていけば良いというものでもなく、とても難しいです*3。特徴量間に強い関係があれば、特徴例をけずるたびに、重要度も大きく変わってきます。良い方法があったら(というか研究されていると思いますが)、ぜひ教えてもらえるとうれしいです。一応、自分のアイディアレベルのネタとしては、
    • (2)(3)であれば、複数の特徴量をシャッフルする手はあるでしょう。ただ、(OOBデータ個数)*(木の本数)*(特徴量の数 ^ シャッフル個数)*(特徴量1つあたりのシャッフル予測の回数 nPerm)回となるので、特徴量が多いときは、シャッフル個数2個にするのすら遅すぎて難しいかもしれません。
    • 探索枝仮り・焼きなまし法・ビームサーチといった最適化手法で削っていく方法もあるでしょう。
    • 特徴量を1つ削ったとき、他の重要度がどう変化するかで、どの特徴量間に強い関係があるかというグラフのようなものを作れば、最適化手法を使っていく上で役に立つかもしれません。例えば、「特徴量はx1,x2は片方を削ると、もう片方の重要度が大きく下がるので、たぶんペアで意味がある特徴量だ。焼きなまし法を試すときは、この2つはまとめて追加したり削除したりするようにしよう。」といったことはできそうな気がします。

ランダムフォレストのパラメータ調整がひどいのに気付くきっかけ

最初にあげた図では、SUBJIDという特徴量は、重要度IncMSEで最小、重要度IncNodePurityで最大になっています(SUBJIDは、本当に無意味なIDで、実際はIncMSEが正しかったです)。
このように、IncMSEの大きさの順位とIncNodePurityの大きさの順位が異なる場合、経験上ですが、予測精度も悪くなる傾向にあります()。恐らく原因としては、

  • 元の特徴量が少なすぎる
  • 木を深いところまで作りすぎている(nodeSizeが大きい or maxLevelが大きい)

のあたりでしょう。
ideoneのコードの645行目

	rf->train(trainingFeatures, trainingAnswers, 5, 999, 2, 5, 0.66f, numTrees);

のnodeSizeとmaxLevelのパラメータを変化させ、IncMSE・IncNodePurity・totalErrorの変化をみてみると良いと思います。


以上です。ご感想・ご意見お待ちしてます。

*1:今回は詳細は述べませんが、詳しく知りたい方は、RのソースコードをimportanceSD, impSDで検索したり、importance.Rからソースを追ってみてください。SDという名前で標準偏差に見えますが、なぜか標準誤差なので注意してください(ネットでも多数の間違いが見られます…)

*2:「とても良い」の可能性もありますが、OOBデータ vs クロスバリデーションの議論は奥深いようなので踏み込みません

*3:あくまで自分のイメージですが、「特徴量の追加削除」を近傍、「重要度」を関数の傾きと置き換えると、近傍にいっただけで関数の傾きが大きく変化したりする関数を最適化しないといけないようなものなので、やばいかんじしかしません…

ランダムフォレストと他の機械学習(or統計)を組み合わせて使う

もしかしたら、プロにとっては当たり前のテクニックかもしれませんが、自分は初めて見たので書きたいと思います。また、おそらく大きい効果を出すのが難しいテクニックだと思われるので、まずは基本的なことを先にやったあとに試したほうがいいでしょう。

追記 3/13 このテクニックは「stacking」と呼ばれるそうです。Twitterで教えていただきました。ありがとうございます!
Stacking, Blending and Stacked Generalization | Garbled Notes
Ensemble learning - Wikipedia, the free encyclopedia



ここから本題です。

前々回のTopCoder機械学習マラソンマッチTripSafetyは、飛行機のフライトデータから、違反フライトを予測するという問題でした。問題やデータはこちら

f:id:shindannin:20150313095715p:plain

予測精度をよくするために特徴量を追加したりするのですが、自分を含めて多くの人がやった手法は、特徴量(x24,x25,x26,x27,x28)だけを使って、新たな特徴量(x29)を追加するものでした。x24-x28は重混雑・中混雑・軽混雑と似たような指標なので、重み付けをして1つの特徴量にまとめるという発想は自然にでてくると思います。自分はx29 = 32*x25 + 8*x26 + 2*x27 + x28としました。
f:id:shindannin:20150313095726p:plain

しかし、優勝したPsyhoさんは、元の特徴量(x24,x25,x26,x27,x28)と目的変数yに対してを使って、ここで正則化つきの線形回帰を使って、その結果を新たな特徴量(x29)として追加しています。新たな特徴量を追加するのにxだけではなくyも使っているのが最大のポイントです。*1
f:id:shindannin:20150313095739p:plain

全体の流れは以下のようになります。

  1. 上記のように、訓練データのx24~x28とyを入力して、線形回帰して、x29=α24x24+...+α28x28+βの、α,βを求める。*2
  2. 訓練データのx1~x29(x29はx24~x28,α24~α28,βから求める)とyを入力して、ランダムフォレストで学習させる
  3. テストデータのx1~x29(同じくx29はx24~x28,α24~α28,βから求める)を学習させたランダムフォレストに入力して、yを予測する

線形回帰での予測結果x29を、y=x29としてそのまま答えにすることもできるのですが、線形回帰とランダムフォレストの2段構えになっています。

なぜ、これがうまいのかというと、

  • 線形回帰の予想精度が良ければ、その結果をほぼそのまま返すことができる。
  • 線形回帰の予想精度が悪かったとしても、ムダな特徴量が1つランダムフォレストに追加されただけなので、全体としての予想精度の低下はほとんどない。(ランダムフォレストでの寄与度がとても低い変数が増えただけ)
  • 線形回帰自体は表現力が高くないけど、その弱点をランダムフォレストでカバーできる。今回の場合でいうと、x29とyの値は近くないけど(特にy=1のとき)、この後にまだランダムフォレストがあるので、そこで精度を上げられる可能性がある。

参考までに、実際に100ケースで、自分のPCでテストしたところ、得点は285363点→288294点と約3000点増えています。あまり大きくはないですが、本番で2位だった112atn323さんとの点差も約3000点差と僅差だったなので、他の部分も十分詰めたあとなら試す価値はあると思います。*3

今回の問題では線形回帰でしたが、簡単にいえば機械学習(or統計)の手法の多くは、2つ組み合わせることができるということです。Psyhoさんのコードを見る限りでは、ロジスティック回帰+ランダムフォレストも試そうとしていたようです。

もちろん、複数の機械学習手法を試す人は多くいると思いますが、以下のどちらかが多いです。

  • 複数の手法を試して、その中で一番良かった手法を使う
  • 問題をいくつかのケースに分けて、ケースごとに一番良いものを使う(例:データが大きいケースはランダムフォレスト、小さいケースはニューラルネット)

でも、今回のPsyhoさんの解のように、2つの機械学習をあわせて使うのは、自分は初めて見ました。

ただ、まだ分かっていないのが、「どの特徴量に対して使うのが効果的か?」。これは難しいです。今回の問題だと、混雑度が大きければ大きいほど、違反フライトは単調増加しそうなので、線形回帰のように表現力低めの手法でも合いそうな気がしますが…。今のところは色々試す、そしてそれを自動化するぐらいしかないですね…。今後の課題です。

*1:yを使わずxだけを使って、線形回帰してxを追加するというのは、以前の記事で少しだけ述べたように普通に使われるテクニックです

*2:実際にはx1~x23を使って、さらにいろんな特徴量を追加してますが、ここでは説明を簡単にするために省きました。

*3:ただ、得点の評価については、本当はもっともっと試さないと正しく評価できないので、これより効果がもっと大きい可能性も小さい可能性もあります。