examp.m

这是在网上下的一个东东

%
% Reset to a known state.
%
clear
rand('state',0);
randn('state',0);
%
% Global variables for observed Theta and h values, etc.
%
global H;
global TM;
global SIGMA;
global Q;
global D;
%
% The unknown/estimated parameters are S=p(1) and T=p(2).
%
%
% The data.  Head is measured to the nearest centimeter.  
%
H=[0.72 0.49 0.30 0.20 0.16 0.12];
TM=[5.0 10.0 20.0 30.0 40.0 50.0];
%
% Fixed parameter values.
%
D=60;
Q=50;
%
% We'll use sigma=1cm.  
%
for i=1:length(TM)
  SIGMA(i)=0.01;
end
%
% Solve the least squares problem with LM.
%
p0=[0.001; 1.0];
[pest,itercnt]=lm('fun','jac',p0,1.0e-12,100);
%
% Check the chi^2 value.
%
chi2=norm(fun(pest))^2
pvalue=1-chi2cdf(chi2,length(H)-2)
%
% Compute the covariance matrix, correlation matrix, and confidence intervals.
%
Jest=jac(pest);
C=inv(Jest'*Jest);
m=max(size(C));
for i=1:m
  for j=1:m
    Corm(i,j)=C(i,j)/(sqrt(C(i,i))*sqrt(C(j,j)));
  end
end
fprintf('S=%.5f +- %.5f\n',[pest(1) 1.96*sqrt(C(1,1))]);
fprintf('T=%.5f +- %.5f\n',[pest(2) 1.96*sqrt(C(2,2))]);
fprintf('Covariance matrix is \n');
C
fprintf('Correlation matrix is \n');
Corm
%
% Now, plot the data and fitted curve.
%
figure(1);
clf;
bookfonts;
tfit=0.001:1.0:51;
for i=1:length(tfit),
  hfit(i)=Q*exp(-D^2*pest(1)/(4*pest(2)*tfit(i)))/(4*pi*pest(2)*tfit(i));
end
plot(tfit,hfit,'k-');
hold on
errorbar(TM,H,SIGMA,'ko');
xlabel('Time (hours)');
ylabel('Head (m)');


%print -deps slugplot.eps
%
% Produce a contour plot of the chi^2 function.
%
figure(2);
clf;
bookfonts;
[X,Y]=meshgrid(0.0001:0.0001:0.01,0.01:0.01:2.0);
Z=zeros(size(X));
[m,n]=size(X);
for i=1:m
  for j=1:n
    Z(i,j)=norm(fun([X(i,j); Y(i,j)]))^2;
  end
end  
[CS,HHH]=contour(X,Y,Z,[10 100 500  1000 2000]);
for i=1:length(HHH),
  set(HHH(i),'EdgeColor','k');
end;
%
% These labels don't look to good, so we'll label by hand.
%
HH=xlabel('S');
set(HH,'FontSize',18);
HH=ylabel('T (m^2/hr)');
set(HH,'FontSize',18);


%print -deps slugcontour.eps
%
% Plot a second contour plot much closer to the fitted values.
%
figure(3);
clf;
bookfonts;
[X,Y]=meshgrid(0.0017:0.000025:0.0025,0.50:0.001:0.70);
Z=zeros(size(X));
[m,n]=size(X);
for i=1:m
  for j=1:n
    Z(i,j)=norm(fun([X(i,j); Y(i,j)]))^2;
  end
end  
[CS,HHH]=contour(X,Y,Z,[chi2+chi2inv(0.95,2) chi2+chi2inv(0.95,2)]);
set(HHH(1),'EdgeColor','k');
HH=xlabel('S');
set(HH,'FontSize',18);
HH=ylabel('T (m^2/hr)');
set(HH,'FontSize',18);
%
% Put in the individual 95% CI's.
%
hold on
plot([pest(1)-1.96*sqrt(C(1,1)) pest(1)-1.96*sqrt(C(1,1))],[0.5 0.7],'k--');
plot([pest(1)+1.96*sqrt(C(1,1)) pest(1)+1.96*sqrt(C(1,1))],[0.5 0.7],'k--');
%
% Comment these out to get the 95% CI for the other parameter. 
%
%plot([0.0017 0.0025],[pest(2)-1.96*sqrt(C(2,2)) pest(2)-1.96*sqrt(C(2,2))],'k--');
%plot([0.0017 0.0025],[pest(2)+1.96*sqrt(C(2,2)) pest(2)+1.96*sqrt(C(2,2))],'k--');
%print -deps slugcontour2.eps