This page was created by the IDL library routine
mk_html_help. For more information on
this routine, refer to the IDL Online Help Navigator
or type:
? mk_html_help
at the IDL command line prompt.
Last modified: Tue Jul 11 20:43:59 2006.
NAME:
AMOEBA
PURPOSE:
Multidimensional minimization of a function FUNC(X), where
X is an N-dimensional vector, using the downhill simplex
method of Nelder and Mead, 1965, Computer Journal, Vol 7, pp 308-313.
This routine is modified from IDL's amoeba routine, which
in turn is based on the AMOEBA routine, Numerical
Recipes in C: The Art of Scientific Computing (Second Edition), Page
411.
CATEGORY:
Function minimization/maximization. Simplex method.
CALLING SEQUENCE:
Result = AMOEBAFIT(Ftol, data, ....)
INPUTS:
FTOL: the fractional tolerance to be achieved in the function
value. e.g. the fractional decrease in the function value in the
terminating step. This should never be less than the
machine's single or double precision.
DATA: The data values
KEYWORD PARAMETERS:
FUNCTION_NAME: a string containing the name of the function to
be minimized. If omitted, the function FUNC is minimized.
This function must accept an Ndim vector as its only parameter and
return a scalar single or double precision floating point value as its
result.
*************************************************************
This function is now the same form as for curvefit, eg
function_name(coords, param, values, deriv)
*************************************************************
COORD: The dependent variable, [0,1,2,3...] if not passed
WEIGHT: The weights for the data, [1,1,1,1,] if not passed
FUNCTION_VALUE: (output) on exit, an Ndim+1 element vector
containing the function values at the simplex points. The first
element contains the function minimum.
NCALLS: (output) the of times the function was evaluated.
NMAX: the maximum number of function evaluations allowed
before terminating. Default = 5000.
P0: Initial starting point, an Ndim element vector. The starting
point must be specified using either the keyword SIMPLEX, or P0 and
SCALE. P0 may be either single or double precision floating.
For example, in a 3-dimensional problem, if the initial guess
is the point [0,0,0], and it is known that the function's
minimum value occurs in the interval: -10 <
X(0) < 10, -100 < X(1) < 100, -200 < X(2) < 200, specify: P0=[0,0,0],
SCALE=[10, 100, 200].
SCALE: a scalar or Ndim element vector contaiing the problem's
characteristic length scale for each dimension.
SCALE is used with P0 to form an initial (Ndim+1) point simplex.
If all dimensions have the same scale, specify a scalar.
SIMPLEX: (output and/or optional input) On input, if P0 and SCALE
are not set, SIMPLEX contains the Ndim+1 vertices, each of
Ndim elements, of starting simplex, in either single or double
precision floating point, in an (Ndim, Ndim+1) array. On output,
SIMPLEX contains the simplex, of dimensions (Ndim, Ndim+1), enclosing
the function minimum. The first point, Simplex(*,0), corresponds to
the function's minimum.
METHOD: What to minimize. Choices are
'sumsqres' to minimize the sum of (data-model)^2,
'sumabres' to minimize the sum of |data-model|,
'maxabdev' to minimize the maximum of |data-model|
OUTPUTS:
Result: If the minimum is found, an Ndim vector, corresponding to
the Function's minimum value is returned. If a function minimum
within the given tolerance, is NOT found in the given number of
evaluations, a scalar value of -1 is returned.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
PROCEDURE:
This procedure implements the Simplex method, described in
Numerical Recipes, Section 10.4. See also the POWELL procedure.
Advantages: requires only function evaluations, not
derivatives, may be more reliable than the POWELL method.
Disadvantages: not as efficient as Powell's method, and usually
requires more function evaluations.
Results are performed in the mode (single or double precision)
returned by the user-supplied function. The mode of the inputs P0,
SCALE, or SIMPLEX, should match that returned by the function. The
mode of the input vector supplied to the user-written function, is
determined by P0, SCALE, or SIMPLEX.
EXAMPLE:
Use Amoeba to find the slope and intercept of a straight line fitting
a given set of points minimizing the maximum error:
The function to be minimized returns the maximum error,
given p(0) = intercept, and p(1) = slope:
FUNCTION FUNC, x, p, y
RETURN, (p[0] + p[1] * x)
END
Define the data points:
x = findgen(17)*5
y = [ 12.0, 24.3, 39.6, 51.0, 66.5, 78.4, 92.7, 107.8, 120.0, $
135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]
Call the function. Fractional tolerance = 1 part in 10^5,
Initial guess = [0,0], and the minimum should be found within
a distance of 100 of that point:
r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_VALUE=fval, METHOD='maxabres')
Check for convergence:
if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
Print results.
print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
Intercept, Slope: 11.4100 2.72800
Function value: 1.33000
MODIFICATION HISTORY:
LAY, Oct, 2002. Modified Amoeba.
(See math/amoebafit.pro)
NAME:
cubic
PURPOSE: (one line)
evaluate a cubic given values and derivatives at two points
DESCRIPTION:
evaluate a cubic given values and derivatives at two points
CATEGORY:
Math
CALLING SEQUENCE:
y = cubic(x0, x1, y0, y1, dy0, dy1, x, dy, dx0, dx1, a)
INPUTS:
x0, x1 - the fixed points
y0, y1 - values at the fixed points
dy0, dy1 - derivatives at the fixed points, same
units as (y1-y0)/(x1-x0)
x - new point at which to evaluate the cubic
OPTIONAL INPUT PARAMETERS:
none
KEYWORD INPUT PARAMETERS:
none
KEYWORD OUTPUT PARAMETERS:
none
OUTPUTS:
returns y = value of the cubic at x
OPTIONAL OUTPUT PARAMETERS:
dy = derivative of the cubic at x
dx0 = (x-x0)/(x1-x0) and dx1 = (x-x1)/(x1-x0) - for debug only
a = coefficients for y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
(for debug only)
COMMON BLOCKS:
None
SIDE EFFECTS:
RESTRICTIONS:
None
PROCEDURE:
Given a cubic function with the following endpoints
y(x0) = y0
y(x1) = y1
dy(x0)/dx = dy0
dy(x1)/dx = dy0
calculate the coefficients so that
y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
dx0 = (x-x0)/(x1-x0)
dx1 = (x-x1)/(x1-x0)
MODIFICATION HISTORY:
Written 2000 October, by Leslie Young, SwRI
Nov 2000. Allow non-scalar argument, LAY
11 Mar 2006. LAY. Moved to $idl/layoung/math, check nparams
before calculating dy
(See math/cubic.pro)
NAME:
CURVEFITLAY
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:
Sigma: 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 math/curvefitlay.pro)
NAME:
CURVEFITLAYTEST
PURPOSE:
tester/driver for curvefitlay
CATEGORY:
E2 - Curve and Surface Fitting.
CALLING SEQUENCE:
CURVEFITLAYTEST
INPUTS:
KEYWORDS:
OUTPUTS:
OPTIONAL OUTPUT PARAMETERS:
COMMON BLOCKS:
SIDE EFFECTS:
RESTRICTIONS:
PROCEDURE:
MODIFICATION HISTORY:
Leslie Young 2004 Apr 30.
(See math/curvefitlaytest.pro)
NAME:
dbrent
PURPOSE:
minimum finding with derivatives using Brent's method
DESCRIPTION:
Given a function f and its derivative function df, and given a
bracketing triplet of abscissas ax, bx, cx [such that bx is between ax
and cx, and f(bx) is less than both f(ax) and f(cx)], this routine
isolates the minimum to a fractional precision of about tol using a modi
cation of Brent's method that uses derivatives. The abscissa of the
minimum is returned as xmin, and the minimum function value is returned
as dbrent, the returned function value.
CALLING SEQUENCE:
min = dbrent(func, ax,bx,cx, xmin, tol=tol, extras=extras)
INPUT PARAMETERS
func: a string containing the name of the function to be minimized.
func, x, f, df, OR func,x,f,df, p
ax, bx, cx: a bracketing triplet of abscissas ax, bx, cx
[such that bx is between ax and cx, and f(bx) is less than both f(ax) and f(cx)]
OPTIONAL INPUT PARAMETERS
tol - tolerance
funcp - optional other parameters to be passed to f
OUTPUT PARAMETERS
xmin abscissa of the minimum
dbrent - minimum function value, returned
PROCEDURE:
Based on dbrent in Numerical Recipes
REVISION HISTORY
(See math/dbrent.pro)
NAME: expm1 PURPOSE: (one line) Return exp(x)-1, even for small x. DESCRIPTION: Return exp(x)-1, even for small x. CATEGORY: Math CALLING SEQUENCE: y = expm1(x) INPUTS: x - exponent OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: exp(x)-1 COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2000 October, by Leslie Young, SwRI Nov 2000. Allow non-scalar argument, LAY
(See math/expm1.pro)
NAME: fMaxLoc PURPOSE: (one line) Returns the location of the maximum in a 1-D array DESCRIPTION: This function returns finds the maximum point in a 1-D array, then fits a parabola to the maximum and its left and right neighbors. The center of the parabola is assumed to be the fractional location of the maximum. NOTICE that this function is useless if the FWHM is small compared to a pixel. CATEGORY: Spectral extraction CALLING SEQUENCE: fMax = fMaxLoc(v) INPUTS: myVec - Input 1-D array. OPTIONAL INPUT PARAMETERS: KEYWORD PARAMETERS: OUTPUTS: fMax - the location of the maximum COMMON BLOCKS: SIDE EFFECTS: RESTRICTIONS: PROCEDURE: MODIFICATION HISTORY: Written June 6, 1996, Eliot Young, NASA Ames Research Center Modified June 22, 1999, Jason Cook & Leslie Young, BU. Found y1=0, y3=0, and arrWidth bugs. Modified Dec 21, 1999, Leslie Young, SwRI. Found bracketing bug
(See math/fmaxloc.pro)
NAME: gcf PURPOSE: (one line) Returns the incomplete gamma function Q(a; x) DESCRIPTION: Returns the incomplete gamma function Q(a; x) evaluated by its continued fraction representation as gammcf. Also returns ln Gamma CATEGORY: Math CALLING SEQUENCE: y = gcf(a,x, gln=gln) INPUTS: a - first parameter for Q x - second parameter for Q OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: gln = ln(Gamma(a)) OUTPUTS: Q(a; x) COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2003 November, by Leslie Young, SwRI based on Numerical Recipes gcf, Section 6.2
(See math/gcf.pro)
NAME:
madline
PURPOSE: (one line)
Fits y = a + bx by the criterion of minimum absolute deviations (MAD).
DESCRIPTION:
Fits y = a + bx by the criterion of least absolute deviations. The arrays x[0..ndata-1] and
y[0..ndata-1] are the input experimental points. The fitted parameters a and b are output,
along with abdev, which is the mean absolute deviation (in y) of the experimental points from
the fitted line.
CATEGORY:
Statistics
CALLING SEQUENCE:
[a,b] = madline(x, y, yfit, abdev)
INPUTS:
x - An n-element vector of independent variables.
y - A vector of dependent variables, the same length as X.
OPTIONAL INPUT PARAMETERS:
INPUT KEYWORD PARAMETERS:
OUTPUT KEYWORD PARAMETERS:
OUTPUTS:
returns [a,b], the fitted parameters for y = a + b * x
Yfit - A named variable that will contain the vector of calculated Y values.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
RESTRICTIONS:
None.
PROCEDURE:
MODIFICATION HISTORY:
Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
based on Numerical Recipes medline.
(See math/madline.pro)
NAME:
madlinew
PURPOSE: (one line)
Fits y = a + bx by the criterion of minimum absolute deviations (MAD).
DESCRIPTION:
Fits y = a + bx by the criterion of least absolute deviations. The arrays x[0..ndata-1] and
y[0..ndata-1] are the input experimental points. The fitted parameters a and b are output,
which minimize the sum of (w * |y - a - bx| ).
CATEGORY:
Statistics
CALLING SEQUENCE:
[a,b] = madlinew(x, y, w, yfit)
INPUTS:
x - An n-element vector of independent variables.
y - A vector of dependent variables, the same length as X.
w - weights
OPTIONAL INPUT PARAMETERS:
INPUT KEYWORD PARAMETERS:
OUTPUT KEYWORD PARAMETERS:
OUTPUTS:
returns [a,b], the fitted parameters for y = a + b * x
Yfit - A named variable that will contain the vector of calculated Y values.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
RESTRICTIONS:
None.
PROCEDURE:
MODIFICATION HISTORY:
Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
based on Numerical Recipes medfit.
(See math/madlinew.pro)
NAME: madscale PURPOSE: (one line) minimize the average deviation (MAD) for y = a * x. DESCRIPTION: This routine estimates the scaling factor, a, such that the average deviation, total(abs(y-a*x) ), is minimized. The statistics are robust in that the result is insensitive to outliers. Note that for y[i] = a * x[i], x can be thought of as either a coordinate, or a template (e.g., a line profile). CATEGORY: Statistics CALLING SEQUENCE: a = madscale(x, y) INPUTS: x - dependent variable. y - data, proportional to x. OPTIONAL INPUT PARAMETERS: INPUT KEYWORD PARAMETERS: OUTPUT KEYWORD PARAMETERS: OUTPUTS: COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: None. PROCEDURE: The goal is to minimize f = total(abs(y-a*x)) = total( abs(x) * abs(y/x - a) ) This is equivalent to finding the weighted median of y/x, with abs(x) as the weights. MODIFICATION HISTORY: Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.
(See math/madscale.pro)
NAME: madscalew PURPOSE: (one line) minimize the average deviation (MAD) for y = a * x. DESCRIPTION: This routine estimates the scaling factor, a, such that the average deviation, total(w * abs(y-a*x) ), is minimized. The statistics are robust in that the result is insensitive to outliers. Note that for y[i] = a * x[i], x can be thought of as either a coordinate, or a template (e.g., a line profile). CATEGORY: Statistics CALLING SEQUENCE: a = madscale(x, y) INPUTS: x - dependent variable. y - data, proportional to x. OPTIONAL INPUT PARAMETERS: w - weights for y INPUT KEYWORD PARAMETERS: OUTPUT KEYWORD PARAMETERS: OUTPUTS: COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: None. PROCEDURE: The goal is to minimize f = total(w * abs(y-a*x)) = total( w * abs(x) * abs(y/x - a) ) This is equivalent to finding the weighted median of y/x, with abs(x) as the weights. MODIFICATION HISTORY: Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.
(See math/madscalew.pro)
NAME:
medianw
PURPOSE: (one line)
minimize the weighted average deviation.
DESCRIPTION:
This routine estimates the average data value, xmed, such that
the average deviation, total(weight * abs(x-xmed) ), is minimized.
For equally-weighted points, this is the median.
The statistics are robust in that the result is insensitive
to outliers.
CATEGORY:
Statistics
CALLING SEQUENCE:
mean = medianw(x, w)
INPUTS:
x - Input data to be analyzed.
OPTIONAL INPUT PARAMETERS:
w - weights. Assumed equally weighted if not passed.
INPUT KEYWORD PARAMETERS:
OUTPUT KEYWORD PARAMETERS:
OUTPUTS:
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
RESTRICTIONS:
None.
PROCEDURE:
Minimizing f(xmed) = total(weight * abs(x-xmed)) is equivalent to finding
df(xmed)/d xmed = 0 = -total(weight * sgn(x-xmed)).
If x is sorted, and xmed=x[imed], this can be written as
0 = total(weight[0:imed-1] * sgn(x[0:imed-1]-xmed))
+ total(weight[imed+1:n-1] * sgn(x[imed+1:n-1]-xmed))
or
0 = - total(weight[0:imed-1]) + total(weight[imed+1:n-1])
or
total(weight[0:imed-1]) = total(weight[imed+1:n-1])
MODIFICATION HISTORY:
Written by Leslie A. Young, Soutwest Research Institute, 2002 Jul 23.
Modified LAY Oct 2002. Better medianwTEST;
Avoid eq for comparing floats in test for one median value
(equivalent to an odd number of equally weighted values)
Average two neighbors if xmed is not one of the listed x's
(equivalent to an even number of equally weighted values)
Modified LAY Oct 2002. Changed name to medianw;
(See math/medianw.pro)
NAME: normv PURPOSE: (one line) Return a normalized vector DESCRIPTION: Return a normalized vector CATEGORY: Math CALLING SEQUENCE: n = normv(v) INPUTS: v - a vector of arbitrary length OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: v/vabs(v) COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2000, by Leslie Young, SwRI
(See math/normv.pro)
NAME:
PIECECUB_INTERP
PURPOSE:
interpolate 1-d data and derivatives with piece-wise cubic functions
EXPLANATION:
Interpolation into a function with continuous values and
derivatives, with tabulated values. Between grid points, the
function is defined by a cubic polynomial. Unlike cubic
splines, the second derivative is not necessarily continuous.
Calculations are done in double precision.
CALLING SEQUENCE:
piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint
INPUT PARAMETERS:
Xtab - Vector containing the current independent variable grid.
Must be monotonic increasing or decreasing
Ytab - Vector containing the current dependent variable values at
the XTAB grid points.
dYtab - Vector containing the current derivative of the dependent
variable values at the XTAB grid points.
Xint - Scalar or vector containing the new independent variable grid
points for which interpolated value(s) of the dependent
variable are sought. Note that -- due to a limitation of the
intrinsic INTERPOLATE() function -- Xint is always converted to
floating point internally.
OUTPUT PARAMETERS:
Yint - Scalar or vector with the interpolated value(s) of the
dependent variable at the XINT grid points.
YINT is double precision if Ytab is, float otherwise.
dYint - Scalar or vector with the interpolated value(s) of the
derivative at the XINT grid points
DYINT is double precision if DYtab is, float otherwise.
OPTIONAL INPUT KEYWORD:
EXAMPLE:
PROCEDURES CALLED:
TABINV, ZPARCHECK
PROCEEDURE
________ DEFINITIONS ___________________
Start with some definitions:
x the independent variable
x0 x at node 0 (or node n)
x1 x at node 1 (or node n+1)
y the dependent variable
y0 y at node 0
dy0 dy/dx at node 0
y1 y at node 0
dy1 dy/dx at node 0
Let's define some things to make our lives easier:
d0 = x-x0
d1 = x-x1
del = x1-x0
so d1 and d0 are functions of x.
d0(x0) = 0 d0(x1) = del
d1(x0) = -del d1(x1) = 0
Using d0, d1, and del, define four functions:
P0(x) = (3*d0-d1) * d1^2 / del^3
P1(x) = -(3*d1-d0) * d0^2 / del^3
Q0(x) = d0 * d1^2 / del^2
Q1(x) = d1 * d0^2 / del^2
with the following derivatives
dP0(x)/dx = 6*d0*d1/del^3
dP1(x)/dx = -6*d0*d1/del^3
dQ0(x)/dx = d1*(d1+2*d0)/del^2
dQ1(x)/dx = d0*(d0+2*d1)/del^2
These functions have very tidy values and derivatives at the nodes:
P0(x0) = 1 P0(x1) = 0 dP0(x0)/dx = 0 P0(x1)/dx = 0
P1(x0) = 0 P1(x1) = 1 dP1(x0)/dx = 0 P1(x1)/dx = 0
Q0(x0) = 0 Q0(x1) = 0 dQ0(x0)/dx = 1 Q0(x1)/dx = 0
Q1(x0) = 0 Q1(x1) = 0 dQ1(x0)/dx = 0 Q1(x1)/dx = 1
________ CALCULATING THE CUBIC FUNCTION ___________________
This makes it trivial to define the cubic function that
matches the nodes:
y(x) = y0 * P0(x) + y1 * P1(x) + dy0 * Q0(x) + dy1 * Q1(x) [CUBIC]
MODIFICATION HISTORY:
(See math/piececub_interp.pro)
NAME:
robofit
PURPOSE: (one line)
Least-squares fitting with the identification of outliers.
DESCRIPTION:
Robust 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.
This routine differs from curvefit in that it assumes that
the distribution may be non-normal -- having an "outerlier"
tail. The function is first fit to the data minimizing
the sum of absolute deviations. From this solution,
outliers are identified, and expunged from the final
least-squares fit.
To maintain compatibility with CURVEFIT, the calling
sequence is the same (with the addition of THRESH).
CATEGORY:
E2 - Curve and Surface Fitting.
CALLING SEQUENCE:
Result = ROBOFIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
INPUTS:
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:
Sigma: A vector of standard deviations for the parameters in A.
COMMON BLOCKS:
None.
SIDE EFFECTS:
None.
RESTRICTIONS:
The function to be fit is defined identically to that used by CURVEFIT.
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:
Call AMOEBA
MODIFICATION HISTORY:
Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
based on Numerical Recipes medfit.
(See math/robofit.pro)
NAME: roboline PURPOSE: (one line) Robust fitting of a line DESCRIPTION: Robust fitting of y = a = b * x CATEGORY: Statistics CALLING SEQUENCE: [a,b] = roboline(x, y, sigma, thresh, bplist, GDINDEX = gdindex, coeferr = coefErr) INPUTS: x - the indices y - the data array, same length as x sigma - estimated noise for each pixel in y thresh - threshold for outliers OPTIONAL INPUT PARAMETERS: bplist - optional bad pixel list (-1 for all OK) INPUT KEYWORD PARAMETERS: OUTPUT KEYWORD PARAMETERS: gdindex - named variable that will contain an array of indices used OUTPUTS: returns scaleFac such that y ~= a + b * x COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: None. PROCEDURE: First iteration finds minimum absolute deviation. Outliers are identified after the first iteration. The second iteration minimizes least squares MODIFICATION HISTORY: Written July 22, 2002 Eliot Young, SwRI Rewritten by Leslie A. Young, SwRI, 2002 Oct 29. - first iteration minimizes absolute deviation - exchanged order of x and y in argument list - only two iterations - optionally returns the error and the list of good pixel indices
(See math/roboline.pro)
NAME:
robomarg
PURPOSE: (one line)
Return the robust averages of each row or each column.
DESCRIPTION:
This routine calls robomean on a row-by-row or col-by-col basis
to get the robust mean of each row or each column. The /GetRowMeans
keyword returns a column containing the robust average of all of
the rows. The /GetColMeans keyword returns a row with the robust
averages of all the columns. /GetRowMeans is the default.
CATEGORY:
spectral extraction
CALLING SEQUENCE:
means = robomarg(array)
INPUTS:
array - Input array to be scanned.
OPTIONAL INPUT PARAMETERS:
KEYWORD INPUT PARAMETERS:
GetRowMeans - FLAG, return a column which contains the
robust means of each row.
GetColMeans - FLAG, return a row which contains the
robust means of each column.
bad - bad pixel mask (1 = bad, 0 = good)
OUTPUTS:
means - The 1-D column or row with averages
KEYWORD OUTPUT PARAMETERS:
COMMON BLOCKS:
SIDE EFFECTS:
RESTRICTIONS:
PROCEDURE:
MODIFICATION HISTORY:
96/06/06, EFY, NASA Ames Research Center
96/10/17, EFY, changed to a function that returns means
01/11/02, LAY, add bad pixel mask option
(See math/robomarg.pro)
NAME: roboscale PURPOSE: (one line) Robust scaling of a template to match an array. DESCRIPTION: Robust fitting of y = scaleFac * x CATEGORY: Statistics CALLING SEQUENCE: scaleFac = roboscale(x, y, sigma, thresh, bplist, ) INPUTS: x - the template y - the data array, same length as x sigma - estimated noise for each pixel in y thesh - threshold for outliers bplist - optional bad pixel list (-1 for all OK) OPTIONAL INPUT PARAMETERS: INPUT KEYWORD PARAMETERS: OUTPUT KEYWORD PARAMETERS: gdindex - named variable that will contain an array of indices used OUTPUTS: returns scaleFac such that y ~= scaleFac * x COMMON BLOCKS: None. SIDE EFFECTS: None. RESTRICTIONS: None. PROCEDURE: First iteration finds minimum absolute deviation. Outliers are identified after the first iteration. The second iteration minimizes least squares MODIFICATION HISTORY: Written July 22, 2002 Eliot Young, SwRI Rewritten by Leslie A. Young, SwRI, 2002 Oct 29. - first iteration minimizes absolute deviation - exchanged order of x and y in argument list - only two iterations - optionally returns the error and the list of good pixel indices
(See math/roboscale.pro)
NAME: sgn PURPOSE: (one line) return the sign of a number or list DESCRIPTION: CATEGORY: Math CALLING SEQUENCE: y = sgn(x) INPUTS: x - argument OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: exp(x)-1 COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2002 August, by Leslie Young, SwRI
(See math/sgn.pro)
NAME: vabs PURPOSE: (one line) Return the length of a vector DESCRIPTION: Return the length of a vector CATEGORY: Math CALLING SEQUENCE: len = abs(v) INPUTS: v - a vector of arbitrary length OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: sqrt(total(v^2)) COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2000, by Leslie Young, SwRI 2006 Jan 12 LAY change to double
(See math/vabs.pro)
NAME: vang PURPOSE: (one line) Return the angle between two vectors DESCRIPTION: Return the angle (in radians) between two vectors CATEGORY: Math CALLING SEQUENCE: d = vang(v1,v2) INPUTS: v1, v2 - two vector of arbitrary length OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: angle between v1 and v2 COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2000, by Leslie Young, SwRI
(See math/vang.pro)
NAME: vdot PURPOSE: (one line) Return the dot product of two vectors DESCRIPTION: Return the dot product of two vectors CATEGORY: Math CALLING SEQUENCE: d = vdot(v1,v2) INPUTS: v1, v2 - two vector of arbitrary length OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: none KEYWORD OUTPUT PARAMETERS: none OUTPUTS: v1 . v2 COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Written 2000, by Leslie Young, SwRI
(See math/vdot.pro)
NAME: wmean PURPOSE: (one line) weighted mean. DESCRIPTION: Return weighted mean and the error of the mean. CATEGORY: Math CALLING SEQUENCE: meanx = wmean(x, xerr, meanerr) INPUTS: x - array of values xerr - array of errors for x OPTIONAL INPUT PARAMETERS: none KEYWORD INPUT PARAMETERS: exterr - calculate the "external error" using scatter KEYWORD OUTPUT PARAMETERS: none OUTPUTS: return weighted mean. manxerr - error in the mean COMMON BLOCKS: None SIDE EFFECTS: RESTRICTIONS: None PROCEDURE: MODIFICATION HISTORY: Added to layoung/math 2004 Feb, by Leslie Young, SwRI fixed typo 2005 Aug 4 LAY
(See math/wmean.pro)