龚光鲁,钱敏平著应用随机过程教程及其在算法与智能计算中的应用 清华大学出版社,2003 第2章随机样本生成法 1一维随机数 随机变量(或随机向量)的样本简称为随机数.由于在统计中常用的是独立样本列, 所以我们不妨假设随机数之间都是独立的.生成随机数的方法,也称为随机数的取样法 (Sampling ) 1.1均匀随机变量的计算机模拟 定义2.1在,1上均匀分布的随机变量的独立样本称为均匀随机数(U[0,随机数) 在计算机上产生的称之为”伪随机数”的数列,是一种具有非常长周期的,且能通过数理 统计中的独立性与均匀性假设检验的数列.实践证明伪随机数是均匀随机数的一种可行的 近似这种伪随机数虽然并不是独立同分布的U[O]随机变量的样本,而是在[0,1中取值 的周期数列,但是由于它可以像均匀随机数一样地通过数理统计中的独立性与均匀性假设 检验,而且它的周期非常长,以至在计算机实际运算过程中不会出现重复,所以在实际计算 中它能很好地替代均匀随机数 最普遍用以产生伪随机数的方法是同余法.典型的例子如下 ),y=1,xn=yn2-36(周期约为21010) n+1 5y, (mo ),y=1, (周期约为102) yn2=yn+yn(mod2),y=0,y1=1,x,n=yn2-(周期约为104) (关于伪随机数,可参见:现代数学手册,随机数学卷,第10篇,孙嘉阳,石坚,丛树 铮,徐映波编蒙特卡罗法.华中科技大学出版社,2000年). 1.2分布函数F(x)的随机数 命题2.2(反函数方法)分布函数为F(x)的独立随机变量列的样本,称为F(x)随机 数若F(x)严格单调递增,ξ是均匀随机数,则F()是F(x)随机数,其中F-为F的反 函数 证明P(F(5)≤x)=P(5≤F(x)=F(x) 命题2.3设随机变量ξ只取有限个值,其分布为 把[0,1 分为n个不交子区间,使第i个区间J的长度为P,任取均匀随机数U,则
33 龚光鲁，钱敏平著 应用随机过程教程及其在算法与智能计算中的应用 清华大学出版社, 2003 第２章 随机样本生成法 1 一维随机数 随机变量（或随机向量）的样本简称为随机数. 由于在统计中常用的是独立样本列， 所以我们不妨假设随机数之间都是独立的．生成随机数的方法，也称为随机数的取样法 (Sampling）． 1. 1 均匀随机变量的计算机模拟 定义２.１ 在[0,1]上均匀分布的随机变量的独立样本称为均匀随机数(U[0,1] 随机数). 在计算机上产生的称之为”伪随机数”的数列, 是一种具有非常长周期的, 且能通过数理 统计中的独立性与均匀性假设检验的数列. 实践证明伪随机数是均匀随机数的一种可行的 近似. 这种伪随机数虽然并不是独立同分布的U[0,1] 随机变量的样本, 而是在[0,1]中取值 的周期数列, 但是由于它可以像均匀随机数一样地通过数理统计中的独立性与均匀性假设 检验, 而且它的周期非常长, 以至在计算机实际运算过程中不会出现重复, 所以在实际计算 中它能很好地替代均匀随机数. 最普遍用以产生伪随机数的方法是同余法. 典型的例子如下: yn+1 = 5 13 yn (mod 2 36 ), y0 = 1, xn = yn × 2 -36 (周期约为2 ×1010）； yn+1 = 5 17 yn (mod 2 42 ), y0 = 1, xn = yn × 2 -42 (周期约为1012）； 2 1 44 0 1 44 周期约为 1014） 4 1 = + (mod 2 ), = 0, = 1, = × 2 ( × - n+ n + n n n y y y y y x y ． （关于伪随机数， 可参见：现代数学手册，随机数学卷，第 10 篇，孙嘉阳， 石坚，丛树 铮，徐映波编 蒙特卡罗法．华中科技大学出版社，2000 年）． 1. 2 分布函数F( x) 的随机数 命题２.２（反函数方法）分布函数为F( x) 的独立随机变量列的样本，称为 F( x) 随机 数．若F( x) 严格单调递增, x 是均匀随机数, 则 (x ) -1 F 是 F( x) 随机数, 其中 -1 F 为 F 的反 函数. 证明 ( ( ) ) ( ( )) ( ) 1 P F £ x = P £ F x = F x - x x . ? 命题２．３ 设随机变量x 只取有限个值，其分布为 ÷ ÷ ø ö ç ç è æ n n p p x x L L 1 1 x ~ . 把[0,1] 分为n 个不交子区间, 使第i 个区间 i J 的长度为 pi . 任取均匀随机数U , 则
1(U) 就是一个5随机数(它的意思是:如果U落入J,就取ξ=x1) 在统计再抽样中的应用 在样本组中再抽样,或者由样本作的参数估计代替分布中的未知参数后,所得到的分布 的随机取样,统一称为 Bootstrap方法.具体地说有如下两种方法 (1)非参数 Bootstrap方法.设自一个未知方差的分布取样X1…,x(不是计算 机模拟取样,而是人工取样)如果要作方差的区间估计,就需要知道方差估计σ的方差 lar().一般Vamr(σ)很不好求,需要对它用再抽样进行估计.为此可将样本分布 X 1作为离散随机变量的分布,独立地取样N次,每次独立地取样m个设从 第k次的m个样本值得到方差的估计σk(k≤N),将此N个的平均记为G,最后用 N-12(0k-0)2估计ar(a) 此法可以用于一般未知参数的方差估计 (2)参数 Bootstrap方法.设自一个带有未知参数(91,…9)的分布p(x,912…9) 的样本X1…,Xn(不是计算机模拟取样,而是人工取样)得到未知参数的估计(91…,9) 后,对分布Pp(x,91,……,91)用计算机模拟取样.独立地取样N次,每次独立地取样m个.其 它与(1)相同 注意,计算机模拟取样只能对已知的分布施行,对于含未知参数的分布,只能作普通的 人工取样.以上的两种再抽样方法,补充了人工取样采样量的限制.因为计算机模拟取样既 快速又经济 1.3正态随机数 N(,1)随机数称为标准正态随机数.生成标准正态随机数有一个比反函数的方法更为 简单的实践方法,就是利用中心极限定理设n,…,n2为均匀随机数(它们是独立的),由 中心极限定理,可以认为=n1+…+n12-6≈N(01),即用5=n1+…+n2-6近似地 作为标准正态随机数.在实际计算中n,(1≤i≤12)们还应该用伪随机数代替 命题2.4(生成标准正态随机数的 Box-Muller方法)取两个独立的均匀随机数
34 å= n i xi I J Ui 1 ( ) 就是一个x 随机数（它的意思是：如果U 落入 i J ，就取 i x = x ）． 在统计再抽样中的应用 在样本组中再抽样，或者由样本作的参数估计代替分布中的未知参数后，所得到的分布 的随机取样，统一称为 Bootstrap 方法．具体地说有如下两种方法 （１）非参数 Bootstrap 方法．设自一个未知方差的分布取样 X Xn , , 1 L （不是计算 机模拟取样，而是人工取样）如果要作方差的区间估计，就需要知道方差估计 2 Ù s 的方差 ( ) 2 Ù Var s ．一般 ( ) 2 Ù Var s 很不好求，需要对它用再抽样进行估计．为此可将样本分布 ÷ ÷ ø ö ç ç è æ n n X Xn 1 1 1 L L 作为离散随机变量的分布, 独立地取样N 次，每次独立地取样m 个．设从 第k 次的m 个样本值得到方差的估计 ( ) 2 k k £ N Ù s ，将此 N 个的平均记为 2 Ù s ，最后用 2 2 2 1 2 ( ) 1 1 Ù Ù = Ù - - = å Ù s s k s N N k 估计 ( ) 2 Ù Var s ． 此法可以用于一般未知参数的方差估计． （２）参数 Bootstrap 方法． 设自一个带有未知参数( , ) J1 LJl 的分布 ( , , ) 1 l p x J LJ 的样本 X Xn , , 1 L （不是计算机模拟取样，而是人工取样）得到未知参数的估计( l ) Ù Ù J１，L，J 后，对分布 p x, l ) Ù Ù （ J１，L，J 用计算机模拟取样．独立地取样 N 次，每次独立地取样 m 个．其 它与（１）相同． 注意，计算机模拟取样只能对已知的分布施行，对于含未知参数的分布，只能作普通的 人工取样．以上的两种再抽样方法，补充了人工取样采样量的限制．因为计算机模拟取样既 快速又经济． 1. 3 正态随机数 N(0,1) 随机数称为标准正态随机数． 生成标准正态随机数有一个比反函数的方法更为 简单的实践方法, 就是利用中心极限定理. 设 1 12 h ,L,h 为均匀随机数（它们是独立的）, 由 中心极限定理，可以认为 6 (0,1) x =h1 +L+h12 - » N , 即用x =h1 +L+h12 -6 近似地 作为标准正态随机数. 在实际计算中hi (1 £ i £ 12) 们还应该用伪随机数代替. 命题２.４ (生成标准正态随机数的 Box-Muller 方法) 取两个独立的均匀随机数
n, n2 2In n, cos(2r,) 52=y-2hn1sn(2x) 则51,2为相互独立的标准正态随机数 证明令1,n2~U[0且独立,则1-1,72也是独立的U,]随机变量.于是由命题 22,p=F(-n1)=√-2hn1~分布函数F(r)=1-e2,9=2n2~U[0,2],且 相互独立.而51=Pco9,52=pSin9,易见它们是独立的标准正态随机数 命题2.5(生成标准正态随机数的 Marsal ia方法)设(X,F)为单位圆上的均匀 随机数.则 In(X X+y 8)0) (提示将直角坐标(x,Y)转换为极坐标(R,9)) 一般正态随机数的生成若为标准正态随机数,则显见G2+为N(μ,a2)随机数 1.4 Poisson随机数 下述结论给出了利用伪随机数生成 Poisson随机数的方法。 命题2.6设U1,U2,…是相互独立的,1均匀随机数.若 Uk<es∏U,则定义N=n·那么N~ poisson 证明令T=-hUk~eXp1在指数流与 Poisson过程的关系(参见第3章)中取参 数为1,取时间t为λ即得 1.5混合分布随机数 对于权重为P12…,Pn(和为1的n个正数)的混合分布随机数,我们有 命题2.7设U~U[0110=10<…<tn=1,t1-1-1=P1(i≤n),F(x)i≤n)为分 布函数,那么 ∑F U--)l-(U)~分布函数为∑pF的混合分布 证明
35 1 2 h ,h , 令 2ln cos(2 ) 1 h1 ph2 x = - , 2ln sin( 2 ) 2 h1 ph2 x = - . 则 1 2 x ,x 为相互独立的标准正态随机数. 证明 令 , ~ [0,1] h1 h2 U 且独立, 则 1 2 1-h ,h 也是独立的U[0,1] 随机变量. 于是由命题 ２.2, (1 1 ) 2 ln 1 ~ 1 r = -h = - h - D F 分布函数 2 2 ( ) 1 r F r e - = - , 2 ~ [0 2 ] J ph2 U ，p D = , 且 相互独立. 而x1 = r cosJ,x2 = r sin J ，易见它们是独立的标准正态随机数. 命题２.５ （生成标准正态随机数的 Marsaglia 方法） 设(X,Y ) 为单位圆上的均匀 随机数. 则 ÷ ÷ ø ö ç ç è æ + - + =÷ ÷ ø ö ç ç è æ D Y X X Y X Y 2 2 2 2 2ln( ) h x ÷ ÷ ø ö ç ç è æ ÷ ÷ ø ö ç ç è æ ÷ ÷ ø ö ç ç è æ 0 1 1 0 , 0 0 ~ N . （提示 将直角坐标(X,Y ) 转换为极坐标(R,J) ）. 一般正态随机数的生成 若x 为标准正态随机数, 则显见sx + m 为 ( , ) 2 N m s 随机数. 1. 4 Poisson 随机数 下述结论给出了利用伪随机数生成 Poisson 随机数的方法。 命题２.６ 设U1 ,U2 ,L是相互独立的[0,1]均匀随机数． 若 Õ Õ + = = -l < £ 1 1 1 n k n k U k e Uk ， 则定义 N = n ． 那么 N ~ Poissonl . 证明 令 1 ln ~ exp Tk =- Uk D . 在指数流与 Poisson 过程的关系 ( 参见第 3 章) 中取参 数为 1, 取时间t 为l 即得. 1. 5 混合分布随机数 对于权重为 p pn , , 1 L (和为 1 的n 个正数) 的混合分布随机数, 我们有 命题２.７ 设 ~ [0,1],0 1, ( ) U U = t 0 <L < t n = t i - t i -1 = pi i £ n , F (x)(i n) i £ 为分 布函数, 那么 å= - - - - n i t t i i i I U p U t F i i 1 ( , ] 1 1 ( ) ( ) ~ 1 分布函数为 i i n i å p F =1 的混合分布. 证明
P1-41(0)≤x) P(F )1(U)≤x,U∈(1-,1]) Pi ∑P(-1<U≤1+pF(x)=∑pF(x) 在实际计算中,应该用伪随机数来取代均匀随机数U,如果取到的伪随机数落在(t-1,l1 中,则取F(U-b ),这样得到的数就是∑PF混合分布随机数 这个方法常用在排队论中.在那里的典型情形是混合指数分布(有的书上称为超指数分 布)即F(x)=1-c4,此时简单地有F(y)=-1-m(1-y),于是计算变得非常简单 入 而有效(当然也可利用命题22通过反函数来得到混合指数分布随机数.但是计算量会增加 很多,因为这个反函数并不简单 1.6 Von Neuman取舍原则 Von Neuman取舍原则 假定我们要生成密度为p(x)的随机数.为此取一个参考分布密度P0(x),使它满足 (1)p0(x)随机数容易生成,例如P0(x)为正态密度,均匀密度,指数密度,及它们的 混合密度; (2)po(x)与p(x)的取值范围差不多,且存在C,使p(x)≤C·P0(x) 那么,我们有 命题2.8设随机变量η具有密度P0(x),而随机变量U~U/[0,且与n独立,则 P(n≤x p() ≥U)=p(vh CP0(7) 证明对的取值用推广了的全概率公式(P(4)=「P(4|n=xg(x)dk),得到 Pn≤x,PO ≥U)「PU≤ pl) Cp()%0()d p(y)dy 左 ))∫PUsp) p() poG 取舍原则( Rejection Principle)的具体作法是 (1)独立地生成n个独立的P0(x)随机数n,…n与n个与之独立的U/[0,随机数
36 £ = - å - = - - ( ( ) ( ) ) ( , ] 1 1 1 1 I U x p U t P F i i t t n i i i i ( ( ) ( ) , ( , ]) ( , ] 1 1 1 1 1 t t i i i i i n i I U x U t t p U t P F i i - - - = £ Î - å - ( ( )) ( ) 1 1 P t 1 U t 1 p F x p F x n i i i n i å i i i i å = = = - < £ - + = ． ? 在实际计算中, 应该用伪随机数来取代均匀随机数U , 如果取到的伪随机数落在 ( , ] i 1 i t t - 中, 则取 ( ) 1 1 i i i p U t F - - - , 这样得到的数就是å= n i piFi 1 混合分布随机数. 这个方法常用在排队论中. 在那里的典型情形是混合指数分布(有的书上称为超指数分 布), 即 x i i F x e -l ( ) = 1- , 此时简单地有 ln(1 ) 1 ( ) 1 F y y i i - l = - - , 于是计算变得非常简单 而有效 (当然也可利用命题 2.2 通过反函数来得到混合指数分布随机数. 但是计算量会增加 很多, 因为这个反函数并不简单). 1. 6 Von Neuman 取舍原则 Von Neuman 取舍原则： 假定我们要生成密度为 p( x) 的随机数. 为此取一个参考分布密度 ( ) 0 p x , 使它满足: (1) ( ) 0 p x 随机数容易生成, 例如 ( ) 0 p x 为正态密度, 均匀密度, 指数密度, 及它们的 混合密度； (2) ( ) 0 p x 与 p( x) 的取值范围差不多, 且存在C ，使 ( ) ( ) 0 p x £ C × p x ． 那么，我们有 命题 2. ８ 设随机变量h 具有密度 ( ) 0 p x , 而随机变量U ~ U[0,1] 且与h 独立, 则 ò -¥ £ ³ = x U p v dv Cp p P x ) ( ) ( ) ( ) ( | 0 h h h . 证明 对h 的取值用推广了的全概率公式（ ò P(A) = P(A |h = x)g(x)dx ），得到 左 = = 右 £ £ = ³ £ ³ = ò ò ò ò ¥ -¥ -¥ ¥ -¥ -¥ p y dy C p y dy C p y dy Cp y p y P U p y dy Cp y p y P U U Cp p P U Cp p P x x x ( ) 1 ( ) 1 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( , 0 0 0 0 0 0 h h h h h ．? 取舍原则(Rejection Principle)的具体作法是: (1) 独立地生成n 个独立的 ( ) 0 p x 随机数h hn , , 1 L 与n 个与之独立的U[0,1] 随机数 U U n , , 1 L
(2)对于1=12,…,如果有P(n) CPo(n,) ≥U1,就保留n1,否则就舍弃n, 由命题2.8,所有这样保留下来的那些n,们就成为一系列独立的p(x)随机数(当然个数会 比n小很多)这种取舍方法称为 Von Neuman取舍原则 取舍原则可以改良为 如果p(x)=y·h(x),只要存在C,使h(x)≤C·p0(x),那么我们可以在取舍原则中 用h(x)代替p(x),得到p(x)随机数.具体为:独立地生成n个独立的P0(x)随机数 h,…,n与n个与之独立的独立UO随机数U1,…,Un,如果n)2U,则保留 0() 否则舍弃n,那么所有保留下的是相互独立的p(x)随机数 证明只需注意到这时有p()y:C.p(x),并且一P(x)=如()即可 yCPo(x) Cp(x) 注1一般地γ需通过 1=[h(x)dx计算,其中的积分不易计算但是上面的事实说 明不必计算γ,即可以忽视这个常数因子.这就使取舍原则变得非常好用.取舍原则的优 点是简单易行,可以适用于非常复杂的分布 注2如果p(x)只在有限区间[a,b上不等于零,而且有界,那么P0(x)就可取均匀分 布U[a,b];如果p(x)只在右半直线不等于零,那么指数分布就可以是p0(x)的一个选择 如果p(x)在实直线上不等于零,且分布密度的尾部不大,则正态分布就可以是P0(x)的一 个选择,如果p(x)在实直线上不等于零,且分布密度的尾部较大(重尾分布),则t分布 就可以是P(x)的一个选择,如果p(x)具有多个峰,则混合正态分布或混合指数分布就可 以是p0(x)的一个选择可见适当精心地选取P0(x)是使计算省时的关键 注3原则上取舍原则也适用于离散分布和多维密度,但是在多维密度的情形,Po(x) 的选取并不容易 注4样本生成的一个重要应用,就是对于函数的一些积分,用 Monte Carlo方法(随机 模拟算法)给出估计·直观地可以猜测,采用的随机数的密度P。(x)的形状与被积函数越像, 则估计的方差会越小,即效果越好.这种取样法称为重要度采样( Impotance Samling)(确切 的定义与方差的最小性证明可参见第8章)
37 (2) 对于i = 1,2,L, 如果有 i i i U Cp p ³ ( ) ( ) 0 h h , 就保留hi , 否则就舍弃hi 。 由命题 2.８, 所有这样保留下来的那些hi 们就成为一系列独立的 p( x) 随机数(当然个数会 比n 小很多). 这种取舍方法称为 Von Neuman 取舍原则. 取舍原则可以改良为： 如果 p(x) = g × h(x) , 只要存在C ，使 ( ) ( ) h x £ C × p0 x , 那么我们可以在取舍原则中 用 h(x) 代替 p( x) ，得到 p( x) 随机数． 具体为：独立地生成n 个独立的 ( ) 0 p x 随机数 h hn , , 1 L 与n 个与之独立的独立U[0,1] 随机数U U n , , 1 L ，如果 i i i U Cp h ³ ( ) ( ) 0 h h , 则保留hi , 否则舍弃hi , 那么所有保留下的是相互独立的 p( x) 随机数. 证明 只需注意到这时有 ( ) ( ) p x £ g ×C × p0 x , 并且 = × ( ) ( ) 0 Cp x p x g ( ) ( ) 0 Cp x h x 即可． 注１ 一般地g 需通过 ò = h(x)dx 1 g 计算, 其中的积分不易计算. 但是上面的事实说 明不必计算g , 即可以忽视这个常数因子．这就使取舍原则变得非常好用．取舍原则的优 点是简单易行, 可以适用于非常复杂的分布. 注２ 如果 p( x) 只在有限区间[a, b]上不等于零, 而且有界, 那么 ( ) 0 p x 就可取均匀分 布 U[a,b] ; 如果 p( x) 只在右半直线不等于零, 那么指数分布就可以是 ( ) 0 p x 的一个选择; 如果 p( x) 在实直线上不等于零，且分布密度的尾部不大, 则正态分布就可以是 ( ) 0 p x 的一 个选择; 如果 p( x) 在实直线上不等于零，且分布密度的尾部较大（重尾分布）, 则t 分布 就可以是 ( ) 0 p x 的一个选择; 如果 p( x) 具有多个峰, 则混合正态分布或混合指数分布就可 以是 ( ) 0 p x 的一个选择. 可见适当精心地选取 ( ) 0 p x 是使计算省时的关键. 注３ 原则上取舍原则也适用于离散分布和多维密度，但是在多维密度的情形， ( ) 0 p x 的选取并不容易． 注４ 样本生成的一个重要应用, 就是对于函数的一些积分, 用 Monte Carlo 方法(随机 模拟算法)给出估计． 直观地可以猜测, 采用的随机数的密度 ( ) 0 p x 的形状与被积函数越像, 则估计的方差会越小, 即效果越好. 这种取样法称为重要度采样 (Impotance Samling)(确切 的定义与方差的最小性证明可参见第 8 章)
1.7 Gamma随机数与Beta随机数的生成 在a≥1时,设[a]是正数a的整部,可以证明r(a,A)分布的密度与r([a]4)分布的密 度的比例是有界函数.由此可以设计一个由指数分布随机数,通过 Von neuman取舍原则得 到r(a,)分布的随机数的方法 在a,B≥1时,设m,2分别为独立的r(a,,I(B)随机数,则可以证明5=n1+n2 为B(a,B)随机数 关于0<a<1时r(a,)分布的取样方法,可以参见文献[LK] 2.多维随机数 2.1连续型多维随机数 对于已知的分布密度,我们可以利用条件密度,把生成多维随机数归结为生成一系列 一维随机数.设随机向量(X1…,Xa)的密度为∫(x12x2…,x).那么,我们有表达式 ∫(x1,x2,…,x4)=fx1(x1)(x2|x1)…f(x|x1,…,x), 其中∫x为X1的边缘密度,f(x4|x1,…,x-1)为在已知X1=x1,…,Xk=x1条件下 X的条件密度于是可以先取一个∫x随机数x1;然后,在x1固定的情形下,生成 ∫(·|x1)随机数x2;再在x1,x2固定的情形下,生成一个f(|x1,x2)随机数x3;…最后 在x1…,x1固定的情形下,生成一个∫(|x12…,x)随机数x,这样得到的 (x,…x)就是随机向量(x1,…,Xd)的一个随机数.当然在生成各个条件密度的随机数 时,仍然可以使用 Von neuman取舍原则.注意在用取舍原则于一维分布取样时,可以忽 略一个因子 然而,在d较大时,更为常用的是Gbs取样法( Gibbs Sampler),它是一种基本的动态 Monte carlo方法,即 Markov链 Monte Carlo方法,在本书第8章中我们将介绍这种方法 在很多实际情形中,多维密度常常并不知道,更无法知道各个条件密度.这时最自然而 粗糙的想法是用条件频率来近似条件密度 2.2离散型多维随机数 只要用概率函数p(x1,…,x)=P(X1=x1,…,Xd=x)代替密度函数,2.1段中的办法 就自然适用 2.3多维正态随机数
38 1. 7 Gamma 随机数与 Beta 随机数的生成 在a ³1时, 设[a]是正数a 的整部, 可以证明G(a,l) 分布的密度与G([a],l) 分布的密 度的比例是有界函数. 由此可以设计一个由指数分布随机数, 通过 Von Neuman 取舍原则得 到G(a,l) 分布的随机数的方法. 在a,b ³ 1时，设 1 2 h ,h 分别为独立的G(a,1),G(b ,1) 随机数, 则可以证明 1 2 1 h h h V + = 为 B(a, b ) 随机数. 关于0 <a <1时G (a,l) 分布的取样方法，可以参见文献 [LK]. 2. 多维随机数 2. 1 连续型多维随机数 对于已知的分布密度, 我们可以利用条件密度, 把生成多维随机数归结为生成一系列 一维随机数. 设随机向量( , , ) X1 L Xd 的密度为 ( , , , ) 1 2 d f x x L x . 那么，我们有表达式 ( , , , ) ( ) ( | ) ( | , , ) 1 2 d = X1 1 2 1 d 1 d -1 f x x L x f x f x x L f x x L x , 其中 X1 f 为 X1的边缘密度, ( | , , ) k 1 k-1 f x x L x 为在已知 1 1 1 1 , , = k - = k- X x L X x 条件下, Xk 的条件密度. 于是可以先取一个 X1 f 随机数 1 x ; 然后，在 1 x 固定的情形下, 生成一个 ( | )1 f × x 随机数 2 x ; 再在 1 2 x , x 固定的情形下, 生成一个 ( | , ) 1 2 f × x x 随机数 3 x ; … 最后, 在 1 1 , , d - x L x 固定的情形下 , 生成一个 ( | , , ) 1 -1 × d f x L x 随机数 d x . 这样得到的 ( , , ) 1 d x L x 就是随机向量 ( , , ) X1 L Xd 的一个随机数. 当然在生成各个条件密度的随机数 时, 仍然可以使用 Von Neuman 取舍原则. 注意在用取舍原则于一维分布取样时，可以忽 略一个因子． 然而，在 d 较大时, 更为常用的是 Gibbs 取样法 (Gibbs Sampler), 它是一种基本的动态 Monte Carlo 方法, 即 Markov 链 Monte Carlo 方法, 在本书第 8 章中我们将介绍这种方法. 在很多实际情形中, 多维密度常常并不知道, 更无法知道各个条件密度. 这时最自然而 粗糙的想法是用条件频率来近似条件密度. 2. 2 离散型多维随机数 只要用概率函数 ( , , ) ( , , ) 1 d 1 1 d d p x x =P X = x X = x D L L 代替密度函数, 2. 1段中的办法 就自然适用. 2. 3 多维正态随机数
-(x-u)∑-(x-p) 设(X1,…,X4)服从N(μΣ),其密度函数为p(x)= (2x)2|∑|2 其中∑=G 为对称正定矩阵.由线性代数我们知道,可以找到下三角形矩阵C cn=0.V1<j),使Σ=CC作独立N(01)随机数51;…,54令5=(51…5),则 C5+~N(4,Σ).可见C5+μ为N(μ,Σ)随机数.于是求多维正态随机数的快速生成 问题,就归结为快速计算下三角形矩阵C的问题 2.4多维Beta随机数 (Dirichlet随机数)的生成 在a1…,ak≥1时,设n1,…,nk分别为独立的r(α1,1),…,r(αk1)-随机数,则可以 证明 nu :+ 就是一个 Dirichlet(ax1,…axk)随机数 附录一用 Matalab生成随机数 Matlab语言的简单提示 Maab提供了强大功能的软件包,其用法可以参阅 Matlab6.l(或 Matlab6.5)的Hep为 了使用方便,我们在本段列出一些基本的记号与最常见的语句 Matlab中的矩阵 矩阵A的分量用方括号括起来.同行的分量用逗号或空格分开,不同行的分量用分 号或回车( Enter)分开.空的方括号[]表示没有分量的矩阵,即空集 2.A(36)表示矩阵A的(3,6)分量,反之给出了它的值也可以调用A(在A的行列不 足3×4时将自动增加行列,即除了在该分量处赋要求的值以外,将其它分量增补为0) 3.矩阵运算除+,一外,还有:*(乘),’(转置)^(幂 (左除,A\B即AB),/(右除,A/B即BA) 常数视其情形可以看成任意阶的矩阵,如A+a(A的所有分量都加常数a) 4.可以用冒号对矩阵进行剪裁如A[3,]表示取A的第3行,A[,3]表示取A的第 列,A[1:3:]表示取A的1,2,3行,A1:23,]表示取A的1,3两行.又用A[BJ=[]表示 删除A的第3行 5.两个矩阵可以拼接成一个 6.特殊的矩阵(在相应的语句后加(m,n)表示其阶)有: zeros(0矩阵),ones(分量全 1的矩阵,eyes(对角为1,其它为0的矩阵),rand(0,1均匀随机数的矩阵)等 二 Matlab中的数组 1.与矩阵的差别是数组不加括号.用冒号表示等差数组,即m:d:n表示公差为d
39 设( , , ) X1 L Xd 服从 N(m,S), 其密度函数为 ( ) ( ) 2 1 2 1 2 1 (2 ) | | 1 ( ) m S m p S - - - - = x x d T p x e ， 其中 ( ) i j d ij £ = , S s 为对称正定矩阵. 由线性代数我们知道, 可以找到下三角形矩阵C , (c 0, i j) ij = " < , 使 T S = CC . 作独立 N(0,1) 随机数 d x , ,x 1 L . 令 T d ( , , ) 1 x x L x D = , 则 Cx + m ~ N(m,S) . 可见Cx + m 为 N(m,S) 随机数. 于是求多维正态随机数的快速生成 问题, 就归结为快速计算下三角形矩阵C 的问题. 2. 4 多维 Beta 随机数 (Dirichlet 随机数) 的生成 在a1 ,L,ak ³ 1时, 设h hk , , 1 L 分别为独立的G(a1 ,1),L,G(ak ,1) - 随机数, 则可以 证明 ( , , 1 1 L L k h h h + + ) 1 k k h h h +L+ 就是一个 Dirichlet ( , , ) a1 L ak 随机数. 3．附录 － 用 Matalab 生成随机数 3．1 Matlab 语言的简单提示 Matlab 提供了强大功能的软件包，其用法可以参阅 Matlab 6.1(或 Matlab 6.5) 的 Help．为 了使用方便，我们在本段列出一些基本的记号与最常见的语句． 一．Matlab 中的矩阵: 1．矩阵 A 的分量用方括号括起来．同行的分量用逗号或空格分开, 不同行的分量用分 号或回车(Enter)分开. 空的方括号[ ]表示没有分量的矩阵，即空集． 2．A(3,6)表示矩阵 A 的(3,6)分量．反之, 给出了它的值也可以调用 A (在 A 的行列不 足3´ 4时将自动增加行列, 即除了在该分量处赋要求的值以外，将其它分量增补为 0)． ３．矩阵运算除＋, －外, 还有：*（乘）, ' (转置), ^ (幂), \ (左除, A\ B即 A B -1 ), / (右除， A/ B 即 -1 BA )． 常数视其情形可以看成任意阶的矩阵，如 A+a (A 的所有分量都加常数 a）． ４．可以用冒号对矩阵进行剪裁, 如 A[3,:]表示取 A 的第 3 行, A[:,3]表示取 A 的第 3 列，A[1:3,:] 表示取 A 的 1,2,3 行, A[1:2:3,:] 表示取 A 的 1,3 两行．又用 A[3,:]=[ ]表示 删除 A 的第 3 行． ５．两个矩阵可以拼接成一个． ６．特殊的矩阵（在相应的语句后加(m, n) 表示其阶）有：zeros(0 矩阵），ones（分量全 1 的矩阵），eyes（对角为 1，其它为 0 的矩阵），rand ([0，1]均匀随机数的矩阵)等． 二 Matlab 中的数组 １. 与矩阵的差别是数组不加括号．用冒号表示等差数组，即 m : d : n 表示公差为 d
起始于m,终止于n的等差数组,而将m:1:n简记为m 2.数组间在运算符号前加一个小点表示按分量运算,有: *(乘),^(幂),/除(即右除)),\(反除(即左除),即除的倒数 3.几个数组可以用方括号拼接成矩阵 单引号括起来的字符串可以作为数组 又如语句 linspace(a,b,m)生成起始于a,终止于b的n个等差数组 逻辑关系 1.一个关系的逻辑值只取两个值:0(表示"错误”),或1(表示"正确”).所以 关系ⅹ==0的逻辑值,就是ⅹ=0的示性函数 2.主要的逻辑关系与逻辑运算有: (逻辑相等),>=(不小于),<=(不大于),~=(不等于),&(与),|(或), 四.语句,函数与操作命令 1.一般的语句的形式是 变量=表达式 (如用分号结束,则表示只运算不显示,如用回车结束则显示它的计算结果.如只有表达式, 则自动生成变量ans,并显示ans=) 2.两个语句间用逗号或分号 3.特殊变量有pi(圆周率丌),i或j(虚数单位),Inf(无穷大), NaN(不定值:0/0),eps(主机中的最小的浮点数,有时用它代替0可以防止溢出,即 用表达式(x=0)eps代替0,它表示在x=0时用eps代替x) 4.函数的变量置于圆括号之中.常见的函数除对数,指数,三角函数外,还有: atan(反正切),abs(绝对值),sqt(开方), round4舍5入取整,foor(负无穷方向取 整),cei(正无穷方向取整,fix(0方向取整,mag(虚部),rats(最近有理数) length(数组的长度),prod(乘),sor(数组的分量按递增排序), trace(求迹),tru(矩阵取 上三角),til(矩阵取下三角)size(矩阵的大小),norm(模),inv(取逆),eg(求特征 值),poly(求特征多项式),expm(求矩阵的指数函数),cond(条件的数目),lu(矩阵的LU分解) gr(矩阵的正交分解),svd矩阵的奇异值分解), feval(赋值语句,如 feval('atan',x)等价于atan(x) 5.用于循环的常见语句有 ause(暂停), all while end if elseif else. end switch. case. otherwise. end break(跳出循环) 6.操作命令(很多与Unⅸ命令相同) help, lookfor, format short(显示到小数点后4位),whos(显示当前工作区的变量) disp(x)显示x的内容 clear,ave,load, diary filename(建立新文件), diary off(停止记录), path, dir, type, delete, quit
40 起始于m ，终止于n 的等差数组, 而将m :1: n 简记为m : n ２. 数组间在运算符号前加一个小点表示按分量运算, 有：. *（乘）, .^（幂）, ./ (除（即右除)）, .\ (反除（即左除），即除的倒数)． ３. 几个数组可以用方括号拼接成矩阵． ４. 单引号括起来的字符串可以作为数组． 又如语句linspace(a, b, n) 生成起始于a ，终止于b 的n 个等差数组． 三．逻辑关系 １．一个关系的逻辑值只取两个值：0（表示＂错误＂），或 1（表示＂正确＂)．所以， 关系 x==0 的逻辑值，就是 x=0 的示性函数． ２．主要的逻辑关系与逻辑运算有： ＝＝（逻辑相等）），> =（不小于），< =（不大于）, ~ =（不等于）, ＆（与），｜（或）， ～（非）． 四．语句，函数与操作命令 １．一般的语句的形式是 变量=表达式 (如用分号结束，则表示只运算不显示．如用回车结束则显示它的计算结果．如只有表达式， 则自动生成变量 ans，并显示 ans =)． ２．两个语句间用逗号或分号． ３．特殊变量有 pi (圆周率p ), i 或 j (虚数单位），Inf （无穷大）， NaN （不定值：0／0）, eps (主机中的最小的浮点数, 有时用它代替０可以防止溢出，即 用表达式(x==0)*eps 代替 0, 它表示在 x=0 时用 eps 代替 x)． ４．函数的变量置于圆括号之中．常见的函数除对数，指数，三角函数外，还有： atan（反正切）, abs（绝对值）, sqrt（开方）, round(4 舍 5 入取整), floor(负无穷方向取 整), ceil(正无穷方向取整), fix(0 方向取整), imag（虚部）, rats(最近有理数) length（数组的长度）, prod（乘）, sort(数组的分量按递增排序)，trace(求迹)，triu (矩阵取 上三角), tril (矩阵取下三角), size（矩阵的大小）, norm（模）, inv（取逆）, eig（求特征 值）, poly(求特征多项式), expm(求矩阵的指数函数), cond(条件的数目), lu(矩阵的 LU 分解), gr(矩阵的正交分解), svd(矩阵的奇异值分解)， feval（赋值语句，如 feval('atan',x) 等价于 atan(x)． ５．用于循环的常见语句有： pause（暂停）， all, any, for, end while, end, if, elseif, else, end switch, case, otherwise, end break(跳出循环) ． ６．操作命令（很多与 Unix 命令相同）： help, lookfor, format short（显示到小数点后４位）, whos（显示当前工作区的变量）, disp(x)(显示 x 的内容), clear, save, load, diary filename（建立新文件）, diary off（停止记录）, path, dir, type, delete, quit 等．
7.画图命令有 ltxy)(一个图), plot(x,y1,xy2)(两个图),或 fplot(图形名 xmin xmax ymin ymax) 或在图上加: hold on,第2个图的语句, hold off 线图: :(点线),-,(虚点线),一-(波折) 点图: 颜色:y(黄),r(红),g(绿),b(兰),w(白),k(黑),m(紫),c(青cyan) 例如:plot(xy1,b:xy2,c") 控制坐标范围用: axis auto, axis off 在图上标字:用ttle字串),或用 gtext("字串)再用鼠标点 分区画多个图:用 subplot(m,nk)(在m×n个区中的第k个区中画图) 曲面:[X,Y]= meshgrid(xy):(把(x,y)改写为数组[X,Y],以便用于表达式中) Z=(X,Y)的表达式;mesh(X,Y,Z) 空间曲线图例如:plot3(sin(t),cos(t)t) 等高线图: contour(X,Y,Zm)(m条等高线) 3.2 Matlab生成随机数的语句 第一种方法是用 random语句,其一般形式为 y= random(分布的英文名,A1,A2,A3m,n), 表示生成m行n列的m×n个参数为(A1,A2,A3)的该分布的随机数。例如 (1)R= random(Noma",0,224):生成期望为0,标准差为1的2行4列)2×4个正态随机数 (2)R= random( Poisson,1:6,1,6)依次生成参数为1到6的(1行6列)6个 Poisson随机数 第二种方法是针对特殊的分布的语句: 几何分布随机数(下面的P,m都可以是矩阵 R= ground(P)(生成参数为P的几何随机数) R= geornd(P,m)(生成参数为P的1×m个几何随机数) R= geornd(Pm,n)(生成参数为P的m行n列的m×n个几何随机数) 例如 (1)R= geornd(l1/2(1:6)(生成参数依次为12,1/22到1/2^6的6个几何随机数) (2)R= geornd(0.01,15])(生成参数为001的(1行5列)5个几何随机数) 二.Beta分布随机数 R= betarnd(A,B(生成参数为AB的Beta随机数) R= betarnd(A,Bm)(生成1×m个数为AB的Bea随机数) R= betarnd(A,Bm,n)(生成m行n列的m×n个数为AB的Beta随机数) 三.正态随机数 R= normand(MU, SIGMA)(生成均值为MU,标准差为 SIGMA的正态随机数) R= normand(MU, SIGMA,m)(生成1×m个正态随机数) R= normand(MU,SGMA,m,n)(生成m行n列的m×n个正态随机数) 例如 (1)R= normand(O,1,15生成5个正态(O,1)随机数 (2)R= normand(123456]0.1,2,3)生成期望依次为1234,56],方差为0.1的6个正态随
41 ７．画图命令有： plot(x,y)（一个图）, plot(x,y1,x,y2)（两个图）, 或 fplot('图形名',[xmin xmax ymin ymax]) 或在图上加：hold on, 第２个图的语句, hold off 线图： -， ：(点线), -．(虚点线)， －－(波折) 点图： ．， +，* 颜色： y（黄），r（红），g（绿），b（兰），w（白），k(黑)，m(紫)，c(青 cyan) 例如：plot(x,y1,'b:',x,y2,'c-') 控制坐标范围用：axis auto, axis off 在图上标字：用 title('字串'), 或用 gtext('字串')再用鼠标点 分区画多个图：用 subplot(m,n,k) (在m´ n 个区中的第 k 个区中画图） 曲面：［X,Y］＝meshgrid(x,y)；（把(x，y)改写为数组［X,Y］，以便用于表达式中）； Z＝(X,Y)的表达式; mesh(X,Y,Z) 空间曲线图例如：plot3(sin(t),cos(t),t) 等高线图：contour(X,Y,Z,m) (m 条等高线)。 ３．２ Matlab 生成随机数的语句 第一种方法是用 random 语句，其一般形式为 y = random('分布的英文名',A1,A2,A3,m,n)， 表示生成m 行n 列的m´ n 个参数为( , , ) A1 A2 A3 的该分布的随机数。例如: (1) R = random('Normal',0,２,2,4): 生成期望为 0,标准差为 1 的(2 行 4 列)2´ 4 个正态随机数 (2) R = random('Poisson',1:6,1,6): 依次生成参数为 1 到 6 的(1 行 6 列)6 个 Poisson 随机数 第二种方法是针对特殊的分布的语句： 一． 几何分布随机数 （下面的 P，m 都可以是矩阵） R = geornd(P) （生成参数为 P 的几何随机数） R = geornd(P,m) （生成参数为 P 的１´ m 个几何随机数） R = geornd(P,m,n) （生成参数为 P 的m 行n 列的m´ n 个几何随机数） 例如 (1) R = geornd(1 ./ 2.^(1:6)) (生成参数依次为 1/2,1/2^2,到 1/2^6 的 6 个几何随机数) (2) R = geornd(0.01,[1 5]) (生成参数为 0.01 的（１行５列）5 个几何随机数). 二．Beta 分布随机数 R = betarnd(A,B) （生成参数为 A,B的 Beta 随机数） R = betarnd(A,B,m) （生成１´ m 个数为 A,B的 Beta 随机数） R = betarnd(A,B,m,n) （生成m 行n 列的m´ n 个数为 A,B的 Beta 随机数）. 三．正态随机数 R = normrnd(MU，SIGMA) （生成均值为 MU，标准差为 SIGMA 的正态随机数） R = normrnd(MU，SIGMA,m) （生成１´ m 个正态随机数） R = normrnd(MU，SIGMA,m,n) （生成m 行n 列的m´ n 个正态随机数） 例如 (1) R = normrnd(0,1,[1 5]) 生成 5 个正态(0,1)随机数 (2) R = normrnd([1 2 3;4 5 6],0.1,2,3) 生成期望依次为[1,2,3;4,5,6],方差为 0.1的 6 个正态随
机数 四.二项随机数:类似地有 R= binornd(N, P) R=binornd(N, P,m) R-binornd(N, p, m, n) 例如 n=10:1060,rl= binornd(n,1/m)或r2= binornd(nh[16)(都生成参数分别为 10,),…,(60,一)的6个二项随机数 五.自由度为V的x2随机数 R=chi2rnd(v) R= chi2rnd( v, m) R=chi2rnd(V,m, n) 六.期望为MU的指数随机数(即Exp1随机数) R= expand(MU) R= expand(MU,m) R= expr rnd(MU,m, n 七.自由度为V1,V2的F分布随机数: R=find(vl, v2) R=find(vi, v2, m) R= frnd(vl, V2, m, n) 八.I(A,λ)随机数: R= gamrnd(A, lambda) R= gamrnd(A, lambda, m)R= gamrnd (A,lambda, m, n) 九.超几何分布随机数: R= hygernd(N, K, M) R= hygernd(N, K, M, m) R=hygernd(N, K, M,m, n) 十.对数正态分布随机数 R= lognrnd(MU, SIGMA)R= lognrnd(MU, SIGMA, m)R=lognrnd(MU, SIGMA,m, n) 十一.负二项随机数: R=nbinrnd(r, p) R=nbinrnd(r, p,m R=nbinrnd(r, p, m, n) 十二. Poisson随机数: R= poissrnd(lambda) R=poisson srnd( lambda, m, n) 例如,以下3种表达有相同的含义: lambda=2,R= poissrmd(lambda,1,10) (Er=poissrnd(lambda, [1 10) Er= poissrnd(lambda( ones(1, 10))) 十三. Rayleigh随机数: R=raylrmd(B) raylrnd(B R= raylrnd(B, m, n) 十四.V个自由度的t分布的随机数: R=trnd(v) R=trnd(v,m) R=trnd(v, m, n)
42 机数． 四．二项随机数：类似地有 R = binornd(N,P) R = binornd(N,P,m) R = binornd(N,p,m,n) 例如 n = 10:10:60; r1 = binornd(n,1./n) 或 r2 = binornd(n,1./n,[1 6]) （都生成参数分别为 ) 60 1 ), ,(60, 10 1 (10, L 的６个二项随机数． 五．自由度为 V 的 2 c 随机数： R = chi2rnd(V) R = chi2rnd(V,m) R = chi2rnd(V,m,n) 六．期望为 MU 的指数随机数（即 MU Exp 1 随机数）： R = exprnd(MU) R = exprnd(MU,m) R = exprnd(MU,m,n) 七．自由度为 V1，V2 的 F 分布随机数： R = frnd(V1，V2) R = frnd(V1，V2,m) R = frnd(V1，V2,m,n) 八．G(A,l) 随机数： R = gamrnd（A,lambda） R = gamrnd（A,lambda,m) R = gamrnd（A,lambda,m,n) 九．超几何分布随机数： R = hygernd(N,K,M) R = hygernd(N,K,M,m) R = hygernd(N,K,M,m,n) 十．对数正态分布随机数 R = lognrnd(MU，SIGMA) R = lognrnd(MU，SIGMA,m) R = lognrnd(MU，SIGMA,m,n) 十一．负二项随机数： R = nbinrnd(r,p) R = nbinrnd(r,p,m) R = nbinrnd(r,p,m,n) 十二．Poisson 随机数： R = poissrnd(lambda) R = poissrnd(lambda,m) R = poissrnd(lambda,m,n) 例如，以下 3 种表达有相同的含义：lambda = 2; R = poissrnd(lambda,1,10) （或 R = poissrnd(lambda,[1 10]) 或 R = poissrnd(lambda(ones(1,10))) 十三．Rayleigh 随机数： R = raylrnd(B) R = raylrnd(B,m) R = raylrnd(B,m,n) 十四．V 个自由度的 t 分布的随机数： R = trnd(V) R = trnd(V,m) R = trnd(V,m,n)
