bp2dim2u3a3a1.m

这是一个很好的关于BP神经网络的程序,用于两类数据的分类.

%%%% 2 classes BP algorithm in non-linear-classfing


% Produce sample set and original classfying line
% Produce samples: R0,R1
close all
clear all

% neuro number in each layers

N=2;        % Input vector dimension 

K=3;        % First layer neuros
J=3;        % Second layer neuros
I=1;        % Third layer neuros


M=10000;                         % Maxmum allowably ierative times
s=10;                              
S=8*s;                          % Nomber of total samples. Learning set, with total 2*S samples, is composed of both S samples class  
alf=0.1;                       % Step length
u0=0.3;                        % Sigmoid function f=1./(1+exp(-u/u0))
clr=['r  ','g  ','y  ','b  ','c  ','m  ','r: ','g: ','y: ','b: ','c: ','m: ','r-.','g-.','y-.','b-.','c-.','m-.'];

            figure(1)
            hold on
            whitebg(1,'k')
            figure(2)
            hold on
            whitebg(2,'k')
            
        % Sigmoid function
            figure(3)
            u=-10:0.1:10;
            plot(u,1./(1+exp(-u./u0)));
            clear u
            

        % Produce and plot sample set X0 and X1, learning samiple set R=R0+R1
            X0=randn(S,N);                        % S X0 input vector, in line vector with N dimension
            X0=[X0,ones(S,1)];                                  % Expand the pattern vector X0 with an 1

                R=4;               
                X11=randn(s,N)+[R*ones(s,1),R*ones(s,1)];           %
                X12=randn(s,N)+[-R*ones(s,1),R*ones(s,1)];          %  
                X13=randn(s,N)+[-R*ones(s,1),-R*ones(s,1)];         % 
                X14=randn(s,N)+[R*ones(s,1),-R*ones(s,1)];          % 
                X15=randn(s,N)+[R*ones(s,1),zeros(s,1)];            %
                X16=randn(s,N)+[zeros(s,1),R*ones(s,1)];            %  
                X17=randn(s,N)+[-R*ones(s,1),zeros(s,1)];           % 
                X18=randn(s,N)+[zeros(s,1),-R*ones(s,1)];           % X1 composed of 8 part, each in a scale of s and 3*s=S

                X1=[X11;X12;X13;X14;X15;X16;X17;X18];                                   %
                X1=[X1,ones(S,1)];                                  % Expand the pattern vector X1 with an 1

            D0=0.9*ones(S,1);               % Expected output for X0; 
            D1=0.1*ones(S,1);               % Expected output for X1;
            
            mX=[X0;X1];             % Sample matrix with dimension S*(N+1)
            D=[D0;D1];              % Expected output matrix with dimension of S*I
            
            figure(1)
            plot(X0(:,1),X0(:,2),'b*')
            plot(X1(:,1),X1(:,2),'gx')

        % Plot axies
                        Xmax=max(mX(:,1));
                        Xmin=min(mX(:,1));
                        Ymax=max(mX(:,2));
                        Ymin=min(mX(:,2));
                        Xmax=ceil(Xmax);
                        Xmin=floor(Xmin);
                        Ymax=ceil(Ymax);
                        Ymin=floor(Ymin);
            plot([Xmin,Xmax],[0,0],'m-.')
            plot([0,0],[Ymin,Ymax],'m-.')
            text(Xmin+0.3,Ymax-0.5,'X0: *');
            text(Xmin+0.3,Ymax-0.9,'X1: x');
            
                        
            clear X11 X12 X13 nn X0 X1 D0 D1 s
            
   aw=1;
   Vn=0;    % Times
   
   disp([0,fix(clock)])
  
   while aw~=0,
        aw=input('Continue(1) or stop (0)?')
        if aw==0,
            break;
        end
            
        %Produce original matrixes W3,W2,W1
            W1=rand(K,N+1);
            W2=rand(J,K+1);
            W3=rand(I,J+1);
            p=1;    % Generations
            
            while p<M ,
                Ep(p)=0;                                            % Square error
                dW1=0;dW2=0;dW3=0;                                  % Weights correctors
                
                    for s=1:2*S
 
                        X=mX(s,:);                          % Take a sample
                        
                        % out of first layer
                        Net1=X*W1';                         % Pure input, K dimension vector;
                        u=exp(-Net1./u0);                   % K'd
                        O1=1./(1+u);                        % Output vectors of layer 1
                        O1=[O1,1];                          % Expand the output vectors
                        F1=diag(u./(u0*(1+u).^2));          % Differential coefficient, K by K dimension matrix
                
                        % out of second layer
                        Net2=O1*W2';                        % Pure input, J dimension vector;
                        u=exp(-Net2./u0);                   % J'd
                        O2=1./(1+u);                        % Output vectors of layer 2
                        O2=[O2,1];                          % Expand the output vectors
                        F2=diag(u./(u0*(1+u).^2));          % Differential coefficient, J by J dimension matrix
                
                        % out of second layer
                        Net3=O2*W3';                        % Pure input, I dimension vector;
                        u=exp(-Net3./u0);                   % I'd
                        O3=1./(1+u);                        % Output vectors of layer 3
                        F3=diag(u./(u0*(1+u).^2));          % Differential coefficient, I by I dimension matrix

                        % calculate errors and amend weights
                        Es=D(s)-O3;                         % Errors
                        Ep(p)=Ep(p)+Es*Es';                 % Summation of square error
                        
                        d3=Es*F3;                           % Error component for layer 3
                        d2=d3*W3(:,1:J)*F2;                 % Error component for layer 2
                        d1=d2*W2(:,1:K)*F1;                 % Error component for layer 1
                      
                        dW3=dW3+alf*d3'*O2; dW2=dW2+alf*d2'*O1;  dW1=dW1+alf*d1'*X;                    %
                        
                    end % of for
                        W1=W1+dW1/(2*S); W2=W2+dW2/(2*S); W3=W3+dW3/(2*S);    % Amending with mean correctors
                        
                        p=p+1;
                        %if mod(p,100)==0,
                        %    disp(p)
                        %end
                
             end % while
             
                % Plot error evaluation line
                figure(2)                   % Plot error evaluation line
                cln=3*(mod(Vn,18))+1;       % 
                clrnum=clr(cln:cln+2);      % Color number 'r  ' or 'y: ' or 'g-.' etc. each with 3 characters
                plot(Ep,clrnum);
                clear p s u dW1 dW2 cln
                
                % Plot the borderline of 2 classes
                figure(1)
                n=1; Lx=[];Ly=[];
                m=1; Rx=[];Ry=[];
                for y=Ymin:0.1:Ymax
                    NL='New Line';               % A new line start
                    for x=Xmin:0.1:Xmax
                        X=[x,y,1];
                        O1=[1./(1+exp(-X*W1'./u0)),1];
                        O2=[1./(1+exp(-O1*W2'./u0)),1];
                        O3=1./(1+exp(-O2*W3'./u0));
                        
                        if NL~='New Line'    
                            if (pr<=0.5)&(O3>=0.5),
                                Lx(n)=x;
                                Ly(n)=y;
                                n=n+1;
                            end
                            if (pr>=0.5)&(O3<=0.5),
                                Rx(m)=x;
                                Ry(m)=y;
                                m=m+1;
                            end
                        end
                        pr=O3; NL='0';          % Storage last value
                    end
                end
                
                % Plot borderline
                if size(Lx)>0,
                    plot(Lx,Ly,clrnum);             % Plot borderline
                    text(Lx(1),Ly(1),num2str(Vn));                          % 
                    text(Lx(fix(length(Lx)/2)),Ly(fix(length(Lx)/2)),num2str(Vn));    % Write bordline number
                    text(Lx(length(Lx)),Ly(length(Lx)),num2str(Vn));        % Repeat
                end
                if size(Rx)>0,    
                    plot(Rx,Ry,clrnum);             % Plot borderline
                    text(Rx(1),Ry(1),num2str(Vn));                          % Write bordline number
                    text(Rx(fix(length(Rx)/2)),Ry(fix(length(Rx)/2)),num2str(Vn));    % Write bordline number
                    text(Rx(length(Rx)),Ry(length(Rx)),num2str(Vn));        % Repeat 
                end
                
                Vn=Vn+1;
                clear  Lx Ly Rx Ry Ep NL pr
                %aw=aw+1
                disp([aw,fix(clock)])
     end % while aw