From cd67da53a407d71d688417e187dfe6076e9befea Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Mon, 4 Mar 2013 06:08:38 +0000 Subject: [PATCH] Patching linear Gaussian estimators to give empirical statistical significance given a certain number of observations. Changed implementation of our conditional MI calculators Kraskov via algorithm 2 to my current interpretation of how this should be done. Added a utility for lagged covariance in MatrixUtils --- ...ualInfoCalculatorMultiVariateGaussian.java | 46 ++++++++++- ...ualInfoCalculatorMultiVariateGaussian.java | 12 ++- ...tualInfoCalculatorMultiVariateKraskov.java | 4 + ...ualInfoCalculatorMultiVariateKraskov2.java | 77 ++++++++++--------- ...ulatorMultiVariateWithDiscreteKraskov.java | 41 ++++++---- .../infodynamics/utils/MatrixUtils.java | 34 +++++++- .../infodynamics/utils/OctaveFileReader.java | 10 ++- 7 files changed, 159 insertions(+), 65 deletions(-) diff --git a/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java index 2c32a75..3c8cea9 100755 --- a/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java +++ b/java/source/infodynamics/measures/continuous/gaussian/ConditionalMutualInfoCalculatorMultiVariateGaussian.java @@ -4,6 +4,7 @@ import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVari import infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon; import infodynamics.utils.AnalyticNullDistributionComputer; import infodynamics.utils.ChiSquareMeasurementDistribution; +import infodynamics.utils.EmpiricalMeasurementDistribution; import infodynamics.utils.MatrixUtils; /** @@ -145,11 +146,14 @@ public class ConditionalMutualInfoCalculatorMultiVariateGaussian * * @param covariance covariance matrix of var1, var2, conditional * variables, considered together. + * @param numObservations the number of observations that the covariance + * was determined from. This is used for later significance calculations * @throws Exception for covariance matrix not matching the expected dimensions, * being non-square, asymmetric or non-positive definite */ - public void setCovariance(double[][] covariance) throws Exception { + public void setCovariance(double[][] covariance, int numObservations) throws Exception { setCovariance(covariance, false); + totalObservations = numObservations; } /** @@ -222,10 +226,14 @@ public class ConditionalMutualInfoCalculatorMultiVariateGaussian * variables, considered together. * @param means mean of var1, var2 and conditional variables (as per * covariance) + * @param numObservations the number of observations that the mean and covariance + * were determined from. This is used for later significance calculations */ - public void setCovarianceAndMeans(double[][] covariance, double[] means) throws Exception { + public void setCovarianceAndMeans(double[][] covariance, double[] means, + int numObservations) throws Exception { + this.means = means; - setCovariance(covariance); + setCovariance(covariance, numObservations); } /** @@ -309,10 +317,40 @@ public class ConditionalMutualInfoCalculatorMultiVariateGaussian // Number of extra parameters in the model incorporating the // extra variable is independent of the number of variables // in the conditional: - return new ChiSquareMeasurementDistribution(2*totalObservations*lastAverage, + return new ChiSquareMeasurementDistribution(2.0*((double)totalObservations)*lastAverage, dimensionsVar1 * dimensionsVar2); } + + + /* (non-Javadoc) + * @see infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon#computeSignificance(int, int) + */ + @Override + public EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int numPermutationsToCheck) throws Exception { + if (var2Observations == null) { + throw new Exception("Cannot compute empirical statistical significance " + + "if user passed in covariance matrix rather than observations."); + } + + return super.computeSignificance(variableToReorder, numPermutationsToCheck); + } + + /* (non-Javadoc) + * @see infodynamics.measures.continuous.ConditionalMutualInfoMultiVariateCommon#computeSignificance(int, int[][]) + */ + @Override + public EmpiricalMeasurementDistribution computeSignificance( + int variableToReorder, int[][] newOrderings) throws Exception { + if (var2Observations == null) { + throw new Exception("Cannot compute empirical statistical significance " + + "if user passed in covariance matrix rather than observations."); + } + + return super.computeSignificance(variableToReorder, newOrderings); + } + /** * @return the number of previously supplied observations for which * the conditional mutual information will be / was computed. diff --git a/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java index d2b5c55..5b7c148 100755 --- a/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java +++ b/java/source/infodynamics/measures/continuous/gaussian/MutualInfoCalculatorMultiVariateGaussian.java @@ -147,11 +147,14 @@ public class MutualInfoCalculatorMultiVariateGaussian * matrix of the source observations, C_dd is the covariance matrix * of the destination observations, and C_sd and C_ds are the covariances * of source to destination and destination to source observations. + * @param numObservations the number of observations that the covariance + * was determined from. This is used for later significance calculations * @throws Exception for covariance matrix not matching the expected dimensions, * being non-square, asymmetric or non-positive definite */ - public void setCovariance(double[][] covariance) throws Exception { + public void setCovariance(double[][] covariance, int numObservations) throws Exception { setCovariance(covariance, false); + totalObservations = numObservations; } /** @@ -218,10 +221,13 @@ public class MutualInfoCalculatorMultiVariateGaussian * and continue into the destination). * @param means mean of the source and destination variables (as per * covariance) + * @param numObservations the number of observations that the mean and covariance + * were determined from. This is used for later significance calculations */ - public void setCovarianceAndMeans(double[][] covariance, double[] means) throws Exception { + public void setCovarianceAndMeans(double[][] covariance, double[] means, + int numObservations) throws Exception { this.means = means; - setCovariance(covariance); + setCovariance(covariance, numObservations); } /** diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java index 1b99c94..75726c6 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov.java @@ -23,6 +23,10 @@ import infodynamics.utils.RandomGenerator; */ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov implements ConditionalMutualInfoCalculatorMultiVariate { + + // TODO - have this class inhereit from ConditionalMutualInfoMultiVariateCommon + // and remove the overlapping code. This will also solve a lot of the + // "Not implemented yet" exceptions that we're throwing here. /** * we compute distances to the kth neighbour in the joint space diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java index d285a9c..cae3296 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java @@ -8,14 +8,16 @@ import infodynamics.utils.MatrixUtils; *

Compute the Conditional Mutual Info using the Kraskov estimation method, * as extended by Frenzel and Pompe.

*

- * Uses the second algorithm (defined at start of p.3 of the Kraskov paper) - + * Attempts to use the second algorithm (defined at start of p.3 of the Kraskov paper) - * note that Frenzel and Pompe only extended the technique directly for the * first algorithm, though here we define it for the second. - * It is unclear exactly how to do this, since we need to account for - * the 1/k factor, which changes in the space of the second - * MI. I've taken a guess, though use of + * It is unclear exactly how to do this - we assume that we only need + * a single 1/k additive factor in the largest joint space (Kraskov seems to do the + * same for the regular MI, noting on p.5 that other correction factors are + * missing, which would go to zero as N -> \infty. + * The implementation here is something of a guess, so use of * {@link ConditionalMutualInfoCalculatorMultiVariateKraskov1} - * is perhaps recommended instead of this class.

+ * is perhaps recommended instead of this algorithm 2.

* *

Computes this directly looking at the marginal space for each variable, rather than * using the multi-info (or integration) in the marginal spaces. @@ -189,23 +191,22 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 // average: averageDiGammas += MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); - double invN_yz = 1.0 / (double) n_yz; - averageInverseCountInJointYZ += invN_yz; if (debug) { // Only tracking this for debugging purposes: double invN_xz = 1.0/(double) n_xz; averageInverseCountInJointXZ += invN_xz; + double invN_yz = 1.0 / (double) n_yz; + averageInverseCountInJointYZ += invN_yz; double localCondMi = MathsUtils.digamma(k) - 1.0 / (double) k + - invN_yz + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, 1/n_yz=%.3f, 1/n_xz=%.3f, local=%.4f\n", t, n_xz, n_yz, n_z, invN_yz, invN_xz, localCondMi); } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; condMi = MathsUtils.digamma(k) - 1.0 / (double) k + - + averageInverseCountInJointYZ + averageDiGammas; + + averageDiGammas; condMiComputed = true; if (debug) { @@ -213,11 +214,12 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 avNyz /= (double)N; avNz /= (double)N; averageInverseCountInJointXZ /= (double) N; + averageInverseCountInJointYZ /= (double) N; System.out.printf("=%.3f, =%.3f, =%.3f\n", avNxz, avNyz, avNz); - System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f + <1/n_yz>=%.3f = %.3f (<1/n_xz>=%.3f)\n", - MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, averageInverseCountInJointYZ, - condMi, averageInverseCountInJointXZ); + System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f = %.3f (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n", + MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, + condMi, averageInverseCountInJointYZ, averageInverseCountInJointXZ); } return condMi; @@ -312,35 +314,35 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 // average: averageDiGammas += MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); - double invN_yz = 1.0 / (double) n_yz; - averageInverseCountInJointYZ += invN_yz; if (debug) { // Only tracking this for debugging purposes: double invN_xz = 1.0/(double) n_xz; averageInverseCountInJointXZ += invN_xz; + double invN_yz = 1.0 / (double) n_yz; + averageInverseCountInJointYZ += invN_yz; double localCondMi = MathsUtils.digamma(k) - 1.0 / (double) k + - invN_yz + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, 1/n_yz=%.3f, 1/n_xz=%.3f, local=%.4f\n", t, n_xz, n_yz, n_z, invN_yz, invN_xz, localCondMi); } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; condMi = MathsUtils.digamma(k) - 1.0 / (double) k + - averageInverseCountInJointYZ + averageDiGammas; + averageDiGammas; condMiComputed = true; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double)N; + averageInverseCountInJointYZ /= (double) N; averageInverseCountInJointXZ /= (double) N; System.out.printf("=%.3f, =%.3f, =%.3f\n", avNxz, avNyz, avNz); - System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f + <1/n_yz>=%.3f = %.3f (<1/n_xz>=%.3f)\n", - MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, averageInverseCountInJointYZ, - condMi, averageInverseCountInJointXZ); + System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f = %.3f (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n", + MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, + condMi, averageInverseCountInJointYZ, averageInverseCountInJointXZ); } return condMi; @@ -437,35 +439,35 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 // average: averageDiGammas += MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); - double invN_yz = 1.0 / (double) n_yz; - averageInverseCountInJointYZ += invN_yz; if (debug) { // Only tracking this for debugging purposes: double invN_xz = 1.0/(double) n_xz; averageInverseCountInJointXZ += invN_xz; + double invN_yz = 1.0 / (double) n_yz; + averageInverseCountInJointYZ += invN_yz; double localCondMi = MathsUtils.digamma(k) - 1.0 / (double) k + - invN_yz + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_z) - MathsUtils.digamma(n_xz) - MathsUtils.digamma(n_yz); System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, 1/n_yz=%.3f, 1/n_xz=%.3f, local=%.4f\n", t, n_xz, n_yz, n_z, invN_yz, invN_xz, localCondMi); } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; condMi = MathsUtils.digamma(k) - 1.0 / (double) k + - averageInverseCountInJointYZ + averageDiGammas; + averageDiGammas; condMiComputed = true; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double)N; + averageInverseCountInJointYZ /= (double) N; averageInverseCountInJointXZ /= (double) N; System.out.printf("=%.3f, =%.3f, =%.3f\n", avNxz, avNyz, avNz); - System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f + <1/n_yz>=%.3f = %.3f (<1/n_xz>=%.3f)\n", - MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, averageInverseCountInJointYZ, - condMi, averageInverseCountInJointXZ); + System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f = %.3f (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n", + MathsUtils.digamma(k), averageDiGammas, 1.0 / (double) k, + condMi, averageInverseCountInJointYZ, averageInverseCountInJointXZ); } return condMi; @@ -559,36 +561,35 @@ public class ConditionalMutualInfoCalculatorMultiVariateKraskov2 double digammaNyz = MathsUtils.digamma(n_yz); double digammaNz = MathsUtils.digamma(n_z); - double invN_yz = 1.0/(double) n_yz; localCondMi[t] = digammaK - digammaNxz - digammaNyz + digammaNz - - invK + invN_yz; + - invK; averageDiGammas += digammaNz - digammaNxz - digammaNyz; - averageInverseCountInJointYZ += invN_yz; if (debug) { // Only tracking this for debugging purposes: + double invN_yz = 1.0/(double) n_yz; double invN_xz = 1.0/(double) n_xz; + averageInverseCountInJointYZ += invN_yz; averageInverseCountInJointXZ += invN_xz; System.out.printf("t=%d, n_xz=%d, n_yz=%d, n_z=%d, 1/n_yz=%.3f, 1/n_xz=%.3f, local=%.4f\n", t, n_xz, n_yz, n_z, invN_yz, invN_xz, localCondMi[t]); } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; - condMi = digammaK + averageDiGammas - invK + - averageInverseCountInJointYZ; + condMi = digammaK + averageDiGammas - invK; condMiComputed = true; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double)N; + averageInverseCountInJointYZ /= (double) N; averageInverseCountInJointXZ /= (double) N; System.out.printf("=%.3f, =%.3f, =%.3f\n", avNxz, avNyz, avNz); - System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f + <1/n_yz>=%.3f = %.3f (<1/n_xz>=%.3f)\n", - digammaK, averageDiGammas, invK, averageInverseCountInJointYZ, - condMi, averageInverseCountInJointXZ); + System.out.printf("Av = digamma(k)=%.3f + =%.3f - 1/k=%.3f = %.3f (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n", + digammaK, averageDiGammas, invK, + condMi, averageInverseCountInJointYZ, averageInverseCountInJointXZ); } return localCondMi; diff --git a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java index bed9d0d..74aa467 100755 --- a/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java +++ b/java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java @@ -14,6 +14,9 @@ import infodynamics.utils.RandomGenerator; *

Uses Kraskov method type 2, since type 1 only looks at points with * distances strictly less than the kth variable, which won't work for one marginal * being discrete.

+ *

The actual equation which should be used for type 2 with conditional MI is + * not completely clear - this calculator makes an assumption as per {#link ConditionalMutualInfoCalculatorMultiVariateKraskov2}. + *

* * @see "Estimating mutual information", Kraskov, A., Stogbauer, H., Grassberger, P., Physical Review E 69, (2004) 066138 * @see http://dx.doi.org/10.1103/PhysRevE.69.066138 @@ -225,10 +228,11 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl // Count the average number of points within eps_x and eps_y double averageDiGammas = 0; - double averageInverseCountInJointYZ = 0; double avNxz = 0; double avNyz = 0; double avNz = 0; + @SuppressWarnings("unused") // The following var is only used in debug mode: + double averageInverseCountInJointYZ = 0; for (int t = 0; t < N; t++) { // Compute eps_x and eps_z for this time step: @@ -284,20 +288,21 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl // average: averageDiGammas += MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_yz) - MathsUtils.digamma(n_z); - averageInverseCountInJointYZ += 1.0 / (double) n_yz; + if (debug) { + averageInverseCountInJointYZ += 1.0 / (double) n_yz; + } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double)N; + averageInverseCountInJointYZ /= (double) N; System.out.println(String.format("Average n_xz=%.3f, Average n_yz=%.3f, Average n_z=%.3f", avNxz, avNyz, avNz)); } - condMi = MathsUtils.digamma(k) - 1.0/(double)k + - averageInverseCountInJointYZ - averageDiGammas; + condMi = MathsUtils.digamma(k) - 1.0/(double)k - averageDiGammas; miComputed = true; return condMi; } @@ -373,25 +378,26 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl // average: averageDiGammas += MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_yz) - MathsUtils.digamma(n_z); - averageInverseCountInJointYZ += 1.0 / (double) n_yz; + if (debug) { + averageInverseCountInJointYZ += 1.0 / (double) n_yz; + } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double) N; + averageInverseCountInJointYZ /= (double) N; System.out.printf("Average n_xz=%.3f (-> digam=%.3f %.3f), Average n_yz=%.3f (-> digam=%.3f)", avNxz, MathsUtils.digamma((int) avNxz), MathsUtils.digamma((int) avNxz - 1), avNyz, MathsUtils.digamma((int) avNyz)); System.out.printf(", Average n_z=%.3f (-> digam=%.3f)\n", avNz, MathsUtils.digamma((int) avNz)); System.out.printf("Independent average num in joint box is %.3f\n", (avNxz * avNyz / (double) N)); - System.out.println(String.format("digamma(k)=%.3f - 1/k=%.3f + <1/n_yz>=%.3f - averageDiGammas=%.3f\n", - MathsUtils.digamma(k), 1.0/(double)k, averageInverseCountInJointYZ, - averageDiGammas)); + System.out.printf("digamma(k)=%.3f - 1/k=%.3f - averageDiGammas=%.3f (<1/n_yz>=%.3f)\n", + MathsUtils.digamma(k), 1.0/(double)k, + averageDiGammas, averageInverseCountInJointYZ); } - condMi = MathsUtils.digamma(k) - 1.0/(double)k + - averageInverseCountInJointYZ - averageDiGammas; + condMi = MathsUtils.digamma(k) - 1.0/(double)k - averageDiGammas; miComputed = true; return condMi; } @@ -410,10 +416,11 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl // Count the average number of points within eps_x and eps_y double averageDiGammas = 0; - double averageInverseCountInJointYZ = 0; double avNxz = 0; double avNyz = 0; double avNz = 0; + @SuppressWarnings("unused") // The following var is only used in debug mode: + double averageInverseCountInJointYZ = 0; for (int t = 0; t < N; t++) { // Compute eps_* for this time step: @@ -471,19 +478,21 @@ public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov impl // average: averageDiGammas += MathsUtils.digamma(n_xz) + MathsUtils.digamma(n_yz) - MathsUtils.digamma(n_z); + if (debug) { + averageInverseCountInJointYZ += 1.0 / (double) n_yz; + } } averageDiGammas /= (double) N; - averageInverseCountInJointYZ /= (double) N; if (debug) { avNxz /= (double)N; avNyz /= (double)N; avNz /= (double)N; + averageInverseCountInJointYZ /= (double) N; System.out.printf("Average n_xz=%.3f, Average n_yz=%.3f, Average n_z=%.3f\n", avNxz, avNyz, avNz); } - condMi = MathsUtils.digamma(k) - 1.0/(double)k + - averageInverseCountInJointYZ - averageDiGammas; + condMi = MathsUtils.digamma(k) - 1.0/(double)k - averageDiGammas; miComputed = true; return condMi; } diff --git a/java/source/infodynamics/utils/MatrixUtils.java b/java/source/infodynamics/utils/MatrixUtils.java index a820795..92b2aa7 100755 --- a/java/source/infodynamics/utils/MatrixUtils.java +++ b/java/source/infodynamics/utils/MatrixUtils.java @@ -2389,6 +2389,33 @@ public class MatrixUtils { return c / (double) (x.length - 1); // -1 for sample covariance } + /** + *

Returns the covariance between the two arrays of data, with + * a given lag between the first and second.

+ *

See - Mathworld + *

+ * + * @param x time series 1 + * @param y time series 2 + * @param delay delay >= 0 to compute the covariance across (from first to second time series) + * @return the covariance + */ + public static double covariance(double[] x, double[] y, int delay) { + double meanX = 0, meanY = 0; + // No error checking if y is same length as x + for (int n = 0; n < x.length - delay; n++) { + meanX += x[n]; + meanY += y[n + delay]; + } + meanX /= (double) (x.length - delay); + meanY /= (double) (x.length - delay); + double c = 0; + for (int t = 0; t < x.length - delay; t++) { + c += (x[t] - meanX)*(y[t + delay]-meanY); + } + return c / (double) (x.length - delay - 1); // -1 for sample covariance + } + /** *

Returns the covariance between the first two columns of data.

* @@ -3093,6 +3120,7 @@ public class MatrixUtils { * @see {@link http://mathworld.wolfram.com/CholeskyDecomposition.html} * @see {@link http://en.wikipedia.org/wiki/Positive-definite_matrix} * @see {@link http://math.nist.gov/javanumerics/jama/} + * @see {@link http://www2.gsu.edu/~mkteer/npdmatri.html} */ public static double[][] CholeskyDecomposition(double[][] A) throws Exception { int n = A.length; @@ -3119,7 +3147,11 @@ public class MatrixUtils { d = A[j][j] - d; // Check the positive definite condition: if (d <= 0.0) { - throw new Exception("CholeskyDecomposition is only performed on positive-definite matrices"); + // Throw an error with some suggestions. The last suggestion is from my observations + // from a simple test with Matlab - I should find a reference for this ... + throw new Exception("CholeskyDecomposition is only performed on positive-definite matrices. " + + "Some reasons for non-positive-definite matrix are listed at http://www2.gsu.edu/~mkteer/npdmatri.html - " + + "note: a correlation matrix is non-positive-definite if you have more variables than observations"); } L[j][j] = Math.sqrt(d); // Set the upper triangular part to all zeros: diff --git a/java/source/infodynamics/utils/OctaveFileReader.java b/java/source/infodynamics/utils/OctaveFileReader.java index 2f1ebcd..7eeb464 100755 --- a/java/source/infodynamics/utils/OctaveFileReader.java +++ b/java/source/infodynamics/utils/OctaveFileReader.java @@ -25,6 +25,7 @@ public class OctaveFileReader { static final String GLOBAL_SCALAR_HEADER = "# type: scalar"; static final String MATRIX_HEADER = "# type: matrix"; static final String GLOBAL_MATRIX_HEADER = "# type: global matrix"; + static final String BOOL_MATRIX_HEADER = "# type: bool matrix"; static final String OCTAVE_DELIMITER = "[ \t]"; public OctaveFileReader(String octaveFilename) throws FileNotFoundException { @@ -123,7 +124,8 @@ public class OctaveFileReader { // We have the header line for this stored variable // Make sure the next line indicates a matrix line = br.readLine(); - if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER)) { + if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER) + && !line.equals(BOOL_MATRIX_HEADER)) { // The variable is not a matrix - it is not compatible with the // int [] type throw new Exception("Stored Octave variable " + name + " is not of a type compatible with int[] - not a matrix"); @@ -180,7 +182,8 @@ public class OctaveFileReader { // We have the header line for this stored variable // Make sure the next line indicates a matrix line = br.readLine(); - if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER)) { + if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER) + && !line.equals(BOOL_MATRIX_HEADER)) { // The variable is not a matrix - it is not compatible with the // int [] type throw new Exception("Stored Octave variable " + name + " is not of a type compatible with long[] - not a matrix"); @@ -294,7 +297,8 @@ public class OctaveFileReader { // We have the header line for this stored variable // Make sure the next line indicates a matirx line = br.readLine(); - if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER)) { + if (!line.equals(MATRIX_HEADER) && !line.equals(GLOBAL_MATRIX_HEADER) + && !line.equals(BOOL_MATRIX_HEADER)) { // The variable is not a matrix - it is not compatible with the // int [][] type throw new Exception("Stored Octave variable " + name + " is not of a type compatible with int[][]");