Hi!
I am doing a polynomial fit to weighted data (measurements with
uncertainties) and need to obtain error estimates for the derived fit
(not error estimates for the fit parameters, which are commonly
returned by many fitting procedures, but error estimates for each
y-point at 1-sigma level). The standard function POLYFITW was found to
be the only one providing this quantity in the output keyword YBAND, in
fact, it's exactly what I need.
BUT the obtained from it YBAND seem to be kind of absurd.
When I expected the code I noticed that in the YBAND calculation, which
is done from the covariance matrix (weighted matrix) and the standard
deviation, the latter is not weighted, while in my opinion it should.
Namely, here is an extract from the POLYFITW code:
VAR = TOTAL((YFIT-Y)^2 )/(N-M) ; VARIANCE ESTIMATE, UNBIASED
SIGMA=SQRT(VAR)
YBAND = FLTARR(N) + A[0,0]
Z=FLTARR(N)+1.
FOR P=1,2*NDEGREE DO BEGIN ; COMPUTE CORRELATED ERROR ESTIMATES ON Y
Z = Z*X ; Z IS NOW X^P
SUM = 0.
FOR J=0 > (P - NDEGREE), NDEGREE < P DO SUM = SUM + A[J,P-J]
YBAND = YBAND + SUM * Z ; ADD IN ALL THE ERROR SOURCES
ENDFOR ; END OF P LOOP
YBAND = YBAND*VAR
YBAND = SQRT( YBAND )
What I did - I just weighted the variance VAR = TOTAL(WEIGHT*(YFIT-Y)^2
)/(N-M) and left the rest.
Have anyone used POLYFITW for weighted data? Any comments? Is is a bug
in the function or in my head?
I appreciate any links to the mathematical recipes of YBAND
calculation, I could not find any.
Yours sincerely,
Maria Loukitcheva
|
