### Computation of the Frequency Function of a Quadratic Form in Random Normal VariablesComputation of the Frequency Function of a Quadratic Form in Random Normal Variables Access Restriction
Subscribed

 Author Freiberger, Walter F. ♦ Jones, Richard H. Source ACM Digital Library Content type Text Publisher Association for Computing Machinery (ACM) File Format PDF Copyright Year ©1960 Language English
 Subject Domain (in DDC) Computer science, information & general works ♦ Data processing & computer science Abstract One of the outstanding problems in the theory of time series analysis is the distribution problem in spectral analysis for small samples. When the number of observations is sufficiently large for the Central Limit Theorem to be applicable, the normal approximation can be used to advantage. In recent years, work has aimed at the discovery of more generally applicable approximate methods and of more rational criteria for the sample size at which the large sample theory becomes useful; the state of the art is summarized in two recent papers by Grenander, Pollak and Slepian  and Freiberger and Grenander .Consider a sample $\textit{x}$ = $(\textit{x}1,$ $\textit{x}2,$ ··· , $\textit{xn})$ of successive values taken from a discrete-valued stationary time series ··· , $\textit{y}-1$ , $\textit{y}0$ , $\textit{y}1$ , ··· , of normally distributed random variables with mean zero and covariance matrix $\textit{R}$ with elements $\textit{r&ngr;μ}$ = $\textit{E}(\textit{x&ngr;xμ}).$ For stationarity we have: $\textit{ri,i}+\textit{&ngr;}≡$ $\textit{r&ngr;}$ = $\textit{E}(\textit{yiyi}+\textit{&ngr;}),$ i, &ngr; = 0, ±1, ±2, ···. The Fourier-Stieltjes representation of the covariances is $\textit{r&ngr;}$ = 1/2π ∫ππ $\textit{e}\textit{i}λ\textit{&ngr;}$ $\textit{dF}(λ)$ (1) where the spectral distribution function $\textit{F}(\textit{λ})$ is bounded and nondecreasing and can, if it is absolutely continuous, be expressed in terms of a spectral density ƒ(λ): $\textit{F}(λ)$ = ∫λπƒ(λ) (1) The practical determination of ƒ(λ) from observations of the process is effected by the introduction of a quadratic form $\textit{Q}$ = $∑\textit{n}\textit{&ngr;},λ=1\textit{w}\textit{&ngr;}-\textit{μ}\textit{x&ngr;xμ}$ which is taken as an estimate for ƒ(λ); the coefficients $\textit{w&ngr;}$ can be written in terms of a spectral weight function or “spectral window” $\textit{w}(λ):$ $\textit{w&ngr;}$ = 1/2π ∫π-π $\textit{x}(λ)\textit{e}\textit{i&ngr;}λ$ $\textit{d}λ.$ (4) For details, see the cited references. In order to obtain confidence limits for these estimates $\textit{Q}$ of the spectral density ƒ(λ), it is important to have methods for computing the distribution of $\textit{Q}.$ This paper deals with one such method which has proved very efficient for digital computer application.It is well known and shown, for instance, by H. Cramèr [3, p. 118], that the characteristic function of (3), where $\textit{w&ngr;μ}$ ≡ $\textit{w&ngr;}-\textit{μ}$ are elements of a non-negative definite symmetric matrix $\textit{W},$ is given by $\textit{φ}(\textit{z})$ = $\textit{EeizQ}$ = | $\textit{I}$ - $2\textit{izRW}$ | -1/2 = $∏n\textit{j}=1$ (1 - $2\textit{iλjz})-1/2$ (5) where the $λ\textit{j}$ are the $\textit{n}$ eigenvalues of the matrix product $\textit{RW}.$ Although $\textit{RW}$ is not necessarily a symmetric matrix, both $\textit{R}$ and $\textit{W}$ are symmetric and non-negative definite, so that all the $λ\textit{j}$ are real and non-negative. The frequency function $\textit{g}(\textit{x})$ of $\textit{Q}$ then follows from Fourier's inversion formula: $\textit{g}(\textit{x})$ = 1/2π ∫∞∞ $\textit{e}-\textit{ixz}$ $∫\textit{n}\textit{j}=1$ (1 - $2\textit{iz}λ\textit{j})-1/2$ $\textit{dz}.$ (6) Several ways have been suggested for evaluating $\textit{g}$ $(\textit{x}).Slepian$  obtains a sum of finite integrals by deforming the contour of integration into a set of circles enclosing pairs of the branch points $\textit{zj}$ = $-\textit{i}/(2λ\textit{j})$ of (6), and collapsing the circles. This method, which was also used in , works well when the eigenvalues cluster toward zero, but not otherwise.A method of repeated convolution of the frequency functions of its individual terms to obtain the frequency function of a quadratic form was developed in  and programmed for the IBM 650. It was found to be slowly convergent in most cases.Taking the logarithmic derivative of (5) and applying the inverse Fourier transform yields the following singular integral equation for $\textit{g}(\textit{x}):$ $\textit{xg}(\textit{x})$ = ∫x0 $g}(\textit{x}$ - $\textit{y})\textit{h}(\textit{y})$ $\textit{dy}$ (7) where $\textit{h}(\textit{x})$ = 1/2 $∑\textit{n}\textit{&ngr;}=1$ $\textit{e}-\textit{i}/2λ\textit{&ngr;}.$ This observation forms the basis of an ingeneous method  for a computational scheme to derive the distribution of $\textit{Q}.$ The difficulties of obtaining initial values, associated with the high order zero of the frequency function $\textit{g}(\textit{x})$ at $\textit{x}$ = 0, are solved and the solution of (7) is discussed in detail in . This method is suitable only for large-scale computers of the order of the IBM 704, and becomes, as does the method in , slowly convergent when the eigenvalues are densely spaced. Gurland  expanded the frequency function in terms of a Laguerre series and () presented a convergent series for the distribution fucntion using Laguerre polynomials. This method, which was found to converge slowly, formed the basis for the following procedure which has proved fast, accurate and reliable for a variety of problems. Essentially, the Laguerre expansion is now taken around Rice's approximation [8, p. 99], which is a type III distribution with appropriately chosen parameters.The following definition of the Laguerre polynomials $\textit{Ln}(α)$ $(\textit{x})$ = $∑\textit{n}\textit{&ngr;}=0$ $(\textit{n}+\textit{α}\textit{n}-\textit{&ngr;})(-\textit{x)&ngr;}/\textit{&ngr;}!$ (8) satisfies the orthogonality relations $∫∞0\textit{e}-\textit{x}$ $\textit{x}\textit{α}\textit{Ln}(\textit{α})$ $(\textit{x})$ $\textit{Lm}(\textit{α})$ $(\textit{x})$ $\textit{dx}$ = ${Γ(\textit{α}$ $=1)(\textit{n}+\textit{αn});$ $\textit{m}$ = $\textit{n}\textit{m}$ ≠ $\textit{n}$ (See Szegö ). Replacing $\textit{x}$ by $λ\textit{x}$ gives ∫∞0 $\textit{e}-λ\textit{x}\textit{x}\textit{α}\textit{L}(\textit{α})\textit{n}(Γ\textit{α})$ $\textit{dx}={Γ(\textit{α}+1/λ\textit{α}+1)$ $(\textit{n}+\textit{αn});\textit{m}$ = $\textit{n}\textit{m}$ ≠ $\textit{n}$ (9) which has the desired weight. The frequency function may now be expanded in a modified Laguerre series $\textit{g}(\textit{x})$ = $\textit{Kx}\textit{α}\textit{e}-λ\textit{x}[\textit{c}0$ + $\textit{c}1\textit{L}1(\textit{α})$ $(λ\textit{x})$ + $\textit{c}2\textit{L}2(α)$ $(λ\textit{x})$ + …] (10) where $\textit{K}$ is chosen so that the weight integrates to one: $∫∞0\textit{Kx}α\textit{e}-λ\textit{x}$ $\textit{dx}$ = 1, ∴ $\textit{K}$ = λα+1/Γ(α + 1). (11) Multiplying both sides of (10) by $\textit{Ln}(α)$ $(λ\textit{x})$ and integrating from 0 to ∞, gives $\textit{cn}$ = $1/(\textit{n}$ + $\textit{αn})$ $∫∞0\textit{L}(\textit{α})\textit{n}(λ\textit{x})\textit{g}(\textit{x})$ $\textit{dx}.$ (12) Using (8), $\textit{cn}$ = $∑\textit{n}\textit{&ngr;}=0(\textit{n}$ + $\textit{α}\textit{n}$ - $\textit{&ngr;})(\textit{n}$ + $\textit{αn}(-λ)\textit{&ngr;}\textit{&ngr;}!∫∞0\textit{x&ngr;9}(\textit{x})\textit{dx}$ $\textit{cn}$ = $∑∞\textit{&ngr;}=0Γ(\textit{α}+1)/Γ(\textit{α}+\textit{&ngr;}+1(\textit{n&ngr;})(-λ)\textit{&ngr;}\textit{α&ngr;}$ (13) where $α\textit{&ngr;}$ is the $\textit{&ngr;}th$ moment of the distribution about zero. Taking the logarithm of the characteristic function (5), log $[\textit{EeizQ}]$ = - 1/2 $∑\textit{n}\textit{j}=1$ log (1 - $2\textit{i}λ\textit{j}\textit{z}),$ and expanding in powers of $\textit{iz},$ one obtains for the cumulants of the distribution (Cramèr [3, p. 186]) $\textit{Xn}$ = $(\textit{n}$ - $1)!2\textit{n}-1$ $∑\textit{j}λ\textit{j}\textit{n}.$ (14) The cumulants are related to the moments about zero by the relation 1 + $\textit{α}1(\textit{iz})$ + $\textit{α}2/2!(\textit{iz})2$ + $\textit{α}3/3!(\textit{iz})3$ + ··· = exp [ $\textit{&khgr;}1(\textit{iz})$ + $\textit{&khgr;}2/2!$ $(\textit{iz})2$ + ··· ]. (15) The mean and standard deviation of the frequency function are arrived at by equating powers of $\textit{iz}:$ $\textit{m}$ = &khgr;1; $\textit{&sgr;}2$ = &khgr;2. (16) The weight function will have the same mean and standard deviation if $\textit{α}$ = $\textit{m}2/\textit{&sgr;}2$ - 1; λ = $\textit{m}/\textit{&sgr;}2.$ (17) Using (13) the remembering that $\textit{α}1$ = $\textit{m}$ and $\textit{α}2$ = $\textit{&sgr;}2$ + $\textit{m}2,$ one obtains $\textit{c}0$ = 1, $\textit{c}1$ = 0, $\textit{c}2$ = 0. (18) Now (10) becomes $\textit{g}(\textit{x})$ = $λ\textit{α}+1/Γ(\textit{α}$ + $1)\textit{x}\textit{α}\textit{e}-λ\textit{x}$ [1 + $\textit{c}3\textit{L}(\textit{α})3(λ\textit{x})$ + $\textit{c}4\textit{L}(\textit{α})4(λ\textit{x})$ + ···] (19) which together with (17), (15), (14), (13) and (8), gives the frequency function in terms of the eigenvalues.Szegö  showed that a sufficient condition for convergence, when expanding in the form $\textit{g}(\textit{x})$ = $∑∞\textit{&ngr;}$ = 0 $\textit{a}&ngr;\textit{L}(\textit{α})&ngr;$ $(\textit{x}),$ is $\textit{g}(\textit{x})$ = $\textit{O}(\textit{e}\textit{x}/2\textit{x}-\textit{α}/2-1/4-\textit{δ}),$ $\textit{δ}$ > 0, $\textit{x}$ → ∞. In our case, this reduces to $ƒ(\textit{x})$ = $\textit{O}(\textit{e}-\textit{x}/2\textit{x}\textit{α}/2-1/4-\textit{δ}),$ $\textit{δ}$ > 0, $\textit{x}$ → ∞. But, (Gurland ) $ƒ(\textit{x})$ = $\textit{O}(\textit{e}-\textit{x}/2λ\textit{m}\textit{x}\textit{n}/2-1,$ $λ\textit{m}$ = max $λ\textit{j},$ $\textit{x}$ → ∞. Therefore, the series will converge if $λλ\textit{m}$ < 1, which is the same restriction imposed by Gurland.For equal eigenvalues, the result is a $\textit{&khgr;}2-distribution$ with $\textit{n}$ degrees of freedom, which serves as an estimate of the computational error.1It may be remarked that Toeplitz theory  gives the asymptotic distribution of the eigenvalues $λ\textit{&ngr;}$ and for large enough sample sizes (see ) it is sufficient to use these results instead of the exact eigenvalue distribution. In these cases, the present method is particularly efficient.The method of computation proceeds as follows. The desired number of cumulants divided by $\textit{n}!$ is calculated from equation (14), to give $&khgr;\textit{n}/\textit{n}!$ = $2\textit{n}-1/\textit{n}$ $∑\textit{j}λ\textit{j}\textit{n};$ after substitution in equation (15), powers of $\textit{iz}$ are compared to give the moments $\textit{α}\textit{&ngr;}.$ The constant $\textit{K}$ follows from (11), with the gamma function calculated by fitting an eighth order polynomial to Γ $(\textit{α}$ + 1) in the range 0 ≦ $\textit{α}$ ≦ 1 and using the recurrence relation Γ $(\textit{α}$ + 1) = $\textit{α}&Ggr(\textit{&agr;α})$ to extend the range to any $\textit{α}$ > - 1. Next, $\textit{α}$ and λ are computed and the series tested for convergence. In all applications to spectral analysis tried so far, the product $λλ\textit{m}$ has been between 0.5 (for the case of equal eigenvalues) and 0.75, so that there has been no trouble with convergence. Successively higher approximations are now calculated, using intermediate quantities $\textit{Jn}$ = $Γ(\textit{α}$ + $1)/Γ(\textit{α}$ = $\textit{n}$ + 1) $(-λ)\textit{n}$ = $(-λ)\textit{n}/\textit{α}$ + $1)(\textit{α}$ + 2) ··· $(\textit{α}$ + $\textit{n})′$ $\textit{Kn}$ = $\textit{Jn}\textit{αn},$ and, with the binomial coefficients $(\textit{n}$ $\textit{&ngr;}),$ $\textit{cn}$ = $∑∞\textit{&ngr;}=0$ $\textit{K}\textit{n}$ (n &ngr;). The Laguerre polynomials may now be expressed in the form $\textit{Ln}(\textit{α})$ $(λ\textit{x})$ = ∑ $\textit{n}\textit{&ngr;}$ = 0 $\textit{E}(\textit{n})&ngr;\textit{x}&ngr;$ where $\textit{E}(\textit{n})\textit{&ngr;}$ = $\textit{J&ngr;}$ $(\textit{&ngr;}$ + $1)(\textit{α}$ + 2) ··· $(\textit{α}$ + $\textit{n}/\textit{n}!,$ and the final series (19) emerges term by term. The speed of convergence depends how close $λλ\textit{m}$ is to 0.5.As a numerical example we present, in the figure, curves for the frequency function (10) computed for a spectral density ƒ(λ) corresponding to white noise (or bandlimited noise with a correlation function sin $π\textit{t}/π\textit{t}$ sampled at the zeros), a rectangular spectral window $\textit{w}(λ)$ of bandwidth π/10 and sample sizes $\textit{n}$ = 20, 30, 40, 50, 70, 100, 150. The computation time for each curve on the IBM 650 was approximately 9 minutes, except for $\textit{n}$ = 70, 100 and 150 where the Toeplitz approximation led to equal eigenvalues and a computation time of 1 minute each. The eigenvalues of $\textit{RW},$ where both $\textit{R}$ and $\textit{W}$ are symmetric but not their product, were computed by Jacobi's method; for orders 20, 30, 40 and 50 all eigenvalues were found in 0.8, 2.7, 6.4 and 12.5 hours respectively.The authors are indebted to Professor Ulf Grenander, formerly of Brown University, now of the University of Stockholm, for stimulating advice. ISSN 00045411 Age Range 18 to 22 years ♦ above 22 year Educational Use Research Education Level UG and PG Learning Resource Type Article Publisher Date 1960-07-01 Publisher Place New York e-ISSN 1557735X Journal Journal of the ACM (JACM) Volume Number 7 Issue Number 3 Page Count 6 Starting Page 245 Ending Page 250

#### Open content in new tab

Source: ACM Digital Library