我在执行模型的MCMC分析方面很难。我相信这与我在模型中具有不完整的伽马功能有关。
我试图最大程度地减少高斯对数的样式，但似乎步行者陷入了井中，而不是试图最大程度地减少可能性功能。这在下图中证明了这一点，其中y轴是模型的参数，x轴是步长。该图显示了步行者如何没有探索参数空间。我添加了另一个图像，以演示参数空间的正确探索。
不正确探索参数空间和参数空间的正确探索
我在下面添加了一些代码，以演示我在做什么，其中x，y和yerr是大约4000点的数组。该代码适用于其他模型，但仅在此模型上断开，因此它必须是其他人没有的固有的。其他模型最明显的变化是添加不完整的伽马函数，否则它与其他模型具有非常相似的功能形式。
我拟合的模型具有：
def singinhomobremss(freq,S_norm,alpha,p,freq_peak): # Single inhomogeneous free-free emission model
return S_norm*(p+1)*(np.power((freq/freq_peak),(2.1*(p+1)-alpha))*special.gammainc((p+1),(freq/freq_peak)**(-2.1)))*special.gamma(p+1)
请注意，我正在使用Python软件包（我会发布链接，但据说我没有足够的声誉...）。我真的不明白为什么步行者为其他模型拒绝为此模型拒绝"步行"。任何帮助都非常感谢，但我知道这是合理的利基区域。
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
import scipy.special as special # For access to the incomplete gamma function.
import emcee
import triangle
import inspect
def singinhomobremss(freq,S_norm,alpha,p,freq_peak): # Single inhomogeneous free-free emission model
return S_norm*(p+1)*(np.power((freq/freq_peak),(2.1*(p+1)-alpha))*special.gammainc((p+1),(freq/freq_peak)**(-2.1)))*special.gamma(p+1)
# initial guess for a fit
p0guess = [7, 0.6, -0.5, 3.5]
# Defining log-likelihood function
def lnlike(theta,x,y,yerr):
S_norm,alpha,p,freq_peak = theta
model = singinhomobremss(x,S_norm,alpha,p,freq_peak)
inv_sigma = 1.0/(yerr**2)
return -0.5*(np.sum((y-model)**2*inv_sigma - np.log(inv_sigma)))
# Use the scipy.opt model to find the optimum of this likelihood function
nll = lambda *args: -lnlike(*args)
result = opt.fmin(nll,p0guess,args=(x,y,yerr),full_output='true')
S_norm_ml,alpha_ml,p_ml,freq_peak_ml = result[0]
# Defining priors
def lnprior(theta):
S_norm,alpha,p,freq_peak = theta
if S_norm_ml/100. < S_norm < S_norm_ml/100. and alpha_ml/100. < alpha < alpha_ml*100. and p_ml/100. < p < p_ml*100. and freq_peak_ml/100. < freq_peak < freq_peak_ml*100:
return 0.00
return -np.inf
# Combining this prior with the definition of the likelihood function, the probablity fucntion is:
def lnprob(theta, x, y, yerr):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, x, y, yerr)
# Now implement emcee
ndim, nwalkers, nsteps = len(inspect.getargspec(singinhomobremss)[0])-1, 200, 2000
# Initialising the walkers in a Gaussian ball around maximum likelihood result
#pos = [result['x'] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
pos = [result[0] + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args = (x,y,yerr))
sampler.run_mcmc(pos, nsteps) # This is the workhorse step.
# Now plotting the walks of the walkers with each step, for each parameter. If they have converged on a good value they should have clumped together.
fig = plt.figure(2,figsize=(10, 10))
fig.clf()
for j in range(ndim):
ax = fig.add_subplot(ndim,1,j+1)
ax.plot(np.array([sampler.chain[:,i,j] for i in range(nsteps)]),"k", alpha = 0.3)
ax.set_ylabel((r'$S_{norm}$',r'$alpha$',r'$p$',r'$nu_{peak}$')[j], fontsize = 15)
plt.xlabel('Steps', fontsize = 15)
fig.show()
# To me it looks like the burn in period is well and truly over by 400 steps. So I will exclude those.
print 'The burnin applied was 400. Make sure the walkers have converged after that many steps.'
samples = sampler.chain[:,400:,:].reshape((-1,ndim))
# Plotting the histograms of the fit.
trifig = triangle.corner(samples, labels = [r'$S_{norm}$',r'$alpha$',r'$p$',r'$nu_{peak}$'])
# Finally to get the 1 sigma final uncertainties you do
S_norm_singinhomobremss_mcmc, alpha_singinhomobremss_mcmc, p_singinhomobremss_mcmc, freq_peak_singinhomobremss_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples,[16,50,84], axis = 0))) # Uncertainites based on the 16th, 50th and 84th percentile.
plt.figure()
plt.clf()
plt.errorbar(nu, flux, flux_err, marker = '.', color = 'gray', linestyle='none', label = 'Data',alpha=0.2)
plt.loglog(nu_peaklong,singinhomobremss(nu_peaklong, *poptsinginhomobremss), 'saddlebrown',label="Best fit inhomogeneous model from least-square")
plt.loglog(nu_peaklong,singinhomobremss(nu_peaklong,S_norm_singinhomobremss_mcmc[0], alpha_singinhomobremss_mcmc[0], p_singinhomobremss_mcmc[0], freq_peak_singinhomobremss_mcmc[0]),color = 'r', linestyle='-', label="Best fit inhomogeneous free-free model.")
plt.title('PKS 1718-649 Epoch - '+filename, fontsize = 15)
minnu = np.array(min(nu))-0.05*np.array(min(nu))
plt.legend(loc='lower center', fontsize=10) # make a legend in the best location
plt.xlabel('Frequency (GHz)', fontsize = 15)
plt.ylabel('Flux Density (Jy)', fontsize = 15)
plt.axis([min(nu)-0.1*min(nu), max(nu)+0.1*max(nu), min(flux)-0.1*min(flux), max(flux)+0.1*max(flux)])
plt.rc('xtick',labelsize=15)
plt.rc('ytick',labelsize=15)
plt.xticks([2,4,6,8,10]) #Setting grid line positions.
plt.yticks([3,4,5])
plt.grid(True)
plt.show()
我解决了问题！如果将来其他人偶然发现了类似情况，请在此处编写解决方案以备将来参考。步行者没有走路，因为lnprob是-inf。该解决方案实际上是补救措施，与我的编码无关。
我将其隔离到lnprior，因为它正在输出-infs。事实证明，参数P为负。因此，条件：
p_ml/100。＆lt;P＆LT;P_ML*100。
永远不会发生。此语句应阅读
p_ml/100。> p> p_ml*100。
将上述代码起作用后。非常愚蠢的疏忽！