📄 rbf18.m
字号:
%为预测做的新程序，这里许多for循环用矩阵运算来代替，希望能够减少运算量,这个程序在逐步
%地被完善，争取做到最最完善
%加入积分,求正交基展开的系数函数名为cross
%这个程序的改动之处是DotNum与s不再相同
%修改了权值调整方程
%数据变化了，采用混沌序列,
%close 是用来关闭已生成图片的命令，其中close只关闭当前图片，close all关闭搜、所有图片
clc;clear;
tic
f=1;%结果图形标号
m=12;%过程神经元隐层的神经元个数
q=8;%一般神经元隐层的神经元个数
T=5000;%计算中心向量时需要的最大循环次数
SamNum=300;%从样本曲线中取的离散点的个数
DotNum=5;%用于预测下一点的样本的个数
s=6;%正交基展开的多项式个数
h=1.2;%h为陡度因子
MaxEpoch=3000;%训练权值的最大循环次数
CC=42;VV=42;UU=42;OO=42;StepL=0.01;E0=0.1;%===================================
rand('state',CC);
C=rand(s,m,T);
rand('state',VV);
V=rand(m,q); % V为第一隐层和第二隐层间的连接权值，权值初始化
rand('state',UU);
U=rand(1,q);
rand('state',OO);
O=rand(1,q);%==================================================
rand('state',42);O1=rand;
fprintf('RBF18.m\nm=%g,q=%g,MaxEpoch=%g,SamNum=%g,DotNum=%g,正交基展开次数s=%g,C%g,V%g,U%g,O%g,步长StepL=%g,网络误差E0=%g\n',m,q,MaxEpoch,SamNum,DotNum,s,CC,VV,UU,OO,StepL,E0);
z=[0.5495 0.5836 0.6283 0.6792 0.7326 0.7858 0.8369 0.8846 0.9283 0.9674...
1.0017 1.0309 1.0549 1.0670 1.0705 1.0696 1.0687 1.0712 1.0794 1.0944...
1.1159 1.1420 1.1693 1.1932 1.2091 1.2136 1.2056 1.1863 1.1580 1.1236...
1.0881 1.0545 1.0244 0.9973 0.9718 0.9458 0.9172 0.8847 0.8477 0.8075...
0.7657 0.7248 0.6867 0.6536 0.6272 0.6097 0.6028 0.6075 0.6232 0.6486...
0.6815 0.7204 0.7643 0.8121 0.8626 0.9135 0.9618 1.0046 1.0400 1.0674...
1.0876 1.1017 1.1116 1.1194 1.1274 1.1372 1.1502 1.1669 1.1871 1.2093...
1.2306 1.2465 1.2514 1.2419 1.2182 1.1837 1.1425 1.0983 1.0536 1.0099...
0.9679 0.9274 0.8878 0.8485 0.8087 0.7681 0.7269 0.6858 0.6460 0.6091...
0.5769 0.5515 0.5350 0.5297 0.5377 0.5602 0.5966 0.6446 0.7006 0.7607...
0.8208 0.8775 0.9281 0.9706 1.0043 1.0290 1.0454 1.0548 1.0589 1.0598...
1.0596 1.0608 1.0659 1.0769 1.0947 1.1191 1.1472 1.1740 1.1928 1.1988...
1.1911 1.1722 1.1459 1.1160 1.0855 1.0562 1.0293 1.0050 0.9824 0.9601...
0.9360 0.9082 0.8755 0.8381 0.7976 0.7569 0.7186 0.6854 0.6592 0.6418...
0.6342 0.6369 0.6494 0.6705 0.6986 0.7323 0.7706 0.8128 0.8584 0.9062...
0.9540 0.9986 1.0374 1.0691 1.0938 1.1124 1.1266 1.1383 1.1493 1.1614...
1.1755 1.1922 1.2113 1.2317 1.2510 1.2655 1.2704 1.2617 1.2384 1.2030...
1.1595 1.1117 1.0625 1.0136 0.9660 0.9198 0.8748 0.8307 0.7870 0.7437...
0.7009 0.6591 0.6191 0.5822 0.5497 0.5234 0.5056 0.4988 0.5055 0.5277...
0.5656 0.6171 0.6777 0.7420 0.8052 0.8631 0.9132 0.9541 0.9854 1.0074...
1.0212 1.0281 1.0296 1.0279 1.0250 1.0236 1.0264 1.0356 1.0527 1.0778...
1.1082 1.1388 1.1629 1.1755 1.1754 1.1647 1.1470 1.1259 1.1041 1.0838...
1.0659 1.0508 1.0375 1.0244 1.0091 0.9890 0.9621 0.9277 0.8877 0.8453...
0.8038 0.7662 0.7344 0.7098 0.6932 0.6846 0.6832 0.6881 0.6979 0.7116...
0.7288 0.7502 0.7768 0.8104 0.8513 0.8977 0.9456 0.9910 1.0311 1.0651...
1.0933 1.1168 1.1370 1.1551 1.1723 1.1892 1.2063 1.2238 1.2417 1.2596...
1.2764 1.2893 1.2939 1.2857 1.2630 1.2276 1.1833 1.1335 1.0810 1.0276...
0.9743 0.9219 0.8705 0.8205 0.7720 0.7251 0.6799 0.6367 0.5959 0.5584...
0.5251 0.4977 0.4778 0.4681 0.4713 0.4903 0.5265 0.5787 0.6422 0.7105...
0.7769 0.8369 0.8875 0.9279 0.9580 0.9787 0.9909 0.9962 0.9960 0.9923...
0.9872 0.9831 0.9828 0.9890 1.0039 1.0280 1.0596 1.0941 1.1248 1.1464...
1.1567 1.1570 1.1504 1.1400 1.1286 1.1184 1.1104 1.1050 1.1015 1.0985...
1.0935 1.0836 1.0657 1.0380 1.0012 0.9584 0.9136 0.8703 0.8309 0.7967...
0.7683 0.7456 0.7279 0.7142 0.7033 0.6946 0.6876 0.6828 0.6817 0.6868...
0.7013 0.7280 0.7672 0.8157 0.8682 0.9199 0.9676 1.0100 1.0469 1.0787...
1.1062 1.1298 1.1502 1.1677 1.1831 1.1968 1.2100 1.2240 1.2398 1.2573];
y=z(1:SamNum);
x=1:SamNum;
%randn('state',42)
yy=y;%yy=y+0.05*randn(1,SamNum);没有加入噪声的样本曲线
[X,minx,maxx,Y,miny,maxy]=premnmx(x,yy); %对x,y进行归一化得到X,Y
for i=1:(SamNum-DotNum)
SamIn1(:,i)=Y(i:i+DotNum-1); %173个训练样本输入
end
SamOut=Y(DotNum+1:SamNum); %173个期望样本输出
%SamIn=SamIn1;
SamIn=cross(SamIn1,s); %正交基展开
figure(f)%============================样本曲线图
hold on
grid
subplot(2,2,1)
plot(X,yy);%==========================
%利用软竞争学习算法求得中心向量及径向基函数的宽度参数
%=================================================
e1=1e-4;P=zeros(m,(SamNum-DotNum));t=1;D=zeros(1,m);sumE1History=[];
BHistory=[];c=zeros(s,m); % P=zeros(m,100);sumQ=0;可以不用定义
while t<=T %迭代次数的循环
B(t)=2.5+(T-2.5)*t/T;%B(t)=B+(B(T)-B)*t/T;
BHistory=[BHistory B(t)];
for j=1:m
Q(j,:)=exp(-B(t)*sum((SamIn-C(:,j,t)*ones(1,SamNum-DotNum)).^2));%得到所有样本的Q(j),(SamIn-C(:,j,t)*ones(1,SamNum-DotNum)).*(SamIn-C(:,j,t)*ones(1,SamNum-DotNum))=(SamIn-C(:,j,t)*ones(1,SamNum-DotNum)).^2
end
P=Q./(ones(m,1)*sum(Q));
for j=1:m
R=(SamIn-repmat(C(:,j,t),1,(SamNum-DotNum)))*P(j,:)';
C(:,j,t+1)=C(:,j,t)+0.04*R/sum(P(j,:),2);
end
E1=abs(sum(C(:,:,t+1)-C(:,:,t))); % 用矩阵按向量求和的方式代替for循环，简化计算
sumE1=sum(E1);
sumE1History=[sumE1History sumE1];
if sumE1>=e1
t=t+1;
else
for j=1:m
%D(j)=sum(P(j,:).*sum((SamIn-C(:,j,t)*ones(1,SamNum-DotNum)).*(SamIn-C(:,j,t)*ones(1,SamNum-DotNum))),2)/sum(P(j,:),2);
D(j)=sum((SamIn-C(:,j,t)*ones(1,SamNum-DotNum)).*(SamIn-C(:,j,t)*ones(1,SamNum-DotNum)))*P(j,:)'/sum(P(j,:),2);
end
break;
end
end
%计算中心向量和宽度系数结束
%==========================================================================
c(:,:)=C(:,:,t);IrU=StepL;IrV=StepL;IrO=StepL; % U为输出层的权值向量，StepL，个学习步长
VGHistory=[];OGHistory=[];MSEHistory=[];MinSSEHistory=[];MinMSEHistory=[];
UHistory=[];VHistory=[];OHistory=[];O1History=[];ErrHistory=[];%定义个参数历史数据矩阵
for epoch=1:MaxEpoch
%y1是第一隐层输出
y1=radbas(dist(SamIn',c)./repmat(sqrt(2*D),SamNum-DotNum,1));%dist(SamIn',c)得到46*15矩阵，即每一列是一个样本的第一隐层输出,y1(46*15).这一步骤正确.这里必须用点除，才能是矩阵对应的每个元素相除
x1=y1*V+repmat(O,SamNum-DotNum,1);%O是第二隐层神经元的阈值，总共有q个 x1=y1*V-repmat(O,SamNum-DotNum,1);
%y2是第二隐层输出
y2=logsig(x1);%得到的y2为SamNum-DotNum*q y2=1./(1+exp(-x1/h))
[U,O1]=solvelin(y2',SamOut);
%YY是网络输出
YY=U*y2'+O1*ones(1,SamNum-DotNum);%得到的YY为1*SamNum-DotNum
Error=YY-SamOut; % 这里y是网络的输出矩阵，1*100，是每个输入样本对应的网络输出
SSE=sumsqr(Error)/2; % SSE及网络误差，是平方和的一半
ErrHistory=[ErrHistory SSE];
UHistory=[UHistory U];
VHistory=[VHistory V];
OHistory=[OHistory O];%O没有进行调整
O1History=[O1History O1];
if SSE<E0
break % 误差达到要求，跳出for循环
end
%============================================
%对权值进行调整
for k=1:q %先对V进行调整，因为U和O的值都是和V有关的
%VGrad=0;
for j=1:m
VGrad=0;
for p=1:(SamNum-DotNum)
VGrad=VGrad+(YY(p)-SamOut(p))*U(k)*y2(p,k)*(1-y2(p,k))*y1(p,j);
end
V(j,k)=V(j,k)-IrV*VGrad;%V(j,k)=V(j,k)-2*IrV*VGrad;
end
OGrad=0;
for p=1:P
OGrad=OGrad-(YY(p)-SamOut(p))*U(k)*y(p,k)*(1-y2(p,k));
%OGrad=OGrad-(YY(p)+SamOut(p))/(1+exp(y1(p,:)*V(:,k)-O(q)))*(1-1/(1+exp(y1(p,:)*V(:,k)-O(q))));
%OGrad=OGrad-(YY(p)-SamOut(p))/(1+exp(y1(p,:)*V(:,k)-O(q)))*(1-1/(1+exp(y1(p,:)*V(:,k)-O(q)))); %原有的
end
O(k)=O(k)-IrO*OGrad;%O(k)=O(k)+2*IrO*OGrad;
%UGrad=Error*y2(:,k);
%U(k)=U(k)-IrU*UGrad;
end
%权值调整结束
%=================================================
MSE=mean(abs((YY-SamOut)./SamOut));
MSEHistory=[MSEHistory MSE];
end
[Ypost]=postmnmx(YY,miny,maxy);%对网络输出的训练数据进行反归一化
%figure(a)
subplot(2,2,2);
plot(X(1:(SamNum-DotNum)),SamOut,'b',X(1:(SamNum-DotNum)),YY,'r-.')
xlabel('样本输入值');ylabel('归一化后的输出');
subplot(2,2,3);
plot(x(1:(SamNum-DotNum)),yy(DotNum+1:SamNum)','b',x(1:(SamNum-DotNum)),Ypost,'r')%这里将网络输出进行反归一化
xlabel('样本输入值');ylabel('反归一化后的输出'); %与最初的y值进行比较
subplot(2,2,4);
plot(x(1:(SamNum-DotNum)),Error)
xlabel('样本输入值');ylabel('网络输出与期望输出的误差');
figure
subplot(1,2,1)
plot(1:MaxEpoch,ErrHistory,1:MaxEpoch,MSEHistory,'r')
%==================================================================
%按照误差平方和最小，找到使网络误差最小的一组权值
[xx,Num] = size(ErrHistory);
MinSSE=min(ErrHistory);
subscriptSSE=findSmallColumn(ErrHistory);
MinMSE=min(MSEHistory);
subscriptMSE=findSmallColumn(MSEHistory);
fprintf('训练结果如下：\nthe min SSE=%5.4f,and the subscript is %g.\n',MinSSE,subscriptSSE);
fprintf('the min MSE=%5.4f,and the subscript is %g.\n',MinMSE,subscriptMSE);
MinSSEHistory=[MinSSEHistory MinSSE];
MinMSEHistory=[MinMSEHistory MinMSE];
U=UHistory(:,(subscriptSSE-1)*q+1:subscriptSSE*q);
V=VHistory(:,(subscriptSSE-1)*q+1:subscriptSSE*q);
O=OHistory(:,(subscriptSSE-1)*q+1:subscriptSSE*q);
O1=O1History(subscriptSSE);
%按照平均相对误差找到的一组权值，两者可进行比较看哪一组预测效果更好，大多数是前者更好
%U=UHistory(:,(subscriptMSE-1)*q+1:subscriptMSE*q);
%V=VHistory(:,(subscriptMSE-1)*q+1:subscriptMSE*q);
%O=OHistory(:,(subscriptMSE-1)*q+1:subscriptMSE*q);
%O1=O1History(subscriptMSE);
%==========================================================================
%网络训练结束
%==========================================================================
%一下步骤为做预测
%获取测试样本
TestNum=60;
a=1:TestNum;
b=z(360-TestNum+1:360);
%b=1-exp(-a).*sin(2*pi*a);
%对a,b进行归一化得到A,B [A,minx,maxx,B,minb,maxb]=premnmx(a,b);
%A=tramnmx(a,minx,maxx);
%B=tramnmx(b,miny,miny);
randn('state',42)
bb=b;%bb=b+0.05*randn(1,TestNum);同样没有加入噪声
[A,mina,maxa,B,minb,maxb]=premnmx(a,bb);
for i=1:(TestNum-DotNum)
TestIn1(:,i)=B(i:i+DotNum-1); %143个训练样本输入
end
TestOut=B(DotNum+1:TestNum);
TestIn=cross(TestIn1,s);
%TestIn=TestIn1;
%计算网络输出
y1=radbas(dist(TestIn',c)./repmat(sqrt(2*D),TestNum-DotNum,1));%dist(TestIn',c)得到46*15矩阵，即每一列是一个样本的第一隐层输出,y1(46*15).这一步骤正确.这里必须用点除，才能是矩阵对应的每个元素相除
x1=y1*V+repmat(O,TestNum-DotNum,1);%x1=y1*V-repmat(O,TestNum-DotNum,1);
y2=logsig(x1);%得到的y2为46*14 y2=1./(1+exp(-x1/h))
YY=U*y2'+O1*ones(1,TestNum-DotNum);%得到的YY为1*46,原来的YY加上了许多噪声数据
TError=YY-TestOut; % 这里y是网络的输出矩阵，1*100，是每个输入样本对应的网络输出
TSSE=1/2*sumsqr(TError); % TSSE及网络误差，是平方和的一半
TMSE=mean(abs((YY-TestOut)./TestOut));%平均相对误差
fprintf('测试结果：TSSE=%5.4f,TMSE=%5.4f\n\n',TSSE,TMSE)
[Bpost]=postmnmx(YY,minb,maxb);
%绘制误差图
subplot(1,2,2)
plot(a(1:(TestNum-DotNum)),bb(DotNum+1:TestNum),'b',a(1:(TestNum-DotNum)),Bpost,'r')
f=f+2;
%end
%end
toc
%当循环次数过多时，会导致O值不合适，出现除数为0的情况
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -