Common Lisp Package: SAPA

README:

FUNCTION

Public

2D-MATRIX->LINEARIZED-UPPER-TRANGULAR (A &KEY (RESULT (LET ((N (CAR (ARRAY-DIMENSIONS A)))) (MAKE-ARRAY (/ (* N (1+ N)) 2)))))

given [1] A (required) ==> a 2d hermitian matrix [2] result (keyword; new vector of appropriate size) <== a vector filled linearly with elements of A converts 2d hermitian array into linearized form and returns [1] result, a vector with 2d hermitian array stored in linearized form

2D-MATRIX-MOVE! (FROM-THIS TO-THIS)

given [1] from-this (required) ==> a 2d matrix [2] to-this (required) <== another 2d matrix transfer contents of from-this to corresponding locations in to-this, and returns [1] to-this

A*X (A X)

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers returns [1] a vector of numbers obtained by multiplying each element of x by a ---- Note: corresponds to SSCAL in LINPACK's bla.

A*X! (A X &KEY (RESULT X))

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers [2] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by multiplying each element of x by a ---- Note: corresponds to SSCAL in LINPACK's bla.

A*X+B (A X B)

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers [3] b (required) ==> a number returns [1] a vector of numbers obtained by multiplying each element of x by a and adding b

A*X+B! (A X B &KEY (RESULT X))

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers [3] b (required) ==> a number [4] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by multiplying each element of x by a and adding b

A*X+Y (A X Y)

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers [3] y (required) ==> a sequence of numbers returns [1] a vector of numbers obtained by multiplying each element of x by a and adding the corresponding element of y ---- Note: corresponds to SAXPY in LINPACK's bla.

A*X+Y! (A X Y &KEY (RESULT Y))

given [1] a (required) ==> a number [2] x (required) ==> a sequence of numbers [3] y (required) ==> a sequence of numbers [4] result (keyword; y) <== a sequence of numbers returns [1] a sequence of numbers obtained by multiplying each element of x by a and adding the corresponding element of y ---- Note: corresponds to SAXPY in LINPACK's bla.

ACVS (TIME-SERIES &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA-P T) (ACS-P NIL) (RESULT (MAKE-ARRAY (- END START))))

given [1] time-series (required) ==> a vector containing the time series [2] start (keyword; 0) ==> start index of time-series to be used [3] end (keyword; length of time-series) ==> 1 + end index of time-series to be used [4] center-data-p (keyword; t) ==> if t, subtract sample mean of time series prior to computing the acvs; if nil, do not subtract the sample mean [5] acs-p (keyword; nil) ==> return autocorrelation sequence rather than autocovariance sequence [6] result (keyword; vector of length n) <== vector to hold acvs (or acs) calculates autocovariance (or autocorrelation) sequence using fft's and returns [1] result, a vector with the sample acvs (or acs) [2] sample variance of time series (with mean treated as specified by center-data-p keyword) --- Note: see Equations (191a) and (190b) of the SAPA book

ACVS-FOR-TIME-SERIES-SIMULATED-FROM-SDF (N-SERIES SDF &KEY (SAMPLING-TIME 1.0) (N-TOTAL (* 4 (NEXT-POWER-OF-2 N-SERIES))) (RESULT (MAKE-ARRAY N-SERIES)))

given [1] n-series (required) ==> length of time series simulated using simulate-time-series-from-sdf (must be a power of 2) [2] sdf (required) ==> spectral density function (sdf) defined for 0 <= f <= 1/(2.0 sampling-time) -- the sdf is assumed to be two-sided and symmetric about f = 0 [3] sampling-time (keyword; 1.0) ==> the assumed sampling time (delta t) [4] n-total (keyword; 4 * next power of 2 for n-series) ==> a power of 2 controlling degree of accuracy of the approximation (the larger, the better -- see Percival, 1992, for details) [5] result (keyword; vector of length n-series) <== a vector to contain simulated time series returns [1] a vector with the theoretical acvs from lag 0 to n-series - 1 for a time series generated via a call to simulate-time-series-from-sdf with n-series, sdf, sampling-time and n-total set to the same values

ADD-SEQUENCES (&REST SEQS)

given [1] a arbitrary number of sequences of numbers, all of the same size returns [1] a new sequence formed by adding the sequences together on an element by element basis

AR-COEFFS->ACVS (AR-COEFFS VARIANCE MAX-LAG &KEY (PROCESS-VARIANCE-P T) (LIST-OF-LOWER-ORDER-PHI NIL) (LIST-OF-LOWER-ORDER-PEV NIL) (RESULT (MAKE-ARRAY (1+ MAX-LAG))))

given [1] AR-coeffs (required) ==> vector of length p with real or complex-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [2] variance (required) ==> process variance (if process-variance-p is true) or innovations variance (if nil) [3] max-lag (required) ==> acvs to be computed out to this lag [4] process-variance-p (keyword; t) ==> if true, required variance parameter is taken to be process variance; otherwise, it is taken to be the innovations variance (both variance and process-variance-p are ignored if the next two keyword parameters are supplied) [5] list-of-lower-order-phi (keyword; nil) ==> if supplied, bypasses call to step-down-Levinson-Durbin-recursions, in which case list-of-lower-order-pev must be supplied also [6] list-of-lower-order-pev (keyword; nil) ==> if supplied, bypasses call to step-down-Levinson-Durbin-recursions [7] result (keyword; vector of size max-lag + 1) <== vector to hold acvs for lags 0, 1, ..., max-lag returns [1] result, i.e., the acvs --- Note: see Section 9.4 of the SAPA book

AR-COEFFS->PREWHITENING-FILTER (AR-COEFFS &KEY (RESULT (MAKE-ARRAY (1+ (LENGTH AR-COEFFS)))))

given [1] AR-coeffs (required) ==> a sequence of AR coefficients phi_{1,p}, ..., phi_{p,p} [2] result (keyword; vector of size (1+ (length AR-coeffs))) <== sequence to hold coefficients for prewhitening filter 1, -phi_{1,p}, ..., -phi_{p,p} returns [1] result, the prewhitening filter --- Note: this filter can be used to compute forward prediction errors; to obtain the filter needed to compute backward prediction errors, reverse the elements in result; i.e., evaluate (reverse result)

AR-COEFFS->REFLECTION-COEFFS (AR-COEFFS &KEY (LIST-OF-LOWER-ORDER-PHI (STEP-DOWN-LEVINSON-DURBIN-RECURSIONS AR-COEFFS 1.0)) (RESULT (MAKE-ARRAY (LENGTH AR-COEFFS))))

given [1] AR-coeffs (required) ==> vector of length p with real or complex-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [2] list-of-lower-order-phi (keyword; calls step-down-Levinson-Durbin-recursions) ==> must contain the list that is returned by the function step-down-Levinson-Durbin-recursions; this keyword parameter is useful if the result of calling that function is already available [3] results (keyword; vector of length p) <== vector of size p into which the reflection coefficients phi_{j,j} of orders j = 1, 2, ..., p are placed returns [1] results, i.e., the reflection coefficients --- Note: see Section 9.4 of the SAPA book

AR-COEFFS->SDF (AR-COEFFS INNOVATIONS-VARIANCE &KEY (N-NONZERO-FREQS 256) (SAMPLING-TIME 1.0) (RETURN-EST-FOR-0-FREQ-P T) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (IF RETURN-EST-FOR-0-FREQ-P (1+ N-NONZERO-FREQS) N-NONZERO-FREQS))) (RETURN-FREQUENCIES-P NIL) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (IF RETURN-EST-FOR-0-FREQ-P (1+ N-NONZERO-FREQS) N-NONZERO-FREQS)))))

given [1] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [2] innovations-variance (required) ==> innovations variance of the process [3] N-nonzero-freqs (keyword; 256) ==> the number of nonzero frequencies at which the AR sdf is to be computed (need NOT be a power of 2); the first nonzero frequency (and also spacing between frequencies) is given by 1/(2 * N-nonzero-freqs * sampling time) = Nyquist frequency/N-nonzero-freqs; the last nonzero frequency is the Nyquist frequency [4] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [5] return-est-for-0-freq-p (keyword; t) ==> if t, sdf is computed at zero frequency; otherwise, it is not computed. [6] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [7] result-sdf (keyword; vector of correct length) <== vector of length N-nonzero-freqs (if return-est-for-0-freq-p is nil) or N-nonzero-freqs +1 (if return-est-for-0-freq-p is t) into which the properly transformed sdf is placed [8] return-frequencies-p (keyword; nil) ==> if t, the frequencies associated with the transfer function are computed and returned in result-freq [9] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [10] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length N-nonzero-freqs (if return-est-for-0-freq-p is nil) or N-nonzero-freqs +1 (if return-est-for-0-freq-p is t) into which the frequencies associated with the values in result-sdf are placed returns [1] result-sdf, a vector holding the properly transformed sdf [2] nil (if return-frequencies-p is nil) or result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-sdf [3] the length of the vector result-sdf --- Note: see Equation (392b) of the SAPA book

AR-COEFFS->SDF-AT-SINGLE-FREQ (AR-COEFFS INNOVATIONS-VARIANCE FREQ &KEY (SAMPLING-TIME 1.0) (SDF-TRANSFORMATION #'CONVERT-TO-DB))

given [1] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [2] innovations-variance (required) ==> innovations variance of the process [3] freq (required) ==> the single frequency at which the sdf is to be computed [4] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [6] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform the value of the sdf at freq returns [1] properly transformed AR sdf evaluated at freq [2] nil (if return-frequencies-p is nil) or result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-sdf [3] the length of the vector result-sdf --- Note: see Equation (392b) of the SAPA book

AR-COEFFS->VARIANCE (AR-COEFFS INNOVATIONS-VARIANCE)

given [1] AR-coeffs (required) ==> a sequence of AR coefficients phi_{1,p}, ..., phi_{p,p} [2] innovations-variance (required) ==> innovations variance of the process returns [1] variance of the process

BANDWIDTH&CONFIDENCE-INTERVALS-FOR-SDF-DB (SDF-DB FREQ SAMPLE-SIZE &KEY (CONFIDENCE-LEVEL 0.95) (LAG-WINDOW-FUNCTION NIL) (MAX-LAG (1- SAMPLE-SIZE)) (SAMPLING-TIME 1.0) (C_H 1.0))

given [1] sdf-dB (required) ==> an sdf estimate (in decibels) for a single frequency [2] freq (required) ==> the frequency associated with sdf-dB [3] sampling-size (required) ==> the sample size of the time series from which sdf-dB was computed [4] confidence-level (keyword; 0.95) ==> the level of the confidence interval (e.g., 0.95 yields a 95% confidence interval) [5] lag-window-function (keyword; nil) ==> the lag window function used to create sdf-dB (nil means a lag window function was NOT used) [6] max-lag (keyword; sample-size - 1) ==> the maximum lag used in conjunction with lag-window-function to create sdf-dB (not used in lag-window-function is nil) [7] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [8] C_h (keyword; 1.0) ==> the variance inflation factor due to tapering (see Table 248 of the SAPA book); the default value of 1.0 is appropriate for a rectangular data taper returns [1] left-hand size of interval describing the bandwidth of the spectral estimate (interval is centered on freq) [2] right-hand size of bandwidth interval [3] lower limit (in decibels) of confidence interval for true sdf based upon sdf-dB [4] upper limit (in decibels) of confidence interval --- Note: see Sections 6.7 and 6.10 of the SAPA book

BARTLETT-BANDWIDTH->M (B_W &KEY (SAMPLING-TIME 1.0))

given desired smoothing window bandwidth B_W and sampling time, returns [1] window parameter m required to approximately achieve B_W using the Bartlett lag window [2] actual B_W achieved --- Note: see Table 269 of the SAPA book

BARTLETT-LAG-WINDOW (TAU M)

given the lag tau and window parameter m, returns the value of the Bartlett lag window --- Note: see Equation (260) of the SAPA book or Priestley, page 439, Equation (6.2.65)

BARTLETT-M->BANDWIDTH (M &KEY (SAMPLING-TIME 1.0))

given window parameter m and sampling time, returns bandwidth B_W for the Bartlett smoothing window --- Note: see Table 269 of the SAPA book

BARTLETT-N-M->DEGREES-OF-FREEDOM (N M &KEY (C_H 1.0))

given sample size N, window parameter m and variance inflation factor C_h, returns equivalent degrees of freedom nu for Bartlett lag window --- Note: see Table 269 of the SAPA book

BIASED-ACVS->UNBIASED-ACVS (BIASED-ACVS)

given a biased estimate of the acvs, returns the corresponding unbiased estimate

BISECTION-WITH-NEWTON-RAPHSON (F F-PRIME X-LEFT X-RIGHT &KEY (ACCURACY (* 10.0 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 100))

given [1] f (required) ==> a function with a single argument [2] f-prime (required) ==> another function with a single argument, this one being the first derivative of f [3] x-left (required) ==> left-hand bracket for the desired root; i.e., left-hand bracket <= desired root [4] x-right (required) ==> right-hand bracket for the desired root; i.e., desired root <= right-hand bracket [5] accuracy (keyword; (* 10.0 single-float-epsilon)) ==> desired relative accuracy for computed rood [6] maximum-number-of-iterations (keyword; 100) ==> maximum number of iterations returns [1] a root of f in [x-left, x-right]; i.e., a value x in that interval such that f(x) = 0 [2] the number of iterations required --- Note: this function is based loosely on rtsafe, Section 9.4, Numerical Recipes, Second Edition

BOX-PLOT (A-SEQ &KEY (SEQ-SORTED-P NIL) (GET-BOTTOM-LINE #'(LAMBDA (V) (LET* ((LOWER-QUARTILE (QUANTILE-OF-ORDERED-SEQ V 0.25)) (IQR (- (QUANTILE-OF-ORDERED-SEQ V 0.75) LOWER-QUARTILE)) (LOWER-LIMIT (- LOWER-QUARTILE (* 1.5 IQR))) (TEMP (BINARY-SEARCH LOWER-LIMIT V)) (LOWER-INDEX (IF TEMP (1+ TEMP) 0)) (LOWER-ADJACENT-VALUE (ELT V LOWER-INDEX))) (VALUES LOWER-ADJACENT-VALUE LOWER-INDEX)))) (GET-BOTTOM-OF-BOX #'(LAMBDA (V) (QUANTILE-OF-ORDERED-SEQ V 0.25))) (GET-LINE-THROUGH-BOX #'(LAMBDA (V) (QUANTILE-OF-ORDERED-SEQ V 0.5))) (GET-TOP-OF-BOX #'(LAMBDA (V) (QUANTILE-OF-ORDERED-SEQ V 0.75))) (GET-TOP-LINE #'(LAMBDA (V) (LET* ((UPPER-QUARTILE (QUANTILE-OF-ORDERED-SEQ V 0.75)) (IQR (- UPPER-QUARTILE (QUANTILE-OF-ORDERED-SEQ V 0.25))) (UPPER-LIMIT (+ UPPER-QUARTILE (* 1.5 IQR))) (TEMP (BINARY-SEARCH UPPER-LIMIT V)) (UPPER-INDEX (IF TEMP TEMP (1- (LENGTH V)))) (UPPER-ADJACENT-VALUE (ELT V UPPER-INDEX))) (VALUES UPPER-ADJACENT-VALUE UPPER-INDEX)))))

given [1] a-seq (required) ==> a sequence of real-valued numbers [2] seq-sorted-p (keyword; nil) ==> if t, a-seq is already sorted; if nil, a copy of a-seq is sorted for use by the function [3] get-bottom-line (keyword; computes lower adjacent value) ==> a function for calculating the line below the box (from sorted data) [4] get-bottom-of-box (keyword; computes lower quartile) ==> a function for calculating bottom of the box [5] get-line-through-box (keyword; computes median) ==> a function for calculating the line through the middle of the box [6] get-top-of-box (keyword; computes upper quartile) ==> a function for calculating top of the box [7] get-top-line (keyword; computes upper adjacent value) ==> a function for calculating the line above the box returns [1] bottom line below box [2] bottom line of box [3] line through middle of box [4] top line of box [5] top line above box [6] a vector containing outside values (note: this vector can be of length 0) --- Note: the values that this function returns can be used to contruct a box plot as defined in Section 2.5 of ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, Tukey

BURG-ALGORITHM (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)) (APPROXIMATE-MSES (MAKE-ARRAY (1+ P))) (REFLECTION-COEFFS (MAKE-ARRAY P)) (FORWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))) (BACKWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))) (EXACT-MSES (MAKE-ARRAY (1+ P))))

given [1] time-series (required) ==> a vector containing a real-valued or complex-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} [7] approximate-MSEs (keyword; a vector of length p+1) <== estimates of innovations variance (mean square errors) for models of order 0, 1, ..., p (the innovations variance for the zeroth order model is the sample variance) [8] reflection-coeffs (keyword; a vector of length p) <== estimates of reflection coefficients phi_{1,1}, ..., phi_{p,p} [9] forward-prediction-errors (keyword; a vector of length end - start - p) <== computed forward prediction errors [10] backward-prediction-errors (keyword; a vector of length end - start - p) <== computed backward prediction errors [11] exact-MSEs (keyword; a vector of length p+1) <== another set of estimates of innovations variance for models of order 0, 1, ..., p; these estimates are based on Equation (419a) of the SAPA book uses Burg's algorithm to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's; this number is also given in (elt approximate-MSEs p) [3] approximate-MSEs [4] reflection-coeffs [5] forward-prediction-errors [6] backward-prediction-errors [7] exact-MSEs --- Note: see Section 9.5 of the SAPA book

CAREFUL-CONVERT-TO-DB (NUM &OPTIONAL (ZERO-MAPPING -100.0))

given [1] num (required) ==> a number [1] zero-mapping (required) ==> a number returns [1] 10 log_10(num) if num > 0 or zero-mapping if num <= 0

CENTER&PREWHITEN&TAPER-TIME-SERIES (TIME-SERIES &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (PREWHITENING-FILTER NIL) (RECENTER-AFTER-PREWHITENING-P T) (DATA-TAPER NIL) (DATA-TAPER-PARAMETERS) (RECENTER-AFTER-TAPERING-P T) (RESTORE-POWER-OPTION-P T) (RESULT (MAKE-ARRAY (IF PREWHITENING-FILTER (- END START (1- (LENGTH PREWHITENING-FILTER))) (- END START)))))

