Transcript 正规化的解谱法 - 清华高能物理中心

粒子物理与核物理实验中的
数据分析
杨振伟
清华大学
第十一讲：解谱法（开拆法）
本讲要点
数学公式,反应函数(矩阵)
 求反应矩阵的逆
 修正因子
 正规化的解谱法

a) Tikhonov 规则
b) MaxEnt 规则
估计量的方差与偏置
 正规化参数的选择
 举例

2
图像还原问题
一个常见的问题：由于实验仪器的原因而出现
图像变形，例如

＝
探测器响应
实验观测分布
真实分布
如果已知探测器响应（可通过探测器模拟
得到探测器响应的形式），
能否还原出不受实验仪
器影响的真实分布？
Unfolding(解谱法)
3
解谱问题的表述(1)
考虑随机变量y，目标：寻找概率密度函数f(y)

1)如果已知参数化形式 f ( y; ) ，
ˆ
最大似然法


f ( y; )
2)不知概率密度函数的形式，可以构造直方图
可以构造直方图(M个区间)，概率密度函数离散化
第j个区间的概率为
pj  
bin j
f ( y )dy
j  1,..., M
 j  tot p j
“真实的直方
图”
为j(或pj)构造估计量
离散化的p.d.f.
(参数的数目=区间的数目M)
4
解谱问题的表述(2)
问题：y 的测量不可能没有误差
1）第i个区间的y在第j个区间测量到
（测量分辨率）
2）第i个区间的某些事例没有测量到
（测量效率）
3）某些区间的事例完全测量不到
（接收度）
测
量
后果：f(y) 模糊化，峰展宽，
甚至分布形状完全变形。
后果严重时需要考虑解谱法还原真实分布
5
响应矩阵
y : 真值； x : 观测值
测量误差的影响可用积分方程表示:
f meas ( x )   R( x | y ) f true ( y )dy
R( x | y) ：响应函数，表示真值为y观测值为x的概率
对y求积分即可得到测量值为x的概率
离散化
观测直方图
(期待值)
M
 i   Rij  j ,
真实直方图
i  1, ..., N
j 1
响应矩阵 Rij=P(观测值在第 i 区|真实值在第 j 区)
6
响应矩阵
Rij=P(观测值在第 i 区|真实值在第 j 区)
数据： n  (n1 ,..., nN )
这里  i  E[ni ]
真实直方图
离散化的p.d.f.

 
注意： , 是常数，n 会受
到统计涨落的影响。
观测数据
和数据的期待值
7
效率、本底
有时，事例可能会不被探测到: 效率
N
N
 R   P (观测值在第 i 区 | 真实值在第 j 区)
i 1
ij
i 1
 P (观测值在全范围 | 真实值在第j区)
  j (效率)
真值

效率
真实直方图第
j 区探测效率
测量值
有时，无真实事例发生，也有事例被观测到：本底
M
 i   Rij  j   i
j 1
i 是在观测直方图上预期的本底数目。
8
各关键量汇总
M
tot    j
“真实”直方图:  ( 1 ,...,  M ),
j 1
概率: p  ( p1 ,..., pM )   / tot
观测直方图的期待值:   ( 1 ,..., N )
观测直方图: n  (n1 ,..., nN )
M 个区间
N 个区间
响应矩阵: Rij  P(观测值在第 i 区| 真实值在第 j 区)
N
效率:  j   Rij
i 1
预期的本底:   ( 1 ,...,  N )
E[n]    R  
一般通过构造 logL 或 2寻求  的估计量，这需要相关的
概率理论，例如:泊松分布(各区间独立)
n
P ( ni ; i ) 
i
i
e  i
ni !
或关联矩阵(各区间不独立)Vij=cov[ni,nj]
9
为什么要用解谱法？
一般而言，我们并不需要解谱法，例如当比较现有理
论的预期值时，最好是将探测器相应叠加到理论中去。即在

n 相比较。
预期值中包含探测器效应并与未修正的原始数据
但是，不将实验数据进行解谱处理，结果发表后，有
关反应矩阵的知识将不在保留。而且，解谱后的分布可以直
接与各种理论的预言相比较，也可以与别的实验经过解谱以
后的分布相比较。
通常解谱后的结果更为有用，否则当反应矩阵不可恢复
时，即使对结果又有新的理论解释，也很难进行理论检验。
在粒子物理研究中，解谱法常用的领域为：
• 强子结构函数
• 的谱函数(也就是强子不变质量谱)
• 强子事例形状分布
• 粒子多重数分布
•．．．
10
为什么用解谱法(举例)
中微子与铁原子核相互作
用，产生电子，测量电子
的能谱。实验上只可能测
量电子在铁中的射程。
能量与平均射程的关系
解谱法
因子修正法：能谱(圆圈)与真值(直方图)
11
响应矩阵的逆
假设   R   的逆存在:   R1 (   )
若数据是泊松分布

n
P ( ni ; i ) 
则有
i
i
ni !
e
 i

N
log L(  )   ( ni log i   i )
i 1
最大似然法的估计量为
ˆ  n
ˆ  R1 (n   )
若Ｒ的非对角元太大，
即区间宽度比分辨率要
小时，会导致上式有很
大的方差，以及在相邻
区间产生很强的负关联。
n
ˆ
？
12
错误的原因(I)
考虑一个简单的例子
 x1 
 b1 
ˆ
A x  b 其中x   , b   ，
 x2 
 b2 
1 1   1   
ˆ
 , 〈
A  
0 〈1
2 1   1   
  1 : 完美探测器
  0 : 分辨率很差
ˆA1  1  1    1   
2   1   1   
1 1 1 1  1  1

    
2 1 1 2   1 1
b1  b2 1 b1  b2  1 
1
ˆ
  
 
x A b
2 1
2   1
很小时，第2项起决定作用
尤其当b1 ~ b2时，考虑到统计涨落，b1  b2实际上是一个随机数。
这个随机数被1 / 放大后，结果中有用的信息完全被非物理振荡湮没。
所以，通常情况下，直接求反应矩阵获得真值x的方法，尽管理论上
是严格的，是无偏有效的，但结果并没有物理意义。
解决办法是进行平滑处理，消除无意义的统计涨落。
但平滑会带来偏向性，需要在涨落与偏向性之间找到平衡。 13
错误的原因(II)
考虑周期分布函数
a0 
f ( x)    (a cosx b sinx), 0  x  2
2  1
如果每一点的测量都按照偏差为 的高斯分布弥散，那么测量到的分布为

 ( y  x) 2 
g ( y)  
exp  
 f ( x)dx (卷积)
2
2
2 
2


将g(y)展开为 cos x和 sin x的函数：
g(y) 
0
2
1

  ( cos x  sin x)
 1
实际上，非物理的剧烈涨落
主要是有高频部分造成的。
平滑处理，主要是消除或压
低高频部分无意义的涨落。
  2 2 
  2 2 
其中 a  exp 
  , b  exp 
 
 2 
 2 
如果可以准确知道g(y)的系数 和 ，则可严格求出真实分布f ( x)的系数a 和b ，
但 和 都有一定的统计误差，这个统计误差被exp( 2 2 / 2)放大。
尤其是对于一般函数，  时, a , b  0,  和 中的统计误差被指数放大后
将完全占据统治地位。
14
错误的原因(III)

假设  真的有精细结构


R作用后(理论上)


得到   R

大部分精细结构被抹平，
但会隐藏精细结构的信息。


 应完全恢复精细结构：
1
R 
R-1 作用回


但我们没有 只有 n（存在涨落）
R-1“认为”这是隐藏的精细结构信息，不遗余力地恢复
这种精细结构，从而造成剧烈的振荡效应。
15
重新研究最大似然法的解(1)
估计量的平均值
1
ˆ
E[ ]  R ( E[n]   )  
无偏!
计算估计量的方差
U ij  cov[ ˆ i , ˆ j ] 
N
1
1
(
R
)
(
R
) jl cov[nk , nl ]

ik
k , l 1
N
  ( R 1 )ik ( R 1 ) jk k
k 1
假设ni是独立的泊松变量
时，cov[nk , nl ]   kl k
16
重新研究最大似然法的解(2)
利用 RCF 边界做无偏估计量
2


log L  N Rik Ril
1
(U )kl   E 

 k l  i 1  i
倒数后得到
  R1 (   )
N
log L(  )   ( ni log i   i )
i 1
N
U ij   ( R1 )ik ( R1 ) jk k
i 1
与刚才得到的结果一致。
ML估计量在所有无偏估计中给出的方差最小。
但这个方差太大了！
为了减小方差，必须引入一些偏置量
策略：接受小的偏置量(系统误差)以换取大幅
减小方差(统计误差)。
17
简单方法：修正因子法
对 , 做相同的分区，并取
ˆ i  Ci (ni   i )(相当于R-1取为对角矩阵)

Ci 

MC
i
MC
i
(R-1)ii=Ci
(修正因子)
 iMC与 iMC 用蒙特卡罗模拟得到（无本底）。
Uij  cov[ˆ i , ˆ j ]  C cov[ni , nj ]
2
i
通常 Ci  O(1) ，因此方差不会被放大。
18
修正因子法中的偏置问题
修正因子法的偏置
bi  E[ˆ i ]  i
 iMC i  sig
bi   MC  sig  i , 其中 isig   i   i =  i i。
i 
 i
除非模拟采用的模型无误，使得 iMC  i ，否则
上式不为零，需要考虑对应的系统误差。
注意：该偏置量存在把 ˆ 拉向  MC 的倾向，造
成模型检验的困难。
1）如果分区宽度  几倍的分辨率，结果不会太坏。
2）实际应用中，中。该方法常用于事例形状变量的分布研究
19
例子：脉冲形状的还原
将理论(真实)的直方图除以受实验仪器影响
的直方图得到修正因子

＝
将观测直方图乘以修正因子直方图得到理论
(真实)的直方图

＝
20
正规化的解谱法
 log L
考虑“合理的”估计量，对选定的
满足
log L(  )  log Lmax   log L
 log L描述了数据 n 与期
待值 之间的“ 距离” 。
 估计量满足该不等式并且最光滑，等价于
将下式求最大值
(  )   log L(  )  S (  )
S (  )  正则化函数(光滑性的量度)
  正则化参量(其选择与给定的  log L对应)
21
正规化的解谱法（续一）
另外，要求解谱后对总事例数的估计为无偏的
N
N
   R 
i 1
i
ij
j
 ntot
因  R    ,
所以是 的函数
i, j
在约束情况下将下式求最大值
N


 (  ,  )   log L(  )  S (  )    ntot   i 
i 1


显然，
 :拉格朗日乘子
 /   0

N

i 1
i
 ntot
22
正规化的解谱法（续二）
N


 (  ,  )   log L(  )  S (  )    ntot   i 
i 1


  0 给出最光滑的解(与数据无关)
   给出最大似然解(方差太大)
需要正规化函数 S (  ) 以及选取  值的方案。
a) Tikhonov 规则
b) MaxEnt 规则
所得到的估计量的好坏由它们的偏置和方差
来判断。
23
Tikhonov 规则
取光滑度等于第 k 阶导数平方的均值，有
2
 d f true ( y ) 
S  f true ( y )    
 dy , 其中k  1, 2...
k
dy


k
通常取 k=2，使得 S 约等于曲率平方的平均
值。对直方图而言，也就是
M 2
S (  )    (  i  2 i 1  i  2 )2
Sov. Math.5(1963)1035
i 1
注意：2 阶导数对直方图的第一和最后的区
间没有很好的定义。
24
Tikhonov 规则（续）
如果在 log L  
1 2

2
下，采用Tikhonov(k=2)规则
 (  ,  )   log L(  )  S (  )
 2
   ( )  S( )
2
是 i 的二次项。对 求偏微分，给出线性方
程，得到 估计值与方差。
在高能物理界现有好几个现成的程序：
RUN，Blobel, SVD, Höcker,…
25
最大熵(MaxEnt)规则
另一种表征光滑度的方法基于熵。对于一组概
率而言，熵定义为
M
H    pi log pi
i 1
i
pi 
tot
Ann. Rev. Astron. Astrophys.24 (1986)127
所有 pi 相等意味着熵最大(最光滑)
有一个 pi  1 ，其它为零，则意味着熵最小。
26
最大熵(MaxEnt)规则（续）
用熵作为正规化函数，
i
i
S (  )  H (  )  
log
tot
i 1 tot
 log( tot填入M个区间中各种可能的总数)
