%本程序实现了把协同模糊聚类算法和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);