马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 牛小贱 于 2015-3-6 11:10 编辑
- function [CRBF width rule e RMSE ]=GDFNN(p,t,width0,parameter)
- %this is GD-FNN training program
- %input P inout data r*q de matrix
- % T outpue data s*q matrix
- % q is the number of sample data
- %parameter(1)=kdmax (max of accommodation criterion)
- %parameter(2)=kdmin (min of accommodation criterion)
- %parameter(3)=gama (decay constant)
- %parameter(4)=emax (max of output error)
- %parameter(5)=emin (min of output error)
- %parameter(6)=beta (convergence constant)
- %parameter(7)=width0 (the width of the frist rule)
- %parameter(8)=k (overlap factor of RBF units)
- %parameter(9)=ks (width updating factor)
- %parameter(10)=kerr (signnificance of a rule)
- %output:
- %CRBF is the centers of the RBF units which is u*r matrix
- %width is the widths of RBF units which is a r by u matrix
- %rule is the number of rules for each iteration
- %e is the output error for each iteration
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化系统%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if nargin<4 ;error('not enough input aruments');end
- if size(p,2)~=size(t,2);error('the input data are not correct');end
- [r,q]=size(p);[s2,q]=size(t); %三个输入3*199 个样本,1*199个输出
- pp=p.';p1=min(pp);p2=max(pp);rang=[p1' p2'];scope=(rang(:,2)-rang(:,1)); %系统的输入最大值与最小值。
- %width0=scope*ones(1,2)/2;
- kdmax=parameter(1);kdmin=parameter(2);gama=parameter(3);emax=parameter(4);%参考动态模糊神经网络这本书
- emin=parameter(5);beta=parameter(6);k=parameter(7);km=parameter(8);ks=parameter(9);kerr=parameter(10);
- ALLIN=[];ALLOUT=[];CRBF=[];width=[];
- %when frist sample data coming
- ALLIN=p(:,1);ALLOUT=t(:,1);
- %--------------seting up the intial FNN---------------------%
- diffc=abs(rang-p(:,1)*ones(1,2));
- for j=1:r
- cvar=abs(diffc(j,:));
- [cdmin,nmin]=min(cvar);
- if cdmin<=km
- CRBF(j,1)=rang(j,nmin);
- width(j,1)=width0(j,nmin);
- else
- [cdb,nb]=max(cvar);
- CRBF(j,1)=p(j,1);
- width(j,1)=cdb/0.8;
- end
- end
- w1=CRBF.';rule(1)=1;
- %cacluating the first out error
- a0=exp(-mdist(ALLIN,CRBF,width));
- a01=[a0 p(:,1).'];
- w2=ALLOUT/a01.';
- a02=w2*a01.';
- sse(1)=sumsqr(ALLOUT-a02)/s2;rmse(1)=sqrt(sse(1));
- %%%%%%%%%%%%%%%%%%%%主循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- for i=2:q
- IN=p(:,i);OUT=t(:,i);ALLOUT=[ALLOUT OUT];ALLIN=[ALLIN IN];
- [r,N]=size(ALLIN);[r,s]=size(CRBF);
- dd=mdist(IN,CRBF,width);
- md=sqrt(dd);
- [d ind]=min(md);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %--------------计算实际输出误差----------%
- ai=exp(-dd);
- ai1=transf(ai,IN);
- ai2=w2*ai1;
- errout=t(:,i)-ai2;
- errout1=errout.*errout;
- errout2=sum(errout1)/s2;
- e(i)=sqrt(errout2);
- %-------------参数选择-----------------%
- if i<=q/3
- ke=emax;
- kd=kdmax;
- ks=0.8;
- elseif i>q/3 && i<2*q/3
- ke=max(emax*beta.^(i-q/3),emin);
- kd=max(kdmax*gama.^(i-q/3),kdmin);
- ks=0.9;
- else
- ke=emin;kd=kdmin;ks=0.95;
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%马氏距离大于给定值的处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if d>kd %假如最小马氏距离的开根号,大于预先给定的Kd ,与模糊规则的完备性有关
- if e(i)>ke
- %计算宽度和中心值
- extc=[CRBF rang]; %因半闭模糊集没有讨论在其论域的起点或终点是否存在两个隶属度函数
- extw=[width width0];
- diffc=extc-IN*ones(1,s+2); %计算输入数据域边界集的欧式距离
- for J=1:r %代表r
- cvar=abs(diffc(J,:)); %计算欧距离
- [cdmin nmin]=min(cvar);
- if cdmin<=km %预定义的参数,该常数控制相邻隶属函数的相似性
- CRBF(J,s+1)=extc(J,nmin); %因而某个输入变量的隶属函数可能会重复
- width(J,s+1)=extw(J,nmin);
- else
- CRBF(J,s+1)=IN(J);
- width(J,s+1)=k*cdmin; %给定新的宽度值
- end
- end
- %%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %-----------------------完全可以删去----------------------------------%
- %但留下来,最初始阶段的近似有较好的效果
- %比较compare the CRBF delect the repeating ones
- oldCRBF=CRBF;oldwidth=width;
- CCRBF=[];newwd=[];
- [r,s1]=size(CRBF);
- while s1>1
- redc=CRBF;
- redc(:,1)=[];
- for j2=1:s1-1
- dcc=CRBF(:,1)-redc(:,j2);
- if all(dcc==0);break;end
- end
- if any(dcc~=0)
- CCRBF=[CCRBF CRBF(:,1)];
- newwd=[newwd width(:,1)];
- end
- CRBF(:,1)=[];width(:,1)=[];
- s1=s1-1;
- end
- CCRBF=[CCRBF CRBF(:,1);];
- newwd=[newwd width(:,1)];
- CRBF=CCRBF;
- width=newwd;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%
-
-
- w1=CRBF.';[u r]=size(w1);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算所有规则的ERR%%%%%%%%%%%%%%%%%%%%%%%%%
- A=exp(-mdist(ALLIN,CRBF,width));
- A0=transf(A,ALLIN);
- if u*(r+1)<=N
- %-------------计算误差衰减率--%
- err=cperr(A0,ALLOUT);
- errT=err.';
- err1=zeros(u,s2*(r+1));
- err1(:)=errT;
- err21=err1.'; err22=sum(err21.*err21)/(s2*(r+1));
- err23=sqrt(err22); %err23 代表各条规则的最终误差的减少率
- No=find(err23<kerr);
- %%%%%%%%%%%%%%%%%%%%%%%%%规则是否需要剔除%%%%%%%%%%%%%%%%%%%%%%%%%
- if ~isempty(No)
- % delc=CRBF(:,No);
- % delw=width(:,No);
- CRBF(:,No)=[];
- w1(No,:)=[];
- width(:,No)=[];
- % rCRBF=CRBF;
- % rwidth=width;
- % err21(:,No)=[];
- [u1 r]=size(w1);
- % ------------结果参数调整与总体误差计算----------%
- A=exp(-mdist(ALLIN,CRBF,width));
- A0=transf(A,ALLIN);
- w2=ALLOUT/A0;
- A1=w2*A0;
- derr=ALLOUT-A1;
- sse=sumsqr(derr)/(i*s2);
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=u1;
- else
- %-------在 最小马氏距离大于给定参数且误差大于也大于给定误差限,后经调整后。误差率值大于给定的值
- %此时该规则不用提出。因而只需参数调整和计算误差。
- w2=ALLOUT/A0; %结果参数调整
- A1=w2*A0;
- derr=ALLOUT-A1;
- sse=sumsqr(derr)/(s2*i);
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=u;
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- else
- w2=ALLOUT/A0;
- A1=w2*A0;
- derr=ALLOUT-A1;
- sse=sumsqr(derr)/(s2*i);
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=u;
- end
- %马氏距离大于给定值且误差小于给定值,直接调整结果参数
- else
- [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
- derr=ALLOUT-A1;
- sse=sumsqr(derr)/(i*s2);
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=s;
- end
- %%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%小于最小马氏距离的处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %------------------或者小于最小马氏距离的开根号------------------------%
- else
- %%
- %高斯宽度修正
- if e(i)>ke % %在最小马氏距离小于给定值情况下,如果系统误差大于系统给定的阀值
- if s*(r+1)<=N %在规则数满足该条件下才运行本模块,N代表,此时具有的样本数,是一个不断递增的数
- aa2=exp(-mdist(ALLIN,CRBF,width));
- aa3=transf(aa2,ALLIN);
- werr=cperr(aa3,ALLOUT) %计算误差下降率
- werrT=werr.';
- werr1=zeros(s,s2*(r+1)); %s stanf for the number of rule and se stand for output vatieble
- werr1(:)=werrT;
- werr21=werr1.'; %矩阵的全新排列有1* (r+1)*u==r+1 *u
- err31=reshape(werr21(:,ind),r+1,s2);
- err32=abs(err31.');
- if s2>1
- err33=sum(err32);
- else
- err33=err32;
- end
- err3=sum(err33);
- err30=err33(1);
- rate0=err30/err3;
- err33(1)=[]; %在该条件下,虽然Xk 能被相邻的模糊规则包含,但该规则并没有重要到该RBF单元能
- err3r=err33; %很好的覆盖Xk 或者反映系统的输出在这个椭球区域有急剧变化。因此该椭球区域需要修正
- rater=err3r./(err3-err30);
- nm=find(rater<=1/r);
- width(nm,ind)=ks.*width(nm,ind); %ind 代表,对一个新的样本,找到以马氏距离最靠近该样本的那条规则
- [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
- derr=ALLOUT-A1; %A1,代表神经网络的输出
- sse=sumsqr(derr)/(i*s2); %输出项全部对应的误差的平方和,再除以总共的项
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=s;
- end
- else %在最小马氏距离小于给定值情况下,如果系统误差小于系统给定的阀值
- [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
- derr=ALLOUT-A1; %什么都不做,直接进去下一个循环。
- sse=sumsqr(derr)/(i*s2); %输出项全部对应的误差的平方和,再除以总共的项
- rmse=sqrt(sse);RMSE(i)=rmse;
- rule(i)=s;
- end
- %%
- end
- end
- end
- function y=mdist(p ,c ,b)
- %计算马氏距离,相当于第三层网络的输入。
- %p= r*q 矩阵 是采样输入矩阵
- %c=r*s 是中心矩阵
- %b =r*s 矩阵是基宽矩阵
- [r,q]=size(p); [r,s]=size(c);y=zeros(s,q);
- if r==1
- for i=1:s %s 代表规则数,r 代表输入的量
- x=c(:,i)*ones(1,q);
- d=abs(p-x);
- xx=b(:,i)*ones(1,q);
- dd=d./xx;
- y(i,:)=dd.*dd; %同一条规则中计算每个样本的马氏距离
- end
- else
- for i=1:s
- x=c(:,i)*ones(1,q); %将C的第i列复制q倍,即样本数
- d=abs(p-x); %使每个输入都有相应的规则计算距离
- xx=b(:,i)*ones(1,q);
- dd=d./xx; %同一条规则中计算每个样本的马氏距离
- %假设 dd=[ 0 0.0022 0.0044 0.0067
- % 0 0.0022 0.0044 0.0067]
- %代表有两个输入,四个样本
- y(i,:)=sum(dd.*dd);
- end
- end
- end
- function [w y a]=mweight(in,c,wb,out)
- %调整TSK FNNS 权值
- %[w2 A1 A]=mweight(p ,CRBF, width, t);
- %in 为采样数据样本,r个输入,n de sample
- %C the centered with u rule
- %out s2 as output number
- [r,n]=size(in);[r u]=size(c);
- [s2 n]=size(out);
- a=exp(-mdist(in,c,wb));%输出结果为马氏距离,后成第三层的输出=规则数*样本数
- a1=transf(a,in); %输出结果为
- w=out/a1;
- y=w*a1;
- end
- function BB=transf(A ,P)
- %计算结果参数
- %A output of RBF
- %P sample input
- [u N]=size(A);[r N]=size(P);
- for j=1:N %样本数
- for i=1:r %输入的变量
- PA((i-1)*u+1:i*u,j)=P(i,j)*A(:,j); %TSK模型,是输入变量的线性组合。各个输入变量需乘以各个规则中的响应的取值。
- %每个输入样本对应的规则的数值应该是不一样的。
- end
- end
- BB=[A;PA];
- end
- function y=cperr(a,b)
- % werr=cperr(aa3,ALLOUT);
- %b 代表系统的期望输出误差
- %a 代表TSK模型中,最后计算的动态变化部分,(包含了系统的输入量)
- %计算误差衰减率
- tT=b';PAT=a';
- [WW AW]=orthogonalize(PAT);
- WSSW=sum(WW.*WW).';
- WSStT=sum(tT.*tT).';
- y=((WW.'*tT).'.^2)./(WSStT*WSSW.'); %源于线性回归的计算。。需差线性回归的计算方法
- end
- function [w a]=orthogonalize(p)
- %该变换是每个W的基向量各自期望输出能量的贡献成为可能
- %p=W*a ,其中p and w have the same column,a is squared matrix
- [u v]=size(p);w(:,1)=p(:,1);a=eye(v);
- for k=2:v
- b=zeros(u,1);
- for i=1:k-1
- a(i,k)=w(:,i).'*p(:,k)/(w(:,i).'*w(:,i));
- b=b+a(i,k)*w(:,i);
- end
- w(:,k)=p(:,k)-b;
- end
- end
复制代码
|