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