given [1] time-series (required) ==> a sequence of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [3] start (keyword; 0) ==> start index of sequence to be used [4] end (keyword; length of time-series) ==> 1 + end index of sequence to be used [5] prewhitening-filter (keyword; nil) ==> vector with coeffients of prewhitening filter or nil (if no prewhitening is to be done) [6] recenter-after-prewhitening-p (keyword; t) ==> if t, prewhitened series is centered using sample mean (not used if prewhitening-filter is nil) [7] data-taper (keyword; nil) ==> nil or a tapering function [8] data-taper-parameters (keyword) ==> parameters for tapering function (not used if data-taper is nil) [9] recenter-after-tapering-p (keyword; t) ==> if t and data-taper is a function, centers tapered series by subtracting off its sample mean [10] restore-power-option-p (keyword; t) ==> if t and data-taper is a function, normalizes tapered series to have same sum of squares as before tapering [11] result (keyword; vector of size (- end start)) <== sequence to hold centered, prewhitened and tapered time series returns [1] result, the sequence containing the centered, prewhitened and tapered time series (time-series are unaltered unless it is bound to result) [2] the number used to center the time series (this is nil if center-data is nil) [3] C_h, the variance inflation factor due to the data taper --- Note: see also center&taper-time-series in tapers.lisp

CENTER&TAPER-TIME-SERIES (TIME-SERIES &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (DATA-TAPER NIL) (DATA-TAPER-PARAMETERS) (RECENTER-AFTER-TAPERING-P T) (RESTORE-POWER-OPTION-P T) (RESULT (MAKE-ARRAY (- END START))))

given [1] time-series (required) ==> a sequence of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [3] start (keyword; 0) ==> start index of sequence to be used [4] end (keyword; length of time-series) ==> 1 + end index of sequence to be used [5] data-taper (keyword; nil) ==> nil or a tapering function [6] data-taper-parameters (keyword) ==> parameters for tapering function (not used if data-taper is nil) [7] recenter-after-tapering-p (keyword; t) ==> if t and data-taper is a function, centers tapered series by subtracting off its sample mean [8] restore-power-option-p (keyword; t) ==> if t and data-taper is a function, normalizes tapered series to have same sum of squares as before tapering [9] result (keyword; vector of size (- end start)) <== sequence to hold centered time series returns [1] result, the sequence containing the centered and/or tapered time series (the contents of time-series are unaltered unless it is also bound to result) [2] the number used to center the time series (this is nil if center-data is nil) [3] C_h, the variance inflation factor due to the data taper --- Note: see also center&prewhiten&taper-time-series in filtering.lisp

CEPSTRUM->I_M (CEPSTRUM &KEY (SAMPLING-TIME 1.0) (N-PRIME (- (* 2 (LENGTH CEPSTRUM)) 2)) (LAG-WINDOW-FUNCTION #'PARZEN-LAG-WINDOW) (LAMBDA-FACTOR 1.0) (RESULT (MAKE-ARRAY (LENGTH CEPSTRUM))))

given [1] cepstrum (required) ==> a vector of length N_U+1 with values of the cepstrum from lag 0 to N_U [2] sampling-time (keyword; 1.0) ==> the sample time (i.e., delta t) [3] N-prime (keyword; 2 * N_U) ==> effective sample size, as discussed in Section 6.15 of the SAPA book (default value ok if N-prime is even, but must be replaced if N-prime is odd) [4] lag-window-function (keyword; #'parzen-lag-window) ==> lag window function of 2 parameters, tau (the lag) and m (the window parameter) [5] lambda-factor (keyword; 1.0) ==> lambda, as discussed in Section 6.15 of the SAPA book [6] result (keyword; vector of same length as cepstrum) <== vector to hold I_m returns [1] the value m at which I_m is minimized [2] result, i.e., I_m for m = 0 to N_U [3] the minimum value of I_m for m = 0 to N_U --- Note: see Equation (283b) of the SAPA book

CHECK-ORTHONORMALITY (LIST-OF-VECTORS)

given a list of vectors, computes and prints their sum of squares and pairwise dot products, and returns the maximum absolute deviation from 0 of the pairwise dot products

CHOLESKY! (A B &KEY (EPS SINGLE-FLOAT-EPSILON))

given [1] A (required) ==> vector of length N*(N+1)/2 representing a two dimensional N by N positive definite matrix with elements stored columnwise; i.e, elements with indices 0 1 2 3 4 5 ... correspond to matrix positions 1,1 1,2 2,2 1,3 2,3 3,3 ... [2] b (required) <=> right-hand vector on input; solution X on output [3] eps (keyword; single-float-epsilon) number used to test for a computationally singular matrix solves A X = b using the Cholesky decomposition method, and returns [1] X, the solution vector (stored in b)

CIRCULAR-SHIFT-SEQUENCE (A-SEQ &KEY (RESULT (MAKE-ARRAY (LENGTH A-SEQ))))

given [1] a-seq (required) ==> a sequence [2] result (keyword; vector of same length as a-seq) <== a sequence returns [1] a `circularly' left-shifted sequence; i.e., maps x(0), x(1), x(2), ..., x(n-2), x(n-1) to x(n-1), x(0), x(1), ..., x(n-3), x(n-2) ---- Note: result can be the same sequence as result

COMPARE-SEQS (ONE-SEQ ANOTHER-SEQ)

given [1] one-seq (required) ==> any sequence [2] another-seq (required) ==> any sequence (must be same length as one-seq) returns [1] maximum absolute difference between corresponding elements of the two sequences [2] average absolute difference --- Note: useful for comparing results that in theory should be identical

COMPOSE-SYMMETRIC-FILTERS (&REST FILTERS)

given any number of symmetric filters (each of odd length and each represented by a vector of coefficients), returns the composite filter that will produce an output identical to that obtained by cascading the individual filters

CONVERT-FROM-DB (NUM)

given [1] num (required) ==> a number returns [1] 10^(x/10)

CONVERT-TO-DB (NUM)

given [1] num (required) ==> a number returns [1] 10 log_10(num) --- Note: num is intended to be a positive number, but this is not checked

COPY-VECTOR (FROM-VECTOR TO-VECTOR &KEY (START 0) (END (LENGTH FROM-VECTOR)))

given [1] from-vector (required) ==> a vector (actually, works with any sequence) [2] to-vector (required) ==> a vector (actually, works with any sequence) [3] start (keyword; 0) ==> start index of from-vector to be transferred [4] end (keyword; length of from-vector) ==> 1 + end index of from-vector to be transferred returns [1] to-vector, with elements start to end - 1 from-vector copied into corresponding elements 0 to end - start - 1 of to-vector --- Note: if to-vector is being created on the spot, might consider using CL's copy-seq instead.

COSINE-DATA-TAPER! (A-TIME-SERIES &KEY (TAPER-PARAMETER 1.0) (NORMALIZATION UNITY) (RESULT A-TIME-SERIES))

given [1] a-time-series (required) <=> a vector containing a time series; this vector is modified UNLESS keyword result is bound to a different vector [2] taper-parameter (keyword; 1.0, i.e., Hanning data taper) ==> a number >= 0.0 and <= 1.0 that specifies the degree of tapering (parameter p of Equation (209)). [3] normalization (keyword; :unity) ==> if :unity, taper normalized such that its sum of squares is equal to unity (Equation (208a)); if :N, taper normalized such that its sum of squares is equal to the number of points in the time series [4] result (keyword; a-time-series) <=> a vector to contain time series multiplied by cosine data taper returns [1] a vector containing the tapered time series [2] C_h, the variance inflation factor computed using Equation (251b) in the SAPA book

CREATE-CI-FOR-AMT-SDF-ESTIMATE (SDF-DB DOFS &KEY (CONFIDENCE-LEVEL 0.95))

given [1] sdf-dB (required) ==> vector containing an adaptive multitaper spectral estimate expressed in decibels [2] dofs (required) ==> vector containing the degrees of freedom associated with the values in sdf-dB [3] confidence-level (keyword; 0.95) ==> the level of the confidence intervals to be created returns [1] a vector containing the upper confidence interval [2] a vector containing the lower confidence interval --- Note: see Section 7.4 of the SAPA book

CREATE-DPSS-LOW-PASS-FILTER (FILTER-LENGTH DELTA W &KEY (NYQUIST-FREQUENCY 0.5) (RESULT (MAKE-ARRAY FILTER-LENGTH)))

given [1] filter-length (required) ==> an odd positive integer = 2L +1 [2] delta (required) ==> ``W'' parameter for dpss in user units (see page 182 of the SAPA book) [2] W (required) ==> a cutoff frequency greater than 0 and less than the Nyquist frequency [3] Nyquist-frequency (keyword; 0.5) ==> the Nyquist frequency [4] result (keyword; vector of length filter-length) <== vector of length filter-length into which filter coefficients are placed (returned by the function) uses a dpss as convergence factors in least squares approximation to a low-pass filter and returns [1] a symmetric low-pass filter of length 2L+1 and a gain of unity at zero frequency; i.e., result(0) = result(2L) result(1) = result(2L-1) etc., and (sum result) = 1.0 --- Note: see Section 5.9 of the SAPA book; (aref result 0) corresponds to element -L of the filter; (aref result L) corresponds to element 0 of the filter; (aref result (* 2 L)) corresponds to element L of the filter

CREATE-LEAST-SQUARES-LOW-PASS-FILTER (FILTER-LENGTH W &KEY (CONVERGENCE-FACTORS NIL) (NYQUIST-FREQUENCY 0.5) (RESULT (MAKE-ARRAY FILTER-LENGTH)))

given [1] filter-length (required) ==> an odd positive integer = 2L +1 [2] W (required) ==> a cutoff frequency greater than 0 and less than the Nyquist frequency [3] convergence-factors (keyword; nil) ==> a one-argument function that maps an integer to a convergence factor; nil is also acceptable, in which case no convergence factors are used [4] Nyquist-frequency (keyword; 0.5) ==> the Nyquist frequency [5] result (keyword; vector of length filter-length) <== vector of length filter-length into which filter coefficients are placed (returned by the function) uses a least squares approximation to a low-pass filter and returns [1] a symmetric low-pass filter of length 2L+1 and a gain of unity at zero frequency; i.e., result(0) = result(2L) result(1) = result(2L-1) etc., and (sum result) = 1.0 --- Note: see Section 5.8 of the SAPA book; (aref result 0) corresponds to element -L of the filter; (aref result L) corresponds to element 0 of the filter; (aref result (* 2 L)) corresponds to element L of the filter

CUMULATIVE-PERIODOGRAM (TIME-SERIES &KEY (SAMPLING-TIME 1.0) (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (CUMULATIVE-PERIODOGRAM-TEST-P T) (SIGNIFICANCE-LEVEL 0.95) (SCRATCH-DFT (MAKE-ARRAY (- END START))) (RETURN-FREQUENCIES-P T) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (TRUNCATE (1- (- END START)) 2)))) (RESULT (MAKE-ARRAY (TRUNCATE (1- (- END START)) 2))))

given [1] time-series (required) ==> a vector containing time series values [2] sampling-time (keyword; 1.0) ==> the sample time (i.e., delta t) [3] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [4] start (keyword; 0) ==> start index of vector to be used [5] end (keyword; length of time-series) ==> 1 + end index of vector to be used [6] cumulative-periodogram-test-p (keyword; t) ==> if t, returns values associated with Kolmogorov test statistic [7] significance-level (keyword; 0.95) ==> level of significance at which Kolmogorov test is performed [8] scratch-dft (keyword; vector of size (- end start)) ==> scratch vector used for in-place dft calculations [9] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the cumulative periodogram are computed and returned in result-freq [10] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of correct length into which the frequencies associated with the values in result are placed [11] result (keyword; vector of correct length) <== vector to hold normalized cumulative periodogram returns [1] vector with normalized cumulative periodogram [2] vector with associated frequencies and -- if cumulative-periodogram-test-p is true -- [3] Kolmogorov test statistic [4] index of maximum deviation [5] frequency of maximum deviation [6] quantile of Kolmogorov test statistic (under null hypothesis) [7] either :reject or :fail-to-reject, depending on whether or not we reject or fail to reject the null hypothesis [8] slope of upper and lower lines [9] intercept of upper line [10] intercept of lower line --- Note: see Section 6.6 of the SAPA book

CUMULATIVE-SUMS (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)) (RESULT (MAKE-ARRAY (- END START))))

given [1] the-seq (required) ==> a sequence of numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sequence with cumulative sum of elements in specified subsequence

DANIELL-BANDWIDTH->M (B_W &KEY (SAMPLING-TIME 1.0))

given desired smoothing window bandwidth B_W and sampling time, returns [1] window parameter m required to approximately achieve B_W using the Daneill lag window [2] actual B_W achieved (in fact, this is always equal to B_W, but we return it anyway for consistency with other lag windows) --- Note: see Table 269 of the SAPA book

DANIELL-LAG-WINDOW (TAU M)

given the lag tau and window parameter m, returns the value of the Daniell lag window --- Note: see equation between Equations (264a) and (264b) of the SAPA book or Priestley, page 441, Equation (6.2.73)

DANIELL-M->BANDWIDTH (M &KEY (SAMPLING-TIME 1.0))

given window parameter m and sampling time, returns bandwidth B_W for the Daniell smoothing window --- Note: see Table 269 of the SAPA book

DANIELL-N-M->DEGREES-OF-FREEDOM (N M &KEY (C_H 1.0))

given sample size N, window parameter m and variance inflation factor C_h, returns equivalent degrees of freedom nu for Daniell lag window --- Note: see Table 269 of the SAPA book

DFT! (X &KEY (N (LENGTH X)) (SAMPLING-TIME 1.0))

given: [1] x (required) <=> a vector of real or complex-valued numbers [2] N (keyword; length of x) ==> number of points in dft [3] sampling-time (keyword; 1.0) returns: [1] x, with contents replaced by the discrete Fourier transform of x, namely, N-1 X(n) = sampling-time SUM x(t) exp(-i 2 pi n t/N) t=0 --- Note: see Equation (110a) of the SAPA book

DFT-CHIRP! (COMPLEX-VECTOR &KEY (N (LENGTH COMPLEX-VECTOR)))

given: [1] complex-vector (required) <=> a vector of real or complex-valued numbers [2] N (keyword; length of complex-vector) ==> number of points (must be a power of 2) computes the discrete Fourier transform of complex-vector using a chirp transform algorithm and returns: [1] complex-vector, with contents replaced by the discrete Fourier transform of the input, namely, N-1 X(n) = SUM x(t) exp(-i 2 pi n t/N) t=0 --- Note: see Equation (110a) of the SAPA book with the sampling time set to unity and also Section 3.10 for a discussion on the chirp transform algorithm

DIFFERENCE (SEQUENCE-OF-NUMBERS &KEY (START 0) (END (LENGTH SEQUENCE-OF-NUMBERS)) (LAG 1) (RESULT (MAKE-ARRAY (- END START LAG))))

given [1] sequence-of-numbers (required) ==> a sequence of numbers [2] start (keyword; 0) ==> start index [3] end (keyword; length of sequence-of-numbers) ==> end index plus one [4] lag (keyword; 1) ==> lag for differencing [5] result (keyword; vector of length (- end start lag)) <== results return [1] a sequence of length (- end start lag) with lag-th order difference of sequence from index start to index end-1 --- Note: See Newton, 1988, page 31.

DIGAMMA (X &KEY (X-RECURSION 8.5) (X-SMALL 1.e-5))

given [1] x (required) ==> a positive number [2] x-recursion (keyword; 8.5) ==> for noninteger x, if x-small < x < x-recursion, recursive formula is used [3] x-small (keyword; 1.0e-5) ==> if x <= x-small, small x approximation used (default ) returns [1] value of digamma (psi, derivative of log(gamma)) function at x --- Note: see Abramowitz and Stegun's equation 6.3.2; expansion 6.3.18 plus recurrence 6.3.5; the small x formula is Equation (5) from Algorithm AS 103 -- Psi (Digamma) Function -- by J. M. Bernardo in Applied Statistics, Vol. 25, No. 3, 1976, pp. 315--317; note that better accuracy can be obtained by increasing x-recursion --- 10.0 to 100.0 are probably useful values

DIRECT-SPECTRAL-ESTIMATE (TIME-SERIES &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (N-NONZERO-FREQS HALF-NEXT-POWER-OF-2) (RETURN-EST-FOR-0-FREQ-P NIL) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS (- END START)))) (DATA-TAPER NIL) (DATA-TAPER-PARAMETERS) (RECENTER-AFTER-TAPERING-P T) (RESTORE-POWER-OPTION-P T) (RETURN-ACVS-P NIL) (RESULT-ACVS (MAKE-ARRAY (- END START))) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P)))))

given [1] time-series (required) ==> a vector of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of time-series) ==> 1 + end index of vector to be used [5] N-nonzero-freqs (keyword; :half-next-power-of-2) ==> specifies at how many nonzero frequencies direct spectral estimate is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [6] return-est-for-0-freq-p (keyword; nil) ==> if t, sdf is computed at zero frequency; otherwise, it is not computed. [7] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [8] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [9] data-taper (keyword; nil) ==> nil or a tapering function [10] data-taper-parameters (keyword) ==> parameters for tapering function (not used if data-taper is nil) [11] recenter-after-tapering-p (keyword; t) ==> if t and data-taper is a function, centers tapered series by subtracting off its sample mean [12] restore-power-option-p (keyword; t) ==> if t and data-taper is a function, normalizes tapered series to have same sum of squares as before tapering [13] return-acvs-p (keyword; nil) ==> if t, computes acvs corresponding to direct spectral estimate [14] result-acvs (keyword; vector of correct length) <== vector into which acvs values are placed (not used if return-acvs-estimate-p is nil) [15] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [16] result-sdf (keyword; vector of correct length) <== vector into which direct spectral estimates are placed; it must be exactly of the length dictated by N-nonzero-freqs and return-est-for-0-freq-p [17] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral estimate are computed and returned in result-freq [18] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [19] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length dictated by N-nonzero-freqs and return-est-for-0-freq-p into which the frequencies associated with the values in result-sdf are placed returns [1] result-sdf, a vector holding the properly transformed sdf [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-sdf -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-sdf [4] C_h, the variance inflation factor due to the data taper [5] result-acvs (if return-acvs-p is t), a vector holding the acvs corresponding to the direct spectral estimate -- or -- nil (if return-acvs-p is nil) --- Note: see Section 6.3 of the SAPA book

DOT-PRODUCT (X Y)

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers returns [1] the dot product of x and y ---- Note: corresponds to SDOT in LINPACK's bla

DPSS-DATA-TAPER! (A-TIME-SERIES &KEY (TAPER-PARAMETER 4.0) (NORMALIZATION UNITY) (RESULT A-TIME-SERIES))

given [1] a-time-series (required) <=> a vector containing a time series; this vector is modified UNLESS keyword result is bound to a different vector [2] taper-parameter (keyword; 4.0) ==> a number > 0.0 that specifies NW, the product of the sample size and the half-width of the concentration interval (see SAPA, page 211) [3] normalization (keyword; :unity) ==> if :unity, taper normalized such that its sum of squares is equal to unity (Equation (208a)); if :N, taper normalized such that its sum of squares is equal to the number of points in the time series [4] result (keyword; a-time-series) <=> a vector to contain time series multiplied by cosine data taper returns [1] a vector containing the tapered time series [2] C_h, the variance inflation factor computed using Equation (251b) in the SAPA book

DPSS-TAPERS-INVERSE-ITERATION (N NUMBER-OF-TAPERS &KEY (TAPER-PARAMETER 4.0) (PRINT-PROGRESS-P NIL) (EPS 5.e-5))

given [1] N (required) the sample size [2] number-of-tapers (required) number of orthonormal dpss data tapers to be computed [3] taper-parameter (keyword; 4.0) NW, the duration-half-bandwidth product (must be such that 0 < NW/N < 1/2) [4] print-progress-p (keyword; nil) if t, prints a dot after each eigenvalue and eigenvector has been computed [5] eps (keyword; 0.5e-4) controls accuracy returns [1] a list of length number-of-tapers of N dimensional vectors with orthonormal dpss's of orders 0, 1, ..., number-of-tapers - 1; [2] a vector of length number-of-tapers with the corresponding eigenvalues --- Note: computes the dpss tapers using inverse iteration (see Section 8.1 of the SAPA book)

DPSS-TAPERS-THOMSON-APPROX (N NUMBER-OF-TAPERS &KEY (TAPER-PARAMETER 4.0) (PRINT-PROGRESS-P NIL) (COMPUTE-TRUE-EIGENVALUES-P NIL) (ABSCISSAS (LET () (DECLARE (SPECIAL *ABSCISSAS-32-POINT*)) *ABSCISSAS-32-POINT*)) (WEIGHTS (LET () (DECLARE (SPECIAL *WEIGHTS-32-POINT*)) *WEIGHTS-32-POINT*)))

given [1] N (required) the sample size [2] number-of-tapers (required) number of orthonormal dpss data tapers to be computed [3] taper-parameter (keyword; 4.0) NW, the duration-half-bandwidth product (must be such that 0 < NW/N < 1/2) [4] print-progress-p (keyword; nil) if t, prints a dot after each eigenvalue and eigenvector has been computed [5] compute-true-eigenvalues-p (keyword; nil) if t, returns eigenvalues for eigenproblem of Equation (378) of the SAPA book; if nil, returns nil [6] abscissas (keyword; *abscissas-32-point*) a vector of abscissas points used in Gauss-Legendre quadrature [7] weights (keyword; *abscissas-32-point*) a vector of weights used in Gauss-Legendre quadrature returns [1] a list of length number-of-tapers of N dimensional vectors with orthonormal dpss's of orders 0, 1, ..., number-of-tapers - 1; [2] a vector of length number-of-tapers with eigenvalues if compute-true-eigenvalues-p is t; nil if if compute-true-eigenvalues-p is nil --- Note: computes the dpss tapers using Thomson's numerical integration scheme (see Section 8.2 of the SAPA book)

DPSS-TAPERS-TRI-DIAG (N NUMBER-OF-TAPERS &KEY (TAPER-PARAMETER 4.0) (PRINT-PROGRESS-P NIL) (COMPUTE-TRUE-EIGENVALUES-P NIL))

given [1] N (required) the sample size [2] number-of-tapers (required) number of orthonormal dpss data tapers to be computed [3] taper-parameter (keyword; 4.0) NW, the duration-half-bandwidth product (must be such that 0 < NW/N < 1/2) [4] print-progress-p (keyword; nil) if t, prints a dot after each eigenvalue and eigenvector has been computed [5] compute-true-eigenvalues-p (keyword; nil) if t, returns eigenvalues for eigenproblem of Equation (378) of the SAPA book; if nil, returns eigenvalues for tridiagonal formulation returns [1] a list of length number-of-tapers of N dimensional vectors with orthonormal dpss's of orders 0, 1, ..., number-of-tapers - 1; [2] a vector of length number-of-tapers with eigenvalues as specified by compute-true-eigenvalues-p --- Note: computes the dpss tapers using the tridiagonal formulation (see Section 8.3 of the SAPA book)

DURBIN-WATSON-TEST-STATISTIC (RESIDUALS)

given a set of residuals in a seqeunce, returns the Durbin-Watson test statistic

EIGENSPECTRA->ADAPTIVE-MULTITAPER-SPECTRAL-ESTIMATE (LIST-OF-EIGENSPECTRA EIGENVALUES VARIANCE-OF-TIME-SERIES &KEY (SAMPLING-TIME 1.0) (N-EIGENSPECTRA (LENGTH LIST-OF-EIGENSPECTRA)) (MAXIMUM-NUMBER-OF-ITERATIONS 100) (RESULT-DOF (MAKE-ARRAY (LENGTH (CAR LIST-OF-EIGENSPECTRA)))) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (LENGTH (CAR LIST-OF-EIGENSPECTRA)))))

given [1] list-of-eigenspectra (required) ==> list of eigenspectra (such as optionally returned by multitaper-spectral-estimate); each eigenspectrum is assumed to be untransformed (e.g., not expressed in dB) [2] eigenvalues (required) ==> vector of eigenvalues corresponding to the dpss's used to create the eigenspectra (the length of eigenvalues should be at least as large as the length specified by N-eigenspectra) [3] variance-of-time-series (required) ==> variance of time series from which eigenspectra were computed [4] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [5] N-eigenspectra (keyword; length of list-of-eigenspectra) ==> number of eigenspectra to be used (must be less than or equal to the length of list-of-eigenspectra) [6] maximum-number-of-iterations (keyword; 100) ==> maximum number of iterations [7] result-dof (keyword; vector of correct length) <== vector into which degrees of freedom for each value in the adaptive multitaper spectral estimate is placed; it must be exactly the same length as each of the vectors in list-of-eigenspectra [8] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [9] result-sdf (keyword; vector of correct length) <== vector into which adaptive multitaper spectral estimate is placed; it must be exactly the same length as each of the vectors in list-of-eigenspectra returns [1] result-sdf, a vector holding the properly transformed adaptive multitaper spectral estimate [2] result-dof, a vector holding corresponding degrees of freedom [3] the maximum number of iterations required to reach convergence for all of the values in result-sdf --- Note: see Section 7.4 of the SAPA book

EIGENSPECTRA->MULTITAPER-SPECTRAL-ESTIMATE (LIST-OF-EIGENSPECTRA &KEY (N-EIGENSPECTRA (LENGTH LIST-OF-EIGENSPECTRA)) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (LENGTH (CAR LIST-OF-EIGENSPECTRA)))))

given [1] list-of-eigenspectra (required) ==> list of eigenspectra (such as optionally returned by multitaper-spectral-estimate); each eigenspectrum is assumed to be untransformed (e.g., not expressed in dB) and represented by a vector [2] N-eigenspectra (keyword; length of list-of-eigenspectra) ==> number of eigenspectra to be used (must be less than or equal to the length of list-of-eigenspectra) [3] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [4] result-sdf (keyword; vector of correct length) <== vector into which multitaper spectral estimate is placed; it must be exactly the same length as each of the vectors in list-of-eigenspectra returns [1] result-sdf, a vector holding the properly transformed multitaper spectral estimate --- Note: see Section 7.1 of the SAPA book

EQUIVALENT-DEGREES-OF-FREEDOM (LAG-WINDOW-FUNCTION MAX-LAG &KEY (SAMPLING-TIME 1.0) (SAMPLE-SIZE (1+ MAX-LAG)) (C_H 1.0))

given a lag window function, a maximum lag, the sampling time, the sample size and the variance inflation factor C_h, returns the corresponding equivalent degrees of freedom --- Note: see Equation (255a) of the SAPA book

EUCLIDEAN-NORM (X)

given [1] x (required) ==> a sequence of numbers returns [1] Euclidean norm of x ---- Note: corresponds to LINPACK's SNRM2, but this is POOR implementation!

EVALUATE-POLYNOMIAL (POLYNOMIAL Z &KEY (DEGREE-OF-POLYNOMIAL (1- (LENGTH POLYNOMIAL))) (NUMBER-OF-DERIVATIVES 0) (COMPLEX-VALUES (MAKE-ARRAY (1+ NUMBER-OF-DERIVATIVES))) (BOUNDS-TO-BE-COMPUTED NIL) (BOUNDS (IF BOUNDS-TO-BE-COMPUTED (MAKE-ARRAY (1+ NUMBER-OF-DERIVATIVES)))))

given [1] polynomial (required) ==> sequence of length n+1 with real or complex-valued numbers giving the n+1 coefficients of a polynomial of nth order; (elt polynomial i) = coefficient of z^(n-i) [2] z (required) ==> complex point at which the polynomial to be evaluated [3] degree-of-polynomial (keyword; (1- (length polynomial))) ==> degree of polynomial [4] number-of-derivatives (keyword; 0) ==> number of derivatives to be evaluated [5] complex-values (keyword; array of length (1+ number-of-derivatives)) <== sequence with values of polynomial and derivatives [6] bounds-to-be-computed (keyword; nil) ==> if t, error bounds are computed for values in complex-values [7] bounds (keyword; array of length (1+ number-of-derivatives)) <== sequence with error bounds returns [1] the sequence complex-values containing the value of the polynomial and its derivatives [2] bounds, a sequence containing optional error bounds --- Note: this is a Lisp version of cmlib routine cpevl; d1 was specified in cpevl as (expt 2 (- 1 (i1mach 11))), where (i1mach 11) ==> THE NUMBER OF BASE-2 DIGITS (SINGLE PRECISION). I have taken this to be equivalent to machine epsilon. If this is in fact NOT the case, then the optional error bounds might not be computed correctly.

EXPONENTIAL-SMOOTHING (TIME-SERIES ALPHA &KEY (INITIAL-PREDICTION (ELT TIME-SERIES 0)) (RETURN-SMOOTHED-VALUES-P T) (N-PREDICTIONS 0) (RESULT-SMOOTHED (IF RETURN-SMOOTHED-VALUES-P (MAKE-ARRAY (LENGTH TIME-SERIES)))) (RESULT-PREDICTIONS (IF (PLUSP N-PREDICTIONS) (MAKE-ARRAY N-PREDICTIONS))))

given [1] time-series (required) ==> a sequence of numbers [2] alpha (required) ==> parameter controlling the degree of exponential smoothing (usually 0 < alpha < 1) [3] initial-prediction (keyword; first element of time-series) ==> to get the smoother going, we need a ``prediction'' for the first element in the sequence time-series -- two common hacks are the first element itself (the default) or 0 (if the time series has a zero sample mean) [4] return-smoothed-values-p (keyword; t) ==> if t, returns smoothed (i.e., filtered) time series [5] n-predictions (keyword; 0) ==> number of predictions desired [6] result-smoothed (keyword; vector of length of time-series or nil) <== a sequence to hold smoothed (filtered) time series -- used only if return-smoothed-values-p is true [7] result-predictions (keyword; vector of length n-predictions or nil) <== a sequence to hold predicted values of the time series -- used only if n-predictions is positive returns [1] sum-of-squares of prediction errors [2] either one step ahead prediction (if n-predictions is 0) or vector of length n-predictions if n-predictions > 0) [3] vector with smoothed values (nil if return-smoothed-values-p is true) --- Note: see Section 7.3 of ``Time Series: A Biostatistical Introduction'' by Diggle, 1990

FACTORIAL (K)

given an integer k, returns k!

FAST-FORWARD-BACKWARD-LS (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)) (PSEUDO-REFLECTION-COEFFS (MAKE-ARRAY P)))

given [1] time-series (required) ==> a vector containing a real-valued or complex-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} [7] pseudo-reflection-coeffs (keyword; a vector of length p) <== sequence of reflection coefficients computed stepwise by the algorithm, but these do not correspond to the estimated phi_{k,p} (even if the estimated AR model is stationary) uses the forward/backward least squares method to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's [3] pseudo-reflection-coeffs --- Note: see Section 9.7 of the SAPA book; this function is an adaptation of the Fortran routine modcovar, pp. 258-60, ``Digital Spectral Analysis with Applications'' by Marple, 1987

FFT! (COMPLEX-VECTOR &KEY (N (LENGTH COMPLEX-VECTOR)))

given: [1] complex-vector (required) <=> a vector of real or complex-valued numbers [2] N (keyword; length of complex-vector) ==> number of points (must be a power of 2) computes the discrete Fourier transform of complex-vector using a fast Fourier transform algorithm and returns: [1] complex-vector, with contents replaced by the discrete Fourier transform of the input, namely, N-1 X(n) = SUM x(t) exp(-i 2 pi n t/N) t=0 --- Note: see Equation (110a) of the SAPA book with the sampling time set to unity

FILTER-TIME-SERIES (TIME-SERIES THE-FILTER &KEY (START 0) (END (LENGTH TIME-SERIES)) (TECHNIQUE FASTEST) (RESULT (MAKE-ARRAY (- (- END START) (1- (LENGTH THE-FILTER))))))

given [1] time-series (required) ==> a vector containing a time series x_0, x_1, ..., x_{N-1} [2] the-filter (required) ==> a vector containing the filter coefficients g_0, g_1, ..., x_{K-1} [3] start (keyword; 0) ==> start index of time-series to be used [4] end (keyword; length of time-series) ==> 1 + end index of time-series to be used [5] technique (keyword; :fastest) ==> if :fastest, tries to pick fastest method; if :fft, uses fft-based method; if :direct, uses direct method [6] result (keyword; vector of appropriate length) <== vector to contain filtered time series K-1 y_t = SUM g_k x_{t+K-1-k}, t = 0, ..., N-K+1 k=0 returns [1] result, a vector containing the filtered time series [2] the number of values in the filtered time series --- Note: result can be the same as time-series

FILTER-TIME-SERIES-DIRECT (TIME-SERIES THE-FILTER &KEY (START 0) (END (LENGTH TIME-SERIES)) (RESULT (MAKE-ARRAY (- (- END START) (1- (LENGTH THE-FILTER))))))

given [1] time-series (required) ==> a vector containing a time series x_0, x_1, ..., x_{N-1} [2] the-filter (required) ==> a vector containing the filter coefficients g_0, g_1, ..., x_{K-1} [3] start (keyword; 0) ==> start index of time-series to be used [4] end (keyword; length of time-series) ==> 1 + end index of time-series to be used [5] result (keyword; vector of appropriate length) <== vector to contain filtered time series K-1 y_t = SUM g_k x_{t+K-1-k}, t = 0, ..., N-K+1 k=0 returns [1] result, a vector containing the filtered time series [2] the number of values in the filtered time series --- Note: result can be the same as time-series

FILTER-TIME-SERIES-FFT (TIME-SERIES THE-FILTER &KEY (START 0) (END (LENGTH TIME-SERIES)) (FFT-SIZE (* 4 (NEXT-POWER-OF-2 (LENGTH THE-FILTER)))) (VERBOSE-P NIL) (RESULT (MAKE-ARRAY (- (- END START) (1- (LENGTH THE-FILTER))) INITIAL-ELEMENT 0.0) RESULT-SUPPLIED-P))

given [1] time-series (required) ==> a vector containing a time series x_0, x_1, ..., x_{N-1} [2] the-filter (required) ==> a vector containing the filter coefficients g_0, g_1, ..., x_{K-1} [3] start (keyword; 0) ==> start index of time-series to be used [4] end (keyword; length of time-series) ==> 1 + end index of time-series to be used [5] fft-size (keyword; 4 * power of 2 that is ceiling for filter length) ==> size of fft's to be used (must be a power of 2) [6] verbose-p (keyword; nil) ==> if t, prints line after each block of data is processed; if nil, prints nothing [7] result (keyword; vector of appropriate length) <== vector to contain filtered time series K-1 y_t = SUM g_k x_{t+K-1-k}, t = 0, ..., N-K+1 k=0 returns [1] result, a vector containing the filtered time series [2] the number of values in the filtered time series --- Note: result can be the same as time-series

FIND-PEAK-OF-AR-SDF-USING-QUADRATIC-APPROX (F_1 F_2 F_3 AR-COEFFS &KEY (SAMPLING-TIME 1.0) (DB-TOLERANCE 0.1))

given [1] f_1 (required) ==> frequency f_1 such that f_1 < f_2 and S(f_1) < S(f_2), where S(.) is the sdf associated with AR-coeffs [2] f_2 (required) ==> middle frequency [3] f_3 (required) ==> frequency f_3 such that f_3 > f_2 and S(f_3) < S(f_2) [4] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [5] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [6] dB-tolerance (keyword; 0.1) ==> convergence criterion returns [1] estimate of frequency of peak/valley in sdf --- Note: for details on this method, see pages 20-21 of Adby and Dempster, ``Introduction to Optimization Methods''; note in particular that we search the inverse of the AR portion of the spectrum instead of the spectrum itself in the belief that the quadratic approximation is better with this reformulation

FIND-PEAK-OR-VALLEY-OF-SDF-USING-BISECTION+NEWTON-RAPHSON (F-LEFT F-RIGHT ERSATZ-ACVS &KEY (SAMPLING-TIME 1.0) (N (LENGTH ERSATZ-ACVS)) (ACCURACY (* 10.0 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 20))

given [1] f-left (required) ==> left bracket for frequency of peak/valley [2] f-right (required) ==> right bracket for frequency of peak/valley [3] ersatz-acvs (required) ==> vector of length N with acvs (or acs) s_0 to s_{N-1} corresponding to sdf (only elements 1 to N-1 are used) [4] sampling-time (keyword; 1.0) ==> positive number [5] N (keyword; length of ersatz-acvs) ==> positive integer [6] accuracy (keyword; 10.0 * single-float-epsilon) ==> passed to bisection-with-Newton-Raphson [7] maximum-number-of-iterations (keyword; 20) ==> passed to bisection-with-Newton-Raphson returns [1] estimate of frequency of peak/valley in sdf [2] number of iterations required --- Note: searches for peak using a combination of bisection and Newton-Raphson; see pages 479--80 and 524--5 of the SAPA book

FISHER-G-STATISTIC (TIME-SERIES &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (ALPHA 0.05))

given [1] time-series (required) ==> a vector of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of time-series) ==> 1 + end index of vector to be used [5] alpha (keyword; 0.05) ==> critical level at which Fisher's g test is performed returns [1] Fisher's g statistic [2] approximation to critical level of Fisher's g test from Equation (491c) [3] either :reject or :fail-to-reject, depending on whether or not we reject or fail to reject the null hypothesis of white noise --- Note: see Section 10.9 of the SAPA book

FORWARD-BACKWARD-LS (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)))

given [1] time-series (required) ==> a vector containing a real-valued or complex-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} uses the forward/backward least squares method to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's --- Note: see Section 9.7 of the SAPA book; this function is an adaptation of the Fortran routine covmcov, pp. 262-4, ``Modern Spectrum Estimation: Theory and Application'' by Kay, 1988

FORWARD-LS (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)))

given [1] time-series (required) ==> a vector containing a real-valued or complex-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} uses the forward least squares method to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's --- Note: see Section 9.7 of the SAPA book; this function is an adaptation of the Fortran routine covmcov, pp. 262-4, ``Modern Spectrum Estimation: Theory and Application'' by Kay, 1988

GAUSS-LEGENDRE-QUADRATURE (LOWER-LIMIT UPPER-LIMIT N)

given [1] lower-limit (required) ==> lower limit of integration [2] upper-limit (required) ==> upper limit of integration [3] N (required) ==> number of points to be computed in Gauss-Legendre quadrature returns [1] an N-dimensional vector of abscissas points [2] an N-dimensional vector of weights --- Note: this function is based on gauleg, Section 4.5, Numerical Recipes, Second Edition

GENERATE-AR-TIME-SERIES (COEFFS VARIANCE SAMPLE-SIZE &KEY (PROCESS-VARIANCE-P T) (LIST-OF-LOWER-ORDER-PHI NIL) (LIST-OF-LOWER-ORDER-PEV NIL) (RESULT (MAKE-ARRAY SAMPLE-SIZE INITIAL-ELEMENT 0.0) RESULT-SUPPLIED-P))

given [1] coeffs (required) ==> sequence of length p with AR coefficients: x_t = coeffs_0*x_{t-1} + coeffs_1*x_{t-2} + ... + coeffs_{p-1}*x_{t-p} + e_t (see Equation (392a) in the SAPA book; the coefficients can be real or complex-valued) [2] variance (required) ==> process variance or innovations variance (see keyword process-variance-p) [3] sample-size (required) ==> length of generated time series [4] process-variance-p (keyword; t) ==> if t, variance is taken to be process variance if nil, variance is taken to be innovations variance [5] list-of-lower-order-phi (keyword; nil) ==> to bypass call to step-down-Levinson-Durbin-recursions (useful for multiple realizations) [6] list-of-lower-order-pev (keyword; nil) ==> to bypass call to step-down-Levinson-Durbin-recursions [7] result (keyword; vector of length sample-size) <== vector with simulated series generates realization of zero mean normally distributed AR(p) process and returns [1] result, vector of length sample-size with realization [2] list-of-lower-order-phi (can be used on subsequent calls to this function to bypass call to step-down-Levinson-Durbin-recursions) [3] list-of-lower-order-pev (can be also used on subsequent calls) --- Note: this function generates the proper stationary initial conditions

GENERATE-BACKWARD-INNOVATIONS (TIME-SERIES AR-COEFFS &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (RESULT (MAKE-ARRAY (- (- END START) (LENGTH AR-COEFFS)))))

given [1] time-series (required) ==> a vector containing a real-valued time series x_t [2] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series (the vector time-series is not altered) [6] result (keyword; vector of appropriate length) <== vector to contain backward innovations p b_t = x_t - SUM phi_{j,p} x_{t+j}, t = 1, ..., N-p j=1 returns [1] result, i.e., the backward innovations

GENERATE-FORWARD-INNOVATIONS (TIME-SERIES AR-COEFFS &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (RESULT (MAKE-ARRAY (- (- END START) (LENGTH AR-COEFFS)))))

given [1] time-series (required) ==> a vector containing a real-valued time series x_t [2] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series (the vector time-series is not altered) [6] result (keyword; vector of appropriate length) <== vector to contain forward innovations p f_t = x_{t+p} - SUM phi_{j,p} x_{t-j}, t = 1, ..., N - p j=1 returns [1] result, i.e., the fprward innovations

GENERATE-MA-TIME-SERIES (COEFFS VARIANCE SAMPLE-SIZE &KEY (PROCESS-VARIANCE-P T) (RESULT (MAKE-ARRAY SAMPLE-SIZE)))

given [1] coeffs (required) ==> sequence of length q with MA coefficients: x_t = e_t - coeffs_0*e_{t-1} - coeffs_1*e_{t-2} - ... - coeffs_{q-1}*e_{t-q} (see Equation (43a) in the SAPA book) [2] variance (required) ==> process variance or innovations variance (see keyword process-variance-p) [3] sample-size (required) ==> length of generated time series [4] process-variance-p (keyword; t) ==> if t, variance is taken to be process variance if nil, variance is taken to be innovations variance [5] result (keyword; vector of length sample-size) <== vector with simulated series generates realization of zero mean normally distributed MA(q) process and returns [1] vector of length sample-size with realization

GENERATE-WHITE-NOISE (N &KEY (DISTRIBUTION NORMAL) (RESULT (MAKE-ARRAY N)))

given [1] n (required) ==> a sample size [2] distribution (keyword; :normal) ==> either a keyword or a function with no arguments that returns a random deviate from a particular distribution. Choices are :binary :cauchy :chi-square-2 :double-exponential :exponential :extreme-value :Gaussian (same as :normal) :logistic :lognormal :normal :uniform [3] result (keyword; vector of length n) <== a sequence in which results are to be stored return [1] result, containing n samples from a white noise process with a distribution specified by the keyword distribution

GRENANDER-SMOOTHING-WINDOW-BANDWIDTH (LAG-WINDOW-FUNCTION MAX-LAG &KEY (SAMPLING-TIME 1.0))

given a lag window function, a maximum lag and the sampling time, returns Grenander's measure of smoothing window bandwidth --- Note: see Equations (241c) and (241b) of the SAPA book; unfortunately, because of the square root operation, this measure can be complex-valued for certain lag windows

HANNING-DATA-TAPER! (A-TIME-SERIES &KEY TAPER-PARAMETER (NORMALIZATION UNITY) (RESULT A-TIME-SERIES))

calls cosine-data-taper! with taper-parameter set to 1.0

HERMITIAN-TRANSPOSE (A-MATRIX &KEY (RESULT (MAKE-ARRAY (REVERSE (ARRAY-DIMENSIONS A-MATRIX)))))

given [1] A (required) ==> a 2d matrix [2] result (keyword; new 2d array of appropriate size) <== a 2d matrix to contain Hermitian transpose of a-matrix returns [1] Hermitian transpose of a-matrix (placed in result)

HISTOGRAM (SEQUENCE-OF-DEVIATES &KEY (NUMBER-OF-BINS 10) (LEFT-OF-FIRST-BIN (MIN-OF-SEQ SEQUENCE-OF-DEVIATES)) (RIGHT-OF-LAST-BIN (MAX-OF-SEQ SEQUENCE-OF-DEVIATES)) (SCALE-AS-DENSITY-P T) (RESULT (MAKE-ARRAY NUMBER-OF-BINS)))

given [1] sequence-of-deviates (required) ==> a sequence of real-valued numbers [2] number-of-bins (keyword; 10) ==> number of bins in histogram [3] left-of-first-bin (keyword; min of sequence-of-deviates) ==> leftmost value of first bin [4] right-of-last-bin (keyword; max of sequence-of-deviates) ==> rightmost value of last bin [5] scale-as-density-p (keyword; t) ==> if t, histogram is scaled as a density; otherwise, cell counts are returned [6] result (keyword; vector of size number-of-bins) <== vector to hold the histogram returns [1] the vector result, which contains the values for the histogram for sequence-of-deviates (one value for each bin); note that, if a bin has boundaries a and b, all points x in sequence-of-deviates for which a<=x<b are counted as being in that bin EXCEPT for the last bin, where the rule is a<=x<=b. [2] maximum value in the histogram (i.e., result) [3] leftmost value of first bin [4] bin-width of bins [5] the number of unbinned elements in sequence-of-deviates less than left-of-first-bin [6] the number of unbinned elements in sequence-of-deviates greater than right-of-last-bin --- Note: the value in (svref result 0) is the histogram value for the interval [leftmost value of first bin, leftmost value of first bin + bin-width of bins); the value in (svref result 1) is the histogram value for the interval [leftmost value of first bin + bin-width of bins, leftmost value of first bin + 2*bin-width of bins); etc.

HUBER-WEIGHT-FUNCTION (X &KEY (PARAMETER 1.5))

given [1] x (required) ==> value at which to evaluate Huber's weight function [2] parameter (keyword; 1.5) ==> tuning parameter returns [1] value of Huber's weight function at x

IDEAL-BAND-PASS-FILTER-IRS (K W-LOW W-HIGH)

given [1] k (required) ==> index of member of impulse response sequence (irs) to be calculated (must be an integer) [2] W-low (required) ==> the low frequency cutoff (in standardized units) [3] W-high (required) ==> the high frequency cutoff (in standardized units so that 0 <= W-low < W-high <= 0.5, the assumed Nyquist frequency). returns [1] k-th member of the impulse response sequence for an ideal band-pass filter

IDEAL-HIGH-PASS-FILTER-IRS (K W)

given [1] k (required) ==> index of member of impulse response sequence (irs) to be calculated (must be an integer) [2] W (required) ==> the cutoff frequency, standardized such that 0 < W < 0.5 = Nyquist frequency returns [1] kth member of the impulse response sequence for an ideal low-pass filter with cutoff frequency W

IDEAL-LOW-PASS-FILTER-IRS (K W)

given [1] k (required) ==> index of member of impulse response sequence (irs) to be calculated (must be an integer) [2] W (required) ==> the cutoff frequency, standardized such that 0 < W < 0.5 = Nyquist frequency returns [1] kth member of the impulse response sequence for an ideal low-pass filter with cutoff frequency W --- Note: see Section 5.8 of the SAPA book

INTEGRATE-SDF (SDF FREQS &KEY (FIDDLE-FACTOR 2.0) (RETURN-INTEGRATED-SDF-P NIL) (RESULT (WHEN RETURN-INTEGRATED-SDF-P (MAKE-ARRAY (LENGTH SDF)))))

given [1] sdf (required) ==> vector with sdf for a real or complex-valued process computed over a grid of frequencies (the grid need not be uniform) -- note that the sdf values must be untransformed (e.g., not in decibels) [2] freqs (required) ==> vector with frequencies over which sdf has been computed (must be same size as vector with sdf values) [3] fiddle-factor (keyword; 2.0) ==> values returned are all multiplied by this factor; the default of 2 is useful for real-valued processes since these have symmetric two-sided sdf's. [4] return-integrated-sdf-p (keyword; nil) ==> if t, integrated sdf as a function of frequency is computed and returned in result [5] result (keyword; vector of length sdf) <== vector with integrated sdf (if return-integrated-sdf-p is t) or nil (if return-integrated-sdf-p is nil) returns [1] the integral of the sdf * fiddle-factor [2] result, i.e., the integrated sdf if return-integrated-sdf-p is t

INTERQUARTILE-RANGE (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)) (THE-SEQ-IS-ORDERED-P NIL))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used [4] the-seq-is-ordered-p (keyword; nil) ==> if t, the-seq is assumed to be sorted from smallest to largest value; if nil, the-seq will be sorted returns [1] interquartile range; i.e, 0.75 quantile minus 0.25 quantile

INVERSE-DFT! (X &KEY (N (LENGTH X)) (SAMPLING-TIME 1.0))

given: [1] X (required) <=> a vector of real or complex-valued numbers [2] N (keyword; length of ) ==> number of points in inverse dft [3] sampling-time (keyword; 1.0) returns: [1] X, with contents replaced by the inverse discrete Fourier transform of X, namely, N-1 x(t) = (N sampling-time)^{-1} SUM X(n) exp(i 2 pi n t/N) n=0 --- Note: see Equation (111a) of the SAPA book

INVERSE-FFT! (COMPLEX-VECTOR &KEY (N (LENGTH COMPLEX-VECTOR)))

given: [1] complex-vector (required) <=> a vector of real or complex-valued numbers [2] N (keyword; length of ) ==> number of points in inverse dft (must be a power of 2) computes the inverse discrete Fourier transform of complex-vector using a fast Fourier transform algorithm and returns: [1] complex-vector, with contents replaced by the inverse discrete Fourier transform of the input X(n), namely, N-1 x(t) = 1/N SUM X(n) exp(i 2 pi n t/N) n=0 --- Note: see Equation (111a) of the SAPA book with the sampling time set to unity

IOTA (START END &KEY (TYPE 'VECTOR) (FLOAT-P NIL) (SKIP 1) (RESULT (MAKE-SEQUENCE TYPE (1+ (/ (- END START) SKIP)))))

given [1] start (required) ==> an integer [2] end (required) ==> an integer >= start [3] type (keyword; 'vector) ==> type of sequence desired [4] float-p (keyword; nil) ==> if true, sequence consists of floats [5] skip (keyword; 1) ==> desired increment between integers in sequence [6] result (keyword; sequence of length (1+ (- start end))) ==> storage space for result returns [1] a sequence of length (1+ (/ (- end start) skip)) with numbers start, start + skip, ..., end - skip, end (these are stored in whatever result points to)

JENKINS-SMOOTHING-WINDOW-BANDWIDTH (LAG-WINDOW-FUNCTION MAX-LAG &KEY (SAMPLING-TIME 1.0))

given a lag window function, a maximum lag and the sampling time, returns Jenkins' measure of smoothing window bandwidth --- Note: see Equation (242c) of the SAPA book

KERNEL-PDF-ESTIMATION (A-SEQ START-X INCREMENT-X N-PDF &KEY (WINDOW-WIDTH (WINDOW-WIDTH-FROM-SILVERMAN A-SEQ)) (RESULT-PDF (MAKE-ARRAY N-PDF INITIAL-ELEMENT 0.0) RESULT-PDF-SUPPLIED-P) (RESULT-X (MAKE-ARRAY N-PDF)))

given [1] a-seq (required) ==> a sequence of real-valued numbers [2] start-x (required) ==> first point at which pdf is to be estimated [3] increment-x (required) ==> increment between points at which pdf is to be estimated (i.e., second point is at start-x + increment-x, third point is at start-x + 2*increment-x, etc) [4] n-pdf (required) ==> number of points at which pdf is to be estimated [5] window-width (keyword; Silverman's formula) ==> window width of kernel pdf estimator [6] result-pdf (keyword; vector of length n-pdf) <== vector to hold pdf estimate [7] result-x (keyword; vector of length n-pdf) <== vector to hold points at which pdf was estimated computes a pdf estimate using the normal (Gaussian) kernel and returns [1] the pdf estimate (in result-pdf) [2] the points associated with the pdf estimates (in result-x) [3] the window width

LAG-WINDOW-SPECTRAL-ESTIMATE (ACVS LAG-WINDOW-FUNCTION &KEY (MAX-LAG (1- (LENGTH ACVS))) (N-TS (LENGTH ACVS)) (N-NONZERO-FREQS HALF-NEXT-POWER-OF-2) (RETURN-EST-FOR-0-FREQ-P NIL) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS (1+ MAX-LAG)))) (C_H 1.0) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (1+ MAX-LAG) RETURN-EST-FOR-0-FREQ-P))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (1+ MAX-LAG) RETURN-EST-FOR-0-FREQ-P)))))

given [1] acvs (required) ==> vector containing autocovariance sequence [2] lag-window-function (required) ==> function of a single variable that computes the value of the lag window for a given lag [3] max-lag (keyword; (1- (length acvs))) ==> maximum lag in acvs to be used [4] N-ts (keyword; length of acvs) ==> length of the time series from which acvs was constructed; this is needed to compute equivalent degrees of freedom [5] N-nonzero-freqs (keyword; :half-next-power-of-2) ==> specifies at how many nonzero frequencies direct spectral estimate is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [6] return-est-for-0-freq-p (keyword; nil) ==> if t, sdf is computed at zero frequency; otherwise, it is not computed. [7] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [8] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [9] C_h (keyword; 1.0) ==> variance inflation factor due to tapering [10] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [11] result-sdf (keyword; vector of correct length) <== vector into which lag window spectral estimates are placed; it must be exactly of the length dictated by N-nonzero-freqs and return-est-for-0-freq-p [12] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral estimate are computed and returned in result-freq [13] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [14] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length dictated by N-nonzero-freqs and return-est-for-0-freq-p into which the frequencies associated with the values in result-sdf are placed returns [1] result-sdf, a vector holding the properly transformed sdf [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-sdf -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-sdf [4] the equivalent degrees of freedom [5] the smoothing window bandwidth --- Note: see Section 6.7 of the SAPA book

LINEAR-INTERPOLATION! (A-SEQUENCE &KEY (INTERPOLATION-PREDICATE #'ZEROP))

given [1] a-sequence (required) <=> a sequence of numbers [2] interpolation-predicate (keyword; #'zerop) ==> a unitary predicate which, when true, indicates that a value in a-sequence is to be replaced by a linear interpolation returns [1] a-sequence, modified with linear interpolates over specified values

LINEARIZED-UPPER-TRANGULAR->2D-MATRIX (A &KEY (RESULT (LET ((N (ROUND (/ (1- (SQRT (1+ (* 8 (LENGTH A))))) 2)))) (MAKE-ARRAY `(,N ,N)))))

given [1] A (required) ==> a vector with an upper triangular matrix stored in linearized form [2] result (keyword; new 2d array of appropriate size) <== a 2d hermitian matrix filled according to contents of A converts linearized form A into 2d hermitian array and returns [1] result, a 2d hermitian matrix

LOG-OF-GAMMA (XX)

given xx (a real or complex valued number whose real part is greater than 0), returns the log (base e) of gamma(xx) --- Note: based upon the discussion in Section 6.1, Numerical Recipes, Second Edition

LOG10 (NUM)

given [1] num (required) ==> a number returns [1] log_10(num) --- Note: shorthand for (log num 10)

LOWER-TRIANGULAR-SOLVE! (A B &KEY (N (LENGTH B)))

given [1] A (required) ==> lower n by n triangular matrix (note: upper superdiagonal part of this matrix is ignored) [2] b (required) <=> on input, the right-hand side vector; on output, the solution X to A X = b [3] n (keyword; length of b) ==> order of matrix A returns [4] X, the solution to A X = b, where X is the same as vector which contained b (note that A is left unchanged)

M-LOCATION-ESTIMATE (THE-SEQ &KEY (WEIGHTING-FUNCTIONS `(,#'(LAMBDA (X) (HUBER-WEIGHT-FUNCTION X PARAMETER 2.0)) ,#'(LAMBDA (X) (THOMSON-WEIGHT-FUNCTION X BETA (/ (- (LENGTH THE-SEQ) 0.5) (LENGTH THE-SEQ)) SYMMETRIC-VERSION T)))) (ITERATIONS '(10 0)))

given [1] the-seq (required) ==> a sequence of numbers [2] weighting-functions (keyword; Huber and Thomson) ==> a list of two weighting functions [3] iterations (keyword; '(10 0)) ==> a list of nonegative integers giving the number of iterations for each weighting function returns [1] m-location estimate for sequence [2] MAD scale estimate (based upon deviations from m-location estimate) --- Note: see article by Hogg in Launer and Wilkinson

MAX-OF-SEQ (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] maximum value in sequence [2] index of maximum value

MEDIAN-ABSOLUTE-DEVIATION (THE-SEQ &KEY (LOCATION-ESTIMATE (SAMPLE-MEDIAN THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] location-estimate (keyword; median of the-seq) ==> robust estimate of location returns [1] the median absolute deviation about location-estimate [2] location-estimate (median of the-seq by default)

MIN-AND-MAX-OF-SEQ (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] minimum value in sequence [2] maximum value in sequence

MIN-OF-SEQ (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] minimum value in sequence [2] index of minimum value

MODIFIED-GRAM-SCHMIDT! (A-MATRIX)

given [1] A (required) <=> a matrix of real or complex-valued numbers with size m by n with rank N computes the factorization A = QR, where Q is of size m by n and has orthonormal columns and R is of size n by n and is upper triangular, and returns [1] Q [2] R --- Note: see Algorithm 5.2.5, p. 219, of Golub and Van Loan, Second Edition, with simple changes to handle matrices with real or complex-valued numbers

MULTIPLY-2-POLYNOMIALS (POLYNOMIAL-1 POLYNOMIAL-2 &KEY (DEGREE-OF-1 (1- (LENGTH POLYNOMIAL-1))) (DEGREE-OF-2 (1- (LENGTH POLYNOMIAL-2))) (PRODUCT-POLYNOMIAL (MAKE-ARRAY (+ DEGREE-OF-1 DEGREE-OF-2 1))))

given [1] polynomial-1 (required) ==> sequence of length n+1 with real or complex-valued numbers giving the n+1 coefficients of a polynomial of nth order; (elt polynomial-1 i) = coefficient of z^(n-i) [2] polynomial-2 (required) ==> another sequence representing another polynomial [3] degree-of-1 (keyword; (1- (length polynomial-1))) ==> degree of first polynomial [4] degree-of-2 (keyword; (1- (length polynomial-2))) ==> degree of second polynomial [5] product-polynomial (keyword; array of length (+ degree-of-1 degree-of-2 1)) <== a sequence with product of two polynomials returns [1] product-polynomial, an sequence with coefficients of polynomial given by the product of polynomial-1 and polynomial-2 --- Note: this routine works for a polynomial p represented either by p(0) + p(1)*z + ... + p(degp)*z^degp or by p(0)*z**degp + p(1)*z**(degp-1) + ... + p(degp)

MULTIPLY-MATRIX-AND-SCALAR (A-MATRIX SCALAR &KEY (RESULT (MAKE-ARRAY (ARRAY-DIMENSIONS A-MATRIX))))

given [1] a-matrix (required) ==> a 2d matrix [2] scalar (required) ==> an arbitrary number [3] result (keyword; new matrix of same size as a-matrix) <== a matrix to contain product of a-matrix and scalar returns [1] product of a-matrix and scalar (placed in result)

MULTIPLY-MATRIX-AND-VECTOR (A-MATRIX B-VECTOR &KEY (RESULT (MAKE-ARRAY (NTH 0 (ARRAY-DIMENSIONS A-MATRIX)))))

given [1] a-matrix (required) ==> a 2d matrix [2] b-vector (required) ==> a vector, with dimensions such that the product of a-matrix and b-vector is defined [3] result (keyword; new vector of appropriate size) <== a vector to contain product of a-matrix and b-vector returns [1] product of a-matrix and b-vector (placed in result)

MULTIPLY-POLYNOMIALS (&REST POLYS)

given [1] sequences representing any number of polynomials (required) returns [1] a sequence representing their product --- Note: this function was written by Andrew G. Bruce

MULTIPLY-SEQUENCES (&REST SEQS)

given [1] a arbitrary number of sequences of numbers, all of the same size returns [1] a new sequence formed by multiplying the sequences together on an element by element basis

MULTIPLY-TWO-MATRICES (A-MATRIX B-MATRIX &KEY (RESULT (MAKE-ARRAY (LIST (NTH 0 (ARRAY-DIMENSIONS A-MATRIX)) (NTH 1 (ARRAY-DIMENSIONS B-MATRIX))))))

given [1] a-matrix (required) ==> a 2d matrix [2] b-matrix (required) ==> another 2d matrix, with dimensions such that the product of a-matrix and b-matrix is defined [3] result (keyword; new 2d array of appropriate size) <== a 2d matrix to contain product of two matrices returns [1] product of two matrices (placed in result)

MULTITAPER-SPECTRAL-ESTIMATE (TIME-SERIES LIST-OF-DATA-TAPERS &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (N-TAPERS (LENGTH LIST-OF-DATA-TAPERS)) (N-NONZERO-FREQS HALF-NEXT-POWER-OF-2) (RETURN-EST-FOR-0-FREQ-P NIL) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS (- END START)))) (RECENTER-AFTER-TAPERING-P T) (RESTORE-POWER-OPTION-P T) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P))) (RETURN-EIGENPSPECTRA-P T) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P)))))

given [1] time-series (required) ==> a sequence of time series values [2] list-of-data-tapers (required) ==> a list of orthonormal data tapers, each of length (- end start) [3] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [4] start (keyword; 0) ==> start index of sequence to be used [5] end (keyword; length of time-series) ==> 1 + end index of sequence to be used [6] N-tapers (keyword; length of list-of-data-tapers) ==> number of data tapers to be used in list-of-data-tapers [7] N-nonzero-freqs (keyword; :half-next-power-of-2) ==> specifies at how many nonzero frequencies multitaper spectral estimate is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [8] return-est-for-0-freq-p (keyword; nil) ==> if t, sdf is computed at zero frequency; otherwise, it is not computed. [9] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [10] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [11] recenter-after-tapering-p (keyword; t) ==> if t and data-taper is a function, centers tapered series by subtracting off its sample mean [12] restore-power-option-p (keyword; t) ==> if t and data-taper is a function, normalizes tapered series to have same sum of squares as before tapering [13] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [14] result-sdf (keyword; vector of correct length) <== vector into which multitaper spectral estimate is placed; it must be exactly of the length dictated by N-nonzero-freqs and return-est-for-0-freq-p [15] return-eigenpspectra-p (keyword; t) ==> if t, individual eigenspectra are returned in a list (each eigenspectrum in the list is associated with the corresponding taper in list-of-data-tapers) [16] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral estimate are computed and returned in result-freq [17] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [18] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length N-nonzero-freqs (if return-est-for-0-freq-p is nil) or N-nonzero-freqs +1 (if return-est-for-0-freq-p is t) into which the frequencies associated with the values in result-sdf are placed returns [1] result-sdf, a vector holding the properly transformed multitaper spectral estimate [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-sdf, or nil (if return-frequencies-p is nil) [3] the length of the vector result-sdf [4] a list of untransformed eigenspectra (return-eigenpspectra-p is t), or nil (if return-eigenpspectra-p is nil) --- Note: see Section 7.1 of the SAPA book

N-APPLICATIONS-OF-THREE-POINT-SMOOTHER (A-SEQ N)

given [1] a-seq (required) ==> a sequence of length 2*n + 1 or greater [2] n (required) ==> positive integer indicating the number of times the three-point smoother is to be applied returns [1] the result of repetitively smoothing a-seq using the function three-point-smoother

NEWTON-RAPHSON (F F-PRIME X-LEFT X-RIGHT &KEY (ACCURACY (* 10.0 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 20) (PREVENT-BRACKET-JUMPING-P T))

given [1] f (required) ==> a function with a single argument [2] f-prime (required) ==> another function with a single argument, this one being the first derivative of f [3] x-left (required) ==> left-hand bracket for the desired root; i.e., left-hand bracket <= desired root [4] x-right (required) ==> right-hand bracket for the desired root; i.e., desired root <= right-hand bracket [5] accuracy (keyword; (* 10.0 single-float-epsilon)) ==> desired relative accuracy for computed rood [6] maximum-number-of-iterations (keyword; 20) ==> maximum number of iterations [7] prevent-bracket-jumping-p (keyword; t) ==> if t, allows Newton-Raphson to continue if it jumps out of the interval [x-left, x-right]; if nil, jumping out of the interval causes an error to be signaled returns [1] a root of f in [x-left, x-right]; i.e., a value x in that interval such that f(x) = 0 [2] the number of iterations required --- Note: this function is based loosely on rtnewt, Section 9.4, Numerical Recipes, Second Edition

NEXT-POWER-OF-2 (N)

given [1] n (required) ==> an integer returns [1] an integer power of 2 greater than or equal to n

NYQUIST-FREQUENCY->SAMPLING-TIME (NYQUIST-FREQUENCY)

given a Nyquist frequency, returns the associated sampling-time --- Note: see page 98 of the SAPA book

ONE-SIDED-FREQ->TWO-SIDED-FREQ (ONE-SIDED-FREQ &KEY (RESULT (MAKE-ARRAY (1- (* 2 (LENGTH ONE-SIDED-FREQ))))))

given a vector of N frequencies 0, f_1, ..., f_{N-1}, returns the vector of 2N-1 frequencies -f_{N-1}, ..., -f_1, 0, f_1, ..., f_{N-1}

ONE-SIDED-SDF->TWO-SIDED-SDF (ONE-SIDED-SDF &KEY (RESULT (MAKE-ARRAY (1- (* 2 (LENGTH ONE-SIDED-SDF))))))

given a one-sided sdf S(0), S(f_1), ..., S_(f_{N-1}) of length N, returns the two-sized sdf S_(-f_{N-1}), ..., S(-f_1), S(0), S(f_1), ..., S_(f_{N-1}) of length 2N-1, where S_(-f_k) = S_(f_k)

ORDINARY-LEAST-SQUARES-CHOLESKY (DEPENDENT-VARIABLE LIST-OF-INDEPENDENT-VARIABLES &KEY (COMPUTE-RESIDUALS-P NIL) (RESULT (IF COMPUTE-RESIDUALS-P (MAKE-ARRAY (LENGTH DEPENDENT-VARIABLE)))))

given [1] dependent-variable (required) ==> a vector [2] list-of-independent-variables (required) ==> a list of either vectors or functions (must be ALL vectors or ALL functions) [3] compute-residuals-p (keyword; nil) ==> if t, residuals are computed [4] result (keyword; new array if compute-residuals-p is t; otherwise nil) <== storage space for residuals (not used unless compute-residuals-p is true) returns [1] vector with estimated parameters [2] residuals --- Note: uses Lisp versions of Cholesky factorization routines from linpack

ORDINARY-LEAST-SQUARES-Q-R (DEPENDENT-VARIABLE LIST-OF-INDEPENDENT-VARIABLES &KEY (COMPUTE-RESIDUALS-P NIL) (COMPUTE-STANDARD-ERRORS-P NIL) (COMPUTE-CORRELATION-MATRIX-P NIL) (RESIDUALS (IF (OR COMPUTE-RESIDUALS-P COMPUTE-STANDARD-ERRORS-P COMPUTE-CORRELATION-MATRIX-P) (MAKE-ARRAY (LENGTH DEPENDENT-VARIABLE)))))

given [1] dependent-variable (required) ==> a vector [2] list-of-independent-variables (required) ==> a list of either vectors, functions or mixture thereof [3] compute-residuals-p (keyword; nil) ==> if t, residuals are computed [4] compute-standard-errors-p (keyword; nil) ==> if t, standard errors for parameter estimates are computed, along with residuals (even if compute-residuals-p is nil) [5] compute-correlation-matrix-p (keyword; nil) ==> if t, correlation matrix for parameter estimates is computed, along with residuals and standard errors (even if either of their associated keywords is nil) [6] residuals (keyword; nil or vector) <== storage space for residuals (nil unless either compute-residuals-p, compute-standard-errors-p or compute-correlation-matrix-p is true) returns [1] vector with parameter estimates [2] vector with residuals [2] standard error of residuals [3] vector with standard errors of parameter estimates [4] correlation matrix --- Note: this routine is based on the discussion in ``Nonlinear Regression Analysis and Its Applications'' by Bates and Watts, Wiley, 1988.

PAPOULIS-BANDWIDTH->M (B_W &KEY (SAMPLING-TIME 1.0))

given desired smoothing window bandwidth B_W and sampling time, returns [1] window parameter m required to approximately achieve B_W using the Papoulis lag window [2] actual B_W achieved --- Note: see Table 269 of the SAPA book

PAPOULIS-LAG-WINDOW (TAU M)

given the lag tau and window parameter m, returns the value of the Papoulis lag window --- Note: see equation near bottom of page 266 of the SAPA book

PAPOULIS-M->BANDWIDTH (M &KEY (SAMPLING-TIME 1.0))

given window parameter m and sampling time, returns bandwidth B_W for the Papoulis smoothing window --- Note: see Table 269 of the SAPA book

PAPOULIS-N-M->DEGREES-OF-FREEDOM (N M &KEY (C_H 1.0))

given sample size N, window parameter m and variance inflation factor C_h, returns equivalent degrees of freedom nu for Papoulis lag window --- Note: see Table 269 of the SAPA book

PARZEN-BANDWIDTH->M (B_W &KEY (SAMPLING-TIME 1.0))

given desired smoothing window bandwidth B_W and sampling time, returns [1] window parameter m required to approximately achieve B_W using the Parzen lag window [2] actual B_W achieved --- Note: see Table 269 of the SAPA book

PARZEN-LAG-WINDOW (TAU M)

given the lag tau and window parameter m, returns the value of the Parzen lag window --- Note: see the equation on page 265 of the SAPA book or Priestley, page 443, Equation (6.2.82)

PARZEN-M->BANDWIDTH (M &KEY (SAMPLING-TIME 1.0))

given window parameter m and sampling time, returns bandwidth B_W for the Parzen smoothing window --- Note: see Table 269 of the SAPA book

PARZEN-N-M->DEGREES-OF-FREEDOM (N M &KEY (C_H 1.0))

given sample size N, window parameter m and variance inflation factor C_h, returns equivalent degrees of freedom nu for Parzen lag window --- Note: see Table 269 of the SAPA book

PERIODOGRAM (TIME-SERIES &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (N-NONZERO-FREQS HALF-NEXT-POWER-OF-2) (RETURN-EST-FOR-0-FREQ-P NIL) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS (- END START)))) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS (- END START) RETURN-EST-FOR-0-FREQ-P)))))

given [1] time-series (required) ==> a vector of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of time-series) ==> 1 + end index of vector to be used [5] N-nonzero-freqs (keyword; :half-next-power-of-2) ==> specifies at how many nonzero frequencies periodogram is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [6] return-est-for-0-freq-p (keyword; nil) ==> if t, periodogram is computed at zero frequency; otherwise, it is not computed. [7] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [8] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [9] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [10] result-sdf (keyword; vector of correct length) <== vector into which periodogram is placed; it must be exactly of the length dictated by N-nonzero-freqs and return-est-for-0-freq-p [11] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the periodogram are computed and returned in result-freq [12] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [13] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length dictated by N-nonzero-freqs and return-est-for-0-freq-p into which the frequencies associated with the values in result-sdf are placed returns [1] result-sdf, a vector holding the properly transformed periodogram [2] result-freq (if return-frequencies-p is t), a vector holding the properly transformed frequencies associated with values in result-sdf -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-sdf --- Note: see Section 6.3 of the SAPA book

PERIODOGRAM-AT-ONE-FREQUENCY (TIME-SERIES FREQ &KEY (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (SAMPLING-TIME 1.0))

given [1] time-series (required) ==> a vector of time series values [2] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of time-series) ==> 1 + end index of vector to be used [5] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) returns [1] value of periodogram at freq [2] approximate conditional least squares estimate for A, the amplitude of cosine term in Equation (461a) of the SAPA book [3] approximate conditional least squares estimate for B the amplitude of sine term --- Note: see Section 10.2 of the SAPA book

POSTCOLOR-SPECTRAL-ESTIMATE (PRECOLORED-SDF-ESTIMATE PREWHITENING-FILTER SAMPLE-SIZE &KEY (N-NONZERO-FREQS HALF-NEXT-POWER-OF-2) (INCLUDES-EST-FOR-0-FREQ-P NIL) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS SAMPLE-SIZE))) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (LENGTH PRECOLORED-SDF-ESTIMATE))))

given [1] precolored-sdf-estimate (required) ==> vector containing sdf estimate to be postcolored; note that this estimate is assumed to be untransformed (i.e., not expressed in dB, etc) [2] prewhitening-filter (required) ==> vector with coefficients of prewhitening filter [3] sample-size (required) ==> length of the time series from which precolored-sdf-estimate was constructed; this is needed to get corresponding grid of frequencies [4] N-nonzero-freqs (keyword; :half-next-power-of-2) ==> specifies at how many nonzero frequencies precolored-sdf-estimate was computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [5] includes-est-for-0-freq-p (keyword; nil) ==> if t, first element of precolored-sdf-estimate corresponds to zero frequency [6] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [7] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [8] result-sdf (keyword; vector of correct length) <== vector into which postcolored sdf estimate is placed; it must be exactly the same length as precolored-sdf-estimate returns [1] result-sdf, a vector holding the postcolored sdf estimate --- Note: see Equation (438) of the SAPA book

POWER-OF-2 (N)

given [1] n (required) ==> an integer returns [1] nil if n is not a power of 2 -- or -- m if n = 2^m for a positive integer m

PREDICTED-Y-AT-X (X INTERCEPT SLOPE)

given [1] x (required) ==> a number [2] intercept (required) ==> intercept of a line [3] slope (required) ==> slope of a line returns [1] intercept + slope*x

Q-Q-PLOT (Y-VALS X-VALS &KEY (N (IF (TYPEP Y-VALS 'SEQUENCE) (LENGTH Y-VALS) 256)) (X-VALS-SORTED-P (FUNCTIONP X-VALS)) (Y-VALS-SORTED-P (FUNCTIONP Y-VALS)))

given [1] y-vals (required) ==> a sequence or function [2] x-vals (required) ==> a sequence or function [3] n (keyword; 256 if y-vals is a function, length of y-vals otherwise) ==> positive integer [4] x-vals-sorted-p (keyword; t if x-vals is a function, nil otherwise) ==> a flag indicating whether x-vals is already sorted [5] y-vals-sorted-p (keyword; t if y-vals is a function, nil otherwise) ==> a flag indicating whether y-vals is already sorted return [1] y values needed to create a q-q-plot for y-vals [2] x values --- Note: For details, see Chapter 6 of ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner and Tukey

Q-R! (A-MATRIX)

same as modified-Gram-Schmidt!

QUANTILE-OF-CHI-SQUARE-2-DISTRIBUTION (P)

given [1] p (required) ==> a number > 0 and < 1 returns [1] p-th quantile of chi-square distribution with 2 degrees of freedom --- Note: Equation (31), P. 640 of Chave, Thomson, and Ander (1987)

QUANTILE-OF-CHI-SQUARE-DISTRIBUTION (NU P &KEY (ACCURACY-LEVEL QUICK))

given [1] nu (required) ==> degrees of freedom [2] p (required) ==> a number > 0 and < 1 [3] accuracy-level (keyword; :quick) if :quick, gives quick approximation; if :accurate, give accurate, but computationally expensive, approximation returns [1] p-th quantile of chi-square distribution with nu degrees of freedom --- Note: for :quick, see Table 6.4, page 225, ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, & Tukey; if :accurate, see AS 91, Applied Statistics

QUANTILE-OF-EXPONENTIAL-DISTRIBUTION (P)

given [1] p (required) ==> a number > 0 and < 1 returns [1] p-th quantile of exponential distribution with mean 1 --- Note: Table 6.4, page 225, ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, & Tukey

QUANTILE-OF-GAMMA-DISTRIBUTION (ALPHA P)

given [1] alpha (required) ==> shape parameter for gamma distribution [2] p (required) ==> a number > 0 and < 1 returns [1] p-th quantile of standard gamma distribution --- Note: Table 6.4, page 225, ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, & Tukey

QUANTILE-OF-KOLMOGOROV-TEST-STATISTIC (SAMPLE-SIZE &KEY (SIGNIFICANCE-LEVEL 0.95))

given [1] sample-size (required) ==> a positive integer [2] significance-level (keyword; 0.95) ==> 0.99, 0.98, 0.95, 0.90, or 0.80; 0.75 is also ok if sample-size > 40 returns [1] quantile of two-sided Kolmogorov test statistic --- Note: see Table 14, page 462, of Conover (2nd edition, 1980); Stephens (1974); and Diggle (1990, page 55); significance-level is currently limited to one of these values: 0.99, 0.98, 0.95, 0.90, or 0.80; if sample-size is greater than 40, 0.75 is also ok

QUANTILE-OF-NORMAL-DISTRIBUTION (P &KEY (ACCURACY-LEVEL QUICK))

given [1] p (required) ==> a number > 0 and < 1 [2] accuracy-level (keyword; :quick) ==> if :quick, gives quick approximation; if :better, gives slower, but better, approximation; if :accurate, uses slowest, but more accurate, approximation; returns [1] quantile of standard normal distribution; --- Note: for :quick, see Table 6.5, page 227, ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, & Tukey; for :better, see Sections 26.2.22 and 26.2.23 of Abramowitz and Stegun; for :accurate, see AS 111 by Beasley and Springer, Applied Statistics, 1977, vol. 26, p.118

QUANTILE-OF-ORDERED-SEQ (THE-SEQ P)

given [1] the-seq (required) ==> a sequence of real-valued numbers, ordered from smallest to largest [2] p (required) ==> a percentile; i.e., 0 <= p <= 1 returns [1] the sample quantile for p --- Note: see Section 2.2 of ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, and Tukey

QUANTILE-PLOT (A-SEQ &KEY (A-SEQ-SORTED-P NIL))

given [1] a-seq (required) ==> any sequence of real-valued numbers [2] a-seq-sorted-p (keyword; nil) ==> if t, a-seq is already sorted; if nil, a copy of a-seq is sorted for use by the function returns [1] sequence with sample quantiles for a-seq [2] sequence with corresponding percentiles --- Note: see Section 2.2 of ``Graphical Methods for Data Analysis'' by Chambers, Cleveland, Kleiner, and Tukey

RANORM

returns a random deviate from a normal distribution with zero mean and unit variance

RANORMS (N &KEY (RESULT (MAKE-ARRAY N)))

given: [1] n (required) ==> number of normal random deviates to be generated [2] result (keyword; vector of length n) <== vector to hold random deviates returns: [1] result, a vector with n normal random deviates (i.e., a realization of length n of a white noise process from a normal distribution)

REFLECTION-COEFFS->AR-COEFFS (REFLECTION-COEFFS &KEY (P (LENGTH REFLECTION-COEFFS)) (SCRATCH (MAKE-ARRAY P)) (RESULT (MAKE-ARRAY P)))

given [1] reflection-coeffs (required) ==> vector of length p with real-valued reflection coefficients phi_{1,1}, phi_{2,2}, ... , phi_{p,p} [2] p (keyword; length of reflection-coeffs) ==> AR model order [3] scratch (keyword; vector of length p) ==> vector of size p used for intermediate results [4] results (keyword; vector of length p) <== vector of size p into which are placed the AR coefficients phi_{j,p} for model x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t returns [1] results, i.e., the AR coefficients --- Note: see Section 9.4 of the SAPA book

REFLECTION-COEFFS->VARIANCE (REFLECTION-COEFFS INNOVATIONS-VARIANCE &KEY (P (LENGTH REFLECTION-COEFFS)))

given [1] reflection-coeffs (required) ==> vector of length p with real-valued reflection coefficients phi_{1,1}, phi_{2,2}, ... , phi_{p,p} [2] innovations-variance (required) ==> innovations variance of the process [3] p (keyword; length of reflection-coeffs) ==> AR model order returns [1] variance of the process --- Note: this is a recursive application of Equation (404c) of the SAPA book

REGRESSION-WITH-AR-ERRORS (TIME-SERIES IND-VARS P &KEY (AR-COEFF-ESTIMATOR #'BURG-ALGORITHM) (MAXIMUM-NUMBER-OF-ITERATIONS 5) (EPS 1.e-7) (TRACE-AR-PARAMETER-ESTIMATES-P NIL) (TRACE-LS-PARAMETER-ESTIMATES-P NIL))

given [1] time-series (required) ==> a vector with dependent variables [2] ind-vars (required) ==> a list of vector(s) with independent variable(s) [3] p (required) ==> autoregressive order [4] AR-coeff-estimator (keyword; #'burg-algorithm) ==> AR estimation procedure [5] maximum-number-of-iterations (keyword; 5) ==> ... like I said ... [6] eps (keyword; 1.0e-7) ==> a small number -- used to test for convergence in estimates of both regression and AR coefficients [7] trace-AR-parameter-estimates-p (keyword; nil) ==> if t, results of each iteration on AR parameters are printed [8] trace-LS-parameter-estimates-p (keyword; nil) ==> if t, results of each iteration on least squares parameters are printed returns [1] estimates of least squares coefficients [2] estimates of AR coefficients [3] variance of residuals [4] either number of iterations required for convergence or nil if maximum number of iterations completed but convergence not yet achieved --- Note: this function is based on Newton, TIMESLAB, page 241; it has NOT been thoroughly tested (particularly with regard to the convergence criterion), so a heavy dose of the usual ``caveat emptor'' is appropriate here!!!

RUNNING-MEDIAN (TIME-SERIES K &KEY (RESULT (MAKE-ARRAY (- (LENGTH TIME-SERIES) (1- K)))))

given [1] time-series (required) ==> a sequence of numbers [2] K (required) ==> positive integer giving the number of points in the running median [3] result (keyword; vector of appropriate length) <== a sequence to hold running medians of time series return [1] result, the vector of running medians of K consecutive points

SAMPLE-CEPSTRUM (LOG-SPECTRUM &KEY (SAMPLING-TIME 1.0) (ZERO-FREQUENCY-P NIL) (NYQUIST-FREQUENCY-P T) (RESULT (MAKE-ARRAY (IF ZERO-FREQUENCY-P (LENGTH LOG-SPECTRUM) (1+ (LENGTH LOG-SPECTRUM))))))

given [1] log-spectrum (required) ==> a vector containing log spectrum [2] sampling-time (keyword; 1.0) ==> the sample time (i.e., delta t) [3] zero-frequency-p (keyword; nil) ==> if t, the first element of log-spectrum is associated with 0 frequency; if nil, first element goes with lowest nonzero frequency [4] Nyquist-frequency-p (keyword; nil) ==> if t, the last element of log-spectrum is associated with Nyquist frequency; if nil, last element goes with a frequency just less than Nyquist [5] result (keyword; vector of appropriate length) <== vector to hold sample cepstrum returns [1] result, i.e., the sample cepstrum from lag 0 to lag N_U where N_U = length of log-spectrum - 1 if zero-frequency-p is t or = length of log-spectrum if zero-frequency-p is nil --- Note: see Equation 282 of the SAPA book

SAMPLE-CORRELATION-COEFFICIENT (SEQ-1 SEQ-2 &KEY (START 0) (END (LENGTH SEQ-1)))

given [1] seq-1 (required) ==> a sequence of real-valued numbers [2] seq-2 (required) ==> another sequence of real-valued numbers; should have the same length as seq-1 [3] start (keyword; 0) ==> start index of sequence to be used [4] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sample correlation coefficient between two sequences

SAMPLE-FROM-A-FUNCTION (A-FUNCTION &KEY (SEQ-OF-X-VALUES 'NIL) (X0 0.0) (DELTA-X 1.0) (N (LENGTH SEQ-OF-X-VALUES)))

given [1] a-function (required) ==> a function to be sampled from [2] seq-of-x-values (keyword; '()) ==> sequence of values at which a-function is to be sampled (default '()) [3] x0 (keyword; 0.0) ==> point of first value to be sampled [4] delta-x (keyword; 1.0) ==> increment between points [5] n (keyword; length of seq-of-x-values) ==> number of samples returns [1] an array with the values of a-function at a specified set of points; and [2] an array with the specified set of points. --- Note that the specified set of points either is in seq-of-x-values -- if it is non-nil -- or is given by x0, x0 + delta-x, x0 + 2*delta-x, ..., x0 + (n-1)*delta-x

SAMPLE-MEAN (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real or complex-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sample mean of the specified numbers in the-seq

SAMPLE-MEAN-AND-VARIANCE (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real or complex-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sample mean of the specified numbers in the-seq [2] sample variance of the specified numbers in the-seq

SAMPLE-MEDIAN (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)) (THE-SEQ-IS-ORDERED-P NIL))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used [4] the-seq-is-ordered-p (keyword; nil) ==> if t, the-seq is assumed to be sorted from smallest to largest value; if nil, the-seq will be sorted returns [1] sample median [2] minimum value in sequence [3] maximum value [4] sorted sequence

SAMPLE-SKEWNESS-AND-KURTOSIS (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sample skewness of the specified numbers in the-seq [2] sample kurtosis [3] sample mean [4] sample variance

SAMPLE-VARIANCE (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real or complex-valued numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sample variance of the specified numbers in the-seq

SAMPLE-VARIOGRAM (TIME-SERIES TIMES)

given: [1] time-series (required) ==> a vector containing the time series [2] times (required) ==> a vector containing the associated times returns [1] the sample variogram [2] associated lags --- Note: see Diggle's book

SAMPLING-TIME->NYQUIST-FREQUENCY (SAMPLING-TIME)

given a sampling time, returns the associated Nyquist frequency --- Note: see page 98 of the SAPA book

SDF->SDF-WITH-ACCURATE-PEAKS (SDF FREQS AR-COEFFS INNOVATIONS-VARIANCE &KEY (SAMPLING-TIME 1.0) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (SEARCH-METHOD QUADRATIC-APPROX) (DB-TOLERANCE 0.1) (ACCURACY (* 10.0 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 20))

given [1] sdf (required) ==> vector with sdf for a real-valued process computed over a grid of frequencies (the grid need not be uniform) -- note that the sdf values can be either untransformed or transformed via an order preserving transformation (such as decibels) [2] freqs (required) ==> vector with frequencies corresponding to values in sdf (must be same size as vector with sdf values) [3] AR-coeffs (required) ==> vector with real-valued AR coefficients phi_{j,p} for an AR(p) model of order p: x_t = phi_{1,p} * x_{t-1} + phi_{2,p} * x_{t-2} + ... + phi_{p,p} * x_{t-p} + e_t [4] innovations-variance (required) ==> innovations variance of the process [5] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [6] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform any newly computed values of the sdf [7] search-method (keyword; :quadratic-approx) ==> selects between two methods to search for the location of peaks; if set to :quadratic-approx, a quadratic approximation is used; if set to :bisection+Newton-Raphson, an interval search using a combination of bisection and Newton-Raphson is used [8] dB-tolerance (keyword; 0.1) ==> convergence criterion (used if search-method is :quadratic-approx) [9] accuracy (keyword; 10 * single-float-epsilon) ==> convergence criterion (used if search-method is :bisection+Newton-Raphson) [10] maximum-number-of-iterations (keyword; 20) ==> controls the number of iterations (used if search-method is :bisection+Newton-Raphson) returns [1] a vector with values of sdf merged with an additional set of values chosen to represent peaks accurately [2] the corresponding augmented vector of frequencies [3] the length of the two vector returned --- Note: see pages 524--5 of the SAPA book

SECANT-METHOD (F X-LEFT X-RIGHT &KEY (ACCURACY (* 10.0 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 50))

given [1] f (required) ==> a function with a single argument [2] x-left (required) ==> left-hand bracket for the desired root; i.e., left-hand bracket <= desired root [3] x-right (required) ==> right-hand bracket for the desired root; i.e., desired root <= right-hand bracket [4] accuracy (keyword; (* 10.0 single-float-epsilon)) ==> desired relative accuracy for computed rood [5] maximum-number-of-iterations (keyword; 50) ==> maximum number of iterations returns [1] a root of f in [x-left, x-right]; i.e., a value x in that interval such that f(x) = 0 [2] the number of iterations required --- Note: this function is based loosely on rtsec, Section 9.2, Numerical Recipes, Second Edition

SIGN (A1 A2)

implementation of Fortran SIGN function

SIMPLE-NUMERICAL-INTEGRATION (F A B &KEY (ACCURACY 1.e-6) (MAXIMUM-NUMBER-OF-ITERATIONS 20))

given [1] f (required) ==> a function with a single argument [2] a (required) ==> left-hand limit for numerical integration [3] b (required) ==> right-hand limit for numerical integration i.e., a < b [4] accuracy (keyword; 1.0e-6) ==> desired relative accuracy for computed rood [5] maximum-number-of-iterations (keyword; 20) ==> maximum number of iterations returns [1] the integral of f over the interval [a, b] [2] the number of iterations required --- Note: this function is based on qtrap, Section 4.2, Numerical Recipes, Second Edition

SIMULATE-TIME-SERIES-FROM-SDF (N-SERIES SDF &KEY (SAMPLING-TIME 1.0) (N-TOTAL (* 4 (NEXT-POWER-OF-2 N-SERIES))) (RESULT (MAKE-ARRAY N-SERIES)))

given [1] n-series (required) ==> length of time series to be simulated [2] sdf (required) ==> spectral density function (sdf) defined for 0 <= f <= 1/(2.0 sampling-time) -- the sdf is assumed to be two-sided and symmetric about f = 0 [3] sampling-time (keyword; 1.0) ==> the assumed sampling time (delta t) [4] n-total (keyword; 4 * next power of 2 for n-series) ==> a power of 2 controlling degree of accuracy of the approximation (the larger, the better -- see Percival, 1992, for details) [5] result (keyword; vector of length n-series) <== a vector to contain simulated time series returns [1] a vector of length n-series generated from sdf --- Note: the method used here is an approximate frequency domain technique; in particular, it will NOT simulate the DC component correctly, so beware of using this function in studies where the process mean is of important

SMOOTHING-WINDOW-FOR-LAG-WINDOW-SPECTRAL-ESTIMATE (SAMPLE-SIZE LAG-WINDOW-FUNCTION &KEY (N-NONZERO-FREQS TWICE-NEXT-POWER-OF-2) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS SAMPLE-SIZE))) (SMOOTH-WIND-TRANSFORMATION #'CAREFUL-CONVERT-TO-DB) (RESULT-SMOOTH-WIND (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T)))))

given [1] sample-size (required) ==> sample size for which smoothing window is to be computed [2] lag-window-function (required) ==> function of a single variable that computes the value of the lag window for a given lag [3] N-nonzero-freqs (keyword; :twice-next-power-of-2) ==> specifies at how many nonzero frequencies smoothing window is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [4] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [5] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [6] smooth-wind-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-smooth-wind [7] result-smooth-wind (keyword; vector of correct length) <== vector into which smoothing window values are placed; it must be exactly of the length dictated by N-nonzero-freqs [8] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the smoothing window are computed and returned in result-freq [9] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [10] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of N-nonzero-freqs +1 into which the frequencies associated with the values in result-smooth-wind are placed returns [1] result-smooth-wind, a vector holding the properly transformed smoothing window [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-smooth-wind -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-smooth-wind --- Note: see equation below Equation (237c) of the SAPA book

SPECTRAL-WINDOW-FOR-DIRECT-SPECTRAL-ESTIMATE (SAMPLE-SIZE &KEY (N-NONZERO-FREQS TWICE-NEXT-POWER-OF-2) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS SAMPLE-SIZE))) (DATA-TAPER NIL) (DATA-TAPER-PARAMETERS) (SPEC-WIND-TRANSFORMATION #'CAREFUL-CONVERT-TO-DB) (RESULT-SPEC-WIND (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T)))))

given [1] sample-size (required) ==> sample size for which spectral window is to be computed [2] N-nonzero-freqs (keyword; :twice-next-power-of-2) ==> specifies at how many nonzero frequencies spectral window is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [3] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [4] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [5] data-taper (keyword; nil) ==> nil or a tapering function [6] data-taper-parameters (keyword) ==> parameters for tapering function (not used if data-taper is nil) [7] spec-wind-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-spec-wind [8] result-spec-wind (keyword; vector of correct length) <== vector into which spectral window values are placed; it must be exactly of the length dictated by N-nonzero-freqs [9] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral window are computed and returned in result-freq [10] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq [11] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of N-nonzero-freqs +1 into which the frequencies associated with the values in result-spec-wind are placed returns [1] result-spec-wind, a vector holding the properly transformed spectral window [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-spec-wind -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-spec-wind --- Note: see equation below Equation (207a) of the SAPA book

SPECTRAL-WINDOW-FOR-LAG-WINDOW-SPECTRAL-ESTIMATE (SAMPLE-SIZE LAG-WINDOW-FUNCTION &KEY (N-NONZERO-FREQS TWICE-NEXT-POWER-OF-2) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (GET-N-DFT N-NONZERO-FREQS SAMPLE-SIZE))) (DATA-TAPER NIL) (DATA-TAPER-PARAMETERS) (SPEC-WIND-TRANSFORMATION #'CAREFUL-CONVERT-TO-DB) (RESULT-SPEC-WIND (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T))) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (GET-N-FREQS N-NONZERO-FREQS SAMPLE-SIZE T)))))

given [1] sample-size (required) ==> sample size for which spectral window is to be computed [2] lag-window-function (required) ==> function of a single variable that computes the value of the lag window for a given lag [3] N-nonzero-freqs (keyword; :twice-next-power-of-2) ==> specifies at how many nonzero frequencies spectral window is to be computed -- choices are: :half-next-power-of-2 ==> 1/2 * power of two >= sample size; :next-power-of-2 ==> power of two >= sample size; :twice-next-power-of-2 ==> 2 * power of two >= sample size; :Fourier ==> just at Fourier frequencies -- or -- any power of 2 >= 1/2 * [power of two >= sample size] [4] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [5] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [6] data-taper (keyword; nil) ==> nil or a tapering function [7] data-taper-parameters (keyword) ==> parameters for tapering function (not used if data-taper is nil) [8] spec-wind-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-spec-wind [9] result-spec-wind (keyword; vector of correct length) <== vector into which spectral window values are placed; it must be exactly of the length dictated by N-nonzero-freqs [10] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral window are computed and returned in result-freq [11] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [12] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of N-nonzero-freqs +1 into which the frequencies associated with the values in result-spec-wind are placed returns [1] result-spec-wind, a vector holding the properly transformed spectral window [2] result-freq (if return-frequencies-p is t), where result-freq is a vector holding the properly transformed frequencies associated with values in result-spec-wind -- or -- nil (if return-frequencies-p is nil) [3] the length of the vector result-spec-wind --- Note: see equation below Equation (244a) of the SAPA book

SPOFA! (A &KEY (N (CAR (ARRAY-DIMENSIONS A))))

given [1] A (required) <=> real symmetric positive definite matrix (this gets trashed) [2] n (keyword; (array-dimensions a)) ==> size of matrix to be factored calculates the upper triangular matrix U of the Cholesky decomposition A = LU, where L = U^T, and returns [1] U, stored in upper triangular part of A

SPOSL! (A B &KEY (N (LENGTH B)))

given [1] A (required) ==> real symmetric positive definite matrix AFTER it has been crunched by spofa! [2] b (required) <=> the right-hand side vector on input; on output, this gets replaced by the solution X [3] N (keyword; length of b) ==> order of matrix A (default is length of b) returns [1] X, the solution to A X = b; note that X is the same as the vector which contained b

STANDARD-NORMAL-PDF (X)

given x, returns value of standard normal pdf at x

STEP-DOWN-LEVINSON-DURBIN-RECURSIONS (COEFFS VARIANCE &KEY (PROCESS-VARIANCE-P T))

given [1] coeffs (required) ==> sequence of length p with AR coefficients: x_t = coeffs_0*x_{t-1} + coeffs_1*x_{t-2} + ... + coeffs_{p-1}*x_{t-p} + e_t (see Equation (392a) in the SAPA book; the coefficients can be real or complex-valued) [2] variance (required) ==> process variance or innovations variance (see keyword process-variance-p) [3] process-variance-p (keyword; t) ==> if t, variance is taken to be process variance if nil, variance is taken to be innovations variance computes best linear prediction coefficients of orders 1, 2, ..., p-1 and prediction error variances of orders 0 (process variance), 1, ..., p and returns [1] list of vectors with best linear prediction coefficients going from order 1 to order p; [2] list of prediction error variances going from order 0 to order p --- Note: see item [4] of the Comments and Extensions to Section 9.4 of the SAPA book. The values returned by this function can be used to set the keyword parameters list-of-lower-order-phi and list-of-lower-order-pev in the function generate-ar-time-series

SUBTRACT-TWO-MATRICES (A-MATRIX B-MATRIX &KEY (RESULT (MAKE-ARRAY (ARRAY-DIMENSIONS A-MATRIX))))

given [1] a-matrix (required) ==> a 2d matrix [2] b-matrix (required) ==> a 2d matrix, with dimensions the same as a-matrix [3] result (keyword; new vector of appropriate size) <== a matrix to contain result of subtracting b-matrix from a-matrix returns [1] a-matrix minus b-matrix (placed in result)

SUM (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sum of elements in specified subsequence

SUM-OF-SQUARES (THE-SEQ &KEY (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of numbers [2] start (keyword; 0) ==> start index of sequence to be used [3] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] sum of squares of elements in specified subsequence

SUPPLIED-DATA-TAPER! (A-TIME-SERIES &KEY (TAPER-PARAMETER (HANNING-DATA-TAPER! (MAKE-ARRAY (LENGTH A-TIME-SERIES) INITIAL-ELEMENT 1.0))) (NORMALIZATION UNITY) (RESULT A-TIME-SERIES))

given [1] a-time-series (required) <=> a vector containing a time series; this vector is modified UNLESS keyword result is bound to a different vector [2] taper-parameter (keyword; vector with Hanning taper) ==> a vector containing the supplied data taper (must be the same length as a-time-series); unmodified upon return [3] normalization (keyword; :unity) ==> if :unity, taper normalized such that its sum of squares is equal to unity (Equation (208a)); if :N, taper normalized such that its sum of squares is equal to the number of points in the time series [4] result (keyword; a-time-series) <=> a vector to contain time series multiplied by the supplied data taper returns [1] a vector containing the tapered time series [2] C_h, the variance inflation factor computed using Equation (251b) in the SAPA book

SYMMETRY-PLOT (A-SEQ &KEY (A-SEQ-SORTED-P NIL))

given [1] a-seq (required) ==> any sequence of real-valued numbers [2] a-seq-sorted-p (keyword; nil) ==> if t, a-seq is already sorted; if nil, a copy of a-seq is sorted for use by the function returns [1] y values for a symmetry plot for a-seq [2] x values for a symmetry plot for a-seq

TAIL-AREA-OF-NORMAL-DISTRIBUTION (Q &KEY (Q-TO-INFINITY-P NIL))

given [1] q (required) ==> a quantile [2] q-to-infinity-p (keyword; nil) ==> if t, return integral from q to +infinity; if nil, return integral from -infinity to q return tail area of the standard norm from either q to +infinity (if q-to-infinity-p is true) or -infinity to q (if q-to-infinity-p is nil) --- see Algorithm AS-66, Applied Statistics, 1973, vol. 22, no. 3

THOMSON-WEIGHT-FUNCTION (X &KEY (BETA 1.0) (SYMMETRIC-VERSION T))

given [1] x (required) ==> value at which to evaluate Thomson's weight function [2] beta (keyword; 1.0) ==> tuning parameter [3] symmetric-version (keyword; t) ==> one-sided or two-sided down weightings returns [1] value of Thomson's weight function at x --- Note: Equation (27), p. 637 of Chave, Thomson, and Ander (1987)

THREE-POINT-SMOOTHER (A-SEQ)

given a sequence of length N, returns a sequence of length N-2 formed by filtering the input sequence with a three-point filter with coefficients 1/4, 1/2 and 1/4 --- Note: see Section 5.7 of the SAPA book

TIME-SERIES-BANDWIDTH (ACVS &KEY (SAMPLING-TIME 1.0))

given [1] acvs (required) ==> a vector of length N with values of the acvs (or acs) from lag 0 to N-1 [2] sampling-time (keyword; 1.0) ==> the sampling time returns [1] unbiased estimate of trace bandwidth (tilde B_T in Equation (280)) [2] biased estimate of trace bandwidth (hat B_T in equation just above Equation (28))) --- Note: see Section 6.14 of the SAPA book

TRACE-MATRIX (A-MATRIX)

given a square matrix, returns its trace

TRANSFER-FUNCTION-FOR-FILTER (THE-FILTER &KEY (TF-TRANSFORMATION #'(LAMBDA (X) (CAREFUL-CONVERT-TO-DB (EXPT (ABS X) 2) -100.0))) (ACTUAL-INDEX-OF-1ST-FILTER-COEFF 0) (N-FFT (* 4 (NEXT-POWER-OF-2 (LENGTH THE-FILTER)))) (RETURN-FREQUENCIES-P NIL) (NYQUIST-FREQUENCY 0.5) (RESULT-TF (MAKE-ARRAY (1+ (/ N-FFT 2)))) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (1+ (/ N-FFT 2))))))

given [1] the-filter (required) ==> a vector of filter coefficients [2] tf-transformation (keyword; mod squared in dB with 0 mapped to -100 dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all of the elements of the transfer function [3] actual-index-of-1st-filter-coeff (keyword; 0) ==> the filter coefficient in (aref the-filter 0) is assumed to have an index equal to whatever is given here (this is only needed if the phase of the transfer function is needed) [4] N-fft (keyword; 4 * (next-power-of-2 (length the-filter))) ==> the length of vector to be used in the fast Fourier transform -- the larger this is, the finer the grid of frequencies over which the transfer function is computed; this number MUST be a power of 2. [5] return-frequencies-p (keyword; nil) ==> if t, the frequencies associated with the transfer function are computed and returned in result-freq [6] Nyquist-frequency (keyword; 0.5) ==> the Nyquist frequency [7] result-tf (keyword; vector of length (1+ (/ N-fft 2))) <== vector of length (1+ (/ N-fft 2)) into which the properly transformed transfer function is placed (returned by the function) [8] result-freq (keyword; vector of length (1+ (/ N-fft 2)) if return-frequencies-p t) <== vector of length (1+ (/ N-fft 2)) into which the frequencies associated with the values in result-tf are placed if return-frequencies-p is true (returned by the function) returns [1] result-tf, a vector holding the properly transformed transfer function [2] nil (if return-frequencies-p is nil) or result-freq (if return-frequencies-p is t), where result-freq is a vector holding the frequencies associated with values in result-tf --- Note: see Section 5.3 of the SAPA book

TRANSFORM-A-SEQUENCE (A-FUNCTION THE-SEQ)

given [1] a-function (required) ==> a function with one argument [2] the-seq (required) ==> a sequence of numbers returns [1] a sequence of numbers obtained by applying a-function to each element of the-seq

TRANSFORM-A-SEQUENCE! (A-FUNCTION THE-SEQ &KEY (RESULT THE-SEQ))

given [1] a-function (required) ==> a function with one argument [2] the-seq (required) ==> a sequence of numbers [3] result (keyword; the-seq) <== a sequence of numbers returns [1] a sequence of numbers obtained by applying a-function to each element of the-seq

TRANSPOSE (A-MATRIX &KEY (RESULT (MAKE-ARRAY (REVERSE (ARRAY-DIMENSIONS A-MATRIX)))))

given [1] A (required) ==> a 2d matrix [2] result (keyword; new 2d array of appropriate size) <== a 2d matrix to contain transpose of a-matrix returns [1] transpose of a-matrix (placed in result)

TRIG-PROLATE-TAPERS (N NUMBER-OF-TAPERS &KEY (TAPER-PARAMETER 4) (PRINT-PROGRESS-P NIL) (COMPUTE-TRUE-EIGENVALUES-P NIL))

given [1] N (required) the sample size [2] number-of-tapers (required) number of orthonormal trig prolate tapers to be computed; currently restricted to one of the following maximum values: 2 if taper-parameter is 2; 4 if taper-parameter is 3; 5 if taper-parameter is 4; 7 if taper-parameter is 5 [3] taper-parameter (keyword; 4.0) NW, the duration-half-bandwidth product; currently restricted to one of the following integers: 2, 3, 4, 5 (must also be such that 0 < NW/N < 1/2) [4] print-progress-p (keyword; nil) if t, prints a dot after each taper has been computed [5] compute-true-eigenvalues-p (keyword; nil) if t, returns eigenvalues for eigenproblem of Equation (378) of the SAPA book; if nil, returns nil returns [1] a list of length number-of-tapers of N dimensional vectors with orthonormal trig prolate data tapers of orders 0, 1, ..., number-of-tapers - 1; [2] a vector of length number-of-tapers with eigenvalues if compute-true-eigenvalues-p is t; nil if if compute-true-eigenvalues-p is nil --- Note: computes the trig prolate approximation to the dpss tapers (see Section 8.4 of the SAPA book)

TRIGAMMA (X &KEY (X-RECURSION 2.0))

given [1] x (required) ==> a positive number [2] x-recursion (keyword; 2.0) ==> if x < x-recursion, recursive formula is used returns [1] value of trigamma function at x --- Note: expansion 6.4.12 plus recurrence 6.4.6 of Abramowitz and Stegun; better accuracy can be obtained by increasing x-recursion ---10.0 to 30.0 is probably a useful upper limit

TWO-STEP-BURG-ALGORITHM (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)) (APPROXIMATE-MSES (MAKE-ARRAY (1+ P))) (REFLECTION-COEFFS (MAKE-ARRAY P)) (FORWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))) (BACKWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))) (EXACT-MSES (MAKE-ARRAY (1+ P))))

given [1] time-series (required) ==> a vector containing a real-valued or complex-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an EVEN integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} [7] approximate-MSEs (keyword; a vector of length p+1) <== estimates of innovations variance (mean square errors) for models of order 0, 1, ..., p (the innovations variance for the zeroth order model is the sample variance) [8] reflection-coeffs (keyword; a vector of length p) <== estimates of reflection coefficients phi_{1,1}, ..., phi_{p,p} [9] forward-prediction-errors (keyword; a vector of length end - start - p) <== computed forward prediction errors [10] backward-prediction-errors (keyword; a vector of length end - start - p) <== computed backward prediction errors [11] exact-MSEs (keyword; a vector of length p+1) <== another set of estimates of innovations variance for models of order 0, 1, ..., p; these estimates are based on Equation (419a) of the SAPA book uses the two step Burg algorithm to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's; this number is also given in (elt approximate-MSEs p) [3] approximate-MSEs [4] reflection-coeffs [5] forward-prediction-errors [6] backward-prediction-errors [7] exact-MSEs --- Note: see Section 9.5 of the SAPA book

UPPER-TRIANGULAR-SOLVE! (A B &KEY (N (LENGTH B)))

given [1] A (required) ==> upper n by n triangular matrix (note: lower subdiagonal part of this matrix is ignored) [2] b (required) <=> on input, the right-hand side vector; on output, the solution X to A X = b [3] n (keyword; length of b) ==> order of matrix A returns [4] X, the solution to A X = b, where X is the same as vector which contained b (note that A is left unchanged)

VAR-PREDICTED-MEAN-AT-X (X COVARIANCE-MATRIX)

given [1] x (required) ==> a number [2] covariance-matrix (required) ==> covariance matrix for estimated intercept and slope as returned by weighted-linear-least-squares returns [1] the variance of the predicted mean of y at x --- Note: see page 56, Draper and Smith, 1966

VAR-PREDICTED-Y-AT-X (X RESIDUAL-VARIANCE COVARIANCE-MATRIX)

given [1] x (required) ==> a number [2] residual-variance (required) ==> the residual variance as returned by weighted-linear-least-squares [3] covariance-matrix (required) ==> covariance matrix for estimated intercept and slope as returned by weighted-linear-least-squares returns [1] the variance of the predicted value of y at x --- Note: see Equation (1.4.8), p. 24, Draper and Smith, 1966

WEIGHTED-LINEAR-LEAST-SQUARES (Y &KEY (VERSUS #'(LAMBDA (I) (FLOAT (- I (/ (1- (LENGTH Y)) 2))))) (WEIGHTS NIL) (COMPUTE-RESIDUALS-P T) (COMPUTE-COVARIANCE-MATRIX-P T) (RESIDUALS (IF COMPUTE-RESIDUALS-P (MAKE-ARRAY (LENGTH Y)))))

given [1] y (required) ==> a sequence with values of y_k's, the dependent variable. [2] versus (keyword; equally spaced & centered values) ==> a sequence or function giving values of x_k's, the independent variable [3] weights (keyword; nil) ==> if vector supplied, these are used as weights on the observations (nil implies that all y_k's are equally weighted) [4] compute-residuals-p (keyword; nil) ==> if t, residuals are computed [5] compute-covariance-matrix-p (keyword; nil) ==> if t, calculates and returns 2x2 var/covar matrix [6] residuals (keyword; nil or vector) <== a sequence to be stuffed with the residuals (not used unless compute-residuals-p is true) fits model y_k = alpha + beta * x_k, and returns [1] estimate of intercept alpha [2] estimate of slope beta [3] estimate of residual variance [4] residuals (if compute-residuals-p is true) [5] covariance matrix for parameter estimates (if compute-covariance-matrix-p is true) --- Note: see Chapter 2 of Draper and Smith, 1966

WEIGHTED-MEAN (THE-SEQ &KEY (WEIGHTS NIL) (START 0) (END (LENGTH THE-SEQ)))

given [1] the-seq (required) ==> a sequence of real or complex-valued numbers [2] weights (keyword; nil) ==> a sequence with (- end start) values to be used as weights; if set to nil, equal weights for all numbers is assumed [3] start (keyword; 0) ==> start index of sequence to be used [4] end (keyword; length of the-seq) ==> 1 + end index of sequence to be used returns [1] the weighted average of the specified numbers in the-seq

WINDOW-WIDTH-FROM-SILVERMAN (A-SEQ)

given a-seq, returns the window width for kernel pdf estimation given on page 48, Equation (3.31), of Silverman, 1986

WOSA-SPECTRAL-ESTIMATE (TIME-SERIES BLOCK-SIZE &KEY (PROPORTION-OF-OVERLAP 0.5) (OVERSAMPLING-FACTOR 1) (CENTER-DATA T) (START 0) (END (LENGTH TIME-SERIES)) (RETURN-EST-FOR-0-FREQ-P T) (SAMPLING-TIME 1.0) (SCRATCH-DFT (MAKE-ARRAY (* OVERSAMPLING-FACTOR BLOCK-SIZE))) (DATA-TAPER #'HANNING-DATA-TAPER!) (DATA-TAPER-PARAMETERS NIL) (RESTORE-POWER-OPTION-P T) (SDF-TRANSFORMATION #'CONVERT-TO-DB) (RESULT-SDF (MAKE-ARRAY (WOSA-GET-N-FREQS BLOCK-SIZE OVERSAMPLING-FACTOR RETURN-EST-FOR-0-FREQ-P))) (RETURN-SDF-ESTIMATES-FOR-EACH-BLOCK-P T) (RETURN-FREQUENCIES-P T) (FREQ-TRANSFORMATION NIL) (RESULT-FREQ (IF RETURN-FREQUENCIES-P (MAKE-ARRAY (WOSA-GET-N-FREQS BLOCK-SIZE OVERSAMPLING-FACTOR RETURN-EST-FOR-0-FREQ-P)))))

given [1] time-series (required) ==> a vector of real-valued numbers [2] block-size (required) ==> a power of two [3] proportion-of-overlap (keyword; 0.5) ==> number greater than 0 and less than 1 [4] oversampling-factor (keyword; 1) ==> a factor that controls the number of frequencies at which the wosa spectral estimate is computed; this factor should be an integer power of two such as 1, 2, 4, etc; for example, 1 yields Fourier frequencies for block-size; 2 yields grid twice as fine as Fourier frequencies; 4 yields griid 4 times as fine as Fourier frequencies; etc. [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, time-series is not centered [6] start (keyword; 0) ==> start index of vector to be used [7] end (keyword; length of time-series) ==> 1 + end index of vector to be used [8] return-est-for-0-freq-p (keyword; nil) ==> if t, sdf is computed at zero frequency; otherwise, it is not computed. [9] sampling-time (keyword; 1.0) ==> sampling time (called delta t in the SAPA book) [10] scratch-dft (keyword; vector of correct length) ==> vector in which the in-place dft is done [11] data-taper (keyword; #'Hanning-data-taper!) ==> a tapering function or nil [12] data-taper-parameters (keyword; nil) ==> parameters for tapering function (not used if data-taper is nil); the default of nil is appropriate for the Hanning data taper because it does not have any parameters [13] restore-power-option-p (keyword; t) ==> if t and data-taper is non-nil, normalizes tapered series to have same sum of squares as before tapering [14] sdf-transformation (keyword; #'convert-to-dB) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-sdf [15] result-sdf (keyword; vector of correct length) <== vector into which wosa sdf estimate is placed; it must be EXACTLY of the length dictated by block-size, oversampling-factor and return-est-for-0-freq-p [16] return-sdf-estimates-for-each-block-p (keyword; t) ==> if t, individual spectra for each block are returned in a list; note that these spectra are untransformed (i.e., the option sdf-transformation applies only to the final wosa estimate) [17] return-frequencies-p (keyword; t) ==> if t, the frequencies associated with the spectral estimate are computed and returned in result-freq [18] freq-transformation (keyword; nil) ==> a function of one argument or nil; if bound to a function, the function is used to transform all elements of result-freq (ignored unless return-frequencies-p is true) [19] result-freq (keyword; nil or vector of correct length) <== not used if return-frequencies-p nil; otherwise, vector of length dictated by block-size, oversampling-factor and return-est-for-0-freq-p into which the frequencies associated with the values in result-sdf are placed returns [1] wosa spectral estimate [2] associated frequencies [3] equivalent degrees of freedom [4] list of individual direct spectral estimates

X*Y (X Y)

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers returns [1] a vector of numbers obtained by multiplying corresponding elements of x and y

X*Y! (X Y &KEY (RESULT X))

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers [3] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by multiplying corresponding elements of x and y

X+B (X B)

given [1] x (required) ==> a sequence of numbers [2] b (required) ==> a number returns [1] a vector of numbers obtained by adding b to each element of x

X+B! (X B &KEY (RESULT X))

given [1] x (required) ==> a sequence of numbers [2] b (required) ==> a number [3] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by adding b to each element of x

X+Y (X Y)

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers returns [1] a vector of numbers obtained by adding x and y together element by element

X+Y! (X Y &KEY (RESULT X))

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers [3] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by adding x and y together element by element

X-Y (X Y)

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers returns [1] a vector of numbers obtained by subtracting elements of y from corresponding elements of x

X-Y! (X Y &KEY (RESULT X))

given [1] x (required) ==> a sequence of numbers [2] y (required) ==> a sequence of numbers [3] result (keyword; x) <== a sequence of numbers returns [1] a sequence of numbers obtained by subtracting elements of y from corresponding elements of x

YULE-WALKER-ALGORITHM-GIVEN-ACVS (ACVS P &KEY (AR-COEFFS (MAKE-ARRAY P)) (APPROXIMATE-MSES (MAKE-ARRAY (1+ P))) (REFLECTION-COEFFS (MAKE-ARRAY P)))

given [1] acvs (required) ==> a vector containing values of the acvs for a real-valued time series x_t from lag 0 up to (and including) p [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} [4] approximate-MSEs (keyword; a vector of length p+1) <== estimates of innovations variance (mean square errors) for models of order 0, 1, ..., p (the innovations variance for the zeroth order model is the sample variance) [5] reflection-coeffs (keyword; a vector of length p) <== estimates of reflection coefficients phi_{1,1}, ..., phi_{p,p} uses the Yule-Walker method to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's; this number is also given in (elt approximate-MSEs p) [3] approximate-MSEs [4] reflection-coeffs --- Note: see Sections 9.3 and 9.4 of the SAPA book

YULE-WALKER-ALGORITHM-GIVEN-DATA (TIME-SERIES P &KEY (START 0) (END (LENGTH TIME-SERIES)) (CENTER-DATA T) (AR-COEFFS (MAKE-ARRAY P)) (APPROXIMATE-MSES (MAKE-ARRAY (1+ P))) (REFLECTION-COEFFS (MAKE-ARRAY P)) (FORWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))) (BACKWARD-PREDICTION-ERRORS (MAKE-ARRAY (- END START P))))

given [1] time-series (required) ==> a vector containing a real-valued time series x_t [2] p (required) ==> autoregressive model order; p should be an integer > 0 and < end - start [3] start (keyword; 0) ==> start index of vector to be used [4] end (keyword; length of the-seq) ==> 1 + end index of vector to be used [5] center-data (keyword; t) ==> if t, subtract sample mean from time series; if a number, subtracts that number from time series; if nil, the effect of this function is to return a copy of the relevant portion of time-series [6] AR-coeffs (keyword; a vector of length p) <== estimates of AR coefficients phi_{1,p}, ..., phi_{p,p} [7] approximate-MSEs (keyword; a vector of length p+1) <== estimates of innovations variance (mean square errors) for models of order 0, 1, ..., p (the innovations variance for the zeroth order model is the sample variance) [8] reflection-coeffs (keyword; a vector of length p) <== estimates of reflection coefficients phi_{1,1}, ..., phi_{p,p} [9] forward-prediction-errors (keyword; a vector of length end - start - p) <== computed forward prediction errors [10] backward-prediction-errors (keyword; a vector of length end - start - p) <== computed backward prediction errors uses the Yule-Walker method to estimate the phi's in the model x_t = phi_{1,p} x_{t-1} + ... + phi_{p,p} x_{t-p} + e_t and returns [1] AR-coeffs, a vector containing phi_{1,p}, ..., phi_{p,p} [2] estimate of the innovations variance, i.e., the variance of the e_t's; this number is also given in (elt approximate-MSEs p) [3] approximate-MSEs [4] reflection-coeffs [5] forward-prediction-errors [6] backward-prediction-errors --- Note: see p. 420 of the SAPA book

ZERO-STRICT-LOWER-DIAGONAL! (A-MATRIX)

given a square matrix, zeros its lower diagonal and returns the modified square matrix

ZEROS-OF-POLYNOMIAL (POLYNOMIAL &KEY (DEGREE-OF-POLYNOMIAL (1- (LENGTH POLYNOMIAL))) (THE-ROOTS (MAKE-ARRAY DEGREE-OF-POLYNOMIAL)) (MAXIMUM-NUMBER-OF-ITERATIONS 25))

given [1] polynomial (required) ==> sequence with coefficients of a polynomial; (elt polynomial 0) must be nonzero [2] degree-of-polynomial (keyword; (1- (length polynomial))) ==> degree of polynomial [3] the-roots (keyword; array of length degree-of-polynomial) <== number of derivatives to be evaluated [4] maximum-number-of-iterations (keyword; 25) ==> maximum number of iterations returns [1] t or nil, where t indicates that all went well, whereas nil indicates that convergence did not occur at end of specificed number of iterations [2] the-roots, a vector with the required roots of the polynomial

Undocumented

TRIANGULAR-CONVERGENCE-FACTORS (K FILTER-LENGTH)

Private

AS-111 (P)

adapted from AS 111 --- FORTRAN routine ppnd

DPSS->EIGENVALUE (DPSS NW)

given dpss (a vector of length N) and NW, computes the corresponding eigenvalue using the method of Exercise [8.1], page 390, of the SAPA book

EIGENVALUES-AND-VECTORS-OF-SYM-TRIDIAG-MATRIX (DIAG OFF-DIAG &KEY (MAXIMUM-NUMBER-OF-ITERATIONS 30) (RETURN-EIGENVECTORS-P T) (Z (IF RETURN-EIGENVECTORS-P (LET* ((N (LENGTH DIAG)) (TEMP (MAKE-ARRAY `(,N ,N) INITIAL-ELEMENT 0.0))) (DOTIMES (I N TEMP) (SETF (AREF TEMP I I) 1.0))))))

