Added Schreiber heart-breath rate example using Kraskov estimator.

This commit is contained in:
joseph.lizier 2014-01-29 05:33:08 +00:00
parent 2b8ab0e939
commit cee9b3dbd3
1 changed files with 73 additions and 0 deletions

View File

@ -0,0 +1,73 @@
% function [teHeartToBreath, teBreathToHeart] = runHeartBreathRateKraskov(knns)
%
% runHeartBreathRateKraskov
% Version 1.0
% Joseph Lizier
% 22/1/2014
%
% Used to explore information transfer in the heart rate / breath rate example of Schreiber --
% but estimates TE using Kraskov-Grassberger estimation.
%
% Runtime is ~3 minutes per knn value (as there is no fast nearest neighbour search implemented yet in the toolkit)
%
% Inputs
% - knns - a scalar specifying a single, or vector specifying multiple, value of K nearest neighbours to evaluate TE (Kraskov) with
% 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] = runHeartBreathRate(knns)
starttime = time;
% 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');
fprintf('TE for heart rate <-> breath rate for Kraskov estimation:\n');
data = load('../../data/SFI-heartRate_breathVol_bloodOx.txt');
% Restrict to the samples that Schreiber mentions:
heartRate = data(2350:3550,:);
% Separate the data from each column:
heart = data(:,1);
chestVol = data(:,2);
bloodOx = data(:,3);
timeSteps = length(heart);
% Using a single conditional mutual information calculator is the least biased way to run this:
teCalc=javaObject('infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov1');
for knnIndex = 1:length(knns)
knn = knns(knnIndex);
% Compute a TE value for knn nearest neighbours
% Perform calculation for heart -> breath (lag 1)
teCalc.initialise(1,1,1); % Use history length 1 (Schreiber k)
teCalc.setProperty('k', sprintf('%d',knn));
teCalc.setObservations(octaveToJavaDoubleMatrix(heart(1:timeSteps-1)), ...
octaveToJavaDoubleMatrix(chestVol(2:timeSteps)), ...
octaveToJavaDoubleMatrix(chestVol(1:timeSteps-1)));
teHeartToBreath(knnIndex) = teCalc.computeAverageLocalOfObservations();
% Perform calculation for breath -> heart (lag 1)
teCalc.initialise(1,1,1); % Use history length 1 (Schreiber k)
teCalc.setProperty('k', sprintf('%d',knn));
teCalc.setObservations(octaveToJavaDoubleMatrix(chestVol(1:timeSteps-1)), ...
octaveToJavaDoubleMatrix(heart(2:timeSteps)), ...
octaveToJavaDoubleMatrix(heart(1:timeSteps-1)));
teBreathToHeart(knnIndex) = teCalc.computeAverageLocalOfObservations();
fprintf('TE(k=%d): heart->breath = %.3f, breath->heart = %.3f\n', knn, teHeartToBreath(knnIndex), teBreathToHeart(knnIndex));
end
endtime = time;
printf("Total runtime was %.1f sec\n", endtime - starttime);
end