木拉提·吐爾德,胡錫健
(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,烏魯木齊 830046)
擬蒙特卡洛方法中Halton序列的隨機(jī)化及其改進(jìn)
木拉提·吐爾德,胡錫健
(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,烏魯木齊 830046)
Halton序列是簡(jiǎn)單并且容易實(shí)現(xiàn)的低差異序列。Wang和Hickrnell[1]提出的隨機(jī)化Halton序列在高維的情形下會(huì)產(chǎn)生高度相關(guān)情況,并且在高維情形下估計(jì)精度和收斂速度都很差。文章結(jié)合隨機(jī)化和置亂化方法改善了Halton序列這些缺點(diǎn)。模擬結(jié)果表明,改進(jìn)的隨機(jī)化Halton序列在精度和相關(guān)性方面都好于隨機(jī)化Halton序列。
—Halton序列;低差異序列;置亂化;隨機(jī)化;數(shù)字積分
自1960年Halton提出Halton序列以來,低差異序列在擬蒙特卡洛(QMC)方法中的應(yīng)用非常廣泛。Halton序列是簡(jiǎn)單的非常容易實(shí)現(xiàn)的低差異序列,但是最大的缺點(diǎn)是在高維的情況下會(huì)產(chǎn)生高度相關(guān)的情況,使得模擬的精度和收斂速度都非常差,為了改進(jìn)Halton序列這個(gè)缺點(diǎn),人們提出了一系列的改進(jìn)方法。1979年Braaten和weller[2]根據(jù)一維偏差率最小的原則提出了的置亂化(scrambling)Halton序列,1995年Owen[3]提出了線性置亂化Halton序列,2000年Wang和Hickrnell[1]提出的隨機(jī)化Halton序列,2004年Atanassovet[4]提出了可以使模擬的誤差邊界最小的Halton序列等。
本文擬結(jié)合隨機(jī)化方法和置亂化方法使對(duì)Wang和Hickrnell提出的隨機(jī)化Halton序列在高度相關(guān)和收斂速度方面進(jìn)行改進(jìn),并且通過實(shí)際模擬來說明改進(jìn)的隨機(jī)化Halton序列在高維的情況下估計(jì)精度和收斂速度都比文獻(xiàn)[1]的隨機(jī)化Halton序列要好。
首先簡(jiǎn)單地介紹經(jīng)典van der corput序列,一維的Halton序列實(shí)際上就是Van der corput序列。其次來介紹文獻(xiàn)[1]提出的隨機(jī)化Halton序列。
如果b≥2的整數(shù),對(duì)于任意的n≥0可以寫成以基為b的展開式
Halton序列是Van der corput序列的擴(kuò)展,一維的Halton序列是實(shí)際上就是Van der corput序列,多維的Halton序列則是由不同的素?cái)?shù)為基的Van der corput序列構(gòu)成??杀硎緸椋?/p>
隨機(jī)化Halton序列[1]要從φb(n)得到φb(n+1),我們只需要給φb(n)加1 b得到φb(n+1),這個(gè)加法不是普通的加減法這叫做向右進(jìn)位加法⊕(rightward carry addition)。
Van Neumann-Kakutani變換定義:若b≥2的整數(shù),對(duì)于x∈[0,1),我們把x寫成如下形式:
定義向右進(jìn)位加法(rightward carry addition)⊕,為:
因?yàn)镠alton序列隨維數(shù)的增加會(huì)產(chǎn)生高度相關(guān)的點(diǎn)列。為了消除這個(gè)高度相關(guān)現(xiàn)象,人們提出了置亂化(scrambling)方法。第一次正式的提出這個(gè)方法的是Braaten和weller[5]
像(1)式中基逆函數(shù)φb(n)一樣定義了置亂化的(scrambled)基逆函數(shù)Sb(n):
這里πb是在數(shù)集(0,1,…,b-1)中的置換,則確定的置亂化Halton序列記作:
改進(jìn)的目的是使隨機(jī)化Halton序列在高維的情形下相關(guān)性得到改善并且在使得隨機(jī)化Halton序列分布更均勻更隨機(jī)。改進(jìn)方法如下:
第一步:假設(shè)在一維的情況下,基為b。對(duì)于任意x0∈[0,1]s,則按照公式(4)得到T(x0)。
第二步:對(duì)得到T(x0)進(jìn)行置亂化變換,由公式(3)可知uk是取自于集合(0,1,…,b-1)的數(shù)。我們提出的置換如果uk≠0則把uk用b-uk來代替。即:
則對(duì)于任意xn∈[0,1]s一維改進(jìn)的隨機(jī)化halton序列是
如果對(duì)于任意xn∈[0,1]s基為b1,b2,…bs的s維多維改進(jìn)隨機(jī)化halton序列是
由定理1可知序列{T(xn)}∞0也是在[0,1]s上均勻分布的隨機(jī)向量。并且很容易證明序列Sn(T(xn))的誤差收斂率也是O(N-2(logN)2s)。
如圖1(b)所示,隨機(jī)化Halton序列在高維的情況下會(huì)產(chǎn)生相關(guān)性,這樣均勻性被破壞了并且估計(jì)方差也會(huì)很大。改進(jìn)隨機(jī)化避免了這種相關(guān)性,見圖1(a)。
圖1
下面考慮在s維[0,1]s下函數(shù) f(x)的多重積分
表1 隨機(jī)化Halton序列和改進(jìn)隨機(jī)化Halton序列模擬結(jié)果
表1是隨機(jī)化Halton序列和改進(jìn)隨機(jī)化Halton序列在維數(shù)s=10,s=20,s=30,s=40,s=50時(shí),并且所取的點(diǎn)數(shù)N=5000,N=10000,N=20000時(shí)多重積分 I(f)估計(jì)值和估計(jì)的樣本方差。
隨機(jī)化halton序列記作:rhalton,改進(jìn)隨機(jī)化halton序列記作irhalton。
從表1可以看出,在維數(shù)從10維到50維變化時(shí)改進(jìn)隨機(jī)化Halton序列的估計(jì)樣本標(biāo)準(zhǔn)差比隨機(jī)化Halton序列要小,并且40維以后隨機(jī)化Halton序列估計(jì)值明顯比改進(jìn)隨機(jī)化Halton序列要差。這時(shí)兩種方法都要求計(jì)算的函數(shù)值個(gè)數(shù)N要足夠大,才能達(dá)到滿意的精度。在估計(jì)的同時(shí)發(fā)現(xiàn)隨著計(jì)算函數(shù)值的個(gè)數(shù)N的增大隨機(jī)化Halton序列計(jì)算估計(jì)值的時(shí)間也非常長(zhǎng),而改進(jìn)隨機(jī)化Halton序列則比隨機(jī)化要短很多。
[1]Wang,X.,Hickernell,F.J.Randomized Halton Sequences[J].Math.Comput.Modeling,2000(,32).
[2]Braaten,E.,Weller,G.An Improved Llow-Discrepancy Sequence for Multidimensional Quasi-Monte Carlo Integrationp[J].Journal of Com?putational Physics,1979(,33).
[3]Owen,A.B.Randomly Permuted(t,m,s)-nets and(t,s)-sequences.In Monte Carlo and Quasi-Monte Carlo Methods in Scientific Comput?ing,H.Niederreiter and P.J.-S.Shiue,Eds.Number 106 in Lecture Notes in Statistics[M].New York:Springer-Verlag,1995.
[4]Atanassov,E.,Durchova,M.Generating and Testing the Modified Halton Sequences[C].In Fifth International Conference on Numerical Methods and Applications,Borovets Springer-Verlag,Ed.Lecture Notes in Computer Science,2002.
[5]Wang,X.,Sloan,I.H.Why are High-dimensional Finance Problems of low Effective Dimension?[J].SIAM Journal on Scientific Computing,2005(,27).
O211
A
1002-6487(2012)24-0015-03
新疆大學(xué)科學(xué)基金資助項(xiàng)目(07020428008)
木拉提·吐爾德(1985-),男,新疆人,碩士,講師,研究方向:概率統(tǒng)計(jì)。
(責(zé)任編輯/亦 民)