Adding common parent implementation of Conditional MI calculator from continuous to discrete variable, plus a linear-Gaussian implementation of it

This commit is contained in:
joseph.lizier 2013-03-18 06:08:12 +00:00
parent b17ba9b05c
commit 75c36dc13a
2 changed files with 456 additions and 0 deletions

View File

@ -0,0 +1,237 @@
package infodynamics.measures.continuous;
import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource;
import infodynamics.utils.MatrixUtils;
import infodynamics.utils.EmpiricalMeasurementDistribution;
import infodynamics.utils.RandomGenerator;
/**
* <p>This is an abstract calculator to
* compute the Conditional Mutual Information I(X;D|Z)
* between a discrete variable D and a
* vector of continuous variables X,
* conditioned on another vector of continuous variables Z.
* It implements common methods for child classes to build on</p>
*
* @author Joseph Lizier
*/
public abstract class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon
implements ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource {
/**
* we compute distances to the kth neighbour
*/
protected double[][] continuousDataX;
protected double[][] conditionedDataZ;
protected int[] discreteData;
protected int[] counts;
protected int base;
protected int dimensionsContinuous;
protected int dimensionsConditional;
protected double[] meansX;
protected double[] stdsX;
protected double[] meansZ;
protected double[] stdsZ;
protected boolean debug;
protected double condMi;
protected boolean miComputed;
protected int totalObservations;
public static final String PROP_NORMALISE = "NORMALISE";
protected boolean normalise = true;
public ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon() {
}
/**
* Initialise the calculator.
*
* @param dimensions number of joint continuous variables
* @param base number of discrete states
* @param dimensionsCond the number of joint continuous variables
* to condition on
*/
public void initialise(int dimensions, int base, int dimensionsCond) {
condMi = 0.0;
miComputed = false;
continuousDataX = null;
discreteData = null;
this.base = base;
dimensionsContinuous = dimensions;
dimensionsConditional = dimensionsCond;
}
/**
* Sets properties for the calculator.
* Valid properties include:
* <ul>
* <li>{@link #PROP_NORMALISE} - whether to normalise the individual
* variables (true by default)</li>
* </ul>
*
* @param propertyName
* @param propertyValue
*/
public void setProperty(String propertyName, String propertyValue) {
if (propertyName.equalsIgnoreCase(PROP_NORMALISE)) {
normalise = Boolean.parseBoolean(propertyValue);
}
}
public void startAddObservations() {
throw new RuntimeException("Not implemented yet");
}
public void finaliseAddObservations() throws Exception {
meansX = MatrixUtils.means(continuousDataX);
stdsX = MatrixUtils.stdDevs(continuousDataX, meansX);
meansZ = MatrixUtils.means(conditionedDataZ);
stdsZ = MatrixUtils.stdDevs(conditionedDataZ, meansZ);
if (normalise) {
// Take a copy since we're going to normalise it
continuousDataX = MatrixUtils.normaliseIntoNewArray(continuousDataX);
conditionedDataZ = MatrixUtils.normaliseIntoNewArray(conditionedDataZ);
}
// count the discrete states:
counts = new int[base];
for (int t = 0; t < discreteData.length; t++) {
counts[discreteData[t]]++;
}
totalObservations = discreteData.length;
}
public void addObservations(double[][] continuousObservations,
int[] discreteObservations, double[][] conditionedObservations)
throws Exception {
throw new RuntimeException("Not implemented yet");
}
public void setObservations(double[][] continuousObservations,
int[] discreteObservations, double[][] conditionedObservations)
throws Exception {
if ((continuousObservations.length != discreteObservations.length) ||
(continuousObservations.length != conditionedObservations.length)) {
throw new Exception("Time steps for observations2 " +
discreteObservations.length + " does not match the length " +
"of observations1 " + continuousObservations.length +
" and of conditionedObservations " + conditionedObservations.length);
}
if (continuousObservations[0].length == 0) {
throw new Exception("Computing MI with a null set of data");
}
if (conditionedObservations[0].length == 0) {
throw new Exception("Computing MI with a null set of conditioned data");
}
continuousDataX = continuousObservations;
discreteData = discreteObservations;
conditionedDataZ = conditionedObservations;
finaliseAddObservations();
}
/**
* Compute the significance of the mutual information of the previously supplied observations.
* We destroy the p(x,d|z) correlations, while retaining the p(x|z) marginals, to check how
* significant this conditional mutual information actually was.
*
* This is in the spirit of Chavez et. al., "Statistical assessment of nonlinear causality:
* application to epileptic EEG signals", Journal of Neuroscience Methods 124 (2003) 113-128
* which was performed for Transfer entropy.
*
* @param reorderDiscreteVariable boolean for whether to reorder the discrete variable (true)
* or the continuous variable (false)
* @param numPermutationsToCheck
* @return the proportion of conditional MI scores from the distribution which have higher or equal MIs to ours.
*/
public synchronized EmpiricalMeasurementDistribution computeSignificance(
boolean reorderDiscreteVariable, int numPermutationsToCheck) throws Exception {
// Generate the re-ordered indices:
RandomGenerator rg = new RandomGenerator();
// Use continuousDataX length (all variables have same length) even though
// we may be randomising the other variable:
int[][] newOrderings = rg.generateDistinctRandomPerturbations(
continuousDataX.length, numPermutationsToCheck);
return computeSignificance(reorderDiscreteVariable, newOrderings);
}
/**
* <p>As per {@link #computeSignificance(boolean, int)} but supplies
* the re-orderings of the observations of the named variable.</p>
*
* <p>We provide a simple implementation which would be suitable for
* any child class, though the child class may prefer to make its
* own implementation to make class-specific optimisations.
* Child classes must implement {@link java.lang.Cloneable}
* for this method to be callable for them, and indeed implement
* the clone() method in a way that protects their structure
* from alteration by surrogate data being supplied to it.</p>
*
* @param reorderDiscreteVariable boolean for whether to reorder the discrete variable (true)
* or the continuous variable (false)
* @param newOrderings the specific new orderings to use
* @return the proportion of conditional MI scores from the distribution which have higher or equal MIs to ours.
*/
public EmpiricalMeasurementDistribution computeSignificance(
boolean reorderDiscreteVariable, int[][] newOrderings) throws Exception {
int numPermutationsToCheck = newOrderings.length;
if (!miComputed) {
computeAverageLocalOfObservations();
}
// Take a clone of the object to compute the MI of the surrogates:
// (this is a shallow copy, it doesn't make new copies of all
// the arrays - child classes should override this)
ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource miSurrogateCalculator =
(ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSource) this.clone();
double[] surrogateMeasurements = new double[numPermutationsToCheck];
// Now compute the MI for each set of shuffled data:
for (int i = 0; i < numPermutationsToCheck; i++) {
if (reorderDiscreteVariable) {
// Reorder the discrete variable
int[] shuffledData =
MatrixUtils.extractSelectedTimePoints(discreteData, newOrderings[i]);
// Perform new initialisations
miSurrogateCalculator.initialise(
dimensionsContinuous, base, dimensionsConditional);
// Set new observations
miSurrogateCalculator.setObservations(continuousDataX,
shuffledData, conditionedDataZ);
} else {
// Reorder the continuous variable
double[][] shuffledData =
MatrixUtils.extractSelectedTimePointsReusingArrays(
continuousDataX, newOrderings[i]);
// Perform new initialisations
miSurrogateCalculator.initialise(
dimensionsContinuous, base, dimensionsConditional);
// Set new observations
miSurrogateCalculator.setObservations(shuffledData,
discreteData, conditionedDataZ);
}
// Compute the MI
surrogateMeasurements[i] = miSurrogateCalculator.computeAverageLocalOfObservations();
if (debug){
System.out.println("New MI was " + surrogateMeasurements[i]);
}
}
return new EmpiricalMeasurementDistribution(surrogateMeasurements, condMi);
}
public void setDebug(boolean debug) {
this.debug = debug;
}
public double getLastAverage() {
return condMi;
}
public int getNumObservations() {
return totalObservations;
}
}

View File

@ -0,0 +1,219 @@
package infodynamics.measures.continuous.gaussian;
import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon;
import infodynamics.utils.MatrixUtils;
/**
* <p>Computes the differential conditional mutual information of a given multivariate set of
* observations with a discrete variable, conditioned on another multivariate set
* of observations,
* assuming that the probability distribution function for these continuous observations is
* a multivariate Gaussian distribution.</p>
*
* <p>This is done by examining the conditional probability distribution for
* the multivariate continuous variable C (given the discrete
* variable D) against the probability distribution for C, all conditioned
* on another continuous variable Z:
* MI(C;D|Z) := H(C|Z) - H(C|D,Z).</p>
*
* <p><b>CAVEAT EMPTOR</b>: The real question this type of calculation asks
* is to what extent does knowing the value of the discrete variable reduce
* variance in the continuous variable(s), given the other known continuous variable.
* Indeed, it may not to behave (I should test this--)
* as we would normally expect a mutual information calculation: if we add more
* continuous variables in, it may increase in spite of redundancy between these
* variables. TODO Further exploration should take place here ...</p>
*
* <p>
* Usage:
* <ol>
* <li>Construct {@link #ConditionalMutualnfoCalculatorMultiVariateWithDiscreteSourceGaussian()}</li>
* <li>{@link #initialise(int, int, int)}</li>
* <li>Set properties using {@link #setProperty(String, String)}</li>
* <li>Provide the observations to the calculator using:
* {@link #setObservations(double[][], int[], double[][])}, or
* a sequence of:
* {@link #startAddObservations()},
* multiple calls to {@link #addObservations(double[][], int[], double[][])}
* and then
* {@link #finaliseAddObservations()}.</li>
* <li>Compute the required information-theoretic results, primarily:
* {@link #computeAverageLocalOfObservations()} to return the average differential
* entropy based on the variance of
* the supplied observations; or other calls to compute
* local values or statistical significance.</li>
* </ol>
* </p>
*
* @see <a href="http://mathworld.wolfram.com/DifferentialEntropy.html">Differential entropy for Gaussian random variables at Mathworld</a>
* @see <a href="http://en.wikipedia.org/wiki/Differential_entropy">Differential entropy for Gaussian random variables at Wikipedia</a>
* @see <a href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">Multivariate normal distribution on Wikipedia</a>
* @author Joseph Lizier joseph.lizier_at_gmail.com
*
*/
public class ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian
extends ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon
implements Cloneable {
/**
* Entropy calculator applied to the whole set of conditional data Z
*/
protected EntropyCalculatorMultiVariateGaussian entCalcZ;
/**
* Entropy calculator applied to the whole set of continuous data C
* and conditional data Z
*
*/
protected EntropyCalculatorMultiVariateGaussian entCalcCZ;
/**
* Entropy calculators applied to the set of conditional data Z
* associated with each discrete value
*/
protected EntropyCalculatorMultiVariateGaussian[] entCalcZForEachDiscrete;
/**
* Entropy calculators applied to the set of continuous data C
* and conditional data Z
* associated with each discrete value
*/
protected EntropyCalculatorMultiVariateGaussian[] entCalcCZForEachDiscrete;
public ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian() {
super();
entCalcZ = new EntropyCalculatorMultiVariateGaussian();
entCalcZForEachDiscrete = null;
entCalcCZ = new EntropyCalculatorMultiVariateGaussian();
entCalcCZForEachDiscrete = null;
}
/* (non-Javadoc)
* @see infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon#initialise(int, int, int)
*/
@Override
public void initialise(int dimensions, int base, int dimensionsCond) {
super.initialise(dimensions, base, dimensionsCond);
entCalcZ.initialise(dimensionsCond);
entCalcCZ.initialise(dimensions + dimensionsCond);
entCalcZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base];
entCalcCZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base];
for (int b = 0; b < base; b++) {
entCalcZForEachDiscrete[b] = new EntropyCalculatorMultiVariateGaussian();
entCalcCZForEachDiscrete[b] = new EntropyCalculatorMultiVariateGaussian();
// If any properties relevant for these calculators were set in
// setProperty then we should set them here
entCalcZForEachDiscrete[b].initialise(dimensionsCond);
entCalcCZForEachDiscrete[b].initialise(dimensions + dimensionsCond);
}
}
/* (non-Javadoc)
* @see infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceCommon#finaliseAddObservations()
*/
@Override
public void finaliseAddObservations() throws Exception {
super.finaliseAddObservations();
// Set the complete set of observations:
// (this will pick up any errors in the dimensions of the continuous
// observations)
entCalcZ.setObservations(conditionedDataZ);
double[][] joinedCZ = MatrixUtils.appendColumns(continuousDataX, conditionedDataZ);
entCalcCZ.setObservations(joinedCZ);
// Set the observations corresponding to each discrete value:
int totalNumberOfSuppliedObservations = 0;
for (int b = 0; b < base; b++) {
// Extract the observations for when this base value occurs:
double[][] obsZForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition(
conditionedDataZ, discreteData, b);
double[][] obsCZForThisDiscValue = MatrixUtils.extractSelectedPointsMatchingCondition(
joinedCZ, discreteData, b);
// Set the observations for each discrete value:
entCalcZForEachDiscrete[b].setObservations(obsZForThisDiscValue);
entCalcCZForEachDiscrete[b].setObservations(obsCZForThisDiscValue);
totalNumberOfSuppliedObservations += obsZForThisDiscValue.length;
}
// Check that all of the supplied observations were extracted corresponding
// to one of the allowed discrete values
if (totalNumberOfSuppliedObservations != discreteData.length) {
throw new Exception("Some values in discreteObservations were not in the range 0..base-1");
}
}
public double computeAverageLocalOfObservations() throws Exception {
// The average mutual information can be expressed
// as a difference between the conditional entropy of the continuous observations
// and the conditional entropy of the continuous given the
// discrete observations, both given the conditional continuous observations:
// I(C;D|Z) = H(C|Z) - H(C|D,Z)
// = H(C,Z) - H(Z) - H(C,Z|D) + H(C|D)
double meanConditionalEntropyCZ = 0;
double meanConditionalEntropyZ = 0;
for (int b = 0; b < base; b++) {
double pOfB = (double) counts[b] / (double) totalObservations;
meanConditionalEntropyZ += pOfB *
entCalcZForEachDiscrete[b].computeAverageLocalOfObservations();
meanConditionalEntropyCZ += pOfB *
entCalcCZForEachDiscrete[b].computeAverageLocalOfObservations();
}
double entCZ = entCalcCZ.computeAverageLocalOfObservations();
double entZ = entCalcZ.computeAverageLocalOfObservations();
condMi = entCZ - entZ
- meanConditionalEntropyCZ
+ meanConditionalEntropyZ;
if (debug) {
System.out.printf("H(C,Z)=%.4f - H(Z)=%.4f - H(C,Z|D)=%.4f + H(Z|D)=%.4f = %.4f\n",
entCZ, entZ, meanConditionalEntropyCZ, meanConditionalEntropyZ, condMi);
}
return condMi;
}
public double[] computeLocalOfPreviousObservations() throws Exception {
throw new RuntimeException("Not implemented yet");
}
public double[] computeLocalUsingPreviousObservations(
double[][] contStates, int[] discreteStates,
double[][] conditionedStates) throws Exception {
throw new RuntimeException("Not implemented yet");
}
/**
* Clone the object - note: while it does create new cloned instances of
* the {@link EntropyCalculatorMultiVariateGaussian} objects,
* I think these only
* have shallow copies to the data.
* This is enough though to maintain the structure across
* various {@link #computeSignificance(boolean, int)} calls.
*
* @see java.lang.Object#clone()
*/
@Override
protected Object clone() throws CloneNotSupportedException {
ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian theClone =
(ConditionalMutualInfoCalculatorMultiVariateWithDiscreteSourceGaussian) super.clone();
// Now assign clones of the EntropyCalculatorMultiVariateGaussian objects:
// First clone those for the conditional variable:
theClone.entCalcZ =
(EntropyCalculatorMultiVariateGaussian) entCalcZ.clone();
if (entCalcZForEachDiscrete != null) {
theClone.entCalcZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base];
for (int b = 0; b < base; b++) {
theClone.entCalcZForEachDiscrete[b] =
(EntropyCalculatorMultiVariateGaussian)
entCalcZForEachDiscrete[b].clone();
}
}
// Next clone those for the conditional variable and the continuous source variable
theClone.entCalcCZ =
(EntropyCalculatorMultiVariateGaussian) entCalcCZ.clone();
if (entCalcCZForEachDiscrete != null) {
theClone.entCalcCZForEachDiscrete = new EntropyCalculatorMultiVariateGaussian[base];
for (int b = 0; b < base; b++) {
theClone.entCalcCZForEachDiscrete[b] =
(EntropyCalculatorMultiVariateGaussian)
entCalcCZForEachDiscrete[b].clone();
}
}
return theClone;
}
}