Hokkaido Univ.& Hitachi 1st New-concept Computing Contest 2017の反省

AtCoderで初の本格マラソンマッチ。楽しかったです。13位/302。最近のマラソンマッチ不振を考えれば良い結果だったけど、反省点も多かったです…。
（方針の公開がOKになったので、公開します。2ndがあるので、ちょっと迷いましたが、北大日立さんにとって良い結果が出たほうが、今後さらなるマラソン開催にもつながるので、公開します。）

問題

- Hokkaido Univ.& Hitachi 1st New-concept Computing Contest 2017 | AtCoder

反省会生放送

TopCoderでプログラムしてみた 第2365回(北大日立マラソン第１回反省会生放送）

ファーストインプレッション(2017/11/16)

• どういうデータなのか？
  • 辺の数とか、平均でいくらぐらいなんでしょう？ →辺は思ってたより多い。100本以上辺があるのは普通。

• 最適解がわかるケースを逆につくれば便利。最適解を知っていれば、焼きなましをやってみたときに、うまくいってるかわかる。
  • King's graphから作る
  • 一定条件を満たせば、そこから辺を足しても最適解をキープすることは可能？より一般的＆難しい問題が作れる
    • 三角不等式を満たすものは除去してもいい？
  • このときに意図的に最適解をつくる操作が、逆の操作が近傍として有力な可能性はある。

• より簡単なケースで考えてみる。
  • King's Graphではなくグリッドグラフにマッピングする場合は？二部グラフに持ち込める？
    • 最大流かカット系で楽にできないか？
    • 使えなくても、初期解として有効な可能性があるので、無駄にはならないと思う。

• 過去の反省
  • 類問の解法に固執しすぎ。
    • 今回のは、TCO17 MM Round 1 TCO17 Marathon Round 1 (GraphDrawing) の反省 - じじいのプログラミングに似てる気もする。そうすると優勝したchokudaiさん解の重力モデルか焼きなましの中間っぽいのが有効なのか？
    • ただ、ちょっと条件がちがっただけで、解法が全然違うことがある。まっさらでスタートしたほうがいいのでは。

• 方針
  • ベタに行けばこの問題は、ただの焼きなまし。そして、最後に焼きなましの解をダメ押しでGreedyするの忘れずに。
  • 粗いグリッドで焼きなましから、徐々に細かいグリッドにして焼きなまし。
    • ただし、先日tomerunさんと似た解法だったのにもかかわらず、自分だけ大敗したので、ちゃんと敗因を確認すること。特に、細かいグリッドになったときに調整できるかどうかもポイント。

• ランダムグラフの頂点間距離は使える？
• 高得点辺だけを使ったのGreedy・焼きなまし
  • 1点や0点の辺は無視してもいいくらい。というか0点は最初から削っていいのでは。

• 何を焼きなまし？（「近傍の種類の多さが重要」と昔Komakiさんが言ってた）
  • 頂点ベース
  • 頂点ベース（位置を気にする。近傍は近くの点）
  • 辺ベース
  • 15点→0点で、辺をじょじょに開放
  • 場所を後から、配置する解法はないか？よく、「置く場所を決める焼きなまし」→「何を置くか決める焼きなまし」のように、２段階に分けると良いこともある。

Submission 1 (2017/11/18) 52935点

同じ番号の頂点をマッピング。同点のみなさん多数。

その後、多忙で無事乗り越えたものの、休日に冷房病になってしまい、計１週間休み。すごく警戒してたのに…。

Submission 2 (2017/11/25) 129838点

場所2箇所の交換（頂点あるなしは気にしない）を100万回する山登り法。

Submission 3 (2017/11/25) 148673点

場所2箇所の交換（頂点あるなしは気にしない）を、9秒間する焼きなまし法。

Submission 4 (2017/11/26) 152871点

焼きなまし法の、温度を5度→1度にする。わりと効果があった。

ビジュアライザ

このタイミングでビジュアライザも追加する。 Siv3Dを今回も使いました。感謝→ Siv3D · GitHub）ビジュアライザは、

• 押すボタンごとに、いろんな近傍への遷移をリアルタイムで試せる。
• 頂点ごとに最高点と現在の得点を追加。ちなみに、最高点はベスト8の総和。

ただ、ビジュアライザをつくっても、なぜ問題が起きてるのかにちゃんと着目した情報を表示して、本質を見抜けないと意味がないので、気を付ける。
このchokudaiさんの記事はおすすめ（というか、悪い例にでてくるモデルは、たぶんわしじゃ…。）
chokudai.hatenablog.com

焼きなまし法の温度と得点、統計情報（どの近傍で、どれだけ上がったか）等も、グラフプロットに追加すればよかった。

Submission 5 (2017/11/26) 159335点

動かす頂点Aと、元グラフでつながっている頂点Bの隣に移動。上位の人が多く使っていた有効な近傍。

セカンドインプレッション(2017/12/26)

• 焼きなまし温度調整必須
• 15点の辺だけ使って、辺を少なくして、橋・関節点で分けて、
• 浮いている点の15とかは重要？
• 15,15,15,15,15,15,15,15,15,15,15,15,15,15,5,5,5 みたいなときは代替があるので、あとで決めていい
• 塊ごと移動できない？
• 番兵で、範囲チェックのifを減らせ
• 範囲広すぎだと、ダメかもしれない。60x60ではなく、範囲を狭くしたら。
• 範囲なしにして、あとからくっつけたほうがいいかも。

Submission 6 (2017/11/27) 159913点

• 元グラフでつながっている頂点Bは、高得点ベスト８の隣接頂点から選ぶ。かつ、焼きなまし法を2倍高速化

• この時点で、焼きなまし法のコツ Ver. 1.2 - じじいのプログラミングに載っている、

★★★★ 高速化がどれだけ効果的か検討する。
のも重要です。ためしに10秒だけやってる計算を、あえて100秒とかに伸ばしてみて、結果がどうなるか見てみましょう。

• 時間をかければかけるほどスコアが伸びる
  • 高速化しよう＆スコアが素早く伸びるように次の状態を決めよう。山登り法に変えるのも良い。
• 時間をかけてもスコアが伸びない & スコアが高い
  • 最適解がでてる可能性が高い。計算を打ち切って、余った時間を他に使う
• 時間をかけてもスコアが伸びない & スコアが低い
  • 次の状態の決め方、スコア関数の見直しを再検討する。焼きなまし法を使わない法がよいときもあるのも忘れずに。

を試した。なんと計算時間を10倍にすると、30ケースあたり平均+5000増える！！（実際には30ケース*10セット）1位が狙える位置に行くことがわかる。

でも、10倍高速化するのは無理だろうなという考えと、まだ試してないアイディアがいっぱいあるので、まずそちらからやることにした。これが敗因の１つ…。
10倍で効果があるのだから、せめて、計算時間を5倍にしたら、計算時間を2倍にしたら、どれだけ増えるのか試して、それをみて高速化にどれだけ作業するか考えるか判断すれば良かった。

失敗したアイディア集

ここから点数がさっぱり増えなかった、ボツアイディア集。反省点の１つとして、結果の確認とまとめを、もうちょっとちゃんとしてれば良かった。失敗したアイディアの中では、点がほとんど変わらなかったものもあるけど、これは下手すると最終盤でちょっと点数を上げたいときに有効な場合もあるので。

初期値

初期値がわるいとダメなケースもある。特に、今回のは、わりとバラバラに頂点が飛んで非効率に見えたので、初期位置を中央に固めた。でも、すぐに結局バラバラに飛んで行って、スコアが上がらなかった。

バラバラにならないようにする。

中央に行くほどボーナス点をつけたり、KingGraphの隣に頂点があるときにボーナス点をつけた。バラバラになりにくくなったけど、スコアはだいぶ下がった。バラバラシャッフルっぽい動きは有効だったっぽい。

頂点Aを意図的に選ぶ。

頂点得点の最大値（頂点の辺ベスト8の合計）をあらかじめ計算しておき、交換する頂点を、以下のような基準で選ぶけど、これもだいぶ点が下がった。

• 低得点率（現在の得点 / 最大値）
• 得点差（最大値 - 現在の得点 ）

yowaさんのなんか他の人の解法をみるに、全てを均等回数を選んだほうがいいという議論もあるらしい。

同じく最初にシャッフルしたあと順番に舐めてた。たしかTopCoderマラソンの↓この回(KnightsAttacks)の感想戦で「焼きなましで変更加える箇所、ランダムじゃなくて順番でいいんじゃね？」みたいな話が出てたのを思い出した

| TCO17 Lighting MM - Togetter - https://t.co/jAG3dXsMvd

— yowa (@yowa) 2017年11月30日

次回は要検討。

隣接頂点Bの選び方

Submission 6「元グラフでつながっている頂点Bは、高得点ベスト８の隣接頂点から選ぶ」の部分が、あまりに雑だったので、ここはスコアは伸びるかなと思ったけど、伸びなかった。辺得点x,y,zだったら、選択確率は、pow(x,t)/(pow(x,t)+pow(y,t)+pow(z,t))のようにして、tを調整した。tが0なら全部等確率1/3, tが∞なら最高辺得点のだけが確率1になる。でも、スコアは伸びなかった。

さまざまな近傍

「近傍をたくさん作ったほうが勝つ」と、Komakiさんが昔言ってたので、いろいろ足した。

• 2-Swap
  • 隣り合う2点を交換。A→B, B→A
• 3-Rotate
  • 三角の形で隣接する3点を、A→B, B→C, C→Aのように移動。ある頂点Aを固定すると、三角形の選び方は12パターンある。
• 4-Rotate
  • 四角の形で隣接する4点を、A→B, B→C, C→D, D→Aのように移動。
• スライド
  • 縦 or 横 or 斜めに一列にならんだx点（2点～4点）を一方向に押してスライドさせる。A,B,C,Dと並んだ場所で、A,B,Cが頂点があって、Dが空きのとき、C→D, B→C, A→Bと移動させる。

それぞれに対して、全てのマスが頂点埋まっているバージョンと、そうでないバージョンを試した。また、近傍については確率で選んで調整する（辺隣接ワープ 95%, 2-swap 4%, 3-Rotate 1%とか、そんなかんじ）

• 空きマスが一部ある場合の結果は、概してよくなかった。
• 2swapについては、低得点の際には効果があるように見えたが、マッチの終盤にはほぼ意味がなくなったので、最終的には、これらの近傍は全てボツとなった。

また、本当なら、別な動きをする近傍をいろいろ作ったほうがいい。今回の3-Rotate, 4-Rotateは、2-swapで2回,3回で代用可能かつそれなりの確率で起こりうるので、あまりよくない近傍（ただ、それでも1回は試すべき。）

レプリカ交換法

ほぼやけくそ。N個の焼きなまし法を同時に走らせる方法。それぞれ別な温度を固定でもつ。一度やってみたかったので、やった。
レプリカ交換法 - Wikipedia
得点差に応じて、温度を入れ替える（状態を入れ替えるより簡単。）
これは、あまり有効でない。中途半端な得点の解がたくさんできるだけだった。特に、今回のように時間をかければスコアが伸びるケースでは、時間が1/Nになるのはもったいない。１つに集中して時間をかけたほうがいい。ただ、並列計算とかできる環境であれば、有効なケースもあるのかも。

プロファイラで高速化。

この時点で、最終日で、残り4時間ぐらいになってしまった…。Visual Studioの標準機能のプロファイラを使って、高速化する。しかし、時間が足りず、結局、以下の２つは高速化せずに終わってしまった。

• 番兵で範囲チェックをなくす
  • ファーストインプレッションにも書いたのに。時間がなくなってしまった。番兵など、最終盤で、高速化するにはバグを埋め込むリスクが高いところは、もっと早めに高速化やったほうがいい。
• 焼きなましのexpの式
  • あわててテーラー展開したやつを書いたけど、点が下がったので、怖くなったので、やめた
  • 焼きなましのexpの式なんて毎回使うのに、高速化を最終盤にしようとして、おざなりになるのを繰り返してるのは、本当にダメである…。効果が少なくても、1度作れば一生使えるのは、暇なときにしっかり作りましょう。

自分のは、7000万回程度しか回らなじゃった。他の人は2億～4億ぐらい回っている。その中でもっとも怪しいのが、ケチらずintとかdoubleとか大きい方を使っていたこと。一般に、プロファイラはボトルネックにはすぐにきづくが、全体的に遅くなってる要素には気づきにくいので、要注意。ただ小さい型を使うのが、これだけ早くなるのは意外。キャッシュミスあたりが予想されるけど、今回のフィールドは小さいので、過小評価してた。自分がなんか誤解している可能性も高い。要調査。

その他

• お尻に火がついたのが、いつもより少し早かった。
• 若くないので、体調管理は大事。冷房病は徹底した温度管理で避けよう。
• 昔は、夜遅くなっても次の日なんともなかったけど、今はリバウンドがくるので、ちゃんと早めに寝たほうがいい。
• AtCoderさん、ありがとう！マラソンマッチでこんな問題を受注できるとは、本当にすばらしい。

TCO17 Marathon Round 1 (GraphDrawing) の反省

問題と1位の方（chokudaiさん）の解法

TCO2017R1

結果

• 力学モデル（バネ）で、長さ比に応じて力を加える。その後、貪欲法（一番悪い辺の点を、最も高得点になる場所へ移動）で解を改善。
• 45位/146位 652,515.28点。30位ぐらいの人ともかなりの点差が開いていて、敗北感ある…。

敗北までの流れ

ファーストインプレッション

• 方針
  • 力学モデル
  • 焼きなまし
  • よさげな頂点同士をマージしながら構築

• 問題の穴
  • 木の足の部分はどうでもよい。
  • 距離が近くになる点どうしは、まとめて動かしてもいいかも。
  • 三角形などサイクルの形で、辺が絶対につながらず満点が取れないケースもあるので、こういうケースでは理論上最高得点を求めて、最高得点にあわせて、目標長さの条件も緩めたほうがいいのでは。
  • MAX_RATIO / MIN_RATIOなのに注意。全部言われた長さにせよとは言ってない。全部スケールしても得点は変わらない？
  • 三角形って何個あるんでしょう。位置関係が確定する気がする。
  • 長さが整数値になることを使ったHackはなさそう。

ざっくり焼きなまし

• 近傍
  • ランダム移動、だんだん移動量を少なく
• 評価値
  • 最小値だと、差分を使いづらく、雑に書くとO(E)、良くてO(logE)ぐらいな気がしたので、とりあえず重み付きの和で対応

Submission 2 得点 248563.14
上位の得点が、70万点近くなので、恐らくハズレ方針だと思い、あきらめる。

力学モデル（バネ）

メモ

• 微分せよ。
• ズレ0のやつで、位置があうか比較する。
• インラインアセンブラがききやすい
• 論文（良い初期位置とか）
• Logをつけてみては？
• 引っ張る力を比ベース
• たぶんランダムな力とかは筋がわるそう。

試したこと

• 普通のバネで引っ張る
  • バネモデル、意外とバネ定数により、発散しやすい。
  • 枠外でも力を加えて、枠に入れたほうが高得点になる。
  • 今回は目標の長さとの比が得点にかかわるので、長さ比に比例した力を加える（普通のバネは差に比例した力を加えるので、）
  • バネの力により簡単に発散するので、バネ定数の調整と、pow()をつけた調整を試す。
    • powやlogを付けた調整は発散しやすかった。
  • 距離の短い辺が、強い力でずっと動き続けて足をひっぱる。
• 序盤は遊びを入れて、自由に動けるようにする。目標得点を0からスタートし、目標票得点内に収まるのであれば力は加えない。そして徐々に目標得点を増やしていく（拘束も強まる）
  • →これは極値にハマる可能性を大幅に減らせる有力な方針だと思ったけど、そうでもなかった。結局同じような場所でハマる。
• エネルギー最小化
  • https://cs.brown.edu/~rt/gdhandbook/chapters/force-directed.pdf のKamada and Kawaiのアルゴリズム。
  • バネエネルギーは積分すればいいので、

差ベース
// energy = integrate k*(x-a) dx from x=a to b ====> k * (a-b)^2 / 2

比ベース
// energy = integrate k*(x/a-1) dx from x=a to b ====> k * (a-b)^2 / (2*a)
// energy = integrate k*(a/x-1) dx from x=b to a ====> k * (a * log(a/b) + (b-a) )

• • 本当だったら、高速化のためKamada and Kawai多変数のニュートン法で = 0系の連立方程式を解くのだけど、まずは遅くてもいいので、エネルギー最小になる点を全探索で試すことにした。バネより結果が悪いのであきらめる。
• このタイミングでビジュアライザも作成。正しくは動いているっぽい。
• ズレなしで満点を取れるテストケースを試す
  • それでも満点がとれない。
  • どうやら初期位置が悪いと簡単なケースでは、ダメなことが分かる。
  • 逆にチートして、初期位置を知っているものとしてやると、バネでも高得点が取れる。
• 簡単なわりには低得点になるケースを探す
  • 変な三角形とか、田んぼの田とか、いろいろと手作業でいろいろと作るも見つからず。

Submission 3 得点 475556.16

恐らくバネ方針も良くないのだと思い、あきらめる。

強引な貪欲法

もう一番得点が低い辺の頂点の片方を、700*700の移動先を全部試して一番良いところに移動させる。というか、なぜこれを最初に試さない。

• 高速化する。得点ハイトマップを作成して、関数形を見る限り、そんなに山の数が多いかんじではないので、焼きなましで探していいのでは？（ドーナツを重ねてくり）

• • 250回ぐらいの焼きなましループで、最適解（=700*700）を90%ぐらいの正解率
• 得点は伸びたが、上位の70万点には程遠い。
• 貪欲→バネ→貪欲→バネを繰り返せば、極大値にハマらないのではないか？
  • ハマる

Submission 4 得点 606403.67

焼きなまし

貪欲でいいのなら、焼きなましもやっぱり行けるでしょうと思い、焼きなましの戻る

メモ

評価値

• 最小値（問題のスコアと一緒）
• 足し算
• 重み付き足し算（ペナルティはV倍とかにしないとダメ）
• 分散

近傍

• 頂点どうしの位置交換
• 悪い頂点のとなり
• 辺対称
• 大ランダム移動（たぶん大きい移動が入った時点で負け）・小近傍移動
• 速度つき近傍（慣性）
• 粗いグリッド⇒細かいグリッド

実際にやったこと

評価値

• 最小値
• 足し算
• 重み付き足し算（に比が1から離れていれば離れるほど、さらに悪くなる）
• 分散
  • O(1)で計算できるし、２番目に悪い点、３番目に悪い点も実質考慮できるので、よいかなと思ったけど、だめだった。

近傍

• 粗いグリッドに分けて、ざっくりした位置を決めたのち、徐々に細かいグリッドにしていく。2*2→4*4→8*8→…
  • 32*32を最高にスコアが落ち始めるので、あきらめる。←これはもったいなかった。32以降は別方針にするとか、2段階の解法で高得点いけた

Submission 2と同レベルの得点しかとれないので、やめる。

Beam search系

迷走する。書くほどでもないか…

反省点

部分最適（極値にハマる）理由を考察するのを第一に考える。

シンプルで統一的な解法にこだわりがち。

• シンプルで統一的な解法でないと高得点が取れないと、その後複雑な解法にした際も、最終的には点が伸ばせないと考える癖がある気がする。
• それ自体はいいのだけど、どうもシンプルな解法をいろいろ試すものの、それぞれの解法への考察が雑になり、点が伸ばせる解法を、謝って枝刈りしてしまっている気がする。
• 悪いところを考察したうえであれば、例外的な対処も入れて良い。部分的対応をどんどん入れていってよい。
  • 力学解法でも、焼きなましでも、高得点は取れてたっぽい。
  • 悪そうな解法でも、部分的にでも良いところを探すようにする。
  • ちゃんとした物理に従う必要はないですよ。
    • バネモデルで、たまにワープさせたって構わない
    • バネモデルで、全部のやつに同時に力を加えなくたって構わない。（一番悪いやつだけとか）

評価関数を状況に応じて変える

（分かってるつもりなのですが…）

• 例えば、焼きなましのグリッド解法は、高得点を取ってる人はいるので、良い方針だったはずなのだけど、2*2のときも、700*700のときも、似たような評価関数を使っているのは、大変良くない。
• そもそも、段階ごとの目的に意識がいってない気がする。粗いグリッドのときは、「正しいセルに振り分ける」のが目的、細かいグリッドのときは「高得点取れるようにする」なのに、評価関数が1つというのは良くない。それぞれの段階での目的をはっきりさせて、評価関数をつくろう。
• これは焼きなましに限らない。Beam Searchの前半と後半で評価値を変えたりもするし、目的ベースで、柔軟にかえる

頭の悪さは、ビジュアライザとランダムテストで補おう

• 極値にハマる、単純なケースを、自分でテストケースを作るものの、見つけられなかった。
  • ランダムで単純なテストケースを大量に作って、試せばいいだけ
  • せっかくビジュアライザもあるんだから、考察も簡単にできるのだし

スコアに大きな影響を与えるものを、必ずしもアルゴリズムの最初に確定する必要はない。

• 今回は短い辺が、低得点の原因になるけど、その場所を先に確定させると高得点にはならない。
• 短い辺は確かにボトルネックになるけど、アルゴリズムの後半でGreedyとかやれば必ず治せる。
• 似たようなテクニックとして、「後で決めても辻褄が合っていればよい。(Peg Jumpingとか)」というのもある。

このブログは、マラソン開催中に非公開で書くべし

• 絶対にマラソン中に書いたほうが効率がよい。
• 他の人の解答をまだあまり検討していない⇒復習としては最悪！

ビームスタックサーチについての勉強メモ

マラソンマッチ界隈でよく知られている通称chokudaiサーチ、それと似ていると言われてるビームスタックサーチ(beam-stack search)について、ちょっと調べてみました。論文はこちらです。
R Zhou, EA Hansen (2005) Beam-Stack Search: Integrating Backtracking with Beam Search, 15th International Conference on Automated Planning and Scheduling

でも、いまだ分からないところも多いので、まとめました。誤解してそうなので、間違いを指摘してもらえると助かります！論文と照らし合わせながら読んでもらえると良いと思います。

この記事を読むのには以下の知識が必要です。

• ビームサーチ
• 全探索（深さ優先探索など）
• A*アルゴリズム等で出てくる、ヒューリスティック関数の意味

ビームスタックサーチの概要

前提

• 許容的なヒューリスティック関数、つまり、最小値を求めるときは、現在の状態～未来で取りうるコストの下限が必要。最大値を求めるときは上限。（Admissible heuristics)
• 最適解を求めるのが目的。枝刈りは行われるが、基本は全部のパスを通ろうとするので、最適解があれば、必ず見つけられる。(Completeness)

