nnmodel.m

人工神经网络：MATLAB源程序用于训练测试

	w2trunc([i mdef+1]) = w2ref([i mdef+1]);	% are set to 0

	% Estimation of training responses using only one hidden node in the model
	[ytrunch1,ytrunch2] = lmeval(topo,w1ref,w2trunc,xtr);

	% Inverse-scaling of responses estimated using one hidden node only
	ytrunchhat = invrange(ytrunch2,0.2,0.8,tableout);
	partyh(i,:) = ytrunchhat;

	% Determination of variance of the predicted responses
	varytrunch0(i) = cov(ytrunchhat);	
end

% The sensitivity of each hidden node is the scaled variance of the responses
% estimated using this variable only
varytrunch = 100*(varytrunch0./sum(varytrunch0));


%%%%%%%%%% INFLUENCE OF REMOVAL OF EACH HIDDEN NODE ON TRAINING RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:mdef
	w2trunc = w2ref;		
	w2trunc(i) = 0;	

	% Estimation of training responses when hidden node i is set to 0
	[ytrunch1,ytrunch2] = lmeval(topo,w1ref,w2trunc,xtr);

	% Inverse-scaling of responses estimated when hidden node i is set to 0 
	ytrunchhat = invrange(ytrunch2,0.2,0.8,tableout);
	partyhh(i,:) = ytrunchhat;	
end
		

