相關變量隨機數序列產生方法?

馬續波 劉佳藝 徐佳意 魯凡 陳義學

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

相關變量隨機數序列產生方法?

馬續波?劉佳藝 徐佳意 魯凡 陳義學

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

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

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

蒙特卡羅方法,相關變量,隨機數產生方法,抽樣

1 引 言

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

2 相關隨機變量模擬理論

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

則就可以利用矩陣Σ產生多維相關正態分布隨機變量.

設Xm×n=(X1,X2,...,Xn)m為需要產生的m組n維相關正態分布隨機變量序列,其均值為u=(u1,u2,...,un),協方差矩陣為A,表示如下:

對于核工程中截面的不確定性分析,Ai,j為分群截面之間的協方差矩陣.產生n維獨立標準正態分布隨機變量抽樣序列Ym×n=(Y1,Y2,...,Yn),每個變量抽樣序列數m,并且Xn,Yn和Σ滿足公式:

并且X的協方差矩陣為

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

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

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

其中Σr由相對協方差矩陣Ar分解得到.由于相對協方差矩陣矩陣元和協方差矩陣Aij滿足關系

則多維相關正態分布隨機變量序列Xi的產生公式為

其中I表示單位矩陣.

通常情況下,能群截面與能群截面之間的相關性處理采用協方差矩陣A表示外,有時還用相關系數矩陣表示,這時用σi表示Xi的標準偏差,計算公式如(10)式所示.

相關系數矩陣中的矩陣元系數

對相關系數矩陣進行分解得

則多維相關正態分布隨機變量序列Xi的產生公式為

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

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

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

3 Jacobi方法產生相關隨機變量驗證

3.1 相關隨機變量的模擬抽樣示例

采用文獻[8]中的例子,要求產生隨機變量X1,X2,X3,其中X1,X2∈ [?1,1],X3∈ [0,4].三個隨機變量具有如下相關系數矩陣R,

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

表1 分別利用兩種方法產生的樣本計算的相關系數對比Tab le 1.CoMparing of two Methods to p roduce correlation coeffi cients for saMp le calcu lating.

圖1 X1,X2和X3兩兩散點圖 (a)X1和X3之間的散點圖;(b)X1和X2之間的散點圖;(c)X3和X2之間的散點圖;(d)變量X1的示例數滿足高斯分布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 大亞灣中微子能譜樣本產生驗證

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

3.3 69群協方差數據庫分解

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

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

圖3 (網刊彩色)利用Jacobi矩陣分解方法抽得樣本分布(a),分別利用Cholesky方法和Jacobi方法抽樣500計算得到每個能量區間的相對誤差(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輻射俘獲截面協方差矩陣分解過程中,由于對本征值進行置零修正而導致的協方差矩陣的變化最大約為1.2%.

如果還采用Cholesky分解,

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

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

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

圖5 238U輻射俘獲截面協方差矩陣采用Jacobi方法分解,由于對本征值進行置零修正而導致的協方差矩陣的變化,最大為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 (網刊彩色)238U輻射俘獲截面協方差矩陣采用Cholesky方法分解,由于對本征值進行置零修正而導致的協方差矩陣的變化,最大為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 結 論

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

[1]Pei L C 1989 CoMpu ter Stochastic SiMu lation(Changsha:Hunan Science and Technology Press)p1(in Chinese)[裴鹿成 1989計算機隨機模擬(長沙:湖南科學出版社)第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蒙特卡羅方法在實驗核物理中的應用 (第二版)(北京:原子能出版社)第1頁]

[3]Zhu Y S 2016 Statistic Analysis in High Energy Physics(Beijing:Science Press)p1(in Chinese)[朱永生 2016高能物理實驗統計分析(北京:科學出版社責任有限公司)第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物理學報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物理學報61 220204]

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

[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原子能科學技術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原子能科學技術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

?國家自然科學基金(批準號:11390383)和中央高校基本科研業務費(批準號:2015ZZD 12)資助的課題.

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

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