Repository: pints-team/pints

Issue:
Add score function 'RRMSE'?
Sometime when I compare different output where the absolute magnitude does not really matter, I prefer to have a 'relative' RMSE score, which is simply `RMSE divided by RMSD of observed data to zero`. Maybe worth adding to PINTS?

Discussion:

Patch:
diff --git a/docs/source/error_measures.rst b/docs/source/error_measures.rst
index dc47ced9..6e998949 100644
--- a/docs/source/error_measures.rst
+++ b/docs/source/error_measures.rst
@@ -17,6 +17,8 @@ Example::
 
 .. autoclass:: MeanSquaredError
 
+.. autoclass:: NormalisedRootMeanSquaredError
+
 .. autoclass:: ProbabilityBasedError
 
 .. autoclass:: ProblemErrorMeasure
diff --git a/pints/__init__.py b/pints/__init__.py
index 47292639..d06f37a5 100644
--- a/pints/__init__.py
+++ b/pints/__init__.py
@@ -135,11 +135,12 @@ from ._boundaries import (
 #
 from ._error_measures import (
     ErrorMeasure,
-    ProblemErrorMeasure,
-    ProbabilityBasedError,
-    SumOfErrors,
     MeanSquaredError,
+    NormalisedRootMeanSquaredError,
+    ProbabilityBasedError,
+    ProblemErrorMeasure,
     RootMeanSquaredError,
+    SumOfErrors,
     SumOfSquaresError,
 )
 
diff --git a/pints/_error_measures.py b/pints/_error_measures.py
index dae77adc..56a4ba3d 100644
--- a/pints/_error_measures.py
+++ b/pints/_error_measures.py
@@ -64,6 +64,81 @@ class ProblemErrorMeasure(ErrorMeasure):
         return self._n_parameters
 
 
+class MeanSquaredError(ProblemErrorMeasure):
+    """
+    Calculates the mean square error: ``f = sum( (x[i] - y[i])**2 ) / n``,
+    where ``n`` is the product of the number of times in the time series and
+    the number of outputs of the problem.
+
+    Extends :class:`ProblemErrorMeasure`.
+
+    Parameters
+    ----------
+    problem
+        A :class:`pints.SingleOutputProblem` or
+        :class:`pints.MultiOutputProblem`.
+    weights
+        An optional sequence of (float) weights, exactly one per problem
+        output. If no weights are specified all sums will be weighted equally.
+    """
+    def __init__(self, problem, weights=None):
+        super(MeanSquaredError, self).__init__(problem)
+        self._ninv = 1.0 / np.product(self._values.shape)
+
+        if weights is None:
+            weights = [1] * self._n_outputs
+        elif self._n_outputs != len(weights):
+            raise ValueError(
+                'Number of weights must match number of problem outputs.')
+        # Check weights
+        self._weights = np.asarray([float(w) for w in weights])
+
+    def __call__(self, x):
+        return np.sum(((np.sum((self._problem.evaluate(x) - self._values)**2,
+                               axis=0) * self._weights) * self._ninv),
+                      axis=0)
+
+    def evaluateS1(self, x):
+        """ See :meth:`ErrorMeasure.evaluateS1()`. """
+        y, dy = self._problem.evaluateS1(x)
+        dy = dy.reshape((self._n_times, self._n_outputs, self._n_parameters))
+        r = y - self._values
+        e = self._ninv * np.sum(np.sum(r**2, axis=0) * self._weights, axis=0)
+        de = 2 * self._ninv * np.sum(np.sum((r.T * dy.T), axis=2) *
+                                     self._weights, axis=1)
+        return e, de
+
+
+class NormalisedRootMeanSquaredError(ProblemErrorMeasure):
+    """
+    Calculates a normalised root mean squared error. This is similar to the
+    (unnormalised) :class:`RootMeanSquaredError`, which is defined as
+    ``f = sqrt( sum( (x[i] - y[i])**2 / n) ) / C``
+    with a normalisation constant ``C`` given by
+    ``C = sqrt( sum((y[i])**2 / n) )``.
+
+    Extends :class:`ProblemErrorMeasure`.
+
+    Parameters
+    ----------
+    problem
+        A :class:`pints.SingleOutputProblem`
+    """
+    def __init__(self, problem):
+        super(NormalisedRootMeanSquaredError, self).__init__(problem)
+
+        if not isinstance(problem, pints.SingleOutputProblem):
+            raise ValueError(
+                'This measure is only defined for single output problems.')
+
+        self._ninv = 1.0 / len(self._values)
+        self._norm = 1.0 / np.sqrt(self._ninv * np.sum(self._values**2))
+
+    def __call__(self, x):
+        return self._norm * np.sqrt(self._ninv * np.sum(
+            (self._problem.evaluate(x) - self._values)**2))
+
+
 class ProbabilityBasedError(ErrorMeasure):
     """
     Changes the sign of a :class:`LogPDF` to use it as an error. Minimising
@@ -101,6 +176,32 @@ class ProbabilityBasedError(ErrorMeasure):
         return self._log_pdf.n_parameters()
 
 
+class RootMeanSquaredError(ProblemErrorMeasure):
+    """
+    Calculates a root mean squared error (RMSE):
+    ``f = sqrt( sum( (x[i] - y[i])**2 / n) )``
+
+    Extends :class:`ProblemErrorMeasure`.
+
+    Parameters
+    ----------
+    problem
+        A :class:`pints.SingleOutputProblem`
+    """
+    def __init__(self, problem):
+        super(RootMeanSquaredError, self).__init__(problem)
+
+        if not isinstance(problem, pints.SingleOutputProblem):
+            raise ValueError(
+                'This measure is only defined for single output problems.')
+
+        self._ninv = 1.0 / len(self._values)
+
+    def __call__(self, x):
+        return np.sqrt(self._ninv * np.sum(
+            (self._problem.evaluate(x) - self._values)**2))
+
+
 class SumOfErrors(ErrorMeasure):
     """
     Calculates a sum of :class:`ErrorMeasure` objects, all defined on the same
@@ -201,77 +302,6 @@ class SumOfErrors(ErrorMeasure):
         return self._n_parameters
 
 
-class MeanSquaredError(ProblemErrorMeasure):
-    """
-    Calculates the mean square error: ``f = sum( (x[i] - y[i])**2 ) / n``,
-    where ``n`` is the product of the number of times in the time series and
-    the number of outputs of the problem.
-
-    Extends :class:`ProblemErrorMeasure`.
-
-    Parameters
-    ----------
-    problem
-        A :class:`pints.SingleOutputProblem` or
-        :class:`pints.MultiOutputProblem`.
-    weights
-        An optional sequence of (float) weights, exactly one per problem
-        output. If no weights are specified all sums will be weighted equally.
-    """
-    def __init__(self, problem, weights=None):
-        super(MeanSquaredError, self).__init__(problem)
-        self._ninv = 1.0 / np.product(self._values.shape)
-
-        if weights is None:
-            weights = [1] * self._n_outputs
-        elif self._n_outputs != len(weights):
-            raise ValueError(
-                'Number of weights must match number of problem outputs.')
-        # Check weights
-        self._weights = np.asarray([float(w) for w in weights])
-
-    def __call__(self, x):
-        return np.sum(((np.sum((self._problem.evaluate(x) - self._values)**2,
-                               axis=0) * self._weights) * self._ninv),
-                      axis=0)
-
-    def evaluateS1(self, x):
-        """ See :meth:`ErrorMeasure.evaluateS1()`. """
-        y, dy = self._problem.evaluateS1(x)
-        dy = dy.reshape((self._n_times, self._n_outputs, self._n_parameters))
-        r = y - self._values
-        e = self._ninv * np.sum(np.sum(r**2, axis=0) * self._weights, axis=0)
-        de = 2 * self._ninv * np.sum(np.sum((r.T * dy.T), axis=2) *
-                                     self._weights, axis=1)
-        return e, de
-
-
-class RootMeanSquaredError(ProblemErrorMeasure):
-    """
-    Calculates a root mean squared error (RMSE):
-    ``f = sqrt( sum( (x[i] - y[i])**2 / n) )``
-
-    Extends :class:`ProblemErrorMeasure`.
-
-    Parameters
-    ----------
-    problem
-        A :class:`pints.SingleOutputProblem`
-    """
-    def __init__(self, problem):
-        super(RootMeanSquaredError, self).__init__(problem)
-
-        if not isinstance(problem, pints.SingleOutputProblem):
-            raise ValueError(
-                'This measure is only defined for single output problems.')
-
-        self._ninv = 1.0 / len(self._values)
-
-    def __call__(self, x):
-        return np.sqrt(self._ninv * np.sum(
-            (self._problem.evaluate(x) - self._values)**2))
-
-
 class SumOfSquaresError(ProblemErrorMeasure):
     """
     Calculates a sum-of-squares error: ``f = sum( (x[i] - y[i])**2 )``
diff --git a/pints/_log_likelihoods.py b/pints/_log_likelihoods.py
index 1d8a2f4a..6a3eaa50 100644
--- a/pints/_log_likelihoods.py
+++ b/pints/_log_likelihoods.py
@@ -137,6 +137,66 @@ class ARMA11LogLikelihood(pints.ProblemLogLikelihood):
                       - np.sum(autocorr_error**2, axis=0) / (2 * sigma**2))
 
 
+class CauchyLogLikelihood(pints.ProblemLogLikelihood):
+    r"""
+    Calculates a log-likelihood assuming independent Cauchy-distributed noise
+    at each time point, and adds one parameter: the scale (``sigma``).
+
+    For a noise characterised by ``sigma``, the log-likelihood is of the form:
+
+    .. math::
+        \log{L(\theta, \sigma)} =
+              -N\log \pi - N\log \sigma
+              -\sum_{i=1}^N\log(1 +
+            \frac{x_i - f(\theta)}{\sigma}^2)
+
+    Extends :class:`ProblemLogLikelihood`.
+
+    Parameters
+    ----------
+    problem
+        A :class:`SingleOutputProblem` or :class:`MultiOutputProblem`. For a
+        single-output problem one parameter is added ``sigma``, where
+        ``sigma`` is scale, for a multi-output problem ``n_outputs``
+        parameters are added.
+    """
+
+    def __init__(self, problem):
+        super(CauchyLogLikelihood, self).__init__(problem)
+
+        # Get number of times, number of outputs
+        self._nt = len(self._times)
+        self._no = problem.n_outputs()
+
+        # Add parameters to problem (one for each output)
+        self._n_parameters = problem.n_parameters() + self._no
+
+        # Pre-calculate
+        self._n = len(self._times)
+        self._n_log_pi = self._n * np.log(np.pi)
+
+    def __call__(self, x):
+        # For multiparameter problems the parameters are stored as
+        # (model_params_1, model_params_2, ..., model_params_k,
+        # sigma_1, sigma_2,...)
+        n = self._n
+        m = self._no
+
+        # problem parameters
+        problem_parameters = x[:-m]
+        error = self._values - self._problem.evaluate(problem_parameters)
+
+        # Distribution parameters
+        sigma = np.asarray(x[-m:])
+
+        # Calculate
+        return np.sum(
+            - self._n_log_pi
+            - n * np.log(sigma)
+            - np.sum(np.log(1 + (error / sigma)**2), axis=0)
+        )
+
+
 class GaussianIntegratedUniformLogLikelihood(pints.ProblemLogLikelihood):
     r"""
     Calculates a log-likelihood assuming independent Gaussian-distributed noise
@@ -278,66 +338,6 @@ class GaussianIntegratedUniformLogLikelihood(pints.ProblemLogLikelihood):
         )
 
 
