-
Notifications
You must be signed in to change notification settings - Fork 0
/
SI_Gauss.m
68 lines (54 loc) · 1.65 KB
/
SI_Gauss.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
function [SI, I, H_cond] = SI_Gauss(Cov_X,Cov_XY,Cov_Y,Z)
%-----------------------------------------------------------------------
% FUNCTION: SI_Gauss.m
% PURPOSE: calculate stochastic interaction given covariance of data
% (Ay, 2001; Ay, 2015, Entropy; Barrett & Seth, 2011, PLoS Comp Biol)
%
% INPUTS:
% Cov_X: covariance of data X (past, t-tau)
% Cov_XY: cross-covariance of X (past, t-tau) and Y (present, t)
% Cov_Y: covariance of data Y (present, t)
% Z: partition of each channel (default: atomic partition)
%
% OUTPUT:
% SI: stochastic interaction SI(Y|X)
% I: mutual information I(X;Y)
% H: entropy, H(Y)
%
% Masafumi Oizumi, 2017
%-----------------------------------------------------------------------
N = size(Cov_X,1); % number of channels
if nargin < 3
Cov_Y = Cov_X;
end
if nargin < 4
Z = 1: 1: N;
end
Cov_Y_X = Cov_cond(Cov_Y,Cov_XY',Cov_X); % conditional covariance matrix
H_cond = H_gauss(Cov_Y_X);
if isinf(H_cond) == 1
fprintf('Alert: Infinite Entropy\n')
end
if isreal(H_cond) == 0
fprintf('Alert: Complex Entropy\n')
end
H = H_gauss(Cov_Y); % entropy
I = H - H_cond; % mutual information
%%
N_c = max(Z); % number of clusters
H_cond_p = zeros(N_c,1);
M_cell = cell(N_c,1);
for i=1: N_c
M_cell{i} = find(Z==i);
end
for i=1: N_c
M = M_cell{i};
Cov_X_p = Cov_X(M,M);
Cov_Y_p = Cov_Y(M,M);
Cov_XY_p = Cov_XY(M,M);
Cov_Y_X_p = Cov_cond(Cov_Y_p,Cov_XY_p',Cov_X_p);
H_cond_p(i) = H_gauss(Cov_Y_X_p);
end
%% stochastic interaction
SI = sum(H_cond_p) - H_cond;
end