CURVEFIT

The CURVEFIT function uses a gradient-expansion algorithm to compute a non-linear least squares fit to a user-supplied function with an arbitrary number of parameters. The user-supplied function may be any non-linear function where the partial derivatives are known or can be approximated. Iterations are performed until the chi square changes by a specified amount, or until a maximum number of iterations have been performed. The CURVEVIT function returns a vector of values for the dependent variables, as fitted by the function fit.

This routine is written in the IDL language. Its source code can be found in the file


curvefit.pro

in the

lib

subdirectory of the IDL distribution.

Arguments

X

An n -element vector of independent variables.

Y

A vector of dependent variables, the same length as X .

Weights

For instrumental (Gaussian) weighting, set Weights _i = 1.0/standard_deviation( Y _i )² . For statistical (Poisson) weighting, Weights _i = 1.0/ Y _i . For no weighting, set Weights _i = 1.0.

A

A vector with as many elements as the number of terms in the user-supplied function, containing the initial estimate for each parameter. On return, the vector A contains the fitted model parameters. If A is double-precision, calculations are performed in double-precision arithmetic, otherwise they are performed in single-precision arithmetic.

Sigma

A named variable that will contain a vector of standard deviations for the elements of the output vector A .

Keywords

CHISQ

Set this keyword equal to a named variable that will contain the value of chi-squared.

FUNCTION_NAME

Use this keyword to specify the name of the function to fit. If this keyword is omitted, CURVEFIT assumes that the IDL procedure FUNCT is to be used. If FUNCT is not already compiled, IDL compiles the function from the file funct.pro , located in the lib subdirectory of the IDL distribution. FUNCT evaluates the sum of a Gaussian and a second-order polynomial.

The function to be fit must be written as an IDL procedure and compiled prior to calling CURVEFIT. The procedure must accept values of X (the independent variable), and A (the fitted function's initial parameter values). It must return values for F (the function's value at X ), and optionally PDER (a 2D array of partial derivatives).

See the Example section below for an example function.

ITER

Set this keyword equal to a named variable that will contain the actual number of iterations performed.

ITMAX

Set this keyword to specify the maximum number of iterations. The default value is 20.

NODERIVATIVE

If this keyword is set, the routine specified by the FUNCTION_NAME keyword will not be requested to provide partial derivatives. The partial derivatives will be estimated by CURVEFIT using forward differences. If analytical derivatives are available they should always be used.

TOL

Use this keyword to specify the desired convergence tolerance. The routine returns when the relative decrease in chi-squared is less than TOL in one iteration. The default value is 1.0 ¥ 10^-3 .

Example

Fit a function of the form F(x) = a * exp(b*x) + c to sample pairs contained in arrays X and Y . The partial derivatives are easily computed symbolically:

df/da = exp(b*x)

df/db = a * x * exp(b*x)

df/dc = 1.0

First, define a procedure to return F(x) and the partial derivatives, given X . Note that A is an array containing the values a , b , and c .

PRO gfunct, X, A, F, pder

bx = EXP(A[1] * X)

F = A[0] * bx + A[2]

IF N_PARAMS() GE 4 THEN $ ; If the procedure is called with four parameters, calculate the partial derivatives.

pder= [[bx], [A[0] * X * bx], [replicate(1.0, N_ELEMENTS(X))]][

END

Compute the fit to the function we have just defined. First, define the independent and dependent variables:

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]

weights = 1.0/Y ; Define a vector of weights.

A = [10.0,-0.1,2.0] ; Provide an initial guess of the function's parameters.

yfit = CURVEFIT(X, Y, weights, A, SIGMA, FUNCTION_NAME='gfunct')
; Compute the parameters.

PRINT, 'Function parameters: ', A ; Print the parameters returned in A.

IDL prints:

Function parameters: 9.91120 -0.100883 2.07773

Thus, the function that best fits the data is:

f ( x ) = 9.91120( e ^-0.100883x ) + 2.07773