I am writing a mcmc python code using the python package emcee. The structure goes as follows `
### MCMC Parameters
# initial guesses for the free parameters
initial = np.float64([1.e-13, -2])
# set no of walkers
nwalkers = 100
# set no of steps
nsteps = 3000
# defining the model
def model(theta, ..):
# here I am calling 1 external codes using subprocess
# writes an input file depending on the freeparameters
write_numdens(theta, rrr)
# calling the first external code
line = "radmc3d image imolspec 1 iline 1 incl {} posang {} widthkms {} vkms {} nodust nostar".format(inc,PA,6,6.41)
subprocess.run([line],shell=True)
# output of the radmc3d is being used by vis_sample
model_vis = vis_sample(imagefile='image.out', uvfile=uvfile, src_distance=dist, verbose=True)
# returns 2d array with complex values
return model_vis
def lnlike(theta,rrr,y_obs):
y_model = model(theta, rrr)
LnLike = -0.5*np.sum(weights*(((y_obs.real-y_model.real))**2 + ((y_obs.imag-y_model.imag))**2))
return LnLike
### prior
def lnprior(theta):
X_0, alpha = theta
# make sure priors are reasonable for source
if 1e-20 < X_0 < 1e-8 and -3. < alpha < 2.:
return 0.0
return -np.inf
### log posterior prob
def lnprob(theta, rrr, y_obs):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, rrr, y_obs) #recall if lp not -inf, its 0, so this just returns likelihood
### The MCMC run
# A list of extra positional arguments which will be called for lnprob along with theta
data = (rrr, y_obs)
# no of parameters
ndim = len(initial)
# this is a multivariate Gaussian centered on each theta, with a small sigma.
p0 = [np.array(initial) + 1e-3 * np.random.randn(ndim) for i in range(nwalkers)]
def main(p0,nwalkers,nsteps,ndim,lnprob,data):
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=data)
print("Running burn-in...")
p0, _, _ = sampler.run_mcmc(p0, 100, progress=True)
sampler.reset()
print("Running production...")
pos, prob, state = sampler.run_mcmc(p0, nsteps, progress=True)
return sampler, pos, prob, state
sampler, pos, prob, state = main(p0, nwalkers, nsteps, ndim, lnprob, data)
The problem I am facing is that while running the code, it does not seem to access the model. Even when I am not giving the necessary inputs required to run the model, the code is running successfully. Anything I am printing inside model is not getting printed. The code is supposed to take days to finish, whereas it is finishing in seconds. Do you have any suggestion to solve this situation?
I have written mcmc codes with similar code structure before, but this is the first time I am facing this kind problem. I have no clue what to do.