精选实验数据处理方法DataAnalysis-ch08课件.ppt

1、第八章第八章 从概率分布函数的抽样从概率分布函数的抽样(Sampling from Probability Distribution Functions)Monte Carlo算法的一个重要组成部分：算法的一个重要组成部分：描述所要模拟的物理系统的一些概率密度函数（描述所要模拟的物理系统的一些概率密度函数（PDF)描述整个系统在空间、能量、时间或多维相空间中的描述整个系统在空间、能量、时间或多维相空间中的发展和演化；发展和演化；Monte Carlo模拟的主要任务模拟的主要任务:通过对这些概率密度函数的随机抽样来模拟物理系统的状通过对这些概率密度函数的随机抽样来模拟物理系统的状态态;为描述系统

2、的演化所必需的一些附加运算为描述系统的演化所必需的一些附加运算.物理过程的描述物理过程的描述从描述物理系统的从描述物理系统的pdf出发，随机抽取系统的出发，随机抽取系统的可能状态。可能状态。本章介绍从一个任意的本章介绍从一个任意的pdf获取样本的抽样方法。获取样本的抽样方法。1.1.直接抽样法（反函数法）直接抽样法（反函数法）2.2.变换抽样法变换抽样法直接抽样法的一般形式直接抽样法的一般形式3.舍选抽样法舍选抽样法4.复合分布的抽样方法复合分布的抽样方法5.混合抽样法混合抽样法6.近似抽样法（列表法近似抽样法（列表法)7.多维分布的抽样多维分布的抽样8.1 等价的连续概率密度函数等价的连续概

3、率密度函数 随机变量：连续型、分离型随机变量：连续型、分离型 概率密度函数：连续分布、分离分布概率密度函数：连续分布、分离分布 利用利用 函数，可将分离型的函数，可将分离型的pdf用连续型的用连续型的pdf 描述描述 用同样的方式来讨论分离型和连续型随机变量的抽样方法用同样的方式来讨论分离型和连续型随机变量的抽样方法Niiixxpxf1)()(NiiiNiiiNiiipxdxxxxpdxxxpxdxxxfXE111)()()()(NiiiNiiiNiiipXExdxxxpXExdxxxpXExdxxfXExXV1212122)()()()()()()()(已知分离型已知分离型pdf:Pi 分离

4、型随机变量分离型随机变量X的取值为的取值为xi的概率的概率 定义一个等价的连续型定义一个等价的连续型pdf:)()()(1)(iiixfdxxxxfdxxx 利用与连续型随机变量相同的方式计算分离型随机变量的期望值利用与连续型随机变量相同的方式计算分离型随机变量的期望值和方差：和方差：8.2 pdf的变换的变换x:连续型的随机变量连续型的随机变量,PDF:f(x)y=y(x)：x的函数的函数,也是随机变量也是随机变量.求求y(x)的概率密度函数的概率密度函数g(x)1、若随机变量、若随机变量x和和y是一一对应的是一一对应的:2、若随机变量若随机变量x和和y不是一一对应的不是一一对应的:x,x+

5、dxy,y+dyX的取值在的取值在x,x+dx的概率的概率=Y的取值的取值在在y,y+dy的概率：的概率：dydxxfyg)()(取绝对值是为了保证取绝对值是为了保证g(y)是非负的是非负的f(x)dx=g(y)dy 即有即有n个区间个区间x,x+dxy,y+dy dydxxfyg)()(需要对这需要对这n个区间求和个区间求和3、推广到、推广到n个随机变量的情况个随机变量的情况:Jacobian行列式行列式nixyyyyyyxxxxiinn,2,1),(),(),(2121nnnnnyxyxyxyxyxyxnnyyyxxxJJxfyg2112111,),(),()()(21214、特例：如果、

6、特例：如果y(x)是是x的累积分布函数的累积分布函数(cdf)(1)(11)()()(xfxFdxdydydxxdxfxFxyyx10,1)(yyg即：即：y在在0,1区间上均匀分布区间上均匀分布 不管不管f(x)取何种形式取何种形式,累积分布累积分布函数总是在函数总是在0,1区间上均匀分布区间上均匀分布8.3 8.3 直接抽样法（反函数法）直接抽样法（反函数法）（Sampling via Inversion of the cdf）)(1yFx原理原理（注意：注意：pdf f(x)必须是归一化的必须是归一化的）：）：设设y=F(x)为随机变量为随机变量x的累积分布函数的累积分布函数 x和和y是

7、一一对应的是一一对应的 先随机抽取先随机抽取y，然后通过求，然后通过求F(x)的反函数的反函数F-1(y)得到随机变量得到随机变量x的值的值 随机变量随机变量y在区间在区间0,1上均匀分布上均匀分布 利用利用0,1区间上均匀分布区间上均匀分布随机数产生器抽取随机数产生器抽取方法：方法：U0,1:0,1区间上均匀分布的随机数区间上均匀分布的随机数1.从从U0,1抽取随机数抽取随机数；2.令令F(x)=;3.解方程得解方程得x：)(1 Fx注：需要知道累积分布函数的解析表达式，且累积分布函数的注：需要知道累积分布函数的解析表达式，且累积分布函数的反函数存在反函数存在分离型随机变量的抽样分离型随机变

8、量的抽样直接抽样法适应于分离型的随机变量直接抽样法适应于分离型的随机变量KiixNiiixpxdxxpxdxfxFy11)()()(Niiixxpxf1)()(方法：方法：1.计算计算yk=yk-1+pk，k=2,3,N,y1=p12.2.从从U0,1U0,1抽取随机数抽取随机数；3.3.求满足求满足yk-1 yk 的的K K值值;4.4.随机变量的第随机变量的第k k个取值即为欲抽取的值。个取值即为欲抽取的值。0 1p apbpc pdp3=0.2b3+c3p2=0.3b2+c2p1=0.5b1+c1a13.05.03.05.05.05.00332211cbcbcb例例1、粒子衰变末态的随机

9、抽样、粒子衰变末态的随机抽样设粒子设粒子a有三种衰变方式，其分支比如下有三种衰变方式，其分支比如下随机选取每次衰变的衰变方式（衰变道）随机选取每次衰变的衰变方式（衰变道）直接抽样法直接抽样法 U0,1rnnrnrnrnnrpprnpnrBrnr)!(!,2,1,0,)1(),;(例例2、二项式分布的抽样、二项式分布的抽样方法方法1：利用上面介绍的直接抽样法，需计算累积分布函数，利用上面介绍的直接抽样法，需计算累积分布函数，当当n很大时，求和计算困难很大时，求和计算困难;方法方法2：利用二项式分布的定义利用二项式分布的定义1.产生产生n个个 i U0,1;2.统计满足条件统计满足条件 i p（表

10、示成功）的（表示成功）的 i的数目的数目r，则，则r表示在表示在n次实验中成功的次数次实验中成功的次数r即为二项式分布的即为二项式分布的抽样值抽样值0,1,2,=r,e r!1=)P(r;-rnppn,0,例例3、泊松分布的抽样、泊松分布的抽样方法方法1：利用直接抽样法，但计算累积分布函数时非常复杂利用直接抽样法，但计算累积分布函数时非常复杂方法方法2：利用泊松分布的定义：二项式分布的极限形式利用泊松分布的定义：二项式分布的极限形式1.选取足够大的选取足够大的n，使，使p=/n相当小，例如，相当小，例如，p=0.12.产生产生n个个 i U0,1;3.统计满足条件统计满足条件 i Divide

11、(1,2);TH1F*h1=new TH1F(h1,h1,100,-5.0,5.0);TH1F*h2=new TH1F(h2,h2,100,-5.0,5.0);SetSeed(9,11)Const Pi=3.1415926;for(int i=0;i Fill(y1);h2-Fill(y2);c1-cd(1);h1-Draw();c1-cd(2);h2-Draw();8.5 舍选抽样法舍选抽样法（acceptance-rejection sampling)直接抽样法的困难：直接抽样法的困难：许多随机变量的累积分布函数无法用解析函数给出许多随机变量的累积分布函数无法用解析函数给出;有些随机变量的

12、累积分布函数的反函数不存在或难以求出；有些随机变量的累积分布函数的反函数不存在或难以求出；即使反函数存在，但计算困难即使反函数存在，但计算困难舍选抽样法：舍选抽样法：抽取随机变量抽取随机变量x的一个随机序列的一个随机序列xi,i=1,2,按一定的舍选规则按一定的舍选规则从中选出一个子序列，使其满足给定的概率分布从中选出一个子序列，使其满足给定的概率分布.假定：假定：随机变量随机变量x的值域为的值域为a,b;x的概率密度函数：的概率密度函数：f(x)=P*(x)/Z,（其中（其中Z为归一化因子）为归一化因子）难以直接抽样难以直接抽样;Q(x)=Q*(x)/ZQ 是另外一个较为简单的函数是另外一个

13、较为简单的函数(ZQ为归一化因子）为归一化因子）可用简单的方法进行抽样可用简单的方法进行抽样;在在x的整个取值范围内：的整个取值范围内：cQ*(x)P*(x),其中其中c为一常数为一常数抽样方法：抽样方法：1.产生两个随机数产生两个随机数 从从Q(x)分布抽取分布抽取 ，a,b;由由0,cQ*()区间上的均匀分布抽取区间上的均匀分布抽取u，u=cQ*(),U0,1.2.接收或舍弃取样值接收或舍弃取样值 x.如果如果 u P*(x)，舍弃，返回到，舍弃，返回到1,重复上述过程重复上述过程;否则，接受否则，接受;几何解释几何解释:在二维图上，随机选取位于曲线在二维图上，随机选取位于曲线cQ*(x)

14、下的点下的点x,u；选取位于曲线选取位于曲线P*(x)下的那些点，则这些点将服从概率密度为下的那些点，则这些点将服从概率密度为f(x)的分布的分布常数常数c的选取的选取 常数常数c应尽可能地小应尽可能地小,因为抽样效率与因为抽样效率与c成反比成反比;c=maxP*(x)/Q*(x),x a,b特例特例:如果取如果取Q(x)=1,x a,b,即均匀分布即均匀分布,则则 X的抽样的抽样:x=(b-a)r+a,r U0,1 c可取为可取为f(x)在在a,b区间上的极大值区间上的极大值abxf(x)c抽取抽取r1,r2 U0,1=a+(b-a)r1r2 f()/cx=2/221)(*xexp2111)

15、(*xxQ52.12)1(2max)(*/)(*max2/22eexxQxPcx)arctan(arctanarctan21)()(11arctan21)(*/)(*)(2axaxdxQxFxadxxQxQxQxaQaa例例1:标准正态分布的抽样标准正态分布的抽样,x-a,a无法用直接抽样法无法用直接抽样法,累积分布函数无解析累积分布函数无解析表达式表达式Breit-wigner or Cauchy分布分布)arctan(arctan2tan 1,0aaxUQuUxxcQQ 1,01152.1)(*2由由Q(x)抽取抽取x 直接抽样法直接抽样法抽取抽取u计算计算P*(x),如果如果u=P*(x

16、),接受接受xfloat gaussian_reject(double a)const float c=1.52;while(true)float eta=randac();float x=tan(eta*2.0*atan(a)+atan(-a);float q=c*1/3.1415926*1.0/(1+x*x);float ksi=randac();float u=ksi*q;float p=1/sqrt(2*3.1415926)*exp(-x*x/2.0);if(u Divide(1,2);TH1F*h1=new TH1F(h1,h1,100,-5.0,5.0);for(int i=0;i

17、 Fill(x);c1-cd(2);h1-Draw();2cos2sin2sin2sin2coscos22SCAB/222222222222sin,cos2sin,2cosBAABSBABACBABBAA例例2：利用舍选法产生随机数：利用舍选法产生随机数C=cos,S=sin,其中其中 为为0,2 区间内区间内均匀分布的随机数均匀分布的随机数方法方法1：先产生：先产生0,2 间均匀分布的随机数：间均匀分布的随机数：=2 r，r U0,1,然后直接计算然后直接计算C和和S 因需要计算三角函数，因需要计算三角函数，故此方法运算速度慢故此方法运算速度慢方法方法2 2：利用舍选法可避免三角函数运算：利

18、用舍选法可避免三角函数运算令令A和和B为单位圆内直角三角形的两个边，则有为单位圆内直角三角形的两个边，则有因此，只要产生单位圆内的随机坐标因此，只要产生单位圆内的随机坐标A和和B,就可用代数运算就可用代数运算求出求出C和和S,算法为算法为1.产生两个产生两个0,1区间上均匀分布的随机数区间上均匀分布的随机数u1和和u2;2.令令v1=2u1-1,v2=u2,则则v1 U-1,1,v2 U0,1;3.计算计算r2=v12+v22,如果如果r2 1,转到转到1，重新产生；，重新产生；4.令令A=v1,B=v2,计算计算C和和S8.6 复合分布的抽样方法复合分布的抽样方法（composition m

19、ethod)1)(;1;0)()(dxxfppxfpxfkkkkkkkkiikiiprpUr111;1,01961年由年由Marsaglia提出的方法提出的方法设随机变量设随机变量X的概率密度函数的概率密度函数f(x)可写成一些可写成一些PDF的线性叠加：的线性叠加：抽样方法：抽样方法：1.利用离散型的随机变量的抽样方法抽取序号利用离散型的随机变量的抽样方法抽取序号k;2.由由fk(x)抽取抽取x)0()()0()(5.0001)(001)()(5.0)(5.0)(5.0)(:2121)0)0,(),0)0,(xexfxexfppxxIxxIxIexIexfexfPDFxxxxx其他其他，例：

20、用复合法产生双指数分布随机数例：用复合法产生双指数分布随机数 产生两个产生两个0,1区间均匀分布的随机数区间均匀分布的随机数 r1和和r2;如果如果r1 0.5,按按f1(x)抽样；抽样；如果如果r1 0.5,按按f2(x)抽样；抽样；f1(x)和和f2(x)都可用直接抽样法都可用直接抽样法8.7 混合抽样法混合抽样法（mixed method)混合法：混合法：直接抽样法直接抽样法舍选法舍选法亦称乘抽样法亦称乘抽样法适用的抽样场合：适用的抽样场合：直接抽样法不可用、舍选抽样法效率低直接抽样法不可用、舍选抽样法效率低概率密度函数难以积分概率密度函数难以积分无累积分布函数的解析表达式；无累积分布函

21、数的解析表达式；累积分布函数的反函数的解析表达式不存在；累积分布函数的反函数的解析表达式不存在；概率密度函数存在尖峰（概率密度函数存在尖峰（spiky);)()()(xgxfxp假定：概率密度函数可写成下面的因子化形式假定：概率密度函数可写成下面的因子化形式其中：其中：f(x)包含了包含了p(x)的峰值部分且可用直接抽样法进行抽样的峰值部分且可用直接抽样法进行抽样 g(x)是一个相对变化平缓的函数，包含了是一个相对变化平缓的函数，包含了p(x)函数的大函数的大部分的数学复杂性；部分的数学复杂性；)(1)(;)(xfCxfCdxxffbaf,1)();(1)(baxxgxgMxgg抽样方法抽样方

22、法：1.将将f(x)归一化：归一化：2.令令Mg为为g(x)在区间在区间a,b上的极大值上的极大值3.利用直接抽样法由利用直接抽样法由 抽取抽取x;)(xf4.抽取抽取0,1区间均匀分布的随机数区间均匀分布的随机数r,如果如果 ，则接受，则接受x,否则，返回到否则，返回到3重新抽样重新抽样rxg)()()()(xgxfxpiiiiiiiiiiixpCCxgxfCxp)()()()()(xpi)(xpi设概率密度函数可写成如下的形式：设概率密度函数可写成如下的形式：推广的形式推广的形式：复合抽样法复合抽样法+混合抽样法混合抽样法=乘加抽样法乘加抽样法抽样方法：抽样方法：1.采用复合抽样法，先确定

23、采用复合抽样法，先确定p(x)的随机数应由哪一个的随机数应由哪一个 抽取；抽取；2.在按在按 抽样时，采用上面的混合方法抽样时，采用上面的混合方法1)()()()(;)()(iibaiiiibaiiiCdxxpCxgxfxpdxxgxfC1,2;)cos1(;1sin11),(00220min1102201102220220010EcmcmEEEcmcmEEEEEcmrnXEEeeeee例：例：compton散射散射微分截面：微分截面：如何抽取散射光子的能量？如何抽取散射光子的能量？=乘加抽样法乘加抽样法1)(0:1,sin11)()(1)()()(;2)1(ln1)()(1)(;ln1ln)

