ysakasin/README.md

## README.md

      
    Raw
  

              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)
	// 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(a67108864.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)
	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
	10: 8378
	20: 8056
	25: 7920

	3904.1: false
	4847.51: false
	5330.18: false
	# カイ 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)