参数估计与置信区间课件

1、参数估计与拟合陈少敏清华大学工程物理系6/28/20081*第1页，共75页。8/5/20222主要内容参数估计最大似然法最小二乘法矩方法MINUIT的使用6/28/20082*第2页，共75页。粒子物理实验研究的目的粒子物理实验的一个重要目的是确定粒子的属性质量、寿命（宽度）、分支比、自旋、宇称高斯布莱-魏格纳指数（时间）布莱-魏格纳（宽度）角动量守恒 J=L+S球谐函数或 d 函数 （角分布）克莱布施 戈丹系数衰变或散射振幅6/28/20083*第3页，共75页。参数估计理论数据数据理论概率微积分给定含有参数 的理论预言分布，对数据能下何种结论？需要一个步骤来从数据 D 中估计参数 最常用

2、的步骤是“拟合”。统计分析6/28/20084*第4页，共75页。什么是估计量？一个估计量对应于这样一个步骤，它能从实际数据测量值中对一个参数或一个分布的属性给出定量的结果6/28/20085*第5页，共75页。如何评判一个参数估计量的好坏？符合程度(一致性)偏置大小(无偏性)方差大小(有效性)6/28/20086*第6页，共75页。参数估计与概率大小的关系如果假设(包括 的取值)为真可以预料会使观测结果具有高的概率。如果假设的 取值远离真值会使观测结果具有低的概率。6/28/20087*第7页，共75页。似然函数在经典统计理论里，L( )并不是 的概率密度函数。根据参数好坏与概率大小的关系，

3、可以认为真实的 应使得下式定义的似然函数有大的数值。6/28/20088*第8页，共75页。最大似然估计量取最大值有时候 L( ) 可以有好几个极大值注意，1）方法利用了所有信息，与如何划分数据分布区间无关； 2）定义的最大似然估计量并不保证它们总是最优的。需要对诸如无偏性，有效性等问题进行研究大样本情况下，最大似然法大都能给出了期待的好结果。小样本的情况，虽然不总是最优，但也能给出最好的实用解。6/28/20089*第9页，共75页。最大似然估计量的唯一性考虑 的最大似然估计值是下列方程的解选用等价参数 h( ) 因为因此，h 的最大似然估计值与参数选取无关，具有唯一性。6/28/20081

4、0*第10页，共75页。指数概率密度函数的参数估计考虑指数概率密度函数设有数据样本 t1,tn。为方便起见采用对数形式(对同样的参数值，该定义不改变最大值的位置)。是平均寿命的最大似然估计6/28/200811*第11页，共75页。指数型最大似然估计是无偏的6/28/200812*第12页，共75页。估计量的方差：数值方法指数分布平均值的估计量为： 6/28/200813*第13页，共75页。估计量的方差：蒙特卡罗方法Nexp=1000蒙特卡罗模拟可给出标准偏差6/28/200814*第14页，共75页。估计量的方差： RCF边界法任何估计量(不仅仅是最大似然法)的方差下界为也称为 Rao-C

5、ramr-Frechet 不等式。通常假设上述结论为真，利用RCF边界估计( b为偏置)最大似然估计量对大的样本统计量 n 几乎总是有效的。6/28/200815*第15页，共75页。指数型函数估计量的 RCF边界已知估计是无偏的 b = 0，所以由RCF边界可得=真的方差对指数型函数求二阶导数6/28/200816*第16页，共75页。RCF边界与HESSE矩阵对于只有一个参数的情况，可以得到求 logL 的最大值可通过数值计算来完成，二阶导数的矩阵(Hessian矩阵)是通过有限差值来估计。调用 CERN 的 MINUIT 软件包中的 HESSE 程序6/28/200817*第17页，共7

6、5页。例子：估计实验所需的统计量质子与反质子弹性散射实验，观测量为散射角 x=cos，服从 f (x;)=0.5(1+x)，其中 是反映反质子极化的参数。目前测量值为0.100.02 ，要想在统计上将相对误差减少到 5%，总共需要多少个事例？由信息不等式，任何估计量的方差下界为 对于本问题，b=0， 6/28/200818*第18页，共75页。例子：实验所需的统计量（续）由信息不等式，带入不等式可以求得 n 1.2 105。6/28/200819*第19页，共75页。估计量的方差: 图解法也就是6/28/200820*第20页，共75页。推广的最大似然法21如果考虑了样本大小 n 也是泊松分布

