mirror of https://github.com/jlizier/jidt
Adding Schreiber example 3 extension to Kraskov code and sample figure for AIS and TE. Also updated main readme file and added readme for the Schreiber demo.
This commit is contained in:
joseph.lizier
parent
8cacb8c660
commit
ebd0c3df17
Binary file not shown.
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 18 KiB |
|
|
@ -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
|
||||
|
||||
=============
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue