diff --git a/pints/_util.py b/pints/_util.py
index 4230243f..7114f12f 100644
--- a/pints/_util.py
+++ b/pints/_util.py
@@ -91,7 +91,8 @@ def vector(x):
     if x.ndim != 1:
         n = np.max(x.shape)
         if np.prod(x.shape) != n:
-            raise ValueError('Unable to convert to 1d vector of scalar values')
+            raise ValueError(
+                'Unable to convert to 1d vector of scalar values.')
         x = x.reshape((n,))
     return x
 
diff --git a/pints/toy/_constant_model.py b/pints/toy/_constant_model.py
index 03df1d3f..22f114c6 100644
--- a/pints/toy/_constant_model.py
+++ b/pints/toy/_constant_model.py
@@ -107,6 +107,7 @@ class ConstantModel(pints.ForwardModelS1):
             #   [df2/dp1, df2/dp2]]
             # i.e.
             #  [[1, 0],
-            #   [0, 1]]
-            dy = np.tile(np.diag(np.ones(self._n)), (len(times), 1, 1))
+            #   [0, 2]]
+            dy = np.tile(
+                np.diag(np.arange(1, self._n + 1)), (len(times), 1, 1))
         return (y, dy)

Add tests for different input shapes for every likelihood / error measure
They should all be tested (and work for)
- Single-value data of shape ``(n_times,)`` (or just a list)
- A NumPy array of shape ``(n_times,1)``
- A NumPy array of shape ``(n_times, n_outputs)`` where n_outputs > 1

See #604, #605 
diff --git a/pints/tests/test_error_measures.py b/pints/tests/test_error_measures.py
index a1e6861a..cfe9108c 100755
--- a/pints/tests/test_error_measures.py
+++ b/pints/tests/test_error_measures.py
@@ -13,6 +13,13 @@ import unittest
 import numpy as np
 
 
+# Unit testing in Python 2 and 3
+try:
+    unittest.TestCase.assertRaisesRegex
+except AttributeError:
+    unittest.TestCase.assertRaisesRegex = unittest.TestCase.assertRaisesRegexp
+
+
 class MiniProblem(pints.SingleOutputProblem):
     def __init__(self):
         self._t = pints.vector([1, 2, 3])
@@ -95,381 +102,512 @@ class MiniLogPDF(pints.LogPDF):
         return 10, np.array([1, 2, 3])
 
 
-class TestErrorMeasures(unittest.TestCase):
-    """
-    Tests the ErrorMeasure classes
-    """
-    def __init__(self, name):
-        super(TestErrorMeasures, self).__init__(name)
+class TestMeanSquaredError(unittest.TestCase):
 
-    def test_mean_squared_error_single(self):
-        # Tests :class:`pints.MeanSquaredError` with a single output.
+    @classmethod
+    def setUpClass(cls):
+        # Create test single output test model
+        cls.model_single = pints.toy.ConstantModel(1)
+        cls.model_multi = pints.toy.ConstantModel(2)
 
-        # Set up problem
-        model = pints.toy.ConstantModel(1)
-        times = [1, 2, 3]
-        values = [1, 1, 1]
-        p = pints.SingleOutputProblem(model, times, values)
-
-        # Test
-        e = pints.MeanSquaredError(p)
-        self.assertEqual(e.n_parameters(), 1)
-        float(e([1]))
-        self.assertEqual(e([1]), 0)
-        self.assertEqual(e([2]), 1)
-        self.assertEqual(e([0]), 1)
-        self.assertEqual(e([3]), 4)
-
-        # Derivatives
-        for x in [1, 2, 3, 4]:
-            y, dy = e.evaluateS1([x])
-            r = x - 1
-            self.assertEqual(y, e([x]))
-            self.assertEqual(dy.shape, (1, ))
-            self.assertTrue(np.all(dy == 2 * r))
-
-    def test_mean_squared_error_multi(self):
-        # Tests :class:`pints.MeanSquaredError` with multiple outputs.
-
-        # Set up problem
-        model = pints.toy.ConstantModel(2)
-        times = [1, 2, 3]
-        values = [[1, 4], [1, 4], [1, 4]]
-        p = pints.MultiOutputProblem(model, times, values)
-
-        # Test
-        e = pints.MeanSquaredError(p)
-        self.assertEqual(e.n_parameters(), 2)
-        float(e([1, 2]))
-        self.assertEqual(e([1, 2]), 0)      # 0
-        self.assertEqual(e([2, 2]), 0.5)    # (3*(1^2+0^2)) / 6 = (1^2+0^2) / 2
-        self.assertEqual(e([2, 3]), 2.5)    # (1^2+2^2) / 2 = 2.5
-        self.assertEqual(e([3, 4]), 10)     # (2^2+4^2) / 2 = 10
-
-        # Derivatives
-        values = np.array([[1, 4], [2, 7], [3, 10]])
-        p = pints.MultiOutputProblem(model, times, values)
-        e = pints.MeanSquaredError(p)
-        x = [1, 2]
-
-        # Model outputs are 3 times [1, 4]
-        # 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, 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)
-        self.assertTrue(np.all(rx == np.array(
-            [[-0, -0],
-             [-1, -3],
-             [-2, -6]]
-        )))
-        self.assertAlmostEqual(e(x), np.sum(rx**2) / 6)
-
-        # Now with derivatives
-        ex, dex = e.evaluateS1(x)
-
-        # Check error
-        self.assertAlmostEqual(ex, e(x))
-
-        # Check derivatives. Shape is (parameters, )
-        self.assertEqual(dex.shape, (2, ))
-
-        # Residuals are: [[0, 0], [-1, -3], [-2, -6]]
-        # 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) * 1 = (1 / 3) * -9 * 1 = -3
-        self.assertEqual(dex[0], -1)
-        self.assertEqual(dex[1], -3)
-
-    def test_mean_squared_error_weighted(self):
-        # Tests :class:`pints.MeanSquaredError` with weighted outputs.
-
-        # Set up problem
-        model = pints.toy.ConstantModel(2)
-        times = [1, 2, 3]
-        values = [[1, 4], [1, 4], [1, 4]]
-        p = pints.MultiOutputProblem(model, times, values)
+        # Generate test data
+        cls.times = [1, 2, 3]
+        cls.n_times = len(cls.times)
+        cls.data_single = np.array([1, 1, 1])
+        cls.data_multi = np.array([
+            [1, 4],
+            [1, 4],
+            [1, 4]])
+
+    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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = mean(input - 1) ** 2
+        self.assertEqual(score, 4)
+
+    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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = mean(input - 1) ** 2
+        self.assertEqual(score, 4)
+
+    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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = mean(input - 1) ** 2
+        self.assertEqual(score, 4)
+
+    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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # exp = (mean(input[0] - 1) ** 2 + mean(2 * input[1] - 4) ** 2) / 2
+        self.assertEqual(score, 10)
+
+    def test_call_two_dim_array_multi_weighted(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Create error measure with weighted input
+        weights = [1, 2]
+        error = pints.MeanSquaredError(problem, weights=weights)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # exp = (weight[0] * mean(input[0] - 1) ** 2 +
+        # weight[1] * mean(2 * input[1] - 4) ** 2) / 2
+        self.assertEqual(score, 18)
+
+    def test_evaluateS1_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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
+
+        # Check that partials are correct
+        # Expected = 2 * mean(input - 1)
+        self.assertEqual(deriv, 2 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
+
+        # Check that partials are correct
+        # Expected = 2 * mean(input - 1)
+        self.assertEqual(deriv, 2 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
 
-        # Test
-        e = pints.MeanSquaredError(p, weights=[1, 2])
+        # Check that partials are correct
+        # Expected = 2 * mean(input - 1)
+        self.assertEqual(deriv, 2 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.MeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (2,))
+
+        # Check that partials are correct
+        # Expectation = [mean(input[0] - 1), 2 * mean(2 * input[1] - 4)]
+        self.assertEqual(deriv[0], test_parameters[0] - 1)
+        self.assertEqual(deriv[1], 2 * (2 * test_parameters[1] - 4))
+
+    def test_evaluateS1_two_dim_array_multi_weighted(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Create error measure with weighted inputs
+        weights = [1, 2]
+        error = pints.MeanSquaredError(problem, weights=weights)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (2,))
+
+        # Check that partials are correct
+        # expectation = [weight [0] * mean(input[0] - 1),
+        # weight[1] * 2 * mean(2 * input[1] - 4)]
+        self.assertEqual(deriv[0], weights[0] * (test_parameters[0] - 1))
+        self.assertEqual(
+            deriv[1], weights[1] * 2 * (2 * test_parameters[1] - 4))
+
+    def test_bad_constructor(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_single, self.times, self.data_single)
+
+        # Test invalid weight shape
+        weights = [1, 2, 3]
         self.assertRaisesRegex(
-            ValueError, 'Number of weights',
-            pints.MeanSquaredError, p, weights=[1, 2, 3])
-        self.assertEqual(e.n_parameters(), 2)
-        float(e([1, 2]))
-        self.assertEqual(e([1, 2]), 0)      # 0
-        self.assertEqual(e([2, 2]), 0.5)    # (3*(1^2*1+0^2*2)) / 6 = 1 / 2
-        self.assertEqual(e([2, 3]), 4.5)    # (1^2*1+2^2*2) / 2 = 4.5
-        self.assertEqual(e([3, 4]), 18)     # (2^2*1+4^2*2) / 2 = 18
-
-        # Derivatives
-        values = np.array([[1, 4], [2, 7], [3, 10]])
-        p = pints.MultiOutputProblem(model, times, values)
-        w = np.array([1, 2])
-        e = pints.MeanSquaredError(p, weights=w)
-        x = [1, 2]
-
-        # Model outputs are 3 times [1, 4]
-        # 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, 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)
-        self.assertTrue(np.all(rx == np.array(
-            [[-0, -0],
-             [-1, -3],
-             [-2, -6]]
-        )))
-        self.assertAlmostEqual(e(x), np.sum(np.sum(rx**2, axis=0) * w) / 6)
-
-        # Now with derivatives
-        ex, dex = e.evaluateS1(x)
-
-        # Check error
-        self.assertAlmostEqual(ex, e(x))
-
-        # Check derivatives. Shape is (parameters, )
-        self.assertEqual(dex.shape, (2, ))
-
-        # Residuals are: [[0, 0], [-1, -3], [-2, -6]]
-        # Derivatives are: [[1, 0], [0, 2]]
-        # dex1 is: (2 / nt / no) * (0 - 1 - 2) * 1 * 1
-        #        = (1 / 3) * -3 * 1 * 1
-        #        = -1
-        # 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], -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)
+            ValueError,
+            'Number of weights must match number of problem outputs.',
+            pints.MeanSquaredError, problem, weights)
+
+
+class TestNormalisedRootMeanSquaredError(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(2)
+
+        # Generate test data
+        cls.times = [1, 2, 3]
+        cls.n_times = len(cls.times)
+        cls.data_single = np.array([2, 2, 2])
+        cls.data_multi = np.array([
+            [1, 4],
+            [1, 4],
+            [1, 4]])
+
+    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 error measure
+        error = pints.NormalisedRootMeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2)) / sqrt(mean(2 ** 2))
+        self.assertEqual(score, 0.5)
 
-        p = MultiMiniProblem()
+    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 error measure
+        error = pints.NormalisedRootMeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2)) / sqrt(mean(2 ** 2))
+        self.assertEqual(score, 0.5)
+
+    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 error measure
+        error = pints.NormalisedRootMeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2)) / sqrt(mean(2 ** 2))
+        self.assertEqual(score, 0.5)
+
+    def test_not_implemented_error(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 error measure
+        error = pints.NormalisedRootMeanSquaredError(problem)
+
+        # Check that not implemented error is raised for evaluateS1
+        self.assertRaisesRegex(
+            NotImplementedError,
+            '',
+            error.evaluateS1, 1)
+
+    def test_bad_constructor(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Check that an error is raised for multi-output problems
         self.assertRaisesRegex(
             ValueError,
             'This measure is only defined for single output problems.',
-            pints.NormalisedRootMeanSquaredError, p)
+            pints.NormalisedRootMeanSquaredError, problem)
 
-    def test_probability_based_error(self):
-        # Tests :class:`pints.ProbabilityBasedError`.
 
-        p = MiniLogPDF()
-        e = pints.ProbabilityBasedError(p)
-        self.assertEqual(e.n_parameters(), 3)
-        self.assertEqual(e([1, 2, 3]), -10)
-        p = MiniProblem()
-        self.assertRaises(ValueError, pints.ProbabilityBasedError, p)
+class TestProbabilityBasedError(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        # Create probability based problem
+        cls.problem = MiniLogPDF()
 
-        # Test derivatives
-        x = [1, 2, 3]
-        y, dy = e.evaluateS1(x)
-        self.assertEqual(y, e(x))
-        self.assertEqual(dy.shape, (3, ))
-        self.assertTrue(np.all(dy == [-1, -2, -3]))
+    def test_call(self):
+        # Create error measure
+        error = pints.ProbabilityBasedError(self.problem)
 
-    def test_root_mean_squared_error(self):
-        # Tests :class:`pints.RootMeanSquaredError`.
+        # Evaluate likelihood for test parameters
+        test_parameters = [1, 2, 3]
+        score = error(test_parameters)
 
-        p = MiniProblem()
-        e = pints.RootMeanSquaredError(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 = np.sqrt((1 + 4 + 9) / 3)
-        self.assertAlmostEqual(e(x), y)
-        x = [1, 1, 1]
-        y = np.sqrt((4 + 1 + 4) / 3)
-        self.assertEqual(e(x), y)
+        # Check that error returns expected value
+        self.assertEqual(score, -10)
+
+    def test_evaluateS1(self):
+        # Create error measure
+        error = pints.ProbabilityBasedError(self.problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [1, 2, 3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that error returns expected value
+        self.assertEqual(score, error(test_parameters))
+
+        # Check dimension of partial derivatives
+        self.assertEqual(deriv.shape, (3,))
 
-        p = MultiMiniProblem()
+        # Check that partials are computed correctly
+        self.assertEqual(deriv[0], -1)
+        self.assertEqual(deriv[1], -2)
+        self.assertEqual(deriv[2], -3)
+
+    def test_n_parameters(self):
+        # Create error measure
+        error = pints.ProbabilityBasedError(self.problem)
+
+        # Get number of parameters
+        n_parameters = error.n_parameters()
+
+        # Check number of parameters
+        self.assertEqual(n_parameters, 3)
+
+    def test_bad_constructor(self):
+        # Check that an error is raised for multi-output problems
         self.assertRaisesRegex(
             ValueError,
-            'This measure is only defined for single output problems.',
-            pints.RootMeanSquaredError, p)
+            'Given log_pdf must be an instance of pints.LogPDF.',
+            pints.ProbabilityBasedError, MiniProblem())
 
-    def test_sum_of_squares_error_single(self):
-        # Tests :class:`pints.MeanSquaredError` with a single output.
 
-        # Set up problem
-        model = pints.toy.ConstantModel(1)
-        times = [1, 2, 3]
-        values = [1, 1, 1]
-        p = pints.SingleOutputProblem(model, times, values)
-
-        # Test
-        e = pints.SumOfSquaresError(p)
-        self.assertEqual(e.n_parameters(), 1)
-        float(e([1]))
-        self.assertEqual(e([1]), 0)
-        self.assertEqual(e([2]), 3)
-        self.assertEqual(e([0]), 3)
-        self.assertEqual(e([3]), 12)
-
-        # Derivatives
-        for x in [1, 2, 3, 4]:
-            ex, dex = e.evaluateS1([x])
-            r = x - 1
-            self.assertEqual(ex, e([x]))
-            self.assertEqual(dex.shape, (1, ))
-            self.assertEqual(dex[0], 2 * 3 * r)
-
-    def test_sum_of_squares_error_multi(self):
-        # Tests :class:`pints.MeanSquaredError` with multiple outputs.
-
-        # Set up problem
-        model = pints.toy.ConstantModel(2)
-        times = [1, 2, 3]
-        values = [[1, 4], [1, 4], [1, 4]]
-        p = pints.MultiOutputProblem(model, times, values)
-
-        # Test
-        e = pints.SumOfSquaresError(p)
-        self.assertEqual(e.n_parameters(), 2)
-        float(e([1, 2]))
-        self.assertEqual(e([1, 2]), 0)      # 0
-        self.assertEqual(e([2, 2]), 3)      # 3*(1^2+0^2) = 3
-        self.assertEqual(e([2, 3]), 15)     # 3*(1^2+2^2) = 15
-        self.assertEqual(e([3, 4]), 60)     # 3*(2^2+4^2) = 60
-
-        # Derivatives
-        values = np.array([[1, 4], [2, 7], [3, 10]])
-        p = pints.MultiOutputProblem(model, times, values)
-        e = pints.SumOfSquaresError(p)
-        x = [1, 2]
-
-        # Model outputs are 3 times [1,4]
-        # 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, 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)
-        self.assertTrue(np.all(rx == np.array(
-            [[-0, -0],
-             [-1, -3],
-             [-2, -6]]
-        )))
-        self.assertAlmostEqual(e(x), np.sum(rx**2))
-
-        # Now with derivatives
-        ex, dex = e.evaluateS1(x)
-
-        # Check error
-        self.assertTrue(np.all(ex == e(x)))
-
-        # Check derivatives. Shape is (parameters, )
-        self.assertEqual(dex.shape, (2, ))
-
-        # Residuals are: [[0, 0], [-1, -3], [-2, -6]]
-        # 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 * 1 = -18
-        self.assertEqual(dex[0], -6)
-        self.assertEqual(dex[1], -18)
-
-    def test_sum_of_squares_error_weighted(self):
-        # Tests :class:`pints.MeanSquaredError` with weighted outputs.
-
-        # Set up problem
-        model = pints.toy.ConstantModel(2)
-        times = [1, 2, 3]
-        values = [[1, 4], [1, 4], [1, 4]]
-        p = pints.MultiOutputProblem(model, times, values)
+class TestRootMeanSquaredError(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(2)
+
+        # Generate test data
+        cls.times = [1, 2, 3]
+        cls.n_times = len(cls.times)
+        cls.data_single = np.array([1, 1, 1])
+        cls.data_multi = np.array([
+            [1, 4],
+            [1, 4],
+            [1, 4]])
+
+    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 error measure
+        error = pints.RootMeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2))
+        self.assertEqual(score, 2)
+
+    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 error measure
+        error = pints.RootMeanSquaredError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2))
+        self.assertEqual(score, 2)
+
+    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 error measure
+        error = pints.RootMeanSquaredError(problem)
 
-        # Test
-        e = pints.SumOfSquaresError(p, weights=[1, 2])
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sqrt(mean((input - 1) ** 2))
+        self.assertEqual(score, 2)
+
+    def test_not_implemented_error(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 error measure
+        error = pints.RootMeanSquaredError(problem)
+
+        # Check that not implemented error is raised for evaluateS1
         self.assertRaisesRegex(
-            ValueError, 'Number of weights',
-            pints.SumOfSquaresError, p, weights=[1, 2, 3])
-        self.assertEqual(e.n_parameters(), 2)
-        float(e([1, 2]))
-        self.assertEqual(e([1, 2]), 0)      # 0
-        self.assertEqual(e([2, 2]), 3)      # 3*(1^2*1+0^2*2) = 3
-        self.assertEqual(e([2, 3]), 27)     # 3*(1^2*1+2^2*2) = 27
-        self.assertEqual(e([3, 4]), 108)    # 3*(2^2*1+4^2*2) = 108
-
-        # Derivatives
-        values = np.array([[1, 4], [2, 7], [3, 10]])
-        p = pints.MultiOutputProblem(model, times, values)
-        w = np.array([1, 2])
-        e = pints.SumOfSquaresError(p, weights=w)
-        x = [1, 2]
-
-        # Model outputs are 3 times [1, 4]
-        # 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, 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)
-        self.assertTrue(np.all(rx == np.array(
-            [[-0, -0],
-             [-1, -3],
-             [-2, -6]]
-        )))
-        self.assertAlmostEqual(e(x), np.sum(np.sum(rx**2, axis=0) * w))
-
-        # Now with derivatives
-        ex, dex = e.evaluateS1(x)
-
-        # Check error
-        self.assertAlmostEqual(ex, e(x))
-
-        # Check derivatives. Shape is (parameters, )
-        self.assertEqual(dex.shape, (2, ))
-
-        # Residuals are: [[0, 0], [-1, -3], [-2, -6]]
-        # 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 * 1
-        #        = 2 * -9 * 2 * 1
-        #        = -36
-        self.assertEqual(dex[0], -6)
-        self.assertEqual(dex[1], -36)
+            NotImplementedError,
+            '',
+            error.evaluateS1, 1)
+
+    def test_bad_constructor(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Check that an error is raised for multi-output problems
+        self.assertRaisesRegex(
+            ValueError,
+            'This measure is only defined for single output problems.',
+            pints.RootMeanSquaredError, problem)
+
+
+class TestSumOfErrors(unittest.TestCase):
 
     def test_sum_of_errors(self):
         # Tests :class:`pints.SumOfErrors`.
@@ -576,5 +714,251 @@ class TestErrorMeasures(unittest.TestCase):
         self.assertTrue(np.all(dy == dy1 + 2 * dy2))
 
 
+class TestSumOfSquaresError(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(2)
+
+        # Generate test data
+        cls.times = [1, 2, 3]
+        cls.n_times = len(cls.times)
+        cls.data_single = np.array([1, 1, 1])
+        cls.data_multi = np.array([
+            [1, 4],
+            [1, 4],
+            [1, 4]])
+
+    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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sum((input - 1) ** 2)
+        self.assertEqual(score, 12)
+
+    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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sum((input - 1) ** 2)
+        self.assertEqual(score, 12)
+
+    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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Expected = sum((input - 1) ** 2)
+        self.assertEqual(score, 12)
+
+    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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Exp = sum((input[0] - 1) ** 2) + sum((2 * input[1] - 4) ** 2)
+        self.assertEqual(score, 60)
+
+    def test_call_two_dim_array_multi_weighted(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Create error measure with weighted input
+        weights = [1, 2]
+        error = pints.SumOfSquaresError(problem, weights=weights)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score = error(test_parameters)
+
+        # Check that error returns expected value
+        # Exp = weight[0] * sum((input[0] - 1) ** 2) +
+        # weight[1] * sum((2 * input[1] - 4) ** 2)
+        self.assertEqual(score, 108)
+
+    def test_evaluateS1_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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
+
+        # Check that partials are correct
+        # Expected = 2 * sum(input - 1)
+        self.assertEqual(deriv, 2 * 3 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
+
+        # Check that partials are correct
+        # Expected = 2 * sum(input - 1)
+        self.assertEqual(deriv, 2 * 3 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (1,))
+
+        # Check that partials are correct
+        # Expected = 2 * sum(input - 1)
+        self.assertEqual(deriv, 2 * 3 * (test_parameters[0] - 1))
+
+    def test_evaluateS1_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 error measure
+        error = pints.SumOfSquaresError(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (2,))
+
+        # Check that partials are correct
+        # Expectation = [2 * sum(input[0] - 1), 4 * sum(2 * input[1] - 4)]
+        self.assertEqual(deriv[0], 2 * 3 * (test_parameters[0] - 1))
+        self.assertEqual(deriv[1], 4 * 3 * (2 * test_parameters[1] - 4))
+
+    def test_evaluateS1_two_dim_array_multi_weighted(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Create error measure with weighted inputs
+        weights = [1, 2]
+        error = pints.SumOfSquaresError(problem, weights=weights)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [3, 4]
+        score, deriv = error.evaluateS1(test_parameters)
+
+        # Check that returned error is correct
+        self.assertEqual(score, error(test_parameters))
+
+        # Check that partial derivatives are returned for each parameter
+        self.assertEqual(deriv.shape, (2,))
+
+        # Check that partials are correct
+        # Expectation = [weight [0] * 2 * sum(input[0] - 1),
+        # weight[1] * 4 * sum(2 * input[1] - 4)]
+        self.assertEqual(
+            deriv[0], weights[0] * 2 * 3 * (test_parameters[0] - 1))
+        self.assertEqual(
+            deriv[1], weights[1] * 4 * 3 * (2 * test_parameters[1] - 4))
+
+    def test_bad_constructor(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_single, self.times, self.data_single)
+
+        # Test invalid weight shape
+        weights = [1, 2, 3]
+        self.assertRaisesRegex(
+            ValueError,
+            'Number of weights must match number of problem outputs.',
+            pints.SumOfSquaresError, problem, weights)
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/pints/tests/test_log_likelihoods.py b/pints/tests/test_log_likelihoods.py
index 81c48e4a..fdeae8b7 100755
--- a/pints/tests/test_log_likelihoods.py
+++ b/pints/tests/test_log_likelihoods.py
@@ -12,499 +12,875 @@ import pints.toy
 import numpy as np
 
 
-class TestLogLikelihood(unittest.TestCase):
+class TestAR1LogLikelihood(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.5, 7.6, 8.5, 3.4],
+            [1.1, -10.3, 15.6, 5.5],
+            [-10, -30.5, -5, 7.6]])
+
+    def test_call_list(self):
+        # Convert data to list
+        values = self.data_single.tolist()
 
-    def test_scaled_log_likelihood(self):
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
 
-        model = pints.toy.LogisticModel()
-        real_parameters = [0.015, 500]
-        test_parameters = [0.014, 501]
-        sigma = 0.001
-        times = np.linspace(0, 1000, 100)
-        values = model.simulate(real_parameters, times)
+        # Create log_likelihood
+        log_likelihood = pints.AR1LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.5, 5]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -19.706737485492436)
+
+    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(model, times, values)
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
+
+        # Create log_likelihood
+        log_likelihood = pints.AR1LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.5, 5]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -19.706737485492436)
+
+    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.AR1LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.5, 5]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -19.706737485492436)
+
+    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.AR1LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [
+            0, 0, 0, 0, 0.5, 1.0, -0.25, 3.0, 0.9, 10.0, 0.0, 2.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -179.22342804581092)
+
+
+class TestARMA11LogLikelihood(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, 4])
+        cls.n_times = len(cls.times)
+        cls.data_single = np.asarray([3, -4.5, 10.5, 0.3])
+        cls.data_multi = np.asarray([
+            [3.5, 7.6, 8.5, 3.4],
+            [1.1, -10.3, 15.6, 5.5],
+            [-10, -30.5, -5, 7.6],
+            [-12, -10.1, -4, 2.3]])
+
+    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.ARMA11LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.9, -0.4, 1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -171.53031588534171)
+
+    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.ARMA11LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.9, -0.4, 1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -171.53031588534171)
+
+    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.ARMA11LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0.9, -0.4, 1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -171.53031588534171)
+
+    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.ARMA11LogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [
+            0, 0, 0, 0, 0.5, 0.34, 1.0, -0.25, 0.1, 3.0, 0.9, 0.0, 10.0, 0.0,
+            0.9, 2.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -214.17034137601107)
+
+
+class TestCauchyLogLikelihood(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.5, 7.6, 8.5, 3.4],
+            [1.1, -10.3, 15.6, 5.5],
+            [-10, -30.5, -5, 7.6]])
+
+    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.CauchyLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -12.339498654173603)
+
+    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.CauchyLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -12.339498654173603)
+
+    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.CauchyLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -12.339498654173603)
+
+    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.CauchyLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [
+            0, 0, 0, 0, 13, 8, 13.5, 10.5]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -49.51182454195375)
+
+
+class TestGaussianIntegratedUniformLogLikelihood(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.GaussianIntegratedUniformLogLikelihood(
+            problem, 2, 4)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -20.441037907121299)
+
+    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.GaussianIntegratedUniformLogLikelihood(
+            problem, 2, 4)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -20.441037907121299)
+
+    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.GaussianIntegratedUniformLogLikelihood(
+            problem, 2, 4)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -20.441037907121299)
+
+    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.GaussianIntegratedUniformLogLikelihood(
+            problem, 2, 4)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -75.443307614807225)
+
+    def test_call_two_dim_array_multi_non_equal_priors(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.GaussianIntegratedUniformLogLikelihood(
+            problem, [1, 0, 5, 2], [2, 4, 7, 8])
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -71.62076263891457)
+
+    def test_bad_constructor_single(self):
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, self.data_single)
+
+        # Check negative bound
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, -1, 2)
+
+        # Check vanishing interval width
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, 0, 0)
+
+        # Check higher lower bound than upper bound
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, 2, 1)
+
+        # Check wrong prior dimensions
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, [1, 2], 2)
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, 1, [2, 3])
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, [1, 2], [2, 3])
+
+    def test_bad_constructor_multi(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_multi, self.times, self.data_multi)
+
+        # Check wrong prior dimensions
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, 2, 2)
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, [1, 2, 3], [2, 4])
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, [1, 2], [2, 4, 5])
+        self.assertRaises(ValueError,
+                          pints.GaussianIntegratedUniformLogLikelihood,
+                          problem, [1, 3], [2, 2])
+
+
+class TestGaussianKnownSigmaLikelihood(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(3)
+
+        # Generate test data
+        cls.times = [1, 2, 3, 4]
+        cls.n_times = len(cls.times)
+        cls.data_single = np.arange(1, 5) / 5.0
+        cls.data_multi = np.array([
+            [10.7, 3.5, 3.8],
+            [1.1, 3.2, -1.4],
+            [9.3, 0.0, 4.5],
+            [1.2, -3, -10]])
+
+    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.GaussianKnownSigmaLogLikelihood(problem, 1.5)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [-1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -7.3420590096957925)
+
+    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.GaussianKnownSigmaLogLikelihood(problem, 1.5)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [-1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -7.3420590096957925)
+
+    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.GaussianKnownSigmaLogLikelihood(problem, 1.5)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [-1]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -7.3420590096957925)
+
+    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.GaussianKnownSigmaLogLikelihood(problem, 1)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -196.9122623984561)
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
-        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
-            problem, sigma)
-        log_likelihood_scaled = pints.ScaledLogLikelihood(
-            log_likelihood_not_scaled)
+        log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, 1.5)
 
-        eval_not_scaled = log_likelihood_not_scaled(test_parameters)
-        eval_scaled = log_likelihood_scaled(test_parameters)
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [3]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
 
-        self.assertEqual(int(eval_not_scaled), -20959169232)
-        self.assertAlmostEqual(eval_scaled * len(times), eval_not_scaled)
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
 
-        # Test bad constructor
-        self.assertRaises(ValueError, pints.ScaledLogLikelihood, model)
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (1, ))
 
-        # Test single-output derivatives
-        y1, dy1 = log_likelihood_not_scaled.evaluateS1(test_parameters)
-        y2, dy2 = log_likelihood_scaled.evaluateS1(test_parameters)
-        self.assertEqual(y1, log_likelihood_not_scaled(test_parameters))
-        self.assertEqual(dy1.shape, (2, ))
-        self.assertEqual(y2, log_likelihood_scaled(test_parameters))
-        self.assertEqual(dy2.shape, (2, ))
-        dy3 = dy2 * len(times)
-        self.assertAlmostEqual(dy1[0] / dy3[0], 1)
-        self.assertAlmostEqual(dy1[1] / dy3[1], 1)
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -4.444444444444445)
 
-        # Test on multi-output problem
-        model = pints.toy.FitzhughNagumoModel()
-        nt = 10
-        no = model.n_outputs()
-        times = np.linspace(0, 100, nt)
-        values = model.simulate([0.5, 0.5, 0.5], times)
-        problem = pints.MultiOutputProblem(model, times, values)
-        unscaled = pints.GaussianKnownSigmaLogLikelihood(problem, 1)
-        scaled = pints.ScaledLogLikelihood(unscaled)
-        p = [0.1, 0.1, 0.1]
-        x = unscaled(p)
-        y = scaled(p)
-        self.assertAlmostEqual(y, x / nt / no)
+    def test_evaluateS1_one_dim_array(self):
+        # Convert data to array of shape (n_times,)
+        values = np.reshape(self.data_single, (self.n_times,))
 
-        # Test multi-output derivatives
-        y1, dy1 = unscaled.evaluateS1(p)
-        y2, dy2 = scaled.evaluateS1(p)
-        self.assertAlmostEqual(y1, unscaled(p), places=6)
-        self.assertEqual(dy1.shape, (3, ))
-        self.assertAlmostEqual(y2, scaled(p))
-        self.assertEqual(dy2.shape, (3, ))
-        dy3 = dy2 * nt * no
-        self.assertAlmostEqual(dy1[0] / dy3[0], 1)
-        self.assertAlmostEqual(dy1[1] / dy3[1], 1)
-
-        # test values of log-likelihood and derivatives
-        model = pints.toy.ConstantModel(3)
-        times = [1, 2, 3, 4]
-        parameters = [0, 0, 0]
-        org_values = [[10.7, 3.5, 3.8],
-                      [1.1, 3.2, -1.4],
-                      [9.3, 0.0, 4.5],
-                      [1.2, -3, -10]]
-        problem = pints.MultiOutputProblem(model, times, org_values)
-        f2 = pints.GaussianKnownSigmaLogLikelihood(problem, [3.5, 1, 12])
-        log_likelihood = pints.ScaledLogLikelihood(f2)
-        # Test Gaussian_logpdf((10.7, 1.1, 9.3, 1.2)|mean=0, sigma=3.5) +
-        #      Gaussian_logpdf((3.5, 3.2, 0.0, -3)|mean=0, sigma=1) +
-        #      Gaussian_logpdf((3.8, -1.4, 4.5, -10)|mean=0, sigma=12)
-        #      = -50.5088...
-        self.assertAlmostEqual(
-            log_likelihood(parameters),
-            -50.508848609684783 / 12.0
-        )
-        l, dl = log_likelihood.evaluateS1(parameters)
-        self.assertAlmostEqual(l, -50.508848609684783 / 12.0)
-        self.assertAlmostEqual(dl[0], 1.820408163265306 / 12.0)
-        self.assertAlmostEqual(dl[1], 3.7000000000000002 / 12.0)
-        self.assertAlmostEqual(dl[2], -0.021527777777777774 / 12.0)
-
-    def test_gaussian_log_likelihoods_single_output(self):
-        # Single-output test for known/unknown noise log-likelihood methods
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, values)
 
-        model = pints.toy.LogisticModel()
-        parameters = [0.015, 500]
-        sigma = 0.1
-        times = np.linspace(0, 1000, 100)
-        values = model.simulate(parameters, times)
-        values += np.random.normal(0, sigma, values.shape)
-        problem = pints.SingleOutputProblem(model, times, values)
+        # Create a scaled and not scaled log_likelihood
+        log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, 1.5)
 
-        # Test if known/unknown give same result
-        l1 = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        l2 = pints.GaussianLogLikelihood(problem)
-        self.assertAlmostEqual(l1(parameters), l2(parameters + [sigma]))
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [3]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (1, ))
+
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -4.444444444444445)
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
+        log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, 1.5)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [3]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
 
-        # Test invalid constructors
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (1, ))
+
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -4.444444444444445)
+
+    def test_evaluateS1_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.GaussianKnownSigmaLogLikelihood(problem, 1)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [0, 0, 0]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (3, ))
+
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], 22.3)
+        self.assertAlmostEqual(deriv[1], 2 * 3.7000000000000002)
+        self.assertAlmostEqual(deriv[2], -9.3)
+
+    def test_deprecated_alias(self):
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, self.data_single)
+
+        # Create deprecated alias
+        log_likelihood = pints.KnownNoiseLogLikelihood(problem, 0.1)
+
+        # Check inheritance from current class
+        self.assertIsInstance(
+            log_likelihood, pints.GaussianKnownSigmaLogLikelihood)
+
+    def test_bad_constructor(self):
+        # Create an object with links to the model and time series
+        problem = pints.MultiOutputProblem(
+            self.model_single, self.times, self.data_single)
+
+        # Check wrong prior dimensions
         self.assertRaises(
             ValueError, pints.GaussianKnownSigmaLogLikelihood, problem, 0)
+        self.assertRaises(
+            ValueError,
+            pints.GaussianKnownSigmaLogLikelihood, problem, [0.1, 0.2])
         self.assertRaises(
             ValueError, pints.GaussianKnownSigmaLogLikelihood, problem, -1)
 
-        # known noise value checks
-        model = pints.toy.ConstantModel(1)
-        times = np.linspace(0, 10, 10)
-        values = model.simulate([2], times)
-        org_values = np.arange(10) / 5.0
-        problem = pints.SingleOutputProblem(model, times, org_values)
-        log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, 1.5)
-        self.assertAlmostEqual(log_likelihood([-1]), -21.999591968683927)
-        l, dl = log_likelihood.evaluateS1([3])
-        self.assertAlmostEqual(l, -23.777369746461702)
-        self.assertAlmostEqual(dl[0], -9.3333333333333321)
-        self.assertEqual(len(dl), 1)
 
-        # unknown noise value checks
+class TestGaussianLogLikelihood(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(3)
+
+        # Generate test data
+        cls.times = np.array([1, 2, 3, 4])
+        cls.n_times = len(cls.times)
+        cls.data_single = np.array([1, 2, 3, 4]) / 5.0
+        cls.data_multi = np.array([
+            [10.7, 3.5, 3.8],
+            [1.1, 3.2, -1.4],
+            [9.3, 0.0, 4.5],
+            [1.2, -3, -10]])
+
+        # Add noise to the data
+        np.random.seed(42)
+        cls.sigma = 0.1
+        cls.data_single += np.random.normal(
+            0, cls.sigma, cls.data_single.shape)
+        cls.data_multi += np.random.normal(0, cls.sigma, cls.data_multi.shape)
+
+    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_likelihoods with known and unknown sigma
         log_likelihood = pints.GaussianLogLikelihood(problem)
-        self.assertAlmostEqual(log_likelihood([-3, 1.5]), -47.777369746461702)
-
-        # unknown noise check sensitivity
-        model = pints.toy.ConstantModel(1)
-        times = np.linspace(0, 10, 10)
-        values = model.simulate([2], times)
-        org_values = np.arange(10) / 5.0
-        problem = pints.SingleOutputProblem(model, times, org_values)
+        log_likelihood_known = pints.GaussianKnownSigmaLogLikelihood(
+            problem, self.sigma)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2, self.sigma]
+        score = log_likelihood(test_parameters)
+
+        # Check that score between known and unknown sigma likelihoods agree
+        self.assertAlmostEqual(
+            log_likelihood(test_parameters),
+            log_likelihood_known(test_parameters[:-1]))
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -421.8952711914118)
+
+    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_likelihoods with known and unknown sigma
+        log_likelihood = pints.GaussianLogLikelihood(problem)
+        log_likelihood_known = pints.GaussianKnownSigmaLogLikelihood(
+            problem, self.sigma)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2, self.sigma]
+        score = log_likelihood(test_parameters)
+
+        # Check that score between known and unknown sigma likelihoods agree
+        self.assertAlmostEqual(
+            log_likelihood(test_parameters),
+            log_likelihood_known(test_parameters[:-1]))
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -421.8952711914118)
+
+    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_likelihoods with known and unknown sigma
         log_likelihood = pints.GaussianLogLikelihood(problem)
-        l, dl = log_likelihood.evaluateS1([7, 2.0])
-        self.assertAlmostEqual(l, -63.04585713764618)
-        self.assertAlmostEqual(dl[0], -15.25)
-        self.assertAlmostEqual(dl[1], 41.925000000000004)
-
-        # Test deprecated aliases
-        l1 = pints.KnownNoiseLogLikelihood(problem, sigma)
-        self.assertIsInstance(l1, pints.GaussianKnownSigmaLogLikelihood)
-
-        l2 = pints.UnknownNoiseLogLikelihood(problem)
-        self.assertIsInstance(l2, pints.GaussianLogLikelihood)
-
-        # test multiple output unknown noise
-        model = pints.toy.ConstantModel(3)
-        parameters = [0, 0, 0]
-        times = [1, 2, 3, 4]
-        values = model.simulate([0, 0, 0], times)
-        org_values = [[10.7, 3.5, 3.8],
-                      [1.1, 3.2, -1.4],
-                      [9.3, 0.0, 4.5],
-                      [1.2, -3, -10]]
-        problem = pints.MultiOutputProblem(model, times, org_values)
+        log_likelihood_known = pints.GaussianKnownSigmaLogLikelihood(
+            problem, self.sigma)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2, self.sigma]
+        score = log_likelihood(test_parameters)
+
+        # Check that score between known and unknown sigma likelihoods agree
+        self.assertAlmostEqual(
+            log_likelihood(test_parameters),
+            log_likelihood_known(test_parameters[:-1]))
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -421.8952711914118)
+
+    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_likelihoods with known and unknown sigma
         log_likelihood = pints.GaussianLogLikelihood(problem)
-        # Test Gaussian_logpdf((10.7, 1.1, 9.3, 1.2)|mean=0, sigma=3.5) +
-        #      Gaussian_logpdf((3.5, 3.2, 0.0, -3)|mean=0, sigma=1) +
-        #      Gaussian_logpdf((3.8, -1.4, 4.5, -10)|mean=0, sigma=12)
-        #      = -50.5088...
+        log_likelihood_known = pints.GaussianKnownSigmaLogLikelihood(
+            problem, self.sigma)
+
+        # Check that score between known and unknown sigma likelihoods agree
         self.assertAlmostEqual(
-            log_likelihood(parameters + [3.5, 1, 12]),
-            -50.508848609684783
-        )
-        l, dl = log_likelihood.evaluateS1(parameters + [3.5, 1, 12])
-        self.assertAlmostEqual(l, -50.508848609684783)
-        self.assertAlmostEqual(dl[0], 1.820408163265306)
-        self.assertAlmostEqual(dl[1], 3.7000000000000002)
-        self.assertAlmostEqual(dl[2], -0.021527777777777774)
-        self.assertAlmostEqual(dl[3], 3.6065306122448981)
-        self.assertAlmostEqual(dl[4], 27.490000000000002)
-        self.assertAlmostEqual(dl[5], -0.25425347222222222)
-
-        # test multiple output model dimensions of sensitivities
-        d = 20
-        model = pints.toy.ConstantModel(d)
-        parameters = [0 for i in range(d)]
-        times = [1, 2, 3, 4]
-        values = model.simulate(parameters, times)
-        org_values = np.ones((len(times), d))
-        extra_params = np.ones(d).tolist()
-        problem = pints.MultiOutputProblem(model, times, org_values)
+            log_likelihood([0, 0, 0, 0.1, 0.1, 0.1]),
+            log_likelihood_known([0, 0, 0]))
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 3.5, 1, 12]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -50.75425117450455)
+
+    def test_evaluateS1_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_likelihoods with known and unknown sigma
+        log_likelihood = pints.GaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [7, 2.0]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (2, ))
+
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -6.436770793841281)
+        self.assertAlmostEqual(deriv[1], 18.75242861278283)
+
+    def test_evaluateS1_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_likelihoods with known and unknown sigma
+        log_likelihood = pints.GaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [7, 2.0]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (2, ))
+
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -6.436770793841281)
+        self.assertAlmostEqual(deriv[1], 18.75242861278283)
+
+    def test_evaluateS1_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_likelihoods with known and unknown sigma
         log_likelihood = pints.GaussianLogLikelihood(problem)
-        l = log_likelihood(parameters + extra_params)
-        l1, dl = log_likelihood.evaluateS1(parameters + extra_params)
-        self.assertTrue(np.array_equal(len(dl),
-                                       len(parameters + extra_params)))
-        self.assertEqual(l, l1)
-
-    def test_gaussian_integrated_uniform_log_likelihood_single(self):
-        # Tests GaussianIntegratedUniformLogLikelihood with single output
-        # problem
-        model = pints.toy.ConstantModel(1)
-        parameters = [0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = np.asarray([1.0, -10.7, 15.5])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.GaussianIntegratedUniformLogLikelihood(
-            problem, 2, 4)
-        self.assertAlmostEqual(log_likelihood([0]), -20.441037907121299)
 
-        # test incorrect constructors
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, -1, 2)
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, 0, 0)
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, 2, 1)
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, [1, 2], [2, 3])
+        # Evaluate likelihood for test parameters
+        test_parameters = [7, 2.0]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
 
-    def test_gaussian_integrated_uniform_log_likelihood_multi(self):
-        # Tests GaussianIntegratedUniformLogLikelihood with multi output
-        # problem
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = 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]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.GaussianIntegratedUniformLogLikelihood(
-            problem, 2, 4)
-        self.assertAlmostEqual(log_likelihood(parameters), -75.443307614807225)
-
-        # test non-equal prior limits
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = 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]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.GaussianIntegratedUniformLogLikelihood(
-            problem, [1, 0, 5, 2], [2, 4, 7, 8])
-        self.assertAlmostEqual(log_likelihood(parameters), -71.62076263891457)
-
-        # test incorrect constructors
-        model = pints.toy.ConstantModel(2)
-        parameters = [0, 0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = [[1, 2],
-                  [3, 4],
-                  [5, 6]]
-        problem = pints.MultiOutputProblem(model, times, values)
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, 2, 2)
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, [1, 2, 3], [2, 4])
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, [1, 2], [2, 4, 5])
-        self.assertRaises(ValueError,
-                          pints.GaussianIntegratedUniformLogLikelihood,
-                          problem, [1, 3], [2, 2])
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
 
-    def test_known_noise_gaussian_single_S1(self):
-        # Simple tests for single known noise Gaussian log-likelihood with
-        # sensitivities.
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (2, ))
 
-        model = pints.toy.LogisticModel()
-        x = [0.015, 500]
-        sigma = 0.1
-        times = np.linspace(0, 1000, 100)
-        values = model.simulate(x, times)
-        values += np.random.normal(0, sigma, values.shape)
-        problem = pints.SingleOutputProblem(model, times, values)
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], -6.436770793841281)
+        self.assertAlmostEqual(deriv[1], 18.75242861278283)
 
-        # Test if values are correct
-        f = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        L1 = f(x)
-        L2, dL = f.evaluateS1(x)
-        self.assertEqual(L1, L2)
-        self.assertEqual(dL.shape, (2,))
+    def test_evaluateS1_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)
 
-        # Test with MultiOutputProblem
-        problem = pints.MultiOutputProblem(model, times, values)
-        f2 = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        L3 = f2(x)
-        L4, dL = f2.evaluateS1(x)
-        self.assertEqual(L3, L4)
-        self.assertEqual(L1, L3)
-        self.assertEqual(dL.shape, (2,))
-
-        # Test without noise
-        values = model.simulate(x, times)
-        problem = pints.SingleOutputProblem(model, times, values)
-        f = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        L1 = f(x)
-        L2, dL = f.evaluateS1(x)
-        self.assertEqual(L1, L2)
-        self.assertEqual(dL.shape, (2,))
-
-        # Test if zero at optimum
-        self.assertTrue(np.all(dL == 0))
-
-        # Test if positive to the left, negative to the right
-        L, dL = f.evaluateS1(x + np.array([-1e-9, 0]))
-        self.assertTrue(dL[0] > 0)
-        L, dL = f.evaluateS1(x + np.array([1e-9, 0]))
-        self.assertTrue(dL[0] < 0)
-
-        # Test if positive to the left, negative to the right
-        L, dL = f.evaluateS1(x + np.array([0, -1e-9]))
-        self.assertTrue(dL[1] > 0)
-        L, dL = f.evaluateS1(x + np.array([0, 1e-9]))
-        self.assertTrue(dL[1] < 0)
-
-        # Plot derivatives
-        if False:
-            import matplotlib.pyplot as plt
-            plt.figure()
-            r = np.linspace(x[0] * 0.95, x[0] * 1.05, 100)
-            L = []
-            dL1 = []
-            dL2 = []
-            for y in r:
-                a, b = f.evaluateS1([y, x[1]])
-                L.append(a)
-                dL1.append(b[0])
-                dL2.append(b[1])
-            plt.subplot(3, 1, 1)
-            plt.plot(r, L)
-            plt.subplot(3, 1, 2)
-            plt.plot(r, dL1)
-            plt.grid(True)
-            plt.subplot(3, 1, 3)
-            plt.plot(r, dL2)
-            plt.grid(True)
-
-            plt.figure()
-            r = np.linspace(x[1] * 0.95, x[1] * 1.05, 100)
-            L = []
-            dL1 = []
-            dL2 = []
-            for y in r:
-                a, b = f.evaluateS1([x[0], y])
-                L.append(a)
-                dL1.append(b[0])
-                dL2.append(b[1])
-            plt.subplot(3, 1, 1)
-            plt.plot(r, L)
-            plt.subplot(3, 1, 2)
-            plt.plot(r, dL1)
-            plt.grid(True)
-            plt.subplot(3, 1, 3)
-            plt.plot(r, dL2)
-            plt.grid(True)
-
-            plt.show()
-
-        # value-based tests (single output tests are above)
-        # multiple outputs
-        model = pints.toy.ConstantModel(3)
-        parameters = [0, 0, 0]
-        times = [1, 2, 3, 4]
-        values = model.simulate(parameters, times)
-        org_values = [[10.7, 3.5, 3.8],
-                      [1.1, 3.2, -1.4],
-                      [9.3, 0.0, 4.5],
-                      [1.2, -3, -10]]
-        problem = pints.MultiOutputProblem(model, times, org_values)
-        sigma = [3.5, 1, 12]
-        log_likelihood = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        # Test Gaussian_logpdf((10.7, 1.1, 9.3, 1.2)|mean=0, sigma=3.5) +
-        #      Gaussian_logpdf((3.5, 3.2, 0.0, -3)|mean=0, sigma=1) +
-        #      Gaussian_logpdf((3.8, -1.4, 4.5, -10)|mean=0, sigma=12)
-        #      = -50.5088...
-        self.assertAlmostEqual(
-            log_likelihood(parameters),
-            -50.508848609684783
-        )
-        l, dl = log_likelihood.evaluateS1(parameters)
-        self.assertAlmostEqual(l, -50.508848609684783)
-        self.assertAlmostEqual(dl[0], 1.820408163265306)
-        self.assertAlmostEqual(dl[1], 3.7000000000000002)
-        self.assertAlmostEqual(dl[2], -0.021527777777777774)
-
-    def test_student_t_log_likelihood_single(self):
-        # Single-output test for Student-t noise log-likelihood methods
-
-        model = pints.toy.ConstantModel(1)
-        parameters = [0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = np.asarray([1.0, -10.7, 15.5])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.StudentTLogLikelihood(problem)
-        # Test Student-t_logpdf(values|mean=0, df = 3, scale = 10) = -11.74..
-        self.assertAlmostEqual(log_likelihood([0, 3, 10]), -11.74010919785115)
-
-    def test_student_t_log_likelihood_multi(self):
-        # Multi-output test for Student-t noise log-likelihood methods
-
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.arange(1, 4)
-        model.simulate(parameters, times)
-        values = np.asarray([[3.5, 7.6, 8.5, 3.4],
-                             [1.1, -10.3, 15.6, 5.5],
-                             [-10, -30.5, -5, 7.6]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.StudentTLogLikelihood(problem)
-        # Test Student-t_logpdf((3.5,1.1,-10)|mean=0, df=2, scale=13) +
-        #      Student-t_logpdf((7.6,-10.3,-30.5)|mean=0, df=1, scale=8) +
-        #      Student-t_logpdf((8.5,15.6,-5)|mean=0, df=2.5, scale=13.5) +
-        #      Student-t_logpdf((3.4,5.5,7.6)|mean=0, df=3.4, scale=10.5)
-        #      = -47.83....
-        self.assertAlmostEqual(
-            log_likelihood(parameters + [2, 13, 1, 8, 2.5, 13.5, 3.4, 10.5]),
-            -47.83720347766945)
+        # Create log_likelihoods with known and unknown sigma
+        log_likelihood = pints.GaussianLogLikelihood(problem)
 
-    def test_cauchy_log_likelihood_single(self):
-        # Single-output test for Cauchy noise log-likelihood methods
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 3.5, 1, 12]
+        score, deriv = log_likelihood.evaluateS1(test_parameters)
 
-        model = pints.toy.ConstantModel(1)
-        parameters = [0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = np.asarray([1.0, -10.7, 15.5])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.CauchyLogLikelihood(problem)
-        # Test Cauchy_logpdf(values|mean=0, scale = 10) = -12.34..
-        self.assertAlmostEqual(log_likelihood([0, 10]), -12.3394986541736)
-
-    def test_cauchy_log_likelihood_multi(self):
-        # Multi-output test for Cauchy noise log-likelihood methods
-
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.arange(1, 4)
-        model.simulate(parameters, times)
-        values = np.asarray([[3.5, 7.6, 8.5, 3.4],
-                             [1.1, -10.3, 15.6, 5.5],
-                             [-10, -30.5, -5, 7.6]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.CauchyLogLikelihood(problem)
-        # Test Cauchy_logpdf((3.5,1.1,-10)|mean=0, scale=13) +
-        #      Cauchy_logpdf((7.6,-10.3,-30.5)|mean=0, scale=8) +
-        #      Cauchy_logpdf((8.5,15.6,-5)|mean=0, scale=13.5) +
-        #      Cauchy_logpdf((3.4,5.5,7.6)|mean=0, scale=10.5)
-        #      = -49.51....
-        self.assertAlmostEqual(
-            log_likelihood(parameters + [13, 8, 13.5, 10.5]),
-            -49.51182454195375)
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score, log_likelihood(test_parameters))
 
-    def test_gaussian_noise_multi(self):
-        # Multi-output test for known/unknown Gaussian noise log-likelihood
-        # methods.
+        # Check that partials have the correct shape
+        self.assertEqual(deriv.shape, (6, ))
 
-        model = pints.toy.FitzhughNagumoModel()
-        parameters = [0.5, 0.5, 0.5]
-        sigma = 0.1
-        times = np.linspace(0, 100, 100)
-        values = model.simulate(parameters, times)
-        values += np.random.normal(0, sigma, values.shape)
-        problem = pints.MultiOutputProblem(model, times, values)
+        # Check that partials are comuted correctly
+        self.assertAlmostEqual(deriv[0], 1.8053598646282394)
+        self.assertAlmostEqual(deriv[1], 6.821148214206516)
+        self.assertAlmostEqual(deriv[2], -0.06083031749704496)
+        self.assertAlmostEqual(deriv[3], 3.5690214591801)
+        self.assertAlmostEqual(deriv[4], 28.09526594310443)
+        self.assertAlmostEqual(deriv[5], -0.25285170370039783)
 
-        # Test if known/unknown give same result
-        l1 = pints.GaussianKnownSigmaLogLikelihood(problem, sigma)
-        l2 = pints.GaussianKnownSigmaLogLikelihood(problem, [sigma, sigma])
-        l3 = pints.GaussianLogLikelihood(problem)
-        self.assertAlmostEqual(
-            l1(parameters),
-            l2(parameters),
-            l3(parameters + [sigma, sigma]))
+    def test_deprecated_alias(self):
+        # Create an object with links to the model and time series
+        problem = pints.SingleOutputProblem(
+            self.model_single, self.times, self.data_single)
 
-        # Test invalid constructors
-        self.assertRaises(
-            ValueError, pints.GaussianKnownSigmaLogLikelihood, problem, 0)
-        self.assertRaises(
-            ValueError, pints.GaussianKnownSigmaLogLikelihood, problem, -1)
-        self.assertRaises(
-            ValueError, pints.GaussianKnownSigmaLogLikelihood, problem, [1])
-        self.assertRaises(
-            ValueError, pints.GaussianKnownSigmaLogLikelihood, problem,
-            [1, 2, 3, 4])
-        self.assertRaises(
-            ValueError, pints.GaussianKnownSigmaLogLikelihood, problem,
-            [1, 2, -3])
+        # Create deprecated alias
+        log_likelihood = pints.UnknownNoiseLogLikelihood(problem)
+
+        # Check inheritance from current class
+        self.assertIsInstance(
+            log_likelihood, pints.GaussianLogLikelihood)
+
+
+class TestKnownNoiseLogLikelihood(unittest.TestCase):
 
     def test_known_noise_gaussian_single_and_multi(self):
         # Tests the output of single-series against multi-series known noise
@@ -551,6 +927,376 @@ class TestLogLikelihood(unittest.TestCase):
         # Check if we get the right output
         self.assertAlmostEqual(log1(0) + log2(0), log3(0))
 
+
+class TestMultiplicativeGaussianLogLikelihood(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(3)
+
+        # Generate test data
+        cls.times = np.array([1, 2, 3, 4])
+        cls.n_times = len(cls.times)
+        cls.data_single = np.asarray([1.9, 2.1, 1.8, 2.2])
+        cls.data_multi = np.array([
+            [10.7, 3.5, 3.8],
+            [1.1, 3.2, -1.4],
+            [9.3, 0.0, 4.5],
+            [1.2, -3, -10]])
+
+    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_likelihoods with known and unknown sigma
+        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2.0, 2.0, 1.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -9.224056577298253)
+
+    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_likelihoods with known and unknown sigma
+        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2.0, 2.0, 1.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -9.224056577298253)
+
+    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_likelihoods with known and unknown sigma
+        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2.0, 2.0, 1.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -9.224056577298253)
+
+    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_likelihoods with known and unknown sigma
+        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [2.0, 2.0, 2.0, 2.0, 1.0, 2.0, 1.0, 2.0, 1.0]
+        score = log_likelihood(test_parameters)
+
+        # Check that likelihood returns expected value
+        self.assertEqual(score, -46.324126706784014)
+
+
+class TestScaledLogLikelihood(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        # Create test single output test model
+        cls.model_single = pints.toy.LogisticModel()
+        cls.model_multi = pints.toy.ConstantModel(2)
+
+        # Generate test data
+        cls.n_times = 10
+        cls.times = np.linspace(0, 1000, cls.n_times)
+        cls.data_single = cls.model_single.simulate(
+            parameters=[0.015, 500], times=cls.times)
+        cls.data_multi = cls.model_multi.simulate(
+            parameters=[1, 2], times=cls.times)
+
+    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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Evaluate likelihoods for test parameters
+        test_parameters = [0.014, 501]
+        score_not_scaled = log_likelihood_not_scaled(test_parameters)
+        score_scaled = log_likelihood_scaled(test_parameters)
+
+        # Check that unscaled likelihood returns expected value
+        self.assertEqual(int(score_not_scaled), -1897896120)
+
+        # Check that scaled likelihood returns expected value
+        self.assertAlmostEqual(score_scaled * self.n_times, score_not_scaled)
+
+    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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Evaluate likelihoods for test parameters
+        test_parameters = [0.014, 501]
+        score_not_scaled = log_likelihood_not_scaled(test_parameters)
+        score_scaled = log_likelihood_scaled(test_parameters)
+
+        # Check that unscaled likelihood returns expected value
+        self.assertEqual(int(score_not_scaled), -1897896120)
+
+        # Check that scaled likelihood returns expected value
+        self.assertAlmostEqual(score_scaled * self.n_times, score_not_scaled)
+
+    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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Evaluate likelihoods for test parameters
+        test_parameters = [0.014, 501]
+        score_not_scaled = log_likelihood_not_scaled(test_parameters)
+        score_scaled = log_likelihood_scaled(test_parameters)
+
+        # Check that unscaled likelihood returns expected value
+        self.assertEqual(int(score_not_scaled), -1897896120)
+
+        # Check that scaled likelihood returns expected value
+        self.assertAlmostEqual(score_scaled * self.n_times, score_not_scaled)
+
+    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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Evaluate likelihoods for test parameters
+        test_parameters = [2, 1]
+        score_not_scaled = log_likelihood_not_scaled(test_parameters)
+        score_scaled = log_likelihood_scaled(test_parameters)
+
+        # Check that unscaled likelihood returns expected value
+        self.assertEqual(int(score_not_scaled), -24999880)
+
+        # Check that scaled likelihood returns expected value
+        number_model_outputs = self.model_multi.n_outputs()
+        self.assertAlmostEqual(
+            score_scaled * self.n_times * number_model_outputs,
+            score_not_scaled)
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [0.014, 501]
+        score_not_scaled, deriv_not_scaled = \
+            log_likelihood_not_scaled.evaluateS1(test_parameters)
+        score_scaled, deriv_scaled = log_likelihood_scaled.evaluateS1(
+            test_parameters)
+
+        # Check that score is computed correctly
+        self.assertEqual(score_not_scaled, log_likelihood_not_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_not_scaled.shape, (2, ))
+
+        # Check that score is computed correctly
+        self.assertEqual(score_scaled, log_likelihood_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_scaled.shape, (2, ))
+
+        # Check that partials of likelihoods agree
+        unscaled_deriv = deriv_scaled * self.n_times
+        self.assertAlmostEqual(deriv_not_scaled[0], unscaled_deriv[0])
+        self.assertAlmostEqual(deriv_not_scaled[1], unscaled_deriv[1])
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [0.014, 501]
+        score_not_scaled, deriv_not_scaled = \
+            log_likelihood_not_scaled.evaluateS1(test_parameters)
+        score_scaled, deriv_scaled = log_likelihood_scaled.evaluateS1(
+            test_parameters)
+
+        # Check that score is computed correctly
+        self.assertEqual(score_not_scaled, log_likelihood_not_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_not_scaled.shape, (2, ))
+
+        # Check that score is computed correctly
+        self.assertEqual(score_scaled, log_likelihood_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_scaled.shape, (2, ))
+
+        # Check that partials of likelihoods agree
+        unscaled_deriv = deriv_scaled * self.n_times
+        self.assertAlmostEqual(deriv_not_scaled[0], unscaled_deriv[0])
+        self.assertAlmostEqual(deriv_not_scaled[1], unscaled_deriv[1])
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [0.014, 501]
+        score_not_scaled, deriv_not_scaled = \
+            log_likelihood_not_scaled.evaluateS1(test_parameters)
+        score_scaled, deriv_scaled = log_likelihood_scaled.evaluateS1(
+            test_parameters)
+
+        # Check that score is computed correctly
+        self.assertEqual(score_not_scaled, log_likelihood_not_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_not_scaled.shape, (2, ))
+
+        # Check that score is computed correctly
+        self.assertEqual(score_scaled, log_likelihood_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_scaled.shape, (2, ))
+
+        # Check that partials of likelihoods agree
+        unscaled_deriv = deriv_scaled * self.n_times
+        self.assertAlmostEqual(deriv_not_scaled[0], unscaled_deriv[0])
+        self.assertAlmostEqual(deriv_not_scaled[1], unscaled_deriv[1])
+
+    def test_evaluateS1_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 a scaled and not scaled log_likelihood
+        log_likelihood_not_scaled = pints.GaussianKnownSigmaLogLikelihood(
+            problem, sigma=0.001)
+        log_likelihood_scaled = pints.ScaledLogLikelihood(
+            log_likelihood_not_scaled)
+
+        # Compute derivatives for scaled and unscaled likelihood
+        test_parameters = [2, 1]
+        score_not_scaled, deriv_not_scaled = \
+            log_likelihood_not_scaled.evaluateS1(test_parameters)
+        score_scaled, deriv_scaled = log_likelihood_scaled.evaluateS1(
+            test_parameters)
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score_not_scaled, log_likelihood_not_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_not_scaled.shape, (2, ))
+
+        # Check that score is computed correctly
+        self.assertAlmostEqual(score_scaled, log_likelihood_scaled(
+            test_parameters))
+
+        # Check that partials have the correct shape
+        self.assertEqual(deriv_scaled.shape, (2, ))
+
+        # Check that partials of likelihoods agree
+        number_model_outputs = self.model_multi.n_outputs()
+        unscaled_deriv = deriv_scaled * self.n_times * number_model_outputs
+        self.assertAlmostEqual(deriv_not_scaled[0], unscaled_deriv[0])
+        self.assertAlmostEqual(deriv_not_scaled[1], unscaled_deriv[1])
+
+    def test_bad_constructor(self):
+        self.assertRaises(
+            ValueError, pints.ScaledLogLikelihood, self.model_single)
+
+
+class TestSumOfIndependentLogPDFs(unittest.TestCase):
+
     def test_sum_of_independent_log_pdfs(self):
 
         # Test single output
@@ -614,102 +1360,92 @@ class TestLogLikelihood(unittest.TestCase):
         y1, dy1 = l1.evaluateS1(x)
         self.assertTrue(np.all(3 * dy1 == dy))
 
-    def test_ar1(self):
-        # single outputs
-        model = pints.toy.ConstantModel(1)
-        parameters = [0]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = np.asarray([1.0, -10.7, 15.5])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.AR1LogLikelihood(problem)
-        self.assertAlmostEqual(
-            log_likelihood([0, 0.5, 5]), -19.706737485492436)
-
-        # multiple outputs
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.arange(1, 5)
-        model.simulate(parameters, times)
-        values = np.asarray([[3.5, 7.6, 8.5, 3.4],
-                             [1.1, -10.3, 15.6, 5.5],
-                             [-10, -30.5, -5, 7.6],
-                             [-12, -10.1, -4, 2.3]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.AR1LogLikelihood(problem)
-        # Test AR1Logpdf((3.5,1.1,-10, -12)|mean=0, rho=0.5, sigma=1) +
-        #      AR1Logpdf((7.6,-10.3,-30.5, -10.1)|mean=0, rho=-0.25, sigma=3) +
-        #      AR1Logpdf((8.5,15.6,-5, -4)|mean=0, rho=0.9, sigma=10) +
-        #      AR1Logpdf((3.4,5.5,7.6, 2.3)|mean=0, rho=0.0, sigma=2)
-        #      = -109.4752924909364 -93.58199 - 18.3833..
-        #        -16.4988
-        self.assertAlmostEqual(
-            log_likelihood(parameters + [0.5, 1.0,
-                                         -0.25, 3.0,
-                                         0.9, 10.0,
-                                         0.0, 2.0]),
-            -237.93936126949615)
-
-    def test_arma11(self):
-        model = pints.toy.ConstantModel(1)
-        parameters = [0]
-        times = np.asarray([1, 2, 3, 4])
-        model.simulate(parameters, times)
-        values = np.asarray([3, -4.5, 10.5, 0.3])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.ARMA11LogLikelihood(problem)
-        self.assertAlmostEqual(
-            log_likelihood([0, 0.9, -0.4, 1]), -171.53031588534171)
-
-        # multiple outputs
-        model = pints.toy.ConstantModel(4)
-        parameters = [0, 0, 0, 0]
-        times = np.arange(1, 5)
-        model.simulate(parameters, times)
-        values = np.asarray([[3.5, 7.6, 8.5, 3.4],
-                             [1.1, -10.3, 15.6, 5.5],
-                             [-10, -30.5, -5, 7.6],
-                             [-12, -10.1, -4, 2.3]])
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.ARMA11LogLikelihood(problem)
-        # ARMA1Logpdf((3.5,1.1,-10, -12)|mean=0, rho=0.5, phi=0.34 sigma=1) +
-        # ARMA1Logpdf((7.6,-10.3,-30.5, -10.1)|
-        #             mean=0, rho=-0.25, phi=0.1, sigma=3) +
-        # ARMA1Logpdf((8.5,15.6,-5, -4)|mean=0, rho=0.9, phi=0.0, sigma=10) +
-        # ARMA1Logpdf((3.4,5.5,7.6, 2.3)|mean=0, rho=0.0, phi=0.9, sigma=2)
-        #      = -116.009 -74.94 -14.32 -8.88
-        self.assertAlmostEqual(
-            log_likelihood(parameters + [0.5, 0.34, 1.0,
-                                         -0.25, 0.1, 3.0,
-                                         0.9, 0.0, 10.0,
-                                         0.0, 0.9, 2.0]),
-            -214.17034137601107)
 
-    def test_multiplicative_gaussian(self):
-        # Test single output
-        model = pints.toy.ConstantModel(1)
-        parameters = [2]
-        times = np.asarray([1, 2, 3, 4])
-        model.simulate(parameters, times)
-        values = np.asarray([1.9, 2.1, 1.8, 2.2])
-        problem = pints.SingleOutputProblem(model, times, values)
-        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+class TestStudentTLogLikelihood(unittest.TestCase):
 
-        self.assertAlmostEqual(log_likelihood(parameters + [2.0, 1.0]),
-                               -9.224056577298253)
+    @classmethod
+    def setUpClass(cls):
+        # Create test single output test model
+        cls.model_single = pints.toy.ConstantModel(1)
+        cls.model_multi = pints.toy.ConstantModel(4)
 
-        # Test multiple output
-        model = pints.toy.ConstantModel(2)
-        parameters = [1, 2]
-        times = np.asarray([1, 2, 3])
-        model.simulate(parameters, times)
-        values = np.asarray([[1.1, 0.9, 1.5], [1.5, 2.5, 2.0]]).transpose()
-        problem = pints.MultiOutputProblem(model, times, values)
-        log_likelihood = pints.MultiplicativeGaussianLogLikelihood(problem)
+        # 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.5, 7.6, 8.5, 3.4],
+            [1.1, -10.3, 15.6, 5.5],
+            [-10, -30.5, -5, 7.6]])
 
-        self.assertAlmostEqual(
-            log_likelihood(parameters + [1.0, 2.0, 1.0, 1.0]),
-            -12.176330824267543)
+    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.StudentTLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 3, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that scaled likelihood returns expected value
+        self.assertEqual(score, -11.74010919785115)
+
+    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.StudentTLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 3, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that scaled likelihood returns expected value
+        self.assertEqual(score, -11.74010919785115)
+
+    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.StudentTLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 3, 10]
+        score = log_likelihood(test_parameters)
+
+        # Check that scaled likelihood returns expected value
+        self.assertEqual(score, -11.74010919785115)
+
+    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.StudentTLogLikelihood(problem)
+
+        # Evaluate likelihood for test parameters
+        test_parameters = [0, 0, 0, 0, 2, 13, 1, 8, 2.5, 13.5, 3.4, 10.5]
+        score = log_likelihood(test_parameters)
+
+        # Check that scaled likelihood returns expected value
+        self.assertEqual(score, -47.83720347766944)
 
 
 if __name__ == '__main__':
diff --git a/pints/tests/test_toy_constant_model.py b/pints/tests/test_toy_constant_model.py
index 7b8d24b0..ac47e0e2 100755
--- a/pints/tests/test_toy_constant_model.py
+++ b/pints/tests/test_toy_constant_model.py
@@ -168,11 +168,11 @@ class TestConstantModel(unittest.TestCase):
         self.assertEqual(dy.shape, (3, 2, 2))
         dmx = np.array(
             [[[1, 0],   # dx1/dp1 = 1, dx1/dp2 = 0
-              [0, 1]],  # dx2/dp2 = 0, dx2/dp2 = 1
+              [0, 2]],  # dx2/dp2 = 0, dx2/dp2 = 2
              [[1, 0],
-              [0, 1]],
+              [0, 2]],
              [[1, 0],
-              [0, 1]]]
+              [0, 2]]]
         )
         self.assertTrue(np.all(dy == dmx))
 

pip install -e .[dev,docs]
parse_log_pytest
pytest
pytest
apt-get update apt-get install -y gcc
3.9
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.16
anyio==4.9.0
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
babel==2.17.0
beautifulsoup4==4.13.3
bleach==6.2.0
certifi==2025.1.31
cffi==1.17.1
charset-normalizer==3.4.1
cma==4.0.0
comm==0.2.2
contourpy==1.3.0
cycler==0.12.1
debugpy==1.8.13
decorator==5.2.1
defusedxml==0.7.1
docutils==0.21.2
exceptiongroup @ file:///croot/exceptiongroup_1706031385326/work
executing==2.2.0
fastjsonschema==2.21.1
flake8==7.2.0
fonttools==4.56.0
fqdn==1.5.1
h11==0.14.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
imagesize==1.4.1
importlib_metadata==8.6.1
importlib_resources==6.5.2
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==6.29.5
ipython==8.18.1
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.10.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.3.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.13
kiwisolver==1.4.7
MarkupSafe==3.0.2
matplotlib==3.9.4
matplotlib-inline==0.1.7
mccabe==0.7.0
mistune==3.1.3
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nest-asyncio==1.6.0
notebook==7.3.3
notebook_shim==0.2.4
numpy==2.0.2
overrides==7.7.0
packaging @ file:///croot/packaging_1734472117206/work
pandocfilters==1.5.1
parso==0.8.4
pexpect==4.9.0
pillow==11.1.0
-e git+https://github.com/pints-team/pints.git@9440aa00a89805b356e5982256a91302185521a6#egg=pints
platformdirs==4.3.7
pluggy @ file:///croot/pluggy_1733169602837/work
prometheus_client==0.21.1
prompt_toolkit==3.0.50
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
pycodestyle==2.13.0
pycparser==2.22
pyflakes==3.3.1
Pygments==2.19.1
pyparsing==3.2.3
pytest @ file:///croot/pytest_1738938843180/work
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
PyYAML==6.0.2
pyzmq==26.3.0
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scipy==1.13.1
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.6
Sphinx==7.4.7
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
stack-data==0.6.3
tabulate==0.9.0
terminado==0.18.1
tinycss2==1.4.0
tomli @ file:///opt/conda/conda-bld/tomli_1657175507142/work
tornado==6.4.2
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.0
uri-template==1.3.0
urllib3==2.3.0
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.13
zipp==3.21.0

