pdfs: add GammaCPdf and InverseGammaCPdf incl. tests

doc/pdfs.rst

@@ -79,6 +79,14 @@ In this section, variable :math:`c` in math exressions denotes condition.
 
    .. automethod:: __init__
 
+.. autoclass:: GammaCPdf
+
+   .. automethod:: __init__
+
+.. autoclass:: InverseGammaCPdf
+
+   .. automethod:: __init__
+
 .. autoclass:: ProdCPdf
 
    .. automethod:: __init__

pybayes/__init__.pxd

@@ -9,5 +9,5 @@ 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, InverseGammaPdf, AbstractEmpPdf, EmpPdf, MarginalizedEmpPdf, ProdPdf
-from pybayes.pdfs cimport MLinGaussCPdf, LinGaussCPdf, GaussCPdf, ProdCPdf
+from pybayes.pdfs cimport MLinGaussCPdf, LinGaussCPdf, GaussCPdf, GammaCPdf, InverseGammaCPdf, ProdCPdf
 from pybayes.filters cimport Filter, KalmanFilter, ParticleFilter, MarginalizedParticleFilter

pybayes/__init__.py

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

pybayes/pdfs.pxd

@@ -149,6 +149,16 @@ cdef class GaussCPdf(CPdf):
     cdef bint _set_gauss_params(self, np.ndarray cond) except False
 
 
+cdef class GammaCPdf(CPdf):
+    cdef public double gamma
+    cdef public GammaPdf gamma_pdf
+
+
+cdef class InverseGammaCPdf(CPdf):
+    cdef public double gamma
+    cdef public InverseGammaPdf igamma_pdf
+
+
 cdef class ProdCPdf(CPdf):
     cdef readonly np.ndarray factors  # dtype=CPdf
     cdef readonly list in_indeces, out_indeces  # dtype=ndarray of ints

pybayes/pdfs.py

@@ -1213,6 +1213,88 @@ def _set_gauss_params(self, cond):
         return True
 
 
