% gdama13_01
% mixing of endmembers, visualized as vectors in 3D space
clear all;
fprintf('gdama13_01\n');
gdama13_01
% box in 3D (x,y,z) space
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;
zmin = 0;
zmax = 1;
% set up 3D plot
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
% view angle
view([-7,56]);
fprintf("Caption: Samples (black arrows) are linear combinations of end-members (red arrows)\n");
Caption: Samples (black arrows) are linear combinations of end-members (red arrows)
% gdama13_02
% factors visualized as vectors in 3D space
clear all;
fprintf('gdama13_02\n');
gdama13_02
% define 3D region of (x,y,z) space
xmin = 0;
xmax = 1;
Lx = xmax-xmin;
ymin = 0;
ymax = 1;
Ly = ymax-ymin;
zmin = 0;
zmax = 1;
Lz = zmax-zmin;
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
if( 1 )
% use SVD. The minus signs are just to get the vectors
% in quadrants that look good in the plot
S=zeros(4,3);
S(1,:)=s1';
S(2,:)=s2';
S(3,:)=s3';
S(4,:)=s4';
[U,LAMBDA,V] = svd(S);
v1 = -V(:,1);
v2 = -V(:,2);
v3 = -V(:,3);
else
% just use the mean vector, and two vectors perpendicular
% to it, one of which is in the F1-F2 plane
v1 = f1+f2;
norm = sqrt(v1'*v1);
v1=v1/norm;
v3 = cross(f1,f2);
norm = sqrt(v3'*v3);
v3=v3/norm;
v2 = cross(v1,v3);
norm = sqrt(v2'*v2);
v2=v2/norm;
end
gda_arrow3(v1,'b-',3);
gda_arrow3(v2,'b-',3);
gda_arrow3(v3,'b:',3);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
% set view
view([-7,56]);
fprintf("Caption: Samples (black arrows) are linear combinations of factors (red arrows)\n");
Caption: Samples (black arrows) are linear combinations of factors (red arrows)
fprintf('Two factors (solid blue) lie in the plane of the samples, one (blue dashed) is normal to it.\n')
Two factors (solid blue) lie in the plane of the samples, one (blue dashed) is normal to it.
% gdama13_03
% visualizatiov of two different sets of factors
% was vectors in 3D space
clear all;
fprintf('gdama13_03\n');
gdama13_03
% 3D box in (x,y,z) space
xmin = 0;
xmax = 1;
Lx = xmax-xmin;
ymin = 0;
ymax = 1;
Ly = ymax-ymin;
zmin = 0;
zmax = 1;
Lz = zmax-zmin;
% plot
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors for Case 1
% factor 1
f1 = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
% factor 2
f2 = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
% samples
s1=0.8*f1+(1-0.8)*f2;
s2=0.6*f1+(1-0.6)*f2;
s3=0.4*f1+(1-0.4)*f2;
s4=0.2*f1+(1-0.2)*f2;
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
view([-7,56]);
fprintf("Caption: Samples (black arrows) and factors (red arrows)\n");
Caption: Samples (black arrows) and factors (red arrows)
% plot
figure(2);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-1, 1, -1, 1, -1, 1]);
% factors for Case 2
f1a = [1, 0.2, 0.8]';
norm = 1.2*sqrt(f1a'*f1a);
f1a=f1a/norm;
f2a = [0.2, 1, 0.8]';
norm = 0.8*sqrt(f2a'*f2a);
f2a=f2a/norm;
f1 = f1a+0.2*f2a;
norm = 1.2*sqrt(f1'*f1);
f1=f1/norm;
f2 = f2a + 0.2*f1a;
norm = 0.8*sqrt(f2'*f2);
f2=f2/norm;
gda_arrow3(f1,'r-',3);
gda_arrow3(f2,'r-',3);
gda_arrow3(s1,'k-',2);
gda_arrow3(s2,'k-',2);
gda_arrow3(s3,'k-',2);
gda_arrow3(s4,'k-',2);
% improvise 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('E_1');
ylabel('E_2');
zlabel('E_3');
view([-7,56]);
fprintf("Caption: Samples (black arrows) and an alternative set of factors (red arrows)\n");
Caption: Samples (black arrows) and an alternative set of factors (red arrows)
% gdama13_04
% factor analysis on Atlantic Rocks dataset
% using singular value decomposition
clear all;
fprintf('gdama13_04\n');
gdama13_04
% load data
D = load('../data/rocks.txt');
sio2 = D(:,1); % SiO2
tio2 = D(:,2); % TiO2
als03 = D(:,3); % Al203
feot = D(:,4); % FeO-total
mgo = D(:,5); % MgO
cao = D(:,6); % CaO
na20 = D(:,7); % Na2O
k20 = D(:,8); % K2O
Ns = size(D);
N = Ns(1);
M = Ns(2);
% compute factors and factor loadings using singular value decompostion
S=D;
[U, LAMBDA, V] = svd(S,0);
lambda = diag(LAMBDA);
Ns = length(lambda);
F = V';
C = U*LAMBDA;
% plot singular values
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot( [1:Ns], lambda, 'k-', 'LineWidth', 2 );
plot( [1:Ns], lambda, 'ko', 'LineWidth', 2 );
% title('singular values, s(i)');
xlabel('index, i');
ylabel('lambda(i)');
fprintf('Caption: Singular-values.\n');
Caption: Singular-values.
% display first five factors
for j = (1:5)
f1=F(j,:);
fprintf('factor %d\n', j);
fprintf('SiO2 %f\n', f1(1));
fprintf('TiO2 %f\n', f1(2));
fprintf('Al203 %f\n', f1(3));
fprintf('FeO-total %f\n', f1(4));
fprintf('MgO %f\n', f1(5));
fprintf('CaO %f\n', f1(6));
fprintf('Na2O %f\n', f1(7));
fprintf(' \n');
end
factor 1
SiO2 -0.908829
TiO2 -0.024638
Al203 -0.275168
FeO-total -0.177851
MgO -0.141341
CaO -0.209989
Na2O -0.044611
factor 2
SiO2 0.007684
TiO2 -0.037474
Al203 -0.301583
FeO-total -0.018421
MgO 0.923193
CaO -0.226917
Na2O -0.058457
factor 3
SiO2 -0.161689
TiO2 -0.126343
Al203 0.567828
FeO-total -0.659205
MgO 0.255748
CaO 0.365682
Na2O -0.041738
factor 4
SiO2 0.209819
TiO2 0.151367
Al203 0.176021
FeO-total -0.427461
MgO -0.118643
CaO -0.780043
Na2O 0.302367
factor 5
SiO2 0.309495
TiO2 -0.100476
Al203 -0.670083
FeO-total -0.585155
MgO -0.195193
CaO 0.207980
Na2O -0.145318
% plot loadings of factors (f2,f3,f4) in 3D (c2,c3,c4) space
cmin=(-50);
cmax=50;
figure(2);
clf;
set(gca, 'LineWidth', 3);
set(gca, 'FontSize', 14);
hold on;
axis( [cmin, cmax, cmin, cmax, cmin, cmax] );
plot3( C(:,2), C(:,3), C(:,4), 'k.', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmin,cmin], [cmin,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmin,cmax], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmax], [cmin,cmin], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmax], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmax,cmax], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmin,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmin,cmin], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmax,cmin], [cmax,cmax], 'k-', 'LineWidth', 2 );
plot3( [cmin,cmin], [cmax,cmax], [cmax,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmax], [cmax,cmin], [cmin,cmin], 'k-', 'LineWidth', 2 );
plot3( [cmax,cmin], [cmax,cmax], [cmin,cmin], 'k-', 'LineWidth', 2 );
xlabel('C(2)');
ylabel('C(3)');
zlabel('C(4)');
view(3);
fprintf('Csption: Scatter plot of the loadings of factors 2, 3, and 4.\n');
Csption: Scatter plot of the loadings of factors 2, 3, and 4.
% gdama13_05
% factor analysis on Atlantic Rocks dataset
% using singular value decomposition
% and the varimax procedure
clear all;
fprintf('gdama13_05');
gdama13_05
% load data
D = load('../data/rocks.txt');
sio2 = D(:,1); % SiO2
tio2 = D(:,2); % TiO2
als03 = D(:,3); % Al203
feot = D(:,4); % FeO-total
mgo = D(:,5); % MgO
cao = D(:,6); % CaO
na20 = D(:,7); % Na2O
k20 = D(:,8); % K2O
Ns = size(D);
N = Ns(1);
M = Ns(2);
% compute factors and factor loadings using singular value decompostion
[U, S, V] = svd(D,0);
% keep only first five singular values
P=5;
F = V(:,1:P)';
C = U(:,1:P)*S(1:P,1:P);
% initial rotated factor matrix, FP, abd rotation matrix, MR
MR=eye(P,P);
FP=F;
% spike these factors using the varimax procedure
k = [2, 3, 4, 5]';
Nk = length(k);
% varimax is an iterative procedure, but converges very rapidly,
% so do only 3 iterations. Within each iteration, one applies
% the procedure to every pair of factors.
for iter = (1:3) % iterations
for ii = (1:Nk) % pairs of factors
for jj = (ii+1:Nk)
% spike factors i and j
i=k(ii);
j=k(jj);
% copy factors from matrix to vectors
fA = FP(i,:)';
fB = FP(j,:)';
% standard varimax procedure to determine rotation angle q
u = fA.^2 - fB.^2;
v = 2* fA .* fB;
A = 2*M*u'*v;
B = sum(u)*sum(v);
top = A - B;
CC = M*(u'*u-v'*v);
DD = (sum(u)^2)-(sum(v)^2);
bot = CC - DD;
q = 0.25 * atan2(top,bot);
cq = cos(q);
sq = sin(q);
fAp = cq*fA + sq*fB;
fBp = -sq*fA + cq*fB;
% put back into factor matrix, FP
FP(i,:) = fAp';
FP(j,:) = fBp';
% accumulate rotation
mA = MR(i,:)';
mB = MR(j,:)';
mAp = cq*mA + sq*mB;
mBp = -sq*mA + cq*mB;
MR(i,:) = mAp';
MR(j,:) = mBp';
end
end
end
% new factor loadings
CP=C*MR';
% display first five factors
for j = (1:5)
f1=FP(j,:);
fprintf('factor %d\n', j);
fprintf('SiO2 %f\n', f1(1));
fprintf('TiO2 %f\n', f1(2));
fprintf('Al203 %f\n', f1(3));
fprintf('FeO-total %f\n', f1(4));
fprintf('MgO %f\n', f1(5));
fprintf('CaO %f\n', f1(6));
fprintf('Na2O %f\n', f1(7));
fprintf('K2O %f\n', f1(8));
fprintf(' \n');
end
factor 1
SiO2 -0.908829
TiO2 -0.024638
Al203 -0.275168
FeO-total -0.177851
MgO -0.141341
CaO -0.209989
Na2O -0.044611
K2O -0.003430
factor 2
SiO2 -0.131621
TiO2 -0.070685
Al203 -0.014311
FeO-total -0.011084
MgO 0.984305
CaO -0.038976
Na2O -0.080453
K2O -0.022437
factor 3
SiO2 0.178651
TiO2 -0.077399
Al203 0.029953
FeO-total -0.979373
MgO 0.010557
CaO 0.014822
Na2O 0.016830
K2O 0.037555
factor 4
SiO2 0.182928
TiO2 0.195179
Al203 0.004503
FeO-total 0.012028
MgO 0.028210
CaO -0.913888
Na2O 0.297483
K2O 0.061593
factor 5
SiO2 0.288635
TiO2 -0.035952
Al203 -0.944592
FeO-total 0.024113
MgO 0.010209
CaO -0.002672
Na2O -0.149835
K2O 0.000397
gda_draw(' ', abs(F(2,:))', 'caption f2', abs(F(3,:))', 'caption f3', abs(F(4,:))', 'caption f4', abs(F(5,:))', 'caption f5', ...
' ', abs(FP(2,:))', 'caption f2p', abs(FP(3,:))', 'caption f3p', abs(FP(4,:))', 'caption f4p', abs(FP(5,:))', 'caption f5p');
fprintf('Caption: (Left) Factors 2-5 before varimax rotation. (Right) After rotation.\n');
Caption: (Left) Factors 2-5 before varimax rotation. (Right) After rotation.
% This is a test that the rotated matrices are correct
% (that is, that they reproduce the original sample matrix)
e = C*F-CP*FP;
E = sqrt(sum(sum(e.^2))) / sqrt(sum(sum((C*F).^2)));
fprintf('Fractional error in reconstruction: %.3f percent\n', 100*E );
Fractional error in reconstruction: 0.000 percent
% gdama13_06
% factor analysis with simple shapes
% that allows these shapes to be classified
clear all;
fprintf('gdama13_06\n');
gdama13_06
% figure 1 is plot of samples
figure(1);
clf;
% define samples. These are simple shapes
% that might represent cross-sections through
% a mountain range
N=15;
M=11;
S=zeros(N,M);
S(1,:) = [0, 1, 2, 3, 4, 5, 4, 3, 2, 1, 0];
S(2,:) = [0, 2, 3, 4, 5, 6, 5, 4, 3, 1, 0];
S(3,:) = [0, 1, 2, 4, 5, 6, 5, 4, 2, 1, 0];
S(4,:) = [0, 1, 3, 5, 5, 4.5, 4, 3.5, 3, 2.5, 0];
S(5,:) = fliplr(S(4,:));
S(6,:) = [0, 6, 6, 6, 5, 5, 5, 3.5, 2, 1, 0];
S(7,:) = fliplr(S(6,:));
S(8,:) = [0, 0.25, 0.5, 1, 3, 6, 3, 1, 0.5, 0.25, 0];
S(9,:) = [0, 1, 2.5, 3, 4, 6, 4, 3, 2.5, 1, 0];
S(10,:) = [0, 1, 2, 4, 6, 4, 2, 1, 1, 0.5, 0];
S(11,:) = fliplr(S(10,:));
S(12,:) = [0, 0, 0.5, 0, 1, 1, 2, 4, 6, 3, 0];
S(13,:) = fliplr(S(12,:));
S(14,:) = [0, 0, 0.5, 1, 3, 6, 5.5, 5.5, 4.5, 3.5, 0];
S(15,:) = fliplr(S(14,:));
% plot samples
for i=(1:N)
subplot(5,3,i);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, M-1, 0, 7] );
plot( [0:M-1], S(i,:), 'k-', 'LineWidth', 2 );
plot( [0:M-1], S(i,:), 'k.', 'LineWidth', 2 );
xlabel('i');
ylabel(sprintf('s_{%d}',i));
end
fprintf('Caption: Mountain profiles 1 through 15.\n');
Caption: Mountain profiles 1 through 15.
% factor analysis
[U, LAMBDA, V] = svd(S);
C = U*LAMBDA;
F = V';
fprintf('First 3 singular values: %f %f %f\n', LAMBDA(1,1), LAMBDA(2,2), LAMBDA(3,3) );
First 3 singular values: 38.356675 12.688680 7.490978
% figure 2 is plot of the first three factors
figure();
clf;
% plot first 3 factors. SIgn flip is just for aesthetic plot
for i=(1:3)
subplot(1,3,i);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 10, -1, 1] );
plot( [0:M-1], -F(i,:), 'k-', 'LineWidth', 2 );
plot( [0:M-1], -F(i,:), 'k.', 'LineWidth', 2 );
plot( [0:M-1], zeros(1,M), 'k:', 'LineWidth', 2 );
xlabel('i');
ylabel(sprintf('F_{%d}',i));
end
fprintf('Caption: First three EOFs (factors).\n');
Caption: First three EOFs (factors).
% figure 3 is a plot of the profiles arranged by loadings 2 and 3
% figure 2 is plot of the first three factors
figure(3);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-8, 8, -8, 8] );
plot( [-8, 8], [0, 0], 'k:', 'LineWidth', 2);
plot( [0, 0], [-8, 8], 'k:', 'LineWidth', 2);
for i=(1:N)
% loadings
c2 = C(i,2);
c3 = C(i,3);
% demean the sample and rescale
scale1 = 0.4;
x = [0:M-1];
s = S(i,:);
s = s - mean(s);
s = s*scale1;
[smax, ismax] = max(s);
xmax = x(ismax);
xbar = mean(x);
x = x-xbar;
xmax = xmax-xbar;
x = x*scale1;
xmax = xmax*scale1;
% offset by scaled version of loadings
x = x+c2;
s = s+c3;
fill( x, s, [0, 0.5, 1.0], 'FaceAlpha', 0.5 );
plot( x, s, 'k-', 'Linewidth', 2);
plot( x, s, 'k.', 'Linewidth', 2);
plot( [x(1), x(M)], [s(1), s(M)], 'k-', 'LineWidth', 2 );
text( c2+xmax, c3, sprintf('%d',i), 'HorizontalAlignment', 'center', 'FontSize', 20 );
end
xlabel('Factor loading, C_{i2}','FontName','Cambria Math');
ylabel('Factor loading, C_{i3}','FontName','Cambria Math');
fprintf('Caption: zMountsin profiles (blue) arranged according to loadings of\n');
Caption: zMountsin profiles (blue) arranged according to loadings of
fprintf('EOFs (factors) 2 and 3');
EOFs (factors) 2 and 3
% gdama13_07
% Synthetic EOF example using 2D images
% that are a mis of spatial patterns with
% time varying amplitudes
clear all;
fprintf('gdama13_07\n');
gdama13_07
% time
Nt=25;
dt=1.0;
t = dt*[0:Nt-1]';
% space Nx by Nx grid
Nx=20;
Nx2=Nx*Nx;
dx=1.0;
L = dx*(Nx-1);
x = dx*[0:Nx-1]';
% cook up data, Nt images, each Nx by Nx
% Each image contains three spatial modes
% (patterns) of variation
% define the modes
m0 = zeros(Nx,Nx);
m1 = zeros(Nx,Nx);
m2 = zeros(Nx,Nx);
mn = 0.1 * random('norm',0,1,Nx,Nx);
for p = (1:Nx)
for q = (1:Nx)
m0(p,q) = sin(pi*x(p)/L).*sin(pi*x(q)/L);
m1(p,q) = sin(pi*x(p)/L).*sin(2*pi*x(q)/L);
m2(p,q) = sin(2*pi*x(p)/L).*sin(3*pi*x(q)/L);
end
end
% build data with modes with these amplitudes at different times
% These loadings were hand-tuned to look aesthetic
c0 = [1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1 ]';
c1 = [1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5 ]';
c2 = [0, 0, 0, 1, 2, 3, 2, 1, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 0, 1, 2, 3, 2, 1 ]';
% data is sum of modes with prescribed coefficients, plus random noise
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
colormap('jet');
hold on;
D = zeros( Nx, Nx, Nt );
for p = (1:Nt)
D(:,:,p) = c0(p)*m0 + c1(p)*m1 + c2(p)*m2 + 0.7*random('norm',0,1,Nx,Nx);
subplot(5,5,p);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0 20 0 20] );
axis off;
hold on;
imagesc(D(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
fprintf('Caption: Set of 25 sequential images (time shown below eacch image).')
Caption: Set of 25 sequential images (time shown below eacch image).
% rearrange Nx by Nx image Nx*Nx as row vector
S=zeros(Nt,Nx2);
for r = (1:Nt)
for p = (1:Nx)
for q = (1:Nx)
s = Nx*(p-1)+q;
S(r,s) = D(p,q,r);
end
end
end
% test of code to rebuild Nx by Nx image froom Nx*Nx row vector
if( 0 )
D2 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D2(p,q,r) = S(r,s);
end
end
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0 20 0 20] );
axis off;
hold on;
imagesc(D2(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
end % end if test
% basic SVD analysis
[U,LAMBDA,V] = svd(S);
C = U*LAMBDA;
S2 = C*(V');
evalue = spdiags(LAMBDA,0);
% plot singular-values
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot([1:Nt], evalue, 'k-', 'LineWidth',3);
plot([1:Nt], evalue, 'ko', 'LineWidth',3);
fprintf('Caption: Singular-values.\n');
Caption: Singular-values.
% test of code to rebuild D from S2=U*LAMBDA*V'=C*V'
if( 0 )
D3 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D3(p,q,r) = S2(r,s);
end
end
figure();
clf;
set(gcf,'pos',[10 10 700 700]);
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(D3(:,:,p));
text(0,-1.5,sprintf('%d',p));
hold on;
end
end % end if test
% select largest eigenvalues, reconstruct data from them
Nf=3;
Ua=U(:,1:Nf);
LAMBDAa=LAMBDA(1:Nf,1:Nf);
Ca = Ua*LAMBDAa;
Va = V(:,1:Nf);
Sa = Ca * (Va');
% retrieve images of modes from V matrix
M = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
M(p,q,r) = V(s,r);
end
end
% plot reconstructed data
D4 = zeros( Nx, Nx, Nt );
for r = (1:Nt)
for s = (1:Nx2)
p = floor( (s-1)/Nx+0.01 ) + 1;
q = s - Nx*(p-1);
D4(p,q,r) = Sa(r,s);
end
end
figure(5);
set(gcf,'pos',[10 10 700 700]);
clf;
for p = (1:Nt)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(D4(:,:,p));
text(0,-1.5,sprintf('%d',p));
end
fprintf('Caption: Reconstruction of set of 25 sequential images using three EOFs\n');
Caption: Reconstruction of set of 25 sequential images using three EOFs
% plot first Nf modes
figure();
set(gcf,'pos',[10 10 700 700]);
clf;
for p = (1:Nf)
subplot(5,5,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
colormap('jet');
axis( [0 20 0 20] );
axis off;
imagesc(M(:,:,p));
text(0,-1.5,sprintf('%d',p));
end
fprintf('Caption: First 3 EOFs (factors).\n');
Caption: First 3 EOFs (factors).
% plot timeseries of coefficients of modes
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
for p = (1:Nf)
subplot(Nf,1,p);
set(gca,'LineWidth',3);
hold on;
set(gca,'FontSize',14);
plot( [1:Nt], Ca(:,p), 'k-', 'LineWidth', 3 );
plot( [1:Nt], Ca(:,p), 'ko', 'LineWidth', 2 );
xlabel('time');
ylabel(sprintf('C%d',p));
end
fprintf('Caption: Loadings C of first 3 EOFs (factors) as time series.\n');
Caption: Loadings C of first 3 EOFs (factors) as time series.