1、粒子物理与核物理实验中的数据分析,杨振伟清华大学第十一讲：解谱法（开拆法）,2,本讲要点,数学公式,反应函数(矩阵)求反应矩阵的逆修正因子正规化的解谱法估计量的方差与偏置正规化参数的选择举例,Tikhonov 规则MaxEnt 规则,3,图像还原问题,一个常见的问题：由于实验仪器的原因而出现图像变形，例如,如果已知探测器响应（可通过探测器模拟得到探测器响应的形式），,探测器响应,实验观测分布,能否还原出不受实验仪器影响的真实分布？,真实分布,4,解谱问题的表述(1),考虑随机变量y，目标：寻找概率密度函数f(y),5,解谱问题的表述(2),问题：y 的测量不可能没有误差,1）第i个区间的y在第
2、j个区间测量到 （测量分辨率）2）第i个区间的某些事例没有测量到 （测量效率）3）某些区间的事例完全测量不到 （接收度）,后果：f(y) 模糊化，峰展宽， 甚至分布形状完全变形。,后果严重时需要考虑解谱法还原真实分布,6,响应矩阵,测量误差的影响可用积分方程表示:,响应矩阵 Rij=P(观测值在第 i 区|真实值在第 j 区),7,响应矩阵,Rij=P(观测值在第 i 区|真实值在第 j 区),注意： 是常数， 会受到统计涨落的影响。,真实直方图离散化的p.d.f.,观测数据和数据的期待值,8,效率、本底,有时，事例可能会不被探测到: 效率,真实直方图第 j 区探测效率,9,各关键量汇总,“真
3、实”直方图:,概率:,观测直方图的期待值:,观测直方图:,响应矩阵:,效率:,预期的本底:,M 个区间,N 个区间,10,为什么要用解谱法？,但是，不将实验数据进行解谱处理，结果发表后，有关反应矩阵的知识将不在保留。而且，解谱后的分布可以直接与各种理论的预言相比较，也可以与别的实验经过解谱以后的分布相比较。,一般而言，我们并不需要解谱法，例如当比较现有理论的预期值时，最好是将探测器相应叠加到理论中去。即在预期值中包含探测器效应并与未修正的原始数据 相比较。,通常解谱后的结果更为有用，否则当反应矩阵不可恢复时，即使对结果又有新的理论解释，也很难进行理论检验。,在粒子物理研究中，解谱法常用的领域为
4、：,强子结构函数 的谱函数(也就是强子不变质量谱) 强子事例形状分布 粒子多重数分布,11,为什么用解谱法(举例),中微子与铁原子核相互作用，产生电子，测量电子的能谱。实验上只可能测量电子在铁中的射程。,能量与平均射程的关系,因子修正法：能谱(圆圈)与真值(直方图),解谱法,12,响应矩阵的逆,若数据是泊松分布,最大似然法的估计量为,假设 的逆存在:,则有,若的非对角元太大，即区间宽度比分辨率要小时，会导致上式有很大的方差，以及在相邻区间产生很强的负关联。,？,13,错误的原因(I),考虑一个简单的例子,解决办法是进行平滑处理，消除无意义的统计涨落。但平滑会带来偏向性，需要在涨落与偏向性之间找
5、到平衡。,14,错误的原因(II),考虑周期分布函数,实际上，非物理的剧烈涨落主要是有高频部分造成的。平滑处理，主要是消除或压低高频部分无意义的涨落。,15,错误的原因(III),假设 真的有精细结构,大部分精细结构被抹平，但会隐藏精细结构的信息。,R-1 作用回 应完全恢复精细结构：,16,重新研究最大似然法的解(1),无偏!,计算估计量的方差,估计量的平均值,17,重新研究最大似然法的解(2),利用 RCF 边界做无偏估计量,倒数后得到,为了减小方差，必须引入一些偏置量,策略：接受小的偏置量(系统误差)以换取大幅 减小方差(统计误差)。,与刚才得到的结果一致。,18,简单方法：修正因子法,
6、对 做相同的分区，并取,与 用蒙特卡罗模拟得到（无本底）。,通常 ，因此方差不会被放大。,(相当于R-1取为对角矩阵) (R-1)ii=Ci,19,修正因子法中的偏置问题,注意：该偏置量存在把 拉向 的倾向，造 成模型检验的困难。,1）如果分区宽度 几倍的分辨率，结果不会太坏。 2）实际应用中，中。该方法常用于事例形状变量的分布研究,除非模拟采用的模型无误，使得 ，否则上式不为零，需要考虑对应的系统误差。,修正因子法的偏置,20,例子：脉冲形状的还原,将理论(真实)的直方图除以受实验仪器影响 的直方图得到修正因子,将观测直方图乘以修正因子直方图得到理论 (真实)的直方图,21,正规化的解谱法,
7、考虑“合理的”估计量，对选定的 满足,正则化函数(光滑性的量度),正则化参量(其选择与给定的 对应),估计量满足该不等式并且最光滑，等价于将下式求最大值,22,正规化的解谱法（续一）,另外，要求解谱后对总事例数的估计为无偏的,:拉格朗日乘子,在约束情况下将下式求最大值,显然，,23,正规化的解谱法（续二）,需要正规化函数 以及选取 值的方案。,所得到的估计量的好坏由它们的偏置和方差来判断。,Tikhonov 规则MaxEnt 规则,24,Tikhonov 规则,取光滑度等于第 k 阶导数平方的均值，有,通常取 k=2，使得 S 约等于曲率平方的平均值。对直方图而言，也就是,注意：2 阶导数对直
8、方图的第一和最后的区 间没有很好的定义。,Sov. Math.5(1963)1035,25,Tikhonov 规则（续）,如果在 下，采用Tikhonov(k=2)规则,是 i 的二次项。对 求偏微分，给出线性方程，得到 估计值与方差。,在高能物理界现有好几个现成的程序：RUN，Blobel, SVD, Hcker,26,最大熵(MaxEnt)规则,另一种表征光滑度的方法基于熵。对于一组概率而言，熵定义为,所有 相等意味着熵最大(最光滑),有一个 ，其它为零，则意味着熵最小。,Ann. Rev. Astron. Astrophys.24 (1986)127,27,最大熵(MaxEnt)规则（续
9、）,用熵作为正规化函数，,有时侯，根据贝叶斯统计,这里，我们仍然采用经典近似: 估计量的好坏由偏置，方差来判断。,注意：熵并不取决于区间的顺序。,的先验概率密度函数(?),28,的方差与偏置,一般来说，决定 的方程是非线性的。当得到正规函数后，将 在 附近展开,G. Cowan, StatisticalData Analysis, OxfordUniversity Press(1998),29,的方差与偏置(续),利用误差传递得到协方差,以及对偏置的估计量，,此处 而且通常情况下,30,正规化参数 的选取, 决定了置于数据的权重大小以便能与光滑度相比较， = 0 给出最大的光滑估计值，并与数据
10、无关。因此虽然方差为零，但有明显的偏置。而取大的 ，则回到高度振荡无偏的最大似然解。,为了在偏置与方差之间达到最大平衡:选择 使均值误差的平方（MSE）最小,或要求偏置不大于它自身的估计方差 。可以找到 的值使得,31,正规化参数选择的重要性,32,迭代法解谱(Bayes方法),应用贝叶斯定理,33,RooUnfold: SVD方法解谱,解谱需要先得到响应矩阵R，并进行求逆。对于性质不太好的矩阵，SVD方法更可靠一些。,34,Singular Value Decomposition(SVD),35,RooUnfold的下载、编译和使用,RooUnfold中提供了SVD以及Bayes的解谱法。下
11、载和编译使用方法见网页：http:/hepunx.rl.ac.uk/adye/software/unfold/RooUnfold.html在training服务器上，已经编译好放到ROOT的库文件目录中，不需要另行下载编译。使用时只要在ROOT脚本文件中加上 gSystem-Load(“libRooUnfold”);,36,RooUnfoldSvd正规化参数的选择,RooUnfoldSvd中正规化函数的选择实际上是Tikhonov(k=2)规则。使用RooUnfoldSvd后自动生成UnfoHist.root文件，其中包含d的直方图，从d的直方图可以读出正规化参数，即kterm=k,其中k之后
12、各个bin的值满足均值为0方差为1的高斯分布。当方差小于(或大于)1时，说明测量值b的误差被低估(或高估)了。不同的分布、分bin数目以及样本大小都会影响正规化参数的选择。训练响应矩阵的MC数据至少应该是测量数据的10倍以上。MC数据的分布与真实分布差别不大时，效果最好。一种做法是，用与真实分布有一定差别但差别不大的分布测试，选定分bin数目以及正规化参数，将该结果直接应用于测量数据。,37,正规化参数选择的例子,在该例子中，kterm选择为12左右。,38,RooUnfoldSvd例子,39,例子:Tikhonov规则(k=2),40,例子:最大熵(MaxEnt)规则,41,一个在图像处理中的最大熵例子,最大熵值方法常用于天文观测图像的重建,与点源的偏置较小,易于推广到两维以上的情况。,42,例子： 的谱函数,=,Eur. Phys. J. C11: 599, 1999,Tikhonov 规则,修正因子方法,由于探测器对两者影响各不相同，因此，需要用开拆法求出“真实”分布。,43,小结,数学上的原理求反应矩阵的逆修正因子正规的解谱过程估计量的方差与偏置正规化参数的选择RooUnfold以及解谱法例子,
