-
Notifications
You must be signed in to change notification settings - Fork 0
/
fMIfromDensities.m
38 lines (30 loc) · 1.23 KB
/
fMIfromDensities.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
function [I,Hr_s,Hr] = fMIfromDensities(Fcond,q)
%FMIFROMDENSITIES Calculates the mutual information from given
%signal-response conditional probability density matrix
%
% Fcond is be a N x N cell array of conditional probability densities
% the cell i,j has the probability to get the observations seen in i given
% stimulus j.
% q is the probability for inputs
Nq = numel(q);
assert(size(Fcond,1)==size(Fcond,2),'F should be square')
assert(size(Fcond,1)==Nq,'Fcond should be size of Nq x Nq');
%% First calculate the non-conditional probability densitites
Q = repmat(q(:)',Nq,1);
FtimesQ = cellfun(@(f,q) f(:)*q,Fcond,num2cell(Q),'uniformoutput',0);
F = cell(Nq,1);
for i=1:Nq
F{i} = sum(cat(2,FtimesQ{i,:}),2);
end
Fself = Fcond(eye(Nq)==1);
%% calculate Entropies for conditional case (easy)
fH = @(F) -sum(log2(F(F>eps)))/nnz(F>eps);
Hr_s = cellfun(fH,Fself);
%% create weights for a weighted average of the contribution of the non-conditional densities
Wq = cellfun(@(q,f) ones(size(f))*q/length(ones(size(f))),num2cell(q),F,'uniformoutput',0);
Wq = cat(1,Wq{:});
Fall = cat(1,F{:});
Fall(Fall<=eps)=nan;
%Wq = Wq./nansum(Wq);
Hr = -nansum(log2(Fall).*Wq);
I = Hr-sum(q.*Hr_s);