%% Matlab script for explaining the covariance decomposition
%% Written by Philippe Lucidarme
%% http://www.lucidarme.me
close all;
clear all;
clc;

%% Parameters tuing

% Number of sample 
N=100;

% Standard deviation along X and Y
StdX=3;
StdY=1;

% Rotation of sample (correlation)
alpha=pi/3;


%% Create and display random sample
X=randn(2,N);
X=[StdX,0;0,StdY]*X;
X=[cos(alpha) , -sin(alpha) ; sin(alpha) , cos(alpha) ]*X;

plot (X(1,:),X(2,:),'.k');
grid on;
axis square equal;

%% Compute covariance 
Sigma=cov(X');


%% Perform singular value decomposition
[U,S,D]=svd(Sigma);


%% Display results (an ellipse with 3x standart deviation)
%% The ellipse must contain 99% of the samples
a=[0:0.1:2*pi];
Xu=[cos(a);sin(a)];
Xu= (U*sqrt(S)) * Xu;
patch (3*Xu(1,:),3*Xu(2,:),'g','FaceAlpha',0.4);

% Display Stadard deviations:
disp ('Standard deviations :');
sqrt(S)