|
楼主 |
发表于 2006-11-13 20:55
|
显示全部楼层
这是我的程序,帮忙看看吧,谢谢!
%KPCA用于过程监控
clear all;
close all;
t=0;%采样周期
Variances=0.040;%高斯核参数
k=0;%主元数
h=0;%特征空间维数
Vsum=0;
%获取实验数据
for i=1:1:100
t=t+0.01;
x(i,1)=t;
x(i,2)=t^2-3*t;
x(i,3)=-t^3+3*t^2;
end
randn('seed',0);
e=[0.1*randn(100,1),0.1*randn(100,1),0.1*randn(100,1)];
xe=x+e;
%归一化实验数据
modelXe_normalization=zscore(xe);
Vmean=mean(xe(1:100,:));
Vstd=std(xe(1:100,:));
%计算核矩阵(使用高斯核)
for m=1:1:100
for n=m:1:100
mediaVector=modelXe_normalization(m,:)-modelXe_normalization(n,:);
KernelMatrix(m,n)=exp(-norm(mediaVector)^2/(2*Variances^2));
KernelMatrix(n,m)=KernelMatrix(m,n);
end
end
%中心化核矩阵
ell=size(KernelMatrix,1);
In=ones(ell,ell)./ell;
centralKernelMatrix=KernelMatrix-KernelMatrix*In-In*KernelMatrix+In*KernelMatrix*In;
%奇异值分解核矩阵
[U,S] = svd(centralKernelMatrix);
while(h<=ell)&&(Vsum<99) %凭经验计算特征空间维数
h=h+1;
Vsum=Vsum+S(h,h)/sum(diag(S)) * 100;
end
%通过交叉检验决定主元数,确定主元模型
k=23;
V=U(:,1:k);
L=S(1:k,1:k);
inverseL=diag(1./diag(L));
sqrtL=diag(sqrt(diag(L)));
invesqrtL=diag(1./diag(sqrtL));
KernelMatrixFeature=invesqrtL*V'*centralKernelMatrix;
allKernelMatrixFeature=diag(1./(sqrt(diag(S))))*U'*centralKernelMatrix;
for i=1:1:100
time(i)=i;
modelT2(i)=KernelMatrixFeature(:,i)'*inverseL*KernelMatrixFeature(:,i)*10;
end
%模型数据曲线输出
plot(time,limitT2,'b',time,modelT2,'r');
xlabel('time(s)');ylabel('T2'); |
|