demo_anderson.m

很好的matlab模式识别工具箱

   maxerr = model.err;
  
%   [nalpha,ntheta,solution,nt,alpha1,alpha2]=...
%       eanders(sets.MI,sets.SIGMA,sets.I,iter,precision/100,t,alpha1,alpha2);
%    if sum(nalpha)==0,
%      solution=-1;
%      nalpha=alpha;
%      ntheta=theta;
%   end
% text=sprintf('Step t=%d\nLine [%f , %f]*x=%f',nt,nalpha(1),nalpha(2),ntheta);
  text=sprintf(['Interation(s) t=%d\nLinear rule q(x)=sgn([%f , %f]*x %+f)\n'...
                 'Minimal r = %.8f\nClassification error = %.4f%%'],...
      nt,nalpha(1),nalpha(2),-ntheta,minr,maxerr*100);
case 2
   % Original Anderson`s solution
   distr.Mean = sets.MI;
   distr.Cov = reshape(sets.SIGMA,2,2,2);

   options.eps = precision;
   options.tmax = t+iter;
   init_model.gamma=lambda;
   init_model.t = t;

   model = androrig( distr, options, init_model);

   nalpha = model.W;
   ntheta = -model.b;
   solution = model.exitflag;
   nt = model.t;
   lambda = model.gamma;
   maxerr = model.err;
   minr = min([model.r1 model.r2]);
  
%   [nalpha,ntheta,solution,nt,lambda,ni,maxerr]=...
%      oanders(sets.MI,sets.SIGMA,sets.I,iter,precision,t,lambda);

  text=sprintf(['Interation(s) t=%d\nLinear rule q(x)=sgn([%f , %f]*x %+f)\n'...
                 'Minimal r = %.8f\nClassification error = %.4f%%'],...
      nt,nalpha(1),nalpha(2),-ntheta,minr,maxerr*100);

%  text=sprintf(...
%      'Step t=%d\nLine [%f , %f]*x=%f\nNi = %f, (1-Lambda)/Lambda = %f, Max error = %f%%',...
%      nt,nalpha(1),nalpha(2),ntheta,ni,(1-lambda)/lambda,maxerr*100);
case 3
   distr = sets;
   options.eps = precision;
   options.tmax = t+iter;
   init_model.W=alpha;
   init_model.b=-theta;
   init_model.t = t;

   % e-Optimal solution
   model = ggradandr( distr, options, init_model );

   nalpha = model.W;
   ntheta = -model.b;
   solution = model.exitflag;
   minr = model.r;
   nt = model.t;
   maxerr = model.err;
  
%   [nalpha,ntheta,solution,nt,alpha1,alpha2]=...
%       eanders(sets.MI,sets.SIGMA,sets.I,iter,precision/100,t,alpha1,alpha2);
%    if sum(nalpha)==0,
%      solution=-1;
%      nalpha=alpha;
%      ntheta=theta;
%   end
% text=sprintf('Step t=%d\nLine [%f , %f]*x=%f',nt,nalpha(1),nalpha(2),ntheta);
  text=sprintf(['Interation(s) t=%d\nLinear rule q(x)=sgn([%f , %f]*x %+f)\n'...
                 'Minimal r = %.8f\nClassification error = %.4f%%'],...
      nt,nalpha(1),nalpha(2),-ntheta,minr,maxerr*100);
 
end

if solution==-1,
    text=sprintf('Solution does not exist.\n');
    play=-1;
    return;
elseif solution==1,
   text=strvcat(text,sprintf('Solution was found in %d step(s)',nt));
   play=0;
else
   play=1;
end

% store new values
h.line.t = nt;
h.line.alpha = nalpha;
h.line.alpha1 = alpha1;
h.line.alpha2 = alpha2;
h.line.lambda = lambda;
h.line.theta = ntheta;

return

%==============================================================
function [handler]=pandr2d(MI,SIGMA,I,alpha1,theta1,handler,anim,alpha2,theta2)
% PANDR2D displays solution of Generalized Anderson's task in 2D.
% [handler]=pandr2d(MI,SIGMA,I,alpha1,theta1,handler,anim,alpha2,theta2)
%
% PANDR2D plots given solution of the Generalized Anderson`s task (GAT) in
%  2-dimensional feature space. This function plots separation line 
%  (2D case of the GAT) and input classes which are symbolized by sets 
%  of ellipsoids.
%
%  Input arguments MI, SIGMA and I describe input mixture of normal 
%  distributions. The pair of input arguments {alpha,theta} describes 
%  separation line which is particular solution of the GAT. 
%
%  For information on the meaning of the arguments refer to help of 
%  functions solving the GAT (ganders,ganders2,eanders etc.).
%
%  When the quadruple of input parameters handler, anim, alpha2 and theta2
%  enter function then a change of the solution is depicted too. The pair
%  alpha1, theta1 is old solution and alpha2, theta2 is new solution. The
%  argument handler contains information about graphics abjects used 
%  in the last call of the function and returned in output variable handle. 
%  When the argument anim=1 then the change is animated.
%
% See also PANDR2D, DEMO_ANDERSON, GANDERS, GANDERS2, EANDERS, OANDERS.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 23.11.1999, 11.5.2000
% Modifications
% 24. 6.00 V. Hlavac, comments polished.
% 2.8.00 V.Franc, comments changed

