Output from imlproj.sas
0 Graphs
|
|
Output from imlproj.sas |
Source
0 Graphs |
|---|
NOTE: Capture of log output started.
NOTE: %INCLUDE (level 1) file n:\psy6140\examples\iml\imlproj.sas is file
n:\psy6140\examples\iml\imlproj.sas.
803 +title 'IMLPROJ: Projections'; 804 +proc iml;
IML Ready
805 + *-- Define a function for a projection of y on X; 806 +start proj(y,X) global(P); 807 + reset noprint; 808 + if nrow(X)=1 then X = t(X); 809 + XPX = t(X) * X; 810 + P = X * ginv(XPX) * t(X); 811 + reset print; 812 + return( P * y ); 813 + finish;
NOTE: Module PROJ defined.
814 + 815 + *-- Define a length function (LEN); 816 +start len(X); 817 + return (sqrt( X[##,] )); 818 + finish;
NOTE: Module LEN defined.
819 +
820 + reset print log fuzz fw=6;
821 +*PROJECTIONS (BOCK: EX.2.6-1);
822 +X= {1 1, 1 1, 1 -1, 1 -1};
X 4 rows 2 cols (numeric)
1 1
1 1
1 -1
1 -1
823 +Y= t(1:4);
Y 4 rows 1 col (numeric)
1
2
3
4
824 +* 1. Project y on X[,1] and X[,2] separately; 825 +r = proj(y,X[,1]);
WARNING: Temporary variable passed as result to module PROJ, argument X.
R 4 rows 1 col (numeric)
2.5
2.5
2.5
2.5
826 +r = proj(y,X[,2]);
WARNING: Temporary variable passed as result to module PROJ, argument X.
R 4 rows 1 col (numeric)
-1
-1
1
1
827 + 828 +* 2. Resolve y into 2 orthogonal components, u and v such that 829 + y = u + v and u * v = 0; 830 +u = proj(y,X);
U 4 rows 1 col (numeric)
1.5
1.5
3.5
3.5
831 +* 3. the projection matrix is P: ; 832 +print P;
P
0.5 0.5 0 0
0.5 0.5 0 0
0 0 0.5 0.5
0 0 0.5 0.5
833 +V=(I(4) - P) * Y;
V 4 rows 1 col (numeric)
-0.5
0.5
-0.5
0.5
834 +* 4. u and v are orthogonal, and sum to Y: ; 835 +r = t(u) * v;
R 1 row 1 col (numeric)
0
836 +y = u + v;
Y 4 rows 1 col (numeric)
1
2
3
4
837 +* 5. the sq-lengths sum to y's; 838 +r = len(Y)##2;
R 1 row 1 col (numeric)
30
839 +r = len(U)##2;
R 1 row 1 col (numeric)
29
840 +r = len(V)##2;
R 1 row 1 col (numeric)
1
841 +; 842 +* 6. P and (I-P) are idempotent matrices; 843 +* --> P = P * P; 844 +* --> det(P)= det(I-P)=0 (not of full rank); 845 +* --> nrow(P) = rank(P) + rank(I-P); 846 +r = P * P;
R 4 rows 4 cols (numeric)
0.5 0.5 0 0
0.5 0.5 0 0
0 0 0.5 0.5
0 0 0.5 0.5
847 +r = DET(P);
R 1 row 1 col (numeric)
0
848 +r = echelon(P);
R 4 rows 4 cols (numeric)
1 1 0 0
0 0 1 1
0 0 0 0
0 0 0 0
849 +r = echelon(I(4) - P);
R 4 rows 4 cols (numeric)
1 -1 0 0
0 0 1 -1
0 0 0 0
0 0 0 0
851 +/*------------------------------------*
852 + | The Gram-Schmidt process |
853 + *------------------------------------*/
854 +*-- MAKES SUCCESSIVE COLUMNS OF A MATRIX ORTHOGONAL.;
855 +X= {1 2 1 6, -1 0 1 2, 1 4 3 3, 0 1 1 0, 1 2 1 0};
X 5 rows 4 cols (numeric)
1 2 1 6
-1 0 1 2
1 4 3 3
0 1 1 0
1 2 1 0
856 +*-- get lengths of each column; 857 +r = len(X);
R 1 row 4 cols (numeric)
2 5 3.6056 7
858 +XPX=t(X) * X;
XPX 4 rows 4 cols (numeric)
4 8 4 7
8 25 17 24
4 17 13 17
7 24 17 49
859 +*-- The diagonal elements of X'X are the squared column lengths; 860 +r = vecdiag(XPX);
R 4 rows 1 col (numeric)
4
25
13
49
861 +*-- but X is singular - det(XPX) = 0; 862 +r = det(XPX);
R 1 row 1 col (numeric)
0
863 +; 864 +* the Gram-Schmidt process finds new matrix, Z, which spans same space.; 865 +Z =J(nrow(X), ncol(X), 0);
Z 5 rows 4 cols (numeric)
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
866 +Z[,1] = X[,1];
Z 5 rows 4 cols (numeric)
1 0 0 0
-1 0 0 0
1 0 0 0
0 0 0 0
1 0 0 0
867 +Z[,2] = X[,2] - proj(X[,2], Z[,1]);
WARNING: Temporary variable passed as result to module PROJ, argument X.
Z 5 rows 4 cols (numeric)
1 0 0 0
-1 2 0 0
1 2 0 0
0 1 0 0
1 0 0 0
868 +* the next col is linearly dependent since orthog. proj = 0; 869 +Z[,3] = X[,3] - proj(X[,3], Z[,1]) - proj(X[,3], Z[,2]);
WARNING: Temporary variable passed as result to module PROJ, argument X. WARNING: Temporary variable passed as result to module PROJ, argument X.
Z 5 rows 4 cols (numeric)
1 0 0 0
-1 2 0 0
1 2 0 0
0 1 0 0
1 0 0 0
870 +Z[,4] = X[,4] - proj(X[,4], Z[,1]) - proj(X[,4], Z[,2]);
WARNING: Temporary variable passed as result to module PROJ, argument X. WARNING: Temporary variable passed as result to module PROJ, argument X.
Z 5 rows 4 cols (numeric)
1 0 0 4.25
-1 2 0 1.5278
1 2 0 -0.972
0 1 0 -1.111
1 0 0 -1.75
871 + 872 + *-- now, all cols of Z are orthogonal; 873 +ZPZ = t(Z) * Z;
ZPZ 4 rows 4 cols (numeric)
4 0 0 0
0 9 0 0
0 0 0 0
0 0 0 25.639
874 + *-- but they differ in length; 875 +r = len(Z);
R 1 row 4 cols (numeric)
2 3 0 5.0635
876 + *-- delete column 3 (0 vector);
877 +Z = Z[,{1 2 4}];
Z 5 rows 3 cols (numeric)
1 0 4.25
-1 2 1.5278
1 2 -0.972
0 1 -1.111
1 0 -1.75
878 + *-- and normalize the columns to make lengths = 1; 879 +Z = Z * diag( 1 / len(Z) ) ;
Z 5 rows 3 cols (numeric)
0.5 0 0.8393
-0.5 0.6667 0.3017
0.5 0.6667 -0.192
0 0.3333 -0.219
0.5 0 -0.346
880 + *--the new columns of y span the same space, 881 + but are orthogonal and = length; 882 +r = t(Z) * Z;
R 3 rows 3 cols (numeric)
1 0 0
0 1 0
0 0 1
883 + 884 +*-- The steps above are carried out by the GSORTH routine; 885 +* Z is orthonormal matrix, same size as X; 886 +* T is an upper-triangular matrix, such that X = Z*T; 887 +* flag is 0 if columns of X are linearly independent, else 1; 888 +* rank deficiency -> 0 columns of Z and 0 rows of T; 889 + 890 + call gsorth(Z,T,flag,X);
Z 5 rows 4 cols (numeric)
0.5 0 0 0.8393
-0.5 0.6667 0 0.3017
0.5 0.6667 0 -0.192
0 0.3333 0 -0.219
0.5 0 0 -0.346
T 4 rows 4 cols (numeric)
2 4 2 3.5
0 3 3 3.3333
0 0 0 0
0 0 0 5.0635
FLAG 1 row 1 col (numeric)
1
891 +
892 +*-- Where orthogonal polynomials come from;
893 + C = {1 1 1, 1 2 4, 1 3 9, 1 4 16};
C 4 rows 3 cols (numeric)
1 1 1
1 2 4
1 3 9
1 4 16
894 + call gsorth(Z,T,flag,C);
Z 4 rows 3 cols (numeric)
0.5 -0.671 0.5
0.5 -0.224 -0.5
0.5 0.2236 -0.5
0.5 0.6708 0.5
T 3 rows 3 cols (numeric)
2 5 15
0 2.2361 11.18
0 0 2
FLAG 1 row 1 col (numeric)
0
895 +quit;
Exiting IML.
NOTE: The PROCEDURE IML used 0.39 seconds.
896 +
NOTE: %INCLUDE (level 1) ending.