二群間の勝利確率の推定について

id:yeomanさんのところで面白い分析をやっていた。
箱ひげ図による可視化はセンスが良いと思ったんだけど、後半のボルトvsガトリンの勝利確率の導出はいろいろと怪しい部分があったように思う。

「諸分野に明るい方ご意見を宜しくお願い致します」というコメントもあるし（こういう姿勢で色々なことに挑戦しているのでこの人のブログは好きだ）、考え始めたらドツボにはまっていって結構な分量になったのでせっかくだから書き出しておく。
一応なるべく分かりやすく書いたつもりだけど、統計統計した話なので興味が無い方は回れ右。

※ブコメで効率のよい解き方のご指摘をいただいたので、結論から知りたい方は末尾の追記を読んでください

問題設定

元記事ではボルトとガトリンのタイムはそれぞれ正規分布に従うと仮定している。
箱ひげ図を見た感じでもほぼ平均に対して対称なので、この仮定はこの二者の比較においては妥当と言ってもいいと思う。

パラメータは書かれていないが、表と対応を取って推測するとボルトが平均9.84、標準偏差0.20。ガトリンが平均9.94、標準偏差0.145。
元記事の裏で用いられているパラメータとそんなに隔たってはいないはず。

後半部分で行われていたのは、上記仮定のもとで、「ガトリンがボルトに勝てる可能性は何パーセントだったのか？」という分析。

ここで行わないといけないのは、「それぞれ異なる正規分布に従う母集団A,Bがあり、そこからそれぞれ1つずつ無作為標本a,bを得た時に、a<bである確率はどの程度か」という問題を解くこと。
問題設定としては単純なんだけど、実際に解くとなるとそんなに自明ではない気がする。

趣味でソシャゲの解析とかやってるとわりとこの手の計算はするんだけど、具体的にこれがどう呼ばれているのかは知らない。
それらしい解説を行っているサイトも見当たらないので、あまり筋のいい名前ではないが適当に「二群間の勝利確率の推定」と名づけて検討してみる。

当方統計の専門家ではないので何かおかしなことを言っているかもしれない。ご指摘があればお願いしたい。

解法1：モンテカルロを回す

いちばん簡単な解法として、実際に二つの正規分布から無数に無作為標本を取って比較して、勝利しているケースの割合を数えるというものがある。
これはベイズ推定で使われるモンテカルロ法と同じ考え方で、十分な数のサンプルからヒストグラムを作ると、ヒストグラム中のパラメータは求めたい分布の真のパラメータに収束していく。
言ってみれば「実際に十分な回数やってみて集計するだけ」なのでコンピュータに任せるぶんには非常に楽である。普段この手の分析はこれで処理している。
100万程度サンプルを取ればまず問題なく十分な精度の近似解になるはずなので、R言語でそういうプログラムを書いて回してみた。

set.seed(0)
b.mean=9.84
b.sd=0.20
g.mean=9.94
g.sd=0.145

# それぞれ100万個サンプルを作ってボルト-ガトリンの差を取る
N = 1000000
b = rnorm(N, b.mean, b.sd)
g = rnorm(N,g.mean, g.sd)
diff = b-g

# 差のうち、値がゼロより大きいものの数を数えて100万で割る
length(diff[diff>0])/N

実行結果は0.342921となった。実際に100万回試したら34万回勝っていたということ。約34%の確率でガトリンはボルトに勝てるということになる。

補足：ほんとに3割も勝てるのか？

実際に回してみた結果そうなるのでもちろんそうなのだが、これは感覚的にもそんなにおかしな数字ではない。あえて大雑把な議論をすると、以下のようなことがいえる。

正規分布は左右対称な分布なので、箱ひげ図の中央の太線（中央値、対称な分布なので平均値でもある）の上側が生じる確率と下側が生じる確率は等しい。コインを投げて表が出るか、裏が出るかが図の上側に入るか下側に入るかに対応していると考えてもいい。

つまり、ボルトが上側を取り、ガトリンが下側を取るケースが25%の確率であるわけだ。このケースだけ考えてみれば、結構な確率で勝ちうることが直感的にわかると思う。他のケースについても若干なりガトリンの勝ち目はあるので、それらを足し合わせれば3割程度にはなりそうだ、ということがお分かりいただけるのではないだろうか。

この結果を導いているのは双方のタイムの分散の大きさによるところが大きいと思う。もう少しお互いの調子の波が小さければ、ガトリンの勝ち目はより薄くなることだろう。

