Some Questions of Programming Style for the APL Statistical Library
Vector 7:1, page 90
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 90
Introduction¶
I should like to discuss some issues of programming style in the context of plans to produce an APL Statistical Library [1]. I shall consider for the sake of example a user who wishes to investigate some method of multivariate analysis and needs to produce simulated data matrices whose columns are standard Normal variates from a population with a known correlation matrix. Briefly, what he needs are independent Normal variates and a method of factorising the correlation matrix into the form \(AA^T\). For my own purposes I usually use the three functions UNIF, BOXMUL and CHOLESKY as shown below:
∇ Z←UNIF N
[1] ⍝ Result is independent U(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter.
[2] Z←(?N⍴2147483647)÷2147483647
∇
∇ Z←BOXMUL N;M
[1] ⍝ Result is independent N(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter.
[2] M←⌈0.5××/,N
[3] Z←N⍴((2,M)⍴¯2×⍟UNIF M)*0.5)×1 2∘.○○2×⌽UNIF M
∇
∇ L←CHOLESKY M;P1;P2;J;N
[1] ⍝ Result is a lower triangular matrix L such that L+.×⍉L is M, the symmetric positive-definite matrix given as the right-hand parameter.
[2] L←(⍴M)⍴0 ⋄ N←1↑⍴Z ⋄ J←1
[3] NEWCOL: P1←J+⍳N-J ⋄ P2←⍳J-1
[4] L[J;J]←(M[J;J]-+/L[J;P2]×L[J;P2])*0.5
[5] L[P1;J]←(M[P1;J]-+/L[P1;P2]×((N-J),(J-1))⍴L[J;P2])÷L[J;J]
[6] →(N≥J←J+1)/NEWCOL
∇
If CORM is the \(p \times p\) population correlation matrix, then the simulated data matrix (of shape \(n \times p\)) is given by:
… or a function DATAMAT may be written, such as:
∇ Z←CORM DATAMAT N
[1] ⍝ Result is a simulated data matrix of N(0,1) variates from a population whose correlation matrix is given by the left-hand parameter and the number of whose rows is given by the right-hand parameter.
[2] Z←⍉(CHOLESKY CORM)+.×BOXMUL(1↑⍴CORM),N
∇
The theory behind this can be obtained from Newman & Odell [2], though there is a misprint in the formula and they refer to the factorisation as Crout’s method, which is incorrect [3].
The Problem¶
A user may find this entirely satisfactory, but he has perhaps unwittingly been forced to accept the following:
(i) The Lewis, Goodman and Miller [4] random number generator
\(x_{n+1} = 16807\,x_n \pmod{2^{31} - 1}\)
when, for the sake of compatibility with other programs, he may prefer the NAG [5] generator
\(x_{n+1} = 13^{13}\,x_n \pmod{2^{59}}\)
(ii) The Box-Müller [6] method of transforming to standard Normal deviates (with a slight peculiarity of which he may be unaware, and which we shall discuss later)
(iii) Cholesky decomposition of the correlation matrix, when he may prefer, or indeed have available, the spectral decomposition, that is he may know the eigenvalues and eigenvectors of the population correlation matrix.
In other words, this programming style, where one function internally calls another, reduces the user’s choice and may not be appropriate for library functions. For let us now suppose that the user wishes to adapt the APL random number generator, either by shuffling the sequence in some way, or by using the sum, modulo unity, of two of its uniform variates, or taking some other precaution against some possible imperfections in it. He may therefore wish to use UNIF2, and instead of Box-Müller he may wish to use HASTINGS which programs a Hastings [7] approximation to the inverse of the Normal distribution function, a procedure which preserves the ranking of the original random number generator:
∇ Z←UNIF2 N
[1] ⍝ Result is independent U(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter.
[2] Z←1|(UNIF N)+(UNIF N)
∇
∇ Z←HASTINGS N;X;H
[1] ⍝ Result is independent approximate N(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter.
[2] H←0.5<X←UNIF2 N
[3] Z←(¯2×⍟|X-H)*0.5
[4] X←(2.515517+Z×0.802853+Z×0.010328)÷1+Z×1.432788+Z×0.189269+Z×0.001308
[5] Z←(1-2×H)×X-Z
∇
He is entitled to believe that, trusting his uniform variates, he should transform them as simply as possible to obtain his Normal variates. We now have a new version of DATAMAT: DATAMAT2. But this is only one of the combinations which may be chosen by the user, and each possibility results in another DATAMAT function.
∇ Z←CORM DATAMAT2 N
[1] ⍝ Result is a simulated data matrix of N(0,1) variates from a population whose correlation matrix is given by the left-hand parameter and the number of whose rows is given by the right-hand parameter.
[2] Z←⍉(CHOLESKY CORM)+.×HASTINGS(1↑⍴CORM),N
∇
The ‘peculiarity’ of BOXMUL referred to is that the uniform variates whose sine and cosine are calculated are treated in reverse order to those which are logarithmically transformed. This is an attempt to avoid the phenomenon described by Neave [8], but a full justification for it has never been carried out. Notice that the use of both sine and cosine is substantially more efficient than the use only of cosine [9], a procedure requiring the generation of twice as many uniform variates and the calculation of twice as many square roots.
Possible Solutions¶
Apart from providing all of the possible combinations, there appear to be two other solutions, neither of which is wholly satisfactory. Firstly, we can prevent functions from calling other functions and merely pass parameters. For example, HASTINGS would not act on the shape of the result required but on the rectangular variates; but this is a little less straightforward for the Box-Müller method, because if an odd number of Normal deviates is required an extra uniform variate must be used, and for acceptance-rejection methods the idea is rather absurd, although if the rejection rate is known it is possible to program such methods efficiently in APL [10].
An alternative is to pass function names. This is more convenient in a second generation APL, and specimen functions are shown below:
∇ Z←UFN ∆BOXMUL N;M
[1] ⍝ Result is independent N(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter; the U(0,1) variates required are supplied by the function whose name is given as the left-hand parameter.
[2] M←⌈0.5××/,N
[3] UFN←UFN,' ',⍕M
[4] Z←N⍴((2,M)⍴(¯2×⍟⍎UFN)*0.5)× 1 2 ∘.○○2×⌽⍎UFN
∇
∇ Z←UFN ∆HASTINGS N;X;H
[1] ⍝ Result is independent approximate N(0,1) variates in a vector or higher dimensional array whose shape is given by the right-hand parameter; the U(0,1) variates required are supplied by the function whose name is given as the left-hand parameter.
[2] H←0.5<X←⍎UFN,' ',⍕N
[3] Z←(¯2×⍟|X-H)*0.5
[4] X←(2.515517+Z×0.802853+Z×0.010328)÷1+Z×1.432788+Z×0.189269+Z×0.001308
[5] Z←(1-2×H)×X-Z
∇
∇ Z←FNS ∆DATAMAT NCORM;TXT1;TXT2
[1] ⍝ Result is a simulated data matrix of correlated N(0,1) variates; the right-hand parameter is a vector consisting of the number of rows required and the correlation matrix; the left-hand parameter is a vector of vectors containing the names respectively of the U(0,1) generator, the method of transforming to N(0,1) variates, and the method of factorising the correlation matrix.
[2] TXT1←(⊃FNS[3]),' ⊃NCORM[2]'
[3] TXT2←'(⊃FNS[1]) ',(⊃FNS[2]),' ',⍕(1↑⍴NCORM[2]),NCORM[1]
[4] Z←⍉(⍎TXT1)+.×(⍎TXT2)
∇
With such functions and a 5×5 correlation matrix, the call:
… gives 100 rows of the required data matrix, using the basic random number generator, Box-Müller, and Cholesky; whereas:
… uses the alternative random number generator and method of transformation. However the user has been given his choice of methods at some cost, both to himself in terms of what he is required to know and do, and to the programmer because the simplicity of DATAMAT as we first encountered it has disappeared entirely.
More importantly, this procedure may be highly inefficient because many successive uses of ∆DATAMAT may require the same factorisation of the same correlation matrix, or the user may wish to retain the original uniform variates for successive runs with a view to investigating some form of variance reduction technique [11].
Discussion¶
The point is that we cannot cater for every conceivable, or even likely, use of such a fundamental system as a statistical library, but if we are to cater for as many as possible, as simply as possible, some difficult decisions have to be made about the style of functions we need.
I have little doubt that in many cases we shall have functions like CHOLESKY where it is fairly clear that a single parameter is operated on and a single result given - though that is to ignore the important question of whether another result indicating that the data are appropriate (e.g. that the matrix is symmetric) and that the calculations have been satisfactorily completed (e.g. that the matrix is positive definite) should also be provided.
However, I am equally sure that to give the potential users control over their analyses we must provide many different functions which carry out essentially the same operation. Of course the decision must be based on who the potential users of the system are. Some may not wish to become involved in details and will require simple-to-use prescriptive functions, but if the APL Statistical Library is to be a unique contribution it must not simply duplicate the NAG library in an APL environment; it must cater for the statistician who wishes to construct new methodologies, and it must therefore provide him with the wherewithal to do so.
Dr. D.R. Appleton,
Division of Medical Statistics, Medical School,
University of Newcastle upon Tyne, NE3 3BU,
England.
References¶
[1] O’Hagan, A. (1990) ASL after Gregynog. Vector 6 (3) 63-65.
[2] Newman, TG & Odell, PL (1971) The generation of random variates. Griffin, London, 38-39.
[3] Ralston, A & Rabinowitz, P. (1965) A first course in numerical analysis. (2nd edition) McGraw-Hill, New York, 421-425.
[4] Lewis, PA, Goodman, AS & Miller, JM (1969) A pseudo-random number generator for the System/360. IBM Systems Journal 8 136-146.
[5] NAGFLIB:1443/0:Mk6:May77.
[6] Box, GEP & Müller, ME (1958) A note on the generation of random normal deviates. Annals of Mathematical Statistics 29 610-611.
[7] Abramowitz, M & Stegun, IA (eds) (1964) Handbook of Mathematical Functions. National Bureau of Standards, Washington, DC.
[8] Neave, HR (1973) On using the Box-Müller transformation with multiplicative congruential pseudo-random generators. Applied Statistics 22 92-97.
[9] Thomson, N. (1989) APL programs for the mathematics classroom. Springer-Verlag, New York, 106-107.
[10] Appleton, DR. (1976) Discussion on paper by Atkinson, AC & Pearce MC "The computer generation of beta, gamma and normal random variables". Journal of the Royal Statistical Society A 139 449-451.
[11] Kleijnen, JPC. (1974) Statistical techniques in simulation. Dekker, New York, 105-239.