最近在写KECA的程序，它是对KPCA的一种改进方法，利用Renyi熵来选取主元个数，并将其应用于TE过程进行故障检测，选取TE过程正常训练数据500*52，测试数据960*52，
但是编程结果一直不对：高斯核参数大一些（100-128左右），计算出来的Renyi熵第一项远远大于其他项，导致主元计算不正确，只有1，但是整体的统计图走向趋势是对的；高斯核函数小一些（3.5-4.5），计算出保留的主元数较为合适（十几到几十不等），但是统计图的趋势完全不对；若是设定大的高斯核参数，自己手动设定合适的保留的主元数（15-27左右），统计图监测效果还行。
所以怎么才能让设定合适的核参数时计算出合适的主元数？我现在的程序到底哪步出了问题？
KECA基本原理
程序的整个流程：
1、导入训练数据X（500*52），计算均值，方差来标准化
2、计算高斯核函数K（500*500），核函数自己设定为128
c=128;
for i = 1 : X_row
for j = 1 : X_row
K(i,j) = exp(-(norm(Xtrain(i,:) - Xtrain(j,:)))^2/c^2);%求核矩阵，采用径向基核函数，参数c
end
end
3、对K直接进行特征值分解（eig）,求解出对应的特征值D和特征向量V，第一个特征值496.8437远远大于第二大的特征向量0.3995。
- 做特征值分解的时候是对K做，还是对K/N做？
- 求出来的特征向量需要单位化吗？怎么单位化?
[V,D]=eig(K);%求特征值和特征向量
lambda1=diag(D); %对角矩阵 求特征值
[lambda,index]=sort(real(lambda1),'descend'); %sort 'descend'每行按从大到小排序，lamda为排序后结果，index为索引 按照降序排列
L1=diag(lambda); %从大到小特征值对角矩阵形式
VVV=V(:,index); %按照特征值降序排列的特征向量 N*N
VV=zeros(X_row,X_row);
for i=1:X_row
VV(:,i) = VVV(:,i)/norm(VVV(:,i)); %特征向量归一化
end
4、计算每一项的Renyi熵值（按照熵大小排列，取前npc个特征值特征向量，熵是平方的形式还是不平方的形式？不平方有负值）
n0=ones(1,X_row);%1*N 元素全为1
entropy=zeros(1,X_row);
for i = 1 : X_row
entropy(i) = lambda(i)*(VV(:,i)'*n0')^2;%求二次熵
%entropy(i) =sqrt(lambda(i))*VV(:,i)'*n0';%求熵
end
5、按照Renyi熵值的大小降序排列对应的特征值和特征向量
[sorted_entropy,r]=sort(entropy,'descend');%将熵值按降序排列，sorted_entropy是排列后的数组，r是序号
lamda=lambda(r);%对应于熵降序排列的特征值
L=diag(lamda);%对应于熵降序排列的特征值对角阵
P=VV(:,r);%对应于熵降序排列的特征向量
6、计算Renyi熵值从大到小累计贡献99%的对应的主元个数
npc=1;
while sum(sorted_entropy(1:npc))/sum(sorted_entropy)<0.99
npc=npc+1; %求保留的主元数
end
一般情况下npc=1可以占据99.99%甚至更高的贡献率，但是这种结果是不对的，但是不知道那里不对
7、利用训练数据计算T2和SPE统计量的控制限
训练数据主元向量T怎么计算？
还是
计算出训练数据T2和SPE统计量使用核密度估计怎么计算控制限？
公式怎么使用？x是代表一个样本的T2统计量值吗？
%% 计算T2 SPE统计量
Z_1=inv(sqrt(L))*P'*K'; %N1*N2
Z=Z_1';%N2*N1
SCORE=K*P;
Z_ECA_1 =inv(sqrt(L(1:npc,1:npc)))*P(:,1:npc)'*K'; %npc*N2
Z_ECA=Z_ECA_1'; %N2*npc
%T2统计量
for i=1:X_row
T2_1(i)=Z_ECA(i,:)*inv(L(1:npc,1:npc))*Z_ECA(i,:)'; % 1*1
end
%T22=diag(Z_ECA*inv(L(1:npc,1:npc))*Z_ECA'); %整体矩阵计算
%Q统计量
Q_1= sum(Z.^2,2)-sum(Z_ECA.^2 ,2);
8、CS统计量怎么计算（训练数据的，以及测试数据的）
CS统计量公式
%训练数据CS统计量及控制限计算计算
Phi_ECA1=sqrt(L(1:npc,1:npc))*P(:,1:npc)'; %npc*N
Phi_ECA=Phi_ECA1';
SM_ECA=0;
for i=1:X_row
SM_ECA=SM_ECA+Phi_ECA1(:,i); %npc*1
end
M_ECA=SM_ECA/X_row; %npc*1
CS1=zeros(1,X_row);
for i=1:X_row
CS1(i)=1-M_ECA'*Phi_ECA1(:,i)/(norm(M_ECA')*norm(Phi_ECA1(:,i)));
end
%利用核密度估计计算控制限
alpha = 0.99;
[f,xi] = ksdensity(CS1);
delta=xi(2)-xi(1);
sum_bb=0;
block=0;
for i=100:-1:1
sum_bb=sum_bb+delta*f(i);
if sum_bb>(1-alpha) & block==0;
bbb=i;
block=1;
end
end
CSUCL=xi(1)+delta*bbb;
%测试数据CS统计量
[K_t] = computeKM(Xtest,Xtrain,c);%求核矩阵，采用径向基核函数，参数c N2*N1
Zt_1=inv(sqrt(L))*P'*K_t'; %N1*N2
Zt=Zt_1';%N2*N1
Zt_ECA1 =inv(sqrt(L(1:npc,1:npc)))*P(:,1:npc)'*K_t'; %npc*N2
Zt_ECA=Zt_ECA1'; %N2*npc
CS2=zeros(1,Xt_row);
for i=1:Xt_row
CS2(i)=1-M_ECA'*Zt_ECA1(:,i)/(norm(M_ECA')*norm(Zt_ECA1(:,i)));
end
