File: mpproperr.pro

package info (click to toggle)
mpfit 1.85+2017.01.03-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 788 kB
  • sloc: python: 126; makefile: 16; sh: 13
file content (250 lines) | stat: -rw-r--r-- 10,061 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
;+
; 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