1234.m

C-means一个简单完整的程序

function [m,g,F,pl,BetP,J,N] = ECM(x,K,version,alpha,beta,delta,init);

[n,d]=size(x);
if nargin<2
    error('ECM needs two arguments');
end

if nargin==2
version=0;
beta=2;alpha = 1;
delta = 100 ;  init=0;
end

if nargin==3
beta=2;alpha = 1;
delta = 100 ;  init=0;
end

if nargin==4
beta=2;
delta = 100 ;  init=0;
end

if nargin==5
delta = 100 ; init=0;
end

if nargin==6
 init=0;
end

if isempty(alpha) alpha=1;end
if isempty(beta) beta=2;end
if isempty(delta) delta=100;end
if isempty(version) version=0;end
if isempty(init) init=0;end
fprintf('-------------------------------------------\n');
fprintf('Evidential c-means\n');
fprintf('-------------------------------------------\n');
fprintf('Number of objects = %5i\n',n);
fprintf('alpha= %5.2f\n',alpha);
fprintf('beta = %5.2f\n',beta);
fprintf('delta = %5.2f\n',delta);

%--------------- construction of the focal set matrix  ---------
ii=1:2^K;
F=zeros(length(ii),K);
for i=1:K,
      F(:,i)=bitget(ii'-1,i);
end;


if init ==0
     fprintf('Random initialization of the centers\n');
g = randn(K,d);
else
    fprintf('Initialization of the centers using FCM\n');
    g = FCM(x,K);
end

pl=[];
BetP=[];
histJ=[];
%------------------------ iterations--------------------------------
fprintf('-----------------------------------------\n');
fprintf('Optimization\n');
pasfini=1;Jold = inf; 
while pasfini

gplus=[];
for i=2:nbFoc
    fi = F(i,:);
    truc = repmat(fi',1,d);
    gplus = [gplus;sum(g.*truc)./sum(truc)];
end

    % calculation of distances to centers 
    D=[];
    for j=1:nbFoc-1
        A = (x - ones(n,1)*gplus(j,:));
        B = diag(A*A');
        D = [D B];
    end
       
    % Calculation of masses
    
    m = zeros(n,nbFoc-1);
    for i=1:n
        for j=1:nbFoc-1
            vect1 = D(i,:);
            vect1 = ((D(i,j)*ones(1,nbFoc-1))./vect1).^(1/(beta-1));
            vect2 =  ((c(j)^(alpha/(beta-1)))*ones(1,nbFoc-1))./(c.^(alpha/(beta-1)));
            vect3 = vect1.*vect2 ;
            m(i,j)=1/(sum(vect3)+((c(j)^(alpha/(beta-1)))*D(i,j)/delta)^(1/(beta-1)));
        end
    end
   
   
        mi = repmat((c(indices).^(alpha-1)),n,1).*(m(:,indices).^beta);
        s = sum(mi,2);
        mats = repmat(s,1,d);
        xim = x.*mats ;
        blq = sum(xim);
        B=[B;blq];
    end
    

    g=A\B;
    
    mvide = ones(n,1)-sum(m,2) ; 
    J = sum(sum((m.^beta).*D.*repmat(c.^alpha,n,1)))+ delta*sum(mvide.^beta);
    fprintf('Objective function = %6.3f\n',J);
    pasfini =abs(J-Jold)>1e-3;
    Jold = J ; 
    histJ=[histJ;J];
end
%--------------------------- end of iterations ----------------------------
m = [ones(n,1)-sum(m,2) m]; 

if version==1
    if K>3
        mm=zeros(n,2^K);
        ii=1:2^K;
        F=zeros(length(ii),K);
        for i=1:K,
            F(:,i)=bitget(ii'-1,i);
        end;
        truc = sum(F,2);
        ind = find(truc<=2);ind=[ind;2^K];
        for j=1:length(ind)
            mm(:,ind(j))=m(:,j);
        end
    else
        mm=m;
    end
    P=[];
    for i=1:n
        pp = mtopl(mm(i,:));
        P=[P;pp];
        bet = betp(mm(i,:));
        BetP = [BetP;bet];
    end
    F=FF;
else
P=[];
for i=1:n
    pp = mtopl(m(i,:));
    P=[P;pp];
    bet = betp(m(i,:));
    BetP = [BetP;bet];
end
end

truc = sum(F,2);
singletons = find(truc==1);
pl = P(:,singletons);

% ----------- validity index ---------------------------
truc = truc';
mat = repmat(truc,n,1);
mat(:,1)=K*ones(n,1);
N = sum(sum(log(mat).*m))/log(K)/n;

fprintf('End of optimization\n');
fprintf('-------------------------------------------\n');


function [out] = betp(in) 
% computing BetP on ?from the m vector (in) 
% out = BetP vector: order a,b,c,... 
% beware: not optimalize, so can be slow for >10 atoms 
 
lm = length(in); 
natoms = round(log2(lm)); 		 
if 2^natoms == lm  
	if in(1) == 1 
		out = 	ones(1,natoms)./natoms; 
	else 
		betp = zeros(1,natoms); 
		for i = 2:lm 
			x = bitget(i-1,1:natoms); % x contains the 1 and 0 of i-1, for a,b,c... 
			betp = betp + in(i)/sum(x).*x; 
		end 
		out = betp./(1-in(1)); 
	end 
else 
	'ACCIDENT in betp: length of input vector not OK: should be a power of 2' 
end 


function [out] = mtopl(in) 
% computing FMT from m to pl 
% in = m vector 
% out = pl vector 
 
in = mtob(in); 
out = btopl(in); 

function [out] = mtob(in) 
% computing FMT from m to b.  
% in = m vector 
% out = b vector 
 
lm = length(in); 
natoms = round(log2(lm)); 		 
if 2^natoms == lm  
for step = 1:natoms  
	i124 = 2^(step-1); 			 
	i842 = 2^(natoms+1-step); 	 
	i421 = 2^(natoms - step); 	 
	in = reshape(in,i124,i842); 
	in(:,(1:i421)*2) = in(:,(1:i421)*2) + in(:,(1:i421)*2-1); 
end	 
out = reshape(in,1,lm); 
else 
	'ACCIDENT in mtob: length of input vector not OK: should be a power of 2' 
end 

function [out] = btopl(in) 
% compute pl from b 
% in = b 
% out = pl 
 
lm = length(in); 
natoms = round(log2(lm)); 		 
if 2^natoms == lm  
in = in(lm)-fliplr(in); 
in(1) = 0; 
out = in; 
else 
	'ACCIDENT in btopl: length of input vector not OK: should be a power of 2' 
end