; $Id: ncurvefit.pro,v 1.2 1993/10/26 23:44:03 doug Exp $ function ncurvefit, x, y, w, a, sigmaa, Function_Name = Function_Name, $ itmax=itmax, iter=iter, tol=tol, chi2=chi2, $ noderivative=noderivative ;+ ; NAME: ; nCURVEFIT ; ; 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 = nCURVEFIT(X, Y, W, A, SIGMAA, FUNCTION_NAME = name, $ ; ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE) ; ; INPUTS: ; X: A row vector of independent variables. ; ; Y: A row vector of dependent variable, the same length as x. ; ; W: A row vector of weights, the same length as x and y. ; For no weighting, ; w(i) = 1.0. ; For instrumental weighting, ; w(i) = 1.0/y(i), etc. ; ; 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. ; ; KEYWORDS: ; 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 nCURVEFIT using forward ; differences. If analytical derivatives are available they ; should always be used. ; ; OUTPUTS: ; Returns a vector of calculated values. ; A: A vector of parameters containing fit. ; ; OPTIONAL OUTPUT PARAMETERS: ; Sigmaa: A vector of standard deviations for the parameters in A. ; ; 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 = Vector of NPOINT independent variables, input. ; A = Vector of NTERMS function parameters, input. ; F = Vector of NPOINT values of function, y(i) = funct(x(i)), 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 nCURVEFIT 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: ; pro gfunct, x, a, f, pder ; f=a(0)*exp(a(1)*x)+a(2) ; pder=findgen(10, 3) ; end ; ; pro fit_curve ; x=float(indgen(10)) ; 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 ; a=[10.0,-0.1,2.0] ; yfit=ncurvefit(x,y,w,a,sigmaa,function_name='gfunct') ; print, yfit ; end ; ; MODIFICATION 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, 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. ; ;- on_error,2 ;Return to caller if error ;Name of function to fit if n_elements(function_name) le 0 then function_name = "FUNCT" ;Convergence tolerance if n_elements(tol) eq 0 then tol = 1.e-3 ;Maximum number of iterations if n_elements(itmax) eq 0 then itmax = 20 a = 1.*a ; Make params floating or double ; Set flag if calculations are to be done in double precision t = size(a) if (t(n_elements(t)-2) eq 5) then double=1 else double=0 ; If we will be estimating partial derivatives then compute machine ; precision if keyword_set(NODERIVATIVE) then begin if (double) then res = nr_machar(/DOUBLE) else res=nr_machar() eps = sqrt(res.eps) endif nterms = n_elements(a) ; # of parameters nfree = (n_elements(y)