given [1] diag (required) ==> sequence of diagonal elements (length of, say, n) [2] off-diag (required) ==> sequence of off-diagonal elements (length of n-1) [3] maximum-number-of-iterations (keyword; 30) ==> maximum number of iterations [4] return-eigenvectors-p (keyword; t) ==> if true, returns eigenvalues and eigenvectors; if nil, just do the eigenvalues (this is faster) [5] z (keyword; n by n diagonal matrix) ==> usually bound to output from Householder-reduction-of-real-sym-matrix returns [1] vector with n eigenvalues (NOT necessarily ordered by size) [2] if return-eigenvectors-p is true, an n by n array whose jth column is the eigenvector corresponding to the jth eigenvalue; if return-eigenvectors-p is nil, that's what you get! --- This is based upon tqli in Numerical Recipes

EQUIVALENT-DOF-FOR-WOSA (N N_S N_B VECTOR-WITH-DATA-TAPER)

given the sample size N, the block size N_S, the number of blocks N_B, and the values of the data taper, returns the equivalent degrees of freedom for wosa using Equation (292b) of the SAPA book

EQUIVALENT-DOF-FOR-WOSA-STANDARD-CASE (N_S)

given number of blocks N_S, returns the equivalent degrees of freedom for wosa with 50% overlap and Hanning data taper using Equation (294) of the SAPA book

