（概率论与数理统计专业论文）一类复杂密度函数抽样的研究.pdf

一类复杂密度函数抽样的研究 中文摘要 一直以来，抽样在统计学中发挥着重要的作用目前，抽样已被广泛地运用到科 学计算中，尤其是统计计算领域从简单的抽样法，到比较复杂的m o n t ec a r l o 方法， 再到基于m o n t ec a r l o 的各种方法，抽样不断更新发展 w r g i l k s 和p w i i d 在1 9 9 2 年提出了适应拒绝抽样来处理任意单变量对数凸概 率密度函数，以函数在某些点的切线构成一个分段包络函数进行抽样本文也是运用 包络的思想，讨论了一类复杂的密度函数，主要是指多峰形式的密度函数，并将重 尾分布的思想溶入到获取函数包络的过程中所得的包络函数形式比较简单，并易于 得到样本，提高了抽样效率 关键词：抽样，复杂密度函数j 包络 作 者：龚焱华 指导老师：汪四水( 副教授) s a m p l i n gf r o mak i n do fc o m p l e xd e n s i t yf u n c t i o n s a b s t r a c t s a m p l i n gh a sb e e np l a y i n ga ni m p o r t a n tr o l ei n s t a t i s t i c s a n dn o wi th a sb e e n w i d e l ya p p l i e d t os c i e n t i f i cc o m p u t i n g ，e s p e c i a l l yt os t a t i s t i c a lc o m p u t i n g f r o m s i m p l es a m p l i n gm e t h o d st oc o m p l e xm o n t ec a r l ow a y s ，a n de v e nt om e t h o d sb a s e d o nm o n t ec a r l o s a m p l i n gi sg o i n go nd e v e l o p i n ga n du p d a t i n g ， w r g i l k sa n de w i l dp r o p o s e dam e t h o dc a l l e da d a p t i v er e j e c t i o ns a m p l i n gt od e a l w i t ha n yu n i v a r i a t el o g c o n c a v ep r o b a b i l i t yd e n s i t yf u n c t i o ni n19 9 2 t og e ts a m p l e s ， t h e yu s et a n g e n t sa ts o m ep o i n t sa se n v e l o p ef u n c t i o nw h i c hi sp i e c e w i s e t h i s a r t i c l e w ew i l ld i s c u s sak i n do fc o m p l e xd e n s i t yf u n c t i o n s ，w h i c ha r eu s u a l l ym u l t i m o d a l a n dw ea l s ou s et h ei d e ao fe n v e l o p ef u n c t i o n ，a n dt a k ea d v a n t a g e so ft h e h e a v y - t a i l e df u n c t i o n st oo b t a i ne n v e l o p e t h ee n v e l o p et h a tw eg e th a sar e l a t i v e l y s i m p l ef o r m ，h e l p i n g u sg e ts a m p l e sm o r ee a s i l y , a n dt h e ni m p r o v i n gt h ee f f i c i e n c y k e y w o r d s ：s a m p l i n g ，c o m p l e xd e n s i t yf u n c t i o n ，e n v e l o p e w r i t t e nb yg o n gy a n h u a s u p e r v i s e db ya s s o c i a t ep r o f w a n gs i s h u i 苏州大学学位论文独创性声明及使用授权声明 学位论文独创性声明 本人郑重声明：所提交的学位论文是本人在导师的指导下，独立 进行研究工作所取得的成果除文中已经注明引用的内容外，本论文 不含其他个人或集体已经发表或撰写过的研究成果，也不含为获得苏 州大学或其它教育机构的学位证书而使用过的材料。对本文的研究作 出重要贡献的个人和集体，均已在文中以明确方式标明。本人承担本 声明的法律责任。 研究生签名：垒美聋日期：兰2 三2 ：竺呈 学位论文使用授权声明 苏州大学、中国科学技术信息研究所、国家图书馆、清华大学论 文合作部、中国社科院文献信息情报中心有权保留本人所送交学位论 文的复印件和电子文档，可以采用影印、缩印或其他复制手段保存论 文。本人电子文档的内容和纸质论文的内容相一致。除在保密期内的 保密论文外，允许论文被查阅和借阅，可以公布( 包括刊登) 论文的 全部或部分内容。论文的公布( 包括刊登) 授权苏州大学学位办办理。 研究生签名：丝基登 日 期：兰1 2 ：竺呈 第一章随机变量的常用抽样方法 1 1 引言 随机模拟方法，又称蒙特卡罗( m o n t ec a r l o ) 方法，是利用计算机进行数值计 算的一类特殊方法它适用的范围非常广泛，既能求解确定性的问题，也能求解随机 性的问题以及科学研究中的理论性问题 在求解数学、物理、工程技术、随机服务系统及备受关注的经济、管理等问题 方面，随机模拟方法通过构造一个概率模型，使所求问题转化为该模型的参数、特 征量或与之有关的量，然后由模拟，得出模型参数或特征量的估计值，从而得出所 求问题的近似解 用随机模拟方法求解实际问题时，首先要解决的是随机数产生的方法对于所涉 及随机现缘形式各异的分布规律，这就要求产生对应于该分布规律的随机数只有得 到了相应分布的随机数后，才能进行模拟计算，称产生某个随机变量的随机数这一 过程为对随机变量进行抽样，并且称产生这个随机变量的随机数的方法为抽样法 对于随机数产生方法的研究已有较长的历史，最早方法称为手工方法，即采用 抽签、掷骰子、摇号等方法，很多电视节目中的抽奖至今仍采用这种方法由于对随 机数的大量需求，手工法已不能满足需要，各类随机数表便应运而生在很长一段时 间里，人们就是利用这些随机数表进行统计模拟计算的随着计算机和模拟方法的广 泛应用，用计算机产生随机数成为新的课题人们以概率统计理论为主要基础，得 到了很多简单易行的抽样方法，并在这些抽样法的基础上，不断改进、更新、推广， 促使了抽样方法的不断发展，在提高抽样效率和确保准确性的前提下，使其使用范 围不断延伸，用于处理更为复杂的的问题，更好地服务于科学发展的各个领域 本文着眼于解决一类比较复杂的分布的抽样问题，这种“复杂性 主要表现为 其密度函数通常有多个众数( 或称为多峰的) 对于这类多峰问题，已有相当一部分 的文献给出了很好的解决，本文主要是在拒绝抽样法的基础上，寻找一个简单易行 的包络函数，使得抽样过程更为简便 1 2 随机变量的常用抽样方法 1 直接抽样法( 也称反函数法或逆变换法) 直接抽样法是依据如下引理得到的 引理1 1 设f ( x ) 是连续且严格单调上升的分布函数，它的反函数存在，记为 f - 1 ( 石) ，即f f _ 1 ( x ) 】= x ，则有如下两个结论： 若随机变量x 的分布函数为f ( 砖，则f ( x ) u ( o ，1 ) ； 若随机变量r u ( o ，1 ) ；则f 1 ( 尺) 的分布函数为f ( x ) 由引理1 1 知，如果随机变量r u ( o ，1 ) ，f - 1 ( ) 是分布函数f ( x ) 的反函数，则 x = f 。( r ) ，( 功 ( 1 ：1 ) 进一步的，若已知随机变量x 的密度函数f ( x ) ，由引理1 1 知 r ：f zf ( t ) d t ( 1 2 ) = i ( 1 ) 利用公式( 1 1 ) 或( 1 2 ) ，进行直接抽样法： 产生r u ( o ，1 ) ； 代入f - 1 ( ) ，即有彳= f - 1 俾) f ( x ) ( 或代入( 1 2 ) 式，解关于x 的方程， 便有x ( 力) 2 变换抽样法 变换抽样法在概率论中有如下的理论依据： 引理1 2 设随机变量x 具有密度函数厂( z ) ，y = g ( x ) 是随机变量x 的函数，又 设工= g 1 ( j ，) 会j j l ( j ，) 存在且有一阶连续导数y = g ( x ) 的密度函数为 p ( y ) = f c h ( y ) 瞰y ) | 根据引理1 2 ，若尺u ( o ，1 ) ，随机变量x = g ( r ) 的分布函数为f ( x ) ，变换抽样 法的抽样步骤为： 产生尺u ( o ，1 ) ； 计算g ( r ) ，x = g ( r ) f ( x ) 实际上，直接抽样法是一种特殊的变换抽样法 2 表i 给出了随机变量r u ( o ，1 ) 的一些函数变换公式 表1 均匀分布随机变量尺的变换公式 变换公式x = g ( r )x 的密度函数p ( 功彳的取值范围 z = a r + b 南( 均匀分布) b ，b + a a o b + a ，b a 0 ，指数分布)( 0 ，o o ) 允 z = 够( 月一j 1 )五丽1 ( 柯西分布) ( a o ，0 0 ) 3 合成法 合成抽样的想法是：如果随机变量z 的密度函数f ( x ) 难于抽样，而z 关于】，的 条件密度函数( 工1 ) ，) 以及y 的密度函数g ( y ) 均易于抽样，则x 的随机数可由如下过 程产生： 由】，的分布g ( y ) 抽取y ； 由条件分布厂( x i y ) 抽取z 则由上述两步所得的x 服从密度f ( x ) 例1 1 设x 的密度函数厂( z ) = q z ( 工) ，其中各 o ，= l ，z ( x ) 是密 i = l i = l 度函数令= 0 ，由合成法获取x 的随机数的过程如下： 产生灭一u ( o ，1 ) ； 卜ij 确定f ，使得吩 “q ； j = oj = o 由z ( 工) 抽取x 如，( 工) = ( 1 + 2 x ) 6 ，0 x 2 ，对( 工) 进行如下分解： 厂( x ) = ；1 三+ 詈主，o j 2 ， 其中石( 工) = 五1 ，五( 工) = 主，= j i ，= 詈，结合逆变换法，具体的抽样步骤为： u ( o ，1 ) 独立地抽取，u 2 5 计算x = 【2 2 u 2 “, ：u ，i “ 11 1 3 , 3 4 拒绝抽样 ( 1 ) 抽样过程 拒绝抽样( r e j e c t i o ns a m p l i n g ) 是从某个分布函数或密度函数抽取独立样本的一 般方法设，( x ) = c x ( x ) 是可计算的，万( 力为概率分布函数或密度函数，c 是未知的 如果可以找到一个抽样分布g ( x ) 和一个“覆盖常数 肘，使得对所有的z 满足包络 性：m g c x ) ，( 工) ，则可以运用下述过程进行抽样 拒绝抽样： 产生，u ( o ，1 ) ，从g ( x ) 抽取工，与x 独立，计算比率 后= 器( 1 ) ；坛( x ) 、一 如果，后，接受x ，_ ) c 万( _ x ) ；否则，返回 引理1 3 由上述抽样过程所得样本工满足：z “万( z ) 证明：令指示函数i ： ，f l ，x g ( x ) 被接受 。一1 0 ，z g ( x ) 未被接受 则有： 川1 ) = 脚= 1 盼枇肛，器咖胁云 因此， p 卜1 ) = 器酬州：1 m 4 ( 2 ) 抽样效率 由于接受一个样本所需操作的期望次数为膨，关键就是寻找到一个较好的分布 g ( x ) 使得m 尽可能的小，保证较高的效率 5 马尔可夫链蒙特卡罗方法( m c m c ) m c m c ( m a r k o vc h a i nm o n t ec a r l o ) 方法是最近发展起来的一种简单且行之有 效的b a y e s 计算方法 m c m c 方法的基本思想是通过建立一个平稳分布为万( 曲的马尔可夫链( m a r k o v c h a i n ) 来得到万( x ) 的样本 定义1 1 随机变量序列仁价，z 【n ，z ( 舢，) 满足：在任一时刻f ( f o ) ，序列中下 一时刻f + 1 处的状态彳h 1 一p ( x l x ) ，它只依赖于时刻f 处的当前状态x ( n ，而与时 刻f 以前的状态仁( 们，j ( i ) ，x ( f - 1 ) ) 无关，则称这样的随机变量序列为马尔可夫链 定义1 2 设留m l ：。为z 上的马尔可夫链，称一步转移概率函数p ( ，) ： p ( x ，x ) p ( x x ) = p ( x = 工l x = x ) ( 离散) 或 p 一8 ) = l p o ，x ) 出( 连续) 为马尔可夫链伍( f ) l ；。的转移核，p ( 工，x 7 ) 表示由当前状态石转移到下一个状态一的概 率 类似地，有z 步转移概率函数： p ( t ；x ，x ) 垒p ( x o + 盯= ) c i x 。= x ) ， 即为当前状态x 经过t 步转移到状态x 的概率 定义1 3 如果转移核p ( ，) 与时间f 无关，则称与之相应的马尔可夫链是时间齐 次的 定义1 4 设x o 的分布为( x ) = p ( x o = z ) ，定义经过t 步后x ( 的边际分布为 ( x ) = p ( x 7 = 工) 定义1 5 如果x ( x ) 满足 ip o ，x ) x ( x ) d x = x ( x ) ，v x z ， 则称x ( x ) 为转移核p ( - ，) 的平稳分布 引理1 4 如果经过t 步后x ”的边际分布为万( x ) ，则马尔可夫链在时刻t 2 后的 任何状态的边际分布都是万o ) 至此，m c m c 方法可以概括为以下三步： 在z 上选择一个“合适”的马尔可夫链，使其转移核为p ( ，) ，这里的“合 适”主要指x ( x ) 应是其相应的平稳分布； 由z 中某一点x 出发，用中的马尔可夫链产生点序列：五，五。： 对某个m 和大的n ，任意函数f ( x ) 的期望估计如下： 岛厂5 志，萎。f ( x ( 0 m ) 以一，= = ：， 6 m e t r o p l i s - h a s t i n g s 算法( 以下简称m h 算法) m e t r o p l i s 等人在1 9 5 3 年提出了一种构造转移核的方法，h a s t i n g s 随后对之加以 推广，形成了m e t r o p l i s h a s t i n g s 方法，其思路如下 引理1 5 任意的不可约转移概率g ( ，- ) 以及一个函数口( ，) ，0 口( ，) 1 ，对任一 组合 ，x ) ( x x ) ，定义 p ( x ，工7 ) = q ( x ，x ) 口( 石，x 7 ) ，x j 7 ， 则p ( x ，工) 构成一个转移核 设平稳分布为万o ) ，时刻f 处于状态x ，即x 。= 工，m e t r o p l i s - h a s t i n g s 算法的 步骤： 产生，g ( i x ) ，这里的g ( l 工) 称为建议分布； 产生r u ( o ，1 ) ，如果r a ( x ，z ) ，则令x ( 川1 = ，否则，z ( 川) = x 在有了目( ，) 后，为了使相应的p ( x ，x ) 以万( x ) 为平稳分布，常取： 口c x ，x ，= m i n t ，弓罢主亨孑乏等 ， 相应的p ( x ，x ) 为： fg ( x ，工) ，x ( x 7 ) g ( x ，x ) x ( x ) q ( x ，x ) 烈毛，卜卜，曲等州冰如拈) 。 6 第二章等能量抽样 2 1 等能量抽样( e q u - e n e r g ys a m p le r ) 等能量抽样( 以下简称e es a m p l e r ) ，从能量函数角度出发，利用温度和能量的 对偶性，可用于处理多维或多峰分布的抽样问题 在具体说明它的抽样过程之前，先对一些概念和符号加以说明 设万( 石) 为目标分布，h ( x ) 为相应的能量函数，有石( x ) o ce x p ( 一办( 工) ) 设能量水平 序列为： 风 且 马 圾 h 川= ， 其中风_ i n f x h ( x ) ，相应的温度序列为： l = 瓦 五 o ； ( i i i ) 对于所有的f ，歹能量集概率= 乞( x d a o 则x 订是遍历的，且以乃作为它的平稳状态分布。 证明：利用向后归纳法来证明 对于f = k 。x k x 足是不可约的和非周期的，下面的证明就是基于这个事实的 现假设对于第( f + 1 ) 阶链结论成立，也就是z ( 州是遍历的且具有平稳状态分布 一- 现在要证明结论对x o 也成立根据e e 抽样的构造过程，如果第以步有砖o = 工， 则从转移核t 滋( 五) 中抽取x ( + o 的概率为1 一儿，e ! ：! l 以概率儿等于一个来自于 o l i + 的y 的概率为 只球加南幽 描灿嘲 因此，对于任意的可测集a 。有条件概率 9 尸( 砖2 彳 砖d = z ，x ( 州) = ( 1 一儿) 磁( 工，彳) + 以网1 州i ( y a ) r a i n ( 1 , 嬲， 训1 一网1 ，毛m 砸1 ，裟糍m ， 利用归纳法对z 似1 遍历性的假设，又高阶链似不受低阶链z ( 订的影响，当刀j 时，有 p ( 、x 训 0 么l 曩o = x ) = j p ( 墨3 么陋) - x ，x ( i + i ) a e ( x 似1 ) | 霹) = z ) = i p ( x y 彳f 曩o = x ，x ( m ) 护州) 一( 1 一儿) 磁( 毛一)( 2 1 ) + 死去k 啪删，署高薏播，砂 圳卜赤k ，啪m ( 1 ，鬻劣川地创 类似的，当玎一0 0 时， p ( 雹尝彳j 雹。- - x ，x ；_ i ，耐。) 一尸( 霹2 彳l 雹。= z ) 一o 现在，定义一个新的转移核s d “) ，它以概率l 一儿发挥转移核磁0 ，) 的作用， 以概率儿发挥具有建议密度q ( z ，少) = i 啬氟。( 力，( 少4 ( ，) ) 的m h 转移的作用注意 尸j “l 到( 2 1 ) 式右边对应于转移核s 0 ，) ，因此，在归纳假设下，z ( 订也渐近的等同于一 个由s f ) ( x ，) 控制的马尔可夫序列 又因为核磁( 石，) 可以连接相邻的能量集，r 建i , 3 ( q ( x ，j ，) 可以连接同一能量集上 的点，它满足c h a p m a n - k o l m o g o r o v 且0 o ) = l ， i ；if ( x ) 0 f 一( x ) = l f ( x ) ，工 0 全体这样的分布族为f 定义3 1 称分布f 属于重尾分布族r ，若对v 兄 0 ，e e 盯= f p 舡卵( x ) = 称旷= f 茁为轻尾分布族 2 对数正态 如果随机变量x 一( ，仃2 ) ，则，的分布就是对数正态分布根据引理2 2 ，可 得对数正态分布的密度函数： 一坠竺) 2 p ( x ) 2 了蒜i e2 r ，石 0 下面证明对数正态分布是重尾分布 证明：对于v r 0 ，计算蹦，x 为对数正态分布 硝：f 。耋一p 一1 一p “出 , i o 4 2 x a x ：志e 一呼抛+ f 志e 一呼抛 垒厶+ 厶厶 下面考虑厶的敛散性： 厶2 f 忐，可抛 m z ：ff ”当p 可p r , d t 5 2 2 2 。一j o 4 2 ；r o 丝e 五杀p 一万p 班垒c 厂( z 胁 因为l i mf ( z ) = ，所以f 二厂( z ) 出发散，即l 发散，从而有歧= ，即对数 z o 一户 正态分布为重尾分布 定义3 2 称q = 茁：e x 声= f x s d f ( 工) 0 为轻度重尾分布族称 k = 盯甄为重度重尾分布族 定义3 3 称分布疋控制分布只，若s u p 露( 石) 瓦 ) 】 0 0 ，记作露 瓦( 或 相应地x ， x ：) 若霞 瓦“茸，则称一与最弱等价若! 驸瓦( 工) 瓦( 硼= 1 ， 则称e 与疋等价 便于表述，记 d r - c = f ：对明茁。，露 o 成立 1 4 引理3 2 对数正态分布属于d 证明：l i m e z f f ( x ) = l i mp 缸( 1 一，o ) ) ：l i r a 1 - f ：( 一x ) x - - 0 0 p “ ：l i m 坐! - j - - o o t e 7 “ 1 ( i n x - p ) 三e 一可 ；l i m 生遮盟一 ” 2 e 。4 4 一l 等e = = = = 一 口 i n x ：f l i m 堂边丝；- 一 。f 。旯p 一埘 = 姆丽1e x p 一警叫 由d x 。的定义，得对数正态分布可以控制所有轻尾分布 当然对于轻度重尾分布，根据d 两的定义及( 3 2 ) 式也可以找到分布，( x ) 来控 制它们也可以通过在对数正态分布上乘以一个常数来控制它们 本文运用了对数正态分布的重尾性质，在处理尾部较高的分布的抽样时，以对 数正态分布作为其尾部的包络，操作起来更加方便 图3 是n ( o ，1 ) 和标准对数正态分布的图像，可以清楚地看到标准对数正态分布 的尾部明显高于n ( o ，1 ) 的尾部 ! 0 6 f j 0 j - 口4 7 心、| 、 d 3 ，0 2 一 层 7 一 图3n ( o ，1 ) 和标准对数正态分布的图像 6 下面分两节来讨论具体的抽样过程 设目标函数为f ( x ) ，其定义域为d ，分别从抽样区间是有限和无限出发，分情 况加以讨论 3 2 1 抽样过程 3 。2 有限区间上的抽样 d 为有限区间，不妨设为陋，6 】，a ，b ，则对于多峰目标函数，( x ) ，为了运用 拒绝抽样法进行抽样，关键是要找至：l jf ( x ) 的包络函数 首先，求出f ( x ) 在区间 a , b 】上的所有极大值点，通常厂o ) 有几个极大值点就表 明它有几个峰不妨设f ( x ) 在区间【a , b 】上有k 个极大值点，按从4 , n 大依次为： x 2 x k ：接着，以直线依次连接各极大值点，得到一个线性分段函数g ( x ) ， 1 6 g ( z ) = 丝丝( x 一工1 ) + m ，工【口，x 2 屯一五 三丝二监伍一再) + 乃，x e x , ，而+ i 】，2 i 0 ，一3 x - 2 ，f 严格增 = 0 ，x = 一2 ，厂( 一2 ) 达极大值 0 ，一2 x 0 ，0 x 2 ，厂严格增 = 0 ，石= 2 ，厂( 2 ) 达极大值 0 ，2 工3 ，俨格降 得到，( z ) 在 一3 ，3 的两个极大值点： j c l = - 2 ，x 2 = 2 ， 且有，( 2 ) = f ( - 2 ) = 2 7 0 9 7 ，计算g ( x ) ，得： g ( 石) = 芴7 0 ，x 【- 3 ，3 】 好胁肛1 _ ，4 0 铽2 菇7 1 帅 3 】蝴均鼢布 取m = 【百1 4 0 】+ l = 2 ，抽样过程如下： 抽取t i u ( o ，1 ) ，r 2 一u ( o ，1 ) 且，i 与r 2 相互独立，有变换法 ： x 垒6 5 3 u ( - 3 ，3 ) ； 计算g + ( z ) ，若吒f ( x ) 2 9 + 0 ) ，则x + f ( x ) ，否则拒绝z 重复上述两步，直至得到所求的样本数 显然，例3 1 运用本文的方法进行抽样时效率是比较高的m = 2 ，平均抽样两 次，就能获得所需样本，且操作过程简单易行，下面再来看无限区间上的情况 3 3 1 抽样过程 3 3 无限区间上的抽样 设目标分布为厂( 功，定义域d 为无限区间，不妨设d 为( - o o ，o o ) 首先，和有限区间上的抽样样，求出f ( x ) 的所有极大值点仍假设f ( x ) 有k 个 极大值点，按从小到大依次为：五 x 2 五 对于区间【五，】，和有限区间上的情况类似，将f ( x ) 的各极大值点依次以直线 相连就得到线性分段函数g ：( 功，使得9 2 ( x ) ，( 工) ，工h ，稚】接着，选取两个点而 和k + i ( 硼，+ i o o ) ，使得x o ，9 2 ( x o ) f ( x o ) ，9 2 ( x , + i ) f ( x k + 1 ) ， 这样，就将9 2 ( x ) 扩充到了区间，吒+ 。】上下面再分别对区间( 哪，：c o 】和 小0 0 ) 加以 处理 先来看右端的尾部：k + ，0 0 ) 直接用对数正态分布来覆盖它取对数正态分布中 相应的 x k + o 。) 部分作为f ( x ) 右端尾部的包络岛( x ) 记对数正态分布的密度函数为 p ( 功，r f l 1 jf ( x ) 右尾 垓小0 0 ) 的包络为：g ，( x ) = c p ( x ) ，其中c 为某个常数，通常取c 1 再来看厂o ) 左端的尾部：( 咖，】由于p ( x ) 的取值范围为( 0 ，0 0 ) ，因此不能直 接用其尾部覆盖f ( x ) 左端尾部作一个简单的变换：求p ( x ) 关于少轴的对称函数， 记为p ( 功，再用p ( 石) 来覆盖f ( x ) 左端的尾部则有f ( x ) 左端尾部的包络： g j ( x ) = c ( z ) ，其中c 为某个常数，通常取c 1 如此，就得到了一个分段函数g ( x ) ： f & ( z ) ，z ( ，而】 g ( 石) = 岛( x ) ，x 【x o ，x t + i 】 【岛( x ) ，z 五+ i ，) 通常情况下，由于r 。g ( x ) 出1 ，不能直接作为抽样分布，需要找一个抽样分 布g + ( z ) 和“覆盖常数m ，使得m g ( z ) f ( x ) ，x d ，然后再进行抽样 分别计算： 二g ( 工) a x 垒p l ，r g ： ) a x 垒p ：，e 9 3 ) a x 兰p 3 ， k l 算p 2 = l p l 一见令，9 2 。( x ) = 旦g ：( 工) 皇r e , 2 ( x ) ，这样就得到一个抽样分布g ( 工) ： 2 f g i ( x ) ，x x k “ 只要取m = 三】+ l ，就可以保证t _ m g + ) f ( x ) ，工d 了 令a 。= p ，a ：= p l + 五，a ，= l ，下面运用拒绝抽样便可得到所需样本了 抽样过程： 抽取吒u ( o ，1 ) ，u ( o ，1 ) 且，i 与相互独立； i 蜀 ) ，r 2 a 1 产生， 9 2 ( 功，a ， 7 m = 【三】+ 1 = 3 ， 口l = p l = 0 0 2 5 8 ，a 2 = p i + p ：= 0 9 7 4 2 ，a 3 = 1 抽样过程： 抽取巧u ( o ，1 ) ，r z u ( 0 ，1 ) 且，i 与眨相互独立； i g l ( 工) ，吒4 l 产生f 9 2 ( x ) ，q r 2 a 2 ，计算g ( ? ) ； 【9 3 ( z ) ，口2 吒1 若巧f ( x ) 3 g ( 工) ，则工f ( x ) ，否则拒绝工 在例3 2 中，m = 3 ，效率还是比较高的g f l 3 ：j ：p , = p 3 = o 0 2 5 8 比较小，抽样 区间绝大部分落在区间 - 7 ，7 】上，i 7 i i - 7 ，7 】的抽样密度g ：( z ) 是一个分段线性函数， 显然抽样时更简单，从而整体上简化了抽样过程 3 4 结论 可以看到本文的方法可以快速有效地从一类复杂分布，尤其是多峰分布中抽取 样本运用拒绝抽样时，采用的包络不仅形式简单，还易于抽样，同时也将“覆盖函 数 m 的取值，控制在最小的范围，以保证效率 文中主要考虑了分布为连续情形下的函数的抽样，事实上，在离散和奇异情形 下这种方法同样适用也是先求出“极大值点”，只是这时的“极大值点 ，与连续情 形下定义的“极大值点 有所差别，方便起见，仍称为“极大值点” 仍然是以依次连接各极大值点得到的分段线性函数作为有限区间上的包络，尾 部也是采用对数正态，及对数正态关于y 轴的对称函数，过程类似 尽管从构造过程来看本方法比较简单，但如果目标函数存在某个或某几个比较 大的峰值时，所得9 2 ( z ) 就会偏离x 轴较远，为覆盖住目标函数，所取的m 的值也会 越大，从而降低效率另外，当目标分布厂0 ) 为重度重尾分布时，和处理轻度重尾 目标分布一样，也可以通过对对数正态分布乘以一个常数来处理，只是需乘上一个 较大的常数才能控制目标分布的尾部，从而影响抽样效率但是，多数情况下，本文 的方法都能以较高的效率获得所需的样本 第四章讨论 本文给出了一个处理复杂密度函数尤其是多峰分布函数的抽样方法，通过运用 特殊的包络函数，不仅简化了抽样过程，降低了复杂性，还提高了效率，而且操作 起来也极为方便，比e es a m p l e r 简单得多 特别是当目标分布的峰比较多时，e es a m p l e r 通常需要划分更多的能量层，并 对各能量层分别采用m h 法进行抽样，随着峰的个数的增加，计算量也随之大幅增 加由于e es a m p l e r 对各能量层内分别使用m h 算法，且各层内的目标函数通常是 不一样的，又因为m h 算法通常需要很长一段时间才能使得样本收敛，大大降低了 效率而本文的方法在有限区间情形下，仅用一个分段线性函数作为包络就可以了， 显然要比m h 法简单得多 此外，e es a m p l e r 在其使用过程时有很多经验的地方，例如，能量层划分和温 度序列的选取，等能量跳跃概率p 。的选取等对这些量的不同取法，常会影响抽样 的结果本文的方法不需要作任何设定，除了目标分布本身外，不存在任何其它参数 而且本文的方法获得的是独立样本，而e es a m p l e r 所得的却不是独立样本 对于本方法的效率，前文中曾提到，如果目标函数存在某个或某几个比较大的 峰值时，所得g ：( 曲就会偏离x 轴较远，为覆盖住目标函数，所取的膨的值也会越大， 从而影响抽样效率这时，可对包络函数的取法作适当调整目标函数的尾部仍用对 数正态分布来覆盖，对中间部分的区间求出拐点，划分出凹凸区间，凸区间用适应 拒绝抽样，凹区间以连接拐点的直线段作为包络，这样处理效率会大大提高 另外，对于多元情形，本方法同样适用，只是包络函数的形式比较复杂寻求更 有效、更简便的包络函数形式，还有待进一步的探讨同样，在处理重度重尾目标分 布时，如果仍用对数正态分布作为目标分布尾部的包络，很有可能会影响抽样效率， 为提高效率，可考虑利用某个重度重尾分布作为包络函数，这也是一项需要进一步 研究的内容 参考文献 【1 】 b rj a ns c a f f o ，j a m e s0 b o o t ha n d 焘c d a v i s o n ，e m p i r i c a ls u p r e m u mr e j e c ti o ns a m p li n g ，b i o m e t r i k a ，2 0 0 2 ，8 9 ( 4 ) ：7 4 5 7 5 4 【2 】e a s e ll a ，ga n dr o b e r t ，c p ，p o s t p r o c e s si n ga c c e p t r e j e c ts a m p l e s ： r e c y c l i n ga n dr e s c a l i n g ，s t a t i s t ，1 9 9 8 ，7 ：1 3 9 1 5 7 3 】3 d u a n e ，s ，k e n n e d y ，a d ，p e n d l e t o n ，b j a n dr o w e t h ，d ，h y b r i dm o n t e c a r l o ，p h y s l e t t b ，1 9 8 7 ，1 9 5 ：2 1 6 2 2 2 【4 】e v a n s ，m a n dt 。s w a r t z ，r a n d o mv a t i a b l eg e n e r a t i o nu s i n gc o n c a v i t y p r o p e r t i e so ft r a n s f o r m e dd e n s i t i e s ，j o u r n a lo fc o m p u t a t l e n a la n d g r a p h i c a ls t a t is t i c s ，1 9 9 8 ，7 ( 4 ) ：5 1 4 5 2 8 5 】5j h a h r e n sa n dk d k o h r t ，c o m p u t e rm e t h o d sf o rs a m p li n gf r o ml a r g e l y a r b i t r a r ys t a t is t i c a ld is t r i b u t i o n ，c o m p u t i n g ，1 9 8 1 ，2 6 ：1 9 3 1 【6 】6j o s e fl e y d o l d ，w o l f g a n gh 6 r m a n n ，o n i v e r s a la i g o r i t h m sa sa na l t e r n a t i v ef o rg e n e r a t i n gn o n u n i f o r mc e n t i n u o u sr a n d o mv a r i a t e s ，h t t p ： s t a t m a t h w u - w i e n a c a t ，2 0 0 2 7 】l u ns l i u ，m o n t ec a r l os t r a t e g i e si ns c l e n t i f i cc o m p u t i n g ，s p r i n g e r ， 2 0 0 1 8 】l i u ，j s ，m e t r o p li z e di n d e p e n d e n c es a m p li n g w it h c o m p a r is o n s t o r e j e c t i o ns a m p li n ga n di m p o r t a n c es a m p li n g ，s t a t is t ，1 9 9 6 ，6 ：13 1 1 9 【9 】9 m i r a ，a ，m o l l e r ，j a n dr o b e r t s ，g ，p e r f e c ts 1 i c es a m p l e r ，j r s t a t s e e s e t bs t a t m e t h o d 0 1 ，2 0 0 1 ，6 3 ，5 9 3 6 0 6 【10 】n e a l ，r m ，s 1 i c es a m p l i n g ，a n n s t a t i s t ，2 0 0 3 ，3 1 ：7 0 5 - 7 6 7 【1 1 】r o b e r t s ，g a n dr o s e n t h a i ，j s ，c o n v e r g e n c eo f s 1i c es a m p l e rm a r k o v c h a i n ，j r s t a t s o c s e r bs t a r m e t o d o l ，1 9 9 9 ，6 1 ：6 4 3 6 6 0 【1 2 】s c k o u ，q i n g z h o ua n dw i n gh u n g w o n g ，e q u i e n e r g ys a m p l e rw it h a p p l i c a t i o n s i n ss t a t is t i c a li n f e r e n c