M
有时侯，根据贝叶斯统计
S (  )   的先验概率密度函数(?)
这里，我们仍然采用经典近似: 估计量的好坏
由偏置，方差来判断。
注意：熵并不取决于区间的顺序。
27
ˆ
 的方差与偏置
ˆ 
一般来说，决定  (n) 的方程是非线性的。当得

ˆ 
到正规函数后，将  (n)在 nobs 附近展开
ˆ ( n )  ˆ obs  A1 B( n  nobs ),
  2
, i , j  1, ..., M ,

  i  j
  2
 1，i  1, ..., M , j  M  1,
Aij  
  i 
  2
 2  0, i  M  1, j  M  1,
 
  2
, i  1, ..., M , j  1, ..., N ,

  i n j
Bij   2
    1, i  M  1, j  1, ..., N .
  n j

为非正规的似然函数
G. Cowan, Statistical
Data Analysis, Oxford
University Press(1998)
28
ˆ
 的方差与偏置(续)
利用误差传递得到协方差
U  CVC
T
Uij  cov[ˆ i , ˆ j ],
1
C  A B,
其中
Vij  cov[ni , n j ]
以及对偏置的估计量，bi  E[ˆ i ]  i ,
ˆ i
ˆb  C (ˆ  n ) 
(ˆ j  n j ),


i
ij
j
j
j 1
j 1 n j
N
N
ˆ
ˆ 
此处  R   .而且通常情况下ˆ  n
29
正规化参数  的选取
 决定了置于数据的权重大小以便能与光滑度
相比较， = 0 给出最大的光滑估计值，并与数据无
关。因此虽然方差为零，但有明显的偏置。而取大
的 ，则回到高度振荡无偏的最大似然解。
为了在偏置与方差之间达到最大平衡:选择
使均值误差的平方（MSE）最小
M
ˆ2
U

b
1 M
1
ii
i
ˆ 2 , 或带权重的 MSE 
MSE 
U

b


M i 1 ii i
M i 1 ˆ i


或要求偏置不大于它自身的估计方差 Wˆ ii 。可以找
到  的值使得
ˆ2
b
 b2   i  M
ˆ
i 1 W
M
ii
这里
Wˆ ij  cov[bˆi , bˆ j ].
30
正规化参数选择的重要性
  0 : 光滑，完全MC，
与数据无关
   : 完全数据，
剧烈振荡，
无物理意义
31
迭代法解谱(Bayes方法)
Rij  P(测量值在i区间 | 真值在j区间)
利用Bayes定理，可以得到估计量
ˆ i 
1
i
N
 P(真值在i区间 | 测量值在j区间)n
j 1
j
P(测量值在j区间 | 真值在i区间) pi
 
nj
 i j 1  P(测量值在j区间 | 真值在k区间) pk
1
N
应用贝叶斯
定理
k

1
i
R ji pi
N

j 1
R
jk
pk
nj
k
选取初始pi(i  1, ,N)，之后用pi  μˆ i /μtot 迭代。
用测试数据确定迭代次数，之后用于实验数据。
当迭代次数很大时，结果趋于剧烈振荡。
32
RooUnfold: SVD方法解谱
解谱需要先得到响应矩阵R，并进行求逆。
对于性质不太好的矩阵，SVD方法更可靠一些。
R可以通过MC获得：
x : MC真值
b : 模拟测量结果
作x  b散点分布(h 2设x和b都分2个bin)
实际上获得了一个2  2的数组(或矩阵) A
A(1,1) : 真值为x1 , 测量值为b1的事例数
 R11
A(1, 2) : 真值为x1 , 测量值为b2的事例数
 R21
A(2,1) : 真值为x2 , 测量值为b1的事例数
 R12
A(2, 2) : 真值为x2 , 测量值为b2的事例数
 R22
Rij 
1
A( j , i ) ,
 A( j, i)
