Harjoitus 4, tehtävä 5
======================




%Harjoitus 4, tehtävä 5: matriisin Sigma_1 valkaisu

Sigma_1 = [1.75 -1.299;-1.299 3.25];
mu = [4 0];
X = mvnrnd(mu,Sigma_1,500);
A = valkaise(Sigma_1);
Y = A*X';
Y = Y';
plot(X(:,1),X(:,2),'o',Y(:,1),Y(:,2),'x');



-----------------------------------------



%Harjoitus 4, tehtävä 5: matriisien Sigma_1 ja Sigma_2 valkaisu

Sigma_1 = [1.75 -1.299;-1.299 3.25];
Sigma_2 = [5.625 3.375;3.375 5.625];
mu_1 = [4 0];
mu_2 = [-5 0];
X1 = mvnrnd(mu_1,Sigma_1,500);
X2 = mvnrnd(mu_2,Sigma_2,500);
X = [X1' X2'];
A = valkaise_2(Sigma_1,Sigma_2);
Y = A*X;
X = X';
Y = Y';
m = max([max(abs(X)),max(abs(Y))])+1;
figure(1)
plot(X(:,1),X(:,2),'o');
axis([-m m -m m])
title('Harjoitus 4, tehtävä 5, osa 2  alkuperäiset jakaumat')
figure(2)
plot(Y(:,1),Y(:,2),'x');
axis([-m m -m m])
title('Harjoitus 4, tehtävä 5, osa 2  muunnetut jakaumat')


--------------------------------------


function A = valkaise(C)
%
%Symmetrisen ja positiivisesti definiitin matriisin C valkaisu.  
%Jos C on satunnaisvektorin X kovarianssimatriisi,
%on satunnaisvektorin Y = A*X kovarianssimatriisi yksikkömatriisi.
%
%   function A = valkaise(C)
%

[U,lam] = eig(C);
%Varmistetaan, että U:n pystyrivit ovat normalisoitu ykkösen pituisiksi
U = normc(U);
A = lam^(-1/2)*U';


----------------------------------------


function A = valkaise_2(C1,C2)
%
%Valkaisee kaksi symmetristä matriisia samanaikaisesti. C1:n tulee olla
%positiivisesti definiitti
%
%   function A = valkaise_2(C1,C2)
%

[U,lam] = eig(C1);
%Varmistetaan, että U:n pystyrivit ovat normalisoitu ykkösen pituisiksi
U = normc(U);

%Muodostetaan C1:n potenssi -1/2
B = U*lam^(-1/2)*U';
%Muodostetaan toinen diagonalisoitava matriisi
B = B*C2*B;
[V,D] = eig(B);

%Varmistetaan, että V:n pystyrivit ovat normalisoitu ykkösen pituisiksi
V = normc(V);

A = V'*U*lam^(-1/2)*U';