January 26, 2017 21:47 · Jan 26, 2017
diff --git a/pca_animation.m b/pca_animation.m
@@ -0,0 +1,113 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Matlab code to produce PCA animations shown here:
+% http://stats.stackexchange.com/questions/2691
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Static image
+
+clear all
+rng(42)
+
+X = randn(100,2);
+X = X*chol([1 0.6; 0.6 0.6]);
+X = bsxfun(@minus, X, mean(X));
+
+[a,b] = eig(cov(X));
+
+fig = figure('Position', [100 100 1000 400]);
+set(gcf,'color','w');
+axis([-3.5 3.5 -3.5 3.5])
+axis square
+hold on
+scatter(X(:,1), X(:,2), 'b', 'filled')
+
+%% Rotating animation
+
+fig = figure('Position', [100 100 1000 400]);
+set(gcf,'color','w');
+axis([-3.5 3.5 -3.5 3.5])
+axis square
+hold on
+
+for alpha = 1:1:179
+    w = [cosd(alpha) sind(alpha)];
+    z = X*w'*w;
+
+    cla
+    for i=1:100
+        plot([X(i,1) z(i,1)], [X(i,2) z(i,2)], 'r')
+    end
+
+    plot(w(1)*3.5*[-1 1], w(2)*3.5*[-1 1], 'k')
+    plot(-w(2)*2*[-1 1], w(1)*2*[-1 1], 'Color', [.6 .6 .6])
+
+    scatter(z(:,1), z(:,2), 'r', 'filled')
+    scatter(X(:,1), X(:,2), 'b', 'filled')
+    scatter(0,0,65,'k','filled', 'LineWidth', 2, 'MarkerFaceColor', [1 1 1], 'MarkerEdgeColor', [0 0 0])
+
+    a1 = 3.5;
+    a2 = 4.5;
+    plot(a(1,2)*[-a2 -a1], a(2,2)*[-a2 -a1], 'm', 'LineWidth', 2)
+    plot(a(1,2)*[ a1  a2], a(2,2)*[ a1  a2], 'm', 'LineWidth', 2)
+    drawnow
+
+    frame = getframe(fig);
+    if alpha == 1
+        [imind,map] = rgb2ind(frame.cdata,16,'nodither');
+    else
+        imind(:,:,1,alpha) = rgb2ind(frame.cdata,map,'nodither');
+    end
+end
+imwrite(imind,map, 'animation_pca.gif', 'DelayTime', 0, 'LoopCount', inf)
+
+%% Pendulum animation
+
+fig = figure('Position', [100 100 1000 400]);
+set(gcf,'color','w');
+axis([-3.5 3.5 -3.5 3.5])
+axis square
+hold on
+
+alpha = -45;
+omega = 0;
+for t = 1:1:1000
+    w = [cosd(alpha) sind(alpha)];
+    z = X*w'*w;
+
+    M = sum(sum((z * [0 1; -1 0]) .* (X-z), 2));
+    omega = omega + M;
+    omega = omega * 0.93;
+    alpha = alpha + omega/40;
+    if(abs(omega)<1 && abs(alpha-abs(atand(a(2,2)/a(1,2)))) < 1)
+        break
+    end
+
+    cla
+    for i=1:100
+        plot([X(i,1) z(i,1)], [X(i,2) z(i,2)], 'r')
+    end
+
+    plot(w(1)*3.5*[-1 1], w(2)*3.5*[-1 1], 'k')
+    plot(-w(2)*2*[-1 1], w(1)*2*[-1 1], 'Color', [.6 .6 .6])
+
+    scatter(z(:,1), z(:,2), 'r', 'filled')
+    scatter(X(:,1), X(:,2), 'b', 'filled')
+    scatter(0,0,65,'k','filled', 'LineWidth', 2, 'MarkerFaceColor', [1 1 1], 'MarkerEdgeColor', [0 0 0])
+
+    a1 = 3.5;
+    a2 = 4.5;
+    plot(a(1,2)*[-a2 -a1], a(2,2)*[-a2 -a1], 'm', 'LineWidth', 2)
+    plot(a(1,2)*[ a1  a2], a(2,2)*[ a1  a2], 'm', 'LineWidth', 2)
+
+    pause(0.01)
+
+    drawnow
+
+    frame = getframe(fig);
+    if t == 1
+        [imind,map] = rgb2ind(frame.cdata,16,'nodither');
+    else
+        imind(:,:,1,t) = rgb2ind(frame.cdata,map,'nodither');
+    end
+end
+imwrite(imind,map, 'animation_pca_pendulum.gif', 'DelayTime', 0, 'LoopCount', inf)