function [c,m,err_c,err_m,r_squared] = linear_regression(x,y)


%performs linear regression on dataset [x,y], fitting to curve
%
%		y=c+mx
%
%The standard deviations of m and c are output together with
%the correlation coefficient (r-squared value)

sz=size(x);
ndata=sz(2);
sx=0.0;
sy=0.0;
st2=0.0;
m=0.0;

sx=sum(x);
sy=sum(y);
ss=1.0*ndata;
m=0.0;
c=0.0;

x_mean=mean(x);
for i=1:ndata,
  t=x(i)-x_mean;
  st2=st2+t*t;
  m=m+t*y(i);
end

m=m/st2;
c=(sy-sx*m)/ss;
sigc=sqrt((1.+sx*sx/(ss*st2))/ss);
sigm=sqrt(1./st2);
chi2=0.0;


for i=1:ndata,
  chi2=chi2+(y(i)-c-m*x(i))^2;
end

q=1.0;
sigdat=sqrt(chi2/(1.0*ndata-2.0));
err_c=sigc*sigdat;
err_m=sigm*sigdat;

y_mean=mean(y);
r_squared=1.-chi2/sum((y-y_mean).^2); 

disp('Slope')
m
disp('error in slope')
err_m
disp('offset')
c
disp('error in offset')
err_c
disp('R squared value for fit')
r_squared





