examp.m

这是在网上下的一个东东

%
% First, setup some global variables for later use.
%
global A;
global b;
global w;
global alpha;
global K;
%
% Reset the random number generator.
%
rand('state',0);
%
% Now, create the problem.  We use the Shaw problem with n=20.
%
A=shaw(20);
%
% A random background level.
%
xtrue=0.2*rand(20,1)+0.4*ones(20,1);
%
% A spike of 10 at index 10.  A second spike of 4 at index 15.
%
xtrue(7)=1000;
xtrue(15)=50;
%
% Calculate the right hand side.
%
btrue=A*xtrue;
b=btrue+1.0*randn(20,1);
%
% We'll use weights of 1.0;
%
w=ones(20,1);
%
% Pick a reasonable bound on ||Ax-b||
%
delta=1.0*sqrt(20);
%
% A useful starting point.
%
x0=norm(b)/norm(A,1)*ones(20,1);

%
% Solve the sample problem using our MaxEnt implementation.
%
%
% A penalty parameter of 1.0e12 seems to work well.
%
K=1.0e8;
%
% First, make an L-curve.
%
alphas=0.15:0.05:0.5;
misfits=zeros(length(alphas),1);
ents=zeros(length(alphas),1);
for i=1:length(alphas)
  alpha=alphas(i);
%
% The toolbox maxent() command isn't very reliable- it fails to get a
% good solution for small values of alpha.  
%
% xmaxent=maxent(A,b,alpha);
%
% This line uses bfgs to do the minimization.  Much more robust, but a
% bit slower.
%
  xmaxent=bfgs('maxentobj','maxentgrd',x0,1.0e-10,2500);
  
  misfits(i)=norm(A*xmaxent-b); 
  ents(i)=(maxentobj(xmaxent)-norm(A*xmaxent-b))/alpha^2;
end
%
% Plot the L-curve.
%
figure(1);
clf;
bookfonts;
plot(misfits,ents,'ko-');
%print -deps maxentl.eps
%
% It appears that the regularization parameter alpha=0.2 will yield 
% a misfit of just about sqrt(20), which is what we want given the
% noise in the data.
%
alpha=0.2;
xmaxent=bfgs('maxentobj','maxentgrd',x0,1.0e-10,2000);
%
% Now, plot the MaxEnt solution.
%
figure(2);
clf;
bookfonts
plotconst(xmaxent,-pi/2,pi/2);
axis([-2 2 -100 1200]);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps maxentsol.eps
%
% Next, we'll try lsqnonneg.
%
xnnls=lsqnonneg(A,b);
figure(3);
clf;
bookfonts
plotconst(xnnls,-pi/2,pi/2);
axis([-2 2 -100 1200]);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps nnlssol.eps
%
% Next, tikhcstr.
%
%
% First, make an L-curve.
%
alphas=0.005:0.005:0.05;
misfits=zeros(length(alphas),1);
normx=zeros(length(alphas),1);
G=eye(20);
d=zeros(20,1);
for i=1:length(alphas)
  alpha=alphas(i);
  xtik=tikhcstr(A,b,G,d,eye(20),alpha);
  misfits(i)=norm(A*xtik-b); 
  normx(i)=norm(xtik);
end
figure(4);
clf;
bookfonts
plot(misfits,normx,'ko-');

%print -deps tikhl.eps
%
% alpha=0.007 satisfies the discrepancy principle.
%
xtik=tikhcstr(A,b,G,d,eye(20),0.007);
figure(5);
clf;
bookfonts
plotconst(xtik,-pi/2,pi/2);
axis([-2 2 -100 1200]);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps tiksol.eps
%
% Plot the true solution.
%
figure(6);
clf;
bookfonts
plotconst(xtrue,-pi/2,pi/2);
axis([-2 2 -100 1200]);
xlabel('\theta (radians)')
ylabel('Intensity')
%print -deps truesol.eps