examp.m

这是在网上下的一个东东

%
% Reinitialize everything.
%
clear
rand('seed',0);
randn('seed',0);
%
% Load in the Shaw problem.
%
load shawexamp.mat
%
% First, plot the l_curve and find its corner.
%
[U,S,V]=svd(G);
figure(1);
bookfonts;
[reg_corner,rho,eta,reg_param]=l_curve(U,diag(S),dspiken,'tikh');
%print -deps c5flcurve.eps
%
% get the solution.
%
figure(2);
bookfonts;
[mtik,rho,eta]=tikhonov(U,diag(S),V,dspiken,reg_corner);
plotconst(mtik,-pi/2,pi/2);
xlabel('\theta');
ylabel('Intensity');
%print -deps c5fmtik.eps

%
% Use the discrepancy principle to get a second solution.
%
[mdisc,alpha]=discrep(U,diag(S),V,dspiken,1.0e-6*sqrt(20));
alpha
figure(3);
bookfonts;
plotconst(mdisc,-pi/2,pi/2);
xlabel('\theta');
ylabel('Intensity');
%print -deps c5fmdisc.eps

%
% Get the Picard plot.
%
figure(4);
bookfonts;
eta=picard(U,diag(S),dspiken);
%
% We manually moved the legend in this plot to make it look better, so
% this is commented out so that we don't have to fix the legend each time.
%
%%print -deps c5fpicard

%
% Now, get at the resolution.
%
figure(5);
bookfonts;
rdspike=tikhonov(U,diag(S),V,dspike,alpha);
plotconst(rdspike,-pi/2,pi/2);
xlabel('\theta');
ylabel('Intensity');
%print -deps c5frdspike.eps
%
% Compute the resolution matrix too.
%
Ginv=inv(G'*G+alpha^2*eye(20))*G';
R=Ginv*G;
%
% Now, produce a plot of the discrepancy principle solution with "error
% bars" versus reality.  
%
figure(6);
bookfonts;
covmdisc=Ginv*1.0e-6^2*eye(20)*Ginv';
H1=plotconstc(spike,-pi/2,pi/2,'k-');
hold on
H2=plotconstc(mdisc,-pi/2,pi/2,'k--');
H3=plotconstc(mdisc+1.96*sqrt(diag(covmdisc)),-pi/2,pi/2,'k:');
H4=plotconstc(mdisc-1.96*sqrt(diag(covmdisc)),-pi/2,pi/2,'k:');
xlabel('\theta');
ylabel('Intensity');
legend([H1; H2; H3],'m_{true}','m_{disc}','95% CI');
%print -deps c5fbias.eps
%
% Tikhonov's theorem calculations
%
w=(G')\spike;
norm(w)
%
% This is hopelessly large, so we'll give up.
%
%
% Next, a smoother model.
%
w=[ones(3,1); zeros(17,1)];
smoothmod=G'*w;
smoothmodd=G*smoothmod;
smoothmoddn=smoothmodd+1.0e-6*randn(20,1);
p=1;
gamma=1/2;
delta=sqrt(20)*1.0e-6;
Delta=delta/norm(w);
alphahat=(Delta/2*gamma*1)^(1/(1+p));
bound=gamma*(p+1)*alphahat^p
malphahat=tikhonov(U,diag(S),V,smoothmoddn,alphahat)
norm(malphahat-smoothmod)
figure(7);
bookfonts;
plotconst(smoothmod,-pi/2,pi/2);
xlabel('\theta');
ylabel('Intensity');
%print -deps c5fsmoothmod.eps

figure(8);
bookfonts;
plotconst(malphahat,-pi/2,pi/2);
xlabel('\theta');
ylabel('Intensity');
%print -deps c5fsmoothmodre.eps