ga1.m

多元线性回归：MATLAB源程序多元线性回归

%#									
%#  function [pv,fv,pop]=ga1(xc,yc,xt,yt,cc,chrom,maxvar,psel,...	
%#			       mate,pcros,pmut,gen,cvgp,valeurlim,...	
%#			       backf,figf);				
%#									
%#  AIM:        Main routine of the GA.					
%#		Subroutine of ga0.m.				
%#		Performs the selection of variables with xc and yc.	
%#		Predicts xt and yt with the output subsets.		
%#									
%#  PRINCIPLE:  Based on the flow diagram (Figure 1) in			
%#		Analytical Chemistry, 67 (1996), 4295.			
%#									
%#  INPUT:      ALL INPUTS ARE EXPLAINED IN THE ALGORITHM IN GA0.M	
%#									
%#	xc, yc	: (row x col) Calibration set, on which the GA will run	
%#	xt, yt	: Test set, on which the GA is validated		
%#		(x: matrix of independent variables (object x variable)	
%#	 	 y: dependent vector)					
%#	cc	: Number of outer cycles, to perform the external	
%#	          validation						
%#	chrom	: Number of strings subsets of variables) in the 	
%#		  population						
%#	maxvar	: Maximum number of variables for selection		
%#	psel	: Probability of selection in the initiation		
%#	      	  (psel=maxvar/col)					
%#	mate	: Number of matings at each generation			
%#	pcros	: Probability of performing cross-over (advised: 0.5)	
%#	pmut	: Probability of performing mutation (advised: 1/col)	
%#	gen	: Number of generations in each cycle			
%#	cvgp	: Number of deletion groups for cross-validation	
%#	valeurlim: Minimal acceptable response (1/RMSECV)		
%#	backf	: Frequency for backward elimination			
%#	figf	: Frequency for plotting the response			
%#									
%#  OUTPUT:								
%#	pv	: (chrom x 1) Vector containing the cross-validated	
%#		  response of the corresponding strings in the
%#		  population						
%#	fv	: (chrom x 1) Vector contained the response, after	
%#		  prediction of the "test set", with the corresponding	
%#		  string in the population				
%#	pop	: (chrom x col) Subsets of variables selected in each	
%#		  string (on the rows). A binary coding is used:	
%#		    0: the corresponding variable is not selected	
%#		    1: the corresponding variable is selected		
%#		  To obtain directly the subset number i:		
%#		    subset_i=find(pop(i,:)==1)				
%#  SUBROUTINES:							
%#	checkson.m : Checking of the children strings			
%#	checksub.m : Comparison of the children strings to other strings
%#		     of the population, before and after entering them	
%#		     in this population.				
%#	gainit.m   : Random initiation of the strings at the beginning	
%#	   	     of the GA						
%#	gamut.m    : Mutation on each child string			
%#	mlr.m      : Computation of the MLR model			
%#	mlrcv.m    : MLR model with cross-validation to compute the	
%#	             the prediction error				
%#	rms.m	   : Calculates the RMSEP between observed and		
%#		     predicted Y vectors.			
%#	stepwise.m : Backward elimination				
%#									
%#  AUTHOR:     Delphine Jouan-Rimbaud					
%#              Copyright(c) 1997 for ChemoAC				
%#              FABI, Vrije Universiteit Brussel		
%#              Laarbeeklaan 103 1090 Jette				
%#									
%# VERSION: 1.1 (28/02/1998)						
%#									
%#  TEST:	Frederic Despagne					
%#									

function [pv,fv,pop]=ga1(xc,yc,xt,yt,cc,chrom,maxvar,psel,mate,...
			 pcros,pmut,gen,cvgp,valeurlim,backf,figf);
clc

[row,col]=size(xc);
chromstep=[];
minrep=[0];	% minimal response of the population
maxrep=[0];	% maximal response of the population
meanrep=[0];	% mean response of the population

% *** 1) Initiation of population ***

[pop,rep]=gainit(xc,yc,chrom,maxvar,psel,cvgp);

  % pop : Initial population, containing 'chrom' strings
  % each string having a maximum of maxvar variables.
  % rep : Vector with the response of each string
%********************************************************


% *** 2) Start the GA for "gen" generations ***

