遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

モナドを使うとは〜Haskellの確率的プログラミング

Haskell 確率的プログラミング
  • 例えば、ベルヌーイ事象を繰り返して、その結果がどうなるか、というのをHaskellでやるとする
  • ベルヌーイ事象の確率質量分布をモナドで作っていれば、２回繰り返しは、モナドを作って、それを>>=で次に渡して、もうひとつのモナドと合体させてモナドを作ればよい
  • それだけだったら、別に、モナドでなくても良さそうだ
*Main Lib Numeric.Probability.Distribution> (choose 0.4 0 1) >>= (\x -> choose 0.4 (x, 0) (x, 1))
fromFreqs [((1,1),0.36),((0,1),0.24),((1,0),0.24),((0,0),0.16000000000000003)]
  • ２回繰り返しは以下のようになる
    • 長くてちょっといやかも
    • とくに、>>= (\x -> choose 0.4 (x, 0) (x, 1))を２回やっているところが、いかにも無駄
*Main Lib Numeric.Probability.Distribution> (choose 0.4 0 1) >>= (\x -> choose 0.4 (x, 0) (x, 1)) >>= (\x -> choose 0.4 (x,0) (x,1))
fromFreqs [(((1,1),1),0.216),(((0,1),1),0.144),(((1,0),1),0.144),(((1,1),0),0.144),(((0,0),1),9.600000000000002e-2),(((0,1),0),9.6e-2),(((1,0),0),9.6e-2),(((0,0),0),6.400000000000002e-2)]
  • ここで、モナドであることを使うと、モナド処理には色々なものが準備されているので、ありがたい。
    • モナドの処理を使うためにControl.Monadをimportして、同じモナド処理　>>= (\x -> choose 0.4 (x, 0) (x, 1))　を繰り返す関数 replicateMを適用する
    • どうしてこんなことができるかと言えば、Monadの連続処理には、結合法則様のものが保証されているから、処理順をやりやすいように切り替えてたたむことが可能
*Main Lib Numeric.Probability.Distribution> import Control.Monad
*Main Lib Numeric.Probability.Distribution Control.Monad> replicateM 2 (choose 0.4 0 1)
fromFreqs [([1,1],0.36),([0,1],0.24),([1,0],0.24),([0,0],0.16000000000000003)]
*Main Lib Numeric.Probability.Distribution Control.Monad> replicateM 3 (choose 0.4 0 1)
fromFreqs [([1,1,1],0.216),([0,1,1],0.144),([1,0,1],0.144),([1,1,0],0.144),([1,0,0],9.600000000000002e-2),([0,0,1],9.6e-2),([0,1,0],9.6e-2),([0,0,0],6.400000000000002e-2)]
  • ほらできた、ということらしい
  • 使用例
    • k <- [1..n] １からnまでのそれぞれの場合kについて、リストを作る
    • dという確率事象の単純繰り返しをk回行う
    • D.map maximum で、k回試行の中の最大値をとり、その値についての確率質量分布にする
    • そのD.expected (期待値を取る)
import Numeric.Probability.Distribution as D
let expmaxseq n d = [D.expected $ D.map maximum $ replicateM k d | k <- [1..n]]
expmaxseq 7 (D.uniform [1..6])
[3.5,4.472222222222223,4.958333333333335,5.2445987654321105,5.4309413580246915,5.560292352536542,5.65411736969124]
    • 少し変えてみる。これはただの目の総和の期待値
*Main Lib Numeric.Probability.Distribution Numeric.Probability.Shape Control.Monad D> expsumseq 7 (D.uniform [1..6])
[3.5,7.0,10.499999999999995,14.000000000000036,17.49999999999999,21.000000000000497,24.499999999998735]