From 75c36dc13a3555cc9f9c109d8f09e1cc88fec22e Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Mon, 18 Mar 2013 06:08:12 +0000 Subject: [PATCH] Adding common parent implementation of Conditional MI calculator from continuous to discrete variable, plus a linear-Gaussian implementation of it --- ...rMultiVariateWithDiscreteSourceCommon.java | 237 ++++++++++++++++++ ...ultiVariateWithDiscreteSourceGaussian.java | 219 ++++++++++++++++ 2 files changed, 456 insertions(+) create mode 100755 java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon.java create mode 100755 java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian.java diff --git a/java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon.java b/java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon.java new file mode 100755 index 0000000..3d41deb --- /dev/null +++ b/java/source/infodynamics/measures/continuous/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon.java @@ -0,0 +1,237 @@ +package infodynamics.measures.continuous; + +import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.EmpiricalMeasurementDistribution; +import infodynamics.utils.RandomGenerator; + +/** + *

This is an abstract calculator to + * compute the Conditional Mutual Information I(X;D|Z) + * between a discrete variable D and a + * vector of continuous variables X, + * conditioned on another vector of continuous variables Z. + * It implements common methods for child classes to build on

+ * + * @author Joseph Lizier + */ +public abstract class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon + implements ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource { + + /** + * we compute distances to the kth neighbour + */ + protected double[][] continuousDataX; + protected double[][] conditionedDataZ; + protected int[] discreteData; + protected int[] counts; + protected int base; + protected int dimensionsContinuous; + protected int dimensionsConditional; + protected double[] meansX; + protected double[] stdsX; + protected double[] meansZ; + protected double[] stdsZ; + protected boolean debug; + protected double condMi; + protected boolean miComputed; + protected int totalObservations; + + public static final String PROP_NORMALISE = "NORMALISE"; + protected boolean normalise = true; + + public ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon() { + } + + /** + * Initialise the calculator. + * + * @param dimensions number of joint continuous variables + * @param base number of discrete states + * @param dimensionsCond the number of joint continuous variables + * to condition on + */ + public void initialise(int dimensions, int base, int dimensionsCond) { + condMi = 0.0; + miComputed = false; + continuousDataX = null; + discreteData = null; + this.base = base; + dimensionsContinuous = dimensions; + dimensionsConditional = dimensionsCond; + } + + /** + * Sets properties for the calculator. + * Valid properties include: + * + * + * @param propertyName + * @param propertyValue + */ + public void setProperty(String propertyName, String propertyValue) { + if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) { + normalise = Boolean.parseBoolean(propertyValue); + } + } + + public void startAddObservations() { + throw new RuntimeException("Not implemented yet"); + } + + public void finaliseAddObservations() throws Exception { + meansX = MatrixUtils.means(continuousDataX); + stdsX = MatrixUtils.stdDevs(continuousDataX, meansX); + meansZ = MatrixUtils.means(conditionedDataZ); + stdsZ = MatrixUtils.stdDevs(conditionedDataZ, meansZ); + if (normalise) { + // Take a copy since we're going to normalise it + continuousDataX = MatrixUtils.normaliseIntoNewArray(continuousDataX); + conditionedDataZ = MatrixUtils.normaliseIntoNewArray(conditionedDataZ); + } + // count the discrete states: + counts = new int[base]; + for (int t = 0; t < discreteData.length; t++) { + counts[discreteData[t]]++; + } + totalObservations = discreteData.length; + } + + public void addObservations(double[][] continuousObservations, + int[] discreteObservations, double[][] conditionedObservations) + throws Exception { + throw new RuntimeException("Not implemented yet"); + } + + public void setObservations(double[][] continuousObservations, + int[] discreteObservations, double[][] conditionedObservations) + throws Exception { + + if ((continuousObservations.length != discreteObservations.length) || + (continuousObservations.length != conditionedObservations.length)) { + throw new Exception("Time steps for observations2 " + + discreteObservations.length + " does not match the length " + + "of observations1 " + continuousObservations.length + + " and of conditionedObservations " + conditionedObservations.length); + } + if (continuousObservations[0].length == 0) { + throw new Exception("Computing MI with a null set of data"); + } + if (conditionedObservations[0].length == 0) { + throw new Exception("Computing MI with a null set of conditioned data"); + } + continuousDataX = continuousObservations; + discreteData = discreteObservations; + conditionedDataZ = conditionedObservations; + + finaliseAddObservations(); + + } + + /** + * Compute the significance of the mutual information of the previously supplied observations. + * We destroy the p(x,d|z) correlations, while retaining the p(x|z) marginals, to check how + * significant this conditional mutual information actually was. + * + * This is in the spirit of Chavez et. al., "Statistical assessment of nonlinear causality: + * application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128 + * which was performed for Transfer entropy. + * + * @param reorderDiscreteVariable boolean for whether to reorder the discrete variable (true) + * or the continuous variable (false) + * @param numPermutationsToCheck + * @return the proportion of conditional MI scores from the distribution which have higher or equal MIs to ours. + */ + public synchronized EmpiricalMeasurementDistribution computeSignificance( + boolean reorderDiscreteVariable, int numPermutationsToCheck) throws Exception { + // Generate the re-ordered indices: + RandomGenerator rg = new RandomGenerator(); + // Use continuousDataX length (all variables have same length) even though + // we may be randomising the other variable: + int[][] newOrderings = rg.generateDistinctRandomPerturbations( + continuousDataX.length, numPermutationsToCheck); + return computeSignificance(reorderDiscreteVariable, newOrderings); + } + + /** + *

As per {@link #computeSignificance(boolean, int)} but supplies + * the re-orderings of the observations of the named variable.

+ * + *

We provide a simple implementation which would be suitable for + * any child class, though the child class may prefer to make its + * own implementation to make class-specific optimisations. + * Child classes must implement {@link java.lang.Cloneable} + * for this method to be callable for them, and indeed implement + * the clone() method in a way that protects their structure + * from alteration by surrogate data being supplied to it.

+ * + * @param reorderDiscreteVariable boolean for whether to reorder the discrete variable (true) + * or the continuous variable (false) + * @param newOrderings the specific new orderings to use + * @return the proportion of conditional MI scores from the distribution which have higher or equal MIs to ours. + */ + public EmpiricalMeasurementDistribution computeSignificance( + boolean reorderDiscreteVariable, int[][] newOrderings) throws Exception { + + int numPermutationsToCheck = newOrderings.length; + if (!miComputed) { + computeAverageLocalOfObservations(); + } + + // Take a clone of the object to compute the MI of the surrogates: + // (this is a shallow copy, it doesn't make new copies of all + // the arrays - child classes should override this) + ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource miSurrogateCalculator = + (ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource) this.clone(); + + double[] surrogateMeasurements = new double[numPermutationsToCheck]; + + // Now compute the MI for each set of shuffled data: + for (int i = 0; i < numPermutationsToCheck; i++) { + if (reorderDiscreteVariable) { + // Reorder the discrete variable + int[] shuffledData = + MatrixUtils.extractSelectedTimePoints(discreteData, newOrderings[i]); + // Perform new initialisations + miSurrogateCalculator.initialise( + dimensionsContinuous, base, dimensionsConditional); + // Set new observations + miSurrogateCalculator.setObservations(continuousDataX, + shuffledData, conditionedDataZ); + } else { + // Reorder the continuous variable + double[][] shuffledData = + MatrixUtils.extractSelectedTimePointsReusingArrays( + continuousDataX, newOrderings[i]); + // Perform new initialisations + miSurrogateCalculator.initialise( + dimensionsContinuous, base, dimensionsConditional); + // Set new observations + miSurrogateCalculator.setObservations(shuffledData, + discreteData, conditionedDataZ); + } + // Compute the MI + surrogateMeasurements[i] = miSurrogateCalculator.computeAverageLocalOfObservations(); + if (debug){ + System.out.println("New MI was " + surrogateMeasurements[i]); + } + } + + return new EmpiricalMeasurementDistribution(surrogateMeasurements, condMi); + } + + public void setDebug(boolean debug) { + this.debug = debug; + } + + public double getLastAverage() { + return condMi; + } + + public int getNumObservations() { + return totalObservations; + } +} diff --git a/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian.java new file mode 100755 index 0000000..69d8f28 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian.java @@ -0,0 +1,219 @@ +package infodynamics.measures.continuous.gaussian; + +import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon; +import infodynamics.utils.MatrixUtils; + +/** + *

Computes the differential conditional mutual information of a given multivariate set of + * observations with a discrete variable, conditioned on another multivariate set + * of observations, + * assuming that the probability distribution function for these continuous observations is + * a multivariate Gaussian distribution.

+ * + *

This is done by examining the conditional probability distribution for + * the multivariate continuous variable C (given the discrete + * variable D) against the probability distribution for C, all conditioned + * on another continuous variable Z: + * MI(C;D|Z) := H(C|Z) - H(C|D,Z).

+ * + *

CAVEAT EMPTOR: The real question this type of calculation asks + * is to what extent does knowing the value of the discrete variable reduce + * variance in the continuous variable(s), given the other known continuous variable. + * Indeed, it may not to behave (I should test this--) + * as we would normally expect a mutual information calculation: if we add more + * continuous variables in, it may increase in spite of redundancy between these + * variables. TODO Further exploration should take place here ...

+ * + *

+ * Usage: + *

    + *
  1. Construct {@link #ConditionalMutualnfoCalculatorMultiVariateWithDiscreteSourceGaussian()}
    2. + *
  2. {@link #initialise(int, int, int)}
    3. + *
  3. Set properties using {@link #setProperty(String, String)}
    4. + *
  4. Provide the observations to the calculator using: + * {@link #setObservations(double[][], int[], double[][])}, or + * a sequence of: + * {@link #startAddObservations()}, + * multiple calls to {@link #addObservations(double[][], int[], double[][])} + * and then + * {@link #finaliseAddObservations()}.
    5. + *
  5. Compute the required information-theoretic results, primarily: + * {@link #computeAverageLocalOfObservations()} to return the average differential + * entropy based on the variance of + * the supplied observations; or other calls to compute + * local values or statistical significance.
    6. + *
+ *

+ * + * @see Differential entropy for Gaussian random variables at Mathworld + * @see Differential entropy for Gaussian random variables at Wikipedia + * @see Multivariate normal distribution on Wikipedia + * @author Joseph Lizier joseph.lizier_at_gmail.com + * + */ +public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian + extends ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon + implements Cloneable { + + /** + * Entropy calculator applied to the whole set of conditional data Z + */ + protected EntropyCalculatorMultiVariateGaussian entCalcZ; + /** + * Entropy calculator applied to the whole set of continuous data C + * and conditional data Z + * + */ + protected EntropyCalculatorMultiVariateGaussian entCalcCZ; + /** + * Entropy calculators applied to the set of conditional data Z + * associated with each discrete value + */ + protected EntropyCalculatorMultiVariateGaussian[] entCalcZForEachDiscrete; + /** + * Entropy calculators applied to the set of continuous data C + * and conditional data Z + * associated with each discrete value + */ + protected EntropyCalculatorMultiVariateGaussian[] entCalcCZForEachDiscrete; + + public ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian() { + super(); + entCalcZ = new EntropyCalculatorMultiVariateGaussian(); + entCalcZForEachDiscrete = null; + entCalcCZ = new EntropyCalculatorMultiVariateGaussian(); + entCalcCZForEachDiscrete = null; + } + + /* (non-Javadoc) + * @see infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon#initialise(int, int, int) + */ + @Override + public void initialise(int dimensions, int base, int dimensionsCond) { + super.initialise(dimensions, base, dimensionsCond); + entCalcZ.initialise(dimensionsCond); + entCalcCZ.initialise(dimensions + dimensionsCond); + entCalcZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + entCalcCZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + for (int b = 0; b < base; b++) { + entCalcZForEachDiscrete[b] = new EntropyCalculatorMultiVariateGaussian(); + entCalcCZForEachDiscrete[b] = new EntropyCalculatorMultiVariateGaussian(); + // If any properties relevant for these calculators were set in + // setProperty then we should set them here + entCalcZForEachDiscrete[b].initialise(dimensionsCond); + entCalcCZForEachDiscrete[b].initialise(dimensions + dimensionsCond); + } + } + + /* (non-Javadoc) + * @see infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon#finaliseAddObservations() + */ + @Override + public void finaliseAddObservations() throws Exception { + super.finaliseAddObservations(); + + // Set the complete set of observations: + // (this will pick up any errors in the dimensions of the continuous + // observations) + entCalcZ.setObservations(conditionedDataZ); + double[][] joinedCZ = MatrixUtils.appendColumns(continuousDataX, conditionedDataZ); + entCalcCZ.setObservations(joinedCZ); + + // Set the observations corresponding to each discrete value: + int totalNumberOfSuppliedObservations = 0; + for (int b = 0; b < base; b++) { + // Extract the observations for when this base value occurs: + double[][] obsZForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition( + conditionedDataZ, discreteData, b); + double[][] obsCZForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition( + joinedCZ, discreteData, b); + // Set the observations for each discrete value: + entCalcZForEachDiscrete[b].setObservations(obsZForThisDiscValue); + entCalcCZForEachDiscrete[b].setObservations(obsCZForThisDiscValue); + totalNumberOfSuppliedObservations += obsZForThisDiscValue.length; + } + // Check that all of the supplied observations were extracted corresponding + // to one of the allowed discrete values + if (totalNumberOfSuppliedObservations != discreteData.length) { + throw new Exception("Some values in discreteObservations were not in the range 0..base-1"); + } + } + + public double computeAverageLocalOfObservations() throws Exception { + // The average mutual information can be expressed + // as a difference between the conditional entropy of the continuous observations + // and the conditional entropy of the continuous given the + // discrete observations, both given the conditional continuous observations: + // I(C;D|Z) = H(C|Z) - H(C|D,Z) + // = H(C,Z) - H(Z) - H(C,Z|D) + H(C|D) + double meanConditionalEntropyCZ = 0; + double meanConditionalEntropyZ = 0; + for (int b = 0; b < base; b++) { + double pOfB = (double) counts[b] / (double) totalObservations; + meanConditionalEntropyZ += pOfB * + entCalcZForEachDiscrete[b].computeAverageLocalOfObservations(); + meanConditionalEntropyCZ += pOfB * + entCalcCZForEachDiscrete[b].computeAverageLocalOfObservations(); + } + double entCZ = entCalcCZ.computeAverageLocalOfObservations(); + double entZ = entCalcZ.computeAverageLocalOfObservations(); + condMi = entCZ - entZ + - meanConditionalEntropyCZ + + meanConditionalEntropyZ; + if (debug) { + System.out.printf("H(C,Z)=%.4f - H(Z)=%.4f - H(C,Z|D)=%.4f + H(Z|D)=%.4f = %.4f\n", + entCZ, entZ, meanConditionalEntropyCZ, meanConditionalEntropyZ, condMi); + } + return condMi; + } + + public double[] computeLocalOfPreviousObservations() throws Exception { + throw new RuntimeException("Not implemented yet"); + } + + public double[] computeLocalUsingPreviousObservations( + double[][] contStates, int[] discreteStates, + double[][] conditionedStates) throws Exception { + throw new RuntimeException("Not implemented yet"); + } + + /** + * Clone the object - note: while it does create new cloned instances of + * the {@link EntropyCalculatorMultiVariateGaussian} objects, + * I think these only + * have shallow copies to the data. + * This is enough though to maintain the structure across + * various {@link #computeSignificance(boolean, int)} calls. + * + * @see java.lang.Object#clone() + */ + @Override + protected Object clone() throws CloneNotSupportedException { + ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian theClone = + (ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian) super.clone(); + // Now assign clones of the EntropyCalculatorMultiVariateGaussian objects: + // First clone those for the conditional variable: + theClone.entCalcZ = + (EntropyCalculatorMultiVariateGaussian) entCalcZ.clone(); + if (entCalcZForEachDiscrete != null) { + theClone.entCalcZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + for (int b = 0; b < base; b++) { + theClone.entCalcZForEachDiscrete[b] = + (EntropyCalculatorMultiVariateGaussian) + entCalcZForEachDiscrete[b].clone(); + } + } + // Next clone those for the conditional variable and the continuous source variable + theClone.entCalcCZ = + (EntropyCalculatorMultiVariateGaussian) entCalcCZ.clone(); + if (entCalcCZForEachDiscrete != null) { + theClone.entCalcCZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + for (int b = 0; b < base; b++) { + theClone.entCalcCZForEachDiscrete[b] = + (EntropyCalculatorMultiVariateGaussian) + entCalcCZForEachDiscrete[b].clone(); + } + } + return theClone; + } +}