（计算数学专业论文）关于蒙特卡罗及拟蒙特卡罗方法的若干研究.pdf

摘要 此论文主要阐述荣特卡罗和拟蒙特卡罗方法在积分、优化和模拟方面 的应用的若干主题。 第1 章和第2 章是关于蒙特卡罗和拟蒙特卡罗方法的预备知识。第1 章介 绍了蒙特卡罗和拟蒙特卡罗积分的误差估计并阐述了拟蒙特卡罗方法的优 势，同时介绍了拟蒙特卡罗的标准优化方法，最后介绍了蒙特卡罗方法的 起源- - m e t r o p o l i s 模拟方法。因为随机数发生器是蒙特卡罗和拟蒙特卡罗方 法的核心之一，所以第2 章介绍了伪随机数和拟随机数发生器。此章还介绍 了如何生成服从一定概率分布函数的随机数。 第3 章提出了b 样条光滑拒绝抽样方法。第2 章介绍的标准拒绝抽样方法 其实跟特征函数的蒙特卡罗积分有密切的关系。而由于特征函数的不连续 性，蒙特卡罗积分应有的误差精度就达不到，拒绝抽样的效果也就受到影 响。我们用b 样条磨光技术在不改变积分值的前提下磨光特征函数，用可微 的权重函数代替特征函数，提高了采样的质量。将b 样条光滑拒绝方法用于 重要抽样估计数值例子显示拟蒙特卡罗积分的精度重新达到7 o ( n - 1 ) 的 阶，而对于蒙特卡罗积分，采用b 样条光滑重要抽样，其精度也比标准积分 的精度o ( 一1 好。由此间接证明了b 样条光滑拒绝抽样比标准拒绝抽样有 效得多。 第4 章是关于蒙特卡罗积分韵，得出了甩于多重积分的精细对偶变数 蒙特卡罗( f i n ea n t l t h e t i cv a r i a b l e sm o n t ec a r l o ，简称f a m c ) 方法的误差估 计式。对维数是8 的= 阶导数连续的函数来说，f a m c 方法理论误差的阶 是0 ( 一( + ；) ) 。我们同时也讨论了对偶变数蒙特卡罗积分( a n t i t h e t i cv a r i - a b l e m o n t e c a r l o ，简称a m c ) 方法。对次数不高于2 的多变量函数，a m c 方 法其误差阶是d ( 一1 ，但其系数比原始蒙特卡罗积分( m c ) 方法的误差阶 的系数小。我们用c 语言实现了并行计算程序，数值实验结果与理论结果吻 合得很好。 第5 章介绍了自适应拟蒙特卡罗全局优化( a q m c ) 方法。a q m c 算法在 不可微函数的蒙特卡罗优化算法上大大前进了一步。首先改进了局部搜 索算法，使得搜索方向，搜索半径和搜索所计算的函数值个数根据已有的 搜索结果而自适应调整。其次引进了种群演化的概念，根据演化度产生新 个体，保证搜索的全局性。因为在搜索过程中会根据搜索中间结果进行调 整，所以此方法不仅可以加快搜索速度。而且可以平衡全局和局部搜索需 求。 第6 章结合遗传程序设计方法和自适应拟蒙特卡罗优化方法用于预测问 题。在现实生活中存在许多随时间而变化的复杂系统和现象，人们通常需 要根据动态系统的观测数据建立合理的微分方程模型( 动力学模型) ，为系统 分析、设计和未来状态的预报提供依据。我们用遗传程序设计方法优化常 微分方程右端的函数，用自适应拟蒙特卡罗优化方法优化函数中的系数。 在预测杭翔市地区全社会用电量的实际应用中，所编写的程序取得了很好 的效果。 第7 章是蒙特卡罗模拟和优化的结合。首先我们介绍了光在组织中传播 的蒙特卡罗模拟的完整过程，解释了如何利用第2 章中介绍的随机数生成方 法根据实际问题产生随机数。然后用自适应拟蒙特卡罗优化方法来解决光 的传播的逆问题。在这一章我们进一步探讨了如何根据实际问题的性质来 平衡全局和局部搜索需求。 附录给出了相关程序的c 语言代码。 a b s t r a c t ih a v es t u d i e da c o u p l e o ft o p i c so nm o n t ec a r l oa n d q u a s i - m o n t ec a r l o m e t h o d s t h i sd i s s e r t a t i o nc o v e r si t sa p p l i c a t i o n si ni n t e g r a t i o n ，o p t i m i z e t i o na n ds i m u l a t i o n c h a p t e r1a n d2a r et h eb a s i ck n o w l e d g e t ou s em o n t ec a r l oa n d q u a s i m o n t ec a r l om e t h o d s c h a p t e r1p r e s e n t st h ee r r o rb o u n d so fm o n t ec a r l o a n dq u a s i m o n t ec a r l oi n t e g r a t i o nm e t h o d s b yc o m p a r i n gt h e s et w om e t h - o d s w es h o wt h ea d v a n t a g e so fq u a s i - m o n t ec a r l om e t h o d w ba l s oi n t r o - d u c et h es t a n d a r dm o a t ec a r l or a n d o ms e a r c hf o ro p t i m i z a t i o n t h el a s t b u tn o tl e a s ta p p l i c a t i o ni s m e t r o p o l i sa l g o r i t h m sw h i c hi s t h eo r i g i no f m o n t ec a r l om e t h o d b e c a u s et h er a n d o mn u m b e r sg e n e r a t o r sa r et h ek e y o fm o n t ec a r l om e t h o d sa n dq u a s i m o n t ec a r l om e t h o d s c h a p t e r2d e - s c r i b e st h ep s e u d o r a n d o mn u m b e rg e n e r a t o r sa n dq u a s i r a n d o mn u m b e r g e n e r a t o r s h o wt og e n e r a t en o n - u n i f o r mr a n d o mn u m b e rf r o mi t s d i s t r i b u t e df u n c t i o ni sa l s oi n t r o d u c e d c h a p t e r3i n t r o d u c e sb s p l i n es m o o t h e dr e j e c t i o ns a m p f i n gm e t h o d t h es t a n d a r dr e j e c t i o ns a m p l i n gm e t h o dw h i d li si n t r o d u c e di nc h a p t e r2i s c l o s e l yr e l a t e dt ot h ep r o b l e m o fq u a s i - m o n t ec a r l oi n t e g r a t i o no fc h a r a c t e r i s t i ef u n c t i o n s ，w h o s e8 c c u 如哟m a yb el o s td u et ot h ed i s c o n t i n u i t yo ft h e c h a r a c t e r i s t i cf u n c t i o n s w eu s eb - s p l i n es m o o t h i n gt e c h n i q u et os m o o t h t h ec h a r a c t e r i s t i cf u n e t i o nw i t h o u tc h a n g i n gt h ei n t e g r a lq u a n t i t ya n dg e ta d i f i e r o n t i a b l ew e i g h tf u n c t i o n t h em e t h o dc o n s i d e r a b l yi r e p r o v e st h eq u a l - i t yo fs a m p l i n gp o i n t s w ja p p l yt h eb - s p l i n es m o o t h e dr o j e c t i o ns a m p l i n g m e t h o dt oi m p o r t a n c es a m p l i n g n u m e r i c a le x p e r i m e n t ss h o wt h a tt h e e r r o r s i z eo ( n 一1 ) i sr e g a i n e db yu s i n gt h eb s p l i n es m o o t h e dr e j e c t i o nm e t h o d f o rq u a s i m o n t ec a r l oe s t i m a t e t h ee r r o rb o u n do fm o n t ec a r l om e t h o d u s i n gb - s p l i n es m o o t h e di m p o r t a n c es a m p l i n g i sa l s ob e t t e rt h a nt h a to ft h e 曼垒垒! ! ! 堕 s t a n d a r dm o n t ec a r l om e t h o d s ot h e b s p l i n es m o o t h e dr q i e c t i o ns a m p l i n g m e t h o di s i n d i r e c t l yp r o v e dt ob es u p e r i o rt ot h es t a n d a r dr e j e c t i o ns g l r l p l i n gm e t h o d 。 c h a p t e r4i sa b o u tt h em o n t ec a r l oi n t e g r a t i o r l w eg e tat h e o r e t i c a le r r o ro ft h ef i n ea n t i t h e t i ev a r i a b t e sm o n t ec a r l o ( f a m c lm e t h o df o r m 山t i d i m e n s i o n a li n t e g r a t i o n t h ee r r o rs i z eo ff a m cj s o ( n 一 + ；) 1 f e r f u n e t i o n sh a v i n gs e c o n dc o n t i n u o u sd e r i v a t i v e 。w h e r esi st h ed i m e n s i o no f t h e i n t e g r a n d w ea l s o 西v et h et h e o r e t i c a le r r o rr e s u l to f a n t i t h e t i cv a r i a b l e m o n t ec a r l o ( a m c lm e t h o df o rm u l t i v a r i a b l ef u n c t i o n sw h o s ed e g r e ei sn o m o r et h a nt w o ，t h ec o n s t a n tb e f o r eo f n 一1i sl e s st h a nt h a to ft h em c m e t h o d 。w er e a l i z et h ep a r a l l e la l g o r i t h mi ncf o rt h ef & m ca n da m c m e t h o d s t h er e s u l 招o ft h en u i n o a i c a le x p e r i m e n t sc o i n c i d ew i t ht h et h e o - r e t i e a lr e s u l t sv a r yw e l l 。 c h a p t e r5i n t r o d u c e sa d a p t i v em o n t e c a r l om e t h o d ( a q m c ) f o r g l o b a l o p t i m i z a t i o n a q m ca l g o r i t h mp r o g r e s s e si nt h en o n d i f f e r e n t i a b l eo p t i - m i z a t i o n 。f i r s t 、w ed e v e l o pt h el o c a ls e a r c hs u c ht h a tt h es e a r d ld i r e c t i o n ， s e a r c hr a d i u sa n da n m b mo fs e a r c hp o i n t sa r ea d j u s t e da c c o r d i n gt o 糠e p r e v i o u ss e a r c hr e s u l t s e c o n d 啪i n t r o d u c et h ei d e a so fp o p u l a t i o na n d g e n e r a t en e w i n d i v i d u a l sa c c o r d i n gt 0p o p u l a t i o ne v o l u t i o nd e g r e e b e c a u s e t h es e a r c hp r o c e d u r ew i t lb ea d j u s t e da c o r d i n g 蜘t h ep r e v i o u sr e s u l t ，t h e m e t h o dn o to n l ys p e e d su pt h er a n d o ms e a r c hb u ta l s ob a l a n c t h eg l o b a l a n dl o c a ld e m a n d s ( a d a p t i v ee q u a l l z a t i o n l c h a p t e r 6c o m b i n e st h eg e n e t i cp r o g r a m m i n g w i t ha q m c o p t i m i z a t i o n m e t h o dt os o l v e 娃l ep r e d i c t i o np r o b l e m s t h e r e a t em a n yc o m p l e xs y s t e m s i nr e a ll i f e 。i no r d e rt oa n a l y z e ，d e s i g na n dp r e d i c tt h es y s t e m ，w eo f t e n w a n tt om o d e lt h ed y n a m i cs y s t e m so fo r d i n a r yd i f f e r e n t i 8 1e q u a t i o n sa c - c o r d i n g t ot h eo b s e r 唰id a t a w eu s eg e n e t i cp r o g r a m m i n gt oo p t i m i z et h e t h er i g h th a n df u n c t i o n so ft h eo r d i n a r yd i f i e r e n t i 虹e q u a t i o n s a d a p t i v e q u a s i m o n t ec 8 r i oo p t i m i z a t i o nm e t h o d sa r eu s e dt oo p t i m i z et h ec o e f f i 。 c i e n t so ft h ef u n c t i o n s t h ep r o g r a mf o rt h ep r e d i c t i o no fe l e c t r i c a lp o w e r c o n s u m p t i o no fh a n g z h o uc i t ys h o w s t h a tt h eh y b r i dm e t h o di sp o w e r f u l - i nc h a p t e r7 ，w ec o m b i n et h em o n t ec a r l os i m u l a t i o na n do p t i m i z a - t i o n w ef i r s ti n t r o d u e et h em o n t ec a r l os i m u l a t i o no fl i g h tt r a n s p o r ti n t i s s u e ，e x p l a i nh o w t og e n e r a t et h er a n d o mn u m b e r a c c o r d i n g _ ot h ep r a c - t i c a lp r o b l e m su s i n gt h et r a n s f o r mm e t h o di n t r o d u c e di nc h a p t e r2 t h e n u s et h ea q m cm e t h o dt os o l v et h ei n v e r s ep r o b l e mo fl i g h tt r a n s p o r t ，w e f u r t h e rd i s c u s sh o w 协b a l a n c et h eg o b a la n dl o e a ls e a r c hd e m a n d s i np r a c - t i c a lp r o b l e m s t h ecc o d e so fs o m ep r o g r a m sa r eg i y e ni nt h ea p p e n d i x 序言 簸特卡罗方法是研究多体藏杂系统以及不确定健过程的强有力的工 其。0 ，d 6 h g a r r s 辆f 。s u l l i v a n 1 l 在杂志c o m p u t i n g i ns c i e n c e & e n g i n e e r - i n g 上发衷文章谭述t - 十世纪对科学和工程计算的发展和爨践最具影响力 的十太算法【2 j ，漾特卡罗方法榜上有名。 裁特卡罗方法的前身m e t r o p o l i s 算法是出j y o nn e u m m m ，s u l a m 和n m e t r o p o l i s 在第= 次谶界大战末为了研究餮交物质中子的扩散丽提出懿。 “蒙特卡罗”是在1 9 4 7 年出n m e t r o p o l i s 命名，并于1 9 4 9 年在n m e t r o p o - i i s 和s u h m 合作的一篇文章3 1 磁式使用。这些科学家提出蒙特卡罗方法的 初衷楚用予攻竞专门豹秘建鼗壤模掇润蘧。起糖入稍黯蘩特专爹方法豹兴 趣不大，但痿来这一方法缀快就坡应用于各个令人惊讶的方向，包括函数 值极小化，计算几何，组含计数等。时至今日，在深刻的理论支持下，瓢 物理模掇到诗算复杂链麴蒸萋，与m e t r o p o l i s 援关魏主题基经涯簸了识箕科 学的整个领域。而随着拟髓帆序列的出现，蒙特卡罗方法也已经发展到拟 浆特卡罗方法。 我们经常要面临一些有非常高维数的向题，或有很多分支路径( p a t h ) 而每一个分支路径都按一定静概率发生翡过程。因舞蒙祷专罗方法楚慰鞫 题进特采撵，所以它不是严格的、糖确的解法。但是我们能用相对于问题 的维数而富相翔小的样本得到近似糟确解。事实上，虽然荣特卡罗和撅蒙 特卡罗方法胃辘毙较耗鞋，毽窀稻憨解决凑维阏露瞧一韬实可行豹蠢法。 幸运的是，照麓现代计算机技术的飞速发展，我们可以计算越米越艇杂的 问题。所以蒙特卡罗和拟蒙特卡罗算法在新的2 l 毽鳃会有更光鞠的前景。 签予照，我选撵了蒙转专罗帮撅蘩酶每罗方法朗羞于主题揍为碜 宠方向。 序言 如果您有任何意见和建议，请联系我。 露桎媛 2 0 0 3 年4 月于杭州 致谢 酋先感谢我的导师王兴华教授，是饿指引我进入蒙特专罗和拟蒙特卡罗研 究镁域，并藏在我的研究除羧给了我缀多帮韵。 嗣辩我要戆谗下述裁骞豹计冀数学教瓣塑熬成费。扶鼗稍组定期举行躲避 论班上我学到了很多知识。 王袋华，李 孛，篷金生，箨圭明，程囊受，韩丹夫，吴疾标，羚隽籀，黄缆选，豫明飞， 叶嫩德，朱建新，杨士俊，粱克维，李大嘲，粱仙红，壬何宇，獬湘江，盒中秋，崔峰和 谳聪聪 另外还有根多人给了我许多帮助，尤其是那些让我有机会将所学的内容 应用于实践的人。金剑歌套绍我参加了关干预测髓编程工作，应嗣遗 传程序设计与自适应拟蒙特卡罗结合的方浃取褥了稚好的效莱。a n d e r s k a r l s s o n 教授与我讨论了光在组织中传播的逆问题，使我得以将自己的优化 方法用子瓣捷这一实薅润瑟。 最后，我袋感谢我的家人和朋友。他们付出的关心和盘持让我实现了求学 豹梦想，簸终褥予成为群磷芏终者。我豹丈失舞汪警鞲士一巍鼗精我数奋 学习，在我准备第一篇论文时我们谯他的实验室一越工作到深夜，那是一 段令太难簿戆鞋毙。在这麓嚣学粳论文豹准餐中 彀也绘了缀多帮助，还 教我画了其中的一魑图。 第一章蒙特卡罗和拟蒙特卡罗方法的基 本知识 蒙特卡罗和拟蒙特卡罗方法在数值方法中应用得很广泛，人们研究这 些方法也已经很多年了。其中最常用的应用是蒙特卡罗积分，蒙特卡罗和 拟蒙特卡罗方法同样广泛地应用于优化和模拟。因为蒙特卡罗方法简单直 接，所以很容易应用。这些方法同样很实用，因为方法的精度仅依赖于问 题复杂性最粗糙的刻划。这章主要介绍蒙特卡罗和拟蒙特卡罗方法的的一 些基本知识。 1 1 蒙特卡罗积分 关于蒙特卡罗积分，c a f l i s c h 4 1 在1 9 9 8 年给出了一个综述。勒贝格可积 函数，( x ) 的积分可以表示成函数，在随机点上函数值的平均值或期望值。考 虑在8 维单位超立方体( u n i tc u b e ) l = 【0 ，1 1 4 上的积分，刚有 j = e 【，( x ) 1 _ z ，( x ) d x = ， ( 1 1 ) 这里x 是在单位超立方体上的均匀分布的随机数向量。 原始蒙特卡罗估计( c r u d em o n t ec a r t o ( m c ) e s t i m 哟r ) 是基于对积分的 概率解释。若有一服从均匀分布的序列靠，则有经验近似估计式： m = 专薹玳t ， 根据强大数定律( s t r o n gl a w o fl a r g en u m b e r s ) j 5 ， 也就是： 溉m 。i ( 1 2 ) 此估计以概率1 收敛 ( 1 3 ) 2 蒙特卡罗和揪雾l 特卡罗方法的基本知识 另外，此近似式是无偏的，这就意昧着对任意，肘的平均值就等予，。落 就是： e l m j = j 这里是在选择的点& 上取平均德 一般来说，藩特卡罗积分误差定义成： ( 1 4 ) ：i 一材( 1 5 ) 这样偏蓑( b i a s ) 是e e j ，均方撤误差( r 打醛0 定义痿： e 陋纠2 ( 1 砩 1 1 1 误羞估计 中心极限定理( c l t ) 【5 】描述了蒙特卡罗积分的误差的大小和统计性质。 定理1 1 。1 对任意太的， a n l 2 pl 。7 这里y 是标准正态分布琏瓤变重f n ( o ，1 ) ，常数扩= 萨f 嗣是，的根方差( s q u a r e r o o t 醪t h ev a r i a n c e ) ，也姥是税： 蚶吲z 。地) 叫2 d x ) 疆8 ) 更精确的表迭式是 脚婚 o l ：= p 卜r o b ( a 广 盼v 即b ) 帆( 1 _ 9 ) 这就是说裴特卡罗积分的误差照以0 ( _ 1 2 ) 为阶，系数是被积函 数，的根方蓑。另外，误差的统计分布大致是一个正态分布变量。与通常酶 数值分耩楚果相比较，这是一个统计结果。它并不提供一个绝对的误菱上 界，而认为误差是概率意义上的估计。另一方面，此估计是一个等式，所 以又是严格的。 1 2 拟蒙特卡罗积分 3 1 1 2 蒙特卡罗方法与网格方法的比较 很多人第一次接触蒙特卡罗方法时会惊讶于此方法的可行性。为什么 随机排列会比网格( g r i d b a s e dm e t h o d ) 更好? 可以从几个方面来回答这个问 题。首先，可以比较蒙特卡罗和网格方法如辛普森方法( s i m p s o n sr u l e ) 的 收敛阶。对于8 维问题的阶网格求积方法的收敛阶是o f - k s 1 。而蒙特卡 罗积分的收敛阶是d ( _ 1 2 ) ，与维数没有关系。所以在高维积分时，当维 数s 满足k s 1 2 时，蒙特卡罗积分方法比网格积分方法精度高。 然而，对于周期解析函数( a n a l y t i cf u n c t i o n ) ，k 的值是无穷的，也就是 说上述简单解释不成立。对蒙特卡罗方法的实用性的更强有力的解释是实 际应用中很难对高维问题进行网格化。对8 问题最简单的超立方网格至少 需要2 s 个点。对8 = 2 0 这样在实际问题中并不很大的维数。就要有1 百万多 个点。另外，要对高维问题进行网格细化就更不现实了，因为每一次细分 网格都要增加2 5 倍的点。相比于网格方法在高维问题中碰到的这些困难， 蒙特卡罗方法的精度几乎与维数无关，而且点的增加都会提高误差精度。 在实际应用时，取个点，精度o ( n - 1 仰并不一定能达到。但经验显示， 对于中等维数( 比如8 = 2 0 ) 的中等复杂的问题，只要中等规模的值，误差 阶o ( n _ 1 2 卜般都能达到。 1 2 拟蒙特卡罗积分 回想蒙特卡罗积分，用个随机点，其积分误差阶的平均数量级 是o ( n 一1 2 ) 。显然，肯定存在这样的个点，使得误差的绝对值不太于平 均值。如果我们能构造这样的点集，就可以对原有的方法进行较大的改 进。拟蒙特卡罗积分方法就是为了发展数值积分而提出的。它致力于构 造其积分误差比平均误差显著要好的那种点集。其积分形式与蒙特卡罗 积分一致，只不过所用的随机数不一样。比如，对于单位超立方体积分区 域，。= 【0 ，1 1 s ，我们定义拟蒙特卡罗估计如下： ， 小x ) d 州击三瓶) ( 1 1 o ) 此式看起来象蒙特卡罗估计，但这里用的是确定性( d e t e r m i 【1 i 8 t i c ) 的点x 1 ， ，) 哪i s 。我们可以明智地选择这些点使得式( 1 1 0 ) 中的误差得以变 小。 4 蒙特卡罗和拟蒙特卡罗方法的蕊本知识 因此，拟蒙特卡罗积分方法的基本思想是用精选的确定性点来取代 嫠特卡罗双分中的随机数点。选择点的标准取决于数值汹题本身。剥于 数值积分而言，选择的的标准穰容易找到，而越由诧；l 进了均匀分布窿 歹l j ( u n i f o r m l yd i s t r i b u t e ds e q u e n c e ) 和偏差( d i s c r e p a n c y ) 的概念。 1 2 1 偏差 绱茇罨殴霉残是对均匀努毒豹偏离豹量铑囊量。 记p 是包含点2 l ，x 1 8 的点集。对任意瓣于五豹子集b ，我们定 义 a ( 县；p ) = x 学( x a ， ( 1 。1 1 ) ；l 这里x 聍是b 的特征函数( c h a r a c t e r i s t i cf u n c t i o n ) 。这样a ；p ) 便是计算满 足酝b ( 1s 锋兰) 的那些熹豹个数的函敬。撼募器是滏于j 5 的勒哭掇测 度子集的非空集合，则点集p 的偏差表示成； 珊= 潞l 掣屯跳 1 1 2 ) 这里沁是8 维羲哭揍测瘦，d n 拶；玛总商osd n ( b ；p ) l e 作集会8 作适 当的专门化，我们有三种最重器的偏差概念。我们记p = 【o 1 ) 。 定义1 , 2 1 ( 曼偏差) 点集p 的燕偏差( s t a r d i s c r e p a n c y ) 可以定义为d ；：，( p ) 一 域红l ，) 嗡) = 玩( 了4 ；p ) ，这里+ 是p 的共有形式n 釜l 【0 ，钺) 的所有 子区阏的集合。 定义1 2 。2 ( 极证偏差) 点褰p 的摄也倍差愀m 8d i s c 御a n c y ) 警找定义 为d ( p ) = d n ( x 1 ，x ) = d ( 了；p ) ，这里了是产的具有彤式n 釜1 陬，毽) 的所宥子区问的集合。 定义1 2 3 ( 筹方性偏藏) 点粜p 的等方性偏差似o t r o p i cd i s c r e p a n c y ) 可以 定叉两括尹) = 西( x l ，x ) = 确( c ；玟，这里c 射4 妁所密插手渠的 集奋。 偏差的往质在【6 1 中有阐述。下蟊的误差 古诗磐辑氇来自这本书。 1 2 拟蒙特卡罗积势5 1 2 2 k o k s m a - h l a w k a 不等式 我们讨论关予式( 1 1 0 ) 所计算的叛蒙特卡罗积努的误差估计的重要结 果。我们从一维情形开始，其中一个经典的结果就是下面的k o k s m a 不等 式。 定理1 2 。1 如慕热数，在撼两【o ，l 】生毒毒赛变分如随睡。嘲v ( d ，粼赌任 意十点1 ，【0 ，1 】，有 ni 1 三n 蚤苁一五施 双力璐犯，i ( 1 1 3 ) 定理1 , 2 2 如暴蕊教，在聪问如，1 l 土是连续的，则对任：卷今点z 一，z 【0 ，1 】，有 i 斋苫n m 沪j ( 1 删删 0 ，即状态扎是 比状态m 能量高的状态，此时系统以概率加p 。转到状态n 。对于正则系 统( c a n o n i c a le n s e m b l e ) ，i p n v t 恒定的系统 堕：e x p ( 一卢6 。) p m 其中口是玻尔兹曼常数。所以产生一个随机数，如果 e x p ( 一础峙。) ， 则系统转到状态n ，否则系统维持在原来的状态。综上所述，系统从状 态m 转到n 的概率是：r a i n ( 1 ，e x p ( 一卢6 k 。) ) 。系统从初始状态一直重复下面 两步： 步骤1 ：试走步，即计算下一步的能量： 步骤2 ：计算走到下一步的概率m i l l ( 1 ，e x p ( 一卢6 k 。) ) 并决定系统下一 步的真正状态。 这是一个马尔科夫( m a r k o vc h a m ) 过程，最终系统会达到稳定状态。 其实m e t r o p o l i s 蒙特卡罗模拟的关键是计算系统转移到下一状态的概 率。我们将在第7 章介绍光的传播的蒙特卡罗模拟方法，同样强调了如何计 算各事件的概率。 1 5 阅读资料 有很多书本介绍蒙特卡罗和拟鼗特卡罗方法。其中 1 3 1 4 1 1 1 5 覆盖积 分、优化、模拟三个方面主题，另外【l6 】注重于随机数发生器，还有【1 _ 7 】则 讨论蒙特卡罗方法在生物领域的应用。 第二章随机数发生器 随机采样是豢特卡罗方法的核心。蒙特卡罗方法酌成功当然取决于随 机摸型的构造，但缀犬理发上也取决卡摸型计冀中照枫数熟性震。这一孝 着重干随机数发生器。 般来说随机数势残均匀醚机数( u n i f o r mr a n d o mn u m b e r s ) 鄹非均匀随 机数( n o n u n i f o r mr a n d o mn u m b e r s ) 。均匀随机数的目标分布( t a r g e td i s t r i b u t i o n ) 楚j 上的均匀分布( u n i f o r md i s t r i b u t i o n ) u 。产生菲溺匈醚梳数逶常蹙先 生成均匀随机数，然后再转换成服从舄标分在f u 的随机数。我们先介 绍伪随机数( p s e u d o r a n d o mn u m b e r s ) 和拟随机数( c l u a s i r a n d o mn u m b e r s ) 发生 器，然嚣奔绍转换方法。 2 1 伪隧视数 在使用蒙特卡罗方法的早期，入们就酷经认识到在现实生活中不可能 骞襄惩的4 夔瓤”数。于是人们求助于能遇过一定冀法及几个参数恩计算 机产生的伪随机数( p m , t ) 。 一令缴典豹且现在扔然广泛使用的产生均匀镄隧枫数豹方法是线性 同余方法( 1 i n e a rc o n g r u e n t i a lm e t h o d ) ，此肖法澄初l i l e h m e r 1 8 提出。线 性同余方法有三个参羧，其中解是值穰穴的征整簸，窃燕满楚is 吐 p ，我 f 拜f 递归公式 撬：龟毽一l 国盘豫一2 0 辔嘞旗呻o b ，妒l 2 t ) 这里。表示二进制按位异或( e x c l u 8 i v e o r ) 。对于 k ，对等的递归公式怒 m 产2 c l 帆102 2 c 2 挑一2o 02 p 邻m l - p m i 呻 ( 2 5 ) 举个铜子，取褥单多项式 对应静谦爨公式楚 茹4 + $ + 1 l l j , i 。8 m l 一3o1 6 i r k 一 o 豫一 1 2 随机数发生器 如果我们取m 的初值为m 1 = 1 ，2 = 1 ，3 = 3 和m 4 _ - 1 3 ，则确 m 5 = 8 0 1 6 0 1 = 0 1 0 呻2o1 0 0 0 0 20 0 0 0 0 12 = 1 1 0 0 1 2 = 2 5 则可以得多j s o b o l 序列的第i 个数 以= b l v l o b 啦o b o 这里6 3 b 2 b l 是珀q 二进制表示形式。 a n t o n o v 和s “e e v 【2 4 】提出了另一种产生s o b 0 1 序列的简单方法。此方法 用下面的递归公式 筑= g i r l o 船t 1 2 0 船地o ( 2 6 ) 其中9 3 9 2 9 1 是以i 为变量的格i t 马( g r a y c o d e ) 的二进制表示。格雷码g ( i ) 是 以非负整数沩变量的函数，使得g ( t ) 和g 0 + 1 ) 的二进制表示只有一个数 位不同，即在g ( ) 的最右一个零位上不同。a n t o n b o v 和s a l e e v 使用的二进 制表示的格雷码是利用下式产生的。 9 3 9 2 9 l = 6 3 6 2 6 1o 6 4 6 3 6 2 ( 这是最常用的格雷码，能产生一系列函数值0 ，i ，3 ，2 ，6 ，7 ，5 ，4 ，) 于是产 生s o b o l 序列的( 2 6 ) 式可以用下列式子代替。 z t = z i 一1 0 r 这里柏自选择是使k 是 一1 的二进制表示的最右零位。 b r a t l e y 和f b 】【在文1 2 5 l 中讨论了起始m l ，m 2 ，点的取值( 此章例子中所 取的简单多项式对应的m 的起始值满足那些标准) 。 图2 1 是一个两维s o b o l ，序列，图2 2 是s o b o l 序列与伪随机数序列的 比较。函数r 伽2 是书【2 6 l 中一个伪随机数发生器。书【2 6 l 中的s o b o l 停列 的c 语言代码经过修改，放在附录中。 2 2 拟随机数 1 3 兰一 。4 4 絮酬气嚣气 阳r f m l 1 i 强2 1 ： 掰维s o b o l ，醚祝数，裁产生的相继的点分奄在先翦产生的点的空隙 中。 d 2口 0v _ 喇。州蛐c 霜嚼坤“h 1 1 0 1 i h 。 蕊 睽 岛 a j o 跏b 时_ 懈h 1 图2 ，2 ：拟随机数序列比伪随机数序列有更好的均匀性。 1 4 隧枧数发生援 2 2 3 n i d e r r e i t e r 的，ms ) 网格和( t ，s ) 序列 这一节我们介绍( t ，m ，$ ) 网格和( t ，s ) 序列。这种序列是非常规贝唾分布的 序到。 为了定义此序列。我们记维数s l ，整数b 2 ，并称p 的子区问e e 兰【啦6 哦，( m + 1 ) 6 一孟) i = l 为以6 为蒸的基本f h - j ( e l e m e n t a r yi n t e r v a l ) ，其中对任意l i s 有如，盔 z 且满足呔0 ，0 蔓如 t ，满足舳佗( k + 1 ) b m 的点x n 组 成戆点拳矮是鞋焉基蝣抟，m ，8 ) 建格，酗爵垮列秘，x l ，是$ ) 摩列。 2 3从均匀随机数到毒# 均匀随机数的转换方法 从均匀随机数通过转抉能产生服从一定概率分布函数的非均匀分布随 帆数。j 鞋：蘩奔缓滤方法是骛逶戆，隳为可以焉寒产擞绝大多数分毒的夔枫 数。 2 3 1 桊辍分韶函数( c d f ) 求邀方法 如果茹是有连续的概率累积分布函数( e d f ) r 的随机标量，则由下式产 生越箍撬爱重有u ( o ，1 ) 分程a = 逸( 嚣) = 艄蠡 ( 2 f 7 ) 这一常谈霹越让我们在均鼙努毒涟祝变量c ，和睫羲燮藿髫之间建立一个篱萃 关系，那就是： g 一巧1 ( 娃 2 3 从均匀随机数到非均匀随机数的转换方法1 5 我们称这一直接了当的转换方式为c d f 求逆技术。因为分布函数求逆容易 计算，所以此方法是很好的产生随机变量的方法。然而，因为计算某些分 布函数的逆函数并不容易，c d f 求逆方法不象预想得那样使用广泛。 2 3 2 拒绝方法( r e j e c t i o nm e t h o d ) 令p ( x ) 是定义在j 。上的概率密度函数。标准的拒绝采样算法如下描述 1 选择7 使得7 s u p x j 。p ( x ) 2 重复下列步骤直至个点被接受： ( a ) 从u ( i o ，1 1 叶1 ) 采样随机向量( x t ，轨) ， ( b ) 如果有y t 7 - 1 p ( ) 则接受点x ： 否则，拒绝这些点。 跟其它产生非均匀随机鼓方法类似，拒绝采样方法也依赖于所使用的均匀 随机数的均匀性。 第三章b 样条光滑拒绝采样方法及其在 重要抽样中的应用 拒绝采样方法是蒙特卡罗方法中最流行的采样算法之一。拒绝采样方 法其实跟特征函数( c h a r a c t e r i s t i cf u n c t i o n s ) 的蒙特卡罗积分有密切的关系。 而由于特征函数的不连续性，蒙特卡罗积分应有的误差精度会达不到，拒 绝采样的效果也就受到影响。我们在文 2 7 1 9 提出了b 样条光滑拒绝采样方 法，通过用b 样条磨光技术在不改变积分值的前提下磨光特征函数。我们 将b 样条光滑拒绝采样方法用于重要性抽样方法中，数值例子显示拟蒙特卡 罗积分的精度重新达到了o ( n - i ) 的阶，而对于蒙特卡罗积分，采用b 样条 光滑重要抽样，其精度也比标准积分的精度0 ( 一 ) 好。由此证明，b 样条 光滑拒绝抽样比标准拒绝抽样有效得多。 3 1引言 标准拒绝采样算法在诸如决策过程，服从密度分布函数的采样等拟蒙 特卡罗实际应用中有非常重要的地位。但是它并不如理论结果所期望的那 么有效。 拒绝采样算法可以解释成对特征函数的积分。回想第1 章所介绍的拟蒙 特卡罗积分，其积分误差j 瞩k o k s m a - h l a w k a 不等式描述，误差精度的阶 - - e 为o ( n 一1 ) 。但是除非函数定义域是边平行于坐标轴的长方形，否则特 征函数的变分是无穷的。这样积分误差的上界就不能用k o k s m a - h l a w k a 不 等式来估计，o ( n “】的理论误差也就无法达到。 我们用b 样条