【高能物理.ROOT】一个Bin的信号强度上限的假设检验(上)

一个Bin的信号强度上限的假设检验

如有错误，欢迎指正！

问题

给定预计的信号个数$S$,本底个数$B$,求出给定数据$N_{obs}$下,信号强度$\mu$的上限.
解释:
一个Bin内真实情况应该有$\mu S+B$个,数据$N_{exp}$应该服从以$\mu S+B$为均值,标准差$\sqrt{\mu S+B}$的泊松分布.但是实际数据$N_{obs}$可能产生了一个值,通过泊松分布来就可以判断哪些$\mu$值是不可能的。
泊松分布的概率密度函数为:
$$f(n|\lambda )=\frac{{\lambda }^n{e}^{-\lambda }}{n!}$$
泊松分布的期望值为$\lambda$,标准差为$\sqrt{\lambda }$.
例如:
假设有一组数据$N_{obs}=10$,预计的（或者说模型的）信号个数$S=4$,本底个数$B=5$,显然当$\mu ==2.5$的时候,
$$f(10|2.5\ast 4+5)=\frac{{\lambda }^n{e}^{-\lambda }}{n!}=0.0486$$
这个概率很小说明$\mu$不太可能在这个值。扫描$f(10|\mu \ast 4+5)>0.05$就可以得到一个$\mu$的区间.但是这个检验方法足够好吗？后面将用假设检验来获得上限。

理论

最大似然函数

实际判选会更加复杂，不过基本结论是一致的,判断在这个$\mu$时的发生的概率应该大于0.05.(正常说法是[拒绝$\mu$当在此$\mu$时给定数据发生的概率小于0.05])
实际中往往用最大似然函数来描述这个概率
$$L(\mu ;\Theta )=L_{else}(\mu ;\Theta )\prod _{Bin=1}^{N_{Bin}}\frac{(\mu S(\Theta )+B(\Theta ){)}^{n_{i,obs}}{e}^{-\mu S(\Theta )-B(\Theta )}}{n_{i,obs}!}$$
其中$\mu ,\Theta$是模型参数,比如说有$N_{Bin}$个Bin,每个Bin的信号个数$S(\Theta )$和本底个数$B(\Theta )$是不同的.观测的结果$n_{i,obs}$也是不同的.
其中$L_{else}(\mu ;\Theta )$是模型一种先验分布或者不关心的区域测量的最大似然函数,

比如说想要限制$\mu$到$\mu =1.0\pm 0.1$,那么$L_{else}(\mu ;\Theta )=\exp [-\frac{(1.0-\mu {)}^2}{2(0.1{)}^2}]$,即$\mu$的限制(或者说constrain)是均值为1.0,方差为0.1的正态分布.这里因为最大似然函数可以差任意一个倍数所以归一化没写。

比如有一个区域没有信号，但是这块位置的拟合可以给信号区域的拟合提供帮助，那么$L_{else}(\mu ;\Theta )=\prod \frac{B^{n_i}{e}^{-B}}{n_i!}$,即边带区域的似然函数(S=0)。

Asimov数据¹

Asimov数据是指用当前的模型产生一个数据,这个数据具有一些好的性质:

• 它是模型中最有可能出现的结果
• 它可以用来估计模型参数的不确定性
• 它可以用来获得假设检验的预期结果及不确定性

Asimov满足(定义):
$$\frac{\partial \log L(n_{i,Asimov};\theta )}{\partial \theta }=0$$
其中$\theta$是模型参数($\mu$,$\Theta$),而$n_{i,Asimov}$是Asimov数据每个Bin的“观测值”.

参数估计及区间估计

way 1

最大似然估计法可以获得最大似然估计值，
$$\hat{\theta }\to s.t.\max L(n;\theta )$$
区间估计从最大似然函数曲线与一条$y=\max L-1/2$的直线的交点得到。
$$({\theta }_{min},{\theta }_{max})\to s.t.L(n;{\theta }_{min})=L(n;{\theta }_{max})=\max L-s^2/2$$
其中$s$表示几倍的$\sigma$值,比如说$s=1.0$表示68%置信区间.

way 2

利用似然比可以得到置信区间
$$\lambda (\theta )\equiv \frac{L(n;\theta )}{L(n;\hat{\theta })},\hat{\theta }\to s.t.\max L(n;\theta )$$
由Wilks定理,$-2ln(\lambda (\theta ))\sim {\chi }^2(m)$,所以可以得到
$$-2ln(\lambda ({\theta }_{min}))=-2ln(\lambda ({\theta }_{max}))=F_{{\chi }^2}^{-1}(1-\alpha ,m)$$
协方差估计:
$$\hat{V}({\theta }_i)=(-\frac{{\partial }^{2}ln(L(n_i;\theta ))}{{\partial }^2{\theta }_i}{|}_{\hat{\theta }}{)}^{-1}$$

