PRO CFIT_ERROR, X, Y, A_NOM, FIT, SIGMAA, WEIGHTS=WEIGHTS, $ SFIT=SFIT, CHI2=CHISQR, DOUBLE=DOUBLE, ERRMSG=ERRMSG ;+ ; Project : SOHO - CDS ; ; Name : CFIT_ERROR ; ; Purpose : Calculate parameter errors from CFIT structure. ; ; Category : Fitting, Class3 ; ; Explanation : This procedure bypasses the CURVEFIT error analysis code. ; Although the CFIT program has the option of returning SIGMAA ; values, there is a possibility that the values returned by ; CURVEFIT may be incorrect if the fit did not improve the ; parameters. In particular, this can happen if CFIT is called ; with parameters that have already converged. ; ; Parameters which are held constant in the fit are treated as ; true constants, and are not used in determining the errors in ; fitted parameters. If a parameter which was held constant ; should be considered as a variable in the error analysis, then ; the FIT structure should be modified before this routine is ; called. ; ; Syntax : CFIT_ERROR, X, Y, A_NOM, FIT, SIGMAA ; ; Examples : YFIT = CFIT(X, Y, A_NOM, FIT, WEIGHTS=WT) ; CFIT_ERROR, X, Y, A_NOM, FIT, SIGMAA, WEIGHTS=WT, CHI2=CHI2 ; ; Inputs : X, Y = Data to be fitted ; FIT = Fit structure as used within CFIT ; ; Opt. Inputs : A_NOM = Parameters of the fit. If not defined, then the ; parameters are derived from the FIT structure. ; However, this is not recommended, since the ; parameters in the FIT structure may not represent the ; best fit to the data. ; ; Outputs : SIGMAA = Errors for each of the parameters in A_NOM. If this ; procedure fails, then SIGMAA will be set to -1. ; ; Opt. Outputs: A_NOM = Only changed if not initially defined. ; ; Keywords : WEIGHTS = The weights for each of the data points used in the ; fit. If not passed, then constant weighting is ; assumed, and the parameter errors are calculated ; based on the value of chi-squared. ; ; CHI2 = Returns the chi-squared value. ; ; SFIT = Structure from MAKE_SFIT_STC. Passing this structure ; speeds up execution when the same structure is used ; for a series of spectra. The best way to use this ; keyword is to pass it in undefined on the first call, ; and then reuse the structure on subsequent calls. ; ; DOUBLE = Set to force double precision calculation. Highly ; recommended. ; ; ERRMSG = If defined and passed, then any error messages will ; be returned to the user in this parameter rather than ; depending on the MESSAGE routine in IDL. If no ; errors are encountered, then a null string is ; returned. In order to use this feature, ERRMSG must ; be defined first, e.g. ; ; ERRMSG = '' ; CFIT_ERROR, ERRMSG=ERRMSG, ... ; IF ERRMSG NE '' THEN ... ; ; ; Calls : MAKE_SFIT_STC, EVAL_SFIT ; ; Common : None. ; ; Restrictions: The data must first be fit using either CFIT or XCFIT. ; ; Side effects: None. ; ; Prev. Hist. : Some of the code is taken from the routine LSTSQR. ; ; History : Version 1, 02-Dec-1999, William Thompson, GSFC ; Version 2, 15-Feb-2000, William Thompson, GSFC ; Use STATUS parameter from INVERT rather than DETERM ; function to catch singular matrices--more forgiving. ; Version 3, 26-Nov-2001, William Thompson, GSFC ; Fixed bug when only one parameter is fitted. ; Return SIGMAA=-1 if procedure fails. ; Added keyword DOUBLE ; ; Contact : WTHOMPSON ;- ; ON_ERROR, 2 ; IF N_PARAMS() LT 5 THEN BEGIN MESSAGE = 'Syntax: CFIT_ERROR, X, Y, A_NOM, FIT, SIGMAA' GOTO, HANDLE_ERROR ENDIF ; NDATA = N_ELEMENTS(X) IF N_ELEMENTS(Y) NE NDATA THEN BEGIN MESSAGE = 'Arrays X and Y must have the same number of points.' GOTO, HANDLE_ERROR ENDIF ; ; Determine whether WEIGHTS were passed. If not, the default is 1. ; WT_PASSED = N_ELEMENTS(WEIGHTS) NE 0 IF WT_PASSED THEN BEGIN IF N_ELEMENTS(WEIGHTS) NE NDATA THEN BEGIN MESSAGE = 'WEIGHTS must be the same size as X and Y' GOTO, HANDLE_ERROR ENDIF WT = WEIGHTS END ELSE BEGIN MESSAGE, /INFORMATIONAL, $ 'WARNING: No weights supplied - constant weights used' WT = REPLICATE(1.D0, NDATA) ENDELSE ; ; Extract the parameters, the linear transformation coefficients, and whether ; each parameter is a constant. ; FOR I = 0,N_TAGS(FIT)-1 DO BEGIN PP = FIT.(I).PARAM P = PP.VALUE A = PP.TRANS_A B = PP.TRANS_B C = PP.CONST IF I EQ 0 THEN BEGIN PARAM = P CONST = C TRANS_A = A TRANS_B = B END ELSE BEGIN PARAM = [PARAM, P] CONST = [CONST, C] TRANS_A = [TRANS_A, A] TRANS_B = [TRANS_B, B] ENDELSE ENDFOR ; ; If A_NOM was not passed, then set it to PARAM. Apply the linear ; transformation coefficients. ; IF N_ELEMENTS(A_NOM) EQ 0 THEN A_NOM = PARAM ELSE PARAM = A_NOM PARAM = PARAM * TRANS_A + TRANS_B WPAR = WHERE(CONST EQ 0, NPAR) ; ; Evaluate the function. Separate out those derivatives which are for ; non-constant parameters. ; IF N_ELEMENTS(SFIT) EQ 0 THEN SFIT = MAKE_SFIT_STC(FIT,DOUBLE=DOUBLE) EVAL_SFIT, X, PARAM, F, FDER, PRIVATE=SFIT FDER = FDER(*,WPAR) ; ; Calculate the matrix of partial derivatives. Extract the diagonal elements. ; ARRAY = TRANSPOSE(FDER) # ((WT # REPLICATE(1,NPAR)) * FDER) DIAG = INDGEN(NPAR)*(NPAR+1) ;Subscripts of diagonal elements. VEC = SQRT(ARRAY(DIAG)) W = WHERE(VEC EQ 0, N_FOUND) IF N_FOUND GT 0 THEN BEGIN MESSAGE = 'Zero partial derivative detected' GOTO, HANDLE_ERROR ENDIF ; ; Renormalize the correlation array based on the diagonal elements. Invert ; the array. ; ARRAY = ARRAY / (VEC # VEC) IF N_ELEMENTS(ARRAY) EQ 1 THEN BEGIN AINV = ARRAY * VEC^2 IF AINV(0) EQ 0 THEN BEGIN MESSAGE = 'A singular matrix was encountered' GOTO, HANDLE_ERROR ENDIF AINV = 1 / AINV END ELSE BEGIN AINV = INVERT(ARRAY, STATUS) / (VEC # VEC) IF STATUS EQ 1 THEN BEGIN MESSAGE = 'A singular matrix was encountered' GOTO, HANDLE_ERROR ENDIF ENDELSE ; ; Calculate chi-squared, and the errors in the parameters. ; NFREE = (NDATA - TOTAL(WT EQ 0) - NPAR) > 1 CHISQR = TOTAL(((F-Y)^2)*WT) / NFREE SIGSQR = ((AINV # TRANSPOSE(FDER))^2) # WT IF NOT WT_PASSED THEN SIGSQR = SIGSQR * CHISQR SIGMAA = 0*PARAM SIGMAA(WPAR) = SQRT(SIGSQR) / TRANS_A(WPAR) RETURN ; ; Error handling point. ; HANDLE_ERROR: IF N_ELEMENTS(PARAM) EQ 0 THEN SIGMAA = -1 ELSE BEGIN SIGMAA = 0*PARAM SIGMAA(WPAR) = -1 ENDELSE IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = MESSAGE ELSE MESSAGE, MESSAGE ; END