Skip to content

Commit

Permalink
Merge pull request #728 from pints-team/i719-log-likelihood-tests
Browse files Browse the repository at this point in the history 
Value-based tests for log-likelihoods
  • Loading branch information
@MichaelClerx
MichaelClerx authored Feb 22, 2019
2 parents c1218dd + 0206af7 commit 9c73c55
Show file tree
Hide file tree
Showing 5 changed files with 217 additions and 60 deletions.
55 changes: 34 additions & 21 deletions pints/_log_likelihoods.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,19 @@ class AR1LogLikelihood(pints.ProblemLogLikelihood):
Calculates a log-likelihood assuming AR1 noise model
.. math::
\log{L(\\theta, \sigma|\\boldsymbol{x})} =
\log{L(\\theta, \sigma'|\\boldsymbol{x})} =
-\\frac{N}{2}\log{2\pi}
-N\log{\sigma}
-\\frac{1}{2\sigma^2}
\sum_{i=1}^N{(\\epsilon_i x_i - \\rho \\epsilon_{i-1} )^2}
-N\log{\sigma'}
-\\frac{1}{2\sigma'^2}
\sum_{i=2}^N{(\\epsilon_i x_i - \\rho \\epsilon_{i-1} )^2}
where
.. math::
\\epsilon_i = x_i - f_i(\\theta)
and :math:`sigma' = \\frac{sigma} \\sqrt{1-\\rho^2}`.
Arguments:
Expand All @@ -54,8 +56,11 @@ def __init__(self, problem):
self._logn = 0.5 * (self._nt) * np.log(2 * np.pi)

def __call__(self, x):
rho = np.asarray(x[-2 * self._no:-self._no])
sigma = np.asarray(x[-self._no:]) * np.sqrt(1 - rho**2)
m = 2 * self._no
parameters = x[-m:]
rho = np.asarray(parameters[0::2])
sigma = np.asarray(parameters[1::2])
sigma = np.asarray(sigma) * np.sqrt(1 - rho**2)
error = self._values - self._problem.evaluate(x[:-2 * self._no])
autocorr_error = error[1:] - rho * error[:-1]
return np.sum(- self._logn - self._nt * np.log(sigma)
Expand All @@ -64,21 +69,29 @@ def __call__(self, x):

class ARMA11LogLikelihood(pints.ProblemLogLikelihood):
"""
Calculates a log-likelihood assuming AR1 noise model
Calculates a log-likelihood assuming ARMA(1,1) noise model.
.. math::
\log{L(\\theta, \sigma|\\boldsymbol{x})} =
-\\frac{N}{2}\log{2\pi}
-N\log{\sigma}
-\\frac{1}{2\sigma^2}
\sum_{i=1}^N{(\\epsilon_i x_i - \\rho \\epsilon_{i-1} -
\\phi \\nu(t-1))^2}
\sum_{i=3}^N{(\\nu_i - \\phi \\nu_{i-1})^2}
where
.. math::
\\nu_i = \\epsilon_i - \\rho \\epsilon_{i-1}
and
..math::
\\epsilon_i = x_i - f_i(\\theta)
and
.. math::
\\sigma = \\sigma\\sqrt{\\frac{1-\\rho^2}{1 + 2\\phi\\rho + \\phi^2}}`
Arguments:
Expand All @@ -104,13 +117,16 @@ def __init__(self, problem):
self._logn = 0.5 * (self._nt) * np.log(2 * np.pi)

def __call__(self, x):
rho = np.asarray(x[-3 * self._no:-2 * self._no])
phi = np.asarray(x[-2 * self._no:-self._no])
m = 3 * self._no
parameters = x[-m:]
rho = np.asarray(parameters[0::3])
phi = np.asarray(parameters[1::3])
sigma = np.asarray(parameters[2::3])
sigma = (
np.asarray(x[-self._no:]) *
sigma *
np.sqrt((1.0 - rho**2) / (1.0 + 2.0 * phi * rho + phi**2))
)
error = self._values - self._problem.evaluate(x[:-3 * self._no])
error = self._values - self._problem.evaluate(x[:-m])
v = error[1:] - rho * error[:-1]
autocorr_error = v[1:] - phi * v[:-1]
return np.sum(- self._logn - self._nt * np.log(sigma)
Expand Down Expand Up @@ -320,29 +336,27 @@ def __call__(self, x):

def evaluateS1(self, x):
""" See :meth:`LogPDF.evaluateS1()`. """
sigma = float(np.asarray(x[-self._no:]))
sigma = np.asarray(x[-self._no:])

# Evaluate, and get residuals
y, dy = self._problem.evaluateS1(x[:-self._no])

# Reshape dy, in case we're working with a single-output problem
dy = dy.reshape(self._nt, self._no, self._n_parameters - 1)
dy = dy.reshape(self._nt, self._no, self._n_parameters - self._no)

# Note: Must be (data - simulation), sign now matters!
r = self._values - y

# Calculate log-likelihood
L = np.sum(-self._logn - self._nt * np.log(sigma)
- (1.0 / (2 * sigma**2)) * np.sum(r**2, axis=0))
L = self.__call__(x)

# Calculate derivatives in the model parameters
dL = np.sum(
(sigma**(-2.0) * np.sum((r.T * dy.T).T, axis=0).T).T, axis=0)

# Calculate derivative wrt sigma
dsigma = np.sum(-self._nt / sigma +
sigma**(-3.0) * np.sum(r**2, axis=0))
dL = np.concatenate((dL, np.array([dsigma])))
dsigma = -self._nt / sigma + sigma**(-3.0) * np.sum(r**2, axis=0)
dL = np.concatenate((dL, np.array(list(dsigma))))

# Return
return L, dL
Expand Down Expand Up @@ -490,4 +504,3 @@ def __init__(self, problem):
'The class `pints.KnownNoiseLogLikelihood` is deprecated.'
' Please use `pints.GaussianLogLikelihood` instead.')
super(UnknownNoiseLogLikelihood, self).__init__(problem)

62 changes: 31 additions & 31 deletions pints/tests/test_error_measures.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,15 +154,15 @@ def test_mean_squared_error_multi(self):
x = [1, 2]

# Model outputs are 3 times [1, 4]
# Model derivatives are 3 times [[1, 0], [0, 2]]
# Model derivatives are 3 times [[1, 0], [0, 1]]
y, dy = p.evaluateS1(x)
self.assertTrue(np.all(y == p.evaluate(x)))
self.assertTrue(np.all(y[0, :] == [1, 4]))
self.assertTrue(np.all(y[1, :] == [1, 4]))
self.assertTrue(np.all(y[2, :] == [1, 4]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 1]]))

# Check residuals
rx = y - np.array(values)
Expand All @@ -183,11 +183,11 @@ def test_mean_squared_error_multi(self):
self.assertEqual(dex.shape, (2, ))

# Residuals are: [[0, 0], [-1, -3], [-2, -6]]
# Derivatives are: [[1, 0], [0, 2]]
# Derivatives are: [[1, 0], [0, 1]]
# dex1 is: (2 / nt / no) * (0 - 1 - 2) * 1 = (1 / 3) * -3 * 1 = -1
# dex2 is: (2 / nt / no) * (0 - 3 - 6) * 2 = (1 / 3) * -9 * 2 = -6
# dex2 is: (2 / nt / no) * (0 - 3 - 6) * 1 = (1 / 3) * -9 * 1 = -3
self.assertEqual(dex[0], -1)
self.assertEqual(dex[1], -6)
self.assertEqual(dex[1], -3)

def test_mean_squared_error_weighted(self):
""" Tests :class:`pints.MeanSquaredError` with weighted outputs. """
Expand Down Expand Up @@ -218,15 +218,15 @@ def test_mean_squared_error_weighted(self):
x = [1, 2]

# Model outputs are 3 times [1, 4]
# Model derivatives are 3 times [[1, 0], [0, 2]]
# Model derivatives are 3 times [[1, 0], [0, 1]]
y, dy = p.evaluateS1(x)
self.assertTrue(np.all(y == p.evaluate(x)))
self.assertTrue(np.all(y[0, :] == [1, 4]))
self.assertTrue(np.all(y[1, :] == [1, 4]))
self.assertTrue(np.all(y[2, :] == [1, 4]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 1]]))

# Check residuals
rx = y - np.array(values)
Expand All @@ -251,11 +251,11 @@ def test_mean_squared_error_weighted(self):
# dex1 is: (2 / nt / no) * (0 - 1 - 2) * 1 * 1
# = (1 / 3) * -3 * 1 * 1
# = -1
# dex2 is: (2 / nt / no) * (0 - 3 - 6) * 2 * 2
# = (1 / 3) * -9 * 2 * 2
# = -12
# dex2 is: (2 / nt / no) * (0 - 3 - 6) * 1 * 2
# = (1 / 3) * -9 * 1 * 2
# = -6
self.assertEqual(dex[0], -1)
self.assertEqual(dex[1], -12)
self.assertEqual(dex[1], -6)

def test_probability_based_error(self):
""" Tests :class:`pints.ProbabilityBasedError`. """
Expand Down Expand Up @@ -344,15 +344,15 @@ def test_sum_of_squares_error_multi(self):
x = [1, 2]

# Model outputs are 3 times [1,4]
# Model derivatives are 3 times [[1, 0], [0, 2]]
# Model derivatives are 3 times [[1, 0], [0, 1]]
y, dy = p.evaluateS1(x)
self.assertTrue(np.all(y == p.evaluate(x)))
self.assertTrue(np.all(y[0, :] == [1, 4]))
self.assertTrue(np.all(y[1, :] == [1, 4]))
self.assertTrue(np.all(y[2, :] == [1, 4]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 1]]))

# Check residuals
rx = y - np.array(values)
Expand All @@ -373,11 +373,11 @@ def test_sum_of_squares_error_multi(self):
self.assertEqual(dex.shape, (2, ))

# Residuals are: [[0, 0], [-1, -3], [-2, -6]]
# Derivatives are: [[1, 0], [0, 2]]
# Derivatives are: [[1, 0], [0, 1]]
# dex1 is: 2 * (0 - 1 - 2) * 1 = 2 * -3 * 1 = -6
# dex2 is: 2 * (0 - 3 - 6) * 2 = 2 * -9 * 2 = -36
# dex2 is: 2 * (0 - 3 - 6) * 2 = 2 * -9 * 1 = -18
self.assertEqual(dex[0], -6)
self.assertEqual(dex[1], -36)
self.assertEqual(dex[1], -18)

def test_sum_of_squares_error_weighted(self):
""" Tests :class:`pints.MeanSquaredError` with weighted outputs. """
Expand Down Expand Up @@ -408,15 +408,15 @@ def test_sum_of_squares_error_weighted(self):
x = [1, 2]

# Model outputs are 3 times [1, 4]
# Model derivatives are 3 times [[1, 0], [0, 2]]
# Model derivatives are 3 times [[1, 0], [0, 1]]
y, dy = p.evaluateS1(x)
self.assertTrue(np.all(y == p.evaluate(x)))
self.assertTrue(np.all(y[0, :] == [1, 4]))
self.assertTrue(np.all(y[1, :] == [1, 4]))
self.assertTrue(np.all(y[2, :] == [1, 4]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 2]]))
self.assertTrue(np.all(dy[0, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[1, :] == [[1, 0], [0, 1]]))
self.assertTrue(np.all(dy[2, :] == [[1, 0], [0, 1]]))

# Check residuals
rx = y - np.array(values)
Expand All @@ -437,15 +437,15 @@ def test_sum_of_squares_error_weighted(self):
self.assertEqual(dex.shape, (2, ))

# Residuals are: [[0, 0], [-1, -3], [-2, -6]]
# Derivatives are: [[1, 0], [0, 2]]
# Derivatives are: [[1, 0], [0, 1]]
# dex1 is: 2 * (0 - 1 - 2) * 1 * 1
# = 2 * -3 * 1 * 1
# = -6
# dex2 is: 2 * (0 - 3 - 6) * 2 * 2
# = 2 * -9 * 2 * 2
# = -72
# dex2 is: 2 * (0 - 3 - 6) * 2 * 1
# = 2 * -9 * 2 * 1
# = -36
self.assertEqual(dex[0], -6)
self.assertEqual(dex[1], -72)
self.assertEqual(dex[1], -36)

def test_sum_of_errors(self):
""" Tests :class:`pints.SumOfErrors`. """
Expand Down
Loading

0 comments on commit 9c73c55

Please sign in to comment.