HOUSEHOLDER-REDUCTION-OF-REAL-SYM-MATRIX! (REAL-SYM-MATRIX &KEY (PREPARE-FOR-EIGENVECTORS-P T))

given [1] real-sym-matrix (required) ==> real, symmetric matrix (n by n) [2] prepare-for-eigenvectors-p (keyword; t) ==> if true, crunch real-sym-matrix for use with eigenvalues-and-vectors-of-sym-tridiag-matrix returns [1] n-dimensional vector with diagonal elements of associated tridiagonal matrix [2] (n-1)-dimensional vector with off-diagonal elements of associated tridiagonal matrix --- See tred2 in Numerical Recipes

INCOMPLETE-GAMMA (A X)

arguments: a --- parameter value > 0 x --- variable >= 0 returns three values: value of incomplete gamma function at a and x; the numerator (little gamma of a and x); the denominator (gamma of a); see Section 6.2 of NR

INCOMPLETE-GAMMA-Q (A X)

arguments: a --- parameter value > 0 x --- variable >= 0 returns three values: value of incomplete gamma function at a and x; the numerator (big gamma of a and x); the denominator (gamma of a); see Section 6.2 of NR

SYMMETRIC-TRIDIAGONAL-SOLVE! (DIAG OFF-DIAG B)