for gene=1:gen	% gene: generation number

home, disp(['cycle ',int2str(cc),', generation ',int2str(gene)])

% 2.1) Plot of the minimal, mean, and maximal response in
% the actual population

if rem(gene,figf)==0 % a plot should be made at generation gene
   minrep=[minrep,min(rep)];	% updating of minrep
   maxrep=[maxrep,max(rep)];	% updating of maxrep
   meanrep=[meanrep,mean(rep)];	% updating of meanrep

   figure(1)
   set(1,'units','normalized','position',[0 0 1 0.8],...
   'menubar','none','name','Response',...
   'numbertitle','off')
   plot(0:figf:gene,minrep,'r',0:figf:gene,maxrep,'c',0:figf:gene,meanrep,'g');
   xlabel('Number of generations')
   ylabel('Cross-validated evaluation function')
   title('Minimal (red), mean (green) and maximal (blue) response (1/RMSECV)');
   figure(gcf)
end


% 2.2) Selection - Cross-over - Mutation

prob=rep/sum(rep); 
sumprob=cumsum(prob);

% The probability of selection of a string depends on its response


for mating=1:mate	% two parents strings are selected "mate" times. 
   a=rand;		% random number between 0 and 1
   [b,index]=sort([a;sumprob]);
   p1=find(index==1);	% p1 : index of the first selected parent
			% corresponds to the position of "a"
   p2=p1;
   while p2==p1		% the two parents strings must be different
	a=rand;
	[b,index]=sort([a;sumprob]);
	p2=find(index==1);	% p2 : index of the second selected parent
   end


   parent1=pop(p1,:);
   parent2=pop(p2,:);
   clear a b index p1 p2

% Cross-over
   child=[];		% matrix with children strings on the rows
   childrep=[];		% response associated to the children strings
   diff=parent1-parent2;
   a=find(diff~=0); 	% The cross-over is done on bits where the value 
			% in each string is different
   son1=parent1;
   son2=parent2;
   for i=1:length(a)
	r=rand;		% Random number between 0 and 1
	if r<=pcros	% Swapping of bits
	   son1(a(i))=parent2(a(i));
	   son2(a(i))=parent1(a(i));
	end
   end
   clear a r i

% Mutation
   son=[son1;son2];
   newson=gamut(son,pmut);	% mutations done in GAMUT.M
				% newson: matrix containing the sons after mutation
				% newson.
   clear son1 son2 son
%*************************************************************************

% 2.3) Check of children strings, and evaluation

   [newson]=checkson(newson,pop,chrom,xc,yc,cvgp,maxvar);

% CHECKSON.M checks that the children strings obey the
% constraints of GA. The output "newson" contains the children strings
% obeying the GA constraints.

   for nbson=1:size(newson,1)	% Evaluation of each child string
	[rr]=mlrcv(xc(:,newson(nbson,:)),yc,1,cvgp);
				% rr=RMSECV of the tested new son
	evson(nbson,1)=1/rr;	% fitness of the new son = 1/RMSECV
   end

   child=[child;newson];	% updating of child
   childrep=[childrep;evson];	% updating of childrep

end %for mating=1:mate

% After 'mate' matings, the matrix child contains a maximum
% of 2 x mate children strings (those which obey the constraints
% of GA), and the associated responses are in the vector childrep.


for lch=1:size(child,1)

   [pop,rep]=checksub(child(lch,:),childrep(lch),pop,rep,chrom);

% CHECKSUB.M checks whether the children strings can enter
% the population, and replace worse solutions.

end

repzero=find(rep==0);
% In order to explore other regions of the search space,
% all the null strings of the output population of checksub
% are now replaced by randomly generated solutions.
% (modification of Leardi's original algorithm).

if repzero~=[]
   while repzero~=[]
	le=length(repzero);	% Number of null vectors in the population
	compteur=1;
	while compteur <= le
	   newstring=(rand(1,col)<=psel); % Random generation of a new string

	   if sum(newstring)>0 & sum(newstring)<=maxvar
% The string meets the requirements to enter the population

		[rr]=mlrcv(xc(:,newstring),yc,1,cvgp);	% evaluation of the string
		newrep=1/rr;

		[pop,rep]=checksub(newstring,newrep,pop,rep,chrom);
% check that the new string can replace a worse solution in the population
		compteur=compteur+1;
	   end
	end
	repzero=find(rep==0);
   end
end
clear newstring newrep 
%*************************************************************************


% Backward elimination
% performed either on the best string, or on the second, third,... best string,
% if better strings, or very similar strings, have already undergone once a backward elimination.

if backf~=0			% backward elimination is allowed in this GA
   if rem(gene,backf)==0,	% A backward elimination should be performed at this actual generation
	if gene==backf		% First cycle of backward elimination
	   chrom1=pop(1,:);
	   rc1=rep(1);
	   chromstep=[chromstep;chrom1];	% all the strings which have gone once through the
						% backward elimination are placed into a matrix
						% chromstep; they will not undergo another backward
						% elimination later, neither will a string too similar too
						% one of them.
	end;

	if gene~=backf				% If not first cycle of backward elimination,
						% one checks whether the best string (or a very similar
						% one) has not already undergone the backward elimination.
						% Otherwise the second best string is candidate to backward
						% elimination, and is also checked in turn, and so on...
	   i=1; 
	   while i<=chrom
		pp=pop(i,:);
		j=1;
		jj=0;	% initialisation of jj parameter
			% jj=0 if string i can undergo the backward elimination.
			% jj=1 if string i cannot undergo the backward elimination
		while j<= size(chromstep,1)
		   compar=chromstep(j,:)+pp;
		   if length(find(compar==1))>=round(0.05*col),

			% in the vector find(compar==1), all elements equal to 1
			% correspond to the variables present in one of the two
			% strings only.
			% Two strings are considered very similar if fewer
			% than 10 % of the selected variables are different.

			j=j+1;	% string j of chromstep is NOT similar to string i of pop.
				% test string j+1 of chromstep.
		   else
			i=i+1;	% string j of chromstep IS similar to string i of pop;
				% therefore, string i of pop cannot undergo the backward elimination.
			j=size(chromstep,1)+2;	% end of while j loop
			jj=1;			% string i cannot undergo the backward elimination
		   end
		end
		if jj==0	% string i of pop can undergo backward elimination
		   indchrom=i;
		   chrom1=pp;
		   rc1=rep(i);
		   chromstep=[chromstep;chrom1];
		   i=chrom+1;			% end of the while i loop
		   j=size(chromstep,1)+2;	% end of the while j loop
		end                                     
	   end %while
	end %if gene~=

	% Once the best candidate to backward elimination was found
	% the backward elimination can be performed.
	% This is done in STEPWISE.M

	[pop,rep]=stepwise(chrom1,rc1,pop,rep,xc,yc,cvgp,chrom,maxvar);
   end %if find


% At the end of STEPWISE.M, CHECKSUB.M is performed, so
% some null vectors may be present in the new 
% population (pop); these should be replaced by randomly
% generated solutions.

repzero=find(rep==0);if repzero~=[]
   while repzero~=[]
	le=length(repzero);			% Number of null strings in the population
	compteur=1;ccc=0;
	while compteur <= le
	   newstring=(rand(1,col)<=psel);	% Random generation of a new string
	   ccc=ccc+1;
	   if sum(newstring)>0 & sum(newstring)<=maxvar
	   % the string meets the requirements to enter the population

		[rr]=mlrcv(xc(:,newstring),yc,1,cvgp);
		newrep=1/rr;
		%evaluation of the new string

		[pop,rep]=checksub(newstring,newrep,pop,rep,chrom);
		% check that the new string can replace a worse solution in the population

		compteur=compteur+1;
	   end
	end
	repzero=find(rep==0);
   end
end
clear newstring newrep 

end %if backf

end     % for gene=1:max_gene

% Now, all generations have been performed. The strings
% which are still in the matrix pop are the final solutions
% found by the GA. Each of these solutions is tested on
% the "external" test set (xt, yt).

pv=rep; 
for i=1:chrom
   [b0,b]=mlr(xc(:,pop(i,:)),yc,1);
   ytp=b0+xt(:,pop(i,:))*b;
   rept=rms(yt,ytp);
   fv(i,1)=1/rept;
end