LDS(1)
超一様分布列/準乱数/Low-discrepancy sequence(LDS)というものを最近知った。low discrepancy sequenceよりもhigh uniformity sequenceとかいう名前のほうがいいと思うけど
金融工学の数値シミュレーションとかで使われてるらしいけど、動機としては経路積分の計算に使うと面白いんじゃないか的なところ。既に誰かやっている人はいそうだけど、化学周りだと、そもそも経路積分で計算すること自体マニアックなので、まあ調べてみる価値はあるかもしれない
LDSの発端は、円周上の無理数回転に関する、Weylの一様分布定理にあるらしい。モンテカルロ法では乱数を使うけれども、(測度に関して)一様に分布してくれるなら、別に乱数である必要はなく(例えば、台形則とかは等間隔に分割しているけども、これも分割数を増やせば、一様分布に近づく)、実際にLDSを使うと収束が速いらしい。
被積分関数について何の情報もなければ、台形則以上の積分法もないように見える。台形則の問題は、要するに数列じゃないというところ。N分割してダメだから(N+1)分割とかすると、点集合は取り直しになるので、被積分関数の計算を再度(N+1)回やらないといけない。これは結構まぬけ。で、誰でも思いつくのが、2分割してダメなら4分割、それでもダメなら8分割、etc...という風にしてけばいいんじゃね？という手段。
区間[0,1)上で考えるとして(普通LDSは単位区間or高次元の場合は超立方体内で考える)、数列にするには並べ方はいくつかある気がするけど(そして、どれも最終的に一様分布に近づくけど)、van der Corput列というのは、
1/2,1/4,3/4,1/8,5/8,3/8,7/8,.....
とかいう順で並べるらしい。別に2分割にこだわる理由もない(計算する時は、ビットシフトで計算できるとかあるけど)ので、一般のb分割の時にも、同じように定義できる。
乱数にしろ、van der Corputにしろ、十分沢山の点を取れば、一様分布に近づくことには変わりないけど、一様分布に収束する速さは異なる。一様性(からのずれ)の測り方として、star-discrepancyというものが定義されている(star discrepancyとは別に普通のdiscrepancyという量も存在するけど、面倒なので以下単にstar discrepancyをdiscrepancyと書く)。discrepancyは、任意の次元で定義できるけど、一次元の場合だと、数列の最初のN点の内、区間[0,a)に入るものの割合と、区間の幅aとの差の上界(0
from math import cos,acos,pi
from random import Random
import time
class RandomSequence:
def __init__(self):
self.rand = Random(time.time())
def __iter__(self):
return self
def next(self):
return self.rand.random()
class vanderCorputSequence:
def __init__(self , b):
self.base = b
self.n = 0
def __iter__(self):
return self
def __getitem__(self,n):
cs = []
ret = 0.0
n = self.n
b = self.base
while n > 0:
cs.append(n%b)
n = n/b
for a in cs[::-1]:
ret = ret/b+a
return ret/b
def next(self):
self.n += 1
return self.__getitem__(self.n)
class WeylSequence:
def __init__(self, _theta):
self.prev = 1
self.x = cos(_theta)
self.cur = self.x
self.mirror = False
self.reflect = False
self.theta = _theta/pi
assert(0 < self.theta and self.theta < 1)
def __iter__(self):
return self
def next(self):
val = acos(self.cur)/pi
ret = 1-val if self.mirror else val
#Chebyshev polynomial
next = 2*(self.x)*(self.cur) - (self.prev)
self.prev = self.cur
self.cur = next
if val < self.theta or 1-val < self.theta:
if self.reflect:
self.reflect = False
else:
self.mirror = not(self.mirror)
self.reflect = True
return ret
def monteCarloIntegral(fun , seq , N):
val = 0.0
for (n , v) in enumerate(seq):
if n==N:break
val += fun(v)
return val/N
if __name__=="__main__":
from math import exp
#integral of exp(-x) over [0,1]
exactval = 1-exp(-1.0)
ns = [10, 100 , 1000 , 10000 , 100000]
for N in ns:
names = ["random" ,
"Weyl" ,
"van der Corput(base 2)" ,
"van der Corput(base 3)"]
seqs = [RandomSequence() ,
WeylSequence(0.46) ,
vanderCorputSequence(2) ,
vanderCorputSequence(3)]
print "N = " , N
for (name,seq) in zip(names,seqs):
val = monteCarloIntegral(lambda x:exp(-x) , seq , N)
print abs(val-exactval) , nameWeylSequenceで与える_thetaは、コンピュータでは基本的に有理数しか扱えないので、(0,1)上の無理数thetaを直接指定する代わりに、pi*thetaを与える。つまり、_theta/piは無理数なので、これを回転角とする無理数回転を考えている
実行結果
N = 10 0.0771923782405 random 0.0560457382959 Weyl 0.0299213553304 van der Corput(base 2) 0.0343910479541 van der Corput(base 3) N = 100 0.0199984031739 random 0.00927106895376 Weyl 0.00599774041195 van der Corput(base 2) 0.00648812817281 van der Corput(base 3) N = 1000 0.000432703075172 random 0.000289339893333 Weyl 0.000684694029749 van der Corput(base 2) 0.00090853938505 van der Corput(base 3) N = 10000 0.00074993581975 random 0.000158183597009 Weyl 0.000104894241615 van der Corput(base 2) 0.000148584365591 van der Corput(base 3) N = 100000 0.000565642700084 random 3.99734242706e-06 Weyl 1.30808546088e-05 van der Corput(base 2) 1.63554735195e-05 van der Corput(base 3)
乱数はランダムなので、運がいいとLDSよりいい結果を与えることもあるが、大体はWeyl列やvan der Corput列より悪い結果となる
実用的な場面でよく使われているLDSは、一般化Faure列とか一般化Sobol列とか一般化Niederreiter列とか呼ばれるものたちらしい。Niederreiterは1987年に、(t,s)-sequenceというクラスのLDSを定義して、(t,s)-sequenceを組織的に構成する手順を与えた。Halton列やWeyl列を除けば、現在知られているLDSは、大抵(t,s)-sequenceになっている。このへんの話は、まだちゃんと理解できてない
