キケン！乱数発生！
衝動的に乱数発生させたくなったので，確率統計の教科書を引っ張りだし，やってみました．
２回生の頃に，乱数を発生させるプログラムをCでやってたのでCのソースが残っていたので，アルゴリズムを考えずにすみました．
乱数を発生させる手続き
結構有名な乱数を発生させる手続きとして，以下のようなものがあります．
rn+1≡a * rn + c(mod m)
これを，中心極限定理，ボックス・ミュラー法で，いい感じに正規乱数にするよ(っていう課題が出たんだよ)．
ちなみに，今回はa = 1229，c = 351750，m = 1664501，r = 2314としました(a,b,mは最も一般的な数値，rは適当です)．
平均と標準偏差
結果を比較しないといけないから，まずは配列の要素の平均と標準偏差を求めるスクリプトを．
# avg.rb
class Array
def avg
sum = 0.0
self.each{|e|
sum += e
}
sum/self.size
end
def stddev
variance = 0.0 #分散
avg = self.avg
self.each{|e|
variance += (e - avg) * (e - avg) / self.size;
}
Math.sqrt(variance)
end
end
中心極限定理
発生した正規分布な乱数r0，r1，r2・・・は0以上でm以下の整数なので，xn = rn/mとすると，0≤x<1の少数からなる数列{x0,x1,・・・}ができます．このxを順番に12個足し合わせ，そこから6を引くと，なんだかいい感じに平均0，標準偏差1に近い，-6〜6までの正規乱数ができるんだってさ．
それを，Rubyでやってみる．
require 'avg'
MAX=100000
def randomizer
a = 1229
c = 351750
m = 1664501
r = 2314
duly = Array.new(MAX){0.0} #正規乱数
x = 0.0
duly.each_index{|i|
(0 ... 12).each{
r = ( r * a + c ) % m
x = r.to_f / m.to_f
duly[i] += x
}
duly[i] -= 6
}
duly
end
duly = randomizer
puts "avg:#{duly.avg}" #平均
puts "stddev:#{duly.stddev}" #標準偏差
=>
$ ruby randomizer1.rb
avg:-0.00179716379263214
stddev:0.998569517897709一般に，正規乱数は多ければ多いほど良いのですが，Rubyでこんな10万個も正規乱数を作ろうとすると，なかなか帰ってこなくなるので不安になります(実際，不安になりました)．Memoizeとかしようがないから，この辺はスクリプト言語の限界ですね．
ボックス・ミュラー法
もう一つの方法，ボックス・ミュラー法．正直な事言いますと，ボクは確率はもともと得意でないので，ボックス・ミュラー法はそんな悟るほど知っている，という訳ではないので，説明しかねます．がんばってググって頂くか，こちらを読んで頂くかして，理解を深めてくださいな．
require 'avg'
MAX = 100000
def randomizer
include Math
a = 1229
c = 351750
m = 1664501
r = 2314
duly = Array.new(MAX){0.0} #正規乱数
x = 0.0 #一様乱数
y = 0.0 #一様乱数
duly.each_index{|i|
next if i % 2 != 0 #iが奇数なら飛ばす
r = ( r * a + c ) % m
x = r.to_f / m.to_f
r = ( r * a + c ) % m
y = r.to_f / m.to_f
duly[i] = sqrt(-2.0 * log(x)) * cos(2.0 * PI * y)
duly[i + 1] = sqrt(-2.0 * log(x)) * sin(2.0 * PI * y)
}
duly
end
duly = randomizer
puts "avg:#{duly.avg}" #平均
puts "stddev:#{duly.stddev}" #標準偏差
=>
$ ruby randomizer2.rb
avg:0.00119758054523837
stddev:0.997899329784094こっちは，わりかし早いです．繰り返しが圧倒的に少ないですからね．
