主持人步行者烧毁，但随后保持不变

Question

我在使用 emcee 时遇到问题。它是一个足够简单的 3 参数拟合，但偶尔（尽管使用很多，但仅在两种情况下发生）我的助行器烧得很好，但然后不动（见下图）。报告的接受分数为 0。

以前有没有其他人遇到过这个问题？我尝试改变我的初始条件以及步行者和迭代的数量等。这段代码在非常相似的数据集上运行良好。它不是一个具有挑战性的参数空间，步行者似乎不太可能“卡住”。

有什么想法吗？我被难住了……我的步行者似乎很懒……

下面的示例代码（和sample data file）。这个代码+数据文件在我运行时失败。

`
import numpy as n
import math
import pylab as py    
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
from scipy import ndimage
import pyfits
from scipy import stats
import emcee
import corner
import scipy.optimize as op
import matplotlib.pyplot as pl
from matplotlib.ticker import MaxNLocator

def sersic(x, B,r_s,m):
      return B * n.exp(-1.0 * (1.9992*m - 0.3271) * ( (x/r_s)**(1.0/m) - 1.))

def lnprior(theta):
      B,r_s,m, lnf = theta
      if  0.0 < B < 500.0 and 0.5 < m < 10. and r_s > 0. and -10.0 < lnf < 1.0: 
        return 0.0
      return -n.inf

def lnlike(theta, x, y, yerr): #"least squares"
      B,r_s,m, lnf = theta
      model = sersic(x,B, r_s, m)
      inv_sigma2 = 1.0/(yerr**2 + model**2*n.exp(2*lnf))
      return -0.5*(n.sum((y-model)**2*inv_sigma2 - n.log(inv_sigma2)))

def lnprob(theta, x, y, yerr):#kills based on priors
    lp = lnprior(theta)
    if not n.isfinite(lp):
        return -n.inf
    return lp + lnlike(theta, x, y, yerr)


profile=open("testprofile.dat",'r')  #read in the data file    
profilelines=profile.readlines()
x=n.empty(len(profilelines))
y=n.empty(len(profilelines))
yerr=n.empty(len(profilelines))

for i,line in enumerate(profilelines):
    col=line.split()
    x[i]=col[0]
    y[i]=col[1]
    yerr[i]=col[2]

# Find the maximum likelihood value.
chi2 = lambda *args: -2 * lnlike(*args)
result = op.minimize(chi2, [50,2.0,0.5,0.5], args=(x, y, yerr))
B_ml, rs_ml,m_ml, lnf_ml = result["x"]
print("""Maximum likelihood result:
          B   = {0}
          r_s = {1}
          m   = {2}
    """.format(B_ml, rs_ml,m_ml))

# Set up the sampler.
ndim, nwalkers = 4, 4000
pos = [result["x"] + 1e-4*n.random.randn(ndim) for i in     range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
# Clear and run the production chain.
print("Running MCMC...")
Niter = 2000 #2000
sampler.run_mcmc(pos, Niter, rstate0=n.random.get_state())
print("Done.")
# Print out the mean acceptance fraction. 
af = sampler.acceptance_fraction
print "Mean acceptance fraction:", n.mean(af)

# Plot sampler chain
pl.clf()
fig, axes = pl.subplots(3, 1, sharex=True, figsize=(8, 9))
axes[0].plot(sampler.chain[:, :, 0].T, color="k", alpha=0.4)
axes[0].yaxis.set_major_locator(MaxNLocator(5))
axes[0].set_ylabel("$B$")

axes[1].plot(sampler.chain[:, :, 1].T, color="k", alpha=0.4)
axes[1].yaxis.set_major_locator(MaxNLocator(5))
axes[1].set_ylabel("$r_s$")

axes[2].plot(n.exp(sampler.chain[:, :, 2]).T, color="k", alpha=0.4)
axes[2].yaxis.set_major_locator(MaxNLocator(5))
axes[2].set_xlabel("step number")

fig.tight_layout(h_pad=0.0)
fig.savefig("line-time_test.png")

# plot MCMC fit
burnin = 100
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim-1))

B_mcmc, r_s_mcmc, m_mcmc = map(lambda v: (v[0]),
                         zip(*n.percentile(samples, [50],
                                            axis=0)))
print("""MCMC  result:
         B   = {0}
         r_s = {1}
         m   = {2}
    """.format(B_mcmc, r_s_mcmc, m_mcmc))

pl.close()

# Make the triangle plot.
burnin = 50
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim-1))

fig = corner.corner(samples, labels=["$B$", "$r_s$", "$m$"])
fig.savefig("line-triangle_test.png")

Answer 1

这是一个更好的结果。我使随机初始样本不太接近最大似然值，并以更少的步行者/链运行链更多的步骤。请注意，我正在绘制 m 参数而不是它的指数，就像你所做的那样。

平均接受分数约为 0.48，在我的笔记本电脑上运行大约需要 1 分钟。您当然可以添加更多步骤并获得更好的结果。

import numpy as n
import emcee
import corner
import scipy.optimize as op
import matplotlib.pyplot as pl
from matplotlib.ticker import MaxNLocator


def sersic(x, B, r_s, m):
    return B * n.exp(
        -1.0 * (1.9992 * m - 0.3271) * ((x / r_s)**(1.0 / m) - 1.))


def lnprior(theta):
    B, r_s, m, lnf = theta
    if 0.0 < B < 500.0 and 0.5 < m < 10. and r_s > 0. and -10.0 < lnf < 1.0:
        return 0.0
    return -n.inf


def lnlike(theta, x, y, yerr):  # "least squares"
    B, r_s, m, lnf = theta
    model = sersic(x, B, r_s, m)
    inv_sigma2 = 1.0 / (yerr**2 + model**2 * n.exp(2 * lnf))
    return -0.5 * (n.sum((y - model)**2 * inv_sigma2 - n.log(inv_sigma2)))


def lnprob(theta, x, y, yerr):  # kills based on priors
    lp = lnprior(theta)
    if not n.isfinite(lp):
        return -n.inf
    return lp + lnlike(theta, x, y, yerr)


profile = open("testprofile.dat", 'r')  # read in the data file
profilelines = profile.readlines()
x = n.empty(len(profilelines))
y = n.empty(len(profilelines))
yerr = n.empty(len(profilelines))

for i, line in enumerate(profilelines):
    col = line.split()
    x[i] = col[0]
    y[i] = col[1]
    yerr[i] = col[2]

# Find the maximum likelihood value.
chi2 = lambda *args: -2 * lnlike(*args)
result = op.minimize(chi2, [50, 2.0, 0.5, 0.5], args=(x, y, yerr))
B_ml, rs_ml, m_ml, lnf_ml = result["x"]
print("""Maximum likelihood result:
          B   = {0}
          r_s = {1}
          m   = {2}
          lnf = {3}
    """.format(B_ml, rs_ml, m_ml, lnf_ml))

# Set up the sampler.
ndim, nwalkers = 4, 10
pos = [result["x"] + 1e-1 * n.random.randn(ndim) for i in range(nwalkers)]

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(x, y, yerr))
# Clear and run the production chain.
print("Running MCMC...")
Niter = 50000
sampler.run_mcmc(pos, Niter, rstate0=n.random.get_state())
print("Done.")
# Print out the mean acceptance fraction.
af = sampler.acceptance_fraction
print("Mean acceptance fraction:", n.mean(af))

# Plot sampler chain
pl.clf()
fig, axes = pl.subplots(3, 1, sharex=True, figsize=(8, 9))
axes[0].plot(sampler.chain[:, :, 0].T, color="k", alpha=0.4)
axes[0].yaxis.set_major_locator(MaxNLocator(5))
axes[0].set_ylabel("$B$")

axes[1].plot(sampler.chain[:, :, 1].T, color="k", alpha=0.4)
axes[1].yaxis.set_major_locator(MaxNLocator(5))
axes[1].set_ylabel("$r_s$")

# axes[2].plot(n.exp(sampler.chain[:, :, 2]).T, color="k", alpha=0.4)
axes[2].plot(sampler.chain[:, :, 2].T, color="k", alpha=0.4)
axes[2].yaxis.set_major_locator(MaxNLocator(5))
axes[2].set_ylabel("$m$")
axes[2].set_xlabel("step number")

fig.tight_layout(h_pad=0.0)
fig.savefig("line-time_test.png")

# plot MCMC fit
burnin = 10000
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim - 1))

B_mcmc, r_s_mcmc, m_mcmc = map(
    lambda v: (v[0]), zip(*n.percentile(samples, [50], axis=0)))
print("""MCMC  result:
         B   = {0}
         r_s = {1}
         m   = {2}
    """.format(B_mcmc, r_s_mcmc, m_mcmc))

pl.close()

# Make the triangle plot.
burnin = 50
samples = sampler.chain[:, burnin:, :3].reshape((-1, ndim - 1))

fig = corner.corner(samples, labels=["$B$", "$r_s$", "$m$"])
fig.savefig("line-triangle_test.png")