% total_cov.m
% 1D and 2D examples illustrating the relation ST = SW + SB
% Dr. David Pan, UAH

clear all;
% 1D case (covariance become variance)
N = 100000;
X = randn(1, N);

% Split into two groups randomly
N1 = floor(N*rand);
N2 = N - N1;

X1 = X(1:N1);
X2 = X(N1+1: N);

m = mean(X);
m1 = mean(X1);
m2 = mean(X2);

ST = sum((X - m).^2);

SW = sum((X1 - m1).^2) + sum((X2 - m2).^2);
SB = N1*(m1-m)^2 + N2*(m2-m)^2;

% ST = SW + SB ?
abs(ST - (SW+SB)) 

% 2D case
m_model = [4, 0];
C_model = [9,4; 4,9];
%rng default  % For reproducibility
Y = mvnrnd(m_model,C_model,N);

% Split Y into three groups randomly
N1 = floor(N*rand);
Diff = N - N1;
N2 = floor(Diff*rand);
N3 = N - (N1 + N2);

Y1 = Y(1:N1,:);
Y2 = Y(N1+1: N1+N2, :);
Y3 = Y(N1+N2+1: N, :);

my = mean(Y);
my1 = mean(Y1);
my2 = mean(Y2);
my3 = mean(Y3);

% Note the definition of cov() in Matlab, need to multiply by (N-1)
STy = (N-1)*cov(Y);
SWy = (N1-1)*cov(Y1) + (N2-1)*cov(Y2) + (N3-1)*cov(Y3);
SBy = N1*(my1-my)'*(my1-my) + N2*(my2-my)'*(my2-my) + N3*(my3-my)'*(my3-my);

abs(STy - (SWy + SBy)) 

