MakeResponseSurface.m
From charlesreid1
The following is the Matlab code for the function that creates the response surfaces and plots/prints useful information about them.
M file is available here: http://files.charlesmartinreid.com/ExperimentalDesign/MakeResponseSurface.m
function response_surface = MakeResponseSurface( X, Y, logscale, names, deg, file_prefix, plotXindex, plotYindex )
% ---------------- Initialization stuff ----------------------
% Ensure X and Y are similar lengths
if( length(X) ~= length(Y) )
error('Mismatch in input lengths for X and Y.\n\n');
end
% Ensure X and Y are oriented correctly
% size(X,1) = nubmer of input variables
% size(X,2) = number of samples
% size(Y,1) = number of output variables
% size(Y,2) = number of samples
if( size(X,1) > size(X,2) )
X=X';
end
if( size(Y,1) > size(Y,2) )
Y=Y';
end
N_xvars = size(X,1);
N_yvars = size(Y,1);
% Y should not be multivariate
% (you can't plot multivariate surf. so it's not worth dealing with)
if( size(Y,1) > 1 )
fprintf('\n\n*******WARNING: Y is multivariate (%0d dimensions). \n\n',size(Y,1));
error('Y cannot be multivariate.\n\n');
end
% If no filename specified, it's OK - just don't save
if nargin < 6
fprintf('Not saving response surface to file. \n\n');
do_save = false;
else
fprintf('Saving response surface to file %s.mat \n',file_prefix);
fprintf('Saving plot to file %s.png \n\n',file_prefix);
do_save = true;
end
% If no plot indices specified, use 1 and 2
if nargin < 7
plotXindex = 1;
plotYindex = 2;
end
% If no degree specified, use quadratic (or other)
if nargin < 5
deg_is_matrix = false;
deg = 2;
fprintf('You did not specify a polynomial degree. Using max. degree %0d. \n',deg);
% If deg is a matrix, using a custom polynomial
elseif size(deg,2) > 1
deg_is_matrix = true;
fprintf('Creating response surface with specified polynomial form. \n');
else
deg_is_matrix = false;
fprintf('Creating response surface of degree %0d. \n',deg);
end
fprintf('\n');
% ------------ Response surface supplementary info ---------
alph = zeros(N_xvars,1); % mins
beta = zeros(N_xvars,1); % maxes
med = zeros(N_xvars,1); % medians
% alph and beta: min/max for each var
% (determines domain for plotting response surface)
for i=1:N_xvars
alph(i) = min(X(i,:)');
beta(i) = max(X(i,:)');
end
% Create an empty response surface
Npoints = 50; % number of points in each direction
% (for plotting response surface)
Xrs = zeros(N_xvars,Npoints);
% Populate x-values for RS
for i=1:N_xvars
%%% if( logscale(i) == true )
%%% Xrs(i,:) = linspace(alph(i),beta(i),Npoints);
%%% else
%%% Xrs(i,:) = logspace(log10(alph(i)),log10(beta(i)),Npoints);
%%% end
Xrs(i,:) = linspace(alph(i),beta(i),Npoints);
med(i) = median(Xrs(i,:));
end
% ---------------- Regression --------------------
fprintf('Regressing data... ');
% Determine poly model form
PolyModelPowers = 0;
if( deg_is_matrix )
fprintf('Setting polynomial model degrees to specified set.\n\n');
PolyModelPowers = deg;
else
fprintf('Setting polynomial model degree to %0.0f.\n\n',deg);
nterms = length(allVL1(N_xvars, deg, '<='));
if( nterms > length(X) )
fprintf('\n');
fprintf('NOTICE: The polynomial of degree %0d has %0d terms; you have %0d function samples.\n', deg, nterms, length(X) );
fprintf('The polynomial is underdetermined, so try reducing the specified polynomial degree.\n');
fprintf('Attempting to reduce the number of terms...\n');
temp_eye = eye(N_xvars)*deg;
PolyModelPowers = [allVL1(N_xvars, deg-1, '<='); temp_eye];
fprintf('Reduced polynomial to %0d terms by eliminating (less important) higher-order interaction terms. \n', length(PolyModelPowers) );
else
PolyModelPowers = allVL1(N_xvars, deg, '<=');
end
end
%X = [X(plotXindex,:); X(plotYindex,:)];
response_surface = regstats(Y',X',PolyModelPowers,{'beta','covb','yhat','mse','r','rsquare','adjrsquare','tstat','fstat'});
b = response_surface.beta;
mse = response_surface.mse;
rsq = response_surface.rsquare;
arsq = response_surface.adjrsquare;
resid = response_surface.r;
fprintf('Done.\n\n');
% ------------ Populate response surface matrix ---------
% Populate y-values for RS automatically:
count = 1;
Yrs = zeros(Npoints,Npoints);
Ysp = zeros(1,Npoints*Npoints);
Xsp = zeros(2,Npoints*Npoints);
for ii=1:Npoints
for jj=1:Npoints
Ysp(count) = 0;
Yrs(ii,jj) = 0;
for kk=1:size(PolyModelPowers,1)
polyprod = 1;
for nn=1:size(PolyModelPowers,2)
if( nn == plotXindex )
polyprod = polyprod*Xrs(nn,ii)^(PolyModelPowers(kk,nn));
elseif ( nn == plotYindex )
polyprod = polyprod*Xrs(nn,jj)^(PolyModelPowers(kk,nn));
else
if( logscale(nn) == true )
polyprod = polyprod*( 10^( mean([log10(alph(nn)) log10(beta(nn)) ]) ) )^(PolyModelPowers(kk,nn));
else
polyprod = polyprod*(mean([alph(nn) beta(nn)]))^(PolyModelPowers(kk,nn));
end
end
end
Ysp(count) = Ysp(count) + b(kk)*polyprod;
Yrs(ii,jj) = Yrs(ii,jj) + b(kk)*polyprod;
end
Xrs2(1,ii) = Xrs(plotXindex,ii);
Xrs2(2,jj) = Xrs(plotYindex,jj);
Xsp(1,count) = Xrs(plotXindex,ii);
Xsp(2,count) = Xrs(plotYindex,jj);
count = count + 1;
end
end
% That whole thing can probably be replaced with x2fx(X,deg)
% where X = [x1_1 x2_1 x3_1 x4_1 etc...; x1_2 x2_2 ....]
% ------------ Plot response surface ---------
% plot that shiznizzle
figure()
surf(Xrs2(1,:),Xrs2(2,:),Yrs','EdgeColor','None')
% ---------------------
% Set plot attributes:
view(3);
grid on;
% If variable is log scale, change axis scale
if( logscale(plotXindex) == true )
set(gca,'XScale','log');
end
if( logscale(plotYindex) == true )
set(gca,'YScale','log');
end
% Set title
if( deg_is_matrix )
title('Polynomial Response Surface');
else
title(['Degree ',num2str(deg),' Polynomial Response Surface']);
end
% Axis labels
xlabel(names{plotXindex});
ylabel(names{plotYindex});
zlabel('Y_p');
% ------------ Plot original responses (scatterplot) ---------
% Option 1: scatter plot in new figure
%figure()
% Option 2 (preferred): scatter plot on top of surf plot
hold on;
% plot that keeznacha
scatter3(X(plotXindex,:),X(plotYindex,:),Y,'ko');
%stem3(X(plotXindex,:),X(plotYindex,:),Y,'Marker','o','MarkerSize',4,'Color','black');
if( do_save )
filename = [file_prefix,'.png'];
saveas(gcf,filename);
end
%%% % -------------- Plot residuals --------------------------
%%%
%%% figure()
%%% stem3(X(plotXindex,:),X(plotYindex,:),abs(resid),'filled');
%%%
%%% % If variable is log scale, change axis scale
%%% if( logscale(plotXindex) == true )
%%% set(gca,'XScale','log');
%%% end
%%% if( logscale(plotYindex) == true )
%%% set(gca,'YScale','log');
%%% end
%%%
%%% % the title
%%% title('Response Surface Regression Residuals');
%%%
%%% if( do_save )
%%% filename = [file_prefix,'_residuals.png'];
%%% saveas(gcf,filename);
%%% end
% ------------ Print goodness-of-fit info ---------
% If saving everything to files, write the statistics to a file
if( do_save )
statistics_filename = [file_prefix,'_statistics.txt'];
fid=fopen(statistics_filename,'w');
fprintf(fid,'---------------------------------------------------\n');
fprintf(fid,'Response surface summary of information:\n');
fprintf(fid,'Number of variables in response surface is %0d. \n',N_xvars);
fprintf(fid,'Number of terms in polynomial is %0d. \n',size(PolyModelPowers,1));
fprintf(fid,'Degree of response surface is ');
if( deg_is_matrix )
fprintf(fid,'varied, deg is a matrix. Max degree = %0d.\n',max(max(deg)));
else
fprintf(fid,'%0d.\n',deg);
end
fprintf(fid,'\n');
fprintf(fid,'MSE =\t\t\t %0.8f \n',mse);
fprintf(fid,'MSE DoF = \t\t %0.0f \n', sum(resid.^2)/mse);
fprintf(fid,'\n');
fprintf(fid,'L-inf norm resid = \t %0.8f \n',norm(resid,Inf));
fprintf(fid,'\n');
fprintf(fid,'R^2 =\t\t\t %0.8f \n',rsq);
fprintf(fid,'adjusted R^2 =\t\t %0.8f \n',arsq);
fprintf(fid,'---------------------------------------------------\n');
fprintf(fid,'\n');
end
% Now print the statistics (whether or not we're saving to files)
fprintf('---------------------------------------------------\n');
fprintf('Response surface summary of information:\n');
fprintf('Number of variables in response surface is %0d. \n',N_xvars);
fprintf('Number of terms in polynomial is %0d. \n',size(PolyModelPowers,1));
fprintf('Degree of response surface is ');
if( deg_is_matrix )
fprintf('varied, deg is a matrix. Max degree = %0d.\n',max(max(deg)));
else
fprintf('%0d.\n',deg);
end
fprintf('\n');
fprintf('MSE =\t\t\t %0.8f \n',mse);
fprintf('MSE DoF = \t\t %0.0f \n', sum(resid.^2)/mse);
fprintf('\n');
fprintf('L-inf norm resid = \t %0.8f \n',norm(resid,Inf));
fprintf('\n');
fprintf('R^2 =\t\t\t %0.8f \n',rsq);
fprintf('adjusted R^2 =\t\t %0.8f \n',arsq);
fprintf('---------------------------------------------------\n');
fprintf('\n');
if( do_save )
model = PolyModelPowers;
filename = [file_prefix,'.mat'];
save(filename,'response_surface','model');
end