特殊なケース

• ビーム幅が十分に大きいときは、分枝限定法-幅優先探索と同じ。バックトラッキングは行われない。
• ビーム幅が1のときは、分枝限定法-深さ優先探索と同じ。（ちなみに、ビームサーチだと、ビーム幅1は貪欲法と同じ。）

疑似コード

上記の論文に載っているビームスタックサーチの疑似コードの引用である。24～28行目の変数relayの部分は、拡張バージョンdivide and conquer beam-stack searchの部分なので、今回は無視する。

ビームスタックサーチの流れ

ビーム幅2で、ビームスタックサーチを使って、最小値を求める場合を考える。（注意：まだ不明な点があるので、グラフをわざと簡単な木にしてます）

各ノードnを状態とし、それぞれのノード内の数字は評価値を表す（f-costと呼ぶ。ノード番号では無い！）。評価関数f(n) は
f(n) = g(n) + h(n)

• g(n)は、スタートノードからノードnまでのコスト
• h(n)は、許容的なヒューリスティック関数、ノードnからゴールノードまでのコストの下限（絶対に過大評価しない値）

スタートノード（f-costが0）から、最小値を求めていく。

Layerごとに、

• f-costでソートする。
• fmin・fmaxという範囲を示すペアの値を持つ。初期値は fmin=0, fmax=U。(Uは現在までの最適解。Uの初期値は+∞で良い)。

