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

勾配ブースティング落穂拾い - マルチクラス分類について

機械学習 XGBoost

勾配ブースティングに特有の問題という訳ではないですが、それはマルチクラス分類をone-vs-allで解いていることに留意しましょうというのが今回の記事です。

2017/04/02追記 LightGBMではマルチクラス分類の際にsoftmaxとone-vs-allのどちらで解くか選択できるようです。 また、XGBoostではsoftmaxのみがサポートされています。 この記事では、マルチクラス分類のためにクラス数と同じ分類器を使っていることを指してone-vs-allとしてしまっています。ご注意ください。

マルチクラス分類について

ここではマルチクラス分類の問題を対象に、アンサンブル学習として代表的なランダムフォレストと勾配ブースティングとを比較します。

ランダムフォレスト

ランダムフォレストによるマルチクラス分類では学習器に分類木を用いています。 scikit-learnではDecisionTreeClassifierです。

コード

分類木はマルチクラスを自然に扱えるので、n_estimatorsと同じ数だけ木を作れば求めるモデルが構築できます。

勾配ブースティング

勾配ブースティングでは弱学習器に回帰木を用いているので事情が変わってきます。 マルチアウトプットの回帰木は一般的ではありませんので、問題をone-vs-allに変換して解くことになります。 そのため勾配ブースティングでは指定したn_estimatorsとクラス数の積の数だけ木を作る必要があります。

コード

勾配ブースティングは速度と精度を高いレベルで両立していますが、マルチクラス分類では

  • クラス数と比例して計算時間が増加する
  • Dark knowledge的なものを利用できない

という欠点が出てくるので、勾配ブースティングのメリットがスポイルされていないか留意する必要がありそうです。

数値例

勾配ブースティング(というかone-vs-all)が苦手そうな例を考えてみました。 セルの中は(クラスAの割合, クラスBの割合)を示しています。

sex/height
(5%, 10%) (15%, 10%) (20%, 20%)
(15%, 15%) (5%, 15%) (20%, 30%)
(20%, 25%) (20%, 25%)

この例にone-vs-allを使うと、クラスBについてsexで分割する以上のものが出てきません。 一方で分類木を使った場合は、クラスBについてsexで分割した後でクラスAの濃い部分を抽出するheightでの分割を発見できます。

マルチクラス分類をXGBoostのみで取組むのは微妙な感じがしますね。 あまり業務でスタッキングしたくないのですが、それも避けられないように思われます。 もしくはDeep Learningするか。

参考

ブースティングやDeep Learningが一般化した時代におけるマルチクラス分類の解き方についてはOtto Group Product Classification ChallengeのForumで盛んな議論があったようです。 実戦的なテクニックも数多いようなので、関心のあるかたはそちらで調査されるのが良いと思われます。

勾配ブースティング落穂拾い - 木の構築について

機械学習 XGBoost

このシリーズについて

XGBoost芸人を自称してちょこちょこ活動をしてきたのですが、最近になって自分の理解の甘さを痛感するようになりました。

気になった箇所を散発的に復習しているのですが、その成果を備忘録として残しておこうと思います。 今のところ体系的にまとめるつもりはないので、これを読んでも勾配ブースティングの全体像はつかめませんので悪しからず。

今回のテーマ以外にはマルチクラス分類の際の挙動等に関心を持っています。

木の構築について

勾配ブースティングでは  M+1 回目のイテレーションで誤差

 \displaystyle
\sum_{i=1}^N L \left( y_i - \sum_{j=1}^M \eta_j f_j \left( x_i \right) \right)
= \sum_{i=1}^N L \left( y_i - f^M \left( x_i \right) \right)

の勾配を上手く表現した木を構築します。

この部分の処理についてscikit-learnとXGBoostでの違いを確認します。

scikit-learn

カステラ本に準拠した処理になっています。 勾配の計算は

 \displaystyle
\frac{\partial}{\partial f^M \left( x_i \right)} L \left( y_i - f^M \left( x_i \right) \right)

となり、これを各サンプルのラベル扱いにして DecisionTreeRegressor に投げます。

コード

このとき当然 DecisionTreeRegressor は2乗誤差*1を最小にするよう葉に属するサンプルにスコアを付けます。つまり平均値となっています。

ロス関数がRMSEのケースではそれで問題ないのですが、他の関数ではズレが生じます。 そのため、木の分岐は DecisionTreeRegressor の結果を使うのですが、最終的なスコアには補正をかけて対応します。

例えば、ロス関数がMAEの場合は、葉に属するサンプルの中央値をスコアとして採用するのが、木の構造を所与としたときの最適となります。

コード

XGBoost

XGBoostでは正則化項が追加されているので最適化の目的関数が異なります。式の導出は別に譲りますが以下のようになります。

 \displaystyle
\sum_{j=1}^M \left( G_j w_j + \alpha \left| w_j \right| + \frac{1}{2} \left( H_j + \lambda \right) w_j^2 \right) + \gamma T

木の構造を所与としたとき最適な  w_j^*

 \displaystyle
w_j^* = \frac{l \left( G_j , \alpha \right)}{H_j + \lambda} \\
\mathrm{where} \quad l(x, y) =
\begin{array}{l}
x - y \quad \mathrm{if} \ x \geq y \\
x + y \quad \mathrm{else}
\end{array}

となります。これは最適化の1階条件

 \displaystyle
G_j + \alpha \mathrm{sign} \left( w_j^* \right) + \left( H_j + \lambda \right) w_j^* = 0 \\
\displaystyle
w_j^* = -\frac{G_j + \alpha \mathrm{sign} \left( w_j^* \right)}{H_j + \lambda}

によります。ここで  w_j^* \leq 0 となるためには  G_j \geq \alpha でなければなりません。

これを代入すれば目的関数から  w_j が消去できます。 したがってXGBoostでは、ある分割についてはじめから最適なスコアを付ける前提で、最適な分岐を探索していることになります。

ポエム

XGBoostは通常の勾配ブースティングより精度が良いとされており、その原因は主としてL1/L2正則化にあるとされています。 確かに通常の勾配ブースティングでもインサンプルの誤差を0にするのは簡単なので、汎化性能の高さは正則化によるものと考えるのが自然です。

しかしながら、今回見てきたようにXGBoostは通常のものより精度の高い近似計算をしています。 個人的には正則化がそれほど効いていると思えないので、最適化の工夫によって高い汎化性能が実現しているのではないかと疑っています。 XGBoostの収束が速い分だけ、同水準のインサンプル誤差で比較したときに森の複雑性がより小さいと期待できるのではないでしょうか。

参考

*1:オプションで絶対誤差も指定可能

4教科で珍しく良い点取った生徒のカンニングを疑ってはいけない、「竜王を弁護する」について

統計

はじめに

将棋ソフト不正疑惑に関する事件について新たな記事が出ていたので、また計算機を叩いてみました。

三浦九段不正疑惑について、渡辺明竜王を弁護する

2017/02/27追記 当該記事は削除されてしまっています。一時、ハフィントンポストが全文を掲載してくれていたのですが、そちらも削除されています。ご覧なりたい方は自身で魚拓を漁ってください。

三浦九段については晴れて「シロ」判定がされて、これから名誉回復をどうしていくかという段で、このような記事が出てくるのには違和感があります。

ネット上でもあまり評判良くないようなので、今更内容に噛み付く必要もないのですが、データジャーナリズムの手習いとして公開してみます。

前提

ここでは、記事中にある一致率の推定値を真の値として扱います。

一致率については、将棋ソフトや局面の違いによって値に大きな差異が出る可能性が指摘されています。 そのような数字を基に組み立てたロジックは無意味だと一蹴することもできるのでしょうが、一致率を受け入れた上でも反論が可能と私は考えました。

設定

