疑似乱数生成

先日、springerの本がAmazonで80%オフになっていたので、おすすめされた本のいくつかを購入しました。これから読んだところの内容を少しずつまとめていきたいともいます。
初回はLCGについての実装です。検索したら割と記事はいろいろ出てきたのですが、LCGが生成する乱数には自己相関があるといった、話はでてこなかったので、そこら辺の実験をメインで行いました。なお、この記事は参考文献のComputational Physicsをもとに作成しています。

さて、タイトルにもある通りLCGは古い世代のポケモンの乱数生成に使われていたそうです。LCGのパラメータの調整がうまくいかず、いろいろ調べていたときに、急にポケモンの話がでてきたのでびっくりしました。ちなみに、その時に読んだポケモンの記事はこちらです。

では、前置きが少し長くなりましたがそろそろ本編スタートです。

$$
$$
$$
$$

筆者「乱数調整？なにそれ。我が家では、コードフリー。。。げふんげふん。」

Linear Congruentioal Generators

日本語だと線形合同法というみたいです。式は漸化式の形になっていて、パラメータと初期値を決めるとあとは順に計算していくだけです。

$$$$

I_{n+1} = [aI_{n} + b] \mod m

a は(奇数)の正整数
m は$2^{32}$ くらい大きな正整数
b は m と互いに素である正整数

$$$$

パラメータの決め方はいろいろあるみたいなので、気になるかたは是非調べてみてください。この記事には参考文献に記載されているものを軽くつまんで載せました。
この漸化式は $[0, m)$ の範囲の値を生成するものになります。 $[0,1)$ の範囲の値が欲しい場合は I_n を m で割ったものを使用します。但し、$[0,1)$ の範囲の値を取得する場合、下記のパラメータの条件で行うとうまく乱数が生成できなかったです。

さて、この式をもう少し詳しく見ていきましょう。mod m であることから、周期の最大値は m です。例えば、mod 3 を考えたときに、I_n の候補は 0 or 1 or 2 の3パターンになります。すなわち m 個ですね。ここで注意したいのは最大周期が m だからと言って必ず周期が最大周期になるとは限らないことです。

Q 初期値 = 3, a = 7, b = 3, m = 16 として下記のLCGの計算式に当てはめてみてください。周期はどうなるでしょうか？

[追記]
@fujitanozomu(藤田望) 様がパラメータについてご教示してくださいました。ありがとうございます。

短所
下位ビットのランダム性が低い。Mの値によっては、最下位に近いビットはほとんどランダムでなく、完全に規則性を持っていることすらある。たとえばMが偶数の場合（コンピュータでの実装が楽であるために、Mに2の冪を採用した場合はこれに当たる）、最下位ビットは、同じものが出続けるか、0と1が交互にでるかのどちらかである。すなわち、偶数ばかりが出る、奇数ばかりが出る、偶数と奇数が交互に出る、のどれかになるということである（最大周期ならば偶数と奇数が交互に出る）。

線形合同法 - Wikipedia より

では、ここらへんで実際に実装してみましょう！

実装！

やってきました。Gistの中身読むだけで、この記事は読まなくていいかもしれないです。

mod 部分の計算が漸化式に当たる部分ですね。

lcg.jl

function LCG(initval, a, b, m, n=1)
    res = zeros(Int, n)
    res[1] = mod(a*initval + b, m)
    n == 1 && return res
    for i in 2:n
        res[i] = mod(a * res[i-1] + b, m)
    end
    res
end

乱数を生成させてみようか

実際に乱数を生成させてプロットしてみましょう。1枚目がLCGで、2枚目がJuliaの組み込み関数を使ったものです。ほぼいい感じじゃないでしょうか！

棒グラフでも見てみましょう。乱数なら一様に出現しているはず。結果は1枚目のLCGが少し1に近い値が少ないのかな？といった印象です。まぁ、許容範囲でしょう！

最後は冒頭にも話したように自己相関を見ていきます。まず、横軸はラグです。すなわち、いくつ隣の要素と比べているかということです。縦軸は相関関係なので $[-1,1]$ です。ラグ0は、同じ要素同士を比べているので相関は1になります。注意してください。

では実際にグラフを見ましょう。1枚目のLCGのほうは自己相関がありそうですね。2枚目もラグ4が怪しいですが、全体的に0に近い値になっています。このように、LCGを単純に使うと自己相関が現れることがあります。

自己相関を消そう！

参考文献では、この自己相関を消すための簡単な方法も紹介されていたので、こちらを紹介して終わりたいと思います。

アルゴリズムの流れを説明すると、

$[0, 1)$ におけるLCGの乱数のリストを作る。但し、リストの長さは100よりも大きな素数とする
$y = LCG \in [0, 1) $ を計算する
$j = int(y * リストの長さ)$ を計算する。但し、intは繰り上げで丸めるとする。
j 番目のリストの値を出力とする。
出力後、j 番目のリストの値を $LCG \in [0, 1) $ で更新し、2に戻る。

となります。実装も一緒に見たほうがわかりやすいと思うので、対応する行に番号を振っておきました。

lcg_trick.jl

function LCG(trick::ErasingTrick, initval, a, b, m, n=1)
    res = zeros(Float64, n)
    rlist = LCG(Float64, initval, a, b, m, trick.Z) # 1
    res[1] = calctrick(trick, (a, b, m), rlist)    
    n == 1 && return res
    for i in 2:n
        res[i] = calctrick(trick, (a, b, m), rlist) # 2 ~ 5
    end
    res
end

function calctrick(trick::ErasingTrick, param_abm, rlist)
    y = LCG(Float64, param_abm...)[begin]          # 2
    j = 1 + Int(floor(y*trick.Z))                  # 3
    y = rlist[j]                                   # 4
    rlist[j] = LCG(Float64, param_abm...)[begin]   # 5
    return y
end

これらを使ってプロットしてみた結果がこちらです。重要なのは3枚目です。見事に自己相関が消えていることがわかります。

まとめ

実際に乱数を使用するアプリケーションによっては自己相関なんか気にならないよといったことがあるかもしれません。その場合は、LCGは実装も簡単ですしメモリも少量で済みます。何事も使いようですね。

参考文献

Q の解答

周期は 4 となります。実装したLCGで確かめてみると 8 → 11 → 0 → 3 とループしているのがわかります。