def MCMC_run(x,y, uncert, params, pmin, pmax, positions = None, num_iter = 10000, return_ = False, name = "mcmc"):
"""_summary_
Args:
x (_type_): _description_
y (_type_): _description_
uncert (_type_): _description_
params (_type_): _description_
pmin (_type_): _description_
pmax (_type_): _description_
positions (_type_): _description_
num_iter (int, optional): _description_. Defaults to 10000.
return_ (bool, optional): _description_. Defaults to False.
name (str, optional): _description_. Defaults to "mcmc".
Returns:
_type_: _description_
"""
if positions is None:
positions = np.linspace(0,len(y)-1,len(y))
names = ["ax","ay", "bx", "by", "ww", "m", "b", "phi"]
funcs = logfuncs(pmin,pmax)
from scipy.optimize import minimize
nll = lambda *args: -1*funcs.log_likelihood(*args)
initial = np.array(params) + 0.1 * np.random.randn(8)
soln = minimize(nll, initial, args=(x, y, uncert, positions))
print(soln.x)
import emcee
from multiprocessing import Pool
pos = np.asarray(soln.x) + 1e-5 * np.random.randn(17, 8)
nwalkers, ndim = pos.shape
filename = str(name) + ".h5"
backend = emcee.backends.HDFBackend(filename)
backend.reset(nwalkers, ndim)
with Pool() as pool:
sampler = emcee.EnsembleSampler(
nwalkers, ndim, funcs.log_probability, args=(x, y, uncert, positions), backend = backend, pool = pool
)
sampler.run_mcmc(pos, num_iter, progress=True)
#, moves=[(emcee.moves.DEMove(), 0.8),(emcee.moves.DESnookerMove(), 0.2)]
samples = sampler.get_chain()
import pickle
with open("emcee_res.pkl", "wb") as f:
pickle.dump(samples, f)
if return_ is True:
return sampler