pro lmcoeff,x,y,a,funcs=funcs,sig=sig,ifit=ifit,alpha=alpha,bvec=bvec,\$ chisq=chisq,yfunc=yfunc,poisson=poisson, _extra=e ;+ ;procedure lmcoeff ; used by LEVMARQ (formerly MRQMIN.F) to evaluate the linearized ; fitting matrix ALPHA and vector BVEC (BETA in Numerical Recipes) ; and calculate chi^2 ; ;syntax ; lmcoeff,x,y,a,funcs=funcs,sig=sig,ifit=ifit,alpha=alpha,bvec=bvec,\$ ; chisq=chisq,yfunc=yfunc,poisson=poisson ; ;parameters ; x [INPUT; required] data points ; y [INPUT; required] Y(X) ; * sizes of X and Y must match ; a [INPUT; required] parameters for user-supplied procedure ; ;keywords ; funcs [INPUT] name of user-defined procedure that takes as input X ; and A, and returns YMODEL(X;A), and the partial derivatives ; of the model function wrt A. Any procedure that was written ; to work with CURVEFIT or GHRS' WFIT will do. ; * default is set to X3MODEL ; * why is it a procedure and not a function? ; well, ask the person who wrote CURVEFIT. ; sig [INPUT; default=sqrt(abs(Y)+0.75)+1] standard deviations on Y ; ifit [INPUT] integer array of same size as A, with 0's indicating ; frozen parameters and 1's indicating thawed parameters ; * if IFIT does not match A for any reason, all parameters ; are assumed to be thawed ; alpha [OUTPUT] fitting matrix ; bvec [OUTPUT] fitting vector ; chisq [OUTPUT] chi^2, or ln(p(D|M)) if POISSON is set ; yfunc [OUTPUT] Y(X;A) -- function values ; poisson [INPUT] if set, returns the log of the poisson likelihood, ; p(D|M) instead of the chi-sq ; * do _not_ use this keyword for background-subtracted data!! ; _extra [JUNK] here only to prevent crashing the program ; ;subroutines ; LNPOISSON ; ;history ; (C) Copr. 1986-92 Numerical Recipes Software =\$j!]Y'1,). ; translated from MRQCOF.F to IDL by Vinay Kashyap (Oct98) ; now returns correct chisq if all params are frozen (VK; JanMM) ; added keyword POISSON and call to LNPOISSON (VK; Aug01) ;- ; initialize ok='ok' nx=n_elements(x) & ny=n_elements(y) & na=n_elements(a) ; usage if nx eq 0 then ok='X not defined' else \$ if ny eq 0 then ok='Y(X) not defined' else \$ if na eq 0 then ok='model parameters not defined' else \$ if nx ne ny then ok='Y and X do not match' if ok ne 'ok' then begin print,'Usage: lmcoeff,x,y,a,funcs=funcs,sig=sig,ifit=ifit,\$' print,' alpha=alpha,bvec=bvec,chisq=chisq,yfunc=yfunc' print,' subroutine of LEVMARQ to evaluate linearized fitting matrix' return endif ; keywords ; if not keyword_set(funcs) then funcs='x3model' ; ns=n_elements(sig) if ns ne ny then ysig=sqrt(abs(y)+0.75)+1. else ysig=sig sig2i=0.*ysig+1. oo=where(ysig gt 0,moo) & if moo gt 0 then sig2i(oo)=1./ysig(oo)^2 ; ni=n_elements(ifit) & kfit=intarr(na)+1 if ni ne na then message,\$ 'thawing all parameters because IFIT has gone bad',/info else kfit=ifit ; initialize ofit=where(kfit ne 0,nfit) if nfit eq 0 then begin message,'all parameters are frozen',/info ; get chisq and return chisq=0. call_procedure,funcs,x,a,yfunc,dyda, _extra=e if keyword_set(poisson) then begin tmp=fltarr(ny) for i=0L,ny-1L do tmp[i]=lnpoisson(y[i],yfunc[i],/chilike) chisq=total(tmp) endif else begin dy=y-yfunc & chisq=total(dy^2*sig2i) endelse return endif alpha=fltarr(nfit,nfit) & bvec=fltarr(nfit) & chisq=0. ; call the external function call_procedure,funcs,x,a,yfunc,dyda, _extra=e dy=y-yfunc ; generate the output matrices by summing over each data point for j=0,nfit-1 do begin jj=ofit(j) wt=dyda(*,jj)*sig2i for k=0,j do alpha(j,k)=total(wt*dyda(*,ofit(k))) ;fill in ALPHA bvec(j)=total(dy*wt) ;fill in BVEC if finite(bvec(j)) eq 0 then bvec(j)=0. endfor ; fill in symmetric side of ALPHA for j=1,nfit-1 do for k=0,j-1 do alpha(k,j)=alpha(j,k) ; get CHISQ if keyword_set(poisson) then begin tmp=fltarr(ny) for i=0L,ny-1L do tmp[i]=lnpoisson(y[i],yfunc[i],/chilike) chisq=total(tmp) endif else chisq=total(dy^2*sig2i) return end