fmin・fmaxのペアをレイヤーの小さい順からスタックにpushしていったものを、ビームスタック(beam stack)と呼ぶ。

スタートノードから、ビームサーチをしていく（8行目）。注意点は、f-costが[fmin,fmax)に入るノードだけを探索するところ21行目?。ビーム幅が2なので、上位2つ（低いコスト）しか入れない。ビーム幅から漏れたノードの評価値が新たにfmaxとなる（3行目）。つまり範囲が狭まる。

ゴールノードにたどりついたので、最小値を更新し、U=6となる。もし、fmax>Uとなった場合は、そのLayerはEmpty Layerとなる。Empty Layerが、ビームスタックのtopに来たら除去される。つまり深いLayerから探索から外れていく（46～48行目）。

（ここが全く分かってない可能性あり！）最深層までたどり着いたら、最深層のfminとfmaxが更新される。fminが今までのfmax、fmaxがUになる（50～51行目）。

ここはfmaxはそのまま。評価値5のノードは辿っておらず、ビーム幅2からあぶれたわけではない。

深さ優先探索順に新しい訪問先を辿っていく

別のゴールノードにたどり着いたので、最小値を更新し、U=5となる。常に後で見つかった解のほうが、良い最適解になるのに注意する(44～45行目)。これでビームスタックが空になったので、終了（49行目）。

