参数估计与拟合.ppt

6 28 2008 BESIII暑期讲习班 1 参数估计与拟合 陈少敏清华大学工程物理系 6 28 2008 BESIII暑期讲习班 2 2020 01 16 2 主要内容 参数估计最大似然法最小二乘法矩方法MINUIT的使用 6 28 2008 BESIII暑期讲习班 3 粒子物理实验研究的目的 粒子物理实验的一个重要目的是确定粒子的属性 6 28 2008 BESIII暑期讲习班 4 参数估计 理论 数据 数据 理论 概率 微积分 给定含有参数 的理论预言分布 对数据能下何种结论 需要一个步骤来从数据D中估计参数 最常用的步骤是 拟合 统计分析 6 28 2008 BESIII暑期讲习班 5 什么是估计量 一个估计量对应于这样一个步骤 它能从实际数据测量值中对一个参数或一个分布的属性给出定量的结果 6 28 2008 BESIII暑期讲习班 6 如何评判一个参数估计量的好坏 符合程度 一致性 偏置大小 无偏性 方差大小 有效性 6 28 2008 BESIII暑期讲习班 7 参数估计与概率大小的关系 如果假设 包括 的取值 为真 可以预料会使观测结果具有高的概率 如果假设的 取值远离真值 会使观测结果具有低的概率 6 28 2008 BESIII暑期讲习班 8 似然函数 在经典统计理论里 L 并不是 的概率密度函数 根据参数好坏与概率大小的关系 可以认为真实的 应使得下式定义的似然函数 有大的数值 6 28 2008 BESIII暑期讲习班 9 最大似然估计量 取最大值 有时候L 可以有好几个极大值 注意 1 方法利用了所有信息 与如何划分数据分布区间无关 2 定义的最大似然估计量并不保证它们总是最优的 需要对诸如无偏性 有效性等问题进行研究 大样本情况下 最大似然法大都能给出了期待的好结果 小样本的情况 虽然不总是最优 但也能给出最好的实用解 6 28 2008 BESIII暑期讲习班 10 最大似然估计量的唯一性 考虑 的最大似然估计值是下列方程的解 选用等价参数h 因为 因此 h的最大似然估计值与参数选取无关 具有唯一性 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 BESIII暑期讲习班 15 估计量的方差 RCF边界法 任何估计量 不仅仅是最大似然法 的方差下界为 也称为Rao Cram r Frechet不等式 通常假设上述结论为真 利用RCF边界估计 b为偏置 最大似然估计量对大的样本统计量n几乎总是有效的 6 28 2008 BESIII暑期讲习班 16 指数型函数估计量的RCF边界 已知估计是无偏的b 0 所以由RCF边界可得 真的方差 对指数型函数求二阶导数 6 28 2008 BESIII暑期讲习班 17 RCF边界与HESSE矩阵 对于只有一个参数的情况 可以得到 求logL的最大值可通过数值计算来完成 二阶导数的矩阵 Hessian矩阵 是通过有限差值来估计 调用CERN的MINUIT软件包中的HESSE程序 6 28 2008 BESIII暑期讲习班 18 例子 估计实验所需的统计量 质子与反质子弹性散射实验 观测量为散射角x cos 服从f x 0 5 1 x 其中 是反映反质子极化的参数 目前测量值为0 10 0 02 要想在统计上将相对误差减少到5 总共需要多少个事例 由信息不等式 任何估计量的方差下界为 对于本问题 b 0 6 28 2008 BESIII暑期讲习班 19 例子 实验所需的统计量 续 由信息不等式 带入不等式可以求得n 1 2 105 6 28 2008 BESIII暑期讲习班 20 估计量的方差 图解法 也就是 6 28 2008 BESIII暑期讲习班 21 推广的最大似然法 21 如果考虑了样本大小n也是泊松分布的随机变量 平均值为 那么 6 28 2008 BESIII暑期讲习班 22 BES上的 质量测量 22 能量点i PRL69 3021 1992 6 28 2008 BESIII暑期讲习班 23 推广的最大似然法 独立 23 可以将其分解成单独求 与 估计值的问题 例如 是信号与本底分量的叠加 6 28 2008 BESIII暑期讲习班 24 推广的最大似然法 独立 续 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 分区处理数据可能出现的问题 27 分区处理的数据满足 对函数变化影响较小 因此 对直方图拟合 一定要确认区间的大小对结果无明显影响 注意 区间无穷小时 与点估计结果一致 6 28 2008 BESIII暑期讲习班 28 似然函数与贝叶斯定理 28 贝叶斯方法 对假设 采用主观概率来描述实验前 对过程的了解由 归纳 先验概率密度函数 利用贝叶斯定理 根据数据来改进先验概率密度函数 6 28 2008 BESIII暑期讲习班 29 似然法估计量与贝叶斯估计量 29 对很纯粹的贝叶斯论者 实际应用中 会碰到如何得到 的问题 对很实际的贝叶斯论者 无金科玉律 主观的 通常假定它在全区域是均匀分布 非贝叶斯方法目前在统计学中占主要地位 但是 如果被测定参数是随机变量 而且它的验前分布已知 就应该用贝叶斯方法 6 28 2008 BESIII暑期讲习班 30 最小二乘法与最大似然法 设有高斯随机变量 yi i 1 N 均值为 对应的对数似然函数 去掉与 无关的项 为 对于独立的高斯变量yi 联合概率密度函数为 6 28 2008 BESIII暑期讲习班 31 最小二乘估计量的定义 如果yi是一多维高斯变量 协方差矩阵为V 满足 那么其对数似然函数为 也就是说 我们应求下式的最小值 它的最小值定义了最小二乘法的估计量 即使yi不是高斯变量 该定义依然适用 6 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暑期讲习班 35 最小二乘估计量的方差 等效地 可以利用下式来计算 如果yi是高斯变量时 其与RCF边界一致 6 28 2008 BESIII暑期讲习班 36 最小二乘估计量的方差 续 6 28 2008 BESIII暑期讲习班 37 约束情况下的最小二乘法拟合 实际问题中会遇到测量量本身要受到某些物理定律的约束 求解可采用拉格朗日乘子法 对每一个约束引入修正因子 i 例如 能动量守恒 衰变顶点约束等等 对一个事例有m个观测量 无参数的最小二乘问题变为 6 28 2008 BESIII暑期讲习班 38 约束情况下的最小二乘法 续一 为了找到最小值 可以通过求微商方法 而n 1次迭代后 设经过n次迭代以后 找到一组解 得到函数的值 在上对 n 进行线性展开 并略去高阶项 得到 6 28 2008 BESIII暑期讲习班 39 约束情况下的最小二乘法 续二 两式联立消掉项 可以得到 因此 可以得到第n 1次迭代的l个拉格朗日乘子取值 以及第n 1次迭代的m个测量量的预期值 6 28 2008 BESIII暑期讲习班 40 约束情况下的最小二乘法 续三 当经过n 1次迭代以后 满足下式时即可终止 实验中 为了提高测量精度而采用的四动量守恒约束拟合 4 Cfit 顶点或质量约束拟合 1 Cfit 大都采用该方式来进行 此时的 2值应满足自由度为 m l 的 2分布 6 28 2008 BESIII暑期讲习班 41 例子 粒子动量分辨的改进 例如 实验观测衰变 通常情况下 探测器对光子探测的能量分辨率较差 从而影响到 0粒子动量重建的精度 已知 6 28 2008 BESIII暑期讲习班 42 例子 粒子动量分辨的改进 续 因此 每一个衰变事例的观测量期待值为 对应于每个观测量有误差估计 而且已知相互间不相关 则无参数的最小二乘问题可写为为 利用一个约束条件下 改进的光子动量观测值进行 重建研究 从 的不变质量谱可以看出光子的动量得到了明显的改进 6 28 2008 BESIII暑期讲习班 43 检验最小二乘法的拟合优度 那么 2min服从N m自由度的最小二乘概率密度函数分布 据此来计算P 值 如在五个数据点的双参数拟合 也就是说 重复实验多次 有26 3 的值将大于 2min 进行1000次蒙特卡罗实验 6 28 2008 BESIII暑期讲习班 44 最小二乘法处理分区数据 最小二乘法拟合使下式有最小值 把yi看做泊松分布的随机变量 方差为 改进的最小二乘法虽方便了计算 但对于有些区间频数太少时 2min不再服从最小二乘的概率密度分布函数 或无定义 6 28 2008 BESIII暑期讲习班 45 最小二乘法的归一化常数问题 例如n 400次 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 BESIII暑期讲习班 47 线性独立函数的方差矩阵 参数 1 m估计值的协方差矩阵的无偏估计 它可以与样本平均值的协方差矩阵相联系 即 6 28 2008 BESIII暑期讲习班 48 线性独立函数均值的方差矩阵 根据线性独立函数均值和其估计值的定义 可以有 6 28 2008 BESIII暑期讲习班 49 参数估计值的方差矩阵 由于待定参数 是e的函数 由误差传递 因此待定参数 的协方差矩阵估计值也可以确定 而根据线性独立函数的均值估计值表达式 可知参数值的估计值可通过求解m个方程组来确定 参数估计任务完成 6 28 2008 BESIII暑期讲习班 50 矩方法估计参数 计算出cos 二阶代数矩的理论期待值 则参数 与二阶代数矩的关系为 只要函数是可积的 采用矩方法原则上就可以测定参数 在实验中 理论预言角分布为 6 28 2008 BESIII暑期讲习班 51 最大似然法 最小二乘法和矩 充分性 估计量应包含观测值对于未知参数的全部信息 一致性 样本容量增大时 估计值收敛于真值 有效性 估计量的分布对其期望值具有最小方差 无偏性 无论样本容量多大 估计值与真值无系统偏差 6 28 2008 BESIII暑期讲习班 52 拟合与MINUIT软件包 52 在粒子物理与核物理实验的参数估计中 利用计算机求极值的数值解越来越普遍 通用的求函数最小值程序是 在使用MINUIT时 为了对结果有正确的诠释 必须对相关输出信息进行理解 MINUIT软件包 求logL最大值等效于求 logL的最小值 在MINUIT框架下单独使用 适于data driven模式 在PAW环境下互动调用MINUIT 基于Fortran 在ROOT环境下互动调用MINUIT 基于C 6 28 2008 BESIII暑期讲习班 53 ROOT平台下的MINUIT调用 53 constintnpoints 2000 Double tx npoints Double tangle cut 0 95 voidminuit fit Getdatapointsget input data Prepareforfitconstintnpar 2 TMinuit gMinuit newTMinuit npar gMinuit SetFCN fcn Double targlist 10 Int tierflg 0 arglist 0 1 gMinuit mnexcm SETERR arglist 1 ierflg SetstartingvaluesandstepsizesforparametersDouble tvstart npar 0 5 0 5 Double tstep npar 0 1 0 1 gMinuit mnparm 0 alpha vstart 0 step 0 0 0 ierflg gMinuit mnparm 1 beta vstart 1 step 1 0 0 ierflg Nowreadyforminimizationsteparglist 0 500 arglist 1 1 gMinuit mnexcm MIGRAD arglist 0 ierflg gMinuit mnexcm HESSE arglist 0 ierflg gMinuit mnexcm MINOS arglist 0 ierflg PrintresultsDouble tfmin fedm errdef covmat npar npar Double talpha alpha err beta beta err Int tnvpar nparx icstat gMinuit mnstat fmin fedm errdef nvpar nparx icstat gMinuit GetParameter 0 alpha alpha err gMinuit GetParameter 1 beta beta err gMinuit mnemat root xminuit fit C 6 28 2008 BESIII暑期讲习班 54 使用MINUIT应注意的问题 54 求极值问题时 不可能自动找到合理的初值 原因 似然函数存在多个极值 提供较好的参数不确定范围可以有助于找到真值 太大的范围会使碰巧在远离真值处找到一个区域最小值 6 28 2008 BESIII暑期讲习班 55 通常在MINUIT的解决方案 55 利用MINUIT函数MIGRAD找 logL的最小值利用MINUIT函数HESSE计算参数的误差还可利用MINUIT函数MINOS做更好的误差估计 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 函数HESSE 续三 62 6 28 2008 BESIII暑期讲习班 63 MINUIT函数MINOS 63 用途 更好的误差计算方法注意 有太多参数时可能会耗费非常多的CPU时间调用 在ROOT或PAW使用 E 选项 6 28 2008 BESIII暑期讲习班 64 经常会碰到的问题 不收敛 64 时常会出现参数拟合过程不收敛 例如MIGRAD不能找到最小值HESSE得到负的二阶导数很多情况下是参数拟合函数有误 程序书写错误或理论模型不正确 数值计算中出现精度问题数值计算中出现稳定性问题HESSE相关系数矩阵是很好的调查起点 表明两个参数几乎100 关联 求最小值过程无法建立 6 28 2008 BESIII暑期讲习班 65 消除影响稳定性的因素 相关性 65 策略一 更恰当地选取拟合参数例如 对于采用类似的双高斯宽度拟合 事例数 HESSE相关系数矩阵 双高斯宽度的大小与它们所占的份额有很强的关联 6 28 2008 BESIII暑期讲习班 66 影响稳定性因素 相关性 续 66 可以采用不同参数化过程来解决 第二个高斯分布宽度与份额比的相关度从0 92减小至0 68完全是参数化过程自身的问题 策略二 固定所有高度相关的参数而只保留一个如果参数高度相关表明它们中有些是多余的 不会对拟合模型的自由度带来贡献 6 28 2008 BESIII暑期讲习班 67 消除影响稳定性的因素 多项式 67 注意 普通多项式a0 a1x a2x2 a3x3的参数化过程通常会导致在拟合参数a之间引入很强的关联 例如 导致拟合稳定性问题而不能找到高阶系数的解 解决方案 采用数学上的正交多项式 例如 勒让德多项式 第一类契贝谢夫多项式等等 6 28 2008 BESIII暑期讲习班 68 参数拟合问题 参数的限界 68 有时候需要对拟合参数进行区间的限制例如 份额比参量只能在 0 1 区间 三角函数值只能在 1 1 之间 等等 但是 可能会导致MINUIT不能正确估计误差 因此 必须小心使用参数的拟合区间限制 参考的解决方案 如果区间限制的引入可使拟合稳定 最好考虑寻找有否其它可替代的无区间限制的参数化方案 如果区间限制的引入可以避免出现 非物理 但仍为合理统计解时 应考虑不要强加区间限制 并对结果采用统计意义上的解释 一般是同时给出物理解与 非物理 解 前者是使结果易于物理解释 后者是使结果可与其它实验结果进行统计意义上的正确合并 6 28 2008 BESIII暑期讲习班 69 举例 加速器中微子振荡实验 Phys Rev D74 072003 2006 6 28 2008 BESIII暑期讲习班 70 如何确认拟合结果的有效性 70 所谓的结果有效性是指 拟合结果无偏拟合误差符合统计不确定性 如果对结果的有效性有疑虑 可以考虑 进行P 值计算 估计结果是否为极端情况 采用