mirror of https://github.com/jlizier/jidt
Amending Kraskov (KSG) MI estimators to have an experimental method to provide conditional entropy of the first variable given the second. Works by removing the Kozachenko Leonenko entropy of variable 1 from the MI, using the same kNN radii as the MI estimator. Includes Unit tests to provide some initial validation.
This commit is contained in:
jlizier
parent
25159009ca
commit
2714650cfa
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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};
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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};
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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];
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue