Ходунки не "ходят" в модели подходят с помощью MCMC

Question

У меня возникли трудности с выполнением 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)

Заметьте, что я использую пакет emcee из пакета 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()

Answer 1

Я решил проблему! Просто напишите здесь решение для дальнейшего использования, если кто-то в будущем столкнется с подобной ситуацией. Ходоки не гуляли, потому что у меня был -инф. Решение на самом деле действительно исправительно и не имеет ничего общего с моим кодированием.

Я изолировал его до lnprior, потому что он выводил -infs. Оказывается, параметр p отрицателен. Поэтому условие:

p_ml / 100.

никогда не происходит. Это утверждение следует читать

p_ml/100. > p > p_ml*100.

Как только это редактирование выполнено, вышеуказанный код работает. Довольно глупый недосмотр!