核主成分分析法代码

clear;

clc;

load normal.mat

X = normal; % X训练数据集

[Xrow, Xcol] = size(X); % Xrow：样本个数 Xcol：样本属性个数

%% 数据预处理，进行标准化出理，处理后均值为0方差为1

Xc = mean(X); % 求原始数据的均值

Xe = std(X); % 求原始数据的方差

X0 = (X-ones(Xrow,1)*Xc) ./ (ones(Xrow,1)*Xe); % 标准阵X0,标准化为均值0，方差1;

c = 20000; %此参数可调

%% 求核矩阵

for i = 1 : Xrow

for j = 1 : Xrow

K(i,j) = exp(-(norm(X0(i,:) - X0(j,:)))^2/c);%求核矩阵，采用径向基核函数，参数c end

end

%% 中心化矩阵

n1 = ones(Xrow, Xrow);

N1 = (1/Xrow) * n1;

Kp = K - N1*K - K*N1 + N1*K*N1; % 中心化矩阵

%% 特征值分解

[V, D] = eig(Kp); % 求协方差矩阵的特征向量（V）和特征值（D）

lmda = real(diag(D)); % 将主对角线上为特征值的对角阵变换成特征值列向量

[Yt, index] = sort(lmda, 'descend'); % 特征值按降序排列，t是排列后的数组，index是序号

%% 确定主元贡献率记下累计贡献率大于85%的特征值的序号放入 mianD中

rate = Yt / sum(Yt); % 计算各特征值的贡献率

sumrate = 0; % 累计贡献率

mpIndex = []; % 记录主元所在特征值向量中的序号

for k = 1 : length(Yt) % 特征值个数

sumrate = sumrate + rate(k); % 计算累计贡献率

mpIndex(k) = index(k); % 保存主元序号

if sumrate > 0.85

break;

end

end

npc = length(mpIndex); % 主元个数

%% 计算负荷向量

for i = 1 : npc

zhuyuan_vector(i) = lmda(mpIndex(i)); % 主元向量

P(:, i) = V(:, mpIndex(i)); % 主元所对应的特征向量（负荷向量）end

zhuyuan_vector2 = diag(zhuyuan_vector); % 构建主元对角阵