Chuck
unread,Jun 6, 2009, 4:26:08 PM6/6/09Sign in to reply to author
Sign in to forward
You do not have permission to delete messages in this group
Sign in to report message
Either email addresses are anonymous for this group or you need the view member email addresses permission to view the original message
to Sysquake
% linfit(x,y,flag) returns the coefficient vector of a linear fit
% to data points specified by column vectors x and y and the
% correlation coefficient of the fit. If flag is true, a plot
% of the data points and the least-squares line are generated
% with error bars from the data points to the line. If optional
% parameter flag is false (the default), no plot is generated.
function (coef,r) = linfit(x,y,flag)
if (nargin < 3)
flag=false;
end
a=[x x.^0];
b=y;
% find coefficients (c(1)*x + c(2))
ata=a'*a;
coef=ata^-1*(a'*b);
n=length(x);
% calculate the correlation coefficient
r=coef(1)*sqrt(det(ata))/sqrt(n*y'*y-(sum(y))^2);
yfit=polyval(coef,x);
% calculate the root-mean-square error
%error=y-yfit;
%rms_error=sqrt(mean(error.^2));
fprintf('y = %f x + %f\n',coef(1),coef(2))
fprintf('r = %f\n', r)
if flag
% plot data points
clf
plot(x,y,'bo')
label('x','y')
title('y as a function of x')
% plot least-squares line in 100 steps
xp=[x(1):(x(n)-x(1))/(100-1):x(n)];
yp=polyval(coef,xp);
plot(xp,yp)
% plot error bars from data points to line
for k=1:n
d=[x(k) x(k)];
e=[y(k) yfit(k)];
plot(d,e,'r-')
end
end