% constants
BORDER=0.95;
POINT_SIZE=8;
POINT_COLOR='k';
POINT_WIDTH=2;
LINE_WIDTH=1;
LINE_COLOR='k';
ELLIPSE_WIDTH=1;
INTERPOL=50;
ANIM_DIF=20;


% handles input arguments
if nargin < 7,
   anim=0;
end
if nargin < 6,
   handler=-1;
end

%%handler


alpha1=alpha1(:);     % alpha1 will column

% number of ellipses
N=size(MI,2);
% dimension, must be equal to 2
DIM=size(MI,1);

if handler==-1,
   % proportions of a separation line
   window=axis*BORDER;
else
   % proportions of a separation line
   window=getaxis(handler(2*N+3))*BORDER;
end;

if anim==0,

   % computes minimal distance among the line alpha*x=theta and
   % the elipses (x-MI)'*inv(SIGMA)*(x-MI)
    [R,inx]=min( abs(alpha1'*MI-theta1)...
        ./sqrt( reshape(alpha1'*SIGMA,DIM,N)'*alpha1 )' );


   if handler==-1,
      % first painting

      % ellipses
      for i=1:N,
         mi=MI(:,i);
         sigma=SIGMA(:,(i-1)*DIM+1:i*DIM);
%%%         isg=inv(sigma);

%%%         [x,y]=ellipse(isg,INTERPOL,R,mi);
         [x,y]=ellips(mi,inv(sigma),R,INTERPOL);

         handler(i)=plot(x,y,'Color',marker_color(I(i)),...
            'EraseMode','xor',...
            'LineWidth',ELLIPSE_WIDTH,...
            'UserData',sigma);
      end
      % line
      [x1,y1,x2,y2]=cliplin1(alpha1,theta1,window);
      handler(N*2+1)=line([x1 x2],[y1 y2],...
         'LineWidth',LINE_WIDTH,...
         'Color',LINE_COLOR,...
         'EraseMode','xor');

%%      get(handler(N*2+1))

      % pull point
      mi=MI(:,inx);
      sigma=SIGMA(:,(inx-1)*DIM+1:inx*DIM);
      x0=mi-(alpha1'*mi-theta1)*sigma*alpha1/(alpha1'*sigma*alpha1);
      handler(N*2+2)=line(x0(1),x0(2),...
         'LineStyle','none',...
         'Color',POINT_COLOR,...
         'MarkerSize',POINT_SIZE,...
         'LineWidth',POINT_WIDTH,...
         'Marker','x',...
         'EraseMode','xor');

       handler(N*2+3)=gca;

   else % if handler==-1,
      % at least second painting

      % hide all objects
%      set(handler,'Visible','off');

      % ellipses
      for i=1:N,
         mi=MI(:,i);
         sigma=get(handler(i),'UserData');
         [x,y]=ellips(mi,inv(sigma),R,INTERPOL);
         set(handler(i),'XData',x,'YData',y,'Visible','on');
      end

      % line
      [x1,y1,x2,y2]=cliplin1(alpha1,theta1,window);
%%      N*2+1
%      get(handler(N*2+1))
%      handler
%      N
 %     x1
%      x2
      set(handler(N*2+1),'XData',[x1 x2],'YData',[y1 y2],'Visible','on');

      % pull point
      mi=MI(:,inx);
      sigma=SIGMA(:,(inx-1)*2+1:inx*2);
      x0=mi-(alpha1'*mi-theta1)*sigma*alpha1/(alpha1'*sigma*alpha1);
      set(handler(N*2+2),'XData',x0(1),'YData',x0(2),'Visible','on');

   end % if handler==-1,

else % if anim==0,

   alpha2=alpha2(:);   % alpha2 will column

   % computes number of the animation steps
   [x11,y11,x12,y12]=cliplin1(alpha1,theta1,window);
   [x21,y21,x22,y22]=cliplin1(alpha2,theta2,window);

   difa=max([sqrt((x11-x21)^2+(y11-y21)^2),sqrt((x11-x21)^2+(y11-y21)^2)]);
   difw=min([window(4)-window(3), window(2)-window(1)]);
   nsteps=max([1,ceil(ANIM_DIF*(difa/difw))]);

   % play
   for i=1:nsteps,
      k=i/nsteps;
      alpha=(1-k)*alpha2+k*alpha1;  % smooth transition of alpha
      theta=(1-k)*theta2+k*theta1;  % --//-- theta

      % recursion
      handler=pandr2d(MI,SIGMA,I,alpha,theta,handler,0);
   end

end % if anim==0,


%=========================================================
function [x1,y1,x2,y2,in]=cliplin1(alpha,theta,window)
% [x1,y1,x2,y2,in]=cliplin1(alpha,theta,window)
%
% CLIPLIN1 clips the line given by the equation alpha*x=theta along
%   the window. It returns two points on the border of the window.
%   If the line is in the window then the argument is equal to 1
%   else it returns 0.
%

minx=window(1);
maxx=window(2);
miny=window(3);
maxy=window(4);

x=zeros(4,1);
y=zeros(4,1);

if alpha(1)==0,
   if alpha(2)~=0,
      x1=minx;
      y1=theta/alpha(2);
      x2=maxx;
      y2=y1;
      in=1;
   else
      % if alpha == 0 then it means the bad input.
      x1=0;
      y1=0;
      x2=0;
      y2=0;
      in=0;
   end
elseif alpha(2)==0,
   x1=theta/alpha(1);
   y1=miny;
   x2=x1;
   y2=maxy;
   in=1;
else
   y(1)=maxy;
   x(1)=(theta-alpha(2)*y(1))/alpha(1);
   y(2)=miny;
   x(2)=(theta-alpha(2)*y(2))/alpha(1);

   x(3)=maxx;
   y(3)=(theta-alpha(1)*x(3))/alpha(2);
   x(4)=minx;
   y(4)=(theta-alpha(1)*x(4))/alpha(2);

   j=0;
   for i=1:4,
      if x(i) <= maxx & x(i) >= minx & y(i) <= maxy & y(i) >= miny,
         if j==0,
            j=j+1;
            x1=x(i);
            y1=y(i);
         elseif j==1,
            j=j+1;
            x2=x(i);
            y2=y(i);
         end
      end
   end

   if j<2,
      x1=0;
      y1=0;
      x2=0;
      y2=0;
      in=0;
   else
      in=1;
   end
end % elseif alpha(2)==0

function []=clrchild(handle)

delete(get(handle,'Children'));

return;

function [win]=cmpwin(mins,maxs,xborder,yborder)

dx=max( (maxs(1)-mins(1)), 1 )*xborder;
dy=max( (maxs(2)-mins(2)), 1 )*yborder;
x1=(mins(1)-dx);
x2=(maxs(1)+dx);
y1=(mins(2)-dy);
y2=(maxs(2)+dx);
win=[x1 x2 y1 y2];

%========================================
function [rect]=getaxis(handle)

rect=[get(handle,'XLim'),get(handle,'YLim'),get(handle,'ZLim')];

return;

function []=setaxis(handle,rect)

set(handle,'XLim',rect(1:2));
set(handle,'YLim',rect(3:4));

if size(rect,2)>=6,
   set(handle,'ZLim',rect(5:6));
end

return;