使用scikit-hep/hepstats计算CLs上限时如何纳入系统不确定性？

在hepstats中，系统不确定性的处理核心是通过**nuisance参数（额外约束参数）**建模，将不确定度转化为可拟合的参数并施加合理的先验约束，最终在CLs计算时通过轮廓化似然（profile likelihood）自动纳入这些参数的影响。结合你的场景（信号高斯分布、背景多项式拟合、类似H→γγ分析），具体实现步骤如下：

1. 信号系统不确定性的建模 #

针对信号的不同不确定来源，分别构建nuisance参数：

• 归一化不确定（如截面、探测效率）：给信号预期事件数引入一个比例型nuisance参数，用高斯或对数正态约束。比如假设信号预期事件数为N_s，相对不确定度为10%，则引入参数alpha_s，约束为alpha_s ~ N(1, 0.1²)，实际信号事件数为N_s * alpha_s。
• 形状不确定（如高斯峰位、宽度偏移）：将高斯分布的均值/宽度设为带约束的nuisance参数。比如原本固定的峰位mu0改为mu = mu0 + delta_mu，delta_mu ~ N(0, sigma_mu²)，其中sigma_mu是峰位的不确定度。

2. 背景系统不确定性的建模 #

你的背景由数据拟合得到，需覆盖两类不确定：

• 背景拟合的统计不确定：将多项式的系数设为nuisance参数，用拟合得到的协方差矩阵施加多元高斯约束。如果是基于Asimov数据集拟合的背景，要确保将拟合参数的误差传递到CLs计算中。
• 背景的系统不确定（如归一化、形状偏移）：类似信号的处理方式，引入比例型或偏移型nuisance参数并施加约束。比如背景归一化不确定度5%，则引入alpha_b ~ N(1, 0.05²)，背景事件数为N_b * alpha_b。

3. 代码层面的实现步骤 #

3.1 构建带nuisance参数的模型 #

from hepstats.model import Model
from hepstats.parameters import Param
import numpy as np
from scipy.stats import norm, poisson

# 定义信号和背景PDF
def signal_pdf(x, mu, sigma, alpha_s):
    return norm.pdf(x, loc=mu, scale=sigma) * alpha_s  # alpha_s控制信号归一化

def background_pdf(x, coeffs, alpha_b):
    poly = np.polyval(coeffs, x)
    return poly * alpha_b  # alpha_b控制背景归一化

# 初始化模型
model = Model()

# 添加信号的nuisance参数（归一化不确定10%）
alpha_s = Param("alpha_s", value=1.0, bounds=(0.5, 1.5))
model.add_param(alpha_s, prior="normal", mean=1.0, sigma=0.1)

# 添加背景的nuisance参数（归一化不确定5%）
alpha_b = Param("alpha_b", value=1.0, bounds=(0.5, 1.5))
model.add_param(alpha_b, prior="normal", mean=1.0, sigma=0.05)

# 添加其他参数（如高斯mu、sigma，多项式系数）
# ... 根据你的实际拟合需求添加固定或浮动参数

3.2 构建似然函数并计算CLs #

hepstats的cls模块会自动处理nuisance参数的轮廓化：

from hepstats.cls import compute_cls

# 假设你已经有数据data，以及信号/背景的预期参数
# 构建似然函数（泊松似然，适用于计数类数据）
def likelihood(params):
    mu = params["mu"]  # 信号强度参数（用于计算截面上限）
    alpha_s_val = params["alpha_s"]
    alpha_b_val = params["alpha_b"]
    # 计算信号和背景的预期事件数
    s = signal_pdf(data, mu=125, sigma=1.0, alpha_s=alpha_s_val).sum() * bin_width
    b = background_pdf(data, coeffs=[0.1, -0.01], alpha_b=alpha_b_val).sum() * bin_width
    # 泊松似然
    return poisson.logpmf(data_counts, s + b).sum()

# 计算CLs上限
cls_results = compute_cls(
    likelihood,
    poi_param="mu",  # 感兴趣的参数（信号强度，对应截面）
    poi_bounds=(0, 5),
    ntoys=1000,  # 玩具数，根据精度需求调整
    model=model
)

# 提取95% CLs上限
upper_limit = cls_results["upper_limit"]

4. 关键注意事项 #

• 对于相关的系统不确定，需使用多元高斯约束（而非独立单变量约束）来建模参数间的相关性，确保不确定度传递准确。
• 选择合适的约束分布：比例类不确定（如效率、归一化）优先用对数正态约束，偏移类不确定（如峰位偏移）用高斯约束。
• 若背景拟合的统计不确定占主导，可通过生成多个背景拟合的伪实验，将其结果纳入CLs的统计涨落分析中。

内容的提问来源于stack exchange，提问作者Diptaparna Biswas