7、的随机变量，平均值为，那么6/28/200821*第21页，共75页。BES上的 质量测量22 能量点 iPRL 69, 3021 (1992) 6/28/200822*第22页，共75页。推广的最大似然法： 独立23可以将其分解成单独求 与 估计值的问题。例如：是信号与本底分量的叠加。6/28/200823*第23页，共75页。推广的最大似然法： 独立(续)24根据概率的定义可知并非所有的 i 独立，因此有在推广的最大似然法中 在总事例数 n 中，第 i 部分事例数为 i = i 如果联合概率密度函数可以表示为6/28/200824*第24页，共75页。可能出现的负值问题25假设有两类事例：

8、信号 (s) 与本底 (b) 问题：如果出现负值，应该如何报告结果？6/28/200825*第25页，共75页。最大似然法处理分区数据26如果样本的联合概率密度函数( ntot 为常数)为通常称为“对直方图拟合”。在某种假设下，有期待值6/28/200826*第26页，共75页。分区处理数据可能出现的问题27分区处理的数据满足 对函数变化影响较小。因此，对直方图拟合，一定要确认区间的大小对结果无明显影响。注意：区间无穷小时，与点估计结果一致。6/28/200827*第27页，共75页。似然函数与贝叶斯定理28贝叶斯方法对假设 ( ) 采用主观概率来描述实验前，对过程的了解由 归纳(先验概率密度

9、函数)利用贝叶斯定理，根据数据来改进先验概率密度函数6/28/200828*第28页，共75页。似然法估计量与贝叶斯估计量29对很纯粹的贝叶斯论者实际应用中，会碰到如何得到 的问题?对很实际的贝叶斯论者无金科玉律(主观的!)，通常假定它在全区域是均匀分布。非贝叶斯方法目前在统计学中占主要地位。但是，如果被测定参数是随机变量，而且它的验前分布已知，就应该用贝叶斯方法。6/28/200829*第29页，共75页。最小二乘法与最大似然法设有高斯随机变量: yi, i=1,N ,均值为对应的对数似然函数(去掉与 无关的项)为对于独立的高斯变量 yi，联合概率密度函数为6/28/200830*第30页，

10、共75页。最小二乘估计量的定义如果 yi 是一多维高斯变量，协方差矩阵为V ，满足那么其对数似然函数为也就是说，我们应求下式的最小值它的最小值定义了最小二乘法的估计量 ，即使 yi 不是高斯变量，该定义依然适用。6/28/200831*第31页，共75页。两种情况下的最小二乘参数估计对参数的估计可以根据理论预期值中所含参数的具体特征而采用不同的参数估计处理方法线性情况：非线性情况：6/28/200832*第32页，共75页。线性情况的最小二乘法估计这里 aj(x) 是 x 的任意线性独立函数。 用矩阵来表示时，令 Aij=aj(xi)，有 对 i 求偏微分，并令结果等于零，有 解方程得到最小二

11、乘法的估计量 6/28/200833*第33页，共75页。非线性情况的最小二乘法估计如果采用牛顿法求上式的最小值，第 n+1次迭代公式可采用6/28/200834*第34页，共75页。最小二乘估计量的方差等效地，可以利用下式来计算如果 yi是高斯变量时，其与RCF边界一致。6/28/200835*第35页，共75页。最小二乘估计量的方差(续)6/28/200836*第36页，共75页。约束情况下的最小二乘法拟合实际问题中会遇到测量量本身要受到某些物理定律的约束。求解可采用拉格朗日乘子法，对每一个约束引入修正因子i，例如，能动量守恒，衰变顶点约束等等。对一个事例有m个观测量，无参数的最小二乘问题

12、变为6/28/200837*第37页，共75页。约束情况下的最小二乘法(续一)为了找到最小值，可以通过求微商方法而 n+1 次迭代后设经过 n 次迭代以后，找到一组解 ，得到函数 的值。在 上对 (n) 进行线性展开，并略去高阶项，得到6/28/200838*第38页，共75页。约束情况下的最小二乘法(续二)两式联立消掉 项，可以得到因此，可以得到第 n+1 次迭代的 l 个拉格朗日乘子取值以及第 n+1 次迭代的 m 个测量量的预期值 6/28/200839*第39页，共75页。约束情况下的最小二乘法(续三)当经过 n+1 次迭代以后，满足下式时即可终止实验中，为了提高测量精度而采用的四动量

