Surrogate Type Extension for MI Kraskov Estimator

demos/python/example10MIContinuousDataKraskovRotatedSurrogates.py

@@ -0,0 +1,159 @@
+##
+##  Java Information Dynamics Toolkit (JIDT)
+##  Copyright (C) 2012, Joseph T. Lizier
+##  
+##  This program is free software: you can redistribute it and/or modify
+##  it under the terms of the GNU General Public License as published by
+##  the Free Software Foundation, either version 3 of the License, or
+##  (at your option) any later version.
+##  
+##  This program is distributed in the hope that it will be useful,
+##  but WITHOUT ANY WARRANTY; without even the implied warranty of
+##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+##  GNU General Public License for more details.
+##  
+##  You should have received a copy of the GNU General Public License
+##  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+##
+
+# = Example 10 - Comparison of Rotated and Shuffled Surrogates for MI using Kraskov estimator  =
+
+# This example by Donovan Rynne, 2024
+
+# Autocorrelated samples result in a bias that is lost in surrogates if created through permutation.
+
+# This example generates random AR(1) data where the source and destination are uncorelated, i.e. generated under the null distribution.
+# It is therefore expected that the p-values are uniformly distributed between 0-1.
+# This demo tests that assumption using the two types surrogate distributions and plots the cumulative probability distribution
+
+from jpype import *
+import random
+import math
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+from tqdm import tqdm
+
+def generate_ar1_data(phi, n, initial_value = random.gauss(0, 1)):
+
+    # Initialize the time series with the initial value
+    ar1_data = [initial_value]
+
+    # Generate AR(1) data
+    for _ in range(1, n):
+        # Generate white noise (epsilon_t)
+        epsilon_t = random.gauss(0, 1)
+        # Calculate the next value in the time series
+        next_value = phi * ar1_data[-1] + epsilon_t
+        ar1_data.append(next_value)
+
+    return ar1_data
+
+def plot_cdf(pvalue):
+
+    pvalue_sorted = np.sort(pvalue)
+    cdf = np.arange(1, len(pvalue_sorted) + 1) / len(pvalue_sorted)
+
+    plt.plot(pvalue_sorted, cdf)
+    plt.axline([0, 0], slope=1, color='red', linestyle='--')  # Add diagonal line for reference
+
+# Change location of jar to match yours (we assume script is called from demos/python):
+jarLocation = os.path.join(os.getcwd(), '..', '..', "infodynamics.jar");
+if (not(os.path.isfile(jarLocation))):
+	exit("infodynamics.jar not found (expected at " + os.path.abspath(jarLocation) + ") - are you running from demos/python?")
+# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
+startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
+
+# Generate some autocorrelated data.
+numObservations = 1000
+# AR(1) correlation parameter, phi
+phi = 0.9
+# Number of realisations
+numRealisations = 1000
+# Number of surrogates
+numSurrogates = 100
+
+# Create the MI Calculator
+calcClass = JPackage("infodynamics.measures.continuous.kraskov").MutualInfoCalculatorMultiVariateKraskov1
+# calcClass = infodynamics_package.measures.continuous.kraskov.MutualInfoCalculatorMultiVariateKraskov1
+calc = calcClass()
+
+# Initialise pvalues array for shuffled and rotated surrogates
+shuffled_pvalues = []
+rotated_pvalues = []
+
+print(f"Generating P-Values using shuffled surrogates for {numRealisations} realisations:")
+for realisation in tqdm(range(numRealisations)):
+
+    # Source array of random autocorrelated data:
+    sourceArray = generate_ar1_data(phi, numObservations)
+    # Destination array of random autocorrelated data
+    destArray = generate_ar1_data(phi, numObservations)
+
+    # Set the surrogate type property to shuffle (default, not neccesary)
+    calc.setProperty("SURROGATE_TYPE", "SHUFFLE")
+    # initialise the calculator
+    calc.initialise()
+
+    #0. load data
+    source = JArray(JDouble, 1)(sourceArray)
+    destination = JArray(JDouble, 1)(destArray)
+
+    # 4. Supply the sample data:
+    calc.setObservations(source, destination)
+    # 5. Compute the estimate:
+    result = calc.computeAverageLocalOfObservations()
+    # 6. Compute the (statistical significance via) null distribution empirically:
+    measDist = calc.computeSignificance(numSurrogates)
+    shuffled_pvalues.append(measDist.pValue)
+
+# Repeat for rotated surrogates
+print(f"Generating P-Values using rotated surrogates for {numRealisations} realisations:")
+for realisation in tqdm(range(numRealisations)):
+    # Source array of random autocorrelated data:
+    sourceArray = generate_ar1_data(phi, numObservations)
+    # Destination array of random autocorrelated data
+    destArray = generate_ar1_data(phi, numObservations)
+
+    # Set the surrogate type property to rotate
+    calc.setProperty("SURROGATE_TYPE", "ROTATE") 
+    #initialise the calculator
+    calc.initialise()
+
+    #0. load data
+    source = JArray(JDouble, 1)(sourceArray)
+    destination = JArray(JDouble, 1)(destArray)
+
+    # 4. Supply the sample data:
+    calc.setObservations(source, destination)
+    # 5. Compute the estimate:
+    result = calc.computeAverageLocalOfObservations()
+    # 6. Compute the (statistical significance via) null distribution empirically:
+    measDist = calc.computeSignificance(numSurrogates)
+    rotated_pvalues.append(measDist.pValue)
+
+# Plot CDFs of shuffled and rotated p-values
+plt.figure(figsize=(12, 6))
+
+# Plot shuffled p-values CDF
+plt.subplot(1, 2, 1)
+plot_cdf(shuffled_pvalues)
+plt.axis([0,1,0,1])
+plt.gca().set_aspect('equal', adjustable='box')
+plt.title('CDF of Shuffled Surrogates P-values')
+plt.xlabel('P-value')
+plt.ylabel('CDF')
+
+# Plot rotated p-values CDF
+plt.subplot(1, 2, 2)
+plot_cdf(rotated_pvalues)
+plt.axis([0,1,0,1])
+plt.gca().set_aspect('equal', adjustable='box')
+plt.title('CDF of Rotated Surrogates P-values')
+plt.xlabel('P-value')
+plt.ylabel('CDF')
+
+plt.tight_layout()
+plt.show()
+
+

java/source/infodynamics/demos/autoanalysis/AutoAnalyserMI.java

@@ -101,15 +101,24 @@ protected void makeSpecificInitialisations() {
 		abstractContinuousClass = MutualInfoCalculatorMultiVariate.class;
 		// Common properties for all continuous calcs:
 		commonContPropertyNames = new String[] {
-				MutualInfoCalculatorMultiVariate.PROP_TIME_DIFF
+				MutualInfoCalculatorMultiVariate.PROP_TIME_DIFF,
+				MutualInfoCalculatorMultiVariate.PROP_SURROGATE_TYPE,
+				MutualInfoCalculatorMultiVariate.PROP_DYN_CORR_EXCL_TIME,
 		};
 		commonContPropertiesFieldNames = new String[] {
-				"PROP_TIME_DIFF"
+				"PROP_TIME_DIFF",
+				"PROP_SURROGATE_TYPE",
+				"PROP_DYN_CORR_EXCL_TIME"
 		};
 		commonContPropertyDescriptions = new String[] {
-				"Time-lag from source to dest to consider MI across; must be >= 0 (0 for standard MI)"
+				"Time-lag from source to dest to consider MI across; must be >= 0 (0 for standard MI)",
+				"Which strategy type to choose for selecting surrogates, default is " + MutualInfoCalculatorMultiVariate.PROP_SHUFFLE,
+				"Dynamic correlation exclusion time or <br/>Theiler window (see Kantz and Schreiber); " +
+					"0 (default) means no dynamic exclusion window. Only used for rotated surrogate selection for Gaussian Estimator",
 		};
 		commonContPropertyValueChoices = new String[][] {
+				null,
+				MutualInfoCalculatorMultiVariateKraskov.VALID_SURROGATE_TYPES,
 				null
 		};
 		// Gaussian properties:
@@ -130,21 +139,17 @@ protected void makeSpecificInitialisations() {
 		// Kernel:
 		kernelProperties = new String[] {
 				MutualInfoCalculatorMultiVariateKernel.KERNEL_WIDTH_PROP_NAME,
-				MutualInfoCalculatorMultiVariateKernel.DYN_CORR_EXCL_TIME_NAME,
 				MutualInfoCalculatorMultiVariateKernel.NORMALISE_PROP_NAME,			
 		};
 		kernelPropertiesFieldNames = new String[] {
 				"KERNEL_WIDTH_PROP_NAME",
-				"DYN_CORR_EXCL_TIME_NAME",
 				"NORMALISE_PROP_NAME"			
 		};
 		kernelPropertyDescriptions = new String[] {
 				"Kernel width to be used in the calculation. <br/>If the property " +
 						MutualInfoCalculatorMultiVariateKernel.NORMALISE_PROP_NAME +
 						" is set, then this is a number of standard deviations; " +
 						"otherwise it is an absolute value.",
-				"Dynamic correlation exclusion time or <br/>Theiler window (see Kantz and Schreiber); " +
-						"0 (default) means no dynamic exclusion window",
 				"(boolean) whether to normalise <br/>each incoming time-series to mean 0, standard deviation 1, or not  (recommended)",
 		};
 		kernelPropertyValueChoices = new String[][] {
@@ -157,7 +162,6 @@ protected void makeSpecificInitialisations() {
 				MutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE,
 				MutualInfoCalculatorMultiVariateKraskov.PROP_K,
 				MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE,
-				MutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME,
 				MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE,
 				MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS,
 				MutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU,
@@ -166,7 +170,6 @@ protected void makeSpecificInitialisations() {
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE",
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_K",
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE",
-				"MutualInfoCalculatorMultiVariateKraskov.PROP_DYN_CORR_EXCL_TIME",
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_NORM_TYPE",
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_NUM_THREADS",
 				"MutualInfoCalculatorMultiVariateKraskov.PROP_USE_GPU"
@@ -177,8 +180,6 @@ protected void makeSpecificInitialisations() {
 				"Standard deviation for an amount <br/>of random Gaussian noise to add to each variable, " +
 						"to avoid having neighbourhoods with artificially large counts. <br/>" +
 						"(\"false\" may be used to indicate \"0\".). The amount is added in after any normalisation.",
-				"Dynamic correlation exclusion time or <br/>Theiler window (see Kantz and Schreiber); " +
-						"0 (default) means no dynamic exclusion window",
 				"<br/>Norm type to use in KSG algorithm between the points in each marginal space. <br/>Options are: " +
 						"\"MAX_NORM\" (default), otherwise \"EUCLIDEAN\" or \"EUCLIDEAN_SQUARED\" (both equivalent here)",
 				"Number of parallel threads to use <br/>in computation: an integer > 0 or \"USE_ALL\" " +
@@ -189,7 +190,6 @@ protected void makeSpecificInitialisations() {
 				{"true", "false"},
 				null,
 				null,
-				null,
 				{"MAX_NORM", "EUCLIDEAN", "EUCLIDEAN_SQUARED"},
 				null,
 				{"true", "false"},

java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariate.java

@@ -18,6 +18,10 @@
 
 package infodynamics.measures.continuous;
 
+import java.util.Arrays;
+import java.util.HashSet;
+import java.util.Set;
+
 /**
  * <p>Interface for implementations of the <b>mutual information</b>,
  * which may be applied to either multivariate or merely univariate
@@ -84,6 +88,28 @@ public interface MutualInfoCalculatorMultiVariate
 	 * mean 0, standard deviation 1 (default true)
 	 */
 	public static final String PROP_NORMALISE = "NORMALISE";
+	/**
+	 * Property name for surrogate type strategy
+	 */
+	public static final String PROP_SURROGATE_TYPE = "SURROGATE_TYPE";
+	/**
+	 * Property name for a dynamics exclusion time window 
+	 * otherwise known as Theiler window (see Kantz and Schreiber).
+	 * Default is 0 which means no dynamic exclusion window.
+	 */
+	public static final String PROP_DYN_CORR_EXCL_TIME = "DYN_CORR_EXCL"; 
+	/**
+	 * Name of shuffle value surrogate type
+	 */
+	public static final String PROP_SHUFFLE = "SHUFFLE";
+	 /**
+	 * Name of rotate value surrogate type
+	 */
+	public static final String PROP_ROTATE = "ROTATE";
+	/**
+	 * Valid options for surrogate type
+	 */
+	public static final String[] VALID_SURROGATE_TYPES = {PROP_SHUFFLE, PROP_ROTATE};
 	/**
 	 * Property name for the std deviation of random Gaussian noise to be
 	 *  added to the data (default is 0, except for Kraskov/KSG estimator