clc
clear
close all
[num,txt,raw] = xlsread('a.xlsx') ;
txt(1,:) = [];
tr = {'正常','低能放电','高能放电','中低温过热','高温过热'};
for ii = 1:length(txt)
for jj = 1:length(tr)
if strcmp(txt{ii,6},tr{jj})
num1(ii,1) = jj;
end
end
end
%% parameters
N=300; % 训练集样本数
Nts=1000; % 测试集/预测集样本数
%% data generation and display
% [x,t]=datagen(N);
[coeff,score,latent,tsquared,explained,mu] = pca(num(1:387,:));
x = score(1:300,1);
t = num1(1:300,1);
figure(1);
hold on;
plot(t,'k*');
%% RVM
K=RBFfun(x,x); %生成K矩阵：N*(N+1)
% 随机初始化 alpha（系数w的先验方差倒数） 和
% beta（样本点误差的方差倒数，p(t|w,x,1/beta）=N(t|y(x),beta)=N(t|w*K,1/beta)
m=size(K,2);
alp=rand(1,m);
beta=rand();
ee = zeros(100,1) ;
best = inf;
for ii=1:100
sig=pinv(diag(alp)+beta*(K'*K)); % 系数w的后验方差矩阵Sigma
mu=sig*K'*t*beta; % 系数w的后验均值mu/u
gamma=1-alp.*diag(sig)';
alp_old=alp;
beta_old=beta;
idx=abs(alp)<1e3; % 部分alp会趋向于无穷大，对应的mu会趋向于会向于0，对于这部分alp不再更新
alp(idx)=gamma(idx)./(mu(idx)'.^2);
beta=(N-sum(gamma))/((t-K*mu)'*(t-K*mu));
% 判断收敛则退出循环
tmp_err=max(abs(alp(idx)-alp_old(idx))./abs(alp(idx)))+abs(beta-beta_old)/abs(beta);
if tmp_err<0.1,
break;
end;
if best>tmp_err
best = tmp_err;
end
ee(ii) = best;
end;
% 计算并呈现训练集的预测结果，注意，这里在计算预测结果时仅使用了有效的mu值，即不为零的mu值
tpred=K(:,idx)*mu(idx);
for ii = 1:length(tpred)
if tpred(ii)>4.1
tpred(ii) = 5;
else
tpred(ii) = round(tpred(ii));
end
end
plot(tpred,'ro');
xts = score(301:387,1);
tts = num1(301:387,1);
figure (2);
plot(tts,'*');
hold on;
xtsK=RBFfun(xts,x);
ttspred=xtsK(:,idx)*mu(idx);
for ii = 1:length(ttspred)
if ttspred(ii)>4.1
ttspred(ii) = 5;
elseif ttspred(ii)>2.5
ttspred(ii) = 3;
elseif ttspred(ii)>1.8
ttspred(ii) = 2;
else
ttspred(ii) = 1;
end
end
plot(ttspred,'o');
hold off
title('测试图');
xlabel('x');
ylabel('y');
legend('期望数据','预测数据');
figure (3)
plot(ee,'k--')
xlabel('迭代次数')
ylabel('误差')
