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,282 @@
1
+function [P] = cca(D, P, epochs, Mdist, alpha0, lambda0)
2
+
3
+%CCA Projects data vectors using Curvilinear Component Analysis.
4
+%
5
+% P = cca(D, P, epochs, [Dist], [alpha0], [lambda0])
6
+%
7
+%  P = cca(D,2,10);           % projects the given data to a plane
8
+%  P = cca(D,pcaproj(D,2),5); % same, but with PCA initialization
9
+%  P = cca(D, 2, 10, Dist);   % same, but the given distance matrix is used
10
+%  
11
+%  Input and output arguments ([]'s are optional):
12
+%   D          (matrix) the data matrix, size dlen x dim
13
+%              (struct) data or map struct            
14
+%   P          (scalar) output dimension
15
+%              (matrix) size dlen x odim, the initial projection
16
+%   epochs     (scalar) training length
17
+%   [Dist]     (matrix) pairwise distance matrix, size dlen x dlen.
18
+%                       If the distances in the input space should
19
+%                       be calculated otherwise than as euclidian
20
+%                       distances, the distance from each vector
21
+%                       to each other vector can be given here,
22
+%                       size dlen x dlen. For example PDIST
23
+%                       function can be used to calculate the
24
+%                       distances: Dist = squareform(pdist(D,'mahal'));
25
+%   [alpha0]   (scalar) initial step size, 0.5 by default
26
+%   [lambda0]  (scalar) initial radius of influence, 3*max(std(D)) by default
27
+%  
28
+%   P          (matrix) size dlen x odim, the projections
29
+%
30
+% Unknown values (NaN's) in the data: projections of vectors with
31
+% unknown components tend to drift towards the center of the
32
+% projection distribution. Projections of totally unknown vectors are
33
+% set to unknown (NaN).
34
+%
35
+% See also SAMMON, PCAPROJ. 
36
+
37
+% Reference: Demartines, P., Herault, J., "Curvilinear Component
38
+%   Analysis: a Self-Organizing Neural Network for Nonlinear
39
+%   Mapping of Data Sets", IEEE Transactions on Neural Networks, 
40
+%   vol 8, no 1, 1997, pp. 148-154.
41
+
42
+% Contributed to SOM Toolbox 2.0, February 2nd, 2000 by Juha Vesanto
43
+% Copyright (c) by Juha Vesanto
44
+% http://www.cis.hut.fi/projects/somtoolbox/
45
+
46
+% juuso 171297 040100
47
+
48
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
+%% Check arguments 
50
+
51
+error(nargchk(3, 6, nargin)); % check the number of input arguments
52
+
53
+% input data
54
+if isstruct(D), 
55
+  if strcmp(D.type,'som_map'), D = D.codebook; else D = D.data; end
56
+end
57
+[noc dim] = size(D);
58
+noc_x_1  = ones(noc, 1); % used frequently
59
+me = zeros(1,dim); st = zeros(1,dim);
60
+for i=1:dim,
61
+  me(i) = mean(D(find(isfinite(D(:,i))),i));
62
+  st(i) = std(D(find(isfinite(D(:,i))),i));
63
+end
64
+
65
+% initial projection
66
+if prod(size(P))==1, 
67
+  P = (2*rand(noc,P)-1).*st(noc_x_1,1:P) + me(noc_x_1,1:P); 
68
+else
69
+  % replace unknown projections with known values
70
+  inds = find(isnan(P)); P(inds) = rand(size(inds));
71
+end
72
+[dummy odim] = size(P);
73
+odim_x_1  = ones(odim, 1); % this is used frequently
74
+
75
+% training length
76
+train_len = epochs*noc;
77
+
78
+% random sample order
79
+rand('state',sum(100*clock));
80
+sample_inds = ceil(noc*rand(train_len,1));
81
+
82
+% mutual distances
83
+if nargin<4 | isempty(Mdist) | all(isnan(Mdist(:))),
84
+  fprintf(2, 'computing mutual distances\r');
85
+  dim_x_1 = ones(dim,1);
86
+  for i = 1:noc,
87
+    x = D(i,:); 
88
+    Diff = D - x(noc_x_1,:);
89
+    N = isnan(Diff);
90
+    Diff(find(N)) = 0; 
91
+    Mdist(:,i) = sqrt((Diff.^2)*dim_x_1);
92
+    N = find(sum(N')==dim); %mutual distance unknown
93
+    if ~isempty(N), Mdist(N,i) = NaN; end
94
+  end
95
+else
96
+  % if the distance matrix is output from PDIST function
97
+  if size(Mdist,1)==1, Mdist = squareform(Mdist); end
98
+  if size(Mdist,1)~=noc, 
99
+    error('Mutual distance matrix size and data set size do not match'); 
100
+  end
101
+end
102
+
103
+% alpha and lambda
104
+if nargin<5 | isempty(alpha0) | isnan(alpha0), alpha0 = 0.5; end
105
+alpha = potency_curve(alpha0,alpha0/100,train_len);
106
+
107
+if nargin<6 | isempty(lambda0) | isnan(lambda0), lambda0 = max(st)*3; end
108
+lambda = potency_curve(lambda0,0.01,train_len);
109
+
110
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111
+%% Action
112
+
113
+k=0; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
114
+
115
+for i=1:train_len, 
116
+  
117
+  ind = sample_inds(i);     % sample index
118
+  dx = Mdist(:,ind);        % mutual distances in input space
119
+  known = find(~isnan(dx)); % known distances
120
+
121
+  if ~isempty(known),
122
+    % sample vector's projection
123
+    y = P(ind,:);                 
124
+
125
+    % distances in output space
126
+    Dy = P(known,:) - y(noc_x_1(known),:); 
127
+    dy = sqrt((Dy.^2)*odim_x_1);           
128
+  
129
+    % relative effect
130
+    dy(find(dy==0)) = 1;        % to get rid of div-by-zero's
131
+    fy = exp(-dy/lambda(i)) .* (dx(known) ./ dy - 1);
132
+
133
+    % Note that the function F here is e^(-dy/lambda)) 
134
+    % instead of the bubble function 1(lambda-dy) used in the 
135
+    % paper.
136
+    
137
+    % Note that here a simplification has been made: the derivatives of the
138
+    % F function have been ignored in calculating the gradient of error
139
+    % function w.r.t. to changes in dy.
140
+    
141
+    % update
142
+    P(known,:) = P(known,:) + alpha(i)*fy(:,odim_x_1).*Dy;
143
+  end
144
+
145
+  % track
146
+  if rem(i,noc)==0, 
147
+    k=k+1; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
148
+  end
149
+
150
+end
151
+
152
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
+%% clear up
154
+
155
+% calculate error
156
+error = cca_error(P,Mdist,lambda(train_len));
157
+fprintf(2,'%d iterations, error %f          \n', epochs, error);
158
+
159
+% set projections of totally unknown vectors as unknown
160
+unknown = find(sum(isnan(D)')==dim);
161
+P(unknown,:) = NaN;
162
+
163
+return;
164
+
165
+
166
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
+%% tips
168
+
169
+% to plot the results, use the code below
170
+
171
+%subplot(2,1,1), 
172
+%switch(odim), 
173
+%  case 1, plot(P(:,1),ones(dlen,1),'x')
174
+%  case 2, plot(P(:,1),P(:,2),'x'); 
175
+%  otherwise, plot3(P(:,1),P(:,2),P(:,3),'x'); rotate3d on
176
+%end
177
+%subplot(2,1,2), dydxplot(P,Mdist);
178
+
179
+% to a project a new point x in the input space to the output space
180
+% do the following:
181
+
182
+% Diff = D - x(noc_x_1,:); Diff(find(isnan(Diff))) = 0; 
183
+% dx = sqrt((Diff.^2)*dim_x_1);
184
+% p = project_point(P,x,dx); % this function can be found from below
185
+% tlen = size(p,1);
186
+% plot(P(:,1),P(:,2),'bx',p(tlen,1),p(tlen,2),'ro',p(:,1),p(:,2),'r-')
187
+
188
+% similar trick can be made to the other direction
189
+
190
+
191
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
192
+%% subfunctions
193
+
194
+function vals = potency_curve(v0,vn,l)
195
+
196
+  % curve that decreases from v0 to vn with a rate that is 
197
+  % somewhere between linear and 1/t
198
+  vals = v0 * (vn/v0).^([0:(l-1)]/(l-1));
199
+
200
+
201
+function error = cca_error(P,Mdist,lambda)
202
+
203
+  [noc odim] = size(P);
204
+  noc_x_1 = ones(noc,1);
205
+  odim_x_1 = ones(odim,1);
206
+
207
+  error = 0;
208
+  for i=1:noc,
209
+    known = find(~isnan(Mdist(:,i)));
210
+    if ~isempty(known),   
211
+      y = P(i,:);                 
212
+      Dy = P(known,:) - y(noc_x_1(known),:);
213
+      dy = sqrt((Dy.^2)*odim_x_1);
214
+      fy = exp(-dy/lambda);
215
+      error = error + sum(((Mdist(known,i) - dy).^2).*fy);
216
+    end
217
+  end
218
+  error = error/2;
219
+
220
+
221
+function [] = dydxplot(P,Mdist)
222
+
223
+  [noc odim] = size(P);
224
+  noc_x_1 = ones(noc,1);
225
+  odim_x_1 = ones(odim,1);
226
+  Pdist = zeros(noc,noc);
227
+    
228
+  for i=1:noc,
229
+    y = P(i,:);                 
230
+    Dy = P - y(noc_x_1,:);
231
+    Pdist(:,i) = sqrt((Dy.^2)*odim_x_1);
232
+  end
233
+
234
+  Pdist = tril(Pdist,-1); 
235
+  inds = find(Pdist > 0); 
236
+  n = length(inds);
237
+  plot(Pdist(inds),Mdist(inds),'.');
238
+  xlabel('dy'), ylabel('dx')
239
+
240
+
241
+function p = project_point(P,x,dx)
242
+
243
+  [noc odim] = size(P);
244
+  noc_x_1 = ones(noc,1);
245
+  odim_x_1 = ones(odim,1);
246
+
247
+  % initial projection
248
+  [dummy,i] = min(dx);
249
+  y = P(i,:)+rand(1,odim)*norm(P(i,:))/20;
250
+ 
251
+  % lambda 
252
+  lambda = norm(std(P));
253
+
254
+  % termination
255
+  eps = 1e-3; i_max = noc*10;
256
+  
257
+  i=1; p(i,:) = y; 
258
+  ready = 0;
259
+  while ~ready,
260
+
261
+    % mutual distances
262
+    Dy = P - y(noc_x_1,:);        % differences in output space
263
+    dy = sqrt((Dy.^2)*odim_x_1);  % distances in output space
264
+    f = exp(-dy/lambda);
265
+  
266
+    fprintf(2,'iteration %d, error %g \r',i,sum(((dx - dy).^2).*f));
267
+
268
+    % all the other vectors push the projected one
269
+    fy = f .* (dx ./ dy - 1) / sum(f);
270
+  
271
+    % update    
272
+    step = - sum(fy(:,odim_x_1).*Dy);
273
+    y = y + step;
274
+  
275
+    i=i+1;
276
+    p(i,:) = y;   
277
+    ready = (norm(step)/norm(y) < eps | i > i_max);
278
+
279
+  end
280
+  fprintf(2,'\n');
281
+     
282
+  
0 283
\ No newline at end of file