mirror of https://github.com/jlizier/jidt
Added TransferEntropyCalculatorViaCondMutualInfo to allow calculation of TE via any implementation of a conditional MI calculator. Added concrete implementations using this for Gaussian estimators (Gaussian TE estimator is equivalent to Granger causality). Replaced existing Kraskov TE calculator (which used two MI calculators) with one which uses one conditional MI estimator via this new method. Also patched TransferEntropyCommon.computeStartAndEndTimePairs
This commit is contained in:
joseph.lizier
parent
c5402e8c48
commit
fd9bf116f9
|
|
@ -22,7 +22,10 @@ public interface ActiveInfoStorageCalculator {
|
|||
*/
|
||||
public static final String K_PROP_NAME = "k_HISTORY";
|
||||
/**
|
||||
* Embedding delay for the past history vector
|
||||
* Embedding delay for the past history vector.
|
||||
* The delay exists between points in the past k-vector
|
||||
* but there is always only one time step between the
|
||||
* past k-vector and the next observation.
|
||||
*/
|
||||
public static final String TAU_PROP_NAME = "TAU";
|
||||
|
||||
|
|
@ -39,6 +42,14 @@ public interface ActiveInfoStorageCalculator {
|
|||
*/
|
||||
public void initialise(int k) throws Exception;
|
||||
|
||||
/**
|
||||
* Initialise the calculator
|
||||
*
|
||||
* @param k Length of past history to consider
|
||||
* @param tau embedding delay to consider
|
||||
*/
|
||||
public void initialise(int k, int tau) throws Exception;
|
||||
|
||||
/**
|
||||
* Allows the user to set properties for the underlying calculator implementation
|
||||
*
|
||||
|
|
@ -86,7 +97,10 @@ public interface ActiveInfoStorageCalculator {
|
|||
* Sets the observations to compute the PDFs from.
|
||||
* Cannot be called in conjunction with start/add/finaliseAddObservations.
|
||||
* valid is a time series (with time indices the same as destination)
|
||||
* indicating whether the observation at that point is valid.
|
||||
* indicating whether the observation at that point is valid;
|
||||
* we only take tuples to add to the observation set where
|
||||
* all points in the time series (even between points in
|
||||
* the embedded k-vector with embedding delays) are valid.
|
||||
*
|
||||
* @param source observations for the source variable
|
||||
* @param destValid
|
||||
|
|
|
|||
|
|
@ -260,7 +260,7 @@ public class ActiveInfoStorageCalculatorViaMutualInfo implements
|
|||
*/
|
||||
public double[] computeLocalOfPreviousObservations() throws Exception {
|
||||
double[] local = miCalc.computeLocalOfPreviousObservations();
|
||||
if (miCalc.getAddedMoreThanOneObservationSet()) {
|
||||
if (!miCalc.getAddedMoreThanOneObservationSet()) {
|
||||
double[] localsToReturn = new double[local.length + (k-1)*tau + 1];
|
||||
System.arraycopy(local, 0, localsToReturn, (k-1)*tau + 1, local.length);
|
||||
return localsToReturn;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,474 @@
|
|||
package infodynamics.measures.continuous;
|
||||
|
||||
import infodynamics.utils.EmpiricalMeasurementDistribution;
|
||||
import infodynamics.utils.MatrixUtils;
|
||||
|
||||
import java.util.Vector;
|
||||
|
||||
|
||||
/**
|
||||
*
|
||||
* <p>Transfer entropy calculator which is implemented using a
|
||||
* Conditional Mutual Information calculator.
|
||||
* </p>
|
||||
*
|
||||
* TODO Delete TransferEntropyCalculatorCommon once we've switched everything over to use this?
|
||||
* Might be useful to leave it after all, and move common functionality from here to there.
|
||||
*
|
||||
* TODO Switch the properties like l etc into TransferEntropyCalculator
|
||||
* once all TE calculators follow this class.
|
||||
*
|
||||
* @see "Schreiber, Physical Review Letters 85 (2) pp.461-464, 2000;
|
||||
* <a href='http://dx.doi.org/10.1103/PhysRevLett.85.461'>download</a>
|
||||
* (for definition of transfer entropy)"
|
||||
* @see "Lizier, Prokopenko and Zomaya, Physical Review E 77, 026110, 2008;
|
||||
* <a href='http://dx.doi.org/10.1103/PhysRevE.77.026110'>download</a>
|
||||
* (for definition of <i>local</i> transfer entropy and qualification
|
||||
* of naming it as <i>apparent</i> transfer entropy)"
|
||||
*
|
||||
* @author Joseph Lizier, <a href="joseph.lizier at gmail.com">email</a>,
|
||||
* <a href="http://lizier.me/joseph/">www</a>
|
||||
*/
|
||||
public abstract class TransferEntropyCalculatorViaCondMutualInfo implements
|
||||
TransferEntropyCalculator {
|
||||
|
||||
/**
|
||||
* Underlying conditional mutual information calculator
|
||||
*/
|
||||
protected ConditionalMutualInfoCalculatorMultiVariate condMiCalc;
|
||||
/**
|
||||
* Length of past destination history to consider (embedding length)
|
||||
*/
|
||||
protected int k = 1;
|
||||
/**
|
||||
* Embedding delay to use between elements of the destination embeding vector.
|
||||
* We're hard-coding a delay of 1 between the history vector and the next
|
||||
* observation however.
|
||||
*/
|
||||
protected int k_tau = 1;
|
||||
/**
|
||||
* Length of past source history to consider (embedding length)
|
||||
*/
|
||||
protected int l = 1;
|
||||
/**
|
||||
* Embedding delay to use between elements of the source embeding vector.
|
||||
*/
|
||||
protected int l_tau = 1;
|
||||
/**
|
||||
* Source-destination next observation delay
|
||||
*/
|
||||
protected int delay = 1;
|
||||
/**
|
||||
* Time index of the last point in the destination embedding of the first
|
||||
* (destination past, source past, destination next) tuple than can be
|
||||
* taken from any set of time-series observations.
|
||||
*/
|
||||
protected int startTimeForFirstDestEmbedding;
|
||||
|
||||
protected boolean debug = false;
|
||||
|
||||
// TODO Move these properties up the TransferEntropyCalculator once all the calcs use them
|
||||
|
||||
/**
|
||||
* Embedding delay for the destination past history vector
|
||||
*/
|
||||
public static final String K_TAU_PROP_NAME = "k_TAU";
|
||||
/**
|
||||
* Embedding delay for the destination past history vector
|
||||
*/
|
||||
public static final String L_PROP_NAME = "l_HISTORY";
|
||||
/**
|
||||
* Embedding delay for the source past history vector
|
||||
*/
|
||||
public static final String L_TAU_PROP_NAME = "l_TAU";
|
||||
/**
|
||||
* Source-destination delay
|
||||
*/
|
||||
public static final String DELAY_PROP_NAME = "DELAY";
|
||||
|
||||
|
||||
public TransferEntropyCalculatorViaCondMutualInfo(String condMiCalculatorClassName)
|
||||
throws InstantiationException, IllegalAccessException, ClassNotFoundException {
|
||||
@SuppressWarnings("unchecked")
|
||||
Class<ConditionalMutualInfoCalculatorMultiVariate> condMiClass =
|
||||
(Class<ConditionalMutualInfoCalculatorMultiVariate>) Class.forName(condMiCalculatorClassName);
|
||||
ConditionalMutualInfoCalculatorMultiVariate condMiCalc = condMiClass.newInstance();
|
||||
construct(condMiCalc);
|
||||
}
|
||||
|
||||
public TransferEntropyCalculatorViaCondMutualInfo(Class<ConditionalMutualInfoCalculatorMultiVariate> condMiCalcClass)
|
||||
throws InstantiationException, IllegalAccessException, ClassNotFoundException {
|
||||
ConditionalMutualInfoCalculatorMultiVariate condMiCalc = condMiCalcClass.newInstance();
|
||||
construct(condMiCalc);
|
||||
}
|
||||
|
||||
/**
|
||||
* Construct this calculator by passing in a constructed but not initialised
|
||||
* underlying Conditional Mutual information calculator.
|
||||
*
|
||||
* @param condMiCalc
|
||||
*/
|
||||
public TransferEntropyCalculatorViaCondMutualInfo(ConditionalMutualInfoCalculatorMultiVariate condMiCalc) {
|
||||
construct(condMiCalc);
|
||||
}
|
||||
|
||||
protected void construct(ConditionalMutualInfoCalculatorMultiVariate condMiCalc) {
|
||||
this.condMiCalc = condMiCalc;
|
||||
}
|
||||
|
||||
public void initialise() throws Exception {
|
||||
initialise(k, k_tau, l, l_tau, delay);
|
||||
}
|
||||
|
||||
public void initialise(int k) throws Exception {
|
||||
initialise(k, k_tau, l, l_tau, delay);
|
||||
}
|
||||
|
||||
public void initialise(int k, int l) throws Exception {
|
||||
initialise(k, k_tau, l, l_tau, delay);
|
||||
}
|
||||
|
||||
public void initialise(int k, int l, int delay) throws Exception {
|
||||
initialise(k, k_tau, l, l_tau, delay);
|
||||
}
|
||||
|
||||
/**
|
||||
* Initialise the calculator
|
||||
*
|
||||
* @param k Length of destination past history to consider
|
||||
* @param k_tau embedding delay for the destination variable
|
||||
* @param l length of source past history to consider
|
||||
* @param l_tau embedding delay for the source variable
|
||||
* @param delay time lag between last element of source and destination next value
|
||||
*/
|
||||
public void initialise(int k, int k_tau, int l, int l_tau, int delay) throws Exception {
|
||||
if (delay < 0) {
|
||||
throw new Exception("Cannot compute TE with source-destination delay < 0");
|
||||
}
|
||||
this.k = k;
|
||||
this.k_tau = k_tau;
|
||||
this.l = l;
|
||||
this.l_tau = l_tau;
|
||||
this.delay = delay;
|
||||
|
||||
// Now check which point we can start taking observations from in any
|
||||
// addObservations call. These two integers represent the last
|
||||
// point of the destination embedding, in the cases where the destination
|
||||
// embedding itself determines where we can start taking observations, or
|
||||
// the case where the source embedding plus delay is longer and so determines
|
||||
// where we can start taking observations.
|
||||
int startTimeBasedOnDestPast = (k-1)*k_tau;
|
||||
int startTimeBasedOnSourcePast = (l-1)*l_tau + delay - 1;
|
||||
startTimeForFirstDestEmbedding = Math.max(startTimeBasedOnDestPast, startTimeBasedOnSourcePast);
|
||||
|
||||
condMiCalc.initialise(l, 1, k);
|
||||
}
|
||||
|
||||
/**
|
||||
* <p>Set the given property to the given value.
|
||||
* These can include:
|
||||
* <ul>
|
||||
* <li>{@link #K_PROP_NAME}</li>
|
||||
* <li>{@link #K_TAU_PROP_NAME}</li>
|
||||
* <li>{@link #L_TAU_PROP_NAME}</li>
|
||||
* <li>{@link #DELAY_PROP_NAME}</li>
|
||||
* </ul>
|
||||
* Or else it is assumed the property
|
||||
* is for the underlying {@link ConditionalMutualInfoCalculatorMultiVariate#setProperty(String, String)} implementation.
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @param propertyValue value of the property.
|
||||
* @throws Exception if there is a problem with the supplied value
|
||||
*/
|
||||
public void setProperty(String propertyName, String propertyValue) throws Exception {
|
||||
boolean propertySet = true;
|
||||
if (propertyName.equalsIgnoreCase(K_PROP_NAME)) {
|
||||
k = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(K_TAU_PROP_NAME)) {
|
||||
k_tau = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(L_PROP_NAME)) {
|
||||
l = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(L_TAU_PROP_NAME)) {
|
||||
l_tau = Integer.parseInt(propertyValue);
|
||||
} else if (propertyName.equalsIgnoreCase(DELAY_PROP_NAME)) {
|
||||
delay = Integer.parseInt(propertyValue);
|
||||
} else {
|
||||
// No property was set on this class, assume it is for the underlying
|
||||
// conditional MI calculator
|
||||
condMiCalc.setProperty(propertyName, propertyValue);
|
||||
propertySet = false;
|
||||
}
|
||||
if (debug && propertySet) {
|
||||
System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
|
||||
" to " + propertyValue);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* <p>Sets the single set of observations to compute the PDFs from.
|
||||
* Cannot be called in conjunction with
|
||||
* {@link #startAddObservations()}/{@link #addObservations(double[], double[])} /
|
||||
* {@link #finaliseAddObservations()}.</p>
|
||||
*
|
||||
* @param source observations for the source variable
|
||||
* @param destination observations for the destination variable
|
||||
* @throws Exception
|
||||
*/
|
||||
public void setObservations(double[] source, double[] destination) throws Exception {
|
||||
|
||||
if (source.length != destination.length) {
|
||||
throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
|
||||
source.length, destination.length));
|
||||
}
|
||||
if (source.length < startTimeForFirstDestEmbedding + 2) {
|
||||
// There are no observations to add here, the time series is too short
|
||||
throw new Exception("Not enough observations to set here given k, k_tau, l, l_tau and delay parameters");
|
||||
}
|
||||
double[][] currentDestPastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(destination, k, k_tau,
|
||||
startTimeForFirstDestEmbedding,
|
||||
destination.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] currentDestNextVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(destination, 1,
|
||||
startTimeForFirstDestEmbedding + 1,
|
||||
destination.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] currentSourcePastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(source, l, l_tau,
|
||||
startTimeForFirstDestEmbedding + 1 - delay,
|
||||
source.length - startTimeForFirstDestEmbedding - 1);
|
||||
condMiCalc.setObservations(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors);
|
||||
}
|
||||
|
||||
public void startAddObservations() {
|
||||
condMiCalc.startAddObservations();
|
||||
}
|
||||
|
||||
public void addObservations(double[] source, double[] destination) throws Exception {
|
||||
if (source.length != destination.length) {
|
||||
throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
|
||||
source.length, destination.length));
|
||||
}
|
||||
if (source.length < startTimeForFirstDestEmbedding + 2) {
|
||||
// There are no observations to add here, the time series is too short
|
||||
// Don't throw an exception, do nothing since more observations
|
||||
// can be added later.
|
||||
return;
|
||||
}
|
||||
double[][] currentDestPastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(destination, k, k_tau,
|
||||
startTimeForFirstDestEmbedding,
|
||||
destination.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] currentDestNextVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(destination, 1,
|
||||
startTimeForFirstDestEmbedding + 1,
|
||||
destination.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] currentSourcePastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(source, l, l_tau,
|
||||
startTimeForFirstDestEmbedding + 1 - delay,
|
||||
source.length - startTimeForFirstDestEmbedding - 1);
|
||||
condMiCalc.addObservations(currentSourcePastVectors, currentDestNextVectors, currentDestPastVectors);
|
||||
}
|
||||
|
||||
public void addObservations(double[] source, double[] destination,
|
||||
int startTime, int numTimeSteps) throws Exception {
|
||||
if (source.length != destination.length) {
|
||||
throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
|
||||
source.length, destination.length));
|
||||
}
|
||||
if (source.length < startTime + numTimeSteps) {
|
||||
// There are not enough observations given the arguments here
|
||||
throw new Exception("Not enough observations to set here given startTime and numTimeSteps parameters");
|
||||
}
|
||||
addObservations(MatrixUtils.select(source, startTime, numTimeSteps),
|
||||
MatrixUtils.select(destination, startTime, numTimeSteps));
|
||||
}
|
||||
|
||||
public void finaliseAddObservations() throws Exception {
|
||||
condMiCalc.finaliseAddObservations();
|
||||
}
|
||||
|
||||
public void setObservations(double[] source, double[] destination,
|
||||
boolean[] sourceValid, boolean[] destValid) throws Exception {
|
||||
|
||||
Vector<int[]> startAndEndTimePairs = computeStartAndEndTimePairs(sourceValid, destValid);
|
||||
|
||||
// We've found the set of start and end times for this pair
|
||||
startAddObservations();
|
||||
for (int[] timePair : startAndEndTimePairs) {
|
||||
int startTime = timePair[0];
|
||||
int endTime = timePair[1];
|
||||
addObservations(source, destination, startTime, endTime - startTime + 1);
|
||||
}
|
||||
finaliseAddObservations();
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute a vector of start and end pairs of time points, between which we have
|
||||
* valid series of both source and destinations.
|
||||
*
|
||||
* Made public so it can be used if one wants to compute the number of
|
||||
* observations prior to setting the observations.
|
||||
*
|
||||
* @param sourceValid
|
||||
* @param destValid
|
||||
* @return
|
||||
* @throws Exception
|
||||
*/
|
||||
public Vector<int[]> computeStartAndEndTimePairs(boolean[] sourceValid, boolean[] destValid) throws Exception {
|
||||
|
||||
if (sourceValid.length != destValid.length) {
|
||||
throw new Exception("Validity arrays must be of same length");
|
||||
}
|
||||
|
||||
int lengthOfDestPastRequired = (k-1)*k_tau + 1;
|
||||
int lengthOfSourcePastRequired = (l-1)*l_tau + 1;
|
||||
// int numSourcePointsBeforeDestStart = delay - 1 + lengthOfSourcePastRequired
|
||||
// - lengthOfDestPastRequired;
|
||||
|
||||
// Scan along the data avoiding invalid values
|
||||
int startTime = 0;
|
||||
Vector<int[]> startAndEndTimePairs = new Vector<int[]>();
|
||||
|
||||
// Simple solution -- this takes more complexity in time, but is
|
||||
// much faster to code:
|
||||
boolean previousWasOk = false;
|
||||
for (int t = startTimeForFirstDestEmbedding; t < destValid.length - 1; t++) {
|
||||
// Check the tuple with the history vector starting from
|
||||
// t and running backwards
|
||||
if (previousWasOk) {
|
||||
// Just check the very next values of each:
|
||||
if (destValid[t + 1] && sourceValid[t + 1 - delay]) {
|
||||
// We can continue adding to this sequence
|
||||
continue;
|
||||
} else {
|
||||
// We need to shut down this sequence now
|
||||
previousWasOk = false;
|
||||
int[] timePair = new int[2];
|
||||
timePair[0] = startTime;
|
||||
timePair[1] = t; // Previous time step was last valid one
|
||||
startAndEndTimePairs.add(timePair);
|
||||
continue;
|
||||
}
|
||||
}
|
||||
// Otherwise we're trying to start a new sequence, so check all values
|
||||
if (!destValid[t + 1]) {
|
||||
continue;
|
||||
}
|
||||
boolean allOk = true;
|
||||
for (int tBack = 0; tBack < lengthOfDestPastRequired; tBack++) {
|
||||
if (!destValid[t - tBack]) {
|
||||
allOk = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!allOk) {
|
||||
continue;
|
||||
}
|
||||
allOk = true;
|
||||
for (int tBack = delay - 1; tBack < delay - 1 + lengthOfSourcePastRequired; tBack++) {
|
||||
if (!sourceValid[t - tBack]) {
|
||||
allOk = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!allOk) {
|
||||
continue;
|
||||
}
|
||||
// Postcondition: We've got a first valid tuple:
|
||||
startTime = t - startTimeForFirstDestEmbedding;
|
||||
previousWasOk = true;
|
||||
}
|
||||
// Now check if we were running a sequence and terminate it:
|
||||
if (previousWasOk) {
|
||||
// We need to shut down this sequence now
|
||||
previousWasOk = false;
|
||||
int[] timePair = new int[2];
|
||||
timePair[0] = startTime;
|
||||
timePair[1] = destValid.length - 1;
|
||||
startAndEndTimePairs.add(timePair);
|
||||
}
|
||||
|
||||
return startAndEndTimePairs;
|
||||
}
|
||||
|
||||
public double computeAverageLocalOfObservations() throws Exception {
|
||||
return condMiCalc.computeAverageLocalOfObservations();
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns a time series of local TE values.
|
||||
* Pads the first (k-1)*tau + 1 elements with zeros (since AIS is undefined here)
|
||||
* if only one time series of observations was used.
|
||||
* Otherwise, local values for all separate series are concatenated, and without
|
||||
* padding of zeros at the start.
|
||||
*
|
||||
* @return an array of local TE values of the previously submitted observations.
|
||||
*/
|
||||
public double[] computeLocalOfPreviousObservations() throws Exception {
|
||||
double[] local = condMiCalc.computeLocalOfPreviousObservations();
|
||||
if (!condMiCalc.getAddedMoreThanOneObservationSet()) {
|
||||
double[] localsToReturn = new double[local.length + startTimeForFirstDestEmbedding + 1];
|
||||
System.arraycopy(local, 0, localsToReturn, startTimeForFirstDestEmbedding + 1, local.length);
|
||||
return localsToReturn;
|
||||
} else {
|
||||
return local;
|
||||
}
|
||||
}
|
||||
|
||||
public double[] computeLocalUsingPreviousObservations(double[] newSourceObservations, double[] newDestObservations) throws Exception {
|
||||
if (newSourceObservations.length != newDestObservations.length) {
|
||||
throw new Exception(String.format("Source and destination lengths (%d and %d) must match!",
|
||||
newSourceObservations.length, newDestObservations.length));
|
||||
}
|
||||
if (newDestObservations.length < startTimeForFirstDestEmbedding + 2) {
|
||||
// There are no observations to compute for here
|
||||
return new double[newDestObservations.length];
|
||||
}
|
||||
double[][] newDestPastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(newDestObservations, k, k_tau,
|
||||
startTimeForFirstDestEmbedding,
|
||||
newDestObservations.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] newDestNextVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(newDestObservations, 1,
|
||||
startTimeForFirstDestEmbedding + 1,
|
||||
newDestObservations.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[][] newSourcePastVectors =
|
||||
MatrixUtils.makeDelayEmbeddingVector(newSourceObservations, l, l_tau,
|
||||
startTimeForFirstDestEmbedding + 1 - delay,
|
||||
newSourceObservations.length - startTimeForFirstDestEmbedding - 1);
|
||||
double[] local = condMiCalc.computeLocalUsingPreviousObservations(
|
||||
newSourcePastVectors, newDestNextVectors, newDestPastVectors);
|
||||
// Pad the front of the array with zeros where local TE isn't defined:
|
||||
double[] localsToReturn = new double[local.length + startTimeForFirstDestEmbedding + 1];
|
||||
System.arraycopy(local, 0, localsToReturn, startTimeForFirstDestEmbedding + 1, local.length);
|
||||
return localsToReturn;
|
||||
|
||||
}
|
||||
|
||||
public EmpiricalMeasurementDistribution computeSignificance(
|
||||
int numPermutationsToCheck) throws Exception {
|
||||
return condMiCalc.computeSignificance(1, numPermutationsToCheck); // Reorder the source
|
||||
}
|
||||
|
||||
public EmpiricalMeasurementDistribution computeSignificance(
|
||||
int[][] newOrderings) throws Exception {
|
||||
return condMiCalc.computeSignificance(1, newOrderings); // Reorder the source
|
||||
}
|
||||
|
||||
public double getLastAverage() {
|
||||
return condMiCalc.getLastAverage();
|
||||
}
|
||||
|
||||
public int getNumObservations() throws Exception {
|
||||
return condMiCalc.getNumObservations();
|
||||
}
|
||||
|
||||
public boolean getAddedMoreThanOneObservationSet() {
|
||||
return condMiCalc.getAddedMoreThanOneObservationSet();
|
||||
}
|
||||
|
||||
public void setDebug(boolean debug) {
|
||||
this.debug = debug;
|
||||
condMiCalc.setDebug(debug);
|
||||
}
|
||||
}
|
||||
|
|
@ -273,7 +273,17 @@ public abstract class TransferEntropyCommon implements
|
|||
startAndEndTimePairs.add(timePair);
|
||||
// System.out.printf("t_s=%d, t_e=%d\n", startTime, endTime);
|
||||
lookingForStart = true;
|
||||
if (!destValid[t]) {
|
||||
// The current destination observation broke our chain;
|
||||
// so we need to start looking all over again:
|
||||
startTime = t + 1;
|
||||
} else {
|
||||
// The current source broke our chain (or we're at the
|
||||
// end of the series anyway, so this doesn't matter);
|
||||
// so we can keep the good destination history
|
||||
// that we've built up here:
|
||||
startTime = t - k + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,57 @@
|
|||
package infodynamics.measures.continuous.gaussian;
|
||||
|
||||
import infodynamics.measures.continuous.TransferEntropyCalculatorViaCondMutualInfo;
|
||||
|
||||
/**
|
||||
*
|
||||
* <p>
|
||||
* Implements a transfer entropy calculator using model of
|
||||
* Gaussian variables with linear interactions.
|
||||
* This is equivalent (up to a multiplicative constant) to
|
||||
* Granger causality (see Barnett et al., below).
|
||||
* This is achieved by plugging in {@link ConditionalMutualInfoCalculatorMultiVariateGaussian}
|
||||
* as the calculator into {@link TransferEntropyCalculatorViaCondMutualInfo}.
|
||||
* </p>
|
||||
*
|
||||
* <p>
|
||||
* Usage:
|
||||
* <ol>
|
||||
* <li>Construct: {@link #TransferEntropyCalculatorGaussian()}</li>
|
||||
* <li>Set properties: {@link #setProperty(String, String)} for each relevant property, including those
|
||||
* of either {@link TransferEntropyCalculatorViaCondMutualInfo#setProperty(String, String)}
|
||||
* or {@link ConditionalMutualInfoCalculatorMultiVariateGaussian#setProperty(String, String)}.</li>
|
||||
* <li>Initialise: by calling one of {@link #initialise()} etc.</li>
|
||||
* <li>Add observations to construct the PDFs: {@link #setObservations(double[])}, or [{@link #startAddObservations()},
|
||||
* {@link #addObservations(double[])}*, {@link #finaliseAddObservations()}]
|
||||
* Note: If not using setObservations(), the results from computeLocal
|
||||
* will be concatenated directly, and getSignificance will mix up observations
|
||||
* from separate trials (added in separate {@link #addObservations(double[])} calls.</li>
|
||||
* <li>Compute measures: e.g. {@link #computeAverageLocalOfObservations()} or
|
||||
* {@link #computeLocalOfPreviousObservations()} etc </li>
|
||||
* </ol>
|
||||
* </p>
|
||||
*
|
||||
* @author Joseph Lizier
|
||||
* @see "Lionel Barnett, Adam B. Barrett, Anil K. Seth, Physical Review Letters 103 (23) 238701, 2009;
|
||||
* <a href='http://dx.doi.org/10.1103/physrevlett.103.238701'>download</a>
|
||||
* (for direct relation between transfer entropy and Granger causality)"
|
||||
*
|
||||
* @see TransferEntropyCalculator
|
||||
*
|
||||
*/
|
||||
public class TransferEntropyCalculatorGaussian
|
||||
extends TransferEntropyCalculatorViaCondMutualInfo {
|
||||
|
||||
public static final String COND_MI_CALCULATOR_GAUSSIAN = ConditionalMutualInfoCalculatorMultiVariateGaussian.class.getName();
|
||||
|
||||
/**
|
||||
* Creates a new instance of the Gaussian-estimate style transfer entropy calculator
|
||||
* @throws ClassNotFoundException
|
||||
* @throws IllegalAccessException
|
||||
* @throws InstantiationException
|
||||
*
|
||||
*/
|
||||
public TransferEntropyCalculatorGaussian() throws InstantiationException, IllegalAccessException, ClassNotFoundException {
|
||||
super(COND_MI_CALCULATOR_GAUSSIAN);
|
||||
}
|
||||
}
|
||||
|
|
@ -36,7 +36,14 @@ public abstract class ConditionalMutualInfoCalculatorMultiVariateKraskov
|
|||
public static boolean tryKeepAllPairsNorms = true;
|
||||
public static int MAX_DATA_SIZE_FOR_KEEP_ALL_PAIRS_NORM = 2000;
|
||||
|
||||
/**
|
||||
* Property name for the number of nearest neighbours k to use in the Kraskov algorithm in
|
||||
* the full joint space.
|
||||
*/
|
||||
public final static String PROP_K = "k";
|
||||
/**
|
||||
* Normalisation to apply to the marginal spaces.
|
||||
*/
|
||||
public final static String PROP_NORM_TYPE = "NORM_TYPE";
|
||||
|
||||
public ConditionalMutualInfoCalculatorMultiVariateKraskov() {
|
||||
|
|
|
|||
|
|
@ -1,36 +1,171 @@
|
|||
package infodynamics.measures.continuous.kraskov;
|
||||
|
||||
import java.util.Hashtable;
|
||||
|
||||
import infodynamics.measures.continuous.ConditionalMutualInfoCalculatorMultiVariate;
|
||||
import infodynamics.measures.continuous.TransferEntropyCalculator;
|
||||
import infodynamics.measures.continuous.TransferEntropyCalculatorViaCondMutualInfo;
|
||||
import infodynamics.measures.continuous.gaussian.ConditionalMutualInfoCalculatorMultiVariateGaussian;
|
||||
|
||||
/**
|
||||
* <p>Compute the Transfer Entropy using the Kraskov estimation method for underlying
|
||||
* mutual information calculation, using direct calculation of the mutual information
|
||||
* rather than breaking it down into component multi-informations.<p>
|
||||
*
|
||||
* @see TransferEntropyCalculatorKraskovByMulti
|
||||
* @see MutualInfoCalculatorMultiVariateDirect
|
||||
* @see MutualInfoCalculatorMultiVariateDirect1
|
||||
* @see MutualInfoCalculatorMultiVariateDirect2
|
||||
* @author Joseph Lizier joseph.lizier at gmail.com
|
||||
* <p>
|
||||
* Implements a transfer entropy calculator using a conditional MI calculator
|
||||
* implementing the Kraskov-Grassberger estimator.
|
||||
* This is achieved by plugging in a {@link ConditionalMutualInfoCalculatorMultiVariateKraskov}
|
||||
* as the calculator into {@link TransferEntropyCalculatorViaCondMutualInfo}.
|
||||
* </p>
|
||||
*
|
||||
* <p>
|
||||
* Usage:
|
||||
* <ol>
|
||||
* <li>Construct: {@link #TransferEntropyCalculatorKraskov()}</li>
|
||||
* <li>Set properties: {@link #setProperty(String, String)} for each relevant property, including those
|
||||
* of either {@link TransferEntropyCalculatorViaCondMutualInfo#setProperty(String, String)}
|
||||
* or {@link ConditionalMutualInfoCalculatorMultiVariateGaussian#setProperty(String, String)}.</li>
|
||||
* <li>Initialise: by calling one of {@link #initialise()} etc.</li>
|
||||
* <li>Add observations to construct the PDFs: {@link #setObservations(double[])}, or [{@link #startAddObservations()},
|
||||
* {@link #addObservations(double[])}*, {@link #finaliseAddObservations()}]
|
||||
* Note: If not using setObservations(), the results from computeLocal
|
||||
* will be concatenated directly, and getSignificance will mix up observations
|
||||
* from separate trials (added in separate {@link #addObservations(double[])} calls.</li>
|
||||
* <li>Compute measures: e.g. {@link #computeAverageLocalOfObservations()} or
|
||||
* {@link #computeLocalOfPreviousObservations()} etc </li>
|
||||
* </ol>
|
||||
* </p>
|
||||
*
|
||||
* @author Joseph Lizier
|
||||
* @see "Kraskov, A., Stoegbauer, H., Grassberger, P., Physical Review E 69, (2004) 066138;
|
||||
* <a href='http://dx.doi.org/10.1103/PhysRevE.69.066138'>download</a>
|
||||
* (for introduction of Kraskov-Grassberger method for MI)"
|
||||
* @see "G. Gomez-Herrero, W. Wu, K. Rutanen, M. C. Soriano, G. Pipa, and R. Vicente,
|
||||
* arXiv:1008.0539, 2010;
|
||||
* <a href='http://arxiv.org/abs/1008.0539'>download</a>
|
||||
* (for introduction of Kraskov-Grassberger technique to transfer entropy)"
|
||||
* @see TransferEntropyCalculator
|
||||
*
|
||||
*/
|
||||
public class TransferEntropyCalculatorKraskov
|
||||
extends TransferEntropyCalculatorKraskovByMulti {
|
||||
extends TransferEntropyCalculatorViaCondMutualInfo {
|
||||
|
||||
public static final String COND_MI_CALCULATOR_KRASKOV1 = ConditionalMutualInfoCalculatorMultiVariateKraskov1.class.getName();
|
||||
public static final String COND_MI_CALCULATOR_KRASKOV2 = ConditionalMutualInfoCalculatorMultiVariateKraskov2.class.getName();
|
||||
|
||||
/**
|
||||
* Create the underlying Kraskov MI calculators using direct multi-variate MI calculators
|
||||
* Property for setting which underlying Kraskov-Grassberger algorithm to use.
|
||||
* Will only be applied at the next initialisation.
|
||||
*/
|
||||
public final static String PROP_KRASKOV_ALG_NUM = "ALG_NUM";
|
||||
|
||||
protected int kraskovAlgorithmNumber = 1;
|
||||
protected boolean algChanged = false;
|
||||
/**
|
||||
* Storage for the properties ready to pass onto the underlying conditional MI calculators should they change
|
||||
*/
|
||||
protected Hashtable<String,String> props;
|
||||
|
||||
/**
|
||||
* Creates a new instance of the Kraskov-estimate style transfer entropy calculator
|
||||
*
|
||||
* Uses algorithm 1 by default, as per Gomez-Herro et al.
|
||||
*
|
||||
* @throws ClassNotFoundException
|
||||
* @throws IllegalAccessException
|
||||
* @throws InstantiationException
|
||||
*
|
||||
*/
|
||||
protected void createKraskovMiCalculators() {
|
||||
if (kraskovAlgorithmNumber == 1) {
|
||||
mickPastToSource = new MutualInfoCalculatorMultiVariateKraskov1();
|
||||
mickNextPastToSource = new MutualInfoCalculatorMultiVariateKraskov1();
|
||||
} else {
|
||||
// Algorithm 2
|
||||
mickPastToSource = new MutualInfoCalculatorMultiVariateKraskov2();
|
||||
mickNextPastToSource = new MutualInfoCalculatorMultiVariateKraskov2();
|
||||
}
|
||||
public TransferEntropyCalculatorKraskov() throws InstantiationException, IllegalAccessException, ClassNotFoundException {
|
||||
super(COND_MI_CALCULATOR_KRASKOV1);
|
||||
kraskovAlgorithmNumber = 1;
|
||||
props = new Hashtable<String,String>();
|
||||
}
|
||||
|
||||
protected void shareDataBetweenUnderlyingCalculators() {
|
||||
// Nothing to do here
|
||||
/**
|
||||
* Creates a new instance of the Kraskov-Grassberger style transfer entropy calculator,
|
||||
* with the supplied conditional MI calculator name
|
||||
*
|
||||
* @param calculatorName fully qualified name of the underlying MI class.
|
||||
* Must be {@link #COND_MI_CALCULATOR_KRASKOV1} or {@link #COND_MI_CALCULATOR_KRASKOV2}
|
||||
* @throws ClassNotFoundException
|
||||
* @throws IllegalAccessException
|
||||
* @throws InstantiationException
|
||||
*
|
||||
*/
|
||||
public TransferEntropyCalculatorKraskov(String calculatorName) throws InstantiationException, IllegalAccessException, ClassNotFoundException {
|
||||
super(calculatorName);
|
||||
// Now check that it was one of our Kraskov-Grassberger calculators:
|
||||
if (calculatorName.equalsIgnoreCase(COND_MI_CALCULATOR_KRASKOV1)) {
|
||||
kraskovAlgorithmNumber = 1;
|
||||
} else if (calculatorName.equalsIgnoreCase(COND_MI_CALCULATOR_KRASKOV2)) {
|
||||
kraskovAlgorithmNumber = 2;
|
||||
} else {
|
||||
throw new ClassNotFoundException("Must be an underlying Kraskov-Grassberger conditional MI calculator");
|
||||
}
|
||||
props = new Hashtable<String,String>();
|
||||
}
|
||||
|
||||
/* (non-Javadoc)
|
||||
* @see infodynamics.measures.continuous.TransferEntropyCalculatorViaCondMutualInfo#initialise(int, int, int, int, int)
|
||||
*/
|
||||
@Override
|
||||
public void initialise(int k, int k_tau, int l, int l_tau, int delay)
|
||||
throws Exception {
|
||||
if (algChanged) {
|
||||
// The algorithm number was changed in a setProperties call:
|
||||
String newCalcName = COND_MI_CALCULATOR_KRASKOV1;
|
||||
if (kraskovAlgorithmNumber == 2) {
|
||||
newCalcName = COND_MI_CALCULATOR_KRASKOV2;
|
||||
}
|
||||
@SuppressWarnings("unchecked")
|
||||
Class<ConditionalMutualInfoCalculatorMultiVariate> condMiClass =
|
||||
(Class<ConditionalMutualInfoCalculatorMultiVariate>) Class.forName(newCalcName);
|
||||
ConditionalMutualInfoCalculatorMultiVariate newCondMiCalc = condMiClass.newInstance();
|
||||
construct(newCondMiCalc);
|
||||
// Set the properties for the Kraskov MI calculators (may pass in properties for our super class
|
||||
// as well, but they should be ignored)
|
||||
for (String key : props.keySet()) {
|
||||
newCondMiCalc.setProperty(key, props.get(key));
|
||||
}
|
||||
algChanged = false;
|
||||
}
|
||||
|
||||
super.initialise(k, k_tau, l, l_tau, delay);
|
||||
}
|
||||
|
||||
/**
|
||||
* Sets properties for the calculator.
|
||||
* Valid properties include:
|
||||
* <ul>
|
||||
* <li>{@link #PROP_KRASKOV_ALG_NUM} - which Kraskov algorithm number to use (1 or 2). Will only be applied at the next initialisation.</li>
|
||||
* <li>Any valid properties for {@link TransferEntropyCalculatorViaCondMutualInfo#setProperty(String, String)}</li>
|
||||
* <li>Any valid properties for {@link ConditionalMutualInfoCalculatorMultiVariateKraskov#setProperty(String, String)}</li>
|
||||
* </ul>
|
||||
* One should set {@link ConditionalMutualInfoCalculatorMultiVariateKraskov#PROP_K} here, the number
|
||||
* of neighbouring points one should count up to in determining the joint kernel size.
|
||||
*
|
||||
* @param propertyName name of the property
|
||||
* @param propertyValue value of the property (as a string)
|
||||
*/
|
||||
public void setProperty(String propertyName, String propertyValue)
|
||||
throws Exception {
|
||||
if (propertyName.equalsIgnoreCase(PROP_KRASKOV_ALG_NUM)) {
|
||||
int previousAlgNumber = kraskovAlgorithmNumber;
|
||||
kraskovAlgorithmNumber = Integer.parseInt(propertyValue);
|
||||
if ((kraskovAlgorithmNumber != 1) && (kraskovAlgorithmNumber != 2)) {
|
||||
throw new Exception("Kraskov algorithm number (" + kraskovAlgorithmNumber
|
||||
+ ") must be either 1 or 2");
|
||||
}
|
||||
if (kraskovAlgorithmNumber != previousAlgNumber) {
|
||||
algChanged = true;
|
||||
}
|
||||
if (debug) {
|
||||
System.out.println(this.getClass().getSimpleName() + ": Set property " + propertyName +
|
||||
" to " + propertyValue);
|
||||
}
|
||||
} else {
|
||||
// Assume it was a property for the parent class or underlying conditional MI calculator
|
||||
super.setProperty(propertyName, propertyValue);
|
||||
props.put(propertyName, propertyValue); // This will keep properties for the super class as well as the cond MI calculator, but this is ok
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue