最近「ゼロからできる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