RNCOR
Generates a pseudorandom orthogonal matrix or a correlation matrix.
Required Arguments
A — N by N random orthogonal matrix. (Output, if IOPT = 0; workspace if IOPT = 1; Input/Output, if IOPT = 2. If IOPT = 2, A is destroyed.)
Optional Arguments
N — The order of the matrices to be generated. (Input)
N must be at least two.
Default: N = size (A,2).
IOPT — Option indicator. (Input)
Default: IOPT = 0.
|
IOPT |
Action |
|
0 |
A random orthogonal matrix is generated in A. |
|
1 |
A random correlation matrix is generated in COR. (A is used as workspace.) |
|
2 |
A random correlation matrix is generated in COR using the orthogonal matrix input in A. |
EV — If IOPT = 1 or 2, a vector of length N containing the eigenvalues of the correlation matrix to be generated. (Input, if IOPT = 1 or 2; not used otherwise.)
The elements of EV must be positive, they must sum to N, and they cannot all be equal.
LDA — Leading dimension of A exactly as specified in the dimension statement in the calling program. (Input)
Default: LDA = size (A,1).
COR — N by N random correlation matrix. (Output, if IOPT = 1 or 2; not used otherwise.)
LDCOR — Leading dimension of COR exactly as specified in the dimension statement in the calling program. (Input)
Default: LDCOR=size(COR, 1)
FORTRAN 90 Interface
Generic: CALL RNCOR (A [, …])
Specific: The specific interface names are S_RNCOR and D_RNCOR.
FORTRAN 77 Interface
Single: CALL RNCOR (N, IOPT, EV, A, LDA, COR, LDCOR)
Double: The double precision name is DRNCOR.
Description
Routine RNCOR generates a pseudorandom orthogonal matrix A from the invariant Haar measure. For each column of A, a random vector from a uniform distribution on a hypersphere is selected and then is projected onto the orthogonal complement of the columns of A already formed. The method is described by Heiberger (1978). (See also Tanner and Thisted 1982.)
A correlation matrix is formed by applying a sequence of planar rotations to the matrix AT DA, where D = diag(EV(1), …, EV(N)), so as to yield ones along the diagonal. The planar rotations are applied in such an order that in the two by two matrix that determines the rotation, one diagonal element is less than 1.0 and one is greater than 1.0. This method is discussed by Bendel and Mickey (1978) and by Lin and Bendel (1985).
The distribution of the correlation matrices produced by this method is not known. Bendel and Mickey (1978) and Johnson and Welch (1980) discuss the distribution.
For larger matrices, rounding can become severe; and the double precision results may differ significantly from single precision results.
Comments
1. Workspace may be explicitly provided, if desired, by use of R2COR/DR2COR. The reference is:
CALL R2COR (N, IOPT, EV, A, LDA, COR, LDCOR, IWK, WK)
The additional arguments are as follows:
IWK — Work vector of length 3 * N.
WK — Work vector of length N.
2. Informational error
|
Type |
Code |
Description |
|
3 |
1 |
Considerable loss of precision occurred in the rotations used to form the correlation matrix. Some of the diagonals of COR differ from 1.0 by more than the machine epsilon. |
3. The routine RNSET can be used to initialize the seed of the random number generator. The routine RNOPT can be used to select the form of the generator.
Example
In this example, RNCOR is used to generate a 4 by 4 pseudorandom correlation matrix with eigenvalues in the ratio 1:2:3:4. (Note that the eigenvalues must sum to 4.) Routines MXTXF (IMSL MATH/LIBRARY) and EVCSF (IMSL MATH/LIBRARY) are used to check the output.
USE IMSL_LIBRARIES
IMPLICIT NONE
INTEGER I, IOPT, ISEED, J, LDA, LDCOR, N, NOUT
REAL A(4,4), COR(4,4), EV(4), EVAL(4), EVEC(4,4), FLOAT, &
SUM, XID(4,4)
INTRINSIC FLOAT
!
CALL UMACH (2, NOUT)
N = 4
LDA = 4
LDCOR = 4
EV(1) = 1.0
EV(2) = 2.0
EV(3) = 3.0
EV(4) = 4.0
! Scale the eigenvalues to sum to N.
SUM = SSUM(N,EV,1)
CALL SSCAL (N, FLOAT(N)/SUM, EV, 1)
ISEED = 123457
CALL RNSET (ISEED)
! Generate an orthogonal matrix.
CALL RNCOR (A)
WRITE (NOUT,99996) ((A(I,J),J=1,N),I=1,N)
99996 FORMAT (' A random orthogonal matrix: ', /, (5X,4F8.4))
! Check it for orthogonality.
CALL MXTXF (A, XID)
WRITE (NOUT,99997) ((XID(I,J),J=1,N),I=1,N)
99997 FORMAT (' The identity matrix?: ', /, (5X,4F8.4))
!
! Now get a correlation matrix using
! the orthogonal matrix in A, which
! will be destroyed.
IOPT = 2
CALL RNCOR (A,IOPT=IOPT, EV=EV, COR=COR)
WRITE (NOUT,99998) ((COR(I,J),J=1,N),I=1,N)
99998 FORMAT (' A random correlation matrix: ', /, (5X,4F8.4))
! Check the eigenvalues.
CALL EVCSF (COR, EVAL, EVEC)
WRITE (NOUT,99999) (EVAL(I),I=1,N)
99999 FORMAT (' The computed eigenvalues:', 4F8.4)
END
Output
A random orthogonal matrix:
-0.8804 -0.2417 0.4065 -0.0351
0.3088 -0.3002 0.5520 0.7141
-0.3500 0.5256 -0.3874 0.6717
-0.0841 -0.7584 -0.6165 0.1941
The identity matrix?:
1.0000 0.0000 0.0000 0.0000
0.0000 1.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000
0.0000 0.0000 0.0000 1.0000
A random correlation matrix:
1.0000 -0.2358 -0.3258 -0.1101
-0.2358 1.0000 0.1906 -0.0172
-0.3258 0.1906 1.0000 -0.4353
-0.1101 -0.0172 -0.4353 1.0000
The computed eigenvalues: 1.6000 1.2000 0.8000 0.4000
