For detail information please refer to the IDL reference manual.
COMFIT - Six common types of gradient-expansion least-square fit
EXPONENTIAL Y = a0 * a1^x + a2
GEOMETRIC Y = a0 * x^a1 + a2
GOMPERTZ Y = a0 * a1^(a2*x) + a3
HYPERBOLIC Y = 1./(a0 + a1*x)
LOGISTIC Y = 1./(a0 * a1^x + a2)
LOGSQUARE Y = a0 + a1*alog10(x) + a2 * alog10(x)^2
CURVEFIT - Least-square fit to an arbitrary non-linear function
GAUSSFIT - Least-square fit to
F(x) = A0*EXP(-((x-A1)/A2)^2/2) + A3 + A4*x + A5*x^2
LADFIT - Least-absolute-deviation fit to Y = A + Bx
LINFIT - Minimize-chi-square error fit to Y = A + Bx
LMFIT - Non-linear least-square fit to n arbitrary non-linear function
POLYFITW - Weighted least-square polynomial fit
POLY_FIT - Least-square polynomial fit
REGRESS - Multiple linear regression fit
SVDFIT - Least-square fit to user-supplied/built-in/legendre polynomial
NAME:
COMFIT
PURPOSE:
This function fits the paired data {X(i), Y(i)} to one of six common
types of appoximating models using a gradient-expansion least-squares
method. The result is a vector containing the model parameters.
CATEGORY:
Statistics.
CALLING SEQUENCE:
Result = COMFIT(X, Y, A)
INPUTS:
X: An n-element vector of type integer, float or double.
Y: An n-element vector of type integer, float or double.
A: A vector of initial estimates for each model parameter. The length
of this vector depends upon the type of model selected.
KEYWORD PARAMETERS:
EXPONENTIAL: If set to a non-zero value, the parameters of the
exponential model are computed. Y = a0 * a1^x + a2.
GEOMETRIC: If set to a non-zero value, the parameters of the
geometric model are computed. Y = a0 * x^a1 + a2.
GOMPERTZ: If set to a non-zero value, the parameters of the
Gompertz model are computed. Y = a0 * a1^(a2*x) + a3.
HYPERBOLIC: If set to a non-zero value, the parameters of the
hyperbolic model are computed. Y = 1./(a0 + a1*x)
LOGISTIC: If set to a non-zero value, the parameters of the
logistic model are computed. Y = 1./(a0 * a1^x + a2)
LOGSQUARE: If set to a non-zero value, the parameters of the
logsquare model are computed.
Y = a0 + a1*alog10(x) + a2 * alog10(x)^2
SIGMA: Use this keyword to specify a named variable which
returns a vector of standard deviations for the computed
model parameters.
WEIGHTS: An n-element vector of weights. If the WEIGHTS vector
is not specified, an n-element vector of 1.0s is used.
YFIT: Use this keyword to specify a named variable which
returns n-element vector of y-data corresponding to the
computed model parameters.
EXAMPLE:
Define two n-element vectors of paired data.
x = [2.27, 15.01, 34.74, 36.01, 43.65, 50.02, 53.84, 58.30, 62.12, $
64.66, 71.66, 79.94, 85.67, 114.95]
y = [5.16, 22.63, 34.36, 34.92, 37.98, 40.22, 41.46, 42.81, 43.91, $
44.62, 46.44, 48.43, 49.70, 55.31]
Define a 3-element vector of initial estimates for the logsquare model.
a = [1.5, 1.5, 1.5]
Compute the model parameters of the logsquare model, a(0), a(1), & a(2).
result = comfit(x, y, a, sigma = sigma, yfit = yfit, /logsquare)
The result should be the 3-element vector:
[1.42494, 7.21900, 9.18794]
REFERENCE:
APPLIED STATISTICS (third edition)
J. Neter, W. Wasserman, G.A. Whitmore
ISBN 0-205-10328-6
MODIFICATION HISTORY:
Written by: GGS, RSI, September 1994
(See /usr/local/rsi/idl_5/lib/comfit.pro)
NAME:
CURVEFIT
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(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
INPUTS:
X: A row vector of independent variables. This routine does
not manipulate or use values in X, it simply passes X
to the user-written function.
Y: A row vector containing the dependent variable.
Weights: A row vector of weights, the same length as Y.
For no weighting,
Weights(i) = 1.0.
For instrumental (Gaussian) weighting,
Weights(i)=1.0/sigma(i)^2
For statistical (Poisson) weighting,
Weights(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. Fitted parameters
are returned in A.
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 (obselete)
CHISQ: The value of reduced 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.
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 = Variable passed into CURVEFIT. It is the job of the user-written
function to interpret this variable.
A = Vector of NTERMS function parameters, input.
F = Vector of NPOINT values of function, y(i) = funct(x), 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 is adapted from:
Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear
Parameters", J. Soc. Ind. Appl. Math., Vol 11, no. 2, pp. 431-441,
June, 1963.
"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(f))]]
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]
Weights=1.0/y ;Weights
a=[10.0,-0.1,2.0] ;Initial guess
yfit=curvefit(x,y,Weights,a,sigma,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.
Steve Penton, RSI, June 1996.
- Changed SIGMAA to SIGMA to be consistant with other fitting
routines.
- Changed CHI2 to CHISQ to be consistant with other fitting
routines.
- Changed W to Weights to be consistant with other fitting
routines.
_ Updated docs regarding weighing.
(See /usr/local/rsi/idl_5/lib/curvefit.pro)
NAME: GAUSSFIT PURPOSE: Fit the equation y=f(x) where: F(x) = A0*EXP(-z^2/2) + A3 + A4*x + A5*x^2 and z=(x-A1)/A2 A0 = height of exp, A1 = center of exp, A2 = sigma (the width). A3 = constant term, A4 = linear term, A5 = quadratic term. Terms A3, A4, and A5 are optional. The parameters A0, A1, A2, A3 are estimated and then CURVEFIT is called. CATEGORY: ?? - fitting CALLING SEQUENCE: Result = GAUSSFIT(X, Y [, A]) INPUTS: X: The independent variable. X must be a vector. Y: The dependent variable. Y must have the same number of points as X. KEYWORD INPUTS: KEYWORD INPUTS: ESTIMATES = optional starting estimates for the parameters of the equation. Should contain NTERMS (6 if NTERMS is not provided) elements. NTERMS = Set NTERMS to 3 to compute the fit: F(x) = A0*EXP(-z^2/2). Set it to 4 to fit: F(x) = A0*EXP(-z^2/2) + A3 Set it to 5 to fit: F(x) = A0*EXP(-z^2/2) + A3 + A4*x OUTPUTS: The fitted function is returned. OPTIONAL OUTPUT PARAMETERS: A: The coefficients of the fit. A is a three to six element vector as described under PURPOSE. COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: The peak or minimum of the Gaussian must be the largest or smallest point in the Y vector. PROCEDURE: The initial estimates are either calculated by the below procedure or passed in by the caller. Then the function CURVEFIT is called to find the least-square fit of the gaussian to the data. Initial estimate calculation: If the (MAX-AVG) of Y is larger than (AVG-MIN) then it is assumed that the line is an emission line, otherwise it is assumed there is an absorbtion line. The estimated center is the MAX or MIN element. The height is (MAX-AVG) or (AVG-MIN) respectively. The width is found by searching out from the extrema until a point is found less than the 1/e value. MODIFICATION HISTORY: DMS, RSI, Dec, 1983. DMS, RSI, Jun, 1995, Added NTERMS keyword. Result is now float if Y is not double. DMS, RSI, Added ESTIMATES keyword.
(See /usr/local/rsi/idl_5/lib/gaussfit.pro)
NAME:
LADFIT
PURPOSE:
This function fits the paired data {X(i), Y(i)} to the linear model,
y = A + Bx, using a "robust" least absolute deviation method. The
result is a two-element vector containing the model parameters, A
and B.
CATEGORY:
Statistics.
CALLING SEQUENCE:
Result = LADFIT(X, Y)
INPUTS:
X: An n-element vector of type integer, float or double.
Y: An n-element vector of type integer, float or double.
KEYWORD PARAMETERS:
ABSDEV: Use this keyword to specify a named variable which returns the
mean absolute deviation for each data-point in the y-direction.
DOUBLE: If set to a non-zero value, computations are done in double
precision arithmetic.
EXAMPLE:
Define two n-element vectors of paired data.
x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $
2.95, 2.18, 3.72, 5.26]
y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $
-0.78, -2.61, 0.31, 1.74]
Compute the model parameters, A and B.
result = ladfit(x, y, absdev = absdev)
The result should be the two-element vector:
[-3.15301, 0.930440]
The keyword parameter should be returned as:
absdev = 0.636851
REFERENCE:
Numerical Recipes, The Art of Scientific Computing (Second Edition)
Cambridge University Press
ISBN 0-521-43108-5
MODIFICATION HISTORY:
Written by: GGS, RSI, September 1994
Modified: GGS, RSI, July 1995
Corrected an infinite loop condition that occured when
the X input parameter contained mostly negative data.
Modified: GGS, RSI, October 1996
If least-absolute-deviation convergence condition is not
satisfied, the algorithm switches to a chi-squared model.
Modified keyword checking and use of double precision.
Modified: GGS, RSI, November 1996
Fixed an error in the computation of the median with
even-length input data. See EVEN keyword to MEDIAN.
(See /usr/local/rsi/idl_5/lib/ladfit.pro)
NAME:
LINFIT
PURPOSE:
This function fits the paired data {X(i), Y(i)} to the linear model,
y = A + Bx, by minimizing the chi-square error statistic. The result
is a two-element vector containing the model parameters,[A,B].
CATEGORY:
Statistics.
CALLING SEQUENCE:
Result = LINFIT(X, Y)
INPUTS:
X: An n-element vector of type integer, float or double.
Y: An n-element vector of type integer, float or double.
KEYWORD PARAMETERS:
CHISQ: Use this keyword to specify a named variable which returns the
chi-square error statistic as the sum of squared errors between
Y(i) and A + BX(i). If individual standard deviations are
supplied, then the chi-square error statistic is computed as
the sum of squared errors divided by the standard deviations.
DOUBLE: If set to a non-zero value, computations are done in double
precision arithmetic.
PROB: Use this keyword to specify a named variable which returns the
probability that the computed fit would have a value of CHISQR
or greater. If PROB is greater than 0.1, the model parameters
are "believable". If PROB is less than 0.1, the accuracy of the
model parameters is questionable.
SDEV: An n-element vector of type integer, float or double that
specifies the individual standard deviations for {X(i), Y(i)}.
SIGMA: Use this keyword to specify a named variable which returns a
two-element vector of probable uncertainties for the model par-
ameters, [SIG_A,SIG_B].
EXAMPLE:
Define two n-element vectors of paired data.
x = [-3.20, 4.49, -1.66, 0.64, -2.43, -0.89, -0.12, 1.41, $
2.95, 2.18, 3.72, 5.26]
y = [-7.14, -1.30, -4.26, -1.90, -6.19, -3.98, -2.87, -1.66, $
-0.78, -2.61, 0.31, 1.74]
Define a vector of standard deviations with a constant value of 0.85
sdev = replicate(0.85, n_elements(x))
Compute the model parameters, A and B.
result = linfit(x, y, chisq = chisq, prob = prob, sdev = sdev)
The result should be the two-element vector:
[-3.44596, 0.867329]
The keyword parameters should be returned as:
chisq = 11.4998, prob = 0.319925
REFERENCE:
Numerical Recipes, The Art of Scientific Computing (Second Edition)
Cambridge University Press
ISBN 0-521-43108-5
MODIFICATION HISTORY:
Written by: GGS, RSI, September 1994
LINFIT is based on the routines: fit.c, gammq.c, gser.c,
and gcf.c described in section 15.2 of Numerical Recipes,
The Art of Scientific Computing (Second Edition), and is
used by permission.
Modified: SVP, RSI, June 1996
Changed SIG_AB to SIGMA to be consistant with the other
fitting functions. Changed CHISQR to CHISQ in the docs
for the same reason. Note that the chisqr and the SIG_AB
keywords are left for backwards compatibility.
Modified: GGS, RSI, October 1996
Modified keyword checking and use of double precision.
Added DOUBLE keyword.
(See /usr/local/rsi/idl_5/lib/linfit.pro)
NAME:
LMFIT
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 = LMFIT(X, Y, A, FITA=FITA, FUNCTION_NAME = name, ITER=ITER,$
ITMAX=ITMAX, SIGMA=SIGMA, TOL=TOL, YFIT=YFIT, WEIGHTS=WEIGHTS)
CONVERGENCE=CONVERGENCE
INPUTS:
X: A row vector of independent variables. This routine does
not manipulate or use values in X, it simply passes X
to the user-written function.
Y: A row vector containing the dependent variable.
A: A vector that contains the initial estimate for each parameter.
WEIGHTS:
Set this keyword equal to a vector of fitting weights for Y(i). This
vector must be the same length as X and Y. If instrumental (Gaussian)
weighting is desired, that is the measurement errors or standard deviations
of Y are known (SIGMA), then WEIGHTS should be set to 1/(SIGMA^2.)
Or if statistical (Poisson) weighting is appropriate, set WEIGHTS=1/Y.
If WEIGHTS is not specified then, no weighting is assumed (WEIGHTS=1.0).
FITA:
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.
KEYWORDS:
FUNCTION_NAME: The name of the function (actually, a procedure) to
fit. If omitted, "LMFUNCT" is used. The procedure must be written as
described under RESTRICTIONS, below.
ALPHA: The value of the Curvature matrix upon exit.
CHISQ: The value of chi-squared on exit
CONVERGENCE: Returns 1 if the fit converges, 0 if it does
not meet the convergence criteria in ITMAX iterations,
or -1 if a singular matrix is encountered.
COVAR: The value of the Covariance matrix upon exit.
DOUBLE: Set this keyword to force the computations to be performed
in double precision.
ITMAX: Maximum number of iterations. Default = 20.
ITER: The actual number of iterations which were performed
SIGMA: A vector of standard deviations for the returned parameters.
TOL: The convergence tolerance. The routine returns when the
relative decrease in chi-squared is less than TOL in an
interation. Default = 1.e-6.
YFIT: The fitted function evaluated at the input X values.
OUTPUTS:
Returns a vector of calculated parameters, which are improvements
upon the initial guesses supplied in A.
SIDE EFFECTS:
None.
RESTRICTIONS:
The function to be fit must be defined and called LMFUNCT,
unless the FUNCTION_NAME keyword is supplied. This function,
must accept a single value of X (the independent variable), and A
(the fitted function's parameter values), and return an
array whose first (zeroth) element is the evalutated function
value, and next n_elements(A) elements are the partial derivatives
with respect to each parameter in A.
If X is passed in as a double, the returned vector MUST be of
type double as well. Likewise, if X is a float, the returned
vector must also be of type float.
For example, here is the default LMFUNCT in the IDL User's Libaray.
which is called as : out_array = LMFUNCT( X, A )
function lmfunct,x,a
;Return a vector appropriate for LMFIT
;
;The function being fit is of the following form:
; F(x) = A(0) * exp( A(1) * X) + A(2) = bx+A(2)
;
;dF/dA(0) is dF(x)/dA(0) = exp(A(1)*X)
;dF/dA(1) is dF(x)/dA(1) = A(0)*X*exp(A(1)*X) = bx * X
;dF/dA(2) is dF(x)/dA(2) = 1.0
;
;return,[[F(x)],[dF/dA(0)],[dF/dA(1)],[dF/dA(2)]]
;
;Note: returning the required function in this manner
; ensures that if X is double the returned vector
; is also of type double. Other methods, such as
; evaluating size(x) are also valid.
bx=A(0)*exp(A(1)*X)
return,[ [bx+A(2)], [exp(A(1)*X)], [bx*X], [1.0] ]
end
PROCEDURE:
Based upon "MRQMIN", least squares fit to a non-linear
function, pages 683-688, Numerical Recipies in C, 2nd Edition,
Press, Teukolsky, Vettering, and Flannery, 1992.
"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 three consequtive iterations fail
to chang the chi square changes by greater than TOL, or until
ITMAX, but at least ITMIN, 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.
The function may fail to converge, or it can encounter
a singular matrix. If this happens, the routine will fail
with the Numerical Recipes error message:
EXAMPLE:
Fit a function of the form:
f(x)=a(0) * exp(a(1)*x) + a(2) + a(3) * sin(x)
Define a lmfit return function:
function myfunct,x,a
;Return a vector appropriate for LMFIT
;The function being fit is of the following form:
; F(x) = A(0) * exp( A(1) * X) + A(2) + A(3) * sin(x)
; dF(x)/dA(0) = exp(A(1)*X)
; dF(x)/dA(1) = A(0)*X*exp(A(1)*X) = bx * X
; dF(x)/dA(2) = 1.0
; dF(x)/dA(3) = sin(x)
bx=A(0)*exp(A(1)*X)
return,[[bx+A(2)+A(3)*sin(x)],[exp(A(1)*X)],[bx*X],[1.0],[sin(x)]]
end
pro run_lmfunct
x=findgen(40)/20. ;Define indep & dep variables.
y=8.8 * exp( -9.9 * X) + 11.11 + 4.9 * sin(x)
sig=0.05 * y
a=[10.0,-7.0,9.0,4.0] ;Initial guess
fita=[1,1,1,1]
ploterr,x,y,sig
yfit=lmfit(x,y,a,WEIGHTS=(1/sig^2.),FITA=FITA,$
SIGMA=SIGMA,FUNCTION_NAME='myfunct')
oplot,x,yfit
for i=0,3 do print,i,a(i),format='("A (",i1,")= ",F6.2)'
end
MODIFICATION HISTORY:
Written, SVP, RSI, June 1996.
(See /usr/local/rsi/idl_5/lib/lmfit.pro)
NAME:
POLYFITW
PURPOSE:
Perform a least-square polynomial fit with optional error estimates.
CATEGORY:
Curve fitting.
CALLING SEQUENCE:
Result = POLYFITW(X, Y, W, NDegree [, Yfit, Yband, Sigma, A])
INPUTS:
X: The independent variable vector.
Y: The dependent variable vector. This vector should be the same
length as X.
W: The vector of weights. This vector should be same length as
X and Y.
NDegree: The degree of polynomial to fit.
OUTPUTS:
POLYFITW returns a vector of coefficients of length NDegree+1.
OPTIONAL OUTPUT PARAMETERS:
Yfit: The vector of calculated Y's. Has an error of + or - Yband.
Yband: Error estimate for each point = 1 sigma.
Sigma: The standard deviation in Y units.
A: Correlation matrix of the coefficients.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
MODIFICATION HISTORY:
Written by: George Lawrence, LASP, University of Colorado,
December, 1981.
Adapted to VAX IDL by: David Stern, Jan, 1982.
Weights added, April, 1987, G. Lawrence
(See /usr/local/rsi/idl_5/lib/polyfitw.pro)
NAME:
POLY_FIT
PURPOSE:
Perform a least-square polynomial fit with optional error estimates.
This routine uses matrix inversion. A newer version of this routine,
SVDFIT, uses Singular Value Decomposition. The SVD technique is more
flexible, but slower.
Another version of this routine, POLYFITW, performs a weighted
least square fit.
CATEGORY:
Curve fitting.
CALLING SEQUENCE:
Result = POLY_FIT(X, Y, NDegree [,Yfit, Yband, Sigma, CORRM] )
INPUTS:
X: The independent variable vector.
Y: The dependent variable vector, should be same length as x.
NDegree: The degree of the polynomial to fit.
OUTPUTS:
POLY_FIT returns a vector of coefficients with a length of NDegree+1.
OPTIONAL OUTPUT PARAMETERS:
Yfit: The vector of calculated Y's. These values have an error
of + or - Yband.
Yband: Error estimate for each point = 1 sigma
Sigma: The standard deviations of the returned coefficients.
Corrm: Correlation matrix of the coefficients.
Keyword Parameters:
DOUBLE = if set, force computations to be in double precision.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
MODIFICATION HISTORY:
Written by: George Lawrence, LASP, University of Colorado,
December, 1981.
Adapted to VAX IDL by: David Stern, Jan, 1982.
Modified: GGS, RSI, March 1996
Corrected a condition which explicitly converted all
internal variables to single-precision float.
Added support for double-precision inputs.
Added a check for singular array inversion.
SVP, RSI, June 1996
Changed A to Corrm to match IDL5.0 docs.
(See /usr/local/rsi/idl_5/lib/poly_fit.pro)
NAME:
REGRESS
PURPOSE:
Perform a multiple linear regression fit.
REGRESS fits the function:
Y[i] = Const + A[0]*X[0,i] + A[1]*X[1,i] + ... +
A[Nterms-1]*X[Nterms-1,i]
CATEGORY:
G2 - Correlation and regression analysis.
CALLING SEQUENCE:
Result = REGRESS(X, Y, Weights, Yfit, Const, Sigma, Ftest, R, Rmul, Chisq)
INPUTS:
X: The array of independent variable data. X must
be dimensioned as an array of Nterms by Npoints, where
there are Nterms coefficients (independent variables) to be
found and Npoints of samples.
Y: The vector of dependent variable points. Y must have Npoints
elements.
Weights: The vector of weights for each equation. Weights must be a vector
of Npoints elements. For instrumental (Gaussian) weighting,
W[i] = 1/standard_deviation(Y[i])^2. For statistical (Poisson)
weighting, w[i] = 1./Y[i]. For no weighting, set w[i]=1,
and also set the RELATIVE_WEIGHT keyword.
OUTPUTS:
REGRESS returns a column vector of coefficients that has Nterms
elements.
OPTIONAL OUTPUT PARAMETERS:
Yfit: Vector of calculated values of Y with Npoints elements.
Const: Constant term. (A0)
Sigma: Vector of standard deviations for coefficients.
Ftest: The value of F for test of fit.
Rmul: The multiple linear correlation coefficient.
R: Vector of linear correlation coefficients.
Rmul: The multiple linear correlation coefficient.
Chisq: Reduced, weighted chi squared.
Status: A named variable to receive the status of the INVERT
(array inversion) computation. A value of 0 indicates
a successful computation. A value of 1 indicates the
inversion is invalid due to a singular array. A value
of 2 indicates the possibility of an inaccurate result
due to the use of a small pivot element.
KEYWORDS:
RELATIVE_WEIGHT: If this keyword is set, the input weights
(W vector) are assumed to be relative values, and not based
on known uncertainties in the Y vector. Set this keyword in
the case of no weighting, W[*] = 1.
PROCEDURE:
Adapted from the program REGRES, Page 172,
Bevington, Data Reduction and Error Analysis for the
Physical Sciences, 1969.
MODIFICATION HISTORY:
Written, DMS, RSI, September, 1982.
Added RELATIVE_WEIGHT keyword W. Landsman August 1991
Fixed bug in invert Bobby Candey 1991 April 22
Added STATUS argument. GGS, RSI, August 1996
(See /usr/local/rsi/idl_5/lib/regress.pro)
NAME:
SVDFIT
PURPOSE:
Perform a general least squares fit with optional error estimates.
This version uses the Numerical Recipies (2nd Edition) function
SVDFIT. A user-supplied function or a built-in polynomial or
legendre polynomial is fit to the data.
CATEGORY:
Curve fitting.
CALLING SEQUENCE:
Result = SVDFIT(X, Y, [M])
INPUTS:
X: A vector representing the independent variable.
Y: Dependent variable vector. This vector must be same length
as X.
OPTIONAL INPUTS:
M: The number of coefficients in the fitting function. For
polynomials, M is equal to the degree of the polynomial + 1.
If not specified and the keyword A is set, THEN
M = N_ELEMENTS(A).
INPUT KEYWORDS:
A: The inital estimates of the desired coefficients. If M
is specified, THEN A must be a vector of M elements.
If A is specified, THEN the input M can be omitted and
M=N_ELEMENTS(A). If not specified, the initial value
of each coefficient is taken to be 1.0. If both M and A are
specified, them must agree as to the number of paramaters.
DOUBLE: Set this keyword to force double precision computations. This
is helpful in reducing roundoff errors and improves the chances
of function convergence.
WEIGHTS: A vector of weights for Y[i]. This vector must be the same
length as X and Y. If this parameter is ommitted, 1's
(No weighting) are assumed. The error for each term is
weighted by Weight[i] when computing the fit. Gaussian or
instrumental uncertianties should be weighted as
Weight = 1/Sigma where Sigma is the measurement
error or standard deviations of Y. For Poisson or statistical
weighting use Weight=1/Y, since Sigma=sqrt(Y).
FUNCTION_NAME:
A string that contains the name of an optional user-supplied
basis function with M coefficients. If omitted, polynomials
are used.
The function is called: R=SVDFUNCT(X,M)
where X and M are scalar values, and the function value is an
M element vector evaluated at X with the M basis functions.
M is the degree of the polynomial +1 if the basis functions are
polynomials. For example, see the function SVDFUNCT or SVDLEG,
in the IDL User Library:
For more examples, see Numerical Recipes in C, second Edition,
page 676-681.
The basis function for polynomials, is R[j] = x)^j.
The function must be able to return R as a FLOAT vector or
a DOUBLE vector depending on the input type of X.
LEGENDRE: Set this keyword to use the IDL function SVDLEG in the lib
directory to fit the data to an M element legendre polynomial.
This keyword overrides the FUNCTION_NAME keyword.
OUTPUTS:
SVDFIT returns a vector of the M coefficients fitted to the
supplied function.
OPTIONAL OUTPUT PARAMETERS:
CHISQ: Sum of squared errors multiplied by weights if weights
are specified.
COVAR: Covariance matrix of the coefficients.
VARIANCE: Sigma squared in estimate of each coeff(M).
That is sqrt(VARIANCE) equals the 1 sigma deviations
of the returned coefficients.
SIGMA: The 1-sigma error estimates of the returned parameters,
SIGMA=SQRT(VARIANCE).
SINGULAR: The number of singular values returned. This value should
be 0. If not, the basis functions do not accurately
characterize the data.
YFIT: Vector of calculated Y's.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
MODIFICATION HISTORY:
Adapted from SVDFIT, from the book Numerical Recipes, Press,
et. al., Page 518.
minor error corrected April, 1992 (J.Murthy)
Completely rewritten to use the actual Numerical Recipes routines
of the 2nd Edition (V.2.06). Added the DOUBLE, SIGMA, A, and
LEGENDRE keywords. Also changed Weight to Weights to match the
other fitting routines.
EXAMPLE:
(See /usr/local/rsi/idl_5/lib/svdfit.pro)