Jacobians should not be computed only for transforms

This commit is contained in:
James Hensman 2015-06-03 21:23:27 +01:00
parent d67e37fcb6
commit c4752c1564

View file

@ -433,35 +433,34 @@ class Indexable(Nameable, Updateable):
if self.priors.size == 0: if self.priors.size == 0:
return 0. return 0.
x = self.param_array x = self.param_array
#evaluate the prior log densities #evaluate the prior log densities
log_p = reduce(lambda a, b: a + b, (p.lnpdf(x[ind]).sum() for p, ind in self.priors.items()), 0) log_p = reduce(lambda a, b: a + b, (p.lnpdf(x[ind]).sum() for p, ind in self.priors.items()), 0)
#account for the transformation by evaluating the log Jacobian (where things are transformed) #account for the transformation by evaluating the log Jacobian (where things are transformed)
log_j = 0. log_j = 0.
priored_indexes = np.hstack([i for p, i in self.priors.items()]) priored_indexes = np.hstack([i for p, i in self.priors.items()])
for c,j in self.constraints.items(): for c,j in self.constraints.items():
if c is 'fixed':continue if not isinstance(c, Transformation):continue
for jj in j: for jj in j:
if jj in priored_indexes: if jj in priored_indexes:
log_j += c.log_jacobian(x[jj]) log_j += c.log_jacobian(x[jj])
return log_p + log_j return log_p + log_j
def _log_prior_gradients(self): def _log_prior_gradients(self):
"""evaluate the gradients of the priors""" """evaluate the gradients of the priors"""
if self.priors.size == 0: if self.priors.size == 0:
return 0. return 0.
x = self.param_array x = self.param_array
ret = np.zeros(x.size) ret = np.zeros(x.size)
#compute derivate of prior density #compute derivate of prior density
[np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.items()] [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.items()]
#add in jacobian derivatives if transformed #add in jacobian derivatives if transformed
priored_indexes = np.hstack([i for p, i in self.priors.items()]) priored_indexes = np.hstack([i for p, i in self.priors.items()])
for c,j in self.constraints.items(): for c,j in self.constraints.items():
for jj in j: for jj in j:
if jj in priored_indexes: if jj in priored_indexes:
ret[jj] += c.log_jacobian_grad(x[jj]) ret[jj] += c.log_jacobian_grad(x[jj])
return ret return ret
#=========================================================================== #===========================================================================
# Tie parameters together # Tie parameters together