Refactored parametrised tests (#5)

* Refactored the tests

* bump version to 0.1.1

* Fixed pep compliance

* Fixed

README.rst

@@ -3,7 +3,7 @@ lsbi: Linear Simulation Based Inference
 =======================================
 :lsbi: Linear Simulation Based Inference
 :Author: Will Handley
-:Version: 0.1.0
+:Version: 0.1.1
 :Homepage: https://github.com/handley-lab/lsbi
 :Documentation: http://lsbi.readthedocs.io/
 

lsbi/_version.py

@@ -1 +1 @@
-__version__ = '0.1.0'
+__version__ = '0.1.1'

lsbi/model.py

@@ -2,19 +2,27 @@
 import numpy as np
 from functools import cached_property
 from scipy.stats import multivariate_normal
+from numpy.linalg import inv, slogdet
+
+
+def logdet(A):
+    """log(abs(det(A)))."""
+    return slogdet(A)[1]
 
 
 class LinearModel(object):
     """A linear model.
 
+    Model M:  D = m + M theta +/- sqrt(C)
+
     Defined by:
     - Parameters: theta (n,)
     - Data: D (d,)
     - Prior mean: mu (n,)
     - Prior covariance: Sigma (n, n)
     - Data mean: m (d,)
     - Data covariance: C (d, d)
-    - Model M: D = m + M theta +/- sqrt(C)
+
 
     Parameters
     ----------
@@ -87,12 +95,12 @@ def d(self):
     @cached_property
     def invSigma(self):
         """Inverse of prior covariance."""
-        return np.linalg.inv(self.Sigma)
+        return inv(self.Sigma)
 
     @cached_property
     def invC(self):
         """Inverse of data covariance."""
-        return np.linalg.inv(self.C)
+        return inv(self.C)
 
     def likelihood(self, theta):
         """P(D|theta) as a scipy distribution object."""
@@ -104,9 +112,8 @@ def prior(self):
 
     def posterior(self, D):
         """P(theta|D) as a scipy distribution object."""
-        Sigma = np.linalg.inv(self.invSigma + self.M.T @ self.invC @ self.M)
-        mu = Sigma @ (self.invSigma @ self.mu
-                      + self.M.T @ self.invC @ (D-self.m))
+        Sigma = inv(self.invSigma + self.M.T @ self.invC @ self.M)
+        mu = self.mu + Sigma @ self.M.T @ self.invC @ (D-self.m)
         return multivariate_normal(mu, Sigma)
 
     def evidence(self):
@@ -132,7 +139,6 @@ def DKL(self, D):
         cov_q = self.prior().cov
         mu_p = self.posterior(D).mean
         mu_q = self.prior().mean
-        return 0.5 * (- np.linalg.slogdet(cov_p)[1]
-                      + np.linalg.slogdet(cov_q)[1]
-                      + np.trace(np.linalg.inv(cov_q) @ cov_p - 1)
-                      + (mu_q - mu_p) @ np.linalg.inv(cov_q) @ (mu_q - mu_p))
+        return 0.5 * (- logdet(cov_p) + logdet(cov_q)
+                      + np.trace(inv(cov_q) @ cov_p - 1)
+                      + (mu_q - mu_p) @ inv(cov_q) @ (mu_q - mu_p))

tests/test_model.py

@@ -5,136 +5,128 @@
 import pytest
 
 
-def _test_shape(model, d, n):
-    assert model.n == n
-    assert model.d == d
-    assert model.M.shape == (d, n)
-    assert model.m.shape == (d,)
-    assert model.C.shape == (d, d)
-    assert model.mu.shape == (n,)
-    assert model.Sigma.shape == (n, n)
-
-
-def test_M():
-    model = LinearModel(M=np.random.rand())
-    _test_shape(model, 1, 1)
-
-    model = LinearModel(M=np.random.rand(1))
-    _test_shape(model, 1, 1)
-
-    model = LinearModel(M=np.random.rand(1, 5))
-    _test_shape(model, 1, 5)
-
-    model = LinearModel(M=np.random.rand(3, 1))
-    _test_shape(model, 3, 1)
-
-    model = LinearModel(M=np.random.rand(3, 5))
-    _test_shape(model, 3, 5)
-
-
-def test_m_mu():
-    model = LinearModel(m=np.random.rand(), mu=np.random.rand())
-    _test_shape(model, 1, 1)
-
-    model = LinearModel(m=np.random.rand(1), mu=np.random.rand(1))
-    _test_shape(model, 1, 1)
-
-    model = LinearModel(m=np.random.rand(1), mu=np.random.rand(5))
-    _test_shape(model, 1, 5)
-
-    model = LinearModel(m=np.random.rand(3), mu=np.random.rand(1))
-    _test_shape(model, 3, 1)
-
-    model = LinearModel(m=np.random.rand(3), mu=np.random.rand(5))
-    _test_shape(model, 3, 5)
-
-
-def test_failure():
-    with pytest.raises(ValueError) as excinfo:
-        LinearModel(m=np.random.rand(5))
-    assert "Unable to determine number of parameters n" in str(excinfo.value)
-
-    with pytest.raises(ValueError) as excinfo:
-        LinearModel(mu=np.random.rand(3))
-    assert "Unable to determine data dimensions d" in str(excinfo.value)
-
-
-def random_model(d, n):
-    M = np.random.rand(d, n)
-    m = np.random.rand(d)
-    C = scipy.stats.wishart(scale=np.eye(d)).rvs()
-    mu = np.random.rand(n)
-    Sigma = scipy.stats.wishart(scale=np.eye(n)).rvs()
-    return LinearModel(M=M, m=m, C=C, mu=mu, Sigma=Sigma)
-
-
-def test_joint():
-    d = 5
-    n = 3
-    N = 100
-    model = random_model(d, n)
-    prior = model.prior()
-    evidence = model.evidence()
-    joint = model.joint()
-
-    samples_1 = prior.rvs(N)
-    samples_2 = joint.rvs(N)[:, -n:]
-
-    for i in range(n):
-        p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
+@pytest.mark.parametrize("n", np.arange(1, 6))
+@pytest.mark.parametrize("d", np.arange(1, 6))
+class TestLinearModel(object):
+
+    def random_model(self, d, n):
+        M = np.random.rand(d, n)
+        m = np.random.rand(d)
+        C = scipy.stats.wishart(scale=np.eye(d)).rvs()
+        mu = np.random.rand(n)
+        Sigma = scipy.stats.wishart(scale=np.eye(n)).rvs()
+        return LinearModel(M=M, m=m, C=C, mu=mu, Sigma=Sigma)
+
+    def _test_shape(self, model, d, n):
+        assert model.n == n
+        assert model.d == d
+        assert model.M.shape == (d, n)
+        assert model.m.shape == (d,)
+        assert model.C.shape == (d, d)
+        assert model.mu.shape == (n,)
+        assert model.Sigma.shape == (n, n)
+
+    def test_M(self, d, n):
+        model = LinearModel(M=np.random.rand())
+        self._test_shape(model, 1, 1)
+
+        model = LinearModel(M=np.random.rand(n))
+        self._test_shape(model, 1, n)
+
+        model = LinearModel(M=np.random.rand(d, n))
+        self._test_shape(model, d, n)
+
+    def test_m_mu(self, d, n):
+        model = LinearModel(m=np.random.rand(), mu=np.random.rand())
+        self._test_shape(model, 1, 1)
+
+        model = LinearModel(m=np.random.rand(), mu=np.random.rand(n))
+        self._test_shape(model, 1, n)
+
+        model = LinearModel(m=np.random.rand(d), mu=np.random.rand())
+        self._test_shape(model, d, 1)
+
+        model = LinearModel(m=np.random.rand(d), mu=np.random.rand(n))
+        self._test_shape(model, d, n)
+
+    def test_failure(self, d, n):
+        with pytest.raises(ValueError) as excinfo:
+            LinearModel(m=np.random.rand(5))
+        string = "Unable to determine number of parameters n"
+        assert string in str(excinfo.value)
+
+        with pytest.raises(ValueError) as excinfo:
+            LinearModel(mu=np.random.rand(3))
+        string = "Unable to determine data dimensions d"
+        assert string in str(excinfo.value)
+
+    def test_joint(self, d, n):
+        N = 100
+        model = self.random_model(d, n)
+        prior = model.prior()
+        evidence = model.evidence()
+        joint = model.joint()
+
+        samples_1 = prior.rvs(N)
+        samples_2 = joint.rvs(N)[:, -n:]
+
+        if n == 1:
+            samples_1 = np.atleast_2d(samples_1).T
+
+        for i in range(n):
+            p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
+            assert p > 1e-5
+
+        p = scipy.stats.kstest(prior.logpdf(samples_2),
+                               prior.logpdf(samples_1)).pvalue
         assert p > 1e-5
 
-    p = scipy.stats.kstest(prior.logpdf(samples_2),
-                           prior.logpdf(samples_1)).pvalue
-    assert p > 1e-5
+        samples_1 = evidence.rvs(N)
+        samples_2 = joint.rvs(N)[:, :d]
 
-    samples_1 = evidence.rvs(N)
-    samples_2 = joint.rvs(N)[:, :d]
+        if d == 1:
+            samples_1 = np.atleast_2d(samples_1).T
 
-    for i in range(d):
-        p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
-        assert p > 1e-5
+        for i in range(d):
+            p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
+            assert p > 1e-5
 
-    p = scipy.stats.kstest(evidence.logpdf(samples_2),
-                           evidence.logpdf(samples_1)).pvalue
-    assert p > 1e-5
-
-
-def test_likelihood_posterior():
-    d = 5
-    n = 3
-    N = 1000
-    model = random_model(d, n)
-    joint = model.joint()
-
-    samples = []
-    theta = model.prior().rvs()
-    for _ in range(N):
-        data = model.likelihood(theta).rvs()
-        theta = model.posterior(data).rvs()
-        samples.append(np.concatenate([data, theta])[:])
-    samples_1 = np.array(samples)[::100]
-    samples_2 = joint.rvs(len(samples_1))
-
-    for i in range(n+d):
-        p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
+        p = scipy.stats.kstest(evidence.logpdf(samples_2),
+                               evidence.logpdf(samples_1)).pvalue
         assert p > 1e-5
 
-    p = scipy.stats.kstest(joint.logpdf(samples_2),
-                           joint.logpdf(samples_1)).pvalue
-    assert p > 1e-5
-
+    def test_likelihood_posterior(self, d, n):
+        N = 100
+        model = self.random_model(d, n)
+        joint = model.joint()
+
+        samples = []
+        model.prior()
+        theta = np.atleast_1d(model.prior().rvs())
+        for _ in range(N):
+            data = np.atleast_1d(model.likelihood(theta).rvs())
+            theta = np.atleast_1d(model.posterior(data).rvs())
+            samples.append(np.concatenate([data, theta])[:])
+        samples_1 = np.array(samples)[::10]
+        samples_2 = joint.rvs(len(samples_1))
+
+        for i in range(n+d):
+            p = scipy.stats.kstest(samples_1[:, i], samples_2[:, i]).pvalue
+            assert p > 1e-5
+
+        p = scipy.stats.kstest(joint.logpdf(samples_2),
+                               joint.logpdf(samples_1)).pvalue
+        assert p > 1e-5
 
-def test_DKL():
-    d = 5
-    n = 3
-    N = 1000
-    model = random_model(d, n)
+    def test_DKL(self, d, n):
+        N = 1000
+        model = self.random_model(d, n)
 
-    data = model.evidence().rvs()
-    posterior = model.posterior(data)
-    prior = model.prior()
+        data = model.evidence().rvs()
+        posterior = model.posterior(data)
+        prior = model.prior()
 
-    samples = posterior.rvs(N)
-    Info = (posterior.logpdf(samples) - prior.logpdf(samples))
-    assert_allclose(Info.mean(), model.DKL(data), atol=5*Info.std()/np.sqrt(N))
+        samples = posterior.rvs(N)
+        Info = (posterior.logpdf(samples) - prior.logpdf(samples))
+        assert_allclose(Info.mean(), model.DKL(data),
+                        atol=5*Info.std()/np.sqrt(N))