-
Notifications
You must be signed in to change notification settings - Fork 15
/
Iccs_mvn.m
executable file
·148 lines (123 loc) · 3.97 KB
/
Iccs_mvn.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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
function Iccs = Iccs_mvn_P2(A, Cfull, varsizes)
% calculate redundancy between a set of Gaussian sources
% from pointwise common change in surprise
% using pairwise marginal maxent solution
%
% A is cell array of sources
% Cfull is full covariance of system
% varsizes specifies the number of variables in each X_i and S (S last)
if sum(varsizes) ~= size(Cfull,1)
error('wrong number of variables specified')
end
if length(varsizes)~=3
error('only 2 variables supported')
end
NA = length(A);
NVs = varsizes(end);
Nx = length(varsizes-1);
NVx = varsizes(1:end-1);
varstart = cumsum(varsizes)+1;
varstart = [1 varstart(1:end-1)];
uniquevars = unique([A{:}]);
sidx = varstart(end):(varstart(end)+NVs-1);
Cs = Cfull(sidx,sidx);
% build Cax for each source
AC = [];
for ai=1:NA
thsA = A{ai};
aidxfull = {};
aidx = {};
thsvstart = 1;
for vi=1:length(thsA)
aidxfull{vi} = varstart(thsA(vi)):(varstart(thsA(vi))+NVx(thsA(vi))-1);
thsL = length(aidxfull{vi});
aidx{vi} = thsvstart:(thsvstart+thsL-1);
thsvstart = thsvstart+thsL;
end
thsNv = length(cell2mat(aidx));
Cas = zeros(thsNv+NVs);
Ca = zeros(thsNv);
% fill in blocks
% diagonal
for vi=1:length(thsA)
Ca(aidx{vi},aidx{vi}) = Cfull(aidxfull{vi},aidxfull{vi});
end
% off diagonal
for vi=1:length(thsA)
for vj=1:length(thsA)
if vi==vj
continue
end
Ca(aidx{vi},aidx{vj}) = Cfull(aidxfull{vi},aidxfull{vj});
end
end
Cas(1:thsNv,1:thsNv) = Ca;
% joint with S
% diagonal
thssidx = thsNv+1:thsNv+NVs;
Cas(thssidx,thssidx) = Cs;
% off diagonal
for vi=1:length(thsA)
Cas(aidx{vi},thssidx) = Cfull(aidxfull{vi},sidx);
Cas(thssidx,aidx{vi}) = Cfull(sidx,aidxfull{vi});
end
Casoff = Cas(1:thsNv,thssidx);
CXYY1 = Casoff * pinv(Cs);
Cacs = Ca - CXYY1*Cas(thssidx,1:thsNv);
MacsF = CXYY1;
AC(ai).Ca = Ca;
AC(ai).Cas = Cas;
AC(ai).Cacs = Cacs;
AC(ai).Casoff = Casoff;
AC(ai).MacsF = CXYY1;
AC(ai).Nv = thsNv;
end
if NA==1
% use closed form expression
chA = chol(AC(1).Ca);
chS = chol(Cs);
chAS = chol(AC(1).Cas);
% normalisations cancel for information
HA = sum(log(diag(chA))); % + 0.5*Nvarx*log(2*pi*exp(1));
HS = sum(log(diag(chS))); % + 0.5*Nvary*log(2*pi*exp(1));
HAS = sum(log(diag(chAS))); % + 0.5*(Nvarx+Nvary)*log(2*pi*exp(1));
Iccs = (HA + HS - HAS) / log(2);
end
if NA==2
% Covariance for Pind(A1,A2)
thsNv = AC(1).Nv + AC(2).Nv + NVs;
ANv = AC(1).Nv + AC(2).Nv;
a1idx = 1:AC(1).Nv;
a2idx = AC(1).Nv+1:AC(1).Nv+AC(2).Nv;
a12idx = 1:AC(1).Nv + AC(2).Nv;
% P2 == full gaussian covariance
C = Cfull(a12idx,a12idx);
Caasoff = Cfull(a12idx,sidx);
CXYY1 = Caasoff * pinv(Cs);
Caacs = C - CXYY1*Cfull(sidx,a12idx);
MaacsF = CXYY1;
thssidx = AC(1).Nv+AC(2).Nv+1:AC(1).Nv+AC(2).Nv+NVs;
% Monte Carlo Integration for Iccs
% 100,000 works well for 5d space.
% 10,000 might be ok but some variance
Nmc = 100000;
% integrate over full space
intD = AC(1).Nv + AC(2).Nv + NVs;
% sample from multivariate gaussian we want expectation over.
mcx = mvnrnd(zeros(1,ANv+NVs), Cfull, Nmc);
px = logmvnpdf(mcx(:,a1idx),zeros(1,AC(1).Nv),AC(1).Ca);
pxcs = logmvnpdf(mcx(:,a1idx), (AC(1).MacsF*mcx(:,thssidx)')', AC(1).Cacs);
py = logmvnpdf(mcx(:,a2idx),zeros(1,AC(2).Nv),AC(2).Ca);
pycs = logmvnpdf(mcx(:,a2idx), (AC(2).MacsF*mcx(:,thssidx)')', AC(2).Cacs);
pxy = logmvnpdf(mcx(:,a12idx), zeros(1,ANv), C);
pxycs = logmvnpdf(mcx(:,a12idx), (MaacsF*mcx(:,thssidx)')', Caacs);
dhx = pxcs - px;
dhy = pycs - py;
dhxy = pxycs - pxy;
lnii = dhx + dhy - dhxy;
keep = sign(dhx) == sign(lnii) & sign(dhx) == sign(dhy) & sign(dhx) == sign(dhxy);
lnii(~keep) = 0;
Iccs = nanmean(lnii);
% convert to bits
Iccs = Iccs ./ log(2);
end