Merge 9588020 into 93082f8

demos/mcorr_dcorr_mantel.ipynb

@@ -435,7 +435,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.4"
+   "version": "3.6.6"
   }
  },
  "nbformat": 4,

mgcpy/benchmarks/simulations.py

@@ -147,8 +147,12 @@ def joint_sim(num_samp, num_dim, noise=0.5):
 
     samp = (np.random.multivariate_normal(cov=sig, mean=np.zeros(2*num_dim),
                                           size=num_samp))
-    y = samp[:, num_dim: (2*num_dim)] + kappa*noise*gauss_noise
-    x = samp[:, 0:num_dim]
+    if num_dim == 1:
+        y = samp[:, (num_dim):(2*num_dim)] + kappa*noise*gauss_noise
+        x = samp[:, 0:num_dim]
+    else:      
+        y = samp[:, (num_dim+1):(2*num_dim)] + kappa*noise*gauss_noise
+        x = samp[:, 0:num_dim]
 
     return x, y
 
@@ -256,27 +260,18 @@ def spiral_sim(num_samp, num_dim, noise=0.4, low=0, high=5):
 
     :return: the data matrix and a response array
     """
-    if num_dim > 1:
-        kappa = 1
-    else:
-        kappa = 0
-    x = gen_x_unif(num_samp, num_dim, low=low, high=high)
-    rx = gen_x_unif(num_samp, num_dim, low=low, high=high)
-    ry = rx
-    z = rx
-    sig = np.diag(np.ones(shape=(num_dim)))
-    gauss_noise = (np.random.multivariate_normal(cov=sig,
-                                                 mean=np.zeros(num_dim),
-                                                 size=num_samp))
+    uniform_dist = gen_x_unif(num_samp, num_dim=1, low=low, high=high)
+    the_x = np.array(np.cos(np.pi * uniform_dist)).reshape(num_samp, 1)
+    y = uniform_dist * np.sin(np.pi * uniform_dist)
+    x = np.zeros(shape=(num_samp, num_dim))
 
-    ry = np.ones((num_samp, num_dim))
-    x[:, 0] = np.cos(z[:, 0].reshape((num_samp)) * np.pi)
-    for i in range(num_dim - 1):
-        x[:, i+1] = (x[:, i].reshape((num_samp)) * np.cos(z[:, i+1].reshape((num_samp)) * np.pi))
-        x[:, i] = (x[:, i].reshape((num_samp)) * np.sin(z[:, i+1].reshape((num_samp)) * np.pi))
-    x = rx * x
+    if num_dim > 1:
+        for i in range(num_dim - 1):
+            x[:, i] = np.squeeze((y * np.power(the_x, i)))
+    x[:, num_dim-1] = np.squeeze(uniform_dist * the_x)
 
-    y = ry * np.sin(z[:, 0].reshape((num_samp, 1)) * np.pi) + kappa*noise*ry*gauss_noise
+    gauss_noise = np.random.normal(loc=0, scale=1, size=(x.shape[0], 1))
+    y = y + noise*num_dim*gauss_noise
 
     return x, y
 

mgcpy/benchmarks/unit_tests/simulations_test.py

@@ -1,9 +1,10 @@
 import numpy as np
+import pytest
 from mgcpy.benchmarks import simulations as sims
 
 
 def test_simulations():
-    num_samps = 100
+    num_samps = 1000
     num_dim1 = 1
     num_dim2 = 300
     independent = True

mgcpy/independence_tests/abstract_class.py

@@ -7,7 +7,7 @@
 
 import numpy as np
 from scipy.spatial.distance import pdist, squareform
-from scipy.stats import t
+from scipy.stats import kendalltau, pearsonr, spearmanr, t
 
 
 def EUCLIDEAN_DISTANCE(x): return squareform(pdist(x, metric='euclidean'))
@@ -129,6 +129,15 @@ def p_value(self, matrix_X, matrix_Y, replication_factor=1000):
                                 "p_local_correlation_matrix": p_local_correlation_matrix,
                                 "local_correlation_matrix": local_correlation_matrix,
                                 "optimal_scale": independence_test_metadata["optimal_scale"]}
+        elif self.get_name() == "kendall":
+            p_value = kendalltau(matrix_X, matrix_Y)[1]
+            p_value_metadata = {}
+        elif self.get_name() == "spearman":
+            p_value = spearmanr(matrix_X, matrix_Y)[1]
+            p_value_metadata = {}
+        elif self.get_name() == "pearson":
+            p_value = pearsonr(matrix_X, matrix_Y)[1]
+            p_value_metadata = {}
         else:
             # estimate the null by a permutation test
             test_stats_null = np.zeros(replication_factor)

mgcpy/independence_tests/hhg.py

@@ -9,7 +9,7 @@ def __init__(self, compute_distance_matrix=None):
         :param compute_distance_matrix: a function to compute the pairwise distance matrix, given a data matrix
         :type compute_distance_matrix: FunctionType or callable()
         """
-        IndependenceTest.__init__(compute_distance_matrix)
+        IndependenceTest.__init__(self, compute_distance_matrix)
         self.which_test = "hhg"
 
     def test_statistic(self, matrix_X, matrix_Y):

mgcpy/independence_tests/kendall_spearman.py

