fcollgkbefore.m

协同模糊聚类建模通过特征选择和协同模糊聚类的模糊建模方法构建T-S模型

%本程序实现了把协同模糊聚类算法和G-K算法相结合,构建T-S模型
%并用该模型对数据进行测试
%输入数据:
%ytrain:训练数据的实际输出,是一个列向量
%xtrain:训练数据矩阵,分为两组,每组的每行代表一个特征.每组特征不同
%ytest:测试数据的实际输出,是一组列向量
%xtest:测试数据矩阵,其每行代表一个特征,并且与训练数据矩阵中的第一组相同
%输出数据:
%trainRMSE:模型对训练数据的均方根误差
%testRMSE:模型对测试数据的均方根误差
%最后输出模型对测试数据的拟合图
%说明:
%由于交叉验证的原因,对不用的测试数据分组不同,所以n也有不同
clear all
fid=fopen('MPG.txt','r');
original=fscanf(fid,'%f',[8,392]);

%--------------------------------------------
%用于训练的数据
ytrain=original(1,1:2:392)';
xtrain(:,:,1)=original([2 5 7],1:2:392);
xtrain(:,:,2)=original([3 4 5],1:2:392);

%--------------------------------------------
%用于测试的数据
ytest=original(1,2:2:262)';
xtest=original([2 5 7],2:2:262);
status=fclose(fid);

%--------------------------------------------
%初始化变量
c=3;                             %规则数
r=0.01;            
n=168;
dis=1;
ntest=size(xtest,2)
datasize=size(xtrain,1);
RMSESUM=0;
coll=0.2;                        %协同系数
L=7                              %用于交叉验证方法,将测试数据分组数

