線形合同法をJuliaで実装する

線形合同法をJuliaで実装する

最近「ゼロからできるMCMC」を読んでいます。Chapter2の乱数と擬似乱数の違い(2.1.3)で、擬似乱数の一例として、線形合同法が登場します。

実際にJuliaで実装してみたので、記録として残しておきます。

線形合同法とは

疑似乱数を発生させる方法の一つに、線形合同法があります。以下の漸化式で与えられます。

$$
x_{n+1} = ax_{n}+ b \pmod M
$$

1つ目に、初項となる整数 を決定し、上の漸化式に投入します。modの計算をして出てきた値が1つめの乱数値になります。2つ目は、その乱数値を上の漸化式に投入し、modの計算をすることで得られます。3つ目以降も同様です。

線形合同法を実装する

実装してみるとこんな感じです。漸化式とmodの関数を定義し、それをもとにfor文で計算するということを行っています。今回は、過去の値をチェックしたかったので、リストに入れています。

## 初期設定(Wikipediaより)
a = 3
b = 5
M = 13

# 関数を定義する
function f(x)
    return  mod(a*x + b, M) 
end

# 実行する
# リストの生成
list1 = [8]

# 繰り返して、list1に乱数列を入れる
for i in 1:100
    push!(list1,  f(list1[i]))
end
ja.wikipedia.org