MCMC를 사용하여 워커가 모델에 맞게 걸을 수 없음

Question

모델의 MCMC 분석을 수행하는 데 어려움이 있습니다. 나는 그것이 모델 내에서 불완전한 감마 함수를 가지고 있다는 사실과 관련이 있다고 생각한다.

가우시안 로그 가능성을 최소화하려고하지만 워커가 우물에 머물러 있고 우도 함수를 최소화하려고 시도하지 않는 것으로 보입니다. 이것은 아래 그림에서 보여집니다. 여기서 y 축은 모델의 매개 변수이고 x 축은 단계 번호입니다. 이미지는 보행자가 매개 변수 공간을 탐색하지 않는 방법을 보여줍니다. 나는 매개 변수 공간의 올바른 탐색이 어떻게 보이는지 보여주기 위해 또 다른 이미지를 추가했다.

Incorrect exploration of parameter space 내가 X, Y 및 yerr 약 4000 점의 배열 어디에 내가 뭐하는 거지 보여주기 위해 아래의 몇 가지 코드를 추가 한

Correct exploration of parameter space. 이 코드는 다른 모델에서는 작동하지만이 모델에서만 작동하므로 다른 모델에는없는 본질적인 것이어야합니다. 다른 모델에 대한 가장 분명한 변화는 불완전한 감마 함수의 추가입니다. 그렇지 않으면 다른 모델과 매우 유사한 기능적 형태를 갖습니다. 나는 (... 내가 링크를 게시 것이지만, 아마도 내가 충분히 명성을 가지고 있지 않음) 파이썬 패키지 사회자를 사용하고

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) 

참고 : 나는 피팅하고

모델은 형태를 갖는다. 나는 다른 모델을 위해 할 때이 모델을 보행자가 '걸을'것을 거부하는 이유를 실제로 알지 못합니다. 어떤 도움이 많이 감사하지만 그것이 합리적으로 틈새 영역임을 이해합니다.

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

나는이 문제를 해결했습니다! 미래에 다른 사람이 비슷한 상황에 걸리면 나중에 참조 할 수 있도록 여기에 해결책을 작성하십시오. lnprob이 -inf 이었기 때문에 도보 자들이 걷지 않고 있었다. 이 솔루션은 실제로 실제로 개선되고 코딩과 관련이 없습니다.

-infs를 출력했기 때문에 lnprior로 격리했습니다. 매개 변수 p가 음수라는 것이 밝혀졌습니다. 따라서, 조건 :

p_ml/100. <p <p_ml * 100.

이 발생하지 않습니다. 이 문장은

p_ml/100으로 표시되어야합니다. > p> p_ml * 100.

편집이 완료되면 위의 코드가 작동합니다. 꽤 어리석은 감시!