NAME:
fitexy
PURPOSE:
Linear Least-squares approximation in one-dimension (y = a + b*x),
when both x and y data have errors (e.g. Gaussian noise).
CALLING EXAMPLE:
fitexy, x, y, A, B, X_SIG=sigx, Y_SIG=sigy, [sigma_A_B, chi_sq, q ]
INPUTS:
x = array of values for independent variable.
y = array of data values assumed to be linearly dependent on x.
REQUIRED INPUT KEYWORDS:
X_SIGMA = scalar or array specifying the standard deviation of x data.
Y_SIGMA = scalar or array specifying the standard deviation of y data.
OPTIONAL INPUT KEYWORD:
TOLERANCE = desired accuracy of minimum & zero location, default=1.e-3.
OUTPUTS:
A_intercept = constant parameter result of linear fit,
B_slope = slope parameter, so that:
( A_intercept + B_slope * x ) approximates the y data.
OPTIONAL OUTPUT:
sigma_A_B = two element array giving standard deviation of
A_intercept and B_slope parameters, respectively.
The standard deviations are not meaningful if (i) the
fit is poor (see parameter q), or (ii) b is so large that
the data are consistent with a vertical (infinite b) line.
If the data are consistent with *all* values of b, then
sigma_A_B = [1e33,e33]
chi_sq = resulting minimum Chi-Square of Linear fit, scalar
q - chi-sq probability, scalar (0-1) giving the probability that
a correct model would give a value equal or larger than the
observed chi squared. A small value of q indicates a poor
fit, perhaps because the errors are underestimated.
COMMON:
common fitexy, communicates the data for computation of chi-square.
CALLS:
function chisq_fitexy
pro minf_bracket
pro minf_parabolic
function zbrent
function chi_sqr1 (from Statistics library)
PROCEDURE:
From "Numerical Recipes" column by Press and Teukolsky:
in "Computer in Physics", May, 1992 Vol.6 No.3
MODIFICATION HISTORY:
Written, Frank Varosi NASA/GSFC September 1992.
Now returns q rather than 1-q W. Landsman December 1992
(See /usr/local/idl/lib/zastron/math/fitexy.pro)
NAME GAMMLN PURPOSE: Return the natural log of the gamma function. Adapted from the algorithm GAMMLN in Section 6.1 in "Numerical Recipes " (Second Edition) by Press et al. 1992. This function became obsolete in IDL 2.4.0 when the equivalent LNGAMMA intrinsic function was introduced. CALLING SEQUENCE: result = gammln ( xx ) INPUTS: xx - numeric scalar or vector for which the log of the gamma function will be evaluated. Must be > 0 OUTPUT: result = alog ( gamma(xx) ). The result will double precision if xx is double, otherwise the result is floating point. NOTES: IDL has an intrinsic gamma function, GAMMA, but overflow occurs in GAMMA for X > 34.5. By computing the log of the gamma function, one can deal with much larger input values. GAMMLN also allows double precision computation, not available with GAMMA. EXAMPLE: Compare the output of GAMMA with GAMMLN IDL> x = findgen(15)+0.5 IDL> print, alog(gamma(x)) IDL> print,gammln(x) METHOD: Uses the expansion of Lanczos as described in Press et al. (1986) REVISION HISTORY: Written, W. Landsman July, 1992 Double Precision update October, 1992
(See /usr/local/idl/lib/zastron/math/gammln.pro)
NAME: GAUSSIAN PURPOSE: Compute the 1-D Gaussian function and optionally the derivative at an array of points. CALLING: y = gaussian( xi, parms,[ pderiv ]) INPUTS: xi = array, independent variable of Gaussian function. parms = parameters of Gaussian, 2 or 3 element array: parms(0) = maximum value (factor) of Gaussian, parms(1) = mean value (center) of Gaussian, parms(2) = standard deviation (sigma) of Gaussian. (if parms has only 2 elements then sigma taken from common). OPTIONAL OUTPUT: pderiv = optional output of partial derivatives, computed only if parameter is present in call. pderiv(*,i) = partial derivative at all xi absisca values with respect to parms(i), i=0,1,2. Function returns array of Gaussian evaluated at xi. EXAMPLE: Evaulate a Gaussian centered at x=0, with sigma=1, and a peak value of 10 at the points 0.5 and 1.5. Also compute the derivative IDL> f = gaussian( [0.5,1.5], [10,0,1], DERIV ) ==> f= [8.825,3.25]. DERIV will be a 2 x 3 array containing the numerical derivative at the two point with respect to the 3 parameters. COMMON BLOCKS: common gaussian, sigma HISTORY: Written, Frank Varosi NASA/GSFC 1992.
(See /usr/local/idl/lib/zastron/math/gaussian.pro)
NAME:
KSONE
PURPOSE:
Return the Kolmogorov-Smirnov statistic and associated probability for
for an array of data values and a user-supplied cumulative distribution
function (CDF) of a single variable. Algorithm from the procedure of
the same name in "Numerical Recipes" by Press et al. 2nd edition (1992)
CALLING SEQUENCE:
ksone, data, func_name, D, prob, [ /PLOT ]
INPUT PARAMATERS:
data - vector of data values, must contain at least 4 elements for the
K-S statistic to be meaningful
func_name - scalar string giving the name of the cumulative distribution
function
OUTPUT PARAMETERS:
D - floating scalar giving the Kolmogorov-Smirnov statistic. It
specified the maximum deviation between the cumulative
distribution of the data and the supplied function
prob - floating scalar between 0 and 1 giving the significance level of
the K-S statistic. Small values of PROB show that the
cumulative distribution function of DATA is significantly
different from FUNC_NAME.
OPTIONAL INPUT KEYWORD:
PLOT - If this keyword is set and non-zero, then KSONE will display a
plot of the CDF of the data with the supplied function
superposed. The data value where the K-S statistic is
computed (i.e. at the maximum difference between the data CDF
and the function) is indicated by a vertical line.
EXAMPLE:
Determine if a vector created by the RANDOMN function is really
consistent with a gaussian distribution.
The CDF of a gaussian is the error function except that a factor
of 2 is included in the error function. So we must create a special
function:
function gauss_cdf
return, errorf( x/sqrt(2) )
end
IDL> data = randomn(seed, 50) ;create data array to be tested
IDL> ksone, abs(data), 'gauss_cdf', D, prob, /PLOT ;Use K-S test
PROB gives the probability that DATA came from a gaussian distribution
NOTES:
Note that the 2nd (1992) edition of Numerical Recipes includes
a more accurate computation of the K-S significance for small
values of N.
PROCEDURE CALLS
procedure prob_ks - computes significance of K-S distribution
REVISION HISTORY:
Written W. Landsman August, 1992
(See /usr/local/idl/lib/zastron/math/ksone.pro)
NAME: KSTWO PURPOSE: Return the Kolmogorov-Smirnov statistic and associated probability that two arrays of data values are drawn from the same distribution Algorithm taken from procedure of the same name in "Numerical Recipes" by Press et al., 2nd edition (1992) CALLING SEQUENCE: kstwo, data1, data2, D, prob INPUT PARAMATERS: data1 - vector of data values, at least 4 data values must be included for the K-S statistic to be meaningful data2 - second set of data values, does not need have the same number of elements as data1 OUTPUT PARAMETERS: D - floating scalar giving the Kolmogorov-Smirnov statistic. It specifies the maximum deviation between the cumulative distribution of the data and the supplied function prob - floating scalar between 0 and 1 giving the significance level of the K-S statistic. Small values of PROB show that the cumulative distribution function of DATA1 is significantly different from DATA2 EXAMPLE: Test whether two vectors created by the RANDOMN function likely came from the same distribution IDL> data1 = randomn(seed,40) ;Create data vectors to be IDL> data2 = randomn(seed,70) ;compared IDL> kstwo, data1, data2, D, prob & print,D,prob PROCEDURE CALLS procedure prob_ks - computes significance of K-S distribution REVISION HISTORY: Written W. Landsman August, 1992
(See /usr/local/idl/lib/zastron/math/kstwo.pro)
NAME: LINTERP PURPOSE: Linearly interpolate tabulated data from one data grid to another. Use QUADTERP for quadratic interpolation. CALLING SEQUENCE: LINTERP, Xtab, Ytab, Xint, Yint, [MISSING = ] 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. Xint - Scalar or vector containing the new independent variable grid points for which interpolated value(s) of the dependent variable are sought. 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 XTAB or YTAB are double, otherwise YINT is REAL*4 OPTIONAL INPUT KEYWORD: MISSING - Scalar specifying YINT value(s) to be assigned, when Xint value(s) are outside of the range of Xtab. Default is to truncate the out of range YINT value(s) to the nearest value of YTAB. See the help for the INTERPOLATE function. EXAMPLE: To linearly interpolate from an IUE spectrum wavelength grid WAVE, FLUX to another grid defined as: WGRID = [1540., 1541., 1542., 1543., 1544, 1545.] IDL> LINTERP, WAVE, FLUX, WGRID, FGRID FGRID will be a 6 element vector containing the values of FLUX linearly interpolated onto the WGRID wavelength scale PROCEDURE: Uses TABINV to calculate the effective index of the values in Xint in the table Xtab. The resulting index is used with the intrinsic INTERPOLATE function to find the corresponding Yint value in Ytab. Unless the MISSING keyword is supplied, out of range Yint values are truncated to the nearest value of Ytab. NOTES: Users with IDL versions before V2.2.2 need to replace the call to INTERPOLATE with the commented lines at the end of the procedure PROCEDURES CALLED: TABINV, ZPARCHECK MODIFICATION HISTORY: Adapted from the IUE RDAF, W. Landsman October, 1988 Modified to use the new INTERPOLATE function June, 1992 Modified to always return REAL*4 October, 1992 Added MISSING keyword August, 1993
(See /usr/local/idl/lib/zastron/math/linterp.pro)
NAME: minF_bracket PURPOSE: Bracket a local minimum of a 1-D function with 3 points, thus ensuring that a minimum exists somewhere in the interval. This routine assumes that the function has a minimum somewhere.... Routine can also be applied to a scalar function of many variables, for such case the local minimum in a specified direction is bracketed, This routine is called by minF_conj_grad, to bracket minimum in the direction of the conjugate gradient of function of many variables CALLING EXAMPLE: xa=0 & xb=1 minF_bracket, xa,xb,xc, fa,fb,fc, FUNC_NAME="name" ;for 1-D func. or: minF_bracket, xa,xb,xc, fa,fb,fc, FUNC="name", $ POINT=[0,1,1], $ DIRECTION=[2,1,1] ;for 3-D func. INPUTS: xa = scalar, guess for point bracketing location of minimum. xb = scalar, second guess for point bracketing location of minimum. KEYWORDS: FUNC_NAME = function name (string) Calling mechanism should be: F = func_name( px ) where: px = scalar or vector of independent variables, input. F = scalar value of function at px. POINT_NDIM = when working with function of N variables, use this keyword to specify the starting point in N-dim space. Default = 0, which assumes function is 1-D. DIRECTION = when working with function of N variables, use this keyword to specify the direction in N-dim space along which to bracket the local minimum, (default=1 for 1-D). (xa,xb,xc) are then relative distances from POINT_NDIM. OUTPUTS: xa,xb,xc = scalars, 3 points which bracket location of minimum, that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists. When working with function of N variables (xa,xb,xc) are then relative distances from POINT_NDIM, in the direction specified by keyword DIRECTION, with scale factor given by magnitude of DIRECTION. OPTIONAL OUTPUT: fa,fb,fc = value of function at 3 points which bracket the minimum, again note that fb < fa and fb < fc if minimum exists. PROCEDURE: algorithm from Numerical Recipes (by Press, et al.), sec.10.1 (p.281). MODIFICATION HISTORY: Written, Frank Varosi NASA/GSFC 1992.
(See /usr/local/idl/lib/zastron/math/minf_bracket.pro)
NAME:
minF_conj_grad
PURPOSE:
Find the local minimum of a scalar function of several variables using
the Conjugate Gradient method (Fletcher-Reeves-Polak-Ribiere algorithm).
Function may be anything with computable partial derivatives.
Each call to minF_conj_grad performs one iteration of algorithm,
and returns an N-dim point closer to the local minimum of function.
CALLING EXAMPLE:
p_min = replicate( 1, N_dim )
minF_conj_grad, p_min, f_min, conv_factor, FUNC_NAME="name",/INITIALIZE
while (conv_factor GT 0) do begin
minF_conj_grad, p_min, f_min, conv_factor, FUNC_NAME="name"
endwhile
INPUTS:
p_min = vector of independent variables, location of minimum point
obtained from previous call to minF_conj_grad, (or first guess).
KEYWORDS:
FUNC_NAME = function name (string)
Calling mechanism should be: F = func_name( px, gradient )
where:
F = scalar value of function at px.
px = vector of independent variables, input.
gradient = vector of partial derivatives of the function
with respect to independent variables, evaluated at px.
This is an optional output parameter:
gradient should not be calculated if parameter is not
supplied in call (Unless you want to waste some time).
/INIT must be specified on first call (whenever p_min is a guess),
to initialize the iteration scheme of algorithm.
/USE_DERIV causes the directional derivative of function to be used
in the 1-D minimization part of algorithm
(default is not to use directional derivative).
TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).
/QUADRATIC runs simpler version which works only for quadratic function.
OUTPUTS:
p_min = vector giving improved solution for location of minimum point.
f_min = value of function at p_min.
conv_factor = gives the current rate of convergence (change in value),
iteration should be stopped when rate gets near zero.
EXTERNAL CALLS:
pro minF_bracket, to find 3 points which bracket the minimum in 1-D.
pro minF_parabolic, to find minimum point in 1-D.
pro minF_parabol_D, to find minimum point in 1-D, using derivatives.
COMMON BLOCKS:
common minf_conj_grad, grad_conj, grad_save, gs_norm
(to keep conjugate gradient, gradient and norm from previous iteration)
PROCEDURE:
Algorithm adapted from Numerical Recipes, sec.10.6 (p.305).
Conjugate gradient is computed from gradient, which then gives
the best direction (in N-dim space) in which to proceed to find
the minimum point. The function is then minimized along
this direction of conjugate gradient (a 1-D minimization).
The algorithm is repeated starting at the new point by calling again.
MODIFICATION HISTORY:
Written, Frank Varosi NASA/GSFC 1992.
(See /usr/local/idl/lib/zastron/math/minf_conj_grad.pro)
NAME: minF_parabolic PURPOSE: Find a local minimum of a 1-D function up to specified tolerance. This routine assumes that the function has a minimum nearby. (recommend first calling minF_bracket, xa,xb,xc, to bracket minimum). Routine can also be applied to a scalar function of many variables, for such case the local minimum in a specified direction is found, This routine is called by minF_conj_grad, to locate minimum in the direction of the conjugate gradient of function of many variables. CALLING EXAMPLES: minF_parabolic, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func. or: minF_parabolic, xa,xb,xc, xmin, fmin, FUNC="name", $ POINT=[0,1,1], $ DIRECTION=[2,1,1] ;for 3-D func. INPUTS: xa,xb,xc = scalars, 3 points which bracket location of minimum, that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists. When working with function of N variables (xa,xb,xc) are then relative distances from POINT_NDIM, in the direction specified by keyword DIRECTION, with scale factor given by magnitude of DIRECTION. KEYWORDS: FUNC_NAME = function name (string) Calling mechanism should be: F = func_name( px ) where: px = scalar or vector of independent variables, input. F = scalar value of function at px. POINT_NDIM = when working with function of N variables, use this keyword to specify the starting point in N-dim space. Default = 0, which assumes function is 1-D. DIRECTION = when working with function of N variables, use this keyword to specify the direction in N-dim space along which to bracket the local minimum, (default=1 for 1-D). (xa, xb, xc, x_min are then relative distances from POINT_NDIM) MAX_ITER = maximum allowed number iterations, default=100. TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7). OUTPUTS: xmin = estimated location of minimum. When working with function of N variables, xmin is the relative distance from POINT_NDIM, in the direction specified by keyword DIRECTION, with scale factor given by magnitude of DIRECTION, so that min. Loc. Pmin = Point_Ndim + xmin * Direction. fmin = value of function at xmin (or Pmin). PROCEDURE: Brent's method to minimize a function by using parabolic interpolation, from Numerical Recipes (by Press, et al.), sec.10.2 (p.285). MODIFICATION HISTORY: Written, Frank Varosi NASA/GSFC 1992.
(See /usr/local/idl/lib/zastron/math/minf_parabolic.pro)
NAME: minF_parabol_D PURPOSE: Find a local minimum of a 1-D function up to specified tolerance, using the first derivative of function in the algorithm. This routine assumes that the function has a minimum nearby. (recommend first calling minF_bracket, xa,xb,xc, to bracket minimum). Routine can also be applied to a scalar function of many variables, for such case the local minimum in a specified direction is found, This routine is called by minF_conj_grad, to locate minimum in the direction of the conjugate gradient of function of many variables. CALLING EXAMPLES: minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC_NAME="name" ;for 1-D func. or: minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC="name", $ POINT=[0,1,1], $ DIRECTION=[2,1,1] ;for 3-D func. INPUTS: xa,xb,xc = scalars, 3 points which bracket location of minimum, that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists. When working with function of N variables (xa,xb,xc) are then relative distances from POINT_NDIM, in the direction specified by keyword DIRECTION, with scale factor given by magnitude of DIRECTION. KEYWORDS: FUNC_NAME = function name (string) Calling mechanism should be: F = func_name( px, gradient ) where: px = scalar or vector of independent variables, input. F = scalar value of function at px. gradient = derivative of function, a scalar if 1-D, a gradient vector if N-D, (should only be computed if arg. is present). POINT_NDIM = when working with function of N variables, use this keyword to specify the starting point in N-dim space. Default = 0, which assumes function is 1-D. DIRECTION = when working with function of N variables, use this keyword to specify the direction in N-dim space along which to bracket the local minimum, (default=1 for 1-D). (xa, xb, xc, x_min are then relative distances from POINT_NDIM) MAX_ITER = maximum allowed number iterations, default=100. TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7). OUTPUTS: xmin = estimated location of minimum. When working with function of N variables, xmin is the relative distance from POINT_NDIM, in the direction specified by keyword DIRECTION, with scale factor given by magnitude of DIRECTION, so that min. Loc. Pmin = Point_Ndim + xmin * Direction. fmin = value of function at xmin (or Pmin). PROCEDURE: Brent's method to minimize a function by using parabolic interpolation and using first derivative of function, from Numerical Recipes (by Press, et al.), sec.10.3 (p.287), MODIFICATION HISTORY: Written, Frank Varosi NASA/GSFC 1992.
(See /usr/local/idl/lib/zastron/math/minf_parabol_d.pro)
PCA PURPOSE: Carry out a Principal Components Analysis (Karhunen-Loeve Transform) Results can be directed to the screen, a file, or output variables CALLING SEQUENCE: PCA, data, eigenval, eigenvect, percentages, proj_obj, proj_atr, [MATRIX =, TEXTOUT = ,/COVARIANCE, /SSQ, /SILENT ] INPUT PARAMETERS: data - 2-d data matrix, data(i,j) contains the jth attribute value for the ith object in the sample. If N_OBJ is the total number of objects (rows) in the sample, and N_ATTRIB is the total number of attributes (columns) then data should be dimensioned N_OBJ x N_ATTRIB. OPTIONAL INPUT KEYWORD PARAMETERS: /COVARIANCE - if this keyword set, then the PCA will be carried out on the covariance matrix (rare), the default is to use the correlation matrix /SSQ - if this keyword is set, then the PCA will be carried out on on the sums-of-squares & cross-products matrix (rare) TEXTOUT - Controls print output device, defaults to !TEXTOUT textout=1 TERMINAL using /more option textout=2 TERMINAL without /more option textout=3.prt textout=4 laser.tmp textout=5 user must open file textout = filename (default extension of .prt) OPTIONAL OUTPUT PARAMETERS: eigenval - N_ATTRIB element vector containing the sorted eigenvalues eigenvect - N_ATRRIB x N_ATTRIB matrix containing the corresponding eigenvectors percentages - N_ATTRIB element containing the cumulative percentage variances associated with the principal components proj_obj - N_OBJ by N_ATTRIB matrix containing the projections of the objects on the principal components proj_atr - N_ATTRIB by N_ATTRIB matrix containing the projections of the attributes on the principal components OPTIONAL OUTPUT PARAMETER MATRIX = analysed matrix, either the covariance matrix if /COVARIANCE is set, the "sum of squares and cross-products" matrix is /SSQ is set, or the (by default) correlation matrix. Matrix will have dimensions N_ATTRIB x N_ATTRIB NOTES: This procedure performs Principal Components Analysis (Karhunen-Loeve Transform) according to the method described in "Multivariate Data Analysis" by Murtagh & Heck [Reidel : Dordrecht 1987], pp. 33-48. Keywords /COVARIANCE and /SSQ are mutually exclusive. The printout contains only (at most) the first seven principle eigenvectors. However, the output variables EIGENVECT contain all the eigenvectors Different authors scale the covariance matrix in different ways. The eigenvalues output by PCA may have to be scaled by 1/N_OBJ or 1/(N_OBJ-1) to agree with other calculations when /COVAR is set. PCA uses the non-standard system variables !TEXTOUT and !TEXTUNIT. These can be added to one's session using the procedure ASTROLIB. Users of V3.5.0 or later could change the calls to TQLI and TRED2 to NR_TQLI and NR_TRED2 EXAMPLE: Perform a PCA analysis on the covariance matrix of a data matrix, DATA, and write the results to a file IDL> PCA, data, /COVAR, t = 'pca.dat' Perform a PCA analysis on the correlation matrix. Suppress all printing, and save the eigenvectors and eigenvalues in output varaibles IDL> PCA, data, eigenval, eigenvect, /SILENT PROCEDURES CALLED: Procedures TEXTOPEN, TEXTCLOSE COPYRIGHT: Copyright 1993, Hughes STX Corporation, Lanham MD 20706. REVISION HISTORY: Immanuel Freedman (after Murtagh F. and Heck A.). December 1993 Wayne Landsman, modified I/O December 1993
(See /usr/local/idl/lib/zastron/math/pca.pro)
NAME
POIDEV
PURPOSE:
Return an integer random deviate drawn from a Poisson distribution with
a specified mean. Adapted from procedure of the same name in
"Numerical Recipes" by Press et al. (1986), Section 7.3
CALLING SEQUENCE:
result = POIDEV( xm, [ SEED = ] )
INPUTS:
xm - integer scalar or vector, specifying mean of the Poisson
distribution
OUTPUT:
result - integer scalar or vector, same size as xm
OPTIONAL KEYWORD INPUT-OUTPUT:
SEED - floating point scalar, used as the seed for the random
distribution. This keyword can be used to have POIDEV
give identical results on consecutive runs.
EXAMPLE:
(1) Add Poisson noise to an integral image array, im
IDL> imnoise = POIDEV( im)
(2) Verify the expected mean and sigma for an input value of 81
IDL> p = POIDEV( intarr(10000) + 81) ;Test for 10,000 points
IDL> print,avg(p),sigma(p)
Average and sigma of the 10000 points should be close to 81 and 9
METHOD:
For small values (< 20) independent exponential deviates are generated
until their sum exceeds the specfied mean, the number of events
required is returned as the Poisson deviate. For large (> 20) values,
uniform random variates are compared with a Lorentzian distribution
function.
NOTES:
Negative values in the input array will be returned as zeros.
PROCEDURES CALLED:
GAMMLN - returns log of the Gamma function. Users with IDL V3.0.0
or later may wish to replace this call with the intrinisc
function LNGAMMA which was introduced in IDL 2.4.0.
REVISION HISTORY:
Version 1 Wayne Landsman July 1992
Added SEED keyword September 1992
(See /usr/local/idl/lib/zastron/math/poidev.pro)
NAME:
POLINT
PURPOSE:
Interpolate a set of N points by fitting a polynomial of degree N-1
Adapted from algorithm in Numerical Recipes, Press et al. (1986)
Section 3.1.
CALLING SEQUENCE
POLINT, xa, ya, x, y, [ dy ]
INPUTS:
XA - X Numeric vector, all values must be distinct. The number
of values in XA should rarely exceed 10 (i.e. a 9th order
polynomial)
YA - Y Numeric vector, same number of elements
X - Numeric scalar specifying value to be interpolated
OUTPUT:
Y - Scalar, interpolated value in (XA,YA) corresponding to X
OPTIONAL OUTPUT
DY - Error estimate on Y, scalar
EXAMPLE:
Find sin(2.5) by polynomial interpolation on sin(indgen(10))
IDL> xa = indgen(10)
IDL> ya = sin( xa )
IDL> polint, xa, ya, 2.5, y ,dy
The above method gives y = .5988 & dy = 3.1e-4 a close
approximation to the actual sin(2.5) = .5985
METHOD:
Uses Neville's algorithm to iteratively build up the correct
polynomial, with each iteration containing one higher order.
REVISION HISTORY:
Written W. Landsman January, 1992
(See /usr/local/idl/lib/zastron/math/polint.pro)
NAME: POLY_SMOOTH PURPOSE: Reduce noise in 1-D data (e.g. time-series, spectrum) but retain dynamic range of variations in the data by applying a least squares smoothing polynomial filter, also called the Savitzky-Golay smoothing filter, cf. Numerical Recipes (Press et al. 1992, Sec.14.8) The low-pass filter coefficients are computed by effectively least-squares fitting a polynomial in moving window, centered on each data point, so the new value will be the zero-th coefficient of the polynomial. Approximate first derivates of the data can be computed by using first degree coefficient of each polynomial, and so on. The filter coefficients for a specified polynomial degree and window width are computed independent of any data, and stored in a common block. The filter is then convolved with the data array to result in smoothed data with reduced noise, but retaining higher order variations (better than SMOOTH). CALLING SEQUENCE: spectrum = poly_smooth( data, [ width, DEGREE = , NLEFT = , NRIGHT = DERIV_ORDER = ,COEFF = ] INPUTS: data = 1-D array, such as a spectrum or time-series. width = total number of data points to use in filter convolution, (default = 5, using 2 past and 2 future data points), must be larger than DEGREE of polynomials, and a guideline is make WIDTH between 1 and 2 times the FWHM of desired features. KEYWORDS: DEGREE = degree of polynomials to use in designing the filter via least squares fits, (default DEGREE = 2), and the higher degrees will preserve sharper features. NLEFT = # of past data points to use in filter convolution, excluding current point, overrides width parameter, so that width = NLEFT + NRIGHT + 1. (default = NRIGHT) NRIGHT = # of future data points to use (default = NLEFT). DERIV_ORDER = order of derivative desired (default = 0, no derivative). COEFFICIENTS = optional output of the filter coefficients applied, but they are all stored in common block for reuse, anyway. RESULTS: Function returns the data convolved with polynomial filter coefs. EXAMPLE: Given a wavelength - flux spectrum (w,f), apply a 31 point quadratic smoothing filter and plot IDL> plot, w, poly_smooth(f,31) COMMON BLOCKS: common poly_smooth, degc, nlc, nrc, coefs, ordermax PROCEDURE: As described in Numerical Recipies, 2nd edition sec.14.8, Savitsky-Golay filter. Matrix of normal eqs. is formed by starting with small terms and then adding progressively larger terms (powers). The filter coefficients of up to derivative ordermax are stored in common, until the specifications change, then recompute coefficients. Coefficients are stored in convolution order, zero lag in the middle. MODIFICATION HISTORY: Written, Frank Varosi NASA/GSFC 1993.
(See /usr/local/idl/lib/zastron/math/poly_smooth.pro)
NAME: PROB_KS PURPOSE: Returns the significance level of an observed value of the Kolmogorov-Smirnov statistic D for an effective number of data points N_eff. Called by KSONE and KSTWO CALLING SEQUENCE: prob_ks, D, N_eff, probks INPUT PARAMATERS: D - Kolmogorov statistic, floating scalar, always non-negative N_eff - Effective number of data points, scalar. For a 2 sided test this is given by (N1*N2)/(N1+N2) where N1 and N2 are the number of points in each data set. OUTPUT PARAMETERS: probks - floating scalar between 0 and 1 giving the significance level of the K-S statistic. Small values of PROB suggest that the distribution being tested are not the same REVISION HISTORY: Written W. Landsman August, 1992
(See /usr/local/idl/lib/zastron/math/prob_ks.pro)
NAME:
QSIMP
PURPOSE:
Integrate a function to specified accuracy using the extended
trapezoidal rule. Adapted from algorithm in Numerical Recipes,
by Press et al. (1986) p. 111. This procedure became mostly
obsolete in IDL V3.5 with the introduction of the intrinsic
procedure NR_QSIMP.
CALLING SEQUENCE:
QSIMP, func, A, B, S, [ EPS = , MAX_ITER = ]
INPUTS:
func - scalar string giving name of function of one variable to
be integrated
A,B - numeric scalars giving the lower and upper bound of the
integration
OUTPUTS:
S - Scalar giving the approximation to the integral of the specified
function between A and B.
OPTIONAL KEYWORD PARAMETERS:
EPS - scalar specify the fractional accuracy before ending the
iteration. Default = 1E-6
MAX_ITER - Integer specifying the total number iterations at which
QSIMP will terminate even if the specified accuracy has not yet
been met. The maximum number of function evaluations will be
2^(MAX_ITER). Default value is MAX_ITER = 20
NOTES:
The function QTRAP is robust way of doing integrals that are not very
smooth. However, if the function has a continuous 3rd derivative then
QSIMP will likely be more efficient at performing the integral.
EXAMPLE:
Compute the integral of sin(x) from 0 to !PI/3.
IDL> QSIMP, 'sin', 0, !PI/3, S & print, S
The value obtained should be cos(!PI/3) = 0.5
PROCEDURES CALLED:
TRAPZD, ZPARCHECK
REVISION HISTORY:
W. Landsman ST Systems Co. August, 1991
(See /usr/local/idl/lib/zastron/math/qsimp.pro)
NAME:
QTRAP
PURPOSE:
Integrate a function to specified accuracy using the extended trapezoidal
rule. Adapted from Numerical Recipes, p. 111.
CALLING SEQUENCE:
QTRAP, func, A, B, S, [EPS = , MAX_ITER = ]
INPUTS:
func - scalar string giving name of function of one variable to
be integrated
A,B - numeric scalars giving the lower and upper bound of the
integration
OUTPUTS:
S - Scalar giving the approximation to the integral of the specified
function between A and B.
OPTIONAL KEYWORD PARAMETERS:
EPS - scalar specify the fractional accuracy before ending the
iteration. Default = 1E-6
MAX_ITER - Integer specifying the total number iterations at which
QTRAP will terminate even if the specified accuracy has not yet
been met. The maximum number of function evaluations will
be 2^(MAX_ITER). Default value is MAX_ITER = 20
NOTES:
QTRAP is robust way of doing integrals that are not very smooth. If the
function has a continuous 3rd derivative then the function QSIMP will
likely be more efficient at performing the integral.
EXAMPLE:
Compute the integral of sin(x) from 0 to !PI/3.
IDL> QTRAP, 'sin', 0, !PI/3, S & print,S
The value obtained should be cos(!PI/3) = 0.5
PROCEDURES CALLED:
TRAPZD, ZPARCHECK
REVISION HISTORY:
W. Landsman ST Systems Co. August, 1991
(See /usr/local/idl/lib/zastron/math/qtrap.pro)
NAME:
QUADTERP
PURPOSE:
Quadratically interpolate (3 point Lagrangian) a function Y = f(X)
at specified grid points. Use LINTERP for linear interpolation
CALLING SEQUENCE:
QUADTERP, Xtab, Ytab, Xint, Yint, [ MISSING = ]
INPUT:
Xtab - Vector (X TABle) containing the current independent variable
Must be either monotonic increasing or decreasing
Ytab - Vector (Y TABle) containing the dependent variable defined
at each of the points of XTAB.
Xint - Scalar or vector giving the values of X for which interpolated
Y values are sought
OUTPUT:
Yint - Interpolated value(s) of Y, same number of points as Xint
OPTIONAL INPUT KEYWORD:
MISSING - Scalar specifying Yint value(s) to be assigned, when Xint
value(s) are outside of the range of Xtab. Default is to
truncate the out of range Yint value(s) to the nearest value
of Ytab. See the help for the INTERPOLATE function.
METHOD:
3-point Lagrangian interpolation. The average of the two quadratics
derived from the four nearest points is returned in YTAB
TABINV is used to locate center point of the interpolation.
RESTRICTIONS:
Unless MISSING keyword is set, points outside the range of Xtab in
which two valid quadratics can be computed are returned at the value
of the next to last point of Ytab (i.e. Ytab(1) and Ytab(NPTS-2) ).
EXAMPLE:
A spectrum has been defined using a wavelength vector WAVE and a
flux vector FLUX. Interpolate onto a new wavelength grid, e.g.
IDL> wgrid = [1540.,1541.,1542.,1543.,1544.,1545.]
IDL> quadterp, wave, flux, wgrid, fgrid
FGRID will be a 5 element vector containing the quadratically
interpolated values of FLUX at the wavelengths given in WGRID.
EXTERNAL ROUTINES:
TABINV, ZPARCHECK, DATATYPE
REVISION HISTORY:
31 October 1986 by B. Boothman, adapted from the IUE RDAF
12 December 1988 J. Murthy, corrected error in Xint
September 1992, W. Landsman, fixed problem with double precision
August 1993, W. Landsman, added MISSING keyword
(See /usr/local/idl/lib/zastron/math/quadterp.pro)
NAME: SIGMA PURPOSE: Calculate the standard deviation value of an array, or calculate the standard deviation over one dimension of an array as a function of all the other dimensions. CALLING SEQUENCE: RESULT = SIGMA(ARRAY) RESULT = SIGMA(ARRAY,N_PAR) RESULT = SIGMA(ARRAY,N_PAR,DIM) INPUTS: ARRAY = Input array. May be any type except string. OPTIONAL INPUT PARAMETERS: N_PAR = Number of parameters. Default value is zero. The number of degrees of freedom is N_FREE = N_POINTS - N_PAR, where N_POINTS is either N_ELEMENTS(ARRAY) or the size of dimension DIM. The value of SIGMA varies as 1 / SQRT(N_FREE). If N_PAR is negative, then the absolute value of N_PAR is used, and an additional factor of 1 / SQRT(N_POINTS) is included, yielding the reduced sigma. DIM = Optional dimension to do standard deviation over. OUTPUTS: The standard deviation value of the array when called with one parameter. If DIM is passed, then the result is an array with all the dimensions of the input array except for the dimension specified, each element of which is the standard deviation of the corresponding vector in the input array. For example, if A is an array with dimensions of (3,4,5), then the command B = SIGMA(A,N,1) is equivalent to B = FLTARR(3,5) FOR J = 0,4 DO BEGIN FOR I = 0,2 DO BEGIN B(I,J) = SIGMA( A(I,*,J) , N ) ENDFOR ENDFOR RESTRICTIONS: Dimension specified must be valid for the array passed; otherwise the input array is returned as the output array. PROCEDURE: Uses the function AVG. When DIM is passed, then the function SUM is used. MODIFICATION HISTORY: William Thompson Applied Research Corporation July, 1986 8201 Corporate Drive Landover, MD 20785
(See /usr/local/idl/lib/zastron/math/sigma.pro)
NAME: SIXLIN PURPOSE: Compute linear regression coefficients by six different methods. Adapted from the FORTRAN program (Rev. 1.1) supplied by Isobe, Feigelson, Akritas, and Babu Ap. J. Vol. 364, p. 104 (1990). Suggested when there is no understanding about the nature of the scatter about a linear relation, and NOT when the errors in the variable are calculable. CALLING SEQUENCE: SIXLIN, xx, yy, a, siga, b, sigb INPUTS: XX - vector of X values YY - vector of Y values, same number of elements as XX OUTPUTS: A - Vector of 6 Y intercept coefficients SIGA - Vector of standard deviations of 6 Y intercepts B - Vector of 6 slope coefficients SIGB - Vector of standard deviations of slope coefficients The output variables are computed using linear regression for each of the following 6 cases: (0) Ordinary Least Squares (OLS) Y vs. X (1) Ordinary Least Squares X vs. Y (2) Ordinary Least Squares Bisector (3) Orthogonal Reduced Major Axis (4) Reduced Major-Axis (5) Mean ordinary Least Squares NOTES: Isobe et al. make the following recommendations (1) If the different linear regression methods yield similar results then quoting OLS(Y|X) is probably the most familiar. (2) If the linear relation is to be used to predict Y vs. X then OLS(Y|X) should be used. (3) If the goal is to determine the functional relationship between X and Y then the OLS bisector is recommended. REVISION HISTORY: Written Wayne Landsman February, 1991 Corrected sigma calculations February, 1992
(See /usr/local/idl/lib/zastron/math/sixlin.pro)
SPLINE_SMOOTH
PURPOSE:
Construct cubic smoothing spline (or give regression solution) to given
data with minimum "roughness" (measured by the energy in the second
derivatives) whilst restricting the weighted mean square distance
of the approximation from the data. The results may be written to
the screen or a file or both and are optionally returned in the
parameters. The results may be optionally displayed graphically.
CATEGORY:
Math and Statistics [PRO.MATH].
CALLING SEQUENCE:
SPLINE_SMOOTH,X,Y,Yerr,distance,[coefficients,smoothness,xplot,yplot XTITLE=xtitle,YTITLE=ytitle, INTERP=interp, TEXTOUT=,/SILENT,/PLOT,/ERRBAR]
INPUT PARAMETERS:
X - N_POINT element vector containing the data abcissae
Y - N_POINT element vector containing the data ordinates
Yerr - estimated uncertainty in ordinates ( positive)
distance - upper bound on the weighted mean square distance
of the approximation from the data
OPTIONAL INPUT PARAMETERS
xplot - vector of spline evaluation abcissae
OPTIONAL INPUT KEYWORD PARAMETERS:
TEXTOUT - Controls print output device, defaults to !TEXTOUT
textout=1 TERMINAL using /more option
textout=2 TERMINAL without /more option
textout=3 .prt
textout=4 laser.tmp
textout=5 user must open file
textout = filename (default extension of .prt)
OPTIONAL OUTPUT PARAMETERS:
coefficients - N_POINT x 4 element array containing the sequence of
spline coefficients including the smoothed ordinates.
smoothness - N_POINT element vector containing the energy in second
derivatives of approximated function.
yplot - vector of evaluated spline ordinates.
OPTIONAL OUTPUT KEYWORD PARAMETERS
/SILENT - suppress all printing.
/PLOT - display smooth curve, data ordinates and error bars
/ERRBAR - display error bars
XTITLE - optional title for X-axis
YTITLE - optional title for Y-axis
INTERP - optionally returned interpolated smooth spline
NOTES:
This procedure constructs a smoothing spline according to the method
described in "Fundamentals of Image Processing" by A. Jain [Prentice-
Hall : New Jersey 1989].
If the distance parameter is sufficiently large a linear regression
is performed, otherwise a cubic smoothing spline is constructed.
This procedure assumes regular sampling and independent identically
distributed normal errors without missing data. The data are sorted.
SPLINE_SMOOTH uses the non-standard system variables !TEXTOUT and
!TEXTUNIT.
These can be added to one's session using the procedure ASTROLIB.
COMMON BLOCKS:
None.
EXAMPLE:
Obtain coefficients of a univariate smoothing spline fitted to data
X,Y assuming normally distributed errors Yerr and write the results to
a file.
IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, smoothness,
t='spline.dat'
Fit a smoothing spline to observational data. Suppress all printing
and save the smoothed ordinates in output variables. Display results.
IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, /SILENT, /PLOT
PROCEDURES CALLED:
Procedures TEXTOPEN, TEXTCLOSE, PLOT, PLOTERR
COPYRIGHT:
Copyright 1993, Hughes STX Corporation, Lanham MD 20706.
AUTHOR:
Immanuel Freedman (after A. Jain). December, 1993
REVISIONS
January 12, 1994 I. Freedman (HSTX) Adjusted formats
March 14, 1994 I. Freedman (HSTX) Improved convergence
March 15, 1994 I. Freedman (HSTX) User-specified interpolates
(See /usr/local/idl/lib/zastron/math/spline_smooth.pro)
NAME:
TABINV
PURPOSE:
To find the effective index of a function value in
an ordered vector.
CALLING SEQUENCE:
TABINV, XARR, X, IEFF
INPUTS:
XARR - the vector array to be searched, must be monotonic
increasing or decreasing
X - the function value(s) whose effective
index is sought (scalar or vector)
OUTPUT:
IEFF - the effective index or indices of X in XARR
real or double precision, same # of elements as X
RESTRICTIONS:
TABINV will abort if XARR is not monotonic. (Equality of
neighboring values in XARR is allowed but results may not be
unique.) This requirement may mean that input vectors with padded
zeroes could cause routine to abort.
PROCEDURE:
A binary search is used to find the values XARR(I)
and XARR(I+1) where XARR(I) < X < XARR(I+1).
IEFF is then computed using linear interpolation
between I and I+1.
IEFF = I + (X-XARR(I)) / (XARR(I+1)-XARR(I))
Let N = number of elements in XARR
if x < XARR(0) then IEFF is set to 0
if x > XARR(N-1) then IEFF is set to N-1
EXAMPLE:
Set all flux values of a spectrum (WAVE vs FLUX) to zero
for wavelengths less than 1150 Angstroms.
IDL> tabinv, wave, 1150.0, I
IDL> flux( 0:fix(I) ) = 0.
FUNCTIONS CALLED:
ISARRAY
REVISION HISTORY:
Adapted from the IUE RDAF January, 1988
More elegant code W. Landsman August, 1989
Mod to work on 2 element decreasing vector August, 1992
(See /usr/local/idl/lib/zastron/math/tabinv.pro)
NAME TRAPZD PURPOSE: Compute the nth stage of refinement of an extended trapezoidal rule. This procedure is called by QSIMP and QTRAP. Algorithm from Numerical Recipes, Section 4.2. TRAPZD is meant to be called iteratively from a higher level procedure. CALLING SEQUENCE: TRAPZD, func, A, B, S, step INPUTS: func - scalar string giving name of function to be integrated. This must be a function of one variable. A,B - scalars giving the limits of the integration INPUT-OUTPUT: S - scalar giving the total sum from the previous interations on input and the refined sum after the current iteration on output. step - LONG scalar giving the number of points at which to compute the function for the current iteration. If step is not defined on input, then S is intialized using the average of the endpoints of limits of integration. NOTES: TRAPZD will check for math errors when computing the function at the endpoints, but not on subsequent iterations. REVISION HISTORY: Written W. Landsman August, 1991
(See /usr/local/idl/lib/zastron/math/trapzd.pro)
NAME:
TSUM
PURPOSE:
Trapezoidal summation of the area under a curve. This procedure
was formerly know as INTEG
CALLING SEQUENCE:
Result = TSUM(y)
or
Result = TSUM( x, y, [ imin, imax ] )
INPUTS:
x = array containing independent variable. If omitted, then
x is assumed to contain the index of the y variable.
x = indgen( N_elements(y) ).
y = array containing dependent variable y = f(x)
OPTIONAL INPUTS:
imin = index of x array at which to begin the integration, integer
scalar. If omitted, then summation starts at x(0).
imax = index of x value at which to end the integration, integer
scalar. If omitted then the integration ends at x(npts).
OUTPUTS:
result = area under the curve y=f(x) between x(imin) and x(imax).
PROCEDURE:
The area is determined of indivdual trapezoids defined by x(i),
x(i+1), y(i) and y(i+1).
MODIFICATION HISTORY:
Written, W.B. Landsman, STI Corp. May 1986
Modified so X is not altered in a one parameter call Jan 1990
(See /usr/local/idl/lib/zastron/math/tsum.pro)
NAME:
ZBRENT
PURPOSE:
Find the zero of a 1-D function up to specified tolerance.
This routine assumes that the function is known to have a zero.
CALLING:
x_zero = ZBRENT( x1, x2, FUNC_NAME="name" )
INPUTS:
x1, x2 = scalars, 2 points which bracket location of function zero,
that is, F(x1) < 0 < F(x2).
Note: computations are performed with
same precision (single/double) as the inputs and user supplied function.
REQUIRED INPUT KEYWORD:
FUNC_NAME = function name (string)
Calling mechanism should be: F = func_name( px )
where: px = scalar independent variable, input.
F = scalar value of function at px,
should be same precision (single/double) as input.
OPTIONAL INPUT KEYWORD:
MAX_ITER = maximum allowed number iterations, default=100.
TOLERANCE = desired accuracy of minimum location, default = 1.e-3.
OUTPUTS:
Returns the location of zero, with accuracy of specified tolerance.
PROCEDURE:
Brent's method to find zero of a function by using bracketing,
bisection, and inverse quadratic interpolation,
from Numerical Recipes (by Press, et al.), sec.9.3 (p.251).
EXAMPLE:
Find the root of the COSINE function between 1. and 2. radians
IDL> print, zbrent( 1, 2, FUNC = 'COS')
and the result will be !PI/2 within the specified tolerance
MODIFICATION HISTORY:
Written, Frank Varosi NASA/GSFC 1992.
FV.1994, mod to check for single/double prec. and set zeps accordingly.
(See /usr/local/idl/lib/zastron/math/zbrent.pro)