diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java index e774fd8..62dbbd7 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java @@ -466,7 +466,7 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov returnValues = new double[numTimePointsToComputeFor]; } else { // We're computing average MI - returnValues = new double[MiKraskovThreadRunner.RETURN_ARRAY_LENGTH]; + returnValues = new double[MiKraskovThreadRunner.RETURN_ARRAY_LENGTH_MI]; } // Distribute the observations to the threads for the parallel processing @@ -504,7 +504,7 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov returnValues, runners[t].myStartTimePoint, runners[t].numberOfTimePoints); } else { // We're computing the average MI, keep the running sums of digammas and counts - MatrixUtils.addInPlace(returnValues, runners[t].getReturnValues()); + MatrixUtils.addInPlace(returnValues, runners[t].getReturnValues(), MiKraskovThreadRunner.RETURN_ARRAY_LENGTH_MI); } } } @@ -533,6 +533,57 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov } } + /** + * Use the same nearest neighbour searches to return a KL-based + * conditional entropy, sharing the neighbour search radius. + * With reference to the KSG paper, this takes the KL-based estimator for H(X) (eqn 22) and adds it to the negative of MI. + * It should also be noted that if the variables are normalised, this changes the absolute values of the entropies, + * and they become relative to normalised variables. + * Since this is a side-method for this estimator, it is currently is only implemented single threaded on CPU. + * It should also be considered experimental - this will eventually be moved to its own estimator. + * + * @return the average conditional entropy of variable1 given variable2 in nats (not bits!) + */ + public double computeAverageConditionalEntropy() throws Exception { + // Compute the MI + double startClockTime = Calendar.getInstance().getTimeInMillis(); + double conditionalEntropy = 0.0; + conditionalEntropy = computeFromObservations(false, null)[0]; + + int N = sourceObservations.length; // number of observations for the PDFs + + // Single-threaded implementation: + ensureKdTreesConstructed(); + double avDigammaCond = 0; + double avLog2xkNNDist = 0; + for (int t = 0; t < N; t++) { + // For this sample only: + double[] returnValues = partialComputeFromObservations(t, 1, false); + int n_y = (int) returnValues[MiKraskovThreadRunner.INDEX_SUM_NY]; + if (isAlgorithm1) { + avDigammaCond += MathsUtils.digamma(n_y+1); + } else { + avDigammaCond += MathsUtils.digamma(n_y); + } + double doublekNNDist = returnValues[MiKraskovThreadRunner.INDEX_SUM_2X_NEIGH_DIST]; + avLog2xkNNDist += Math.log(doublekNNDist); + } + avDigammaCond /= (double) N; + avLog2xkNNDist *= (double) dimensionsSource / (double) N; + if (isAlgorithm1) { + conditionalEntropy = -digammaK + avDigammaCond + avLog2xkNNDist; + } else { + conditionalEntropy = -digammaK + (1.0 / (double)k) + avDigammaCond + avLog2xkNNDist; + } + + if (debug) { + Calendar rightNow2 = Calendar.getInstance(); + long endTime = rightNow2.getTimeInMillis(); + System.out.println("Calculation time: " + ((endTime - startClockTime)/1000.0) + " sec" ); + } + return conditionalEntropy; + } + /** * Protected method to be used internally for threaded implementations. * This method implements the guts of each Kraskov algorithm, computing the number of @@ -860,7 +911,9 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov public static final int INDEX_SUM_DIGAMMAS = 0; public static final int INDEX_SUM_NX = 1; public static final int INDEX_SUM_NY = 2; - public static final int RETURN_ARRAY_LENGTH = 3; + public static final int INDEX_SUM_2X_NEIGH_DIST = 3; // Not used by GPU, and only read for conditional entropy + public static final int RETURN_ARRAY_LENGTH_MI = 3; + public static final int RETURN_ARRAY_LENGTH = 4; public MiKraskovThreadRunner( MutualInfoCalculatorMultiVariateKraskov miCalc, diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java index 6d44a1c..c68040e 100644 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov1.java @@ -71,6 +71,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 double sumDiGammas = 0; double sumNx = 0; double sumNy = 0; + double sum2xkNNDist = 0; for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) { // Compute eps for this time step by @@ -89,6 +90,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 int n_y = nnSearcherDest.countPointsStrictlyWithinR( t, kthNnData.distance, dynCorrExclTime); + sum2xkNNDist += kthNnData.distance + kthNnData.distance; sumNx += n_x; sumNy += n_y; // And take the digammas: @@ -116,7 +118,7 @@ public class MutualInfoCalculatorMultiVariateKraskov1 if (returnLocals) { return localMi; } else { - return new double[] {sumDiGammas, sumNx, sumNy}; + return new double[] {sumDiGammas, sumNx, sumNy, sum2xkNNDist}; } } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java index 01ed869..5e565c3 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov2.java @@ -73,6 +73,7 @@ public class MutualInfoCalculatorMultiVariateKraskov2 double sumDiGammas = 0; double sumNx = 0; double sumNy = 0; + double sum2xkNNDist = 0; for (int t = startTimePoint; t < startTimePoint + numTimePoints; t++) { // Compute eps_x and eps_y for this time step by @@ -94,6 +95,7 @@ public class MutualInfoCalculatorMultiVariateKraskov2 } } + sum2xkNNDist += eps_x + eps_x; // Count the number of points whose x distance is less // than or equal to eps_x, and whose y distance is less // than or equal to eps_y @@ -126,7 +128,7 @@ public class MutualInfoCalculatorMultiVariateKraskov2 if (returnLocals) { return localMi; } else { - return new double[] {sumDiGammas, sumNx, sumNy}; + return new double[] {sumDiGammas, sumNx, sumNy, sum2xkNNDist}; } } diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index d02f263..c04f1de 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -615,7 +615,8 @@ public class MatrixUtils { } /** - * Adds two arrays together, returning the result in input1 + * Adds two arrays together, returning the result in input1. + * Requires the two arrays to be of equal length. * * @param input1 * @param input2 @@ -624,7 +625,18 @@ public class MatrixUtils { if (input1.length != input2.length) { throw new Exception("Lengths of arrays are not equal"); } - for (int i = 0; i < input1.length; i++) { + addInPlace(input1, input2, input1.length); + } + + /** + * Adds two arrays together, returning the result in input1. + * Adds only the first length items. + * + * @param input1 + * @param input2 + */ + public static void addInPlace(double[] input1, double[] input2, int length) throws Exception { + for (int i = 0; i < length; i++) { input1[i] = input1[i] + input2[i]; } } diff --git a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java index e173a1f..bed843f 100644 --- a/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java +++ b/java/unittests/infodynamics/measures/continuous/kraskov/MutualInfoMultiVariateTester.java @@ -560,5 +560,97 @@ public class MutualInfoMultiVariateTester } } } + + /** + * Test the experimental conditional entropy method + * on random Gaussians + * + * @throws Exception if file not found + * + */ + public void testConditionalEntropyforNoisyIndependentVariablesFromFile() throws Exception { + + // We'll just take the first two columns from this data set + // Works well on 10ColsRandomGaussian-1.txt because it is Gaussian + ArrayFileReader afr = new ArrayFileReader("demos/data/10ColsRandomGaussian-1.txt"); + double[][] data = afr.getDouble2DMatrix(); + + // When normalising the marginals to std dev 1, we know that entropy of the marginal + // is expected to be: + double expected_Hx = 0.5 * Math.log(2.0 * Math.PI * Math.E); + + System.out.println("Kraskov comparison - conditional entropy"); + + MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(1); + + miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests + + int[] numSamplesToCheck = new int[] {1000, 3000, 5000, data.length}; + double conditionalEnt = 0, expected_H_X_given_Y = 0; + for (int ni = 0; ni < numSamplesToCheck.length; ni++) { + miCalc.initialise(1, 1); + miCalc.setObservations( + MatrixUtils.selectRowsAndColumns(data, 0, numSamplesToCheck[ni], 0, 1), + MatrixUtils.selectRowsAndColumns(data, 0, numSamplesToCheck[ni], 1, 1)); + miCalc.setDebug(true); + double mi = miCalc.computeAverageLocalOfObservations(); + miCalc.setDebug(false); + conditionalEnt = miCalc.computeAverageConditionalEntropy(); + expected_H_X_given_Y = -mi + expected_Hx; + + System.out.printf("k=%s: Average MI %.8f; Average H(X|Y) = %.8f (expected %.8f) from %d samples\n", + miCalc.getProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_K), mi, + conditionalEnt, expected_H_X_given_Y, miCalc.getNumObservations()); + } + // Only check the assertion on the full data set: + assertEquals(expected_H_X_given_Y, conditionalEnt, 0.01); + + } + + /** + * Test the experimental conditional entropy method + * on random coupled uniform variables + * + * @throws Exception if file not found + * + */ + public void testConditionalEntropyforNoisyDependentVariablesFromFile() throws Exception { + + // We'll just take the first two columns from this data set + // We'll use 4ColsPairedNoisyDependence where first column is randomly distributed + // on 0..1 and third adds some noise to that. + ArrayFileReader afr = new ArrayFileReader("demos/data/4ColsPairedNoisyDependence-1.txt"); + double[][] data = afr.getDouble2DMatrix(); + + // When the marginal x is uniformly distributed on 0..1, we know that entropy of the marginal + // is expected to be 0 -- when we don't normalise the variables! + double expected_Hx = 0; + + System.out.println("Kraskov comparison - conditional entropy"); + + MutualInfoCalculatorMultiVariateKraskov miCalc = getNewCalc(1); + + miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_ADD_NOISE, "0"); // Need consistency for unit tests + miCalc.setProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_NORMALISE, "false"); // Need to avoid normalising for the expected value to hold + + int[] numSamplesToCheck = new int[] {1000, 2000, data.length}; + double conditionalEnt = 0, expected_H_X_given_Y = 0; + for (int ni = 0; ni < numSamplesToCheck.length; ni++) { + miCalc.initialise(1, 1); + miCalc.setObservations( + MatrixUtils.selectRowsAndColumns(data, 0, numSamplesToCheck[ni], 0, 1), + MatrixUtils.selectRowsAndColumns(data, 0, numSamplesToCheck[ni], 2, 1)); + double mi = miCalc.computeAverageLocalOfObservations(); + conditionalEnt = miCalc.computeAverageConditionalEntropy(); + expected_H_X_given_Y = -mi + expected_Hx; + + System.out.printf("k=%s: Average MI %.8f; Average H(X|Y) = %.8f (expected %.8f) from %d samples\n", + miCalc.getProperty(MutualInfoCalculatorMultiVariateKraskov.PROP_K), mi, + conditionalEnt, expected_H_X_given_Y, miCalc.getNumObservations()); + } + // Only check the assertion on the full data set: + assertEquals(expected_H_X_given_Y, conditionalEnt, 0.02); + + } }