相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法?

馬續(xù)波 劉佳藝 徐佳意 魯凡 陳義學(xué)

(華北電力大學(xué)核科學(xué)與工程學(xué)院,北京 102206)

相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法?

馬續(xù)波?劉佳藝 徐佳意 魯凡 陳義學(xué)

(華北電力大學(xué)核科學(xué)與工程學(xué)院,北京 102206)

(2017年4月17日收到;2017年5月16日收到修改稿)

當(dāng)采用蒙特卡羅方法對很多問題進(jìn)行研究時(shí),有時(shí)需要對多維相關(guān)隨機(jī)變量進(jìn)行抽樣.之前的研究表明:在協(xié)方差矩陣滿足正定條件時(shí),可以采用Cholesky分解方法產(chǎn)生多維相關(guān)隨機(jī)變量.本文首先對產(chǎn)生多維相關(guān)隨機(jī)變量的理論公式進(jìn)行了推導(dǎo),發(fā)現(xiàn)采用Cholesky分解并不是產(chǎn)生多維相關(guān)隨機(jī)變量的唯一方法,其他的矩陣分解方法只要能滿足協(xié)方差矩陣的分解條件,同樣可以用來產(chǎn)生多維相關(guān)隨機(jī)變量.同時(shí)給出了采用協(xié)方差矩陣、相對協(xié)方差矩陣和相關(guān)系數(shù)矩陣產(chǎn)生多維隨機(jī)變量的公式,以方便以后使用.在此基礎(chǔ)上,利用一個(gè)簡單測試題和Jacobi矩陣分解方法對上述理論進(jìn)行了驗(yàn)證.通過對大亞灣中微子能譜進(jìn)行抽樣分析,Jacobi矩陣分解和Cholesky矩陣分解結(jié)果一致.針對核工程中的不確定性分析常用的238U輻射俘獲截面協(xié)方差矩陣進(jìn)行分解時(shí),由于協(xié)方差矩陣的矩陣本征值有負(fù)值,導(dǎo)致很多矩陣分解方法無法使用,在引入置零修正以后發(fā)現(xiàn),與Cholesky對角線置零修正相比,Jacobi負(fù)本征值置零修正的誤差更小.

蒙特卡羅方法,相關(guān)變量,隨機(jī)數(shù)產(chǎn)生方法,抽樣

1 引 言

蒙特卡羅方法在很多領(lǐng)域得到了廣泛的應(yīng)用[1?7].由具有已知分布的隨機(jī)變量進(jìn)行抽樣的方法在蒙特卡羅方法中占有重要的地位,并且已得到廣泛的研究.對于多維隨機(jī)變量的蒙特卡羅模擬問題,通常假定各隨機(jī)變量之間是相互獨(dú)立的,由此可以根據(jù)每個(gè)變量的統(tǒng)計(jì)分布去獨(dú)立產(chǎn)生多變量的樣本.但在某些學(xué)科,比如核工程中的不確定性分析,不僅同一反應(yīng)類型每個(gè)能群與每個(gè)能群的多群截面是相關(guān)的,不同反應(yīng)類型同一反應(yīng)能群的截面也是相關(guān)的,故需要對多維具有相關(guān)性的隨機(jī)變量進(jìn)行模擬.這就需要產(chǎn)生具有相關(guān)系數(shù)矩陣的多個(gè)隨機(jī)數(shù)序列.之前的研究表明,基于協(xié)方差矩陣的Cholesky因子分解的線性變換方法被認(rèn)為是最好的一種方法[8],但Cholesky因子分解的線性變換方法要求協(xié)方差矩陣必須是正定的,而對于很多實(shí)際問題的協(xié)方差矩陣,正定條件并不一定能完全滿足,這時(shí)如果還采用Cholesky因子分解并做近似處理可能出現(xiàn)較大誤差.針對這樣的問題,本文提出了基于Jacobi矩陣分解法進(jìn)行相關(guān)隨機(jī)變量序列產(chǎn)生的方法,相比于Cholesky因子分解,Jacobi矩陣分解法不要求矩陣正定,只要滿足對稱即可,這樣就進(jìn)一步擴(kuò)展了相關(guān)隨機(jī)變量序列的產(chǎn)生方法.并且進(jìn)一步發(fā)現(xiàn),如果協(xié)方差矩陣本征值部分為負(fù)值時(shí),采用不同的置零修正方法,誤差有較大差別.

