APPLICATIONS OF SAS IML: SPH 5421 notes.021
EIGENVECTORS AND EIGENVALUES
AND
SIMULATING MULTIVARIATE NORMAL DISTRIBUTIONS
========================================================================
EIGENVECTORS AND EIGENVALUES
Let A be a matrix.
Vector x is said to be an eigenvector for matrix A if there exists
a constant a and a nonzero vector x = (x1 x2)' such that
Ax = a*x.
The value of a is defined to be the eigenvalue associated with the
eigenvector x.
Example:
| 5 2 |
Find an eigenvector for the matrix A = | | .
| 2 8 |
Solution: Suppose x = (x1 x2)', and Ax = a*x for some constant a.
This implies that:
5 * x1 + 2 * x2 = a * x1, and
2 * x1 + 8 * x2 = a * x2.
This looks like 2 equations in 3 unknowns: x1, x2, and a. In
general you cannot expect to find solutions for such a system.
There is another way to view this. Start with the equation
Ax = a*x.
This can be written as
Ax - a*x = 0, or
(A - a*I)x = 0.
This is again a system of two equations with 3 unknowns. If you
knew what a was, it would be two equations in two unknowns, x1 and x2.
If the matrix A - a*I were nonsingular, the only solution to a
system like this is x = (0 0)'. This is not an eigenvector, however,
because by definition eigenvectors have to be nonzero. The only way to
have a nonzero solution is if A - a*I is a singular matrix.
This means that the determinant of A - a*I must be 0.
| 5 - a 2 |
Note that A - a*I = | |,
| 2 8 - a |
and that the determinant is : (5 - a) * (8 - a) - 2 * 2 .
Hence we need to find the solutions a to the equation
(5 - a) * (8 - a) - 4 = 0.
This reduces to
a^2 - 13*a + 36.
The roots of this quadratic equation are:
a1 = 9 and a2 = 4.
Thus 9 and 4 are the eigenvalues of the matrix A.
To find the corresponding eigenvectors, you need to solve, for
example,
A * x = 9 * x, where x = (x1, x2)' .
This reduces to two equations,
5 * x1 + 2 * x2 = 9 * x1 and
2 * x1 + 8 * x2 = 9 * x2.
There is not a unique solution to this system of equations. Any
solution will satisfy the requirement that
x2 = 2 * x1.
We can therefore specify x1 and find the corresponding x2. For
example, let x1 = 1; then x2 = 2 and a corresponding eigenvector is
x = (1 2)'. However, any other vector of the form (u 2*u)' is
also one of the eigenvectors of a. You can check that the other
eigenvectors are of the form (2*v -v)'.
PROC IML has functions which find eigenvalues and eigenvectors.
The functions are eigval() and eigvec(). The arguments for these
functions must be square numeric matrices. Below is a little SAS program
which finds the eigenvalues and eigenvectors for the matrix A above.
The printout follows the program.
========================================================================
footnote "~john-c/5421/eigen.sas &sysdate &systime" ;
options linesize = 80 ;
proc iml ;
a = {5 2, 2 8} ;
eigvals = eigval(a) ;
eigvecs = eigvec(a) ;
print a eigvals eigvecs ;
finish ;
========================================================================
The SAS System 1
17:39 Tuesday, March 6, 2001
A EIGVALS EIGVECS
5 2 9 0.4472136 0.8944272
2 8 4 0.8944272 -0.447214
~john-c/5421/eigen.sas 06MAR01 17:39
========================================================================
Note that SAS returns eigenvectors which are of unit length; that
is, x1^2 + x2^2 = 1. Thus in the example above,
.447^2 + .894^2 = .2 + .8 = 1.
Such vectors are called 'normalized'.
It is a surprising fact that the eigenvectors for a given symmetric
matrix are orthogonal (or perpendicular) to each other. See Problem 21a.
Two vectors (x1 x2)` and (y1 y2)` are orthogonal provided their dot
product is 0. You may recall from undergrad calculus/linear algebra
that the dot product of such vectors is a scalar, and is defined to be
(x1 x2) * (y1 y2)` = x1 * y1 + x2 * y2.
Again, you can easily check that this dot product is 0 for the two
eigenvectors in the example above.
Some matrices do not have any [real] eigenvectors. For example, matrices
of the form
| cos(x) -sin(x)|
| |
| sin(x) cos(x)|
for any value x which is not an integer multiple of pi, do not have
[real] eigenvectors.
Why not?
That is, what is the underlying geometric reason?
Note that the above matrix is not symmetric, and therefore cannot
be the covariance matrix for two random variables.
All nonsingular 3 x 3 matrices have at least one real eigenvector.
See if you can find (1) an algebraic reason for this, and (2) a
geometric reason.
All matrices do have complex eigenvalues and eigenvectors. That
however is off-topic for what we want to do here related to statistics.
========================================================================
Why are eigenvectors and eigenvalues useful in statistics?
If variables X1, X2, X3, ..., Xp have a multivariate normal
distribution with in general nonzero covariance between Xi and Xj,
it is possible to find Y1, Y2, Y3, ..., Yp such that each Yi is
a linear combination of the X's, and the Y's are independent. That is,
the X's and the Y's are related by a linear transformation, but the
Y's are easier to deal with than the original X's because they are
independent. Finding the linear transformation depends on first
finding the eigenvectors and eigenvalues of the covariance matrix
of the X's. See the next section for more details on this.
========================================================================
SIMULATING MULTIVARIATE NORMAL DISTRIBUTIONS
Suppose (X, Y)' is a bivariate normal random variable with expectation (EX,
EY)' and covariance matrix V. How can you simulate observations which have this
distribution?
The question is easy to answer if V is a diagonal matrix. In this case, X and
Y are independent, and the two diagonal entries of V are var(X) and var(Y). Let Sx
= sqrt(var(X)) and Sy = sqrt(var(Y)). Then you can simulate (X, Y) with the
desired distribution as follows:
X = EX + Sx * rannor(seed) ;
Y = EY + Sy * rannor(seed) ;
The question is somewhat harder if V is not diagonal. It is harder still if
you want to simulate observations with a given covariance matrix in dimensions
higher than 2.
A simple case is the following. Suppose you want a covariance matrix
of the form
| 20 14 14 |
| |
[*****] S = | 14 20 14 | .
| |
| 14 14 20 |
for a multivariate normal distribution.
Thus you want three variables, X1, X2, and X3 such that
var(X1) = var(X2) = var(X3) = 20,
and
covar(X1, X2) = covar(X1, X3) = covar(X2, X3) = 14.
Let Z, W, U, and V be chosen so that (1) all are independent of
each other, and (2) var(Z) = 14, var(W) = var(U) = var(V) = 6.
Define X1 = Z + W,
X2 = Z + U, and
X3 = Z + V.
Note that
covar(X1, X2) = covar(Z + W, Z + U)
= covar(Z, Z) + covar(Z, U) + covar(Z, W) + covar(W, U).
= var(Z) + 0 + 0 + 0 = 14,
because the last three covariance terms here are 0 [independence].
Also,
var(X1) = var(Z) + var(W) = 14 + 6 = 20,
and similarly for X2 and X3.
What this means is that if you want to generate X1, X2, and X3 with
so that their covariance matrix is S as above in equation [*****], then
you need only generate Z ~ N(0, 14), W ~ N(0, 6), U ~ N(0, 6) and V ~ N(0, 6),
and add them as shown above to get X1, X2, and X3.
A general approach to this problem is the following. Given a covariance
matrix A, find a transformation P and a diagonal matrix D with only nonnegative entries, such
that
A = P * D * P`.
Then note that D is the covariance matrix of random variables (U, V)` which are
independent and normal, and var(U) = D11 (the first diagonal entry in D, and
similarly var(V) = D22. Perform the transformation
(X, Y)` = P(U, V)`.
The covariance matrix of (X, Y)` is
cov(X, Y)` = cov(P * (U, V)`)
= P * cov(U, V)` * P`
= P * D * P` = A.
This essentially solves the problem. To simulate (X, Y)`, you first simulate
(U, V)` (independent normal variates, as described above), then apply the
transformation P.
But how do you find P and D ?
PROC IML has a function called 'eigen' which does exactly what is needed for
this. Specifically, the command
call eigen(E, P, A) ;
works as follows. The input matrix A is an n x n symmetric matrix. The output
vector E is a column vector, E = (e1, e2, e3, ..., en)` of eigenvalues of A, and the
matrix P is an orthogonal n x n matrix whose columns are the eigenvectors of A.
What this means is that, if Pi is the i-th column of P, then
A * Pi = Pi * ei, which implies
A * P = P * D,
where D is the diagonal matrix whose diagonal entries are the eigenvalues ei.
'Orthogonality' of P means that P * P' = P' * P = I, the identity matrix.
This implies that P` is the inverse of P. Thus
A = A * I = A * P * P` = P * D * P`,
which is what we want.
Note that not all symmetric matrices A are covariance matrices. They must also
have the property of being positive definite. What this means is that for ANY nonzero
n x 1 column vector X,
X` * A * X > 0.
In general it is hard to tell just by looking at a symmetric matrix that it
is positive definite. But as it happens, if a symmetric matrix A is positive
definite, all of its eigenvalues are positive. Since IML gives you a way of finding
all the eigenvalues, you can always test whether a given symmetric matrix is
positive definite and therefore whether it can be a covariance matrix.
A positive semi-definite n x n matrix A is such that
X` * A * X >= 0,
for any n x 1 vector X.
Fact: if an n x n matrix which has all positive real eigenvalues is
positive definite.
Fact: All covariance matrices are positive semi-definite.
The appended program illustrates the use of IML to simulate a bivariate normal
distribution which has a pre-specified covariance matrix:
=================================================================================
/* Program to simulate data from a multivariate normal distribution */
/* with a prespecified covariance matrix. */
options linesize = 80 ;
footnote "program: /home/walleye/john-c/5421/multi.sas &sysdate &systime";
data imlgen ;
title1 'Use of PROC IML to simulate bivariate normal data with';
title2 'a given covariance matrix.' ;
proc iml ;
A = {3 2, 2 3} ; /* A covariance matrix. The object here is to */
/* simulate (x1, x2)` such that their covariance */
/* matrix is A. */
/* */
/* The IML function 'eigen' finds the eigenvalues and eigenvectors */
/* of the matrix A. The eigenvalues are stored in the vector */
/* DIAG. The eigenvectors are the columns of the matrix P. */
/* */
/* Note that P is orthogonal, and A = P * DIAG * P'. */
/* */
call eigen(DIAG, P, A); /* Note DIAG here is actually a 2 x 1 column matrix. */
print A P DIAG ;
s1 = sqrt(DIAG[1]) ; /* DIAG[1] and DIAG[2] are the eigenvalues */
s2 = sqrt(DIAG[2]) ;
seed = 20000222 ;
e1 = {1, 0} ; e2 = {0, 1} ;
V = j(100, 2, 0) ;
do i = 1 to 100 ;
y1 = s1 * rannor(seed) ; /* Simulate y1, y2 independent and */
y2 = s2 * rannor(seed) ; /* y1 ~ N(0, s1), ys ~ N(0, s2). */
/* Note that the y's have covariance */
/* matrix DIAG. */
V[i, 1] = y1 ; /* Create an n x 2 matrix of the y's */
V[i, 2] = y2 ;
end ;
W = V * P` ; /* Transform the y's ... the result will */
/* have covariance P * DIAG * P` */
/* Note W is also n x 2. */
varnames = 'y1':'y2' ;
create yobs from V [colname = varnames] ; /* Put y's on an output */
append from V ; /* dataset, yobs. */
varnames = 'x1':'x2' ; /* Put x's on an output */
create trans from W [colname = varnames] ; /* dataset, trans. */
append from W ;
quit ;
run ;
proc corr cov data = yobs ; /* Print correlation of the y's */
var y1 y2 ;
title1 'Correlation and covariance of the Ys' ;
title2 '' ;
proc corr cov data = trans ; /* Print correlation of the x's */
var x1 x2 ;
title1 'Correlation and covariance of the Xs' ;
=================================================================================
Use of PROC IML to simulate bivariate normal data with 1
a given covariance matrix.
19:32 Thursday, February 24, 2000
A P DIAG
3 2 0.7071068 0.7071068 5
2 3 0.7071068 -0.707107 1
program: /home/walleye/john-c/5421/multi.sas 24FEB00 19:32
Correlation and covariance of the Ys 2
' 19:32 Thursday, February 24, 2000
Correlation Analysis
2 'VAR' Variables: Y1 Y2
Covariance Matrix DF = 99
Y1 Y2
Y1 5.109249639 0.023975642
Y2 0.023975642 1.010058145
Simple Statistics
Variable N Mean Std Dev Sum Minimum Maximum
Y1 100 -0.14649 2.26036 -14.64895 -5.70825 4.58132
Y2 100 -0.13865 1.00502 -13.86500 -3.29248 1.85551
Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0 / N = 100
Y1 Y2
Y1 1.00000 0.01055
0.0 0.9170
Y2 0.01055 1.00000
0.9170 0.0
program: /home/walleye/john-c/5421/multi.sas 24FEB00 19:32
Correlation and covariance of the Xs 3
19:32 Thursday, February 24, 2000
Correlation Analysis
2 'VAR' Variables: X1 X2
Covariance Matrix DF = 99
X1 X2
X1 3.083629533 2.049595747
X2 2.049595747 3.035678250
Simple Statistics
Variable N Mean Std Dev Sum Minimum Maximum
X1 100 -0.20162 1.75603 -20.16241 -3.98214 3.62798
X2 100 -0.00554 1.74232 -0.55433 -4.54450 3.72547
Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0 / N = 100
X1 X2
X1 1.00000 0.66990
0.0 0.0001
X2 0.66990 1.00000
0.0001 0.0
program: /home/walleye/john-c/5421/multi.sas 24FEB00 19:32
=================================================================================
The program simulates 100 observations. What you expect to get here are the
following:
1. The variables Y1 and Y2 should be approximately independent (covariance
close to 0) and normal with mean (0, 0)` and variances approximately
equal to the eigenvalues of A.
2. The variables X1 and X2 should be approximately bivariate normal
with covariance matrix approximately equal to A.
The printed output indicates approximate agreement with both of these
statements. You don't get exact agreement because this is a finite sample. Note the
use of the 'cov' option in PROC CORR.
Note that in this case the matrix P represents a transformation
consisting of a reflection through the line x = y, followed by a
45-degree rotation in the clockwise direction.
---------------------------------------------------------------------------------
PROBLEM 20
Simulate 256 observations from a multivariate normal distribution which has the
following covariance matrix:
| 2 1 0 |
| |
A = | 1 2 1 | .
| |
| 0 1 2 |
Verify that the multivariate distribution approximately agrees with the
expected distribution by the use of PROC CORR.
----------------------------------------------------------------------------------
PROBLEM 21
Write a program which simulates a bivariate distribution for (X, Y) such that
both of the following conditions are true:
(1) X and Y both have Bernoulli distributions Ber(.5)
(2) the correlation of X and Y is 0.3.
Check that your program works by using PROC CORR.
(Note - this is not a problem where you need to use IML)
----------------------------------------------------------------------------------
PROBLEM 21a.
Let A be a symmetric matrix, and suppose u and v are two
eigenvectors with eigenvalues a and b respectively. Prove that
if a is not equal to b, the u and v must be orthogonal.
----------------------------------------------------------------------------------
/home/walleye/john-c/5421/notes.021 Last update: November 9, 2011.