diff --git a/java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java b/java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java index 40eb210..828b268 100755 --- a/java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java +++ b/java/source/infodynamics/measures/continuous/EntropyCalculatorMultiVariate.java @@ -18,7 +18,7 @@ public interface EntropyCalculatorMultiVariate { public double[] computeLocalUsingPreviousObservations(double states[][]) throws Exception; - public double[] computeLocalOfPreviousObservations(); + public double[] computeLocalOfPreviousObservations() throws Exception; public double getLastAverage(); diff --git a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java index b372713..8d0eabb 100755 --- a/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java +++ b/java/source/infodynamics/measures/continuous/MutualInfoCalculatorMultiVariateWithDiscrete.java @@ -2,27 +2,158 @@ package infodynamics.measures.continuous; import infodynamics.utils.EmpiricalMeasurementDistribution; +/** + *

Interface defining computation of the differential mutual information between a given multivariate set of + * observations + * and a discrete variable. + * This is done by examining the conditional probability distribution (given the discrete + * variable) against the probability distribution for the multivariate set.

+ * + *

+ * Usage intended for the child classes: + *

    + *
  1. Construct the child class
  2. + *
  3. {@link #initialise(int, int)} to tell the class the number of dimensions of the continuous observations, + * and the base (number of available states) of the discrete variable.
  4. + *
  5. Set properties using {@link #setProperty(String, String)}
  6. + *
  7. Provide the observations to the calculator using: + * {@link #setObservations(double[][], int[])}.
  8. + *
  9. Compute the required information-theoretic results, primarily: + * {@link #computeAverageLocalOfObservations()} to return the average + * MI; or other calls to compute + * local values or statistical significance.
  10. + *
+ *

+ * + * @author Joseph Lizier joseph.lizier_at_gmail.com + * + */ public interface MutualInfoCalculatorMultiVariateWithDiscrete { + /** + * Initialise the calculator for use or reuse. + * This clears any previously supplied observations, but + * preserves the supplied properties. + * + * @param dimensions number of joint continuous variables + * @param base number of states in the discrete observations + * @throws Exception + */ public void initialise(int dimensions, int base) throws Exception; + /** + *

Set the required property of the calculator to the given value.

+ * + *

There are no general properties settable on all child classes; + * each child class may define their own properties.

+ * + * @param propertyName name of property + * @param propertyValue value of property + */ public void setProperty(String propertyName, String propertyValue); + /** + * Set the properties from which the mutual information should be computed. + * + * @param continuousObservations observations of the joint continuous variables; + * first index is time or observation number, second index is variable number + * (the number of variables should be equal to the number of dimensions set + * in {@link #initialise(int, int)}.) + * @param discreteObservations observations of the discrete variable; must be + * the same number of observations as supplied for continuousObservations. + * Each observation must lie in the range 0..base-1, where base was set by + * {@link #initialise(int, int)}. + * @throws Exception + */ public void setObservations(double[][] continuousObservations, int[] discreteObservations) throws Exception; + /** + * Compute the average mutual information from the previously supplied + * observations via {@link #setObservations(double[][], int[])} + * + * @return a scalar for the average mutual information + * @throws Exception + */ public double computeAverageLocalOfObservations() throws Exception; + /** + * Compute the local mutual information for each pair of continuous + * and discrete values supplied here, using probability distribution + * functions generated using the observations previously supplied + * to the method {@link #setObservations(double[][], int[])}. + * + * @param contStates observations of the joint continuous variables, + * as specified in {@link #setObservations(double[][], int[])} + * @param discreteStates observations of the discrete variable, + * as specified in {@link #setObservations(double[][], int[])} + * @return a time series of the local mutual information values + * for each supplied observation + * @throws Exception + */ public double[] computeLocalUsingPreviousObservations(double[][] contStates, int[] discreteStates) throws Exception; + /** + *

Compute the empiricial statistical significance of the mutual information + * of the previously supplied observations + * via {@link #setObservations(double[][], int[])}. + * We destroy the p(x,y) correlations, while retaining the p(x), p(y) marginals, to check how + * significant this mutual information actually was. + * Specifically, this method returns an {@link EmpiricalMeasurementDistribution} + * object containing numPermutationsToCheck surrogate measurements + * where the discrete data is shuffled against the continuous data set. + *

+ * + *

This is in the spirit of Chavez et. al. + * which was performed for Transfer entropy. + *

+ * + * @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: + * application to epileptic EEG signals', Journal of Neuroscience Methods 124 (2003) 113-128" + * @throws Exception + */ public EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception; + /** + *

Compute the empiricial statistical significance of the mutual information + * of the previously supplied observations + * via {@link #setObservations(double[][], int[])}. + * This method performs as per {@link #computeSignificance(int)} except + * that the shuffling of the discrete time series is + * specified here by the newOrderings parameter. + *

+ * + * @param newOrderings the specific new orderings to use. First index + * is for the new ordering number; second index is, for that + * particular reordering, which time point of the original series + * to grab at that point. + * @return the proportion of MI scores from the distribution which have higher or equal MIs to ours. + * @see #computeSignificance(int) + */ public EmpiricalMeasurementDistribution computeSignificance(int[][] newOrderings) throws Exception; + /** + * Set whether to print extra debug messages + * + * @param debug whether to print extra debug messages + */ public void setDebug(boolean debug); + /** + * Return the previously computed average + * + * @return the previously computed average + */ public double getLastAverage(); + /** + * Get the number of observations supplied via + * {@link #setObservations(double[][], int[])} + * + * @return the number of supplied observations + */ public int getNumObservations(); } diff --git a/java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java index f23b34d..50a1823 100755 --- a/java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java +++ b/java/source/infodynamics/measures/continuous/gaussian/EntropyCalculatorMultiVariateGaussian.java @@ -25,13 +25,20 @@ import infodynamics.utils.MatrixUtils; * @author Joseph Lizier joseph.lizier_at_gmail.com * */ -public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorMultiVariate { +public class EntropyCalculatorMultiVariateGaussian + implements EntropyCalculatorMultiVariate, Cloneable { /** * Covariance matrix of the most recently supplied observations */ protected double[][] covariance; + /** + * Means of the most recently supplied observations (source variables + * listed first, destination variables second). + */ + protected double[] means; + /** * The set of observations, retained in case the user wants to retrieve the local * entropy values of these @@ -43,6 +50,11 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM */ protected int dimensions; + /** + * Determinant of the covariance matrix; stored to save computation time + */ + protected double detCovariance; + protected double lastAverage; protected boolean debug; @@ -59,8 +71,10 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM */ public void initialise(int dimensions) { covariance = null; + means = null; observations = null; this.dimensions = dimensions; + detCovariance = 0; } /** @@ -71,7 +85,9 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM */ public void setObservations(double[][] observations) { this.observations = observations; - covariance = MatrixUtils.covarianceMatrix(observations); + means = MatrixUtils.means(observations); + covariance = MatrixUtils.covarianceMatrix(observations, means); + detCovariance = 0; // Check that the observations was of the correct number of dimensions: // (done afterwards since the covariance matrix computation checks that // all rows had the right number of columns @@ -81,12 +97,16 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM } /** - * Set the covariance of the distribution for which we will compute the - * entropy. + *

Set the covariance of the distribution for which we will compute the + * entropy.

* + *

Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}.

+ * * @param covariance covariance matrix */ public void setCovariance(double[][] covariance) throws Exception { + detCovariance = 0; observations = null; // Make sure the supplied covariance matrix is square: int rows = covariance.length; @@ -102,6 +122,25 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM this.covariance = covariance; } + /** + *

Set the covariance and mean of the distribution for which we will compute the + * entropy.

+ * + *

Note that without setting any observations, you cannot later + * call {@link #computeLocalOfPreviousObservations()}.

+ * + * @param covariance covariance matrix of the variables + * @param means mean of the variables + * @throws Exception where the dimensions of the covariance or means are not correct + */ + public void setCovarianceAndMeans(double[][] covariance, double[] means) throws Exception { + if (means.length != dimensions) { + throw new Exception("Supplied mean matrix does not match initialised number of dimensions"); + } + this.means = means; + setCovariance(covariance); + } + /** *

The joint entropy for a multivariate Gaussian-distribution of dimension n * with covariance matrix C is 0.5*\log_e{(2*pi*e)^n*|det(C)|}, @@ -116,8 +155,9 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM */ public double computeAverageLocalOfObservations() { try { + detCovariance = MatrixUtils.determinant(covariance); lastAverage = 0.5 * (dimensions* (1 + Math.log(2.0*Math.PI)) + - Math.log(Math.abs(MatrixUtils.determinant(covariance)))); + Math.log(Math.abs(detCovariance))); return lastAverage; } catch (Exception e) { // Should not happen, since we check the validity of the supplied @@ -151,27 +191,81 @@ public class EntropyCalculatorMultiVariateGaussian implements EntropyCalculatorM return lastAverage; } + /** + * Compute the local entropy for each of the given supplied + * observations + * + * @param states joint vectors of observations; first index + * is observation number, second is variable number. + * @return an array of local values + * @see Multivariate normal distribution on Wikipedia + */ public double[] computeLocalUsingPreviousObservations(double[][] states) throws Exception { - // TODO Implement me - throw new RuntimeException("Not implemented yet"); + if (means == null) { + throw new Exception("Cannot compute local values without having means either supplied or computed via setObservations()"); + } + + // Check that the covariance matrix was positive definite: + // it is known (i think because it is symmetric) that the covariance + // matrix will be positive definite unless one variable is an exact + // linear combination of the others (ref: wikipedia, above). + // We can check for linear dependence by ensuring that the determinant + // of the covariance matrix is non-zero: + if (detCovariance == 0) { + // We need to check: + detCovariance = MatrixUtils.determinant(covariance); + if (detCovariance == 0) { + throw new Exception("Covariance matrix is not positive definite"); + } + } + // Now we are clear to take the matrix inverse (via Cholesky decomposition, + // since we have a symmetric positive definite matrix): + double[][] invCovariance = MatrixUtils.invertSymmPosDefMatrix(covariance); + + // If we have a time delay, slide the local values + double[] localValues = new double[states.length]; + for (int t = 0; t < states.length; t++) { + double[] deviationsFromMean = + MatrixUtils.subtract(states[t], means); + // Computing PDF + // (see the PDF defined at the wikipedia page referenced in the method header) + double jointExpArg = MatrixUtils.dotProduct( + MatrixUtils.matrixProduct(deviationsFromMean, + invCovariance), + deviationsFromMean); + double pJoint = Math.pow(2.0 * Math.PI, -(double) dimensions / 2.0) * + Math.exp(-0.5 * jointExpArg) / + Math.sqrt(detCovariance); + localValues[t] = - Math.log(pJoint); + } + + // Don't set average if this was the previously supplied observations, + // since it won't be the same as what would have been computed + // analytically. + + return localValues; } - public double[] computeLocalOfPreviousObservations() { - - // TODO Implement this function - if (true) - throw new RuntimeException("Not implemented yet"); - + public double[] computeLocalOfPreviousObservations() throws Exception { if (observations == null) { - throw new RuntimeException("Cannot compute local values since no observations were supplied"); + throw new Exception("Cannot compute local values since no observations were supplied"); } - double[] localEntropy = new double[observations.length]; - for (int t=0; t < observations.length; t++) { - // Compute the probability for the given observation, based on - // the assumption of a multivariate Gaussian PDF: - - } - return null; + return computeLocalUsingPreviousObservations(observations); + } + + /** + * Provide an implementation of the clone() method. + * This does not deeply copy all of the underlying data, just providing + * a copy of the references to it all. + * This is enough to protect the integrity of the calculator + * however if the clone is supplied different data (though the + * clone should not alter the data). + * + * @see java.lang.Object#clone() + */ + @Override + protected Object clone() throws CloneNotSupportedException { + return super.clone(); } } diff --git a/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java index 92e91fd..91d254b 100755 --- a/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java +++ b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java @@ -34,8 +34,9 @@ import infodynamics.utils.MatrixUtils; * *

* - * @see Differential entropy for Gaussian random variables defined at - * {@link http://mathworld.wolfram.com/DifferentialEntropy.html} + * @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 * */ @@ -147,6 +148,8 @@ public class MutualInfoCalculatorMultiVariateGaussian * @param covariance covariance matrix of the source and destination * variables, considered together (variable indices start with the source * and continue into the destination). + * @param means mean of the source and destination variables (as per + * covariance) */ public void setCovarianceAndMeans(double[][] covariance, double[] means) throws Exception { this.means = means; @@ -306,7 +309,8 @@ public class MutualInfoCalculatorMultiVariateGaussian * If the {@link MutualInfoCalculatorMultiVariate#PROP_TIME_DIFF} * property was set to say k, then the local values align with the * destination value (i.e. after the given delay k). As such, the - * first k values of the array will be zeros. + * first k values of the array will be zeros. + * @see Multivariate normal distribution on Wikipedia * @throws Exception */ protected double[] computeLocalUsingPreviousObservations(double[][] newSourceObs, @@ -360,7 +364,6 @@ public class MutualInfoCalculatorMultiVariateGaussian } // If we have a time delay, slide the local values double[] localValues = new double[lengthOfReturnArray]; - double newAverage = 0.0; for (int t = offset; t < newDestObs.length; t++) { // Computing local values for: // a. sourceObservations[t - offset] @@ -376,6 +379,7 @@ public class MutualInfoCalculatorMultiVariateGaussian destDeviationsFromMean); // Computing PDFs WITHOUT (2*pi)^dim factor, since these will cancel: + // (see the PDFs defined at the wikipedia page referenced in the method header) double sourceExpArg = MatrixUtils.dotProduct( MatrixUtils.matrixProduct(sourceDeviationsFromMean, invSourceCovariance), @@ -398,17 +402,13 @@ public class MutualInfoCalculatorMultiVariateGaussian // Returning results in nats: double localValue = Math.log(adjustedPJoint / (adjustedPSource * adjustedPDest)); - localValues[t] = localValue; - - if (isPreviousObservations) { - newAverage += localValue; - } + localValues[t] = localValue; } - if (isPreviousObservations) { - lastAverage = newAverage / (newDestObs.length - timeDiff); - miComputed = true; - } + // if (isPreviousObservations) { + // Don't store the average value here, since it won't be exactly + // the same as what would have been computed under the analytic expression + // } return localValues; } diff --git a/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateWithDiscreteGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateWithDiscreteGaussian.java new file mode 100755 index 0000000..8d2d466 --- /dev/null +++ b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateWithDiscreteGaussian.java @@ -0,0 +1,308 @@ +package infodynamics.measures.continuous.gaussian; + +import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariateWithDiscrete; +import infodynamics.utils.EmpiricalMeasurementDistribution; +import infodynamics.utils.MatrixUtils; +import infodynamics.utils.RandomGenerator; + +/** + *

Computes the differential mutual information between a given multivariate set of + * observations + * (assuming that the probability distribution function for these observations is + * a multivariate Gaussian distribution) + * and a discrete variable. + * This is done by examining the conditional probability distribution (given the discrete + * variable) against the probability distribution for the mulitvariate set.

+ * + *

+ * Usage: + *

    + *
  1. Construct {@link #MutualInfoCalculatorMultiVariateWithDiscreteGaussian()}
  2. + *
  3. {@link #initialise(int, int)}
  4. + *
  5. Set properties using {@link #setProperty(String, String)}
  6. + *
  7. Provide the observations to the calculator using: + * {@link #setObservations(double[][], int[])}, or + * {@link #setCovariances(double[][], double[][])}.
  8. + *
  9. Compute the required information-theoretic results, primarily: + * {@link #computeAverageLocalOfObservations()} to return the average differential + * entropy based on either the set variance or the variance of + * the supplied observations; or other calls to compute + * local values or statistical significance.
  10. + *
+ *

+ * + * @see Differential entropy for Gaussian random variables defined at + * {@link http://mathworld.wolfram.com/DifferentialEntropy.html} + * @author Joseph Lizier joseph.lizier_at_gmail.com + * + */ +public class MutualInfoCalculatorMultiVariateWithDiscreteGaussian implements + MutualInfoCalculatorMultiVariateWithDiscrete, Cloneable { + + /** + * Entropy calculator applied to the whole set of continuous data + */ + protected EntropyCalculatorMultiVariateGaussian entCalc; + /** + * Entropy calculators applied to the set of continuous data + * associated with each discrete value + */ + protected EntropyCalculatorMultiVariateGaussian[] entCalcForEachDiscrete; + + /** + * Keep a copy of the discrete observations, to enable us to compute the + * statistical significance later. We don't need to keep a copy of + * the continuous observations, since they're kept intact by entCalc. + */ + protected int[] discreteObservations; + + /** + * Number of supplied observations + */ + protected int totalObservations = 0; + + /** + * The number of possible discrete states + */ + protected int base = 0; + + /** + * Whether to print extra debug messages + */ + protected boolean debug = false; + + /** + * The last computed average MI value + */ + protected double lastAverage = 0; + + public MutualInfoCalculatorMultiVariateWithDiscreteGaussian() { + entCalc = new EntropyCalculatorMultiVariateGaussian(); + entCalcForEachDiscrete = null; + } + + public void initialise(int dimensions, int base) throws Exception { + totalObservations = 0; + lastAverage = 0; + discreteObservations = null; + this.base = base; + entCalc.initialise(dimensions); + entCalcForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + for (int b = 0; b < base; b++) { + entCalcForEachDiscrete[b] = new EntropyCalculatorMultiVariateGaussian(); + // If any properties relevant for these calculators were set in + // setProperty then we should set them here + entCalcForEachDiscrete[b].initialise(dimensions); + } + } + + /** + *

Set the required property to the given value.

+ * + *

At this stage, there are no settable properties for this calculator.

+ * + * @param propertyName name of property + * @param propertyValue value of property + */ + public void setProperty(String propertyName, String propertyValue) { + // No properties for this calculator + } + + public void setObservations(double[][] continuousObservations, + int[] discreteObservations) throws Exception { + if (continuousObservations.length != discreteObservations.length) { + throw new Exception("Observations are not of the same length"); + } + totalObservations = continuousObservations.length; + // Set the complete set of observations: + entCalc.setObservations(continuousObservations); + // Set the observations corresponding to each discrete value: + setDiscreteData(continuousObservations, discreteObservations); + } + + protected void setDiscreteData(double continuousObservations[][], + int discreteObservations[]) throws Exception { + int totalNumberOfSuppliedObservations = 0; + for (int b = 0; b < base; b++) { + // Extract the observations for when this base value occurs: + double[][] obsForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition( + continuousObservations, discreteObservations, b); + // Set the observations for each discrete value: + entCalcForEachDiscrete[b].setObservations(obsForThisDiscValue); + totalNumberOfSuppliedObservations += obsForThisDiscValue.length; + } + // Check that all of the supplied observations were extracted corresponding + // to one of the allowed discrete values + if (totalNumberOfSuppliedObservations != discreteObservations.length) { + throw new Exception("Some values in discreteObservations were not in the range 0..base-1"); + } + this.discreteObservations = discreteObservations; + } + + public double computeAverageLocalOfObservations() throws Exception { + // The average mutual information can be expressed + // as a difference between the entropy of the continuous observations + // and the conditional entropy of the continuous given the + // discrete observations: + // I(C;D) = H(C) - H(C|D) + // Subtract that conditional entropy from the + // entropy of all observations: + lastAverage = entCalc.computeAverageLocalOfObservations() + - computeAverageLocalConditionalEntropyOfObservations(); + return lastAverage; + } + + /** + * Compute the average conditional entropy of the continuous data given + * the discrete data (averaged over all discrete values) + * + * @return average conditional entropy + */ + protected double computeAverageLocalConditionalEntropyOfObservations() { + double meanConditionalEntropy = 0; + for (int b = 0; b < base; b++) { + double pOfB = (double) entCalcForEachDiscrete[b].observations.length / + (double) totalObservations; + meanConditionalEntropy += pOfB * + entCalcForEachDiscrete[b].computeAverageLocalOfObservations(); + } + return meanConditionalEntropy; + } + + public double[] computeLocalUsingPreviousObservations( + double[][] contStates, int[] discreteStates) throws Exception { + + // The local mutual information can be expressed + // as a difference between the local + // entropy of the continuous observations + // and the local conditional entropy of the continuous given the + // discrete observations: + // i(C;D) = h(C) - h(C|D) + // First compute the local entropy for the continuous + // observations: + double[] localValues = entCalc.computeLocalUsingPreviousObservations(contStates); + // Next compute the local conditional entropes from each + // conditional entropy calculator: + double[][] localConditionalEntropies = new double[base][]; + for (int b = 0; b < base; b++) { + // Extract the observations for when this base value occurs: + double[][] obsForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition( + contStates, discreteStates, b); + // and compute the local conditional entropies for these time points: + localConditionalEntropies[b] = + entCalcForEachDiscrete[b]. + computeLocalUsingPreviousObservations(obsForThisDiscValue); + } + // Now subtract the correct local conditional entropy from the local + // entropy of the continuous observations at the correct time points + int[] nextTForDiscreteState = new int[base]; + for (int t = 0; t < contStates.length; t++) { + // Find the local conditional entropy value at this time point: + // 1. Grab the current discrete value + int b = discreteStates[t]; + // 2. Pick out the next index in the local conditional entropies + // for this discrete value, nextTForDiscreteState[b], and pull + // out the local value at this index: + double localCondEntropy = + localConditionalEntropies[b][nextTForDiscreteState[b]]; + // 3. Update the next index for this discrete value: + nextTForDiscreteState[b]++; + // Now finalise the local MI at this point: + localValues[t] -= localCondEntropy; + } + + return localValues; + } + + public EmpiricalMeasurementDistribution computeSignificance( + int numPermutationsToCheck) throws Exception { + if (totalObservations == 0) { + throw new Exception("Must have set observations before computing significance"); + } + // Generate the re-ordered indices: + RandomGenerator rg = new RandomGenerator(); + int[][] newOrderings = rg.generateDistinctRandomPerturbations(totalObservations, numPermutationsToCheck); + return computeSignificance(newOrderings); + } + + public EmpiricalMeasurementDistribution computeSignificance( + int[][] newOrderings) throws Exception { + int numPermutationsToCheck = newOrderings.length; + if (lastAverage == 0) { + // This may execute even if the average was computed and found + // to be == 0; this won't hurt, just costs a tiny bit of execution + // time. + 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) + MutualInfoCalculatorMultiVariateWithDiscreteGaussian miSurrogateCalculator = + (MutualInfoCalculatorMultiVariateWithDiscreteGaussian) clone(); + + double[] surrogateMeasurements = new double[numPermutationsToCheck]; + + // Now compute the MI for each set of shuffled data: + for (int i = 0; i < numPermutationsToCheck; i++) { + // Generate a new re-ordered discrete data + int[] shuffledDiscreteData = + MatrixUtils.extractSelectedTimePoints( + discreteObservations, newOrderings[i]); + // Re-initialise the data in the conditional entropy calculators: + // (in theory, we should call initialise and properly call + // setObservations(), but we know we can short circuit that from + // inside this calculator, and avoid recomputing covariances etc + // on the continuous data set) + miSurrogateCalculator.setDiscreteData( + entCalc.observations, shuffledDiscreteData); + // Compute the MI: + surrogateMeasurements[i] = miSurrogateCalculator.entCalc.computeAverageLocalOfObservations() - + miSurrogateCalculator.computeAverageLocalConditionalEntropyOfObservations(); + if (debug){ + System.out.println("New MI was " + surrogateMeasurements[i]); + } + } + + return new EmpiricalMeasurementDistribution(surrogateMeasurements, lastAverage); + } + + public void setDebug(boolean debug) { + this.debug = debug; + } + + public double getLastAverage() { + return lastAverage; + } + + public int getNumObservations() { + return totalObservations; + } + + /** + * 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(int)} calls. + * + * @see java.lang.Object#clone() + */ + @Override + protected Object clone() throws CloneNotSupportedException { + MutualInfoCalculatorMultiVariateWithDiscreteGaussian theClone = + (MutualInfoCalculatorMultiVariateWithDiscreteGaussian) super.clone(); + // Now assign clones of the EntropyCalculatorMultiVariateGaussian objects: + theClone.entCalc = + (EntropyCalculatorMultiVariateGaussian) entCalc.clone(); + if (entCalcForEachDiscrete != null) { + theClone.entCalcForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base]; + for (int b = 0; b < base; b++) { + theClone.entCalcForEachDiscrete[b] = + (EntropyCalculatorMultiVariateGaussian) + entCalcForEachDiscrete[b].clone(); + } + } + return theClone; + } +} diff --git a/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateKernel.java b/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateKernel.java index cfd9a00..517ccd8 100755 --- a/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateKernel.java +++ b/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateKernel.java @@ -608,7 +608,7 @@ public class MutualInfoCalculatorMultiVariateKernel /** * Clone the object - note: while it does create new cloned instances of - * the MulitVariateKernelEstimator objects, I think these only + * the {@link KernelEstimatorMultiVariate} objects, I think these only * have shallow copies to the data. * This is enough though to maintain the structure across * various {@link #computeSignificance(int)} calls. diff --git a/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateWithDiscreteKernel.java b/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateWithDiscreteKernel.java index 73646c2..a173da7 100755 --- a/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateWithDiscreteKernel.java +++ b/java/source/infodynamics/measures/continuous/kernel/MutualInfoCalculatorMultiVariateWithDiscreteKernel.java @@ -10,9 +10,9 @@ import infodynamics.utils.RandomGenerator; public class MutualInfoCalculatorMultiVariateWithDiscreteKernel implements MutualInfoCalculatorMultiVariateWithDiscrete { - KernelEstimatorMultiVariate mvke = null; - KernelEstimatorMultiVariate[] mvkeForEachDiscrete = null; - int base = 0; + protected KernelEstimatorMultiVariate mvke = null; + protected KernelEstimatorMultiVariate[] mvkeForEachDiscrete = null; + protected int base = 0; private int totalObservations = 0; // private int dimensions1 = 0; @@ -207,18 +207,6 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKernel implements return newMI; } - /** - * Compute the significance of the mutual information of the previously supplied observations. - * We destroy the p(x,y) correlations, while retaining the p(x), p(y) marginals, to check how - * significant this 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 numPermutationsToCheck - * @return the proportion of MI scores from the distribution which have higher or equal MIs to ours. - */ public synchronized EmpiricalMeasurementDistribution computeSignificance(int numPermutationsToCheck) throws Exception { // Generate the re-ordered indices: RandomGenerator rg = new RandomGenerator(); @@ -226,18 +214,6 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKernel implements return computeSignificance(newOrderings); } - /** - * Compute the significance of the mutual information of the previously supplied observations. - * We destroy the p(x,y) correlations, while retaining the p(x), p(y) marginals, to check how - * significant this 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 newOrderings the specific new orderings to use - * @return the proportion of MI scores from the distribution which have higher or equal MIs to ours. - */ public EmpiricalMeasurementDistribution computeSignificance(int[][] newOrderings) throws Exception { int numPermutationsToCheck = newOrderings.length; if (!miComputed) { @@ -540,19 +516,20 @@ public class MutualInfoCalculatorMultiVariateWithDiscreteKernel implements } /** - * Set properties for the mutual information calculator. + *

Set properties for the mutual information calculator. * These can include: *

+ *

* - * Note that dynamic correlation exclusion may have unexpected results if multiple + *

Note that dynamic correlation exclusion may have unexpected results if multiple * observation sets have been added. This is because multiple observation sets * are treated as though they are from a single time series, so observations from * near the end of observation set i will be excluded from comparison to * observations near the beginning of observation set (i+1). + *

* * @param propertyName * @param propertyValue diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index 12bad56..31d193b 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -2359,11 +2359,10 @@ public class MatrixUtils { /** *

Returns the covariance between the first two columns of data.

- *

See - Mathworld - *

* * @param data * @return the covariance + * @see Mathworld */ public static double covarianceFirstTwoColumns(double[][] data) { return covarianceTwoColumns(data, 0, 1); @@ -2445,13 +2444,21 @@ public class MatrixUtils { * @return covariance matrix */ public static double[][] covarianceMatrix(double[][] data) { + return covarianceMatrix(data, means(data)); + } + + /** + * Compute the covariance matrix between all column pairs (variables) in the + * multivariate data set + * + * @param data multivariate array of data; first index is time, second is + * variable number + * @param means the mean of each variable (column) in the data + * @return covariance matrix + */ + public static double[][] covarianceMatrix(double[][] data, double[] means) { int numVariables = data[0].length; double[][] covariances = new double[numVariables][numVariables]; - // Compute means of each variable once up front to save time - double[] means = new double[numVariables]; - for (int r = 0; r < numVariables; r++) { - means[r] = mean(data, r); - } for (int r = 0; r < numVariables; r++) { for (int c = r; c < numVariables; c++) { // Compute the covariance between variable r and c: @@ -3017,6 +3024,7 @@ public class MatrixUtils { * @param B matrix with as many rows as A and any number of columns * @return X so that A*X = B * @see {@link http://math.nist.gov/javanumerics/jama/} + * @see #CholeskyDecomposition(double[][]) */ public static double[][] solveViaCholeskyResult(double[][] L, double[][] B) { int aRows = L.length;