2 相關(guān)隨機(jī)變量模擬理論

由于正態(tài)分布具有在線性變換下保持分布規(guī)律不變的獨(dú)特性質(zhì),因此多維相關(guān)正態(tài)分布隨機(jī)變量抽樣序列可以通過方差矩陣進(jìn)行分解,然后對獨(dú)立正態(tài)分布抽樣序列進(jìn)行線性變換來產(chǎn)生[8].之前研究較多的是采用Cholesky分解,下面證明:不僅可以采用Cholesky分解,其他的分解方式只要能滿足把協(xié)方差矩陣A分解成矩陣Σ和ΣT的相乘的形式,即

則就可以利用矩陣Σ產(chǎn)生多維相關(guān)正態(tài)分布隨機(jī)變量.

設(shè)Xm×n=(X1,X2,...,Xn)m為需要產(chǎn)生的m組n維相關(guān)正態(tài)分布隨機(jī)變量序列,其均值為u=(u1,u2,...,un),協(xié)方差矩陣為A,表示如下:

對于核工程中截面的不確定性分析,Ai,j為分群截面之間的協(xié)方差矩陣.產(chǎn)生n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣序列Ym×n=(Y1,Y2,...,Yn),每個(gè)變量抽樣序列數(shù)m,并且Xn,Yn和Σ滿足公式:

并且X的協(xié)方差矩陣為

其中矩陣Σ需要通過對協(xié)方差矩陣A分解得到,

由(4),(5),(6)式可知,只要把協(xié)方差矩陣A分解成矩陣Σ和ΣT的相乘的形式,就可以利用(3)式來產(chǎn)生多維相關(guān)正態(tài)分布隨機(jī)變量序列.

如果問題給出的是相對協(xié)方差矩陣Ar,而不是直接給出協(xié)方差矩A,可以證明,同樣可以把相對協(xié)方差矩陣Ar進(jìn)行分解

其中Σr由相對協(xié)方差矩陣Ar分解得到.由于相對協(xié)方差矩陣矩陣元和協(xié)方差矩陣Aij滿足關(guān)系

則多維相關(guān)正態(tài)分布隨機(jī)變量序列Xi的產(chǎn)生公式為

其中I表示單位矩陣.

通常情況下,能群截面與能群截面之間的相關(guān)性處理采用協(xié)方差矩陣A表示外,有時(shí)還用相關(guān)系數(shù)矩陣表示,這時(shí)用σi表示Xi的標(biāo)準(zhǔn)偏差,計(jì)算公式如(10)式所示.

相關(guān)系數(shù)矩陣中的矩陣元系數(shù)

對相關(guān)系數(shù)矩陣進(jìn)行分解得

則多維相關(guān)正態(tài)分布隨機(jī)變量序列Xi的產(chǎn)生公式為

(13)式與文獻(xiàn)[8]中給出的(5)式的形式是一致的,不同的是文獻(xiàn)[8]中給出的(5)式針對Cholesky分解是成立的,并且也進(jìn)行了驗(yàn)證,但本文中的(13)式并沒有要求一定是Cholesky分解,只要相關(guān)系數(shù)矩陣R能滿足(12)式,則就可以利用(13)式產(chǎn)生多維相關(guān)正態(tài)分布隨機(jī)變量序列Xi.

與Cholesky分解相比,Jacobi矩陣分解并不要求協(xié)方差矩陣一定是正定的,只要求其具有對稱性即可.另外Jacobi方法還可以用來求解奇異矩陣,并且求解出的矩陣的奇異值擁有較高的相對精度,奇異向量的正交性好、有較強(qiáng)的數(shù)值穩(wěn)定性,并且算法實(shí)現(xiàn)簡單[9].Jacobi方法用于求實(shí)對稱陣的全部特征值、特征向量.對于實(shí)對稱矩陣A,必有正交矩陣U,使

其中D是對角陣,其主對角線元是A的特征值.

3 Jacobi方法產(chǎn)生相關(guān)隨機(jī)變量驗(yàn)證

3.1 相關(guān)隨機(jī)變量的模擬抽樣示例

采用文獻(xiàn)[8]中的例子,要求產(chǎn)生隨機(jī)變量X1,X2,X3,其中X1,X2∈ [?1,1],X3∈ [0,4].三個(gè)隨機(jī)變量具有如下相關(guān)系數(shù)矩陣R,

我們利用Fortran2003語言,分別采用Cholesky分解和Jacobi分解法對矩陣R進(jìn)行了計(jì)算,結(jié)果如表1和圖1.

表1 分別利用兩種方法產(chǎn)生的樣本計(jì)算的相關(guān)系數(shù)對比Tab le 1.CoMparing of two Methods to p roduce correlation coeffi cients for saMp le calcu lating.

圖1 X1,X2和X3兩兩散點(diǎn)圖 (a)X1和X3之間的散點(diǎn)圖;(b)X1和X2之間的散點(diǎn)圖;(c)X3和X2之間的散點(diǎn)圖;(d)變量X1的示例數(shù)滿足高斯分布Fig.1.ScattergraMof X1,X2and X3:(a)X1and X3;(b)X1and X2;(c)X1and X2;(d)variab le X1satisfies the Gauss d istribu tion.

3.2 大亞灣中微子能譜樣本產(chǎn)生驗(yàn)證

江門中微子實(shí)驗(yàn)的物理目標(biāo)是大約利用6年左右時(shí)間測量中微子的質(zhì)量順序.反應(yīng)堆中微子能譜是反應(yīng)堆中微子實(shí)驗(yàn)重要的輸入?yún)?shù),其抽樣方法和誤差計(jì)算方法也會對反應(yīng)堆中微子實(shí)驗(yàn)產(chǎn)生重要影響.文獻(xiàn)[10,11]利用世界上最大的中微子樣本測量了最精確的中微子能譜以及中微子能譜的協(xié)方差矩陣,如圖2所示.利用測量的中微子能譜和協(xié)方差矩陣,分別采用Cholesky分解和Jacobi分解兩種矩陣分解方法進(jìn)行抽樣,樣本總數(shù)為500.利用Jacobi方法產(chǎn)生的樣本如圖3(a)所示,計(jì)算得到每個(gè)能量區(qū)間的相對誤差如圖3(b)所示.由圖3可見,兩種計(jì)算方法計(jì)算得到的結(jié)果是一致的,誤差的最大不一致小于5%.

3.3 69群協(xié)方差數(shù)據(jù)庫分解

核工程設(shè)計(jì)中,通過對復(fù)雜核系統(tǒng)進(jìn)行不確定性分析,往往可以減小系統(tǒng)的近似程度,提高系統(tǒng)的經(jīng)濟(jì)性,近來收到國內(nèi)外學(xué)者進(jìn)行的大量研究[6,7,12?19].由于利用蒙特卡羅抽樣方法進(jìn)行不確定性分析具有多種優(yōu)點(diǎn)(比如可以同時(shí)研究多個(gè)輸出變量的不確定性,可以考慮高階不確定度影響等),蒙特卡羅抽樣方法越來越受到國內(nèi)外學(xué)者的重視.在利用蒙特卡羅抽樣方法進(jìn)行不確定性分析時(shí),其中關(guān)鍵的一步就是要對多群截面協(xié)方差數(shù)據(jù)進(jìn)行矩陣分解.多群截面的協(xié)方差矩陣一般由NJOY軟件[20]產(chǎn)生.壓水堆設(shè)計(jì)計(jì)算中,通常采用69群或者172群的能群結(jié)構(gòu),我們采用NJOY軟件產(chǎn)生了69群的235U和238U輻射俘獲截面的協(xié)方差矩陣,如圖4所示.238U輻射俘獲截面協(xié)方差矩陣的本征值如表2所列,可以看到,部分能群的本征值出現(xiàn)了負(fù)值.為了解決這個(gè)問題,國內(nèi)外的通用做法是對負(fù)本征值置零處理.本征值出現(xiàn)為負(fù)值的情況,主要原因是利用NJOY產(chǎn)生協(xié)方差矩陣過程中采用了很多近似,在多種近似的情況下,最終導(dǎo)致協(xié)方差矩陣的本征值出現(xiàn)負(fù)值.這種情況下,無法再用Cholesky分解方法對矩陣進(jìn)行分解,因?yàn)镃holesky分解必須要求協(xié)方差矩陣的本征值為大于等于零的數(shù)值.我們采用對本征值置零的方法后得到的新的本征值矩陣為D′,然后利用新的本征值矩陣和已經(jīng)得到的正交矩陣U求得新的矩陣A′,

