如何为PHOEBE中MCMC自动选择lnprobability截断值？

直接用中位数（50分位数）确实容易过度截断，因为MCMC链的lnprob分布通常是右偏的——大部分样本集中在较低lnprob区域，而有效样本（收敛后的）集中在高lnprob的峰值附近。以下几种自动方法可以解决这个问题：

1. 结合burn-in判断的峰值区间截断 #

先剔除未收敛的burn-in阶段样本，再针对收敛后的样本计算截断值，能避免低概率的初始样本干扰：

• 用emcee的自相关时间计算确定burn-in长度，一般取最大自相关时间的2-3倍
• 对burn-in后的样本，选择90或95分位数作为截断值，既能过滤噪声样本，又不会过度压缩有效样本的分布

示例代码：

import emcee
import numpy as np

# 获取PHOBE输出的lnprobability（形状为nwalkers x nsteps）
lnprob = sol.get_value('lnprobabilities@fastcompute@emcee_solver@emcee_sol@emcee@solution')

# 过滤非有限值，同时保留walker维度
finite_mask = np.isfinite(lnprob)
lnprob_clean = lnprob[finite_mask].reshape(lnprob.shape[0], -1)

# 计算自相关时间，确定burn-in长度
tau = emcee.autocorr.integrated_time(lnprob_clean, quiet=True)
burn_in = int(2 * np.max(tau))

# 保留burn-in后的样本
lnprob_post_burn = lnprob_clean[:, burn_in:]
lnprob_flat_post = lnprob_post_burn.flatten()

# 取90分位数作为截断值
ln_cutoff = np.percentile(lnprob_flat_post, 90)

2. 基于KDE峰值的自适应截断 #

用核密度估计（KDE）拟合lnprob的分布，找到峰值后，取峰值减去1-2倍样本标准差作为截断值，能自适应匹配每个系统的lnprob分布特征：

from scipy.stats import gaussian_kde

# 先完成burn-in处理（同方法1）
lnprob_flat_post = lnprob_post_burn.flatten()

# 拟合KDE曲线
kde = gaussian_kde(lnprob_flat_post)
xvals = np.linspace(lnprob_flat_post.min(), lnprob_flat_post.max(), 1000)
pdf = kde(xvals)

# 找到lnprob的峰值位置
peak_lnprob = xvals[np.argmax(pdf)]
# 计算样本标准差
lnprob_std = np.std(lnprob_flat_post)

# 截断值设为峰值减1倍标准差（可根据实际情况调整为1.5倍）
ln_cutoff = peak_lnprob - 1 * lnprob_std

3. 基于有效样本占比的自适应截断 #

设定一个目标保留的高概率样本占比（比如30%-50%），直接找到对应的lnprob分位数，确保不会过度截断：

# 完成burn-in处理后
lnprob_flat_post = lnprob_post_burn.flatten()

# 设定保留30%的高概率样本
target_fraction = 0.3
# 计算对应分位数（分位数从小到大排序，所以取100*(1-target_fraction)）
ln_cutoff = np.percentile(lnprob_flat_post, 100*(1-target_fraction))

关键注意事项 #

• 先处理burn-in是核心：未收敛的初始样本lnprob普遍偏低，会严重干扰截断值的合理性
• 可针对部分系统做验证：对比截断前后的参数不确定性区间，确保截断后样本能覆盖真实的参数分布
• 若食双星系统存在多模态参数分布，需先通过高斯混合模型等方法检测模态，再针对每个模态单独计算截断值

内容的提问来源于stack exchange，提问作者Valentina Bonilla