%**********Fisher數據的PCA算法************
clear;
clc;
XP=[50 33 14 02 1;64 28 56 22 3;65 28 46 15 2;67 31 56 24 3;63 28 51 15 3;46 34 14 03 1;69 31 51 23 3;62 22 45 15 2;59 32 48 18 2; 46 36 10 02 1; 61 30 46 14 2; 60 27 51 16 2;65 30 52 20 3; 56 25 39 11 2; 65 30 55 18 3; 58 27 51 19 3;68 32 59 23 3; 51 33 17 05 1; 57 28 45 13 2; 62 34 54 23 3;77 38 67 22 3; 63 33 47 16 2; 67 33 57 25 3; 76 30 66 21 3;49 25 45 17 3; 55 35 13 02 1; 67 30 52 23 3; 70 32 47 14 2;64 32 45 15 2; 61 28 40 13 2; 48 31 16 02 1; 59 30 51 18 3;55 24 38 11 2; 63 25 50 19 3; 64 32 53 23 3; 52 34 14 02 1;49 36 14 01 1; 54 30 45 15 2; 79 38 64 20 3; 44 32 13 02 1;67 33 57 21 3; 50 35 16 06 1; 58 26 40 12 2; 44 30 13 02 1;77 28 67 20 3; 63 27 49 18 3; 47 32 16 02 1; 55 26 44 12 2;50 23 33 10 2; 72 32 60 18 3; 48 30 14 03 1; 51 38 16 02 1;61 30 49 18 3; 48 34 19 02 1; 50 30 16 02 1; 50 32 12 02 1;61 26 56 14 3; 64 28 56 21 3; 43 30 11 01 1; 58 40 12 02 1;51 38 19 04 1; 67 31 44 14 2; 62 28 48 18 3; 49 30 14 02 1;51 35 14 02 1; 56 30 45 15 2; 58 27 41 10 2; 50 34 16 04 1;46 32 14 02 1; 60 29 45 15 2; 57 26 35 10 2; 57 44 15 04 1;50 36 14 02 1; 77 30 61 23 3; 63 34 56 24 3; 58 27 51 19 3;57 29 42 13 2; 72 30 58 16 3; 54 34 15 04 1; 52 41 15 01 1;71 30 59 21 3; 64 31 55 18 3; 60 30 48 18 3; 63 29 56 18 3;49 24 33 10 2; 56 27 42 13 2; 57 30 42 12 2; 55 42 14 02 1;49 31 15 02 1; 77 26 69 23 3; 60 22 50 15 3; 54 39 17 04 1;66 29 46 13 2; 52 27 39 14 2; 60 34 45 16 2; 50 34 15 02 1;44 29 14 02 1; 50 20 35 10 2; 55 24 37 10 2; 58 27 39 12 2;47 32 13 02 1; 46 31 15 02 1; 69 32 57 23 3; 62 29 43 13 2;74 28 61 19 3; 59 30 42 15 2; 51 34 15 02 1; 50 35 13 03 1;56 28 49 20 3; 60 22 40 10 2; 73 29 63 18 3; 67 25 58 18 3;49 31 15 01 1; 67 31 47 15 2; 63 23 44 13 2; 54 37 15 02 1;56 30 41 13 2; 63 25 49 15 2; 61 28 47 12 2; 64 29 43 13 2;51 25 30 11 2; 57 28 41 13 2; 65 30 58 22 3; 69 31 54 21 3;54 39 13 04 1; 51 35 14 03 1; 72 36 61 25 3; 65 32 51 20 3;61 29 47 14 2; 56 29 36 13 2; 69 31 49 15 2; 64 27 53 19 3;68 30 55 21 3; 55 25 40 13 2; 48 34 16 02 1; 48 30 14 01 1;45 23 13 03 1; 57 25 50 20 3; 57 38 17 03 1; 51 38 15 03 1;55 23 40 13 2; 66 30 44 14 2; 68 28 48 14 2; 54 34 17 02 1;51 37 15 04 1;52 35 15 02 1;58 28 51 24 3;67 30 50 17 2;63 33 60 25 3;53 37 15 02 1];
X=XP(:,1:4);

[n,m]=size(X);

%取各類的樣本
index = input('選擇第幾類數據作為正常模型（1～3）：');
XClus=[];
for k=1:n
    if XP(k,5)==index
        XClus=[XClus;X(k,:)];
    end
end
%樣本數據標準化
XClusB = XStd(XClus,'標準差標準化');
%零均值化
x0 = XClusB - ones(50,1)*mean(XClusB);
%主元分析
[pc,score,latent,tsquare] = princomp(x0);

percent = input('確定方差貢獻率限(0～1）：');
%主元個數選擇
latent_sum=0;
block=0;
for i=1:m
    latent_sum=latent_sum+latent(i);
    if latent_sum/sum(latent) >= percent & block==0;
        kp=i;
        block=1;
    end
end

%求算統計量閾值
alpha = input('確定統計閾值置信度(0～1）：');
[ts_ctr,spe_ctr] = PCAThrd(x0,latent,alpha,kp);


%待測數據的預處理（步驟同樣本數據）
x1 = (X - ones(n,1)*mean(XClus))./ (ones(n,1)*std(XClus));
x1b= x1*pc;

x1bp=[x1b XP(:,m+1)];
%第一、二主元模式分布圖
figure(1);
for k=1:n;
    if x1bp(k,m+1)==1;
        plot(x1bp(k,1),x1bp(k,2),'b*')
        title('PCA based Fisher-Set Classification')
        hold on
    end
    if  x1bp(k,m+1)==2;
        plot(x1bp(k,1),x1bp(k,2),'rp')
        hold on
    end
    if  x1bp(k,m+1)==3;
        plot(x1bp(k,1),x1bp(k,2),'k.')
        hold on
    end
end
hold off;

%T方統計量圖
figure(2);
for i=1:n;
    ts(i) = x1b(i,1:kp) * diag( 1 ./ latent(1:kp) ) * x1b(i,1:kp)';

    if x1bp(i,m+1)==1;
        stem(i,ts(i),'b*')
        hold on
    end
    if  x1bp(i,m+1)==2;
        stem(i,ts(i),'rp')
        hold on
    end
    if  x1bp(i,m+1)==3;
        stem(i,ts(i),'k.')
        hold on
    end
end
title('Tsquare Control limits');
hold on;
plot(ts_ctr(ones(1,n)),'y');
hold off;

%SPE統計量圖
figure(3);
for i=1:n;
    spe(i) = x1(i,:)*(eye(m)-pc(:,1:kp)*pc(:,1:kp)')*x1(i,:)';
    if x1bp(i,m+1)==1;
        stem(i,spe(i),'b*')

        hold on
    end
    if  x1bp(i,m+1)==2;
        stem(i,spe(i),'rp')
        hold on
    end
    if  x1bp(i,m+1)==3;
        stem(i,spe(i),'k.')
        hold on
    end
end
title('SPE Control limits');
hold on;
plot(spe_ctr(ones(1,n)),'y');
hold off;