乱数の検定

『アルゴリズムC〈第3巻〉グラフ・数理・トピックス』の35章 乱数を読みました。線形合同法，加法的合同法，乱数の検定について紹介がありました。乱数の検定について興味を持ったので，Rubyのrand関数についてχ2乗検定を行ってみました。

Rubyのソースファイルを見てみると，random.cの冒頭で，

This is based on trimmed version of MT19937. To get the original version,
contact http://www.math.sci.hiroshima-u.ac.jp/~m-mat/mt/emt.html.

とあるように，MT19937をベースにしているようです。上記のURLにアクセスしてみると，メルセンヌツイスターという方法のようです。

chi_square_test.rb

0以上100未満の乱数を1000個発生させ，χ2乗検定を行います。この操作を10回繰り返してみました。

# -*- coding: utf-8 -*-

# χ2乗検定(n:乱数の個数, r:0以上r未満の乱数)
def chisquare(n, r)
  time = Time.now
  srand(time.sec ^ time.usec ^ Process.pid)
  f = Array.new(r) { 0 }
  n.times do
    val = rand(r)
    f[val] += 1
  end
  t = f.map! { |x| x ** 2 }.inject(0) { |result, item| result + item }
  return ((r.to_f * t / n) - n)
end

10.times do
  # 0以上100未満の乱数を1000個発生させ，χ2乗検定を行う。
  printf "χ^2 = %6.1f\n", chisquare(1000, 100)
end

実行結果

実行結果は，以下のようになりました。
χ2乗の統計量がr(今回は100)に近ければランダムで，離れていればランダムではないのですが，『アルゴリズムC〈第3巻〉グラフ・数理・トピックス』によると，線形合同法では，100.8，加法的合同法では105.4であることと比較すると，以下の結果はちょっと疑問の残るものとなりました…。

$ ruby -v
ruby 1.9.1p378 (2010-01-10 revision 26273) [i386-darwin10.3.0]
$ ruby chi_square_test.rb
χ^2 = 117.6
χ^2 = 115.6
χ^2 = 113.0
χ^2 = 95.2
χ^2 = 95.0
χ^2 = 100.4
χ^2 = 87.8
χ^2 = 104.0
χ^2 = 86.0
χ^2 = 91.0