examp.m

这是在网上下的一个东东

%
% First, reset the state of the random number generator.
%
rand('state',0);
randn('state',0);
%
% Generate the x values
%
x=(25:3:97)';
n=length(x)
%
% and the y values.
%
ytrue=10*x+0.00*ones(size(x));
%
% and the noisy y values.
%
y=ytrue+0.05*randn(size(x)).*ytrue;
%
% Do the regression.  Note that \ here finds a least squares solution.
%
G=[ones(size(x)) x];
m=G\y;
%
% display the model.
%
disp('the model parameters are');
m
%
% Compute the modeled y's.
%
ymod=G*m;
%
% Compute the residuals and estimate s and C.
%
r=G*m-y;
s=sqrt(norm(r)^2/(n-2))
C=s^2*inv(G'*G)
%
% Compute confidence intervals for the parameters.
%
disp('half widths')
tinv(0.975,n-2)*sqrt(diag(C))
%
% Plots
%  1. Data and fitted line.
%  2. residuals.
%
figure(1)
clf;
bookfonts;
plot(x,y,'ko');
hold on
plot(x,ymod,'k');
xlabel('x');
ylabel('y');
%print -deps fitted.eps
figure(2)
clf;
bookfonts
plot(x,y-ymod,'ko');
xlabel('x');
ylabel('r');
%print -deps resid.eps
%
% Now, rescale the problem.
% 
W=inv(diag(y));
Gw=W*G;
yw=W*y;
mw=Gw\yw;
%
% display the model.
%
disp('the model parameters are');
mw
%
% Compute the modeled y's.
%
ymodw=Gw*mw;
%
% Compute the residuals and estimate s and C.
%
rw=Gw*mw-yw;
sw=sqrt(norm(rw)^2/(n-2))
Cw=sw^2*inv(Gw'*Gw)
%
% Compute confidence intervals for the parameters.
%
disp('half widths')
tinv(0.975,n-2)*sqrt(diag(Cw))
%
% Compute the modeled y's.
%
ymodw=G*mw;
%
% Plots
%  3. Data and fitted line.
%  4. residuals.
%
figure(3)
clf;
bookfonts;
plot(x,y,'ko');
hold on
plot(x,ymodw,'k');
xlabel('x');
ylabel('y');
%print -deps fittedsc.eps
figure(4)
clf;
bookfonts;
plot(x,rw,'ko');
xlabel('x');
ylabel('r');
%print -deps residsc.eps