S_MaxMinVariance.m

clc; clear; close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input parameters

Mu=[.05 .1]';
s=[.1 .1]';
nu=40;
r=-.9;

NumSimulations=10000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate sample
C=[1 r; r 1];
Ones=ones(NumSimulations,1);
Y = Ones*Mu' + (Ones*s').*mvtrnd(C,nu,NumSimulations);
X = exp(Y);
m=mean(X)';
S=cov(X);

% evaluate standard deviation on a one-dim projection (versor)
Theta=[0 : pi/100 : 2*pi];
for n=1:length(Theta)
    th=Theta(n);
    
    e=[cos(th)  % versor
        sin(th)];
    Z=X*e;      % projection
    SDev(n)=std(Z);  % standard deviation
    
    Radius(n)=(e'*inv(S)*e)^(-1/2);   % radius of ellipsoid
end

% compute min and max standard deviation and respective versor
[dd,min_n]=min(SDev);
e_min= [cos(Theta(min_n))
        sin(Theta(min_n))];
s_min=SDev(min_n);

[dd,max_n]=max(SDev);
e_max=[cos(Theta(max_n))
        sin(Theta(max_n))];
s_max=SDev(max_n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plots

figure 
% scatter plot simulations
h=plot(X(:,1),X(:,2),'.');
hold on 

% plot ellipsoid
Scale=2;
PlotEigVectors=1;
PlotSquare=1;
TwoDimEllipsoid(m,S,Scale,PlotEigVectors,PlotSquare)
grid on

% plot special directions defined by the max-min versors
Center=mean(X)'*.7; % de-center plot of special directions for better display
v=Scale*[-1:.1:1]';
Ones=1+0*v;

v_min=Ones*Center'+v*s_min*e_min';
hold on
h=plot(v_min(:,1),v_min(:,2),'r');
v_max=Ones*Center'+v*s_max*e_max';
hold on 
h=plot(v_max(:,1),v_max(:,2),'r');
axis equal
grid on

% plot statistics versus geometry
figure
Scaled_Theta=Theta/(pi/2);
h=plot(Scaled_Theta,SDev);  % plot standard deviation as function of direction
hold on 
h=plot(Scaled_Theta,Radius,'r'); % plot radius of ellipsoid as function of direction
xlim([Scaled_Theta(1) Scaled_Theta(end)])
grid on
xlabel('theta/(pi/2)')
legend('st.dev on projection','radius of ellipsoid')