主持人拒绝正确探索功能

Question

我有一个一维函数，如下所示：

即明显下降到 0 附近的稳定值。该函数写为：

((1. / np.sqrt(1. + x ** 2)) - (1. / np.sqrt(1. + C ** 2))) ** 2

其中

C

是我尝试使用 emcee 探索的参数。问题在于主持人（使用统一的先验）拒绝探索大可能性区域，而是看似随机地在该参数的整个允许范围内徘徊。痕迹如下（完整代码如下）：

其中真实值用红线显示。虽然

scipy.optimize.minimize

能够轻松放大真实值，但主持人显然无法做到这一点。

我做错了什么，还是这个函数无法像我一样使用统一的先验来探索？

---

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()

Answer 1

在对数似然函数中，在分母中包含噪声标准差和因子 2，即将其更改为：

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

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

并且您的采样器运行到：

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y_obs, C_max, noise))

你应该得到你所期望的：