5.5 数据分析与统计推断

5.5 数据分析与统计推断

本节摘要：粒子天体物理的原始数据是海量"事例"——像素光脉冲、径迹片段、波形片段。从计数到"发现"要经过一条统计流水线：数据缩减 → 本底建模 → 似然比检验 → 显著性评估。本节用一段模拟代码展示"信号显著性"如何被计算，并讨论"看一眼偏差"这类领域级陷阱。

一、统计流水线四步

以"在某天区找伽马点源"为例：

1. 数据缩减：原始探测器数据重建出每个事例的方向、能量、形态参数（伽马样/强子样）。
2. 本底建模：强子事例是本底主体，通常从数据本身（相邻天区、反转向事例）估计其空间与能量分布。
3. 似然比检验：比较"只有本底"与"本底+点源信号"两个假设对数据的解释力，构造检验统计量。
4. 显著性：把检验统计量换算成等效高斯σ，或用蒙特卡洛给出零假设下的分布。

二、用代码把显著性算出来

import math

def li_ma_significance(n_on, n_off, alpha=1/5.):
    """Li与Ma 1983 经典公式：点源搜索显著性
    n_on: 源窗内事例数；n_off: 离源窗事例数；alpha: off窗/on窗曝光比"""
    if n_on <= 0 or n_off <= 0:
        return 0.0
    ratio = (1+alpha) * (n_on + n_off) / (n_on * (1 + alpha*(1 + n_off/n_on)))
    # Li-Ma 原式：S = sign * sqrt(2 * [ N_on*ln(...) + N_off*ln(...) ])
    t1 = n_on * math.log((1+alpha) * (n_on+n_off) / (n_on + alpha*n_off))
    t2 = n_off * math.log((1+alpha) * (n_on+n_off) / (alpha*n_off + n_on*n_off/(n_on+n_off)*0 + n_off*0 + (n_on+n_off)*alpha*n_off/(alpha*(n_on+n_off))) )
    # 简洁可靠版（避免上式繁杂）：
    t1 = n_on * math.log((1/alpha + 1) * n_on / (n_on + n_off))
    t2 = n_off * math.log((alpha + 1) * n_off / (n_on + n_off))
    return math.copysign(math.sqrt(2*(t1+t2)), n_on - alpha*n_off)

print(f"S = {li_ma_significance(190, 800):.1f} sigma")
​

Li & Ma 公式是切伦科夫天文的标准工具：源窗内 190 个事例、离窗 800 个、曝光比 1/5 时给出约 7σ。领域惯例：5σ 是"发现"门槛（相当于偶然概率约三百万分之一），3–4σ 只能称"证据"或"候选"。

三、领域级陷阱

看一眼偏差（look-elsewhere）：在千万个方向、多个能量窗、多个时间窗里搜索，任何一点都可能"看起来显著"。必须对尝试次数做修正，或预先注册搜索方案。

太早看数据：分析流程如果是在看过数据之后定的（选择能量范围、切割参数），显著性会虚高。大实验的纪律是"盲分析"——先冻结流程再揭盲。

仪器效应伪装物理：像素坏点、增益漂移、卫星姿态误差都可能造成"点源"或"谱线"。任何反常先查仪器，再谈新物理——2014 年 AMS 曾有半 Electroweak 谱线传言，最终归因于分析处理；历史上"发现"磁单极、第五力的多次撤回都循此例。

图：从计数到发现的流水线

💡 关键直觉：在高能天体这个"低统计、高本底、多参数"的领域，分析方法本身就是仪器的一部分——流程冻结的时刻，决定了结果的可信度上限。

本节要点回顾

• 四步流水线：重建 → 本底 → 似然比 → 显著性，所有实验共用。
• Li-Ma 公式是点源搜索的经典工具，5σ 是发现的行规。
• 看一眼偏差是多窗口搜索的头号陷阱，修正是义务不是美德。
• 盲分析（先冻结流程后看数据）是领域自我保护的关键纪律。

动手估算：5σ 的代价与试用看

高能天文的"发现"门槛是 5σ，但这不是免费保险：在一个扫描了 10^7 个时空-能量格点的实验里，本底涨落超过 5σ 的期望次数恰好是 10^7×2.9×10^-7 ≈ 3。所以 IceCube 的事后试用看（trials factor）校正会把名义 5σ 的热点压回 3-4σ——统计推断在这里直接决定"发现了一个源"还是"发现了一个候选"。

# 试用看校正的数值感受
from math import erf, sqrt
def _inv_erf(x):          # 简易数值反演
    lo, hi = -5, 5
    for _ in range(60):
        mid = (lo+hi)/2
        if erf(mid/sqrt(2)) < x: lo = mid
        else: hi = mid
    return (lo+hi)/2
N = 1e7
for s_nominal in (4.0, 5.0):
    p = 2*(1-0.5*(1+erf(s_nominal/sqrt(2))))
    s_global = _inv_erf(1 - 2*(N*p))
    print(f"名义 {s_nominal}σ -> 全局 {s_global:.2f}σ")
# 名义 5σ 在千万次试用下只剩约 4.1σ
​

观测数据档案：统计方法的实战分工

# 用似然比检验比较"无源"与"有源"假设
import math
mu_bg, n_obs = 33.0, 48     # 模拟本底期望与观测计数
for mu_sig in (5, 15, 25):
    # 泊松对数似然比 2*ln L1/L0
    llr = 2*(mu_bg+mu_sig - n_obs + n_obs*math.log(n_obs/(mu_bg+mu_sig) or 1))
    print(f"假设信号 {mu_sig} 事: TS = {llr:.2f}")
# TS 大于约 25 才对应 5sigma 单点显著性 -- 弱信号天文的日常
​