私見と疑問点

アルゴリズムの観点から

• 疑似コードで、どうやってバックトラッキング（深さ優先探索）を実行しているのか？
  • 普通は深さ優先探索をするには、再帰にするか、stackにpush・popをするか、どちらかだけど、8～37行目を見る限り、popはしていないし、再帰呼び出しもない。
  • 関数searchは毎回スタートから行われているので、新たに訪問するノードの順が、深さ優先探索順になるだけかもしれない。
    • ただ、そうだとすると、毎回スタートから同じノードを何度もたどり、計算量がひどいことになりそうな…。
• 以下の例は大丈夫なのか？f-cost3のノードを訪問する前に、f-cost4と5のノードを訪問するので、fmaxが6になり(3行目)、その後すぐfminが6になるので、f-cost3のノードを訪れずに終わってしまいそう。

• wikipediaでビームサーチを見たら、そこでもヒューリスティック関数の例がでていたし、先行する論文のいくつかもそうだった。もしかして、それが普通？
• ビームスタックサーチを分枝限定法（branch and bound）と呼ぶのは大げさな感じが…。ただの枝刈りでは…。でも、英語版wikipedia見ると、ヒューリスティック関数で枝刈りするだけでも、そう呼べるみたい。

マラソンマッチの観点から

• 許容的なヒューリスティック関数が必要だとすると、使用箇所は限定されそう。
  • 最短距離のようにコストが増えていくだけなら、許容的なヒューリスティック関数を常に0とみなしても、よさげですが。自由に値をとりうる評価関数だと、枝刈りが行われず、ビームスタックサーチは、ただの探索になってしまいます。
