%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Response surface method for finding optimal approximator size p
% and data set G size M for a function approximation problem.
% Program also allow you to keep either p or M constant and
% vary the other one (to generate a 2D response surface).
%
% By: Kevin Passino
% Version: 3/22/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
% Consider the follow range of sizes of the approximator
p=4:2:40; % Notice that we consider only even p for convenience since
% it is easiest to add two parameters at a time for the
% type of approximator that we are using (TS fuzzy system)
% Also, must start at p(1)>=4 due to how the spreads are defined
% in terms of the centers.
% Consider the following range of sizes of the data set
M=100:100;
% Generate a test set
% Define a parameter to control the size of the test set
% relative to sizes of training sets
gamma=2; % Choose the size of the test set to be twice the size
% of the largest training set
[xt,Gamma]=unknownfunction(gamma*max(M));
% Set the number of trials for each point on the RSM
% (to use an average MSE for each point on the RSM)
Ntrials=1;
%Ntrials=100;
% Initialize the variables that hold the training data:
x=0*ones(length(M),length(p),length(M),Ntrials);
G=0*ones(length(M),length(p),length(M),Ntrials);
% Set a variable for the number of times to run the program (for checking
% that the response surface does not change too much for additional trials)
Ntests=1;
Ntests=5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop for generating data for RSM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ll=1:Ntests
for ii=1:length(p)
for jj=1:length(M)
for kk=1:Ntrials
[x(1:M(jj),ii,jj,kk),G(1:M(jj),ii,jj,kk)]=unknownfunction(M(jj)); % Generate training data
% First, form the vector Y
Y=G(1:M(jj),ii,jj,kk);
% Find Phi, which involves processing x through phi (for the training data)
% Find the centers of the membership functions
c(1,:)=-6:(12/(0.5*p(ii)-1)):6; % Produces one center for every two parameters
sigma(1,:)=((c(1,2)-c(1,1))/(2))*ones(1,length(c(1,:))); % Chooses width to be half the distance
% between adjacent centers
% Initialize Phi
for j=1:length(c(1,:))
mu(j,1)=exp(-0.5*((x(1,ii,jj,kk)-c(1,j))/sigma(1,j))^2);
end
denominator(1)=sum(mu(:,1));
for j=1:length(c(1,:))
xi(j,1)=mu(j,1)/denominator(1);
end
Phi=[xi(:,1)', x(1,ii,jj,kk)*xi(:,1)'];
% Form the rest of Phi
for i=2:M(jj)
for j=1:length(c(1,:))
mu(j,i)=exp(-0.5*((x(i,ii,jj,kk)-c(1,j))/sigma(1,j))^2);
end
denominator(i)=sum(mu(:,i));
for j=1:length(c(1,:))
xi(j,i)=mu(j,i)/denominator(i);
end
Phi=[Phi; xi(:,i)', x(i,ii,jj,kk)*xi(:,i)'];
end
% Find the least squares estimate
theta=Phi\Y;
% Compute the approximator values for the test set (xt,Gamma)
for i=1:length(Gamma)
for j=1:length(c(1,:))
mut(j,i)=exp(-0.5*((xt(i)-c(1,j))/sigma(1,j))^2);
end
denominatort(i)=sum(mut(:,i));
for j=1:length(c(1,:))
xit(j,i)=mut(j,i)/denominatort(i);
end
phit=[xit(:,i)', xt(i)*xit(:,i)']';
Fts(i,1)=theta'*phit;
end
% Plot the data (optional - just to gain insight)
% figure(1)
% clf
% plot(xt,Fts,'o')
% ylabel('Approximator output F_t_s')
% xlabel('Test data')
% Tp=num2str(p(ii));
% TM=num2str(M(jj));
% T=strcat('Approximator mapping, p=',Tp,'Training set size, M=',TM);
% title(T)
% pause
% Compute the approximation error at the test data
J(jj,ii,kk)=(1/length(Gamma))*(Gamma-Fts(:,1))'*(Gamma-Fts(:,1)); % Compute the mean-squared error
end % End Ntrials loop
% Need to clear variables as dimensions change in two outer-most loops
clear Y c sigma mu denominator Phi xi theta mut denominatort xit phit
% Compute the average mean squared error
Jsurf(jj,ii,ll)=(1/Ntrials)*sum(J(jj,ii,:));
end % End M loop
end % End p loop
end % Ends the Ntests loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the response surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
clf
% First consider 2D surfaces
if length(M)==1, % Plot across p
for ll=1:Ntests
figure(2)
plot(p,squeeze(Jsurf(1,:,ll)),'k-')
xlabel('Approximator size, p');
ylabel('Average MSE, J');
TM=num2str(M(1));
T=strcat('Response surface, Training set size, M=',TM);
title(T)
% Next line for zooming in on response surface region
%axis([10 max(p) min(Jsurf(1,:,ll)) 0.03])
hold on
end
end
%clear J
% Below, some additional plotting possibilities
%if length(p)==1, % Plot across M
%figure(3)
%clf
%plot(M,Jsurf(:,1),'k-')
%xlabel('Training data set size, M');
%ylabel('Average MSE, J');
%Tp=num2str(p(1));
%T=strcat('Response surface, Approximator size, p=',Tp);
%title(T)
%end
% Next, 3D
%if length(M)>1 & length(p)>1, % If have data for a 3D surface
%figure(4)
%clf
%surf(p,M,Jsurf)
%colormap(white);
%xlabel('Approximator size, p');
%ylabel('Training data set size, M');
%zlabel('Average MSE, J');
%title('Response surface for average MSE');
%rotate3d
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%