i
即 R  AT
33
Singular Value Decomposition(SVD)
任意n  m矩阵A可以分解成
A  USV T
它们通过对角矩阵S联系起来, 对角元素为
其中U  AAT的本征矢构成的n  n矩阵
V  AT A的本征矢构成的m  m矩阵
S  diag(eig ( AA )) n  m对角矩阵
T
U , V满足 UU T  I,
SVD 之后，已知量b变为d，未知量x变为z，
VV T  I
S ii  S jj , 如果i  j。显然A1  VS 1U T
矩阵A的奇异值。
实际上，奇异值很小的部分是结果剧烈
振荡的根源(参见2维例子中  0的情况),
对应于傅立叶分解中高频部分，相应的，
向量d中相应的元素的服从标准正态分布。
等价地说，低频部分是物理的，高频部分
需要截断，从d的元素分布可以找到截断。
Ax  b  USV T x  b  SV T x  U T b
定义 z  V T x， d  U T b
 Sz  d
 z  S 1d  x  Vz  VS 1d  VS 1U T b  A1b
SVD 的结果是把b和x都分解成正交归一向量的线性组合，
其系数分别构成向量d和z，基分别是矩阵U和V的列。
34
RooUnfold的下载、编译和使用

RooUnfold中提供了SVD以及Bayes的解谱法。
下载和编译使用方法见网页：
http://hepunx.rl.ac.uk/~adye/software/
unfold/RooUnfold.html
在training服务器上，已经编译好放到ROOT的
库文件目录中，不需要另行下载编译。使用时
只要在ROOT脚本文件中加上
gSystem->Load(“libRooUnfold”);
35
RooUnfoldSvd正规化参数的选择






RooUnfoldSvd中正规化函数的选择实际上是
Tikhonov(k=2)规则。
使用RooUnfoldSvd后自动生成UnfoHist.root文件，其中
包含d的直方图，从d的直方图可以读出正规化参数，即
kterm=k,其中k之后各个bin的值满足均值为0方差为1的
高斯分布。
当方差小于(或大于)1时，说明测量值b的误差被低估(或高
估)了。
不同的分布、分bin数目以及样本大小都会影响正规化参数
的选择。训练响应矩阵的MC数据至少应该是测量数据的10
倍以上。
MC数据的分布与真实分布差别不大时，效果最好。
一种做法是，用与真实分布有一定差别但差别不大的分布
测试，选定分bin数目以及正规化参数，将该结果直接应用
于测量数据。
36
正规化参数选择的例子
有效部分
非物理振荡部分
在该例子中，kterm选择为12左右。
37
RooUnfoldSvd例子
38
例子:Tikhonov规则(k=2)
39
例子:最大熵(MaxEnt)规则
40
一个在图像处理中的最大熵例子
最大熵值方
法常用于天
文观测图像
的重建,与点
源的偏置较
小,易于推广
到两维以上
的情况。
41
例子：  的谱函数
为了测定奇异夸克质量，实验上可采用比较
  X ( ud )  与  X ( us )  中,X 的质量平方差
ms
Eur. Phys. J. C11: 599, 1999
00  M 2 (ud )  M 2 (us)
由于探测器对两者影响各不相同，因此，需要
用开拆法求出“真实”分布。
—
Tikhonov 规则
=
修正因子方法
42
小结
1. 数学上的原理
真实直方图，数据 n 以及其期待值  ，满足   R   ,目标是构造  的估计量。
2. 求反应矩阵的逆
有很大的振荡行为( 及大的方差) ，但在各种无偏解中具有零偏置与最小的方差。
3. 修正因子
Ci  iMC / iMC , 方法又快又简单。
4. 正规的解谱过程
Tikhonov:从第 k阶导数的均方值中进行光滑处理
MaxEnt:从H   pi log pi 熵中进行光滑处理
5. 估计量的方差与偏置
i
在求解过程中采用了线性近似，因而不是无偏的
6. 正规化参数的选择
无最好的方案，可以采用 2  M(区间总数)的方案。
7. RooUnfold以及解谱法例子
只要探测器的响应可知，就一定可以得到真实的分布
43