-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathanlys_KLD.m
77 lines (59 loc) · 2 KB
/
anlys_KLD.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
function [KLD] = anlys_KLD(currycdr, curryloc, source, Phi)
bglevel = 0.01; % background level
[P, T] = findpeaks(std(Phi, 1, 1));
T = T(P >= max(P)/sqrt(10));
%%
for t = 1 : length(T)
cdr = norms(squeeze(currycdr(4, :, :)), 2, 2).';
% cdr = std(squeeze(currycdr(4, :, :)), 1, 2).';
% cdr = abs(squeeze(currycdr(4, :, T(t))));
if(sum(cdr ~= 0))
cdr(cdr == 0) = bglevel * min(cdr(cdr ~= 0)) / sum(cdr == 0);
%%
ref = zeros(size(cdr));
for i = 1 : length(source.epl)
[~, ind] = find_nvoxel(source.epl{i}, curryloc);
ref(ind) = norms(source.sdm(i, :), 2, 2);
% ref(ind) = std(source.sdm(i, :), 1, 2);
% ref(ind) = abs(source.sdm(i, T(t)));
end
ref(ref == 0) = bglevel * min(ref(ref ~= 0)) / sum(ref == 0);
%%
KLD(t) = (KLDiv(ref, cdr)) / 2;
% KLD = (KLDiv(ref, cdr) + KLDiv(cdr, ref)) / 2;
else
KLD(t) = inf;
end
end
KLD = mean(KLD);
end
function dist=KLDiv(P,Q)
% dist = KLDiv(P,Q) Kullback-Leibler divergence of two discrete probability
% distributions
% P and Q are automatically normalised to have the sum of one on rows
% have the length of one at each
% P = n x nbins
% Q = 1 x nbins or n x nbins(one to one)
% dist = n x 1
% downloaded from Mathworks File Exchange
if size(P,2)~=size(Q,2)
error('the number of columns in P and Q should be the same');
end
if sum(~isfinite(P(:))) + sum(~isfinite(Q(:)))
error('the inputs contain non-finite values!')
end
% normalizing the P and Q
if size(Q,1)==1
Q = Q ./sum(Q);
P = P ./repmat(sum(P,2),[1 size(P,2)]);
temp = P.*log(P./repmat(Q,[size(P,1) 1]));
temp(isnan(temp))=0;% resolving the case when P(i)==0
dist = sum(temp,2);
elseif size(Q,1)==size(P,1)
Q = Q ./repmat(sum(Q,2),[1 size(Q,2)]);
P = P ./repmat(sum(P,2),[1 size(P,2)]);
temp = P.*log(P./Q);
temp(isnan(temp))=0; % resolving the case when P(i)==0
dist = sum(temp,2);
end
end