-class CauchyLogLikelihood(pints.ProblemLogLikelihood):
-    r"""
-    Calculates a log-likelihood assuming independent Cauchy-distributed noise
-    at each time point, and adds one parameter: the scale (``sigma``).
-
-    For a noise characterised by ``sigma``, the log-likelihood is of the form:
-
-    .. math::
-        \log{L(\theta, \sigma)} =
-              -N\log \pi - N\log \sigma
-              -\sum_{i=1}^N\log(1 +
-            \frac{x_i - f(\theta)}{\sigma}^2)
-
-    Extends :class:`ProblemLogLikelihood`.
-
-    Parameters
-    ----------
-    problem
-        A :class:`SingleOutputProblem` or :class:`MultiOutputProblem`. For a
-        single-output problem one parameter is added ``sigma``, where
-        ``sigma`` is scale, for a multi-output problem ``n_outputs``
-        parameters are added.
-    """
-
-    def __init__(self, problem):
-        super(CauchyLogLikelihood, self).__init__(problem)
-
-        # Get number of times, number of outputs
-        self._nt = len(self._times)
-        self._no = problem.n_outputs()
-
-        # Add parameters to problem (one for each output)
-        self._n_parameters = problem.n_parameters() + self._no
-
-        # Pre-calculate
-        self._n = len(self._times)
-        self._n_log_pi = self._n * np.log(np.pi)
-
-    def __call__(self, x):
-        # For multiparameter problems the parameters are stored as
-        # (model_params_1, model_params_2, ..., model_params_k,
-        # sigma_1, sigma_2,...)
-        n = self._n
-        m = self._no
-
-        # problem parameters
-        problem_parameters = x[:-m]
-        error = self._values - self._problem.evaluate(problem_parameters)
-
-        # Distribution parameters
-        sigma = np.asarray(x[-m:])
-
-        # Calculate
-        return np.sum(
-            - self._n_log_pi
-            - n * np.log(sigma)
-            - np.sum(np.log(1 + (error / sigma)**2), axis=0)
-        )
-
-
 class GaussianKnownSigmaLogLikelihood(pints.ProblemLogLikelihood):
     r"""
     Calculates a log-likelihood assuming independent Gaussian noise at each


Test patch:
diff --git a/pints/tests/test_error_measures.py b/pints/tests/test_error_measures.py
index 25511b74..a1e6861a 100755
--- a/pints/tests/test_error_measures.py
+++ b/pints/tests/test_error_measures.py
@@ -256,6 +256,28 @@ class TestErrorMeasures(unittest.TestCase):
         self.assertEqual(dex[0], -1)
         self.assertEqual(dex[1], -6)
 
+    def test_normalised_root_mean_squared_error(self):
+        # Tests :class:`pints.NormalisedRootMeanSquaredError`.
+
+        p = MiniProblem()
+        e = pints.NormalisedRootMeanSquaredError(p)
+        self.assertEqual(e.n_parameters(), 3)
+        float(e([1, 2, 3]))
+        self.assertEqual(e([-1, 2, 3]), 0)
+        self.assertNotEqual(np.all(e([1, 2, 3])), 0)
+        x = [0, 0, 0]
+        y = 1.0
+        self.assertEqual(e(x), y)
+        x = [1, 1, 1]
+        y = np.sqrt((4 + 1 + 4) / 3) / np.sqrt((1 + 4 + 9) / 3)
+        self.assertAlmostEqual(e(x), y)
+
+        p = MultiMiniProblem()
+        self.assertRaisesRegex(
+            ValueError,
+            'This measure is only defined for single output problems.',
+            pints.NormalisedRootMeanSquaredError, p)
+
     def test_probability_based_error(self):
         # Tests :class:`pints.ProbabilityBasedError`.
 
@@ -290,7 +312,10 @@ class TestErrorMeasures(unittest.TestCase):
         self.assertEqual(e(x), y)
 
         p = MultiMiniProblem()
-        self.assertRaises(ValueError, pints.RootMeanSquaredError, p)
+        self.assertRaisesRegex(
+            ValueError,
+            'This measure is only defined for single output problems.',
+            pints.RootMeanSquaredError, p)
 
     def test_sum_of_squares_error_single(self):
         # Tests :class:`pints.MeanSquaredError` with a single output.