Skip to content
Snippets Groups Projects
Commit 52baf97c authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Update from log likelihood to log posterior

parent 3b03a766
No related branches found
No related tags found
1 merge request!842Add a mean-log-likelihood method to improve the ACT estimation
Hide whitespace changes
Inline Side-by-side
bilby/core/sampler/ptemcee.py
+ 48
− 39
View file @ 52baf97c
......@@ -26,7 +26,7 @@ ConvergenceInputs = namedtuple(
"autocorr_tol",
"autocorr_tau",
"gradient_tau",
"gradient_mean_logl",
"gradient_mean_log_posterior",
"Q_tol",
"safety",
"burn_in_nact",
......@@ -78,7 +78,7 @@ class Ptemcee(MCMCSampler):
gradient_tau: float, (0.1)
The maximum (smoothed) local gradient of the ACT estimate to allow.
This ensures the ACT estimate is stable before finishing sampling.
gradient_mean_logl: float, (0.1)
gradient_mean_log_posterior: float, (0.1)
The maximum (smoothed) local gradient of the logliklilhood to allow.
This ensures the ACT estimate is stable before finishing sampling.
Q_tol: float (1.01)
......@@ -159,7 +159,7 @@ class Ptemcee(MCMCSampler):
safety=1,
autocorr_tau=1,
gradient_tau=0.1,
gradient_mean_logl=0.1,
gradient_mean_log_posterior=0.1,
Q_tol=1.02,
min_tau=1,
check_point_deltaT=600,
......@@ -213,7 +213,7 @@ class Ptemcee(MCMCSampler):
mean_logl_frac=mean_logl_frac,
thin_by_nact=thin_by_nact,
gradient_tau=gradient_tau,
gradient_mean_logl=gradient_mean_logl,
gradient_mean_log_posterior=gradient_mean_log_posterior,
Q_tol=Q_tol,
nsamples=nsamples,
ignore_keys_for_tau=ignore_keys_for_tau,
......@@ -382,6 +382,7 @@ class Ptemcee(MCMCSampler):
self.iteration = data["iteration"]
self.chain_array = data["chain_array"]
self.log_likelihood_array = data["log_likelihood_array"]
self.log_posterior_array = data["log_posterior_array"]
self.pos0 = data["pos0"]
self.beta_list = data["beta_list"]
self.sampler._betas = np.array(self.beta_list[-1])
......@@ -424,7 +425,8 @@ class Ptemcee(MCMCSampler):
# Initialize storing results
self.iteration = 0
self.chain_array = self.get_zero_chain_array()
self.log_likelihood_array = self.get_zero_log_likelihood_array()
self.log_likelihood_array = self.get_zero_array()
self.log_posterior_array = self.get_zero_array()
self.beta_list = []
self.tau_list = []
self.tau_list_n = []
......@@ -437,7 +439,7 @@ class Ptemcee(MCMCSampler):
def get_zero_chain_array(self):
return np.zeros((self.nwalkers, self.max_steps, self.ndim))
def get_zero_log_likelihood_array(self):
def get_zero_array(self):
return np.zeros((self.ntemps, self.nwalkers, self.max_steps))
def get_pos0(self):
......@@ -483,14 +485,18 @@ class Ptemcee(MCMCSampler):
self.chain_array = np.concatenate((
self.chain_array, self.get_zero_chain_array()), axis=1)
self.log_likelihood_array = np.concatenate((
self.log_likelihood_array, self.get_zero_log_likelihood_array()),
self.log_likelihood_array, self.get_zero_array()),
axis=2)
self.log_posterior_array = np.concatenate((
self.log_posterior_array, self.get_zero_array()),
axis=2)
self.pos0 = pos0
self.chain_array[:, self.iteration, :] = pos0[0, :, :]
self.log_likelihood_array[:, :, self.iteration] = log_likelihood
self.mean_log_likelihood = np.mean(
self.log_likelihood_array[:, :, :self. iteration], axis=1)
self.log_posterior_array[:, :, self.iteration] = log_posterior
self.mean_log_posterior = np.mean(
self.log_posterior_array[:, :, :self. iteration], axis=1)
# Calculate time per iteration
self.time_per_check.append((datetime.datetime.now() - t0).total_seconds())
......@@ -499,15 +505,16 @@ class Ptemcee(MCMCSampler):
self.iteration += 1
# Calculate minimum iteration step to discard
mean_logl_min_it = get_mean_logl_min_it(
self.mean_log_likelihood,
minimum_iteration = get_minimum_stable_itertion(
self.mean_log_posterior,
frac=self.convergence_inputs.mean_logl_frac
)
logger.debug("Mean logl min it = {}".format(mean_logl_min_it))
logger.debug("Minimum iteration = {}".format(minimum_iteration))
# Calculate the maximum discard number
discard_max = np.max(
[self.convergence_inputs.burn_in_fixed_discard, mean_logl_min_it]
[self.convergence_inputs.burn_in_fixed_discard,
minimum_iteration]
)
if self.iteration > discard_max + self.nwalkers:
......@@ -536,7 +543,7 @@ class Ptemcee(MCMCSampler):
self.tau_list,
self.tau_list_n,
self.Q_list,
self.mean_log_likelihood,
self.mean_log_posterior,
)
if stop:
......@@ -650,8 +657,8 @@ class Ptemcee(MCMCSampler):
logger.info("tau plot failed with exception {}".format(e))
try:
plot_mean_log_likelihood(
self.mean_log_likelihood,
plot_mean_log_posterior(
self.mean_log_posterior,
self.outdir,
self.label,
)
......@@ -659,13 +666,13 @@ class Ptemcee(MCMCSampler):
logger.info("mean_logl plot failed with exception {}".format(e))
def get_mean_logl_min_it(mean_log_likelihood, frac):
nsteps = mean_log_likelihood.shape[1]
if nsteps < 10:
def get_minimum_stable_itertion(mean_array, frac, nsteps_min=10):
nsteps = mean_array.shape[1]
if nsteps < nsteps_min:
return 0
min_it = 0
for x in mean_log_likelihood:
for x in mean_array:
maxl = np.max(x)
fracdiff = (maxl - x) / np.abs(maxl)
idxs = fracdiff < frac
......@@ -685,7 +692,7 @@ def check_iteration(
tau_list,
tau_list_n,
Q_list,
mean_log_likelihood,
mean_log_posterior,
):
""" Per-iteration logic to calculate the convergence check
......@@ -776,28 +783,28 @@ def check_iteration(
)
tau_usable = False
check_mean_log_likelihood = mean_log_likelihood[:, -nsteps_to_check:]
gradient_mean_logl = get_max_gradient(
check_mean_log_likelihood, axis=1, window_length=GRAD_WINDOW_LENGTH,
check_mean_log_posterior = mean_log_posterior[:, -nsteps_to_check:]
gradient_mean_log_posterior = get_max_gradient(
check_mean_log_posterior, axis=1, window_length=GRAD_WINDOW_LENGTH,
smooth=True)
if gradient_mean_logl < ci.gradient_mean_logl:
if gradient_mean_log_posterior < ci.gradient_mean_log_posterior:
logger.debug(
"tau usable as {} < gradient_mean_logl={}"
.format(gradient_mean_logl, ci.gradient_mean_logl)
"tau usable as {} < gradient_mean_log_posterior={}"
.format(gradient_mean_log_posterior, ci.gradient_mean_log_posterior)
)
tau_usable = True
else:
logger.debug(
"tau not usable as {} > gradient_mean_logl={}"
.format(gradient_mean_logl, ci.gradient_mean_logl)
"tau not usable as {} > gradient_mean_log_posterior={}"
.format(gradient_mean_log_posterior, ci.gradient_mean_log_posterior)
)
tau_usable = False
else:
logger.debug("ACT is nan")
gradient_tau = np.nan
gradient_mean_logl = np.nan
gradient_mean_log_posterior = np.nan
tau_usable = False
if nsteps < tau_int * ci.autocorr_tol:
......@@ -816,7 +823,7 @@ def check_iteration(
samples_per_check,
tau_int,
gradient_tau,
gradient_mean_logl,
gradient_mean_log_posterior,
tau_usable,
convergence_inputs,
Q
......@@ -856,7 +863,7 @@ def print_progress(
samples_per_check,
tau_int,
gradient_tau,
gradient_mean_logl,
gradient_mean_log_posterior,
tau_usable,
convergence_inputs,
Q,
......@@ -880,7 +887,9 @@ def print_progress(
sampling_time = datetime.timedelta(seconds=np.sum(time_per_check))
tau_str = "{}(+{:0.2f},+{:0.2f})".format(tau_int, gradient_tau, gradient_mean_logl)
tau_str = "{}(+{:0.2f},+{:0.2f})".format(
tau_int, gradient_tau, gradient_mean_log_posterior
)
if tau_usable:
tau_str = "={}".format(tau_str)
......@@ -1050,21 +1059,21 @@ def plot_tau(
plt.close(fig)
def plot_mean_log_likelihood(mean_log_likelihood, outdir, label):
def plot_mean_log_posterior(mean_log_posterior, outdir, label):
ntemps, nsteps = mean_log_likelihood.shape
ymax = np.max(mean_log_likelihood)
ymin = np.min(mean_log_likelihood[:, -100:])
ntemps, nsteps = mean_log_posterior.shape
ymax = np.max(mean_log_posterior)
ymin = np.min(mean_log_posterior[:, -100:])
ymax += 0.1 * (ymax - ymin)
ymin -= 0.1 * (ymax - ymin)
fig, ax = plt.subplots()
idxs = np.arange(nsteps)
ax.plot(idxs, mean_log_likelihood.T)
ax.plot(idxs, mean_log_posterior.T)
ax.set(xlabel="Iteration", ylabel=r"$\langle\log\mathcal{L}\rangle$",
ylim=(ymin, ymax))
fig.tight_layout()
fig.savefig("{}/{}_checkpoint_meanloglike.png".format(outdir, label))
fig.savefig("{}/{}_checkpoint_meanlogposterior.png".format(outdir, label))
plt.close(fig)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment