求指导：泊松生成数据集的单参数拟合参数反推

数据生成逻辑 #

以下是生成探测器活跃数量频率分布的代码：

import numpy as np

def generator(photons):
    runs = 10000
    sectors = []
    events = np.random.poisson(photons, runs)
    for i in events:
        hit = np.random.randint(low=1, high=9, size=i)
        sectors.append(hit)
    freq = np.zeros((9,), dtype=int)
    for j in sectors:
        act_sec = len(np.unique(j))
        freq[act_sec] += 1
    freq = [ k / runs for k in freq ]

    return freq

该函数逻辑：

• 按泊松分布生成runs次事件的光子数
• 每个光子随机击中1-8号探测器中的一个
• 统计每次事件中活跃探测器数量（即被击中过的不同探测器数）的频率分布，返回长度为9的数组（索引对应活跃探测器数0-8，其中0的频率始终为0）

需求与问题 #

需要通过观测到的活跃探测器频率分布，反推生成数据时的输入参数photons（入射光子数）。要求：

• 拟合时以探测器数量为截断条件（即只使用部分活跃探测器数对应的频率数据）
• 拟合模型仅允许单个参数（即目标参数photons）

尝试过scipy.optimize.curve_fit和lmfit工具，但拟合未成功，寻求具体实现指导。

解决思路与实现步骤 #

1. 推导理论频率模型

拟合的核心是建立活跃探测器数k（k=1~8）与入射光子数μ（即photons参数）的理论概率关系：
对于单个事件，入射光子数服从泊松分布P(n) = e^(-μ) * μ^n / n!；当有n个光子时，用容斥原理计算恰好激活k个探测器的概率：

P(k | n) = C(8, k) * [ (k/8)^n - C(k, k-1)*((k-1)/8)^n + ... + (-1)^(k-1)*C(k,1)*(1/8)^n ]

最终活跃探测器数k的边缘概率为：

P(k | μ) = Σ(n=k到∞) P(n) * P(k | n)

实际计算时，只需求和到μ+5*sqrt(μ)即可覆盖泊松分布99.9%的概率区间。

2. 构建拟合目标函数

定义输入参数μ、输出对应活跃探测器数理论频率的函数：

import numpy as np
from scipy.special import comb, poisson

def theoretical_freq(μ, max_n=None):
    if max_n is None:
        # 取μ+5倍标准差作为最大n，覆盖绝大多数概率
        max_n = int(μ + 5 * np.sqrt(μ)) + 2
    freq = np.zeros(9)
    for k in range(1, 9):
        prob_k = 0.0
        for n in range(k, max_n+1):
            # 泊松概率P(n)
            p_n = poisson.pmf(n, μ)
            # 容斥计算n个光子恰好击中k个探测器的概率
            p_k_n = 0.0
            for i in range(k):
                sign = (-1)**i
                c = comb(k, i)
                term = c * ((k - i)/8)**n
                p_k_n += sign * term
            prob_k += p_n * p_k_n
        freq[k] = prob_k
    return freq

3. 拟合实现（以scipy为例）

使用curve_fit时，需注意用卡方统计量作为拟合依据（更适合频率分布），示例代码：

from scipy.optimize import curve_fit

# 生成观测数据（实际使用时替换为真实观测数据）
observed_freq = generator(3)
# 提取k=1~8的有效频率数据
obs_data = observed_freq[1:]

# 定义拟合用函数：输入μ，输出k=1~8的理论频率
def fit_func(μ):
    return theoretical_freq(μ)[1:]

# 初始猜测值（通过活跃探测器数均值估算：E[k]≈8*(1 - e^(-μ/8))，反推得μ≈-8*ln(1 - E[k]/8)）
initial_guess = [4.0]

# 执行拟合，用频率的标准误作为权重
popt, pcov = curve_fit(fit_func, xdata=[], ydata=obs_data, p0=initial_guess, 
                       sigma=np.sqrt(obs_data/10000),  # runs=10000，标准误为sqrt(freq/runs)
                       absolute_sigma=True)

print(f"拟合得到的光子数μ={popt[0]:.3f}，标准差={np.sqrt(pcov[0][0]):.3f}")

4. 截断条件的处理

如果需要只使用部分活跃探测器数的数据（比如k=2~7），只需修改数据提取和拟合函数：

# 截断到k=2~7的频率数据
obs_data = observed_freq[2:8]

def fit_func_truncated(μ):
    return theoretical_freq(μ)[2:8]

# 重新执行拟合
popt_trunc, pcov_trunc = curve_fit(fit_func_truncated, xdata=[], ydata=obs_data, p0=initial_guess,
                                   sigma=np.sqrt(obs_data/10000), absolute_sigma=True)

5. 拟合失败常见原因排查

• 初始猜测偏差过大：用活跃探测器数均值估算初始μ，避免拟合陷入局部最优
• 理论模型与生成逻辑不一致：检查容斥原理的实现是否正确，或是否遗漏了足够多的n取值
• 权重设置不合理：必须用频率的标准误作为权重，否则低频率数据会干扰拟合结果

内容的提问来源于stack exchange，提问作者No Telling