From cda8689afb99f9a8fa20cb6d0d4abbdb53d19705 Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Wed, 5 Sep 2012 06:23:48 +0000 Subject: [PATCH] Made sure all discrete calculators use log_2 to make answers in bits (some were still using log of base). Merged usual and debug methods for computing average MI in discrete calculator Added setDebug to super InfoMeasureCalculator Cleaned up header of InfoMeasureCalculator Adding TODO comments for Kraskov calculators Bug fix on MutualInfoCalculatorMultiVariateWithDiscreteKraskov - we weren't normalising incoming observations for computing local MI properly, this is fixed now. --- ...nfoCalculatorMultiVariateWithDiscrete.java | 2 +- ...tualInfoCalculatorMultiVariateKraskov.java | 4 +- ...ulatorMultiVariateWithDiscreteKraskov.java | 7 ++ .../kraskov/MultiInfoCalculatorKraskov.java | 3 + ...tualInfoCalculatorMultiVariateKraskov.java | 4 +- ...oCalculatorMultiVariateKraskovByMulti.java | 2 + ...ulatorMultiVariateWithDiscreteKraskov.java | 99 +++++++++++++++++-- .../ApparentTransferEntropyCalculator.java | 2 +- .../discrete/InfoMeasureCalculator.java | 47 +++++++-- .../discrete/MultiInformationCalculator.java | 2 +- .../discrete/MutualInformationCalculator.java | 52 ++-------- .../discrete/SeparableInfoCalculator.java | 2 +- .../infodynamics/utils/MatrixUtils.java | 11 ++- 13 files changed, 171 insertions(+), 66 deletions(-) diff --git a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java index 8d0eabb..a075ba3 100755 --- a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java +++ b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java @@ -110,7 +110,7 @@ public interface MutualInfoCalculatorMultiVariateWithDiscrete { * * @param numPermutationsToCheck the number of permuted surrogates to examine * @return the proportion of MI scores from the distribution which have higher or equal MIs to ours. - * @link "Chavez et. al., 'Statistical assessment of nonlinear causality: + * @see "Chavez et. al., 'Statistical assessment of nonlinear causality: * application to epileptic EEG signals', Journal of Neuroscience Methods 124 (2003) 113-128" * @throws Exception */ diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java index 95b7131..c7959d5 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java @@ -228,7 +228,9 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov { public abstract double[] computeLocalOfPreviousObservations() throws Exception; public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { - // If implemented, will need to incorporate any time difference here. + // If implemented, will need to incorporate any normalisation here + // (normalising the incoming data the same way the previously + // supplied observations were normalised). throw new Exception("Local method not implemented yet"); } diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java index e26c32f..55049ff 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java @@ -541,6 +541,10 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl public double[] computeLocalUsingPreviousObservations(double[][] continuousStates, int[] discreteStates) throws Exception { + // TODO Implement local method. + // Note: will need to keep the means and stds of supplied observations + // if we normalised them (since we'll need to normalise the + // observations supplied here to match them) throw new Exception("Local method not implemented yet"); } @@ -560,6 +564,9 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl double[][] contStates, int[] discreteStates, double[][] conditionedStates) throws Exception { // TODO Auto-generated method stub + // Note: will need to keep the means and stds of supplied observations + // if we normalised them (since we'll need to normalise the + // observations supplied here to match them) throw new Exception("Not implemented yet"); } } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov.java index d1ec679..0c650e8 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MultiInfoCalculatorKraskov.java @@ -244,6 +244,9 @@ public abstract class MultiInfoCalculatorKraskov implements public abstract double[] computeLocalOfPreviousObservations() throws Exception; public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { + // TODO If this is implemented, will need to normalise the incoming + // observations the same way that previously supplied ones were + // normalised (if they were normalised, that is) throw new Exception("Local method not implemented yet"); } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java index 5a92fa2..2fbd91b 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskov.java @@ -229,7 +229,9 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov implements public abstract double[] computeLocalOfPreviousObservations() throws Exception; public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { - // If implemented, will need to incorporate any time difference here. + // TODO If implemented, will need to incorporate any time difference here. + // Will also need to handle normalisation of the incoming data + // appropriately throw new Exception("Local method not implemented yet"); } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskovByMulti.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskovByMulti.java index f885c40..17bbed9 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskovByMulti.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateKraskovByMulti.java @@ -328,6 +328,8 @@ public abstract class MutualInfoCalculatorMultiVariateKraskovByMulti implements } public double[] computeLocalUsingPreviousObservations(double[][] states1, double[][] states2) throws Exception { + // TODO When implemented, we'll need to normalise the data in the same way + // as previously supplied observations (if they were normalised here) throw new Exception("Local method not implemented yet"); } diff --git a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateWithDiscreteKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateWithDiscreteKraskov.java index aa2f750..a118cb5 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateWithDiscreteKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/MutualInfoCalculatorMultiVariateWithDiscreteKraskov.java @@ -35,7 +35,14 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu protected double[][] continuousData; protected int[] discreteData; protected int[] counts; + /** + * Number of possible states of the discrete variable + */ protected int base; + /** + * Number of dimenions of the joint continuous variable + */ + protected int dimensions; protected boolean debug; protected double mi; protected boolean miComputed; @@ -50,7 +57,19 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu public final static String PROP_K = "k"; public final static String PROP_NORM_TYPE = "NORM_TYPE"; public static final String PROP_NORMALISE = "NORMALISE"; - private boolean normalise = true; + + /** + * Track whether we're going to normalise the joint variables individually + */ + protected boolean normalise = true; + /** + * Track the means of the joint variables if we are normalising them + */ + protected double[] means; + /** + * Track the std devs of the joint variables if we are normalising them + */ + protected double[] stds; public MutualInfoCalculatorMultiVariateWithDiscreteKraskov() { super(); @@ -68,8 +87,10 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu miComputed = false; xNorms = null; continuousData = null; + means = null; + stds = null; discreteData = null; - // No need to keep the dimenions here + this.dimensions = dimensions; this.base = base; } @@ -122,11 +143,18 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu if (continuousObservations[0].length == 0) { throw new Exception("Computing MI with a null set of data"); } + if (continuousObservations[0].length != dimensions) { + throw new Exception("The continuous observations do not have the expected number of variables (" + dimensions + ")"); + } continuousData = continuousObservations; discreteData = discreteObservations; if (normalise) { // Take a copy since we're going to normalise it - continuousData = MatrixUtils.normaliseIntoNewArray(continuousObservations); + // And we need to keep the means/stds ready to normalise local values + // that are supplied later: + means = MatrixUtils.means(continuousObservations); + stds = MatrixUtils.stdDevs(continuousObservations, means); + continuousData = MatrixUtils.normaliseIntoNewArray(continuousObservations, means, stds); } // count the discrete states: counts = new int[base]; @@ -237,6 +265,9 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu return computeAverageLocalOfObservationsWhileComputingDistances(); } + // Postcondition: we'll compute the norms before the main loop, + // unless they have already been computed: + if (xNorms == null) { computeNorms(); } @@ -247,6 +278,8 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu double avNx = 0; double avNy = 0; + double testSum = 0.0; // Used for debugging prints + for (int t = 0; t < N; t++) { // Compute eps_x and eps_y for this time step: // using x norms to all neighbours @@ -269,7 +302,20 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu avNy += n_y; // And take the digamma before adding into the // average: - averageDiGammas += MathsUtils.digamma(n_x) + MathsUtils.digamma(n_y); + double localSum = MathsUtils.digamma(n_x) + MathsUtils.digamma(n_y); + averageDiGammas += localSum; + + if (debug) { + double localValue = MathsUtils.digamma(k) - 1.0/(double)k - localSum + MathsUtils.digamma(N); + testSum += localValue; + if (dimensions == 1) { + System.out.printf("t=%d: x=%.3f, eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, continuousData[t][0], eps_x, n_x, n_y, localValue, testSum); + } else { + System.out.printf("t=%d: eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, eps_x, n_x, n_y, localValue, testSum); + } + } } averageDiGammas /= (double) N; if (debug) { @@ -304,6 +350,8 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu double avNx = 0; double avNy = 0; + double testSum = 0.0; // Used for debugging prints + for (int t = 0; t < N; t++) { // Compute eps_x and eps_y for this time step: // First get x norms to all neighbours @@ -335,7 +383,20 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu avNy += n_y; // And take the digamma before adding into the // average: - averageDiGammas += MathsUtils.digamma(n_x) + MathsUtils.digamma(n_y); + double localSum = MathsUtils.digamma(n_x) + MathsUtils.digamma(n_y); + averageDiGammas += localSum; + + if (debug) { + double localValue = MathsUtils.digamma(k) - 1.0/(double)k - localSum + MathsUtils.digamma(N); + testSum += localValue; + if (dimensions == 1) { + System.out.printf("t=%d: x=%.3f, eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, continuousData[t][0], eps_x, n_x, n_y, localValue, testSum); + } else { + System.out.printf("t=%d: eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, eps_x, n_x, n_y, localValue, testSum); + } + } } averageDiGammas /= (double) N; if (debug) { @@ -425,7 +486,7 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu * it should be ignored, since the data point was not counted in its * own counts in the standard version anyway.

* - *

The supplied observations are intended to be new observations, + *

Importantly, the supplied observations are intended to be new observations, * not those fed in to compute the PDFs from. There would be * some subtle changes necessary to accomodate computing locals on * the same data set (e.g. not counting the current point as one @@ -440,11 +501,25 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu public double[] computeLocalUsingPreviousObservations(double[][] continuousNewStates, int[] discreteNewStates) throws Exception { + if (normalise) { + // The stored observations continuousData have been normalised + // according to their stored means and stds; we need to + // normalise the incoming observations the same way before + // comparing them + continuousNewStates = MatrixUtils.normaliseIntoNewArray( + continuousNewStates, means, stds); + } + int N = continuousNewStates.length; // number of observations double[] locals = new double[N]; double fixedPartOfLocals = MathsUtils.digamma(k) - 1.0/(double)k + MathsUtils.digamma(N); + double testSum = 0.0; + if (debug) { + System.out.printf("digamma(k)=%.3f - 1/k=%.3f + digamma(N)=%.3f\n", + MathsUtils.digamma(k), 1.0/(double)k, MathsUtils.digamma(N)); + } double avNx = 0; double avNy = 0; @@ -476,11 +551,21 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKraskov implements Mutu // Now compute the local value: locals[t] = fixedPartOfLocals - MathsUtils.digamma(n_x) - MathsUtils.digamma(n_y); + if (debug) { + testSum += locals[t]; + if (dimensions == 1) { + System.out.printf("t=%d: x=%.3f, eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, continuousNewStates[t][0], eps_x, n_x, n_y, locals[t], testSum); + } else { + System.out.printf("t=%d: eps_x=%.3f, n_x=%d, n_y=%d, local=%.3f, running total = %.5f\n", + t, eps_x, n_x, n_y, locals[t], testSum); + } + } } if (debug) { avNx /= (double)N; avNy /= (double)N; - System.out.println(String.format("Average n_x=%.3f, Average n_y=%.3f", avNx, avNy)); + System.out.printf("Average n_x=%.3f, Average n_y=%.3f\n", avNx, avNy); } return locals; diff --git a/java/source/infodynamics/measures/discrete/ApparentTransferEntropyCalculator.java b/java/source/infodynamics/measures/discrete/ApparentTransferEntropyCalculator.java index 177a05a..4b5dd1c 100755 --- a/java/source/infodynamics/measures/discrete/ApparentTransferEntropyCalculator.java +++ b/java/source/infodynamics/measures/discrete/ApparentTransferEntropyCalculator.java @@ -638,7 +638,7 @@ public class ApparentTransferEntropyCalculator extends ContextOfPastMeasureCalcu double logTerm = (double) nextPastCount[nextVal][prevVal] / (double) pastCount[prevVal] / p_next; - double localValue = Math.log(logTerm) / log_base; + double localValue = Math.log(logTerm) / log_2; activeCont = (nextPastCount[nextVal][prevVal] / (double) observations) * localValue; } else { diff --git a/java/source/infodynamics/measures/discrete/InfoMeasureCalculator.java b/java/source/infodynamics/measures/discrete/InfoMeasureCalculator.java index 17d5b66..2a827d3 100755 --- a/java/source/infodynamics/measures/discrete/InfoMeasureCalculator.java +++ b/java/source/infodynamics/measures/discrete/InfoMeasureCalculator.java @@ -1,15 +1,31 @@ package infodynamics.measures.discrete; /** - * Info theoretic measure calculator base class + *

Info theoretic measure calculator base class, providing common functionality + * for user-level measure classes.

* - * Usage: - * 1. Continuous accumulation of observations before computing : - * Call: a. initialise() - * b. addObservations() several times over - * c. computeLocalFromPreviousObservations() or computeAverageLocalOfObservations() - * 2. Standalone computation from a single set of observations: - * Call: computeLocal() or computeAverageLocal() + *

Usage of child classes is intended to follow this general pattern: + *

    + *
  1. Construct;
    2. + *
  2. {@link #initialise()};
    3. + *
  3. Then, either of the following: + *
      + *
    1. Continuous accumulation of observations before computing, via: + *
        + *
      1. calling "addObservations()" methods of children + * several times over;
        2. + *
      2. Compute required quantities, using + * {@link #computeAverageLocalOfObservations()} or + * "computeLocalUsinPreviousObservations()".
        3. + *
      2. + *
    2. Standalone computation from a single set of observations; call: + * "computeLocal()" or "computeAverageLocal()".
      3. + *
    + *

+ * + *

Note: various functions referred to above (e.g. "addObservations()") + * are not specified here, so that this class can be a superclass + * for both univariate, pairwise and multivariate methods.

* * @author Joseph Lizier * joseph.lizier at gmail.com @@ -120,4 +136,19 @@ public abstract class InfoMeasureCalculator { } return true; } + + /** + * Compute the average value of the measure from the previously supplied + * observations + * + * @return average value + */ + public abstract double computeAverageLocalOfObservations(); + + /** + * @param debug the debug status to set + */ + public void setDebug(boolean debug) { + this.debug = debug; + } } diff --git a/java/source/infodynamics/measures/discrete/MultiInformationCalculator.java b/java/source/infodynamics/measures/discrete/MultiInformationCalculator.java index 6b19397..d77ef1a 100755 --- a/java/source/infodynamics/measures/discrete/MultiInformationCalculator.java +++ b/java/source/infodynamics/measures/discrete/MultiInformationCalculator.java @@ -151,7 +151,7 @@ public class MultiInformationCalculator extends InfoMeasureCalculator { } double jointProb = (double) jointCount[jointValue] / (double) observations; double logValue = jointProb / prodMarginalProbs; - double localValue = Math.log(logValue) / log_base; + double localValue = Math.log(logValue) / log_2; if (jointProb > 0.0) { if (!checkedFirst) { max = localValue; diff --git a/java/source/infodynamics/measures/discrete/MutualInformationCalculator.java b/java/source/infodynamics/measures/discrete/MutualInformationCalculator.java index 3946e48..830b844 100755 --- a/java/source/infodynamics/measures/discrete/MutualInformationCalculator.java +++ b/java/source/infodynamics/measures/discrete/MutualInformationCalculator.java @@ -103,47 +103,9 @@ public class MutualInformationCalculator extends InfoMeasureCalculator max = 0; min = 0; double meanSqLocals = 0; - for (int i = 0; i < base; i++) { - // compute p_i - double probi = (double) iCount[i] / (double) observations; - for (int j = 0; j < base; j++) { - // compute p_j - double probj = (double) jCount[j] / (double) observations; - // compute p(veci=i, vecj=j) - double jointProb = (double) jointCount[i][j] / (double) observations; - // Compute MI contribution: - if (jointProb * probi * probj > 0.0) { - double localValue = Math.log(jointProb / (probi * probj)) / log_base; - miCont = jointProb * localValue; - if (localValue > max) { - max = localValue; - } else if (localValue < min) { - min = localValue; - } - // Add this contribution to the mean - // of the squared local values - meanSqLocals += miCont * localValue; - } else { - miCont = 0.0; - } - mi += miCont; - } + if (debug) { + System.out.println("i\tj\tp_i\tp_j\tp_joint\tlocal"); } - - average = mi; - std = Math.sqrt(meanSqLocals - average * average); - - return mi; - } - - public double debugLocalOfObservations() { - double mi = 0.0; - double miCont = 0.0; - - max = 0; - min = 0; - double meanSqLocals = 0; - System.out.println("i\tj\tp_i\tp_j\tp_joint\tlocal"); for (int i = 0; i < base; i++) { // compute p_i double probi = (double) iCount[i] / (double) observations; @@ -154,10 +116,12 @@ public class MutualInformationCalculator extends InfoMeasureCalculator double jointProb = (double) jointCount[i][j] / (double) observations; // Compute MI contribution: if (jointProb * probi * probj > 0.0) { - double localValue = Math.log(jointProb / (probi * probj)) / log_base; + double localValue = Math.log(jointProb / (probi * probj)) / log_2; miCont = jointProb * localValue; - System.out.printf("%d\t%d\t%.4f\t%.4f\t%.4f\t%.4f\n", - i, j, probi, probj, jointProb, localValue); + if (debug) { + System.out.printf("%d\t%d\t%.4f\t%.4f\t%.4f\t%.4f\n", + i, j, probi, probj, jointProb, localValue); + } if (localValue > max) { max = localValue; } else if (localValue < min) { @@ -273,7 +237,7 @@ public class MutualInformationCalculator extends InfoMeasureCalculator // and we've got two counts on the bottom // but one count on the top: logTerm *= (double) observations; - localMI[r] = Math.log(logTerm) / log_base; + localMI[r] = Math.log(logTerm) / log_2; average += localMI[r]; if (localMI[r] > max) { max = localMI[r]; diff --git a/java/source/infodynamics/measures/discrete/SeparableInfoCalculator.java b/java/source/infodynamics/measures/discrete/SeparableInfoCalculator.java index f1e0bbc..cbc5009 100755 --- a/java/source/infodynamics/measures/discrete/SeparableInfoCalculator.java +++ b/java/source/infodynamics/measures/discrete/SeparableInfoCalculator.java @@ -1618,7 +1618,7 @@ public class SeparableInfoCalculator extends ContextOfPastMeasureCalculator { if (computeMultiInfoCoherence) { miCalc.setDebug(debug); } - this.debug = debug; + super.setDebug(debug); } } diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index 6e8a02b..aa4419c 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -1838,9 +1838,18 @@ public class MatrixUtils { * @param matrix 2D matrix of doubles */ public static double[][] normaliseIntoNewArray(double[][] matrix) { - double[][] newMatrix = new double[matrix.length][matrix[0].length]; double[] means = means(matrix); double[] stds = stdDevs(matrix, means); + return normaliseIntoNewArray(matrix, means, stds); + } + + /** + * Normalises the elements along each column of the matrix + * + * @param matrix 2D matrix of doubles + */ + public static double[][] normaliseIntoNewArray(double[][] matrix, double[] means, double[] stds) { + double[][] newMatrix = new double[matrix.length][matrix[0].length]; for (int r = 0; r < newMatrix.length; r++) { for (int c = 0; c < newMatrix[r].length; c++) { newMatrix[r][c] = matrix[r][c] - means[c];