確率統計メルセンヌ・ツイスタ; Mersenne Twister
概要
メルセンヌ・ツイスタ[1]は、科学研究に広く使われている乱数生成法の一つ。1996年に松本眞氏、西村拓士氏によって発表された。
ソースコードは以下の公式ページに公開されている。
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html
特徴
・長い周期
出力列が一巡するまで\(2^{19937}-1\)と非常に長い周期を有する。実行時間の観点から周期をもつことを無視できるほど長い。
・高次元均等分布
メルセンヌ・ツイスタは623次元均等分布である。\(N\)次元均等分布とは、数列のうち隣り合う\(N+1\)個の数を抜き出し、それを\(N+1\)次元上にプロットしていくと、 点群が幾つかの並行した超平面に並ぶが、\(N\)個の数を抜き出した場合は超平面に並ばず均等に分布している性質のことである。
線形合同法:2次元均等分布
メルセンヌ・ツイスタ:623次元均等分布
・ビット均等性
ビット均等性とはすべてのビットが統計的にランダムである性質のことである。線形合同法は下位から\(n\)ビット目は\(2^n\)の周期しか持たない。 2で剰余すると交互に0，1が出ることになる。 ただしUNIX以外のrand()は下位16ビットを切り捨てることで，2で剰余すると交互に0，1が出る問題は解消している(それでもビット均等性ではない)。
32bitのバイナリプロット 下位ビットから左に並んでいる。 左:線形合同法 右:メルセンヌ・ツイスタ
ソースコード
namespace ExRandom {
public class MT19937 {
UInt32 pos;
UInt32[] state = new UInt32[624];
/// <summary>コンストラクタ 時間による乱数種設定</summary>
public MT19937() {
Initialize();
}
/// <summary>コンストラクタ 乱数種指定</summary>
public MT19937(int seed) {
Initialize(seed);
}
/// <summary>コンストラクタ 乱数種指定配列版</summary>
public MT19937(int[] seeds) {
Initialize(seeds);
}
/// <summary>32bit符号なし整数の乱数値の生成</summary>
private UInt32 Generate() {
//BSDライセンスのためメルセンヌ・ツイスタ公式ページを参照してください
}
/// <summary>時間による乱数種設定</summary>
public void Initialize() {
Initialize(Environment.TickCount);
}
/// <summary>乱数種指定</summary>
public void Initialize(int seed) {
//BSDライセンスのためメルセンヌ・ツイスタ公式ページを参照してください
}
/// <summary>乱数種指定配列版</summary>
public void Initialize(int[] seeds) {
//BSDライセンスのためメルセンヌ・ツイスタ公式ページを参照してください
}
/// <summary>32bit符号なし整数の乱数値の生成</summary>
public UInt32 Next() {
return Generate();
}
/// <summary>32bit符号なし整数の乱数値配列の生成</summary>
public UInt32[] NextArray(int size) {
var array = new UInt32[size];
for(int i = 0; i < array.Length; i++) {
array[i] = Generate();
}
return array;
}
/// <summary>64bit符号なし整数の乱数値配列の生成</summary>
public UInt64 Next64() {
return ((UInt64)Generate()) << 32 | (UInt64)Generate();
}
/// <summary>二値乱数の生成 確率thrで1を返す</summary>
public bool NextBool(double thr = 0.5) {
return NextDouble_OpenInterval1() < thr;
}
/// <summary>閉区間[0,1]の乱数の生成</summary>
public double NextDouble() {
const double INV_2POW32_M1 = 1.0 / (double)UInt32.MaxValue;
return (double)Generate() * INV_2POW32_M1;
}
/// <summary>半開区間(0,1]の乱数の生成</summary>
public double NextDouble_OpenInterval0() {
const double INV_2POW32 = 1.0 / (((double)UInt32.MaxValue) + 1.0);
return (double)Generate() * INV_2POW32 + INV_2POW32;
}
/// <summary>半開区間[0,1)の乱数の生成</summary>
public double NextDouble_OpenInterval1() {
const double INV_2POW32 = 1.0 / (((double)UInt32.MaxValue) + 1.0);
return (double)Generate() * INV_2POW32;
}
/// <summary>開区間(0,1)の乱数の生成</summary>
public double NextDouble_OpenInterval01() {
const double INV_2POW32 = 1.0 / (((double)UInt32.MaxValue) + 1.0);
return (((double)Generate()) + 0.5) * INV_2POW32;
}
}
}
namespace ExRandom.Tests {
[TestClass()]
public class MT19937Tests {
[TestMethod()]
public void MT19937Test() {
int[] seeds = { 0x123, 0x234, 0x345, 0x456 };
MT19937 mt = new MT19937(seeds);
var output = mt.NextArray(1000);
Assert.AreEqual(output[0], 1067595299u);
Assert.AreEqual(output[1], 955945823u);
Assert.AreEqual(output[624], 3768408841u);
Assert.AreEqual(output[999], 3460025646u);
}
}
}
[1] “Mersenne twister: A 623-dimensionally equidistributed uniform pseudorandom number generator”, Makoto Matsumoto, Takuji Nishimura, ACM Transactions on Modeling and Computer Simulations, vol.8, issue.1, pp.3-30, January 1998
