1、6/28/2008,BESIII 暑期讲习班,1,参数估计与拟合,陈少敏清华大学工程物理系,6/28/2008,BESIII 暑期讲习班,2,2018-10-08,2,主要内容,参数估计最大似然法最小二乘法矩方法MINUIT的使用,6/28/2008,BESIII 暑期讲习班,3,粒子物理实验研究的目的,粒子物理实验的一个重要目的是确定粒子的属性,6/28/2008,BESIII 暑期讲习班,4,参数估计,理论,数据,数据,理论,概率,微积分,给定含有参数 的理论预言分布，对数据能下何种结论？,需要一个步骤来从数据 D 中估计参数 ,最常用的步骤是“拟合”。,统计分析,6/28/2008,BE
2、SIII 暑期讲习班,5,什么是估计量？,一个估计量对应于这样一个步骤，它能从实际数据测量值中对一个参数或一个分布的属性给出定量的结果,6/28/2008,BESIII 暑期讲习班,6,如何评判一个参数估计量的好坏？,符合程度(一致性),偏置大小(无偏性),方差大小(有效性),6/28/2008,BESIII 暑期讲习班,7,参数估计与概率大小的关系,如果假设(包括 的取值)为真,可以预料会使观测结果具有高的概率。,如果假设的 取值远离真值,会使观测结果具有低的概率。,6/28/2008,BESIII 暑期讲习班,8,似然函数,在经典统计理论里，L( )并不是 的概率密度函数。,根据参数好坏与
3、概率大小的关系，可以认为真实的 应使得下式定义的似然函数,有大的数值。,6/28/2008,BESIII 暑期讲习班,9,最大似然估计量,取最大值,有时候 L( ) 可以有好几个极大值,注意，1）方法利用了所有信息，与如何划分数据分布区间无关； 2）定义的最大似然估计量并不保证它们总是最优的。,需要对诸如无偏性，有效性等问题进行研究,大样本情况下，最大似然法大都能给出了期待的好结果。,小样本的情况，虽然不总是最优，但也能给出最好的实用解。,6/28/2008,BESIII 暑期讲习班,10,最大似然估计量的唯一性,考虑 的最大似然估计值是下列方程的解,选用等价参数 h( ),因为,因此，h 的
4、最大似然估计值与参数选取无关，具有唯一性。,6/28/2008,BESIII 暑期讲习班,11,指数概率密度函数的参数估计,考虑指数概率密度函数,设有数据样本 t1,tn。为方便起见采用对数形式(对同样的参数值，该定义不改变最大值的位置)。,是平均寿命的最大似然估计,6/28/2008,BESIII 暑期讲习班,12,指数型最大似然估计是无偏的,6/28/2008,BESIII 暑期讲习班,13,估计量的方差：数值方法,指数分布平均值的估计量为：,6/28/2008,BESIII 暑期讲习班,14,估计量的方差：蒙特卡罗方法,Nexp=1000,蒙特卡罗模拟可给出标准偏差,6/28/2008,
5、BESIII 暑期讲习班,15,估计量的方差： RCF边界法,任何估计量(不仅仅是最大似然法)的方差下界为,也称为 Rao-Cramr-Frechet 不等式。,通常假设上述结论为真，利用RCF边界估计,( b为偏置),最大似然估计量对大的样本统计量 n 几乎总是有效的。,6/28/2008,BESIII 暑期讲习班,16,指数型函数估计量的 RCF边界,已知估计是无偏的 b = 0，所以由RCF边界可得,=真的方差,对指数型函数求二阶导数,6/28/2008,BESIII 暑期讲习班,17,RCF边界与HESSE矩阵,对于只有一个参数的情况，可以得到,求 logL 的最大值可通过数值计算来完
6、成，二阶导数的矩阵(Hessian矩阵)是通过有限差值来估计。,调用 CERN 的 MINUIT 软件包中的 HESSE 程序,6/28/2008,BESIII 暑期讲习班,18,例子：估计实验所需的统计量,质子与反质子弹性散射实验，观测量为散射角 x=cos，服从 f (x;)=0.5(1+x)，其中 是反映反质子极化的参数。目前测量值为0.100.02 ，要想在统计上将相对误差减少到 5%，总共需要多少个事例？,由信息不等式，任何估计量的方差下界为,对于本问题，b=0，,6/28/2008,BESIII 暑期讲习班,19,例子：实验所需的统计量（续）,由信息不等式，,带入不等式可以求得 n
7、 1.2 105。,6/28/2008,BESIII 暑期讲习班,20,估计量的方差: 图解法,也就是,6/28/2008,BESIII 暑期讲习班,21,推广的最大似然法,21,如果考虑了样本大小 n 也是泊松分布的随机变量，平均值为，那么,6/28/2008,BESIII 暑期讲习班,22,BES上的 质量测量,22,能量点 i,PRL 69, 3021 (1992),6/28/2008,BESIII 暑期讲习班,23,推广的最大似然法： 独立,23,可以将其分解成单独求 与 估计值的问题。,例如：是信号与本底分量的叠加。,6/28/2008,BESIII 暑期讲习班,24,推广的最大似然
8、法： 独立(续),24,根据概率的定义可知并非所有的 i 独立，因此有,在推广的最大似然法中,在总事例数 n 中，第 i 部分事例数为 i = i,如果联合概率密度函数可以表示为,6/28/2008,BESIII 暑期讲习班,25,可能出现的负值问题,25,假设有两类事例：信号 (s) 与本底 (b),问题：如果出现负值，应该如何报告结果？,6/28/2008,BESIII 暑期讲习班,26,最大似然法处理分区数据,26,如果样本的联合概率密度函数( ntot 为常数)为,通常称为“对直方图拟合”。,在某种假设下，有期待值,6/28/2008,BESIII 暑期讲习班,27,分区处理数据可能出
9、现的问题,27,分区处理的数据满足, 对函数变化影响较小。,因此，对直方图拟合，一定要确认区间的大小对结果无明显影响。注意：区间无穷小时，与点估计结果一致。,6/28/2008,BESIII 暑期讲习班,28,似然函数与贝叶斯定理,28,贝叶斯方法,对假设 ( ) 采用主观概率来描述实验前，对过程的了解由 归纳(先验概率密度函数)利用贝叶斯定理，根据数据来改进先验概率密度函数,6/28/2008,BESIII 暑期讲习班,29,似然法估计量与贝叶斯估计量,29,对很纯粹的贝叶斯论者,实际应用中，会碰到如何得到 的问题?,对很实际的贝叶斯论者,无金科玉律(主观的!)，通常假定它在全区域是均匀分布
10、。,非贝叶斯方法目前在统计学中占主要地位。但是，如果被测定参数是随机变量，而且它的验前分布已知，就应该用贝叶斯方法。,6/28/2008,BESIII 暑期讲习班,30,最小二乘法与最大似然法,设有高斯随机变量: yi, i=1,N ,均值为,对应的对数似然函数(去掉与 无关的项)为,对于独立的高斯变量 yi，联合概率密度函数为,6/28/2008,BESIII 暑期讲习班,31,最小二乘估计量的定义,如果 yi 是一多维高斯变量，协方差矩阵为V ，满足,那么其对数似然函数为,也就是说，我们应求下式的最小值,它的最小值定义了最小二乘法的估计量 ，即使 yi 不是高斯变量，该定义依然适用。,6/
11、28/2008,BESIII 暑期讲习班,32,两种情况下的最小二乘参数估计,对参数的估计可以根据理论预期值中所含参数的具体特征而采用不同的参数估计处理方法,线性情况：,非线性情况：,6/28/2008,BESIII 暑期讲习班,33,线性情况的最小二乘法估计,这里 aj(x) 是 x 的任意线性独立函数。,用矩阵来表示时，令 Aij=aj(xi)，有,对 i 求偏微分，并令结果等于零，有,解方程得到最小二乘法的估计量,6/28/2008,BESIII 暑期讲习班,34,非线性情况的最小二乘法估计,如果采用牛顿法求上式的最小值，第 n+1次迭代公式可采用,6/28/2008,BESIII 暑期
12、讲习班,35,最小二乘估计量的方差,等效地，可以利用下式来计算,如果 yi是高斯变量时，其与RCF边界一致。,6/28/2008,BESIII 暑期讲习班,36,最小二乘估计量的方差(续),6/28/2008,BESIII 暑期讲习班,37,约束情况下的最小二乘法拟合,实际问题中会遇到测量量本身要受到某些物理定律的约束。,求解可采用拉格朗日乘子法，对每一个约束引入修正因子i，,例如，能动量守恒，衰变顶点约束等等。对一个事例有m个观测量，无参数的最小二乘问题变为,6/28/2008,BESIII 暑期讲习班,38,约束情况下的最小二乘法(续一),为了找到最小值，可以通过求微商方法,而 n+1 次
13、迭代后,设经过 n 次迭代以后，找到一组解 ，得到函数 的值。在 上对 (n) 进行线性展开，并略去高阶项，得到,6/28/2008,BESIII 暑期讲习班,39,约束情况下的最小二乘法(续二),两式联立消掉 项，可以得到,因此，可以得到第 n+1 次迭代的 l 个拉格朗日乘子取值,以及第 n+1 次迭代的 m 个测量量的预期值,6/28/2008,BESIII 暑期讲习班,40,约束情况下的最小二乘法(续三),当经过 n+1 次迭代以后，满足下式时即可终止,实验中，为了提高测量精度而采用的四动量守恒约束拟合(4-C fit)，顶点或质量约束拟合(1-C fit)，大都采用该方式来进行。,此
14、时的 2 值应满足自由度为（m - l）的 2 分布。,6/28/2008,BESIII 暑期讲习班,41,例子：粒子动量分辨的改进,例如，实验观测衰变,通常情况下，探测器对光子探测的能量分辨率较差，从而影响到 0 粒子动量重建的精度。,已知：,6/28/2008,BESIII 暑期讲习班,42,例子：粒子动量分辨的改进（续）,因此，每一个衰变事例的观测量期待值为,对应于每个观测量有误差估计，而且已知相互间不相关。则无参数的最小二乘问题可写为为,利用一个约束条件下，改进的光子动量观测值进行 重建研究，从 的不变质量谱可以看出光子的动量得到了明显的改进。,6/28/2008,BESIII 暑期讲
15、习班,43,检验最小二乘法的拟合优度,那么2min 服从 N-m自由度的最小二乘概率密度函数分布。,据此来计算P-值,如在五个数据点的双参数拟合,也就是说，重复实验多次，有 26.3% 的值将大于 2min 。,进行 1000 次蒙特卡罗实验,6/28/2008,BESIII 暑期讲习班,44,最小二乘法处理分区数据,最小二乘法拟合使下式有最小值,把 yi 看做泊松分布的随机变量，方差为,改进的最小二乘法虽方便了计算，但对于有些区间频数太少时2min不再服从最小二乘的概率密度分布函数(或无定义)。,6/28/2008,BESIII 暑期讲习班,45,最小二乘法的归一化常数问题,例如 n = 4
16、00次，N = 20个区间,解决的方法是从数据中直接得到 n ，或者最好是用最大似然法定n。,6/28/2008,BESIII 暑期讲习班,46,矩的一般表达式,假设对随机变量 x 有n 次测量 x1,xn，服从概率密度函数分布 f (x; ) 。其中有 m 个未知参数 1,m。如果可以构造 m 个线性独立函数 ai(x), i=1,m，其均值可写为,为了确定参数，上述独立函数必须进行适当选择使得含参数的函数 ei( ) 可以确定。,此时，函数 ei( ) 可以通过计算无偏的样本平均值来估计,因此，参数值可以通过求解 m 个 ei( ) 方程组来确定。,矩的一般表达式,6/28/2008,BE
17、SIII 暑期讲习班,47,线性独立函数的方差矩阵,参数 1,m 估计值的协方差矩阵的无偏估计,它可以与样本平均值的协方差矩阵相联系, 即,6/28/2008,BESIII 暑期讲习班,48,线性独立函数均值的方差矩阵,根据线性独立函数均值和其估计值的定义，可以有,6/28/2008,BESIII 暑期讲习班,49,参数估计值的方差矩阵,由于待定参数 是 e 的函数，由误差传递,因此待定参数 的协方差矩阵估计值也可以确定。而根据线性独立函数的均值估计值表达式,可知参数值的估计值 可通过求解 m 个 方程组来确定。,参数估计任务完成 ,6/28/2008,BESIII 暑期讲习班,50,矩方法估
18、计参数,计算出 cos 二阶代数矩的理论期待值,则参数 与二阶代数矩的关系为,只要函数是可积的，采用矩方法原则上就可以测定参数。,在实验 中，理论预言角分布为,6/28/2008,BESIII 暑期讲习班,51,最大似然法、最小二乘法和矩,充分性：估计量应包含观测值对于未知参数的全部信息；一致性：样本容量增大时，估计值收敛于真值；有效性：估计量的分布对其期望值具有最小方差；无偏性：无论样本容量多大，估计值与真值无系统偏差。,6/28/2008,BESIII 暑期讲习班,52,拟合与 MINUIT 软件包,52,在粒子物理与核物理实验的参数估计中，利用计算机求极值的数值解越来越普遍。通用的求函数
19、最小值程序是,在使用MINUIT时，为了对结果有正确的诠释，必须对相关输出信息进行理解。,MINUIT 软件包,求 log L 最大值等效于求 log L 的最小值。,在 MINUIT 框架下单独使用(适于data-driven 模式)；在 PAW 环境下互动调用 MINUIT (基于Fortran)；在 ROOT 环境下互动调用 MINUIT (基于C+)。,6/28/2008,BESIII 暑期讲习班,53,ROOT 平台下的 MINUIT 调用,53,const int npoints=2000;Double_t xnpoints;Double_t angle_cut=0.95;void
20、 minuit_fit() / Get data points get_input_data(); / Prepare for fit const int npar=2; TMinuit *gMinuit = new TMinuit(npar); gMinuit-SetFCN(fcn); Double_t arglist10; Int_t ierflg = 0; arglist0 = 1; gMinuit-mnexcm(SET ERR, arglist ,1,ierflg);/ Set starting values and step sizes for parameters Double_t
21、 vstartnpar = 0.5, 0.5 ; Double_t stepnpar = 0.1 , 0.1 ; gMinuit-mnparm(0, alpha, vstart0, step0, 0,0,ierflg); gMinuit-mnparm(1, beta, vstart1, step1, 0,0,ierflg);,/ Now ready for minimization step arglist0 = 500; arglist1 = 1.; gMinuit-mnexcm(MIGRAD, arglist ,0,ierflg); gMinuit-mnexcm(HESSE, arglis
22、t ,0,ierflg); gMinuit-mnexcm(“MINOS, arglist ,0,ierflg);/ Print results Double_t fmin,fedm,errdef,covmatnparnpar; Double_t alpha,alpha_err,beta,beta_err; Int_t nvpar,nparx,icstat; gMinuit-mnstat(fmin,fedm,errdef,nvpar,nparx,icstat); gMinuit-GetParameter(0, alpha, alpha_err ); gMinuit-GetParameter(1,
23、 beta, beta_err ); gMinuit-mnemat( ,root.x minuit_fit.C,6/28/2008,BESIII 暑期讲习班,54,使用 MINUIT 应注意的问题,54,求极值问题时，不可能自动找到合理的初值。,原因：似然函数存在多个极值。,提供较好的参数不确定范围可以有助于找到真值。太大的范围会使碰巧在远离真值处找到一个区域最小值。,6/28/2008,BESIII 暑期讲习班,55,通常在MINUIT的解决方案,55,利用MINUIT函数 MIGRAD 找 log L 的最小值利用MINUIT函数 HESSE 计算参数的误差还可利用MINUIT函数 MIN
24、OS 做更好的误差估计,6/28/2008,BESIII 暑期讲习班,56,MINUIT函数 MIGRAD,56,6/28/2008,BESIII 暑期讲习班,57,函数 MIGRAD（续一）,57,6/28/2008,BESIII 暑期讲习班,58,函数 MIGRAD（续二）,58,6/28/2008,BESIII 暑期讲习班,59,MINUIT函数 HESSE,59,6/28/2008,BESIII 暑期讲习班,60,函数 HESSE （续一）,60,6/28/2008,BESIII 暑期讲习班,61,函数 HESSE （续二）,61,6/28/2008,BESIII 暑期讲习班,62,函
25、数 HESSE （续三）,62,6/28/2008,BESIII 暑期讲习班,63,MINUIT函数 MINOS,63,用途：更好的误差计算方法注意：有太多参数时可能会耗费非常多的CPU时间调用：在ROOT 或 PAW 使用“E” 选项,6/28/2008,BESIII 暑期讲习班,64,经常会碰到的问题：不收敛,64,时常会出现参数拟合过程不收敛，例如MIGRAD 不能找到最小值HESSE 得到负的二阶导数很多情况下是参数拟合函数有误（程序书写错误或理论模型不正确）数值计算中出现精度问题数值计算中出现稳定性问题HESSE 相关系数矩阵是很好的调查起点,表明两个参数几乎100%关联，求最小值过
26、程无法建立。,6/28/2008,BESIII 暑期讲习班,65,消除影响稳定性的因素：相关性,65,策略一：更恰当地选取拟合参数例如，对于采用类似的双高斯宽度拟合,事例数,HESSE 相关系数矩阵,双高斯宽度的大小与它们所占的份额有很强的关联,6/28/2008,BESIII 暑期讲习班,66,影响稳定性因素：相关性（续）,66,可以采用不同参数化过程来解决,第二个高斯分布宽度与份额比的相关度从0.92减小至0.68 完全是参数化过程自身的问题,策略二：固定所有高度相关的参数而只保留一个如果参数高度相关表明它们中有些是多余的，不会对拟合模型的自由度带来贡献。,6/28/2008,BESIII
27、 暑期讲习班,67,消除影响稳定性的因素：多项式,67,注意：普通多项式 a0+a1x+a2x2+a3x3 的参数化过程通常会导致在拟合参数 a 之间引入很强的关联。例如，导致拟合稳定性问题而不能找到高阶系数的解。,解决方案：采用数学上的正交多项式，例如，勒让德多项式，第一类契贝谢夫多项式等等,6/28/2008,BESIII 暑期讲习班,68,参数拟合问题：参数的限界,68,有时候需要对拟合参数进行区间的限制例如，份额比参量只能在【0，1】区间，三角函数值只能在【-1，+1】之间，等等。,但是，可能会导致 MINUIT 不能正确估计误差。因此，必须小心使用参数的拟合区间限制。参考的解决方案：
28、如果区间限制的引入可使拟合稳定，最好考虑寻找有否其它可替代的无区间限制的参数化方案。如果区间限制的引入可以避免出现“非物理”但仍为合理统计解时，应考虑不要强加区间限制，并对结果采用统计意义上的解释。一般是同时给出物理解与“非物理”解。前者是使结果易于物理解释，后者是使结果可与其它实验结果进行统计意义上的正确合并。,6/28/2008,BESIII 暑期讲习班,69,举例：加速器中微子振荡实验,Phys.Rev.D74:072003,2006,6/28/2008,BESIII 暑期讲习班,70,如何确认拟合结果的有效性,70,所谓的结果有效性是指,拟合结果无偏拟合误差符合统计不确定性,如果对结果
29、的有效性有疑虑，可以考虑,进行P-值计算，估计结果是否为极端情况。采用蒙特卡罗样本，重复同样统计量的参数拟合过程至少100次，确认拟合的参数值是否有偏向性。如果有偏向性，例如，发现有偏差量 ab，应考虑对结果进行 a 因子修正，并把 b 当作系统误差。检查蒙特卡罗研究中均值的妳散是否与拟合参数的误差符合，以确定误差的合理性。,6/28/2008,BESIII 暑期讲习班,71,举例：拟合信号数目,71,一维直方图信号的拟合,信号服从高斯分布本底分布已知,研究拟合参数 Nsig,产生 Nsig,Nsig(fit),6/28/2008,BESIII 暑期讲习班,72,举例：拟合结果的“Pull”分
30、析,72,如何判断误差是否合理？,(Nsig),从蒙特卡罗研究中很难进行诠释没有与 Nsig 等效的分布研究误差,解决方案是检查“Pull”分布,定义Pull 的特点本例是一个在统计精度内无偏而且误差正确的拟合。,如果拟合结果无偏则均值为0；如果拟合误差正确则宽度为1。,Pull(Nsig),6/28/2008,BESIII 暑期讲习班,73,低统计下的拟合结果检查,73,对低统计样本或者大样本小信号情况需要特别注意,可能出现的问题,似然评估量不再有效，导致从二阶导数中估计误差不再准确最小二乘法中误差不再服从高斯分布而是泊松分布偏置项在评估量中受统计量影响1/N部分不再可以忽略,一般地，偏置的
31、出现使得拟合误差正确性会有问题,采用点估计而不是分区的最大似然法拟合尽最大可能检查拟合的有效性,6/28/2008,BESIII 暑期讲习班,74,举例：小样本情况的偏向性,74,低统计量情况举例,Nsig=10020， Nbg=200,蒙特卡罗模拟研究,分布在低统计下已经不再是对称分布,Pull 的均值已经不再是 0，有2.3 的偏离,必须对结果进行修正并引入恰当的误差。,6/28/2008,BESIII 暑期讲习班,75,总结,如何比较理论与实验,当理论含待定参数时，如何构造参数的估计量,如何判断参数估计量的好坏：三个标准,参数估计量的好坏与样本概率大小的关系,似然函数、最大似然函数,估计量方差的估计：四种方法,最小二乘法、矩方法,使用 MINUIT 拟合时如何诠释结果,
