From 370fd59f3b4e9e32c4a06cd7810ebfc3830786ef Mon Sep 17 00:00:00 2001 From: William Gurecky Date: Sat, 25 Jul 2020 04:58:25 -0500 Subject: [PATCH 1/2] adds initial support for mixture of two copula --- starvine/bvcopula/copula/__init__.py | 1 + starvine/bvcopula/copula/mixture_copula.py | 91 +++++++++++++++++++ starvine/bvcopula/copula_factory.py | 33 +++++++ .../bvcopula/tests/test_mixture_copula.py | 37 ++++++++ 4 files changed, 162 insertions(+) create mode 100644 starvine/bvcopula/copula/mixture_copula.py create mode 100644 starvine/bvcopula/tests/test_mixture_copula.py 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..10a5e87 --- /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.theta0) + list(copula_b.theta0) + [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..2dec675 --- /dev/null +++ b/starvine/bvcopula/tests/test_mixture_copula.py @@ -0,0 +1,37 @@ +## +# \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 +import matplotlib.pyplot as plt +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(1), 0.5, + GumbelCopula(2), 0.5) + + def testMixCoplulaPdf(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.pdf(u, v) + self.assertTrue(np.all(c_pdf >= 0)) + + 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 From 8ad82552caf2fc87f04a0bc48bdf713bdf7b9928 Mon Sep 17 00:00:00 2001 From: William Gurecky Date: Sat, 25 Jul 2020 05:55:48 -0500 Subject: [PATCH 2/2] fix bug in mixuture copula parameter init --- starvine/bvcopula/copula/mixture_copula.py | 2 +- starvine/bvcopula/tests/test_mixture_copula.py | 15 +++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/starvine/bvcopula/copula/mixture_copula.py b/starvine/bvcopula/copula/mixture_copula.py index 10a5e87..bb6428d 100644 --- a/starvine/bvcopula/copula/mixture_copula.py +++ b/starvine/bvcopula/copula/mixture_copula.py @@ -20,7 +20,7 @@ def __init__(self, copula_a, wt_a, copula_b, 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.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) diff --git a/starvine/bvcopula/tests/test_mixture_copula.py b/starvine/bvcopula/tests/test_mixture_copula.py index 2dec675..b6d503a 100644 --- a/starvine/bvcopula/tests/test_mixture_copula.py +++ b/starvine/bvcopula/tests/test_mixture_copula.py @@ -6,7 +6,7 @@ # COPULA IMPORTS from starvine.bvcopula.copula.mixture_copula import MixtureCopula as mc import numpy as np -import matplotlib.pyplot as plt +from starvine.bvcopula import bv_plot import os pwd_ = os.path.dirname(os.path.abspath(__file__)) from starvine.bvcopula.copula.gumbel_copula import GumbelCopula @@ -17,14 +17,17 @@ class TestMixtureCopula(unittest.TestCase): def setUp(self): - self._mix_copula = mc(GumbelCopula(1), 0.5, - GumbelCopula(2), 0.5) + self._mix_copula = mc(GumbelCopula(2, [2.1]), 0.5, + GumbelCopula(3, [2.7]), 0.5) def testMixCoplulaPdf(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.pdf(u, v) + 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)