記事中では、ある対局で高い一致率の示される確率について2つの数字が提示されています。

  • 10%: 「当方の計算では、出現率でいえば、良くて10%程度」
  • 25%: 「第三者委員会の報告は、(中略)上位4分の1に入るような出来」

はっきりと「ある対局で高い一致率の示される確率」と書かれている訳ではないですが、文脈からそう判断して問題ないと思われます。

私も以前の記事で計算をしたのですが、ある対局で一致率85%以上となる確率を9.79%と概算しています。 不安定な一致率の推定値から計算したものなので、これが正しいと考えてもいませんが、ここでは三浦九段に不利な10%を採用してみます。

確率計算

弁護記事では「そんな出来の良い将棋が、疑いをもたれた短期間に、4局も」あったと指摘しています。「短期間」が何を指すか明らかではありませんが、私の前回記事では某まとめサイト情報を参考に12局としていました。

12局中4局以上で高い一致率が示される確率

r   = 0.1
num = 12
lim = 3
s   = 1.0 - scipy.stats.binom.cdf(lim, num, r)

2.56%となります。 教科書にあるような有意水準5%だと、通常では起こりにくいことが起こっていて疑わしいとなります。

以前の記事では、12局というのが恣意的に切り取った一部なのでシロともクロとも言えないとコメントしていましたが、よくよく考えれば恣意的な切り取りについても数値的に評価できそうなので、その計算を今回追加します。

一致率の高くなる期間を切り取れる確率

三浦九段の2015年度対局数は42局です。 42局の中に一致率の高い4局を含む12局が出現するか、という問題に変わると、この確率は当然ながら2.56%より高くなります。

解析的に解く方法が思いつかないので、モンテカルロ法を使ってみます。エレガントな解法があれば教えてください。

確率10%のクジを42回引いて、結果のrolling_sumの最大値が3を超えるかで判定しています。

num_ext  = 42
num_iter = 10000
res = np.zeros(num_iter)
for ii in range(num_iter):
    samp = np.random.choice([0, 1], num_ext, replace=True, p=[1.0-r, r])
    seq  = pd.Series(samp).rolling(window=num).sum()
    res[ii] = seq.max()
s = (res > lim).mean()

モンテカルロ法なので結果にブレはありますが、17%程度という結果になりました。

トップ棋士は、ベースの一致率が高くなりがちというだけでなく、対局数つまり試行回数が増えてしまうため、一致率による疑いをかけられやすい構造にあるようです。

計算は以下にまとめています。 github.com

主張

サンプルを恣意的に切り取ることで、17%程度発生するありふれた事象を2%程度に見せることができます。 我々が統計の結果を作る・見るときには、このようなテクニックが使われていないか常に注意する必要があります。

いつもソフトを使ってる訳でないので、一致率の高い4局が問題で、12局だとか42局は関係ないという論を張っている人もいるようですが、これは誤りです。 無実の人が同じ試行をしたときに起こりそうもない結果が出たときのみ疑いをかけるのが正しい姿勢です。 羽生さんの対局から一致率の高い1局を引っ張ってきて、いつもはソフト指ししていないけれど、当該局ではソフト指ししていたという主張は馬鹿げているでしょう。

また、弁護記事では試験のカンニングについて言及されていますが、もし教師の方がいらっしゃたら記事のような手続きで生徒を疑うのはやめてください。個人的に疑いを持って物的証拠を探すのはセーフかもしれませんが、疑いを同僚教師に話したりするのはアウトです。 過去の偏差値からは取れそうもない点数を取る程度の事象であれば、「卒業までに受ける試験の数×カンニングが疑われるような素行の生徒数」の回数の試行があれば、簡単に発生してしまうのが確率というものです。

最期に

最近、話題のBIGについて計算したのですが、そちらでは直感通りのつまらない結果になりました。

MXNetを中心としたCustom Loss functionの話

機械学習 Python Deep Learning

