%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Response surface methodology for studying % instinct-learning balance % % By: Kevin Passino % Version: 5/2/01 % % Note: This code can take a relatively long % time to run (of course depending on % what computer you are using). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all % Initialize memory NL=200; % Set the number of steps in the lifetime of the organism time=1:NL; % Time steps during a lifetime Ng=100; % Set the number of generations timeepoch=1:Ng; % Set environmental characteristics y=zeros(NL,Ng); % Initialize/allocate memory for sensed variable from environment z=zeros(NL,Ng); % Initialize sigmaz2=0.5*ones(NL,Ng); % Initialize the variance in the variable we want to estimate (could % change during a lifetime or over generation) xmin=0; % Set limits on size of x values xmax=4; xbar=2; % For case where have an unknown constant in the environment that want to estimate x=xbar*ones(NL,Ng); % Set variable want to estimate as a constant xhat=zeros(NL,Ng); Ndivxhat0=40; % Set the number of points on a xbar=xhat0 axis xhat0=linspace(xmin,xmax,Ndivxhat0); % Initialize the estimator initial conditions Ndivn=20; % Number of divisions on n space, same as the max value considered for n n=1:Ndivn; % Initialize cost function J=zeros(Ndivn,Ndivxhat0,Ng); % Set weighting parameters for cost function % A choice to illustrate some concepts: w1=0.01; w2=1; w3=0.05; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Start loop for generation of data for response surfaces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ell=1:Ng % Loop for testing multiple times... for ii=1:length(xhat0) for jj=1:length(n) clear Y % Dimension changes on the regressor vector so clear it each time it changes Y=xhat0(ii)*ones(n(jj),1); % Set initial condition % Run estimator for lifetime of animal for k=1:NL % Generate sensed signals, including initial values that are instincts z(k,ell)=sigmaz2(k,ell)*randn; % Generate noise y(k,ell)=x(k,ell)+z(k,ell); % Generate signal that is sensed Y=[y(k,ell); Y(1:(n(jj)-1),1)]; % Shift regression vector, load in new value xhat(k,ell)=mean(Y); end J(jj,ii,ell)=w1*n(jj)... % Cost of storage +w2*(1/NL)*(x(:,ell)-xhat(:,ell))'*(x(:,ell)-xhat(:,ell))... % Quality of estimation +w3*exp(-(x(1,ell)-xhat0(ii))^2/((0.1)^2)); % Models cost of instinct end % End ii loop end % End jj loop end % End ell loop % Compute average of performance measure for Ng runs for ii=1:length(xhat0) for jj=1:length(n) Javg(jj,ii)=mean(J(jj,ii,:)); % Compute average of cost measure across the generations end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [val,rowindex]=min(min(Javg)); % Gives the min value (val) and its row index [val,colindex]=min(min(Javg')); % Gives the min value and its column index % Best point bestpoint=[xhat0(rowindex); n(colindex)] figure(1) clf surf(xhat0,n,Javg) view(-162,48) colormap(jet) colormap(white); hold on bestpointplot=plot3(bestpoint(1,1),bestpoint(2,1),Javg(colindex,rowindex),'r.'); set(bestpointplot,'MarkerSize',20); hold off title('Average cost') xlabel('Initial condition for estimate (x=2)') ylabel('n') rotate3d on figure(2) clf for ii=1:length(xhat0) for jj=1:length(n) plot(timeepoch,squeeze(J(jj,ii,:)),'k-') hold on end end hold off title('Cost') xlabel('Generation') ylabel('J') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%