@@ -42,10 +42,6 @@ def test_statistic(self, matrix_X, matrix_Y):
         >>> kendall_spearman = KendallSpearman()
         >>> kendall_spearman_stat = kendall_spearman.test_statistic(X, Y)
         """
-        if matrix_X is None:
-            matrix_X = self.matrix_X
-        if matrix_Y is None:
-            matrix_Y = self.matrix_Y
         assert matrix_X.shape[1] == 1, "Data matrix should be (n, 1) shape"
         assert matrix_Y.shape[1] == 1, "Data matrix should be (n, 1) shape"
 

mgcpy/independence_tests/unit_tests/hhg_test.py

@@ -0,0 +1,28 @@
+import mgcpy.benchmarks.simulations as sims
+import numpy as np
+import pytest
+from mgcpy.independence_tests.hhg import HHG
+
+
+def test_hhg():
+    # Against a randomly defined data set
+    X = np.array([1.1728, 2.4941, 2.4101, 0.1814, 1.1978, 1.5806, 1.2504,
+                  1.9706, 1.8839, 0.8760])[:, np.newaxis]
+    Y = np.array([3.2311, 12.1113, 11.1350, 1.1989, 3.3127, 4.8580, 3.4917,
+                  7.1748, 6.5792, 2.4012])[:, np.newaxis]
+    hhg = HHG()
+    test_stat = hhg.test_statistic(X, Y)[0]
+
+    assert np.round(test_stat, decimals=2) == 411.88
+
+    # Against linear simulations
+    np.random.seed(0)
+    X, Y = sims.linear_sim(100, 1)
+    test_stat = hhg.test_statistic(X, Y)[0]
+
+    assert np.round(test_stat, decimals=2) == 28986.52
+
+    X, Y = sims.linear_sim(100, 1, noise=0)
+    test_stat = hhg.test_statistic(X, Y)[0]
+
+    assert np.round(test_stat, decimals=2) == 950600.00

mgcpy/independence_tests/unit_tests/kendall_spearman_test.py

@@ -0,0 +1,35 @@
+import mgcpy.benchmarks.simulations as sims
+import numpy as np
+import pytest
+from mgcpy.independence_tests.kendall_spearman import KendallSpearman
+
+
+def test_kendall_spearman():
+    # Against a randomly defined data set
+    X = np.array([1.1728, 2.4941, 2.4101, 0.1814, 1.1978, 1.5806, 1.2504,
+                  1.9706, 1.8839, 0.8760])[:, np.newaxis]
+    Y = np.array([3.2311, 12.1113, 11.1350, 1.1989, 3.3127, 4.8580, 3.4917,
+                  7.1748, 6.5792, 2.4012])[:, np.newaxis]
+    kspear = KendallSpearman(None)
+    kspear2 = KendallSpearman(None, 'spearman')
+
+    test_stat1 = kspear.test_statistic(X, Y)[0]
+    test_stat2 = kspear2.test_statistic(X, Y)[0]
+
+    assert np.round(test_stat1, decimals=2) == 1.00
+    assert np.round(test_stat2, decimals=2) == 1.00
+
+    # Against linear simulations
+    np.random.seed(0)
+    X, Y = sims.linear_sim(100, 1)
+    kspear = KendallSpearman(None)
+    kspear2 = KendallSpearman(None, 'spearman')
+
+    assert kspear.get_name() == 'kendall'
+    assert kspear2.get_name() == 'spearman'
+
+    test_stat1 = kspear.test_statistic(X, Y)[0]
+    test_stat2 = kspear2.test_statistic(X, Y)[0]
+
+    assert np.round(test_stat1, decimals=2) == 0.33
+    assert np.round(test_stat2, decimals=2) == 0.48

mgcpy/independence_tests/unit_tests/rv_corr_test.py

@@ -1,9 +1,8 @@
 import mgcpy.benchmarks.simulations as sims
 import numpy as np
+import pytest
 from mgcpy.independence_tests.rv_corr import RVCorr
 
-# from scipy.io import savemat
-
 
 def test_local_corr():
     # Against a randomly defined data set
@@ -18,17 +17,13 @@ def test_local_corr():
     test_stat1 = rvcorr.test_statistic(X, Y)[0]
     test_stat2 = rvcorr2.test_statistic(X, Y)[0]
     test_stat3 = rvcorr3.test_statistic(X, Y)[0]
-
     assert np.round(test_stat1, decimals=2) == 0.90
     assert np.round(test_stat2, decimals=2) == 0.95
     assert np.round(test_stat3, decimals=2) == 0.90
 
-    del X, Y, rvcorr, rvcorr2, rvcorr3, test_stat1, test_stat2, test_stat3
-
     # Against linear simulations
     np.random.seed(0)
     X, Y = sims.linear_sim(100, 1)
-    # savemat('distance matrix data', {'X' : X, 'Y' : Y})
     rvcorr = RVCorr(None)
     rvcorr2 = RVCorr(None, 'pearson')
     rvcorr3 = RVCorr(None, 'cca')
@@ -40,10 +35,6 @@ def test_local_corr():
     test_stat1 = rvcorr.test_statistic(X, Y)[0]
     test_stat2 = rvcorr2.test_statistic(X, Y)[0]
     test_stat3 = rvcorr3.test_statistic(X, Y)[0]
-
     assert np.round(test_stat1, decimals=2) == 0.24
     assert np.round(test_stat2, decimals=2) == 0.49
     assert np.round(test_stat3, decimals=2) == 0.24
-
-
-test_local_corr()