• マラソンマッチの探索上手な人だと、そもそも、幅優先だか深さ優先だか最良優先だかあまり区別のつかないような、自由な探索を書く人も多く、まさにビームスタックサーチは、その中に含まれそう。
• ビームスタックサーチとほぼ同じと言われてたchokudaiサーチとは、アルゴリズムも目的も全然違うっぽい。

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

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

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

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

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

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

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

ボツネタ1 : ID分解

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

まず、データを眺めてください。

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

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

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

まぁ、ただ注意点は、

• 分解前の元の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色。

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

念のため、上記の画像を出力するコードも貼っておきます。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のランダムフォレストには２種類の重要度(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)

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

• 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. ランダムフォレストで、決定木を作るときに訓練に使われなかったデータサンプル（ブートストラップサンプリングで外れたデータ）を使います。これを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の変わりに使うといったこともできるでしょう。
• ただ「重要度を低い特徴量を削る」と簡単に書きましたが、貪欲に１番重要度が低い特徴量を削っていけば良いというものでもなく、とても難しいです*3。特徴量間に強い関係があれば、特徴例をけずるたびに、重要度も大きく変わってきます。良い方法があったら（というか研究されていると思いますが）、ぜひ教えてもらえるとうれしいです。一応、自分のアイディアレベルのネタとしては、
  • (2)(3)であれば、複数の特徴量をシャッフルする手はあるでしょう。ただ、（OOBデータ個数）＊（木の本数）＊（特徴量の数 ^ シャッフル個数）＊（特徴量1つあたりのシャッフル予測の回数 nPerm）回となるので、特徴量が多いときは、シャッフル個数2個にするのすら遅すぎて難しいかもしれません。
  • 探索枝仮り・焼きなまし法・ビームサーチといった最適化手法で削っていく方法もあるでしょう。
  • 特徴量を１つ削ったとき、他の重要度がどう変化するかで、どの特徴量間に強い関係があるかというグラフのようなものを作れば、最適化手法を使っていく上で役に立つかもしれません。例えば、「特徴量はx1,x2は片方を削ると、もう片方の重要度が大きく下がるので、たぶんペアで意味がある特徴量だ。焼きなまし法を試すときは、この２つはまとめて追加したり削除したりするようにしよう。」といったことはできそうな気がします。

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

最初にあげた図では、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は、飛行機のフライトデータから、違反フライトを予測するという問題でした。問題やデータはこちら

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

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

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

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としてそのまま答えにすることもできるのですが、線形回帰とランダムフォレストの２段構えになっています。

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

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

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

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

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

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

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

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

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

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

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

2015年 競技プログラミングの目標

評価

目標：2015年は100点以上

200～	偉大すぎるので、誰かが奢ってくれるはず
150～199	PERFECT
120～149	GREAT
100～119	GOOD
60～99	進歩なし
30～59	怠惰・堕落(FUJIYAMAに乗る)
0～29	人としてダメ(FUJIYAMAに乗る)

(なお、2014年は60点でした)

得点表

SRMに56.6%以上参加*1	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点
プロコン賞金10万円ごとに	10点
プロコンでオンサイト決勝出場	50点
ニコ生放送200枠分（=100時間）	10点
海外でプロコンオフ会をやる	10点
Codeforcesで栗原市ランキングいり	10点

*1:2015/8/2 SRM30回参加にしてましたが、当初予定されていた月4回SRMでなくなってしまったので、割合に変えます。30回/53回(SRM48回+TCO5回)=56.6%