長周期乱数
独立な複数の系列への分割
別のサイトで説明しているL'Ecuyer の方法による乱数は、周期が約
と長い。乱数を用いるシミュレーションの場合、乱数生成に掛かる時間よりも、生成した乱数に基づいて行われる計算時間の方が長いことが多い。また、周期が長いと、乱数生成の系列を適当に分割して、独立した複数の乱数生成系列を容易に得ることができる。このことについては本ウェブサイトの後半で説明する。
L'Ecuyer の方法による乱数の周期の長さをコンピュータ上での実行時間の長さで調べるため以下のPythonプログラムを作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ CheckTime.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from datetime import datetime
from RNGen import RNp2to191 as
MyRN
# A random number
generator class
def get_MicroSec():
"""
Timer in microsecond
"""
now = datetime.now()
return now.microsecond +
1000000.0 * (now.second + 60.0 * (now.minute \
+ 60.0 * now.hour))
rn = MyRN()
t0 = get_MicroSec()
for v in range(2**20):
v = rn.uni()
t1 = get_MicroSec();
print("(2**20): ", t1 - t0,
"microseconds")
s = 1000000.0 * 60.0 * 60.0 * 24.0 *
365.0
v = ((t1 - t0) * (2**191) / (2**20)) /
s
print("(2**191): ", v,
"years")
w = ((t1 - t0) * (2**100) / (2**20)) /
s
print("(2**100): ", w,
"years")
u = ((t1 - t0) * (2**50) / (2**20)) /
s
print("(2**50): ", u,
"years")
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムでは、
回一様乱数rn.uni()を呼び出したときの時間をもとに、周期である回呼び出すときの時間、および
回呼び出すときの時間と
回呼び出すときの時間を求めている。実行結果は図１のようになる。
図１
回の呼び出しに約３秒掛かっている。このとき、回の呼び出しには約
年掛かるという結果である。１周期内の乱数すべて生成することは事実上不可能である。回で約
年、回で約94年である。回で約年掛かるので、１周期個の長さの系列を個ずつの部分系列に分けて使えば、独立な複数の系列として扱える。このとき、最大
個の独立な部分系列を得ることになる。
乱数生成のクラス型のアトリビュートjump(n)は、乱数生成を
個先まで飛ぶものである。いま、３個の乱数生成オブジェクトrn_1、rn_2、rn_3を用意して、図２のように個ずつシフトする。
図２
例として、次のようなコードを作成した。このシフトにより３個の乱数生成オブジェクトrn_1、rn_2、rn_3は、３個の独立な乱数生成オブジェクトとして用いることができる。
from RNGen import RNp2to191 as MyRN
rn_1 = MyRN()
rn_2 = MyRN()
rn_2.jump(100)
rn_3 = MyRN()
rn_3.jump(100)
rn_3.jump(100)
print("{0:<20}{1:20}{2:20}".format('rn_1',
'rn_2', 'rn_3'))
for i in range(10):
print("{0:<20.15f}{1:<20.15f}{2:<20.15f}".format(rn_1.uni(),
rn_2.uni(), rn_3.uni()))
上のコードを実行すると図３のような結果になる。
図３
本ウェブサイトのプログラムは圧縮ファイルRN_Prgs.zipにまとめた。
乱数生成クラスのパッケージファイルRNGen.pyは以下のようである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ RNGen.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
#
#
A progtram purely in Python
#
of
#
Pierre L'Ecuyer's algorithm of random number generation
#
#
Yasuharu Okamoto,2016.09
#
# cf.
R. L. Eubank and A. Kupresanin (2012). "Statistical computing in C++ and
R".
#
class RNp2to191:
c209 = 2**32 - 209
c22853 = 2**32 - 22853
def __init__(self,*t):
if
(len(t) == 0):
self.x10 = 3
self.x11 = 2
self.x12 = 1
self.x20 = 6
self.x21 = 5
self.x22 = 4
else:
if (len(t) != 6):
print("len(t) = %d"%(len(t)))
raise "Number of parameters Error !"
else:
self.x10 = t[0]
self.x11 = t[1]
self.x12 = t[2]
self.x20 = t[3]
self.x21 = t[4]
self.x22 = t[5]
#
# Uniform
distribution, 0 < uni() < 1
#
def uni(self):
tx1
= (1403580 * self.x11 - 810728 * self.x12) % RNp2to191.c209
self.x12 = self.x11
self.x11 = self.x10
self.x10 = tx1
tx2
= (527612 * self.x20 - 1370589 * self.x22) % RNp2to191.c22853
self.x22 = self.x21
self.x21 = self.x20
self.x20 = tx2
z =
(self.x10 - self.x20) % RNp2to191.c209
if (z > 0):
return z / (RNp2to191.c209 + 1.0)
else:
return RNp2t0191.c209 / (RNp2to191.c209 + 1.0)
#
# Jump over
2**n numbers to be generated
#
def jump(self, n):
if
(n > 0):
A = [[0, 1403580, -810728],
[1, 0, 0],
[0, 1, 0]]
B = [[527612, 0, -1370589],
[1, 0, 0],
[0, 1, 0]]
C = [[0, 0, 0],
[0, 0, 0],
[0, 0, 0]]
for i in range(0, n):
for j in range(0, 3):
for k in range(0, 3):
tv = 0
for h in range(0, 3):
tv = tv + A[j][h] * A[h][k]
C[j][k]
= tv
for j in range(0, 3):
for k in range(0, 3):
A[j][k] = C[j][k] % RNp2to191.c209
for j in range(0, 3):
for k in range(0, 3):
tv = 0
for h in range(0, 3):
tv = tv + B[j][h] * B[h][k]
C[j][k] = tv
for j in range(0, 3):
for k in range(0, 3):
B[j][k] = C[j][k] % RNp2to191.c22853
t0 = A[0][0] * self.x10 + A[0][1] * self.x11 + A[0][2] * self.x12
t1 = A[1][0] * self.x10 + A[1][1] * self.x11 + A[1][2] * self.x12
t2 = A[2][0] * self.x10 + A[2][1] * self.x11 + A[2][2] * self.x12
self.x10 = t0 % RNp2to191.c209
self.x11 = t1 % RNp2to191.c209
self.x12 = t2 % RNp2to191.c209
t0 = B[0][0] * self.x20 + B[0][1] * self.x21 + B[0][2] * self.x22
t1 = B[1][0] * self.x20 + B[1][1] * self.x21 + B[1][2] * self.x22
t2 = B[2][0] * self.x20 + B[2][1] * self.x21 + B[2][2] * self.x22
self.x20 = t0 % RNp2to191.c22853
self.x21 = t1 % RNp2to191.c22853
self.x22
= t2 % RNp2to191.c22853
#
#
Standard Normal Distribution
# Rejection
polar method for normal variate
#
def normal(self):
w =
0.0
v1 =
0.0
while True:
v1 = 2.0 * self.uni() - 1.0
v2 = 2.0 * self.uni() - 1.0
w = (v1 ** 2.0) + (v2 ** 2.0)
if ((w < 1.0) and (w > 0.0)):
break
c =
math.sqrt(-2.0 * math.log(w) / w)
return c * v1
#
# Normal
Distribution of mean m and standard
deviation s
#
def normalMS(self, m, s):
return s * self.normal() + m
#
# A Pair of
Independent Random Numbers of
#
the Standard Normal Distribution
#
def normalPair(self):
w =
0.0
v1 =
0.0
v2 =
0.0
while True:
v1 = 2.0 * self.uni() - 1.0
v2 = 2.0 * self.uni() - 1.0
w = (v1 ** 2.0) + (v2 ** 2.0)
if ((w < 1.0) and (w > 0.0)):
break
c =
math.sqrt(-2.0 * math.log(w) / w)
return c * v1, c * v2
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