diff --git a/demos/octave/SchreiberTransferEntropyExamples/README-SchreiberTeDemos.pdf b/demos/octave/SchreiberTransferEntropyExamples/README-SchreiberTeDemos.pdf new file mode 100755 index 0000000..cfdace4 Binary files /dev/null and b/demos/octave/SchreiberTransferEntropyExamples/README-SchreiberTeDemos.pdf differ diff --git a/demos/octave/SchreiberTransferEntropyExamples/activeInfoStorageHeartBreathRatesKraskov.m b/demos/octave/SchreiberTransferEntropyExamples/activeInfoStorageHeartBreathRatesKraskov.m new file mode 100755 index 0000000..cb8a613 --- /dev/null +++ b/demos/octave/SchreiberTransferEntropyExamples/activeInfoStorageHeartBreathRatesKraskov.m @@ -0,0 +1,109 @@ +% function [aisHeart, aisBreath] = activeInfoStorageHeartBreathRatesKraskov(kHistories, knn, numSurrogates) +% +% activeInfoStorageHeartBreathRatesKraskov +% Version 1.0 +% Joseph Lizier +% 04/04/2014 +% +% Used to explore active information storage in the heart rate / breath rate example of Schreiber -- +% estimated using Kraskov-Grassberger estimation. +% +% Usually plateaus of AIS indicate that the correct embedding is found; for Kraskov estimation, +% a peak will indicate this (given that bias correction will pull down values for larger k. +% +% +% Inputs +% - kHistories - a vector of which embedded history lengths to evaluate for both variables +% - knn - a scalar specifying a single value of K nearest neighbours to evaluate AIS (Kraskov) with +% - numSurrogates - a scalar specifying the number of surrogates to evaluate AIS from null distribution +% Outputs +% - aisHeart - active information storage TE (heart -> breath) for each value of k nearest neighbours +% - aisBreath - TE (breath -> heart) for each value of k nearest neighbours + + +function [aisHeart, aisBreath] = activeInfoStorageHeartBreathRatesKraskov(kHistories, knn, numSurrogates) + + tic; + + % Add utilities to the path + addpath('..'); + + % Assumes the jar is two levels up - change this if this is not the case + % Octave is happy to have the path added multiple times; I'm unsure if this is true for matlab + javaaddpath('../../../infodynamics.jar'); + + if (nargin < 3) + numSurrogates = 0; + end + + data = load('../../data/SFI-heartRate_breathVol_bloodOx.txt'); + + % Restrict to the samples that Schreiber mentions: + data = data(2350:3550,:); + + % Separate the data from each column: + heart = data(:,1); + chestVol = data(:,2); + bloodOx = data(:,3); + timeSteps = length(heart); + + fprintf('AIS for heart rate and breath rate for Kraskov estimation with %d samples:\n', timeSteps); + + % Using a single conditional mutual information calculator is the least biased way to run this: + aisCalc=javaObject('infodynamics.measures.continuous.kraskov.ActiveInfoStorageCalculatorKraskov'); + + for kIndex = 1:length(kHistories) + kHistory = kHistories(kIndex); + % Compute an AIS value for this embedding length for each variable: + + % Perform calculation for heart + aisCalc.initialise(kHistory); % Use history length kHistory (Schreiber k) + aisCalc.setProperty('k', sprintf('%d',knn)); + aisCalc.setProperty('NORMALISE', 'true'); + aisCalc.setObservations(octaveToJavaDoubleArray(heart(1:timeSteps))); + aisHeart(kIndex) = aisCalc.computeAverageLocalOfObservations(); + if (numSurrogates > 0) + aisHeartNullDist = aisCalc.computeSignificance(numSurrogates); + aisHeartNullMean = aisHeartNullDist.getMeanOfDistribution(); + aisHeartNullStd = aisHeartNullDist.getStdOfDistribution(); + end + + % Perform calculation for breath + aisCalc.initialise(kHistory); % Use history length kHistory (Schreiber k) + aisCalc.setProperty('k', sprintf('%d',knn)); + aisCalc.setProperty('NORMALISE', 'true'); + aisCalc.setObservations(octaveToJavaDoubleArray(chestVol(1:timeSteps))); + aisBreath(kIndex) = aisCalc.computeAverageLocalOfObservations(); + if (numSurrogates > 0) + aisBreathNullDist = aisCalc.computeSignificance(numSurrogates); + aisBreathNullMean = aisBreathNullDist.getMeanOfDistribution(); + aisBreathNullStd = aisBreathNullDist.getStdOfDistribution(); + end + + fprintf('AIS(k=%d,knns=%d): h = %.3f', kHistory, knn, aisHeart(kIndex)); + if (numSurrogates > 0) + fprintf(' (null = %.3f +/- %.3f)', aisHeartNullMean, aisHeartNullStd); + end + fprintf(', b = %.3f', aisBreath(kIndex)); + if (numSurrogates > 0) + fprintf('(null = %.3f +/- %.3f)\n', aisBreathNullMean, aisBreathNullStd); + else + fprintf('\n'); + end + end + + totaltime = toc; + fprintf('Total runtime was %.1f sec\n', totaltime); + + hold off; + plot(kHistories, aisHeart, 'rx'); + hold on; + plot(kHistories, aisBreath, 'bo'); + hold off; + legend(['AIS(Heart) '; 'AIS(Breath)']) + set (gca,'fontsize',26); + xlabel('embedding history k', 'FontSize', 36, 'FontWeight', 'bold'); + ylabel('AIS(k)', 'FontSize', 36, 'FontWeight', 'bold'); + print('heartBreathResults-kraskovAIS.eps', '-depsc'); +end + diff --git a/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m b/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m index ddf27ce..f9e9f27 100755 --- a/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m +++ b/demos/octave/SchreiberTransferEntropyExamples/runHeartBreathRateKraskov.m @@ -13,22 +13,27 @@ % - kHistory - destination embedding length % - lHistory - source embedding length % - knns - a scalar specifying a single, or vector specifying multiple, value of K nearest neighbours to evaluate TE (Kraskov) with +% - numSurrogates - a scalar specifying the number of surrogates to evaluate AIS from null distribution % Outputs % - teHeartToBreath - TE (heart -> breath) for each value of k nearest neighbours % - teBreathToHeart - TE (breath -> heart) for each value of k nearest neighbours -function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(kHistory, lHistory, knns) +function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(kHistory, lHistory, knns, numSurrogates) starttime = tic; % Add utilities to the path addpath('..'); - % Assumes the jar is two levels up - change this if this is not the case + % Assumes the jar is three levels up - change this if this is not the case % Octave is happy to have the path added multiple times; I'm unsure if this is true for matlab javaaddpath('../../../infodynamics.jar'); + if (nargin < 4) + numSurrogates = 0; + end + data = load('../../data/SFI-heartRate_breathVol_bloodOx.txt'); % Restrict to the samples that Schreiber mentions: @@ -55,6 +60,11 @@ function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(kHistory teCalc.setObservations(octaveToJavaDoubleArray(heart), ... octaveToJavaDoubleArray(chestVol)); teHeartToBreath(knnIndex) = teCalc.computeAverageLocalOfObservations(); + if (numSurrogates > 0) + teHeartToBreathNullDist = teCalc.computeSignificance(numSurrogates); + teHeartToBreathNullMean = teHeartToBreathNullDist.getMeanOfDistribution(); + teHeartToBreathNullStd = teHeartToBreathNullDist.getStdOfDistribution(); + end % Perform calculation for breath -> heart (lag 1) teCalc.initialise(kHistory,1,lHistory,1,1); @@ -62,13 +72,28 @@ function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(kHistory teCalc.setObservations(octaveToJavaDoubleArray(chestVol), ... octaveToJavaDoubleArray(heart)); teBreathToHeart(knnIndex) = teCalc.computeAverageLocalOfObservations(); + if (numSurrogates > 0) + teBreathToHeartNullDist = teCalc.computeSignificance(numSurrogates); + teBreathToHeartNullMean = teBreathToHeartNullDist.getMeanOfDistribution(); + teBreathToHeartNullStd = teBreathToHeartNullDist.getStdOfDistribution(); + end - fprintf('TE(k=%d): heart->breath = %.3f, breath->heart = %.3f nats\n', knn, teHeartToBreath(knnIndex), teBreathToHeart(knnIndex)); + fprintf('TE(k=%d,l=%d,knn=%d): h->b = %.3f', kHistory, lHistory, knn, teHeartToBreath(knnIndex)); + if (numSurrogates > 0) + fprintf(' (null = %.3f +/- %.3f)', teHeartToBreathNullMean, teHeartToBreathNullStd); + end + fprintf(', b->h = %.3f nats', teBreathToHeart(knnIndex)); + if (numSurrogates > 0) + fprintf('(null = %.3f +/- %.3f)\n', teBreathToHeartNullMean, teBreathToHeartNullStd); + else + fprintf('\n'); + end end tElapsed = toc(starttime); fprintf('Total runtime was %.1f sec\n', tElapsed); + hold off; plot(knns, teHeartToBreath, 'rx', 'markersize', 10); hold on; plot(knns, teBreathToHeart, 'mo', 'markersize', 10); @@ -77,6 +102,6 @@ function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(kHistory set (gca,'fontsize',26); xlabel('K nearest neighbours', 'FontSize', 44, 'FontWeight', 'bold'); ylabel('TE', 'FontSize', 44, 'FontWeight', 'bold'); - print('results.eps', '-deps', '-color'); + % print('results.eps', '-deps', '-color'); end diff --git a/demos/octave/SchreiberTransferEntropyExamples/sampleResults/heartBreathResults-kraskovAIS.png b/demos/octave/SchreiberTransferEntropyExamples/sampleResults/heartBreathResults-kraskovAIS.png new file mode 100755 index 0000000..0e20510 Binary files /dev/null and b/demos/octave/SchreiberTransferEntropyExamples/sampleResults/heartBreathResults-kraskovAIS.png differ diff --git a/readme-template.txt b/readme-template.txt index 18ff430..4cad671 100755 --- a/readme-template.txt +++ b/readme-template.txt @@ -57,12 +57,12 @@ That's it. Documentation ============= +A research paper describing the toolkit is included in the top level directory -- "InfoDynamicsToolkit.pdf". + Javadocs for the toolkit are included in the full distribution at javadocs. They can also be generated using "ant javadocs" (useful if you are on an SVN view). Further, they will soon be posted on the web. -A research paper describing the toolkit and its use is in preparation and will be included in the distribution in future. - Further documentation is provided by the Usage examples below. ============= @@ -79,9 +79,11 @@ Several sets of demonstration code are distributed with the toolkit: d. demos/octave/CellularAutomata -- using the Java toolkit to plot local information dynamics profiles in cellular automata; the toolkit is run under Octave or Matlab -- see description at http://code.google.com/p/information-dynamics-toolkit/wiki/CellularAutomataDemos - e. demos/octave/DetectingInteractionLags -- brief examples using the transfer entropy to examine source-delay interaction lags. Documentation to come soon; in the interim, see header comments in the .m files. + e. demos/octave/SchreiberTransferEntropyExamples -- recreates the transfer entropy examples in Schreiber's original paper presenting this measure; shows the correct parameter settings to reproduce these results -- see description at http://code.google.com/p/information-dynamics-toolkit/wiki/SchreiberTeDemos + + f. demos/octave/DetectingInteractionLags -- brief examples using the transfer entropy to examine source-delay interaction lags. Documentation to come soon; in the interim, see header comments in the .m files. - f. java/unittests -- the JUnit test cases for the Java toolkit are included in the distribution -- these case also be browsed to see simple use cases for the various calculators in the toolkit -- see description at http://code.google.com/p/information-dynamics-toolkit/wiki/JUnitTestCases + g. java/unittests -- the JUnit test cases for the Java toolkit are included in the distribution -- these case also be browsed to see simple use cases for the various calculators in the toolkit -- see description at http://code.google.com/p/information-dynamics-toolkit/wiki/JUnitTestCases =============