mykmeans.m 1.56 KB
Newer Older
Stéphane Adjemian's avatar
Stéphane Adjemian committed
1
function [c,SqrtVariance,Weights] = mykmeans(x,g,init,cod)
2

Stéphane Adjemian's avatar
Stéphane Adjemian committed
3
% Copyright (C) 2013-2017 Dynare Team
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

[n,m] = size(x) ;
indold = zeros(1,m) ;
if cod==0
Stéphane Adjemian's avatar
Stéphane Adjemian committed
23
24
25
26
    d = transpose(sum(bsxfun(@power,bsxfun(@minus,x,mean(x)),2)));
    d = sortrows( [transpose(1:m) d],2) ;
    d = d((1+(0:1:g-1))*m/g,1) ;
    c = x(:,d);
27
else
Stéphane Adjemian's avatar
Stéphane Adjemian committed
28
    c = init ;
29
end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
for iter=1:300
    dist = zeros(g,m) ;
    for i=1:g
        dist(i,:) = sum(bsxfun(@power,bsxfun(@minus,x,c(:,i)),2));
    end
    [rien,ind] = min(dist) ;
    if isequal(ind,indold)
        break ;
    end
    indold = ind ;
    for i=1:g
        lin = bsxfun(@eq,ind,i.*ones(1,m)) ;
        h = x(:,lin) ;
        c(:,i) = mean(h,2) ;
    end
end
SqrtVariance = zeros(n,n,g) ;
Weights = zeros(1,g) ;
48
for i=1:g
Stéphane Adjemian's avatar
Stéphane Adjemian committed
49
50
51
52
    temp = x(:,bsxfun(@eq,ind,i*ones(1,m))) ;
    u = bsxfun(@minus,temp,mean(temp,2)); %temp-mean(temp,1)' ;
    SqrtVariance(:,:,i) = chol( (u*u')/size(temp,2) )' ;
    Weights(i) = size(temp,2)/m ;
53
end