diff --git a/test/likelihood_tests.py b/test/likelihood_tests.py
index a27b36df5f4aae3e53a44d9d1295f2019d1e1f67..a3d6da004a6b6690332c957acee07256c6b3d67b 100644
--- a/test/likelihood_tests.py
+++ b/test/likelihood_tests.py
@@ -4,12 +4,127 @@ from tupak.core import prior
 from tupak.core.result import Result
 import unittest
 from mock import MagicMock
+import mock
 import numpy as np
 import inspect
 import os
 import copy
 
 
+class TestLikelihoodBase(unittest.TestCase):
+
+    def setUp(self):
+        self.likelihood = tupak.core.likelihood.Likelihood()
+
+    def tearDown(self):
+        del self.likelihood
+
+    def test_base_log_likelihood(self):
+        self.assertTrue(np.isnan(self.likelihood.log_likelihood()))
+
+    def test_base_noise_log_likelihood(self):
+        self.assertTrue(np.isnan(self.likelihood.noise_log_likelihood()))
+
+    def test_base_log_likelihood_ratio(self):
+        self.assertTrue(np.isnan(self.likelihood.log_likelihood_ratio()))
+
+
+class TestAnalytical1DLikelihood(unittest.TestCase):
+
+    def setUp(self):
+        self.x = np.arange(start=0, stop=100, step=1)
+        self.y = np.arange(start=0, stop=100, step=1)
+
+        def test_func(x, parameter1, parameter2):
+            return parameter1 * x + parameter2
+
+        self.func = test_func
+        self.parameter1_value = 4
+        self.parameter2_value = 7
+        self.analytical_1d_likelihood = tupak.likelihood.Analytical1DLikelihood(x=self.x, y=self.y, func=self.func)
+        self.analytical_1d_likelihood.parameters['parameter1'] = self.parameter1_value
+        self.analytical_1d_likelihood.parameters['parameter2'] = self.parameter2_value
+
+    def tearDown(self):
+        del self.x
+        del self.y
+        del self.func
+        del self.analytical_1d_likelihood
+        del self.parameter1_value
+        del self.parameter2_value
+
+    def test_init_x(self):
+        self.assertTrue(np.array_equal(self.x, self.analytical_1d_likelihood.x))
+
+    def test_set_x_to_array(self):
+        new_x = np.arange(start=0, stop=50, step=2)
+        self.analytical_1d_likelihood.x = new_x
+        self.assertTrue(np.array_equal(new_x, self.analytical_1d_likelihood.x))
+
+    def test_set_x_to_int(self):
+        new_x = 5
+        self.analytical_1d_likelihood.x = new_x
+        expected_x = np.array([new_x])
+        self.assertTrue(np.array_equal(expected_x, self.analytical_1d_likelihood.x))
+
+    def test_set_x_to_float(self):
+        new_x = 5.3
+        self.analytical_1d_likelihood.x = new_x
+        expected_x = np.array([new_x])
+        self.assertTrue(np.array_equal(expected_x, self.analytical_1d_likelihood.x))
+
+    def test_init_y(self):
+        self.assertTrue(np.array_equal(self.y, self.analytical_1d_likelihood.y))
+
+    def test_set_y_to_array(self):
+        new_y = np.arange(start=0, stop=50, step=2)
+        self.analytical_1d_likelihood.y = new_y
+        self.assertTrue(np.array_equal(new_y, self.analytical_1d_likelihood.y))
+
+    def test_set_y_to_int(self):
+        new_y = 5
+        self.analytical_1d_likelihood.y = new_y
+        expected_y = np.array([new_y])
+        self.assertTrue(np.array_equal(expected_y, self.analytical_1d_likelihood.y))
+
+    def test_set_y_to_float(self):
+        new_y = 5.3
+        self.analytical_1d_likelihood.y = new_y
+        expected_y = np.array([new_y])
+        self.assertTrue(np.array_equal(expected_y, self.analytical_1d_likelihood.y))
+
+    def test_init_func(self):
+        self.assertEqual(self.func, self.analytical_1d_likelihood.func)
+
+    def test_set_func(self):
+        def new_func(x):
+            return x
+        with self.assertRaises(AttributeError):
+            # noinspection PyPropertyAccess
+            self.analytical_1d_likelihood.func = new_func
+
+    def test_parameters(self):
+        expected_parameters = dict(parameter1=self.parameter1_value,
+                                   parameter2=self.parameter2_value)
+        self.assertDictEqual(expected_parameters, self.analytical_1d_likelihood.parameters)
+
+    def test_n(self):
+        self.assertEqual(len(self.x), self.analytical_1d_likelihood.n)
+
+    def test_set_n(self):
+        with self.assertRaises(AttributeError):
+            # noinspection PyPropertyAccess
+            self.analytical_1d_likelihood.n = 2
+
+    def test_model_parameters(self):
+        sigma = 5
+        self.analytical_1d_likelihood.sigma = sigma
+        self.analytical_1d_likelihood.parameters['sigma'] = sigma
+        expected_model_parameters = dict(parameter1=self.parameter1_value,
+                                         parameter2=self.parameter2_value)
+        self.assertDictEqual(expected_model_parameters, self.analytical_1d_likelihood.model_parameters)
+
+
 class TestGaussianLikelihood(unittest.TestCase):
 
     def setUp(self):
