occam.m

这是在网上下的一个东东

%
% occam(fun,jac,L,d,m0,delta)
%
function m=occam(fun,jac,L,d,m0,delta)
m=m0;
oldm=zeros(size(m));
iter=0;
while ((norm(oldm-m)/norm(m) >5.0e-3) | (mchi2 > delta^2*1.01))
  oldm=m;
  iter=iter+1;
  if (iter > 30)
    return;
  end
  G=feval(fun,m);
  norm(G-d);
  J=feval(jac,m);
  dhat=d-G+J*m;
  
%
%  This is a simple brute force way to do the line search.  Much more
%  sophisticated methods are available.  Note: we've restricted the line
%  search to the range from 1.0e-20 to 1.  This seems to work well in
%  practice, but might need to be adjusted for a particular problem.
% 
  alphas=logspace(-20,0,100);
  for i=1:100,
    M=J'*J+alphas(i)^2*L'*L;
    if (cond(M) < 1.0e15)
      m=inv(J'*J+alphas(i)^2*L'*L)*J'*dhat;
      chis(i)=norm(feval(fun,m)-d,2)^2;
    else
      chis(i)=+Inf;
    end
  end
  [Y,I]=min(chis);
  if (Y > delta^2)
    disp('Improving Chi^2');
    alpha=alphas(I(1));
  else
    disp('Smoothing m');
    I=find(chis<=delta^2);
    alpha=alphas(max(I));
  end
  m=inv(J'*J+alpha^2*L'*L)*J'*dhat;
  mchi2=norm(feval(fun,m)-d,2)^2;
%  plotconst(m,0,1000);
%  plot(m);
%  shg;
end