エレガントな解法、エレファントな解法 〜モンテカルロ法を添えて〜
コインを100回投げたうちに、表もしくは裏が連続して10回でる確率をモンテカルロ法で計算する。
なにはともかくC++で実装してみた。
template < typename Random >
bool coin_toss( const unsigned int try_count, const unsigned int length, Random & r )
{
unsigned int count{} ;
int prev{} ;
std::uniform_int_distribution<> d( 0, 1 ) ;
for ( unsigned int i = 0 ; i != try_count ; ++i )
{
int result = d(r) ;
if ( result == prev )
{
++count ;
if ( count == length )
return true ;
}
else
{
prev = result ;
count = 1 ;
}
}
return false ;
}
template < typename Random >
double monte_carlo( unsigned int try_count, Random & r )
{
unsigned int count{} ;
for ( unsigned int i = 0 ; i != try_count ; ++i )
{
if ( coin_toss( 100, 10, r ) )
++count ;
}
return double(count) / double(try_count) ;
}
int main()
{
std::array<unsigned int, sizeof(std::mt19937)/sizeof(unsigned int)> c ; ;
std::generate( std::begin(c), std::end(c), std::random_device{} ) ;
std::seed_seq s( std::begin(c), std::end(c) ) ;
std::mt19937 engine(s) ;
for ( unsigned int i = 100 ; i != 1000000 ; i *= 10 )
{
auto r = engine ;
std::cout << i << "\t : " << monte_carlo( i, r ) * 100.0 << "%\n" ;
}
}
C++11から高度な乱数ライブラリが使えるようになったのだが、現在、実行時に適切に初期化する方法に面倒なボイラープレートコードが必要だ。実行時に乱数で初期化してくれる機能が現在提案されている。
ところで、最近Haskellを学んでいるので、このコードをHaskellでも書いてみようと思った。
まずは乱数だ。コインの表裏はBoolで表現するとして、Haskellらしく無限の乱数リストを生成しよう。
coin_seq :: (RandomGen g) => g -> [Bool]
coin_seq gen = randoms gen
あとはこれをtake 100して一回分の試行とする。1000回試行したければ、drop 100したのを再帰的にtake 100すればよい。しかし、もっといい方法がある。要素数100の[Bool]が入ったリスト、つまり[[Bool]]を無限に作ればいい。そのリストをtake 1000すれば1000回の試行分になる。
split_n :: Int -> [a] -> [[a]]
split_n n seq = take n seq : split_n n (drop n seq)
さて、seq_n = take 100 $ split_n 100 (coin_seq gen)のようなリストseq_n :: [[Bool]]は作れた。あとは[Bool]に連続した要素がn個あるかどうかを調べるcoin_toss nと、それをリストのすべての要素に対して行うmonte_carloを書けばいいだけだ。
coin_tossが結果をBoolで返すとすると、とりあえずseq_nにcoin_tossをmapしてTrueだけfilterして、その結果のリストの要素数を数えれば、成功した試行数がわかる。あとは試行数で割って確率を求めればいいだけだ。
monte_carlo :: Int -> [Bool] -> Double
monte_carlo try_count seq = ((fromIntegral n) / (fromIntegral try_count)) * 100.0
where
seq_n = take try_count (split_n 100 seq)
n = length . filter id $ map (coin_toss 10) seq_n
coin_toss len seqはInt -> [Bool] -> Boolな関数で、seqの中にlen個の連続したTrueもしくはFalseがあるかどうかを数えればよい。
いろいろ考えたが、まずgroup seqして連続した同じ要素を分割して[[Bool]]を作り、それぞれの[Bool]の要素数をいれたリスト[Int]をmapし、len以上の要素だけをfilterして、結果のリストに要素があればTrueを返す計算をすればいいのではないか。つまり、
has_contiguous_elems :: (Eq a) => Int -> [a] -> Bool
has_contiguous_elems len seq = not $ null $ filter (>= len) $ map length (group seq)
coin_toss :: Int -> [Bool] -> Bool
coint_toss _ [] = 0
coin_toss len seq = has_contiguous_elems len seq
となる。全体的にはこうだ。
import System.Random
import Data.List
coin_seq :: (RandomGen g) => g -> [Bool]
coin_seq gen = randoms gen
split_n :: Int -> [a] -> [[a]]
split_n n seq = take n seq : split_n n (drop n seq)
has_contiguous_elems :: (Eq a) => Int -> [a] -> Bool
has_contiguous_elems len seq = not $ null $ filter (>= len) $ map length (group seq)
coin_toss :: Int -> [Bool] -> Bool
coint_toss _ [] = 0
coin_toss len seq = has_contiguous_elems len seq
monte_carlo :: Int -> [Bool] -> Double
monte_carlo try_count seq = ((fromIntegral n) / (fromIntegral try_count)) * 100.0
where
seq_n = take try_count (split_n 100 seq)
n = length . filter id $ map (coin_toss 10) seq_n
main = do
gen <- getStdGen
let s = coin_seq gen
in mapM (\n -> putStrLn ((show n) ++ "\t : " ++ (show (monte_carlo n s)) ++ "%") )
[100, 1000, 10000, 100000]
このコードは動いた。ただし、最適化オプションを使っても、C++の10倍遅かった。
やはりlen個の連続した要素が存在するかどうかを調べるhas_contiguous_elemsが遅いのではないかと思い、これを再帰で実装することにした。
has_contiguous_elems' :: (Eq a) => Int -> [a] -> Bool
has_contiguous_elems' _ [] = False
has_contiguous_elems' len (x:xs) =
if count >= len
then True
else has_contiguous_elems' len $ dropWhile (== x) xs
where count = 1 + (length $ takeWhile (== x) xs)
この実装に差し替えたところ、実行速度はC++の9倍遅かった。
すると、他の部分もリストからリストを生成するという記述をやめて、再帰で実装すればもう少し早くなるのだろうか。
そして思うのは、同様のリストからリストを清々するような処理は、所詮リストのメモリーサイズがキャッシュに収まる程度なので、C++で書くとHaskellより早いのではないかとも思ったが、まだ試していない。