%--------------------------------------------
%初始化矩阵
v=zeros(c,datasize,2);           %原形矩阵
sigma=zeros(c,datasize,2);       %协方差矩阵
d=zeros(c,n,2);                  %距离矩阵
u=rand(c,n,2);                   %模糊划分矩阵
p=zeros(1,c,2);                  %规则概率
w=ones(c,2);                     %权值
first=ones(n,1);
[center,u(:,:,1),objFcn] = fcm(xtrain(:,1:n,1)',c);        %模糊划分矩阵的初始化
[center,u(:,:,2),objFcn] = fcm(xtrain(:,1:n,2)',c);        %模糊划分矩阵的初始化

%--------------------------------------------
%交叉验证，分七组,其中六组训练，一组测试
for p=1:L
    AA=xtrain;BB=ytrain;
    pr=(p-1)*28+1;
    for pp=1:28
        BB(pr)=[];
        AA(:,pr,:)=[];
    end
    data=AA;
    y=BB;
    validationx=xtrain(:,(pr:pr+27),1);
    validationy=ytrain(pr:pr+27);
    fi(:,:,1)=data(:,:,1)';
    fi(:,:,2)=data(:,:,2)';
    
%--------------------------------------------
%计算前件参数--原形矩阵v
while(dis>r)
      oldu=u;
      for i=1:c
          for j=1:datasize
              A(i,j,1)=0;
              A(i,j,2)=0;
              B(i,1)=0;
              B(i,2)=0;
              C(i,j,1)=0;
              C(i,j,2)=0;
              D(i,j,1)=0;
              D(i,j,2)=0;
              for k=1:n
                  A(i,j,1)=A(i,j,1)+u(i,k,1)^2*data(j,k,1);
                  A(i,j,2)=A(i,j,2)+u(i,k,2)^2*data(j,k,2);
                  B(i,1)=B(i,1)+u(i,k,1)^2;
                  B(i,2)=B(i,2)+u(i,k,2)^2;
                  C(i,j,1)=C(i,j,1)+(u(i,k,1)-u(i,k,2))^2*data(j,k,1);
                  C(i,j,2)=C(i,j,2)+(u(i,k,2)-u(i,k,1))^2*data(j,k,2);
                  D(i,j,1)=D(i,j,1)+(u(i,k,1)-u(i,k,2))^2;
                  D(i,j,2)=D(i,j,2)+(u(i,k,2)-u(i,k,1))^2;
              end
              v(i,j,1)=(A(i,j,1)+coll*C(i,j,1))/(B(i,1)+coll*D(i,1));
              v(i,j,2)=(A(i,j,2)+coll*C(i,j,2))/(B(i,2)+coll*D(i,2));
          end
      end
      
%--------------------------------------------
%计算前件参数--sigma
      for i=1:c                
          for j=1:datasize
              sumdv1=0;
              sumdv2=0;
              for k=1:n
                  sumdv1=sumdv1+u(i,k,1)*(data(j,k,1)-v(i,j,1))^2;
                  sumdv2=sumdv2+u(i,k,2)*(data(j,k,2)-v(i,j,2))^2;
              end
              sigma(i,j,1)=sumdv1/sum(u(i,:,1));
              sigma(i,j,2)=sumdv2/sum(u(i,:,2));
          end
      end
      
%--------------------------------------------
%计算后件参数--th
      for i=1:c          
          bet1=diag(u(i,:,1));
          bet2=diag(u(i,:,2));
          th(:,i,1)=inv([fi(:,:,1) first]'*bet1*[fi(:,:,1) first])*[fi(:,:,1) first]'*bet1*y;
          th(:,i,2)=pinv([fi(:,:,2) first]'*bet2*[fi(:,:,2) first])*[fi(:,:,2) first]'*bet2*y;
      end
      
%--------------------------------------------
%计算规则的概率p和权值w               
      for i=1:c                    
          p(i,1)=sum(u(i,:,1))/n;
          p(i,2)=sum(u(i,:,2))/n;
          wdown1=1;
          wdown2=1;
          for j=1:datasize
              wdown1=wdown1*(sigma(i,j,1)*2*pi);
              wdown2=wdown2*(sigma(i,j,2)*2*pi);
          end
          if sqrt(wdown1)==0;
             w(i,1)=1;
          else
             w(i,1)=p(i,1)/sqrt(wdown1);
          end
          if sqrt(wdown2)==0;
             w(i,2)=1;
          else
             w(i,2)=p(i,2)/sqrt(wdown2);
          end
      end
      
%--------------------------------------------
%计算样本到原形距离矩阵d
      for i=1:c                               
          for k=1:n
              fexp1=1;
              fexp2=1;
              for j=1:datasize
                  fexp1=fexp1*exp(-(data(j,k,1)-v(i,j,1))^2/(2*sigma(i,j,1)));
                  fexp2=fexp2*exp(-(data(j,k,2)-v(i,j,2))^2/(2*sigma(i,j,2)));
              end
              d(i,k,1)=qz(i,1)*fexp1;
              d(i,k,2)=qz(i,2)*fexp2;
          end
      end
      
%--------------------------------------------
%更新模糊划分矩阵u
      uu(:,:,1)=coll*u(:,:,1);                           
      uu(:,:,2)=coll*u(:,:,2);
      for i=1:c          
          for a=1:n
              if sum(d(:,a,1))<=0
                 u(i,a,1)=0;
              else
                 u(i,a,1)=(coll*uu(i,a,1))/(coll+1)+(d(i,a,1)/sum(d(:,a,1)))*(1-sum(uu(:,a,1))/(1+coll)); 
              end
              if sum(d(:,a,2))<=0
                 u(i,a,2)=0;
              else
                 u(i,a,2)=(coll*uu(i,a,2))/(coll+1)+(d(i,a,2)/sum(d(:,a,2)))*(1-sum(uu(:,a,2))/(1+coll)); 
              end
          end
      end
      dis=sum(sum(oldu(:,:,1)-u(:,:,1)).^2)/n*c;
      dis=sqrt(dis);
end

trainy=testing(validationx,w(:,1)',c,th(:,:,1),sigma(:,:,1),v(:,:,1));          %模型对每组训练数据的输出
TRAINRMSE=sqrt(sum((validationy-trainy').^2)/42);                               %每组训练数据的RMSE
RMSESUM=RMSESUM+TRAINRMSE;                                                      %每组训练数据的RMSE和
end

%--------------------------------------------
%训练及测试RMSE
trainRMSE=RMSESUM/L
testy=testing(xtest,w(:,1)',c,th(:,:,1),sigma(:,:,1),v(:,:,1));                 %模型对测试数据的输出,参数采用第一组特征子集的参数
testRMSE=sqrt((sum((ytest-testy').^2))/ntest)                                   %测试数据的RMSE

%--------------------------------------------
%测试的拟合图
put=1:ntest;
plot(put,ytest,'-',put,testy',':')
ylabel('MPG')
xlabel('样本')
h=legend('实际输出','模型输出',2);