clear all;
close all;

%兩類數據的概率
p1 = 0.6;
p2 = 0.4;
d = 3;%樣本點的維數；
%各樣本點；
load 'D:\fisher\x11.txt';% w1中各樣本點的X坐標
load 'D:\fisher\x12.txt';% w1中各樣本點的Y坐標
load 'D:\fisher\x13.txt';% w1中各樣本點的Z坐標
x11 = reshape(x11,1,[]);  % 把X坐標化為1行N1列的矩陣
x12 = reshape(x12,1,[]);
x13 = reshape(x13,1,[]);
N1 = 36; %1類樣本點有N1個；
X1 = zeros(N1,d);
for i=1:N1
    X1(i,:) = [x11(i),x12(i),x13(i)]; %矩陣X1的每一行代表一個三維的樣本點；
end
%求樣本均值向量 mi
m1 = mean(X1)'; %對X1按列求平均，得到1類樣本均值向量（3*1的列向量）；
load 'D:\fisher\x21.txt'; % w2中各樣本點的X坐標
load 'D:\fisher\x22.txt'; % w2中各樣本點的Y坐標
load 'D:\fisher\x23.txt'; % w2中各樣本點的Z坐標

x21 = reshape(x21,1,[]);
x22 = reshape(x22,1,[]);
x23 = reshape(x23,1,[]);
N2 = 36;%2類樣本點個數；
X2 = zeros(N2,d);
for i=1:N2
    X2(i,:) = [x21(i),x22(i),x23(i)]; %矩陣X2的每一行代表一個三維的樣本點；
end
%dlmwrite('sample.txt',X1);
%dlmwrite('sample.txt',X2,'-append');
dlmwrite ('sample1.txt',X1);%把變換得到的樣本點保存在文件中
dlmwrite ('sample2.txt',X2);
m2 = mean(X2)'; %對X2按列求平均，得到2類樣本均值向量（3*1的列向量）；
figure(2);
plot3(X1(:,1),X1(:,2),X1(:,3),'+r',m1(1),m1(2),m1(3),'.b',X2(:,1),X2(:,2),X2(:,3),'*r',m2(1),m2(2),m2(3),'.g');
grid on;

%樣本類內離散度矩陣Si與總類內離散度矩陣Sw 
S1 = zeros(3); %S1為3*3的矩陣
for i=1:N1
    S1 = S1 + (X1(i,:)'-m1) * (X1(i,:)'-m1)'; %X1(i,:)為1類樣本的第I個點，3維行向量；
end
S2=zeros(3);%S2為3*3的矩陣
for i=1:N2
    S2 = S2 + (X2(i,:)'-m2) * (X2(i,:)'-m2)';
end
Sw = S1 + S2;

W = inv(Sw) * (m1-m2); %最佳W值的確定,W為3*1的列向量
display('最佳的W為：');
display(W);

M1 = W' * m1; %一維Y空間各類樣本均值
M2 = W' * m2;
w0 = -((M1+M2)/2)-(log(p1/p2))/(N1+N2-2) %P(w1)與P(w2)已知時,w0的確定

load  'D:\fisher\test.txt'; %輸入的測試數據text為矩陣，每一行代表一個三維樣本點；
%load  D:\fisher\sample.txt';
%test = sample;
Y = test*W  + w0; %Y為一判決向量
a1 = find (Y>0); %分類
a2 = find (Y<=0);

figure(1)
shad1 = shadow(W,test(a1,:)'); %1類的各點在W上的投影向量
shad2 = shadow(W,test(a2,:)'); %2類的各點在W上的投影向量
shad =  shadow(W,test');
plot3(test(a1,1),test(a1,2),test(a1,3),'*r');%畫出1類點
hold on
plot3(test(a2,1),test(a2,2),test(a2,3),'+b');%畫出2類點
hold on
plot3(shad1(1,:),shad1(2,:),shad1(3,:),'ro');%畫出1類點在W上的投影向量
hold on
plot3(shad2(1,:),shad2(2,:),shad2(3,:),'bo');%畫出2類點在W上的投影向量
hold on
plot3(shad(1,:),shad(2,:),shad(3,:),'y');%畫出向量W
legend('第一類','第二類','一類投影','二類投影','W');
grid on;
hold off;