From ab21895321ebd410e82dd702bf4b5bff30d99b1d Mon Sep 17 00:00:00 2001
From: "joseph.lizier" <joseph.lizier@gmail.com>
Date: Tue, 5 Aug 2014 02:04:56 +0000
Subject: [PATCH] Correcting TE Kernel estimator's bias correction method to
 match Schreiber's approach to including contributions from other correlation
 sums when the full joint space has no neighbours.

---
 .../kernel/KernelEstimatorUniVariate.java     |  4 ++
 .../TransferEntropyCalculatorKernel.java      | 41 ++++++++++++++++++-
 2 files changed, 43 insertions(+), 2 deletions(-)

diff --git a/java/source/infodynamics/measures/continuous/kernel/KernelEstimatorUniVariate.java b/java/source/infodynamics/measures/continuous/kernel/KernelEstimatorUniVariate.java
index dc869c9..ba490c8 100755
--- a/java/source/infodynamics/measures/continuous/kernel/KernelEstimatorUniVariate.java
+++ b/java/source/infodynamics/measures/continuous/kernel/KernelEstimatorUniVariate.java
@@ -258,6 +258,10 @@ public class KernelEstimatorUniVariate {
 			System.out.println("Count after upper bin " + (bin + 1) + " = " + count);
 		}
 
+		// TODO Should we divide this result by the kernel width (in use)?
+		//  Schreiber doesn't do this in the TE paper, but it is done in 
+		//  the Kaiser and Schreiber paper. It won't matter to TE, but
+		//  will make a difference to individual entropy estimates.
 		return (double) count / (double) totalTimePointsCompared;
 	}
 	
diff --git a/java/source/infodynamics/measures/continuous/kernel/TransferEntropyCalculatorKernel.java b/java/source/infodynamics/measures/continuous/kernel/TransferEntropyCalculatorKernel.java
index 0c7b5c1..1040105 100755
--- a/java/source/infodynamics/measures/continuous/kernel/TransferEntropyCalculatorKernel.java
+++ b/java/source/infodynamics/measures/continuous/kernel/TransferEntropyCalculatorKernel.java
@@ -237,6 +237,8 @@ public class TransferEntropyCalculatorKernel
 	 *   though the adverse correction of the bias correction is worse than the correction
 	 *   itself (because the probabilities being multiplied/divided are not independent), 
 	 *   so recommend not to use this method.
+	 *   (Bias correction is implemented properly in the Kraskov et al. 
+	 *   estimators, see {@link infodynamics.measures.continuous.kraskov.TransferEntropyCalculatorKraskov}
 	 * </p>
 	 * <p>It is implemented here for testing purposes only.</p>
 	 * 
@@ -244,25 +246,60 @@ public class TransferEntropyCalculatorKernel
 	public double computeAverageLocalOfObservationsWithCorrection() throws Exception {
 		
 		double te = 0.0;
+		int numNoNeighbours = 0;
+		double contributionsNoNeighbours = 0.0;
 		for (int b = 0; b < totalObservations; b++) {
 			TransferEntropyKernelCounts kernelCounts = teKernelEstimator.getCount(destPastVectors[b],
 					destNextValues[b], sourceValues[b], b);
 			double cont = 0.0;
-			if (kernelCounts.countNextPastSource > 0) {
+			// Original code:
+			/* if (kernelCounts.countNextPastSource > 0) {
 				cont = MathsUtils.digamma(kernelCounts.countNextPastSource) - 
 						MathsUtils.digamma(kernelCounts.countPastSource) -
 						MathsUtils.digamma(kernelCounts.countNextPast) +
 						MathsUtils.digamma(kernelCounts.countPast);
+			} */
+			// But Schreiber confirmed to me that with dynamic correlation 
+			//  exclusion, when you may have no nearest neighbours,
+			//  he was allowing contributions from other groups (e.g. countPastSource)
+			//  to be added even if the full joint count was zero.
+			// Implement it like this:
+			//  TODO Do we need to correct the totalObservations
+			//   divisor for each digamma sum to account for this?
+			//   (Schreiber does not do that; for the moment we'll accept that
+			//    as the right approach)
+			if (kernelCounts.countPastSource > 0) {
+				cont -= MathsUtils.digamma(kernelCounts.countPastSource);
+			}
+			if (kernelCounts.countNextPast > 0) {
+				cont -= MathsUtils.digamma(kernelCounts.countNextPast);
+			}
+			if (kernelCounts.countPast > 0) {
+				cont += MathsUtils.digamma(kernelCounts.countPast);
+			}
+			if (kernelCounts.countNextPastSource > 0) {
+				cont += MathsUtils.digamma(kernelCounts.countNextPastSource);
+			} else {
+				// These contributions are from a set with no neighbours in 
+				//  the full joint space.
+				contributionsNoNeighbours += cont;
+				numNoNeighbours++;
 			}
 			te += cont;
 			/*
 			if (debug) {
-				System.out.println(b + ": " + logTerm + " -> " + (cont/Math.log(2.0)) + " -> sum: " + (te/Math.log(2.0)));
+				System.out.println(b + ": " + cont + " -> " + (cont/Math.log(2.0)) + " -> sum: " + (te/Math.log(2.0)));
 			}
 			*/
 		}
 		// Average it, and convert results to bytes
 		lastAverage = te / (double) totalObservations / Math.log(2.0);
+		if (debug) {
+			System.out.printf("TE=%.4f, with %d contributions from 0 neighbour sets being %.4f\n",
+					lastAverage,
+					numNoNeighbours,
+					contributionsNoNeighbours / (double) totalObservations / Math.log(2.0));
+		}
 		return lastAverage;
 	}