; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; UPDATED VERSIONs can be found on my WEB PAGE:
; Propagate fitted model uncertainties to measured data points
; MAJOR TOPICS:
; Curve and Surface Fitting
; CALLING SEQUENCE:
; YCOVAR = MPPROPERR(BEST_FJAC, PCOVAR, PFREE_INDEX, [/DIAGONAL])
; 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
; 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)
; 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
; 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
; 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.
; ; 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, $
; ; 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
; 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
status = 0
szf = size(fjac)
if szf NE 2 then begin
errmsg = 'ERROR: BEST_FJAC must be an NPOINTxNFREE array'
npoints = szf ;; Number of measurement points
nfree = szf ;; 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'
szc = size(pcovar)
if szc NE 2 then begin
errmsg = 'ERROR: PCOVAR must be an NPARxNPAR array'
if szc NE szc then goto, PCOVAR_BAD_DIMS
npar = szc
if npar LT nfree then begin
errmsg = 'ERROR: size of PCOVAR array is smaller than PFREE_INDEX'
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
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]