given [1] diag (required) <=> diagonal part of symmetric tridiagonal matrix A; trash on output [1] off-diag (required) <=> off-diagonal part of symmetric tridiagonal matrix; trash on output [2] b (required) <=> on input, the right-hand side vector; on output, the solution X to A X = b returns [4] X, the solution to A X = b, where X is the vector that contained b on input --- Note: this is an implementation of Algorithm 4.3.6, p. 156, Golub and Van Loan, 1989, with modifications to avoid divides by zero

TOEPLITZ (R Y &OPTIONAL (X (MAKE-ARRAY (LENGTH Y))))

given [1] r (required) ==> vector of length 2*N-1 representing the Toeplitz matrix R(i,j) with 0 <= i,j <= N-1, via the 2*N-1 values r(0,N-1), r(0,N-2),..., r(0,1), r(0,0), r(1,0),..., r(N-1,0) [2] y (required) ==> right-hand side vector of length N [3] x (optional; vector of same size as y) <== solution vector of length N solves the Toeplitz system sum_{j=0}^{N-1} r_{N-1+i-j} x_j = y_j for i = 0, 1, ..., N-1 and returns [1] x (a vector of length N) or nil (in case of a failure of the Levinson method due to a singular principal minor --- Note: this function is essentially a Lisp version of the Fortran routine toeplz from Numerical Recipes

Undocumented

ABSOLUTE-SUM-OF-VECTOR-ELEMENTS (VECTOR)

AS-91 (P DEGREES-OF-FREEDOM)

AS-91-BLOCK-1 (CH A G C)

AS-91-BLOCK-4 (CH XX C P G)

BESSI0-NR (X)

BINARY-SEARCH-INTERNAL (VALUE ORDERED-SEQ START END LOWER-TEST UPPER-TEST)

CALCULATE-NUMBER-OF-BLOCKS (SAMPLE-SIZE BLOCK-SIZE PROPORTION-OF-OVERLAP)

COMPLEX-OF-ABSOLUTES (Z)

CONTINUED-FRACTION-REPRESENTATION-FOR-INCOMPLETE-GAMMA (A X)

EVAL-POLY-2 (POLYNOMIAL X)

FAST-TRIDIAG-EIGENVALUE->DPSS! (EIGENVALUE ORDER-OF-TAPER DIAG OFF-DIAG &KEY (EPS (* 10 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 25) (B (GENERATE-INITIAL-GUESS-AT-DPSS (LENGTH DIAG) ORDER-OF-TAPER HALF-SIZE-P T)))

FIND-SYMMETRIC-FREQUENCY (F_1 F_2 F_3 F_N)

FLIP-SIGN-OF-COEFFICIENTS! (A-SEQUENCE)

GENERATE-INITIAL-GUESS-AT-DPSS (N ORDER-OF-TAPER &KEY (HALF-SIZE-P NIL) (RESULT (MAKE-ARRAY (IF HALF-SIZE-P (/ (IF (EVENP N) N (1- N)) 2) N))))

GENERATE-TRI-PROLATE-TAPER (N GENERATING-FUNCTION)

GET-N-DFT (N-NONZERO-FREQS SAMPLE-SIZE)

GET-N-FREQS (N-NONZERO-FREQS SAMPLE-SIZE RETURN-EST-FOR-0-FREQ-P)

GET-OFFSET-TO-KTH-BLOCK (SAMPLE-SIZE BLOCK-SIZE NUMBER-OF-BLOCKS K)

H-FUNC (F G X)

INITIAL-ESTIMATES-FOR-ZEROS-OF-POLYNOMIAL (POLYNOMIAL &KEY (DEGREE-OF-POLYNOMIAL (1- (LENGTH POLYNOMIAL))) (INITIAL-ESTIMATES (MAKE-ARRAY DEGREE-OF-POLYNOMIAL)) (SCRATCH (MAKE-ARRAY (* 2 (1+ DEGREE-OF-POLYNOMIAL)))))

KTH-EIGENVALUE-OF-TRIDIAGONAL-MATRIX (DIAG OFF-DIAG K &KEY (SQUARED-OFF-DIAG (MAP 'VECTOR #'(LAMBDA (X) (* X X)) OFF-DIAG)) (MACHEPS SINGLE-FLOAT-EPSILON))

LARGEST-EIGENVALUES-OF-TRIDIAGONAL-MATRIX (DIAG OFF-DIAG NUMBER-OF-TAPERS &KEY (SQUARED-OFF-DIAG (MAP 'VECTOR #'(LAMBDA (X) (* X X)) OFF-DIAG)) (MACHEPS SINGLE-FLOAT-EPSILON) (PRINT-PROGRESS-P NIL))

ONE-ITERATION-OF-TWO-STEP-BURG (K STEP N FK BK PHI LIMIT)

POWER-OF-2+1 (N)

PRE-FFT (N &KEY (COMPLEX-EXP-VECTOR NIL))

PRE-FFT-WITH-CACHE (N)

PRODUCE-UPSILON-FUNC (COEFFICIENTS)

SERIES-REPRESENTATION-FOR-INCOMPLETE-GAMMA (A X)

STURM-SEQUENCE-COUNT (TEST-LAMBDA DIAG OFF-DIAG SQUARED-OFF-DIAG MACHEP &KEY (START 0) (END (LENGTH DIAG)))

SYMMETRIC-TRIDIAGONAL-INFINITY-NORM (DIAG OFF-DIAG)

TRAPEZOIDAL-RULE (F A B S N)

TRIDIAG-EIGENVALUE->DPSS! (EIGENVALUE ORDER-OF-TAPER DIAG OFF-DIAG &KEY (EPS (* 10 SINGLE-FLOAT-EPSILON)) (MAXIMUM-NUMBER-OF-ITERATIONS 25) (B (GENERATE-INITIAL-GUESS-AT-DPSS (LENGTH DIAG) ORDER-OF-TAPER)))

TRIG-PROLATE-SIGN-COSINE-COEFFICIENTS (M K)

WHICH-IS-FASTER? (N-FILTER N-TS)

WOSA-GET-N-FREQS (BLOCK-SIZE OVERSAMPLING-FACTOR RETURN-EST-FOR-0-FREQ-P)

ZEROS-OF-POLYNOMIAL-WITH-INITIAL-ESTIMATES (POLYNOMIAL INITIAL-ESTIMATES &KEY (DEGREE-OF-POLYNOMIAL (1- (LENGTH POLYNOMIAL))) (SCRATCH (MAKE-ARRAY (1+ DEGREE-OF-POLYNOMIAL))) (FINAL-RESULTS INITIAL-ESTIMATES) (MAXIMUM-NUMBER-OF-ITERATIONS 25))

MACRO

Public

DIVF (PLACE &OPTIONAL (DELTA 1) &ENVIRONMENT ENV)

Like <incf>

MULTF (PLACE &OPTIONAL (DELTA 1) &ENVIRONMENT ENV)

Like <incf>

VARIABLE

Private

Undocumented

*ABSCISSAS-32-POINT*

*SAPA-CACHED-PRE-FFT-ARRAYS*

*WEIGHTS-32-POINT*

+COEFFICENTS-FOR-LOG-OF-GAMMA+

+NR-GSER-EPS+

+NR-GSER-ITMAX+

+USE-PRE-FFT-WITH-CACHE-P+

CONSTANT

Private

Undocumented

+EULER-CONSTANT+

+LOG-PI+

+PPCHI2-AA+

+PPCHI2-C1+

+PPCHI2-C10+

+PPCHI2-C11+

+PPCHI2-C12+

+PPCHI2-C13+

+PPCHI2-C14+

+PPCHI2-C15+

+PPCHI2-C16+

+PPCHI2-C17+

+PPCHI2-C18+

+PPCHI2-C19+

+PPCHI2-C2+

+PPCHI2-C20+

+PPCHI2-C21+

+PPCHI2-C22+

+PPCHI2-C23+

+PPCHI2-C24+

+PPCHI2-C25+

+PPCHI2-C26+

+PPCHI2-C27+

+PPCHI2-C28+

+PPCHI2-C29+

+PPCHI2-C3+

+PPCHI2-C30+

+PPCHI2-C31+

+PPCHI2-C32+

+PPCHI2-C33+

+PPCHI2-C34+

+PPCHI2-C35+

+PPCHI2-C36+

+PPCHI2-C37+

+PPCHI2-C38+

+PPCHI2-C4+

+PPCHI2-C5+

+PPCHI2-C6+

+PPCHI2-C7+

+PPCHI2-C8+

+PPCHI2-C9+

+PPCHI2-E+

+PPCHI2-PMAX+

+PPCHI2-PMIN+

+PPND-A0+

+PPND-A1+

+PPND-A2+

+PPND-A3+

+PPND-B1+

+PPND-B2+

+PPND-B3+

+PPND-B4+

+PPND-C0+

+PPND-C1+

+PPND-C2+

+PPND-C3+

+PPND-D1+

+PPND-D2+

+SAPA-10-OVER-LOG-10+

+SAPA-2-PI+

+SAPA-4-TIME-EXP-OF-1-OVER-4+

+SAPA-4-TIME-EXP-OF-MINUS-1-POINT-35+

+SAPA-MINUS-2-PI+

+SAPA-SQRT-8-OVER-E+

+SQRT-OF-TWO-PI+