13、守恒约束拟合(4-C fit)，顶点或质量约束拟合(1-C fit)，大都采用该方式来进行。此时的 2 值应满足自由度为（m - l）的 2 分布。6/28/200840*第40页，共75页。例子：粒子动量分辨的改进例如，实验观测衰变通常情况下，探测器对光子探测的能量分辨率较差，从而影响到 0 粒子动量重建的精度。已知：r6/28/200841*第41页，共75页。例子：粒子动量分辨的改进（续）因此，每一个衰变事例的观测量期待值为对应于每个观测量有误差估计，而且已知相互间不相关。则无参数的最小二乘问题可写为为利用一个约束条件下，改进的光子动量观测值进行 重建研究，从 的不变质量谱可以看出光子的

14、动量得到了明显的改进。 6/28/200842*第42页，共75页。检验最小二乘法的拟合优度那么2min 服从 N-m自由度的最小二乘概率密度函数分布。 据此来计算P-值 如在五个数据点的双参数拟合 也就是说，重复实验多次，有 26.3% 的值将大于 2min 。进行 1000 次蒙特卡罗实验6/28/200843*第43页，共75页。最小二乘法处理分区数据最小二乘法拟合使下式有最小值把 yi 看做泊松分布的随机变量，方差为改进的最小二乘法虽方便了计算，但对于有些区间频数太少时2min不再服从最小二乘的概率密度分布函数(或无定义)。6/28/200844*第44页，共75页。最小二乘法的归一化

15、常数问题例如 n = 400次，N = 20个区间解决的方法是从数据中直接得到 n ，或者最好是用最大似然法定n。6/28/200845*第45页，共75页。矩的一般表达式假设对随机变量 x 有n 次测量 x1,xn，服从概率密度函数分布 f (x; ) 。其中有 m 个未知参数 1,m。如果可以构造 m 个线性独立函数 ai(x), i=1,m，其均值可写为为了确定参数，上述独立函数必须进行适当选择使得含参数的函数 ei( ) 可以确定。此时，函数 ei( ) 可以通过计算无偏的样本平均值来估计因此，参数值可以通过求解 m 个 ei( ) 方程组来确定。矩的一般表达式6/28/200846*

16、第46页，共75页。线性独立函数的方差矩阵参数 1,m 估计值的协方差矩阵的无偏估计它可以与样本平均值的协方差矩阵相联系, 即6/28/200847*第47页，共75页。线性独立函数均值的方差矩阵根据线性独立函数均值和其估计值的定义，可以有6/28/200848*第48页，共75页。参数估计值的方差矩阵由于待定参数 是 e 的函数，由误差传递 因此待定参数 的协方差矩阵估计值也可以确定。而根据线性独立函数的均值估计值表达式 可知参数值的估计值 可通过求解 m 个 方程组来确定。参数估计任务完成 6/28/200849*第49页，共75页。矩方法估计参数计算出 cos 二阶代数矩的理论期待值 则

17、参数 与二阶代数矩的关系为 只要函数是可积的，采用矩方法原则上就可以测定参数。在实验 中，理论预言角分布为6/28/200850*第50页，共75页。最大似然法、最小二乘法和矩矩方法最大似然法最小二乘法数据输入单个事例单个事例直方图多维问题最容易归一化较复杂较难充分性会有信息丢失最具充分性有时与区间大小有关一致性收敛于真值收敛于真值收敛于真值有效性不是最有效通常最有效基本上与似然法一样无偏性渐进无偏渐进无偏渐进无偏拟合优度较难评估较难评估很容易充分性：估计量应包含观测值对于未知参数的全部信息；一致性：样本容量增大时，估计值收敛于真值；有效性：估计量的分布对其期望值具有最小方差；无偏性：无论样本

18、容量多大，估计值与真值无系统偏差。6/28/200851*第51页，共75页。拟合与 MINUIT 软件包52在粒子物理与核物理实验的参数估计中，利用计算机求极值的数值解越来越普遍。通用的求函数最小值程序是在使用MINUIT时，为了对结果有正确的诠释，必须对相关输出信息进行理解。MINUIT 软件包求 log L 最大值等效于求 log L 的最小值。 在 MINUIT 框架下单独使用(适于data-driven 模式)；在 PAW 环境下互动调用 MINUIT (基于Fortran)；在 ROOT 环境下互动调用 MINUIT (基于C+)。6/28/200852*第52页，共75页。ROO

19、T 平台下的 MINUIT 调用53const int npoints=2000;Double_t xnpoints;Double_t angle_cut=0.95;void 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(S

20、ET ERR, arglist ,1,ierflg);/ Set starting values and step sizes for parameters Double_t 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;

21、arglist1 = 1.; gMinuit-mnexcm(MIGRAD, arglist ,0,ierflg); gMinuit-mnexcm(HESSE, arglist ,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,nv

22、par,nparx,icstat); gMinuit-GetParameter(0, alpha, alpha_err ); gMinuit-GetParameter(1, beta, beta_err ); gMinuit-mnemat(&covmat00,npar); Double_t rho=covmat01/(alpha_err*beta_err); cout alpha = alpha +/- alpha_errendl; cout beta = beta +/- beta_errendl; cout covalphabeta= covmat01endl; cout rhoalpha

23、beta= rhoendl; cout log_l = -0.5*fmin.x minuit_fit.C6/28/200853*第53页，共75页。使用 MINUIT 应注意的问题54求极值问题时，不可能自动找到合理的初值。-log L参数值区域最小值真实最小值原因：似然函数存在多个极值。提供较好的参数不确定范围可以有助于找到真值。太大的范围会使碰巧在远离真值处找到一个区域最小值。6/28/200854*第54页，共75页。通常在MINUIT的解决方案55利用MINUIT函数 MIGRAD 找 log L 的最小值利用MINUIT函数 HESSE 计算参数的误差还可利用MINUIT函数 MIN

24、OS 做更好的误差估计6/28/200855*第55页，共75页。MINUIT函数 MIGRAD566/28/200856*第56页，共75页。函数 MIGRAD（续一）576/28/200857*第57页，共75页。函数 MIGRAD（续二）586/28/200858*第58页，共75页。MINUIT函数 HESSE596/28/200859*第59页，共75页。函数 HESSE （续一）606/28/200860*第60页，共75页。函数 HESSE （续二）616/28/200861*第61页，共75页。函数 HESSE （续三）626/28/200862*第62页，共75页。MINUI

25、T函数 MINOS63用途：更好的误差计算方法注意：有太多参数时可能会耗费非常多的CPU时间调用：在ROOT 或 PAW 使用“E” 选项6/28/200863*第63页，共75页。经常会碰到的问题：不收敛64时常会出现参数拟合过程不收敛，例如MIGRAD 不能找到最小值HESSE 得到负的二阶导数很多情况下是参数拟合函数有误（程序书写错误或理论模型不正确）数值计算中出现精度问题数值计算中出现稳定性问题HESSE 相关系数矩阵是很好的调查起点表明两个参数几乎100%关联，求最小值过程无法建立。6/28/200864*第64页，共75页。消除影响稳定性的因素：相关性65策略一：更恰当地选取拟合参

26、数例如，对于采用类似的双高斯宽度拟合事例数HESSE 相关系数矩阵双高斯宽度的大小与它们所占的份额有很强的关联6/28/200865*第65页，共75页。影响稳定性因素：相关性（续）66可以采用不同参数化过程来解决第二个高斯分布宽度与份额比的相关度从0.92减小至0.68 完全是参数化过程自身的问题策略二：固定所有高度相关的参数而只保留一个如果参数高度相关表明它们中有些是多余的，不会对拟合模型的自由度带来贡献。6/28/200866*第66页，共75页。消除影响稳定性的因素：多项式67注意：普通多项式 a0+a1x+a2x2+a3x3 的参数化过程通常会导致在拟合参数 a 之间引入很强的关联。

27、例如，导致拟合稳定性问题而不能找到高阶系数的解。解决方案：采用数学上的正交多项式，例如，勒让德多项式，第一类契贝谢夫多项式等等6/28/200867*第67页，共75页。参数拟合问题：参数的限界68有时候需要对拟合参数进行区间的限制例如，份额比参量只能在【0，1】区间，三角函数值只能在【-1，+1】之间，等等。但是，可能会导致 MINUIT 不能正确估计误差。因此，必须小心使用参数的拟合区间限制。参考的解决方案：如果区间限制的引入可使拟合稳定，最好考虑寻找有否其它可替代的无区间限制的参数化方案。如果区间限制的引入可以避免出现“非物理”但仍为合理统计解时，应考虑不要强加区间限制，并对结果采用统计意义上的解释。一般是同时给出物理解与“非物理”解。前者是使结果易于物理解释，后者是使结果可与其它实验结果进行统计意义上的正确合并。6/28/200868*第68页，共75页。举例：加速器中微子振荡实验Phys.Rev.D74:072003,20066/28/200869*第69页，共75页。如何确认拟合结果的有效性70所谓的结果有效性是指拟合结果无偏拟合误差符合统计不确定性如果对结果的有效性有疑虑，可以考虑进行P-值计算，估计结果是否为极端情况。采用蒙特卡罗样本，重复同样统计量的参数拟合过程至少100次，确认拟