これはDeepLearning Advent Calendar 2016の12月24日分の記事です。

Deep Learning一般の内容ではないので申し訳ないですが、来年はMXNet Advent Calendarやりたいですねという願いもあり、空き枠に滑り込みました。

TensorFlowとChainerが強くてMXNetは肩身が狭いですが、AWSという錦の御旗も得たのでユーザ増加に期待したいです。 MXNetはドキュメントが弱いとよく言われますが、最近はドキュメントサンプルコードも充実してきています。 9月頃はドキュメントがリンク切れしておりサードパーティーのドキュメントを読んでいたくらいなので随分良くなりましたよ。

さて、今回はCustom Loss functionについて紹介しようと思います。機能自体は前からあったようですが今月になって簡単な資料が公開され、一般ユーザにも存在が知られるようになったものです。

ロス関数について

ここでは、モデル出力値とラベルの差を表現する関数で、勾配法によって最小化する対象をロス関数と呼びます。 どのライブラリでも回帰なら2乗誤差、分類ならクロスエントロピーが使われているはずです。

注意が必要なのは、モデル精度の最終的な評価をするための指標とロス関数は別になるケースもあるということです。 分類問題ではクロスエントロピーよりもAccuracyやAUCで精度評価することが多いと思いますが、そういった場合でもロス関数にはクロスエントロピーが選ばれてクロスエントロピーを最小化するのが通常です。

評価指標をAUCに変えたら、それを最適化してくれると思っているユーザは結構多く、私もそういった方への案内誘導をボランティアでやっています。*1

一方で、ロス関数に拘ったからといって報われるとは限らず、絶対誤差で評価されるタスク等でも2乗誤差を最小化する方が上手くいくケースも多いのが実際と思われます。

ごく限定的な状況ですが、標準のロス関数では対応できないような問題のためにロス関数をカスタムするインターフェースが各ライブラリから提供されています。

MXNetでの事情

MXNetのLinearRegressionOutputSoftmaxOutputといった出力層向けのシンボルは、出力とロス関数が一体化しておりパラメータでロス関数を変更できません。 ロス関数を変えたい場合はMakelossというシンボルを利用する必要があります。

回帰の例

LinearRegressionOutputを使った通常の方法は以下のようになります。これはMSEを最小化します。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=1)
net1  = mx.sym.LinearRegressionOutput(fco, label)

同じようにMSEの最小化をCustom Loss functionでやると以下のようになります。

loss = mx.sym.square(
    mx.sym.Reshape(fco, shape=(-1,)) - label
)
net2 = mx.sym.MakeLoss(loss)

mx.sym.squareやelement-wise minusの導関数はMXNetが知っているので自動微分を使って上手く最小化してくれる訳です。 numpy等を使って好き勝手にロス関数を作ることはできません。

注意 MakeLossにmx.sym.squareした結果を渡すとpredictした際にも2乗した結果が出力されるそうです。

custom loss symbol in R/Python · Issue #3368 · dmlc/mxnet · GitHub

確かにどれが出力として求められているかはMXNetで判断しにくいかもしれません。これについては使い易いように修正がされるそうです。

分類の例

MXNetのLogisticRegressionOutputは名前のような挙動をしないことが知られています。

RPubs - MXnet LogisticRegressionOutput test

SoftmaxOutputを使うことが推奨されていますが、これは手前のレイヤーで隠れユニット数をカテゴリ数と同じだけ確保する必要があります。 2値分類のときに隠れユニットを2つ持っておくのは冗長な気がするので、Custom Loss functionでロジスティック回帰を組んでみようというのがこの試みです。

