PHOEBE 中如何自动设置 MCMC 的 lnprobability cutoff？

2026-11-04数据小许

结论：不要用 lnprobability 中位数筛选后验样本

如果目标是估计参数及其后验不确定性，就不应把 np.percentile(lnprob_clean, 50) 作为通用的样本保留阈值。这样做通常会删除约一半的有限样本，其中可能包含正常的后验尾部，从而使区间偏窄。

可以先按迭代位置去掉尚未稳定的初始阶段（burn-in），再确认后续采样是否充分，并保留有效的后验样本。如果采样仍然卡住，或没有充分探索参数空间，应继续采样或重新初始化，再决定是否报告结果。

仅凭 lnprobability 数值，无法设定一个适用于所有食双星系统并保证不确定性正确的自动 cutoff。

为什么这个阈值会截掉正常样本？

MCMC 的参数估计依赖样本在参数空间中的分布。较低的 lnprobability 并不自动意味着样本无效；采样器正常探索后验尾部时，也会产生较低的值。

假设这里的 lnprobability 是采样目标的对数密度，记作 (\ell(\theta))。如果只保留 (\ell(\theta)\ge c) 的样本，得到的分布就变成了：

[ p_c(\theta)\propto p(\theta),\mathbf{1}{\ell(\theta)\ge c}. ]

这相当于截断了原目标分布。此后计算的标准差或分位数区间描述的是截断后的分布，不能直接当作完整后验的不确定性。

flatten() 本身只是改变数组形状，但展平后会失去迭代位置和 walker 的区分。仅凭展平后的数值分布，无法判断某个低值属于初始化阶段，还是稳定采样中的正常波动。

批量处理应自动判断什么？

建议保留迭代和 walker 两个维度，并让每个系统经过以下流程：

检查初始化阶段。同时查看各参数和 lnprobability 的轨迹，寻找持续漂移、长时间停滞或不同区域之间缺少跳转的情况。lnprobability 看起来稳定，并不足以说明所有参数都已稳定。
确定候选 burn-in。在若干候选截断位置之后，比较前后时间窗口的参数分位数和区间宽度。窗口长度和允许变化量应由业务精度要求及代表性系统的诊断结果决定，没有通用数值。
评估剩余采样是否充分。对各参数估计 integrated autocorrelation time，检查估计值是否随采样延长而稳定，以及剩余链是否足够长。emcee 官方将自相关时间作为评估采样效率和 Monte Carlo 误差的重要工具，同时指出这些诊断不能保证收敛。emcee 自相关与收敛说明
将诊断失败的系统转入继续采样或人工检查。不要通过提高 cutoff，把没有采好的结果筛成一个看起来集中的分布。

自动化流程最好返回“采样充分”“需要继续采样”“需要检查初始化或多峰结构”等状态，而不是为每个系统强行生成一个阈值。

代码示例：按 burn-in 截取，再计算参数区间

下面的代码处理已经从 PHOEBE 提取的数组，不假定某个版本的参数样本读取接口。使用前需要确认并整理轴顺序：

chain：(n_steps, n_walkers, n_params)。
log_prob：(n_steps, n_walkers)，与 chain 逐点对应。
burnin：通过前述诊断确定的初始迭代数。

import numpy as np
from emcee.autocorr import integrated_time


def summarize_mcmc(chain, log_prob, burnin):
    chain = np.asarray(chain, dtype=float)
    log_prob = np.asarray(log_prob, dtype=float)

    if chain.ndim != 3:
        raise ValueError("chain must have shape (steps, walkers, parameters)")
    if log_prob.shape != chain.shape[:2]:
        raise ValueError("log_prob must match the first two axes of chain")
    if not isinstance(burnin, (int, np.integer)):
        raise ValueError("burnin must be an integer")
    if not 0 <= burnin < chain.shape[0]:
        raise ValueError("burnin is outside the available step range")

    post_chain = chain[burnin:]
    post_log_prob = log_prob[burnin:]

    # 保留时间结构进行诊断，不先展平或逐点删除样本。
    if not np.all(np.isfinite(post_chain)):
        raise ValueError("Non-finite parameter values require investigation")
    if not np.all(np.isfinite(post_log_prob)):
        raise ValueError("Non-finite log probabilities require investigation")

    # 链太短时可能抛出 AutocorrError，应转入继续采样流程。
    tau = integrated_time(post_chain, tol=50, quiet=False)
    if not np.all(np.isfinite(tau)) or np.any(tau <= 0):
        raise ValueError("Invalid autocorrelation time estimate")

    # 诊断之后再展平；不按 log probability 大小过滤。
    samples = post_chain.reshape(-1, post_chain.shape[-1])
    q16, q50, q84 = np.percentile(samples, [16, 50, 84], axis=0)

    return {
        "median": q50,
        "lower_error": q50 - q16,
        "upper_error": q84 - q50,
        "tau": tau,
        "steps_per_tau": post_chain.shape[0] / tau,
    }

integrated_time 接受上述三维数组结构；tol=50 表示估计所需的最小自相关时间倍数，链太短时会报告问题。通过这个检查也不等于证明收敛，仍需结合轨迹、时间窗口比较和多峰检查。emcee 自相关 API

这里用第 16、50、84 百分位数报告中位数和约 68% 的等尾可信区间，不要求后验服从高斯分布。

PHOEBE 中还需要注意什么？

具体如何关闭或绕过 lnprob_cutoff，需要依据所用 PHOEBE 版本和实际后处理接口确认。无论采用哪个入口，都要检查绘图、导出样本和计算参数区间时是否仍有隐含的阈值过滤。

如果某个 walker 长期处于低概率区域，不应只凭数值较低就删除。需要判断它是初始化失败、尚未混合，还是停留在一个真实的次要模式；这几种情况的处理方式不同。

同一 emcee ensemble 内的 walkers 不是相互独立的多条链，不能直接把它们当作独立链计算 Gelman-Rubin 诊断。需要这类交叉检查时，应使用独立初始化的多个 ensemble。emcee 收敛诊断注意事项