;+ ; NAME: ; MPPROPERR ; ; AUTHOR: ; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 ; craigm@lheamail.gsfc.nasa.gov ; UPDATED VERSIONs can be found on my WEB PAGE: ; http://cow.physics.wisc.edu/~craigm/idl/idl.html ; ; PURPOSE: ; Propagate fitted model uncertainties to measured data points ; ; MAJOR TOPICS: ; Curve and Surface Fitting ; ; CALLING SEQUENCE: ; YCOVAR = MPPROPERR(BEST_FJAC, PCOVAR, PFREE_INDEX, [/DIAGONAL]) ; ; DESCRIPTION: ; ; MPPROPERR propagates the parameter uncertainties of a fitted ; model to provide estimates of the model uncertainty at each ; measurement point. ; ; When fitting a model to data with uncertainties, the parameters ; will have estimated uncertainties. In fact, the parameter ; variance-covariance matrix indicates the estimated uncertainties ; and correlations between parameters. These uncertainties and ; correlations can, in turn, be used to estimate the "error in the ; model" for each measurement point. In a sense, this quantity also ; reflects the sensitivity of the model to each data point. ; ; The algorithm used by MPPROPERR uses standard propagation of error ; techniques, assuming that errors are small. The input values of ; MPPROPERR should be found from the output keywords of MPFIT or ; MPFITFUN, as documented below. ; ; The user has a choice whether to compute the *full* ; variance-covariance matrix or not, depending on the setting of the ; DIAGONAL keyword. The full matrix is large, and indicates the ; correlation the sampled model function between each measurement ; point and every other point. The variance terms lie on the ; diagonal, and the covariance terms are on the off-diagonal. ; ; Usually however, the user will want to set /DIAGONAL, which only ; returns the "diagonal" or variance terms, which represent the ; model "uncertainty" at each measurement point. The /DIAGONAL ; setting only controls the amount of data returned to the user. ; the full *parameter* covariance matrix is always used to compute ; the output regardless of the setting for /DIAGONAL. ; ; When using MPPROPERR, keep in mind the following dimensions of ; the problem: ; NPOINTS - number of measurement points ; NPAR - total number of fit parameters ; NFREE - number of *free* fit parameters ; ; The inputs to this function are: ; BEST_FJAC - the partial derivative matrix, or Jacobian matrix, ; as estimated by MPFIT or MPFITFUN (see below), ; which has dimensions of ARRAY(NPOINTS,NFREE). ; PCOVAR - the parameter covariance matrix, as estimated by MPFIT ; or MPFITFUN (see below), which has dimensions of ; ARRAY(NPAR,NPAR). ; PFREE_INDEX - an index array which describes which of the ; parameter set were variable, as returned by MPFIT ; or MPFITFUN. Of the total parameter set PARMS, ; only PARMS[PFREE_INDEX] were varied by MPFIT. ; ; There are special considerations about the values returned by ; MPPROPERR. First, if a parameter is touching a boundary ; limit when the fit is complete, then it will be marked as having ; variance and covariance of zero. To avoid this situation, one can ; re-run MPFIT or MPFITFUN with MAXITER=0 and boundary limits ; disabled. This will permit MPFIT to estimate variance and ; covariance for all parameters, without allowing them to actually ; vary during the fit. ; ; Also, it is important to have a quality parameter covariance ; matrix PCOVAR. If the matrix is singular or nearly singular, then ; the measurement variances and covariances will not be meaningful. ; It helps to parameterize the problem to minimize parameter ; covariances. Also, consider fitting with double precision ; quantities instead of single precision to minimize the chances of ; round-off error creating a singular covariance matrix. ; ; IMPORTANT NOTE: the quantities returned by this function are the ; *VARIANCE* and covariance. If the user wishes to compute ; estimated standard deviation, then one should compute ; SQRT(VARIANCE). (see example below) ; ; INPUTS: ; ; BEST_FJAC - the Jacobian matrix, as estimated by MPFIT/MPFITFUN ; (returned in keyword BEST_FJAC). This should be an ; array ARRAY(NPOINTS,NFREE) where NFREE is the number ; of free parameters. ; ; PCOVAR - the full parameter covariance matrix, as returned in the ; COVAR keyword of MPFIT/MPFITFUN. This should be an array ; ARRAY(NPAR,NPAR) where NPAR is the *total* number of ; parameters. ; ; RETURNS: ; ; The estimated uncertainty at each measurement point, due to ; propagation of errors. The dimensions depend on the value of the ; DIAGONAL keyword. ; DIAGONAL=1: returned value is ARRAY(NPOINTS) ; corresponding to the *VARIANCE* of the model ; function sampled at each measurment point ; **NOTE**: the propagated standard deviation would ; then be SQRT(RESULT). ; ; DIAGONAL=0: returned value is ARRAY(NPOINTS,NPOINTS) ; corresponding to the variance-covariance matrix of ; the model function, sampled at the measurement ; points. ; ; ; KEYWORD PARAMETERS: ; ; DIAGONAL - if set, then compute only the "diagonal" (variance) ; terms. If not set, then propagate the full covariance ; matrix for each measurement point. ; ; NAN - if set, then ignore NAN values in BEST_FJAC or PCOVAR ; matrices (they would be set to zero). ; ; PFREE_INDEX - index list of free parameters, as returned in the ; PFREE_INDEX keyword of MPFIT/MPFITFUN. This should ; be an integer array ARRAY(NFREE), such that ; parameters PARMS[PFREE_INDEX] were freely varied during ; the fit, and the remaining parameters were not. ; Thus it should also be the case that PFREE_INDEX ; indicates the rows and columns of the parameter ; covariance matrix which were allowed to vary freely. ; Default: All parameters will be considered free. ; ; ; EXAMPLE: ; ; ; First, generate some synthetic data ; npts = 200 ; x = dindgen(npts) * 0.1 - 10. ; Independent variable ; yi = gauss1(x, [2.2D, 1.4, 3000.]) ; "Ideal" Y variable ; y = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise ; sy = sqrt(1000.D + y) ; Poisson errors ; ; ; Now fit a Gaussian to see how well we can recover ; p0 = [1.D, 1., 1000.] ; Initial guess (cent, width, area) ; p = mpfitfun('GAUSS1', x, y, sy, p0, $ ; Fit a function ; best_fjac=best_fjac, pfree_index=pfree_index, /calc_fjac, $ ; covar=pcovar) ; ; Above statement calculates best Jacobian and parameter ; ; covariance matrix ; ; ; Propagate errors from parameter covariance matrix to estimated ; ; measurement uncertainty. The /DIAG call returns only the ; ; "diagonal" (variance) term for each measurement. ; ycovar = mpproperr(best_fjac, pcovar, pfree_index=pfree_index, /diagonal) ; ; sy_prop = sqrt(ycovar) ; Estimated sigma ; ; ; REFERENCES: ; ; MINPACK-1, Jorge More', available from netlib (www.netlib.org). ; "Optimization Software Guide," Jorge More' and Stephen Wright, ; SIAM, *Frontiers in Applied Mathematics*, Number 14. ; ; MODIFICATION HISTORY: ; Written, 2010-10-27, CM ; Updated documentation, 2011-06-26, CM ; ; $Id: mpproperr.pro,v 1.5 2011/12/22 02:08:22 cmarkwar Exp $ ;- ; Copyright (C) 2011, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- function mpproperr, fjac, pcovar, pfree_index=ifree, diagonal=diag, $ nan=nan, status=status, errmsg=errmsg COMPILE_OPT strictarr status = 0 szf = size(fjac) if szf[0] NE 2 then begin errmsg = 'ERROR: BEST_FJAC must be an NPOINTxNFREE array' return, !values.d_nan endif npoints = szf[1] ;; Number of measurement points nfree = szf[2] ;; Number of free parameters nfree1 = n_elements(ifree) if nfree1 EQ 0 then begin ifree1 = lindgen(nfree) endif else if nfree1 NE nfree then begin errmsg = 'ERROR: Dimensions of PFREE_INDEX and BEST_FJAC must match' return, !values.d_nan endif szc = size(pcovar) if szc[0] NE 2 then begin PCOVAR_BAD_DIMS: errmsg = 'ERROR: PCOVAR must be an NPARxNPAR array' return, !values.d_nan endif if szc[1] NE szc[2] then goto, PCOVAR_BAD_DIMS npar = szc[1] if npar LT nfree then begin errmsg = 'ERROR: size of PCOVAR array is smaller than PFREE_INDEX' return, !values.d_nan endif fjac1 = fjac ;; NOTE: if there are parts of the covariance matrix which are zero, ;; that is OK, since they will contribute nothing to the output. pcovar1 = (pcovar[ifree,*])[*,ifree] ;; Check for NAN values and, if requested, set them to zero. if keyword_set(nan) then begin wh = where(finite(pcovar1) EQ 0, ct) if ct GT 0 then pcovar1[wh] = 0 wh = where(finite(fjac1) EQ 0, ct) if ct GT 0 then fjac1[wh] = 0 endif if NOT keyword_set(diag) then begin ;; Pull out the full covariance matrix (using matrix notation) ycovar = (fjac # pcovar1) # transpose(fjac) endif else begin ;; Only pull out the variance (diagonal) terms, and optimize a ;; little so that we don't use all the memory. ycovar = 0 for i = 0, nfree-1 do begin for j = 0, nfree-1 do begin ycovar = ycovar + fjac[*,i]*fjac[*,j]*pcovar1[i,j] endfor endfor endelse return, ycovar end