15. 乱数

15.1 概要

ランダムな値をとる数列を乱数といい, 例えば, 無作為にサンプルを選ぶような場合に使われます. 数値計算の分野で乱数が必要になる例としては, モンテカルロ法, 確率的現象のシミュレーションなどがあげられます.

乱数をコンピュータで生成するとはあるアルゴリズムに基づいたプログラムを使って乱数を得るということなので, 完全にはランダムではなく何らかの規則性を持つことになります. そのためこれを疑似乱数と呼ぶこともあります.

疑似乱数の利点としては, 必要なときにいつでも生成することができ, 結果を検証したいときに完全に再現することができることがあげられます. 単に疑似乱数というと一般的には一様分布に従う乱数のことを指します.

Excel では標準で疑似乱数を生成することができ, VBA では Rnd 関数と Randomize (初期化), ワークシート関数では RAND, RANDARRAY と RANDBETWEEN が提供されています.

15.2 乱数の生成方法

15.2.1 線形合同法

古くから広く使われている代表的な疑似乱数生成法です.

\(0\) と \(m\) の間の整数の初期値 \(X_0\) を選び, 次の漸化式により数列を生成します.
\[
X_{i+1} = (aX_i + c) \space mod \space m
\] ただし, \(0 \leq a < m, 0 \leq c < m\) とします.

パラメータ \(a, c\) および \(m\) を適切に選ぶと乱数列が得られます. 乱数列を繰り返す周期や統計的性質 (ランダム性) はこれらのパラメータの選び方によって決まります.

同じ初期値 \(X_0\) からは同じ乱数列が得られます. 従って, 最大でも \(m\) の周期で同じ数列を繰り返すことがわかります. そのため \(m\) はできるだけ大きい値を選んだ方がよいことになります. 加えて, 漸化式中の \(mod \space m\) の計算をオーバーフローを利用して速く行うために, ワード長 32 ビットであれば \(m = 2^{32}\) とすることがよくあります.

実数の乱数が必要な場合には, 得られた整数乱数を浮動小数演算で \(m\) で割るなどして変換します.

VBA 関数 Rnd では \(m = 2^{24}, a = 16598013, c = 12820163\) が使われていて, 周期は \(2^{24}\) とあまり長くありません. 関数値として \([0, 1)\) の単精度浮動小数を返します.

例として, 文献[1]で推奨されている \(m = 2^{32}, a = 69069, c = 1\) を使って実験してみます. 周期は \(2^{32}\) となります.
\[
X_{i+1} = (69069X_i + 1) \space mod \space 2^{32}
\] VB.NET プログラム例を以下に示します.

Class Test

Shared Function TRand(ByRef IX As Integer) As Integer
    IX = 69069 * IX + 1
    TRand = IX
End Function

Public Shared Sub Main(ByVal args() As String)
    Dim IX As Integer, I As Integer
    IX = 13
    For I = 1 To 10
        Console.WriteLine(TRand(IX))
    Next
End Sub

End Class

注 – 整数オーバーフローチェックを解除する必要があります (vbc -removeintchecks+ としてコンパイルします). Excal VBA ではその方法が不明だったので VB.Net を使用しました.

シードを 13 としたときには次のような乱数列が得られました.

 
897898
1887374819
-1755994680
784918825
-1813853482
-1145091233
1516383764
-1962284539
-1042831614
-735193445

 
線形合同法の欠点としては下位桁がランダムでないことがあげられます. 上の例でも奇数と偶数を繰り返している (最下位ビットが 1 と 0 を繰り返している) ことがわかります.

しかし, この場合でも上位桁のランダム性は良いので, 例えば上位 24 ビットだけを使うなどすれば良いランダム性が得られます. なお, 実数乱数に変換して使う場合には下位ビットの周期性はあまり問題になりません.

15.2.2 メルセンヌツイスター (MT)

1996年頃に日本で開発された比較的新しいものです. General Feedback Shift Register (GFSR) 方式を改良したアルゴリズムを使用しており, 文献[2]のホームページに詳しい説明とプログラムが掲載されています.

非常に長周期 \((2^{19937} – 1)\) かつ統計的性質も良い (よりランダムと言える) 乱数発生器です.

プログラミング言語 Python の乱数生成モジュールでもメルセンヌツイスターが採用されています.

15.3 XLPack を使った乱数の生成方法

Excel の VBA関数 Rnd やワークシート関数 RAND を使えば乱数を生成することができ, 一般的にはこれで問題ありません. しかし, 倍精度浮動小数が必要な場合や特に長周期ですぐれたランダム性を必要とする場合のためにメルセンヌツイスターの VBA 関数 InitGenrand, GenrandInt32, GenrandInt31, GenrandReal53 を提供しています.

InitGenrand は初期化を行います. この引数 (シードという) により生成する乱数列の始点が決まります. GenrandInt32, GenrandInt31, GenrandReal53 はそれぞれ符号付き 32 ビット整数, 符号なし 31 ビット整数, 区間 [0, 1) における 53 ビット実数の疑似乱数を生成します.

例題

はじめに InitGenrand による初期化を行い GenrandInt32 を呼び出して 32 ビット整数乱数を 10 個生成する.

VBA プログラム例を以下に示します.

Sub Start()
    Dim Seed As Long, I As Long
    Seed = 13
    Call InitGenrand(Seed)
    For I = 1 To 10
        Debug.Print GenrandInt32()
    Next
End Sub

次のような結果が得られました.

  
-954760878 
-1686456144 
 1020231754 
-603726320 
-754717978 
-459635870 
-147106060 
 769458329 
-117677332 
-1036873798 

参考文献