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)