解法2：定式化して愚直に手計算して解く

実際に定式化して解こうとするとこれがまた難しい。（結論から言うと解けなかった）

ボルトがある値を出したという条件のもとで
ガトリンがそれを上回る値を出す確率

なので、条件付き確率として考える必要がある。ボルトのタイムをb、ガトリンのタイムをgとすると、 ${ p(g \lt b | b) }$ というのが求めたい条件付きの確率密度関数である。
${ p(g \lt b | b)}$ というのはbより小さい全てのgについて確率を合計したものなので、連続型の確率密度関数である正規分布では積分によって導出する必要がある。

そして、これをボルトが出しうる全ての値について考慮しなくてはいけない。
ボルトがタイムbを出す確率がp(b)だとすると、 ${ p(g \lt b | b) p(b) }$ を、すべてのとりうるbについて合計する必要がある。これも積分計算である。

したがって、

を求めればよいことになる。

累積分布関数 ${F_X(x)=p(X \leq x)=\int_{-\infty}^x f_X(t) dt}$ であるから、内側の積分に関しては正規分布の累積分布関数の定義をそのまま使えばいい。 ……と思ったのだが、そう簡単な話でもなかった。

正規分布の累積分布関数は以下のとおり。

erfってなんだろうと思っても日本語版wikipediaには定義がない。解説もない。（統計関係の日本語版wikipediaはいつもこんな感じであてにならない）仕方がないので英語版を参照して解説を引っ張ってくる。誤差関数というらしい。

誤差関数 - Wikipedia

あまり数学が得意な方ではないので、見た瞬間にうっと思ってしまった。引数が積分の定義域にかかる関数になっている。雲行きが怪しくなってきた。

こいつを計算した上でbについて積分してやる必要があるわけだけど、どうしていいのか見当もつかない。誤差関数の積分の話もwikipediaにあったが、以下のようになっていた。

相補誤差関数の累次積分

これはもう無理ですわ...
解けるかどうかの判断も自分の知識ではおぼつかないので早々に諦めた。

解法3：定式化だけして累積分布関数と積分はプログラムに計算させる

R言語のパッケージにはpnormという関数があって、内部的にどういう計算をしているのかは知らないが正規分布の累積分布関数の値を求めることができる。
数値積分の関数も用意されている。解法2で導いた途中までの式（ ${ \int_{-\infty}^{\infty} p(g \lt b | b) p(b) db }$ ）をプログラムに落とし込んで解いてみることにする。

f = function(b) {
  pnorm(b, g.mean, g.sd) * dnorm(b, b.mean, b.sd)
}

integrate(f, 0,20)

実行結果がこれ。

0.3428103 with absolute error < 5e-06

ありがとうございました。やっぱり34％でした。
誤差があるあたり結局近似してそうに見えるし、なんだか負けた気分だけど仕方がない。一応上記の立式が正しいであろうことはわかった。

まとめ

元々の問題に対する答えとしては3割程度勝てるということでいいと思う。
もちろんあくまで両者のタイムが正規分布からのランダムサンプルに従うという理論上の話で、ボルトがガトリンを相手取った時だけ競争心を燃やすとか、そういう定量化できない条件がつくとどうにもならないんだけど。

むしろ一見単純な問題でも真っ当に向き合うとこれほど闇が広がってるんだなという学びを得た気がする。正規分布はおそろしい。
式展開のほうはちょっと手に負えなかったので、この式もちゃんと手計算で解けるよ！という方がいたらぜひ教えていただきたい。（読んでも自分には理解できない可能性はあるが）

あと、モンテカルロの強力さが身にしみてわかった。乱数発生で片付くものは片付けた方が応用的には手っ取り早い。

追記(2015/12/14 20:20)

二群間の勝利確率の推定について - shigusa_t’s diary
正規分布の差は正規分布になるはずなので、差の分布の0以下の部分の確率の和を求めればよいような
2015/12/14 18:06

ありがとうございます。初めて知りましたが、これが一番妥当な解き方だと思います。再生性という性質があるんですね。

再生性 - Wikipedia

正規分布について、

が成り立つらしいので、単純にパラメータを足し引きした分布について確率の和を計算するだけでいい模様。

pnorm(0, b.mean-g.mean, sqrt(b.sd^2 + g.sd^2), lower=F)

コードはこれ一行で済む。結果も上の二つの手法と同じになった。勉強になりました...