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.

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.

Syntax

Result = CURVEFIT( X, Y, Weights, A [, Sigma] [, CHISQ=variable] [, /DOUBLE] [, FITA=vector] [, FUNCTION_NAME=string] [, ITER=variable] [, ITMAX=value] [, /NODERIVATIVE] [, STATUS={0 | 1 | 2}] [, TOL=value] [, YERROR=variable] )

Return Value

Returns a vector of values for the dependent variables, as fitted by the function fit. If A is double-precision or if the DOUBLE keyword is set, calculations are performed in double-precision arithmetic, otherwise they are performed in single-precision arithmetic.

Arguments

X

An n-element vector containing the independent variable values. X may be of type integer, floating point, or double-precision floating-point.

For multivariate models, X can also be an n x m or an m x n array, where n is the number of values and m is the number of independent variables. In this case your fit function needs to be written to correctly interpret the dimensions of X.

Y

A vector of dependent variables. Y must have the same number of elements as F returned by the user-defined function.

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. If Weights is set to an undefined variable then no weighting will be used. Set Weights_i to 0 to not include it in the fit.

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.

Sigma

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

Note: If Weights is omitted, then you are assuming that your supplied model is the correct model for your data, and therefore, no independent goodness-of-fit test is possible. In this case, the values returned for the Sigma argument are multiplied by SQRT(CHISQ/(N-M)), where N is the number of points in X, and M is the number of coefficients. See Section 15.2 of Numerical Recipes in C (Second Edition) for details.

Keywords

CHISQ

Set this keyword to a named variable that will contain the value of the chi-square goodness-of-fit statistic, weighted by the measurement error:

where D_f = N - M is the degrees of freedom (N is the number of points in X, and M is the number of coefficients), W is the weighting value, and Yfit is the result.

DOUBLE

Set this keyword to force the computation to be done in double-precision arithmetic.

FITA

Set this keyword to a vector, with as many elements as A, which contains a zero for each fixed parameter, and a non-zero value for elements of A to fit. If not supplied, all parameters are taken to be non-fixed.

Note: If the NODERIVATIVE keyword is set then CURVEFIT will compute the partial derivatives by calling your function once for each parameter, using a value of that parameter that has been incremented by a small epsilon. This call will be made for each parameter including the fixed parameters. During the actual fitting procedure, IDL will ensure that the fixed parameters remain constant.

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).

The return value for F must have the same number of elements as Y. The return value for PDER (if supplied) must be a 2D array with dimensions
[N_ELEMENTS(Y), N_ELEMENTS(A)].

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.

STATUS

Set this keyword to a named variable that will contain an integer indicating the status of the computation. Possible return values are:

0	The computation was successful.
1	The computation failed. Chi-square was increasing without bounds.
2	The computation failed to converge in ITMAX iterations.

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 x 10^-3.

YERROR

Set this keyword to a named variable that will contain the standard error between YFIT and Y.

Examples

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 the procedure is called with four parameters, calculate the

;partial derivatives.

IF N_PARAMS() GE 4 THEN $

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]

;Define a vector of weights.

weights = 1.0/Y

;Provide an initial guess of the function’s parameters.

A = [10.0,-0.1,2.0]

;Compute the parameters.

yfit = CURVEFIT(X, Y, weights, A, SIGMA, FUNCTION_NAME='gfunct')

;Print the parameters returned in A.

PRINT, 'Function parameters: ', 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

Here is another example for CURVEFIT, which shows how to use multivariate input. The model function is:

Z = A X^B + C X Y + D Y^E

Click the following program:

PRO mymulticurvefit, XY, A, F, PDER
     X = XY[*,0]
     Y = XY[*,1]
     F = A[0]*X^A[1] + A[2]*X*Y + A[3]*Y^A[4]
     PDER = [[X^A[1]], [A[0]*X^A[1]*ALOG(X)], $
     [X*Y], [Y^A[4]], [A[3]*Y^A[4]*ALOG(Y)]]
END

Then click to run the following commands on the command line to show the plot and overplot output:

seed = 1
n = 1000
X = 2*RANDOMU(seed, n)
Y = 2*RANDOMU(seed, n)
XY = [[X], [Y]]
A = [-3, 2, 5, 3, 3]
noise = RANDOMN(seed, n)
Z = A[0]*X^A[1] + A[2]*X*Y + A[3]*Y^A[4] + noise

;Provide an initial guess of the function's parameters.
A = [1d, 1, 1, 1, 1]

;Compute the parameters.
yfit = CURVEFIT(XY, Z, weights, A, $
FUNCTION_NAME='mymulticurvefit')

PRINT, 'Function parameters: ', A

; Plot the original data as scattered points.
iPlot, X, Y, Z, /SCATTER

; Create a 2D surface of the final fit.
X1 = 2*FINDGEN(101)/100
X = REBIN(X1, 101, 101)
Y = TRANSPOSE(X)

Z = A[0]*X^A[1] + A[2]*X*Y + A[3]*Y^A[4]
iSurface, Z, X1, X1, /OVERPLOT

Version History

Pre 4.0	Introduced
5.6	Added YERROR keyword
6.0	Added FITA and STATUS keywords