%
%Harjoitus 8, tehtävä 5
%

m1 = [0; 0]; m2 = [0; 2];
C1 = [2 1; 1 2]; C2 = [2 -1; -1 2];
n = 5000;
x1 = mynorgen2(m1, C1, n);
x2 = mynorgen2(m2, C2, n);
P1 = 0.5;
P2 = 1 - P1;
axis('square')
axis([-6 6 -6 6])
plot(x1(1,:), x1(2,:), '.', x2(1,:), x2(2,:), 'x')
title('\bullet: luokka 1, x: luokka 2')

pause %paina return jatkaaksesi
x = [x1 x2];
cl = norbay(x, P1, m1, C1, m2, C2);
est = P1 * sum(cl(1:n) ~= 1)/n + P2 * sum(cl(n+1:2*n) ~= 2)/n
% visualisoidaan päätösalueet:
[xx, yy] = meshgrid(-6:0.3:6, -6:0.4:6);
xx = xx(:); yy = yy(:);
x = [xx'; yy'];
cl = norbay(x, P1, m1, C1, m2, C2);
y1 = x(:, (cl == 1)); y2 = x(:, (cl == 2));
plot(y1(1,:), y1(2,:), '.', y2(1,:), y2(2,:), 'x')
title('\bullet: päätösalue 1, x: päätösalue 2')


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

function x = mynorgen2(mu, Sigma, n)
%
%Generoidaan n kpl d-ulotteisia multinormaalisti jakautuneita
%satunnaisvektoreita
%
%mu = odotusarvo (dx1 vektori)
%Sigma =  kovarianssimatriisi (d x d) 
%
% function x = mynorgen2(mu, Sigma, n)
%

d = length(mu);

[U, lam] = eig(Sigma);

%Normalisoidaan ominaisvektorit
t = sqrt(sum(U.*U));
U = U./t(ones(d,1),:);
A = U*sqrt(lam);
x = A*randn(d,n) + mu(:, ones(n,1));


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

function class = norbay(x, P1, m1, C1, m2, C2)
% class = norbay(x, P1, m1, C1, m2, C2)
%
% Kvadraattinen luokitin kahdelle normaalisti jakautuneelle luokalle
% kun ehdolliset luokkien odotusarvot ja kovarianssimatriisit ovat
% m1, m2, C1 ja C2.  P1 on luokan 1 prioritodennäköisyys.
%
% x(:, i)	i's luokiteltava vektori 
% class(i)	vektorin x(:, i) luokka luokittimen mukaan; 1 tai 2

% tarkistetaan dimensiot:
[dim, n] = size(x);
[pm1, qm1] = size(m1);
[pm2, qm2] = size(m2);
[pc1, qc1] = size(C1);
[pc2, qc2] = size(C2);
if (dim ~= pm1 | dim ~= pm2 | qm1 ~= 1 | qm2 ~= 1 | pc1 ~= qc1 | pc2 ~= qc2 ...
	| pc1 ~= dim | pc2 ~= dim)
    error('Epäyhteensopivat dimensiot.')
end

% calculate qf1(i) = (x(:,i) - m1)' inv(C1) (x(:,i) - m1), i = 1..n
X1 = x - m1(:, ones(1, n));
A = inv(C1) * X1;
if (dim == 1), qf1 = X1 .* A;
else, qf1 = sum(X1 .* A); end

% calculate qf2(i) = (x(:,i) - m2)' inv(C2) (x(:,i) - m2), i = 1..n
X1 = x - m2(:, ones(1, n));
A = inv(C2) * X1;
if (dim == 1), qf2 = X1 .* A;
else, qf2 = sum(X1 .* A); end

thresh = 0.5 * log(det(C2) / det(C1)) + log(P1 / (1 - P1));
class = 1 + (0.5 * (qf1 - qf2) > thresh);