- 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 $