特别的,边缘似然比:
$$\lambda ({\theta }_1)\equiv \frac{L(n_i;{\theta }_1,{\hat{\hat{\theta }}}_2)}{L(n_i;\hat{\theta })},{\hat{\hat{\theta }}}_2({\theta }_1)\to 给定{\theta }_1,s.t.\max L(n_i;{\theta }_1,{\theta }_2)$$

假设检验

假设检验一定要有一个检验统计量，对于特定问题有特定的假设检验统计量。
奈曼-皮尔逊引理可以确定$t\equiv -2ln(\lambda (\theta ))$是最优统计量。定义$t_{\mu }\equiv -2ln(\lambda (\mu ))$,其中$\mu$是信号强度，即跟前面一致的$\mu$,这里是边缘似然比,其他参数看不到了。由于实际问题往往有$\mu >0$的限制实际的统计量稍微复杂一些。

现在来做一个上限的假设检验:

做上限的假设检验一般是不能排除$\mu =0$情况下的检验。(也称为灵敏度分析)
对于一个给定的${\mu }_{test}$值做假设检验。

H0: $\mu ={\mu }_{test}$,H1: $\mu =0$

$t_{\mu ,test}$在当模型$\mu ={\mu }^{\prime }$的分布是清楚的，若${\mu }_{test}={\mu }^{\prime }$则服从中心卡方分布${\chi }^2(m)$,m为参数个数(这里只有一个参数$\mu$,m=1)。反之服从非中心卡方分布。
那么其实目前很清晰了，(如果真实世界确实是在当前模型下,一般假定)，若真实世界服从${\mu }^{\prime }$的模型(注意：我们不知道这个数值)，那么我们只需要确定$t_{\mu ,test}$的观测数值是不是比较符合中心卡方分布即可！
如果符合，那么说明当前检验的${\mu }_{test}$有可能是正确的。不能拒绝H0。
如果不太符合(例如用$p_0$量化描述一下不可能程度$p_0<0.05$)，那我们就应该拒绝当前${\mu }_{test}$作为真实信号强度的假设(当前${\mu }_{test}$非常有可能不对)。拒绝H0
通过不断调整${\mu }_{test}$的值，预计能够得到信号强度$\mu$的上限。

CLs,CLsb,CLb

大致方向清楚了，但是细致想来如何量化"符合"，"不符合"呢？

上面一张图描述了$t_{\mu ,test}$的分布情况，上面两个分布一个是真实世界${\mu }^{\prime }={\mu }_{test}$情况下$t_{\mu ,test}$应该服从的分布H0，
另一个是真实世界${\mu }^{\prime }=0$的情况下$t_{\mu .test}$的分布，上面说过,若${\mu }_{test}\ne {\mu }^{\prime }$那么$t_{\mu ,test}$服从非中心卡方分布，即图中的所示分布H1。
一个简单方法衡量图中的$p_0$值($p_{s+b},\alpha$),如果$p_0<0.05$就拒绝H0。
但是如果本底$B$很多(H1和H0分布几乎一致,因为$S$影响不大),实际上如果$p_0<0.05$，很有可能这只是因为本底的波动导致的，不应该拒绝H0.
所以，我们需要一个更加严格的检验方法。

CLs就是一个更严格的检验方法。

• $C{L}_s=\frac{p_0}{p_1}=\frac{p_{s+b}}{1-P_b}=\frac{C{L}_{s+b}}{C{L}_b}=\frac{\alpha }{1-\beta }$
• $C{L}_{s+b}=p_{s+b}=p_0$
• $C{L}_b=p_1=1-P_b$

值的注意的是,检验统计量是可以不同的选择的，不过$C{L}_s$判断方法是一致的。通过找到参数值${\mu }_{test}$,使得在对$t_{\mu ,test}$的检验发现$C{L}_s==0.05$，对应的${\mu }_{test}$就是在目前数据情况下给出的信号强度的上限。

提到数据,自然联想到之前说的Asimov数据,通过Asimov数据能够估计一个信号强度$\mu$的上限，以及其不确定性。例如在纯本底假设下$\mu =0$的Asimov数据，就能得到纯本底假设下预期上限即不确定度。

下一步代码实现。

---

1. https://arxiv.org/abs/1007.1727v3 ↩︎