pdfs: introduce InverseGammaPdf, including sampling and tests

doc/pdfs.rst

@@ -44,6 +44,10 @@ Unconditional Probability Density Functions (pdfs)
 
    .. automethod:: __init__
 
+.. autoclass:: InverseGammaPdf
+
+   .. automethod:: __init__
+
 .. autoclass:: AbstractEmpPdf
 
 .. autoclass:: EmpPdf

pybayes/__init__.pxd

@@ -8,6 +8,6 @@ Cython definition file for PyBayes
 
 # TODO: cython bug(?): Cannot type from pdfs cimport Something
 from pybayes.pdfs cimport RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf
-from pybayes.pdfs cimport GammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
+from pybayes.pdfs cimport GammaPdf, InverseGammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
 from pybayes.pdfs cimport MLinGaussCPdf, LinGaussCPdf, GaussCPdf, ProdCPdf
 from pybayes.filters cimport Filter, KalmanFilter, ParticleFilter, MarginalizedParticleFilter

pybayes/__init__.py

@@ -14,6 +14,6 @@
 """
 
 from pdfs import RVComp, RV, CPdf, Pdf, UniPdf, AbstractGaussPdf, GaussPdf, LogNormPdf
-from pdfs import GammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
+from pdfs import GammaPdf, InverseGammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
 from pdfs import MLinGaussCPdf, LinGaussCPdf, GaussCPdf, ProdCPdf
 from filters import Filter, KalmanFilter, ParticleFilter, MarginalizedParticleFilter

pybayes/pdfs.pxd

@@ -84,6 +84,10 @@ cdef class GammaPdf(Pdf):
     cdef public double k, theta
 
 
+cdef class InverseGammaPdf(Pdf):
+    cdef public double alpha, beta
+
+
 cdef class AbstractEmpPdf(Pdf):
     cdef public np.ndarray weights  # dtype=double, ndims=1
     cdef public np.ndarray particles  # dtype=double, ndims=2

pybayes/pdfs.py

@@ -677,6 +677,50 @@ def sample(self, cond = None):
         return random.gamma(self.k, self.theta, size=(1,))
 
 
+class InverseGammaPdf(Pdf):
+    r"""Inverse gamma distribution with shape parameter :math:`\alpha` and scale parameter
+    :math:`\beta`. Extends :class:`Pdf`.
+
+    If random variable :math:`X \sim \text{Gamma}(k, \theta)` then :math:`Y = 1/X` will
+    have distribution :math:`\text{InverseGamma}(k, 1/\theta)` i.e.
+    :math:`\alpha = k, \beta = 1/\theta`
+
+    .. math:: f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-\alpha-1} e^{\frac{-\beta}{x}}
+    """
+
+    def __init__(self, alpha, beta, rv = None):
+        r"""Initialise Inverse gamma pdf.
+
+        :param double alpha: :math:`\alpha` shape parameter above
+        :param double beta: :math:`\beta` scale parameter above
+        """
+        assert alpha > 0.
+        assert beta > 0.
+        self.alpha = alpha
+        self.beta = beta
+        self._set_rv(1, rv)
+
+    def mean(self, cond = None):
+        if self.alpha <= 1.0:
+            raise NotImplementedError("Indeterminate form")
+        return np.array([self.beta / (self.alpha - 1.0)])
+
+    def variance(self, cond = None):
+        if self.alpha <= 2.0:
+            raise NotImplementedError("Indeterminate form")
+        return np.array([self.beta**2 / ((self.alpha - 2.0)*(self.alpha - 1.0)**2)])
+
+    def eval_log(self, x, cond = None):
+        self._check_x(x)
+        if x[0] <= 0.:
+            return float('-inf')
+        return self.alpha*math.log(self.beta) - math.lgamma(self.alpha) \
+               + (-self.alpha - 1.0)*math.log(x[0]) - self.beta/x[0]
+
+    def sample(self, cond = None):
+        return 1.0 / random.gamma(self.alpha, 1.0/self.beta, size=(1,))
+
+
 class AbstractEmpPdf(Pdf):
     r"""An abstraction of empirical probability density functions that provides common methods such
     as weight normalisation. Extends :class:`Pdf`.

pybayes/tests/test_pdfs.py

@@ -532,6 +532,70 @@ def test_sample(self):
         self.assertTrue(np.all(abs(emp2.variance() - self.gamma2.variance()) <= 1.3))
 
 
+class TestInverseGammaPdf(PbTestCase):
+    """Test Inverse gamma pdf"""
+
+    def setUp(self):
+        self.gamma1 = pb.InverseGammaPdf(2.2, 3.3)
+        self.gamma2 = pb.InverseGammaPdf(4.2, 1.3)
+
+    def test_init(self):
+        self.assertEqual(type(self.gamma1), pb.InverseGammaPdf)
+        self.assertEqual(type(self.gamma2), pb.InverseGammaPdf)
+
+    def test_invalid_init(self):
+        constructor = pb.InverseGammaPdf
+
+        # non-positive alpha:
+        self.assertRaises(AssertionError, constructor, 0., 1.)
+        self.assertRaises(AssertionError, constructor, -0.3, 2.)
+        # non-positive beta:
+        self.assertRaises(AssertionError, constructor, 2.1, 0.)
+        self.assertRaises(AssertionError, constructor, 2.7, -12.3)
+
+    def test_shape(self):
+        self.assertEqual(self.gamma1.shape(), 1)
+        self.assertEqual(self.gamma2.shape(), 1)
+
+    def test_mean(self):
+        self.assertApproxEqual(self.gamma1.mean(), np.array([2.75]))
+        self.assertApproxEqual(self.gamma2.mean(), np.array([0.40625]))
+
+    def test_variance(self):
+        self.assertApproxEqual(self.gamma1.variance(), np.array([37.81249999999996]))
+        self.assertApproxEqual(self.gamma2.variance(), np.array([0.07501775568181818]))
+
+    def test_eval_log(self):
+        exp_results = np.array([
+            [0.                , 0.                ],  # eval in -1
+            [0.                , 0.                ],  # in 0
+            [0.4628658368306967, 0.1057554711028215],  # in 1
+            [0.2622678382649004, 0.0055110998545603],  # in 2
+            [0.1241981059139586, 0.0008311145217714],  # ...
+            [0.0651241395203701, 0.0002075047340156],
+            [0.0376087194223439, 0.0000693945050975]
+        ])
+
+        x = np.zeros(1)
+        for i in range(7):
+            x[0] = i - 1.
+            self.assertApproxEqual(exp(self.gamma1.eval_log(x)), exp_results[i][0])
+            self.assertApproxEqual(exp(self.gamma2.eval_log(x)), exp_results[i][1])
+
+    @stochastic
+    def test_sample(self):
+        """Test GaussPdf.sample() mean and variance."""
+        N = 1000  # number of samples, variance is very sensible here
+        emp1 = pb.EmpPdf(self.gamma1.samples(N))  # Emipirical pdf computes sample mean and variance for us
+        emp2 = pb.EmpPdf(self.gamma2.samples(N))  # Emipirical pdf computes sample mean and variance for us
+
+        self.assertTrue(np.all(abs(emp1.mean() - self.gamma1.mean()) <= 0.4))
+        self.assertTrue(np.all(abs(emp2.mean() - self.gamma2.mean()) <= 0.025))
+
+        self.assertTrue(np.all(abs(emp1.variance() - self.gamma1.variance()) <= 32.0))
+        self.assertTrue(np.all(abs(emp2.variance() - self.gamma2.variance()) <= 0.03))
+
+
 class TestEmpPdf(PbTestCase):
     """Test empirical pdf"""