LoginSignup
0
0

More than 5 years have passed since last update.

VBAHaskellの紹介 その25(乱数生成:メルセンヌ・ツイスタ mt19937)

Last updated at Posted at 2015-10-14

(2016/1/5) 新しく追加した discrete_dist と random_shuffle の説明を追記。


VBAHaskellに乱数列を生成するモジュールを追加した。
擬似乱数列生成器の1つであるメルセンヌ・ツイスタ、正確に言うとC++11の標準ライブラリstd::default_random_engine の実装として今の Visual C++ が採用しているmt19937を使った。
VBAのRnd関数よりも高品質な乱数を高速に得ることを目的としている。

新規ファイルとしてC++側に VBA_Random.cpp、VBA側に misc_random.bas を追加した。
今のところ実装しているのは乱数ジェネレータと5つの分布と2つのユーティリティだけだ。

関数名 機能 引数 返り値
seed_Engine 乱数ジェネレータを初期化 シード値整数(Optional) シード値
uniform_int_dist 一様整数乱数の生成 個数、from、to 乱数値(スカラ/配列)
uniform_real_dist 一様実数乱数の生成 個数、from、to 乱数値(スカラ/配列)
normal_dist 正規分布乱数の生成 個数、平均、標準偏差 乱数値(スカラ/配列)
bernoulli_dist ベルヌーイ乱数の生成 個数、true確率 乱数値(スカラ/配列)
discrete_dist 離散分布の生成 個数、発生比率配列 乱数値(スカラ/配列)
random_iota iotaのランダム版 スタート、エンド ランダムに並んだiota
random_shuffle 配列のランダムな並べ替え 配列 並べ替え後の配列

分布関数は引数の「個数」が1以上の時は配列を返し、1未満の場合は単一の値を返すようにしている。疑似乱数なので seed_Engine に同じシード値を代入すると同じ乱数列を再現することもできる。通常はシード値を省略してデフォルトの(その度に違う)値を使えばよい。
random_iota はスタートからエンドまでの整数がランダムに並んだ配列を返す。
random_shuffle は引数に渡した配列の要素をランダムに並べ替えた配列を返す。

サンプルプログラム

?seed_Engine()     ' デフォルトの初期化(std::random_deviceを使用)
 1889001354        ' 結果的に使われたシード値
N = 1000000
m = normal_dist(N, 0, 1)     ' 標準正規分布列を 1,000,000 個生成(けっこう速い!)
' 平均値
mean = foldl1(p_plus, m) / N 
' 標準偏差
stddev = (foldl1(p_plus, mapF(p_mult, mapF(p_minus(, mean), m))) / N) ^ 0.5
' 累積密度
density = foldl1(p_plus, mapF(p_less(, 1.0), m)) / N
?mean  stddev  density
 2.75740534904507E-04  1.00073839053398  0.84098 

標準正規分布における 1σ(片側) は0.8413447...だが、100万個だとこんなものなのだろう。
乱数列の生成はあっけないほど速い。それに対してVBAHaskellで求めた平均値や標準偏差は非常に遅い。これらを求める関数は別途作ることになるかもしれない。

次にサイコロの出目を作ってみる。

?seed_Engine()     ' デフォルトの初期化(std::random_deviceを使用)
 1593330603        ' 結果的に使われたシード値
' サイコロを100万回振る
m = uniform_int_dist(1000000, 1, 6)
' 出目=1,2,3,4,5,6に等しいかそれぞれ判定
eachN = mapF(p_mapF(p_equal(yield_1), m), Array(1, 2, 3, 4, 5, 6))
' 2次元配列にバラして(unzip)横方向に加算する
printM foldl1(p_plus, unzip(eachN,2), 2)
  166501  166990  166108  167039  166455  166907
?1000000 / 6      ' 確率は1/6
 166666.666666667 

次は discrete_dist を使って離散分布を作り、モールス符号のような文字列を作成してみる。

' (0, 1, 2)が順に出現比率 3:1:2 で離散分布する長さ48の配列
m = discrete_dist(48, Array(3, 1, 2))
printM m
  2  0  0  0  1  2  2  2  0  2  2  2  0  1  2  0  1  2  0  2  1  2  0  0  1  0  2  0  2  2  0  0  0  2  0  2  0  0  1  0  2  0  1  1  2  2  0  0
' m を Array("・"," ","-") から要素を取得するインデックスとして使って
' 0 => "・",  1 => " ",   2 => "-"  にそれぞれマップする
s = mapF(p_getNth(, Array("・"," ","-")), m)
?Join(s, "")
-・・・ ---・---・ -・ -・- -・・ ・-・--・・・-・-・・ ・-・  --・・

次は random_shuffle を使って配列をランダムにシャッフルする。

m = iota(1, 20)
printM m
  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20
m = random_shuffle(m)
printM m
  13  20  4  1  5  6  12  14  16  17  18  7  10  15  11  19  8  3  9  2

VBAHaskellの紹介 その24(yield式の導入~lambdaExprの廃止)
VBAHaskellの紹介 その1 (最初はmapF)
GitHub VBAHaskell

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0