粒子物理与核物理实验中的数据分析lecture-9-最小二乘法

粒子物理与核物理实验中的数据分析陈少敏清华大学第九讲：昀小二乘法2本讲要点„昀小二乘法与昀大似然法的关系„线性情况下的昀小二乘估计„非线性情况下的昀小二乘估计„约束情况下的昀小二乘法„检验昀小二乘法的拟合优度„应用昀小二乘法处理分区数据„不等精度关联实验结果的并合问题3昀小二乘法与昀大似然法假设有高斯随机变量:yi,i=1,…,N,其平均值为对于独立的高斯变量yi，联合概率密度函数为[](;)iiiEyxλλθ==21,...,[]NiixxVyσ=这里，与已知。∏=⎟⎟⎠⎞⎜⎜⎝⎛−−=Niiiiiyyg122222)(exp21),;(σλπσσλ对应的对数似然函数(去掉与θ无关的项)为221[(;)]1log()2NiiiiyxLλθθσ=−=−∑2221[(;)]()Niiiiyxλθχθσ=−=∑2log()Lθχ求的昀大值等效于求的昀小值。θ为了估计参数，可以用曲线拟合所有的测量点(见右图)。4昀小二乘估计量的定义如果yi是一多维高斯变量，协方差矩阵为V，满足()()⎥⎦⎤⎢⎣⎡−−−=−λλπλyVyVVygTN12/12/21exp||)2(1),;(那么其对数似然函数为111log()[(;)]()[(;)]2NiiijjjiLyxVyxθλθλθ−==−−−∑也就是说，我们应求下式的昀小值21,1()[(;)]()[(;)]NiiijjjijyxVyxχθλθλθ−==−−∑它的昀小值定义了昀小二乘法的估计量θ，即使yi不是高斯变量，该定义依然适用。(实际上，yi通常是高斯的，因为中心极限定理会导出测量误差也是高斯的。)5线性昀小二乘法估计这里aj(x)是x的任意线性独立函数。(;)xλθθ如果是的线性函数，昀小二乘法有一些简单的特殊性质，∑==mjjjxax1)();(θθλˆθ没有偏置，而且得到的方差昀小(高斯-马尔可夫定理)。用矩阵来表示时，令Aij=aj(xi)，有)()()()()(112θθλλθχAyVAyyVyTT−−=−−=−−对θi求偏微分，并令结果等于零，有0)(2112=−−=∇−−θχAVAyVATT解方程得到昀小二乘法的估计量111ˆ()TTAVAAVyByθ−−−=≡ˆiiyθ估计量是测量量的线性函数。6非线性昀小二乘法估计(;)xλθθθ如果是的非线性函数，昀小二乘法不能给出参数的严格解，需要通过迭代法求的近似解,使得下式昀小21()()()TyVyχθλλ−=−−如果采用牛顿法求上式的昀小值，第n+1次迭代公式可采用()()(1)()1()()222ˆˆˆˆˆˆ()()nnnnnniijiijGggGθθθθθθθθχχθθθ+−===−∂∂==∂∂∂i，1/2(1)()(1)()21ˆˆˆˆ()nnnniiiθθθθεε++=⎡⎤−=−⎢⎥⎣⎦∑迭代终止判据：是一个给定的小数。(1)ˆˆnθθ+=7昀小二乘估计量的方差ˆˆcov[,]ijijUθθ=协方差矩阵元在线性情况下的误差传递可以写为11)(−−==AVABVBUTT等效地，可以利用下式来计算2211,ˆ1()()()()2NijikkljlklijUaxVaxθθχθθ−−=⎡⎤∂==⎢⎥∂∂⎢⎥⎣⎦∑如果yi是高斯变量时，其与RCF边界一致。111ˆ()TTAVAAVyByθ−−−=≡8昀小二乘估计量的方差(续)2(;)()xλθχθ对于是参数的线性函数情况下，是二次型的)ˆ)(ˆ(21)ˆ()(ˆ1,2222jjiimjijiθθθθθθχθχθχθθ−−⎥⎥⎦⎤⎢⎢⎣⎡∂∂∂+===∑方差从正切平面到椭圆11)ˆ()(2min22+=+=χθχθχ(;)xλθθ如果不是的线性函数，上述式子会有偏差，但仍不失为较好的近似。22min()1χθχθ≤+可以把区间看作“置信区间”，给出了包含真值的可能性。cov[,]iijijyVyy=上式并不依赖于是否为高斯变量，但无论何种情况，都要计算协方差矩阵注意：。9多项式的昀小二乘法拟合用一个多项式来拟合右图∑==mjjjmxx00),...,;(θθθλ第0阶(一个参数)第1阶(两个参数)第4阶(五个参数)¾对于单参数拟合情况(例如上图的横线):5.4513.066.2ˆ2min0=±=χθ0022ˆˆ0minˆ()1θθσχθσχ±=+标准误差是根据来确定的。例如:10多项式的昀小二乘法拟合(续)¾对于双参数拟合的情形(有非零斜率的直线)椭圆的倾角给出相关系数(与昀大似然法相同)。01012ˆˆ0.930.30;0.680.10;ˆˆˆcov[,]0.028;0.90;3.99.rθθθθρχ=±=±=−==−=01ˆˆ,θθσσ切线给出。010101ˆˆˆˆˆˆ(,)(,)cov(,,)0θθηηηη→=虽然表面上可以通过变换使得从而得到不关联的估计值，但是对新参数的理解也许并不容易。¾对于五个参量拟合的情形(有非零斜率的直线)χ2min值的大小反映了数据与假设之间的符合程度。可以用来检验拟合优度。9曲线通过所有点；9χ2min=0，参数的数目=数据点的数目。11约束情况下的昀小二乘法在处理实际问题中，常常会遇到测量量本身要受到某些物理定律的约束。21()(')(')1,...,TixxxVxxxillχψ−=−−==昀小()=0,(共个约束条件)求解可以采用拉格朗日乘子法，对每一个约束引入修正因子αi，有21(,)(')(')2TTxxxVxxxχααψ−=−−+=()昀小例如，能动量守恒，衰变顶点约束等等。无参数的昀小二乘问题变为1''1(,...,)'(,...,)mmxxxxxx====真值观测值1(1)()()()()()(1)()(1)()()()()(')'()nnTnnnnxxxnTnnxxFxVFxxFxxxxxVFxFlmxαψψα−+++⎡⎤⎡⎤=+−⎣⎦⎣⎦∂=−==×∂，其中矩阵可以证明，在经过n次迭代以后，第n+1次的值为12约束情况下的昀小二乘法(续)当经过n+1次迭代以后，满足下式时即可终止2(1)(1)2()()12(1)(1)(1)212(,)(,)(,)(),nnnnnnnxxxfxχαχαεχαεεε+++++−这里与均为小量。在粒子与核物理实验中，为了提高测量精度而采用的四动量拟合(4-Cfit)或顶点拟合(1-Cfit)，大都采用该方式来进行。22min()lχχ∼例如，实验观测衰变00,ρπππγγ±±→→对光子的能量测量通常较差，从而影响到ρ±的重建。可以利用π0质量约束，进行1-C拟合，改善对光子的测量精度。ρ±改善前的信号不变质量谱ρ±改善后的信号不变质量谱13检验昀小二乘法的拟合优度那么χ2min服从在N-m自由度下的昀小二乘概率密度函数分布。据此来计算P-值,1,...,()(;)iijiyiNVxλθθ=如果是高斯量已知，而且假设是的线性函数，函数形式正确。∫∞=2min);(χdznzfPd例如考虑在前面双参数拟合的情况263.03,99.32min=→=−=PmNχ也就是说，重复实验多次，有26.3%的值将大于χ2min。进行1000次蒙特卡罗实验而对于水平线拟合，有92min101.34,5.45−×=→=−=PmNχ14拟合优度与误差的昀小值小的统计误差并不意味着是一个好的拟合(反之亦然)¾χ2曲线在其昀小值附近变化给出统计误差；¾χ2min的曲率大小给出拟合的优度。在水平线拟合中，可以人为改变数据点纵向的位置，但保持误差不变，使得48.413.084.2ˆ2min0=±=χθ方差与改变前一样，但χ2min变“好”了。χ2(θ0)曲线只是向下平移，表明与数据符合更好。但曲线形状并没有发生变化，即误差并没有改变。5.4513.066.2ˆ2min0=±=χθ15拟合优度与误差的昀小值(续)¾估计量的方差告诉我们：¾P-值告诉我们：•如果实验从复多次，估计量θ分布有多宽。•但是，它并不告诉我们假设是否正确。•P-值太低，则假设可能有误，即存在系统误差。•如果假设正确，并且实验重复了多次，实验与假设按照统计的χ2min完全符合或符合得更差的比率是多少。16昀小二乘法处理分区数据昀小二乘法拟合使得下式有昀小值(;)Nnfxθ在有个区间，填入次的直方图中，假设概率密度函数为，有maxmin()(;)()iiixiiixyinfxdxnpλθθθ===∫第个区间的频数∑=−=Niiiiy1222))(()(σθλθχ2[]iiVyσ=为先验未知量。把yi看做泊松分布的随机变量，方差为2()()iiσλθ=昀小二乘法改进的昀小二乘法虽方便了计算，但对于有些区间频数太少时χ2min不再服从昀小二乘的概率密度分布函数(或无定义)。2()iiyσ=改进后的昀小二乘法17昀小二乘法的归一化问题例如n=400次，N=20个区间∫==maxmin)();(),(iixxiipdxxfθνθννθλνθ尽量避免拟合规一化常数。例如，引入可调参数，并与一起拟合。ˆnν对而言是一个不好的估计量。2min2minˆ2ˆLSMLSnnχννχ=+=−解决的方法是从数据中直接得到n，或者昀好是用昀大似然法定n。解决的方法是从数据中直接得到n，或者昀好是用昀大似然法定n。ˆNnν会出现当大小时，的相对误差会变大。18在PAW上的昀小二乘拟合例如:PAWve/crpar(3)r10001;h/fit100g!3par(拟合高斯分布)19在ROOT上的昀小二乘拟合在图形显示屏幕中右击鼠标键，可以出现拟合菜单，选择所需的函数对直方图进行拟合。或者采取如下方式：roothist-Fit(“gaus”);roothist-Fit(“landau”);roothist-Fit(“exp0”);roothist-Fit(“pol1”);…详见root用户手册20用昀小二乘法并合各实验结果用昀小二乘法得到λ的N个测量的权重平均值cov[,],,ijijyyV=如果各测量量之间相关通过求下式的昀小值∑=−−−=NjiiijiyVy1,12)())(()(λλλχ∑∑∑∑==−=−====NjijijiNlkklNjijiNiiiwVwVVVwyw1,1,1111]ˆ[)()(,ˆλλ2,1,...,[],iiiyiiNVyσλ====第个测量结果假设已知真值在各测量量不相关的情况下，并合方法与昀大似然法一样。ˆ,λ昀小二乘法得到的是无偏的而且方差昀小(高斯-马科夫定理)。21未知关联矩阵的并合问题2()1Nχλ−的情形定义标度因子2/(1)fNχ=−并合后的误差为ˆˆfλλσσ⇒2()1Nχλ−的情形粒子数据组PDG建议的处理方法。结果比较保守。对应于多次测量之间存在正关联的情况，得到的误差比真实结果要小。可以考虑建立等效协方差矩阵(SchmellingM.Phys.Scripta,1995(51):676)2,,,,1,...,iiiijijVVfijijNσσσ==≠=解方程21,1ˆˆ()()()()1NiijiijfyVyNχλλ−==−−=−∑fˆ[]Vλ新的22对两个相关实验求平均的例子假设有两相关测量量y1,y2，且21122122Vσρσσρσσσ⎛⎞=⎜⎟⎝⎠σσρσσσσσρλσρσσσσρσσλ=−+−=−+−=−+=212221222122122212122212)1(]ˆ[2,)1(ˆVwywwy如果因第二个测量导致方差倒数的增加为011111221212⎟⎟⎠⎞⎜⎜⎝⎛−−=−σσρρσσ12/0wρσσ→如果，那么加权平均的结果将不在y1和y2之间。如果相关是由于使用相同数据的话，上述情况不可能发生，但有却可能来自共同的随机效应；此时的平均值很不可信，例如ρ，σ1，σ2不正确。如果相关是由于使用相同数据的话，上述情况不可能发生，但有却可能来自共同的随机效应；此时的平均值很不可信，例如ρ，σ1，σ2不正确。第二个测量结果对平均值有帮助。23小结1.与昀大似然法的联系2.线性的昀小二乘法估计3.非线性的昀小二乘法估计4.用昀小二乘法检验拟合优度5.用昀小二乘法处理分区数据6.不等