杨振伟 清华大学 第十一讲：解谱法（开拆法） 粒子物理与核物理实验中的数据分析 杨振伟 清华大学 第十一讲：解谱法（开拆法）

本讲要点 数学公式,反应函数(矩阵) 求反应矩阵的逆 修正因子 正规化的解谱法 估计量的方差与偏置 正规化参数的选择 举例 Tikhonov 规则 MaxEnt 规则

图像还原问题 一个常见的问题：由于实验仪器的原因而出现 图像变形，例如 真实分布 ＝ 真实分布 探测器响应 实验观测分布 如果已知探测器响应（可通过探测器模拟得到探测器响应的形式）， 能否还原出不受实验仪器影响的真实分布？ Unfolding(解谱法)

解谱问题的表述(1) 考虑随机变量y，目标：寻找概率密度函数f(y) 1)如果已知参数化形式 ， 最大似然法 1)如果已知参数化形式　　 ， 最大似然法 2)不知概率密度函数的形式，可以构造直方图 可以构造直方图(M个区间)，概率密度函数离散化 第j个区间的概率为 “真实的直方图” 为j(或pj)构造估计量 (参数的数目=区间的数目M) 离散化的p.d.f.

解谱问题的表述(2) 问题：y 的测量不可能没有误差 后果：f(y) 模糊化，峰展宽， 甚至分布形状完全变形。 1）第i个区间的y在第j个区间测量到 （测量分辨率） 2）第i个区间的某些事例没有测量到 （测量效率） 3）某些区间的事例完全测量不到 （接收度） 后果：f(y) 模糊化，峰展宽， 甚至分布形状完全变形。 后果严重时需要考虑解谱法还原真实分布

响应矩阵 真值； 观测值 测量误差的影响可用积分方程表示: 响应矩阵 Rij=P(观测值在第 i 区|真实值在第 j 区) ：响应函数，表示真值为y观测值为x的概率 对y求积分即可得到测量值为x的概率 离散化 观测直方图 (期待值) 真实直方图 响应矩阵 Rij=P(观测值在第 i 区|真实值在第 j 区)

响应矩阵 Rij=P(观测值在第 i 区|真实值在第 j 区) 真实直方图 观测数据 离散化的p.d.f. 和数据的期待值 数据： 这里 注意： 是常数， 会受 到统计涨落的影响。 真实直方图 离散化的p.d.f. 观测数据 和数据的期待值

效率、本底 有时，事例可能会不被探测到: 效率 有时，无真实事例发生，也有事例被观测到：本底 i 是在观测直方图上预期的本底数目。 真实直方图第 j 区探测效率 真值 效率 测量值 i 是在观测直方图上预期的本底数目。 有时，无真实事例发生，也有事例被观测到：本底

各关键量汇总 “真实”直方图: M 个区间 概率: 观测直方图的期待值: N 个区间 观测直方图: 响应矩阵: 效率: 预期的本底: 或关联矩阵(各区间不独立)Vij=cov[ni,nj] 一般通过构造 logL 或 2寻求 的估计量，这需要相关的概率理论，例如:泊松分布(各区间独立)

为什么要用解谱法？ 一般而言，我们并不需要解谱法，例如当比较现有理论的预期值时，最好是将探测器相应叠加到理论中去。即在预期值中包含探测器效应并与未修正的原始数据 相比较。 但是，不将实验数据进行解谱处理，结果发表后，有关反应矩阵的知识将不在保留。而且，解谱后的分布可以直接与各种理论的预言相比较，也可以与别的实验经过解谱以后的分布相比较。 通常解谱后的结果更为有用，否则当反应矩阵不可恢复时，即使对结果又有新的理论解释，也很难进行理论检验。 在粒子物理研究中，解谱法常用的领域为： 强子结构函数 的谱函数(也就是强子不变质量谱) 强子事例形状分布 粒子多重数分布 ．．．

为什么用解谱法(举例) 中微子与铁原子核相互作用，产生电子，测量电子的能谱。实验上只可能测量电子在铁中的射程。 解谱法 能量与平均射程的关系 因子修正法：能谱(圆圈)与真值(直方图) 解谱法

