最小二乘线性拟合系数标准差

matlab中polyfit似乎没能给出拟合系数的方差？
《Numerical Recipes (Fortran Version)》Willian H. Press, Brain P. Flanny, Saul A. Teukolsky, Willam T. Vetterling有一章叫Modeling of Data，不同的edtion章节号可能会有变化，其中谈到least squares fitting中的a和b的uncertainty,也就是standard deviation标准差。用的是最大似然估计和 chi-square分布求得的。其子程序如下：
SUBROUTINE fit(x,y,a,b,siga,sigb,chi2,q,sig)
! Given a set of data points x(i),y(i), fit them to a straight line y=a+bx by minimizing chi-square.
! Returned are a, b and their respective probable uncertainties siga and sigb, the chi-square chi2.
! and the goodness-of-fit probability q (that the fit would have chi-square this large or larger.
! The sig is an optional input parameter. When the standard deviation of y(i) is unavailable, sig
! is absent and q is returned as 1.0 and the normalization of chi2 is to unit standard deviation on
! all points.
! chi2=( 1 - r**2 ) * sum( y(i) - mean(y) )

USE nrtype; USE nrutil, ONLY : assert_eq
USE nr, ONLY : gammq
IMPLICIT NONE
REAL(SP), DIMENSION(:), INTENT(IN) :: x,y
REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q
REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig
INTEGER(I4B) :: ndata
REAL(SP) :: sigdat,ss,sx,sxoss,sy,st2
REAL(SP), DIMENSION(size(x)), TARGET :: t
REAL(SP), DIMENSION(:), POINTER :: wt
if (present(sig)) then
ndata=assert_eq(size(x),size(y),size(sig),'fit')
wt=>t
wt(:)=1.0_sp/(sig(:)**2)
ss=sum(wt(:))
sx=dot_product(wt,x)
sy=dot_product(wt,y)
else
ndata=assert_eq(size(x),size(y),'fit')
ss=real(size(x),sp)
sx=sum(x)
sy=sum(y)
end if
sxoss=sx/ss
t(:)=x(:)-sxoss
if (present(sig)) then
t(:)=t(:)/sig(:)
b=dot_product(t/sig,y)
else
b=dot_product(t,y)
end if
st2=dot_product(t,t)
b=b/st2
a=(sy-sx*b)/ss
siga=sqrt((1.0_sp+sx*sx/(ss*st2))/ss)
sigb=sqrt(1.0_sp/st2)
t(:)=y(:)-a-b*x(:)
q=1.0
if (present(sig)) then
t(:)=t(:)/sig(:)
chi2=dot_product(t,t)
if (ndata > 2) q=gammq(0.5_sp*(size(x)-2),0.5_sp*chi2)
else
chi2=dot_product(t,t)
sigdat=sqrt(chi2/(size(x)-2))
siga=siga*sigdat
sigb=sigb*sigdat
end if
END SUBROUTINE fit

!==将它的思想转化为Matlab子程序linearfitstd(x,y)
function [siga,sigb]=linearfitstd(x,y)
%This function is aimming for Standard deivations
% for a and b which are got from the least square
% fitting using 'polyfit'.
p=polyfit(x,y,1);
y1=polyval(p,x);
n=length(x);
sx=sum(x);
sxx=sum(x.^2);
dyy=sum((y1-y).^2);
delt=n*sxx-sx^2;
sigy=sqrt( dyy/ (n-2) );
siga=sqrt( sxx/delt )*sigy;
sigb=sqrt(n/delt)*sigy;