OpelのMT19937の検定をしてみる

README.md

OpelのMT19937の検定をしてみる

参考

Algorithms with Python 番外編：擬似乱数の検定

検定内容

• 等確率性の検定 (exec_test01)
  • カイ二乗分布
  • 9500前後なら許容可能
• 系列相関係数による検定 (exex_test02)
  • trueならOK

ファイル

• opal_mt.js
  • OpalのMT19937実装を抜き出して検定を組み込んだもの
• ruby_mt.rb
  • RubyのRandomクラスの検定
• result_js.txt
  • opal_mt.jsの出力
• result_ruby.txt
  • ruby_mt.rbの出力
• result_transpiled_ruby.txt
  • ruby_mt.rbをOpalでトランスパイルして実行した出力

opal_mt.js

// https://github.com/opal/opal/blob/b1074a37150f9b0b53a35f0efa4a1cbddc02ba29/opal/corelib/random/mersenne_twister.rb
// node opal_mt.js
var mt19937 = (function() {
    /* Period parameters */
    var N = 624;
    var M = 397;
    var MATRIX_A = 0x9908b0df;      /* constant vector a */
    var UMASK = 0x80000000;         /* most significant w-r bits */
    var LMASK = 0x7fffffff;         /* least significant r bits */
    var MIXBITS = function(u,v) { return ( ((u) & UMASK) | ((v) & LMASK) ); };
    var TWIST = function(u,v) { return (MIXBITS((u),(v)) >>> 1) ^ ((v & 0x1) ? MATRIX_A : 0x0); };
    function init(s) {
      var mt = {left: 0, next: N, state: new Array(N)};
      init_genrand(mt, s);
      return mt;
    }
    /* initializes mt[N] with a seed */
    function init_genrand(mt, s) {
      var j, i;
      mt.state[0] = s >>> 0;
      for (j=1; j<N; j++) {
        mt.state[j] = (1812433253 * ((mt.state[j-1] ^ (mt.state[j-1] >> 30) >>> 0)) + j);
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
        /* In the previous versions, MSBs of the seed affect   */
        /* only MSBs of the array state[].                     */
        /* 2002/01/09 modified by Makoto Matsumoto             */
        mt.state[j] &= 0xffffffff;  /* for >32 bit machines */
      }
      mt.left = 1;
      mt.next = N;
    }
    /* generate N words at one time */
    function next_state(mt) {
      var p = 0, _p = mt.state;
      var j;
      mt.left = N;
      mt.next = 0;
      for (j=N-M+1; --j; p++)
        _p[p] = _p[p+(M)] ^ TWIST(_p[p+(0)], _p[p+(1)]);
      for (j=M; --j; p++)
        _p[p] = _p[p+(M-N)] ^ TWIST(_p[p+(0)], _p[p+(1)]);
      _p[p] = _p[p+(M-N)] ^ TWIST(_p[p+(0)], _p[0]);
    }
    /* generates a random number on [0,0xffffffff]-interval */
    function genrand_int32(mt) {
      /* mt must be initialized */
      var y;
      if (--mt.left <= 0) next_state(mt);
      y = mt.state[mt.next++];
      /* Tempering */
      y ^= (y >>> 11);
      y ^= (y << 7) & 0x9d2c5680;
      y ^= (y << 15) & 0xefc60000;
      y ^= (y >>> 18);
      return y >>> 0;
    }
    function int_pair_to_real_exclusive(a, b) {
      a >>>= 5;
      b >>>= 6;
      return(a*67108864.0+b)*(1.0/9007199254740992.0);
    }
    // generates a random number on [0,1) with 53-bit resolution
    function genrand_real(mt) {
      /* mt must be initialized */
      var a = genrand_int32(mt), b = genrand_int32(mt);
      return int_pair_to_real_exclusive(a, b);
    }
    return { genrand_real: genrand_real, init: init, genrand_int32: genrand_int32 };
  })()



// カイ 2 乗値を求める
function get_chi2(count, F) {
  let chi2 = 0.0;
  count.forEach(x => chi2 += (x - F) ** 2 / F);

  return chi2;
}

// 度数検定
function test01(prg, n, k) {
  let count = new Array(k)
  count.fill(0)

  for (let i = 0; i < n; i++) {
    count[parseInt(mt19937.genrand_real(prg) * k)] += 1
  }

  return get_chi2(count, n / k)
}

