mirror of https://github.com/jlizier/jidt
Added python examples 7 (ensemble method with TE Kraskov) and 9 (TE Kraskov with auto-embedding Ragwitz)
This commit is contained in:
jlizier
parent
f761a291e8
commit
bea2c77aa3
|
|
@ -0,0 +1,83 @@
|
|||
##
|
||||
## Java Information Dynamics Toolkit (JIDT)
|
||||
## Copyright (C) 2012, Joseph T. Lizier
|
||||
##
|
||||
## This program is free software: you can redistribute it and/or modify
|
||||
## it under the terms of the GNU General Public License as published by
|
||||
## the Free Software Foundation, either version 3 of the License, or
|
||||
## (at your option) any later version.
|
||||
##
|
||||
## This program is distributed in the hope that it will be useful,
|
||||
## but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
## GNU General Public License for more details.
|
||||
##
|
||||
## You should have received a copy of the GNU General Public License
|
||||
## along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
##
|
||||
|
||||
# = Example 7 - Ensemble method with transfer entropy on continuous data using Kraskov estimators =
|
||||
|
||||
# Calculation of transfer entropy (TE) by supplying an ensemble of samples from multiple time series.
|
||||
# We use continuous-valued data using the Kraskov-estimator TE calculator here.
|
||||
|
||||
from jpype import *
|
||||
import random
|
||||
import math
|
||||
|
||||
# Change location of jar to match yours:
|
||||
jarLocation = "../../infodynamics.jar"
|
||||
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
|
||||
startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
|
||||
|
||||
# Generate some random normalised data.
|
||||
numObservations = 1000
|
||||
covariance=0.4
|
||||
numTrials=10
|
||||
kHistoryLength=1
|
||||
|
||||
# Create a TE calculator and run it:
|
||||
teCalcClass = JPackage("infodynamics.measures.continuous.kraskov").TransferEntropyCalculatorKraskov
|
||||
teCalc = teCalcClass()
|
||||
teCalc.setProperty("k", "4") # Use Kraskov parameter K=4 for 4 nearest points
|
||||
teCalc.initialise(kHistoryLength) # Use target history length of kHistoryLength (Schreiber k)
|
||||
teCalc.startAddObservations()
|
||||
|
||||
for trial in range(0,numTrials):
|
||||
# Create a new trial, with destArray correlated to
|
||||
# previous value of sourceArray:
|
||||
sourceArray = [random.normalvariate(0,1) for r in range(numObservations)]
|
||||
destArray = [0] + [sum(pair) for pair in zip([covariance*y for y in sourceArray[0:numObservations-1]], \
|
||||
[(1-covariance)*y for y in [random.normalvariate(0,1) for r in range(numObservations-1)]] ) ]
|
||||
|
||||
# Add observations for this trial:
|
||||
print("Adding samples from trial %d ..." % trial)
|
||||
teCalc.addObservations(JArray(JDouble, 1)(sourceArray), JArray(JDouble, 1)(destArray))
|
||||
|
||||
# We've finished adding trials:
|
||||
print("Finished adding trials")
|
||||
teCalc.finaliseAddObservations()
|
||||
|
||||
# Compute the result:
|
||||
print("Computing TE ...")
|
||||
result = teCalc.computeAverageLocalOfObservations()
|
||||
# Note that the calculation is a random variable (because the generated
|
||||
# data is a set of random variables) - the result will be of the order
|
||||
# of what we expect, but not exactly equal to it; in fact, there will
|
||||
# be some variance around it (smaller than example 4 since we have more samples).
|
||||
print("TE result %.4f nats; expected to be close to %.4f nats for these correlated Gaussians " % \
|
||||
(result, math.log(1.0/(1-math.pow(covariance,2)))))
|
||||
|
||||
# And here's how to pull the local TEs out corresponding to each input time series.
|
||||
# Normally you would need to track how to split these up yourself -- here
|
||||
# it's easy because our input time series are all of the same length
|
||||
localTEs=teCalc.computeLocalOfPreviousObservations()
|
||||
localValuesPerTrial = int(len(localTEs)/numTrials) # Need to convert to int for indices later
|
||||
for trial in range(0,numTrials):
|
||||
startIndex = localValuesPerTrial*trial
|
||||
endIndex = localValuesPerTrial*(trial+1)-1
|
||||
print("Local TEs for trial %d go from array index %d to %d" % (trial, startIndex, endIndex))
|
||||
print(" corresponding to time points %d:%d (indexed from 0) of that trial" % (kHistoryLength, numObservations-1))
|
||||
# Access the local TEs for this trial as:
|
||||
localTEForThisTrial = localTEs[startIndex:endIndex]
|
||||
|
||||
|
|
@ -0,0 +1,93 @@
|
|||
##
|
||||
## Java Information Dynamics Toolkit (JIDT)
|
||||
## Copyright (C) 2012, Joseph T. Lizier
|
||||
##
|
||||
## This program is free software: you can redistribute it and/or modify
|
||||
## it under the terms of the GNU General Public License as published by
|
||||
## the Free Software Foundation, either version 3 of the License, or
|
||||
## (at your option) any later version.
|
||||
##
|
||||
## This program is distributed in the hope that it will be useful,
|
||||
## but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
## GNU General Public License for more details.
|
||||
##
|
||||
## You should have received a copy of the GNU General Public License
|
||||
## along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
##
|
||||
|
||||
# = Example 9 - Transfer entropy on continuous data using Kraskov estimators with auto-embedding =
|
||||
|
||||
# Transfer entropy (TE) calculation on continuous-valued data using the Kraskov-estimator TE calculator,
|
||||
# with automatic selection of embedding parameters
|
||||
|
||||
|
||||
from jpype import *
|
||||
import random
|
||||
import math
|
||||
import numpy
|
||||
import readFloatsFile
|
||||
|
||||
# Change location of jar to match yours:
|
||||
jarLocation = "../../infodynamics.jar"
|
||||
# Start the JVM (add the "-Xmx" option with say 1024M if you get crashes due to not enough memory space)
|
||||
startJVM(getDefaultJVMPath(), "-ea", "-Djava.class.path=" + jarLocation)
|
||||
|
||||
# Examine the heart-breath interaction that Schreiber originally looked at:
|
||||
datafile = '../data/SFI-heartRate_breathVol_bloodOx.txt'
|
||||
data = readFloatsFile.readFloatsFile(datafile)
|
||||
# As numpy array:
|
||||
A = numpy.array(data)
|
||||
# Select data points 2350:3550, pulling out the relevant columns:
|
||||
breathRate = A[2350:3551,1];
|
||||
heartRate = A[2350:3551,0];
|
||||
|
||||
# Create a Kraskov TE calculator:
|
||||
teCalcClass = JPackage("infodynamics.measures.continuous.kraskov").TransferEntropyCalculatorKraskov
|
||||
teCalc = teCalcClass()
|
||||
|
||||
# Set properties for auto-embedding of both source and destination
|
||||
# using the Ragwitz criteria:
|
||||
# a. Auto-embedding method
|
||||
teCalc.setProperty(teCalcClass.PROP_AUTO_EMBED_METHOD,
|
||||
teCalcClass.AUTO_EMBED_METHOD_RAGWITZ)
|
||||
# b. Search range for embedding dimension (k) and delay (tau)
|
||||
teCalc.setProperty(teCalcClass.PROP_K_SEARCH_MAX, "6")
|
||||
teCalc.setProperty(teCalcClass.PROP_TAU_SEARCH_MAX, "6")
|
||||
# Since we're auto-embedding, no need to supply k, l, k_tau, l_tau here:
|
||||
teCalc.initialise()
|
||||
# Compute TE from breath (column 1) to heart (column 0)
|
||||
teCalc.setObservations(breathRate, heartRate)
|
||||
teBreathToHeart = teCalc.computeAverageLocalOfObservations()
|
||||
|
||||
# Check the auto-selected parameters and print out the result:
|
||||
optimisedK = int(teCalc.getProperty(teCalcClass.K_PROP_NAME))
|
||||
optimisedKTau = int(teCalc.getProperty(teCalcClass.K_TAU_PROP_NAME))
|
||||
optimisedL = int(teCalc.getProperty(teCalcClass.L_PROP_NAME))
|
||||
optimisedLTau = int(teCalc.getProperty(teCalcClass.L_TAU_PROP_NAME))
|
||||
print(("TE(breath->heart) was %.3f nats for (heart embedding:) k=%d," + \
|
||||
"k_tau=%d, (breath embedding:) l=%d,l_tau=%d optimised via Ragwitz criteria") % \
|
||||
(teBreathToHeart, optimisedK, optimisedKTau, optimisedL, optimisedLTau))
|
||||
|
||||
# Next, embed the destination only using the Ragwitz criteria:
|
||||
teCalc.setProperty(teCalcClass.PROP_AUTO_EMBED_METHOD,
|
||||
teCalcClass.AUTO_EMBED_METHOD_RAGWITZ_DEST_ONLY)
|
||||
teCalc.setProperty(teCalcClass.PROP_K_SEARCH_MAX, "6")
|
||||
teCalc.setProperty(teCalcClass.PROP_TAU_SEARCH_MAX, "6")
|
||||
# Since we're only auto-embedding the destination, we supply
|
||||
# source embedding here (to overwrite the auto embeddings from above):
|
||||
teCalc.setProperty(teCalcClass.L_PROP_NAME, "1")
|
||||
teCalc.setProperty(teCalcClass.L_TAU_PROP_NAME, "1")
|
||||
# Since we're auto-embedding, no need to supply k and k_tau here:
|
||||
teCalc.initialise()
|
||||
# Compute TE from breath (column 1) to heart (column 0)
|
||||
teCalc.setObservations(breathRate, heartRate)
|
||||
teBreathToHeartDestEmbedding = teCalc.computeAverageLocalOfObservations()
|
||||
|
||||
# Check the auto-selected parameters and print out the result:
|
||||
optimisedK = int(teCalc.getProperty(teCalcClass.K_PROP_NAME))
|
||||
optimisedKTau = int(teCalc.getProperty(teCalcClass.K_TAU_PROP_NAME))
|
||||
print(("TE(breath->heart) was %.3f nats for (heart embedding:) k=%d," + \
|
||||
"k_tau=%d, optimised via Ragwitz criteria, plus (breath embedding:) l=1,l_tau=1") % \
|
||||
(teBreathToHeartDestEmbedding, optimisedK, optimisedKTau))
|
||||
|
||||
Loading…
Reference in New Issue