home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 3.ddi / MMLE3.DI$ / MMLEMARQ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.0 KB  |  69 lines
  1. function [pnew,llf,marq]=mmlemarq(p,pid,p2snam,gg,marq,llf)
  2. % USED BY MMLE.M TO GET NEAR MINIMUM (=>QUADRATIC AREA ) OF COST FUNCTION
  3.  
  4.  % -  USES FUNCTIONS mmlelike, AND USER CREATED FUNCTION p2ss.
  5.  % -  COMPUTES PARAMETER VECTOR "pnew" TO IMPROVE llf FOR TOUGH PROBLEMS.
  6.  % -  USES MARQUARD PARAMETER marq TO FIND BEST GAIN/DIRECTION COMPROMISE OF
  7.  %       GRADIENT AND NEWTON UPDATE TO OBTAIN THE BEST llf IMPROVEMENT
  8.  %            FROM EACH EXPENSIVE GRADIENT/HESSIAN COMPUTATION.
  9.  % -  marq >>1 GIVES SLOW GRADIENT ALGO.
  10.  % -  marq <<1 APPROACHES NEWTON UPDATE.
  11.  % -  marq =0 DEFAULTS TO marq =1.
  12.  % -  marq IS AUTOMATICALLY REDUCED AS FAST AS PROBLEM PERMITS
  13.  %         TO SPEED USE OF CONSTRAINED NEWTON ALGO BY mmle.m
  14.  % -  pnew IS NEW PARAM, llf IS BEST llf, marq IS LAST marq.
  15.  % -  mmlestep, mmlelike DESCRIBE THE OTHER VARIABLES.
  16.  %---------------------------------------------------------------------------
  17.  
  18.  
  19. global hes
  20. global wersum
  21. global grad
  22. global rr
  23.  
  24.  
  25. eyenpid = eye(length(pid));
  26. %----------------------------------------SCALE UNITS FOR UNIT DIAGONAL
  27. scale = diag(exp(-.5*log(diag(hes))));
  28. shes = scale*hes*scale;sgrad = scale*grad;
  29. %
  30. pnew_t=p;pnew=p;if marq == 0,marq = 1;end %   DEFAULT OPTION
  31. trend     = 0; %  MEANINGS =>>>   0:FIRST LOOP, -1:WAS WORSE, +1:WAS BETTER
  32.  
  33. %------------------------------------------------DISPLAY MAX GRADIENT ELEMENT
  34. absgr=abs(grad);mgr=max(absgr);igr=find(absgr==mgr);
  35. ans4=sprintf('%f',mgr);ans5=sprintf('%f',pid(igr));
  36. ans1=[' --- MMLEMARQ ---   max(abs(grad)) = ',ans4(1:8), ...
  37.           ' -- p( ',ans5(1:2),')'];disp(ans1)
  38. %--------------------------------------------------------------START LOOP
  39. for nstep = 1:5,
  40.   ans =-scale*real((shes+marq*eyenpid)\sgrad ); %  PARAMETER STEP
  41.   pnew_t(pid)=p(pid) + ans'; %                     TRIAL PARAMETER VECTOR
  42.   [llf_t,wersum_t,rr_t] = mmlelike( pnew_t, p2snam, gg );%    GET ITS LLF
  43.     if llf_t==NaN;llf_t=8888888;end;
  44.   ans1=sprintf('%f',marq);ans2=sprintf('%f',llf);
  45.   ans3=sprintf('%f',llf_t);ans4=sprintf('%f',wersum_t);
  46.   ans4=[ ' marq = ', ans1(1:7), '   best llf = ',ans2(1:7),...
  47.    '   trial llf = ', ans3(1:7),'   wersum = ',ans4(1:7)];disp(ans4)
  48.  
  49. %-----------------------------------------------------MARQUARDT DECISIONS
  50.   if llf_t > llf, %           LLF IS WORSE
  51.       if trend == 1,return,end %   WAS BETTER, NOW WORSE => RETURN
  52.          if trend == -1,
  53.          if abs(llf_t-llf)<.5,return% BACK WITHIN .5 OF BEST LLF => RETN
  54.          end,end
  55.       marq=marq*10; %         INCREASE MARQ
  56.       trend=-1; %             REMEMBER THIS STEP WAS WORSE
  57.    else
  58.       llf=llf_t;rr=rr_t;%     UPDATE BEST LLF,RR
  59.       wersum=wersum_t; %      WERSUM
  60.       pnew=pnew_t; %          SAVE BEST PARAMETERS FOUND
  61.          if trend == -1, return,end %   MARQ NOW HIGH ENOUGH TO RETURN
  62.       mold=marq;marq = marq/sqrt(10);
  63.       trend=1; %              REMEMBER THIS WAS BETTER
  64.    end
  65. end %                         END OF 'FOR' LOOP
  66. return
  67. %------------------------------------------------------------- end mmlemarq
  68.  
  69.