Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 02-04-2010
さいころをいっぱい振ったときの最大値の 話．再び．
例えば
１：普通の６面のサイコロを４つ振る．一番大きい目が３の確率は？
２：４面，６面，８面のサイコロを２つずつ振る．一番大きい目が５の確率は？
ちなみに答は，１が 2.7%で，２が 16.0%．
この，サイコロ最大値の問題は，計算はやや面倒くさいけれど，考え方もプログラムも極めてシンプル．
サイコロの数を n，それぞれの面の数を m として，最大値が j になる確率は
非常にシンプルなことをめんどくさく書いてみた・・・
（シンプルな）考え方：４面，６面，８面のサイコロを１つずつ振る場合を例にする．
まず最大値が１になる場合を考える．３つのサイコロとも１を出さなくてはいけないので，その確率は
P1 = (1/4)(1/6)(1/8)
次に最大値が２になる場合：３つのサイコロとも１か２を出さなくてはいけないので，その確率は (2/4)(2/6)(2/8) ぽいのだけど，３つとも１だった場合 （確率は P1）を除かなくてはいけないので
P2 = (2/4)(2/6)(2/8) - P1
最大値が３になる場合：３つとも１か２か３を出さなくていけないんだけど，P1+P2 を引かなくてはいけないから
P3 = (3/4)(3/6)(3/8) - ( P1+P2 )
それで，上の式のような 掛け算 - 足し算 になる．最大値が j 以下になる確率は？だったら，マイナス足し算の部分がいらない．
で，min( j, mi ) なんだけど，これはただ単に，例えば「最大値が５」の確率の時は４面のサイコロはどうがんばっても貢献できないので，(5/4) となってしまう代わりに (4/4) にする処置．例１のように目の数が同じサイコロしか出てこない場合は無視してよい．
これを計算するコードはＲなら２行で書ける．無理すれば１行で済むけれど．given.j <- function(j,v){ prod( pmin(v,j)/v) }
diff( sapply(as.list(0:max(v)), given.j, v=v) ) v が面の数．例えば例題２だったら， v <- rep(c(4,6,8),each=2) としてから走らせればよい．