Adding NORMALISE and NOISE_LEVEL_TO_ADD properties to all mutual information calculators for continuous-valued data (extends this capability from only KSG to Gaussian and kernel as well, the latter already had NORMALISE)

This commit is contained in:
Joseph Lizier 2019-04-02 12:38:17 +11:00
parent ecec4d20e3
commit b1437452e4
6 changed files with 176 additions and 110 deletions

View File

@ -79,7 +79,18 @@ public interface MutualInfoCalculatorMultiVariate
* must be >= 0)
*/
public static final String PROP_TIME_DIFF = "TIME_DIFF";
/**
* Property name for whether to normalise the incoming data to
* mean 0, standard deviation 1 (default true)
*/
public static final String PROP_NORMALISE = "NORMALISE";
/**
* Property name for the std deviation of random Gaussian noise to be
* added to the data (default is 1e-8, matching the MILCA toolkit);
* if the data is to be normalised, that will be done before adding this noise.
*/
public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
/**
* Compute the mutual information if the observations of the
* first variable (source)

View File

@ -24,6 +24,7 @@ import infodynamics.utils.MatrixUtils;
import infodynamics.utils.RandomGenerator;
import java.util.Iterator;
import java.util.Random;
import java.util.Vector;
/**
@ -60,6 +61,14 @@ public abstract class MutualInfoMultiVariateCommon implements
* {@link addObservations(double[][], double[][])} functions.
*/
protected double[][] sourceObservations;
/**
* Stored means of source variables before any normalising or noise addition
*/
protected double[] sourceMeansBeforeNorm;
/**
* Stored standard deviations of source variables before any normalising or noise addition
*/
protected double[] sourceStdsBeforeNorm;
/**
* The set of destination observations, retained in case the user wants to retrieve the local
@ -68,6 +77,14 @@ public abstract class MutualInfoMultiVariateCommon implements
* {@link addObservations(double[][], double[][])} functions.
*/
protected double[][] destObservations;
/**
* Stored means of destination variables before any normalising or noise addition
*/
protected double[] destMeansBeforeNorm;
/**
* Stored standard deviations of destination variablesbefore any normalising or noise addition
*/
protected double[] destStdsBeforeNorm;
/**
* Total number of observations supplied.
@ -114,6 +131,18 @@ public abstract class MutualInfoMultiVariateCommon implements
* Whether the user has supplied more than one (disjoint) set of samples
*/
protected boolean addedMoreThanOneObservationSet;
/**
* Whether to normalise the incoming data
*/
protected boolean normalise = true;
/**
* Whether to add an amount of random noise to the incoming data
*/
protected boolean addNoise = false;
/**
* Amount of random Gaussian noise to add to the incoming data
*/
protected double noiseLevel = (double) 0.0;
/* (non-Javadoc)
* @see infodynamics.measures.continuous.ChannelCalculatorCommon#initialise()
@ -130,6 +159,10 @@ public abstract class MutualInfoMultiVariateCommon implements
miComputed = false;
sourceObservations = null;
destObservations = null;
sourceMeansBeforeNorm = null;
sourceStdsBeforeNorm = null;
destMeansBeforeNorm = null;
destStdsBeforeNorm = null;
addedMoreThanOneObservationSet = false;
}
@ -145,6 +178,18 @@ public abstract class MutualInfoMultiVariateCommon implements
* Time difference between source and destination (0 by default).
* Must be >= 0.
* </li>
* <li>{@link #PROP_NORMALISE} -- whether to normalise the incoming individual
* variables to mean 0 and standard deviation 1 (true by default, except for Gaussian calculator)</li>
* <li>{@link #PROP_ADD_NOISE} -- a standard deviation for an amount of
* random Gaussian noise to add to
* each variable, to avoid having neighbourhoods with artificially
* large counts. (We also accept "false" to indicate "0".)
* The amount is added in after any normalisation,
* so can be considered as a number of standard deviations of the data.
* Default is 0 for most estimators; this is strongly recommended by
* by Kraskov for the KSG method though, so for that estimator we
* use 1e-8 to match the MILCA toolkit (though note it adds in
* a random amount of noise in [0,noiseLevel) ).</li>
* </ul>
*
* <p>Unknown property values are ignored.</p>
@ -156,16 +201,23 @@ public abstract class MutualInfoMultiVariateCommon implements
public void setProperty(String propertyName, String propertyValue)
throws Exception {
// TODO Have a NORMALISE property which is true by default,
// except for the linear Gaussian calculator (see conditional
// mutual info calculators)
boolean propertySet = true;
if (propertyName.equalsIgnoreCase(PROP_TIME_DIFF)) {
timeDiff = Integer.parseInt(propertyValue);
if (timeDiff < 0) {
throw new Exception("Time difference must be >= 0. Flip data1 and data2 around if required.");
}
} else if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
normalise = Boolean.parseBoolean(propertyValue);
} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
if (propertyValue.equals("0") ||
propertyValue.equalsIgnoreCase("false")) {
addNoise = false;
noiseLevel = 0;
} else {
addNoise = true;
noiseLevel = Double.parseDouble(propertyValue);
}
} else {
// No property was set here
propertySet = false;
@ -182,6 +234,10 @@ public abstract class MutualInfoMultiVariateCommon implements
if (propertyName.equalsIgnoreCase(PROP_TIME_DIFF)) {
return Integer.toString(timeDiff);
} else if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
return Boolean.toString(normalise);
} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
return Double.toString(noiseLevel);
} else {
// No property was recognised here
return null;
@ -356,6 +412,47 @@ public abstract class MutualInfoMultiVariateCommon implements
// We don't need to keep the vectors of observation sets anymore:
vectorOfSourceObservations = null;
vectorOfDestinationObservations = null;
// Normalise the data if required, and store means/stds any normalising
if (normalise) {
normaliseData();
} else {
sourceMeansBeforeNorm = MatrixUtils.means(sourceObservations);
sourceStdsBeforeNorm = MatrixUtils.stdDevs(sourceObservations, sourceMeansBeforeNorm);
destMeansBeforeNorm = MatrixUtils.means(destObservations);
destStdsBeforeNorm = MatrixUtils.stdDevs(destObservations, destMeansBeforeNorm);
}
// Add Gaussian noise of std dev noiseLevel to the data if required
if (addNoise) {
Random random = new Random();
for (int r = 0; r < sourceObservations.length; r++) {
for (int c = 0; c < dimensionsSource; c++) {
sourceObservations[r][c] +=
random.nextGaussian()*noiseLevel;
}
for (int c = 0; c < dimensionsDest; c++) {
destObservations[r][c] +=
random.nextGaussian()*noiseLevel;
}
}
}
}
/**
* Protected method to normalise the stored data samples for each variable.
* This method can be overriden by children if required to perform
* specific actions for their estimation methods.
*/
protected void normaliseData() {
// We can overwrite these since they're already
// a copy of the users' data.
double[][] stats = MatrixUtils.normalise(sourceObservations);
sourceMeansBeforeNorm = stats[0];
sourceStdsBeforeNorm = stats[1];
stats = MatrixUtils.normalise(destObservations);
destMeansBeforeNorm = stats[0];
destStdsBeforeNorm = stats[1];
}
/**

View File

@ -94,12 +94,6 @@ public class MutualInfoCalculatorMultiVariateGaussian
*/
protected double[][] Ldest;
/**
* Means of the most recently supplied observations (source variables
* listed first, destination variables second).
*/
protected double[] means;
/**
* Cached determinant of the joint covariance matrix
*/
@ -127,7 +121,7 @@ public class MutualInfoCalculatorMultiVariateGaussian
* Construct an instance of the Gaussian MI calculator
*/
public MutualInfoCalculatorMultiVariateGaussian() {
// Nothing to do
normalise = false; // Not much need to have this set for Gaussian, just creating work
}
public void initialise(int sourceDimensions, int destDimensions) {
@ -135,7 +129,6 @@ public class MutualInfoCalculatorMultiVariateGaussian
L = null;
Lsource = null;
Ldest = null;
means = null;
detCovariance = 0;
detSourceCovariance = 0;
detDestCovariance = 0;
@ -214,13 +207,6 @@ public class MutualInfoCalculatorMultiVariateGaussian
// destObservations[][] arrays.
super.finaliseAddObservations();
// Store the means of each variable (useful for local values later)
means = new double[dimensionsSource + dimensionsDest];
double[] sourceMeans = MatrixUtils.means(sourceObservations);
double[] destMeans = MatrixUtils.means(destObservations);
System.arraycopy(sourceMeans, 0, means, 0, dimensionsSource);
System.arraycopy(destMeans, 0, means, dimensionsSource, dimensionsDest);
// Store the covariances of the variables
// Generally, this should not throw an exception, since we checked
// the observations had the correct number of variables
@ -361,8 +347,19 @@ public class MutualInfoCalculatorMultiVariateGaussian
*/
public void setCovarianceAndMeans(double[][] covariance, double[] means,
int numObservations) throws Exception {
this.means = means;
sourceMeansBeforeNorm = MatrixUtils.select(means, 0, dimensionsSource);
destMeansBeforeNorm = MatrixUtils.select(means, dimensionsSource, dimensionsDest);
setCovariance(covariance, numObservations);
// set the std deviations from the covariance:
sourceStdsBeforeNorm = new double[dimensionsSource];
destStdsBeforeNorm = new double[dimensionsDest];
for (int i = 0; i < covariance.length; i++) {
if (i < dimensionsSource) {
sourceStdsBeforeNorm[i] = Math.sqrt(covariance[i][i]);
} else {
destStdsBeforeNorm[i - dimensionsSource] = Math.sqrt(covariance[i][i]);
}
}
}
/**
@ -589,10 +586,16 @@ public class MutualInfoCalculatorMultiVariateGaussian
protected double[] computeLocalUsingPreviousObservations(double[][] newSourceObs,
double[][] newDestObs, boolean isPreviousObservations) throws Exception {
if (means == null) {
if (sourceMeansBeforeNorm == null) {
throw new Exception("Cannot compute local values without having means either supplied or computed via setObservations()");
}
if ((!isPreviousObservations) && normalise) {
// Need to normalise new observations
newSourceObs = MatrixUtils.normaliseIntoNewArray(newSourceObs, sourceMeansBeforeNorm, sourceStdsBeforeNorm);
newDestObs = MatrixUtils.normaliseIntoNewArray(newDestObs, destMeansBeforeNorm, destStdsBeforeNorm);
}
// Check that the covariance matrix was positive definite:
// (this was done earlier in computing the Cholesky decomposition,
// we may still need to compute the determinant)
@ -635,11 +638,14 @@ public class MutualInfoCalculatorMultiVariateGaussian
double[][] invDestCovariance = MatrixUtils.solveViaCholeskyResult(Ldest,
MatrixUtils.identityMatrix(Ldest.length));
// Use the following array to index directly into dimensions of the destination sample vectors.
// We don't need a sourceIndicesSelected because there is no offset
// from zero for sourceIndicesInCovariance (unlike for destIndicesInCovariance)
int[] destIndicesSelected = MatrixUtils.subtract(destIndicesInCovariance, dimensionsSource);
// Now, only use the means from the subsets of linearly independent variables:
// double[] sourceMeans = MatrixUtils.select(means, 0, dimensionsSource);
double[] sourceMeans = MatrixUtils.select(means, sourceIndicesInCovariance);
// double[] destMeans = MatrixUtils.select(means, dimensionsSource, dimensionsDest);
double[] destMeans = MatrixUtils.select(means, destIndicesInCovariance);
double[] sourceMeans = MatrixUtils.select(sourceMeansBeforeNorm, sourceIndicesInCovariance);
double[] destMeans = MatrixUtils.select(destMeansBeforeNorm, destIndicesSelected);
int lengthOfReturnArray, offset;
if (isPreviousObservations && addedMoreThanOneObservationSet) {
@ -654,11 +660,6 @@ public class MutualInfoCalculatorMultiVariateGaussian
offset = timeDiff;
}
// Use the following array to index directly into dimensions of the destination sample vectors.
// We don't need a sourceIndicesSelected because there is no offset
// from zero for sourceIndicesInCovariance (unlike for destIndicesInCovariance)
int[] destIndicesSelected = MatrixUtils.subtract(destIndicesInCovariance, dimensionsSource);
// In case we need this for bias correction:
ChiSquareMeasurementDistribution analyticMeasDist = computeSignificance();

View File

@ -163,6 +163,13 @@ public class MutualInfoCalculatorMultiVariateKernel
}
}
@Override
protected void normaliseData() {
// Here we actually do nothing -- normalising is implemented for
// the kernel estimator by adjusting the kernel widths only, not by
// changing the data we operate on!!
}
/**
* Compute the average MI from the previously supplied observations.
*
@ -624,8 +631,6 @@ public class MutualInfoCalculatorMultiVariateKernel
* kernel width to be used in the calculation. If {@link #normalise} is set,
* then this is a number of standard deviations; otherwise it
* is an absolute value. Default is {@link #DEFAULT_KERNEL_WIDTH}.</li>
* <li>{@link #NORMALISE_PROP_NAME} -- whether to normalise the incoming variables
* to mean 0, standard deviation 1, or not (default false). Sets {@link #normalise}.</li>
* <li>{@link #DYN_CORR_EXCL_TIME_NAME} -- a dynamics exclusion time window (see Kantz and Schreiber),
* default is 0 which means no dynamic exclusion window.</li>
* <li>{@link #FORCE_KERNEL_COMPARE_TO_ALL} -- whether to force the underlying kernel estimators to compare
@ -655,6 +660,8 @@ public class MutualInfoCalculatorMultiVariateKernel
propertyName.equalsIgnoreCase(EPSILON_PROP_NAME)) {
kernelWidth = Double.parseDouble(propertyValue);
} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
// This is nominally a property of the super class, but we
// have extra work to do with this one
normalise = Boolean.parseBoolean(propertyValue);
mvkeSource.setNormalise(normalise);
mvkeDest.setNormalise(normalise);
@ -695,8 +702,6 @@ public class MutualInfoCalculatorMultiVariateKernel
if (propertyName.equalsIgnoreCase(KERNEL_WIDTH_PROP_NAME) ||
propertyName.equalsIgnoreCase(EPSILON_PROP_NAME)) {
return Double.toString(kernelWidth);
} else if (propertyName.equalsIgnoreCase(NORMALISE_PROP_NAME)) {
return Boolean.toString(normalise);
} else if (propertyName.equalsIgnoreCase(DYN_CORR_EXCL_TIME_NAME)) {
return Integer.toString(dynCorrExclTime);
} else if (propertyName.equalsIgnoreCase(FORCE_KERNEL_COMPARE_TO_ALL)) {

View File

@ -20,7 +20,6 @@ package infodynamics.measures.continuous.kraskov;
import java.util.Calendar;
import java.util.PriorityQueue;
import java.util.Random;
import infodynamics.measures.continuous.MutualInfoCalculatorMultiVariate;
import infodynamics.measures.continuous.MutualInfoMultiVariateCommon;
@ -91,16 +90,6 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
* default is {@link EuclideanUtils#NORM_MAX_NORM}.
*/
public final static String PROP_NORM_TYPE = "NORM_TYPE";
/**
* Property name for whether to normalise the incoming data to
* mean 0, standard deviation 1 (default true)
*/
public static final String PROP_NORMALISE = "NORMALISE";
/**
* Property name for an amount of random Gaussian noise to be
* added to the data (default is 1e-8, matching the MILCA toolkit).
*/
public static final String PROP_ADD_NOISE = "NOISE_LEVEL_TO_ADD";
/**
* Property name for a dynamics exclusion time window
* otherwise known as Theiler window (see Kantz and Schreiber).
@ -129,18 +118,6 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
*/
public static final String PROP_GPU_LIBRARY_PATH = "GPU_LIBRARY_PATH";
/**
* Whether to normalise the incoming data
*/
protected boolean normalise = true;
/**
* Whether to add an amount of random noise to the incoming data
*/
protected boolean addNoise = true;
/**
* Amount of random Gaussian noise to add to the incoming data
*/
protected double noiseLevel = (double) 1e-8;
/**
* Whether we use dynamic correlation exclusion
*/
@ -202,6 +179,9 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
*/
public MutualInfoCalculatorMultiVariateKraskov() {
super();
// Switch on adding noise to the data by default for the KSG estimator
addNoise = true;
noiseLevel = (double) 1e-8;
}
@Override
@ -226,20 +206,9 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
* working out the norms between the points in each marginal space.
* Options are defined by {@link KdTree#setNormType(String)} -
* default is {@link EuclideanUtils#NORM_MAX_NORM}.</li>
* <li>{@link #PROP_NORMALISE} -- whether to normalise the incoming individual
* variables to mean 0 and standard deviation 1 (true by default)</li>
* <li>{@link #PROP_DYN_CORR_EXCL_TIME} -- a dynamics exclusion time window,
* also known as Theiler window (see Kantz and Schreiber);
* default is 0 which means no dynamic exclusion window.</li>
* <li>{@link #PROP_ADD_NOISE} -- a standard deviation for an amount of
* random Gaussian noise to add to
* each variable, to avoid having neighbourhoods with artificially
* large counts. (We also accept "false" to indicate "0".)
* The amount is added in after any normalisation,
* so can be considered as a number of standard deviations of the data.
* (Recommended by Kraskov. MILCA uses 1e-8; but adds in
* a random amount of noise in [0,noiseLevel) ).
* Default 1e-8 to match the noise order in MILCA toolkit.</li>
* <li>{@link #PROP_NUM_THREADS} -- the integer number of parallel threads
* to use in the computation. Can be passed as a string "USE_ALL"
* to use all available processors on the machine.
@ -259,20 +228,9 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
k = Integer.parseInt(propertyValue);
} else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) {
normType = KdTree.validateNormType(propertyValue);
} else if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
normalise = Boolean.parseBoolean(propertyValue);
} else if (propertyName.equalsIgnoreCase(PROP_DYN_CORR_EXCL_TIME)) {
dynCorrExclTime = Integer.parseInt(propertyValue);
dynCorrExcl = (dynCorrExclTime > 0);
} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
if (propertyValue.equals("0") ||
propertyValue.equalsIgnoreCase("false")) {
addNoise = false;
noiseLevel = 0;
} else {
addNoise = true;
noiseLevel = Double.parseDouble(propertyValue);
}
} else if (propertyName.equalsIgnoreCase(PROP_NUM_THREADS)) {
if (propertyValue.equalsIgnoreCase(USE_ALL_THREADS)) {
numThreads = Runtime.getRuntime().availableProcessors();
@ -315,12 +273,8 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
return Integer.toString(k);
} else if (propertyName.equalsIgnoreCase(PROP_NORM_TYPE)) {
return KdTree.convertNormTypeToString(normType);
} else if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
return Boolean.toString(normalise);
} else if (propertyName.equalsIgnoreCase(PROP_DYN_CORR_EXCL_TIME)) {
return Integer.toString(dynCorrExclTime);
} else if (propertyName.equalsIgnoreCase(PROP_ADD_NOISE)) {
return Double.toString(noiseLevel);
} else if (propertyName.equalsIgnoreCase(PROP_NUM_THREADS)) {
return Integer.toString(numThreads);
} else if (propertyName.equalsIgnoreCase(PROP_USE_GPU)) {
@ -355,29 +309,6 @@ public abstract class MutualInfoCalculatorMultiVariateKraskov
k + ") and any dynamic correlation exclusion (" + dynCorrExclTime + ")");
}
// Normalise the data if required
if (normalise) {
// We can overwrite these since they're already
// a copy of the users' data.
MatrixUtils.normalise(sourceObservations);
MatrixUtils.normalise(destObservations);
}
if (addNoise) {
Random random = new Random();
// Add Gaussian noise of std dev noiseLevel to the data
for (int r = 0; r < sourceObservations.length; r++) {
for (int c = 0; c < dimensionsSource; c++) {
sourceObservations[r][c] +=
random.nextGaussian()*noiseLevel;
}
for (int c = 0; c < dimensionsDest; c++) {
destObservations[r][c] +=
random.nextGaussian()*noiseLevel;
}
}
}
// Set the constants:
digammaK = MathsUtils.digamma(k);
digammaN = MathsUtils.digamma(totalObservations);

View File

@ -2336,11 +2336,16 @@ public class MatrixUtils {
}
/**
* Normalises the elements along each column of the matrix
* Normalises the elements along each column of the matrix, storing the result
* back into the original matrix.
* If the standard deviation of the column is zero, we just remove the mean
*
* @param matrix 2D matrix of doubles
* @param matrix 2D matrix of doubles -- normalised results are returned in place (not
* in the return value of the method)
* @return a double[][] returnMatrix where returnMatrix[0] is an array of the means
* of each column, and returnMatrix[1] is an array of the std deviations of each column.
*/
public static void normalise(double[][] matrix) {
public static double[][] normalise(double[][] matrix) {
double[] means = means(matrix);
double[] stds = stdDevs(matrix, means);
@ -2357,6 +2362,8 @@ public class MatrixUtils {
} // else we just subtract off the mean
}
}
return new double[][] {means, stds};
}
/**
@ -4272,4 +4279,18 @@ public class MatrixUtils {
}
return array;
}
/**
* Return an array correpsonding to the diagonal of the given matrix
*
* @param matrix
* @return
*/
public static double[] diagonal(double[][] matrix) {
double[] diag = new double[matrix.length];
for (int i = 0; i < matrix.length; i++) {
diag[i] = matrix[i][i];
}
return diag;
}
}