%% 2. Génération d'une variable aléatoire gaussienne
close all;

res = zeros(100, 3);

for N = 1 : 10 : 10001
    k = floor((N - 1) / 100) + 1;
    tab = 3 + 2 * randn(N, 1);
    mu = mean(tab);
    res(k, 1) = N;
    res(k, 2) = mu;
    sigma = var(tab);
    res(k, 3) = sigma;
end;

figure(); plot(res(:,1), res(:,2));
figure(); plot(res(:,1), res(:,3));

%% 3. Génération de vecteurs aléatoires gaussiens
close all;

% 1)
tab = repmat([4 9], 100, 1) + randn(100, 2);

figure(); plot(tab(:,1), tab(:,2), '.');

% 2)
N = 100
tab = randn(N, 2);
m = [4 9];
s = sqrt([1 0; 0 6]);
for k = 1 : N
    tab(k,:) = tab(k,:) * s + m;
end;

mean(tab)
cov(tab)
figure(); plot(tab(:,1), tab(:,2), '.');

% 3) sur papier

% 4)
N = 100
tab = randn(N, 2);
s = [2 2; 2 5];
u = sqrtm(s);
for k = 1 : N
    tab(k,:) = tab(k,:) * u;
end;
mean(tab)
cov(tab)
alpha = atan(4 / -3) / 2
figure(); plot(tab(:,1), tab(:,2), '.');
[V, Sd] = eig(s);
Sd
V * Sd * transpose(V)

% 5)
m1 = [4 9]; m2 = [8.5 7.5]; m3 = [6 3.5];
s1 = [2 2; 2 5]; s2 = [2 -2; -2 5]; s3 = [7 -4; -4 7];

x1 = repmat(m1, 100, 1) + randn(100, 2) * sqrtm(s1);
x2 = repmat(m2, 100, 1) + randn(100, 2) * sqrtm(s2);
x3 = repmat(m3, 100, 1) + randn(100, 2) * sqrtm(s3);

figure(); plot(x1(:, 1), x1(:, 2), '.r', x2(:, 1), x2(:, 2), '.g', x3(:, 1), x3(:, 2), '.b');

%% 4. Courbes d'équidensité

% 1)
x_range = linspace(0.27, 12.5, 57); y_range = linspace(-2, 15, 57);

% 2)
dens1 = zeros(57, 57);

invs1 = inv(s1);

for i = 1 : 57
    for j = 1 : 57
        x = [x_range(i) y_range(j)] - m1;
        dens1(i, j) = exp(-x * invs1 * transpose(x) / 2);
    end;
end;

dens1 = dens1 / (2 * pi * sqrt(det(s1)));

% 3)
figure();
hold on;
contour(x_range, y_range, transpose(dens1));
% On voit que les courbes d'équidensité sont des ellipses centrées en m1

% 4)

dens2 = zeros(57, 57);

invs2 = inv(s2);

for i = 1 : 57
    for j = 1 : 57
        x = [x_range(i) y_range(j)] - m2;
        dens2(i, j) = exp(-x * invs2 * transpose(x) / 2);
    end;
end;

dens2 = dens2 / (2 * pi * sqrt(det(s2)));

contour(x_range, y_range, transpose(dens2));

dens3 = zeros(57, 57);

invs3 = inv(s3);

for i = 1 : 57
    for j = 1 : 57
        x = [x_range(i) y_range(j)] - m3;
        dens3(i, j) = exp(-x * invs3 * transpose(x) / 2);
    end;
end;

dens3 = dens3 / (2 * pi * sqrt(det(s3)));

contour(x_range, y_range, transpose(dens3));
hold off;

figure();
mesh(x_range, y_range, transpose(dens1 + dens2 + dens3));

%% 5. Visualisation des frontières

Z = ones(57, 57);

for i = 1 : 57
    for j = 1 : 57
        [a, Z(i, j)] = max([dens1(i, j) dens2(i, j) dens3(i, j)]);
    end;
end;

figure();
hold on;
contour(x_range, y_range, transpose(Z));
plot(x1(:, 1), x1(:, 2), '.r', x2(:, 1), x2(:, 2), '.g', x3(:, 1), x3(:, 2), '.b');
hold off;

%% 6. Application

load('/cal/homes/refeyton/tps/matlab/bayes/voitures.mat');

tab = [cars(:,1) cars(:, 5)];

class1 = tab(find(cars(:,8) == 1), :);
m1 = mean(class1);
s1 = cov(class1);
class2 = tab(find(cars(:,8) == 2), :);
m2 = mean(class2);
s2 = cov(class2);
class3 = tab(find(cars(:,8) == 3), :);
m3 = mean(class3);
s3 = cov(class3);

figure();
hold on;
plot(class1(:, 1), class1(:, 2), '.r');
plot(class2(:, 1), class2(:, 2), '.g');
plot(class3(:, 1), class3(:, 2), '.b');
hold off;

x_range = linspace(5, 50, 100); y_range = linspace(1500, 5500, 100);

% Calculate all 3 densities

dens1 = zeros(100, 100);

invs1 = inv(s1);

for i = 1 : 100
    for j = 1 : 100
        x = [x_range(i) y_range(j)] - m1;
        dens1(i, j) = exp(-x * invs1 * transpose(x) / 2);
    end;
end;

dens1 = dens1 / (2 * pi * sqrt(det(s1)));

figure();
hold on;
contour(x_range, y_range, transpose(dens1));

dens2 = zeros(100, 100);

invs2 = inv(s2);

for i = 1 : 100
    for j = 1 : 100
        x = [x_range(i) y_range(j)] - m2;
        dens2(i, j) = exp(-x * invs2 * transpose(x) / 2);
    end;
end;

dens2 = dens2 / (2 * pi * sqrt(det(s2)));

contour(x_range, y_range, transpose(dens2));

dens3 = zeros(100, 100);

invs3 = inv(s3);

for i = 1 : 100
    for j = 1 : 100
        x = [x_range(i) y_range(j)] - m3;
        dens3(i, j) = exp(-x * invs3 * transpose(x) / 2);
    end;
end;

dens3 = dens3 / (2 * pi * sqrt(det(s3)));

contour(x_range, y_range, transpose(dens3));
hold off;

Z = ones(100, 100);

for i = 1 : 100
    for j = 1 : 100
        [a, Z(i, j)] = max([dens1(i, j) dens2(i, j) dens3(i, j)]);
    end;
end;

figure();
hold on;
contour(x_range, y_range, transpose(Z));
plot(class1(:, 1), class1(:, 2), '.r', class2(:, 1), class2(:, 2), '.g', class3(:, 1), class3(:, 2), '.b');
hold off;