Commit 8d8c3853 authored by Gregory Ashton's avatar Gregory Ashton
Browse files

Update initial check to prior and posterior

Previously, a check was done to ensure that all initial positions of the
walkers would yield a non -np.inf prior. This extends this to also check
the likelihood evaluations. The method is also simplified partially

Partially resolves issues in #6
parent 6358d069
...@@ -25,9 +25,7 @@ theta_prior = {'F0': {'type': 'unif', ...@@ -25,9 +25,7 @@ theta_prior = {'F0': {'type': 'unif',
'F2': F2, 'F2': F2,
'Alpha': Alpha, 'Alpha': Alpha,
'Delta': Delta, 'Delta': Delta,
'transient_tstart': {'type': 'unif', 'transient_tstart': minStartTime,
'lower': minStartTime,
'upper': maxStartTime},
'transient_duration': {'type': 'halfnorm', 'transient_duration': {'type': 'halfnorm',
'loc': 0.001*Tspan, 'loc': 0.001*Tspan,
'scale': 0.5*Tspan} 'scale': 0.5*Tspan}
...@@ -212,38 +212,39 @@ class MCMCSearch(core.BaseSearchClass): ...@@ -212,38 +212,39 @@ class MCMCSearch(core.BaseSearchClass):
self.theta_symbols = [self.theta_symbols[i] for i in idxs] self.theta_symbols = [self.theta_symbols[i] for i in idxs]
self.theta_keys = [self.theta_keys[i] for i in idxs] self.theta_keys = [self.theta_keys[i] for i in idxs]
def _evaluate_logpost(self, p0vec):
init_logp = np.array([
self.logp(p, self.theta_prior, self.theta_keys,
for p in p0vec])
init_logl = np.array([
for p in p0vec])
return init_logl + init_logp
def _check_initial_points(self, p0): def _check_initial_points(self, p0):
for nt in range(self.ntemps): for nt in range(self.ntemps):'Checking temperature {} chains'.format(nt))'Checking temperature {} chains'.format(nt))
initial_priors = np.array([ num = sum(self._evaluate_logpost(p0[nt]) == -np.inf)
self.logp(p, self.theta_prior, self.theta_keys, if num > 0:
for p in p0[nt]])
number_of_initial_out_of_bounds = sum(initial_priors == -np.inf)
if number_of_initial_out_of_bounds > 0:
logging.warning( logging.warning(
'Of {} initial values, {} are -np.inf due to the prior' 'Of {} initial values, {} are -np.inf due to the prior'
.format(len(initial_priors), .format(len(p0[0]), num))
p0 = self._generate_new_p0_to_fix_initial_points( p0 = self._generate_new_p0_to_fix_initial_points(
p0, nt, initial_priors) p0, nt)
def _generate_new_p0_to_fix_initial_points(self, p0, nt, initial_priors): def _generate_new_p0_to_fix_initial_points(self, p0, nt):'Attempting to correct intial values')'Attempting to correct intial values')
idxs = np.arange(self.nwalkers)[initial_priors == -np.inf] init_logpost = self._evaluate_logpost(p0[nt])
idxs = np.arange(self.nwalkers)[init_logpost == -np.inf]
count = 0 count = 0
while sum(initial_priors == -np.inf) > 0 and count < 100: while sum(init_logpost == -np.inf) > 0 and count < 100:
for j in idxs: for j in idxs:
p0[nt][j] = (p0[nt][np.random.randint(0, self.nwalkers)]*( p0[nt][j] = (p0[nt][np.random.randint(0, self.nwalkers)]*(
1+np.random.normal(0, 1e-10, self.ndim))) 1+np.random.normal(0, 1e-10, self.ndim)))
initial_priors = np.array([ init_logpost = self._evaluate_logpost(p0[nt])
self.logp(p, self.theta_prior, self.theta_keys,
for p in p0[nt]])
count += 1 count += 1
if sum(initial_priors == -np.inf) > 0: if sum(init_logpost == -np.inf) > 0:'Failed to fix initial priors')'Failed to fix initial priors')
else: else:'Suceeded to fix initial priors')'Suceeded to fix initial priors')
Supports Markdown
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