Skip to content

Commit

Permalink
draft of log-normal function working
Browse files Browse the repository at this point in the history
partly addresses #1376
  • Loading branch information
ben18785 committed Aug 6, 2021
1 parent ce67499 commit e852fdc
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 10 deletions.
17 changes: 9 additions & 8 deletions pints/_log_likelihoods.py
Original file line number Diff line number Diff line change
Expand Up @@ -842,8 +842,7 @@ def __init__(self, problem):
# Pre-calculate parts
self._logn = 0.5 * self._nt * np.log(2 * np.pi)
vals = np.asarray(self._values)
print(not np.any(vals))
if not np.any(vals):
if np.any(vals <= 0):
raise ValueError('All data points must exceed zero.')
self._log_values = np.log(self._values)

Expand All @@ -852,13 +851,21 @@ def __call__(self, x):
if any(sigma < 0):
return -np.inf
soln = self._problem.evaluate(x[:-self._no])
if np.any(soln < 0):
return -np.inf
error = np.log(self._values) - np.log(soln)
return np.sum(- self._logn - self._nt * np.log(sigma)
- np.sum(self._log_values, axis=0)
- np.sum(error**2, axis=0) / (2 * sigma**2))

def evaluateS1(self, x):
""" See :meth:`LogPDF.evaluateS1()`. """
# Calculate log-likelihood and when log-prob is infinite, gradients are
# not defined
L = self.__call__(x)
if L == -np.inf:
return L, np.tile(np.nan, self._n_parameters)

sigma = np.asarray(x[-self._no:])

# Evaluate, and get residuals
Expand All @@ -870,12 +877,6 @@ def evaluateS1(self, x):
# Note: Must be (np.log(data) - simulation), sign now matters!
r = np.log(self._values) - np.log(y)

# Calculate log-likelihood and when log-prob is infinite, gradients are
# not defined
L = self.__call__(x)
if L == -np.inf:
return L, np.tile(np.nan, self._n_parameters)

# Calculate derivatives in the model parameters
dL = np.sum(
(sigma**(-2.0) * np.sum((r.T * dy.T / y.T).T, axis=0).T).T, axis=0)
Expand Down
12 changes: 10 additions & 2 deletions pints/tests/test_log_likelihoods.py
Original file line number Diff line number Diff line change
Expand Up @@ -1429,8 +1429,10 @@ def setUpClass(cls):
def test_bad_constructor(self):
# tests that bad data types result in error
data = [0, 4, 5.5, 7.2]
self.assertRaises(ValueError, pints.SingleOutputProblem,
self.model_single, self.times, data)
problem = pints.SingleOutputProblem(
self.model_single, self.times, data)
self.assertRaises(ValueError, pints.LogNormalLogLikelihood,
problem)

def test_call(self):
# test calls of log-likelihood
Expand Down Expand Up @@ -1482,6 +1484,12 @@ def test_evaluateS1(self):
self.assertEqual(y, -np.inf)
self.assertTrue(np.array_equal(dL, np.tile(np.nan, 2), equal_nan=True))

mu = -1
sigma = 1
y, dL = self.log_likelihood.evaluateS1([mu, sigma])
self.assertEqual(y, -np.inf)
self.assertTrue(np.array_equal(dL, np.tile(np.nan, 2), equal_nan=True))

# two dim output problem
mu1 = 1.5
mu2 = 3.4 / 2
Expand Down

0 comments on commit e852fdc

Please sign in to comment.