Mixture copula

starvine/bvcopula/copula/__init__.py

@@ -5,4 +5,5 @@
            'clayton_copula',
            'gumbel_copula',
            'marshall_olkin_copula',
+           'mixture_copula',
            ]

starvine/bvcopula/copula/mixture_copula.py

@@ -0,0 +1,91 @@
+##
+# \brief Mixture of two copula
+from __future__ import print_function, absolute_import, division
+import numpy as np
+from scipy import stats
+# STARVINE IMPORTS
+from starvine.bvcopula.copula.copula_base import CopulaBase
+from starvine.bvcopula.copula.mvtdstpack import mvtdstpack as mvt
+
+
+class MixtureCopula(CopulaBase):
+    """!
+    @brief Gaussian copula model
+    single parameter
+    \f$\theta[0] \in (-1, 1)\f$
+    """
+    def __init__(self, copula_a, wt_a, copula_b, wt_b):
+        wt_a = wt_a / (wt_a + wt_b)
+        wt_b = wt_b / (wt_a + wt_b)
+        self._name = copula_a.name + '-' + copula_b.name
+        self._thetaBounds = tuple(list(copula_a.thetaBounds) + list(copula_b.thetaBounds) + [(0.,1.), (0.,1.)])
+        self._theta0 = tuple(list(copula_a.theta0) + list(copula_b.theta0) + [wt_a, wt_b])
+        self.fittedParams = list(copula_a.fittedParams) + list(copula_b.fittedParams) + [wt_a, wt_b]
+        self._copula_a = copula_a
+        self._copula_b = copula_b
+        self._n_params_a = len(copula_a.thetaBounds)
+        self._n_params_b = len(copula_b.thetaBounds)
+
+    @property
+    def wgts(self):
+        return self._fittedParams[-2:]
+
+    @property
+    def thetaBounds(self):
+        return self._thetaBounds
+
+    @property
+    def theta0(self):
+        return self._theta0
+
+    @property
+    def name(self):
+        return self._name
+
+    @property
+    def rotation(self):
+        # mixture copula has no defined rotation
+        return 0
+
+    @CopulaBase._rotPDF
+    def _pdf(self, u, v, rotation=0, *theta):
+        theta_a = theta[0:self._n_params_a]
+        theta_b = theta[self._n_params_a:-2]
+        wgt_a = theta[-2]
+        wgt_b = theta[-1]
+        return wgt_a * self._copula_a.pdf(u, v, self._copula_a.rotation, *theta_a) + \
+            wgt_b * self._copula_b.pdf(u, v, self._copula_b.rotation, *theta_b)
+
+    @CopulaBase._rotCDF
+    def _cdf(self, u, v, rotation=0, *theta):
+        theta_a = theta[0:self._n_params_a]
+        theta_b = theta[self._n_params_a:-2]
+        wgt_a = theta[-2]
+        wgt_b = theta[-1]
+        return wgt_a * self._copula_a.cdf(u, v, self._copula_a.rotation, *theta_a) + \
+            wgt_b * self._copula_b.cdf(u, v, self._copula_b.rotation, *theta_b)
+
+    @CopulaBase._rotH
+    def _h(self, v, u, rotation=0, *theta):
+        theta_a = theta[0:self._n_params_a]
+        theta_b = theta[self._n_params_a:-2]
+        wgt_a = theta[-2]
+        wgt_b = theta[-1]
+        return wgt_a * self._copula_a.h(u, v, self._copula_a.rotation, *theta_a) + \
+            wgt_b * self._copula_b.h(u, v, self._copula_b.rotation, *theta_b)
+
+    @CopulaBase._rotHinv
+    def _hinv(self, v, u, rotation=0, *theta):
+        theta_a = theta[0:self._n_params_a]
+        theta_b = theta[self._n_params_a:-2]
+        wgt_a = theta[-2]
+        wgt_b = theta[-1]
+        return wgt_a * self._copula_a.hinv(u, v, self._copula_a.rotation, *theta_a) + \
+            wgt_b * self._copula_b.hinv(u, v, self._copula_b.rotation, *theta_b)
+
+    @CopulaBase._rotGen
+    def _gen(self, t, *theta):
+        """!
+        @brief Copula generating function
+        """
+        raise NotImplementedError

starvine/bvcopula/copula_factory.py

@@ -10,12 +10,29 @@ def validateRotation(rotation):
         raise RuntimeError("Invalid rotation.")
 
 
+def validateMixtureParams(mix_list):
+    for tuple_copula_settings in mix_list:
+        # ensure ("string", int, float) tuples in mix list
+        assert len(tuple_copula_settings) == 3
+        assert isinstance(tuple_copula_settings[0], str)  # type
+        assert isinstance(tuple_copula_settings[1], int)  # rotation
+        assert isinstance(tuple_copula_settings[2], float)  # wgt
+
+
 def Copula(copulatype, rotation=0):
     """!
     @brief Factory method. Returns a bivariate copula instance.
     @param copulatype <b>string</b> Copula type.
     @param rotation <b>int</b>  Copula rotation 1==90 deg, 2==180 deg, 3==270 deg
     """
+    if isinstance(copulatype, list):
+        if len(copulatype) == 1:
+            copulatype, rotation = copulatype[0], copulatype[1]
+        elif len(copulatype) == 2:
+            validateMixtureParams(copulatype)
+            return MixCopula(*copulatype[0], *copulatype[1])
+        else:
+            raise RuntimeError("Mixture copula of more than 2 distributions is not implemented")
     validateRotation(rotation)
     if re.match("t", copulatype):
         return t_copula.StudentTCopula(0)
@@ -34,3 +51,19 @@ def Copula(copulatype, rotation=0):
     else:
         # default
         raise RuntimeError("Copula type %s is not available" % str(copulatype))
+
+
+def MixCopula(copulatype_a, rotation_a=0, wgt_a=0.5, copulatype_b=None, rotation_b=0, wgt_b=0.5):
+    """!
+    @brief Helper method to generate a mixture distribution of two copula
+    @param copulatype_a <b>string</b> Copula A type.
+    @param rotation_a <b>int</b>  Copula rotation 1==90 deg, 2==180 deg, 3==270 deg
+    @param copulatype_b <b>string</b> Copula B type.
+    @param rotation_b <b>int</b>  Copula rotation 1==90 deg, 2==180 deg, 3==270 deg
+    @param wgt_a <b>float</b>  Copula A weight in mixture
+    @param wgt_b <b>float</b>  Copula B weight in mixture
+    """
+    wgt_a = wgt_a / (wgt_a + wgt_b)
+    wgt_b = wgt_b / (wgt_a + wgt_b)
+    return mixture_copula.MixtureCopula(Copula(copulatype_a, rotation_a), wgt_a,
+                                        Copula(copulatype_b, rotation_b), wgt_b)

starvine/bvcopula/tests/test_mixture_copula.py

@@ -0,0 +1,40 @@
+##
+# \brief Tests ability to fit and sample a mixture of 2 copula
+from __future__ import print_function, division
+import unittest
+from scipy.stats.mstats import rankdata
+# COPULA IMPORTS
+from starvine.bvcopula.copula.mixture_copula import MixtureCopula as mc
+import numpy as np
+from starvine.bvcopula import bv_plot
+import os
+pwd_ = os.path.dirname(os.path.abspath(__file__))
+from starvine.bvcopula.copula.gumbel_copula import GumbelCopula
+dataDir = pwd_ + "/data/"
+np.random.seed(123)
+tol = 4e-2
+
+
+class TestMixtureCopula(unittest.TestCase):
+    def setUp(self):
+        self._mix_copula = mc(GumbelCopula(2, [2.1]), 0.5,
+                              GumbelCopula(3, [2.7]), 0.5)
+
+    def testMixCoplulaPdf(self):
+        u = np.linspace(6.0e-2, 1.0-6e-2, 50)
+        v = np.linspace(6.0e-2, 1.0-6e-2, 50)
+        uu, vv = np.meshgrid(u, v)
+        c_pdf = self._mix_copula.pdf(uu.flatten(), vv.flatten())
+        self.assertTrue(np.all(c_pdf >= 0))
+        # plot mixture pdf
+        bv_plot.bvContourf(uu.flatten(), vv.flatten(), c_pdf, savefig="mix.png")
+
+    def testMixCoplulaCdf(self):
+        u = np.linspace(1.0e-8, 1.0-1e-8, 50)
+        v = np.linspace(1.0e-8, 1.0-1e-8, 50)
+        c_pdf = self._mix_copula.cdf(u, v)
+        self.assertTrue(np.all(c_pdf >= 0))
+        self.assertTrue(np.all(c_pdf <= 1))
+
+    def testMixCoplulaSample(self):
+        pass