clear
clc
DATA = xlsread('datag.xlsx');
rowrank=randperm(size(DATA,1));
DATA_R = DATA(rowrank,:);
data_train = DATA_R(1:8000,1:2);
Y_train = DATA_R(1:8000,3);
data_test = DATA_R(8001:9223,1:2);
Y_real = DATA_R(8001:9223,3);
hideNums = 18;
maxnumber = 100000;
alpha = 0.1;
yta = 0.5;
train_number = 8000;
test_number = 1223;
C = rand(hideNums,2)*20-10;
d = zeros(hideNums,1); %% RBF宽度
% for i = 1:hideNums
% dmin = Inf;
% for j = 1:hideNums
% D = sqrt((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2);
% if(D < dmin&&i ~= j)
% dmin = D;
% end
% end
% d(i) = 5*dmin;
% end
for i = 1:hideNums
dmax = 0;
for j = 1:hideNums
D = sqrt((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2);
if(D > dmax&&i ~= j)
dmax = D;
end
end
d(i) = dmax/sqrt(2*hideNums);
end
%% 初始化隐层输出矩阵
H=zeros(train_number,hideNums); %%隐含层输出矩阵
for i = 1:train_number
for j = 1:hideNums
H(i,j)=exp(-((data_train(i,1)-C(j,1))^2+(data_train(i,2)-C(j,2))^2)/(2*d(j)^2));
end
end
A = H';
B = A*H;
G = inv(B);
W = G*A*Y_train;
%% 计算测试误差
% 隐含层输出矩阵
dC = zeros(18,2);
dd = zeros(18,1);
dW = zeros(18,1);
for k = 1:maxnumber
test = zeros(train_number,hideNums);
for i = 1:train_number
for j=1:hideNums
test(i,j)=exp( -( (data_train(i,1)-C(j,1))^2+(data_train(i,2)-C(j,2))^2 )/(2*d(j)^2) );
T(i,j) = (data_train(i,1)-C(j,1))^2+(data_train(i,2)-C(j,2))^2;
end
end
Y_pre = test*W;
en = Y_pre - Y_train;
EN = 20*en;
E = (1/length(EN)) *sum((EN).^2)
END(k) = E;
if(E < 0.00001)
break;
end
for j =1 : hideNums
dW(j) = -yta/train_number*en'*test(:,j);
SUM1 = zeros(1,2);
SUM2 = zeros(1,1);
for l = 1:train_number
SUM1 = SUM1 + en(l)*test(l,j)*(data_train(l,:)-C(j,:));
end
dC(j,:)= -yta/(train_number*d(j)^2)*SUM1;
%dC(j,:)= -yta*W(j)/d(j)^2*SUM1;
for l = 1:train_number
SUM2 = SUM2 + en(l)*test(l,j)*T(l,j);
end
dd(j)= -yta/(train_number*d(j)^3)*SUM2;
%dd(j)= -yta*W(j)/d(j)^3*SUM2;
end
W = W + alpha*dW; %更新权重矩阵
C = C + alpha*dC; %更新中心向量
d = d + alpha*dd; %更新宽度
end
for i = 1:test_number
for j=1:hideNums
H_real(i,j)=exp( -( (data_test(i,1)-C(j,1))^2+(data_test(i,2)-C(j,2))^2 )/(2*d(j)^2) );
end
end
Y_test = H_real*W;
en_test = Y_test - Y_real;
EN_test = 20*en_test;
E_test = (1/length(EN_test)) *sum((EN_test).^2)
figure(1)
plot(en,'b-');
axis([1,8000,-0.1,0.1])
figure(2)
plot(END,'r-');
figure(3)
plot(en_test,'y');
