GAUSSFIT

The GAUSSFIT function computes a non-linear least-squares fit to a function f (x) with from three to six unknown parameters. f (x) is a linear combination of a Gaussian and a quadratic; the number of terms is controlled by the keyword parameter NTERMS.

This routine is written in the IDL language. Its source code can be found in the file gaussfit.pro in the lib subdirectory of the IDL distribution.

Example

The following example shows how to use GAUSSFIT to fit to four different functions, with NTERMS=3,4,5,6. To simulate actual data, random noise has been added to each function. Save the code in a file named "ex_gaussfit.pro".

; Define the independent variable.

n = 101

x = (FINDGEN(n)-(n/2))/4

; Define the coefficients.

a = [4.0, 1.0, 2.0, 1.0, 0.25, 0.01]

print, 'Expect: ', a

z = (x - a[1])/a[2] ; Gaussian variable

seed = 123321 ; Pick a starting seed value

for nterms=3,6 do begin

; Define the dependent variables. Start with random noise.

y = 0.4*RANDOMN(seed, n)

; Use a switch statement so we fall through to each term.

switch nterms of

6: y = y + a[5]*x^2

5: y = y + a[4]*x

4: y = y + a[3]

3: y = y + a[0]*exp(-z^2/2)

endswitch

; Fit the data to the function, storing coefficients in

; coeff:

yfit = GAUSSFIT(x, y, coeff, NTERMS=nterms)

print, 'Result: ', coeff[0:nterms-1]

; Plot the original data and the fitted curve:

p1 = PLOT(x, y, TITLE='nterms='+STRTRIM(nterms,2), $

LAYOUT=[4,1,nterms-2], CURRENT=(nterms gt 3), $

DIMENSIONS=[800,200], MARGIN=[0.1,0.2,0.1,0.2])

p2 = PLOT(x, yfit, THICK=2, /OVERPLOT)

endfor

end

When this program is compiled and executed, IDL generates the plots at the top of this topic, and then prints the following results:

IDL> .run ex_gaussfit

Expect: 4.00000 1.00000 2.00000 1.00000 0.250000 0.0100000

Result: 3.91262 1.05290 1.99806

Result: 4.04331 1.03996 2.06383 0.983630

Result: 4.24642 1.04692 1.95604 0.959738 0.242130

Result: 3.74969 0.934981 2.02861 1.10497 0.251669 0.00969010

Syntax

Result = GAUSSFIT( X, Y [, A] [, CHISQ=variable] [, ESTIMATES=array] [, MEASURE_ERRORS=vector] [, NTERMS=integer{3 to 6}] [, SIGMA=variable] [, YERROR=variable])

Return Value

The result is a vector with the same length as X that contains the function f(x) evaluated at each X value. The coefficients of the fit are returned in the argument A.

Arguments

X

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

Y

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

A

A named variable that will contain the coefficients A of the fit.

Keywords

CHISQ

Set this keyword to a named variable that will contain the value of the reduced chi-square goodness-of-fit statistic.

ESTIMATES

Set this keyword equal to an array of starting estimates for the parameters of the equation. If the NTERMS keyword is specified, the ESTIMATES array should have NTERMS elements. If NTERMS is not specified, the ESTIMATES array should have six elements. If the ESTIMATES array is not specified, estimates are calculated by first subtracting a constant (if NTERMS >= 4) or a linear term (if NTERMS >= 5), and then forming a simple estimate of the Gaussian coefficients.

Note: If ESTIMATES[2] (the Gaussian width) is zero, then GAUSSFIT will automatically compute a best guess for the starting value for this term.

MEASURE_ERRORS

Set this keyword to a vector containing standard measurement errors for each point Y[i]. This vector must be the same length as X and Y.

Note: For Gaussian errors (e.g., instrumental uncertainties), MEASURE_ERRORS should be set to the standard deviations of each point in Y. For Poisson or statistical weighting, MEASURE_ERRORS should be set to SQRT(Y).

NTERMS

Set this keyword to an integer value between three and six to specify the function to be used for the fit. The values correspond to the functions shown below. In all cases:

NTERMS=6

NTERMS=5

NTERMS=4

NTERMS=3

For all of these cases, z is defined as,

NTERMS=6 is the default setting. Here, A₀ is the height of the Gaussian, A₁ is the center of the Gaussian, A₂ is the width (the standard deviation) of the Gaussian, A₃ is the constant term, A₄ is the linear term, and A₅ is the quadratic term.

Tip: The full width at half maximum (FWHM) of the Gaussian may be computed as 2*SQRT(2*ALOG(2))*A2.

SIGMA

Set this keyword to a named variable that will contain the 1-sigma error estimates of the returned parameters.

Note: If MEASURE_ERRORS is omitted, then you are assuming that a polynomial is the correct model for your data, and therefore, no independent goodness-of-fit test is possible. In this case, the values returned in SIGMA 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.

YERROR

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

Version History

Original	Introduced
5.6	Added CHISQ, SIGMA, and YERROR keywords
6.0	Added MEASURE_ERRORS keyword