圖2 (網(wǎng)刊彩色)大亞灣中微子實(shí)驗(yàn)測量的中微子能譜(a)和協(xié)方差矩陣(b)Fig.2.(color on line)Neutrino spectrum(a)and covarianceMatrix(b)Measured in Daya Bay neutrino experiMent.

圖3 (網(wǎng)刊彩色)利用Jacobi矩陣分解方法抽得樣本分布(a),分別利用Cholesky方法和Jacobi方法抽樣500計(jì)算得到每個(gè)能量區(qū)間的相對誤差(b)Fig.3.(color on line)SaMp le d istribution extracted by JacobiMatrix decoMposition Method(a),the relative error of each energy interval calcu lated by Cholesky and JacobiMethodsWith 500 saMp les(b).

求得矩陣A′和矩陣A的相對誤差分布如圖5所示,238U輻射俘獲截面協(xié)方差矩陣分解過程中,由于對本征值進(jìn)行置零修正而導(dǎo)致的協(xié)方差矩陣的變化最大約為1.2%.

如果還采用Cholesky分解,

其中L為下三角矩陣矩陣,P為只有對角線上有值的對角陣,如果對矩陣P上的負(fù)值置零,得到的新矩陣為P′,利用新矩陣求得的矩陣A′的相對誤差誤差分布如圖6所示.由圖6可見,得到的矩陣A′的相對誤差的最大值為153%,誤差要遠(yuǎn)大于采用本征值置零的方法.

圖4 (網(wǎng)刊彩色)238U(n,γ)反應(yīng)的協(xié)方差矩陣Fig.4.(color on line)covarianceMatrix of238U(n,γ)reaction.

表2 238U輻射俘獲截面協(xié)方差矩陣的本征值Tab le 2.Eigenvalue of238U radiative cap tu re cross section covariance Matrix.

圖5 238U輻射俘獲截面協(xié)方差矩陣采用Jacobi方法分解,由于對本征值進(jìn)行置零修正而導(dǎo)致的協(xié)方差矩陣的變化,最大為1.2%Fig.5.The covariance Matrix of238U radiative capture cross section is decoMposed by Jacobi Method With zero correction of the eigenvalue,the MaxiMuMrelative error is 1.2%.

圖6 (網(wǎng)刊彩色)238U輻射俘獲截面協(xié)方差矩陣采用Cholesky方法分解,由于對本征值進(jìn)行置零修正而導(dǎo)致的協(xié)方差矩陣的變化,最大為153%Fig.6. (color on line)The covariance Matrix of238U rad iative cap tu re cross section is decoMposed by Cholesky Method With zero correction of the eigenvalue,the MaxiMuMrelative error is 153%.

4 結(jié) 論

通過對產(chǎn)生多維相關(guān)隨機(jī)變量的理論公式進(jìn)行推導(dǎo),發(fā)現(xiàn)采用Cholesky分解并不是產(chǎn)生多維相關(guān)隨機(jī)變量的唯一方法,其他的矩陣分解方法只要能滿足協(xié)方差矩陣的分解條件,同樣可以用來產(chǎn)生多維相關(guān)隨機(jī)變量.同時(shí)給出了采用協(xié)方差矩陣、相對協(xié)方差矩陣和相關(guān)系數(shù)矩陣產(chǎn)生多維隨機(jī)變量的公式.在此基礎(chǔ)上,利用一個(gè)簡單測試題驗(yàn)證了Jacobi矩陣分解方法也同樣適用于相關(guān)隨機(jī)變量的產(chǎn)生.通過對大亞灣中微子能譜進(jìn)行抽樣分析,Jacobi矩陣分解和Cholesky矩陣分解計(jì)算得到的每個(gè)能量區(qū)間的相對誤差結(jié)果是一致的.針對核工程中的不確定性分析常用的238U輻射俘獲截面協(xié)方差矩陣進(jìn)行分解時(shí),由于協(xié)方差矩陣的矩陣本征值有負(fù)值,導(dǎo)致很多矩陣分解方法無法使用,在引入置零修正以后發(fā)現(xiàn),與Cholesky對角線置零修正相比,Jacobi負(fù)本征值置零修正的誤差更小.雖然238U輻射俘獲截面協(xié)方差矩陣采用Jacobi負(fù)本征值置零修正更好,但不能保證該結(jié)論適用于其他例子,因此還需要對以上現(xiàn)象的理論做深入研究,以便找到更好的修正方法.

