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

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.

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);
 	}
 
 	/**

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

java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateKraskov2.java

@@ -8,14 +8,16 @@ import infodynamics.utils.MatrixUtils;
  * <p>Compute the Conditional Mutual Info using the Kraskov estimation method,
  * as extended by Frenzel and Pompe.</p>
  * <p>
- * 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.</p>
+ *  is perhaps recommended instead of this algorithm 2.</p>
  *  
  * <p>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("<n_xz>=%.3f, <n_yz>=%.3f, <n_z>=%.3f\n",
 					avNxz, avNyz, avNz);
-			System.out.printf("Av = digamma(k)=%.3f + <digammas>=%.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 + <digammas>=%.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("<n_xz>=%.3f, <n_yz>=%.3f, <n_z>=%.3f\n",
 					avNxz, avNyz, avNz);
-			System.out.printf("Av = digamma(k)=%.3f + <digammas>=%.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 + <digammas>=%.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("<n_xz>=%.3f, <n_yz>=%.3f, <n_z>=%.3f\n",
 					avNxz, avNyz, avNz);
-			System.out.printf("Av = digamma(k)=%.3f + <digammas>=%.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 + <digammas>=%.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("<n_xz>=%.3f, <n_yz>=%.3f, <n_z>=%.3f\n",
 					avNxz, avNyz, avNz);
-			System.out.printf("Av = digamma(k)=%.3f + <digammas>=%.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 + <digammas>=%.3f - 1/k=%.3f  = %.3f (<1/n_yz>=%.3f, <1/n_xz>=%.3f)\n",
+					digammaK, averageDiGammas, invK,
+					condMi, averageInverseCountInJointYZ, averageInverseCountInJointXZ);
 		}
 		
 		return localCondMi;

java/source/infodynamics/measures/continuous/kraskov/ConditionalMutualInfoCalculatorMultiVariateWithDiscreteKraskov.java

@@ -14,6 +14,9 @@ import infodynamics.utils.RandomGenerator;
  * <p>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.</p>
+ * <p>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}.
+ * </p>
  * 
  * @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;
 	}

java/source/infodynamics/utils/MatrixUtils.java

@@ -2389,6 +2389,33 @@ public class MatrixUtils {
 		return c / (double) (x.length - 1); // -1 for sample covariance
 	}
 
+	/**
+	 * <p>Returns the covariance between the two arrays of data, with
+	 * a given lag between the first and second.</p>
+	 * <p>See - <a href="http://mathworld.wolfram.com/Covariance.html">Mathworld</a>
+	 * </p>
+	 * 
+	 * @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
+	}
+
 	/**
 	 * <p>Returns the covariance between the first two columns of data.</p>
 	 * 
@@ -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:

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[][]");