+class GammaCPdf(CPdf):
+    r"""Conditional pdf based on :class:`GammaPdf` tuned in a way to have mean
+    :math:`\mu` and standard deviation :math:`\gamma \mu`. In other words,
+    :math:`\text{GammaCpdf}(\mu, \gamma) = \text{GammaPdf}\left(
+    k = \gamma^{-2}, \theta = \gamma^2 \mu \right)`
+
+    The :math:`\gamma` parameter is specified in constructor and the :math:`\mu`
+    parameter is the conditioning variable.
+    """
+
+    def __init__(self, gamma, rv = None, cond_rv = None):
+        """Initialise conditional gamma pdf.
+
+        :param float gamma: :math:`\gamma` parameter above
+        """
+        self.gamma = gamma
+        self.gamma_pdf = GammaPdf(self.gamma**(-2.), 1.)  # theta is set later
+        self._set_rvs(1, rv, 1, cond_rv)
+
+    def mean(self, cond = None):
+        self._set_cond(cond)
+        return self.gamma_pdf.mean()
+
+    def variance(self, cond = None):
+        self._set_cond(cond)
+        return self.gamma_pdf.variance()
+
+    def eval_log(self, x, cond = None):
+        self._set_cond(cond)
+        # x is checked in self.igamma_pdf
+        return self.gamma_pdf.eval_log(x)
+
+    def sample(self, cond = None):
+        self._set_cond(cond)
+        return self.gamma_pdf.sample()
+
+    def _set_cond(self, cond):
+        self._check_cond(cond)
+        self.gamma_pdf.theta = self.gamma**2. * cond[0]
+
+
+class InverseGammaCPdf(CPdf):
+    r"""Conditional pdf based on :class:`InverseGammaPdf` tuned in a way to have mean
+    :math:`\mu` and standard deviation :math:`\gamma \mu`. In other words,
+    :math:`\text{InverseGammaCpdf}(\mu, \gamma) = \text{InverseGammaPdf}\left(
+    \alpha = \gamma^{-2} + 2, \beta = \left( \gamma^{-2} + 1 \right) \mu \right)`
+
+    The :math:`\gamma` parameter is specified in constructor and the :math:`\mu`
+    parameter is the conditioning variable.
+    """
+
+    def __init__(self, gamma, rv = None, cond_rv = None):
+        """Initialise conditional inverse gamma pdf.
+
+        :param float gamma: :math:`\gamma` parameter above
+        """
+        self.gamma = gamma
+        self.igamma_pdf = InverseGammaPdf(self.gamma**(-2.) + 2., 1.)  # beta is set later
+        self._set_rvs(1, rv, 1, cond_rv)
+
+    def mean(self, cond = None):
+        self._set_cond(cond)
+        return self.igamma_pdf.mean()
+
+    def variance(self, cond = None):
+        self._set_cond(cond)
+        return self.igamma_pdf.variance()
+
+    def eval_log(self, x, cond = None):
+        self._set_cond(cond)
+        # x is checked in self.igamma_pdf
+        return self.igamma_pdf.eval_log(x)
+
+    def sample(self, cond = None):
+        self._set_cond(cond)
+        return self.igamma_pdf.sample()
+
+    def _set_cond(self, cond):
+        self._check_cond(cond)
+        self.igamma_pdf.beta = (self.gamma**(-2.) + 1.)*cond[0]
+
+
 class ProdCPdf(CPdf):
     r"""Pdf that is formed as a chain rule of multiple conditional pdfs.
     Extends :class:`CPdf`.

pybayes/tests/test_pdfs.py

@@ -1051,6 +1051,76 @@ def test_sample(self):
         self.assertTrue(np.all(abs(emp.variance() - var) <= fuzz))
 
 
+class TestGammaCPdf(PbTestCase):
+    """Test conditional gamma pdf"""
+
+    def setUp(self):
+        self.gamma1 = pb.GammaCPdf(0.2)
+        self.gamma2 = pb.GammaCPdf(1.2)
+
+    def test_mean(self):
+        self.assertApproxEqual(self.gamma1.mean(np.array([5.])), np.array([5.]))
+        self.assertApproxEqual(self.gamma1.mean(np.array([2.])), np.array([2.]))
+        self.assertApproxEqual(self.gamma2.mean(np.array([3.])), np.array([3.]))
+        self.assertApproxEqual(self.gamma2.mean(np.array([4.])), np.array([4.]))
+
+    def test_variance(self):
+        self.assertApproxEqual(self.gamma1.variance(np.array([5.])), np.array([1.0])**2)
+        self.assertApproxEqual(self.gamma1.variance(np.array([2.])), np.array([0.4])**2)
+        self.assertApproxEqual(self.gamma2.variance(np.array([3.])), np.array([3.6])**2)
+        self.assertApproxEqual(self.gamma2.variance(np.array([4.])), np.array([4.8])**2)
+
+    def test_eval_log(self):
+        cond1 = np.array([12.])
+        cond2 = np.array([0.1])
+        equiv_gamma1 = pb.GammaPdf(25., 0.48)  # for gamma = 0.2, mu = 12
+        equiv_gamma2 = pb.GammaPdf(25., 0.004)  # for gamma = 0.2, mu = 0.1
+        for i in range(-2, 20):
+            x = np.array([i * 1.9], dtype=float)
+            # exp because assertApproxEqual has problems with -inf
+            self.assertApproxEqual(np.exp(self.gamma1.eval_log(x, cond1)),
+                                   np.exp(equiv_gamma1.eval_log(x)))
+            self.assertApproxEqual(np.exp(self.gamma1.eval_log(x, cond2)),
+                                   np.exp(equiv_gamma2.eval_log(x)))
+
+    # TODO: test sample()
+
+
+class TestInverseGammaCPdf(PbTestCase):
+    """Test conditional inverse gamma pdf"""
+
+    def setUp(self):
+        self.igamma1 = pb.InverseGammaCPdf(0.2)
+        self.igamma2 = pb.InverseGammaCPdf(1.2)
+
+    def test_mean(self):
+        self.assertApproxEqual(self.igamma1.mean(np.array([5.])), np.array([5.]))
+        self.assertApproxEqual(self.igamma1.mean(np.array([2.])), np.array([2.]))
+        self.assertApproxEqual(self.igamma2.mean(np.array([3.])), np.array([3.]))
+        self.assertApproxEqual(self.igamma2.mean(np.array([4.])), np.array([4.]))
+
+    def test_variance(self):
+        self.assertApproxEqual(self.igamma1.variance(np.array([5.])), np.array([1.0])**2)
+        self.assertApproxEqual(self.igamma1.variance(np.array([2.])), np.array([0.4])**2)
+        self.assertApproxEqual(self.igamma2.variance(np.array([3.])), np.array([3.6])**2)
+        self.assertApproxEqual(self.igamma2.variance(np.array([4.])), np.array([4.8])**2)
+
+    def test_eval_log(self):
+        cond1 = np.array([12.])
+        cond2 = np.array([0.1])
+        equiv_gamma1 = pb.InverseGammaPdf(27., 312.)  # for gamma = 0.2, mu = 12
+        equiv_gamma2 = pb.InverseGammaPdf(27., 2.6)  # for gamma = 0.2, mu = 0.1
+        for i in range(-2, 20):
+            x = np.array([i * 1.9], dtype=float)
+            # exp because assertApproxEqual has problems with -inf
+            self.assertApproxEqual(np.exp(self.igamma1.eval_log(x, cond1)),
+                                   np.exp(equiv_gamma1.eval_log(x)))
+            self.assertApproxEqual(np.exp(self.igamma1.eval_log(x, cond2)),
+                                   np.exp(equiv_gamma2.eval_log(x)))
+
+    # TODO: test sample()
+
+
 class TestProdCPdf(PbTestCase):
     """Test conditional product of pdfs"""