[1]Pei L C 1989 CoMpu ter Stochastic SiMu lation(Changsha:Hunan Science and Technology Press)p1(in Chinese)[裴鹿成 1989計(jì)算機(jī)隨機(jī)模擬(長沙:湖南科學(xué)出版社)第1頁]

[2]Xu SY 2006Mon te Carlo Method and its Application in Nuclear Physics ExperiMen t(2nd Ed.)(Beijing:AtoMic Energy Press)p1(in Chinese)[許淑艷 2006蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用 (第二版)(北京:原子能出版社)第1頁]

[3]Zhu Y S 2016 Statistic Analysis in High Energy Physics(Beijing:Science Press)p1(in Chinese)[朱永生 2016高能物理實(shí)驗(yàn)統(tǒng)計(jì)分析(北京:科學(xué)出版社責(zé)任有限公司)第1頁]

[4]Landau D P,Binder K 2000 A Guide to Mon te Carlo SiMu lations in Statistical Physics(2nd Ed.)(NeWYork:CaMbridge University Press)p1

[5]Ferguson D M,SiepMann J I,Truh lar D G 1999 Mon te Carlo Methods in CheMical Physics(NeWYork:John Wiley&Sons,Inc.)p1

[6]MattheWR B 2011 Ph.D.Dissertation(HaMilton:McMaster university)

[7]Hu Z H,Ye T,Liu X G,Wang J 2017 Acta Phys.Sin.66 012801(in Chinese)[胡澤華,葉濤,劉雄國,王佳2017物理學(xué)報(bào)66 012801]

[8]Wen D Z,Zhuo R H,D ing D J,Zheng H,Cheng J,Li Z H 2012 Acta Phys.Sin.61 220204(in Chinese)[文德智,卓仁鴻,丁大杰,鄭慧,成晶,李正宏2012物理學(xué)報(bào)61 220204]

[9]Guo Q 2011 M.S.D issertation(Jiangsu:SoochoWUniversity)(in Chinese)[郭強(qiáng) 2001碩士學(xué)位論文(江蘇:蘇州大學(xué))]

[10]An F P,Balantekin A B,Band H R,et a l.2017 Chin.Phys.C 41 013002

[11]An F P,Balantekin A B,Band H R,et al.2016 Phys.Rev.Lett.116 061801

[12]Ivanov K,Av raMova M,KaMeroWS,Kodeli I,Sartori E,Ivanov E,Cabellos O 2013 BenchMarks for UncertaintyAnalysis in Modelling(UAM)for the Design,Operation and Safety Analysis of LWRs,VoluMe I:Specification and Support Data for Neutronics Cases(Phase I),NEA/NSC/DOC(2013)7 https://www.oecd-nea.org/science/docs/2013/nsc-doc2013-7 pd f

[13]YaMaMoto A,K inoshita K,Watanabe T,Endo T,KodaMa Y,Ohoka Y,Ushino T,Nagano H 2015 Nucl.Sci.Engineer.181 160

[14]Park H J,ShiMH J,K iMC H 2012 Sci.Technol.Nucl.Insta ll.2012 247

[15]Curtis E M2013 M.S.D issertation (HaMilton:McMaster University)

[16]Dan G C,Mihaela I B 2004 Nucl.Sci.Engineer.147 204

[17]ZwerMann W,Aures A,Gallner L,et al.2014 Nucl.Engineer.Techno l.46 3