SoftmaxOutputを使った通常の方法は以下のようになります。これはLog Lossを最小化します。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('sm_label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=2)
net1  = mx.sym.SoftmaxOutput(fco, name='sm')

ロジスティック回帰らしく書いた場合は以下のようになります。

data  = mx.sym.Variable('data')
label = mx.sym.Variable('label')
fc    = mx.sym.FullyConnected(data, num_hidden=3)
act   = mx.sym.Activation(fc, act_type='relu')
fco   = mx.sym.FullyConnected(act, num_hidden=1)
p     = mx.sym.Activation(fco, act_type='sigmoid')
eps   = 1e-6
p     = mx.sym.minimum(mx.sym.maximum(p, eps), 1.0-eps)
q     = 1.0 - p
lp    = mx.sym.Reshape(mx.sym.log(p), shape=(-1,))
lq    = mx.sym.Reshape(mx.sym.log(q), shape=(-1,))
loss  = label * lp + (1.0 - label) * lq
net2  = mx.sym.MakeLoss(loss)

出力を確率に直すやり方がわからないので、ちゃんと動いているかわからないのですが、 方針自体はこれで良いはず…

中途半端になってしまいましたが、これ以上掘るのは効率が悪いので、本日はここまでとします。 手元のJupyter Notebookはこちらです。

github.com

やっぱり

以前に比べればドキュメントが強化されているとはいえ、最近どういう機能が追加されたか等を個人で調べるのは大変です。 ユーザが増えて互助していければ素敵ですね。

*1:こんなことは働き盛りの人間がやる仕事ではないので、ボット化したいなぁ

LightGBMのPythonパッケージ触ってみた

機械学習 Python

DMLCひとりアドベントカレンダー0日目の記事です。 強い競合が現れたということで、DMLCとは直接関係ないですがLightGBMについて紹介します。

LightGBMとは

勾配ブースティング木の高速な実装としてXGBoostが有名ですが、Microsoftの開発した更に高速な実装がLightGBMです。 実験によるとXGBoostの数倍高速だということです。

パフォーマンス向上の大きな要因は、決定木の分割をexactに探索するのではなく、ヒストグラムのビン単位にして計算を間引いていることにあるようです。 ビン単位にすると精度が落ちそうな気もしますが、汎化性能を考えたら特徴量が粗いくらいでちょうど良いのだろうと腹落ちしています。

同じアプローチでFastBDTが先に発表されましたが開発が滞ってしまっている一方、LightGBMでは最近Pythonパッケージがベータリリースされました。 今日、我々がPythonで勾配ブースティングをする際にはXGBoostかLightGBMの2択*1となります。

導入

ビルドはこちら、インストールはこちらを参考に。

XGBoostとそれほど変わりません。

使い方

Exampleより

xgboost.DMatrix的なものを作って…

# create dataset for lightgbm
lgb_train = lgb.Dataset(X_train, y_train)
lgb_eval  = lgb.Dataset(X_test, y_test, reference=lgb_train)

パラメータをディクショナリに格納して…

# specify your configurations as a dict
params = {
    'task' : 'train',
    'boosting_type' : 'gbdt',
    'objective' : 'regression',
    'metric' : {'l2', 'auc'},
    'num_leaves' : 31,
    'learning_rate' : 0.05,
    'feature_fraction' : 0.9,
    'bagging_fraction' : 0.8,
    'bagging_freq': 5,
    'verbose' : 0
}

trainするだけ。XGBoostと同じ流れですね。

# train
gbm = lgb.train(params,
                lgb_train,
                num_boost_round=20,
                valid_sets=lgb_eval,
                early_stopping_rounds=5)

predictpandas.DataFrameを直接受け入れてくれるのには注意が必要です。

# predict
y_pred = gbm.predict(X_test, num_iteration=gbm.best_iteration)

XGBoostのパラメータとの対応は以下のようになります。

2016/12/18追記 tksさんに指摘を頂いたので記述を修正しました。 github.com

LightGBM XGBoost 候補 備考
objective objective regression, binary, multiclass, ...
boosting_type booster gbdt, dart 線形モデルはなし
max_depth max_depth num_leavesの制約を受けるので注意
learning_rate learning_rate
feature_fraction colsample_bytree
bagging_fraction subsample bagging_freqを1以上にしないと有効にならない
min_gain_to_split gamma
lamnda_l1 alpha
lambda_l2 lambda
metric eval_metrix l1, l2, auc, ...

LightGBMにもmax_depthパラメータが用意されていますが、オフィシャルにはnum_leavesで葉の数を制御します。

また、特徴量をまとめるビンの数をmax_binで設定できます。

詳細は公式のWikiを調べて下さい。

カテゴリ変数のサポート

LightGBMではカテゴリ変数がサポートされています。 One-hot encodingが不要になったので、メモリ的に有利であったり、カラム方向のサンプリングやfeature importanceが直感的になる等のメリットがありそうです。

下の実験からのコード片ですが、traincategorical_featureパラメータにカテゴリ変数として扱うカラムのインデックスをリストで渡してやるだけです。 カテゴリ変数はintのみ対応しているので、文字だったりする場合は0始まりのコードに置き換える必要があります。

obj_lgb = lgb.train(
    prm_lgb, dt_lgb_c, num_boost_round=num_round,
    valid_sets=dv_lgb_c,
    categorical_feature=list(range(len(x.columns.values))))

処理を完全には追っていないのですが、rpartのようなカテゴリ併合をしてくれる訳ではないようで、数値変数であればx <= thresholdで分割されるところが、x is categoryで分割されるようになるだけのようです。

また以下の実験の出力からなのですが、数値変数として扱うと以下のように小数点の閾値で分割されます。

threshold=1.5 3.5 2.5 0.5 0.5 0.5 0.5 0.5 4.5 14.5 0.5 1.5 0.5 0.5 1.5 2.5 0.5 0.5 0.5 0.5 0.5 4.5 0.5 0.5 0.5 0.5 1.5 0.5 2.5 2.5

一方、カテゴリ変数として扱うとデータにある値そのもので分割されます。

threshold=0 0 0 0 1 1 3 2 0 1 0 0 0 1 1 0 0 3 1 1 0 3 0 0 0 0 1 0 1 0

実験

手元でもパフォーマンスを確認してみようと、Ottoのデータで実験してみました。

github.com

モデル MLogLoss 計算時間
XGBoost 0.570366 65.20
LightGBM 0.522378 16.02
LightGBM Cat 0.671286 15.78

この設定ではLightGBMが4倍以上速い結果となりました。精度もLightGBMの方が良好です。

全変数をカテゴリ変数として扱ったLightGBM Catの有り難みがないように見えますが、One-hot encodingによってカラム数が膨らんでいる場合には計算時間の短縮が実感できるはずです。 精度が悪いのは、今回のデータが順序性を持っているためと思われます。

最後に

XGBoost4jのような使い方をすると話は別かもしれませんが、クライアントマシンでせこせこ計算する分にはLightGBMで良さそうです。 開発も活発なので、XGBoostにしかないような機能も続々と取り込まれている上に、コードもスッキリしており改造が比較的容易になっています。

CNTKなんかもパフォーマンス良いようですし、巻き返しに力の入っているMicrosoftさんのLightGBMには将来性も期待できるのではないでしょうか。

次は23日にMXNetのネタを考えています。

*1:scikit-learnのGradientBoostingもMAEで分割ができる等、面白そうではあるんですが…

ピュアオーディオ的な解像度で語る乱数のコク

統計

周回遅れも甚だしいですが、「乱数のコク」について思う所があったので計算をしてみました。

オリジナル

一様乱数と比べて、5つの一様乱数の平均の方がコクがあるとのことです。 シミュレーションで得たヒストグラムは以下のようになりました。

def rand_koku(n, seed):
    np.random.seed(seed)
    u1 = np.random.rand(n)
    u2 = np.random.rand(n)
    u3 = np.random.rand(n)
    u4 = np.random.rand(n)
    u5 = np.random.rand(n)
    return (u1, (u1+u2+u3+u4+u5)/5.0)

f:id:marugari2:20161108002149p:plain

平均を取った側は正規分布のような釣り鐘型になっています。 果たしてこれをコクと呼べるでしょうか。2次モーメントも異なっていますし、分布の形も別物です。 オーディオで言ったらSONYBOSEくらい差があります。

乱数のコクが感じられる例

今回は2種類の正規乱数を用意しました。

1つはBox-Muller法による正規乱数です。 これはnumpyの正規乱数生成にも用いられているベーシックなものです。

def normal_box_muller(n, seed):
    np.random.seed(seed)
    nm = int(n / 2)
    u  = np.random.rand(n)
    u1 = u[:nm]
    u2 = u[nm:n]
    r  = -2.0 * np.log(u1)
    v  = 2.0 * np.pi * u2
    z1 = np.sqrt(r) * np.cos(v)
    z2 = np.sqrt(r) * np.sin(v)
    z  = np.concatenate([z1, z2, -z1, -z2])
    return z

もう1つはMoroの逆関数法による正規乱数です。 正規分布の累積分布関数の逆関数を利用して、一様乱数から正規乱数を得ることができます。

def normal_moro(n, seed):
    np.random.seed(seed)
    a0  = 2.50662823884
    a1  = -18.61500062529
    a2  = 41.39119773534
    a3  = -25.44106049637
    b0  = -8.47351093090
    b1  = 23.08336743743
    b2  = -21.06224101826
    b3  = 3.13082909833
    c0  = 0.3374754822726147
    c1  = 0.9761690190917186
    c2  = 0.1607979714918209
    c3  = 0.0276438810333863
    c4  = 0.0038405729373609
    c5  = 0.0003951896511919
    c6  = 0.0000321767881768
    c7  = 0.0000002888167364
    c8  = 0.0000003960315187
    u   = np.random.rand(n)
    y   = u - 0.5
    f1  = (np.abs(y) < 0.42)
    y1  = y[f1]
    r   = y1 ** 2
    num = ((a3*r+a2)*r+a1)*r+a0
    den = (((b3*r+b2)*r+b1)*r+b0)*r+1.0
    x1  = y1 * num / den
    f2  = np.logical_not(f1)
    u2  = u[f2]
    y2  = y[f2]
    r   = np.where(y2 > 0.0, 1.0 - u2, u2)
    r   = np.log(-np.log(r))
    x2  = c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*(c6+r*(c7+r*c8)))))))
    x2  = np.where(y2 < 0.0, -x2, x2)
    z1  = np.concatenate([x1, x2])
    z   = np.concatenate([z1, -z1])
    return z

