function lstsqr,dat,funs,wt,nfun,rms,chisq,outp,type ; This routine does a weighted linear least-squares fit of the nfun functions ; contained in array funs to the data in array dat. The weights are given ; in array wt. Fit coefficients are returned. If arguments outp and type ; are given, then on return outp contains: ; type = 0 => the fitted function ; type = 1 => residuals around fit,in the sense (data - fit) ; type = 2 => ratio (data/fit) ; if outp is an argument but type is not given, type defaults to 0. ; on return, chisq contains the average of err^2*wt^2, so that wt is implicitly ; taken to be the reciprocal of the expected sigma at each data point. ; The technique used is to construct and solve the normal equations. ; Note that this technique will give garbage for ill-conditioned systems. ; get dimensions of things, make extended arrays for generating normal eqn ; matrix npr=n_params() s=size(dat) if (s(0) ne 1) then begin print,'bad dimension in lstsqr data' return,0. end nx=s(1) wte=rebin(wt,nx,nfun) datw=reform(dat,nx,1) datw=rebin(datw,nx,nfun)*wte funw=funs*wte ; make normal eqn matrix, rhs a=fltarr(nfun,nfun) rhs=rebin(funw*datw,1,nfun) rhs=reform(rhs,nfun) for i=0,nfun-1 do begin for j=0,nfun-1 do begin if(i ge j) then begin prod=rebin(funw(*,i)*funw(*,j),1) a(i,j)=prod a(j,i)=prod end end end ; solve equations by lu decomposition ludcmp,a,index,d lubksb,a,index,rhs ; Make fit function, rms outpt=reform(rhs,1,nfun) outpt=rebin(outpt,nx,nfun)*funs outpt=rebin(outpt,nx)*nfun dif=dat-outpt s=where(wt gt 0) wt2=wt^2 dif2=dif^2 rmst=sqrt(total(dif2(s))/n_elements(s)) ndegfree=n_elements(s)-nfun chisq=total(dif2(s)*wt(s))/ndegfree ; if rms is defined, set its value if (npr ge 5) then rms=rmst if (npr le 5) then return,rhs ; make final outp,depending on type if (npr lt 8) then typ = 0 else typ = type if (typ eq 0) then begin outp=outpt return,rhs end if (typ eq 1) then begin outp=dif return,rhs end if(typ eq 2) then begin outp=dat/outpt return,rhs end end