clear
%input为待分析的原始矩阵
input=[0.0167 0.0397 0.1635 0.8238 0.2830 0.0767 0.0737 0.0052
0.0115 0.0265 0.3939 0.8565 0.2735 0.1355 0.0807 0.0087
0.0157 0.0353 0.2817 0.7746 0.4827 0.0876 0.0514 0.0077
0.0259 0.0324 0.2259 0.8872 0.3518 0.0918 0.0625 0.0048
0.0005 0.2467 0.7612 0.2568 0.0871 0.0204 0.0872 0.0751
0.0229 0.6534 0.5854 0.3270 0.0582 0.0410 0.0456 0.0256
0.0046 0.7312 0.4225 0.1201 0.0765 0.0367 0.0404 0.0787
0.0525 0.5248 0.5801 0.2126 0.0954 0.0238 0.0877 0.0589
0.0050 0.0141 0.1750 0.5829 0.4590 0.1022 0.2024 0.0841
0.2382 0.0372 0.1346 0.7423 0.3067 0.2099 0.1058 0.2867
0.0081 0.0416 0.3731 0.2726 0.4139 0.1879 0.0937 0.0814
0.0556 0.0299 0.1334 0.4826 0.5161 0.1219 0.2035 0.0660
0.1438 0.3418 0.4200 0.1942 0.0950 0.0670 0.0895 0.0626
0.0990 0.2271 0.6179 0.2019 0.0918 0.1184 0.0900 0.0105
0.1352 0.3025 0.4419 0.1826 0.1620 0.0765 0.0527 0.0093
0.2230 0.2776 0.3544 0.2091 0.1181 0.0802 0.0750 0.0058
0.0016 0.0202 0.0633 0.1612 0.4463 0.1493 0.4639 0.0317
0.0074 0.0538 0.0487 0.2053 0.2982 0.3001 0.2426 0.0108
0.0102 0.0600 0.0351 0.0754 0.3919 0.2686 0.2149 0.0333
0.0431 0.0436 0.0482 0.1335 0.4888 0.1742 0.4666 0.0249
0.0870 0.0739 0.1496 0.3291 0.2592 0.0742 0.2470 0.0726
0.1776 0.1951 0.1150 0.4191 0.1731 0.1525 0.1769 0.0758
0.1090 0.2182 0.1696 0.1539 0.2337 0.1365 0.1681 0.0621
0.1180 0.0912 0.1363 0.2600 0.1955 0.0744 0.2308 0.0657
];
[e,f]=size(input);
%求每一列的平均值,并将列平均值保存在向量x2中
x=input;
for j=1:f
x1(j,1)=input(1,j);
for i=2:e
x1(j,1)=x1(j,1)+input(i,j);
end
x2(j,1)=x1(j,1)/e;
end
%对输入矩阵进行标准化计算过程
for j=1:f
s1(j,1)=(input(1,j)-x2(j,1))*(input(1,j)-x2(j,1));
for i=2:e
s1(j,1)=s1(j,1)+(input(i,j)-x2(j,1))*(input(i,j)-x2(j,1));
end
s(j,1)=sqrt(s1(j,1)/(e-1));
end
for i=1:e
for j=1:f
standard(i,j)=(input(i,j)-x2(j,1))/s(j,1);
end
end
%将标准化的矩阵送入data变量进行核主成分分析
data=standard;
kernel='rbf';
q=4;
p1=2;
p2=1;
[m,n]=size(data);
k=zeros(m,m);
%计算核矩阵
for i=1:m
for j=1:m
k(i,j)=svkernel(kernel,data(i,:),data(j,:));
end
end
D=mean(k);
E=mean(D);
J=repmat(D,m,1);
k=k-J-J'+repmat(E,m,m); %将核矩阵中心化
[pcs,evalue]=eigs(k); %核矩阵的特征向量和特征值
%特征向量规范化
sqrtL=diag(sqrt(diag(evalue)));
invsqrtL=diag(1./diag(sqrtL));
pcs=invsqrtL*pcs';
%选择数据投影后的维数 (85%的方差贡献率)
a=sum(evalue);
z=sum(a);
cc=0;
j=0;
b=zeros(1,n);
for i=1:size(a,2)
%b(i)=evalue(i)/a; %方差贡献率
b(i)=a(i)/z;
cc=cc+b(i); % 累计方差贡献率
c(i)=cc;
if cc>0.85 &j==0
j=i;
end
end
b=b';
trainFeat=pcs(1:q,:)*k;%映射到主成分空间
trainFeat=trainFeat';
bj=b(1:q); %主成分的方差贡献率
plot(trainFeat(1:4,1),trainFeat(1:4,2),'b+')%绘制2维分布散点图
hold on
grid
plot(trainFeat(5:8,1),trainFeat(5:8,2),'r*')
plot(trainFeat(9:12,1),trainFeat(9:12,2),'gh')
plot(trainFeat(13:16,1),trainFeat(13:16,2),'kd')
plot(trainFeat(17:20,1),trainFeat(17:20,2),'bx')
plot(trainFeat(21:24,1),trainFeat(21:24,2),'ro')
hold off