アルゴリズムはどちらも Paul Glasserman, Monte Carlo Methods in Financial Engineering を参考にしています。

乱数を200,000,000個発生させて得たヒストグラムは以下のようになりました。

  • 両側0.5シグマ f:id:marugari2:20161108003322p:plain 両側0.5シグマの範囲では差を見て取ることができません。

  • 片側4シグマ f:id:marugari2:20161108003533p:plain 同じ一様乱数系列を用いていますが、正規分布への変換法によって片側4シグマの裾の分布が異なっています。

英語版WikipediaではBox-Muller法の弱点について記載があります。 裾の分布に拘る場合はMoroの逆関数のようなものを使うことが推奨されているようです。

何気なくrandnで正規乱数を得ることも多いですが、そのアルゴリズムの違いによって得られる結果が異なることを見てきました。 こういったところにこそ「乱数のコク」があるのではないでしょうか。

最後に

正規分布は非常に基礎的な確率分布ですが、それに従う乱数を得るためのアルゴリズムは(私にとっては)そんなに簡単ではありません。 正規乱数でさえそうなのです。

ベイズ統計では、複雑な事前分布と複雑な尤度関数の積からなる事後分布に従うサンプルを生成したいケースがあるのですが、そのような複雑な分布関数にはBox-Mullerのような解析的なアプローチもMoroのような精度の良い近似関数も存在しないことが常です。 そのようなマトモでない分布に従う乱数を生成するためにMarcov Chain Monte Carloが必要となっている訳です。

