% J. C. Spall, February 2000. % This code implements the random algorithm B in ISSO. This code assumes % noisy loss measurements. % Evaluates both the averaging and threshold idea. % global sigma p noise n=1000; % no. of loss measurements per run p=10; % problem dimension sigma=1; % noise stand. deviation cases=5; % no. of runs being averaged t=2.0096; % t-value for 2-sided test based on d.o.f.=cases-1 loss_noise='loss4thnoise'; loss='loss4thorder'; s_B=.5; % standard dev. of perturbations in Alg. B avg=20; % no. of loss measurements being averaged threshold=1; % the positive no. expressing the "acceptance threshold" theta_0=ones(10,1); feval(loss,theta_0) theta_lo=-5*ones(p,1); %lower bounds on theta theta_hi=5*ones(p,1); %upper bounds on theta % Algorithm B under averaging meanlossB=0; meanlosssqB=0; meandist=0; truetheta=zeros(10,1); randn('seed',31415927) for i=1:cases theta=theta_0; d(1,i,1)=feval(loss,theta); yold=0; for l=1:avg yold=yold+feval(loss_noise,theta); end yoldavg=yold/avg; for k=1:n/avg -1 thetanew=theta+s_B*randn(p,1); % Two lines below invoke constraints thetanew=min(thetanew,theta_hi); thetanew=max(thetanew,theta_lo); ynew=0; for l=1:avg ynew=ynew+feval(loss_noise,thetanew); end ynewavg=ynew/avg; if ynewavg < yoldavg theta=thetanew; yoldavg=ynewavg; else end d(k+1,i,1)=feval(loss,theta); %records true loss value for ith case, kth iteration end meanlossB=(i-1)*meanlossB/i+feval(loss,theta)/i; meanlosssqB=(i-1)*meanlosssqB/i+feval(loss,theta)^2/i; meandist=(i-1)*meandist/i+norm(theta-truetheta,2)/i; theta end disp('mean and confidence intervals under averaging') meanlossB if cases > 1 sd_B=(cases^(-.5))*((cases/(cases-1))^.5)*(meanlosssqB-meanlossB^2)^.5;sd_B % Confidence intervals [meanlossB-t*sd_B,meanlossB+t*sd_B] else end meandist % Algorithm B with modified acceptance criterion("threshold") meanlossB=0; meanlosssqB=0; meandist=0; randn('seed',11113); for i=1:1 theta=theta_0; yold=feval(loss_noise,theta); for k=1:n-1 thetanew=theta+s_B*randn(p,1); thetanew=min(thetanew,theta_hi); thetanew=max(thetanew,theta_lo); ynew=feval(loss_noise,thetanew); if ynew < yold - threshold theta=thetanew; yold=ynew; else end end meanlossB=(i-1)*meanlossB/i+feval(loss,theta)/i; meanlosssqB=(i-1)*meanlosssqB/i+feval(loss,theta)^2/i; % meandist=(i-1)*meandist/i+norm(theta-[1,0.66667]',2)/i; end disp('mean and confidence intervals under altered threshold') meanlossB if cases > 1 sd_B=(cases^(-.5))*((cases/(cases-1))^.5)*(meanlosssqB-meanlossB^2)^.5;sd_B % Confidence intervals [meanlossB-t*sd_B,meanlossB+t*sd_B] else end meandist