LDS(2)
http://d.hatena.ne.jp/m-a-o/20091212#p3
の続き。
前回のあらすじ。ある領域上で定義された関数の積分を計算するには、その領域で一様分布する数列に対する関数値の平均を取ればいい。一様分布する数列としては乱数がよく使われてモンテカルロ法として知られている。けども、乱数でなくても、一様分布する数列は存在する。存在するだけでなく、乱数より"ずっと速く一様分布に収束する"数列が存在してLow-discrepancy sequence(LDS)と呼ばれている。LDSを使ってモンテカルロ積分は、乱数を使ったモンテカルロ積分よりも収束が速くなる(Koksma-Hlawkaの定理)。
LDSを使ったモンテカルロ法は、準モンテカルロ法と呼ばれるらしい。けど、モンテカルロ法を最初に行ったvon Neumann-Ulamが使った数列は、そもそも乱数ではなかったので、全部モンテカルロでいいじゃんとか思う
積分領域が一次元の場合は、被積分関数のクラスに応じて、よい方法が沢山あるので、よっぽど病的な関数ばっかを扱うとかでない限り、ありがたみはそこまでないと思う。ので、必然的に高次元LDSを考える。高次元LDSを作る方法としては、「互いに相関のない」1次元LDSを並べるというのが、一つある。この方向で、van der Corput列を並べたのがHalton列。有理数体上一次独立な回転角を持つ無理数回転を並べたものを、Kronecker列とか呼ぶらしい。
実用上は、Faure列とかSobol列とかいうものがよく使われるらしい。Owenという人が考えたScramblingという操作を行ったものを、一般化Faure列とか一般化Sobol列とかいうらしいけど、これは、Faure列とか、Sobol列とかは、隣の点は近くにあることが多いので、一般には積分が収束していく時の振動幅が大きくなったり偏りが生じがちだけど、適当に順番入れ替えることで、そこらへんを改善できるという感じなのだと思う。ユーザ的には、一般化Faure列使っとけばいいんじゃね？という感じっぽい。
まあ、scramblingはいいとして、Faure列の定義は分かりやすい。Faure列がLDSであることを示すには、Niederreiterが導入した(t,s)-sequenceの概念を知っていれば比較的簡単(むしろ、(t,s)-sequenceの定義は、Faure列とか見て思いついたんだと思う)。
(t,s)-sequenceの定義は何も知らないで見たら意味不明だけど、要は、一様分布するということは、ある領域に含まれる点の数がその領域の体積に比例していればよいということを言ってるに過ぎない。とはいえ、有限個の点を取ってきた場合、適当な領域を取ったら、1個も点が含まれていないということは起こるので、代わりに基本区間というものを考える。van der Corput列とかの時と同じ意味でbを基数とすると、基本区間は、超直方体(?)で、各辺の座標軸への射影が、区間[a/b^j,(a+1)/b^j]になってるようなもの。で、最初のb^m個の点のうち、この基本区間に入っている点の数が、基本区間の体積に比例していればいいよねとかいうのが、アイデア。
(t,s)-sequence自体は、数列の性質を述べたに過ぎないので、その性質を満たす数列がどのくらい存在するかとか、どうやって構成するかとかは別問題だけど、Niederreiterは、(t,s)-sequenceを構成する体系的な方法を与えている。Faure列とかSobol列も(t,s)-sequenceで、Niederreiterの構成の特殊な場合とみなせる。
試しに7次元領域上の関数
(1*2*3*4*5*6*7)*cos(x1)*cos(2*x2)*cos(3*x3)*cos(4*x4)*cos(5*x5)*cos(6*x6)*cos(7*x7)
の数値積分を、乱数・Halton列・Faure列・Kronecker列で行って、積分誤差を比較してみた。Nは、サンプル点の数で、結果。
N = 1000 乱数: 9.9055483122 Halton: 9.15376775417 Faure: 5.72693149852 Kronecker: 3.04231847824 N = 10000 乱数: 6.46412740778 Halton: 1.11944591219 Faure: 4.12439321412 Kronecker: 1.70737628191 N = 100000 乱数: 0.21314196244 Halton: 0.13139014844 Faure: 0.189523901958 Kronecker: 0.245784242938 N = 500000 乱数: 0.478886945071 Halton: 0.34728002196 Faure: 0.0320830069256 Kronecker: 0.0997440714485
被積分関数がKronecker列に有利だったり不利だったりするんじゃないかという疑惑があるが、この程度の次元数だと、乱数も頑張ってるし、あんまりどの方法でも差がでない。あと、Faure列とかHalton列の生成が結構遅い(実装のせいだけど
以下、実装。
from random import Random import time from math import acos,pi #一次元LDS 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, _x): self.prev = 1 self.x = _x self.cur = _x self.mirror = False self.reflect = False self.theta = acos(_x)/pi def __iter__(self): return self def next(self): val = acos(self.cur)/pi if self.mirror: ret = 1-val else: ret = val #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 #高次元LDS #pure random number class RandomSequence: def __init__(self , dimension): self.generators = [Random(time.time()) for i in xrange(dimension)] def __iter__(self): return self def next(self): return tuple([g.random() for g in self.generators]) #Halton sequence class HaltonSequence: def __init__(self , dimension): p = 2 primes = [] while True: if(len(primes)==dimension):break for i in xrange(2,p): if p%i==0:break else: primes.append(p) p += 1 self.generators = [vanderCorputSequence(p) for p in primes] def __iter__(self): return self def next(self): return tuple([s.next() for s in self.generators]) #Kronecker class KroneckerSequence: def __init__(self , dimension): self.generators = [WeylSequence(0.01-float(2*i+1)/(2*dimension+1)) for i in xrange(dimension)] def __iter__(self): return self def next(self): return tuple([s.next() for s in self.generators]) def combination(n,k): if n==0 or k==0 or n==k:return 1 else: return combination(n-1,k-1)+combination(n-1,k) #Faure class FaureSequence: def __init__(self , dimension): p = dimension+1 while True: for i in xrange(2,p): if p%i==0: p += 1 break else: self.base = p self.n = 0 self.dimension = dimension return None def __iter__(self): return self def __getitem__(self,n): cs = [] ret = [] n = self.n b = self.base dimension = self.dimension while n > 0: cs.append(n%b) n = n/b for i in xrange(dimension): val = 0.0 for a in cs[::-1]: val = val/b + a ret.append( val/b ) tmp = [i for i in cs] cs = [] for i in xrange(len(tmp)): a = 0 for j in xrange(i,len(tmp)): a += combination(j,i)*tmp[j] cs.append( a%b ) return tuple(ret) def next(self): self.n += 1 return self.__getitem__(self.n) def monteCarloIntegral(fun , seq , N): val = 0.0 for (n , v) in enumerate(seq): if n==N:break val += apply(fun , v) return val/N if __name__=="__main__": from math import exp,sin,cos def f(x1,x2,x3,x4,x5,x6,x7): return 5040*cos(x1)*cos(2*x2)*cos(3*x3)*cos(4*x4)*cos(5*x5)*cos(6*x6)*cos(7*x7) D = 7 exactval = sin(1)*sin(2)*sin(3)*sin(4)*sin(5)*sin(6)*sin(7) for N in [1000 , 10000 , 100000 , 500000]: print "N =",N seq = RandomSequence(D) print "乱数:",abs(exactval - monteCarloIntegral(f , seq , N)) seq = HaltonSequence(D) print "Halton:",abs(exactval - monteCarloIntegral(f , seq , N)) seq = FaureSequence(D) print "Faure:",abs(exactval - monteCarloIntegral(f , seq , N)) seq = KroneckerSequence(D) print "Kronecker:",abs(exactval - monteCarloIntegral(f , seq , N)) print ""
