-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathlinreg.m
58 lines (53 loc) · 1.67 KB
/
linreg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
function linreg(x,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LINREG(x,y) computes linear regression of y on x, i.e.,
% y = a + b*x
% and returns coefficients a and b with 95% uncertainity.
% Assumes x,y vectors have no NaN elements, and that x includes
% no measurement error. Uses TINV to estimate student-t value
% for confidence limits. See NLINREG when both x and y are
% considered as random variabiles.
%
% Ref: Draper and Smith, Applied Regression Analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ver 1.0: 11/25/96 (RB)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert to row vectors
[N,M]=size(x);
if M==1
x=x';y=y';
end
% 95% confidence limit student-t distribution
L=length(x);
P=.95;
p=.5.*(1 + P);
st=tinv(p,N-2);
mx=mean(x);
my=mean(y);
dx=x-mx;
dy=y-my;
reg=sum(dx.*dy)./sum(dx.^2);
b=sum(dx.*y)./sum(dx.^2);
a=my-b.*mx;
syx=(sum(dy.^2)-(sum(dx.*dy)^2)./sum(dx.^2))./(N-2);
syx=sqrt(syx);
aerr=st.*syx.*sqrt(((mx.^2)./sum(dx.^2)+1./N));
berr=st.*syx./sqrt(sum(dx.^2));
corr=corrcoef(x,y);
stdy=std(y-a-b*x);
A=[a aerr b berr corr(1,2)];
disp(['mean x = ',num2str(mx),'; mean y = ',num2str(my)])
disp(['std(y-a-bx) = ',num2str(stdy),'; N = ',int2str(N)])
disp(['a aerr = ',num2str(a),' +- ',num2str(aerr)])
disp(['b berr = ',num2str(b),' +- ',num2str(berr)])
disp(['corr coef = ',num2str(corr(1,2))])
disp(['(corr coef)^2 = ',num2str(corr(1,2)^2)])
xr=[min(x) max(x)];
yp=a+b*xr;
ypu=a+aerr + (b+berr)*xr;
ypl=a-aerr + (b-berr).*xr;
plot(x,y,'o',xr,yp,xr,ypu,'r:',xr,ypl,'r:')
xlabel('X')
ylabel('Y')
title('Least-squares fit of Y on X with maximum 95% CIs')
end