24、()()()()(1sin11)(0222022222220201111100122112200ggdfFfdfFfgffgf 抽样方法：抽样方法：r1,r2,r3是三个在是三个在0,1区间上均匀分布的随机数区间上均匀分布的随机数1.确定由哪一个确定由哪一个f函数来抽取函数来抽取：如果：如果r1抽样方法：抽样方法：8.9 多维分布的抽样多维分布的抽样 将一维分布的抽样方法推广到多维分布将一维分布的抽样方法推广到多维分布设设n维概率密度为维概率密度为 ，的取值域为的取值域为n维长方体维长方体 ,)(xf,bax,2121nnbbbbaaaa)(xf的极大值为的极大值为M,则则 的随机数的随机数

25、的舍选法的舍选法)(xf,21n抽取抽取r,ri U0,1,i=1,2,ni=ai+(bi-ai)ri,i=1,2,n Mfr/)(x舍选法用于多维抽样舍选法用于多维抽样混合抽样法用于多维分布的抽样混合抽样法用于多维分布的抽样设多维随机变量概率密度函数设多维随机变量概率密度函数 可表示为可表示为)(xf)()()(xgxHxf式中式中 为任意多维概率密度函数，其随机数为为任意多维概率密度函数，其随机数为)(xg,21nH()0的极大值为的极大值为M，则，则f()的随机数的随机数 可可由下列方式抽取：由下列方式抽取：xx,21n抽取抽取r和和MHr/)(条件密度法条件密度法:利用条件概率密度将多

26、维模拟转化为一维模拟问题利用条件概率密度将多维模拟转化为一维模拟问题.任意任意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联合概联合概率密度率密度nnnnndxxxxfxxxf),(),(211211fn(x1,x2,xn)=f(xn|x1,x2,xn-1)f(xn-1|x1,x2,xn-2)f(x2|x1)f(x1)n n个一维概率密度的个一维概率密度的乘积乘积),.,(),.,(),.,|(121121121nnnnnnxxxfxxxfxxxxf条件密度法条件密度法:fn(x1,x2,xn)=f(xn|x1,x2,xn-1)f(xn-1|x1,x2,xn-2)f(x2|x1)f(x1)抽样方法：抽样方法：由由f(x1)抽样得到抽样得到x1的随机数的随机数 1，将，将 1代入代入f(x2|x1)中得到中得到x2的一维分布的一维分布f(x2|1)；由由f(x2|1)抽样得到抽样得到x2的随机数的随机数 2，