const MAX_INT = Number.MAX_SAFE_INTEGER || Math.pow(2, 53) - 1

// test01 の実行
function exec_test01() {
  [[10, 16.919], [20, 30.144], [25, 36.415]].forEach((setting) => {
    const k = setting[0];
    const limit = setting[1];


    const seed = Math.round(Math.random() * MAX_INT);
    var mt = mt19937.init(seed)
    const n = k * 100   //1 つの区間に 100 個
    let i = 0
    for (let y = 0; y < 10000; y++) {
      if (test01(mt, n, k) < limit) {
        i += 1
      }

    }
    console.log(k, i)
  })
}

exec_test01()
console.log()



// カイ 2 乗分布表 (0.10 きざみ)
const chi2_9 = [
  14.684, 12.242, 10.656, 9.414, 8.343,
   7.357,  6.393,  5.380, 4.168, 0.000,
]

const chi2_19 = [
  27.203, 23.900, 21.689, 19.910, 18.338,
  16.850, 15.352, 13.716, 11.651,  0.000,
]

const chi2_24 = [
  33.196,  29.553,  27.096,  25.106,  23.337,
  21.652,  19.943,  18.062,  15.659,   0.000,
]

// test01 のカイ 2 乗検定
function chi2_test01(k, chi2_table){
  const seed = Math.round(Math.random() * MAX_INT);
  var prg = mt19937.init(seed)

  let count = new Array(10)
  count.fill(0)

  for (let p = 0; p < 2000; p++) {
    let chi2 = test01(prg, 100 * k, k)
    for (let i = 0; i < 10; i++) {
      if (chi2 >= chi2_table[i]) {
        count[i] += 1
        break
      }
    }
  }
  let chi2 = get_chi2(count, 200.0)
  console.log(chi2, chi2 < 16.919)
}

chi2_test01(10, chi2_9)
chi2_test01(20, chi2_19)
chi2_test01(25, chi2_24)

result_js.txt

10 9496
20 9505
25 9471

8.459999999999999 true
4.98 true
1.74 true


see http://www.nct9.ne.jp/m_hiroi/light/pystat04.html

result_ruby.txt

10: 9486
20: 9491
25: 9487

10.93: true
7.47: true
9.81: true

result_transpiled_ruby.txt

10: 8378
20: 8056
25: 7920

3904.1: false
4847.51: false
5330.18: false

ruby_mt.rb

# カイ 2 乗値を求める
def get_chi2(count, f)
  chi2 = 0.0
  count.each do |x|
      chi2 += (x - f) ** 2 / f
  end
  return chi2
end

# 度数検定
def test01(prg, n, k)
  count = Array.new(k, 0)
  n.times do
      count[(prg.rand() * k).to_i] += 1
  end
  return get_chi2(count, n.to_f / k)
end

# test01 の実行
def exec_test01()
  [[10, 16.919], [20, 30.144], [25, 36.415]].each do |k, limit|
      prg = Random.new
      n = k * 100    # 1 つの区間に 100 個
      i = 0
      10000.times do
        if test01(prg, n, k) < limit
          i += 1
        end
      end

      puts "#{k}: #{i}"
  end
end

exec_test01()
puts


# カイ 2 乗分布表 (0.10 きざみ)
chi2_9 = [
  14.684, 12.242, 10.656, 9.414, 8.343,
   7.357,  6.393,  5.380, 4.168, 0.000,
]

chi2_19 = [
  27.203, 23.900, 21.689, 19.910, 18.338,
  16.850, 15.352, 13.716, 11.651,  0.000,
]

chi2_24 = [
  33.196,  29.553,  27.096,  25.106,  23.337,
  21.652,  19.943,  18.062,  15.659,   0.000,
]

# test01 のカイ 2 乗検定
def chi2_test01(k, chi2_table)
  prg = Random.new
  count = Array.new(10, 0)
  2000.times do
    chi2 = test01(prg, 100 * k, k)
    10.times do |i|
      if chi2 >= chi2_table[i]
        count[i] += 1
        break
      end
    end
  end
  chi2 = get_chi2(count, 200.0)
  puts "#{chi2}: #{chi2 < 16.919}"
end

chi2_test01(10, chi2_9)
chi2_test01(20, chi2_19)
chi2_test01(25, chi2_24)