こういった話はRobert and Casella, Monte Carlo Statistical Methodsに感動的にわかりやすく書かれているので、私は強くオススメをしています。

一致率90%以上はどれくらいあり得るのか計算してみる(ガバガバ設定)

統計

将棋界で大騒ぎになっている疑惑について、連盟公式掲示板でのレスバトルをするための材料として手元で計算をしてみました。

計算の前提となっているパラメータがガバガバなので、白か黒かの材料にはなりようがないですが、とりあえず今出ている数字を使ってみたらどうなるかというメモです。

私は

  • 統計検定持ってません。変な計算してたら叱って下さい。
  • 調査設計については無知です。調査法がイケてない点については触れません。
  • 観る将なので、戦法や序盤・終盤でどれくらい手の広さが異なるか知りません。
  • 電王戦は見てるけど、ponaと技巧でどれくらい強さが違うかは知りません。

設定

千田先生のDropboxにある資料の数字から独自に作成します。 「千田棋譜 資料」内に当該ファイルがあります。 www.dropbox.com

  • ソフト指ししないプロ棋士の一手毎のソフト一致率は75%とする。
  • 一致するかどうかは系列相関なしで独立とする。

2点目の仮定について、一本道の変化だと一致率の高い局が多くなるはずなのでシロ派に不利な仮定と考えられます。一方、全部の手をソフト指しする必要はなく、ソフト指し側が一致率をコントロールできるという点ではシロ派に有利な仮定とも考えられます。 これは水掛け論になるので、計算を簡単にするために独立としておきます。