[18]Wan C H,Cao L Z,Wu H C,Zu T J,Shen W2015 Atom.Energy Sci.Technol.49 11(in Chinese)[萬承輝,曹良志,吳宏春,祖鐵軍,沈偉2015原子能科學(xué)技術(shù)49 11]

[19]Wang X Z,Yu H,Wang WM,Hu Y,Yang X Y 2014 Atom.Energy Sci.Technol.48 9(in Chinese)[王新哲,喻宏,王文明,胡赟,楊曉燕2014原子能科學(xué)技術(shù)48 9]

[20]Macfarlanef R,Kah ler A 2010 Nuclear Data Sheets 111 2739

PACS:02.50.Ng DO I:10.7498/aps.66.160201

*Project supported by National Natural Science Foundation of China(G rant No.11390383)and FundaMental Research Funds for the Central Universities of China(G rant No.2015ZZD 12).

?Corresponding author.E-Mail:Maxb@ncepu.edu.cn

G eneration o f correlated pseudorandoMvaribales?

Ma Xu-Bo?Liu Jia-Yi Xu Jia-Yi Lu Fan Chen Yi-Xue
(School of Nuclear Science and Technolgy,North China E lectric Power University,Beijing 102206,China)

17 Ap ril 2017;revised Manuscrip t

16 May 2017)

When Monte Carlo method is used to study many probleMs,it is sometimes necessary to saMp le correlated pseudorandoMvariables.Previous studies have shown that the Cholesky decoMposition Method can be used to generate correlated pseudorandoMvariab les when the covariance Matrix satisfies the positive eigenvalue condition.However,some covariancematrices do not satisfy the condition.In this study,the theoretical formula for generating correlated pseudorandoMvariables is deduced,and it is found that Cholesky decoMposition is not the only way to generatemultidiMensional correlated pseudorandoMvariables.The other Matrix decoMposition Methods can be used to generate multidimensional relevant randoMvariables if the positive eigenvalue condition is satisfied.At the same time,we give the formula for generating themultidiMensional randoMvariab le by using the covarianceMatrix,the relative covariance Matrix and the correlation coeffi cientMatrix to facilitate the later use.In order to verify the above theory,a siMp le test exaMp leWith 3×3 relative covariancematrix is used,and it is found that the correlation coeffi cient results obtained by JacobiMethod are consistent With those froMthe Cholesky Method.The correlation coeffi cients areMore close to the real valuesWith increasing the saMp ling number.A fter that,the antineutrino energy spectra of Daya Bay are generated by using JacobiMatrix decoMposition and Cholesky Matrix decoMposition Method,and their relative errors of each energy bin are in good agreeMent,and the diff erences are less than 5.0%in alMost all the energy bins.The above two tests demonstrate that the theoretical formula for generating correlated pseudorandoMvariab les is corrected.Generating correlated pseudorandoMvariab les is used in nuclear energy to analyze the uncertainty of nuclear data library in reactor simulation,and Many codes have been developed,such as one-,two-and three-diMensional TSUNAMI,SCALE-SS,XSUSA,and SUACL.However,when themethod of generating correlated pseudorandoMvariables is used to decoMpose the238U radiation cross section covarianceMatrix,it is found that the negative eigenvalue appears and p revious study Method cannot be used.In order to deal With the238U radiation cross section covariance Matrix and other siMilar matrices,the zero correction is proposed.When the zero correction is used in Cholesky diagonal correction and Jacobi eigenvalue zero correction,it is found that Jacobi negative eigenvalue zero correction error is sMaller than that With Cholesky diagonal correction.In future,the theory about zero correction Will be studied and itWill focus on ascertaining which correction Method is better for the negative eigenvalueMatrix.

Moto Carlomethod,pseudorandoMnumbers,correlated randoMvariables,samp ling

10.7498/aps.66.160201

?國家自然科學(xué)基金(批準(zhǔn)號:11390383)和中央高校基本科研業(yè)務(wù)費(fèi)(批準(zhǔn)號:2015ZZD 12)資助的課題.

?通信作者.E-Mail:Maxb@ncepu.edu.cn

?2017中國物理學(xué)會C h inese P hysica l Society http://Wu lixb.iphy.ac.cn