求助：emcee实现MCMC拟合对数正态函数未获合理结果

有人能帮我找出代码里的问题吗？我用emcee做MCMC，把模拟数据拟合到对数正态函数，但拟合结果不如curve-fit的合理。翻了好多网上示例都没找到错在哪，麻烦帮忙看看，谢谢！

以下是我的代码：

import matplotlib.pyplot as plt
import numpy as np
import emcee
from scipy.optimize import curve_fit
import corner

def lognormal(x, A, mu, sigma):
     return A * 1/(x * sigma * np.sqrt(2 * np.pi)) * np.exp(-(np.log(x) - mu)**2/(2 * sigma**2))

#MCMC

def mcmc(x_data, y_data, func, initial_guess = [1,1,1], burning_in = 200, thinner = 50, steps = 5000, walkers = 100, dimension = 3):
#Scipy curve fit

    curve_fit_params, cov= curve_fit(func, x_data, y_data, p0 = initial_guess)

    curve_fit_error = np.sqrt(np.diag(cov))


    #MCMCM setup

    # Define the log-likelihood function

    def log_likelihood(theta, x, y, func):
        A, mu, sigma = theta
        y_model = func(x, A, mu, sigma)
        residual = y - y_model
        chi_squared = np.sum(residual**2)
        log_like = -0.5 * chi_squared
        return log_like


    # Define the log-prior function
    def log_prior(theta):
        A, mu, sigma = theta
        # Add any prior constraints here
        # For example, uniform priors for A, mu and sigma between certain ranges
        if 0 < A < 2.*curve_fit_params[0] and 0 < mu < 2.*curve_fit_params[1] and 0 < sigma < 2.*curve_fit_params[2]:
            return 0.0  # log(1) = 0
        return -np.inf  # log(0) = -inf

    # Define the log-posterior function
    def log_posterior(theta, x, y, func):
        log_prior_val = log_prior(theta)
        if not np.isfinite(log_prior_val):
            return -np.inf  # log(0) = -inf
        log_like_val = log_likelihood(theta, x, y, func)
        return log_prior_val + log_like_val

    #MCMC

    # Set up the number of walkers and dimensions
    nwalkers = walkers
    ndim = dimension

    # Initialize the walkers with random positions near the maximum likelihood solution
    pos = [0.1 * np.random.randn(ndim) + curve_fit_params for _ in range(nwalkers)]

    # Set up the emcee sampler
    sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x_data, y_data, func))

    # Run the MCMC sampling
    nsteps = steps  # Number of steps to run the sampler
    sampler.run_mcmc(pos, nsteps, progress = True)

    #tau = sampler.get_autocorr_time()
    #print(tau)

    # Get the chain of samples and flatten it
    samples = sampler.get_chain(discard=burning_in, thin = thinner, flat=True)

    # Get the best-fit parameters (maximum a posteriori or MAP estimate)
    mcmc_fit_mean = np.mean(samples, axis=0)

    # Get the uncertainties in the parameters (credible intervals)
    mcmc_fit_std = np.std(samples, axis=0)

    #corner plot

    labels = ["A", "mu", "sigma"]

    fig = corner.corner(samples, bins = 100, plot_datapoints=False, smooth = True, labels=labels, truths=curve_fit_params, truth_color = 'darkred', color = '#002448', quantiles=[0.16, 0.50, 0.84], hist_kwargs={"density": True, "alpha": 0.4, "histtype": "stepfilled"})
    #plt.savefig('MCMC_test_convergence.png', bbox_inches='tight', dpi=800)
    plt.show()
    plt.clf()
    
    return samples, mcmc_fit_mean, mcmc_fit_std, curve_fit_params, curve_fit_error



#test data

n_sample = 1000

x_test = np.linspace(0.001,250, n_sample)

y_test = lognormal(x_test, 0.8, 4.7, 0.4) + np.random.normal(0, 0.0003, n_sample)


plt.scatter(x_test, y_test, linewidth  = 0.5, s = 1, color = 'darkred')
plt.show()
plt.clf()



result = mcmc(x_test, y_test, lognormal)
print('MCMC results: ', *result[1])
print('curvefit results: ', *result[3])

plt.plot(x_test, lognormal(x_test, *result[1]), linewidth = 0.5, color = 'navy', label = 'MCMC')
plt.plot(x_test, lognormal(x_test, *result[3]), linewidth = 0.5, color = 'orange', label = 'curvefit')
plt.scatter(x_test, y_test, linewidth  = 0.5, s = 1, color = 'darkred')
plt.legend(loc = 'best')
plt.clf()

问题分析与修复方案 #

1. 对数似然函数缺失噪声归一化 #

当前log_likelihood直接计算残差平方和，但未考虑数据噪声的幅度，导致似然函数缩放不合理，MCMC采样的后验分布偏离最优解。

修复后的似然函数（加入噪声标准差估计）：

def log_likelihood(theta, x, y, func):
    A, mu, sigma = theta
    y_model = func(x, A, mu, sigma)
    # 估计噪声标准差
    sigma_noise = np.std(y - y_model)
    residual = y - y_model
    # 归一化残差的卡方统计量
    chi_squared = np.sum((residual / sigma_noise)**2)
    # 加入似然的归一化常数（不影响采样，但更符合统计定义）
    log_like = -0.5 * (chi_squared + len(y) * np.log(2 * np.pi * sigma_noise**2))
    return log_like

2. 先验范围过于狭窄 #

原代码将参数限制在0到2*curve_fit_params之间，若curve-fit结果有偏差，MCMC会被限制在局部区域无法探索全局最优解。建议基于真实数据范围放宽先验：

def log_prior(theta):
    A, mu, sigma = theta
    # 基于模拟数据的真实参数设置合理宽范围
    if 0.1 < A < 2.0 and 2.0 < mu < 6.0 and 0.1 < sigma < 1.0:
        return 0.0
    return -np.inf

3. 采样参数设置不合理 #

• 未验证自相关时间：打开tau = sampler.get_autocorr_time()的注释，确保thinner设置为自相关时间的2-3倍，避免样本自相关。
• 有效采样数不足：可适当增加steps到10000，减少thinner（比如取tau的2倍），保证有效样本量足够。

4. 最优参数选择用中位数而非均值 #

对数正态参数的后验分布通常不对称，中位数（corner图中的0.5分位数）比均值更能代表最优拟合结果：

# 替换原mcmc_fit_mean的计算
mcmc_fit_median = np.percentile(samples, 50, axis=0)

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