Browse code

starting som prediction fine-tuned class-performance visualisation

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@112 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on21/01/2009 16:34:25
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,138 @@
1
+function [adm,admu,tdmu] = som_distortion(sM, D, arg1, arg2)
2
+
3
+%SOM_DISTORTION Calculate distortion measure for the map.
4
+%
5
+% [adm,admu,tdmu] = som_distortion(sMap, D, [radius], ['prob'])
6
+%
7
+%  adm = som_distortion(sMap,D);
8
+%  [adm,admu] = som_distortion(sMap,D);
9
+%  som_show(sMap,'color',admu);
10
+%
11
+%  Input and output arguments: 
12
+%   sMap     (struct) a map struct
13
+%   D        (struct) a data struct
14
+%            (matrix) size dlen x dim, a data matrix
15
+%   [radius] (scalar) neighborhood function radius to be used.
16
+%                     Defaults to the last radius_fin in the 
17
+%                     trainhist field of the map struct, or 1 if
18
+%                     that is missing.
19
+%   ['prob'] (string) If given, this argument forces the 
20
+%                     neigborhood function values for each map
21
+%                     unit to be normalized so that they sum to 1.
22
+%
23
+%   adm      (scalar) average distortion measure (sum(dm)/dlen)
24
+%   admu     (vector) size munits x 1, average distortion in each unit 
25
+%   tdmu     (vector) size munits x 1, total distortion for each unit
26
+%
27
+% The distortion measure is defined as: 
28
+%                                           2
29
+%    E = sum sum h(bmu(i),j) ||m(j) - x(i)|| 
30
+%         i   j    
31
+% 
32
+% where m(i) is the ith prototype vector of SOM, x(j) is the jth data
33
+% vector, and h(.,.) is the neighborhood function. In case of fixed
34
+% neighborhood and discreet data, the distortion measure can be
35
+% interpreted as the energy function of the SOM. Note, though, that
36
+% the learning rule that follows from the distortion measure is
37
+% different from the SOM training rule, so SOM only minimizes the
38
+% distortion measure approximately.
39
+% 
40
+% If the 'prob' argument is given, the distortion measure can be 
41
+% interpreted as an expected quantization error when the neighborhood 
42
+% function values give the likelyhoods of accidentally assigning 
43
+% vector j to unit i. The normal quantization error is a special case 
44
+% of this with zero incorrect assignement likelihood. 
45
+% 
46
+% NOTE: when calculating BMUs and distances, the mask of the given 
47
+%       map is used.
48
+%
49
+% See also SOM_QUALITY, SOM_BMUS, SOM_HITS.
50
+
51
+% Reference: Kohonen, T., "Self-Organizing Map", 2nd ed., 
52
+%    Springer-Verlag, Berlin, 1995, pp. 120-121.
53
+%
54
+%    Graepel, T., Burger, M. and Obermayer, K., 
55
+%    "Phase Transitions in Stochastic Self-Organizing Maps",
56
+%    Physical Review E, Vol 56, No 4, pp. 3876-3890 (1997).
57
+
58
+% Contributed to SOM Toolbox vs2, Feb 3rd, 2000 by Juha Vesanto
59
+% Copyright (c) by Juha Vesanto
60
+% http://www.cis.hut.fi/projects/somtoolbox/
61
+
62
+% Version 2.0beta juuso 030200
63
+
64
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65
+%% check arguments
66
+
67
+% input arguments
68
+if nargin < 2, error('Not enough input arguments.'); end
69
+
70
+% map
71
+M = sM.codebook;
72
+munits = prod(sM.topol.msize);
73
+
74
+% data
75
+if isstruct(D), D = D.data; end
76
+[dlen dim] = size(D);
77
+
78
+% arg1, arg2
79
+rad = NaN;
80
+normalize = 0;
81
+if nargin>2, 
82
+  if isnumeric(arg1), rad = arg1;
83
+  elseif ischar(arg1) & strcmp(arg1,'prob'), normalize = 0;
84
+  end
85
+end
86
+if nargin>3, 
87
+  if isnumeric(arg2), rad = arg2;
88
+  elseif ischar(arg2) & strcmp(arg2,'prob'), normalize = 0;
89
+  end
90
+end
91
+
92
+% neighborhood radius
93
+if isempty(rad) | isnan(rad), 
94
+  if ~isempty(sM.trainhist), rad = sM.trainhist(end).radius_fin;
95
+  else rad = 1; 
96
+  end
97
+end
98
+if rad<eps, rad = eps; end
99
+
100
+% neighborhood  
101
+Ud = som_unit_dists(sM.topol); 
102
+switch sM.neigh, 
103
+ case 'bubble',   H = (Ud <= rad);
104
+ case 'gaussian', H = exp(-(Ud.^2)/(2*rad*rad)); 
105
+ case 'cutgauss', H = exp(-(Ud.^2)/(2*rad*rad)) .* (Ud <= rad);
106
+ case 'ep',       H = (1 - (Ud.^2)/rad) .* (Ud <= rad);
107
+end  
108
+if normalize, 
109
+  for i=1:munits, H(:,i) = H(:,i)/sum(H(:,i)); end
110
+end
111
+
112
+% total distortion measure
113
+mu_x_1 = ones(munits,1);
114
+tdmu = zeros(munits,1);
115
+hits = zeros(munits,1);
116
+for i=1:dlen,
117
+  x = D(i,:);                        % data sample
118
+  known = ~isnan(x);                 % its known components
119
+  Dx = M(:,known) - x(mu_x_1,known); % each map unit minus the vector
120
+  dist2 = (Dx.^2)*sM.mask(known);    % squared distances  
121
+  [qerr bmu] = min(dist2);           % find BMU
122
+  tdmu = tdmu + dist2.*H(:,bmu);     % add to distortion measure
123
+  hits(bmu) = hits(bmu)+1;           % add to hits
124
+end 
125
+
126
+% average distortion per unit
127
+admu = tdmu; 
128
+ind = find(hits>0);
129
+admu(ind) = admu(ind) ./ hits(ind);
130
+  
131
+% average distortion measure
132
+adm = sum(tdmu)/dlen;
133
+
134
+return;
135
+
136
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137
+
138
+