emcee无法有效探索目标函数参数的问题排查求助

我有一个一维函数，形式如下：
((1. / np.sqrt(1. + x ** 2)) - (1. / np.sqrt(1. + C ** 2))) ** 2
该函数在0附近有明显下降并趋于稳定值。我尝试使用emcee对其中的参数C进行探索，但采用均匀先验时，emcee并未探索高似然区域，反而在参数允许范围内随机游走，其采样轨迹如图所示。而scipy.optimize.minimize却能轻松定位到C的真实值。请问是我的操作有误，还是该函数不适合使用均匀先验进行探索？

附完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import emcee


def main():

    # Set true value for the variable
    C_true = 27.

    # Generate synthetic data
    x = np.arange(.1, 100)
    y_true = func(x, C_true)
    noise = 0.01
    y_obs = np.random.normal(y_true, noise)

    # Set up the MCMC
    nwalkers = 4
    ndim = 1
    nburn = 500
    nsteps = 5000
    # Maximum value for the 'C' parameter
    C_max = 5 * C_true
    # Use a 10% STDDEV around the true value for the initial state
    p0 = [np.random.normal(C_true, C_true * .1, nwalkers)]
    p0 = np.array(p0).T

    # Run the MCMC
    print("Running emcee...")
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y_obs, C_max))
    # Burn-in
    state = sampler.run_mcmc(p0, nburn)
    sampler.reset()
    sampler.run_mcmc(state, nsteps)
    samples = sampler.chain.reshape((-1, ndim))

    # Print the median and 1-sigma uncertainty of the parameters
    C_median = np.median(samples)
    C_percnt = np.percentile(samples, [16, 84])
    print(f'C = {C_median:.2f} ({C_percnt[0]:.2f}, {C_percnt[1]:.2f})')

    # Chains
    plt.plot(sampler.chain[:, :, 0].T, c='k', alpha=0.1)
    plt.axhline(C_true, color='r')
    plt.ylabel('C')
    plt.xlabel('Step')
    plt.tight_layout()
    plt.show()

    # Fitted func
    plt.scatter(x, y_obs)
    y_emcee = func(x, C_median)
    plt.scatter(x, y_emcee)
    plt.show()


def func(x, C):
    x_C = ((1. / np.sqrt(1. + x ** 2)) - (1. / np.sqrt(1. + C ** 2))) ** 2
    # Beyond C, the function is fixed to 0
    return np.where(x < C, x_C, 0)


def lnlike(C, x, y_obs):
    model = func(x, C)
    lkl = -np.sum((y_obs - model) ** 2)
    return lkl


def lnprior(C, C_max):
    if 0 < C < C_max:
        return 0.0
    return -np.inf


def lnprob(C, x, y_obs, C_max):
    lp = lnprior(C, C_max)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(C, x, y_obs)


if __name__ == '__main__':
    main()

问题原因分析 #

• 似然函数平坦性：当C大于数据中最大x值（99）时，模型输出固定为0，和C=99的输出完全一致，导致似然函数在C>99区域极度平坦。emcee作为MCMC采样器，在无梯度引导的平坦区域会随机游走；而scipy.minimize是局部优化器，能快速锁定真实值所在的高似然峰值区域。
• 先验范围过大：设置的C_max=135包含了似然无差异的平坦区域，采样器容易陷入这些无信息区域。
• 步行者数量不足：nwalkers=4对于一维问题来说偏少，采样效率低，难以脱离非高似然区域。

修正方案 #

1. 缩小先验范围 #

将C_max设置为数据的最大x值，因为C>=x.max()时模型输出无差异，无需采样：

C_max = x.max()  # 替换原5*C_true

2. 修正似然函数 #

添加噪声项的高斯似然更符合统计模型，能让似然峰值更尖锐，帮助采样器定位：

def lnlike(C, x, y_obs, noise):
    model = func(x, C)
    # 正确的高斯对数似然形式
    return -0.5 * np.sum(((y_obs - model)/noise)**2 + np.log(2*np.pi*noise**2))

同时更新lnprob和sampler初始化：

def lnprob(C, x, y_obs, C_max, noise):
    lp = lnprior(C, C_max)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(C, x, y_obs, noise)

# 初始化sampler时传入noise参数
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y_obs, C_max, noise))

3. 增加步行者数量 #

将nwalkers提升至10（一维问题建议至少2*ndim，这里ndim=1，更多步行者能提升采样稳定性）：

nwalkers = 10

4. 优化初始值（可选） #

可以让初始值覆盖更宽的合理范围，比如在(0.5*C_true, 1.5*C_true)内随机生成，避免初始值过于集中：

p0 = np.random.uniform(0.5*C_true, 1.5*C_true, (nwalkers, ndim))

效果验证 #

调整后，emcee的采样链会收敛到真实值附近，不再出现大范围随机游走，采样结果的中位数和置信区间会更接近真实值。

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