計算結果

一局で一致率が90%以上になる確率

序盤を通過した後で終局まで30手*1指すとします。27手以上一致すれば率が90%以上になります。

p   = 0.75
num = 30
lim = int(num * 0.9) - 1
q   = 1.0 - scipy.stats.binom.cdf(lim, num, p)

3.74%になりました。結構高いような気もします。

実績値と比べて検討もできますが、私にはヘビーなのでこの数字を使うことにします。

一致率の高い対局が複数出る確率

疑惑の7月以降12局*2で一致率90%以上が複数回現れる確率はどうなるでしょう。

r   = q
num = 12
s1  = 1.0 - scipy.stats.binom.cdf(0, num, r)
s2  = 1.0 - scipy.stats.binom.cdf(1, num, r)
s3  = 1.0 - scipy.stats.binom.cdf(2, num, r)

1局以上現れる確率は36.7%、2局以上現れる確率は7.2%、3局以上現れる確率は0.9%です。 三浦久保戦の一致率が高かったのも考慮して、2局あったとすると*3確率は7.2%です。

棋士を十数人集めれば、1人くらいは石に当たる奴が出てくる計算ですね。

90%という区切りが悪かった?

疑惑の4局は一致率が85%以上です。区切りを変えてみましょう。

p   = 0.75
num = 30
lim = int(num * 0.85)
q   = 1.0 - scipy.stats.binom.cdf(lim, num, p)

9.79%です。

r   = q
num = 12
lim = 3
s   = 1.0 - scipy.stats.binom.cdf(lim, num, r)

12局中で4局以上現れる確率は2.4%です。

この設定だと5%は割ってきますが、12局がいいとこ取りなのをどう見るかでしょうか。

所感

最初から分かっていたことですが、一致率でクロ・シロの議論はできそうもありません。 しかし、当初想定していたよりはシロを支持できる数字が出たようにも思えます。

レスバトルで盛り上がる前に計算機を叩くよう心がけたいです。

*1:対局毎に終局までの手数は異なりますが、計算を簡単にするため30で固定しています。

*2:いいとこ取りになるので、この局数で良いのかという議論はあると思います。

*3:ある一局の一致率が高かったからといって、ヒット数を多くするのは適正ではないですが、サンプルも少ないのでとりあえず。