响应矩阵的逆 假设 的逆存在: 若数据是泊松分布 则有 最大似然法的估计量为 ？ 假设 的逆存在: 若数据是泊松分布 ？ 则有 最大似然法的估计量为 若Ｒ的非对角元太大，即区间宽度比分辨率要小时，会导致上式有很大的方差，以及在相邻区间产生很强的负关联。

错误的原因(I) 考虑一个简单的例子 解决办法是进行平滑处理，消除无意义的统计涨落。 但平滑会带来偏向性，需要在涨落与偏向性之间找到平衡。

错误的原因(II) 考虑周期分布函数 实际上，非物理的剧烈涨落主要是有高频部分造成的。平滑处理，主要是消除或压低高频部分无意义的涨落。

错误的原因(III) 假设 真的有精细结构 R-1 作用回 应完全恢复精细结构： 但我们没有 只有 （存在涨落） 假设 真的有精细结构 R作用后(理论上) 得到 大部分精细结构被抹平， 但会隐藏精细结构的信息。 R-1 作用回 应完全恢复精细结构： 但我们没有 只有 （存在涨落） R-1“认为”这是隐藏的精细结构信息，不遗余力地恢复这种精细结构，从而造成剧烈的振荡效应。

重新研究最大似然法的解(1) 估计量的平均值 无偏! 计算估计量的方差

重新研究最大似然法的解(2) 利用 RCF 边界做无偏估计量 倒数后得到 ML估计量在所有无偏估计中给出的方差最小。 但这个方差太大了！ 与刚才得到的结果一致。 ML估计量在所有无偏估计中给出的方差最小。 但这个方差太大了！ 为了减小方差，必须引入一些偏置量 策略：接受小的偏置量(系统误差)以换取大幅 减小方差(统计误差)。

简单方法：修正因子法 对 做相同的分区，并取 通常 ，因此方差不会被放大。 (相当于R-1取为对角矩阵) (R-1)ii=Ci 对 做相同的分区，并取 (相当于R-1取为对角矩阵) (R-1)ii=Ci 与 用蒙特卡罗模拟得到（无本底）。 通常 ，因此方差不会被放大。

修正因子法中的偏置问题 修正因子法的偏置 除非模拟采用的模型无误，使得 ，否则 上式不为零，需要考虑对应的系统误差。 除非模拟采用的模型无误，使得 ，否则 上式不为零，需要考虑对应的系统误差。 注意：该偏置量存在把 拉向 的倾向，造 成模型检验的困难。 1）如果分区宽度  几倍的分辨率，结果不会太坏。 2）实际应用中，中。该方法常用于事例形状变量的分布研究

例子：脉冲形状的还原 将理论(真实)的直方图除以受实验仪器影响 的直方图得到修正因子 将观测直方图乘以修正因子直方图得到理论 ＝ 将观测直方图乘以修正因子直方图得到理论 (真实)的直方图 ＝

正规化的解谱法 考虑“合理的”估计量，对选定的 满足 将下式求最大值 正则化函数(光滑性的量度) 正则化参量(其选择与给定的 对应) 考虑“合理的”估计量，对选定的　　 满足 估计量满足该不等式并且最光滑，等价于 将下式求最大值 正则化函数(光滑性的量度) 正则化参量(其选择与给定的 对应)

正规化的解谱法（续一） 另外，要求解谱后对总事例数的估计为无偏的 在约束情况下将下式求最大值 显然， :拉格朗日乘子

正规化的解谱法（续二） 需要正规化函数 以及选取  值的方案。 所得到的估计量的好坏由它们的偏置和方差 来判断。 给出最光滑的解(与数据无关) 给出最大似然解(方差太大) 需要正规化函数 以及选取  值的方案。 Tikhonov 规则 MaxEnt 规则 所得到的估计量的好坏由它们的偏置和方差 来判断。

Tikhonov 规则 取光滑度等于第 k 阶导数平方的均值，有 通常取 k=2，使得 S 约等于曲率平方的平均 值。对直方图而言，也就是 Sov. Math.5(1963)1035 注意：2 阶导数对直方图的第一和最后的区 间没有很好的定义。

