# 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

# BINARY-SEARCH (VALUE ORDERED-SEQ &KEY (LOWER-TEST #'<=) (UPPER-TEST #'<=))

given
[1] value (required)
==> a real-valued number
[2] ordered-seq (required)
==> an ordered sequence of real-valued numbers
[3] lower-test (keyword; #'<=)
==> a predicate of two arguments
that defines the lower test
[4] upper-test (keyword; #'<=)
==> a predicate of two arguments
that defines the upper test
returns
[1] if lower-test and upper-test are set
to their default values, this function
returns an index i such that v[i] <= value <= v[i+1],
where 0 <= i <= n-2 (if there is no such i, nil is returned);
if nondefault values are supplied for lower-test and upper-test,
then the condition `v[i] <= value <= v[i+1]'
gets modified in an obvious way
---
Note: value can be an arbitrary object, and ordered-seq
can be a list of arbitrary objects if the binary predicates
lower-test and upper-test are appropriate set

# 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

# PRINT-1-OR-2-D-ARRAY (AN-ARRAY &KEY (TAG array) (FORMAT-CONTROL-FOR-ARRAY-ELEMENT ~F))

given
[1] an-array (required)
==> a one-dimensional or two-dimensional array
[2] tag (keyword; `array')
==> an optional tag to be printed
along with array elements
[3] format-control-for-array-element (keyword; `~F')
==> format control for a single array element
prints the elements of an-array and
returns
[1] an-array

# PRINT-ROWS-OF-NXM-MATRIX (AN-NXM-MATRIX &KEY (FORMAT-CONTROL-FOR-ARRAY-ELEMENT ~F ))

given
[1] an-nxm-matrix (required)
==> a two-dimensional array
[3] format-control-for-array-element (keyword; `~F')
==> format control for a single array element
prints the elements of the 2d array and
returns
[1] an-nxm-matrix

# 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>