1、蒙特卡罗方法及应用,第三章由已知分布的随机抽样,3.1 随机抽样及其特点,随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。,由己知分布的随机抽样指的是由己知分布的总体中抽取简单子样。,本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用严格的数学方法，从已知分布中随机抽样产生。,只要所采用的随机数满足独立性和均匀性要求，那么由上述方法产生的简单子样严格服从具有相同的总体分布并相互独立。,对于任意给定的累积分布函数F(x)，要从中抽取随机变量序列X1,X2,Xn ，则可以产生一0,1间均匀分布的随机数序列， 1, 2, n ，按
2、下式来确定Xi值,也就是说:当抽取j后，判断所要取得各点值中，所对应的分布函数值大于j，则其最小值/下界为所取值。,可以证明随机变量序列X1,X2,XN具有相同分布。,3.2.1 离散型分布的直接抽样方法,直接抽样法适应于离散型的随机变量,设离散型随机变量X的可能取值为x1, x2, , xN, 其概率为,累积分布函数：,0,x1,xN-1,xN,x2,xk-1,xk,1,F(x),pk,0,x1,xN-1,xN,p1,p2,pN,x2,pk,xk-1,xk,3.2.1 离散型分布的直接抽样方法,方法：,计算yk = yk-1 + pk，k = 2,3,N, y1 = p1 产生在0,1区间上
3、均匀分布的随机数 ； 求满足yk-1yk 的k值; 随机变量的第k个取值即为欲抽取的值。,0,x1,xN-1,xN,x2,xk-1,xk,1,F(x)=y,pk,0,x1,xN-1,xN,p1,p2,pN,x2,pk,xk-1,xk,3.2.1 离散型分布的直接抽样方法,证明：,0,x1,xN-1,xN,x2,xk-1,xk,1,F(x),pk,0,x1,xN-1,xN,p1,p2,pN,x2,pk,xk-1,xk,即:所产生的随机数的概率密度函数/pdf为pk,3.2.1 几个典型的例子,例1、某核素衰变末态的随机抽样,设某核素X有三种衰变方式，其分支比如下,随机选取每次衰变的衰变方式（衰变
4、道）直接抽样法, = U0,1,3.2.1 几个典型的例子,例2、中子与核的反应类型确定,假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，辐射俘获，裂变，相应的反应截面分别为el ,in,f ,a。则发生每一种反应类型的概率依次为 ：,其中反应总截面t el in f a。,反应类型的确定方法为：产生一个随机数,3.2.1 几个典型的例子,例3、碰撞目标核种类的确定,中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定到底与那种元素发生碰撞。假定入射粒子与介质中每种核的发生反应的宏观总截面分别为1, 2, , n，则中子或光子与每种核碰撞的概率分别为：,其中t12n。碰撞核种
5、类的确定方法为：产生一个随机数 ，如果,则中子或光子与第I种核发生碰撞。,3.2.1 几个典型的例子,例4、掷骰子点数的抽样,掷骰子点数X=n的概率为：,选取随机数 ，如果,则所掷骰子的点数,在等概率的情况下，可使用如下更简单的方法：,其中表示取整数。,3.2.1 几个典型的例子,例5、二项式分布的抽样,方法1：利用上面介绍的直接抽样法，需计算累积分布函数，当n很大时，求和计算困难;,方法2：利用二项式分布的定义,产生n个 iU0,1; 统计满足条件 i p（表示成功的概率）的 i的数目r，则r表示在n次实验中成功的次数,r即为二项式分布的抽样值,3.2.1 几个典型的例子,例6、泊松分布的抽
6、样,方法1：利用直接抽样法，但计算累积分布函数时非常复杂,方法2：利用泊松分布的定义：二项式分布的极限形式,选取足够大的n，使p=/n相当小，例如，p=0.1 产生n个 iU0,1; 统计满足条件 i p（表示成功）的 i的数目k，则k表示在n次实验中成功的次数。,k即为泊松分布的抽样值的近似值， n越大，近似程度越好,3.2.2 连续型分布的直接抽样方法,对于连续型分布，如果分布函数y=F(x) 的反函数x=F1(y)存在，则直接抽样方法是 ：,方法：,产生在0,1区间上均匀分布的随机数 = U0,1 ；,令F(x) = , 解方程得x：,3.2.2 几个典型的例子,例7、在a,b上均匀分布
7、的抽样,则,3.2.2 几个典型的例子,例8、 分布抽样,分布为连续型分布，作为它的一个特例是：,其分布函数为：,则,3.2.2 几个典型的例子,例9、 指数分布抽样（中子在介质中的输运长度服从负指数分布）,指数分布为连续型分布，其一般形式如下：,其分布函数为：,则,因为1也是随机数，可将上式简化为,3.3 挑选抽样方法,连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。,2、分布函数可以给出其解析形式，但是反函数给不出来。,1、分布函数无法用解析形式给出，因而其反函数也无法给出。,3、分布函数即使能够给出反函
8、数，但运算量很大。下面叙述的挑选抽样方法是克服这些困难的比较好的方法。,3.3 挑选抽样方法,由 von Neumann 1951年提出并完成。,抽取随机变量x的一个随机序列xi, i=1,2, 按一定的舍选规则从中选出一个子序列，使其满足给定的概率分布.,3.3 挑选抽样方法,挑选法抽样步骤：,产生a, b区间内均匀分布的随机数x: x = (b-a)r1+a, r1 U0, 1; 产生0,c区间内均匀分布的随机数y: y = cr2, r2 U0,1; 当y f(x)时，接受x为所需的随机数，否则，返回到第一步重新抽取一对(x,y).,设随机变量x的取值区间为xa,b, 其间概率密度函数f
9、(x)有界，即,3.3 挑选抽样方法,a,b,x,f(x),c,几何解释:,在二维图上，随机选取位于矩形abfe内的点(x,y)； 选取位于曲线f(x)下的那些点，则这些点将服从概率密度为f(x)的分布,e,f,3.3 挑选抽样方法,证明:,按挑选抽样法抽出的随机数d的概率：,a,b,x,f(x),c,e,f,x和y的概率密度函数分别为,联合概率密度函数为,即d的概率函数为f(x),d,3.3 挑选抽样方法,抽样效率:,对挑选抽样法：欲产生m个随机变量x的值需产生n对(x,y)，显然，m n,如果选出某特定分布的一个随机数平均地需要n个随机数r1 U0, 1，则抽样效率定义为,z,3.3 挑选
10、抽样方法,改进的挑选抽样法,3.3 挑选抽样方法,抽样方法：,1. 产生两个随机数,产生分布为h(x) 的随机数x，xa,b; 产生0, Mh(x) 区间上均匀分布的随机数y，y= Mh(x) , U0,1.,2. 接收或舍弃取样值 x.,如果 y f(x)，舍弃，返回到1,重复上述过程; 否则，接受x;,3.3 挑选抽样方法,几何解释:,在二维图上，随机选取位于曲线Mh(x)下的点x,y； 选取位于曲线f(x)下的那些点，则这些点将服从概率密度为f(x)的分布,x,3.3 挑选抽样方法,证明:,按挑选抽样法抽出的随机数d的概率：,d,x和y的概率密度函数分别为,联合概率密度函数为,即d的概率
11、函数为f(x),x,3.3 挑选抽样方法,抽样效率:,x,常数M的选取,常数M应尽可能地小,因为抽样效率与M成反比; M=maxf(x)/h(x), x a,b,3.3 几个典型的例子,例10:标准正态分布的抽样,x-a,a,无法用直接抽样法,累积分布函数无解析表达式,Breit-wigner or Cauchy分布,先对h(x)归一化，以保证其累积分布最大值为1,3.3 几个典型的例子,由h(x)抽取x 直接抽样法,抽取u,计算f(x), 如果u= f(x), 接受x,例10:标准正态分布的抽样,x-a,a,3.3 几个典型的例子,例11：利用挑选法产生随机数C=cos, S=sin,其中为
12、0, 2区间内均匀分布的随机数,方法1：先产生0, 2间均匀分布的随机数： = 2 r， rU0,1, 然后直接计算C和S 因需要计算三角函数，故此方法运算速度慢,方法2：利用挑选法可避免三角函数运算,令A和B为单位圆内直角三角形的两个边，则有,3.3 几个典型的例子,因此，只要产生单位圆内的随机坐标A和B,就可用代数运算求出C和S,算法为,产生两个0,1区间上均匀分布的随机数1和2; 令v1=21-1, v2=2,则v1U-1,1,v2U0,1; 计算r2=v12+v22, 如果r2 1,转到1，重新产生； 令A=v1, B=v2, 计算C和S,例11：利用挑选法产生随机数C=cos, S=
13、sin,其中为0, 2区间内均匀分布的随机数,作业：从上面可知，该方法的抽样效率为/4，试列举一种抽样效率更高的方案，并给出详细的抽样原理及步骤。,3.4 复合抽样方法,1961年由Marsaglia提出的方法,设随机变量X的概率密度函数f(x)可写成一些PDF的线性叠加：,抽样方法：,利用离散型的随机变量的抽样方法抽取序号k;,由fk(x)抽取x,3.4 几个典型的例子,例12：用复合法产生双指数分布的随机数,产生两个0,1区间均匀分布的随机数1和2; 如果1 0.5, 按f1(x)抽样； 如果1 0.5, 按f2(x)抽样；,3.4 几个典型的例子,例13：指数函数分布的随机数,产生n+1
14、个0,1区间均匀分布的随机数 1, 2, n+1; 按f1(y)抽取y； 按f2(x/Yf1)抽取x；,3.5 替换抽样方法,为了实现某个复杂的随机变量 y 的抽样，将其表示成若干个简单的随机变量 x1，x2，xn 的函数,得到 x1，x2，xn 的抽样后，即可确定 y 的抽样，这种方法叫作替换法抽样。即,3.5 几个典型的例子,例14：高斯分布的抽样方法,进行变量变换：,由N(0,1)分布抽样得到y,如何抽取服从N(0,1)分布的随机变量？,3.5 几个典型的例子,方法1、利用中心极限定理,设x1,x2,xn是一组n个独立的随机变量，xi的平均值和方差分别为i和i，则当n时，变量,服从标准正
15、态分布N(0,1),设x 是在0, 1之间均匀分布的随机数, 对n个x的取值xi,在n时，y服从N(0,1)，在实际应用时，可取n=12,抽样方法：,1、产生12个U0,1的随机数i,2、,3.5 几个典型的例子,方法2、利用替换抽样法,y1和y2是两个相互独立的、服从标准正态分布的随机数,变换：,反变换：,抽样方法：,1）产生一对0,1区间均匀分布的随机数1和2； 2）,3.5 几个典型的例子,证明用此方法抽取的y1,y2满足上面的联合概率分布,雅可比行列式：,3.5 几个典型的例子,方法3、采用极坐标变换,x和y是两个相互独立的、服从标准正态分布的随机数,变换：,反变换：,则变量和的联合概
16、率密度函数为,即在0,2 区间上均匀分布, 服从 = 1的指数分布,3.5 几个典型的例子,抽样方法：,产生 ：= 2 r，rU0,1 产生 ： = -ln r， rU0,1x= (2)1/2 cos ， y=(2)1/2 sin,3.6 随机抽样的一般方法,3.6.1、加抽样方法 =复合抽样方法 3.6.2、减抽样方法 3.6.3、乘抽样方法 3.6.4、乘加抽样方法 3.6.5、乘减抽样方法 =减抽样方法 3.6.6、对称抽样方法 3.6.7、多维分布抽样方法 3.6.8、积分抽样方法,3.6 随机抽样的一般方法,乘抽样方法,由复合抽样的概率密度函数可以看出,那么接受ff(x)抽样到的x值
17、的概率为,用一函数g(x)代替分布函数 ，得乘抽样概率密度函数,3.6 随机抽样的一般方法,乘抽样方法步骤：,1. 将ff(x)归一化：,2. 令,Mg为g(x)在区间a,b上的极大值即上界,3. 利用直接抽样法由 抽取x;,4 . 抽取0,1区间均匀分布的随机数r, 如果 ，则接受x, 否则，返回到3重新抽样,3.6 随机抽样的一般方法,减抽样方法,变换得到,可以看出上式可以看做乘抽样方法，其中,若取,3.6 随机抽样的一般方法,减抽样方法,则,可知各自的接受概率分布函数为（归一化函数g(x)）,抽样步骤：,由f1(x)抽取x; 2. 抽取U0,1; 3. 若 ，那么 接受 x，否则返回1重
18、复,由f2(x)抽取x; 2. 抽取U0,1; 3. 若 ，那么 接受 x，否则返回1重复,3.6 随机抽样的一般方法,设概率密度函数可写成如下的形式：,乘抽样方法的推广形式,复合抽样法+混合抽样法=乘加抽样法,抽样方法：,1. 采用复合抽样法，先确定p(x)的随机数应由哪一个 抽取； 2. 在按 抽样时，采用上面的乘抽样方法确定从fi(x)中的抽样是否被接受。,3.6 几个典型的例子,例15：compton散射,微分截面：,如何抽取散射光子的能量？=乘加抽样法,3.6 几个典型的例子,3.6 几个典型的例子,抽样方法：r1,r2,r3是三个在0,1区间上均匀分布的随机数,1. 确定由哪一个f
19、函数来抽取：如果r11/(1+2)，选择f1(), 否则选f2()； 2. 根据f1或f2抽取：直接抽样法,3. 计算sin2:,4.计算g(), 如果 接受， 否则返回到第一步,3.6 随机抽样的一般方法,乘减抽样方法,乘减分布的形式为：,其中H1(x) 、H2(x)为非负函数，f1(x)、f2(x) 为任意分布密度函数。,与减抽样方法类似，乘减分布的抽样方法也分为两种。,3.6 随机抽样的一般方法,对称抽样方法,课本上存在原理错误：必须在保证 成立的情况下，原理上才能解释的通。,从乘抽样方法了解到：,1. 先采用直接抽样方法或挑选抽样方法从f1(x)中抽取 ；如采用直接抽样方法： 2. 那
20、么可选取随机数2 ，那么接受 的概率为,3.6 随机抽样的一般方法,对称抽样方法,1. 先采用直接抽样方法或挑选抽样方法从f1(x)中抽取得到 如采用直接抽样方法： 2. 那么可选取随机数2 ，那么接受 的概率为,为了提高抽样效率，已知U0,1为随机数，那么1- U0,1那么做一次实验相当于做两次，用随机数失败后，可以考虑随机数1- 。,3.6 随机抽样的一般方法,对称抽样方法,抽样过程为：,3.6 随机抽样的一般方法,多维分布抽样方法,根据条件概率密度公式，两变量联合概率密度分布函数可表示成,式中,1、由f1(x)抽样得到x的抽样值X； 2、将X代入f2(y|x)中得到y的一维分布f2(y|
21、X)；由f2(y|X)抽样得到y的抽样值Y。,抽样方法：,3.6 随机抽样的一般方法,多维分布抽样方法,利用条件概率密度将多维模拟转化为一维模拟问题。,任意n维概率密度函数fn(x1,x2,xn)可表示为,fn(x1,x2,xn) = f(xn| x1,x2,xn-1) fn-1(x1,x2,xn-1),条件概率,x1,x2,xn-1联合概率密度,fn(x1,x2,xn) = f(xn| x1,x2,xn-1) f(xn-1| x1,x2,xn-2)f(x2| x1)f(x1),n个一维概率密度的乘积,3.6 随机抽样的一般方法,多维分布抽样方法,fn(x1,x2,xn) = f(xn| x1
22、,x2,xn-1) f(xn-1| x1,x2,xn-2)f(x2| x1)f(x1),抽样方法：,1、由f(x1)抽样得到x1的抽样值X1 ； 2、将X1代入f(x2| x1)中得到x2的一维分布f(x2|X1)；由f(x2|X1)抽样得到x2的抽样值X2 ； ,3.6 随机抽样的一般方法,积分抽样方法,称为积分分布密度函数，其中 f0(x，y) 为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样步骤为：,1、由f0(x，y)中采用二维函数抽样方法得到X,Y； 2、如果Y H(X)，那么接受X；否则返回第一步继续执行。,3.7 随机抽样的其它方法,3.7.1、偏倚抽样方法 3.
23、7.2、近似抽样方法 3.7.3、近似-修正抽样方法 3.7.4、指数分布抽样方法,3.7.2 近似抽样法,使用的场合：,概率密度函数的形式非常复杂，在模拟过程中进行计算需花费相当多的CPU时间; 概率密度函数无解析形式，只能用数值或曲线的形式表示。,3.7.2 近似抽样法,基本方法:,设概率密度函数：f(x), xa,b,将区间a,b分成n个子区间，分点为,分点对应的函数值为:,对于每一个小区间，利用一简单的函数fa(x)来近似地表示原概率分布函数，并使fa(x)在该区间内的积分与f(x)在该区间内的积分相等，即俩者的概率相等。,利用f(x)计算每一个子区间的概率值 并计算f(x)的累积分布函数在xi处的值 ,并将F(xi)和xi存入数据表中;,3.7.2 近似抽样法,抽样方法: 随机选择子区间:选取0,1区间内均匀分布的随机数r,找出满足 的子区间xi-1,xi; 根据fa(x)的特点,确定欲抽取的随机数；,几点说明：,分点的选取:,f0, f1, ., fn应能充分反映f(x)的变化状况,即,在 f(x) 变化迅速的区域分点密一点,变化缓慢的区域分点稀一点.,fa(x)的选取:,阶梯近似: 线性近似 二次曲线近似 样条函数近似,
