; $Id: curvefit.pro,v 1.5 1995/03/07 15:29:46 dave Exp $ function curvefit_c, x, y, w, a, sigmaa, Function_Name = Function_Name, $ itmax=itmax, iter=iter, tol=tol, chi2=chi2, $ noderivative=noderivative, $ a_expected=a_expected, a_weight=a_weight ; Copyright (c) 1988-1995, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;+ ; NAME: ; CURVEFIT_C ; ; 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 = CURVEFIT_C(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 CURVEFIT using forward ; differences. If analytical derivatives are available they ; should always be used. ; A_EXPECTED : The expected value of the parameters (INPUT). ; If omitted, initial guesses are used. ; A_WEIGHT : the inverse of the variances of the allowable difference between ; the expected and the final parameters ; ; 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 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(x))]] ; 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 ; ; 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, 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. ; Stein Vidar H. Haugan, Univ. of Oslo, 7 May 1996 ; - Detecting NaN errors and deadlock repetitions inside the ; REPEAT loop. ; Jongchul CHAE, 1996 ; - Introduce a regularizing term to prevent the degeneracy between ; parameters ; - Add new keyword paramters A_EXPECTED and A_weight ; - Eliminate forward difference derivative calulations ; - Vectorize the program so as to fit a set of data points ; at a same time ; ; ;- on_error,2 ;Return to caller if error ;Name of function to fit if n_elements(function_name) le 0 then function_name = "FUNCT" if n_elements(tol) eq 0 then tol = 1.e-2 ;Convergence tolerance if n_elements(itmax) eq 0 then itmax = 20 ;Maximum # iterations if n_elements(a_expected) eq 0 then a_expected=a if n_elements(a_weight) eq 0 then a_weight=a*0. ; Data Input Check ; if not (same_dim(a,a_expected) and same_dim(a, a_weight)) then begin print, 'Either A_EXPECTED or A_WEIGHT does not have the same ' print, 'data structure as A.' return,0 end if not (same_dim(x,w) and same_dim(x, y)) then begin print, 'Either Y or W does not have the same ' print, 'data structure as X.' return,0 end sa = size(a) sx = size(x) if not (sa(0) eq sx(0)) then begin print, 'A is not compatible with X.' return,0 end if sx(0) lt 1 or sx(0) gt 2 then begin print, 'Data should be 1-d or 2-d.' return,0 endif if sx(0) eq 1 then begin nset = 1 nterms = sa(1) ndata = sx(1) endif if sa(0) eq 2 then if sa(2) eq sx(2) then begin nset = sa(2) nterms = sa(1) ndata = sx(1) endif else begin print, ' Number of parameter sets should be equal to that of data sets.' return,0 endelse ; nfree =ndata-nterms ; Degrees of freedom if nfree le 0 then message, 'Curvefit - not enough data points.' flambda = replicate(0.01, nset) ;Initial lambda diag = lindgen(nterms)*(nterms+1) ; Subscripts of diagonal elements ;print, 'ndata=', ndata, ' nterm=', nterms,' nset=', nset b=a ; for iter = 1, itmax do begin ; Iteration loop call_procedure, function_name, x, a, yfit, pder h1 = total(w*(y-yfit)^2, 1)/nfree $ + total(a_weight*(a-a_expected)^2, 1)/nfree for k =0, nset-1 do begin beta = reform((y(*,k)-yfit(*,k))*w(*,k)) # reform(pder(*,*,k)) alpha = (transpose(reform(pder(*,*,k))) # $ (reform(w(*,k)) # (fltarr(nterms)+1)*reform(pder(*,*,k)))) alpha_new =alpha alpha_new(diag) = (alpha_new(diag)+ a_weight(*,k))*(1+flambda(k)) beta_new = beta beta_new =beta_new - a_weight(*,k)*(a(*,k)-a_expected(*,k)) ; ; Invert modified curvature matrix to find new parameters. array = invert(alpha_new) b(*,k) = a(*,k)+array # transpose(beta_new) ; New params endfor ;print, 'iterc =', iter, h1, flambda call_procedure, function_name, x, b, yfit ; Evaluate function h =total(w*(y-yfit)^2, 1)/nfree $ + total(a_weight*(a-a_expected)^2, 1)/nfree good = h lt h1 flambda = flambda*(good/10.+(1-good)*10.) > 1.0e-8 < 1.0e4 ss=where(good, count) if count ge 1 then begin a(*, ss) = b(*,ss) h1(*, ss)=h(*,ss) endif endfor ;iteration loop return,yfit ;return result END