Merge pull request #530 from LinkedEarth/generalized_surrogages

Generalized surrogates

environment.yml

@@ -19,4 +19,5 @@ dependencies:
   - pytest
   - pip:
     - pyhht
+    - stochastic
     - '-e .'

pyleoclim/core/coherence.py

@@ -422,7 +422,7 @@ def plot(self, var='wtc', xlabel=None, ylabel=None, title='auto', figsize=[10, 8
     def dashboard(self, title=None, figsize=[9,12], overlap = True, phase_style = {}, 
                   line_colors = ['tab:blue','tab:orange'], savefig_settings={},
                   ts_plot_kwargs = None, wavelet_plot_kwargs= None):
-         ''' Cross-wavelet dashboard, including the two series, WTC and XWT.
+         ''' Cross-wavelet dashboard, including the two series, their WTC and XWT.
 
              Note: this design balances many considerations, and is not easily customizable.
 
@@ -595,7 +595,7 @@ def signif_test(self, number=200, method='ar1sim', seed=None, qs=[0.95], setting
         
             Number of surrogate series to create for significance testing. The default is 200.
         
-        method : {'ar1sim'}, optional
+        method : {'ar1sim','phaseran'}, optional
         
             Method through which to generate the surrogate series. The default is 'ar1sim'.
        
@@ -685,19 +685,26 @@ def signif_test(self, number=200, method='ar1sim', seed=None, qs=[0.95], setting
         (pseudo)random number at every execution, which may be important for marginal features
         in small ensembles. In general, however, we recommend increasing the
         number of draws to check that features are robust.
+        
+        One can also specifiy a different method to obtain surrogates, e.g. phase randomization:
+            
+        .. jupyter-execute::
 
+            coh.signif_test(method='phaseran').plot() 
         '''
 
         if number == 0:
             return self
-
+        
         new = self.copy()
         surr1 = self.timeseries1.surrogates(
             number=number, seed=seed, method=method, settings=settings
         )
         surr2 = self.timeseries2.surrogates(
             number=number, seed=seed, method=method, settings=settings
         )
+
+        # adjust time axis
 
         wtcs, xwts = [], []
 

pyleoclim/core/ensembleseries.py

@@ -3,6 +3,7 @@
 In addition to a MultipleSeries object, an EnsembleSeries object has the following properties:
 - All series members are assumed to share the same units and other metadata.
 - The class enables ensemble-oriented methods for computation (e.g., quantiles) and visualization (e.g., envelope plot).    
+
 """
 
 from ..utils import plotting
@@ -246,7 +247,8 @@ def quantiles(self, qs=[0.05, 0.5, 0.95], axis = 'value'):
 
         return ens_qs
 
-    def correlation(self, target=None, timespan=None, alpha=0.05, settings=None, fdr_kwargs=None, common_time_kwargs=None, mute_pbar=False, seed=None):
+    def correlation(self, target=None, timespan=None, alpha=0.05, method = 'ttest', statistic = 'pearsonr',
+                    settings=None, fdr_kwargs=None, common_time_kwargs=None, mute_pbar=False, seed=None):
         ''' Calculate the correlation between an EnsembleSeries object to a target.
 
         If the target is not specified, then the 1st member of the ensemble will be the target
@@ -271,6 +273,15 @@ def correlation(self, target=None, timespan=None, alpha=0.05, settings=None, fdr
         alpha : float
 
             The significance level (0.05 by default)
+       
+        method : str, {'ttest','built-in','ar1sim','phaseran'}
+            method for significance testing. Default is 'ttest'
+            
+        statistic : str
+            The name of the statistic used to measure the association, to be chosen from a subset of
+            https://docs.scipy.org/doc/scipy/reference/stats.html#association-correlation-tests
+            Currently supported: ['pearsonr','spearmanr','pointbiserialr','kendalltau','weightedtau']
+            The default is 'pearsonr'.
 
         settings : dict
 
@@ -332,11 +343,15 @@ def correlation(self, target=None, timespan=None, alpha=0.05, settings=None, fdr
                 series_list.append(ts)
 
             ts_ens = pyleo.EnsembleSeries(series_list)
-
-            # set an arbitrary random seed to fix the result
+              
+            # to set an arbitrary random seed to fix the result
             corr_res = ts_ens.correlation(ts, seed=2333)
             print(corr_res)
             
+            # to change the statistic: 
+            corr_res = ts_ens.correlation(ts, statistic='kendalltau', method='phaseran', settings = {'nsim':20})
+            print(corr_res)
+            
         The `print` function tabulates the output, and conveys the p-value according
         to the correlation test applied ("isospec", by default). To plot the result:
         
@@ -365,8 +380,11 @@ def correlation(self, target=None, timespan=None, alpha=0.05, settings=None, fdr
                 value2 = target.value
                 time2 = target.time
 
-            ts2 = Series(time=time2, value=value2, verbose=idx==0)
-            corr_res = ts1.correlation(ts2, timespan=timespan, settings=settings, common_time_kwargs=common_time_kwargs, seed=seed)
+            ts2 = Series(time=time2, value=value2, verbose=idx==0, auto_time_params=False)
+            corr_res = ts1.correlation(ts2, timespan=timespan, method=method,
+                                       statistic=statistic,
+                                       settings=settings, mute_pbar=True,
+                                       common_time_kwargs=common_time_kwargs, seed=seed)
             r_list.append(corr_res.r)
             signif_list.append(corr_res.signif)
             p_list.append(corr_res.p)
@@ -1193,5 +1211,4 @@ def to_array(self, axis='value', labels=True):
         else:
             return vals
 
-
 

pyleoclim/core/multipleseries.py

@@ -59,7 +59,7 @@ class MultipleSeries:
         soi = pyleo.utils.load_dataset('SOI')
         nino = pyleo.utils.load_dataset('NINO3')
         ms = soi & nino
-        ms.name = 'ENSO'
+        ms.label = 'ENSO'
         ms
                 
     '''
@@ -1628,13 +1628,10 @@ def plot(self, figsize=[10, 4],
 
         .. jupyter-execute::
 
-            import pyleoclim as pyleo
-            url = 'http://wiki.linked.earth/wiki/index.php/Special:WTLiPD?op=export&lipdid=MD982176.Stott.2004'
-            data = pyleo.Lipd(usr_path = url)
-            tslist = data.to_LipdSeriesList()
-            tslist = tslist[2:] # drop the first two series which only concerns age and depth
-            ms = pyleo.MultipleSeries(tslist)
-
+            soi = pyleo.utils.load_dataset('SOI')
+            nino = pyleo.utils.load_dataset('NINO3')
+            ms = soi & nino
+            ms.name = 'ENSO'
             fig, ax = ms.plot()
 
         '''