Repository: pints-team/pints (language: python)

Issue:
Create GaussianIntegratedLogUniformLogLikelihood
Similar to https://pints.readthedocs.io/en/stable/log_likelihoods.html#pints.GaussianIntegratedUniformLogLikelihood, but assume a prior that is uniform in log-transformed space instead of the original parameter space, i.e. `p(\sigma) \propto 1 / \sigma`, following Eq. (7)-(9) in https://doi.org/10.1093/gji/ggaa168.
I1483 gaussian integrated log uniform log likelihood
Closes #1483 


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ac3969350..eef1acb71 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@ All notable changes to this project will be documented in this file.
 - [#1479](https://github.com/pints-team/pints/pull/1479) PINTS is no longer tested on Python 3.6.
 
 ### Added
+- [#1484](https://github.com/pints-team/pints/pull/1484) Added a `GaussianIntegratedLogUniformLogLikelihood` class that calculates the log-likelihood with its Gaussian noise integrated out with an uninformative prior in log-space.
 - [#1466](https://github.com/pints-team/pints/pull/1466) Added a `TransformedRectangularBoundaries` class that preserves the `RectangularBoundaries` methods after transformation.
 - [#1462](https://github.com/pints-team/pints/pull/1461) The `OptimisationController` now has a stopping criterion `max_evaluations`.
 - [#1460](https://github.com/pints-team/pints/pull/1460) [#1468](https://github.com/pints-team/pints/pull/1468) Added the `Adam` local optimiser.
diff --git a/docs/source/log_likelihoods.rst b/docs/source/log_likelihoods.rst
index ce9903f04..d2d4dbe6b 100644
--- a/docs/source/log_likelihoods.rst
+++ b/docs/source/log_likelihoods.rst
@@ -20,6 +20,7 @@ Overview:
 - :class:`ARMA11LogLikelihood`
 - :class:`CauchyLogLikelihood`
 - :class:`ConstantAndMultiplicativeGaussianLogLikelihood`
+- :class:`GaussianIntegratedLogUniformLogLikelihood`
 - :class:`GaussianIntegratedUniformLogLikelihood`
 - :class:`GaussianKnownSigmaLogLikelihood`
 - :class:`GaussianLogLikelihood`
@@ -39,6 +40,8 @@ Overview:
 
 .. autoclass:: ConstantAndMultiplicativeGaussianLogLikelihood
 
+.. autoclass:: GaussianIntegratedLogUniformLogLikelihood
+
 .. autoclass:: GaussianIntegratedUniformLogLikelihood
 
 .. autoclass:: GaussianKnownSigmaLogLikelihood
diff --git a/pints/__init__.py b/pints/__init__.py
index 65bbe9268..e6448e64a 100644
--- a/pints/__init__.py
+++ b/pints/__init__.py
@@ -115,6 +115,7 @@ def version(formatted=False):
     ARMA11LogLikelihood,
     CauchyLogLikelihood,
     ConstantAndMultiplicativeGaussianLogLikelihood,
+    GaussianIntegratedLogUniformLogLikelihood,
     GaussianIntegratedUniformLogLikelihood,
     GaussianKnownSigmaLogLikelihood,
     GaussianLogLikelihood,
diff --git a/pints/_log_likelihoods.py b/pints/_log_likelihoods.py
index 49fb22bff..3c7792a7b 100644
--- a/pints/_log_likelihoods.py
+++ b/pints/_log_likelihoods.py
@@ -465,23 +465,106 @@ def evaluateS1(self, parameters):
         return L, dL
 
 
+class GaussianIntegratedLogUniformLogLikelihood(pints.ProblemLogLikelihood):
+    r"""
+    Calculates a log-likelihood assuming independent Gaussian-distributed noise
+    at each time point where :math:`p(\sigma)\propto 1/\sigma` has been
+    integrated out of the joint posterior of :math:`p(\theta,\sigma|X)`,
+
+    .. math::
+        \begin{align} p(\theta|X) &= \int_{0}^{\infty} p(\theta, \sigma|X)
+        \mathrm{d}\sigma\\
+        &\propto \int_{0}^{\infty} p(X|\theta, \sigma) p(\theta, \sigma)
+        \mathrm{d}\sigma.\end{align}
+
+    Note that this is exactly the same statistical model as
+    :class:`pints.GaussianLogLikelihood` with a uniform prior on
+    :math:`\sigma` in the logarithmic transformed space.
+    This is also known as Jeffrey's prior, as the standard deviation is a scale
+    parameter, as discussed here: https://doi.org/10.1093/gji/ggaa168.
+
+    A possible advantage of this log-likelihood compared with using a
+    :class:`pints.GaussianLogLikelihood`, is that it has one fewer parameters
+    (:math:`sigma`) which may speed up convergence to the posterior
+    distribution, especially for multi-output problems which will have
+    ``n_outputs`` fewer parameter dimensions.
+
+    The log-likelihood is given in terms of the sum of squared errors:
+
+    .. math::
+        SSE = \sum_{i=1}^n (f_i(\theta) - y_i)^2,
+
+    and is given up to a normalisation constant by:
+
+    .. math::
+        \begin{align}
+        \text{log} L =
+            & - n / 2 \text{log}(SSE) \\
+            & - (n - 2) / 2 \text{log}(2) \\
+            & + \text{log}\left[\Gamma(n / 2)\right],
+        \end{align}
+
+    where :math:`\Gamma(a)` is the incomplete gamma function.
+
+    This log-likelihood is inherently a Bayesian method since it assumes a
+    uniform prior on :math:`\sigma` _in the logarithmic transformed space.
+    However using this likelihood in optimisation routines should yield the
+    same estimates as the full :class:`pints.GaussianLogLikelihood`.
+
+    Extends :class:`ProblemLogLikelihood`.
+
+    Parameters
+    ----------
+    problem
+        A :class:`SingleOutputProblem` or :class:`MultiOutputProblem`.
+    """
+
+    def __init__(self, problem):
+        super(GaussianIntegratedLogUniformLogLikelihood,
+              self).__init__(problem)
+
+        # Get number of times, number of outputs
+        self._nt = len(self._times)
+        self._no = problem.n_outputs()
+
+        # Add parameters to problem
+        self._n_parameters = problem.n_parameters()
+
+        # Pre-calculate
+        n = self._nt
+        self._n_over_2 = n / 2
+        self._log_gamma = scipy.special.gammaln(self._n_over_2)
+        self._constant_1 = -np.log(2) * (n - 2) / 2
+        self._constant = self._constant_1 + self._log_gamma
+
+    def __call__(self, x):
+        error = self._values - self._problem.evaluate(x)
+        sse = np.sum(error**2, axis=0)
+
+        # Calculate
+        sse = pints.vector(sse)
+        return np.sum(self._constant - self._n_over_2 * np.log(sse))
+
+
 class GaussianIntegratedUniformLogLikelihood(pints.ProblemLogLikelihood):
     r"""
     Calculates a log-likelihood assuming independent Gaussian-distributed noise
     at each time point where :math:`\sigma\sim U(a,b)` has been integrated out
-    of the joint posterior of :math:`p(\theta,\sigma|X)`,
+    of the joint posterior of :math:`p(\theta,\sigma|X)` in the same way as in
+    :class:`pints.GaussianIntegratedLogUniformLogLikelihood`,
 
     .. math::
         \begin{align} p(\theta|X) &= \int_{0}^{\infty} p(\theta, \sigma|X)
         \mathrm{d}\sigma\\
         &\propto \int_{0}^{\infty} p(X|\theta, \sigma) p(\theta, \sigma)
-        \mathrm{d}\sigma,\end{align}
+        \mathrm{d}\sigma.\end{align}
 
     Note that this is exactly the same statistical model as
     :class:`pints.GaussianLogLikelihood` with a uniform prior on
     :math:`\sigma`.
 
-    A possible advantage of this log-likelihood compared with using a
+    Similar to :class:`pints.GaussianIntegratedLogUniformLogLikelihood`, a
+    possible advantage of this log-likelihood compared with using a
     :class:`pints.GaussianLogLikelihood`, is that it has one fewer parameters
     (:math:`sigma`) which may speed up convergence to the posterior
     distribution, especially for multi-output problems which will have
@@ -490,7 +573,7 @@ class GaussianIntegratedUniformLogLikelihood(pints.ProblemLogLikelihood):
     The log-likelihood is given in terms of the sum of squared errors:
 
     .. math::
-        SSE = \sum_{i=1}^n (f_i(\theta) - y_i)^2
+        SSE = \sum_{i=1}^n (f_i(\theta) - y_i)^2,
 
     and is given up to a normalisation constant by:
 
@@ -501,11 +584,11 @@ class GaussianIntegratedUniformLogLikelihood(pints.ProblemLogLikelihood):
             & - \text{log}(2 (b - a) \sqrt(2)) \\
             & + (1 / 2 - n / 2) \text{log}(SSE) \\
             & + \text{log}\left[\Gamma((n - 1) / 2, \frac{SSE}{2 b^2}) -
-                \Gamma((n - 1) / 2, \frac{SSE}{2 a^2}) \right]
+                \Gamma((n - 1) / 2, \frac{SSE}{2 a^2}) \right],
         \end{align}
 
     where :math:`\Gamma(u,v)` is the upper incomplete gamma function as defined
-    here: https://en.wikipedia.org/wiki/Incomplete_gamma_function
+    here: https://en.wikipedia.org/wiki/Incomplete_gamma_function.
 
     This log-likelihood is inherently a Bayesian method since it assumes a
     uniform prior on :math:`\sigma\sim U(a,b)`. However using this likelihood
diff --git a/pints/plot/_pairwise.py b/pints/plot/_pairwise.py
index 83ffa8eec..00433de7b 100644
--- a/pints/plot/_pairwise.py
+++ b/pints/plot/_pairwise.py
@@ -43,7 +43,7 @@ def pairwise(samples,
     heatmap
         Set to ``True`` to plot heatmap for the pairwise plots.
         Cannot be used together with ``kde``.
-    Opacity
+    opacity
         This value can be used to manually set the opacity of the
         points in the scatter plots (when ``kde=False`` and ``heatmap=False``
         only).


Test patch:
diff --git a/pints/tests/test_log_likelihoods.py b/pints/tests/test_log_likelihoods.py
index 9fcf2323c..693097f64 100755
--- a/pints/tests/test_log_likelihoods.py
+++ b/pints/tests/test_log_likelihoods.py
@@ -802,6 +802,97 @@ def test_negative_sigma(self):
             self.assertTrue(np.isnan(dl))
 
 
+class TestGaussianIntegratedLogUniformLogLikelihood(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        # Create test single output test model
+        cls.model_single = pints.toy.ConstantModel(1)
+        cls.model_multi = pints.toy.ConstantModel(4)
+
+        # Generate test data
+        cls.times = np.asarray([1, 2, 3])
+        cls.n_times = len(cls.times)
+        cls.data_single = np.asarray([1.0, -10.7, 15.5])
+        cls.data_multi = np.asarray([
+            [3.4, 4.3, 22.0, -7.3],
+            [11.1, 12.2, 13.9, 5.0],
+            [-0.4, -12.3, -8.3, -1.2]])
+
+    def test_call_list(self):
+        # Convert data to list
+        values = self.data_single.tolist()
+
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
+
+        # Create log_likelihood
+        log_likelihood = pints.GaussianIntegratedLogUniformLogLikelihood(
+            problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertAlmostEqual(score, -9.278656018336216)
+
+    def test_call_one_dim_array(self):
+        # Convert data to array of shape (n_times,)
+        values = np.reshape(self.data_single, (self.n_times,))
+
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
+
+        # Create log_likelihood
+        log_likelihood = pints.GaussianIntegratedLogUniformLogLikelihood(
+            problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertAlmostEqual(score, -9.278656018336216)
+
+    def test_call_two_dim_array_single(self):
+        # Convert data to array of shape (n_times, 1)
+        values = np.reshape(self.data_single, (self.n_times, 1))
+
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
+
+        # Create log_likelihood
+        log_likelihood = pints.GaussianIntegratedLogUniformLogLikelihood(
+            problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertAlmostEqual(score, -9.278656018336216)
+
+    def test_call_two_dim_array_multi(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Create log_likelihood
+        log_likelihood = pints.GaussianIntegratedLogUniformLogLikelihood(
+            problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertAlmostEqual(score, -34.36281460402985)
+
+
 class TestGaussianIntegratedUniformLogLikelihood(unittest.TestCase):
 
     @classmethod