name: pints
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - exceptiongroup=1.2.0=py39h06a4308_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=3.0.16=h5eee18b_0
  - packaging=24.2=py39h06a4308_0
  - pip=25.0=py39h06a4308_0
  - pluggy=1.5.0=py39h06a4308_0
  - pytest=8.3.4=py39h06a4308_0
  - python=3.9.21=he870216_1
  - readline=8.2=h5eee18b_0
  - setuptools=75.8.0=py39h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - tomli=2.0.1=py39h06a4308_0
  - tzdata=2025a=h04d1e81_0
  - wheel=0.45.1=py39h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.16
      - anyio==4.9.0
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.3.0
      - asttokens==3.0.0
      - async-lru==2.0.5
      - attrs==25.3.0
      - babel==2.17.0
      - beautifulsoup4==4.13.3
      - bleach==6.2.0
      - certifi==2025.1.31
      - cffi==1.17.1
      - charset-normalizer==3.4.1
      - cma==4.0.0
      - comm==0.2.2
      - contourpy==1.3.0
      - cycler==0.12.1
      - debugpy==1.8.13
      - decorator==5.2.1
      - defusedxml==0.7.1
      - docutils==0.21.2
      - executing==2.2.0
      - fastjsonschema==2.21.1
      - flake8==7.2.0
      - fonttools==4.56.0
      - fqdn==1.5.1
      - h11==0.14.0
      - httpcore==1.0.7
      - httpx==0.28.1
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==8.6.1
      - importlib-resources==6.5.2
      - ipykernel==6.29.5
      - ipython==8.18.1
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.10.0
      - jsonpointer==3.0.0
      - jsonschema==4.23.0
      - jsonschema-specifications==2024.10.1
      - jupyter==1.1.1
      - jupyter-client==8.6.3
      - jupyter-console==6.6.3
      - jupyter-core==5.7.2
      - jupyter-events==0.12.0
      - jupyter-lsp==2.2.5
      - jupyter-server==2.15.0
      - jupyter-server-terminals==0.5.3
      - jupyterlab==4.3.6
      - jupyterlab-pygments==0.3.0
      - jupyterlab-server==2.27.3
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.7
      - markupsafe==3.0.2
      - matplotlib==3.9.4
      - matplotlib-inline==0.1.7
      - mccabe==0.7.0
      - mistune==3.1.3
      - nbclient==0.10.2
      - nbconvert==7.16.6
      - nbformat==5.10.4
      - nest-asyncio==1.6.0
      - notebook==7.3.3
      - notebook-shim==0.2.4
      - numpy==2.0.2
      - overrides==7.7.0
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pexpect==4.9.0
      - pillow==11.1.0
      - platformdirs==4.3.7
      - prometheus-client==0.21.1
      - prompt-toolkit==3.0.50
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pure-eval==0.2.3
      - pycodestyle==2.13.0
      - pycparser==2.22
      - pyflakes==3.3.1
      - pygments==2.19.1
      - pyparsing==3.2.3
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.3.0
      - pyyaml==6.0.2
      - pyzmq==26.3.0
      - referencing==0.36.2
      - requests==2.32.3
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - rpds-py==0.24.0
      - scipy==1.13.1
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.6
      - sphinx==7.4.7
      - sphinxcontrib-applehelp==2.0.0
      - sphinxcontrib-devhelp==2.0.0
      - sphinxcontrib-htmlhelp==2.1.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==2.0.0
      - sphinxcontrib-serializinghtml==2.0.0
      - stack-data==0.6.3
      - tabulate==0.9.0
      - terminado==0.18.1
      - tinycss2==1.4.0
      - tornado==6.4.2
      - traitlets==5.14.3
      - types-python-dateutil==2.9.0.20241206
      - typing-extensions==4.13.0
      - uri-template==1.3.0
      - urllib3==2.3.0
      - wcwidth==0.2.13
      - webcolors==24.11.1
      - webencodings==0.5.1
      - websocket-client==1.8.0
      - widgetsnbextension==4.0.13
      - zipp==3.21.0
prefix: /opt/conda/envs/pints

pints/tests/test_error_measures.py::TestSumOfSquaresError::test_evaluateS1_two_dim_array_multi pints/tests/test_error_measures.py::TestSumOfSquaresError::test_evaluateS1_two_dim_array_multi_weighted pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_evaluateS1_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_evaluateS1_two_dim_array_multi pints/tests/test_toy_constant_model.py::TestConstantModel::test_derivatives
pints/tests/test_error_measures.py::TestMeanSquaredError::test_bad_constructor pints/tests/test_error_measures.py::TestMeanSquaredError::test_call_list pints/tests/test_error_measures.py::TestMeanSquaredError::test_call_one_dim_array pints/tests/test_error_measures.py::TestMeanSquaredError::test_call_two_dim_array_multi pints/tests/test_error_measures.py::TestMeanSquaredError::test_call_two_dim_array_multi_weighted pints/tests/test_error_measures.py::TestMeanSquaredError::test_call_two_dim_array_single pints/tests/test_error_measures.py::TestMeanSquaredError::test_evaluateS1_list pints/tests/test_error_measures.py::TestMeanSquaredError::test_evaluateS1_one_dim_array pints/tests/test_error_measures.py::TestMeanSquaredError::test_evaluateS1_two_dim_array_multi pints/tests/test_error_measures.py::TestMeanSquaredError::test_evaluateS1_two_dim_array_multi_weighted pints/tests/test_error_measures.py::TestMeanSquaredError::test_evaluateS1_two_dim_array_single pints/tests/test_error_measures.py::TestSumOfErrors::test_sum_of_errors pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_evaluateS1_list pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_evaluateS1_one_dim_array pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_evaluateS1_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_evaluateS1_two_dim_array_single
pints/tests/test_error_measures.py::TestNormalisedRootMeanSquaredError::test_bad_constructor pints/tests/test_error_measures.py::TestNormalisedRootMeanSquaredError::test_call_list pints/tests/test_error_measures.py::TestNormalisedRootMeanSquaredError::test_call_one_dim_array pints/tests/test_error_measures.py::TestNormalisedRootMeanSquaredError::test_call_two_dim_array_single pints/tests/test_error_measures.py::TestNormalisedRootMeanSquaredError::test_not_implemented_error pints/tests/test_error_measures.py::TestProbabilityBasedError::test_bad_constructor pints/tests/test_error_measures.py::TestProbabilityBasedError::test_call pints/tests/test_error_measures.py::TestProbabilityBasedError::test_evaluateS1 pints/tests/test_error_measures.py::TestProbabilityBasedError::test_n_parameters pints/tests/test_error_measures.py::TestRootMeanSquaredError::test_bad_constructor pints/tests/test_error_measures.py::TestRootMeanSquaredError::test_call_list pints/tests/test_error_measures.py::TestRootMeanSquaredError::test_call_one_dim_array pints/tests/test_error_measures.py::TestRootMeanSquaredError::test_call_two_dim_array_single pints/tests/test_error_measures.py::TestRootMeanSquaredError::test_not_implemented_error pints/tests/test_error_measures.py::TestSumOfSquaresError::test_bad_constructor pints/tests/test_error_measures.py::TestSumOfSquaresError::test_call_list pints/tests/test_error_measures.py::TestSumOfSquaresError::test_call_one_dim_array pints/tests/test_error_measures.py::TestSumOfSquaresError::test_call_two_dim_array_multi pints/tests/test_error_measures.py::TestSumOfSquaresError::test_call_two_dim_array_multi_weighted pints/tests/test_error_measures.py::TestSumOfSquaresError::test_call_two_dim_array_single pints/tests/test_error_measures.py::TestSumOfSquaresError::test_evaluateS1_list pints/tests/test_error_measures.py::TestSumOfSquaresError::test_evaluateS1_one_dim_array pints/tests/test_error_measures.py::TestSumOfSquaresError::test_evaluateS1_two_dim_array_single pints/tests/test_log_likelihoods.py::TestAR1LogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestAR1LogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestAR1LogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestAR1LogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestARMA11LogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestARMA11LogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestARMA11LogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestARMA11LogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestCauchyLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestCauchyLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestCauchyLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestCauchyLogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_bad_constructor_multi pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_bad_constructor_single pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_call_two_dim_array_multi_non_equal_priors pints/tests/test_log_likelihoods.py::TestGaussianIntegratedUniformLogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_bad_constructor pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_deprecated_alias pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_evaluateS1_list pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_evaluateS1_one_dim_array pints/tests/test_log_likelihoods.py::TestGaussianKnownSigmaLikelihood::test_evaluateS1_two_dim_array_single pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_deprecated_alias pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_evaluateS1_list pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_evaluateS1_one_dim_array pints/tests/test_log_likelihoods.py::TestGaussianLogLikelihood::test_evaluateS1_two_dim_array_single pints/tests/test_log_likelihoods.py::TestKnownNoiseLogLikelihood::test_known_noise_gaussian_single_and_multi pints/tests/test_log_likelihoods.py::TestMultiplicativeGaussianLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestMultiplicativeGaussianLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestMultiplicativeGaussianLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestMultiplicativeGaussianLogLikelihood::test_call_two_dim_array_single pints/tests/test_log_likelihoods.py::TestScaledLogLikelihood::test_bad_constructor pints/tests/test_log_likelihoods.py::TestSumOfIndependentLogPDFs::test_sum_of_independent_log_pdfs pints/tests/test_log_likelihoods.py::TestStudentTLogLikelihood::test_call_list pints/tests/test_log_likelihoods.py::TestStudentTLogLikelihood::test_call_one_dim_array pints/tests/test_log_likelihoods.py::TestStudentTLogLikelihood::test_call_two_dim_array_multi pints/tests/test_log_likelihoods.py::TestStudentTLogLikelihood::test_call_two_dim_array_single pints/tests/test_toy_constant_model.py::TestConstantModel::test_1d pints/tests/test_toy_constant_model.py::TestConstantModel::test_3d pints/tests/test_toy_constant_model.py::TestConstantModel::test_errors pints/tests/test_toy_constant_model.py::TestConstantModel::test_in_problem pints/tests/test_toy_constant_model.py::TestConstantModel::test_varying_numbers_of_parameters
pints/_util.py pints/toy/_constant_model.py
pints/_util.py pints/toy/_constant_model.py