Tikhonov 规则（续） 如果在 下，采用Tikhonov(k=2)规则 是 i 的二次项。对 求偏微分，给出线性方 程，得到 估计值与方差。 在高能物理界现有好几个现成的程序：RUN，Blobel, SVD, Höcker,…

最大熵(MaxEnt)规则 另一种表征光滑度的方法基于熵。对于一组概 率而言，熵定义为 所有 相等意味着熵最大(最光滑) Ann. Rev. Astron. Astrophys.24 (1986)127 所有 相等意味着熵最大(最光滑) 有一个 ，其它为零，则意味着熵最小。

最大熵(MaxEnt)规则（续） 用熵作为正规化函数， 有时侯，根据贝叶斯统计 这里，我们仍然采用经典近似: 估计量的好坏 的先验概率密度函数(?) 这里，我们仍然采用经典近似: 估计量的好坏 由偏置，方差来判断。 注意：熵并不取决于区间的顺序。

的方差与偏置 一般来说，决定 的方程是非线性的。当得 到正规函数后，将 在 附近展开 G. Cowan, Statistical 一般来说，决定 的方程是非线性的。当得 到正规函数后，将 在 附近展开 G. Cowan, Statistical Data Analysis, Oxford University Press(1998)

的方差与偏置(续) 利用误差传递得到协方差 以及对偏置的估计量， 此处 而且通常情况下

正规化参数  的选取  决定了置于数据的权重大小以便能与光滑度相比较， = 0 给出最大的光滑估计值，并与数据无关。因此虽然方差为零，但有明显的偏置。而取大的 ，则回到高度振荡无偏的最大似然解。 为了在偏置与方差之间达到最大平衡:选择 使均值误差的平方（MSE）最小 或要求偏置不大于它自身的估计方差 。可以找 到  的值使得

正规化参数选择的重要性

迭代法解谱(Bayes方法) 应用贝叶斯定理

RooUnfold: SVD方法解谱 解谱需要先得到响应矩阵R，并进行求逆。 对于性质不太好的矩阵，SVD方法更可靠一些。

Singular Value Decomposition(SVD)

RooUnfold的下载、编译和使用 RooUnfold中提供了SVD以及Bayes的解谱法。 下载和编译使用方法见网页：http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html 在training服务器上，已经编译好放到ROOT的库文件目录中，不需要另行下载编译。使用时只要在ROOT脚本文件中加上 gSystem->Load(“libRooUnfold”);

RooUnfoldSvd正规化参数的选择 RooUnfoldSvd中正规化函数的选择实际上是Tikhonov(k=2)规则。 使用RooUnfoldSvd后自动生成UnfoHist.root文件，其中包含d的直方图，从d的直方图可以读出正规化参数，即kterm=k,其中k之后各个bin的值满足均值为0方差为1的高斯分布。 当方差小于(或大于)1时，说明测量值b的误差被低估(或高估)了。 不同的分布、分bin数目以及样本大小都会影响正规化参数的选择。训练响应矩阵的MC数据至少应该是测量数据的10倍以上。 MC数据的分布与真实分布差别不大时，效果最好。 一种做法是，用与真实分布有一定差别但差别不大的分布测试，选定分bin数目以及正规化参数，将该结果直接应用于测量数据。

正规化参数选择的例子 有效部分 非物理振荡部分 在该例子中，kterm选择为12左右。

RooUnfoldSvd例子

例子:Tikhonov规则(k=2)

例子:最大熵(MaxEnt)规则

一个在图像处理中的最大熵例子 最大熵值方法常用于天文观测图像的重建,与点源的偏置较小,易于推广到两维以上的情况。

例子：  的谱函数 由于探测器对两者影响各不相同，因此，需要 用开拆法求出“真实”分布。 — = Tikhonov 规则 修正因子方法 Eur. Phys. J. C11: 599, 1999 由于探测器对两者影响各不相同，因此，需要 用开拆法求出“真实”分布。 — = Tikhonov 规则 修正因子方法

小结 数学上的原理 求反应矩阵的逆 修正因子 正规的解谱过程 估计量的方差与偏置 正规化参数的选择 RooUnfold以及解谱法例子