@@ -20,6 +135,7 @@ class TestGaussianLikelihood(unittest.TestCase):
 
         def test_function(x, m, c):
             return m * x + c
+
         self.function = test_function
 
     def tearDown(self):
@@ -47,7 +163,7 @@ class TestGaussianLikelihood(unittest.TestCase):
         self.assertTrue(type(likelihood.sigma) == type(sigma_array))
         self.assertTrue(all(likelihood.sigma == sigma_array))
 
-    def test_unknown_float_sigma(self):
+    def test_set_sigma_None(self):
         likelihood = tupak.core.likelihood.GaussianLikelihood(
             self.x, self.y, self.function, sigma=None)
         likelihood.parameters['m'] = 2
@@ -55,37 +171,29 @@ class TestGaussianLikelihood(unittest.TestCase):
         self.assertTrue(likelihood.sigma is None)
         with self.assertRaises(TypeError):
             likelihood.log_likelihood()
+
+    def test_sigma_float(self):
+        likelihood = tupak.core.likelihood.GaussianLikelihood(
+            self.x, self.y, self.function, sigma=None)
+        likelihood.parameters['m'] = 2
+        likelihood.parameters['c'] = 0
         likelihood.parameters['sigma'] = 1
         likelihood.log_likelihood()
         self.assertTrue(likelihood.sigma is None)
 
-    def test_y(self):
-        likelihood = tupak.core.likelihood.GaussianLikelihood(
-            self.x, self.y, self.function, self.sigma)
-        self.assertTrue(all(likelihood.y == self.y))
-
-    def test_x(self):
-        likelihood = tupak.core.likelihood.GaussianLikelihood(
-            self.x, self.y, self.function, self.sigma)
-        self.assertTrue(all(likelihood.x == self.x))
-
-    def test_N(self):
-        likelihood = tupak.core.likelihood.GaussianLikelihood(
-            self.x, self.y, self.function, self.sigma)
-        self.assertTrue(likelihood.N == len(self.x))
-
 
 class TestStudentTLikelihood(unittest.TestCase):
-    
+
     def setUp(self):
         self.N = 100
-        self.nu = self.N-2
+        self.nu = self.N - 2
         self.sigma = 1
         self.x = np.linspace(0, 1, self.N)
         self.y = 2 * self.x + 1 + np.random.normal(0, self.sigma, self.N)
 
         def test_function(x, m, c):
             return m * x + c
+
         self.function = test_function
 
     def tearDown(self):
@@ -103,42 +211,61 @@ class TestStudentTLikelihood(unittest.TestCase):
         likelihood.log_likelihood()
         self.assertEqual(likelihood.sigma, self.sigma)
 
-    def test_unknown_float_nu(self):
+    def test_set_nu_none(self):
         likelihood = tupak.core.likelihood.StudentTLikelihood(
             self.x, self.y, self.function, nu=None)
         likelihood.parameters['m'] = 2
         likelihood.parameters['c'] = 0
         self.assertTrue(likelihood.nu is None)
-        with self.assertRaises((TypeError, ValueError)):
+
+    def test_log_likelihood_nu_none(self):
+        likelihood = tupak.core.likelihood.StudentTLikelihood(
+            self.x, self.y, self.function, nu=None)
+        likelihood.parameters['m'] = 2
+        likelihood.parameters['c'] = 0
+        with self.assertRaises((ValueError, TypeError)):
+            # ValueError in Python2, TypeError in Python3
             likelihood.log_likelihood()
-        likelihood.parameters['nu'] = 98
-        likelihood.log_likelihood()
-        self.assertTrue(likelihood.nu is None)
 
-    def test_y(self):
+    def test_log_likelihood_nu_zero(self):
         likelihood = tupak.core.likelihood.StudentTLikelihood(
-            self.x, self.y, self.function, self.nu, self.sigma)
-        self.assertTrue(all(likelihood.y == self.y))
+            self.x, self.y, self.function, nu=0)
+        likelihood.parameters['m'] = 2
+        likelihood.parameters['c'] = 0
+        with self.assertRaises(ValueError):
+            likelihood.log_likelihood()
 
-    def test_x(self):
+    def test_log_likelihood_nu_negative(self):
         likelihood = tupak.core.likelihood.StudentTLikelihood(
-            self.x, self.y, self.function, self.nu, self.sigma)
-        self.assertTrue(all(likelihood.x == self.x))
+            self.x, self.y, self.function, nu=-1)
+        likelihood.parameters['m'] = 2
+        likelihood.parameters['c'] = 0
+        with self.assertRaises(ValueError):
+            likelihood.log_likelihood()
 
-    def test_N(self):
+    def test_setting_nu_positive_does_not_change_class_attribute(self):
         likelihood = tupak.core.likelihood.StudentTLikelihood(
-            self.x, self.y, self.function, self.nu, self.sigma)
-        self.assertTrue(likelihood.N == len(self.x))
+            self.x, self.y, self.function, nu=None)
+        likelihood.parameters['m'] = 2
+        likelihood.parameters['c'] = 0
+        likelihood.parameters['nu'] = 98
+        self.assertTrue(likelihood.nu is None)
+
+    def test_lam(self):
+        likelihood = tupak.core.likelihood.StudentTLikelihood(
+            self.x, self.y, self.function, nu=0, sigma=0.5)
+
+        self.assertAlmostEqual(4.0, likelihood.lam)
 
 
 class TestPoissonLikelihood(unittest.TestCase):
-    
+
     def setUp(self):
         self.N = 100
         self.mu = 5
         self.x = np.linspace(0, 1, self.N)
         self.y = np.random.poisson(self.mu, self.N)
-        self.yfloat = np.copy(self.y)*1.
+        self.yfloat = np.copy(self.y) * 1.
         self.yneg = np.copy(self.y)
         self.yneg[0] = -1
 
@@ -146,10 +273,11 @@ class TestPoissonLikelihood(unittest.TestCase):
             return c
 
         def test_function_array(x, c):
-            return np.ones(len(x))*c
+            return np.ones(len(x)) * c
 
         self.function = test_function
         self.function_array = test_function_array
+        self.poisson_likelihood = tupak.core.likelihood.PoissonLikelihood(self.x, self.y, self.function)
 
     def tearDown(self):
         del self.N
@@ -160,23 +288,22 @@ class TestPoissonLikelihood(unittest.TestCase):
         del self.yneg
         del self.function
         del self.function_array
+        del self.poisson_likelihood
 
-    def test_non_integer(self):
+    def test_init_y_non_integer(self):
         with self.assertRaises(ValueError):
             tupak.core.likelihood.PoissonLikelihood(
                 self.x, self.yfloat, self.function)
 
-    def test_negative(self):
+    def test_init__y_negative(self):
         with self.assertRaises(ValueError):
             tupak.core.likelihood.PoissonLikelihood(
                 self.x, self.yneg, self.function)
 
     def test_neg_rate(self):
-        likelihood = tupak.core.likelihood.PoissonLikelihood(
-            self.x, self.y, self.function)
-        likelihood.parameters['c'] = -2
+        self.poisson_likelihood.parameters['c'] = -2
         with self.assertRaises(ValueError):
-            likelihood.log_likelihood()
+            self.poisson_likelihood.log_likelihood()
 
     def test_neg_rate_array(self):
         likelihood = tupak.core.likelihood.PoissonLikelihood(
@@ -185,24 +312,53 @@ class TestPoissonLikelihood(unittest.TestCase):
         with self.assertRaises(ValueError):
             likelihood.log_likelihood()
 
-    def test_y(self):
-        likelihood = tupak.core.likelihood.PoissonLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(all(likelihood.y == self.y))
+    def test_init_y(self):
+        self.assertTrue(np.array_equal(self.y, self.poisson_likelihood.y))
 
-    def test_x(self):
-        likelihood = tupak.core.likelihood.PoissonLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(all(likelihood.x == self.x))
+    def test_set_y_to_array(self):
+        new_y = np.arange(start=0, stop=50, step=2)
+        self.poisson_likelihood.y = new_y
+        self.assertTrue(np.array_equal(new_y, self.poisson_likelihood.y))
 
-    def test_N(self):
-        likelihood = tupak.core.likelihood.PoissonLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(likelihood.N == len(self.x))
+    def test_set_y_to_positive_int(self):
+        new_y = 5
+        self.poisson_likelihood.y = new_y
+        expected_y = np.array([new_y])
+        self.assertTrue(np.array_equal(expected_y, self.poisson_likelihood.y))
+
+    def test_set_y_to_negative_int(self):
+        with self.assertRaises(ValueError):
+            self.poisson_likelihood.y = -5
+
+    def test_set_y_to_float(self):
+        with self.assertRaises(ValueError):
+            self.poisson_likelihood.y = 5.3
+
+    def test_log_likelihood_wrong_func_return_type(self):
+        poisson_likelihood = tupak.likelihood.PoissonLikelihood(x=self.x, y=self.y, func=lambda x: 'test')
+        with self.assertRaises(ValueError):
+            poisson_likelihood.log_likelihood()
+
+    def test_log_likelihood_negative_func_return_element(self):
+        poisson_likelihood = tupak.likelihood.PoissonLikelihood(x=self.x, y=self.y, func=lambda x: np.array([3, 6, -2]))
+        with self.assertRaises(ValueError):
+            poisson_likelihood.log_likelihood()
+
+    def test_log_likelihood_zero_func_return_element(self):
+        poisson_likelihood = tupak.likelihood.PoissonLikelihood(x=self.x, y=self.y, func=lambda x: np.array([3, 6, 0]))
+        self.assertEqual(-np.inf, poisson_likelihood.log_likelihood())
+
+    def test_log_likelihood_dummy(self):
+        """ Merely tests if it goes into the right if else bracket """
+        poisson_likelihood = tupak.likelihood.PoissonLikelihood(x=self.x, y=self.y,
+                                                                func=lambda x: np.linspace(1, 100, self.N))
+        with mock.patch('numpy.sum') as m:
+            m.return_value = 1
+            self.assertEqual(0, poisson_likelihood.log_likelihood())
 
 
 class TestExponentialLikelihood(unittest.TestCase):
-    
+
     def setUp(self):
         self.N = 100
         self.mu = 5
@@ -215,10 +371,12 @@ class TestExponentialLikelihood(unittest.TestCase):
             return c
 
         def test_function_array(x, c):
-            return c*np.ones(len(x))
+            return c * np.ones(len(x))
 
         self.function = test_function
         self.function_array = test_function_array
+        self.exponential_likelihood = tupak.core.likelihood.ExponentialLikelihood(x=self.x, y=self.y,
+                                                                                  func=self.function)
 
     def tearDown(self):
         del self.N
@@ -245,20 +403,44 @@ class TestExponentialLikelihood(unittest.TestCase):
         likelihood.parameters['c'] = -1
         self.assertEqual(likelihood.log_likelihood(), -np.inf)
 
-    def test_y(self):
-        likelihood = tupak.core.likelihood.ExponentialLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(all(likelihood.y == self.y))
+    def test_init_y(self):
+        self.assertTrue(np.array_equal(self.y, self.exponential_likelihood.y))
 
-    def test_x(self):
-        likelihood = tupak.core.likelihood.ExponentialLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(all(likelihood.x == self.x))
+    def test_set_y_to_array(self):
+        new_y = np.arange(start=0, stop=50, step=2)
+        self.exponential_likelihood.y = new_y
+        self.assertTrue(np.array_equal(new_y, self.exponential_likelihood.y))
 
-    def test_N(self):
-        likelihood = tupak.core.likelihood.ExponentialLikelihood(
-            self.x, self.y, self.function)
-        self.assertTrue(likelihood.N == len(self.x))
+    def test_set_y_to_positive_int(self):
+        new_y = 5
+        self.exponential_likelihood.y = new_y
+        expected_y = np.array([new_y])
+        self.assertTrue(np.array_equal(expected_y, self.exponential_likelihood.y))
+
+    def test_set_y_to_negative_int(self):
+        with self.assertRaises(ValueError):
+            self.exponential_likelihood.y = -5
+
+    def test_set_y_to_positive_float(self):
+        new_y = 5.3
+        self.exponential_likelihood.y = new_y
+        self.assertTrue(np.array_equal(np.array([5.3]), self.exponential_likelihood.y))
+
+    def test_set_y_to_negative_float(self):
+        with self.assertRaises(ValueError):
+            self.exponential_likelihood.y = -5.3
+
+    def test_set_y_to_nd_array_with_negative_element(self):
+        with self.assertRaises(ValueError):
+            self.exponential_likelihood.y = np.array([4.3, -1.2, 4])
+
+    def test_log_likelihood_default(self):
+        """ Merely tests that it ends up at the right place in the code """
+        exponential_likelihood = tupak.core.likelihood.ExponentialLikelihood(x=self.x, y=self.y,
+                                                                             func=lambda x: np.array([4.2]))
+        with mock.patch('numpy.sum') as m:
+            m.return_value = 3
+            self.assertEqual(-3, exponential_likelihood.log_likelihood())
 
 
 if __name__ == '__main__':
diff --git a/tupak/core/likelihood.py b/tupak/core/likelihood.py
index 98a716bb4430dcb47160c5653dadb31c7f5d8c7f..b4bfaace99e47ed9a01868e9f2dae2bd7ebc3ce8 100644
--- a/tupak/core/likelihood.py
+++ b/tupak/core/likelihood.py
@@ -4,6 +4,7 @@ import inspect
 import numpy as np
 from scipy.special import gammaln
 
+
 class Likelihood(object):
 
     def __init__(self, parameters=None):
@@ -43,8 +44,89 @@ class Likelihood(object):
         return self.log_likelihood() - self.noise_log_likelihood()
 
 
-class GaussianLikelihood(Likelihood):
-    def __init__(self, x, y, function, sigma=None):
+class Analytical1DLikelihood(Likelihood):
+    """
+    A general class for 1D analytical functions. The model
+    parameters are inferred from the arguments of function
+
+    Parameters
+    ----------
+    x, y: array_like
+        The data to analyse
+    func:
+        The python function to fit to the data. Note, this must take the
+        dependent variable as its first argument. The other arguments
+        will require a prior and will be sampled over (unless a fixed
+        value is given).
+    """
+
+    def __init__(self, x, y, func):
+        parameters = self._infer_parameters_from_function(func)
+        Likelihood.__init__(self, dict.fromkeys(parameters))
+        self.x = x
+        self.y = y
+        self.__func = func
+        self.__function_keys = list(self.parameters.keys())
+
+    @property
+    def func(self):
+        """ Make func read-only """
+        return self.__func
+
+    @property
+    def model_parameters(self):
+        """ This sets up the function only parameters (i.e. not sigma for the GaussianLikelihood) """
+        return {key: self.parameters[key] for key in self.function_keys}
+
+    @property
+    def function_keys(self):
+        """ Makes function_keys read_only """
+        return self.__function_keys
+
+    @staticmethod
+    def _infer_parameters_from_function(func):
+        """ Infers the arguments of function (except the first arg which is
+            assumed to be the dep. variable)
+        """
+        parameters = inspect.getargspec(func).args
+        parameters.pop(0)
+        return parameters
+
+    @property
+    def n(self):
+        """ The number of data points """
+        return len(self.x)
+
+    @property
+    def x(self):
+        """ The independent variable. Setter assures that single numbers will be converted to arrays internally """
+        return self.__x
+
+    @x.setter
+    def x(self, x):
+        if isinstance(x, int) or isinstance(x, float):
+            x = np.array([x])
+        self.__x = x
+
+    @property
+    def y(self):
+        """ The dependent variable. Setter assures that single numbers will be converted to arrays internally """
+        return self.__y
+
+    @y.setter
+    def y(self, y):
+        if isinstance(y, int) or isinstance(y, float):
+            y = np.array([y])
+        self.__y = y
+
+    @property
+    def residual(self):
+        """ Residual of the function against the data. """
+        return self.y - self.func(self.x, **self.model_parameters)
+
+
+class GaussianLikelihood(Analytical1DLikelihood):
+    def __init__(self, x, y, func, sigma=None):
         """
         A general Gaussian likelihood for known or unknown noise - the model
         parameters are inferred from the arguments of function
@@ -53,7 +135,7 @@ class GaussianLikelihood(Likelihood):
         ----------
         x, y: array_like
             The data to analyse
-        function:
+        func:
             The python function to fit to the data. Note, this must take the
             dependent variable as its first argument. The other arguments
             will require a prior and will be sampled over (unless a fixed
@@ -66,52 +148,31 @@ class GaussianLikelihood(Likelihood):
             to that for `x` and `y`.
         """
 
-        parameters = self._infer_parameters_from_function(function)
-        Likelihood.__init__(self, dict.fromkeys(parameters))
-
-        self.x = x
-        self.y = y
+        Analytical1DLikelihood.__init__(self, x=x, y=y, func=func)
         self.sigma = sigma
-        self.function = function
 
         # Check if sigma was provided, if not it is a parameter
-        self.function_keys = list(self.parameters.keys())
         if self.sigma is None:
             self.parameters['sigma'] = None
 
-    @staticmethod
-    def _infer_parameters_from_function(func):
-        """ Infers the arguments of function (except the first arg which is
-            assumed to be the dep. variable)
-        """
-        parameters = inspect.getargspec(func).args
-        parameters.pop(0)
-        return parameters
-
-    @property
-    def N(self):
-        """ The number of data points """
-        return len(self.x)
-
     def log_likelihood(self):
-        # This checks if sigma has been set in parameters. If so, that value
-        # will be used. Otherwise, the attribute sigma is used. The logic is
-        # that if sigma is not in parameters the attribute is used which was
-        # given at init (i.e. the known sigma as either a float or array).
-        sigma = self.parameters.get('sigma', self.sigma)
-
-        # This sets up the function only parameters (i.e. not sigma)
-        model_parameters = {k: self.parameters[k] for k in self.function_keys}
+        return self.__summed_log_likelihood(sigma=self.__get_sigma())
 
-        # Calculate the residual
-        res = self.y - self.function(self.x, **model_parameters)
+    def __get_sigma(self):
+        """
+        This checks if sigma has been set in parameters. If so, that value
+        will be used. Otherwise, the attribute sigma is used. The logic is
+        that if sigma is not in parameters the attribute is used which was
+        given at init (i.e. the known sigma as either a float or array).
+        """
+        return self.parameters.get('sigma', self.sigma)
 
-        # Return the summed log likelihood
-        return -0.5 * (np.sum((res / sigma)**2)
-                       + self.N * np.log(2 * np.pi * sigma**2))
+    def __summed_log_likelihood(self, sigma):
+        return -0.5 * (np.sum((self.residual / sigma) ** 2)
+                       + self.n * np.log(2 * np.pi * sigma ** 2))
 
 
-class PoissonLikelihood(Likelihood):
+class PoissonLikelihood(Analytical1DLikelihood):
     def __init__(self, x, y, func):
         """
         A general Poisson likelihood for a rate - the model parameters are
@@ -134,82 +195,37 @@ class PoissonLikelihood(Likelihood):
             fixed value is given).
         """
 
-        parameters = self._infer_parameters_from_function(func)
-        Likelihood.__init__(self, dict.fromkeys(parameters))
-
-        self.x = x           # the dependent variable
-        self.y = y           # the counts
-
-        # check values are non-negative integers
-        if isinstance(self.y, int):
-            # convert to numpy array if passing a single integer
-            self.y = np.array([self.y])
-
-        # check array is an integer array
-        if self.y.dtype.kind not in 'ui':
-            raise ValueError("Data must be non-negative integers")
-
-        # check for non-negative integers
-        if np.any(self.y < 0):
-            raise ValueError("Data must be non-negative integers")
-
-        self.function = func
-
-        self.function_keys = list(self.parameters.keys())
-
-    @staticmethod
-    def _infer_parameters_from_function(func):
-        """ Infers the arguments of function (except the first arg which is
-            assumed to be the dep. variable)
-        """
-        parameters = inspect.getargspec(func).args
-        parameters.pop(0)
-        return parameters
+        Analytical1DLikelihood.__init__(self, x=x, y=y, func=func)
 
     @property
-    def N(self):
-        """ The number of data points """
-        return len(self.y)
+    def y(self):
+        """ Property assures that y-value is a positive integer. """
+        return self.__y
+
+    @y.setter
+    def y(self, y):
+        if not isinstance(y, np.ndarray):
+            y = np.array([y])
+        # check array is a non-negative integer array
+        if y.dtype.kind not in 'ui' or np.any(y < 0):
+            raise ValueError("Data must be non-negative integers")
+        self.__y = y
 
     def log_likelihood(self):
-        # This sets up the function only parameters (i.e. not sigma)
-        model_parameters = {k: self.parameters[k] for k in self.function_keys}
-
-        # Calculate the rate
-        rate = self.function(self.x, **model_parameters)
-
-        # sum of log factorial of counts
-        sumlogfactorial = np.sum(gammaln(self.y + 1))
-
-        # check if rate is a single value
-        if isinstance(rate, float):
-            # check rate is positive
-            if rate < 0.:
-                raise ValueError(("Poisson rate function returns a negative ",
-                                  "value!"))
-
-            if rate == 0.:
-                return -np.inf
-            else:
-                # Return the summed log likelihood
-                return (-self.N*rate + np.sum(self.y*np.log(rate))
-                        -sumlogfactorial)
-        elif isinstance(rate, np.ndarray):
-            # check rates are positive
-            if np.any(rate < 0.):
-                raise ValueError(("Poisson rate function returns a negative",
-                                  " value!"))
-
-            if np.any(rate == 0.):
-                return -np.inf
-            else:
-                return (np.sum(-rate + self.counts*np.log(rate))
-                        -sumlogfactorial)
+        rate = self.func(self.x, **self.model_parameters)
+        if not isinstance(rate, np.ndarray):
+            raise ValueError("Poisson rate function returns wrong value type! "
+                             "Is {} when it should be numpy.ndarray".format(type(rate)))
+        elif np.any(rate < 0.):
+            raise ValueError(("Poisson rate function returns a negative",
+                              " value!"))
+        elif np.any(rate == 0.):
+            return -np.inf
         else:
-            raise ValueError("Poisson rate function returns wrong value type!")
+            return np.sum(-rate + self.y * np.log(rate)) - np.sum(gammaln(self.y + 1))
 
 
-class ExponentialLikelihood(Likelihood):
+class ExponentialLikelihood(Analytical1DLikelihood):
     def __init__(self, x, y, func):
         """
         An exponential likelihood function.
@@ -226,52 +242,29 @@ class ExponentialLikelihood(Likelihood):
             value is given). The model should return the expected mean of
             the exponential distribution for each data point.
         """
-
-        parameters = self._infer_parameters_from_function(func)
-        Likelihood.__init__(self, dict.fromkeys(parameters))
-
-        self.x = x           # the dependent variable
-        self.y = y           # the observed data
-
-        # check for non-negative values
-        if np.any(self.y < 0):
-            raise ValueError("Data must be non-negative")
-
-        self.function = func
-
-        # Check if sigma was provided, if not it is a parameter
-        self.function_keys = list(self.parameters.keys())
-
-    @staticmethod
-    def _infer_parameters_from_function(func):
-        """ Infers the arguments of function (except the first arg which is
-            assumed to be the dep. variable)
-        """
-        parameters = inspect.getargspec(func).args
-        parameters.pop(0)
-        return parameters
+        Analytical1DLikelihood.__init__(self, x=x, y=y, func=func)
 
     @property
-    def N(self):
-        """ The number of data points """
-        return len(self.y)
+    def y(self):
+        """ Property assures that y-value is positive. """
+        return self.__y
+
+    @y.setter
+    def y(self, y):
+        if not isinstance(y, np.ndarray):
+            y = np.array([y])
+        if np.any(y < 0):
+            raise ValueError("Data must be non-negative")
+        self.__y = y
 
     def log_likelihood(self):
-        # This sets up the function only parameters (i.e. not sigma)
-        model_parameters = {k: self.parameters[k] for k in self.function_keys}
-
-        # Calculate the mean of the distribution
-        mu = self.function(self.x, **model_parameters)
-
-        # return -inf if any mean values are negative
+        mu = self.func(self.x, **self.model_parameters)
         if np.any(mu < 0.):
             return -np.inf
-
-        # Return the summed log likelihood
-        return -np.sum(np.log(mu) + (self.y/mu))
+        return -np.sum(np.log(mu) + (self.y / mu))
 
 
-class StudentTLikelihood(Likelihood):
+class StudentTLikelihood(Analytical1DLikelihood):
     def __init__(self, x, y, func, nu=None, sigma=1.):
         """
         A general Student's t-likelihood for known or unknown number of degrees
@@ -301,56 +294,33 @@ class StudentTLikelihood(Likelihood):
             Set the scale of the distribution. If not given then this defaults
             to 1, which specifies a standard (central) Student's t-distribution
         """
+        Analytical1DLikelihood.__init__(self, x=x, y=y, func=func)
 
-        parameters = self._infer_parameters_from_function(func)
-        Likelihood.__init__(self, dict.fromkeys(parameters))
-
-        self.x = x
-        self.y = y
         self.nu = nu
         self.sigma = sigma
-        self.function = func
 
         # Check if nu was provided, if not it is a parameter
-        self.function_keys = list(self.parameters.keys())
         if self.nu is None:
             self.parameters['nu'] = None
 
-    @staticmethod
-    def _infer_parameters_from_function(func):
-        """ Infers the arguments of function (except the first arg which is
-            assumed to be the dep. variable)
-        """
-        parameters = inspect.getargspec(func).args
-        parameters.pop(0)
-        return parameters
-
     @property
-    def N(self):
-        """ The number of data points """
-        return len(self.x)
+    def lam(self):
+        """ Converts 'scale' to 'precision' """
+        return 1. / self.sigma ** 2
 
     def log_likelihood(self):
-        # This checks if nu or sigma have been set in parameters. If so, those
-        # values will be used. Otherwise, the attribute sigma is used. The logic is
-        # that if nu is not in parameters the attribute is used which was
-        # given at init (i.e. the known nu as a float).
-        nu = self.parameters.get('nu', self.nu)
-
-        if nu <= 0.:
+        if self.__get_nu() <= 0.:
             raise ValueError("Number of degrees of freedom for Student's t-likelihood must be positive")
 
-        # This sets up the function only parameters (i.e. not sigma)
-        model_parameters = {k: self.parameters[k] for k in self.function_keys}
-
-        # Calculate the residual
-        res = self.y - self.function(self.x, **model_parameters)
+        return self.__summed_log_likelihood(self.__get_nu())
 
-        # convert "scale" to "precision"
-        lam = 1./self.sigma**2
+    def __get_nu(self):
+        """ This checks if nu or sigma have been set in parameters. If so, those
+        values will be used. Otherwise, the attribute nu is used. The logic is
+        that if nu is not in parameters the attribute is used which was
+        given at init (i.e. the known nu as a float)."""
+        return self.parameters.get('nu', self.nu)
 
-        # Return the summed log likelihood
-        return (self.N*(gammaln((nu + 1.0) / 2.0)
-                     + .5 * np.log(lam / (nu * np.pi))
-                     - gammaln(nu / 2.0))
-                     - (nu + 1.0) / 2.0 * np.sum(np.log1p(lam * res**2 / nu)))
+    def __summed_log_likelihood(self, nu):
+        return self.n * (gammaln((nu + 1.0) / 2.0) + .5 * np.log(self.lam / (nu * np.pi)) - gammaln(nu / 2.0)) \
+               - (nu + 1.0) / 2.0 * np.sum(np.log1p(self.lam * self.residual ** 2 / nu))