diff --git a/starvine/bvcopula/copula/__init__.py b/starvine/bvcopula/copula/__init__.py
index 74289e0..1110a54 100644
--- a/starvine/bvcopula/copula/__init__.py
+++ b/starvine/bvcopula/copula/__init__.py
@@ -5,4 +5,5 @@
'clayton_copula',
'gumbel_copula',
'marshall_olkin_copula',
+ 'mixture_copula',
]
diff --git a/starvine/bvcopula/copula/mixture_copula.py b/starvine/bvcopula/copula/mixture_copula.py
new file mode 100644
index 0000000..bb6428d
--- /dev/null
+++ b/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
diff --git a/starvine/bvcopula/copula_factory.py b/starvine/bvcopula/copula_factory.py
index 1905c54..b511619 100644
--- a/starvine/bvcopula/copula_factory.py
+++ b/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 string Copula type.
@param rotation int 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 string Copula A type.
+ @param rotation_a int Copula rotation 1==90 deg, 2==180 deg, 3==270 deg
+ @param copulatype_b string Copula B type.
+ @param rotation_b int Copula rotation 1==90 deg, 2==180 deg, 3==270 deg
+ @param wgt_a float Copula A weight in mixture
+ @param wgt_b float 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)
diff --git a/starvine/bvcopula/tests/test_mixture_copula.py b/starvine/bvcopula/tests/test_mixture_copula.py
new file mode 100644
index 0000000..b6d503a
--- /dev/null
+++ 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