《最小二乘线性拟合系数差.docx》由会员分享，可在线阅读，更多相关《最小二乘线性拟合系数差.docx（5页珍藏版）》请在三一办公上搜索。
1、最小二乘线性拟合系数差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
2、 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 ch
3、i-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 nr
4、util, ONLY : assert_eqUSE nr, ONLY : gammqIMPLICIT NONEREAL(SP), DIMENSION(:), INTENT(IN) : x,yREAL(SP), INTENT(OUT) : a,b,siga,sigb,chi2,qREAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) : sigINTEGER(I4B) : ndataREAL(SP) : sigdat,ss,sx,sxoss,sy,st2REAL(SP), DIMENSION(size(x), TARGET : tREAL(SP), DIMENS
5、ION(:), POINTER : wtif (present(sig) thenndata=assert_eq(size(x),size(y),size(sig),fit)wt=>twt(:)=1.0_sp/(sig(:)*2)ss=sum(wt(:)sx=dot_product(wt,x)sy=dot_product(wt,y)elsendata=assert_eq(size(x),size(y),fit)ss=real(size(x),sp)sx=sum(x)sy=sum(y)end ifsxoss=sx/sst(:)=x(:)-sxossif (present(sig) thent
6、(:)=t(:)/sig(:)b=dot_product(t/sig,y)elseb=dot_product(t,y)end ifst2=dot_product(t,t)b=b/st2a=(sy-sx*b)/sssiga=sqrt(1.0_sp+sx*sx/(ss*st2)/ss)sigb=sqrt(1.0_sp/st2)t(:)=y(:)-a-b*x(:)q=1.0if (present(sig) thent(:)=t(:)/sig(:)chi2=dot_product(t,t)if (ndata > 2) q=gammq(0.5_sp*(size(x)-2),0.5_sp*chi2)e
7、lsechi2=dot_product(t,t)sigdat=sqrt(chi2/(size(x)-2)siga=siga*sigdatsigb=sigb*sigdatend ifEND 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-sx2;sigy=sqrt( dyy/ (n-2) );siga=sqrt( sxx/delt )*sigy;sigb=sqrt(n/delt)*sigy;
