function [p,FN,par,d]= noisyminimize(par)
% minimized a function "fun" using noisy estimates of the function and its
% gradient. start with the parameter values par(1,:)...par(K,:)
% I had problems with the parameter delta so I fixed it at .513
% apparent minimum is parbest=[1.89   -8.66   -2.29]  corresponding
% to mu=0.06, alpha=13.752, beta=-7.72, delta=0.098
delta=.5;
m=size(par,1); FN=[];GR=[]; delt=.01/m; sx=[]; d=[]; q=.1;  % q governs how fast delt ->0
for i=1:m
    [f,g]=fun2(par(i,:),i);
    FN=[FN;f]; GR=[GR;g];
end
[x,i]=sort(FN);   FN=FN(i,:);   GR=GR(i,:);  % sort the function from worst to best
x=repmat(par(m,:),m,1)-par;
x=(sum(x.*x,2)<delt);  % 0 or one as the parameter is within r=sqrt(delt);
gav=mean(GR(x,:),1); fav=mean(FN(x,:));
FN(m)=fav; GR(m,:)=gav;
while m<=5000;     % <-------change the number of iterations here.............
parnew=par(m,:)-max(-.02,min(.02,delt*GR(m,:)));  % the maximum jump size is .01
m=m+1; par=[par;parnew];
[f,g]=fun2(parnew,m); 
display([' new function value ' num2str(f) ' new parameters ' mat2str(parnew)]);
Rm=1-(1+q)*((GR(m-1,:)*g')>0);   %base the change in step size on the observed new gradient
delt=max(.000001,min(.5,delt*(1-Rm*delt))); % do not let delta go negative or >.5
d=[d delt];
FN=[FN;f]; GR=[GR;g];
% choose all parameter values within radius r=delt^p  of par(m,:) and let
% fav=average function value, gav=average gradient
%Rm=1-cosangle(gav,GR(m-1,:)); % continuous version
x=repmat(par(m,:),m,1)-par;
x=(sum(x.*x,2)<delt^1.8);  % 0 or one as the parameter is within r=sqrt(delt);
sx=[sx sum(x)];
gav=mean(GR(x,:),1); fav=mean(FN(x,:));
FN(m)=fav; GR(m,:)=gav;   
end
plot(FN); title('n versus value of function')
figure;    plot(d.*(1:length(d)));    title('n \delta')
figure;      comet(par(:,1),par(:,2))
figure;    plot(sx);   title('number of points combined')
[m,i]=min(FN);
p=par(i,:);