[AIT logo]

Institut für Astronomie und Astrophysik

Abteilung Astronomie

Sand 1, D-72076 Tübingen, Germany
[Uni logo]


JWCURVEFIT Source code in jwcurvefit.pro

JWCURVEFIT

Name
       JWCURVEFIT
Purpose
       Non-linear least squares fit to a function of an arbitrary
       number of parameters.  The function may be any non-linear
       function.  If available, partial derivatives can be calculated by
       the user function, else this routine will estimate partial derivatives
       with a forward difference approximation.
Category
       E2 - Curve and Surface Fitting.
Calling Sequence
       Result = JWCURVEFIT(X, Y, W, A, SIGMAA, FUNCTION_NAME = name, $
                         ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
Input Parameters
       X:  A row vector of independent variables.  This routine does
               not manipulate or use values in X, it simply passes X
               to the user-written function.
       Y:  A row vector containing the dependent variable.
       W:  A row vector of weights, the same length as Y.
               For no weighting,
               w(i) = 1.0.
               For instrumental weighting,
               w(i) = 1.0/y(i), etc.
            comment jw: for REAL Chi^2 values, set w(i)=1./sigma(i)^2
                                                       =1./variance(i)
       A:  A vector, with as many elements as the number of terms, that
           contains the initial estimate for each parameter.  If A is double-
           precision, calculations are performed in double precision,
           otherwise they are performed in single precision.
Keyword Parameters
       FUNCTION_NAME:  The name of the function (actually, a procedure) to
       fit.  If omitted, "FUNCT" is used. The procedure must be written as
       described under RESTRICTIONS, below.
       ITMAX:  Maximum number of iterations. Default = 20.
       ITER:   The actual number of iterations which were performed
       TOL:    The convergence tolerance. The routine returns when the
               relative decrease in chi-squared is less than TOL in an
               interation. Default = 1.e-3.
       CHI2:   The value of chi-squared on exit
       NODERIVATIVE:   If this keyword is set then the user procedure will not
               be requested to provide partial derivatives. The partial
               derivatives will be estimated in CURVEFIT using forward
               differences. If analytical derivatives are available they
               should always be used.
       BOUNDS: 2D-Array providing lower and upper bounds for the
               parameters. The fitting algorithm will not use any
               parameters outside the bounds.
Output Parameters
       Returns a vector of calculated values.
       A:  A vector of parameters containing fit.
Optional Output
       Sigmaa:  A vector of standard deviations for the parameters in
                A. (JW: DO NOT USE AS VALUES ARE COMPLETELY OFF!!!)
       dof   :  number of degrees of freedom of fit
Common Blocks
       NONE.
Side Effects
       None.
Restrictions
       The function to be fit must be defined and called FUNCT,
       unless the FUNCTION_NAME keyword is supplied.  This function,
       (actually written as a procedure) must accept values of
       X (the independent variable), and A (the fitted function's
       parameter values), and return F (the function's value at
       X), and PDER (a 2D array of partial derivatives).
       For an example, see FUNCT in the IDL User's Libaray.
       A call to FUNCT is entered as:
       FUNCT, X, A, F, PDER
 where:
       X = Variable passed into CURVEFIT.  It is the job of the user-written
               function to interpret this variable.
       A = Vector of NTERMS function parameters, input.
       F = Vector of NPOINT values of function, y(i) = funct(x), output.
       PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
               PDER(I,J) = Derivative of function at ith point with
               respect to jth parameter.  Optional output parameter.
               PDER should not be calculated if the parameter is not
               supplied in call. If the /NODERIVATIVE keyword is set in the
               call to CURVEFIT then the user routine will never need to
               calculate PDER.
Procedure
       Copied from "CURFIT", least squares fit to a non-linear
       function, pages 237-239, Bevington, Data Reduction and Error
       Analysis for the Physical Sciences.
       "This method is the Gradient-expansion algorithm which
       combines the best features of the gradient search with
       the method of linearizing the fitting function."
       Iterations are performed until the chi square changes by
       only TOL or until ITMAX iterations have been performed.
       The initial guess of the parameter values should be
       as close to the actual values as possible or the solution
       may not converge.
 EXAMPLE:  Fit a function of the form f(x) = a * exp(b*x) + c to
       sample pairs contained in x and y.
       In this example, a=a(0), b=a(1) and c=a(2).
       The partials are easily computed symbolicaly:
               df/da = exp(b*x), df/db = a * x * exp(b*x), and df/dc = 1.0
               Here is the user-written procedure to return F(x) and
               the partials, given x:
       pro gfunct, x, a, f, pder       ; Function + partials
         bx = exp(a(1) * x)
         f= a(0) * bx + a(2)           ;Evaluate the function
         if N_PARAMS() ge 4 then $     ;Return partials?
               pder= [[bx], [a(0) * x * bx], [replicate(1.0, N_ELEMENTS(y))]]
       end
         x=findgen(10)                 ;Define indep & dep variables.
         y=[12.0, 11.0,10.2,9.4,8.7,8.1,7.5,6.9,6.5,6.1]
         w=1.0/y                       ;Weights
         a=[10.0,-0.1,2.0]             ;Initial guess
         yfit=curvefit(x,y,w,a,sigmaa,function_name='gfunct')
         print, 'Function parameters: ', a
         print, yfit
       end
Revision History
       Written, DMS, RSI, September, 1982.
       Does not iterate if the first guess is good.  DMS, Oct, 1990.
       Added CALL_PROCEDURE to make the function's name a parameter.
              (Nov 1990)
       12/14/92 - modified to reflect the changes in the 1991
            edition of Bevington (eq. II-27) (jiy-suggested by CreaSo)
       Mark Rivers, U of Chicago, Feb. 12, 1995
           - Added following keywords: ITMAX, ITER, TOL, CHI2, NODERIVATIVE
             These make the routine much more generally useful.
           - Removed Oct. 1990 modification so the routine does one iteration
             even if first guess is good. Required to get meaningful output
             for errors.
           - Added forward difference derivative calculations required for
             NODERIVATIVE keyword.
           - Fixed a bug: PDER was passed to user's procedure on first call,
             but was not defined. Thus, user's procedure might not calculate
             it, but the result was then used.
       Joern Wilms, Univ. Tuebingen, Institute for Astronomy, 1997/1998:
           - Added bounds parameter
           - Added routines to enable the freezing of individual
             parameters (set low and high bounds to same value)
             (necessary for the computation of meaningful errors, see
             subroutine fiterror)
           - made more stable by checking for Nan's (presence of
             which could result in infinite loops).
           - return degrees of freedom to enable us to compute
             meaningful chi^2 instead of reduced chi2 only.
           - added _extra message passing to the fit function
           - better bounds handling (could lead to infinite loops)
 NEW VERSION: wrapper to mpfit -- PLEASE USE MPFIT DIRECTLY IN FUTURE
 PROGRAMS!!!!!!!!!!!
       Katja Pottschmidt, IAAT, 2001/02/11:
           - corrected calculation of dof for new version
            (nof=nof+1 instead of nof=nof-1 for frozen parameters)
   CVS Version 1.5: 2001/03/13 Joern Wilms, IAAT
           bug correction: compatibility code does not require
              function_name keyword anymore
   Version 1.6: handling of _error improved

Last modified by pro2html on 2005 January 04 at 16:36 UTC

[Home Page] [Software, Documentation] [IDL Documentation] [Quick Reference] [Feedback]

Jörn Wilms (wilms@astro.uni-tuebingen.de)
Updated automatically