%%%%%%%%%%%%%%%%%%%%%%%%%%% MENU FOR MODEL INSPECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contflag = 1;			% Flag for menu display	
while contflag			% While summary results and figures are requested, display menu
	men1 = menu('Model inspection plots',...
	'1. Prediction results',...
	'2. Input variables sensitivity',...
	'3. Hidden nodes sensitivity',...
	'4. Partial models - Input variables',...
	'5. Partial models - Hidden nodes',...
	'6. Projection on hidden nodes',... 
	'7. Deviations from linearity',...
	'8. Summary results',...
	'9. Draw the neural network',...
	'10. End');
		
	if men1 == 1	% Prediction results
		continue1 == 2;
		clc;close
		while continue1
			men2 = menu('Summary results',...
			'1. Training set',...
			'2. Monitoring set',...
			'3. Test set',...
			'4. Back to main menu');
			if men2 == 1		% Training results
				clc;close
				rtra = corrplot(ytrain,ytrref',1);	% Correlation plot, training set
				xlabel('Observed values - Training set')
				set(gcf,'Position',setfig);
				men3 = menu('Training results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close
				plot(1:nxtr,ytrain,'+',1:nxtr,ytrref,'co',1:nxtr,ytrain,'y',1:nxtr,ytrref,'c--');
				renstra = sprintf('RMSEC=%9.6f',rmstraref);
				xlabel(renstra); ylabel('+ Experimental     o Predicted');title('Training set, adjustment')
				set(gcf,'Position',setfig);
				men3 = menu('Training results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close;
				plot(1:nxtr,residtrain,1:nxtr,0*residtrain,'c--',1:nxtr,residtrain,'yo')
				xlabel('Sample number'); ylabel('Residuals');title('Training set, residuals')
				set(gcf,'Position',setfig);
				men3 = menu('Training results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close;
				hist(rmstra);
				xlabel('RMSEC values'); ylabel('Cumulative frequencies');title('Training set')
				set(gcf,'Position',setfig);
				men3 = menu('Training results',...
				'1. Print',...
				'2. Back to previous menu');
				if men3 == 1
					print
				end;close
			elseif men2 == 2	% Monitoring results		
				clc;close
				rmon = corrplot(ymon,ymonref',1);	% Correlation plot, monitoring set
				xlabel('Observed values - Monitoring set')
				set(gcf,'Position',setfig);
				men3 = menu('Monitoring results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close
				plot(1:nxm,ymon,'+',1:nxm,ymonref,'co',1:nxm,ymon,'y',1:nxm,ymonref,'c--')		
				rensmon = sprintf('RMSEM=%9.6f',rmsmonref);
				xlabel(rensmon); ylabel('+ Experimental     o Predicted');title('Monitoring set, adjustment')
				set(gcf,'Position',setfig);
				men3 = menu('Monitoring results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close
				plot(1:nxm,residmon,1:nxm,0*residmon,'c--',1:nxm,residmon,'yo')	
				xlabel('Sample number'); ylabel('Residuals');title('Monitoring set, residuals')
				set(gcf,'Position',setfig);
				men3 = menu('Monitoring results',...
				'1. Print',...
				'2. Next plot');
				if men3 == 1
					print
				end;close
				hist(rmsmon);
				xlabel('RMSEM values'); ylabel('Cumulative frequencies');title('Monitoring set')
				set(gcf,'Position',setfig);
				men3 = menu('Monitoring results',...
				'1. Print',...
				'2. Back to previous menu');
				if men3 == 1
					print
				end;close
			elseif men2 == 3	% Test results
				clc;close
				if testflag
				rtest = corrplot(ytest,ytesref',1);	% Correlation plot, test set
				xlabel('Observed values - Test')
				set(gcf,'Position',setfig);
				men3 = menu('Test results',...
					'1. Print',...
					'2. Next plot');
					if men3 == 1
						print
					end;close
					plot(1:nxte,ytest,'+',1:nxte,ytesref,'co',1:nxte,ytest,'y',1:nxte,ytesref,'c--') 
					renstest = sprintf('RMSEP=%9.6f',rmstestref);
					xlabel(renstest); ylabel('+ Experimental     o Predicted');title('Test set, adjustment')
					set(gcf,'Position',setfig);
					men3 = menu('Test results',...
					'1. Print',...
					'2. Next plot');
					if men3 == 1
						print
					end;close
					plot(1:nxte,residtest,1:nxte,0*residtest,'c--',1:nxte,residtest,'yo')
					xlabel('Sample number'); ylabel('Residuals');title('Test set, residuals')
					set(gcf,'Position',setfig);
					men3 = menu('Test results',...
					'1. Print',...
					'2. Next plot');
					if men3 == 1
						print
					end;close
					hist(rmstest);
					xlabel('RMSEP values'); ylabel('Cumulative frequencies');title('Test set')
					set(gcf,'Position',setfig);
					men3 = menu('Test results',...
					'1. Print',...
					'2. Back to previous menu');
					if men3 == 1
						print
					end;close
				else
					men3 = menu('No test set was provided',...
					'Back to previous menu');
				end
			elseif men2 == 4
				continue1 = 0;
			end
		end
	end
	
	if men1 == 2	% Input variables sensitivity
		clc;close;
		bar(varytrunc)		
		ylabel('Sensitivity');title('Partial modeling')
		%axis([0 mxtr+1 0 1])
		set(gcf,'Position',setfig);
		men2 = menu('Input variables sensitivity',...
		'1. Print',...
		'2. Back to main menu');
		if men2 == 1
			print
		end;close
	end

	if men1 == 3	% Hidden nodes sensitivity
		clc;close
		bar(varytrunch)
		set(gcf,'Position',setfig);
		xlabel('Hidden node');ylabel('Sensitivity');title('Partial modeling')
		men2 = menu('Hidden nodes sensitivity',...
		'1. Print',...
		'2. Back to main menu');
		if men2 == 1
			print
		end;close
	end

	if men1 == 4	% Partial models (input variables)
		for i = 1:mxtr
			partyy = concyhtr(i,:)';
			[y_sort,indd] = sort(partyy);		% Sort y-values of partial models
			partyy = partyy(indd);			% Rank y-values in ascending order
			clc;
			[qc] = polyfit(xtrain(indd,i),partyy,2);	% Quadratic polynomial fit
			qcm(i,:) = qc;
		end
		clc;close
		continue1 = 1;
		while continue1
			men2 = menu('Partial models -Input variables',...
			'1. Training set',...
			'2. Monitoring set',...
			'3. Test set',...
			'4. ANOVA lack-of-fit',...
			'5. Back to main menu');
			if men2 == 1		% Partial model, training set
				clc;close;close
				for i = 1:mxtr
					[b0,bb] = mlr(xtrain(:,i),concyhtr(i,:)',1);	
					subplot(211);plot(xtrain(:,i),concyhtr(i,:),'*');
					renstit = sprintf('Training set - Partial model Y=f(X%2.0f)',i);
					title(renstit)
					ylabel('Y predicted')
					hold;clc;ax1 = axis;					% These lines are necessary to avoid that						
					if abs(qcm(i,1)) > 1e-3					% the printer automatically modifies the axes.
						b0 = qcm(i,3)-((qcm(i,2)-bb)^2)/(4*qcm(i,1));
						ylin = bb*ax1+b0;
						line([ax1(1) ax1(2)],[ylin(1) ylin(2)],'color','r');
						lim = (bb-qcm(i,2))/(2*qcm(i,1));
						line([lim lim],[ax1(3) ax1(4)],'color', 'r');
					end;
					line([ax1(1) ax1(2)],[ax1(3) ax1(3)],'color','w');						
					axis(ax1);
					subplot(212);plot(xtrain(:,i),residtrain,'o')
					hold;clc;line([ax1(1) ax1(2)],[0 0],'color','r');
					ax2 = axis;axis([ax1(1) ax1(2) ax2(3) ax2(4)]);
					if abs(qcm(i,1)) > 1e-3
						line([lim lim],[ax2(3) ax2(4)],'color', 'r');
					end	
					rensx = sprintf('X%2.0f  (Sensitivity: %3.2f)',[i varytrunc(i)]);
					xlabel(rensx)
					ylabel('Residuals')
					set(gcf,'Position',setfig);
					men3 = menu('Partial model - Training',...
					'1. Print',...
					'2. Next plot');
					if men3 == 1
						print
					end;clc;close;close
				end
			elseif men2 == 2	% Partial model, monitoring set
				for i = 1:mxtr
					clc;close
					subplot(211);plot(xmon(:,i),concymhat(i,:),'*');
					renstit = sprintf('Monitoring set - Partial model Y=f(X%2.0f)',i);
					title(renstit)
					ylabel('Y predicted')
					hold;clc;ax1 = axis;					% These 2 lines are necessary to avoid that
					line([ax1(1) ax1(2)],[ax1(3) ax1(3)],'color','w');	% the printer automatically modifies the axes.
					axis(ax1);
					subplot(212);plot(xmon(:,i),residmon,'o')
					hold;clc;line([ax1(1) ax1(2)],[0 0],'color','r');
					ax2 = axis;axis([ax1(1) ax1(2) ax2(3) ax2(4)]);
					rensx = sprintf('X%2.0f  (Sensitivity: %3.2f)',[i varytrunc(i)]);
					xlabel(rensx)
					ylabel('Residuals')
					set(gcf,'Position',setfig);
					men3 = menu('Partial model - Monitoring',...
					'1. Print',...
					'2. Next plot');
					if men3 == 1
						print
					end;clc;close
				end
			elseif men2 == 3	% Partial model, test set
				if testflag
					for i = 1:mxtr
						clc;close
						subplot(211);plot(xtest(:,i),concyteshat(i,:),'*');
						renstit = sprintf('Test set - Partial model Y=f(X%2.0f)',i);
						title(renstit)
						ylabel('Y predicted')
						hold;clc;ax1 = axis;					% These 2 lines are necessary to avoid that
						line([ax1(1) ax1(2)],[ax1(3) ax1(3)],'color','w');	% the printer automatically modifies the axes.
						axis(ax1);
						subplot(212);plot(xtest(:,i),residtest,'o')
						hold;clc;line([ax1(1) ax1(2)],[0 0],'color','r');
						ax2 = axis;axis([ax1(1) ax1(2) ax2(3) ax2(4)]);
						rensx = sprintf('X%2.0f  (Sensitivity: %3.2f)',[i varytrunc(i)]);
						xlabel(rensx)
						ylabel('Residuals')
						set(gcf,'Position',setfig);
						men3 = menu('Partial model - Test',...
						'1. Print',...
						'2. Next plot');
						if men3 == 1
							print
						end;clc;close
					end
				else
					men4 = menu('No test set was provided',...
					'Back to previous menu');
				end
			elseif men2 == 4	% ANOVA for lack-of-fit
				noilev = 0.1*max(max(abs(xtr)));	% The noise level for the input matrix is taken as 10% of the maximum absolute value
				for j = 1:32
					% Noise is added to the input matrix for estimation of 32 replicate training responses
					if j == 1
						nois = zeros(size(xtr));	% No noise added for the first of the ten replicate trials (reference trial)