粒子物理与核物理实验中的数据分析lecture-12-开拆法

粒子物理与核物理实验中的数据分析陈少敏清华大学第十二讲：开拆法2本讲要点„数学公式,反应函数(矩阵)„求反应矩阵的逆„修正因子„正规化的开拆法„估计量的方差与偏置„正规化参数的选择„举例a)Tikhonov规则b)MaxEnt规则3图像还原问题⊗＝一个常见的问题：由于实验仪器的原因而出现图像变形，例如如果，已知通过探测器模拟可以给出其影响的形式实验观测分布能否还原出不受实验仪器影响的分布？能否还原出不受实验仪器影响的分布？Unfolding(开拆法)真实分布4开拆问题的表述考虑有随机变量，我们的目标是要找到概率密度函数如果函数可参数化为　　　　，那么确定概率密度函数，等效于若无参数化形式，可通过构造直方图y)(yf);(θGyf最大似然法θˆGbin()1,...,jjjtotjpfydyjMpμμ===∫“真实的直方图”目标：为(或)构造估计量jμjp问题：在测量时不可能没有误差y参数的数目=区间的数目M()fy各区间之间填入的数目互串,导致散开变宽。5反应矩阵测量误差的影响:数据：=y真值；观测值=x∫=dyyfyxRxftruemeas)()|()(写成离散形式NiRMjjiji,...,1,1==∑=μν观测直方图(期待值)反应矩阵真实直方图(|)ijRPij=观测值在第区真实值在第区][),...,(1iiNnEnnn==νG这里注意：是常数，会受到统计涨落的影响。νμGG,nG6效率，本底有时侯，事例可能会不被探测到是在观测直方图上预期的本底数目，并假设它是已知的。有时在无真实事例发生的时候，也有事例被观测到11(|)(|()NNijiijRPijPjε=====∑∑观测值在第区真实值在第区观测值在全范围真实值在第区)效率∑=+=MjijijiR1βμνiβ取决于在第j区的真实直方图7各关键量总汇“真实”直方图:观测直方图的期待值:反应矩阵:效率:∑===MjjtotM11),,...,(μμμμμ概率:观测直方图:GtotMpppμμ/),...,(1GG==),...,(1Nννν=G),...,(1Nnnn=G(|)ijRPij=观测值在第区真实值在第区∑==NiijjR1ε),...,(1Nβββ=预期的本底:GβμνGGGG+==RnE][为了找到的估计量，需要相关的概率理论，例如:iiennPiniiiννν−=!);(泊松分布，或关联矩阵以便构造],cov[jiijnnV=2logχL或μGM个区间N个区间8为什么要用开拆法一般而言，我们并不需要开拆法，例如当比较现有理论的预期值时，最好是将探测器相应叠加到理论中去。即在预期值中包含探测器效应并与未修正的原始数据相比较。但是，不将实验数据进行开拆处理，结果发表后，有关反应矩阵的知识将不在保留。而且，开拆后的分布可以直接与各种理论的预言相比较，也可以与别的实验经过开拆以后的分布相比较。nG通常开拆的结果更为有用，因为当反应矩阵变得不可恢复时，即使对实验结果可能又有了新的理论解释，也很难进行理论检验。在粒子物理研究中，开拆法常用的领域为：•结构函数•τ的谱函数(也就是强子不变质量谱)•强子事例形状分布•粒子多重数分布•．．．9反应矩阵的逆假设的逆存在:若数据是泊松分布βμ最大似然法的估计量为νGGG+=R)(1βνμGGG−=−RiiennPiniiiννν−=!);(∑=−=NiiiinL1)log()(logννμG)(ˆˆ1βμνGGGGG−==−nRn则有μνnˆμ若Ｒ的非对角元太大，即区间宽度比分辨率要小时，会导致上式有很大的方差，以及在相邻区间产生很强的负关联。若Ｒ的非对角元太大，即区间宽度比分辨率要小时，会导致上式有很大的方差，以及在相邻区间产生很强的负关联。？10错误的原因假设真的有精细结构应用R给出观测期待值时，虽然一些结构还能留下，但大部分的精细结构都被抹平了。应用R-1到恢复精细结构：采用观测值时，由于统计涨落的缘故，有不少非物理因素造成的突起。μGμGμνGGR=νGνμGG1−=R但我们没有只有νGnGnGR-1“认为”这是与原来的精细结构有关，导致有振荡效应。1ˆRnμ−=GG11重新研究最大似然法的解是无偏的!计算估计量的方差μβμGGGG=−=−)][(]ˆ[1nERE1111,11ˆˆcov[,]()()cov[,]()()NNijijikjlklikjkkklkURRnnRRμμν−−−−=====∑∑利用RCF边界做无偏估计量∑=−=⎥⎦⎤⎢⎣⎡∂∂∂−=NiiiliklkklRRLEU121log)(νμμ倒数后给出∑=−−=NikjkikijRRU111)()(ν即使最大似然法在各无偏估计中给出的方差最小。但得到的方差可能仍然很大。为了减小方差，必须引入一些偏置量策略：接受小的偏置量(系统误差)以换取大幅减小方差(统计误差)。策略：接受小的偏置量(系统误差)以换取大幅减小方差(统计误差)。,]iklklknnnδν=是独立的泊松变量时，cov[12简单方法：修正因子法对做相同的分区，并取,与是来自无本底情况下的蒙特卡罗模拟结果。通常，因此方差不会被放大。但偏置为νμGG,)(ˆiiiinCβμ−=()MCiiMCiCμν=修正因子MCiνMCiμ2ˆˆcov[,]cov[,]ijijiijUCnnμμ==)1(OCi≈iiiEbμμ−=]ˆ[,MCsigsigiiiiiiiMCsigiibμμνννβνν⎛⎞=−=−⎜⎟⎝⎠这里。注意：该偏置量存在把拉向的倾向，造成模型检验的困难。μˆGMCμG1）如果分区宽度≥几倍的分辨率，结果不会太坏。2）实际应用中，该方法常用于事例形状变量的分布研究中。1）如果分区宽度≥几倍的分辨率，结果不会太坏。2）实际应用中，该方法常用于事例形状变量的分布研究中。除非模拟采用的模型无误，否则上式不为零，需要考虑对应的系统误差。13例子：脉冲形状的还原÷＝¾将理论(真实)的直方图除以受实验仪器影响的直方图得到修正因子¾将观测直方图乘以修正因子直方图得到理论(真实)的直方图＝⊗14正规化的开拆法考虑“合理的”估计量，使得某些　　　　满足LlogΔLLLloglog)(logmaxΔ−≥μG估计量可通过将下式求最大值，选出最“光滑的”一个来构造)()(log)(μμαμGGGSL+=Φ=)(μGS正则化函数(光滑的量度)=α正则化参量(选择给出欲求的)LlogΔ另外，要求开拆后对总事例数的估计为无偏的∑∑===NijitotijijinR1,μν⎥⎦⎤⎢⎣⎡−++=∑=NiitotnSL1)()(log),(νλμμαλμϕGGGλ:拉格朗日乘子在约束情况下将下式求最大值,Rνμβμ=+GGG因所以是的函数logLnνΔGG描述了数据与期待值之间的“距离”。μG15正规化的开拆法(续)totNiin=→=∂∂∑=10/νλϕ∞→=αα0给出最光滑的解(数据无关)给出最大似然解(方差可能太大)()SμG显然，需要正规化函数与如何取α值的方案。所得到的估计量的好坏由它们的偏置和方差来判断。a)Tikhonov规则b)MaxEnt规则16Tikhonov规则取光滑度等于第k阶导数均值的平方，有[]...2,1,)()(2=⎟⎟⎠⎞⎜⎜⎝⎛−=∫kdydyyfdyfSktruektrue这里通常取k=2，使得S约等于曲率平方的平均值。对直方图而言，也就是∑−=++−+−−=21221)2()(MiiiiSμμμμG注意：2阶导数对直方图的第一和最后的区间没有很好的定义。如果在下，采用Tikhonov(k=2)规则，221logχ−=L2(,)()()2iSαϕμλχμμμ=−+GGG是的二次项令的导数为零，给出线性方程。ϕ在高能物理界现有好几个现成的程序：RUN，Blobel,SVD,Höcker,…在高能物理界现有好几个现成的程序：RUN，Blobel,SVD,Höcker,…Sov.Math.5(1963)103517最大熵(MaxEnt)规则另一种表征光滑度的方法基于熵。对于一组概率而言，它表示为∑=−=MiiippH1log所有相等意味着熵最大(最光滑)ip有一个，其它为零，则意味着熵最小1=ip用熵作为正规化函数，1()()loglog(MMiiitottottotSHμμμμμμμ===−∝∑GG填入个区间中各种可能的总数)有时侯，根据贝叶斯统计的先验概率密度函数(?)μμGG→)(S这里，我们坚持采用经典近似:估计量的好坏由偏置，方差来判断。注意：熵并不取决于区间的顺序。注意：熵并不取决于区间的顺序。Ann.Rev.Astron.Astrophys.24(1986)12718的方差与偏置μˆG一般来说，决定的方程是非线性的。在附近展开)(ˆnGGμobsnG)(ˆnGGμ1222222ˆˆ()(),,,1,...,,11,...,,1,0,1,1,,1,...,,1,...,,1,1,1,...,.obsobsijijiijijjnABnnijMAiMjiMjMiMjNnBiMjNnμμϕμμϕμλϕλϕμϕλ−≈−−⎧∂=⎪∂∂⎪⎪∂⎪M==−==+⎨∂∂⎪⎪∂⎪==+=+∂⎪⎩⎧∂==⎪∂∂⎪=⎨∂⎪==+=⎪∂∂⎩GGGGG，G.Cowan,StatisticalDataAnalysis,OxfordUniversityPress(1998)ϕ为非正规的似然函数19的方差与偏置(续)μˆG利用误差传递得到协方差ˆˆcov[,],ijijUμμ=,1BACCVCUT−==这里以及对偏置的估计量，,]ˆ[iiiEbμμ−=∑∑==−∂∂=−=NjNjjjjijjijinnnCb11),ˆ(ˆ)ˆ(ˆνμν此处而且通常情况下.ˆˆβμνGGG+=Rˆnν≠GG20正规化参数α的选取α决定了置于数据的权重大小以便能与光滑度相比较，α=0给出最大的光滑估计值，并与数据无关。因此虽然方差为零，但有明显的偏置。而取大的α，则回到高度振荡无偏的最大似然解。为了在偏置与方差之间达到最大平衡:选择α使均值误差的平方最小()2211ˆ11ˆ,.ˆMMiiiiiiiiiUbMSEUbWeightedMSEMMμ==+=+=∑∑或或要求偏置不大于它自身的估计方差。它可以找到α的值使得221ˆˆˆˆcov[,].ˆMibijijiiibMWbbWχ====∑这里ˆiiWG.Cowan,StatisticalDataAnalysis,OxfordUniversityPress(1998)M。Schmelling,NIMA340(1994)40021例子:Tikhonov规则(k=2)22例子:最大熵(MaxEnt)规则23一个在图像处理中的最大熵例子最大熵值方法常用于天文观测图像的重建,与点源的偏置较小,易于推广到两维以上的情况。最大熵值方法常用于天文观测图像的重建,与点源的偏置较小,易于推广到两维以上的情况。24例子：τ的谱函数()()XudXXusτττντν→→为了测定奇异夸克质量，实验上可采用比较与中,的质量平方差0022()()MudMusΔ=−—=smEur.Phys.J.C11:599,1999Tikhonov规则修正因子方法由于探测器对两者影响各不相同，因此，需要用开拆法求出“真实”分布。25小结1.数学上的原理2.求反应矩阵的逆3.修正因子4.正规的开拆过程5.估计量的方差与偏置6.正规化参数的选择7.例子,nRμνμνβμ=+GGGGGGG真实直方图，数据以及其期待值，满足目标是构造的估计量。有很大的振荡行为(及大的方差)，但在各种无偏解中具有零偏置与最小的方差。/MCMCiiiCμν=,方法又快又简单。Tikhonov:kMaxEnt:logiiiHpp=−∑从第阶导数的均方值中进行光滑处理从熵中进行光滑处理在求解过程中采用了线性近似，因而不是无偏的2M()χ=无最好的方案，可以采用区间总数的方案。只要探测器的响应可知，就一定可以得到真实的分布