clc;
close all;
clear all;

sigma1 = 0.5*0.5;
sigma2 = 2*sigma1;
sigma3 = 0.5*sigma1;

d1 = 0; m1 = 0;
d2 = 1.5; m2 = sqrt(2*d2*d2);
d3 = 5; m3 = sqrt(2*d3*d3);

N = 100;

c{1}.x = sqrt(sigma1)*randn(N, 1);
c{1}.y = sqrt(sigma1)*randn(N, 1);

c{2}.x = sqrt(sigma2)*randn(N, 1) + m2;
c{2}.y = sqrt(sigma2)*randn(N, 1) + m2;

c{3}.x = sqrt(sigma3)*randn(N, 1) + m3;
c{3}.y = sqrt(sigma3)*randn(N, 1) + m3;

N_CLASSES = 3;


% projection of a point on a line defined by points p0 and p1
% from http://cs.nyu.edu/~yap/classes/visual/03s/hw/h2/math.pdf
p0.x = -2; p0.y = -6;
p1.x = 10; p1.y = 6;
A = [p1.x - p0.x, p1.y - p0.y
     p0.y - p1.y, p1.x - p0.x];
 
invA = inv(A);

% projection
pr{1}.x = zeros(N, 1);
pr{1}.y = zeros(N, 1);
pr{2}.x = zeros(N, 1);
pr{2}.y = zeros(N, 1);
pr{3}.x = zeros(N, 1);
pr{3}.y = zeros(N, 1);

for k=1:N
    
    for cl = 1:N_CLASSES
        q.x = c{cl}.x(k);
        q.y = c{cl}.y(k);

        p = invA * [q.x*(p1.x - p0.x) + q.y*(p1.y - p0.y); ...
                    p0.y*(p1.x - p0.x) - p0.x*(p1.y - p0.y)];

        pr{cl}.x(k) = p(1);
        pr{cl}.y(k) = p(2);
    end
    
end


plot(c{1}.x, c{1}.y, 'b*', c{2}.x, c{2}.y, 'r*', c{3}.x, c{3}.y, 'b*');
xlabel('Feature 1');
ylabel('Feature 2');
legend('Class A', 'Class B');
hold on;
plot(pr{1}.x, pr{1}.y, 'b.', pr{2}.x, pr{2}.y, 'r.', pr{3}.x, pr{3}.y, 'b.');
plot([p0.x, p1.x], [p0.y, p1.y], 'k');
axis([-2, 10, -6, 8]);



x = -3:0.01:10;
g1 = 1 / ( sigma1*sqrt(2*pi) ) * exp(-x.*x / ( 2*sigma1*sigma1 ) );
g2 = 1 / ( sigma2*sqrt(2*pi) ) * exp(-(x - d2).*(x - d2) / ( 2*sigma2*sigma2 ) );
g3 = 1 / ( sigma3*sqrt(2*pi) ) * exp(-(x - d3).*(x - d3) / ( 2*sigma3*sigma3 ) );
figure;
plot(x, 0.5*g1, 'b', x, g2, 'r', x, 0.5*g3, 'b');
xlabel('New Feature "x"');
ylabel('probability density function');
legend('Class A', 'Class B');
grid on;



