2. 测量数据直方图与噪声分布拟合排障

为什么直方图拟合能区分“统计噪声”与“非统计异常”

重复测量得到的计数或概率数据，如果噪声来源只有量子投影噪声、探测器散粒噪声等统计成分，那么大量重复结果的频数分布应当收敛到某个理论分布：计数型数据通常接近泊松分布，均值较大（一般 λ≳20）时可用高斯近似；而归一化概率、差分信号等连续量更常用高斯分布描述。排障的关键不是“拟合得像不像”，而是看残差结构：纯统计噪声的残差应在零附近随机散布、无明显形状；一旦出现周期性条纹、单侧堆积、长尾或孤立尖峰，就说明混入了非统计成分（如周期性电磁干扰、基线漂移、离群事件、死时间效应）。

排障流程：从数据到结论

1. 整理数据：把 N 次重复测量结果按分箱统计频数。分箱数可用 Sturges 规则 k ≈ 1 + log2(N)，或 Freedman–Diaconis 规则 k ≈ 2·IQR·N^(-1/3) 的倒数；分箱过少会掩盖结构，过多会放大统计涨落。
2. 选择候选模型：计数型数据先试泊松分布 P(k)=λ^k e^(-λ)/k!；连续量或大均值计数用高斯 N(μ,σ²)。用样本均值/方差估计参数，不要用“目测”调参。
3. 计算拟合优度：对每个分箱计算期望频数 E_i，要求 E_i ≥ 5（不足则合并相邻箱），再算 χ² = Σ(O_i−E_i)²/E_i，自由度 ν = k − 1 − p（p 为估计参数个数，泊松 p=1，高斯 p=2）。
4. 看残差图：标准化残差 r_i = (O_i−E_i)/√E_i，正常应在 ±2 内随机分布；出现系统性形状即判为非统计成分。

可复现的排障脚本（Python）

下面脚本对一组示例计数数据完成分箱、泊松/高斯拟合、χ² 检验与残差输出。示例数据为教学构造，用于演示判读方法，不代表真实实验测量。

import numpy as np
from scipy import stats

rng = np.random.default_rng(7)
# 教学示例：主体为泊松噪声(λ=50)，另叠加少量离群事件
clean = rng.poisson(lam=50, size=4000)
outliers = rng.poisson(lam=90, size=60)      # 非统计成分：离群事件
data = np.concatenate([clean, outliers])

# 1) 分箱
counts, edges = np.histogram(data, bins="auto")
centers = 0.5 * (edges[:-1] + edges[1:])
width = edges[1] - edges[0]
N = data.size

# 2) 候选模型：泊松（用样本均值估计 λ）
lam = data.mean()
exp_pois = stats.poisson.pmf(np.arange(edges[0], edges[-1] + 1), lam) * N
# 按分箱边界聚合期望频数
exp_bins = np.array([
    exp_pois[max(0, int(np.floor(edges[i])) - int(edges[0])):
             max(1, int(np.ceil(edges[i+1])) - int(edges[0]))].sum()
    for i in range(len(counts))
])

# 3) 合并期望频数 < 5 的箱
def merge(counts, exp, min_exp=5):
    c, e = [], []
    acc_c = acc_e = 0.0
    for ci, ei in zip(counts, exp):
        acc_c += ci; acc_e += ei
        if acc_e >= min_exp:
            c.append(acc_c); e.append(acc_e)
            acc_c = acc_e = 0.0
    if acc_e > 0:
        c[-1] += acc_c; e[-1] += acc_e
    return np.array(c), np.array(e)

obs, exp = merge(counts, exp_bins)
chi2 = ((obs - exp) ** 2 / exp).sum()
dof = len(obs) - 1 - 1          # 泊松估计了 1 个参数 λ
pval = stats.chi2.sf(chi2, dof)

# 4) 标准化残差
resid = (obs - exp) / np.sqrt(exp)
print(f"λ={lam:.2f}  χ²={chi2:.2f}  ν={dof}  χ²/ν={chi2/dof:.2f}  p={pval:.3g}")
print("标准化残差:", np.round(resid, 2))
print("|r|>2 的箱数:", int((np.abs(resid) > 2).sum()))

常见误判与对策